text:
#ifndef FUSED_LASSO_B_H
#define FUSED_LASSO_B_H

#include "ALM_Katyusha.h"
#include <string>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <stdio.h>      /* printf */
#include <time.h>
#include <fstream>
#include <algorithm>
#include <iomanip>
#include <ctime>
#include <math.h>

//This class solves problem: min_{x} sum_i 0.5*(Ax- b)^2+ sum_i |x^i- x^{i+1}|+ P(x) by IPALM_Katyusha
// where P(x)=frac{sig2}{2}|x|_2^2 +sig1|x|_1.

// f^i(x)= 0.5*(x- b_i)^2, h^i(x)= |x|, Mx= (x^1- x^2, ..., x_n- x_1)

template<typename L, typename D>
class Fused_lasso_b: public ALM_Katyusha<L, D>
{
private:

	D lambda3;

	Matrix<L,D> my_M;

	Matrix<L,D> my_A;

	D val_lambda_f;

protected:

public:

	D sig1;
	D sig2;

	Fused_lasso_b(const char* Matrix_file, D val_lambda1, D val_lambda2, D val_lambda3)
	:ALM_Katyusha<L,D>(),my_A(Matrix_file)
	{
		my_M.construct_fused_matrix(my_A);
		this->matrix_merge(my_M,my_A);
		sig1=val_lambda1;
		sig2=val_lambda2;
		lambda3=val_lambda3;
		this->L_h= val_lambda3;
		this->mu_g= sig2;
	}

	L get_n(){return my_A.nfeatures;}
	L get_m(){return my_A.nfeatures+ my_A.nsamples;}

	inline void set_matrix_A(){
		this->data_A= my_A;
	}

	inline void set_matrix_M(){
		this->data_M= my_M;
	}

	inline D value_of_f_j(D x, L i){
		return 0.5*(x- my_A.b[i])*(x- my_A.b[i]);
	}

	inline D gradient_of_f_j(D x, L i){
		return x- my_A.b[i];
	}

	inline D value_of_f_star_j(D x, L i){
		return 0.5*x*x+ my_A.b[i]*x;
	}

	inline void rescale(){
		this->lambda1= sig1;
		this->lambda2= sig2+ this->tau_s*this->beta_s;
	}

	inline D value_of_P_j(D x, L j){
		return sig1*fabs(x)+ sig2/2*x*x;
	}

	inline D prox_of_P_j(D x1, D x2, L j){
		if (x1*x2- sig1> 0){
			return (x1*x2- sig1)/(sig2+ x2);
		}
		else if(x1*x2+ sig1< 0){
			return (x1*x2+ sig1)/(sig2+ x2);
		}
		else{
			return 0;
		}
	}

	inline D value_of_h_j(D x, L i){
		return lambda3*fabs(x);
	}

	inline D value_of_h_star_j(D x, L j){
		if (x<= 2*lambda3 && x>= -2*lambda3){
			return 0;
		}
		else{
			cout<< "error in h*(x)"<< endl;
			cout<< "lambda3= "<< lambda3<< "; x= "<< x<< endl;
			system("pause");
			return std::numeric_limits<double>::max();
		}
	}

	inline D compute_one_step(D tau, D u, D x){
		D new_x;
		if(x>tau*(this->lambda1+u))
		new_x=(x-tau*(this->lambda1+u))/(1+this->lambda2*tau);
		else if(x<tau*(u-this->lambda1))
		new_x=(x-tau*(u-this->lambda1))/(1+this->lambda2*tau);
		else
		new_x=0;
		return new_x;
	}

	inline D prox_of_h_star_j(D x1, D x2, L j){
		if (x1> lambda3){
			return lambda3;
		}
		else if (x1< -lambda3){
			return -lambda3;
		}
		else{
			return x1;
		}
	}

	inline D feasible_dual(vector<D> & x){

		if(this->lambda2>0)
		{
			return 1;
		}
		else
		{
			D scal=1;
			L l=x.size();
			for(L i=0;i<l;i++)
			if(fabs(x[i])>this->lambda1)
			scal=min(this->lambda1/fabs(x[i]),scal);
			return scal;
		}
	}

	//We implement Section 6.2 of SPDC paper for 'just-in-time' update
	//t0 is t0+1 in the paper SPDC ; x is the z vector in Katyusha
	inline void compute_just_in_time_prox_grad_without_x_s(D tau, D u, D &x, L t0,L t1, D w, D &y, L j){
		if(t0==t1)
		return;
		else{
			if(this->lambda2>0)
			{
				D theta=1./(1+this->lambda2*tau);
				D hplus=(u+this->lambda1)/this->lambda2;
				D hminus=(u-this->lambda1)/this->lambda2;
				D p=pow(theta,t1-t0);
				if(this->lambda1==0){
					y=compute_aid_2(p,theta, hplus, t1-t0, x, y, w);
					x=p*(x)-(1-p)*u/this->lambda2;
				}else{
					if(x==0){
						if(this->lambda1+u<0){
							y=compute_aid_2(p,theta, hplus, t1-t0, x, y, w);
							x=p*x-(1-p)*hplus;
						}
						else if(u-this->lambda1>0){
							y=compute_aid_2(p,theta, hminus, t1-t0, x, y, w);
							x=p*x-(1-p)*hminus;
						}else{
							x=0;
							D p2=pow(this->theta3,t1-t0);
							D tmp2=(1-p2)/(this->theta1+this->theta2);
							y=(this->theta2*w)*tmp2+p2*y;
						}
					}else if(x>0){
						if(this->lambda1+u<=0){
							y=compute_aid_2(p,theta, hplus, t1-t0, x, y, w);
							x=p*x-(1-p)*hplus;
						}else{
							D t=t0-log(1+this->lambda2*x/(this->lambda1+u))/log(theta);
							if(t<t1){
								L t_tmp=floor(t);
								p=pow(theta,t_tmp-t0);
								y=compute_aid_2(p,theta, hplus, t_tmp-t0, x, y, w);
								x=p*x-(1-p)*hplus;
								x=compute_one_step(tau,u,x);
								y=this->theta1*x+this->theta2*w+this->theta3*y;
								compute_just_in_time_prox_grad_without_x_s(tau, u,  x,  t_tmp+1, t1, w, y, j);
							}
							else{
								y=compute_aid_2(p,theta, hplus, t1-t0, x, y, w);
								x=p*x-(1-p)*hplus;
							}
						}

					}else if(x<0){
						D minusw=-w;
						D minusx=-x;
						D minusy=-y;
						compute_just_in_time_prox_grad_without_x_s(tau, -u,  minusx,  t0, t1,minusw, minusy, j);
						x=-minusx;
						y=-minusy;
					}

				}

			}
			else{

				if(this->lambda1==0){
					y=compute_aid(tau*u, t1-t0, x, y, w);
					x=x-(t1-t0)*tau*u;
				}
				else{
					if(x==0){
						if(this->lambda1+u<0){
							y=compute_aid(tau*(this->lambda1+u), t1-t0, x, y, w);
							x=x-(t1-t0)*tau*(this->lambda1+u);
						}
						else if(u-this->lambda1>0){
							y=compute_aid(tau*(-this->lambda1+u), t1-t0, x, y, w);
							x=x-(t1-t0)*tau*(u-this->lambda1);
						}else{
							D p2=pow(this->theta3,t1-t0);
							D tmp2=(1-p2)/(this->theta1+this->theta2);
							y=(this->theta2*w)*tmp2+p2*y;
							x=0;
						}
					}else if(x>0){
						if(this->lambda1+u<=0){
							y=compute_aid(tau*(this->lambda1+u), t1-t0, x, y, w);
							x=x-(t1-t0)*tau*(this->lambda1+u);
						}else{
							D t=t0+x/(this->lambda1+u);
							if(t<t1){
								L t_tmp=floor(t);
								y=compute_aid(tau*(this->lambda1+u), t_tmp-t0, x, y, w);
								x=x-(t_tmp-t0)*tau*(this->lambda1+u);
								x=compute_one_step(tau,u,x);
								y=this->theta1*x+this->theta2*w+this->theta3*y;
								compute_just_in_time_prox_grad_without_x_s(tau, u,  x,  t_tmp+1, t1, w, y, j);
							}
							else{
								y=compute_aid(tau*(this->lambda1+u), t1-t0, x, y, w);
								x=x-(t1-t0)*tau*(this->lambda1+u);
							}
						}

					}else if(x<0){
						D minusw=-w;
						D minusx=-x;
						D minusy=-y;
						compute_just_in_time_prox_grad_without_x_s(tau, -u,  minusx,  t0, t1,minusw, minusy, j);
						x=-minusx;
						y=-minusy;
					}

				}
			}
		}

	};

	// compute_aid returns the value of y_{k+s} given by: y_{k+1}=theta1*z_{k+1}+theta2*w+(1-theta1-theta2)*y_k; z_{k+i}=z_k-i*h
	D compute_aid(D h, D s, D x, D y, D w){
		D p2=pow(this->theta3,s);
		D tmp2=(1-p2)/(this->theta1+this->theta2);
		D tmp3=(tmp2*this->theta3-p2*s)/(this->theta1+this->theta2);
		y=this->theta1*(tmp2*(x-s*h)+h*tmp3)+this->theta2*w*tmp2+p2*y;
		return y;
	}

	// compute_aid_2 returns the value of y_{k+s} given by: y_{k+1}=theta1*z_{k+1}+theta2*w+(1-theta1-theta2)*y_k; z_{k+i}=q^i(z_k+h)-h
	//p=q^s;
	D compute_aid_2(D p, D q, D h, D s, D x, D y, D w){
		D theta3overq=this->theta3/q;
		D p2=pow(this->theta3,s);
		D p3=pow(theta3overq,s);
		D tmp=(1-p3)/(1-theta3overq);
		D tmp2=(1-p2)/(this->theta1+this->theta2);
		y=this->theta1*p*(x+h)*tmp+(this->theta2*w-this->theta1*h)*tmp2+p2*y;
		return y;
	}

	void Katyusha_solver(D beta_0, D epsilon_0,  D eta , D rho, vector<D> & x0,vector<D> & y0,L val_tau, L max_nb_outer, L p_N_1, L p_N_2,string filename1, string filename2, D time){
		this->ALM_solve_with_L_Katyusha(beta_0, epsilon_0,  eta, rho, x0,y0, val_tau, max_nb_outer,  p_N_1,  p_N_2,  filename1,  filename2, time);
	}

};

#endif


text:
/*
* -----------------------------------------------------------------
*  In Situ Adaptive Tabulation Library --- isat_lib.h
*  Version: 1.6180
*  Date: Feb 25, 2010
* -----------------------------------------------------------------
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010 by Americo Barbosa da Cunha Junior
*
*  This program is free software: you can redistribute it and/or
*  modify it under the terms of the GNU General Public License as
*  published by the Free Software Foundation, either version 3 of
*  the License, or (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  A copy of the GNU General Public License is available in
*  LICENSE.txt or http://www.gnu.org/licenses/.
* -----------------------------------------------------------------
*  This is the header file for a library
*  with the ISAT algorithm.
* -----------------------------------------------------------------
*/

#ifndef __ISAT_LIB_H__
#define __ISAT_LIB_H__

#include <gsl/gsl_blas.h>
#include <gsl/gsl_math.h>

#include "../include/bst_lib.h"

#undef  RTOL
#define RTOL 1.0e-6

#undef  ATOL
#define ATOL 1.0e-9

/*
*------------------------------------------------------------
*   structure of ISAT workspace
*
*   last update: Jul 13, 2010
*------------------------------------------------------------
*/

typedef struct isat_struct
{
    bst_node *root;         /* binary search tree root      */
    unsigned int lf;        /* # of leaves                  */
    unsigned int nd;        /* # of nodes                   */
    unsigned int add;       /* # of adds                    */
    unsigned int grw;       /* # of grows                   */
    unsigned int rtv;       /* # of retrieves               */
    unsigned int dev;       /* # of discarted evaluations   */
    unsigned int hgt;       /* tree height                  */
    unsigned int max_lf;    /* maximum value of tree leaves */
    clock_t time_add;       /* cpu clock spent by additions */
    clock_t time_grw;       /* cpu clock spent by growths   */
    clock_t time_rtv;       /* cpu clock spent by retrieves */
    clock_t time_dev;       /* cpu clock spent by dev       */

} isat_wrk;
/*------------------------------------------------------------*/

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

isat_wrk *isat_alloc();

void isat_free(void **isat_bl);

int isat_set(isat_wrk *isat);

void isat_statistics(isat_wrk *isat);

int isat_input(unsigned int *max_lf,
               double *etol,
               double *n0);

void isat_eoa_mtrx(gsl_matrix *A,
                   double etol,
                   double n0,
                   gsl_matrix *L);

double isat_lerror(gsl_vector *Rphi,
                   gsl_vector *Rphi0,
		           gsl_matrix *A,
                   gsl_vector *phi,
                   gsl_vector *phi0);

int isat4(isat_wrk *isat,
          void *thrm_data,
          void *cvode_mem,
          double etol,
          double n0,
          double t0,
          double delta_t,
          gsl_vector *phi,
          gsl_matrix *A,
          gsl_matrix *L,
          gsl_vector *Rphi);

#endif /* __ISAT_LIB_H__ */

ext: h

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
/* statistics/gsl_statistics_uint.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Jim Davies, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_STATISTICS_UINT_H__
#define __GSL_STATISTICS_UINT_H__

#include <stddef.h>
#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

GSL_EXPORT double gsl_stats_uint_mean (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_uint_variance (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_uint_sd (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_uint_variance_with_fixed_mean (const unsigned int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_uint_sd_with_fixed_mean (const unsigned int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_uint_absdev (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_uint_skew (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_uint_kurtosis (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT double gsl_stats_uint_lag1_autocorrelation (const unsigned int data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_uint_covariance (const unsigned int data1[], const size_t stride1,const unsigned int data2[], const size_t stride2, const size_t n);

GSL_EXPORT double gsl_stats_uint_variance_m (const unsigned int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_uint_sd_m (const unsigned int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_uint_absdev_m (const unsigned int data[], const size_t stride, const size_t n, const double mean);
GSL_EXPORT double gsl_stats_uint_skew_m_sd (const unsigned int data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_uint_kurtosis_m_sd (const unsigned int data[], const size_t stride, const size_t n, const double mean, const double sd);
GSL_EXPORT double gsl_stats_uint_lag1_autocorrelation_m (const unsigned int data[], const size_t stride, const size_t n, const double mean);

GSL_EXPORT double gsl_stats_uint_covariance_m (const unsigned int data1[], const size_t stride1,const unsigned int data2[], const size_t stride2, const size_t n, const double mean1, const double mean2);

GSL_EXPORT double gsl_stats_uint_pvariance (const unsigned int data1[], const size_t stride1, const size_t n1, const unsigned int data2[], const size_t stride2, const size_t n2);
GSL_EXPORT double gsl_stats_uint_ttest (const unsigned int data1[], const size_t stride1, const size_t n1, const unsigned int data2[], const size_t stride2, const size_t n2);

GSL_EXPORT unsigned int gsl_stats_uint_max (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT unsigned int gsl_stats_uint_min (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_uint_minmax (unsigned int * min, unsigned int * max, const unsigned int data[], const size_t stride, const size_t n);

GSL_EXPORT size_t gsl_stats_uint_max_index (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT size_t gsl_stats_uint_min_index (const unsigned int data[], const size_t stride, const size_t n);
GSL_EXPORT void gsl_stats_uint_minmax_index (size_t * min_index, size_t * max_index, const unsigned int data[], const size_t stride, const size_t n);

GSL_EXPORT double gsl_stats_uint_median_from_sorted_data (const unsigned int sorted_data[], const size_t stride, const size_t n) ;
GSL_EXPORT double gsl_stats_uint_quantile_from_sorted_data (const unsigned int sorted_data[], const size_t stride, const size_t n, const double f) ;

__END_DECLS

#endif /* __GSL_STATISTICS_UINT_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/**
 *
 * @file qwrapper_cpotrf.c
 *
 *  PLASMA core_blas quark wrapper
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Jakub Kurzak
 * @date 2010-11-15
 * @generated c Tue Jan  7 11:44:56 2014
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 **/
void QUARK_CORE_cpotrf(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo, int n, int nb,
                       PLASMA_Complex32_t *A, int lda,
                       PLASMA_sequence *sequence, PLASMA_request *request,
                       int iinfo)
{
    DAG_CORE_POTRF;
    QUARK_Insert_Task(quark, CORE_cpotrf_quark, task_flags,
        sizeof(PLASMA_enum),                &uplo,      VALUE,
        sizeof(int),                        &n,         VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,    A,                 INOUT,
        sizeof(int),                        &lda,       VALUE,
        sizeof(PLASMA_sequence*),           &sequence,  VALUE,
        sizeof(PLASMA_request*),            &request,   VALUE,
        sizeof(int),                        &iinfo,     VALUE,
        0);
}

/***************************************************************************//**
 *
 **/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cpotrf_quark = PCORE_cpotrf_quark
#define CORE_cpotrf_quark PCORE_cpotrf_quark
#endif
void CORE_cpotrf_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int n;
    PLASMA_Complex32_t *A;
    int lda;
    PLASMA_sequence *sequence;
    PLASMA_request *request;
    int iinfo;

    int info;

    quark_unpack_args_7(quark, uplo, n, A, lda, sequence, request, iinfo);
    info = LAPACKE_cpotrf_work(
        LAPACK_COL_MAJOR,
        lapack_const(uplo),
        n, A, lda);
    if (sequence->status == PLASMA_SUCCESS && info != 0)
        plasma_sequence_flush(quark, sequence, request, iinfo+info);
}

ext: c


text:
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>

#include "oct.h"

void
print_octave(const gsl_matrix *m, const char *str)
{
  FILE *fp;
  size_t i, j;
  const size_t N = m->size1;
  const size_t M = m->size2;

  if (str == NULL)
    fp = stdout;
  else
    fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %zu\n", M);

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < M; ++j)
        {
          fprintf(fp,
                  "%10.12e%s",
                  gsl_matrix_get(m, i, j),
                  (j < M - 1) ? " " : "\n");
        }
    }

  if (str != NULL)
    fclose(fp);
}

void
printv_octave(const gsl_vector *v, const char *str)
{
  FILE *fp;
  size_t i;
  const size_t N = v->size;

  fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %u\n", 1);

  for (i = 0; i < N; ++i)
    {
      fprintf(fp, "%10.12e\n", gsl_vector_get(v, i));
    }

  fclose(fp);
}

void
printc_octave(const gsl_matrix_complex *m, const char *str)
{
  FILE *fp;
  size_t i, j;
  const size_t N = m->size1;
  const size_t M = m->size2;

  fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: complex matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %zu\n", M);

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < M; ++j)
        {
          gsl_complex z = gsl_matrix_complex_get(m, i, j);
          fprintf(fp,
                  "(%.12e,%.12e)%s",
                  GSL_REAL(z),
                  GSL_IMAG(z),
                  (j < M - 1) ? " " : "\n");
        }
    }

  fclose(fp);
}

/* print symmetric matrix, using upper triangle */
void
printsym_octave(const gsl_matrix *m, const char *str)
{
  FILE *fp;
  size_t i, j;
  const size_t N = m->size1;
  const size_t M = m->size2;

  fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %zu\n", M);

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < M; ++j)
        {
          double z;

          if (j >= i)
            z = gsl_matrix_get(m, i, j);
          else
            z = gsl_matrix_get(m, j, i);

          fprintf(fp,
                  "%.12e%s",
                  z,
                  (j < M - 1) ? " " : "\n");
        }
    }

  fclose(fp);
}

/* print Hermitian matrix, using upper triangle */
void
printherm_octave(const gsl_matrix_complex *m, const char *str)
{
  FILE *fp;
  size_t i, j;
  const size_t N = m->size1;
  const size_t M = m->size2;

  fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: complex matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %zu\n", M);

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < M; ++j)
        {
          gsl_complex z;

          if (j >= i)
            z = gsl_matrix_complex_get(m, i, j);
          else
            z = gsl_complex_conjugate(gsl_matrix_complex_get(m, j, i));

          fprintf(fp,
                  "(%.12e,%.12e)%s",
                  GSL_REAL(z),
                  GSL_IMAG(z),
                  (j < M - 1) ? " " : "\n");
        }
    }

  fclose(fp);
}

void
printcv_octave(const gsl_vector_complex *v, const char *str)
{
  FILE *fp;
  size_t i;
  const size_t N = v->size;

  fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: complex matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %u\n", 1);

  for (i = 0; i < N; ++i)
    {
      gsl_complex z = gsl_vector_complex_get(v, i);
      fprintf(fp, "(%.12e,%.12e)\n", GSL_REAL(z), GSL_IMAG(z));
    }

  fclose(fp);
}

/* print triangular matrix, using upper triangle */
void
printtri_octave(const gsl_matrix *m, const char *str)
{
  FILE *fp;
  size_t i, j;
  const size_t N = m->size1;
  const size_t M = m->size2;

  if (str == NULL)
    fp = stdout;
  else
    fp = fopen(str, "w");

  if (!fp)
    return;

  fprintf(fp, "%% Created by Octave 2.1.73, Tue Aug 01 15:00:27 2006 MDT <blah@blah>\n");
  fprintf(fp, "%% name: %s\n", str);
  fprintf(fp, "%% type: matrix\n");
  fprintf(fp, "%% rows: %zu\n", N);
  fprintf(fp, "%% columns: %zu\n", M);

  for (i = 0; i < N; ++i)
    {
      for (j = 0; j < M; ++j)
        {
          double mij;

          if (j >= i)
            mij = gsl_matrix_get(m, i, j);
          else
            mij = 0.0;

          fprintf(fp,
                  "%10.12e%s",
                  mij,
                  (j < M - 1) ? " " : "\n");
        }
    }

  if (str != NULL)
    fclose(fp);
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#include "../include/paralleltt.h"

#include <stdlib.h>
#include <lapacke.h>
#include <mpi.h>
#include <stdio.h>
#include <string.h>

void SSTT_sketch_to_train(sketch* s, tensor_train* tt, int train_ind)
{
    int head = 0;
    MPI_tensor* ten = s->ten;
    MPI_Comm comm = ten->comm;
    int iscol = s->iscol;
    int flattening = s->flattening;

    flattening_info* fi = flattening_info_init(ten, flattening, iscol, 0);

    // train stuff
    int d = ten->d;
    int n = tt->n[train_ind];
    int r1 = tt->r[train_ind];
    int r2 = tt->r[train_ind + 1];
    double* train = tt->trains[train_ind];
    matrix_tt* train_mat = matrix_tt_wrap(r1 * n, r2, train);
    matrix_tt_fill_zeros(train_mat);

    int* op = s->owner_partition;
    int rank = ten->rank;

    matrix_tt* train_submat = submatrix(train_mat, 0, 0, 0, 0);
    for (int ii = op[rank]; ii < op[rank+1]; ++ii){
        matrix_tt* X_submat = own_submatrix(s, ii, 0);
        matrix_tt* X_subsubmat = submatrix(X_submat, 0, 0, 0, 0);

        flattening_info_f_update(fi, ten, ii);

        int i0_train = fi->f_t_index[0];
        int sz = fi->f_t_sizes[0];

        int jj_max = 1;
        if (train_ind != 0){
            i0_train = i0_train + r1 * fi->f_t_index[1];
            jj_max = fi->f_t_sizes[1];
        }

        for (int jj = 0; jj < jj_max; ++jj){
            submatrix_update(train_submat, i0_train + jj*r1, i0_train + sz + jj*r1, 0, r2);
            submatrix_update(X_subsubmat, jj*sz, (jj+1)*sz, 0, r2);
            matrix_tt_copy_data(train_submat, X_subsubmat);
        }

        free(X_submat);     X_submat = NULL;
        free(X_subsubmat);  X_subsubmat = NULL;
    }
    free(train_submat); train_submat = NULL;

    matrix_tt_reshape(r1*n*r2, 1, train_mat); // Reshape so MPI_Allreduce is only called once
    matrix_tt_allreduce(comm, train_mat);

    flattening_info_free(fi);
    free(train_mat);
    return;
}

MPI_tensor* SSTT_next_ten_init(MPI_tensor* prev_ten, tensor_train* tt, int train_ind)
{
    int r2 = tt->r[train_ind + 1];

    MPI_Comm comm = prev_ten->comm;
    int rank = prev_ten->rank;
    int size = prev_ten->comm_size;
    int is_first = (train_ind == 0);

    int d = (is_first) ? prev_ten->d : prev_ten->d - 1;
    int prev_offset = (is_first) ? 0 : 1;
    int* n = (int*) calloc(d, sizeof(int));
    int* nps = (int*) calloc(d, sizeof(int));
    int Nblocks = 1;
    n[0] = r2;
    nps[0] = 1;

    int** partitions = (int**) malloc(d * sizeof(int*));
    partitions[0] = (int*) calloc(2, sizeof(int));
    partitions[0][1] = r2;
    for (int ii = 1; ii < d; ++ii){
        n[ii] = prev_ten->n[prev_offset + ii];
        nps[ii] = prev_ten->nps[prev_offset + ii];
        Nblocks = Nblocks * nps[ii];

        partitions[ii] = (int*) calloc(nps[ii]+1, sizeof(int));
        int* partition_ii = partitions[ii];
        int* prev_partition_ii = prev_ten->partitions[prev_offset + ii];
        for (int jj = 0; jj < nps[ii]+1; ++jj){
            partition_ii[jj] = prev_partition_ii[jj];
        }
    }

    int** schedule = (int**) calloc(size, sizeof(int*));
    int n_schedule = get_schedule(schedule, nps, d, size);
    int* inverse_schedule = (int*) calloc(Nblocks, sizeof(int));
    for (int ii = 0; ii < size; ++ii){
        int* schedule_ii = schedule[ii];
        for (int jj = 0; jj < n_schedule; ++jj){
            if (schedule_ii[jj] != -1){
                inverse_schedule[schedule_ii[jj]] = ii*n_schedule + jj;
            }
        }
    }

    long* subtensor_sizes = (long*) malloc(n_schedule * sizeof(long));
    int* t_t_block = (int*) malloc(d*sizeof(int));
    long X_size = 0;
    for (int jj = 0; jj < n_schedule; ++jj){

        int t_v_block = schedule[rank][jj];

        if (t_v_block != -1){
            subtensor_sizes[jj] = 1;
            to_tensor_ind(t_t_block, t_v_block, nps, d);

            for (int kk = 0; kk < d; ++kk){
                subtensor_sizes[jj] = subtensor_sizes[jj] * (partitions[kk][t_t_block[kk]+1] - partitions[kk][t_t_block[kk]]);
            }
            X_size = X_size + subtensor_sizes[jj];
        }
        else{
            subtensor_sizes[jj] = 0;
        }
    }
    free(t_t_block);
    double* X = (double*) calloc(X_size, sizeof(double)); // callocing here, because I will probably need it to start at 0

    double* X_subtensor = X;
    double** subtensors = (double**) malloc(n_schedule * sizeof(double*));
    for (int jj = 0; jj < n_schedule; ++jj){
        subtensors[jj] = X_subtensor;
        X_subtensor = X_subtensor + subtensor_sizes[jj];
    }
    int flattening = d;
    void* parameters = p_static_init(subtensors);

    free(subtensor_sizes);

    // Assigning fields
    MPI_tensor* next_ten = (MPI_tensor*) malloc(sizeof(MPI_tensor));

    next_ten->d = d;
    next_ten->n = n;

    next_ten->comm = comm;
    next_ten->rank = rank;
    next_ten->comm_size = size;

    next_ten->schedule = schedule;
    next_ten->n_schedule = n_schedule;
    next_ten->inverse_schedule = inverse_schedule;

    next_ten->partitions = partitions;
    next_ten->nps = nps;

    next_ten->current_part = -1;

    next_ten->f_ten = NULL;
    next_ten->parameters = parameters;

    next_ten->X_size = X_size;
    next_ten->X = X;

    next_ten->ind1 = (int*) malloc(d * sizeof(int));
    next_ten->ind2 = (int*) malloc(d * sizeof(int));
    next_ten->tensor_part = (int*) malloc(d * sizeof(int));
    next_ten->group_ranks = (int*) malloc(size * sizeof(int));
    next_ten->t_kk = (int*) malloc(d * sizeof(int));

    return next_ten;
}

MPI_tensor* SSTT_next_ten(MPI_tensor* prev_ten, tensor_train* tt, int train_ind)
{
    // Get the train matrix. This is what we multiply prev_ten by!
    int n0 = tt->n[train_ind];
    int r1 = tt->r[train_ind];
    int r2 = tt->r[train_ind + 1];
    double* train = tt->trains[train_ind];
    matrix_tt* train_mat = matrix_tt_wrap(r1 * n0, r2, train);

    // Allocate memory for next_ten
    MPI_tensor* next_ten = SSTT_next_ten_init(prev_ten, tt, train_ind);

    // Allocate scratch memory
    long scratch_size = 1;
    for (int ii = 0; ii < next_ten->d; ++ii){
        int* partition_ii = next_ten->partitions[ii];
        int size_ii = 0;
        for (int jj = 0; jj < next_ten->nps[ii]; ++jj){
            int tmp = partition_ii[jj+1] - partition_ii[jj];
            size_ii = (size_ii > tmp) ? size_ii : tmp;
        }
        scratch_size = scratch_size * size_ii;
    }
    double* scratch = (double*) malloc(scratch_size * sizeof(double));

    // Get some indexing stuff
    int iscol = 1;
    int prev_flattening = (train_ind == 0) ? 1 : 2;
    flattening_info* prev_fi = flattening_info_init(prev_ten, prev_flattening, iscol, 0);
    flattening_info* prev_fi_tmp = flattening_info_init(prev_ten, prev_flattening, iscol, 0);

    int next_flattening = 1;
    flattening_info* next_fi = flattening_info_init(next_ten, next_flattening, iscol, 0);
    flattening_info* next_fi_tmp = flattening_info_init(next_ten, next_flattening, iscol, 0);

    // Streaming loop
    int rank = prev_ten->rank;
    int size = prev_ten->comm_size;
    int* prev_schedule_rank = prev_ten->schedule[rank];
    int* next_schedule_rank = next_ten->schedule[rank];

    MPI_Group world_group;
    MPI_Comm_group(next_ten->comm, &world_group);

    matrix_tt* prev_mat = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    matrix_tt* next_mat = (matrix_tt*) calloc(1, sizeof(matrix_tt));
    matrix_tt* share_mat = (matrix_tt*) calloc(1, sizeof(matrix_tt));

    long buf_size = 1;
    for (int ii = 0; ii < next_ten->d; ++ii){
        int sz_ii = 0;
        int* partition_ii = next_ten->partitions[ii];
        for (int jj = 0; jj < next_ten->nps[ii]; ++jj){
            int sz_ii_jj = partition_ii[jj+1] - partition_ii[jj];
            sz_ii = (sz_ii > sz_ii_jj) ? sz_ii : sz_ii_jj;
        }
        buf_size *= sz_ii;
    }
    matrix_tt* buf = matrix_tt_init(buf_size, 1);

    for (int ii = 0; ii < prev_ten->n_schedule; ++ii){
        // Some initialization - figuring out who owns and needs what
        int prev_block_rank = prev_schedule_rank[ii];

        int* next_t_v_blocks = (int*) malloc(size*sizeof(int));
        for (int jj = 0; jj < size; ++jj){
            int prev_block_jj = prev_ten->schedule[jj][ii];

            if (prev_block_jj != -1){
                flattening_info_update(prev_fi_tmp, prev_ten, prev_block_jj);
                int* next_t_t_block = (prev_fi_tmp->t_t_block) + prev_flattening;
                int* next_nps = (next_fi->t_nps) + 1;
                next_t_v_blocks[jj] = to_vec_ind(next_t_t_block, next_nps, next_fi->t_d - 1);
            }
            else{
                next_t_v_blocks[jj] = -1;
            }
        }

        if (prev_block_rank != -1){
            // Get prev_mat
            stream(prev_ten, prev_block_rank);
            flattening_info_update(prev_fi, prev_ten, prev_block_rank);
            matrix_tt_wrap_update(prev_mat, prev_fi->f_N, prev_fi->s_N, get_X(prev_ten));

            // Get next_mat
            int next_t_v_block = next_t_v_blocks[rank];
            flattening_info_update(next_fi, next_ten, next_t_v_block);

            int next_owner;
            int next_epoch;
            MPI_tensor_get_owner(next_ten, next_t_v_block, &next_owner, &next_epoch);

            double beta;

            if (next_owner == rank){
                beta = 1.0;
                stream(next_ten, next_schedule_rank[next_epoch]);
                matrix_tt_wrap_update(next_mat, next_fi->f_N, next_fi->s_N, get_X(next_ten));
            }
            else{
                beta = 0.0;
                matrix_tt_wrap_update(next_mat, next_fi->f_N, next_fi->s_N, scratch);
            }

            // Get train submat
            int i0 = r1 * prev_fi->f_t_index[prev_flattening - 1];
            int i1 = i0 + r1 * prev_fi->f_t_sizes[prev_flattening - 1];
            submatrix_update(train_mat, i0, i1, 0, r2);
            train_mat->transpose = 1;

            // Multiply
            matrix_tt_dgemm(train_mat, prev_mat, next_mat, 1.0, beta);
        }

        // Share
        int* already_shared = (int*) calloc(size, sizeof(int));
        int* group_ranks = prev_ten->group_ranks;
//        printf("Starting share\n");
        for (int jj = 0; jj < size; ++jj){
            int group_size = 0;
            int t_v_block_jj = next_t_v_blocks[jj];
            int owner_jj;
            int epoch_jj;
            MPI_tensor_get_owner(next_ten, t_v_block_jj, &owner_jj, &epoch_jj);
            if (owner_jj == jj){
                already_shared[jj] = 1;
            }
            else if (t_v_block_jj == -1){
                already_shared[jj] = 1;
            }
            else if(already_shared[jj] == 0){
                int owner_ind;
                int rank_ind;
                for (int kk = 0; kk < size; ++kk){
                    if (kk == rank){
                        rank_ind = group_size;
                    }

                    if (kk == owner_jj){
                        group_ranks[group_size] = kk;
                        group_size = group_size + 1;
                    }
                    else if (next_t_v_blocks[kk] == t_v_block_jj){
                        already_shared[kk] = 1;
                        group_ranks[group_size] = kk;
                        group_size = group_size + 1;
                    }
                }

                matrix_tt* reduce_mat = NULL;
                if (owner_jj == rank){
                    flattening_info_update(next_fi_tmp, next_ten, t_v_block_jj);
                    stream(next_ten, t_v_block_jj);
                    matrix_tt_wrap_update(share_mat, next_fi_tmp->t_N, 1, get_X(next_ten));
                    reduce_mat = share_mat;
                }
                else{
                    matrix_tt_reshape(next_mat->n * next_mat->m, 1, next_mat);
                    reduce_mat = next_mat;
                }

                matrix_tt_reshape(reduce_mat->m, reduce_mat->n, buf);
                matrix_tt_group_reduce(next_ten->comm, rank, reduce_mat, buf, owner_jj, group_ranks, group_size);
            }
        }

        free(next_t_v_blocks);
        free(already_shared);
    }
    free(prev_mat);
    free(next_mat);
    free(share_mat);
    matrix_tt_free(buf);

    free(train_mat);
    free(scratch);

    flattening_info_free(prev_fi);
    flattening_info_free(prev_fi_tmp);
    flattening_info_free(next_fi);
    flattening_info_free(next_fi_tmp);

    MPI_Group_free(&world_group);
    return next_ten;
}

void SSTT_copy_final_train(MPI_tensor* ten, tensor_train* tt)
{
    int r1 = tt->r[tt->d - 1];
    int n  = tt->n[tt->d - 1];
    matrix_tt* train_mat = matrix_tt_wrap(r1, n, tt->trains[tt->d - 1]);

    int rank = ten->rank;
    MPI_Comm comm = ten->comm;
    int size = ten->comm_size;

    int flattening = 1;
    int iscol = 0;
    flattening_info* fi = flattening_info_init(ten, flattening, iscol, 0);

    int Nblocks = ten->nps[1];
    int head = 0;
    for (int ii = 0; ii < Nblocks; ++ii){
        int rank_ii;
        int epoch_ii;
        MPI_tensor_get_owner(ten, ii, &rank_ii, &epoch_ii);

        flattening_info_update(fi, ten, ii);
        matrix_tt* train_submat = submatrix(train_mat, 0, r1, fi->f_t_index[0], fi->f_t_index[0] + fi->f_N);

        if (rank_ii == rank){
            stream(ten, ii);
            matrix_tt* ten_mat = matrix_tt_wrap(train_submat->m, train_submat->n, get_X(ten));
            if (rank == head){
                matrix_tt_copy_data(train_submat, ten_mat);
            }
            else{
                matrix_tt_send(comm, ten_mat, head);
            }
            free(ten_mat);
        }
        else if (rank == head){
            matrix_tt_recv(comm, train_submat, rank_ii);
        }

        free(train_submat);
    }

    free(train_mat);
    flattening_info_free(fi);
}

// This frees ten. So just be careful
void SSTT(tensor_train* tt, MPI_tensor* ten){
    int d = ten->d;
    int buf = 2;
    int iscol = 1;

    for (int ii = 0; ii < d-1; ++ii){
        int flattening = (ii == 0) ? 1 : 2;
        sketch* s = sketch_init(ten, flattening, tt->r[ii+1], buf, iscol); // Initialize sketch

        perform_sketch(s); // Sketch

        sketch_qr(s); // QR

        SSTT_sketch_to_train(s, tt, ii); // Get the train

        MPI_tensor* prev_ten = ten;
        ten = SSTT_next_ten(prev_ten, tt, ii); // Perform train^T * ten

        sketch_free(s);
        MPI_tensor_free(prev_ten); prev_ten = NULL;
    }
    SSTT_copy_final_train(ten, tt); // copy ten into trains[d-1]

    MPI_tensor_free(ten);
}


text:
/*
Copyright (C) 2019-2020 JingWeiZhangHuai <jingweizhanghuai@163.com>
Licensed under the Apache License, Version 2.0; you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
*/

// #include <cblas.h>
#include "morn_tensor.h"

struct TensorConnectPara
{
    MLayer *prev;

    int height;
    int width;
    int channel;

    int res_valid;

    float rate;
    float decay;
    float momentum;
};
void *mTensorConnectPara(MList *ini,char *name)
{
    struct TensorConnectPara *para=(struct TensorConnectPara *)mMalloc(sizeof(struct TensorConnectPara));

    para->prev = mNetworkLayer(ini,mINIRead(ini,name,"prev"));
    mException((para->prev == NULL),EXIT,"invalid prev");
    para->res_valid = (strcmp("Input",mLayerType(para->prev))!=0);

    para->height = 1; mINIRead(ini,name,"height" ,"%d",&(para->height ));
    para->width  = 1; mINIRead(ini,name,"width"  ,"%d",&(para->width  ));
    para->channel= 1; mINIRead(ini,name,"channel","%d",&(para->channel));

    para->rate    =0.001;if(mINIRead(ini,name,"rate"    ,"%f",&(para->rate    ))==NULL) mINIRead(ini,"para","rate"    ,"%f",&(para->rate    ));
    para->momentum= 0.9 ;if(mINIRead(ini,name,"momentum","%f",&(para->momentum))==NULL) mINIRead(ini,"para","momentum","%f",&(para->momentum));
    para->decay   = 0.01;if(mINIRead(ini,name,"decay"   ,"%f",&(para->decay   ))==NULL) mINIRead(ini,"para","decay"   ,"%f",&(para->decay   ));
    mException((para->decay<0.0f)||(para->decay>=1.0f),EXIT,"invalid para decay");

    return para;
}

struct HandleTensorConnect
{
    MMemoryBlock *weight;
    MMemoryBlock *update;
};
void endTensorConnect(struct HandleTensorConnect *handle)
{
    if(handle->weight != NULL) mMemoryBlockRelease(handle->weight);
    if(handle->update != NULL) mMemoryBlockRelease(handle->update);
}
#define HASH_TensorConnect 0xb8986c4a

struct HandleTensorConnect *TensorConnectSet(MLayer *layer)
{
    struct TensorConnectPara *para = (struct TensorConnectPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *res= para->prev->res;
    MTensor *out=layer->tns;
    int device = out->device;

    MHandle *hdl=mHandle(out,TensorConnect);
    struct HandleTensorConnect *handle = (struct HandleTensorConnect *)(hdl->handle);
    if(layer->state != DFLT) return handle;

    int weight_height= para->channel*para->height*para->width;
    int weight_width =   in->channel*  in->height*  in->width +1;
    int data_size = weight_height*weight_width;

    mTensorRedefine(out,in->batch,para->channel,para->height,para->width,NULL);
    if(morn_network_flag == MORN_TRAIN)
    {
        mTensorRedefine(res,in->batch,in->channel,in->height,in->width,NULL);

        if(handle->update != NULL) mMemoryBlockRelease(handle->update);
        handle->update =mMemoryBlockCreate(data_size*sizeof(float),device);
        memset(handle->update->data,0,data_size*sizeof(float));
        handle->update->flag = MORN_HOST;
        mMemoryBlockWrite(handle->update);
    }

    if(handle->weight != NULL) mMemoryBlockRelease(handle->weight);
    handle->weight =mMemoryBlockCreate(data_size*sizeof(float),device);
    float *weight_data = (float *)(handle->weight->data);
    if(morn_network_parafile==NULL)
    {
        float scale = sqrt(2.0f/weight_width);
        for(int i=0;i<data_size;i++)
            weight_data[i] = scale*mNormalRand(0.0f,1.0f);//((float)mRand(-16384,16383))/16384.0f;
    }
    else
    {
        mNetworkParaRead(layer,"weight",weight_data,data_size*sizeof(float));
    }
    handle->weight->flag = MORN_HOST;
    mMemoryBlockWrite(handle->weight);

    hdl->valid = 1;
    return handle;
}

void mTensorConnectForward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Connect",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorConnectPara *para = (struct TensorConnectPara *)(layer->para);
    MTensor *in = para->prev->tns;
    MTensor *out= layer->tns;
    int device = out->device;

    struct HandleTensorConnect *handle = TensorConnectSet(layer);

    int weight_height= out->height*out->width*out->channel;
    int weight_width =  in->height* in->width* in->channel +1;
    for(int b=0;b<in->batch;b++)
    {
        mSgemm(device,MORN_NO_TRANS,MORN_NO_TRANS,weight_height,1,weight_width,1.0f,handle->weight,weight_width,mTensorMemory(in,b),1,0.0f,mTensorMemory(out,b),1);
        printf("out_data[8]=%f\n",((float *)(out->data[b]))[8]);
    }

    layer->state = MORN_FORWARD;
}

void mTensorConnectBackward(MLayer *layer)
{
    mException(INVALID_POINTER(layer),EXIT,"invalid input");
    mException(strcmp("Connect",mLayerType(layer)),EXIT,"invalid layer type");
    struct TensorConnectPara *para = (struct TensorConnectPara *)(layer->para);
    MTensor *in  = para->prev->tns;
    MTensor *res = para->prev->res;
    MTensor *out = layer->res;
    int device = out->device;

    MHandle *hdl=mHandle(layer->tns,TensorConnect);
    struct HandleTensorConnect *handle = (struct HandleTensorConnect *)(hdl->handle);
    mException((hdl->valid==0),EXIT,"no forward operate");

    int weight_height= out->height*out->width*out->channel;
    int weight_width =  in->height* in->width* in->channel +1;
    float *update_data = (float*)(handle->update->data);

    mNetworkParaWrite(layer,"weight",handle->weight->data,weight_height*weight_width*sizeof(float));

    float beta = para->momentum;
    for(int b=0;b<in->batch;b++)
    {
        mSgemm(device,MORN_NO_TRANS,MORN_TRANS,weight_height,weight_width,1,1.0f,mTensorMemory(out,b),1,mTensorMemory(in,b),1,beta,handle->update,weight_width);
        printf("update_data[8]=%f\n",((float*)(handle->update->data))[8]);
        beta = 1.0;
    }

    if(para->res_valid) for(int b=0;b<in->batch;b++)
    {
        mSgemm(device,MORN_TRANS,MORN_NO_TRANS,weight_width,1,weight_height,1.0f,handle->weight,weight_width,mTensorMemory(out,b),1,((para->prev->state==MORN_FORWARD)?0.0f:1.0f),mTensorMemory(res,b),1);
        printf("res_data[8]=%f\n",((float *)(res->data[b]))[8]);
    }
    mSaxpby(device,weight_height*weight_width,(0.0f-(para->rate/(float)(in->batch))),handle->update,1,(1.0f-(para->decay*para->rate)),handle->weight,1);

    para->prev->state = MORN_BACKWARD;
}

max_stars_repo_name: ishine/Morn

max_issues_repo_name: ishine/Morn

max_forks_repo_name: ishine/Morn


text:
/*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*
**                                                                                **
** This file forms part of the Underworld geophysics modelling application.         **
**                                                                                **
** For full license and copyright information, please refer to the LICENSE.md file  **
** located at the project root, or contact the authors.                             **
**                                                                                **
**~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*~*/
/*
This provides a couple of routines which allowed the time per iteration
to be stored. This is achieved useing the KSPMonitors. The issue with using
these monitors is the petsc support for these is terrible. We need to be
able to access our user monitor context, but petsc only lets you get the
context from the FIRST monitor, which we cannot guarentee will be the one
defined to store the timing information.

To get around this we store the index into the monitor array which our
monitor gets inersted into. This works, but we have no 'nice' way to
recover the user monitor context except via the monitor_index. Thus when
the monitor is set, we return the value of the monitor index. This number
is required whenever you wish to extract the timing information.

*/
#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif /* _GNU_SOURCE needed for asprintf */
#include <stdio.h>
#include <stdlib.h>

#include <petsc.h>
#include <petscvec.h>
#include <petscmat.h>
#include <petscksp.h>

#include <StGermain/StGermain.h>
#include <StgDomain/StgDomain.h>

#include "common-driver-utils.h"

#include <petscversion.h>
#if ( (PETSC_VERSION_MAJOR >= 3) && (PETSC_VERSION_MINOR >=3) )
  #if (PETSC_VERSION_MINOR >=6)
     #include "petsc/private/kspimpl.h"
  #else
     #include "petsc-private/kspimpl.h"
  #endif
#else
  #include "private/kspimpl.h"
#endif

struct timed_cvg_ctx {
	PetscInt monitor_index;
	PetscInt hist_len, hist_len_max;
	PetscLogDouble *time_log;
	PetscLogDouble time;
	PetscReal *r_log;
};

/* ~~~~~~~~~~~~~~~~ SolverMonitor routines ~~~~~~~~~~~~~~~~~~~~~~ */
/*
I shouldn't need to include the residual log, but the petsc one does strange things, when used
with gmres and fgmres. I know it's related to restarts but I just don't get it. This one will
produce the same data as that which is displayed on the screen when you used -ksp_monitor.
*/

#undef __FUNCT__
#define __FUNCT__ "BSSCR_check_validity_of_timing_monitor"
PetscErrorCode BSSCR_check_validity_of_timing_monitor( KSP ksp, PetscInt monitor_index )
{

	if( monitor_index < 0 ) {
		Stg_SETERRQ1( PETSC_ERR_ARG_WRONG, "Monitor index cannot be negative. You had %D", monitor_index );
	}
	if( monitor_index >= ksp->numbermonitors ) {
		Stg_SETERRQ2( PETSC_ERR_ARG_WRONG, "Monitor index >= number of monitors set (%D). You had %D", ksp->numbermonitors, monitor_index );
	}

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogHistory"
PetscErrorCode BSSCR_KSPLogHistory(struct timed_cvg_ctx *ctx,PetscReal norm)
{
	PetscLogDouble t, time;

	PetscGetTime(&t);
	time = t - ctx->time;
	if (ctx->hist_len==0) {
		/* Get reference time and zero first result */
		PetscGetTime(&ctx->time);
		time = 0.0;
	}

	if (ctx->r_log && ctx->time_log && ctx->hist_len_max > ctx->hist_len) {
		ctx->r_log[ctx->hist_len] = norm;
		ctx->time_log[ctx->hist_len] = time;

		ctx->hist_len++;
	}

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_timed_cvg_test"
PetscErrorCode BSSCR_timed_cvg_test(KSP ksp, int it, PetscReal rnorm, void *mctx)
{
	struct timed_cvg_ctx *ctx = (struct timed_cvg_ctx*)(mctx);

	BSSCR_KSPLogHistory(ctx,rnorm);

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogDestroyMonitor"
PetscErrorCode BSSCR_KSPLogDestroyMonitor(void *_ctx)
{
	struct timed_cvg_ctx *ctx = (struct timed_cvg_ctx*)_ctx;

	PetscFree( ctx->r_log );
	PetscFree( ctx->time_log );
	PetscFree( ctx );

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogSetMonitor"
PetscErrorCode BSSCR_KSPLogSetMonitor(KSP ksp,PetscInt len,PetscInt *monitor_index)
{
	PetscLogDouble *time_log;
	PetscReal *r_log;
	struct timed_cvg_ctx *ctx;

	PetscMalloc( sizeof(PetscLogDouble)*len, &time_log );
	PetscMalloc( sizeof(PetscReal)*len, &r_log );
	PetscMalloc( sizeof(struct timed_cvg_ctx), &ctx );

	ctx->hist_len_max  = len;
	ctx->hist_len      = 0;
	ctx->time_log      = time_log;
	ctx->r_log         = r_log;

	ctx->monitor_index = ksp->numbermonitors;
	*monitor_index = ksp->numbermonitors;

	//KSPSetConvergenceTest(ksp,timed_cvg_test,(void*)ctx );
	KSPMonitorSet( ksp, BSSCR_timed_cvg_test, (void*)ctx, BSSCR_KSPLogDestroyMonitor );

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogGetTimeHistory"
PetscErrorCode BSSCR_KSPLogGetTimeHistory(KSP ksp,PetscInt monitor_index, PetscInt *na,PetscLogDouble **log)
{
	struct timed_cvg_ctx *ctx;
	PetscInt its;

	ctx = (struct timed_cvg_ctx*)ksp->monitorcontext[ monitor_index ];

	KSPGetIterationNumber(ksp,&its);
	its = its + 1;
	*na = PetscMin( (PetscInt)(its), (PetscInt)(ctx->hist_len_max) );
	*log = ctx->time_log;

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogGetResidualHistory"
PetscErrorCode BSSCR_KSPLogGetResidualHistory(KSP ksp,PetscInt monitor_index, PetscInt *na,PetscReal **rlog)
{
	struct timed_cvg_ctx *ctx;
	PetscInt its;

	ctx = (struct timed_cvg_ctx*)ksp->monitorcontext[ monitor_index ];
	KSPGetIterationNumber(ksp,&its);
	its = its + 1;
	*na = PetscMin( (PetscInt)(its), (PetscInt)(ctx->hist_len_max) );
	*rlog = ctx->r_log;

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogGetResidualTimeHistory"
PetscErrorCode BSSCR_KSPLogGetResidualTimeHistory(KSP ksp,PetscInt monitor_index, PetscInt *na,PetscReal **rlog,PetscLogDouble **tlog)
{
	struct timed_cvg_ctx *ctx;
	PetscInt its;

	KSPGetIterationNumber(ksp,&its);
	ctx = (struct timed_cvg_ctx*)ksp->monitorcontext[ monitor_index ];
	its = its + 1;
	*na = PetscMin( (PetscInt)(its), (PetscInt)(ctx->hist_len_max) );
	*rlog = ctx->r_log;
	*tlog = ctx->time_log;

	PetscFunctionReturn(0);
}

KSPNormType  BS_NO_PC=  KSP_NORM_UNPRECONDITIONED;
KSPNormType  BS_W_PC=   KSP_NORM_PRECONDITIONED;

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPGetNormType"
PetscErrorCode BSSCR_KSPGetNormType(KSP ksp, KSPNormType *normtype) {
	*normtype = ksp->normtype;
	return(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_KSPLogSolve"
PetscErrorCode BSSCR_KSPLogSolve(PetscViewer v,PetscInt monitor_index, KSP ksp)
{
	Vec b,u;
	PetscInt i;
	KSPType type;
	PetscReal b_nrm, *rlog;
	PC pc;
	PetscLogDouble *log;
	PetscInt nits;
	KSPNormType ntype;
	PCSide side;
	PetscTruth has_unpc_config, ksp_matches;
	const char *list[]            = { "cg"   , "bicg" , "bcgs" , "gmres" , "fgmres", "minres", "lgmres", 0 };
	const PCSide   pc_side[]      = { PC_LEFT, PC_LEFT, PC_LEFT, PC_RIGHT, PC_RIGHT, PC_LEFT , PC_RIGHT, 0 };
	const KSPNormType norm_type[] = { BS_NO_PC  , BS_NO_PC  , BS_NO_PC  , BS_W_PC    , BS_NO_PC   , BS_NO_PC   , BS_W_PC    , 0 };

	KSPGetRhs( ksp, &b );
	KSPGetSolution( ksp, &u );

	KSPGetType( ksp, &type );
	BSSCR_KSPGetNormType( ksp, &ntype );
	KSPGetPreconditionerSide( ksp, &side );
	has_unpc_config = PETSC_FALSE;

	i = 0;
	while (list[i]) {
		PetscStrcmp( "fgmres", type, &ksp_matches );
		if (ksp_matches) {
			has_unpc_config = PETSC_TRUE;
			break;
		}

		PetscStrcmp( list[i],type,&ksp_matches );
		if (ksp_matches) {
			if( pc_side[i] == side && ntype == norm_type[i] ) {
				has_unpc_config = PETSC_TRUE;
				break;
			}
		}
		i++;
	}

	if (has_unpc_config) {
		PetscViewerASCIIPrintf( v, "[%s] Detected options consistent with unpreconditioned norm: using side=%d, norm_type=%d \n", type, side, ntype );
	}

	BSSCR_KSPLogGetResidualTimeHistory(ksp,monitor_index, &nits,&rlog,&log);
	if (!rlog || !log) {
		Stg_SETERRQ( PETSC_ERR_ORDER, "You need to call BSSCR_KSPLogSetMonitor() before BSSCR_KSPLogGetResidualTimeHistory()" );
	}

	if (!has_unpc_config) {
		Vec pc_b;
		VecDuplicate( u, &pc_b );
		KSPGetPC( ksp, &pc );
		PCApply( pc, b, pc_b );
		VecNorm( pc_b, NORM_2, &b_nrm );

		Stg_VecDestroy(&pc_b );

		PetscViewerASCIIPrintf( v, "# it        |Qr|               |Qr0|              |f|             |Qr|/|Qr0|         log10(Qr/Qr0)       |Qr|/|Qf|          log10(Qr/Qf)           time\n" );
	}
	else if (has_unpc_config) {
		VecNorm( b, NORM_2, &b_nrm );
		PetscViewerASCIIPrintf( v, "# it         |r|               |r0|               |f|              |r|/|r0|          log10(r/r0)          |r|/|f|           log10(r/f)            time\n" );
	}

	for( i=0; i<nits; i++ ) {
		PetscViewerASCIIPrintf( v, "%1.4d %14.12e %14.12e %14.12e %14.12e %14.12e %14.12e %14.12e %14.12e\n",
				i, rlog[i],rlog[0],b_nrm, rlog[i]/rlog[0],log10(rlog[i]/rlog[0]), rlog[i]/b_nrm,log10(rlog[i]/b_nrm), log[i] );
	}

	PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "BSSCR_BSSCR_KSPLogSolveSummary"
PetscErrorCode BSSCR_BSSCR_KSPLogSolveSummary(PetscViewer v,PetscInt monitor_index, KSP ksp)
{
	KSPConvergedReason reason;
	char *reas;
	PetscReal *rlog;
	PetscLogDouble *log;
	PetscInt nits;

	KSPGetConvergedReason( ksp, &reason );
	if( reason == KSP_CONVERGED_RTOL )
		asprintf(&reas, "KSP_CONVERGED_RTOL" );
	else if( reason == KSP_CONVERGED_ATOL )
		asprintf(&reas, "KSP_CONVERGED_ATOL" );
	else if( reason == KSP_CONVERGED_ITS )
		asprintf(&reas, "KSP_CONVERGED_ITS" );
	else if( reason == KSP_DIVERGED_ITS )
		asprintf(&reas, "KSP_DIVERGED_ITS" );
	else if( reason == KSP_DIVERGED_DTOL )
		asprintf(&reas, "KSP_DIVERGED_DTOL" );
	else if( reason == KSP_DIVERGED_INDEFINITE_PC )
		asprintf(&reas, "KSP_DIVERGED_INDEFINITE_PC" );
	else {
		Stg_SETERRQ( PETSC_ERR_SUP, "unknown convergence reason detected");
	}

	BSSCR_KSPLogGetResidualTimeHistory(ksp,monitor_index, &nits,&rlog,&log);
	if (!rlog || !log) {
		Stg_SETERRQ( PETSC_ERR_ORDER, "You need to call BSSCR_KSPLogSetMonitor() before BSSCR_KSPLogGetResidualTimeHistory()" );
	}
	PetscViewerASCIIPrintf( v, "\n");
	PetscViewerASCIIPrintf( v, "# ===========================================================\n");
	PetscViewerASCIIPrintf( v, "#  KSP summary \n" );
	PetscViewerASCIIPrintf( v, "#  time (its) [rnorm] reason \n" );
	PetscViewerASCIIPrintf( v, "#  %.2f (%d) %.2e %s \n", log[nits-1], nits-1, rlog[nits-1], reas );
	PetscViewerASCIIPrintf( v, "# ===========================================================\n");

	PetscFunctionReturn(0);
}


text:
/*
 * BRAINS
 * (B)LR (R)everberation-mapping (A)nalysis (I)n AGNs with (N)ested (S)ampling
 * Yan-Rong Li, liyanrong@ihep.ac.cn
 * Thu, Aug 4, 2016
 */

/*! \file dnest_sa.c
 *  \brief run dnest sampling for sa analysis.
 */

#ifdef SA

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
 #include <stddef.h>
#include <math.h>
#include <float.h>
#include <string.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_interp.h>
#include <mpi.h>

#include "brains.h"

DNestFptrSet *fptrset_sa;

int dnest_sa(int argc, char **argv)
{
  int i;
  double logz_sa;

  set_sa_blr_model();

  num_params_blr = 0;  /* RM parameter numbers */
  num_params_sa_blr = num_params_sa_blr_model + num_params_sa_extpar;
  num_params_blr_tot = num_params_sa_blr;
  num_params_sa = num_params_sa_blr;
  num_params = num_params_sa;

  par_fix = (int *) malloc(num_params * sizeof(int));
  par_fix_val = (double *) malloc(num_params * sizeof(double));
  par_range_model = malloc( num_params * sizeof(double *));
  par_prior_gaussian = malloc( num_params * sizeof(double *));
  for(i=0; i<num_params; i++)
  {
    par_range_model[i] = malloc(2*sizeof(double));
    par_prior_gaussian[i] = malloc(2*sizeof(double));
  }
  par_prior_model = malloc( num_params * sizeof(int));

  fptrset_sa = dnest_malloc_fptrset();
  /* setup functions used for dnest*/
  fptrset_sa->from_prior = from_prior_sa;
  fptrset_sa->print_particle = print_particle_sa;
  fptrset_sa->restart_action = restart_action_sa;
  fptrset_sa->accept_action = accept_action_sa;
  fptrset_sa->kill_action = kill_action_sa;
  fptrset_sa->perturb = perturb_sa;
  fptrset_sa->read_particle = read_particle_sa;

  fptrset_sa->log_likelihoods_cal_initial = log_likelihoods_cal_initial_sa;
  fptrset_sa->log_likelihoods_cal_restart = log_likelihoods_cal_restart_sa;
  fptrset_sa->log_likelihoods_cal = log_likelihoods_cal_sa;

  set_par_range_sa();
  set_par_fix_sa_blrmodel();

  for(i=num_params_sa_blr_model; i<num_params; i++)
  {
    par_fix[i] = 0;
    par_fix_val[i] = -DBL_MAX;
  }

  /* fix DA */
  par_fix[num_params_sa_blr_model] = 1;
  par_fix_val[num_params_sa_blr_model] = log(550.0);

  /* fix FA */
  par_fix[num_params_sa_blr_model+2] = 1;
  par_fix_val[num_params_sa_blr_model+2] = log(1.0);

  print_par_names_sa();

  force_update = parset.flag_force_update;
  if(parset.flag_para_name != 1)
    logz_sa = dnest(argc, argv, fptrset_sa, num_params, dnest_options_file);

  dnest_free_fptrset(fptrset_sa);

  return 0;
}

void set_par_range_sa()
{
  int i;

  /* setup parameter range, BLR parameters first */
  for(i=0; i<num_params_sa_blr_model; i++)
  {
    par_range_model[i][0] = sa_blr_range_model[i][0];
    par_range_model[i][1] = sa_blr_range_model[i][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }

  /* sa extra parameters */
  for(i=num_params_sa_blr_model; i<num_params_sa; i++)
  {
    par_range_model[i][0] = sa_extpar_range[i-num_params_sa_blr_model][0];
    par_range_model[i][1] = sa_extpar_range[i-num_params_sa_blr_model][1];

    par_prior_model[i] = UNIFORM;
    par_prior_gaussian[i][0] = 0.0;
    par_prior_gaussian[i][1] = 0.0;
  }
  return;
}

/*!
 *  print names and prior ranges for parameters
 *
 */
void print_par_names_sa()
{
  if(thistask != roottask)
    return;

  int i, j;
  FILE *fp;
  char fname[BRAINS_MAX_STR_LENGTH], str_fmt[BRAINS_MAX_STR_LENGTH];

  sprintf(fname, "%s/%s", parset.file_dir, "data/para_names_sa.txt");
  fp = fopen(fname, "w");
  if(fp == NULL)
  {
    fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
    exit(0);
  }

  strcpy(str_fmt, "%4d %-15s %10.6f %10.6f %4d %4d %15.6e\n");

  printf("# Print parameter name in %s\n", fname);

  fprintf(fp, "#*************************************************\n");
  fprint_version(fp);
  fprintf(fp, "#*************************************************\n");

  fprintf(fp, "%4s %-15s %10s %10s %4s %4s %15s\n", "#", "Par", "Min", "Max", "Prior", "Fix", "Val");

  i=-1;
  for(j=0; j<num_params_sa_blr_model; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "SA BLR model", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }

  for(j=0; j<num_params_sa_extpar; j++)
  {
    i++;
    fprintf(fp, str_fmt, i, "SA ExtPar", par_range_model[i][0], par_range_model[i][1], par_prior_model[i],
                            par_fix[i], par_fix_val[i]);
  }
  fclose(fp);
  return;
}

/*!
 * this function generates a sample from prior.
 */
void from_prior_sa(void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    if(par_prior_model[i] == GAUSSIAN)
    {
      pm[i] = dnest_randn()*par_prior_gaussian[i][1] + par_prior_gaussian[i][0];
      dnest_wrap(&pm[i], par_range_model[i][0], par_range_model[i][1]);
    }
    else
    {
      pm[i] = par_range_model[i][0] + dnest_rand() * (par_range_model[i][1] - par_range_model[i][0]);
    }
  }

  /* cope with fixed parameters */
  for(i=0; i<num_params; i++)
  {
    if(par_fix[i] == 1)
      pm[i] = par_fix_val[i];
  }

  which_parameter_update = -1;

  return;
}

/*!
 * this function calculates likelihood at initial step.
 */
double log_likelihoods_cal_initial_sa(const void *model)
{
  double logL;
  logL = prob_sa(model);
  return logL;
}

/*!
 * this function calculates likelihood at initial step.
 */
double log_likelihoods_cal_restart_sa(const void *model)
{
  double logL;
  logL = prob_sa(model);
  return logL;
}

/*!
 * this function calculates likelihood.
 */
double log_likelihoods_cal_sa(const void *model)
{
  double logL;
  logL = prob_sa(model);
  return logL;
}

/*!
 * this function prints out parameters.
 */
void print_particle_sa(FILE *fp, const void *model)
{
  int i;
  double *pm = (double *)model;

  for(i=0; i<num_params; i++)
  {
    fprintf(fp, "%e ", pm[i] );
  }
  fprintf(fp, "\n");
  return;
}

/*!
 * This function read the particle from the file.
 */
void read_particle_sa(FILE *fp, void *model)
{
  int j;
  double *psample = (double *)model;

  for(j=0; j < dnest_num_params; j++)
  {
    if(fscanf(fp, "%lf", psample+j) < 1)
    {
      printf("%f\n", *psample);
      fprintf(stderr, "#Error: Cannot read file %s.\n", options.sample_file);
      exit(0);
    }
  }
  return;
}

/*!
 * this function perturbs parameters.
 */
double perturb_sa(void *model)
{
  double *pm = (double *)model;
  double logH = 0.0, limit1, limit2, width;
  int which, which_level, count_saves;

  /*
   * fixed parameters need not to update
   * perturb important parameters more frequently
   */
  do
  {
    which = dnest_rand_int(num_params);
  }while(par_fix[which] == 1);

  which_parameter_update = which;

  /* level-dependent width */
  count_saves = dnest_get_count_saves();
  which_level_update = dnest_get_which_level_update();
  which_level = which_level_update > (size_levels-10)?(size_levels-10):which_level_update;

  if( which_level > 0 && count_saves > 1000)
  {
    limit1 = limits[(which_level-1) * num_params *2 + which *2];
    limit2 = limits[(which_level-1) * num_params *2 + which *2 + 1];
    width = (limit2 - limit1);
    width /= (3*2.35);
  }
  else
  {
    limit1 = par_range_model[which][0];
    limit2 = par_range_model[which][1];
    width = (par_range_model[which][1] - par_range_model[which][0]);
    width /= (2.35);
  }

  if(par_prior_model[which] == GAUSSIAN)
  {
    logH -= (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
    pm[which] += dnest_randh() * width;
    dnest_wrap(&pm[which], par_range_model[which][0], par_range_model[which][1]);
    //dnest_wrap(&pm[which], limit1, limit2);
    logH += (-0.5*pow((pm[which] - par_prior_gaussian[which][0])/par_prior_gaussian[which][1], 2.0) );
  }
  else
  {
    pm[which] += dnest_randh() * width;
    dnest_wrap(&(pm[which]), par_range_model[which][0], par_range_model[which][1]);
    //dnest_wrap(&pm[which], limit1, limit2);
  }

  return logH;
}

/*!
 * this function calculates likelihood.
 */
double log_likelihoods_cal_sa_exam(const void *model)
{
  return 0.0;
}

/*
 * action when a particle's move is accepted.
 * usually store some values with no need to recompute.
 */
void accept_action_sa()
{
  return;
}

/*
 * action when particle i is killed in cdnest sampling.
 * particle i_copy's properties is copyed to particle i.
 */
void kill_action_sa(int i, int i_copy)
{
  return;
}

#endif

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
#pragma once

#include "memcached/types.h"
#include "slabs.h"

#include <gsl/gsl-lite.h>
#include <atomic>
#include <cstddef>
#include <cstring>

/*
 * You should not try to aquire any of the item locks before calling these
 * functions.
 */
typedef struct _hash_item {
    struct _hash_item* next;
    struct _hash_item* prev;
    struct _hash_item* h_next; /* hash chain next */
    /**
     * The unique identifier for this item (it is guaranteed to be unique
     * per key, which means that a two different version of a document
     * cannot have the same CAS value (This is not true after a server
     * restart given that default_bucket is an in-memory bucket).
     */
    uint64_t cas;

    /** least recent access */
    rel_time_t time;

    /** When the item will expire (relative to process startup) */
    rel_time_t exptime;

    /**
     * When the current lock for the object expire. If locktime < "current
     * time" the item isn't locked anymore (timed out). If locktime >=
     * "current time" the object is locked.
     */
    rel_time_t locktime;

    /** The total size of the data (in bytes) */
    uint32_t nbytes;

    /** Flags associated with the item (in network byte order) */
    uint32_t flags;

    /**
     * The number of entities holding a reference to this item object (we
     * operate in a copy'n'write context so it is always safe for all of
     * our clients to share an existing object, but we need the refcount
     * so that we know when we can release the object.
     */
    uint16_t refcount;

    /** Intermal flags used by the engine.*/
    std::atomic<uint8_t> iflag;

    /** which slab class we're in */
    uint8_t slabs_clsid;

    /** to identify the type of the data */
    uint8_t datatype;

    // There is 3 spare bytes due to alignment
} hash_item;

/*
    The structure of the key we hash with.

    This is a combination of the bucket index and the client's key.

    To respect the memcached protocol we support keys > 250, even
    though the current frontend doesn't.

    Keys upto 128 bytes long will be carried wholly on the stack,
    larger keys go on the heap.
*/
typedef struct _hash_key_sized {
    bucket_id_t bucket_index;
    uint8_t client_key[128];
} hash_key_sized;

typedef struct _hash_key_data {
    bucket_id_t bucket_index;
    uint8_t client_key[1];
} hash_key_data;

typedef struct _hash_key_header {
    uint16_t len; /* length of the hash key (bucket_index+client) */
    hash_key_data* full_key; /* points to hash_key::key_storage or a malloc blob*/
} hash_key_header;

typedef struct _hash_key {
    hash_key_header header;
    hash_key_sized key_storage;
} hash_key;

static inline uint8_t* hash_key_get_key(const hash_key* key) {
    return (uint8_t*)key->header.full_key;
}

static inline bucket_id_t hash_key_get_bucket_index(const hash_key* key) {
    return key->header.full_key->bucket_index;
}

static inline void hash_key_set_bucket_index(hash_key* key,
                                             bucket_id_t bucket_index) {
    key->header.full_key->bucket_index = bucket_index;
}

static inline uint16_t hash_key_get_key_len(const hash_key* key) {
    return key->header.len;
}

static inline void hash_key_set_len(hash_key* key, uint16_t len) {
    key->header.len = len;
}

static inline uint8_t* hash_key_get_client_key(const hash_key* key) {
    return key->header.full_key->client_key;
}

static inline uint16_t hash_key_get_client_key_len(const hash_key* key) {
    return hash_key_get_key_len(key) -
           gsl::narrow<uint16_t>(sizeof(key->header.full_key->bucket_index));
}

static inline void hash_key_set_client_key(hash_key* key,
                                           const void* client_key,
                                           const size_t client_key_len) {
    memcpy(key->header.full_key->client_key, client_key, client_key_len);
}

/*
 * return the bytes needed to store the hash_key structure
 * in a single contiguous allocation.
 */
static inline size_t hash_key_get_alloc_size(const hash_key* key) {
    return offsetof(hash_key, key_storage) + hash_key_get_key_len(key);
}

typedef struct {
    unsigned int evicted;
    unsigned int evicted_nonzero;
    rel_time_t evicted_time;
    unsigned int outofmemory;
    unsigned int tailrepairs;
    unsigned int reclaimed;
} itemstats_t;

struct items {
   hash_item *heads[POWER_LARGEST];
   hash_item *tails[POWER_LARGEST];
   itemstats_t itemstats[POWER_LARGEST];
   unsigned int sizes[POWER_LARGEST];
   /*
    * serialise access to the items data
   */
   std::mutex lock;
};

/**
 * Allocate and initialize a new item structure
 * @param engine handle to the storage engine
 * @param key the key for the new item
 * @param nkey the number of bytes in the key
 * @param flags the flags in the new item
 * @param exptime when the object should expire
 * @param nbytes the number of bytes in the body for the item
 * @return a pointer to an item on success NULL otherwise
 */
hash_item *item_alloc(struct default_engine *engine,
                      const void *key, const size_t nkey, int flags,
                      rel_time_t exptime, int nbytes, const void *cookie,
                      uint8_t datatype);

/**
 * Get an item from the cache
 *
 * @param engine handle to the storage engine
 * @param cookie connection cookie
 * @param key the key for the item to get
 * @param nkey the number of bytes in the key
 * @param state Only return documents in this state
 * @return pointer to the item if it exists or NULL otherwise
 */
hash_item* item_get(struct default_engine* engine,
                    const void* cookie,
                    const void* key,
                    const size_t nkey,
                    const DocStateFilter state);

/**
 * Get an item from the cache using a hash_key
 *
 * @param engine handle to the storage engine
 * @param cookie connection cookie
 * @param key to lookup
 * @param state Only return documents in this state
 * @return pointer to the item if it exists or NULL otherwise
 */
hash_item* item_get(struct default_engine* engine,
                    const void* cookie,
                    const hash_key& key,
                    const DocStateFilter state);

/**
 * Get an item from the cache and acquire the lock.
 *
 * @param engine handle to the storage engine
 * @param cookie connection cookie
 * @param where to return the item (if found)
 * @param key the key for the item to get
 * @param nkey the number of bytes in the key
 * @param locktime when the item expire
 * @return ENGINE_SUCCESS for success
 */
ENGINE_ERROR_CODE item_get_locked(struct default_engine* engine,
                                  const void* cookie,
                                  hash_item** it,
                                  const void* key,
                                  const size_t nkey,
                                  rel_time_t locktime);

/**
 * Get and touch an item
 *
 * @param engine handle to the storage engine
 * @param cookie connection cookie
 * @param where to return the item (if found)
 * @param key the key for the item to get
 * @param nkey the number of bytes in the key
 * @param exptime The new expiry time
 * @return ENGINE_SUCCESS for success
 */
ENGINE_ERROR_CODE item_get_and_touch(struct default_engine* engine,
                                     const void* cookie,
                                     hash_item** it,
                                     const void* key,
                                     const size_t nkey,
                                     rel_time_t exptime);

/**
 * Unlock an item in the cache
 *
 * @param engine handle to the storage engine
 * @param cookie connection cookie
 * @param key the key for the item to unlock
 * @param nkey the number of bytes in the key
 * @param cas value for the locked value
 * @return ENGINE_SUCCESS for success
 */
ENGINE_ERROR_CODE item_unlock(struct default_engine* engine,
                              const void* cookie,
                              const void* key,
                              const size_t nkey,
                              uint64_t cas);

/**
 * Reset the item statistics
 * @param engine handle to the storage engine
 */
void item_stats_reset(struct default_engine *engine);

/**
 * Get item statitistics
 * @param engine handle to the storage engine
 * @param add_stat callback provided by the core used to
 *                 push statistics into the response
 * @param cookie cookie provided by the core to identify the client
 */
void item_stats(struct default_engine* engine,
                const AddStatFn& add_stat,
                const void* cookie);

/**
 * Get detaild item statitistics
 * @param engine handle to the storage engine
 * @param add_stat callback provided by the core used to
 *                 push statistics into the response
 * @param cookie cookie provided by the core to identify the client
 */
void item_stats_sizes(struct default_engine* engine,
                      const AddStatFn& add_stat,
                      const void* cookie);

/**
 * Flush expired items from the cache
 * @param engine handle to the storage engine
 */
void  item_flush_expired(struct default_engine *engine);

/**
 * Release our reference to the current item
 * @param engine handle to the storage engine
 * @param it the item to release
 */
void item_release(struct default_engine *engine, hash_item *it);

/**
 * Unlink the item from the hash table (make it inaccessible)
 * @param engine handle to the storage engine
 * @param it the item to unlink
 */
void item_unlink(struct default_engine *engine, hash_item *it);

/**
 * Unlink the item from the hash table (make it inaccessible),
 * but only if the CAS value in the item is the same as the
 * one in the hash table (two different connections may operate
 * on the same objects, so the cas value for the value in the
 * hashtable may be different than the items value. We need
 * to have exclusive access to the hashtable to do the actual
 * unlink)
 *
 * @param engine handle to the storage engine
 * @param it the item to unlink
 */
ENGINE_ERROR_CODE safe_item_unlink(struct default_engine *engine,
                                   hash_item *it);

/**
 * Store an item in the cache
 * @param engine handle to the storage engine
 * @param item the item to store
 * @param cas the cas value (OUT)
 * @param operation what kind of store operation is this (ADD/SET etc)
 * @param document_state the state of the document to store
 * @return ENGINE_SUCCESS on success
 *
 * @todo should we refactor this into hash_item ** and remove the cas
 *       there so that we can get it from the item instead?
 */
ENGINE_ERROR_CODE store_item(struct default_engine *engine,
                             hash_item *item,
                             uint64_t *cas,
                             ENGINE_STORE_OPERATION operation,
                             const void *cookie,
                             const DocumentState document_state);

/**
 * Run a single scrub loop for the engine.
 * @param engine handle to the storage engine
 */
void item_scrubber_main(struct default_engine *engine);

/**
 * Start the item scrubber for the engine
 * @param engine handle to the storage engine
 * @return true if the scrubber has been invoked
 */
bool item_start_scrub(struct default_engine *engine);


text:
/* ///////////////////////////////////////////////////////////////////// */
/*!
  \file
  \brief Contains basic functions for problem initialization.

  The init.c file collects most of the user-supplied functions useful
  for problem configuration.
  It is automatically searched for by the makefile.

  \author A. Mignone (mignone@ph.unito.it)
  \date   March 5, 2017
*/
/* ///////////////////////////////////////////////////////////////////// */
#include "pluto.h"
#include <math.h>
#include "coolit.h"
#include <stdbool.h>
//#include <gsl.h>

double mu, Rcl;
double rho0, P0, cs0, v0, s0, cool_0, tcool0, r_crit;

int q = 2, counter = 0;
double T0 , Mach0, n0, chi, singular;

/* ********************************************************************* */
void Init (double *v, double x1, double x2, double x3)
/*!
 * The Init() function can be used to assign initial conditions as
 * as a function of spatial position.
 *
 * \param [out] v   a pointer to a vector of primitive variables
 * \param [in] x1   coordinate point in the 1st dimension
 * \param [in] x2   coordinate point in the 2nd dimension
 * \param [in] x3   coordinate point in the 3rdt dimension
 *
 * The meaning of x1, x2 and x3 depends on the geometry:
 * \f[ \begin{array}{cccl}
 *    x_1  & x_2    & x_3  & \mathrm{Geometry}    \\ \noalign{\medskip}
 *     \hline
 *    x    &   y    &  z   & \mathrm{Cartesian}   \\ \noalign{\medskip}
 *    R    &   z    &  -   & \mathrm{cylindrical} \\ \noalign{\medskip}
 *    R    & \phi   &  z   & \mathrm{polar}       \\ \noalign{\medskip}
 *    r    & \theta & \phi & \mathrm{spherical}
 *    \end{array}
 *  \f]
 *
 * Variable names are accessed by means of an index v[nv], where
 * nv = RHO is density, nv = PRS is pressure, nv = (VX1, VX2, VX3) are
 * the three components of velocity, and so forth.
 *
 *********************************************************************** */
{

  T0 = g_inputParam[T_CRIT]; //K
  Mach0 = g_inputParam[MACH_CRIT];
  n0 = g_inputParam[N_CRIT]; //cm^-3
  chi = g_inputParam[CHI];
  singular = g_inputParam[SINGULAR];

  mu = MeanMolecularWeight(v);
  rho0 = n0*mu*CONST_mp; //g cm^-3
  P0 = n0*CONST_kB*T0; //dyne cm^-2
  cs0 = sqrt(g_gamma*P0/rho0); //cm s^-1
  v0 = -Mach0*cs0; //cm s^-1
  s0 = P0/pow(rho0,g_gamma); //CGS
  cool_0 = lambda(T0); //CGS
  tcool0 = (1./(g_gamma-1))*CONST_kB*T0/(n0*cool_0); //s

  r_crit = -q*v0*g_gamma*tcool0; //cm
  Rcl = singular*r_crit/UNIT_LENGTH; // code units (pc usually)

  //double steep=20., f=0.007;
  //double shift = steep*atanh(2.*f*Tcl/(Tw-Tcl)-1);
  if (prank==0 && counter==0) printf("T0=%.1e\tM0=%.1f\tn0=%.2f\tchi=%d\tsingular=%.2f\n",T0,Mach0, n0, (int)chi, singular);
  if (prank==0 && counter==0) printf("rho0=%.1e\t\tP0=%.1e\t\tv0=%.2e\t\tr_crit=%.2e pc\tRcl=%.2e pc\t\ttcool0=%.2e\n",rho0/UNIT_DENSITY,P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2)), v0/UNIT_VELOCITY, r_crit/UNIT_LENGTH, Rcl, tcool0/(CONST_PI*1e12));
  if (prank==0 && counter==0) printf("x from %.2f to %.2f\n",g_domBeg[IDIR],g_domEnd[IDIR],x1);
  if (prank==0 && (x1<g_domBeg[IDIR] || x1>g_domEnd[IDIR])) printf("%.2f pc\n",x1);
  counter++;

  double T = (x1/Rcl<1.)?T0/chi:T0; //Tcl + (Tw-Tcl)/2. *(1+ tanh(x1/steep+shift));
  //T = Tcl + (chi-1.)*Tcl*(x1-g_domBeg[IDIR])/(g_domEnd[IDIR]-g_domBeg[IDIR]); //a linear temperature going from cooler to hotter
  v[RHO] = (x1/Rcl<1.)?chi*n0*mu*CONST_mp/UNIT_DENSITY:n0*mu*CONST_mp/UNIT_DENSITY;
  EXPAND(
  v[VX1] = 0.0;,
  v[VX2] = 0.0;,
  v[VX3] = 0.0;);

  //#if HAVE_ENERGY
  v[PRS] = P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
  //#endif
  v[TRC] = 1.0;
  //Pwind = P/(UNIT_DENSITY*UNIT_VELOCITY*UNIT_VELOCITY);
  //Dwind = ncl*mu*mp/chi/UNIT_DENSITY;

  #if PHYSICS == MHD || PHYSICS == RMHD
  v[BX1] = 0.0;
  v[BX2] = 0.0;
  v[BX3] = 0.0;

  v[AX1] = 0.0;
  v[AX2] = 0.0;
  v[AX3] = 0.0;
  #endif
}

/* ********************************************************************* */
void InitDomain (Data *d, Grid *grid)
/*!
 * Assign initial condition by looping over the computational domain.
 * Called after the usual Init() function to assign initial conditions
 * on primitive variables.
 * Value assigned here will overwrite those prescribed during Init().
 *
 *
 *********************************************************************** */
{
}

/* ********************************************************************* */
void Analysis (const Data *d, Grid *grid)
/*!
 *  Perform runtime data analysis.
 *
 * \param [in] d the PLUTO Data structure
 * \param [in] grid   pointer to array of Grid structures
 *
 *********************************************************************** */
{
	/*
	int i, j, k;
	double *dx, *dy, *dz;*/
	/* ---- Parallel data reduction ---- */
	/*
	#ifdef PARALLEL
	int transfer_size = 3;
	int transfer = 0;
	double sendArray [transfer_size], recvArray[transfer_size];*/

	/* ---- Set pointer shortcuts ---- */
	/*
	dx = grid->dx[IDIR];
	dy = grid->dx[JDIR];
	dz = grid->dx[KDIR];*/

}
#if PHYSICS == MHD
/* ********************************************************************* */
void BackgroundField (double x1, double x2, double x3, double *B0)
/*!
 * Define the component of a static, curl-free background
 * magnetic field.
 *
 * \param [in] x1  position in the 1st coordinate direction \f$x_1\f$
 * \param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
 * \param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
 * \param [out] B0 array containing the vector componens of the background
 *                 magnetic field
 *********************************************************************** */
{
   B0[0] = 0.0;
   B0[1] = 0.0;
   B0[2] = 0.0;
}
#endif

/* ********************************************************************* */
void UserDefBoundary (const Data *d, RBox *box, int side, Grid *grid)
/*!
 *  Assign user-defined boundary conditions.
 *
 * \param [in,out] d  pointer to the PLUTO data structure containing
 *                    cell-centered primitive quantities (d->Vc) and
 *                    staggered magnetic fields (d->Vs, when used) to
 *                    be filled.
 * \param [in] box    pointer to a RBox structure containing the lower
 *                    and upper indices of the ghost zone-centers/nodes
 *                    or edges at which data values should be assigned.
 * \param [in] side   specifies the boundary side where ghost zones need
 *                    to be filled. It can assume the following
 *                    pre-definite values: X1_BEG, X1_END,
 *                                         X2_BEG, X2_END,
 *                                         X3_BEG, X3_END.
 *                    The special value side == 0 is used to control
 *                    a region inside the computational domain.
 * \param [in] grid  pointer to an array of Grid structures.
 *
 *********************************************************************** */
{
  int   i, j, k, nv;
  double  *x1, *x2, *x3;

  x1 = grid->x[IDIR];
  x2 = grid->x[JDIR];
  x3 = grid->x[KDIR];

  if (side == X1_BEG) {

    BOX_LOOP(box,k,j,i) {
		d->Vc[VX1][k][j][i] = d->Vc[VX1][k][j][IBEG];   //outflow velocity
		d->Vc[RHO][k][j][i] = 3.*rho0/UNIT_DENSITY;
		d->Vc[PRS][k][j][i] = 1.52*P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
		//d->Vc[VX1][k][j][i] = 0.14*sqrt(g_gamma*d->Vc[PRS][k][j][i]/d->Vc[RHO][k][j][i]);
	}

  }else if (side == X1_END){

    BOX_LOOP(box,k,j,i) {
		d->Vc[VX1][k][j][i] = d->Vc[VX1][k][j][IEND]; //outflow velocity
		d->Vc[RHO][k][j][i] = 0.02*rho0/UNIT_DENSITY;
		d->Vc[PRS][k][j][i] = 1.12*P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
		//d->Vc[VX1][k][j][i] = 0.*sqrt(g_gamma*d->Vc[PRS][k][j][i]/d->Vc[RHO][k][j][i]);
	}

  }

  /*
  if (side == 0) {  */  /* -- check solution inside domain -- */
  /* TOT_LOOP(k,j,i){

    if ((x1[i]*UNIT_LENGTH/r_crit>=(1.-1e-1)) && (x1[i]*UNIT_LENGTH/r_crit<=(1.+1e-1))) {
		d->Vc[RHO][k][j][i] = rho0/UNIT_DENSITY;
		d->Vc[VX1][k][j][i] = v0/UNIT_VELOCITY;
		d->Vc[PRS][k][j][i] = P0/(UNIT_DENSITY*pow(UNIT_VELOCITY,2));
	}/*
	if (x1[i]/Rcl<1.)
		d->Vc[VX1][k][j][i] = 0.; *//*
   }
  }*/
}

#if BODY_FORCE != NO
/* ********************************************************************* */
void BodyForceVector(double *v, double *g, double x1, double x2, double x3)
/*!
 * Prescribe the acceleration vector as a function of the coordinates
 * and the vector of primitive variables *v.
 *
 * \param [in] v  pointer to a cell-centered vector of primitive
 *                variables
 * \param [out] g acceleration vector
 * \param [in] x1  position in the 1st coordinate direction \f$x_1\f$
 * \param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
 * \param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
 *
 *********************************************************************** */
{
  g[IDIR] = 0.0;
  g[JDIR] = 0.0;
  g[KDIR] = 0.0;
}
/* ********************************************************************* */
double BodyForcePotential(double x1, double x2, double x3)
/*!
 * Return the gravitational potential as function of the coordinates.
 *
 * \param [in] x1  position in the 1st coordinate direction \f$x_1\f$
 * \param [in] x2  position in the 2nd coordinate direction \f$x_2\f$
 * \param [in] x3  position in the 3rd coordinate direction \f$x_3\f$
 *
 * \return The body force potential \f$ \Phi(x_1,x_2,x_3) \f$.
 *
 *********************************************************************** */
{
  return 0.0;
}
#endif

max_stars_repo_name: dutta-alankar/cooling-flow-model

max_issues_repo_name: dutta-alankar/cooling-flow-model

max_forks_repo_name: dutta-alankar/cooling-flow-model


text:
/*
 * Copyright (c) 2011-2013  University of Texas at Austin. All rights reserved.
 *
 * $COPYRIGHT$
 *
 * Additional copyrights may follow
 *
 * This file is part of PerfExpert.
 *
 * PerfExpert is free software: you can redistribute it and/or modify it under
 * the terms of the The University of Texas at Austin Research License
 *
 * PerfExpert is distributed in the hope that it will be useful, but WITHOUT ANY
 * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 * A PARTICULAR PURPOSE.
 *
 * Authors: Leonardo Fialho and Ashay Rane
 *
 * $HEADER$
 */

#ifndef HISTOGRAM_H_
#define HISTOGRAM_H_

#include <algorithm>
#include <gsl/gsl_histogram.h>

#include "analysis_defs.h"

bool pair_sort(const pair_t& p1, const pair_t& p2);

int flatten_and_sort_histogram(gsl_histogram*& hist, pair_list_t& pair_list);

int create_histogram_if_null(gsl_histogram*& hist, size_t bins);

#endif  /* HISTOGRAM_H_ */

max_stars_repo_name: roystgnr/perfexpert

max_issues_repo_name: roystgnr/perfexpert

max_forks_repo_name: roystgnr/perfexpert


text:
#include <stdio.h>
#include <gsl/gsl_vector.h>

int
main (void)
{
  int i;
  gsl_vector * v = gsl_vector_alloc (3);

  for (i = 0; i < 3; i++)
    {
      gsl_vector_set (v, i, 1.23 + i);
    }

  for (i = 0; i < 100; i++) /* OUT OF RANGE ERROR */
    {
      printf ("v_%d = %g\n", i, gsl_vector_get (v, i));
    }

  gsl_vector_free (v);
  return 0;
}

max_stars_repo_name: manggoguy/parsec-modified

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
/*
** estimate false discovery rate based on histograms
** input histograms are generated by "vted"
**
** Ref: Lohmann et al (2016) PLoS One, 11(6):e0158185
**
** G.Lohmann, MPI-KYB, Jan 2015
*/
#include <viaio/VImage.h>
#include <viaio/Vlib.h>
#include <viaio/mu.h>
#include <viaio/option.h>

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>

#define SQR(x) ((x)*(x))
#define ABS(x) ((x) > 0 ? (x) : -(x))

#define HMIN 0
#define HMAX 1.001

#define NBINS 10000

size_t CountLines(FILE *fp)
{
  char ch;
  size_t numlines=0;
  while(!feof(fp)) {
    ch = fgetc(fp);
    if(ch == '\n') {
      numlines++;
    }
  }
  rewind(fp);
  return numlines;
}

void ReadTxtHist(char *filename,gsl_histogram *hist)
{
  FILE *fp=fopen(filename,"r");
  if (!fp) VError("err opening %s",filename);

  size_t nbins = CountLines(fp);
  if (nbins != (size_t)NBINS)
    VError(" inconsistent bins: %lu",gsl_histogram_bins (hist));

  size_t i=0;
  double lower=0,hz=0,mx=0;
  for (i=0; i<nbins; i++) {
    if (fscanf(fp,"%lf %lf %lf",&lower,&hz,&mx) != 3) VError(" read error");
    if (lower < HMIN || lower >= HMAX) VError(" value out of histogram range: %f",lower);
    hist->bin[i] = mx;
  }
  fclose(fp);
}

int main(int argc, char *argv[])
{
  static VString realfilename = "";
  static VString nullfilename = "";
  static VString outfilename = "";
  static VFloat  alpha = 0.05;
  static VOptionDescRec  options[] = {
    {"real",VStringRepn,1,(VPointer) &realfilename,VRequiredOpt,NULL,"input real histogram file"},
    {"null",VStringRepn,1,(VPointer) &nullfilename,VRequiredOpt,NULL,"input null histogram file"},
    {"out",VStringRepn,1,(VPointer) &outfilename,VRequiredOpt,NULL,"output txt file"},
    {"alpha",VFloatRepn,1,(VPointer) &alpha,VOptionalOpt,NULL,"alpha level"},
  };
  double lower=0,upper=0,tiny=1.0e-8;
  long i;
  VParseFilterCmd(VNumber(options),options,argc,argv,NULL,NULL);
  gsl_set_error_handler_off();

  /* ini histogram structs */
  size_t nbins = NBINS;
  double hmin = HMIN,hmax = HMAX;
  gsl_histogram *realhist = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (realhist,hmin,hmax);
  gsl_histogram *nullhist = gsl_histogram_alloc (nbins);
  gsl_histogram_set_ranges_uniform (nullhist,hmin,hmax);

  /* read input histograms */
  ReadTxtHist(realfilename,realhist);
  ReadTxtHist(nullfilename,nullhist);

  /* cumulative distribution functions */
  gsl_histogram_pdf *cdfz = gsl_histogram_pdf_alloc(nbins);
  gsl_histogram_pdf *cdf0 = gsl_histogram_pdf_alloc(nbins);
  gsl_histogram_pdf_init (cdfz,realhist);
  gsl_histogram_pdf_init (cdf0,nullhist);

  FILE *fp = fopen(outfilename,"w");
  if (!fp) VError(" err opening file %s",outfilename);
  fprintf(fp,"#     ted            F0            Fz           FDR \n");
  fprintf(fp,"#------------------------------------------------------\n");

  /* tail-area FDR */
  double cutoff=0;
  int flag = 0;
  for (i=0; i<nbins; i++) {
    gsl_histogram_get_range (realhist,i,&lower,&upper);
    double z = lower + 0.5*(upper-lower);
    double Fz = 1.0-cdfz->sum[i];
    double F0 = 1.0-cdf0->sum[i];
    double Fdr = 0.0;
    if (Fz  > tiny) Fdr = F0/Fz;
    if (Fdr > 1.0) Fdr = 1.0;
    if (Fdr < 0.0) Fdr = 0.0;

    if (Fdr < alpha && flag == 0 && i > 10) {
      flag = 1;
      cutoff = z;
    }
    fprintf(fp," %12.8f  %12.8f  %12.8f  %12.8f\n",z,F0,Fz,Fdr);
  }
  fprintf(fp,"# Fdr < %.3f  if ted > %lf\n",alpha,cutoff);
  fprintf(stderr," Fdr < %.3f  if ted > %lf\n",alpha,cutoff);

  fclose(fp);
  exit(0);
}

max_stars_repo_name: zrajna/lipsia

max_issues_repo_name: zrajna/lipsia

max_forks_repo_name: zrajna/lipsia


text:
/* vector/gsl_vector_ulong.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_VECTOR_ULONG_H__
#define __GSL_VECTOR_ULONG_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_inline.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_block_ulong.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  unsigned long *data;
  gsl_block_ulong *block;
  int owner;
}
gsl_vector_ulong;

typedef struct
{
  gsl_vector_ulong vector;
} _gsl_vector_ulong_view;

typedef _gsl_vector_ulong_view gsl_vector_ulong_view;

typedef struct
{
  gsl_vector_ulong vector;
} _gsl_vector_ulong_const_view;

typedef const _gsl_vector_ulong_const_view gsl_vector_ulong_const_view;

/* Allocation */

GSL_FUN gsl_vector_ulong *gsl_vector_ulong_alloc (const size_t n);
GSL_FUN gsl_vector_ulong *gsl_vector_ulong_calloc (const size_t n);

GSL_FUN gsl_vector_ulong *gsl_vector_ulong_alloc_from_block (gsl_block_ulong * b,
                                                     const size_t offset,
                                                     const size_t n,
                                                     const size_t stride);

GSL_FUN gsl_vector_ulong *gsl_vector_ulong_alloc_from_vector (gsl_vector_ulong * v,
                                                      const size_t offset,
                                                      const size_t n,
                                                      const size_t stride);

GSL_FUN void gsl_vector_ulong_free (gsl_vector_ulong * v);

/* Views */

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_view_array (unsigned long *v, size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_view_array_with_stride (unsigned long *base,
                                         size_t stride,
                                         size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_view_array (const unsigned long *v, size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_view_array_with_stride (const unsigned long *base,
                                               size_t stride,
                                               size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_subvector (gsl_vector_ulong *v,
                            size_t i,
                            size_t n);

GSL_FUN _gsl_vector_ulong_view
gsl_vector_ulong_subvector_with_stride (gsl_vector_ulong *v,
                                        size_t i,
                                        size_t stride,
                                        size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_subvector (const gsl_vector_ulong *v,
                                  size_t i,
                                  size_t n);

GSL_FUN _gsl_vector_ulong_const_view
gsl_vector_ulong_const_subvector_with_stride (const gsl_vector_ulong *v,
                                              size_t i,
                                              size_t stride,
                                              size_t n);

/* Operations */

GSL_FUN void gsl_vector_ulong_set_zero (gsl_vector_ulong * v);
GSL_FUN void gsl_vector_ulong_set_all (gsl_vector_ulong * v, unsigned long x);
GSL_FUN int gsl_vector_ulong_set_basis (gsl_vector_ulong * v, size_t i);

GSL_FUN int gsl_vector_ulong_fread (FILE * stream, gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_fwrite (FILE * stream, const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_fscanf (FILE * stream, gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_fprintf (FILE * stream, const gsl_vector_ulong * v,
                              const char *format);

GSL_FUN int gsl_vector_ulong_memcpy (gsl_vector_ulong * dest, const gsl_vector_ulong * src);

GSL_FUN int gsl_vector_ulong_reverse (gsl_vector_ulong * v);

GSL_FUN int gsl_vector_ulong_swap (gsl_vector_ulong * v, gsl_vector_ulong * w);
GSL_FUN int gsl_vector_ulong_swap_elements (gsl_vector_ulong * v, const size_t i, const size_t j);

GSL_FUN unsigned long gsl_vector_ulong_max (const gsl_vector_ulong * v);
GSL_FUN unsigned long gsl_vector_ulong_min (const gsl_vector_ulong * v);
GSL_FUN void gsl_vector_ulong_minmax (const gsl_vector_ulong * v, unsigned long * min_out, unsigned long * max_out);

GSL_FUN size_t gsl_vector_ulong_max_index (const gsl_vector_ulong * v);
GSL_FUN size_t gsl_vector_ulong_min_index (const gsl_vector_ulong * v);
GSL_FUN void gsl_vector_ulong_minmax_index (const gsl_vector_ulong * v, size_t * imin, size_t * imax);

GSL_FUN int gsl_vector_ulong_add (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_sub (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_mul (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_div (gsl_vector_ulong * a, const gsl_vector_ulong * b);
GSL_FUN int gsl_vector_ulong_scale (gsl_vector_ulong * a, const unsigned long x);
GSL_FUN int gsl_vector_ulong_add_constant (gsl_vector_ulong * a, const double x);
GSL_FUN int gsl_vector_ulong_axpby (const unsigned long alpha, const gsl_vector_ulong * x, const unsigned long beta, gsl_vector_ulong * y);
GSL_FUN unsigned long gsl_vector_ulong_sum (const gsl_vector_ulong * a);

GSL_FUN int gsl_vector_ulong_equal (const gsl_vector_ulong * u,
                            const gsl_vector_ulong * v);

GSL_FUN int gsl_vector_ulong_isnull (const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_ispos (const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_isneg (const gsl_vector_ulong * v);
GSL_FUN int gsl_vector_ulong_isnonneg (const gsl_vector_ulong * v);

GSL_FUN INLINE_DECL unsigned long gsl_vector_ulong_get (const gsl_vector_ulong * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_ulong_set (gsl_vector_ulong * v, const size_t i, unsigned long x);
GSL_FUN INLINE_DECL unsigned long * gsl_vector_ulong_ptr (gsl_vector_ulong * v, const size_t i);
GSL_FUN INLINE_DECL const unsigned long * gsl_vector_ulong_const_ptr (const gsl_vector_ulong * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
unsigned long
gsl_vector_ulong_get (const gsl_vector_ulong * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, 0);
    }
#endif
  return v->data[i * v->stride];
}

INLINE_FUN
void
gsl_vector_ulong_set (gsl_vector_ulong * v, const size_t i, unsigned long x)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  v->data[i * v->stride] = x;
}

INLINE_FUN
unsigned long *
gsl_vector_ulong_ptr (gsl_vector_ulong * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (unsigned long *) (v->data + i * v->stride);
}

INLINE_FUN
const unsigned long *
gsl_vector_ulong_const_ptr (const gsl_vector_ulong * v, const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return (const unsigned long *) (v->data + i * v->stride);
}
#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_ULONG_H__ */

max_stars_repo_name: mgreter/astrometrylib

max_issues_repo_name: zzpwahaha/VimbaCamJILA

max_forks_repo_name: zzpwahaha/VimbaCamJILA