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 init.c
 *  \brief initialize the program.
 */

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

#include "brains.h"

/*!
 * This function initialize the program.
 */
void init()
{
  int i, j;

  nq = 1 + parset.flag_trend;
  num_params_trend = nq;
  num_params_drw = 3; /* include systematic error */
  if(parset.flag_dim == 0)
  {
    num_params_resp = 0; /* no A and Ag for 0D */
  }
  else
  {
    num_params_resp = 2; /* response A and Ag */
  }

  if(parset.flag_trend_diff > 0)
  {
    num_params_difftrend = parset.flag_trend_diff; /* differences of the trends between continuum and line */
  }
  else
  {
    num_params_difftrend = 0;
  }

  num_params_var = num_params_drw + num_params_trend + num_params_resp + num_params_difftrend;

  /* number of parameters for narrow line, only valid for 2d RM. */
  num_params_nlr = 0;
  if(parset.flag_narrowline >= 2)
    num_params_nlr = 3;

  /* number of parameters for spectral broadening, only valid for 2d RM */
  num_params_res = 1;
  if(parset.flag_InstRes > 0)
    num_params_res = n_line_data;

  /* number of parameters for line center, only valid for 2d RM */
  num_params_linecenter = 0;
  if(parset.flag_linecenter >  0)
    num_params_linecenter = 1;
  else if(parset.flag_linecenter <  0)
    num_params_linecenter = n_line_data;

  switch(parset.flag_blrmodel)
  {
    case -1:  /* user defined analytical transfer function */
      BLRmodel_size = num_params_MyTransfun2d * sizeof(double);
      set_blr_range_model = set_par_range_mytransfun;
      break;
    case 0:
      BLRmodel_size = num_params_MyBLRmodel2d * sizeof(double);
      set_blr_range_model = set_blr_range_mymodel;
      break;
    case 1:
      BLRmodel_size = sizeof(BLRmodel1);
      set_blr_range_model = set_blr_range_model1;
      break;
    case 2:
      BLRmodel_size = sizeof(BLRmodel2);
      set_blr_range_model = set_blr_range_model2;
      break;
    case 3:
      BLRmodel_size = sizeof(BLRmodel3);
      set_blr_range_model = set_blr_range_model3;
      break;
    case 4:
      BLRmodel_size = sizeof(BLRmodel4);
      set_blr_range_model = set_blr_range_model4;
      break;
    case 5:
      BLRmodel_size = sizeof(BLRmodel5);
      set_blr_range_model = set_blr_range_model5;
      break;
    case 6:
      BLRmodel_size = sizeof(BLRmodel6);
      set_blr_range_model = set_blr_range_model6;
      break;
    case 7:
      BLRmodel_size = sizeof(BLRmodel7);
      set_blr_range_model = set_blr_range_model7;
      break;
    case 8:
      BLRmodel_size = sizeof(BLRmodel8);
      set_blr_range_model = set_blr_range_model8;
      break;
    case 9:
      BLRmodel_size = sizeof(BLRmodel9);
      set_blr_range_model = set_blr_range_model9;
      break;
    default:
      BLRmodel_size = sizeof(BLRmodel1);
      set_blr_range_model = set_blr_range_model1;
      break;
  }

#ifdef SA
  if(parset.flag_dim > 2 || parset.flag_dim < 0)
  {
    num_params_sa_extpar = sizeof(SAExtPar)/sizeof(double);
    switch(parset.flag_sa_blrmodel)
    {
      case 0:
        SABLRmodel_size = num_params_MyBLRmodel2d * sizeof(double);
        set_sa_blr_range_model = set_sa_blr_range_mymodel;
        break;
      case 1:
        SABLRmodel_size = sizeof(SABLRmodel1);
        set_sa_blr_range_model = set_sa_blr_range_model1;
        break;
      case 2:
        SABLRmodel_size = sizeof(SABLRmodel2);
        set_sa_blr_range_model = set_sa_blr_range_model2;
        break;
      case 3:
        SABLRmodel_size = sizeof(SABLRmodel3);
        set_sa_blr_range_model = set_sa_blr_range_model3;
        break;
      case 4:
        SABLRmodel_size = sizeof(SABLRmodel4);
        set_sa_blr_range_model = set_sa_blr_range_model4;
        break;
      case 5:
        SABLRmodel_size = sizeof(SABLRmodel5);
        set_sa_blr_range_model = set_sa_blr_range_model5;
        break;
      case 6:
        SABLRmodel_size = sizeof(SABLRmodel6);
        set_sa_blr_range_model = set_sa_blr_range_model6;
        break;
      case 7:
        SABLRmodel_size = sizeof(SABLRmodel7);
        set_sa_blr_range_model = set_sa_blr_range_model7;
        break;
      case 8:
        SABLRmodel_size = sizeof(SABLRmodel8);
        set_sa_blr_range_model = set_sa_blr_range_model8;
        break;
      case 9:
        SABLRmodel_size = sizeof(SABLRmodel9);
        set_sa_blr_range_model = set_sa_blr_range_model9;
        break;
      default:
        SABLRmodel_size = sizeof(SABLRmodel1);
        set_sa_blr_range_model = set_sa_blr_range_model1;
        break;
    }
  }
#endif

  if(parset.flag_dim != 3)
  {
    /* set maximum continuum point */
    n_con_max = parset.n_con_recon;
    if(parset.flag_dim >=-1)
    {
      if(n_con_data > n_con_max)
        n_con_max = n_con_data;
    }
  }

  allocate_memory();

  /* initialize GSL */
  gsl_T = gsl_rng_default;
  gsl_r = gsl_rng_alloc (gsl_T);

#ifndef Debug
  if(parset.flag_rng_seed != 1)
  {
    gsl_rng_set(gsl_r, time(NULL)+thistask+1350);
  }
  else
  {
    gsl_rng_set(gsl_r, parset.rng_seed+thistask+1350);
  }
#else
  if(parset.flag_rng_seed != 1)
  {
    gsl_rng_set(gsl_r, 6666+thistask+1350);
    printf("# debugging, task %d brains random seed %d.\n", thistask, 6666+thistask+1350);
  }
  else
  {
    gsl_rng_set(gsl_r, parset.rng_seed+thistask+1350);
  }
#endif

  /* default BH mass range */
  mass_range[0] = 1.0;
  mass_range[1] = 1.0e4;

  /* default rcloud_max_set */
  rcloud_max_set = 1.0e6;

  /* not RM */
  if(parset.flag_dim != 3)
  {
    gsl_acc = gsl_interp_accel_alloc();
    gsl_linear = gsl_interp_alloc(gsl_interp_linear, parset.n_con_recon);

    if(parset.flag_dim >=-1)
    {
      /* set Larr_data */
      for(i=0;i<n_con_data;i++)
      {
        Larr_data[i*nq + 0]=1.0;
        for(j=1; j<nq; j++)
          Larr_data[i*nq + j] = pow(Tcon_data[i], j);
      }
    }

    Tcad_data = 1.0;
    Tspan_data = 1.0e4;
    if(parset.flag_dim == 0)
    {
      /* set cadence and time span of data */
      Tspan_data = (Tcon_data[n_con_data -1] - Tcon_data[0]);
      Tspan_data_con = Tspan_data;
      if(Tspan_data < 0.0)
      {
        if(thistask == roottask)
        {
          printf("# Incorrect epochs in continuum, please check the input data.\n");
          exit(0);
        }
      }
      Tcad_data = Tspan_data;
      for(i=1; i< n_con_data; i++)
      {
        if(Tcad_data > Tcon_data[i] - Tcon_data[i-1])
          Tcad_data = Tcon_data[i] - Tcon_data[i-1];
      }
    }

    if(parset.flag_dim > 0 || parset.flag_dim == -1)
    {
      /* set cadence and time span of data */
      Tspan_data_con = (Tcon_data[n_con_data -1] - Tcon_data[0]);
      Tspan_data = (Tline_data[n_line_data -1] - Tcon_data[0]);
      if(Tspan_data < 0.0)
      {
        if(thistask == roottask)
        {
          printf("# Incorrect epochs in continuum and line, please check the input data.\n");
          exit(0);
        }
      }
      Tcad_data = Tspan_data;
      for(i=1; i< n_con_data; i++)
      {
        if(Tcad_data > Tcon_data[i] - Tcon_data[i-1])
          Tcad_data = Tcon_data[i] - Tcon_data[i-1];
      }
      for(i=1; i< n_line_data; i++)
      {
        if(Tcad_data > Tline_data[i] - Tline_data[i-1])
          Tcad_data = Tline_data[i] - Tline_data[i-1];
      }

      /* set mediate time of continuum data */
      Tmed_data = 0.5*(Tcon_data[0] + Tcon_data[n_con_data-1]);

      for(i=0; i<num_params_difftrend; i++)
      {
        pow_Tcon_data[i] = (pow(Tcon_data[n_con_data-1]-Tmed_data, i+2)
                          - pow(Tcon_data[0]-Tmed_data, i+2)) / (i+2) / Tspan_data_con;
      }

      /* set time back for continuum reconstruction */
      time_back_set = Tspan_data_con + (Tcon_data[0] - Tline_data[0]);
      time_back_set = fmax(2.0*Tcad_data, time_back_set);

      /* make rcloud_max and time_back consistent with each other, rcloud_max has a higher priority */
      double DT=Tcon_data[0] - time_back_set;
      if(parset.rcloud_max > 0.0)
      {
        DT = fmax(DT, Tline_data[0] - parset.rcloud_max*2.0);
      }
      else if(parset.time_back > 0.0) /* neglect when parset.rcloud_max is set */
      {
        DT = fmax(DT, Tcon_data[0] - parset.time_back);
      }
      time_back_set = Tcon_data[0] - DT;

      /* set the range of cloud radial distribution */
      rcloud_min_set = 0.0;
      rcloud_max_set = Tspan_data/2.0;

      /* rcloud_max should smaller than  (Tl0 - Tc0)/2 */
      rcloud_max_set = fmin( rcloud_max_set,  (Tline_data[0] - Tcon_data[0] + time_back_set)/2.0 );

      if(parset.rcloud_max > 0.0)
        rcloud_max_set = fmin(rcloud_max_set, parset.rcloud_max);

      /* check whether n_con_recon is appropriate */
      double med_cad, med_cad_recon;
      //med_cad = get_mediate_cad(Tcon_data, n_con_data);
      med_cad = (Tcon_data[n_con_data-1] - Tcon_data[0])/(n_con_data -1);
      med_cad_recon = (Tcon_data[n_con_data-1] - Tcon_data[0] + time_back_set)/(parset.n_con_recon-1);
      if(med_cad_recon > 1.05*med_cad)
      {
        if(thistask == roottask)
        {
          printf("# Too small NConRecon.\n"
               "# Better to change it to %d or set RCloudMax/TimeBack to smaller than %f/%f.\n",
               (int)(parset.n_con_recon * med_cad_recon/med_cad),
               (Tline_data[0] - Tcon_data[0] + time_back_set)/2.0 ,time_back_set);
          printf("\e[1;35m" "# Use '-f' option in command line if want to ignore this check.\n" "\e[0m");
        }
        if(parset.flag_force_run != 1)
        {
          exit(0);
        }
      }
    }

    if(thistask == roottask && parset.flag_dim >=-1)
    {
      printf("rcloud_min_max_set: %f %f\n", rcloud_min_set, rcloud_max_set);
    }

    // response range A and Ag
    resp_range[0][0] = log(0.1);
    resp_range[0][1] = log(10.0);

    resp_range[1][0] = -1.0;
    resp_range[1][1] =  3.0;

    /* set the range of continuum variation  */
    var_range_model[0][0] = log(1.0); /* systematic error in continuum */
    var_range_model[0][1] = log(1.0+10.0);

    var_range_model[1][0] = -15.0; /* log(sigma) */
    var_range_model[1][1] = -1.0;

    var_range_model[2][0] = log(Tcad_data); /* log(tau) */
    var_range_model[2][1] = log(Tspan_data);

    var_range_model[3][0] = -10.0; /* mean value or trend parameter values */
    var_range_model[3][1] =  10.0;

    for(i=0; i<num_params_difftrend; i++)
    {
      var_range_model[4+i][0] = -1.0/pow(Tspan_data, i+1); /* slope of the trend in the differences between contiuum and line */
      var_range_model[4+i][1] =  1.0/pow(Tspan_data, i+1);
    }

    var_range_model[4+num_params_difftrend][0] = -10.0; /* light curve values */
    var_range_model[4+num_params_difftrend][1] =  10.0;

    if(num_params_nlr > 1)
    {
      if(parset.flag_narrowline == 2) /* Gaussian prior */
      {
        nlr_range_model[0][0] = -10.0;
        nlr_range_model[0][1] =  10.0;

        nlr_prior_model[0] = GAUSSIAN;
      }
      else
      {
        nlr_range_model[0][0] = log(1.0e-1); /* logrithmic prior */
        nlr_range_model[0][1] = log(1.0e4);

        nlr_prior_model[0] = UNIFORM;
      }

      nlr_range_model[1][0] = -10.0;
      nlr_range_model[1][1] =  10.0;

      nlr_prior_model[1] = GAUSSIAN;

      nlr_range_model[2][0] = -10.0;
      nlr_range_model[2][1] =  10.0;

      nlr_prior_model[2] = GAUSSIAN;
    }

    // range for systematic error
    sys_err_line_range[0] = log(1.0);
    sys_err_line_range[1] = log(1.0+10.0);

    set_blr_range_model();

  }

#ifdef SA
  /* SA */
  if(parset.flag_dim > 2)
  {
    set_sa_blr_range_model();

    sa_extpar_range[0][0] = log(100.0);  /* DA */
    sa_extpar_range[0][1] = log(1000.0);

    sa_extpar_range[1][0] = 0.0;         /* PA */
    sa_extpar_range[1][1] = 180.0;

    sa_extpar_range[2][0] = log(0.5);    /* FA  */
    sa_extpar_range[2][1] = log(2.0);

    sa_extpar_range[3][0] = -(wave_sa_data[1]-wave_sa_data[0]);    /* CO  */
    sa_extpar_range[3][1] =  (wave_sa_data[1]-wave_sa_data[0]);
  }
#endif
}

/*!
 * This function allocates memory for variables used throughout the code.
 */
void allocate_memory()
{
  int i;

  if(parset.flag_dim != 3)
  {
    Tcon = malloc(parset.n_con_recon * sizeof(double));
    Fcerrs = malloc(parset.n_con_recon * sizeof(double));

    PSmat = malloc(parset.n_con_recon * parset.n_con_recon * sizeof(double));
    //PNmat = malloc(parset.n_con_recon * parset.n_con_recon * sizeof(double));
    USmat = malloc(parset.n_con_recon * n_con_data * sizeof(double));
    PSmat_data = malloc(n_con_data * n_con_data * sizeof(double));
    //PNmat_data = malloc(n_con_data * n_con_data * sizeof(double));
    PCmat_data = malloc(n_con_data * n_con_data * sizeof(double));
    IPCmat_data = malloc(n_con_data * n_con_data * sizeof(double));
    PQmat = malloc(parset.n_con_recon * parset.n_con_recon * sizeof(double));
    PEmat1 = malloc(parset.n_con_recon * n_con_data * sizeof(double));
    PEmat2 = malloc(parset.n_con_recon * parset.n_con_recon * sizeof(double));

    blr_range_model = malloc(BLRmodel_size/sizeof(double) * sizeof(double *));
    for(i=0; i<BLRmodel_size/sizeof(double); i++)
    {
      blr_range_model[i] = malloc(2*sizeof(double));
    }

    workspace = malloc((n_con_data*nq + 7*n_con_max + nq*nq + nq)*sizeof(double));
    workspace_uv = malloc(2*parset.n_con_recon*sizeof(double));
    Larr_data = malloc(n_con_data*nq*sizeof(double));
    Larr_rec = malloc(parset.n_con_recon*nq*sizeof(double));
    pow_Tcon_data = malloc(num_params_difftrend*sizeof(double));

    var_param = malloc(num_params_var * sizeof(double));
    var_param_std = malloc(num_params_var * sizeof(double));

    for(i=0; i<num_params_var; i++)
    {
      var_param[i] = var_param_std[i] = 0.0;
    }
  }

#ifdef SA
  if(parset.flag_dim > 2 || parset.flag_dim < 0)
  {
    sa_extpar_range = malloc(num_params_sa_extpar * sizeof(double *));
    for(i=0; i<num_params_sa_extpar; i++)
    {
      sa_extpar_range[i]=malloc(2*sizeof(double));
    }

    sa_blr_range_model = malloc(SABLRmodel_size/sizeof(double) * sizeof(double *));
    for(i=0; i<SABLRmodel_size/sizeof(double); i++)
    {
      sa_blr_range_model[i] = malloc(2*sizeof(double));
    }

    if(parset.flag_sa_par_mutual != 0)
    {
      idx_sa_par_mutual = malloc(2*sizeof(int));
      idx_rm_par_mutual = malloc(2*sizeof(int));
    }
  }
#endif

  return;
}

/*!
 * This function free memory.
 */
void free_memory()
{
  int i;

  if(parset.flag_dim != 3)
  {
    free(Tcon);
    free(Fcerrs);

    free(PSmat);
    //(PNmat);
    free(USmat);
    free(PSmat_data);
    //free(PNmat_data);
    free(PCmat_data);
    free(IPCmat_data);
    free(PQmat);
    free(PEmat1);
    free(PEmat2);

    for(i=0; i<BLRmodel_size/sizeof(double); i++)
    {
      free(blr_range_model[i]);
    }
    free(blr_range_model);

    free(workspace);
    free(workspace_uv);
    free(Larr_data);
    free(Larr_rec);
    free(pow_Tcon_data);

    free(var_param);
    free(var_param_std);
  }

#ifdef SA
  if(parset.flag_dim > 2 || parset.flag_dim < 0)
  {
    for(i=0; i<num_params_sa_extpar; i++)
    {
      free(sa_extpar_range[i]);
    }
    free(sa_extpar_range);

    for(i=0; i<SABLRmodel_size/sizeof(double); i++)
    {
      free(sa_blr_range_model[i]);
    }
    free(sa_blr_range_model);

    if(parset.flag_sa_par_mutual != 0)
    {
      free(idx_sa_par_mutual);
      free(idx_rm_par_mutual);
    }
  }
#endif
  return;
}

/*!
 * This function normalise the light curves to a scale of unity.
 */
void scale_con_line()
{
  int i, j;
  double ave_con, ave_line;

  con_scale = 1.0;
  line_scale = 1.0;

  ave_con = 0.0;
  ave_line = 0.0;
  for(i=0; i<n_con_data; i++)
  {
    ave_con += Fcon_data[i];
  }

  ave_con /= n_con_data;
  con_scale = 1.0/ave_con;

  for(i=0; i<n_con_data; i++)
  {
    Fcon_data[i] *=con_scale;
    Fcerrs_data[i] *=con_scale;
  }

  if(thistask == roottask)
    printf("con scale: %e\t%e\n", con_scale, ave_con);

  con_error_mean *= con_scale;

  if(parset.flag_dim == 0)
  {
    return;
  }

  for(i=0; i<n_line_data; i++)
  {
    ave_line += Fline_data[i];
  }
  ave_line /=n_line_data;

  line_scale = 1.0/ave_line;

  if(thistask == roottask)
    printf("line scale: %e\t%e\n", line_scale, ave_line);

  for(i=0; i<n_line_data; i++)
  {
    // note mask with error < 0.0
    if(Flerrs_data[i] > 0.0)
    {
      Fline_data[i] *= line_scale;
      Flerrs_data[i] *= line_scale;
    }

    if(parset.flag_dim==2 || parset.flag_dim == -1 || parset.flag_dim == 5)
      for(j=0; j<n_vel_data; j++)
      {
        // note mask with error < 0.0
        if(Flerrs2d_data[i*n_vel_data + j] > 0.0)
        {
          Fline2d_data[i*n_vel_data + j] *= line_scale;
          Flerrs2d_data[i*n_vel_data + j] *= line_scale;
        }
      }
  }
  line_error_mean *= line_scale;

  if( (parset.flag_dim==2 || parset.flag_dim == -1 || parset.flag_dim == 5) && parset.flag_narrowline!=0)
  {
    parset.flux_narrowline *= line_scale;
    parset.flux_narrowline_err *= line_scale;
  }
  return;
}

/*!
 * This function copes with parameter fixing.\n
 * Only fix BLR model parameters.
 */
void set_par_fix_blrmodel()
{
  int i;
  char *pstr;

  npar_fix = 0;

  if(thistask == roottask)
  {
    pstr = parset.str_par_fix_val;
    // set the default value if not provided.
    for(i=strlen(parset.str_par_fix); i<num_params_blr_model; i++)
      parset.str_par_fix[i] = '0';

    for(i=0; i<num_params_blr_model; i++)
    {
      if(parset.str_par_fix[i] == '0')
      {
        par_fix[i] = 0;
        par_fix_val[i] = -DBL_MAX;  /* set to be the smallest value */
      }
      else if(parset.str_par_fix[i] == '1')
      {
        if(pstr == NULL)
        {
          printf("# %d-th BLR parameter value is not provided (counting from 0).\n", i);
          exit(0);
        }
        par_fix[i] = 1;
        sscanf(pstr, "%lf", &par_fix_val[i]);
        npar_fix++;
        printf("# %d-th BLR parameter fixed, value= %f.\n", i, par_fix_val[i]);
        pstr = strchr(pstr, ':'); /* values are separated by ":" */
        if(pstr!=NULL)
        {
          pstr++;
        }
      }
      else   // default value
      {
        par_fix[i] = 0;
        par_fix_val[i] = -DBL_MAX;
      }
    }
  }

  MPI_Bcast(par_fix, num_params_blr_model, MPI_INT, roottask, MPI_COMM_WORLD);
  MPI_Bcast(par_fix_val, num_params_blr_model, MPI_DOUBLE, roottask, MPI_COMM_WORLD);

  return;
}

/*
 * setup BLR model parameter range.
 */

// model 1
void set_blr_range_model1()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 0.999;
  //inc
  blr_range_model[i][0] = 0.0;  // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //lambda
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.5;
  //q
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;

  //blr_range_model[2][1] = fmin(blr_range_model[2][1], log(rcloud_max_set));
  return;
}

// model 2
void set_blr_range_model2()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 0.999;
  //inc
  blr_range_model[i][0] = 0.0;  // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //sigr
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //sigtheta
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;

  return;
}

// model 3
void set_blr_range_model3()
{
  int i;

  i = 0;
  //Rin
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //F
  blr_range_model[i][0] = log(1.0);
  blr_range_model[i++][1] = log(1.0e2);
  //alpha
  blr_range_model[i][0] = -3.0;
  blr_range_model[i++][1] = 3.0;
  //inc
  blr_range_model[i][0] = 0.0;  // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;   // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //q
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;

  //rcloud_max_set = fmax(rcloud_max_set, exp(blr_range_model[3][1] + blr_range_model[4][1]));
  return;
}

// model 4
void set_blr_range_model4()
{
  int i;

  i = 0;
  //Rin
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //F
  blr_range_model[i][0] = log(1.0);
  blr_range_model[i++][1] = log(1.0e2);
  //alpha
  blr_range_model[i][0] = -3.0;
  blr_range_model[i++][1] = 3.0;
  //inc
  blr_range_model[i][0] = 0.0; // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = sqrt(2.0)/2.0;
  //q
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;

  //rcloud_max_set = fmax(rcloud_max_set, exp(blr_range_model[3][1] + blr_range_model[4][1]));
  return;
}

// model 5
void set_blr_range_model5()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //Fin
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //Fout
  blr_range_model[i][0] = log(1.0);
  blr_range_model[i++][1] = log(10.0);
  //alpha
  blr_range_model[i][0] = 1.0;
  blr_range_model[i++][1] = 3.0;
  //inc
  blr_range_model[i][0] = 0.0; // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in degree
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //beta
  blr_range_model[i][0] = 1.0;
  blr_range_model[i++][1] = 5.0;
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //fellip
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //fflow
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //sigr_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //theta_rot
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 90.0;
  //sig_turb
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);

  //rcloud_max_set = fmax(rcloud_max_set, exp(blr_range_model[2][1] + blr_range_model[4][1]));

  return;
}

// model 6
void set_blr_range_model6()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //inc
  blr_range_model[i][0] = 0.0; //in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //gamma
  blr_range_model[i][0] = 1.0;
  blr_range_model[i++][1] = 5.0;
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //fellip
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //fflow
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //sigr_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //theta_rot
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 90.0;
  //sig_turb
  blr_range_model[i][0] = log(0.0001);
  blr_range_model[i++][1] = log(0.1);

  //rcloud_max_set = 44.85;
  //blr_range_model[2][1] = fmin(blr_range_model[2][1], log(rcloud_max_set));
  return;
}

// model 7
void set_blr_range_model7()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //inc
  blr_range_model[i][0] = 0.0;  // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //gamma
  blr_range_model[i][0] = 1.0;
  blr_range_model[i++][1] = 5.0;
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;

  //fsh
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //mu_un
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta_un
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F_un
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //opn_un
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 90.0;

  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);
  //fellip
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //fflow
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //sigr_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_circ
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //sigr_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);
  //sigthe_rad
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(1.0);
  //theta_rot
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 90.0;

  //fellip_un
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //fflow_un
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;

  //sig_turb
  blr_range_model[i][0] = log(0.001);
  blr_range_model[i++][1] = log(0.1);

  return;
}

// model 8
void set_blr_range_model8()
{
  int i;

  i = 0;
  //theta_min
  blr_range_model[i][0] = 20.0;     // in degree
  blr_range_model[i++][1] = 90.0;
  //dtheta_max
  blr_range_model[i][0] = 0.0;      // in degree
  blr_range_model[i++][1] = 90.0;
  //r_min
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.25);
  //fr_max
  blr_range_model[i][0] = 1.0;  //
  blr_range_model[i++][1] = 10.0;
  //gamma
  blr_range_model[i][0] = 0.0;  //
  blr_range_model[i++][1] = 3.0;
  //alpha
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 3.0;
  //lambda
  blr_range_model[i][0] = -3.0;
  blr_range_model[i++][1] = 0.0;
  //k
  blr_range_model[i][0] = -0.5;
  blr_range_model[i++][1] = 0.5;
  //xi
  blr_range_model[i][0] = 0.0;
  blr_range_model[i++][1] = 1.0;
  //Rv
  blr_range_model[i][0] = log(10.0);
  blr_range_model[i++][1] = log(50.0);
  //Rblr
  blr_range_model[i][0] = log(10.0);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //inc
  blr_range_model[i][0] = 0.0;   // cos(inc)
  blr_range_model[i++][1] = 1.0;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);

  return;
}

// model 9
void set_blr_range_model9()
{
  int i;

  i = 0;
  //mu
  blr_range_model[i][0] = log(0.1);
  blr_range_model[i++][1] = log(rcloud_max_set*0.5);
  //beta
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 2.0;
  //F
  blr_range_model[i][0] = 0.001;
  blr_range_model[i++][1] = 0.999;
  //inc
  blr_range_model[i][0] = 0.0;  // in cosine
  blr_range_model[i++][1] = 1.0;
  //opn
  blr_range_model[i][0] = 0.0;  // in rad
  blr_range_model[i++][1] = 90.0;
  //mbh
  blr_range_model[i][0] = log(mass_range[0]);
  blr_range_model[i++][1] = log(mass_range[1]);

  return;
}

/*
 * set 1D BLR model and functions.
 */
void set_blr_model1d()
{
  switch(parset.flag_blrmodel)
  {
    case -1:
      num_params_blr_model = num_params_MyTransfun1d;
      transfun_1d_cal = transfun_1d_cal_mytransfun;
      break;

    case 0:
      num_params_blr_model = num_params_MyBLRmodel1d;
      gen_cloud_sample = gen_cloud_sample_mymodel;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 1:
      num_params_blr_model = 6;
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 2:
      num_params_blr_model = 6;
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 3:
      num_params_blr_model = 6;
      gen_cloud_sample = gen_cloud_sample_model3;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 4:
      num_params_blr_model = 6;
      gen_cloud_sample = gen_cloud_sample_model3;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 5:
      num_params_blr_model = 9;
      gen_cloud_sample = gen_cloud_sample_model5;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 6:
      num_params_blr_model = 8;
      gen_cloud_sample = gen_cloud_sample_model6;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 7:
      num_params_blr_model = 13;
      gen_cloud_sample = gen_cloud_sample_model7;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 8:
      num_params_blr_model = 13;
      gen_cloud_sample = gen_cloud_sample_model8;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    case 9:
      num_params_blr_model = 5;
      gen_cloud_sample = gen_cloud_sample_model9;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;

    default:
      num_params_blr_model = 6;
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_1d_cal = transfun_1d_cal_cloud;
      break;
  }
  return;
}

/*
 * set 2D BLR model and functions.
 */
void set_blr_model2d()
{
  switch(parset.flag_blrmodel)
  {
    case -1:
      num_params_blr_model = num_params_MyTransfun2d;
      transfun_2d_cal = transfun_2d_cal_mytransfun;
      break;

    case 0:
      num_params_blr_model = num_params_MyBLRmodel2d;
      gen_cloud_sample = gen_cloud_sample_mymodel;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 1:
      num_params_blr_model = sizeof(BLRmodel1)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 2:
      num_params_blr_model = sizeof(BLRmodel2)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model2;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 3:
      num_params_blr_model = sizeof(BLRmodel3)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model3;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 4:
      num_params_blr_model = sizeof(BLRmodel4)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model4;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 5:
      num_params_blr_model = sizeof(BLRmodel5)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model5;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 6:
      num_params_blr_model = sizeof(BLRmodel6)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model6;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 7:
      num_params_blr_model = sizeof(BLRmodel7)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model7;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 8:
      num_params_blr_model = sizeof(BLRmodel8)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model8;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    case 9:
      num_params_blr_model = sizeof(BLRmodel9)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model9;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;

    default:
      num_params_blr_model = sizeof(BLRmodel1)/sizeof(double);
      gen_cloud_sample = gen_cloud_sample_model1;
      transfun_2d_cal = transfun_2d_cal_cloud;
      break;
  }
  return;
}

max_stars_repo_name: yzxamos/BRAINS

max_issues_repo_name: yzxamos/BRAINS

max_forks_repo_name: yzxamos/BRAINS


text:
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
#include <math.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_integration.h>

#include "ccl.h"

//
// Macros for replacing relative paths
#define EXPAND_STR(s) STRING(s)
#define STRING(s) #s

const ccl_configuration default_config = {
  ccl_boltzmann_class, ccl_halofit, ccl_nobaryons,
  ccl_tinker10, ccl_duffy2008, ccl_emu_strict};

//Precision parameters
/**
 * Default relative precision if not otherwise specified
 */
#define GSL_EPSREL 1E-4

/**
 * Default number of iterations for integration and root-finding if not otherwise
 * specified
 */
#define GSL_N_ITERATION 1000

/**
 * Default number of Gauss-Kronrod points in QAG integration if not otherwise
 * specified
 */
#define GSL_INTEGRATION_GAUSS_KRONROD_POINTS GSL_INTEG_GAUSS41

/**
 * Relative precision in sigma_R calculations
 */
#define GSL_EPSREL_SIGMAR 1E-5

/**
 * Relative precision in k_NL calculations
 */
#define GSL_EPSREL_KNL 1E-5

/**
 * Relative precision in distance calculations
 */
#define GSL_EPSREL_DIST 1E-6

/**
 * Relative precision in growth calculations
 */
#define GSL_EPSREL_GROWTH 1E-6

/**
 * Relative precision in dNdz calculations
 */
#define GSL_EPSREL_DNDZ 1E-6

const ccl_gsl_params default_gsl_params = {
  GSL_N_ITERATION,                     // N_ITERATION
  GSL_INTEGRATION_GAUSS_KRONROD_POINTS,// INTEGRATION_GAUSS_KRONROD_POINTS
  GSL_EPSREL,                          // INTEGRATION_EPSREL
  GSL_INTEGRATION_GAUSS_KRONROD_POINTS,// INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS
  GSL_EPSREL,                          // INTEGRATION_LIMBER_EPSREL
  GSL_EPSREL_DIST,                     // INTEGRATION_DISTANCE_EPSREL
  GSL_EPSREL_SIGMAR,                   // INTEGRATION_SIGMAR_EPSREL
  GSL_EPSREL_KNL,                      // INTEGRATION_KNL_EPSREL
  GSL_EPSREL,                          // ROOT_EPSREL
  GSL_N_ITERATION,                     // ROOT_N_ITERATION
  GSL_EPSREL_GROWTH,                   // ODE_GROWTH_EPSREL
  1E-6,                                // EPS_SCALEFAC_GROWTH
  1E7,                                 // HM_MMIN
  1E17,                                // HM_MMAX
  0.0,                                 // HM_EPSABS
  1E-4,                                // HM_EPSREL
  1000,                                // HM_LIMIT
  GSL_INTEG_GAUSS41                    // HM_INT_METHOD
  };

#undef GSL_EPSREL
#undef GSL_N_ITERATION
#undef GSL_INTEGRATION_GAUSS_KRONROD_POINTS
#undef GSL_EPSREL_SIGMAR
#undef GSL_EPSREL_KNL
#undef GSL_EPSREL_DIST
#undef GSL_EPSREL_GROWTH
#undef GSL_EPSREL_DNDZ

const ccl_spline_params default_spline_params = {

  // scale factor spline params
  250,  // A_SPLINE_NA
  0.1,  // A_SPLINE_MIN
  0.01,  // A_SPLINE_MINLOG_PK
  0.1,  // A_SPLINE_MIN_PK,
  0.01,  // A_SPLINE_MINLOG_SM,
  0.1,  // A_SPLINE_MIN_SM,
  1.0,  // A_SPLINE_MAX,
  0.0001,  // A_SPLINE_MINLOG,
  250,  // A_SPLINE_NLOG,

  // mass splines
  0.025,  // LOGM_SPLINE_DELTA
  50,  // LOGM_SPLINE_NM
  6,  // LOGM_SPLINE_MIN
  17,  // LOGM_SPLINE_MAX

  // PS a and k spline
  13,  // A_SPLINE_NA_SM
  6,  // A_SPLINE_NLOG_SM
  40,  // A_SPLINE_NA_PK
  11,  // A_SPLINE_NLOG_PK

  // k-splines and integrals
  50,  // K_MAX_SPLINE
  1E3,  // K_MAX
  5E-5,  // K_MIN
  0.025,  // DLOGK_INTEGRATION
  167,  // N_K
  100000,  // N_K_3DCOR

  // correlation function parameters
  0.01,  // ELL_MIN_CORR
  60000,  // ELL_MAX_CORR
  5000,  // N_ELL_CORR

  //Spline types
  NULL,
  NULL,
  NULL,
  NULL,
  NULL,
  NULL,
  NULL
};

ccl_physical_constants ccl_constants = {
  /**
   * Lightspeed / H0 in units of Mpc/h (from CODATA 2014)
   */
  2997.92458,

  /**
   * Newton's gravitational constant in units of m^3/Kg/s^2
   */
  //6.6738e-11,  /(from PDG 2013) in m^3/Kg/s^2
  //6.67428e-11, // CLASS VALUE
  6.67408e-11, // from CODATA 2014

  /**
   * Solar mass in units of kg (from GSL)
   */
  //GSL_CONST_MKSA_SOLAR_MASS,
  //1.9885e30, //(from PDG 2015) in Kg
  1.9884754153381438E+30, //from IAU 2015

  /**
   * Mpc to meters (from PDG 2016 and using M_PI)
   */
  3.085677581491367399198952281E+22,

  /**
   * pc to meters (from PDG 2016 and using M_PI)
   */
  3.085677581491367399198952281E+16,

  /**
   * Rho critical in units of M_sun/h / (Mpc/h)^3
   */
  ((3*100*100)/(8*M_PI*6.67408e-11)) * (1000*1000*3.085677581491367399198952281E+22/1.9884754153381438E+30),

  /**
   * Boltzmann constant in units of J/K
  */
  //GSL_CONST_MKSA_BOLTZMANN,
  1.38064852E-23, //from CODATA 2014

  /**
   * Stefan-Boltzmann constant in units of kg/s^3 / K^4
   */
  //GSL_CONST_MKSA_STEFAN_BOLTZMANN_CONSTANT,
  5.670367E-8, //from CODATA 2014

  /**
   * Planck's constant in units kg m^2 / s
   */
  //GSL_CONST_MKSA_PLANCKS_CONSTANT_H,
  6.626070040E-34, //from CODATA 2014

  /**
   * The speed of light in m/s
   */
  //GSL_CONST_MKSA_SPEED_OF_LIGHT,
  299792458.0, //from CODATA 2014

  /**
   * Electron volt to Joules convestion
   */
  //GSL_CONST_MKSA_ELECTRON_VOLT,
  1.6021766208e-19,  //from CODATA 2014

  /**
   * Temperature of the CMB in K
   */
  2.725,
  //2.7255, // CLASS value

  /**
   * T_ncdm, as taken from CLASS, explanatory.ini
   */
  0.71611,

  /**
   * neutrino mass splitting differences
   * See Lesgourgues and Pastor, 2012 for these values.
   * Adv. High Energy Phys. 2012 (2012) 608515,
   * arXiv:1212.6154, page 13
  */
  7.62E-5,
  2.55E-3,
  -2.43E-3
};

/* ------- ROUTINE: ccl_cosmology_create ------
INPUTS: ccl_parameters params
        ccl_configuration config
TASK: creates the ccl_cosmology struct and passes some values to it
DEFINITIONS:
chi: comoving distance [Mpc]
growth: growth function (density)
fgrowth: logarithmic derivative of the growth (density) (dlnD/da?)
E: E(a)=H(a)/H0
growth0: growth at z=0, defined to be 1
sigma: ?
p_lnl: nonlinear matter power spectrum at z=0?
computed_distances, computed_growth,
computed_power, computed_sigma: store status of the computations
*/
ccl_cosmology * ccl_cosmology_create(ccl_parameters params, ccl_configuration config)
{
  ccl_cosmology * cosmo = malloc(sizeof(ccl_cosmology));
  cosmo->params = params;
  cosmo->config = config;
  cosmo->gsl_params = default_gsl_params;
  cosmo->spline_params = default_spline_params;
  cosmo->spline_params.A_SPLINE_TYPE = gsl_interp_akima;
  cosmo->spline_params.K_SPLINE_TYPE = gsl_interp_akima;
  cosmo->spline_params.M_SPLINE_TYPE = gsl_interp_akima;
  cosmo->spline_params.D_SPLINE_TYPE = gsl_interp_akima;
  cosmo->spline_params.PNL_SPLINE_TYPE = gsl_interp2d_bicubic;
  cosmo->spline_params.PLIN_SPLINE_TYPE = gsl_interp2d_bicubic;
  cosmo->spline_params.CORR_SPLINE_TYPE = gsl_interp_akima;

  cosmo->data.chi = NULL;
  cosmo->data.growth = NULL;
  cosmo->data.fgrowth = NULL;
  cosmo->data.E = NULL;
  cosmo->data.growth0 = 1.;
  cosmo->data.achi = NULL;

  cosmo->data.logsigma = NULL;

  cosmo->data.rsd_splines[0] = NULL;
  cosmo->data.rsd_splines[1] = NULL;
  cosmo->data.rsd_splines[2] = NULL;

  cosmo->computed_distances = false;
  cosmo->computed_growth = false;
  cosmo->computed_sigma = false;
  cosmo->status = 0;
  ccl_cosmology_set_status_message(cosmo, "");

  if(cosmo->spline_params.A_SPLINE_MAX !=1.) {
    cosmo->status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo,
                                     "ccl_core.c: ccl_cosmology_create(): "
                                     "A_SPLINE_MAX needs to be 1.\n");
  }

  return cosmo;
}

/* ------ ROUTINE: ccl_parameters_fill_initial -------
INPUT: ccl_parameters: params
TASK: fill parameters not set by ccl_parameters_create with some initial values
DEFINITIONS:
Omega_g = (Omega_g*h^2)/h^2 is the radiation parameter; "g" is for photons, as in CLASS
T_CMB: CMB temperature in Kelvin
Omega_l: Lambda
A_s: amplitude of the primordial PS, enforced here to initially set to NaN
sigma8: variance in 8 Mpc/h spheres for normalization of matter PS, enforced here to initially set to NaN
z_star: recombination redshift
 */
void ccl_parameters_fill_initial(ccl_parameters * params, int *status)
{
  // Fixed radiation parameters
  // Omega_g * h**2 is known from T_CMB
  params->T_CMB =  ccl_constants.T_CMB;
  // kg / m^3
  double rho_g = 4. * ccl_constants.STBOLTZ / pow(ccl_constants.CLIGHT, 3) * pow(params->T_CMB, 4);
  // kg / m^3
  double rho_crit =
    ccl_constants.RHO_CRITICAL *
    ccl_constants.SOLAR_MASS/pow(ccl_constants.MPC_TO_METER, 3) *
    pow(params->h, 2);
  params->Omega_g = rho_g/rho_crit;

  // Get the N_nu_rel from Neff and N_nu_mass
  params->N_nu_rel = params->Neff - params->N_nu_mass * pow(ccl_constants.TNCDM, 4) / pow(4./11.,4./3.);

  // Temperature of the relativistic neutrinos in K
  double T_nu= (params->T_CMB) * pow(4./11.,1./3.);
  // in kg / m^3
  double rho_nu_rel =
    params->N_nu_rel* 7.0/8.0 * 4. *
    ccl_constants.STBOLTZ / pow(ccl_constants.CLIGHT, 3) *
    pow(T_nu, 4);
  params-> Omega_nu_rel = rho_nu_rel/rho_crit;

  // If non-relativistic neutrinos are present, calculate the phase_space integral.
  if((params->N_nu_mass)>0) {
    params->Omega_nu_mass = ccl_Omeganuh2(
      1.0, params->N_nu_mass, params->m_nu, params->T_CMB, status) / ((params->h)*(params->h));
  }
  else{
    params->Omega_nu_mass = 0.;
  }

  params->Omega_m = params->Omega_b + params-> Omega_c + params->Omega_nu_mass;
  params->Omega_l = 1.0 - params->Omega_m - params->Omega_g - params->Omega_nu_rel - params->Omega_k;
  // Initially undetermined parameters - set to nan to trigger
  // problems if they are mistakenly used.
  if (isfinite(params->A_s)) {params->sigma8 = NAN;}
  if (isfinite(params->sigma8)) {params->A_s = NAN;}
  params->z_star = NAN;

  if(fabs(params->Omega_k)<1E-6)
    params->k_sign=0;
  else if(params->Omega_k>0)
    params->k_sign=-1;
  else
    params->k_sign=1;
  params->sqrtk=sqrt(fabs(params->Omega_k))*params->h/ccl_constants.CLIGHT_HMPC;
}

/* ------ ROUTINE: ccl_parameters_create -------
INPUT: numbers for the basic cosmological parameters needed by CCL
TASK: fill params with some initial values provided by the user
DEFINITIONS:
Omega_c: cold dark matter
Omega_b: baryons
Omega_m: matter
Omega_k: curvature
little omega_x means Omega_x*h^2
Neff : Effective number of neutrino speces
mnu : Pointer to either sum of neutrino masses or list of three masses.
mnu_type : how the neutrino mass(es) should be treated
w0: Dark energy eq of state parameter
wa: Dark energy eq of state parameter, time variation
H0: Hubble's constant in km/s/Mpc.
h: Hubble's constant divided by (100 km/s/Mpc).
A_s: amplitude of the primordial PS
n_s: index of the primordial PS

 */
ccl_parameters ccl_parameters_create(double Omega_c, double Omega_b, double Omega_k,
				     double Neff, double* mnu, int n_mnu,
				     double w0, double wa, double h, double norm_pk,
				     double n_s, double bcm_log10Mc, double bcm_etab,
				     double bcm_ks, double mu_0, double sigma_0,
				     int nz_mgrowth, double *zarr_mgrowth,
				     double *dfarr_mgrowth, int *status)
{
  #ifndef USE_GSL_ERROR
    gsl_set_error_handler_off();
  #endif

  ccl_parameters params;
  // Initialize params
  params.m_nu = NULL;
  params.z_mgrowth=NULL;
  params.df_mgrowth=NULL;
  params.sigma8 = NAN;
  params.A_s = NAN;
  params.Omega_c = Omega_c;
  params.Omega_b = Omega_b;
  params.Omega_k = Omega_k;
  params.Neff = Neff;
  params.m_nu = malloc(n_mnu*sizeof(double));
  params.sum_nu_masses = 0.;
  for(int i = 0; i<n_mnu; i=i+1){
     params.m_nu[i] = mnu[i];
     params.sum_nu_masses = params.sum_nu_masses + mnu[i];
  }

  if(params.sum_nu_masses<1e-15){
    params.N_nu_mass = 0;
  }else{
    params.N_nu_mass = n_mnu;
   }

  // Dark Energy
  params.w0 = w0;
  params.wa = wa;

  // Hubble parameters
  params.h = h;
  params.H0 = h*100;

  // Primordial power spectra
  if(norm_pk<1E-5)
    params.A_s=norm_pk;
  else
    params.sigma8=norm_pk;
  params.n_s = n_s;

  //Baryonic params
  if(bcm_log10Mc<0)
    params.bcm_log10Mc=log10(1.2e14);
  else
    params.bcm_log10Mc=bcm_log10Mc;
  if(bcm_etab<0)
    params.bcm_etab=0.5;
  else
    params.bcm_etab=bcm_etab;
  if(bcm_ks<0)
    params.bcm_ks=55.0;
  else
    params.bcm_ks=bcm_ks;

  // Params of the mu / Sigma parameterisation of MG
  params.mu_0 = mu_0;
  params.sigma_0 = sigma_0;

  // Set remaining standard and easily derived parameters
  ccl_parameters_fill_initial(&params, status);

  //Trigger modified growth function if nz>0
  if(nz_mgrowth>0) {
    params.has_mgrowth=true;
    params.nz_mgrowth=nz_mgrowth;
    params.z_mgrowth=malloc(params.nz_mgrowth*sizeof(double));
    params.df_mgrowth=malloc(params.nz_mgrowth*sizeof(double));
    memcpy(params.z_mgrowth,zarr_mgrowth,params.nz_mgrowth*sizeof(double));
    memcpy(params.df_mgrowth,dfarr_mgrowth,params.nz_mgrowth*sizeof(double));
  }
  else {
    params.has_mgrowth=false;
    params.nz_mgrowth=0;
    params.z_mgrowth=NULL;
    params.df_mgrowth=NULL;
  }

  return params;
}

/**
 * Write a cosmology parameters object to a file in yaml format.
 * @param cosmo Cosmological parameters
 * @param f FILE* pointer opened for reading
 * @return void
 */
void ccl_parameters_write_yaml(ccl_parameters * params, const char * filename, int *status)
{
  FILE * f = fopen(filename, "w");

  if (!f){
    *status = CCL_ERROR_FILE_WRITE;
    return;
  }

#define WRITE_DOUBLE(name) fprintf(f, #name ": %le\n",params->name)
#define WRITE_INT(name) fprintf(f, #name ": %d\n",params->name)

  // Densities: CDM, baryons, total matter, curvature
  WRITE_DOUBLE(Omega_c);
  WRITE_DOUBLE(Omega_b);
  WRITE_DOUBLE(Omega_m);
  WRITE_DOUBLE(Omega_k);
  WRITE_INT(k_sign);

  // Dark Energy
  WRITE_DOUBLE(w0);
  WRITE_DOUBLE(wa);

  // Hubble parameters
  WRITE_DOUBLE(H0);
  WRITE_DOUBLE(h);

  // Neutrino properties
  WRITE_DOUBLE(Neff);
  WRITE_INT(N_nu_mass);
  WRITE_DOUBLE(N_nu_rel);

  if (params->N_nu_mass>0){
    fprintf(f, "m_nu: [");
    for (int i=0; i<params->N_nu_mass; i++){
      fprintf(f, "%le, ", params->m_nu[i]);
    }
    fprintf(f, "]\n");
  }

  WRITE_DOUBLE(sum_nu_masses);
  WRITE_DOUBLE(Omega_nu_mass);
  WRITE_DOUBLE(Omega_nu_rel);

  // Primordial power spectra
  WRITE_DOUBLE(A_s);
  WRITE_DOUBLE(n_s);

  // Radiation parameters
  WRITE_DOUBLE(Omega_g);
  WRITE_DOUBLE(T_CMB);

  // BCM baryonic model parameters
  WRITE_DOUBLE(bcm_log10Mc);
  WRITE_DOUBLE(bcm_etab);
  WRITE_DOUBLE(bcm_ks);

  // Modified gravity parameters
  WRITE_DOUBLE(mu_0);
  WRITE_DOUBLE(sigma_0);

  // Derived parameters
  WRITE_DOUBLE(sigma8);
  WRITE_DOUBLE(Omega_l);
  WRITE_DOUBLE(z_star);

  WRITE_INT(has_mgrowth);
  WRITE_INT(nz_mgrowth);

  if (params->has_mgrowth){
    fprintf(f, "z_mgrowth: [");
    for (int i=0; i<params->nz_mgrowth; i++){
      fprintf(f, "%le, ", params->z_mgrowth[i]);
    }
    fprintf(f, "]\n");

    fprintf(f, "df_mgrowth: [");
    for (int i=0; i<params->nz_mgrowth; i++){
      fprintf(f, "%le, ", params->df_mgrowth[i]);
    }
    fprintf(f, "]\n");
  }

#undef WRITE_DOUBLE
#undef WRITE_INT

  fclose(f);
}

/**
 * Write a cosmology parameters object to a file in yaml format.
 * @param cosmo Cosmological parameters
 * @param f FILE* pointer opened for reading
 * @return void
 */
ccl_parameters ccl_parameters_read_yaml(const char * filename, int *status) {
  FILE * f = fopen(filename, "r");

  if (!f) {
    *status = CCL_ERROR_FILE_READ;
    ccl_parameters bad_params;
    ccl_raise_warning(CCL_ERROR_FILE_READ,
                      "ccl_core.c: ccl_parameters_read_yaml(): "
                      "Failed to read parameters from file.");

    return bad_params;
  }

#define READ_DOUBLE(name) double name; *status |= (0==fscanf(f, #name ": %le\n",&name));
#define READ_INT(name) int name; *status |= (0==fscanf(f, #name ": %d\n",&name))

  // Densities: CDM, baryons, total matter, curvature
  READ_DOUBLE(Omega_c);
  READ_DOUBLE(Omega_b);
  READ_DOUBLE(Omega_m);
  READ_DOUBLE(Omega_k);
  READ_INT(k_sign);

  // Dark Energy
  READ_DOUBLE(w0);
  READ_DOUBLE(wa);

  // Hubble parameters
  READ_DOUBLE(H0);
  READ_DOUBLE(h);

  // Neutrino properties
  READ_DOUBLE(Neff);
  READ_INT(N_nu_mass);
  READ_DOUBLE(N_nu_rel);

  double mnu[3] = {0.0, 0.0, 0.0};
  if (N_nu_mass>0){
    *status |= (0==fscanf(f, "m_nu: ["));
    for (int i=0; i<N_nu_mass; i++){
      *status |= (0==fscanf(f, "%le, ", mnu+i));
    }
    *status |= (0==fscanf(f, "]\n"));
  }

  READ_DOUBLE(sum_nu_masses);
  READ_DOUBLE(Omega_nu_mass);
  READ_DOUBLE(Omega_nu_rel);

  // Primordial power spectra
  READ_DOUBLE(A_s);
  READ_DOUBLE(n_s);

  // Radiation parameters
  READ_DOUBLE(Omega_g);
  READ_DOUBLE(T_CMB);

  // BCM baryonic model parameters
  READ_DOUBLE(bcm_log10Mc);
  READ_DOUBLE(bcm_etab);
  READ_DOUBLE(bcm_ks);

  // Modified gravity parameters
  READ_DOUBLE(mu_0);
  READ_DOUBLE(sigma_0);

  // Derived parameters
  READ_DOUBLE(sigma8);
  READ_DOUBLE(Omega_l);
  READ_DOUBLE(z_star);

  READ_INT(has_mgrowth);
  READ_INT(nz_mgrowth);

  double *z_mgrowth;
  double *df_mgrowth;

  if (has_mgrowth){
    z_mgrowth = malloc(nz_mgrowth*sizeof(double));
    df_mgrowth = malloc(nz_mgrowth*sizeof(double));
    *status |= (0==fscanf(f, "z_mgrowth: ["));
    for (int i=0; i<nz_mgrowth; i++){
      *status |= (0==fscanf(f, "%le, ", z_mgrowth+i));
    }
    *status |= (0==fscanf(f, "]\n"));

    *status |= (0==fscanf(f, "df_mgrowth: ["));
    for (int i=0; i<nz_mgrowth; i++){
      *status |= (0==fscanf(f, "%le, ", df_mgrowth+i));
    }
    *status |= (0==fscanf(f, "]\n"));
  }
  else{
    z_mgrowth = NULL;
    df_mgrowth = NULL;
  }

#undef READ_DOUBLE
#undef READ_INT

  fclose(f);

  if (*status) {
    ccl_raise_warning(
      *status,
      "ccl_core.c: ccl_parameters_read_yaml():"
      "Structure of YAML file incorrect: %s",
      filename);
  }

  double norm_pk;

  if (isnan(A_s)){
    norm_pk = sigma8;
  }
  else{
   norm_pk = A_s;
  }

  ccl_parameters params = ccl_parameters_create(
    Omega_c, Omega_b, Omega_k,
    Neff, mnu, N_nu_mass,
    w0, wa, h, norm_pk,
    n_s, bcm_log10Mc, bcm_etab,
    bcm_ks, mu_0, sigma_0, nz_mgrowth, z_mgrowth,
    df_mgrowth, status);

  if(z_mgrowth) free(z_mgrowth);
  if (df_mgrowth) free(df_mgrowth);

  return params;
}

/* ------- ROUTINE: ccl_data_free --------
INPUT: ccl_data
TASK: free the input data
*/
void ccl_data_free(ccl_data * data) {
  //We cannot assume that all of these have been allocated
  //TODO: it would actually make more sense to do this within ccl_cosmology_free,
  //where we could make use of the flags "computed_distances" etc. to figure out
  //what to free up
  gsl_spline_free(data->chi);
  gsl_spline_free(data->growth);
  gsl_spline_free(data->fgrowth);
  gsl_spline_free(data->E);
  gsl_spline_free(data->achi);
  gsl_spline2d_free(data->logsigma);
  ccl_f1d_t_free(data->rsd_splines[0]);
  ccl_f1d_t_free(data->rsd_splines[1]);
  ccl_f1d_t_free(data->rsd_splines[2]);
}

/* ------- ROUTINE: ccl_cosmology_set_status_message --------
INPUT: ccl_cosmology struct, status_string
TASK: set the status message safely.
*/
void ccl_cosmology_set_status_message(ccl_cosmology * cosmo, const char * message, ...) {
  const int trunc = 480; /* must be < 500 - 4 */

  va_list va;
  va_start(va, message);

  #pragma omp critical
  {
    vsnprintf(cosmo->status_message, trunc, message, va);

    /* if truncation happens, message[trunc - 1] is not NULL, ... will show up. */
    strcpy(&cosmo->status_message[trunc], "...");
  }

  va_end(va);
}

/* ------- ROUTINE: ccl_parameters_free --------
INPUT: ccl_parameters struct
TASK: free allocated quantities in the parameters struct
*/
void ccl_parameters_free(ccl_parameters * params) {
  if (params->m_nu != NULL){
    free(params->m_nu);
    params->m_nu = NULL;
  }
  if (params->z_mgrowth != NULL){
    free(params->z_mgrowth);
    params->z_mgrowth = NULL;
  }
  if (params->df_mgrowth != NULL){
    free(params->df_mgrowth);
    params->df_mgrowth = NULL;
  }
}

/* ------- ROUTINE: ccl_cosmology_free --------
INPUT: ccl_cosmology struct
TASK: free the input data and the cosmology struct
*/
void ccl_cosmology_free(ccl_cosmology * cosmo) {
  if (cosmo != NULL)
    ccl_data_free(&cosmo->data);
  free(cosmo);
}

int ccl_get_pk_spline_na(ccl_cosmology *cosmo) {
  return cosmo->spline_params.A_SPLINE_NA_PK + cosmo->spline_params.A_SPLINE_NLOG_PK - 1;
}

void ccl_get_pk_spline_a_array(ccl_cosmology *cosmo,int ndout,double* doutput,int *status) {
  double *d = NULL;
  if (ndout != ccl_get_pk_spline_na(cosmo))
    *status = CCL_ERROR_INCONSISTENT;
  if (*status == 0) {
    d = ccl_linlog_spacing(cosmo->spline_params.A_SPLINE_MINLOG_PK,
      cosmo->spline_params.A_SPLINE_MIN_PK,
      cosmo->spline_params.A_SPLINE_MAX,
      cosmo->spline_params.A_SPLINE_NLOG_PK,
      cosmo->spline_params.A_SPLINE_NA_PK);
    if (d == NULL)
      *status = CCL_ERROR_MEMORY;
  }
  if(*status==0)
    memcpy(doutput, d, ndout*sizeof(double));

  free(d);
}

int ccl_get_pk_spline_nk(ccl_cosmology *cosmo) {
  double ndecades = log10(cosmo->spline_params.K_MAX) - log10(cosmo->spline_params.K_MIN);
  return (int)ceil(ndecades*cosmo->spline_params.N_K);
}

void ccl_get_pk_spline_lk_array(ccl_cosmology *cosmo,int ndout,double* doutput,int *status) {
  double *d = NULL;
  if (ndout != ccl_get_pk_spline_nk(cosmo))
    *status = CCL_ERROR_INCONSISTENT;
  if (*status == 0) {
    d = ccl_log_spacing(cosmo->spline_params.K_MIN, cosmo->spline_params.K_MAX, ndout);
    if (d == NULL)
      *status = CCL_ERROR_MEMORY;
  }
  if (*status == 0) {
    for(int ii=0; ii < ndout; ii++)
      doutput[ii] = log(d[ii]);
  }
  free(d);
}

ext: c

max_stars_repo_name: borisbolliet/CCL

max_issues_repo_name: borisbolliet/CCL

max_forks_repo_name: borisbolliet/CCL


text:
#include <cblas.h>
#include <stdio.h>

double m[] = {
  3, 1, 3,
  1, 5, 9,
  2, 6, 5
};

double x[] = {
  -1, 3, -3
};

#ifdef __cplusplus
extern "C" {
#endif

     void dgesv_(int *n, int *nrhs,  double *a,  int  *lda,
           int *ipivot, double *b, int *ldb, int *info);

#ifdef __cplusplus
}
#endif

int main(void) {
  int i;
  // blas:
  double A[6] = {1.0, 2.0, 1.0, -3.0, 4.0, -1.0};
  double B[6] = {1.0, 2.0, 1.0, -3.0, 4.0, -1.0};
  double C[9] = {.5, .5, .5, .5, .5, .5, .5, .5, .5};
  cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans,
              3, 3, 2, 1, A, 3, B, 3, 2, C, 3);
  for (i = 0; i < 9; i++)
    printf("%f\n", C[i]);

  // lapack:
  int ipiv[3];
  int j;
  int info;
  int n = 1;
  int nrhs = 1;
  int lda = 3;
  int ldb = 3;
  dgesv_(&n,&nrhs, &m[0], &lda, ipiv, &x[0], &ldb, &info);
  for (i=0; i<3; ++i)
    printf("%5.1f\n", x[i]);

  return 0;
}

max_stars_repo_name: kkauder/spack

max_issues_repo_name: kkauder/spack

max_forks_repo_name: kkauder/spack


text:
/* interpolation/accel.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * 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.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>

gsl_interp_accel *
gsl_interp_accel_alloc (void)
{
  gsl_interp_accel *a = (gsl_interp_accel *) malloc (sizeof (gsl_interp_accel));
  if (a == 0)
    {
      GSL_ERROR_NULL("could not allocate space for gsl_interp_accel", GSL_ENOMEM);
    }

  a->cache = 0;
  a->hit_count = 0;
  a->miss_count = 0;

  return a;
}

int
gsl_interp_accel_reset (gsl_interp_accel * a)
{
  a->cache = 0;
  a->hit_count = 0;
  a->miss_count = 0;

  return GSL_SUCCESS;
}

void
gsl_interp_accel_free (gsl_interp_accel * a)
{
  RETURN_IF_NULL (a);
  free (a);
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#include <qdm.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <hdf5.h>
#include <hdf5_hl.h>

static
void
qdm_mcmc_workspace_init(
    qdm_mcmc *mcmc
)
{
  mcmc->w.xi      = gsl_vector_calloc(mcmc->p.xi->size);
  mcmc->w.xi_acc  = gsl_vector_calloc(mcmc->p.xi->size);
  mcmc->w.xi_p    = gsl_vector_calloc(mcmc->p.xi->size);
  mcmc->w.xi_tune = gsl_vector_calloc(mcmc->p.xi->size);

  gsl_vector_memcpy(mcmc->w.xi,   mcmc->p.xi);
  gsl_vector_memcpy(mcmc->w.xi_p, mcmc->p.xi);

  gsl_vector_set_all(mcmc->w.xi_tune, mcmc->p.xi_tune_sd);

  mcmc->w.ll   = gsl_vector_calloc(mcmc->p.ll->size);
  mcmc->w.ll_p = gsl_vector_calloc(mcmc->p.ll->size);

  gsl_vector_memcpy(mcmc->w.ll,   mcmc->p.ll);
  gsl_vector_memcpy(mcmc->w.ll_p, mcmc->p.ll);

  mcmc->w.tau   = gsl_vector_calloc(mcmc->p.tau->size);
  mcmc->w.tau_p = gsl_vector_calloc(mcmc->p.tau->size);

  gsl_vector_memcpy(mcmc->w.tau,   mcmc->p.tau);
  gsl_vector_memcpy(mcmc->w.tau_p, mcmc->p.tau);

  mcmc->w.theta        = gsl_matrix_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2);
  mcmc->w.theta_acc    = gsl_matrix_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2);
  mcmc->w.theta_p      = gsl_matrix_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2);
  mcmc->w.theta_star   = gsl_matrix_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2);
  mcmc->w.theta_star_p = gsl_matrix_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2);
  mcmc->w.theta_tune   = gsl_matrix_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2);

  gsl_matrix_memcpy(mcmc->w.theta,      mcmc->p.theta);
  gsl_matrix_memcpy(mcmc->w.theta_star, mcmc->p.theta);

  gsl_matrix_set_all(mcmc->w.theta_tune, mcmc->p.theta_tune_sd);
}

static
void
qdm_mcmc_workspace_fini(
    qdm_mcmc *mcmc
)
{
  gsl_vector_free(mcmc->w.xi);
  gsl_vector_free(mcmc->w.xi_acc);
  gsl_vector_free(mcmc->w.xi_p);
  gsl_vector_free(mcmc->w.xi_tune);

  gsl_vector_free(mcmc->w.ll);
  gsl_vector_free(mcmc->w.ll_p);

  gsl_vector_free(mcmc->w.tau);
  gsl_vector_free(mcmc->w.tau_p);

  gsl_matrix_free(mcmc->w.theta);
  gsl_matrix_free(mcmc->w.theta_acc);
  gsl_matrix_free(mcmc->w.theta_p);
  gsl_matrix_free(mcmc->w.theta_star);
  gsl_matrix_free(mcmc->w.theta_star_p);
  gsl_matrix_free(mcmc->w.theta_tune);

  mcmc->w.xi      = NULL;
  mcmc->w.xi_acc  = NULL;
  mcmc->w.xi_p    = NULL;
  mcmc->w.xi_tune = NULL;

  mcmc->w.ll   = NULL;
  mcmc->w.ll_p = NULL;

  mcmc->w.tau   = NULL;
  mcmc->w.tau_p = NULL;

  mcmc->w.theta        = NULL;
  mcmc->w.theta_acc    = NULL;
  mcmc->w.theta_p      = NULL;
  mcmc->w.theta_star   = NULL;
  mcmc->w.theta_star_p = NULL;
  mcmc->w.theta_tune   = NULL;
}

static
void
qdm_mcmc_results_init(
    qdm_mcmc *mcmc
)
{
  mcmc->r.s = mcmc->p.iter / mcmc->p.thin;

  mcmc->r.theta      = qdm_ijk_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2, mcmc->r.s);
  mcmc->r.theta_star = qdm_ijk_calloc(mcmc->p.theta->size1, mcmc->p.theta->size2, mcmc->r.s);

  mcmc->r.ll  = qdm_ijk_calloc(1, mcmc->p.ll->size, mcmc->r.s);
  mcmc->r.tau = qdm_ijk_calloc(1, mcmc->p.tau->size, mcmc->r.s);

  mcmc->r.xi = qdm_ijk_calloc(1, mcmc->p.xi->size, mcmc->r.s);
}

static
void
qdm_mcmc_results_fini(
    qdm_mcmc *mcmc
)
{
  qdm_ijk_free(mcmc->r.theta);
  qdm_ijk_free(mcmc->r.theta_star);

  qdm_ijk_free(mcmc->r.ll);
  qdm_ijk_free(mcmc->r.tau);

  qdm_ijk_free(mcmc->r.xi);

  mcmc->r.theta      = NULL;
  mcmc->r.theta_star = NULL;

  mcmc->r.ll  = NULL;
  mcmc->r.tau = NULL;

  mcmc->r.xi = NULL;
}

qdm_mcmc *
qdm_mcmc_alloc(
    qdm_mcmc_parameters p
)
{
  qdm_mcmc *mcmc = malloc(sizeof(qdm_mcmc));

  mcmc->p = p;

  qdm_mcmc_workspace_init(mcmc);
  qdm_mcmc_results_init(mcmc);

  return mcmc;
}

void
qdm_mcmc_free(
    qdm_mcmc *mcmc
)
{
  if (mcmc == NULL) {
    return;
  }

  qdm_mcmc_results_fini(mcmc);
  qdm_mcmc_workspace_fini(mcmc);

  free(mcmc);
}

int
qdm_mcmc_run(
    qdm_mcmc *mcmc
)
{
  int status = 0;

  for (size_t i = 0; i < mcmc->p.burn; i++) {
    status = qdm_mcmc_next(mcmc);
    if (status != 0) {
      goto cleanup;
    }

    if (i % mcmc->p.acc_check == 0) {
      qdm_mcmc_update_tune(mcmc);
    }
  }

  for (size_t i = 0; i < mcmc->p.iter; i++) {
    status = qdm_mcmc_next(mcmc);
    if (status != 0) {
      goto cleanup;
    }

    if (i % mcmc->p.thin == 0) {
      status = qdm_mcmc_save(mcmc, i / mcmc->p.thin);
      if (status != 0) {
        goto cleanup;
      }
    }
  }

cleanup:
  return status;
}

int
qdm_mcmc_next(
    qdm_mcmc *mcmc
)
{
  int status = 0;

  status = qdm_mcmc_update_theta(mcmc);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_mcmc_update_xi(mcmc);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  return status;
}

void
qdm_mcmc_update_tune(
    qdm_mcmc *mcmc
)
{
  for (size_t p = 0; p < mcmc->w.theta->size1; p++) {
    for (size_t m = 0; m < mcmc->w.theta->size2; m++) {
      if (gsl_matrix_get(mcmc->w.theta_acc, p, m) / (double)mcmc->p.acc_check > 0.5) {
        gsl_matrix_set(mcmc->w.theta_tune, p, m, gsl_min(gsl_matrix_get(mcmc->w.theta_tune, p, m) * 1.2, 5));
      }

      if (gsl_matrix_get(mcmc->w.theta_acc, p, m) / (double)mcmc->p.acc_check < 0.3) {
        gsl_matrix_set(mcmc->w.theta_tune, p, m, gsl_matrix_get(mcmc->w.theta_tune, p, m) * 0.8);
      }
    }
  }

  if (gsl_vector_get(mcmc->w.xi_acc, 0) / mcmc->p.acc_check > 0.5) {
    gsl_vector_set(mcmc->w.xi_tune, 0, gsl_min(gsl_vector_get(mcmc->w.xi_tune, 0) * 1.2, 5));
  }
  if (gsl_vector_get(mcmc->w.xi_acc, 0) / mcmc->p.acc_check < 0.3) {
    gsl_vector_set(mcmc->w.xi_tune, 0, gsl_vector_get(mcmc->w.xi_tune, 0) * 0.8);
  }

  if (gsl_vector_get(mcmc->w.xi_acc, 1) / mcmc->p.acc_check > 0.5) {
    gsl_vector_set(mcmc->w.xi_tune, 1, gsl_min(gsl_vector_get(mcmc->w.xi_tune, 1) * 1.2, 5));
  }
  if (gsl_vector_get(mcmc->w.xi_acc, 1) / mcmc->p.acc_check < 0.3) {
    gsl_vector_set(mcmc->w.xi_tune, 1, gsl_vector_get(mcmc->w.xi_tune, 1) * 0.8);
  }

  /* Reset the acceptance counter. */
  gsl_matrix_set_zero(mcmc->w.theta_acc);
  gsl_vector_set_zero(mcmc->w.xi_acc);
}

int
qdm_mcmc_update_theta(
    qdm_mcmc *mcmc
)
{
  int status = 0;

  /* Update the mth spline. */
  size_t p_max = mcmc->w.theta->size1;
  if (mcmc->p.truncate) {
    p_max = 1;
  }

  for (size_t p = 0; p < p_max; p++) {
    for (size_t m = 0; m < mcmc->w.theta->size2; m++) {
      /* Propose new theta star. */
      status = gsl_matrix_memcpy(mcmc->w.theta_star_p, mcmc->w.theta_star);
      if (status != 0) {
        goto cleanup;
      }

      double proposal =
        gsl_matrix_get(mcmc->w.theta_tune, p, m) *
        gsl_ran_ugaussian(mcmc->p.rng) +
        gsl_matrix_get(mcmc->w.theta_star, p, m);

      gsl_matrix_set(mcmc->w.theta_star_p, p, m, proposal);

      gsl_matrix_memcpy(mcmc->w.theta_p, mcmc->w.theta_star_p);
      qdm_theta_matrix_constrain(mcmc->w.theta_p, mcmc->p.theta_min);

      for (size_t i = 0; i < mcmc->p.y->size; i++) {
        qdm_logl_3(
            gsl_vector_ptr(mcmc->w.ll_p, i),
            gsl_vector_ptr(mcmc->w.tau_p, i),

            gsl_vector_get(mcmc->p.x, i),
            gsl_vector_get(mcmc->p.y, i),

            mcmc->p.t,
            mcmc->w.xi,

            mcmc->w.theta_p
        );
      }

      double ratio = qdm_vector_sum(mcmc->w.ll_p) - qdm_vector_sum(mcmc->w.ll);

      if (log(gsl_rng_uniform(mcmc->p.rng)) < ratio) {
        status = gsl_matrix_memcpy(mcmc->w.theta_star, mcmc->w.theta_star_p);
        if (status != 0) {
          goto cleanup;
        }

        status = gsl_matrix_memcpy(mcmc->w.theta, mcmc->w.theta_p);
        if (status != 0) {
          goto cleanup;
        }

        status = gsl_vector_memcpy(mcmc->w.ll, mcmc->w.ll_p);
        if (status != 0) {
          goto cleanup;
        }

        status = gsl_vector_memcpy(mcmc->w.tau, mcmc->w.tau_p);
        if (status != 0) {
          goto cleanup;
        }

        gsl_matrix_set(mcmc->w.theta_acc, p, m, gsl_matrix_get(mcmc->w.theta_acc, p, m) + 1);
      }
    }
  }

cleanup:
  return status;
}

int
qdm_mcmc_update_xi(
    qdm_mcmc *mcmc
)
{
  int status = 0;

  /* Update xi_low. */
  {
    double xi_low = gsl_vector_get(mcmc->w.xi, 0);
    double xi_low_p = exp(log(xi_low) + gsl_vector_get(mcmc->w.xi_tune, 0) * gsl_ran_ugaussian(mcmc->p.rng));

    gsl_vector_memcpy(mcmc->w.xi_p, mcmc->w.xi);
    gsl_vector_set(mcmc->w.xi_p, 0, xi_low_p);

    for (size_t i = 0; i < mcmc->p.y->size; i++) {
      qdm_logl_3(
          gsl_vector_ptr(mcmc->w.ll_p, i),
          gsl_vector_ptr(mcmc->w.tau_p, i),

          gsl_vector_get(mcmc->p.x, i),
          gsl_vector_get(mcmc->p.y, i),

          mcmc->p.t,
          mcmc->w.xi_p,

          mcmc->w.theta
      );
    }

    double ratio = -0.5 * (1 / mcmc->p.xi_prior_var) *  pow(log(xi_low_p) - mcmc->p.xi_prior_mean, 2) +
                    0.5 * (1 / mcmc->p.xi_prior_var) * (pow(log(xi_low  ) - mcmc->p.xi_prior_mean, 2) + qdm_vector_sum(mcmc->w.ll_p) - qdm_vector_sum(mcmc->w.ll));
    if (log(gsl_rng_uniform(mcmc->p.rng)) < ratio) {
      gsl_vector_set(mcmc->w.xi, 0, xi_low_p);

      status = gsl_vector_memcpy(mcmc->w.ll, mcmc->w.ll_p);
      if (status != 0) {
        goto cleanup;
      }

      gsl_vector_set(mcmc->w.xi_acc, 0, gsl_vector_get(mcmc->w.xi_acc, 0) + 1);
    }
  }

  /* Update xi_high. */
  {
    double xi_high = gsl_vector_get(mcmc->w.xi, 1);
    double xi_high_p = exp(log(xi_high) + gsl_vector_get(mcmc->w.xi_tune, 1) * gsl_ran_ugaussian(mcmc->p.rng));

    gsl_vector_memcpy(mcmc->w.xi_p, mcmc->w.xi);
    gsl_vector_set(mcmc->w.xi_p, 1, xi_high_p);

    for (size_t i = 0; i < mcmc->p.y->size; i++) {
      qdm_logl_3(
          gsl_vector_ptr(mcmc->w.ll_p, i),
          gsl_vector_ptr(mcmc->w.tau_p, i),

          gsl_vector_get(mcmc->p.x, i),
          gsl_vector_get(mcmc->p.y, i),

          mcmc->p.t,
          mcmc->w.xi_p,

          mcmc->w.theta
      );
    }

    double ratio = -0.5 * (1 / mcmc->p.xi_prior_var) *  pow(log(xi_high_p) - mcmc->p.xi_prior_mean, 2) +
                    0.5 * (1 / mcmc->p.xi_prior_var) * (pow(log(xi_high  ) - mcmc->p.xi_prior_mean, 2) + qdm_vector_sum(mcmc->w.ll_p) - qdm_vector_sum(mcmc->w.ll));
    if (log(gsl_rng_uniform(mcmc->p.rng)) < ratio) {
      gsl_vector_set(mcmc->w.xi, 1, xi_high_p);

      status = gsl_vector_memcpy(mcmc->w.ll, mcmc->w.ll_p);
      if (status != 0) {
        goto cleanup;
      }

      gsl_vector_set(mcmc->w.xi_acc, 1, gsl_vector_get(mcmc->w.xi_acc, 1) + 1);
    }
  }

cleanup:
  return status;
}

int
qdm_mcmc_save(
    qdm_mcmc *mcmc,
    size_t k
)
{
  int status = 0;

  {
    gsl_matrix_view m = qdm_ijk_get_ij(mcmc->r.theta, k);

    status = gsl_matrix_memcpy(&m.matrix, mcmc->w.theta);
    if (status != 0) {
      goto cleanup;
    }
  }

  {
    gsl_matrix_view m = qdm_ijk_get_ij(mcmc->r.theta_star, k);

    status = gsl_matrix_memcpy(&m.matrix, mcmc->w.theta_star);
    if (status != 0) {
      goto cleanup;
    }
  }

  {
    gsl_matrix_view m = qdm_ijk_get_ij(mcmc->r.ll, k);

    status = gsl_matrix_set_row(&m.matrix, 0, mcmc->w.ll);
    if (status != 0) {
      goto cleanup;
    }
  }

  {
    gsl_matrix_view m = qdm_ijk_get_ij(mcmc->r.tau, k);

    status = gsl_matrix_set_row(&m.matrix, 0, mcmc->w.tau);
    if (status != 0) {
      goto cleanup;
    }
  }

  {
    gsl_matrix_view m = qdm_ijk_get_ij(mcmc->r.xi, k);

    status = gsl_matrix_set_row(&m.matrix, 0, mcmc->w.xi);
    if (status != 0) {
      goto cleanup;
    }
  }

cleanup:
  return status;
}

int
qdm_mcmc_write(
    hid_t id,
    const qdm_mcmc *mcmc
)
{
  int status = 0;

  status = qdm_ijk_write(id, "theta", mcmc->r.theta);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_write(id, "theta_star", mcmc->r.theta_star);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_write(id, "ll", mcmc->r.ll);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_write(id, "tau", mcmc->r.tau);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_write(id, "xi", mcmc->r.xi);
  if (status != 0) {
    goto cleanup;
  }

cleanup:
  return status;
}

int
qdm_mcmc_read(
    hid_t id,
    qdm_mcmc **mcmc
)
{
  int status = 0;

  *mcmc = malloc(sizeof(qdm_mcmc));

  status = qdm_ijk_read(id, "theta", &(*mcmc)->r.theta);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_read(id, "theta_star", &(*mcmc)->r.theta_star);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_read(id, "ll", &(*mcmc)->r.ll);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_read(id, "tau", &(*mcmc)->r.tau);
  if (status != 0) {
    goto cleanup;
  }

  status = qdm_ijk_read(id, "xi", &(*mcmc)->r.xi);
  if (status != 0) {
    goto cleanup;
  }

  (*mcmc)->r.s = (*mcmc)->r.theta->size3;

cleanup:
  return status;
}

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>

#define RESOLUTION 28  // bmp image 28*28
#define LINE_LEN (RESOLUTION * RESOLUTION * 4) // All the pixel data stored within a line
#define DATA_NUM_WITHOUT_LABEL (RESOLUTION * RESOLUTION) // Data length without label
#define DATA_NUM_WITH_LABEL ((RESOLUTION * RESOLUTION) + 1) // For training, the first number is label, so +1
#define HIDDEN_NODE 250 // node number in hidden layer
#define OUTPUT_NODE 10  // output node number
#define LEARN_RATE 0.1
#define TRAIN_SET_SIZE 5000
#define TEST_SET_SIZE 1000
#define EPOCH 8

#define SIGMOID(x) (1.0 / (1.0 + (pow(M_E, -x))))
#define NOR_DIS_DEV(x) (pow(x, -0.5))  // For generating init weight

FILE* openfile(char* file_name);
void read_parse(FILE* fp, int input_data[], int* label, char* mode);
gsl_matrix* scale_input(int input_data[]);
gsl_matrix* vectorize_target(int label);
gsl_matrix* init_weight(int row, int col);
double get_output_error(gsl_matrix* target, gsl_matrix* output, gsl_matrix* output_error);
void training(gsl_matrix* better_input, gsl_matrix* target, gsl_matrix* input_weight, gsl_matrix* hidden_weight, gsl_matrix* hidden_x, gsl_matrix* hidden_o, gsl_matrix* output_o, gsl_matrix* output_x, gsl_matrix* output_error, gsl_matrix* hidden_error, double* error);
void store_weight(gsl_matrix* input_weight, gsl_matrix* hidden_weight);
void load_weight(gsl_matrix* input_weight, gsl_matrix* hidden_weight);
int recognize(gsl_matrix* better_input, gsl_matrix* input_weight, gsl_matrix* hidden_weight, gsl_matrix* hidden_x, gsl_matrix* hidden_o, gsl_matrix* output_x, gsl_matrix* output_o);

int main(int argc, char **argv)
{
	int label;
	double error;
	double error_acc = 0.0;
	char* file_name;
	char* mode;
	FILE* fp;
	int input_data[DATA_NUM_WITHOUT_LABEL];
	gsl_matrix* better_input;
	gsl_matrix* target;
	gsl_matrix* output_error = gsl_matrix_alloc(OUTPUT_NODE, 1);
	gsl_matrix* hidden_error = gsl_matrix_alloc(HIDDEN_NODE, 1);
	gsl_matrix* output_x = gsl_matrix_alloc(OUTPUT_NODE, 1);
	gsl_matrix* output_o = gsl_matrix_alloc(OUTPUT_NODE, 1);
    gsl_matrix* hidden_x = gsl_matrix_alloc(HIDDEN_NODE, 1);
    gsl_matrix* hidden_o = gsl_matrix_alloc(HIDDEN_NODE, 1);
	gsl_matrix* input_weight;
	gsl_matrix* hidden_weight;

	if(argc < 3){
		fprintf(stderr, "Usage: %s [-t(raining), -r(cognizing), -e(poch)] [file_name]\n", argv[0]);
		exit(EXIT_FAILURE);
	}else{
		mode = argv[1];
		file_name = argv[2];
	}

	fp = openfile(file_name);

	if(!strcmp(mode, "-t")){
		//training
		input_weight = init_weight(HIDDEN_NODE, DATA_NUM_WITHOUT_LABEL); // initiate the input-hidden weight
    	hidden_weight = init_weight(OUTPUT_NODE, HIDDEN_NODE);  // initiate the hidden-output weight

    	int training_times = 0;
    	puts("Training started\n");
        for (int i = 0; i < EPOCH; i++){

			do{
				read_parse(fp, input_data, &label, mode);
				better_input = scale_input(input_data);
				target = vectorize_target(label);
				training(better_input, target, input_weight, hidden_weight, hidden_x, hidden_o, output_o, output_x, output_error, hidden_error, &error);
				training_times++;
                error_acc += error;

			}while(training_times < TRAIN_SET_SIZE);

			rewind(fp);  // back to the start of the input file, start epoch
			if(i == 0){
				printf("In the first training, Accumulative Error = %f\n", error_acc);
				puts("========Start Epoch=========");
			}else{
				printf("Epoch %d, Accumulative Error = %f\n", i, error_acc);
			}
			error_acc = 0.0;
            training_times = 0;

        }

        gsl_matrix_free(better_input), gsl_matrix_free(target), gsl_matrix_free(output_error), gsl_matrix_free(hidden_error), gsl_matrix_free(output_x), gsl_matrix_free(output_o);
        gsl_matrix_free(hidden_x), gsl_matrix_free(hidden_o);

        store_weight(input_weight, hidden_weight);
        puts("Training weight has been stored!");
        gsl_matrix_free(input_weight), gsl_matrix_free(hidden_weight);

	}

	if(!strcmp(mode, "-r")){
		input_weight = gsl_matrix_alloc(HIDDEN_NODE, DATA_NUM_WITHOUT_LABEL);
		hidden_weight = gsl_matrix_alloc(OUTPUT_NODE, HIDDEN_NODE);

		load_weight(input_weight, hidden_weight);

		int test_times = 0;
		int wrong_times = 0;
		int anwser;

		do{
            read_parse(fp, input_data, &label, mode);
			better_input = scale_input(input_data);
			anwser = recognize(better_input, input_weight, hidden_weight, hidden_x, hidden_o, output_x, output_o);

			printf("Label: %d, Guess: %d\n", label, anwser);

			if(anwser != label) wrong_times++;

			test_times++;

		}while(test_times < TEST_SET_SIZE);
		printf("Performance: %2f\n", (1 - ((float)wrong_times/(float)test_times)) * 100);
		gsl_matrix_free(better_input), gsl_matrix_free(target), gsl_matrix_free(output_error), gsl_matrix_free(hidden_error), gsl_matrix_free(output_x), gsl_matrix_free(output_o);
		gsl_matrix_free(hidden_x), gsl_matrix_free(hidden_o);
		gsl_matrix_free(input_weight), gsl_matrix_free(hidden_weight);

	}
        fclose(fp);
		return 0;
}

FILE* openfile(char* file_name)
{
	FILE* fp;
	if(!(fp = fopen(file_name, "r")))
		fprintf(stderr, "Can't open the file\n");
	return fp;
}

void read_parse(FILE* fp, int input_data[], int* label, char* mode)
{
	char line[LINE_LEN];
	char temp;
	int line_c = 0;
	int input_c = 0;
	char* token;

	while((temp = fgetc(fp)) != '\n'){
        if(temp == EOF) break;
        line[line_c++] = temp;
	}
    line[line_c] = '\0';
    line_c = 0;
	if(!strcmp(mode, "-t")){
		token = strtok(line, ",");
		*label = atoi(token);
		while(token != NULL){
			token = strtok(NULL, ",");
			if(token == NULL) break;
			input_data[input_c++] = atoi(token);
		}

	}
	if(!strcmp(mode, "-r")){
		token = strtok(line, ",");
		*label = atoi(token);
		while(token != NULL){
			token = strtok(NULL, ",");
			if(token == NULL) break;
			input_data[input_c++] = atoi(token);
		}
	}
}

gsl_matrix* scale_input(int input_data[])
{
    gsl_matrix* better_input = gsl_matrix_alloc(DATA_NUM_WITHOUT_LABEL, 1);
	for(int i = 0; i < DATA_NUM_WITHOUT_LABEL; i++){
		gsl_matrix_set(better_input, i, 0, ((double)input_data[i]/255.0*0.99)+0.01);
	}
    return better_input;
}

gsl_matrix* vectorize_target(int label)
{
	gsl_matrix* target = gsl_matrix_alloc(10, 1);

	for(int i = 0; i < 10; i++){
		if(i == label){
			gsl_matrix_set(target, i, 0, 0.99);
		}else{
			gsl_matrix_set(target, i, 0, 0.01);
		}
	}

	return target;
}

gsl_matrix* init_weight(int row, int col)
{
    gsl_matrix* input_weight = gsl_matrix_alloc(row, col);

    const gsl_rng_type * T;
    gsl_rng * r;
    double sigma = NOR_DIS_DEV(col);
    r = gsl_rng_alloc (gsl_rng_mt19937);
    for(int i = 0; i < row; i++){

        for(int j = 0; j < col; j++){
            gsl_matrix_set(input_weight, i, j, gsl_ran_gaussian(r, sigma));
        }
    }
    gsl_rng_free(r);
    return input_weight;
}

double get_output_error(gsl_matrix* target, gsl_matrix* output, gsl_matrix* output_error)
{
	double error = 0.0;
	for (int i = 0; i < OUTPUT_NODE; i++){
		gsl_matrix_set(output_error, i, 0, (gsl_matrix_get(target, i, 0) - gsl_matrix_get(output, i, 0)));
	}

	for (int i = 0; i < OUTPUT_NODE; i++){
		error += pow(gsl_matrix_get(output_error, i, 0), 2);
	}

	return error;
}

void training(gsl_matrix* better_input, gsl_matrix* target, gsl_matrix* input_weight, gsl_matrix* hidden_weight, gsl_matrix* hidden_x, gsl_matrix* hidden_o, gsl_matrix* output_o, gsl_matrix* output_x, gsl_matrix* output_error, gsl_matrix* hidden_error, double* error)
{
	gsl_matrix* delta_hidden_weight = gsl_matrix_alloc(OUTPUT_NODE, HIDDEN_NODE);
	gsl_matrix* delta_input_weight = gsl_matrix_alloc(HIDDEN_NODE, DATA_NUM_WITHOUT_LABEL);
	gsl_matrix* temp_hidden = gsl_matrix_alloc(OUTPUT_NODE, 1);
	gsl_matrix* temp_input = gsl_matrix_alloc(HIDDEN_NODE, 1);

	 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, input_weight, better_input, 0.0, hidden_x);  // X = W*I
	 for(int i = 0; i < HIDDEN_NODE; i++){
	 	gsl_matrix_set(hidden_o, i, 0, SIGMOID(gsl_matrix_get(hidden_x, i, 0)));
	 }

	 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, hidden_weight, hidden_o, 0.0, output_x);
	 for(int i = 0; i < OUTPUT_NODE; i++){
	 	gsl_matrix_set(output_o, i, 0, SIGMOID(gsl_matrix_get(output_x, i, 0)));
	 }

	 *error = get_output_error(target, output_o, output_error);
	 gsl_blas_dgemm(CblasTrans, CblasNoTrans, 1.0, hidden_weight, output_error, 0.0, hidden_error); // Error hidden

	 for(int i = 0; i < OUTPUT_NODE; i++){
	 	gsl_matrix_set(temp_hidden, i, 0, (LEARN_RATE * gsl_matrix_get(output_error, i, 0) * gsl_matrix_get(output_o, i, 0) * (1.0 - gsl_matrix_get(output_o, i, 0))));
	 }
	 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, temp_hidden, hidden_o, 0.0, delta_hidden_weight);

	 for(int i = 0; i < HIDDEN_NODE; i++){
	 	gsl_matrix_set(temp_input, i, 0, (LEARN_RATE * gsl_matrix_get(hidden_error, i, 0) * gsl_matrix_get(hidden_o, i, 0) * (1.0 - gsl_matrix_get(hidden_error, i, 0))));
	 }
	 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, temp_input, better_input, 0.0, delta_input_weight);

	 gsl_matrix_add(hidden_weight, delta_hidden_weight);
	 gsl_matrix_add(input_weight, delta_input_weight);

	 gsl_matrix_free(delta_input_weight);
	 gsl_matrix_free(delta_hidden_weight);
	 gsl_matrix_free(temp_input);
	 gsl_matrix_free(temp_hidden);

}

void store_weight(gsl_matrix* input_weight, gsl_matrix* hidden_weight)
{
	FILE* fp = fopen("Good_Weight", "w+");

	for(int i = 0; i < HIDDEN_NODE; i++){
		for(int j = 0; j < DATA_NUM_WITHOUT_LABEL; j++){

			fprintf(fp, "%f,", gsl_matrix_get(input_weight, i, j));
		}
	}

	fseek(fp, -1, SEEK_CUR);
	fputc('\n', fp);

	for(int i = 0; i < OUTPUT_NODE; i++){
		for(int j = 0; j < HIDDEN_NODE; j++){

			fprintf(fp, "%f,", gsl_matrix_get(hidden_weight, i, j));
		}
	}
	fseek(fp, -1, SEEK_CUR);
	fputc('\n', fp);

	fclose(fp);
}

void load_weight(gsl_matrix* input_weight, gsl_matrix* hidden_weight)
{
	FILE* fp = fopen("Good_Weight", "r");
	char num_temp[20];
	double array[HIDDEN_NODE*DATA_NUM_WITHOUT_LABEL];
	char c;
	int j = 0, k;

	if(fp == NULL){
		fprintf(stderr, "No Good_Weight here, Please training first\n");
		exit(EXIT_FAILURE);
	}

	while((c = fgetc(fp)) != '\n'){

		if(c == ','){
			k = 0;
			array[j++] = atof(num_temp);
			bzero(num_temp, sizeof(num_temp));
			continue;
		}
		num_temp[k++] = c;
	}

	array[j] = atof(num_temp);

	j = 0;

	for(int n = 0; n < HIDDEN_NODE; n++){
		for(int m = 0; m < DATA_NUM_WITHOUT_LABEL; m++){
			gsl_matrix_set(input_weight, n, m, array[j++]);
		}
	}

	j = 0;
	k = 0;
	while((c = fgetc(fp)) != '\n'){

		if(c == ','){
			k = 0;
			array[j++] = atof(num_temp);
			bzero(num_temp, sizeof(num_temp));
			continue;
		}
		num_temp[k++] = c;
	}

	array[j] = atof(num_temp);

	j = 0;

	for(int n = 0; n < OUTPUT_NODE; n++){
		for(int m = 0; m < HIDDEN_NODE; m++){
			gsl_matrix_set(hidden_weight, n, m, array[j++]);
		}
	}

	fclose(fp);
}

int recognize(gsl_matrix* better_input, gsl_matrix* input_weight, gsl_matrix* hidden_weight, gsl_matrix* hidden_x, gsl_matrix* hidden_o, gsl_matrix* output_x, gsl_matrix* output_o)
{
	size_t row, col;
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, input_weight, better_input, 0.0, hidden_x);  // X = W*I
	 for(int i = 0; i < HIDDEN_NODE; i++){
	 	gsl_matrix_set(hidden_o, i, 0, SIGMOID(gsl_matrix_get(hidden_x, i, 0)));
	 }

	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, hidden_weight, hidden_o, 0.0, output_x);
	 for(int i = 0; i < OUTPUT_NODE; i++){
	 	gsl_matrix_set(output_o, i, 0, SIGMOID(gsl_matrix_get(output_x, i, 0)));
	 }

	 gsl_matrix_max_index(output_o, &row, &col);

	 return (int)row;
}


text:
/* matrix/gsl_matrix_uchar.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 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 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_MATRIX_UCHAR_H__
#define __GSL_MATRIX_UCHAR_H__

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_uchar.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 size1;
  size_t size2;
  size_t tda;
  unsigned char * data;
  gsl_block_uchar * block;
  int owner;
} gsl_matrix_uchar;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_view;

typedef _gsl_matrix_uchar_view gsl_matrix_uchar_view;

typedef struct
{
  gsl_matrix_uchar matrix;
} _gsl_matrix_uchar_const_view;

typedef const _gsl_matrix_uchar_const_view gsl_matrix_uchar_const_view;

/* Allocation */

GSL_EXPORT
gsl_matrix_uchar *
gsl_matrix_uchar_alloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_uchar *
gsl_matrix_uchar_calloc (const size_t n1, const size_t n2);

GSL_EXPORT
gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_block (gsl_block_uchar * b,
                                   const size_t offset,
                                   const size_t n1,
                                   const size_t n2,
                                   const size_t d2);

GSL_EXPORT
gsl_matrix_uchar *
gsl_matrix_uchar_alloc_from_matrix (gsl_matrix_uchar * m,
                                    const size_t k1,
                                    const size_t k2,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
gsl_vector_uchar *
gsl_vector_uchar_alloc_row_from_matrix (gsl_matrix_uchar * m,
                                        const size_t i);

GSL_EXPORT
gsl_vector_uchar *
gsl_vector_uchar_alloc_col_from_matrix (gsl_matrix_uchar * m,
                                        const size_t j);

GSL_EXPORT void gsl_matrix_uchar_free (gsl_matrix_uchar * m);

/* Views */

GSL_EXPORT
_gsl_matrix_uchar_view
gsl_matrix_uchar_submatrix (gsl_matrix_uchar * m,
                            const size_t i, const size_t j,
                            const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_matrix_uchar_row (gsl_matrix_uchar * m, const size_t i);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_matrix_uchar_column (gsl_matrix_uchar * m, const size_t j);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_matrix_uchar_diagonal (gsl_matrix_uchar * m);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_matrix_uchar_subdiagonal (gsl_matrix_uchar * m, const size_t k);

GSL_EXPORT
_gsl_vector_uchar_view
gsl_matrix_uchar_superdiagonal (gsl_matrix_uchar * m, const size_t k);

GSL_EXPORT
_gsl_matrix_uchar_view
gsl_matrix_uchar_view_array (unsigned char * base,
                             const size_t n1,
                             const size_t n2);

GSL_EXPORT
_gsl_matrix_uchar_view
gsl_matrix_uchar_view_array_with_tda (unsigned char * base,
                                      const size_t n1,
                                      const size_t n2,
                                      const size_t tda);

GSL_EXPORT
_gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector (gsl_vector_uchar * v,
                              const size_t n1,
                              const size_t n2);

GSL_EXPORT
_gsl_matrix_uchar_view
gsl_matrix_uchar_view_vector_with_tda (gsl_vector_uchar * v,
                                       const size_t n1,
                                       const size_t n2,
                                       const size_t tda);

GSL_EXPORT
_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_submatrix (const gsl_matrix_uchar * m,
                                  const size_t i, const size_t j,
                                  const size_t n1, const size_t n2);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_row (const gsl_matrix_uchar * m,
                            const size_t i);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_column (const gsl_matrix_uchar * m,
                               const size_t j);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_diagonal (const gsl_matrix_uchar * m);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_subdiagonal (const gsl_matrix_uchar * m,
                                    const size_t k);

GSL_EXPORT
_gsl_vector_uchar_const_view
gsl_matrix_uchar_const_superdiagonal (const gsl_matrix_uchar * m,
                                      const size_t k);

GSL_EXPORT
_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array (const unsigned char * base,
                                   const size_t n1,
                                   const size_t n2);

GSL_EXPORT
_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_array_with_tda (const unsigned char * base,
                                            const size_t n1,
                                            const size_t n2,
                                            const size_t tda);

GSL_EXPORT
_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector (const gsl_vector_uchar * v,
                                    const size_t n1,
                                    const size_t n2);

GSL_EXPORT
_gsl_matrix_uchar_const_view
gsl_matrix_uchar_const_view_vector_with_tda (const gsl_vector_uchar * v,
                                             const size_t n1,
                                             const size_t n2,
                                             const size_t tda);

/* Operations */

GSL_EXPORT unsigned char   gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_EXPORT void    gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x);

GSL_EXPORT unsigned char * gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_EXPORT const unsigned char * gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j);

GSL_EXPORT void gsl_matrix_uchar_set_zero (gsl_matrix_uchar * m);
GSL_EXPORT void gsl_matrix_uchar_set_identity (gsl_matrix_uchar * m);
GSL_EXPORT void gsl_matrix_uchar_set_all (gsl_matrix_uchar * m, unsigned char x);

GSL_EXPORT int gsl_matrix_uchar_fread (FILE * stream, gsl_matrix_uchar * m) ;
GSL_EXPORT int gsl_matrix_uchar_fwrite (FILE * stream, const gsl_matrix_uchar * m) ;
GSL_EXPORT int gsl_matrix_uchar_fscanf (FILE * stream, gsl_matrix_uchar * m);
GSL_EXPORT int gsl_matrix_uchar_fprintf (FILE * stream, const gsl_matrix_uchar * m, const char * format);

GSL_EXPORT int gsl_matrix_uchar_memcpy(gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);
GSL_EXPORT int gsl_matrix_uchar_swap(gsl_matrix_uchar * m1, gsl_matrix_uchar * m2);

GSL_EXPORT int gsl_matrix_uchar_swap_rows(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_uchar_swap_columns(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_uchar_swap_rowcol(gsl_matrix_uchar * m, const size_t i, const size_t j);
GSL_EXPORT int gsl_matrix_uchar_transpose (gsl_matrix_uchar * m);
GSL_EXPORT int gsl_matrix_uchar_transpose_memcpy (gsl_matrix_uchar * dest, const gsl_matrix_uchar * src);

GSL_EXPORT unsigned char gsl_matrix_uchar_max (const gsl_matrix_uchar * m);
GSL_EXPORT unsigned char gsl_matrix_uchar_min (const gsl_matrix_uchar * m);
GSL_EXPORT void gsl_matrix_uchar_minmax (const gsl_matrix_uchar * m, unsigned char * min_out, unsigned char * max_out);

GSL_EXPORT void gsl_matrix_uchar_max_index (const gsl_matrix_uchar * m, size_t * imax, size_t *jmax);
GSL_EXPORT void gsl_matrix_uchar_min_index (const gsl_matrix_uchar * m, size_t * imin, size_t *jmin);
GSL_EXPORT void gsl_matrix_uchar_minmax_index (const gsl_matrix_uchar * m, size_t * imin, size_t * jmin, size_t * imax, size_t * jmax);

GSL_EXPORT int gsl_matrix_uchar_isnull (const gsl_matrix_uchar * m);

GSL_EXPORT int gsl_matrix_uchar_add (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_EXPORT int gsl_matrix_uchar_sub (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_EXPORT int gsl_matrix_uchar_mul_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_EXPORT int gsl_matrix_uchar_div_elements (gsl_matrix_uchar * a, const gsl_matrix_uchar * b);
GSL_EXPORT int gsl_matrix_uchar_scale (gsl_matrix_uchar * a, const double x);
GSL_EXPORT int gsl_matrix_uchar_add_constant (gsl_matrix_uchar * a, const double x);
GSL_EXPORT int gsl_matrix_uchar_add_diagonal (gsl_matrix_uchar * a, const double x);

/***********************************************************************/
/* The functions below are obsolete                                    */
/***********************************************************************/
GSL_EXPORT int gsl_matrix_uchar_get_row(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t i);
GSL_EXPORT int gsl_matrix_uchar_get_col(gsl_vector_uchar * v, const gsl_matrix_uchar * m, const size_t j);
GSL_EXPORT int gsl_matrix_uchar_set_row(gsl_matrix_uchar * m, const size_t i, const gsl_vector_uchar * v);
GSL_EXPORT int gsl_matrix_uchar_set_col(gsl_matrix_uchar * m, const size_t j, const gsl_vector_uchar * v);

/* inline functions if you are using GCC */

#ifdef HAVE_INLINE
extern inline
unsigned char
gsl_matrix_uchar_get(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VAL("first index out of range", GSL_EINVAL, 0) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VAL("second index out of range", GSL_EINVAL, 0) ;
    }
#endif
  return m->data[i * m->tda + j] ;
}

extern inline
void
gsl_matrix_uchar_set(gsl_matrix_uchar * m, const size_t i, const size_t j, const unsigned char x)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_VOID("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_VOID("second index out of range", GSL_EINVAL) ;
    }
#endif
  m->data[i * m->tda + j] = x ;
}

extern inline
unsigned char *
gsl_matrix_uchar_ptr(gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (unsigned char *) (m->data + (i * m->tda + j)) ;
}

extern inline
const unsigned char *
gsl_matrix_uchar_const_ptr(const gsl_matrix_uchar * m, const size_t i, const size_t j)
{
#if GSL_RANGE_CHECK
  if (i >= m->size1)
    {
      GSL_ERROR_NULL("first index out of range", GSL_EINVAL) ;
    }
  else if (j >= m->size2)
    {
      GSL_ERROR_NULL("second index out of range", GSL_EINVAL) ;
    }
#endif
  return (const unsigned char *) (m->data + (i * m->tda + j)) ;
}

#endif

__END_DECLS

#endif /* __GSL_MATRIX_UCHAR_H__ */

ext: h

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws