text:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

#define DECORRCRITS 1
#define SET 1

int pcl_ran_mvgaussian(gsl_rng *rseed, gsl_vector *mean, gsl_matrix *Sigma, gsl_vector *output)
{
  int i;
  gsl_linalg_cholesky_decomp(Sigma);
  for (i = 0; i < output->size; i++){
    gsl_vector_set (output, i , gsl_ran_gaussian(rseed,1));
  }
  gsl_blas_dtrmv (CblasLower, CblasNoTrans, CblasNonUnit, Sigma, output);
  gsl_blas_daxpy (1.0,mean, output);
  return 0;
}

unsigned long int random_seed()
{

  unsigned int seed;
  FILE *devrandom;
  if ((devrandom = fopen("/dev/random","r")) == NULL) {
    fprintf(stderr,"Cannot open /dev/random, setting seed to 0\n");
    seed = 0;
  } else {
    fread(&seed,sizeof(seed),1,devrandom);
    fclose(devrandom);
  }

  return(seed);

}

double ran_inv_gamma(const gsl_rng * r, double a, double b){
  return(1/gsl_ran_gamma(r, a, 1/b));
}

double ran_trunc_norm_upper(const gsl_rng * r, double mu, double sigma, double a){
  double u = gsl_rng_uniform_pos(r);
  double v = gsl_cdf_ugaussian_P((-a+mu)/sigma);
  return(-sigma * gsl_cdf_ugaussian_Pinv(u*v)+mu);
}

double ran_trunc_norm_lower(const gsl_rng * r, double mu, double sigma, double a){
  double u = gsl_rng_uniform_pos(r);
  double v = gsl_cdf_ugaussian_P((a-mu)/sigma);
  return(sigma * gsl_cdf_ugaussian_Pinv(u*v)+mu);
}

double ran_trunc_norm_both(const gsl_rng * r, double mu, double sigma, double a, double b){
  double low = gsl_cdf_ugaussian_P((a-mu)/sigma);
  double high = gsl_cdf_ugaussian_P((b-mu)/sigma);
  double u = gsl_ran_flat(r,low,high);
  return(sigma * gsl_cdf_ugaussian_Pinv(u)+mu);
}

int main(int argc, char *argv[]){

  unsigned long int seed;
  seed = random_seed();
  int ITERS,iter,i,j,k,K,I,J,nRespCat;
  char *fnEstimates,*fnChains,*fnInput;
  int yij,sij;
  int errs=0;
  char *name;

  double a,b,e,f,sigDecorr,mu_crits,sd_crits;
  double sig2mu,mumu;

  ////This is all data and prior entry code.

  if(argc<9){
    fprintf(stderr,"Arguments needed: #iterations nRespCat sig2_mu a b e f sigDecorr [analysis.name] [input.file]\n");
    fprintf(stderr,"\n****Not enough arguments.\n\n");
    return(1);
  }
  if(argc>9){
    name=argv[9];
  }else{
    name="analysis";
  }

  ITERS=atoi(argv[1]);
  nRespCat=atoi(argv[2]);
  sig2mu=atof(argv[3]);
  a=atof(argv[4]);
  b=atof(argv[5]);
  e=atof(argv[6]);
  f=atof(argv[7]);
  sigDecorr=atof(argv[8]);

  //confirm and check command line arguments.
  printf("Name      : %s\n",name);
  printf("nRespCat = %i\n",nRespCat);
  printf("sig2mu   = %.2f\n",sig2mu);
  printf("a        = %.2f\n",a);
  printf("b        = %.2f\n",b);
  printf("e        = %.2f\n",e);
  printf("f        = %.2f\n",f);
  printf("sigDecorr= %.2f\n",sigDecorr);

  if(nRespCat<=0){printf("\n*****nRespCat must be greater than 0.\n");errs++;}
  if(sig2mu<=0){printf("\n*****sig2mu must be greater than 0.\n");errs++;}
  if(a<=0){printf("\n*****a must be greater than 0.\n");errs++;}
  if(b<=0){printf("\n*****b must be greater than 0.\n");errs++;}
  if(e<=0){printf("\n*****e must be greater than 0.\n");errs++;}
  if(f<=0){printf("\n*****f must be greater than 0.\n");errs++;}
  if(sigDecorr<=0){printf("\n*****sigDecorr must be greater than 0.\n");errs++;}

  if(errs>0){printf("Exiting...\n\n");exit(1);}

  asprintf(&fnEstimates,"%s.est",name);
  asprintf(&fnChains,"%s.chn",name);

  if(argc>10){
    fnInput=argv[10];
  }else{
    asprintf(&fnInput,"%s.dat",name);
  }

  K=nRespCat-1;

  //open files.
  FILE *CHAINS;
  if ((CHAINS = fopen(fnChains,"w")) == NULL) {
    fprintf(stderr,"\n****Cannot open chain output file!\n\n");
    return(1);
  }

  FILE *ESTIMATES;
  if ((ESTIMATES = fopen(fnEstimates,"w")) == NULL) {
    fprintf(stderr,"\n****Cannot open estimate output file!\n\n");
    fclose(CHAINS);
    return(1);
  }

  FILE *INPUT;
  if ((INPUT = fopen(fnInput,"r")) == NULL) {
    fprintf(stderr,"\n****Cannot open input file %s!\n\n",fnInput);
    fclose(CHAINS);
    fclose(ESTIMATES);
    return(1);
  }

  //Now that we have opened the files, get the data.

  i=0;
  int newline=0,spaces=0;
  int got;
  while(!newline&&!feof(INPUT)){
    i++;
    got=fgetc(INPUT);
    //printf("Character %i: %i %c\n",i,got,got);
    if(got==32) spaces++;
    if(got==10) newline=1;
  }
  rewind(INPUT);

  if(!newline){
    printf("Could not determine number of conditions!\n");
    exit(1);
  }
  J=(spaces+1)/2;

  i=0;
  int totsig=0;
  int totnoi=0;

  printf("Reading data...\n");
  while(!feof(INPUT)){
    for(j=0;j<J;j++){
      if(j==(J-1)){
        fscanf(INPUT,"%i %i\n",&yij,&sij);
      }else{
        fscanf(INPUT,"%i %i ",&yij,&sij);
      }
      if((yij>=nRespCat)||(sij!=0&&sij!=1)){
        printf("\n****Invalid data for participant %i, item %i:  (yij=%i,sij=%i)\n\n",i+1,j+1,yij,sij);
        fclose(CHAINS);
        fclose(ESTIMATES);
        fclose(INPUT);
        return(1);
      }
      if(yij!=-1){
	totsig+=sij;
	totnoi+=1-sij;
      }
    }
    i++;
  }
  rewind(INPUT);

  I=i;

  //printf("Read %d subjects and %d items, with %d total signal trials.\n",I,J,totsig);

  //initialize matrices

  printf("Reserving memory for matrices...\n");
  int Ys[totsig],Yn[totnoi],Subs[totsig],Subn[totnoi];
  gsl_vector *Ws,*Wn,*Ts,*Tn,*WsMean,*WnMean,*WsDev,*TsMean,*TnMean;
  gsl_matrix *Xs,*Xn,*Ss,*Sn,*TpSs,*TpSn,*XstXs,*XntXn,*Vs,*Vn;
  Xs = gsl_matrix_calloc(totsig, I+J+1);
  Xn = gsl_matrix_calloc(totnoi, I+J+1);
  XstXs = gsl_matrix_calloc(I+J+1, I+J+1);
  XntXn = gsl_matrix_calloc(I+J+1, I+J+1);
  Vs = gsl_matrix_calloc(I+J+1, I+J+1);
  Vn = gsl_matrix_calloc(I+J+1, I+J+1);
  gsl_matrix *Vn0 = gsl_matrix_calloc(I+J+1, I+J+1);
  gsl_matrix *Vs0 = gsl_matrix_calloc(I+J+1, I+J+1);
  Ws = gsl_vector_calloc(totsig);
  Wn = gsl_vector_calloc(totnoi);
  Ts = gsl_vector_calloc(I+J+1);
  Tn = gsl_vector_calloc(I+J+1);
  WsMean = gsl_vector_calloc(totsig);
  WnMean = gsl_vector_calloc(totnoi);
  WsDev = gsl_vector_calloc(totsig);
  TsMean= gsl_vector_calloc(I+J+1);
  TnMean= gsl_vector_calloc(I+J+1);
  gsl_vector *TsMean0= gsl_vector_calloc(I+J+1);
  gsl_vector *TnMean0= gsl_vector_calloc(I+J+1);

  gsl_vector *isAlph = gsl_vector_calloc(I+J+1);
  gsl_vector *isBeta = gsl_vector_calloc(I+J+1);

  TpSs = gsl_matrix_calloc(I+J+1,I+J+1);
  TpSn = gsl_matrix_calloc(I+J+1,I+J+1);
  gsl_matrix_set_identity(TpSs);
  gsl_matrix_set_identity(TpSn);

  Ss = gsl_matrix_calloc(totsig, totsig);
  gsl_matrix_set_identity(Ss);
  Sn = gsl_matrix_calloc(totnoi, totnoi);
  gsl_matrix_set_identity(Sn);

  gsl_permutation *gslPerm = gsl_permutation_alloc(I+J+1);
  gsl_permutation_init(gslPerm);

  gsl_matrix_set(TpSs,0,0,1.0/sig2mu);
  gsl_matrix_set(TpSn,0,0,1.0/sig2mu);

  //end initialize matrices

  printf("Creating design matrix...\n");
  int y[I][J],s[I][J],sigindex=0,noiindex=0;
  i=0;
  while(!feof(INPUT)){
    for(j=0;j<J;j++){
      if(j==(J-1)){
	fscanf(INPUT,"%i %i\n",&yij,&sij);
      }else{
	fscanf(INPUT,"%i %i ",&yij,&sij);
      }
      y[i][j]=yij;
      s[i][j]=sij;
      if(yij!=-1){
	if(sij){
	  Ys[sigindex] = y[i][j];
	  Subs[sigindex] = i;
	  gsl_matrix_set(Xs,sigindex,0,1);
	  gsl_matrix_set(Xs,sigindex,i+1,1);
	  gsl_matrix_set(Xs,sigindex++,I+1+j,1);
	}else{
	  Yn[noiindex] = y[i][j];
	  Subn[noiindex] = i;
	  gsl_matrix_set(Xn,noiindex,0,1);
	  gsl_matrix_set(Xn,noiindex,i+1,1);
	  gsl_matrix_set(Xn,noiindex++,I+1+j,1);
	}
      }
    }
    i++;
  }

  gsl_blas_dgemm(CblasTrans,CblasNoTrans,1,Xn,Xn,1,XntXn);
  gsl_blas_dgemm(CblasTrans,CblasNoTrans,1,Xs,Xs,1,XstXs);

  //gsl_vector_set(Tn,0,-.6);
  //gsl_vector_set(Ts,0,.6);

  //gsl_matrix_fprintf(stdout,Xs,"%f");

  //exit(0);

  printf("Read %i subjects values for %i items.\n\n",I,J);
  fclose(INPUT);
  ////End data entry.

  //Create starting values
  int keq0=0;
  double crits[I][K],critCand[I][K],w[I][J];
  double maxw[I][K+1];
  double minw[I][K+1];
  double sig2=1;
  double sumWnSqr=0,sumWsSqr=0;
  int nSignal[I],totSignal=0;
  int nRespS[I][K+1];
  int nRespN[I][K+1];
  double sumwN[I];
  double sumwS[I];
  double zDecorr[I][K];
  double sig2N=1;
  double sig2beta0=.2,sig2beta1=.2,sig2alph0=.2,sig2alph1=.2;

  double bDecorr[I];
  int accDecorr[I],badCrits[I];
  double decorrRate=0,sumAlp20=0,sumAlp21=0,sumBet20=0,sumBet21=0;

  for(i=0;i<I;i++){
    accDecorr[i]=0;
    nSignal[i]=0;
    nRespN[i][K]=0;
    nRespS[i][K]=0;
    crits[i][0]=0;
    crits[i][K-1]=SET;
    for(k=0;k<K;k++){
      nRespS[i][k]=0;
      nRespN[i][k]=0;
      if(k!=0 && k<(K-1)){
	crits[i][k]=k*SET/((float)(K-1));
	//printf("i: %i k: %i crit: %f\n",i,k,crits[i][k]);
      }
    }
    for(j=0;j<J;j++){
      nSignal[i]+=s[i][j];
      totSignal+=s[i][j];
      if(y[i][j]!=-1){
	nRespN[i][y[i][j]]+=1-s[i][j];
	nRespS[i][y[i][j]]+=s[i][j];
      }
    }
  }

  //for(k=0;k<(K+1);k++)
  //{
  //  printf("nRespN[0][%d]=%d\n",k,nRespN[0][k]);
  //  printf("nRespS[0][%d]=%d\n",k,nRespS[0][k]);
  //}
  //end starting values

  //Initialize GSL
  const gsl_rng_type * T;
  gsl_rng * r;

  gsl_rng_env_setup();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  //gsl_rng_set(r,seed);
  gsl_rng_set(r,0);

  gsl_blas_dgemv(CblasNoTrans,1,Xs,Ts,0,WsMean);
  gsl_blas_dgemv(CblasNoTrans,1,Xn,Tn,0,WnMean);

  printf("Total signal trials: %i\nTotal noise trials: %i\n\n",totsig,totnoi);

  //begin MCMC loop.
  i=0;
  printf("Starting %i MCMC iterations...\n",ITERS);
  for(iter=0;iter<ITERS;iter++){

    for(i=0;i<I;i++){
      sumwS[i]=0;
      sumwN[i]=0;
      for(k=0;k<(K+1);k++){
	minw[i][k]=GSL_POSINF;
	maxw[i][k]=GSL_NEGINF;
      }
    }
    //sig2=1;

    //sample latent variables

    for(i=0;i<totsig;i++){
      //printf("%d %d %f %f\n",Ys[i],Subs[i],crits[Subs[i]][0],crits[Subs[i]][1]);
      if( Ys[i]==0 ){
	gsl_vector_set(Ws,i,ran_trunc_norm_lower(r, gsl_vector_get(WsMean,i), sqrt(sig2), crits[Subs[i]][0]));
	if(gsl_vector_get(Ws,i)>maxw[Subs[i]][0]) { maxw[Subs[i]][0]=gsl_vector_get(Ws,i);}
      }else if( Ys[i]==K ){
      	gsl_vector_set(Ws,i,ran_trunc_norm_upper(r, gsl_vector_get(WsMean,i), sqrt(sig2), crits[Subs[i]][K-1]));
	if(gsl_vector_get(Ws,i)<minw[Subs[i]][K]) { minw[Subs[i]][K]=gsl_vector_get(Ws,i);}
      }else{
	gsl_vector_set(Ws,i,ran_trunc_norm_both(r, gsl_vector_get(WsMean,i), sqrt(sig2), crits[Subs[i]][Ys[i]-1],crits[Subs[i]][Ys[i]]));
	if(gsl_vector_get(Ws,i)>maxw[Subs[i]][Ys[i]]) { maxw[Subs[i]][Ys[i]]=gsl_vector_get(Ws,i);}
	if(gsl_vector_get(Ws,i)<minw[Subs[i]][Ys[i]]) { minw[Subs[i]][Ys[i]]=gsl_vector_get(Ws,i);}
      }
    }
    for(i=0;i<totnoi;i++){
      //if(i==511) printf("%d %d %f %f\n",Yn[i],Subn[i],crits[Subn[i]][Yn[i]-1],crits[Subs[i]][Yn[i]]);
      if( Yn[i]==0 ){
	gsl_vector_set(Wn,i,ran_trunc_norm_lower(r, gsl_vector_get(WnMean,i), sqrt(sig2N), crits[Subn[i]][0]));
	if(gsl_vector_get(Wn,i)>maxw[Subn[i]][0]) { maxw[Subn[i]][0]=gsl_vector_get(Wn,i);}
      }else if( Yn[i]==K ){
      	gsl_vector_set(Wn,i,ran_trunc_norm_upper(r, gsl_vector_get(WnMean,i), sqrt(sig2N), crits[Subn[i]][K-1]));
	if(gsl_vector_get(Wn,i)<minw[Subn[i]][K]) { minw[Subn[i]][K]=gsl_vector_get(Wn,i);}
      }else{
	gsl_vector_set(Wn,i,ran_trunc_norm_both(r, gsl_vector_get(WnMean,i), sqrt(sig2N), crits[Subn[i]][Yn[i]-1],crits[Subn[i]][Yn[i]]));
	if(gsl_vector_get(Wn,i)>maxw[Subn[i]][Yn[i]]) { maxw[Subn[i]][Yn[i]]=gsl_vector_get(Wn,i);}
	if(gsl_vector_get(Wn,i)<minw[Subn[i]][Yn[i]]) { minw[Subn[i]][Yn[i]]=gsl_vector_get(Wn,i);}
      }
    }
    //printf("***\n");
    //gsl_vector_fprintf(stderr,Wn,"%f");
    //gsl_vector_fprintf(stderr,WnMean,"%f");
    //for(i=0;i<totnoi;i++) printf("%i : %f\n",i,gsl_vector_get(WnMean,i));

    gsl_blas_daxpy(-1,Ws,WsMean);
    gsl_blas_daxpy(-1,Wn,WnMean);

    gsl_blas_ddot(WsMean,WsMean,&sumWsSqr);
    gsl_blas_ddot(WnMean,WnMean,&sumWnSqr);

    //gsl_vector_fprintf(stderr,WnMean,"%f");
    //for(i=0;i<totnoi;i++) printf("%i : %f\n",i,gsl_vector_get(Wn,i));

    //printf("\nSS: %f,%f,%f,%d\n--------\n",sumWnSqr,a,b,totnoi);

    //sample sigma^2
    sig2=ran_inv_gamma(r, a+.5*totsig,b+(.5*sumWsSqr));
    //sig2=1;
    //printf("IG: %f   %f\n",a0+.5*totSignal,(1.0)*b0+(.5*sumWSqr));
    //sig2=1.3;
    fwrite(&sig2,sizeof(double),1,CHAINS);

    sig2N=ran_inv_gamma(r, a+.5*totnoi,b+(.5*sumWnSqr));

    fwrite(&sig2N,sizeof(double),1,CHAINS);

    // set up variances for full conditional on parameter vector
    for(i=0;i<(I+J+1);i++){
      if(i>I){
	gsl_matrix_set(TpSs,i,i,1.0/sig2beta1);
	gsl_matrix_set(TpSn,i,i,1.0/sig2beta0);
      }else if(i>0){
  	gsl_matrix_set(TpSs,i,i,1.0/sig2alph1);
	gsl_matrix_set(TpSn,i,i,1.0/sig2alph0);
      }
      //printf("i: %d  |  %f\n",i,gsl_matrix_get(TpSs,i,i));
    }

    gsl_matrix_memcpy(Vs,TpSs);
    gsl_matrix_memcpy(Vn,TpSn);

    gsl_matrix_scale(Vs,sig2);
    gsl_matrix_scale(Vn,sig2N);
    gsl_matrix_add(Vn,XntXn);
    gsl_matrix_add(Vs,XstXs);
    gsl_matrix_scale(Vs,1.0/sig2);
    gsl_matrix_scale(Vn,1.0/sig2N);

    //for(i=0;i<(I+J+1);i++){
    // for(j=0;j<(I+J+1);j++){
    //	printf("%f ",gsl_matrix_get(Vs,i,j));
    //  }
    //  printf("\n");
    //}

    gsl_linalg_LU_decomp(Vn,gslPerm,&i);
    gsl_linalg_LU_invert(Vn,gslPerm,Vn0);

    gsl_linalg_LU_decomp(Vs,gslPerm,&i);
    gsl_linalg_LU_invert(Vs,gslPerm,Vs0);

    //for(i=0;i<(I+J+1);i++){
    //  for(j=0;j<(I+J+1);j++){
    //	printf("%f ",gsl_matrix_get(Vs0,i,j));
    //  }
    //  printf("\n");
    //}

    //sample parameter vectors

    gsl_blas_dgemv(CblasTrans,1.0/sig2,Xs,Ws,0,TsMean);
    gsl_blas_dgemv(CblasNoTrans,1,Vs0,TsMean,0,TsMean0);

    gsl_blas_dgemv(CblasTrans,1.0/sig2N,Xn,Wn,0,TnMean);
    gsl_blas_dgemv(CblasNoTrans,1,Vn0,TnMean,0,TnMean0);

    pcl_ran_mvgaussian(r, TsMean0, Vs0, Ts);
    pcl_ran_mvgaussian(r, TnMean0, Vn0, Tn);

    gsl_vector_fwrite(CHAINS,Ts);
    gsl_vector_fwrite(CHAINS,Tn);

    //Sample variances of parameters
    sumAlp20=0;
    sumAlp21=0;
    sumBet20=0;
    sumBet21=0;
    for(i=0;i<(I+J+1);i++){
      if(i>I){
	sumBet20+=pow(gsl_vector_get(Tn,i),2);
	sumBet21+=pow(gsl_vector_get(Ts,i),2);
      }else if(i>0){
	sumAlp20+=pow(gsl_vector_get(Tn,i),2);
	sumAlp21+=pow(gsl_vector_get(Ts,i),2);
      }
    }
    //printf("%f\n",sumBet20);

    //sample sigma^2
    sig2alph0=ran_inv_gamma(r, e+.5*I,f+(.5*sumAlp20));
    sig2alph1=ran_inv_gamma(r, e+.5*I,f+(.5*sumAlp21));
    sig2beta0=ran_inv_gamma(r, e+.5*J,f+(.5*sumBet20));
    sig2beta1=ran_inv_gamma(r, e+.5*J,f+(.5*sumBet21));

    //printf("IG: %f   %f\n",a0+.5*totSignal,(1.0)*b0+(.5*sumWSqr));
    //sig2=1.3;
    fwrite(&sig2alph0,sizeof(double),1,CHAINS);
    fwrite(&sig2alph1,sizeof(double),1,CHAINS);
    fwrite(&sig2beta0,sizeof(double),1,CHAINS);
    fwrite(&sig2beta1,sizeof(double),1,CHAINS);

    gsl_blas_dgemv(CblasNoTrans,1,Xs,Ts,0,WsMean);
    gsl_blas_dgemv(CblasNoTrans,1,Xn,Tn,0,WnMean);

    //gsl_vector_fprintf(stderr,WsMean,"%f");
    //printf("\n\n");
    //gsl_vector_fprintf(stderr,WnMean,"%f");
    //printf("\n\n");
    //gsl_vector_fprintf(stderr,WsMean,"%g");

    //printf("\ns2=%f\n%f %f,%f %f\n",sig2,sumwN[i]/(1.0*(J-nSignal[i])),1.0/sqrt(1.0*(J-nSignal[i])),sumwS[i]/(1.0*nSignal[i]),sqrt(sig2/(1.0*(nSignal[i]))));

      //sample criteria
    for(i=0;i<I;i++){
      //zDecorr[i]=gsl_ran_gaussian(r,sigDecorr);
      bDecorr[i]=1;
      badCrits[i]=0;
      crits[i][0]=0; crits[i][K-1]=SET;
      for(k=0;k<K;k++){
	zDecorr[i][k]=gsl_ran_gaussian(r,sigDecorr);
	//printf("i=%d k=%i max=%f min=%f\n",i,k,maxw[i][k],minw[i][k+1]);
	//if((k!=0) && (k!=(K-1))) crits[i][k]=(gsl_rng_uniform_pos(r)*(minw[i][k+1]-maxw[i][k])+maxw[i][k]);
	//crits[i][k]=ran_trunc_norm_both(r,mu_crits,sd_crits,maxw[i][k],minw[i][k+1])*(k!=keq0);
	//printf("newcrit: %f\n",crits[i][k]);
	//for decorrelating
	critCand[i][k]=crits[i][k]+zDecorr[i][k]*(k!=0 && k!=(K-1));
	//bDecorr[i]=bDecorr[i]*exp(-.5*pow((critCand[i][k]-mu_crits)/sd_crits,2)+.5*pow((crits[i][k]-mu_crits)/sd_crits,2));
	//printf("candcrit: %f\n",critCand[i][k]);
	if(k>0 && critCand[i][k]<critCand[i][k-1] ) badCrits[i]=1;
	//printf("CRITS: %d %d d:%f s:%f  :   %f %f\n",iter,k,ds[i],sig2,crits[i][k],critCand[i][k]);
      }
    }

    if(DECORRCRITS){
      for(i=0;i<totsig;i++){
	if(!badCrits[Subs[i]])
	  if(Ys[i]==0){
	    bDecorr[Subs[i]]=bDecorr[Subs[i]]*(
					       gsl_cdf_ugaussian_P((critCand[Subs[i]][0]-gsl_vector_get(WsMean,i))/sqrt(sig2)) /
					       gsl_cdf_ugaussian_P((crits[Subs[i]][0]   -gsl_vector_get(WsMean,i))/sqrt(sig2))
					       );
	  }else if(Ys[i]==K){
	    bDecorr[Subs[i]]=bDecorr[Subs[i]]*(
					       (1-gsl_cdf_ugaussian_P((critCand[Subs[i]][K-1]-gsl_vector_get(WsMean,i))/sqrt(sig2))) /
					       (1-gsl_cdf_ugaussian_P((crits[Subs[i]][K-1]   -gsl_vector_get(WsMean,i))/sqrt(sig2)))
					       );
	  }else{
	    bDecorr[Subs[i]]=bDecorr[Subs[i]]*(
					       (gsl_cdf_ugaussian_P((critCand[Subs[i]][Ys[i]]-gsl_vector_get(WsMean,i))/sqrt(sig2)) - gsl_cdf_ugaussian_P((critCand[Subs[i]][Ys[i]-1]-gsl_vector_get(WsMean,i))/sqrt(sig2))) /
					       (gsl_cdf_ugaussian_P((crits[Subs[i]][Ys[i]]   -gsl_vector_get(WsMean,i))/sqrt(sig2)) - gsl_cdf_ugaussian_P((crits[Subs[i]][Ys[i]-1]   -gsl_vector_get(WsMean,i))/sqrt(sig2)))
					       );
	  }
      }

      for(i=0;i<totnoi;i++){
	if(!badCrits[Subn[i]])
	  if(Yn[i]==0){
	    bDecorr[Subn[i]]=bDecorr[Subn[i]]*(
					       gsl_cdf_ugaussian_P((critCand[Subn[i]][0]-gsl_vector_get(WnMean,i))/sqrt(sig2N)) /
					       gsl_cdf_ugaussian_P((crits[Subn[i]][0]   -gsl_vector_get(WnMean,i))/sqrt(sig2N))
					       );
	  }else if(Yn[i]==K){
	    bDecorr[Subn[i]]=bDecorr[Subn[i]]*(
					       (1-gsl_cdf_ugaussian_P((critCand[Subn[i]][K-1]-gsl_vector_get(WnMean,i))/sqrt(sig2N))) /
					       (1-gsl_cdf_ugaussian_P((crits[Subn[i]][K-1]   -gsl_vector_get(WnMean,i))/sqrt(sig2N)))
					       );
	  }else{
	    bDecorr[Subn[i]]=bDecorr[Subn[i]]*(
					       (gsl_cdf_ugaussian_P((critCand[Subn[i]][Yn[i]]-gsl_vector_get(WnMean,i))/sqrt(sig2N)) - gsl_cdf_ugaussian_P((critCand[Subn[i]][Yn[i]-1]-gsl_vector_get(WnMean,i))/sqrt(sig2N))) /
						(gsl_cdf_ugaussian_P((crits[Subn[i]][Yn[i]]   -gsl_vector_get(WnMean,i))/sqrt(sig2N)) - gsl_cdf_ugaussian_P((crits[Subn[i]][Yn[i]-1]   -gsl_vector_get(WnMean,i))/sqrt(sig2N)))
					       );
	  }
      }

      //for(i=0;i<I;i++)
      //for(k=0;k<K;k++)
      //printf("m: %d, i: %d, k: %d,  z: %f crit: %f  cand: %f, BAD: %d\n",iter,i,k,zDecorr[i][k],crits[i][k],critCand[i][k],badCrits[i]);

      for(i=0;i<I;i++){
	if(!badCrits[i]){

	  accDecorr[i]=gsl_ran_bernoulli(r,(bDecorr[i]>1)?1:bDecorr[i]);
	  //printf("BDECORR: m%d i%d    bdecorr %f\n",iter,i,bDecorr[i]);
	  if(accDecorr[i]){
	    decorrRate+=1/(1.0*I*ITERS);

	    for(k=0;k<K;k++){
	      //printf("m: %d, i: %d, k: %d,  z: %f crit: %f  cand: %f\n",iter,i,k,zDecorr[i][k],crits[i][k],critCand[i][k]);
	      crits[i][k]=critCand[i][k];
	    }
	  }
	}
      }
    }

    fwrite(crits,sizeof(double),I*K,CHAINS);

    if(!((iter+1)%25)) { printf("Iteration %i\n",iter+1); }

  }//End MCMC loop

  printf("\nDone. Decorrelating step acceptance rate: %f\n\n",decorrRate);
  //close files
  fclose(ESTIMATES);
  fclose(CHAINS);
  return 0;
}

max_stars_repo_name: richarddmorey/ConfidenceRatingsSigDet

max_issues_repo_name: richarddmorey/ConfidenceRatingsSigDet

max_forks_repo_name: richarddmorey/ConfidenceRatingsSigDet


text:
/*
 Copyright (C) 2002 M. Marques, A. Castro, A. Rubio, G. Bertsch

 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, 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 _SYMBOLS_H
#define _SYMBOLS_H

#include <gsl/gsl_complex.h>
#include "liboct_parser.h"

typedef enum{
  S_CMPLX, S_STR, S_BLOCK, S_FNCT
} symrec_type;

/* Data type for links in the chain of symbols. */
typedef struct symrec{
  char *name;                  /* name of symbol */
  symrec_type type;            /* type of symbol: complex, string, block, or function */
  int def;                     /* has this symbol been defined */
  int used;                    /* this symbol has been used before */

  int nargs;                   /* if type==S_FNCT contains the number of arguments of the function */

  union {
    gsl_complex c;             /* value of a S_CMPLX */
    char *str;                 /* value of a S_STR */
    sym_block *block;          /* to store blocks */
    gsl_complex (*fnctptr)();  /* value of a S_FNCT */
  } value;

  struct symrec *next;         /* link field */
} symrec;

/* The symbol table: a chain of struct symrec. */
extern symrec *sym_table;
extern char *reserved_symbols[];

symrec *putsym (const char *sym_name, symrec_type sym_type);
symrec *getsym (const char *sym_name);
int      rmsym (const char *sym_name);

void sym_notdef(symrec *sym);
void sym_redef(symrec *sym);
void sym_wrong_arg(symrec *sym);
void sym_init_table(void);
void sym_end_table(void);
void sym_output_table(int only_unused, int mpiv_node);
void str_tolower(char *in);
void sym_mark_table_used();
void sym_print(FILE *f, const symrec *ptr);

#endif

ext: h

max_stars_repo_name: shunsuke-sato/octopus

max_issues_repo_name: shunsuke-sato/octopus

max_forks_repo_name: shunsuke-sato/octopus


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_statistics_double.h>
#include <math.h>
#include "tz_error.h"
#include "tz_constant.h"
#include "tz_voxel_linked_list.h"
#include "tz_stack_sampling.h"
#include "tz_voxel_graphics.h"
#include "tz_neurotrace.h"
#include "tz_stack_math.h"
#include "tz_stack_utils.h"
#include "tz_objdetect.h"
#include "tz_voxeltrans.h"
#include "tz_stack_stat.h"
#include "tz_stack_draw.h"
#include "tz_geo3d_vector.h"
#include "tz_stack_bwmorph.h"
#include "tz_stack_bwdist.h"
#include "tz_stack_threshold.h"

INIT_EXCEPTION_MAIN(e)

int main(int argc, char* argv[])
{
  char neuron_name[100];
  if (argc == 1) {
    strcpy(neuron_name, "fly_neuron");
  } else {
    strcpy(neuron_name, argv[1]);
  }

  char file_path[100];

  sprintf(file_path, "../data/%s/traced.tif", neuron_name);
  Stack *stack = Read_Stack(file_path);

  sprintf(file_path, "../data/%s/soma.tif", neuron_name);
  Stack *label = Read_Stack(file_path);

  Stack *label2 = Copy_Stack(label);

  Rgb_Color color;
  Set_Color(&color, 0, 255, 0);

  Stack_Threshold(label, 2);
  Stack_Binarize(label);

  Stack_Threshold(label2, 1);
  Stack_Binarize(label2);
  Stack_Xor(label, label2, label2);

  sprintf(file_path, "../data/%s.tif", neuron_name);
  Stack *signal = Read_Stack(file_path);
  //label = Stack_Boundary_N(label, NULL, 8);
  //Stack_Label_Bwc(stack, label, color);
  Stack_Label_Color(stack, label, 5.4, 1.0, signal);

  //label2 = Stack_Boundary_N(label2, NULL, 8);
  //Set_Color(&color, 0, 0, 255);
  //Stack_Label_Bwc(stack, label2, color);
  Stack_Label_Color(stack, label2, 3.6, 1.0, signal);

  sprintf(file_path, "../data/%s/label.tif", neuron_name);
  Write_Stack(file_path, stack);

  Kill_Stack(signal);
  Kill_Stack(stack);
  Kill_Stack(label);
  Kill_Stack(label2);

  return 0;
}

ext: c


text:
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>

#ifdef SUBFIND

#include "allvars.h"
#include "proto.h"
#include "subfind.h"

/*! Structure for communication during the density computation. Holds data that is sent to other processors.
 */
static struct linkngbdata_in
{
  MyDouble Pos[3];
  MyFloat DM_Hsml;
  int NodeList[NODELISTLENGTH];
}
 *LinkngbDataIn, *LinkngbDataGet;

static struct linkngbdata_out
{
  int Ngb;
}
 *LinkngbDataResult, *LinkngbDataOut;

void subfind_find_linkngb(void)
{
  long long ntot;
  int i, j, ndone, ndone_flag, npleft, dummy, iter = 0, save_DesNumNgb;
  MyFloat *Left, *Right;
  char *Todo;
  int ngrp, sendTask, recvTask, place, nexport, nimport;
  double dmax1, dmax2, t0, t1;

  if(ThisTask == 0)
    printf("Start find_linkngb (%d particles on task=%d)\n", NumPartGroup, ThisTask);

  save_DesNumNgb = All.DesNumNgb;
  All.DesNumNgb = All.DesLinkNgb;	/* for simplicity, reset this value */

  /* allocate buffers to arrange communication */

  Ngblist = (int *) mymalloc(NumPartGroup * sizeof(int));
  Dist2list = (double *) mymalloc(NumPartGroup * sizeof(double));

  All.BunchSize =
    (int) ((All.BufferSize * 1024 * 1024) / (sizeof(struct data_index) + sizeof(struct data_nodelist) +
					     sizeof(struct linkngbdata_in) + sizeof(struct linkngbdata_out) +
					     sizemax(sizeof(struct linkngbdata_in),
						     sizeof(struct linkngbdata_out))));
  DataIndexTable = (struct data_index *) mymalloc(All.BunchSize * sizeof(struct data_index));
  DataNodeList = (struct data_nodelist *) mymalloc(All.BunchSize * sizeof(struct data_nodelist));

  Left = mymalloc(sizeof(MyFloat) * NumPartGroup);
  Right = mymalloc(sizeof(MyFloat) * NumPartGroup);
  Todo = mymalloc(sizeof(char) * NumPartGroup);

  for(i = 0; i < NumPartGroup; i++)
    {
      Left[i] = Right[i] = 0;
      Todo[i] = 1;
    }

  /* we will repeat the whole thing for those particles where we didn't find enough neighbours */
  do
    {
      t0 = second();

      i = 0;			/* begin with this index */

      do
	{
	  for(j = 0; j < NTask; j++)
	    {
	      Send_count[j] = 0;
	      Exportflag[j] = -1;
	    }

	  /* do local particles and prepare export list */

	  for(nexport = 0; i < NumPartGroup; i++)
	    {
	      if(Todo[i])
		{
		  if(subfind_linkngb_evaluate(i, 0, &nexport, Send_count) < 0)
		    break;
		}
	    }

	  qsort(DataIndexTable, nexport, sizeof(struct data_index), data_index_compare);

	  MPI_Allgather(Send_count, NTask, MPI_INT, Sendcount_matrix, NTask, MPI_INT, MPI_COMM_WORLD);

	  for(j = 0, nimport = 0, Recv_offset[0] = 0, Send_offset[0] = 0; j < NTask; j++)
	    {
	      Recv_count[j] = Sendcount_matrix[j * NTask + ThisTask];
	      nimport += Recv_count[j];

	      if(j > 0)
		{
		  Send_offset[j] = Send_offset[j - 1] + Send_count[j - 1];
		  Recv_offset[j] = Recv_offset[j - 1] + Recv_count[j - 1];
		}
	    }

	  LinkngbDataGet = (struct linkngbdata_in *) mymalloc(nimport * sizeof(struct linkngbdata_in));
	  LinkngbDataIn = (struct linkngbdata_in *) mymalloc(nexport * sizeof(struct linkngbdata_in));

	  /* prepare particle data for export */
	  for(j = 0; j < nexport; j++)
	    {
	      place = DataIndexTable[j].Index;

	      LinkngbDataIn[j].Pos[0] = P[place].Pos[0];
	      LinkngbDataIn[j].Pos[1] = P[place].Pos[1];
	      LinkngbDataIn[j].Pos[2] = P[place].Pos[2];
	      LinkngbDataIn[j].DM_Hsml = P[place].DM_Hsml;

	      memcpy(LinkngbDataIn[j].NodeList,
		     DataNodeList[DataIndexTable[j].IndexGet].NodeList, NODELISTLENGTH * sizeof(int));
	    }

	  /* exchange particle data */
	  for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	    {
	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;

	      if(recvTask < NTask)
		{
		  if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		    {
		      /* get the particles */
		      MPI_Sendrecv(&LinkngbDataIn[Send_offset[recvTask]],
				   Send_count[recvTask] * sizeof(struct linkngbdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A,
				   &LinkngbDataGet[Recv_offset[recvTask]],
				   Recv_count[recvTask] * sizeof(struct linkngbdata_in), MPI_BYTE,
				   recvTask, TAG_DENS_A, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		    }
		}
	    }

	  myfree(LinkngbDataIn);
	  LinkngbDataResult = (struct linkngbdata_out *) mymalloc(nimport * sizeof(struct linkngbdata_out));
	  LinkngbDataOut = (struct linkngbdata_out *) mymalloc(nexport * sizeof(struct linkngbdata_out));

	  /* now do the particles that were sent to us */
	  for(j = 0; j < nimport; j++)
	    subfind_linkngb_evaluate(j, 1, &dummy, &dummy);

	  if(i >= NumPartGroup)
	    ndone_flag = 1;
	  else
	    ndone_flag = 0;

	  MPI_Allreduce(&ndone_flag, &ndone, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

	  /* get the result */
	  for(ngrp = 1; ngrp < (1 << PTask); ngrp++)
	    {
	      sendTask = ThisTask;
	      recvTask = ThisTask ^ ngrp;
	      if(recvTask < NTask)
		{
		  if(Send_count[recvTask] > 0 || Recv_count[recvTask] > 0)
		    {
		      /* send the results */
		      MPI_Sendrecv(&LinkngbDataResult[Recv_offset[recvTask]],
				   Recv_count[recvTask] * sizeof(struct linkngbdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B,
				   &LinkngbDataOut[Send_offset[recvTask]],
				   Send_count[recvTask] * sizeof(struct linkngbdata_out),
				   MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		    }
		}
	    }

	  /* add the result to the local particles */
	  for(j = 0; j < nexport; j++)
	    {
	      place = DataIndexTable[j].Index;

	      P[place].DM_NumNgb += LinkngbDataOut[j].Ngb;
	    }

	  myfree(LinkngbDataOut);
	  myfree(LinkngbDataResult);
	  myfree(LinkngbDataGet);
	}
      while(ndone < NTask);

      /* do final operations on results */
      for(i = 0, npleft = 0; i < NumPartGroup; i++)
	{
	  /* now check whether we had enough neighbours */
	  if(Todo[i])
	    {
	      if(P[i].DM_NumNgb != All.DesLinkNgb &&
		 ((Right[i] - Left[i]) > 1.0e-3 * Left[i] || Left[i] == 0 || Right[i] == 0))
		{
		  /* need to redo this particle */
		  npleft++;

		  if(P[i].DM_NumNgb < All.DesLinkNgb)
		    Left[i] = DMAX(P[i].DM_Hsml, Left[i]);
		  else
		    {
		      if(Right[i] != 0)
			{
			  if(P[i].DM_Hsml < Right[i])
			    Right[i] = P[i].DM_Hsml;
			}
		      else
			Right[i] = P[i].DM_Hsml;
		    }

		  if(iter >= MAXITER - 10)
		    {
		      printf
			("i=%d task=%d ID=%d DM_Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n   pos=(%g|%g|%g)\n",
			 i, ThisTask, (int) P[i].ID, P[i].DM_Hsml, Left[i], Right[i],
			 (double) P[i].DM_NumNgb, Right[i] - Left[i], P[i].Pos[0], P[i].Pos[1], P[i].Pos[2]);
		      fflush(stdout);
		    }

		  if(Right[i] > 0 && Left[i] > 0)
		    P[i].DM_Hsml = pow(0.5 * (pow(Left[i], 3) + pow(Right[i], 3)), 1.0 / 3);
		  else
		    {
		      if(Right[i] == 0 && Left[i] == 0)
			endrun(8189);	/* can't occur */

		      if(Right[i] == 0 && Left[i] > 0)
			P[i].DM_Hsml *= 1.26;

		      if(Right[i] > 0 && Left[i] == 0)
			P[i].DM_Hsml /= 1.26;
		    }
		}
	      else
		Todo[i] = 0;
	    }
	}

      sumup_large_ints(1, &npleft, &ntot);

      t1 = second();

      if(ntot > 0)
	{
	  iter++;

	  if(iter > 0 && ThisTask == 0)
	    {
	      printf("find linkngb iteration %d: need to repeat for %d%09d particles. (took %g sec)\n", iter,
		     (int) (ntot / 1000000000), (int) (ntot % 1000000000), timediff(t0, t1));
	      fflush(stdout);
	    }

	  if(iter > MAXITER)
	    {
	      printf("failed to converge in neighbour iteration in density()\n");
	      fflush(stdout);
	      endrun(1155);
	    }
	}
    }
  while(ntot > 0);

  myfree(Todo);
  myfree(Right);
  myfree(Left);

  myfree(DataNodeList);
  myfree(DataIndexTable);

  myfree(Dist2list);
  myfree(Ngblist);

  All.DesNumNgb = save_DesNumNgb;	/* restore it */
}

/*! This function represents the core of the SPH density computation. The
 *  target particle may either be local, or reside in the communication
 *  buffer.
 */
int subfind_linkngb_evaluate(int target, int mode, int *nexport, int *nsend_local)
{
  int startnode, numngb, ngbs, listindex = 0;
  double hmax;
  double h, h2, hinv, hinv3;
  MyDouble *pos;

  numngb = 0;

  if(mode == 0)
    {
      pos = P[target].Pos;
      h = P[target].DM_Hsml;
    }
  else
    {
      pos = LinkngbDataGet[target].Pos;
      h = LinkngbDataGet[target].DM_Hsml;
    }

  h2 = h * h;
  hinv = 1.0 / h;
  hinv3 = hinv * hinv * hinv;

  if(mode == 0)
    {
      startnode = All.MaxPart;	/* root node */
    }
  else
    {
      startnode = LinkngbDataGet[target].NodeList[0];
      startnode = Nodes[startnode].u.d.nextnode;	/* open it */
    }

  numngb = 0;

  while(startnode >= 0)
    {
      while(startnode >= 0)
	{
	  ngbs = subfind_ngb_treefind_linkngb(pos, h, target, &startnode, mode, &hmax, nexport, nsend_local);

	  if(ngbs < 0)
	    return -1;

	  if(mode == 0 && hmax > 0)
	    P[target].DM_Hsml = hmax;

	  numngb += ngbs;
	}

      if(mode == 1)
	{
	  listindex++;
	  if(listindex < NODELISTLENGTH)
	    {
	      startnode = LinkngbDataGet[target].NodeList[listindex];
	      if(startnode >= 0)
		startnode = Nodes[startnode].u.d.nextnode;	/* open it */
	    }
	}
    }

  if(mode == 0)
    {
      P[target].DM_NumNgb = numngb;
    }
  else
    {
      LinkngbDataResult[target].Ngb = numngb;
    }

  return 0;
}

int subfind_ngb_treefind_linkngb(MyDouble searchcenter[3], double hsml, int target, int *startnode, int mode,
				 double *hmax, int *nexport, int *nsend_local)
{
  int numngb, i, no, p, task, nexport_save, exported = 0;
  struct NODE *current;
  double dx, dy, dz, dist, r2;

#ifdef PERIODIC
  MyDouble xtmp;
#endif

  nexport_save = *nexport;

  *hmax = 0;
  numngb = 0;
  no = *startnode;

  while(no >= 0)
    {
      if(no < All.MaxPart)	/* single particle */
	{
	  p = no;
	  no = Nextnode[no];

#ifdef DENSITY_SPLIT_BY_TYPE
	  if(!((1 << P[p].Type) & (DENSITY_SPLIT_BY_TYPE)))
#else
	  if(!((1 << P[p].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
	    continue;

	  dist = hsml;
	  dx = NGB_PERIODIC_LONG_X(P[p].Pos[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(P[p].Pos[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(P[p].Pos[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  if((r2 = (dx * dx + dy * dy + dz * dz)) > dist * dist)
	    continue;

	  Dist2list[numngb] = r2;
	  Ngblist[numngb++] = p;
	}
      else
	{
	  if(no >= All.MaxPart + MaxNodes)	/* pseudo particle */
	    {
	      if(mode == 1)
		endrun(12312);

	      if(target >= 0)	/* if no target is given, export will not occur */
		{
		  exported = 1;
		  if(Exportflag[task = DomainTask[no - (All.MaxPart + MaxNodes)]] != target)
		    {
		      Exportflag[task] = target;
		      Exportnodecount[task] = NODELISTLENGTH;
		    }

		  if(Exportnodecount[task] == NODELISTLENGTH)
		    {
		      if(*nexport >= All.BunchSize)
			{
			  *nexport = nexport_save;
			  if(nexport_save == 0)
			    endrun(13004);	/* in this case, the buffer is too small to process even a single particle */
			  for(task = 0; task < NTask; task++)
			    nsend_local[task] = 0;
			  for(no = 0; no < nexport_save; no++)
			    nsend_local[DataIndexTable[no].Task]++;
			  return -1;
			}
		      Exportnodecount[task] = 0;
		      Exportindex[task] = *nexport;
		      DataIndexTable[*nexport].Task = task;
		      DataIndexTable[*nexport].Index = target;
		      DataIndexTable[*nexport].IndexGet = *nexport;
		      *nexport = *nexport + 1;
		      nsend_local[task]++;
		    }

		  DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]++] =
		    DomainNodeIndex[no - (All.MaxPart + MaxNodes)];

		  if(Exportnodecount[task] < NODELISTLENGTH)
		    DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]] = -1;
		}

	      no = Nextnode[no - MaxNodes];
	      continue;
	    }

	  current = &Nodes[no];

	  if(mode == 1)
	    {
	      if(current->u.d.bitflags & (1 << BITFLAG_TOPLEVEL))	/* we reached a top-level node again, which means that we are done with the branch */
		{
		  *startnode = -1;
		  return numngb;
		}
	    }

	  no = current->u.d.sibling;	/* in case the node can be discarded */

	  dist = hsml + 0.5 * current->len;;
	  dx = NGB_PERIODIC_LONG_X(current->center[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(current->center[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(current->center[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  /* now test against the minimal sphere enclosing everything */
	  dist += FACT1 * current->len;
	  if(dx * dx + dy * dy + dz * dz > dist * dist)
	    continue;

	  no = current->u.d.nextnode;	/* ok, we need to open the node */
	}
    }

  if(mode == 0)			/* local particle */
    if(exported == 0)		/* completely local */
      if(numngb >= All.DesNumNgb)
	{
	  R2list = mymalloc(sizeof(struct r2data) * numngb);
	  for(i = 0; i < numngb; i++)
	    {
	      R2list[i].index = Ngblist[i];
	      R2list[i].r2 = Dist2list[i];
	    }

	  qsort(R2list, numngb, sizeof(struct r2data), subfind_ngb_compare_dist);

	  *hmax = sqrt(R2list[All.DesNumNgb - 1].r2);
	  numngb = All.DesNumNgb;

	  for(i = 0; i < numngb; i++)
	    {
	      Ngblist[i] = R2list[i].index;
	      Dist2list[i] = R2list[i].r2;
	    }

	  myfree(R2list);
	}

  *startnode = -1;
  return numngb;
}

#ifdef DENSITY_SPLIT_BY_TYPE
/*! This routine finds all neighbours `j' that can interact with
 *  \f$ r_{ij} < h_i \f$  OR if  \f$ r_{ij} < h_j \f$.
 */
int subfind_ngb_treefind_linkpairs(MyDouble searchcenter[3], double hsml, int target, int *startnode,
				   int mode, double *hmax, int *nexport, int *nsend_local)
{
  int numngb, i, no, p, task, nexport_save, exported = 0;
  struct NODE *current;
  double dx, dy, dz, dist, r2, dmax1, dmax2;

#ifdef PERIODIC
  MyDouble xtmp;
#endif

  nexport_save = *nexport;

  *hmax = 0;
  numngb = 0;
  no = *startnode;

  while(no >= 0)
    {
      if(no < All.MaxPart)	/* single particle */
	{
	  p = no;
	  no = Nextnode[no];

#ifdef DENSITY_SPLIT_BY_TYPE
	  if(!((1 << P[p].Type) & (DENSITY_SPLIT_BY_TYPE)))
#else
	  if(!((1 << P[p].Type) & (FOF_PRIMARY_LINK_TYPES)))
#endif
	    continue;

	  dist = DMAX(P[p].DM_Hsml, hsml);
	  dx = NGB_PERIODIC_LONG_X(P[p].Pos[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(P[p].Pos[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(P[p].Pos[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  if((r2 = (dx * dx + dy * dy + dz * dz)) > dist * dist)
	    continue;

	  Dist2list[numngb] = r2;
	  Ngblist[numngb++] = p;
	}
      else
	{
	  if(no >= All.MaxPart + MaxNodes)	/* pseudo particle */
	    {
	      if(mode == 1)
		endrun(12312);

	      if(target >= 0)	/* if no target is given, export will not occur */
		{
		  exported = 1;
		  if(Exportflag[task = DomainTask[no - (All.MaxPart + MaxNodes)]] != target)
		    {
		      Exportflag[task] = target;
		      Exportnodecount[task] = NODELISTLENGTH;
		    }

		  if(Exportnodecount[task] == NODELISTLENGTH)
		    {
		      if(*nexport >= All.BunchSize)
			{
			  *nexport = nexport_save;
			  if(nexport_save == 0)
			    endrun(13004);	/* in this case, the buffer is too small to process even a single particle */
			  for(task = 0; task < NTask; task++)
			    nsend_local[task] = 0;
			  for(no = 0; no < nexport_save; no++)
			    nsend_local[DataIndexTable[no].Task]++;
			  return -1;
			}
		      Exportnodecount[task] = 0;
		      Exportindex[task] = *nexport;
		      DataIndexTable[*nexport].Task = task;
		      DataIndexTable[*nexport].Index = target;
		      DataIndexTable[*nexport].IndexGet = *nexport;
		      *nexport = *nexport + 1;
		      nsend_local[task]++;
		    }

		  DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]++] =
		    DomainNodeIndex[no - (All.MaxPart + MaxNodes)];

		  if(Exportnodecount[task] < NODELISTLENGTH)
		    DataNodeList[Exportindex[task]].NodeList[Exportnodecount[task]] = -1;
		}

	      no = Nextnode[no - MaxNodes];
	      continue;
	    }

	  current = &Nodes[no];

	  if(mode == 1)
	    {
	      if(current->u.d.bitflags & (1 << BITFLAG_TOPLEVEL))	/* we reached a top-level node again, which means that we are done with the branch */
		{
		  *startnode = -1;
		  return numngb;
		}
	    }

	  dist = DMAX(Extnodes[no].hmax, hsml) + 0.5 * current->len;
	  no = current->u.d.sibling;	/* in case the node can be discarded */
	  dx = NGB_PERIODIC_LONG_X(current->center[0] - searchcenter[0]);
	  if(dx > dist)
	    continue;
	  dy = NGB_PERIODIC_LONG_Y(current->center[1] - searchcenter[1]);
	  if(dy > dist)
	    continue;
	  dz = NGB_PERIODIC_LONG_Z(current->center[2] - searchcenter[2]);
	  if(dz > dist)
	    continue;
	  /* now test against the minimal sphere enclosing everything */
	  dist += FACT1 * current->len;
	  if(dx * dx + dy * dy + dz * dz > dist * dist)
	    continue;

	  no = current->u.d.nextnode;	/* ok, we need to open the node */
	}
    }

  if(mode == 0)			/* local particle */
    if(exported == 0)		/* completely local */
      if(numngb >= All.DesNumNgb)
	{
	  R2list = mymalloc(sizeof(struct r2data) * numngb);
	  for(i = 0; i < numngb; i++)
	    {
	      R2list[i].index = Ngblist[i];
	      R2list[i].r2 = Dist2list[i];
	    }

	  qsort(R2list, numngb, sizeof(struct r2data), subfind_ngb_compare_dist);

	  *hmax = sqrt(R2list[All.DesNumNgb - 1].r2);
	  numngb = All.DesNumNgb;

	  for(i = 0; i < numngb; i++)
	    {
	      Ngblist[i] = R2list[i].index;
	      Dist2list[i] = R2list[i].r2;
	    }

	  myfree(R2list);
	}

  *startnode = -1;
  return numngb;
}
#endif

int subfind_ngb_compare_dist(const void *a, const void *b)
{
  if(((struct r2data *) a)->r2 < (((struct r2data *) b)->r2))
    return -1;

  if(((struct r2data *) a)->r2 > (((struct r2data *) b)->r2))
    return +1;

  return 0;
}

#endif

max_stars_repo_name: egpbos/egp

max_issues_repo_name: egpbos/egp

max_forks_repo_name: egpbos/egp


text:
/**
 * @file      EigenAdaptors.h
 *
 * @brief     This header contains adaptor functions for Eigen types
 *
 * @author    Stefan Reinhold
 * @copyright Copyright (C) 2018 Stefan Reinhold  -- All Rights Reserved.
 *            You may use, distribute and modify this code under the terms of
 *            the AFL 3.0 license; see LICENSE for full license details.
 */

#pragma once

#include <Eigen/Core>
#include <gsl/gsl>

#include <array>

namespace CortidQCT {
namespace Internal {

namespace Adaptor {

namespace Detail_ {

template <class T, class I>
constexpr inline auto subscript(T &&x,
                                I i) noexcept(noexcept(std::declval<T>()[i])) {
  return x[i];
}

template <class Derived, class I>
constexpr inline auto
subscript(Eigen::MatrixBase<Derived> const &x,
          I i) noexcept(noexcept(x(std::declval<Eigen::Index>()))) {
  return x(gsl::narrow_cast<Eigen::Index>(i));
}

template <class Derived, class I>
constexpr inline auto
subscript(Eigen::MatrixBase<Derived> &&x,
          I i) noexcept(noexcept(x(std::declval<Eigen::Index>()))) {
  return x(gsl::narrow_cast<Eigen::Index>(i));
}

template <class T, std::size_t N> struct ConstructionHelper {
  template <class P, class... Args>
  inline constexpr T operator()(P &&param, Args &&... args) const
      noexcept(noexcept(ConstructionHelper<T, N - 1>{}(
          std::forward<P>(param), subscript(param, N - 1),
          std::forward<Args>(args)...))) {
    return ConstructionHelper<T, N - 1>{}(std::forward<P>(param),
                                          subscript(param, N - 1),
                                          std::forward<Args>(args)...);
  }
};

template <class T> struct ConstructionHelper<T, 0> {
  template <class P, class... Args>
  inline constexpr T operator()(P &&, Args &&... args) const
      noexcept(noexcept(T{std::forward<Args>(args)...})) {
    return T{std::forward<Args>(args)...};
  }
};

template <class T, std::size_t N>
struct ConstructionHelper<std::array<T, N>, 0> {
  template <class P, class... Args>
  inline constexpr std::array<T, N> operator()(P &&, Args &&... args) const
      noexcept(noexcept(std::array<T, N>{{std::forward<Args>(args)...}})) {
    return std::array<T, N>{{std::forward<Args>(args)...}};
  }
};

} // namespace Detail_

/// Converts an std::array based vector into an Eigen vector
template <class T, std::size_t N>
inline Eigen::Matrix<T, static_cast<Eigen::Index>(N), 1>
vec(std::array<T, N> const &arr) {
  using Result = Eigen::Matrix<T, static_cast<Eigen::Index>(N), 1>;
  return Detail_::ConstructionHelper<Result, N>{}(arr);
}

/// Converts an Eigen vector to a std::array based one
template <class Derived>
inline std::array<typename Derived::Scalar,
                  static_cast<std::size_t>(Derived::RowsAtCompileTime)>
arr(Eigen::MatrixBase<Derived> const &vec) noexcept {
  static_assert(Derived::RowsAtCompileTime != Eigen::Dynamic &&
                    Derived::ColsAtCompileTime == 1,
                "Conversion only available for fix sized vectors");

  using T = typename Derived::Scalar;
  auto constexpr N = static_cast<std::size_t>(Derived::RowsAtCompileTime);
  using Result = std::array<T, N>;
  return Detail_::ConstructionHelper<Result, N>{}(vec);
}

template <class T>
inline Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> const>
map(std::vector<T> const &v) {
  return Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1> const>{
      v.data(), gsl::narrow_cast<Eigen::Index>(v.size())};
}

template <class T>
inline Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1>> map(std::vector<T> &v) {
  return Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, 1>>{
      v.data(), gsl::narrow_cast<Eigen::Index>(v.size())};
}

template <std::size_t N, class T>
inline Eigen::Map<
    Eigen::Matrix<T, static_cast<Eigen::Index>(N), Eigen::Dynamic> const>
map(std::vector<std::array<T, N>> const &m) {
  return Eigen::Map<
      Eigen::Matrix<T, static_cast<Eigen::Index>(N), Eigen::Dynamic> const>{
      reinterpret_cast<T const *>(m.data()), gsl::narrow_cast<Eigen::Index>(N),
      gsl::narrow_cast<Eigen::Index>(m.size())};
}

template <std::size_t N, class T>
inline Eigen::Map<
    Eigen::Matrix<T, static_cast<Eigen::Index>(N), Eigen::Dynamic>>
map(std::vector<std::array<T, N>> &m) {
  return Eigen::Map<
      Eigen::Matrix<T, static_cast<Eigen::Index>(N), Eigen::Dynamic>>{
      reinterpret_cast<T *>(m.data()), gsl::narrow_cast<Eigen::Index>(N),
      gsl::narrow_cast<Eigen::Index>(m.size())};
}

} // namespace Adaptor
} // namespace Internal
} // namespace CortidQCT

max_stars_repo_name: ithron/CortidQCT

max_issues_repo_name: ithron/CortidQCT

max_forks_repo_name: ithron/CortidQCT


text:
/*
 * cube.c -- solutions of a cubic system
 *
 * (c) 2016 Prof Dr Andreas Mueller, Hochschule Rapperswil
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

int	f(double x, const double y[], double f[], void *params) {
	f[0] = (y[0] - 1) * y[0] * (y[0] + 1);
	return GSL_SUCCESS;
}

int	singlecurve(FILE *out, const char *name, double x0, double y0) {
	printf("%s: %f\n", name, y0);
	gsl_odeiv2_system	system = { f, NULL, 1, NULL };
	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			1e-6, 1e-6, 0.0);
	fprintf(out, "path %s;\n", name);
	double	x = x0;
	double	y[1] = { y0 };
	fprintf(out, "%s = (%.3f,%.3f)", name, x, y[0]);
	for (int i = 1; i <= 400; i++) {
		double	xnext = x + 0.01;
		int	status = gsl_odeiv2_driver_apply(driver, &x, xnext, y);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "error: rc = %d\n", status);
			fprintf(out, ";\n");
			return -1;
		}
		if (fabs(y[0]) > 100) {
			goto end;
		}
		fprintf(out, "\n    --(%.3f,%.3f)", x, y[0]);
	}
end:
	fprintf(out, ";\n");
	gsl_odeiv2_driver_free(driver);
	return 1;
}

int	centercurve(FILE *out, const char *name, double y0) {
	gsl_odeiv2_system	system = { f, NULL, 1, NULL };
	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			-1e-6, 1e-6, 0.0);
	double	x = 0;
	double	y[1] = { y0 };
	double	xnext = -2;
	// integrate backwards to the initial value
	int	status = gsl_odeiv2_driver_apply(driver, &x, xnext, y);
	gsl_odeiv2_driver_free(driver);
	driver = gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			1e-6, 1e-6, 0.0);
/*
	if (fabs(y[0]) > 100) {
		return -1;
	}
*/
	singlecurve(out, name, -2, y[0]);
	return 1;
}

int	main(int argc, char *argv[]) {
	FILE	*out = fopen("cube.mp", "w");
	char	suffix = 'a';
	char	pathname[6];

	for (int i = 0; i < 16; i++) {
		double	y0 = -1.5 + i * 0.2;
		snprintf(pathname, sizeof(pathname), "path%c", suffix++);
		centercurve(out, pathname, y0);
	}

	return EXIT_SUCCESS;
}

max_stars_repo_name: MatthiasRubin/SeminarDGL

max_issues_repo_name: MatthiasRubin/SeminarDGL

max_forks_repo_name: MatthiasRubin/SeminarDGL


text:
//
// Created by robin on 2018/9/1.
// Copyright (c) 2018 Robin. All rights reserved.
//

#pragma once

#include <type_traits>
#include <vector>
#include <algorithm>
#include <gsl/gsl>

namespace satz::bytes {

/**
 * @brief Generate n random bytes.
 */
std::vector<uint8_t> make_random_bytes (int n);

/**
 * @brief Given a byte sequence, return the indices of non-zero bits, regarding
 *    the bytes as a little-endian bit array with index starting at 0.
 * @param bytes
 * @post The resultant vector should be in increasing order.
 */
template <typename T>
std::vector<T> make_indices (gsl::span<const uint8_t> bytes) {
  std::vector<T> d;
  d.reserve(bytes.size() * 8);

  for (int i = 0; i < bytes.size(); ++i) {
    for (int j = 0; j < 8; ++j) {
      if ((bytes[i] >> j) & 0x1) d.push_back(8 * i + j);
    }
  }

  return d;
}

/**
 * @brief Given a byte sequence, join the bytes into a single object.
 *    If the number of bytes is greater than the size of T, discard extra
 *    bytes. If the number of bytes is less than the size of T, pad with
 *    leading zeros.
 * @param bytes
 * @example
 *    [0xef, 0xbe, 0xed, 0xfe] ==> 0xfeedbeef
 *    [0xed, 0xfe] ==> 0x0000feed
 *    [0xef, 0xbe, 0xed, 0xfe, 0xee] => 0xfeedbeef
 * @post The relative byte order of the result should be the same as the argument.
 */
template <typename T>
T from_bytes (gsl::span<const uint8_t> bytes) {

  static_assert(std::is_trivially_copyable_v<T>);
  using U = std::make_unsigned_t<T>;
  static_assert(sizeof(U) == sizeof(T));

  U             ret = 0;
  const ssize_t m   = bytes.size();
  const ssize_t n   = sizeof(U);

  std::for_each(bytes.rbegin() + std::max(m - n, ssize_t(0)),
                bytes.rend(),
                [&ret] (uint8_t b) {
                  ret = (ret << 8) | b;
                });

  return static_cast<T>(ret);
}

}

max_stars_repo_name: lie-yan/rabin-fingerprint

max_issues_repo_name: lie-yan/rabin-fingerprint

max_forks_repo_name: lie-yan/rabin-fingerprint


text:
/*
 * C version of Diffusive Nested Sampling (DNest4) by Brendon J. Brewer
 *
 * Yan-Rong Li, liyanrong@mail.ihep.ac.cn
 * Jun 30, 2016
 *
 */
#ifndef _DNESTVARS_H
#define _DNESTVARS_H

#ifdef __cplusplus
extern "C" {
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include <gsl/gsl_rng.h>

#define DNEST_MAJOR_VERSION 0  // Dec 2, 2018
#define DNEST_MINOR_VERSION 1
#define DNEST_PATCH_VERSION 0

#define STR_MAX_LENGTH (100)
#define BUF_MAX_LENGTH (200)
#define LEVEL_NUM_MAX (2000)

/* output files */
extern FILE *fsample, *fsample_info;

/* random number generator */
extern const gsl_rng_type * dnest_gsl_T;
extern gsl_rng * dnest_gsl_r;

typedef struct
{
  double value;
  double tiebreaker;
}LikelihoodType;

typedef struct
{
  LikelihoodType log_likelihood;
  double log_X;
  unsigned long long int visits, exceeds;
  unsigned long long int accepts, tries;
}Level;

// struct for options
typedef struct
{
  unsigned int num_particles;
  unsigned int new_level_interval;
  unsigned int save_interval;
  unsigned int thread_steps;
  unsigned int max_num_levels;
  double lambda, beta;
  unsigned int max_num_saves;
  double max_ptol;

  char sample_file[STR_MAX_LENGTH];
  char sample_info_file[STR_MAX_LENGTH];
  char levels_file[STR_MAX_LENGTH];
  char sampler_state_file[STR_MAX_LENGTH];
  char posterior_sample_file[STR_MAX_LENGTH];
  char posterior_sample_info_file[STR_MAX_LENGTH];
  char limits_file[STR_MAX_LENGTH];
}Options;
extern Options options;
extern char options_file[STR_MAX_LENGTH];

typedef struct
{
  void (*from_prior)(void *model, const void *arg);
  double (*log_likelihoods_cal)(const void *model, const void *arg);
  double (*log_likelihoods_cal_initial)(const void *model, const void *arg);
  double (*log_likelihoods_cal_restart)(const void *model, const void *arg);
  double (*perturb)(void *model, const void *arg);
  void (*print_particle)(FILE *fp, const void *model, const void *arg);
  void (*read_particle)(FILE *fp, void *model);
  void (*restart_action)(int iflag);
  void (*accept_action)();
  void (*kill_action)(int i, int i_copy);
}DNestFptrSet;

extern void *particles;
extern int dnest_size_of_modeltype;
extern int particle_offset_size, particle_offset_double;

// sampler
extern bool save_to_disk;
extern unsigned int num_threads;
extern double compression;
extern unsigned int regularisation;

extern void *particles;
extern LikelihoodType *log_likelihoods;
extern unsigned int *level_assignments;

// number account of unaccepted times
extern unsigned int *account_unaccepts;

extern int size_levels;
extern Level *levels;
extern unsigned int count_saves, num_saves, num_saves_restart;
extern unsigned long long int count_mcmc_steps;
extern LikelihoodType *above;
extern unsigned int size_above;

extern int dnest_flag_restart, dnest_flag_postprc, dnest_flag_sample_info, dnest_flag_limits;
extern double dnest_post_temp;
extern char file_restart[STR_MAX_LENGTH], file_save_restart[STR_MAX_LENGTH];

extern double post_logz;
extern int dnest_num_params;
extern char dnest_sample_postfix[STR_MAX_LENGTH], dnest_sample_tag[STR_MAX_LENGTH], dnest_sample_dir[STR_MAX_LENGTH];

//the limits of parameters for each level;
extern double *limits, *copies_of_limits;

extern int dnest_which_particle_update; // which particle to be updated
extern int dnest_which_level_update;    // which level to be updated;
extern int *dnest_perturb_accept;
extern int dnest_root;

extern void *dnest_arg;
//***********************************************
/*                  functions                  */
extern double mod(double y, double x);
extern void dnest_wrap(double *x, double min, double max);
extern void wrap_limit(double *x, double min, double max);
extern int mod_int(int y, int x);
extern int dnest_cmp(const void *pa, const void *pb);

extern void options_load(int max_num_saves, double ptol);
extern void setup(int argc, char** argv, DNestFptrSet *fptrset, int num_params, char *sample_dir, int max_num_saves, double ptol);
extern void finalise();

extern double dnest(int argc, char **argv, DNestFptrSet *fptrset,  int num_params, char *sample_dir,
             int max_num_saves, double pdff, const void *arg);
extern void dnest_run();
extern void dnest_mcmc_run();
extern void update_particle(unsigned int which);
extern void update_level_assignment(unsigned int which);
extern double log_push(unsigned int which_level);
extern bool enough_levels(Level *l, int size_l);
extern void do_bookkeeping();
extern void save_levels();
extern void save_particle();
extern void save_limits();
extern void kill_lagging_particles();
extern void renormalise_visits();
extern void recalculate_log_X();
extern double dnest_randh();
extern double dnest_rand();
extern double dnest_randn();
extern int dnest_rand_int(int size);
extern void dnest_postprocess(double temperature, int max_num_saves, double ptol);
extern void postprocess(double temperature);
extern void initialize_output_file();
extern void close_output_file();
extern void dnest_save_restart();
extern void dnest_restart();
extern void dnest_restart_action(int iflag);
extern void dnest_accept_action();
extern void dnest_kill_action(int i, int i_copy);
extern void dnest_print_particle(FILE *fp, const void *model, const void *arg);
extern void dnest_read_particle(FILE *fp, void *model);
extern int dnest_get_size_levels();
extern int dnest_get_which_level_update();
extern int dnest_get_which_particle_update();
extern void dnest_get_posterior_sample_file(char *fname);
extern int dnest_check_version(char *verion_str);
extern unsigned int dnest_get_which_num_saves();
extern unsigned int dnest_get_count_saves();
extern unsigned long long int dnest_get_count_mcmc_steps();
extern void dnest_check_fptrset(DNestFptrSet *fptrset);
extern DNestFptrSet * dnest_malloc_fptrset();
extern void dnest_free_fptrset(DNestFptrSet * fptrset);
/*=====================================================*/
// users responsible for following functions
extern void (*print_particle)(FILE *fp, const void *model, const void *arg);
extern void (*read_particle)(FILE *fp, void *model);
extern void (*from_prior)(void *model, const void *arg);
extern double (*log_likelihoods_cal)(const void *model, const void *arg);
extern double (*log_likelihoods_cal_initial)(const void *model, const void *arg);
extern double (*log_likelihoods_cal_restart)(const void *model, const void *arg);
extern double (*perturb)(void *model, const void *arg);
extern void (*restart_action)(int iflag);
extern void (*accept_action)();
extern void (*kill_action)(int i, int i_copy);
/*=====================================================*/

#ifdef __cplusplus
}
#endif

#endif

max_stars_repo_name: LiyrAstroph/PyCALI

max_issues_repo_name: LiyrAstroph/pyCALI

max_forks_repo_name: LiyrAstroph/pyCALI


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 reconstruct_line2d.c
 *  \brief reconstruct 2d line and BLR model.
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <float.h>
#include <gsl/gsl_interp.h>

#include "brains.h"

void *best_model_line2d;      /*!< best model */
void *best_model_std_line2d;  /*!< standard deviation of the best model */

/*!
 * postprocessing.
 */
void postprocess2d()
{
  char posterior_sample_file[BRAINS_MAX_STR_LENGTH];
  int num_ps, i, j, k, nc;
  double *pm, *pmstd;
  double *lag;
  void *posterior_sample, *post_model;
  double mean_lag, mean_lag_std, sum1, sum2;
  int size_of_modeltype = num_params * sizeof(double);

  best_model_line2d = malloc(size_of_modeltype);
  best_model_std_line2d = malloc(size_of_modeltype);

  if(thistask == roottask)
  {
    // initialize smoothing workspace
    smooth_init(n_vel_data_ext, Vline_data_ext);
    char fname[200];
    FILE *fp, *fcon, *fline, *ftran, *fline1d;
    double *Fline1d, dV;

    // velocity grid width, in term of wavelength of Hbeta.
    dV = (Vline_data[n_vel_data-1]-Vline_data[0])/(n_vel_data-1) * parset.linecenter/C_Unit;
    Fline1d = malloc(n_line_data * sizeof(double));

    // get number of lines in posterior sample file
    get_posterior_sample_file(dnest_options_file, posterior_sample_file);

    //file for posterior sample
    fp = fopen(posterior_sample_file, "r");
    if(fp == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", posterior_sample_file);
      exit(0);
    }
    //file for continuum reconstruction
    sprintf(fname, "%s/%s", parset.file_dir, "data/con_rec.txt");
    fcon = fopen(fname, "w");
    if(fcon == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
      exit(0);
    }
    //file for line reconstruction
    sprintf(fname, "%s/%s", parset.file_dir, "data/line2d_rec.txt");
    fline = fopen(fname, "w");
    if(fline == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s.\n", fname);
      exit(0);
    }

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

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

    // read number of lines
    if(fscanf(fp, "# %d", &num_ps) < 1)
    {
      fprintf(stderr, "# Error: Cannot read file %s.\n", posterior_sample_file);
      exit(0);
    }
    printf("# Number of points in posterior sample: %d\n", num_ps);

    lag = malloc(num_ps * sizeof(double));
    post_model = malloc(size_of_modeltype);
    posterior_sample = malloc(num_ps * size_of_modeltype);

    force_update = 1;
    which_parameter_update = -1; // force to update the transfer function
    which_particle_update = 0;
    mean_lag = 0.0;
    nc = 0;

    Fcon_rm = Fcon_rm_particles[which_particle_update];
    TransTau = TransTau_particles[which_particle_update];
    Trans2D_at_veldata = Trans2D_at_veldata_particles[which_particle_update];
    Fline2d_at_data = Fline_at_data_particles[which_particle_update];

    for(i=0; i<num_ps; i++)
    {
      // read lines
      for(j=0; j<num_params; j++)
      {
        if(fscanf(fp, "%lf", (double *)post_model + j) < 1)
        {
          fprintf(stderr, "# Error: Cannot read file %s.\n", posterior_sample_file);
          exit(0);
        }
      }
      fscanf(fp, "\n");

      //store model
      memcpy(posterior_sample+i*size_of_modeltype, post_model, size_of_modeltype);

      //calculate_con_from_model(post_model + num_params_blr *sizeof(double));
      calculate_con_from_model_semiseparable(post_model + num_params_blr *sizeof(double));
      calculate_con_rm(post_model);
      gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);

      transfun_2d_cal(post_model, Vline_data_ext, Trans2D_at_veldata,
                                      n_vel_data_ext, 0);
      calculate_line2d_from_blrmodel(post_model, Tline_data, Vline_data_ext, Trans2D_at_veldata,
                                      Fline2d_at_data, n_line_data, n_vel_data_ext);

      // calculate integrated line fluxes
      for(j = 0; j < n_line_data; j++)
      {
        Fline1d[j] = 0.0;
        for(k=0; k<n_vel_data; k++)
        {
          Fline1d[j] += Fline2d_at_data[j * n_vel_data_ext + (k+n_vel_data_incr)] * dV /line_scale;
        }
      }

      // calculate mean time lag
      sum1 = 0.0;
      sum2 = 0.0;
      for(j=0; j<parset.n_tau; j++)
      {
        for(k=0; k<n_vel_data_ext; k++)
        {
          sum1 += Trans2D_at_veldata[j * n_vel_data_ext + k] * TransTau[j];
          sum2 += Trans2D_at_veldata[j * n_vel_data_ext + k];
        }
      }
      // take care of zero transfer function
      if(sum2 > 0.0)
      {
        lag[i] = sum1/sum2;
        mean_lag += lag[i];
        nc++;
      }
      else
      {
        lag[i] = -DBL_MAX;
      }

      //if( i % (num_ps/10+1) == 0)
      {
        // output continuum
        for(j=0; j<parset.n_con_recon; j++)
        {
          fprintf(fcon, "%e %e\n", Tcon[j]*(1.0+parset.redshift), Fcon[j]/con_scale);
        }
        fprintf(fcon, "\n");

        // output 2d line
        for(j=0; j<n_line_data; j++)
        {
          for(k=0; k<n_vel_data; k++)
          {
            fprintf(fline, "%e ", Fline2d_at_data[j * n_vel_data_ext + (k+n_vel_data_incr)]/line_scale);
          }
          fprintf(fline, "\n");
        }
        fprintf(fline, "\n");

        // output transfer function
        for(j=0; j<parset.n_tau; j++)
        {
          fprintf(ftran, "%e ", TransTau[j]);

          for(k=0; k<n_vel_data; k++)
          {
            fprintf(ftran, "%e ", Trans2D_at_veldata[j * n_vel_data_ext + (k+n_vel_data_incr)]);
          }
          fprintf(ftran, "\n");
        }
        fprintf(ftran, "\n");

        // output 1d line
        for(j = 0; j<n_line_data; j++)
        {
          fprintf(fline1d, "%e %e\n", Tline_data[j]*(1.0+parset.redshift), Fline1d[j]);
        }
        fprintf(fline1d, "\n");
      }
    }

    smooth_end();

    fclose(fp);
    fclose(fcon);
    fclose(fline);
    fclose(ftran);
    fclose(fline1d);

    /* calculate mean lags */
    mean_lag /= (nc);
    mean_lag_std = 0.0;
    for(i=0; i<num_ps; i++)
    {
      if(lag[i] > -DBL_MAX)
        mean_lag_std += (lag[i] - mean_lag) * (lag[i] - mean_lag);
    }
    if(nc > 1)
      mean_lag_std = sqrt(mean_lag_std/(nc -1.0));
    else
      mean_lag_std = 0.0;
    printf("Mean time lag: %f+-%f\n", mean_lag, mean_lag_std);

    pm = (double *)best_model_line2d;
    pmstd = (double *)best_model_std_line2d;
    for(j=0; j<num_params; j++)
    {
      pm[j] = pmstd[j] = 0.0;
    }
    for(i=0; i<num_ps; i++)
    {
      for(j =0; j<num_params; j++)
        pm[j] += *((double *)posterior_sample + i*num_params + j );
    }

    for(j=0; j<num_params; j++)
      pm[j] /= num_ps;

    for(i=0; i<num_ps; i++)
    {
      for(j=0; j<num_params; j++)
        pmstd[j] += pow( *((double *)posterior_sample + i*num_params + j ) - pm[j], 2.0 );
    }

    for(j=0; j<num_params; j++)
    {
      if(num_ps > 1)
        pmstd[j] = sqrt(pmstd[j]/(num_ps-1.0));
      else
        pmstd[j] = 0.0;
    }

    for(j = 0; j<num_params_blr + num_params_var; j++)
      printf("Best params %d %f +- %f\n", j, *((double *)best_model_line2d + j),
                                             *((double *)best_model_std_line2d+j) );

    free(lag);
    free(post_model);
    free(posterior_sample);
    free(Fline1d);
  }
  return;
}

/*!
 * this function run dnest sampleing, reconstruct light curves using the best estimates for parameters.
 */
void reconstruct_line2d()
{
  int i, argc=0;
  char **argv;

  //configure restart of dnest
  argv = malloc(9*sizeof(char *));
  for(i=0; i<9; i++)
  {
    argv[i] = malloc(BRAINS_MAX_STR_LENGTH*sizeof(char));
  }
  //setup argc and argv
  strcpy(argv[argc++], "dnest");
  strcpy(argv[argc++], "-s");
  strcpy(argv[argc], parset.file_dir);
  strcat(argv[argc++], "/data/restart2d_dnest.txt");

  if(parset.flag_restart == 1)
  {
    strcpy(argv[argc++], "-r");
    strcpy(argv[argc], parset.file_dir);
    strcat(argv[argc], "/");
    strcat(argv[argc++], "data/restart2d_dnest.txt");
  }
  if(parset.flag_postprc == 1)
  {
    strcpy(argv[argc++], "-p");
  }
  if(parset.flag_temp == 1)
  {
    sprintf(argv[argc++], "-t%f", parset.temperature);
  }
  if(parset.flag_sample_info == 1)
  {
    strcpy(argv[argc++], "-c");
  }

  //level-dependent sampling
  {
    strcpy(argv[argc++], "-l");
  }

  reconstruct_line2d_init();

  smooth_init(n_vel_data_ext, Vline_data_ext);
  dnest_line2d(argc, argv);
  smooth_end();

  if(parset.flag_exam_prior != 1 && parset.flag_para_name != 1)
  {
    postprocess2d();

    // calculate light curves using the best model
    if(thistask == roottask)
    {
      force_update = 1;
      which_parameter_update = -1;
      which_particle_update = 0;
      Fcon_rm = Fcon_rm_particles[which_particle_update];
      TransTau = TransTau_particles[which_particle_update];
      Trans2D_at_veldata = Trans2D_at_veldata_particles[which_particle_update];
      Fline2d_at_data = Fline_at_data_particles[which_particle_update];

      //calculate_con_from_model(best_model_line2d + num_params_blr *sizeof(double));
      calculate_con_from_model_semiseparable(best_model_line2d + num_params_blr *sizeof(double));
      calculate_con_rm(best_model_line2d);
      gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);

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

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

      for(i=0; i<parset.n_con_recon; i++)
      {
        fprintf(fp, "%e %e %e\n", Tcon[i]*(1.0+parset.redshift), Fcon[i] / con_scale, Fcerrs[i] / con_scale);
      }
      fclose(fp);

      smooth_init(n_vel_data_ext, Vline_data_ext);
      // recovered line2d at data points
      transfun_2d_cal(best_model_line2d, Vline_data_ext, Trans2D_at_veldata,
                                              n_vel_data_ext, parset.flag_save_clouds);
      calculate_line2d_from_blrmodel(best_model_line2d, Tline_data, Vline_data_ext, Trans2D_at_veldata,
                                                       Fline2d_at_data, n_line_data, n_vel_data_ext);

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

      fprintf(fp, "# %d %d\n", n_line_data, n_vel_data);
      for(i=0; i<n_line_data; i++)
      {
        fprintf(fp, "# %f\n", Tline_data[i]*(1.0+parset.redshift));
        for(j=0; j<n_vel_data; j++)
        {
          fprintf(fp, "%e %e\n", Wline_data[j],
                            Fline2d_at_data[i*n_vel_data_ext + (j+n_vel_data_incr)] / line_scale);
        }
        fprintf(fp, "\n");
      }
      fclose(fp);

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

      fprintf(fp, "# %d %d\n", parset.n_tau, n_vel_data);
      for(i=0; i<parset.n_tau; i++)
      {
        fprintf(fp, "# %f\n", TransTau[i]);
        for(j=0; j<n_vel_data; j++)
        {
          fprintf(fp, "%e %e\n", Vline_data[j], Trans2D_at_veldata[i*n_vel_data_ext + (j+n_vel_data_incr)]);
        }
        fprintf(fp, "\n");
      }
      fclose(fp);

      smooth_end();

      // recovered line2d at specified points
      smooth_init(parset.n_vel_recon, TransV);

      which_parameter_update = -1;
      which_particle_update = 0;
      transfun_2d_cal(best_model_line2d, TransV, Trans2D, parset.n_vel_recon, 0);

      /* there is no data for spectral broadening at given specified epoch, using the mean value
       * and set InstRes_err=0.0.
       */
      double *pm = (double *)best_model_line2d;
      if(parset.flag_InstRes > 0)
      {
        parset.flag_InstRes = 0; /* force to be uniform prior */
        double instres_mean = 0.0;
        for(i=0; i<n_line_data; i++)
        {
          instres_mean += instres_epoch[i];
        }
        parset.InstRes = instres_mean/n_line_data;
        parset.InstRes_err = 0.0;

        instres_mean = 0.0;
        for(i=0; i<n_line_data; i++)
        {
          instres_mean += pm[num_params_blr_model+num_params_nlr+i];
        }
        pm[num_params_blr_model + num_params_nlr ] = instres_mean/n_line_data;
      }

      /* similarly, there is no data for line center information at given specified epoch,
       * using the mean value
       */
      if(parset.flag_linecenter < 0)
      {
        parset.flag_linecenter = 1; /* force to be uniform prior, note num_params_linecenter is still unchanged */
        double linecenter_mean = 0.0;
        for(i=0; i<n_line_data; i++)
        {
          linecenter_mean += pm[idx_linecenter + i];
        }
        pm[idx_linecenter] = linecenter_mean/n_line_data;
      }

      calculate_line2d_from_blrmodel(best_model_line2d, Tline, TransV,
          Trans2D, Fline2d, parset.n_line_recon, parset.n_vel_recon);

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

      fprintf(fp, "# %d %d\n", parset.n_line_recon, parset.n_vel_recon);
      for(i=0; i<parset.n_line_recon; i++)
      {
        fprintf(fp, "# %f\n", Tline[i]*(1.0+parset.redshift));
        for(j=0; j<parset.n_vel_recon; j++)
        {
          fprintf(fp, "%e %e\n", (1.0 + TransV[j]/C_Unit) * parset.linecenter * (1.0+parset.redshift),
                Fline2d[i*parset.n_vel_recon + j] / line_scale);
        }

        fprintf(fp, "\n");
      }
      fclose(fp);

      // output 2d transfer function
      sprintf(fname, "%s/%s", parset.file_dir, parset.tran2d_out_file);
      fp = fopen(fname, "w");
      if(fp == NULL)
      {
        fprintf(stderr, "# Error: Cannot open file %s\n", fname);
        exit(-1);
      }
      fprintf(fp, "# %d %d\n", parset.n_tau, parset.n_vel_recon);
      for(i=0; i<parset.n_tau; i++)
      {
        fprintf(fp, "# %f\n", TransTau[i]);
        for(j=0; j<parset.n_vel_recon; j++)
        {
          fprintf(fp, "%e %e\n", TransV[j]*VelUnit, Trans2D[i*parset.n_vel_recon + j]);
        }

        fprintf(fp, "\n");
      }
      fclose(fp);
      smooth_end();
    }
  }

  reconstruct_line2d_end();

  //clear up argv
  for(i=0; i<9; i++)
  {
    free(argv[i]);
  }
  free(argv);

  return;
}

/*!
 * this function initializes 2d reconstruction.
 */
void reconstruct_line2d_init()
{
  int i, j;
  double dT, Tspan;

  Tspan = Tcon_data[n_con_data -1] - Tcon_data[0];

  /* set time grid for continuum */
  /*if(parset.time_back > 0.0)
    Tcon_min = Tcon_data[0] - parset.time_back;
  else
    Tcon_min = Tcon_data[0] - fmax(0.05*Tspan, Tspan/2.0 + (Tcon_data[0] - Tline_data[0]));*/

  Tcon_min = Tcon_data[0] - time_back_set - 10.0;

  Tcon_max = Tcon_data[n_con_data-1] + fmax(0.05*Tspan, 20.0);
  Tcon_max = fmax(Tcon_max, Tline_data[n_line_data -1] + 10.0);  /* The time span should cover that of the emission line data */

  if(thistask == roottask)
    printf("Tcon_min_max: %f %f\n", Tcon_min - Tcon_data[0], Tcon_max - Tcon_data[n_con_data-1]);

  dT = (Tcon_max - Tcon_min)/(parset.n_con_recon -1);

  for(i=0; i<parset.n_con_recon; i++)
  {
    Tcon[i] = Tcon_min + i*dT;
  }

  /* set Larr_rec */
  for(i=0;i<parset.n_con_recon;i++)
  {
    Larr_rec[i*nq + 0]=1.0;
    for(j=1; j<nq; j++)
      Larr_rec[i*nq + j] = pow(Tcon[i], j);
  }

  //TransTau = malloc(parset.n_tau * sizeof(double));
  //Trans2D_at_veldata = malloc(parset.n_tau * n_vel_data * sizeof(double));

  TransV = malloc(parset.n_vel_recon * sizeof(double));
  Trans2D = malloc(parset.n_tau * parset.n_vel_recon * sizeof(double));

  //Fline2d_at_data = malloc(n_line_data * n_vel_data * sizeof(double));

  Tline = malloc(parset.n_line_recon * sizeof(double));
  Fline2d = malloc(parset.n_line_recon * parset.n_vel_recon * sizeof(double));

  Tline_min = Tline_data[0] - fmin(0.1*(Tline_data[n_line_data - 1] - Tline_data[0]), 10.0);
  if(parset.time_back <= 0.0)
    Tline_min = fmax(Tline_min, Tcon_min + time_back_set);

  Tline_max = Tline_data[n_line_data -1] + fmin(0.1*(Tline_data[n_line_data - 1] - Tline_data[0]), 10.0);
  Tline_max = fmin(Tline_max, Tcon_max - 1.0);  /* The time span should be smaller than that of the continuum */

  if(thistask == roottask)
    printf("Tline_min_max: %f %f\n", Tline_min - Tline_data[0], Tline_max - Tline_data[n_line_data -1]);

  dT = (Tline_max - Tline_min)/(parset.n_line_recon - 1);

  for(i=0; i<parset.n_line_recon; i++)
  {
    Tline[i] = Tline_min + i*dT;
  }

  double vel_max_set = Vline_data_ext[n_vel_data_ext -1], vel_min_set = Vline_data_ext[0];
  double dVel = (vel_max_set- vel_min_set)/(parset.n_vel_recon -1.0);

  for(i=0; i<parset.n_vel_recon; i++)
  {
  	TransV[i] = vel_min_set + dVel*i;
  }

  sprintf(dnest_options_file, "%s/%s", parset.file_dir, "src/OPTIONS2D");
  if(thistask == roottask)
  {
    get_num_particles(dnest_options_file);
  }
  MPI_Bcast(&parset.num_particles, 1, MPI_INT, roottask, MPI_COMM_WORLD);

  Fcon = malloc(parset.n_con_recon * sizeof(double *));
  Fcon_rm_particles = malloc(parset.num_particles * sizeof(double *));
  Fcon_rm_particles_perturb = malloc(parset.num_particles * sizeof(double *));
  for(i=0; i<parset.num_particles; i++)
  {
    Fcon_rm_particles[i] = malloc(parset.n_con_recon * sizeof(double));
    Fcon_rm_particles_perturb[i] = malloc(parset.n_con_recon * sizeof(double));
  }

  TransTau_particles = malloc(parset.num_particles * sizeof(double *));
  TransTau_particles_perturb = malloc(parset.num_particles * sizeof(double *));
  for(i=0; i<parset.num_particles; i++)
  {
    TransTau_particles[i] = malloc(parset.n_tau * sizeof(double));
    TransTau_particles_perturb[i] = malloc(parset.n_tau * sizeof(double));
  }

  Trans2D_at_veldata_particles = malloc(parset.num_particles * sizeof(double *));
  Trans2D_at_veldata_particles_perturb = malloc(parset.num_particles * sizeof(double *));
  for(i=0; i<parset.num_particles; i++)
  {
    Trans2D_at_veldata_particles[i] = malloc(parset.n_tau * n_vel_data_ext * sizeof(double));
    Trans2D_at_veldata_particles_perturb[i] = malloc(parset.n_tau * n_vel_data_ext * sizeof(double));
  }

  Fline_at_data_particles = malloc(parset.num_particles * sizeof(double *));
  Fline_at_data_particles_perturb = malloc(parset.num_particles * sizeof(double *));
  for(i=0; i<parset.num_particles; i++)
  {
    Fline_at_data_particles[i] = malloc(n_line_data * n_vel_data_ext * sizeof(double));
    Fline_at_data_particles_perturb[i] = malloc(n_line_data * n_vel_data_ext * sizeof(double));
  }

  clouds_tau = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_weight = malloc(parset.n_cloud_per_task * sizeof(double));
  clouds_vel = malloc(parset.n_cloud_per_task * parset.n_vel_per_cloud * sizeof(double));

  if(parset.flag_save_clouds && thistask == roottask)
  {
    if(parset.n_cloud_per_task <= 10000)
      icr_cloud_save = 1;
    else
      icr_cloud_save = parset.n_cloud_per_task/10000;

    char fname[200];
    sprintf(fname, "%s/%s", parset.file_dir, parset.cloud_out_file);
    fcloud_out = fopen(fname, "w");
    if(fcloud_out == NULL)
    {
      fprintf(stderr, "# Error: Cannot open file %s\n", fname);
      exit(-1);
    }
  }
  return;
}

/*!
 * this function finalizes 2d reconstruction.
 */
void reconstruct_line2d_end()
{
  free(Tline);
  //free(Fline2d_at_data);
  free(Fline2d);

  //free(TransTau);
  free(TransV);
  free(Trans2D);
  //free(Trans2D_at_veldata);
  int i;
  for(i=0; i<parset.num_particles; i++)
  {
    free(Fcon_rm_particles[i]);
    free(Fcon_rm_particles_perturb[i]);
  }
  free(Fcon);
  free(Fcon_rm_particles);
  free(Fcon_rm_particles_perturb);

  free(Fline_at_data);

  free(par_fix);
  free(par_fix_val);

  free(best_model_line2d);
  free(best_model_std_line2d);

  for(i=0; i<parset.num_particles; i++)
  {
    free(Trans2D_at_veldata_particles[i]);
    free(Trans2D_at_veldata_particles_perturb[i]);
    free(TransTau_particles[i]);
    free(TransTau_particles_perturb[i]);
    free(Fline_at_data_particles[i]);
    free(Fline_at_data_particles_perturb[i]);
  }
  free(Trans2D_at_veldata_particles);
  free(Trans2D_at_veldata_particles_perturb);
  free(TransTau_particles);
  free(TransTau_particles_perturb);
  free(Fline_at_data_particles);
  free(Fline_at_data_particles_perturb);

  for(i=0; i<num_params; i++)
  {
    free(par_range_model[i]);
    free(par_prior_gaussian[i]);
  }
  free(par_range_model);
  free(par_prior_gaussian);
  free(par_prior_model);

  free(clouds_tau);
  free(clouds_weight);
  free(clouds_vel);

  if(parset.flag_save_clouds && thistask==roottask)
  {
    fclose(fcloud_out);
  }

  if(thistask == roottask)
  {
    printf("Ends reconstruct_line2d.\n");
  }
  return;
}

/*!
 * this function calculate probability at initial step.
 */
double prob_initial_line2d(const void *model)
{
  double prob_line = 0.0, var2, dy, var2_se;
  int i, j;
  double *pm = (double *)model;

  which_particle_update = dnest_get_which_particle_update();

  Fcon_rm = Fcon_rm_particles[which_particle_update];
  //calculate_con_from_model(model + num_params_blr*sizeof(double));
  calculate_con_from_model_semiseparable(model + num_params_blr*sizeof(double));
  calculate_con_rm(model);
  gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);

  TransTau = TransTau_particles[which_particle_update];
  Trans2D_at_veldata = Trans2D_at_veldata_particles[which_particle_update];
  Fline2d_at_data = Fline_at_data_particles[which_particle_update];
  which_parameter_update = -1;
  transfun_2d_cal(model, Vline_data_ext, Trans2D_at_veldata, n_vel_data_ext, 0);
  calculate_line2d_from_blrmodel(model, Tline_data, Vline_data_ext, Trans2D_at_veldata, Fline2d_at_data, n_line_data, n_vel_data_ext);

  var2_se = (exp(pm[num_params_blr-1])-1.0) * (exp(pm[num_params_blr-1])-1.0) * line_error_mean_sq;
  for(i=0; i<n_line_data; i++)
  {
    for(j=0; j<n_vel_data; j++)
    {
      //note mask with error < 0.0
      if(Flerrs2d_data[i] > 0.0)
      {
        dy = Fline2d_data[i*n_vel_data + j] - Fline2d_at_data[i * n_vel_data_ext + (j+n_vel_data_incr)];
        var2 = Flerrs2d_data[i*n_vel_data+j]*Flerrs2d_data[i*n_vel_data+j] + var2_se;
        prob_line += (-0.5 * (dy*dy)/var2) - 0.5*log(var2 * 2.0*PI);
      }
    }
  }

  return prob_line;
}

/*!
 * this function calculate probability at restart step.
 */
double prob_restart_line2d(const void *model)
{
  double prob_line = 0.0, var2, dy, var2_se;
  int i, j;
  double *pm = (double *)model;

  which_particle_update = dnest_get_which_particle_update();
  Fcon_rm = Fcon_rm_particles[which_particle_update];
  //calculate_con_from_model(model + num_params_blr*sizeof(double));
  calculate_con_from_model_semiseparable(model + num_params_blr*sizeof(double));
  calculate_con_rm(model);
  gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);

  TransTau = TransTau_particles[which_particle_update];
  Trans2D_at_veldata = Trans2D_at_veldata_particles[which_particle_update];
  Fline2d_at_data = Fline_at_data_particles[which_particle_update];
  which_parameter_update = num_params + 1; // so as not to update clouds.
  transfun_2d_cal(model, Vline_data_ext, Trans2D_at_veldata, n_vel_data_ext, 0);
  calculate_line2d_from_blrmodel(model, Tline_data, Vline_data_ext, Trans2D_at_veldata, Fline2d_at_data, n_line_data, n_vel_data_ext);

  var2_se = (exp(pm[num_params_blr-1])-1.0) * (exp(pm[num_params_blr-1])-1.0) * line_error_mean_sq;
  for(i=0; i<n_line_data; i++)
  {
    for(j=0; j<n_vel_data; j++)
    {
      //note mask with error < 0.0
      if(Flerrs2d_data[i] > 0.0)
      {
        dy = Fline2d_data[i*n_vel_data + j] - Fline2d_at_data[i * n_vel_data_ext + (j+n_vel_data_incr)];
        var2 = Flerrs2d_data[i*n_vel_data+j]*Flerrs2d_data[i*n_vel_data+j] + var2_se;
        prob_line += (-0.5 * (dy*dy)/var2) - 0.5*log(var2 * 2.0*PI);
      }
    }
  }

  return prob_line;
}
/*!
 * this function calculate probability.
 *
 * At each MCMC step, only one parameter is updated, which only changes some values; thus,
 * optimization that reuses the unchanged values can improve computation efficiency.
 */
double prob_line2d(const void *model)
{
  double prob_line = 0.0, var2, dy, var2_se;
  int i, j;
  double *pm = (double *)model;

  which_particle_update = dnest_get_which_particle_update();

  // only update continuum reconstruction when the corresponding parameters are updated
  if(which_parameter_update >= num_params_blr)
  {
    Fcon_rm = Fcon_rm_particles_perturb[which_particle_update];
    //calculate_con_from_model(model + num_params_blr*sizeof(double));
    calculate_con_from_model_semiseparable(model + num_params_blr*sizeof(double));
    calculate_con_rm(model);
    gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);
  }
  else /* continuum has no change, use the previous values */
  {
    Fcon_rm = Fcon_rm_particles[which_particle_update];
    gsl_interp_init(gsl_linear, Tcon, Fcon_rm, parset.n_con_recon);
  }

   // only update transfer function when BLR model is changed.
   // or when forced to updated
  if( (which_parameter_update < num_params_blr-1) || force_update == 1)
  {
    TransTau = TransTau_particles_perturb[which_particle_update];
    Trans2D_at_veldata = Trans2D_at_veldata_particles_perturb[which_particle_update];
    transfun_2d_cal(model, Vline_data_ext, Trans2D_at_veldata, n_vel_data_ext, 0);
  }
  else
  {
    TransTau = TransTau_particles[which_particle_update];
    Trans2D_at_veldata = Trans2D_at_veldata_particles[which_particle_update];
  }

  /* no need to calculate line when only systematic error parameter of line are updated.
   * otherwise, always need to calculate line.
   */
  if( which_parameter_update != num_params_blr-1 || force_update == 1 )
  {
    Fline2d_at_data = Fline_at_data_particles_perturb[which_particle_update];
    calculate_line2d_from_blrmodel(model, Tline_data, Vline_data_ext, Trans2D_at_veldata, Fline2d_at_data, n_line_data, n_vel_data_ext);

    var2_se = (exp(pm[num_params_blr-1])-1.0) * (exp(pm[num_params_blr-1])-1.0) * line_error_mean_sq;
    for(i=0; i<n_line_data; i++)
    {
      for(j=0; j<n_vel_data; j++)
      {
        //note mask with error < 0.0
        if(Flerrs2d_data[i] > 0.0)
        {
          dy = Fline2d_data[i*n_vel_data + j] - Fline2d_at_data[i * n_vel_data_ext + (j+n_vel_data_incr)];
          var2 = Flerrs2d_data[i*n_vel_data+j]*Flerrs2d_data[i*n_vel_data+j] + var2_se;
          prob_line += (-0.5 * (dy*dy)/var2) - 0.5*log(var2 * 2.0*PI);
        }
      }
    }
  }
  else
  {
    /* re-point */
    Fline2d_at_data = Fline_at_data_particles[which_particle_update];
    var2_se = (exp(pm[num_params_blr-1])-1.0) * (exp(pm[num_params_blr-1])-1.0) * line_error_mean_sq;
    for(i=0; i<n_line_data; i++)
    {
      for(j=0; j<n_vel_data; j++)
      {
        //note mask with error < 0.0
        if(Flerrs2d_data[i] > 0.0)
        {
          dy = Fline2d_data[i*n_vel_data + j] - Fline2d_at_data[i * n_vel_data_ext + (j+n_vel_data_incr)];
          var2 = Flerrs2d_data[i*n_vel_data+j]*Flerrs2d_data[i*n_vel_data+j] + var2_se;
          prob_line += (-0.5 * (dy*dy)/var2) - 0.5*log(var2 * 2.0*PI);
        }
      }
    }
  }

  return prob_line;
}

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS