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

#include <gsl/gsl_integration.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>

#include <class.h> /* from extern/ */

#include "ccl.h"
#include "ccl_params.h"
#include "ccl_emu17.h"
#include "ccl_emu17_params.h"

/*------ ROUTINE: ccl_cosmology_compute_power_class -----
INPUT: ccl_cosmology * cosmo
*/
static void ccl_free_class_structs(ccl_cosmology *cosmo,
				   struct background *ba,
				   struct thermo *th,
				   struct perturbs *pt,
				   struct transfers *tr,
				   struct primordial *pm,
				   struct spectra *sp,
				   struct nonlinear *nl,
				   struct lensing *le,
				   int *init_arr,
				   int * status)
{
  int i_init=6;
  if(init_arr[i_init--]) {
    if (spectra_free(sp) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS spectra:%s\n", sp->error_message);
      return;
    }
  }

  if(init_arr[i_init--]) {
    if (transfer_free(tr) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS transfer:%s\n", tr->error_message);
      return;
    }
  }

  if(init_arr[i_init--]) {
    if (nonlinear_free(nl) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS nonlinear:%s\n", nl->error_message);
      return;
    }
  }

  if(init_arr[i_init--]) {
    if (primordial_free(pm) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS pm:%s\n", pm->error_message);
      return;
    }
  }

  if(init_arr[i_init--]) {
    if (perturb_free(pt) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS pt:%s\n", pt->error_message);
      return;
    }
  }

  if(init_arr[i_init--]) {
    if (thermodynamics_free(th) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS thermo:%s\n", th->error_message);
      return;
    }
  }

  if(init_arr[i_init--]) {
    if (background_free(ba) == _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_free_class_structs(): Error freeing CLASS bg:%s\n", ba->error_message);
      return;
    }
  }

  return;
}

static void ccl_class_preinit(struct background *ba,
			      struct thermo *th,
			      struct perturbs *pt,
			      struct transfers *tr,
			      struct primordial *pm,
			      struct spectra *sp,
			      struct nonlinear *nl,
			      struct lensing *le)
{
  //pre-initialize all fields that are freed by *_free() routine
  //prevents crashes if *_init()failed and did not initialize all tables freed by *_free()

  //init for background_free
  ba->tau_table = NULL;
  ba->z_table = NULL;
  ba->d2tau_dz2_table = NULL;
  ba->background_table = NULL;
  ba->d2background_dtau2_table = NULL;

  //init for thermodynamics_free
  th->z_table = NULL;
  th->thermodynamics_table = NULL;
  th->d2thermodynamics_dz2_table = NULL;

  //init for perturb_free
  pt->tau_sampling = NULL;
  pt->tp_size = NULL;
  pt->ic_size = NULL;
  pt->k = NULL;
  pt->k_size_cmb = NULL;
  pt->k_size_cl = NULL;
  pt->k_size = NULL;
  pt->sources = NULL;

  //init for primordial_free
  pm->amplitude = NULL;
  pm->tilt = NULL;
  pm->running = NULL;
  pm->lnpk = NULL;
  pm->ddlnpk = NULL;
  pm->is_non_zero = NULL;
  pm->ic_size = NULL;
  pm->ic_ic_size = NULL;
  pm->lnk = NULL;

  //init for nonlinear_free
  nl->k = NULL;
  nl->tau = NULL;
  nl->nl_corr_density = NULL;
  nl->k_nl = NULL;

  //init for transfer_free
  tr->tt_size = NULL;
  tr->l_size_tt = NULL;
  tr->l_size = NULL;
  tr->l = NULL;
  tr->q = NULL;
  tr->k = NULL;
  tr->transfer = NULL;

  //init for spectra_free
  //spectra_free checks all other data fields before freeing
  sp->is_non_zero = NULL;
  sp->ic_size = NULL;
  sp->ic_ic_size = NULL;
}

static void ccl_run_class(ccl_cosmology *cosmo,
			  struct file_content *fc,
			  struct precision* pr,
			  struct background* ba,
			  struct thermo* th,
			  struct perturbs* pt,
			  struct transfers* tr,
			  struct primordial* pm,
			  struct spectra* sp,
			  struct nonlinear* nl,
			  struct lensing* le,
			  struct output* op,
			  int *init_arr,
			  int * status)
{
  ErrorMsg errmsg;            // for error messages
  int i_init=0;
  ccl_class_preinit(ba,th,pt,tr,pm,sp,nl,le);

  if(input_init(fc,pr,ba,th,pt,tr,pm,sp,nl,le,op,errmsg) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS input:%s\n", errmsg);
    return;
  }
  if (background_init(pr,ba) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS background:%s\n", ba->error_message);
    return;
  }
  init_arr[i_init++]=1;
  if (thermodynamics_init(pr,ba,th) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS thermodynamics:%s\n", th->error_message);
    return;
  }
  init_arr[i_init++]=1;
  if (perturb_init(pr,ba,th,pt) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS pertubations:%s\n", pt->error_message);
    return;
  }
  init_arr[i_init++]=1;
  if (primordial_init(pr,pt,pm) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS primordial:%s\n", pm->error_message);
    return;
  }
  init_arr[i_init++]=1;
  if (nonlinear_init(pr,ba,th,pt,pm,nl) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS nonlinear:%s\n", nl->error_message);
    return;
  }
  init_arr[i_init++]=1;
  if (transfer_init(pr,ba,th,pt,nl,tr) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS transfer:%s\n", tr->error_message);
    return;
  }
  init_arr[i_init++]=1;
  if (spectra_init(pr,ba,pt,pm,nl,tr,sp) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS spectra:%s\n", sp->error_message);
    return;
  }
  init_arr[i_init++]=1;
}

static double ccl_get_class_As(ccl_cosmology *cosmo, struct file_content *fc, int position_As,
			       double sigma8, int * status)
{
  //structures for class test run
  struct precision pr;        // for precision parameters
  struct background ba;       // for cosmological background
  struct thermo th;           // for thermodynamics
  struct perturbs pt;         // for source functions
  struct transfers tr;        // for transfer functions
  struct primordial pm;       // for primordial spectra
  struct spectra sp;          // for output spectra
  struct nonlinear nl;        // for non-linear spectra
  struct lensing le;
  struct output op;

  //temporarily overwrite P_k_max_1/Mpc to speed up sigma8 calculation
  double k_max_old = 0.;
  int position_kmax =2;
  double A_s_guess;
  int init_arr[7]={0,0,0,0,0,0,0};

  if (strcmp(fc->name[position_kmax],"P_k_max_1/Mpc")) {
    k_max_old = strtof(fc->value[position_kmax],NULL);
    sprintf(fc->value[position_kmax],"%.15e",10.);
  }
  A_s_guess = 2.43e-9/0.87659*sigma8;
  sprintf(fc->value[position_As],"%.15e",A_s_guess);

  ccl_run_class(cosmo, fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,init_arr,status);
  if (cosmo->status != CCL_ERROR_CLASS) A_s_guess*=pow(sigma8/sp.sigma8,2.);
  ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);

  if (k_max_old >0) {
    sprintf(fc->value[position_kmax],"%.15e",k_max_old);
  }
  return A_s_guess;
}

static void ccl_fill_class_parameters(ccl_cosmology * cosmo, struct file_content * fc,
				      int parser_length, int * status)

{

  // initialize fc fields
  //silences Valgrind's "Conditional jump or move depends on uninitialised value" warning
  for (int i = 0; i< parser_length; i++){
    strcpy(fc->name[i]," ");
    strcpy(fc->value[i]," ");
  }

  strcpy(fc->name[0],"output");
  strcpy(fc->value[0],"mPk");

  strcpy(fc->name[1],"non linear");
  if (cosmo->config.matter_power_spectrum_method == ccl_halofit)
    strcpy(fc->value[1],"Halofit");
  else
    strcpy(fc->value[1],"none");

  strcpy(fc->name[2],"P_k_max_1/Mpc");
  sprintf(fc->value[2],"%.15e",ccl_splines->K_MAX_SPLINE); //in units of 1/Mpc, corroborated with ccl_constants.h

  strcpy(fc->name[3],"z_max_pk");
  sprintf(fc->value[3],"%.15e",1./ccl_splines->A_SPLINE_MINLOG_PK-1.);

  strcpy(fc->name[4],"modes");
  strcpy(fc->value[4],"s");

  strcpy(fc->name[5],"lensing");
  strcpy(fc->value[5],"no");

  // now, copy over cosmology parameters
  strcpy(fc->name[6],"h");
  sprintf(fc->value[6],"%.15e",cosmo->params.h);

  strcpy(fc->name[7],"Omega_cdm");
  sprintf(fc->value[7],"%.15e",cosmo->params.Omega_c);

  strcpy(fc->name[8],"Omega_b");
  sprintf(fc->value[8],"%.15e",cosmo->params.Omega_b);

  strcpy(fc->name[9],"Omega_k");
  sprintf(fc->value[9],"%.15e",cosmo->params.Omega_k);

  strcpy(fc->name[10],"n_s");
  sprintf(fc->value[10],"%.15e",cosmo->params.n_s);

  //cosmological constant?
  // set Omega_Lambda = 0.0 if w !=-1
  if ((cosmo->params.w0 !=-1.0) || (cosmo->params.wa !=0)) {
    strcpy(fc->name[11],"Omega_Lambda");
    sprintf(fc->value[11],"%.15e",0.0);

    strcpy(fc->name[12],"w0_fld");
    sprintf(fc->value[12],"%.15e",cosmo->params.w0);

    strcpy(fc->name[13],"wa_fld");
    sprintf(fc->value[13],"%.15e",cosmo->params.wa);
  }
  //neutrino parameters
  //massless neutrinos
  if (cosmo->params.N_nu_rel > 1.e-4) {
    strcpy(fc->name[14],"N_ur");
    sprintf(fc->value[14],"%.15e",cosmo->params.N_nu_rel);
  }
  else {
    strcpy(fc->name[14],"N_ur");
    sprintf(fc->value[14],"%.15e", 0.);
  }
  if (cosmo->params.N_nu_mass > 0) {
    strcpy(fc->name[15],"N_ncdm");
    sprintf(fc->value[15],"%d",cosmo->params.N_nu_mass);
    strcpy(fc->name[16],"m_ncdm");
    sprintf(fc->value[16],"%f", (cosmo->params.mnu)[0]);
    if (cosmo->params.N_nu_mass >=1){
		for (int i = 1; i < cosmo->params.N_nu_mass; i++) {
			char tmp[20];
			sprintf(tmp,", %f",(cosmo->params.mnu)[i]);
			strcat(fc->value[16],tmp);
		}
	}

  }

  strcpy(fc->name[17],"T_cmb");
  sprintf(fc->value[17],"%.15e",cosmo->params.T_CMB);

  //normalization comes last, so that all other parameters are filled in for determining A_s if sigma8 is specified
  if (isfinite(cosmo->params.sigma8) && isfinite(cosmo->params.A_s)){
      *status = CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: class_parameters(): Error initializing CLASS parameters: both sigma8 and A_s defined\n");
    return;
  }
  if (isfinite(cosmo->params.sigma8)) {
    strcpy(fc->name[parser_length-1],"A_s");
    sprintf(fc->value[parser_length-1],"%.15e",ccl_get_class_As(cosmo,fc,parser_length-1,cosmo->params.sigma8, status));
  }
  else if (isfinite(cosmo->params.A_s)) {
    strcpy(fc->name[parser_length-1],"A_s");
    sprintf(fc->value[parser_length-1],"%.15e",cosmo->params.A_s);
  }
  else {
    *status = CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: class_parameters(): Error initializing CLASS pararmeters: neither sigma8 nor A_s defined\n");
    return;
  }

}

static void ccl_cosmology_compute_power_class(ccl_cosmology * cosmo, int * status)
{
  struct precision pr;        // for precision parameters
  struct background ba;       // for cosmological background
  struct thermo th;           // for thermodynamics
  struct perturbs pt;         // for source functions
  struct transfers tr;        // for transfer functions
  struct primordial pm;       // for primordial spectra
  struct spectra sp;          // for output spectra
  struct nonlinear nl;        // for non-linear spectra
  struct lensing le;
  struct output op;
  struct file_content fc;

  ErrorMsg errmsg; // for error messages
  // generate file_content structure
  // CLASS configuration parameters will be passed through this structure,
  // to avoid writing and reading .ini files for every call
  int parser_length = 20;
  int init_arr[7]={0,0,0,0,0,0,0};
  if (parser_init(&fc,parser_length,"none",errmsg) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): parser init error:%s\n", errmsg);
    return;
  }

  ccl_fill_class_parameters(cosmo,&fc,parser_length, status);

  if (*status != CCL_ERROR_CLASS)
    ccl_run_class(cosmo, &fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,init_arr,status);

  if (*status == CCL_ERROR_CLASS) {
    //printed error message while running CLASS
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }

  if (parser_free(&fc)== _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error freeing CLASS parser\n");
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }

  //These are the limits of the splining range
  cosmo->data.k_min_lin=2*exp(sp.ln_k[0]);
  cosmo->data.k_max_lin=ccl_splines->K_MAX_SPLINE;

  //CLASS calculations done - now allocate CCL splines
  double kmin = cosmo->data.k_min_lin;
  double kmax = ccl_splines->K_MAX_SPLINE;
  //Compute nk from number of decades and N_K = # k per decade
  double ndecades = log10(kmax) - log10(kmin);
  int nk = (int)ceil(ndecades*ccl_splines->N_K);
  double amin = ccl_splines->A_SPLINE_MINLOG_PK;
  double amax = ccl_splines->A_SPLINE_MAX;
  int na = ccl_splines->A_SPLINE_NA_PK+ccl_splines->A_SPLINE_NLOG_PK-1;

  // The x array is initially k, but will later
  // be overwritten with log(k)
  double * x = ccl_log_spacing(kmin, kmax, nk);
  double * a = ccl_linlog_spacing(amin, ccl_splines->A_SPLINE_MIN_PK, amax, ccl_splines->A_SPLINE_NLOG_PK, ccl_splines->A_SPLINE_NA_PK);
  double * y2d_lin = malloc(nk * na * sizeof(double));
  double * y2d_nl = malloc(nk * na * sizeof(double));

  //If error, store status, we will free later
  if (a==NULL|| x==NULL || y2d_lin==NULL || y2d_nl==NULL) {
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): memory allocation error\n");
  }

  //Status flags
  int newstatus=0;
  int pwstatus=0;

  //If not, proceed
  if(!*status){

    // After this loop x will contain log(k)
    // all in Mpc, not Mpc/h units!
    double psout_l,ic;
    for (int i=0; i<nk; i++) {
      for (int j = 0; j < na; j++) {
	//The 2D interpolation routines access the function values y_{k_ia_j} with the following ordering:
	//y_ij = y2d[j*N_k + i]
	//with i = 0,...,N_k-1 and j = 0,...,N_a-1.
	newstatus |= spectra_pk_at_k_and_z(&ba, &pm, &sp,x[i],1./a[j]-1., &psout_l,&ic);
	y2d_lin[j*nk+i] = log(psout_l);
      }
      x[i] = log(x[i]);
    }

    //If error, store status, we will free later
    if(newstatus) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error computing CLASS power spectrum\n");
    }

  }

  //If no error, proceed
  if(!*status) {

    gsl_spline2d * log_power = gsl_spline2d_alloc(PLIN_SPLINE_TYPE, nk,na);
    pwstatus = gsl_spline2d_init(log_power, x, a, y2d_lin,nk,na);

    //If not, proceed
    if(!pwstatus){
      cosmo->data.p_lin = log_power;
    } else {
      gsl_spline2d_free(log_power);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error creating log_power spline\n");
    }

  }

  if(*status){ //Linear power spec failed, so we return without proceeding to nonlinear.
    free(x);
    free(a);
    free(y2d_nl);
    free(y2d_lin);
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }

  //Non-linear power
  //At the moment KMIN can't be less than CLASS's kmin in the nonlinear case.

  //If error, store status, we will free later
  if (kmin<(exp(sp.ln_k[0]))) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): K_MIN is less than CLASS's kmin. Not yet supported for nonlinear P(k).\n");
  }

  //If not, proceed
  if(!*status){

    //These are the limits of the splining range
    cosmo->data.k_min_nl=2*exp(sp.ln_k[0]);
    cosmo->data.k_max_nl=ccl_splines->K_MAX_SPLINE;

    if(cosmo->config.matter_power_spectrum_method==ccl_halofit) {

      double psout_nl;
      for (int i=0; i<nk; i++) {
	for (int j = 0; j < na; j++) {
	  newstatus |= spectra_pk_nl_at_k_and_z(&ba, &pm, &sp,exp(x[i]),1./a[j]-1.,&psout_nl);
	  y2d_nl[j*nk+i] = log(psout_nl);
	}
      }
    }

    if(newstatus){
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error computing CLASS power spectrum\n");
    }

  }

  if(!*status){

    gsl_spline2d * log_power_nl = gsl_spline2d_alloc(PNL_SPLINE_TYPE, nk,na);
    pwstatus = gsl_spline2d_init(log_power_nl, x, a, y2d_nl,nk,na);

    if(!pwstatus){
      cosmo->data.p_nl = log_power_nl;
    } else {
      gsl_spline2d_free(log_power_nl);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error creating log_power_nl spline\n");
    }

  }

  ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
  free(x);
  free(a);
  free(y2d_nl);
  free(y2d_lin);

  return;

}

/* BCM correction
   See Schneider & Teyssier (2015) for details of the model.
   The equations inside this function correspond to those in that
   work for:
   gf -> G(k)
   scomp -> S(k)
*/
double ccl_bcm_model_fka(ccl_cosmology * cosmo, double k, double a, int *status){

  double fka;
  double b0;
  double bfunc, bfunc4;
  double kg;
  double gf,scomp;
  double kh;
  double z;

  z=1./a-1.;
  kh = k/cosmo->params.h; //convert to h/Mpc
  b0 = 0.105*cosmo->params.bcm_log10Mc-1.27; //Eq. 4.4
  bfunc = b0/(1.+pow(z/2.3,2.5)); //Eq. 4.3
  bfunc4 = (1-bfunc)*(1-bfunc)*(1-bfunc)*(1-bfunc); //B^4(z)
  kg = 0.7*bfunc4*pow(cosmo->params.bcm_etab,-1.6); //Eq. 4.3
  gf = bfunc/(1+pow(kh/kg,3.))+1.-bfunc; //Eq. 4.2, k in h/Mpc
  scomp = 1+(kh/cosmo->params.bcm_ks)*(kh/cosmo->params.bcm_ks); //Eq 4.5, k in h/Mpc
  fka = gf*scomp;
  return fka;
}

void ccl_cosmology_write_power_class_z(char *filename, ccl_cosmology * cosmo, double z, int * status)
{
  struct precision pr;        // for precision parameters
  struct background ba;       // for cosmological background
  struct thermo th;           // for thermodynamics
  struct perturbs pt;         // for source functions
  struct transfers tr;        // for transfer functions
  struct primordial pm;       // for primordial spectra
  struct spectra sp;          // for output spectra
  struct nonlinear nl;        // for non-linear spectra
  struct lensing le;
  struct output op;
  struct file_content fc;

  ErrorMsg errmsg; // for error messages
  // generate file_content structure
  // CLASS configuration parameters will be passed through this structure,
  // to avoid writing and reading .ini files for every call
  int parser_length = 20;
  int init_arr[7]={0,0,0,0,0,0,0};
  if (parser_init(&fc,parser_length,"none",errmsg) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): parser init error:%s\n", errmsg);
    return;
  }

  ccl_fill_class_parameters(cosmo,&fc,parser_length, status);

  if (*status != CCL_ERROR_CLASS)
    ccl_run_class(cosmo, &fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,init_arr,status);

  if (*status == CCL_ERROR_CLASS) {
    //printed error message while running CLASS
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }
  if (parser_free(&fc)== _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): Error freeing CLASS parser\n");
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    return;
  }
  FILE *f;
  f = fopen(filename,"w");
  if (!f){
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): Error opening output file\n");
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    fclose(f);
    return;
  }
  double psout_l,ic;
  int s=0;
  for (int i=0; i<sp.ln_k_size; i++) {
    s |= spectra_pk_at_k_and_z(&ba, &pm, &sp,exp(sp.ln_k[i]),z, &psout_l,&ic);
    fprintf(f,"%e %e\n",exp(sp.ln_k[i]),psout_l);
  }
  fclose(f);
  if(s) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: write_power_class_z(): Error writing CLASS power spectrum\n");
  }

  ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
}

typedef struct {
  double rsound;
  double zeq;
  double keq;
  double zdrag;
  double kSilk;
  double rsound_approx;
  double th2p7;
  double alphac;
  double alphab;
  double betac;
  double betab;
  double bnode;
} eh_struct;

static eh_struct *eh_struct_new(ccl_parameters *params)
{
  //////
  // Computes Eisenstein & Hu parameters for
  // P_k and r_sound
  // see astro-ph/9709112 for the relevant equations
  double OMh2,OBh2;
  double th2p7;
  eh_struct *eh=malloc(sizeof(eh_struct));
  if(eh==NULL)
    return NULL;

  OMh2=params->Omega_m*params->h*params->h; //Cosmo params scaled by h^2
  OBh2=params->Omega_b*params->h*params->h;
  th2p7=params->T_CMB/2.7; //This is Theta_{2.7} in E&Hu notation
  eh->th2p7=th2p7; //This is Theta_{2.7} in E&Hu notation
  eh->zeq=2.5E4*OMh2/pow(th2p7,4);// Eq 2
  eh->keq=0.0746*OMh2/(params->h*th2p7*th2p7);//Eq. 3

  //These group corresponds to Eq. 4
  double b1,b2;
  b1=0.313*pow(OMh2,-0.419)*(1+0.607*pow(OMh2,0.674));
  b2=0.238*pow(OMh2,0.223);
  eh->zdrag=1291*pow(OMh2,0.251)*(1+b1*pow(OBh2,b2))/(1+0.659*pow(OMh2,0.828));

  //These are the baryon-to-photon ratios
  //at equality (Req) and drag (Rd) epochs
  //Eq. 5
  double Req,Rd;
  Req=31.5*OBh2*1000./(eh->zeq*pow(th2p7,4));
  Rd=31.5*OBh2*1000./((1+eh->zdrag)*pow(th2p7,4));
  eh->rsound=2/(3*eh->keq)*sqrt(6/Req)*
    log((sqrt(1+Rd)+sqrt(Rd+Req))/(1+sqrt(Req)));

  //This is Eq. 7 (but in h/Mpc)
  eh->kSilk=1.6*pow(OBh2,0.52)*pow(OMh2,0.73)*(1+pow(10.4*OMh2,-0.95))/params->h;

  //These are Eqs. 11
  double a1,a2,b_frac;
  a1=pow(46.9*OMh2,0.670)*(1+pow(32.1*OMh2,-0.532));
  a2=pow(12.0*OMh2,0.424)*(1+pow(45.0*OMh2,-0.582));
  b_frac=OBh2/OMh2;
  eh->alphac=pow(a1,-b_frac)*pow(a2,-b_frac*b_frac*b_frac);

  //These are Eqs. 12
  double bb1,bb2;
  bb1=0.944/(1+pow(458*OMh2,-0.708));
  bb2=pow(0.395*OMh2,-0.0266);
  eh->betac=1/(1+bb1*(pow(1-b_frac,bb2)-1));

  double y=eh->zeq/(1+eh->zdrag);
  double sqy=sqrt(1+y);
  double gy=y*(-6*sqy+(2+3*y)*log((sqy+1)/(sqy-1))); //Eq 15
  //Baryon suppression Eq. 14
  eh->alphab=2.07*eh->keq*eh->rsound*pow(1+Rd,-0.75)*gy;

  //Baryon envelope shift Eq. 24
  eh->betab=0.5+b_frac+(3-2*b_frac)*sqrt(pow(17.2*OMh2,2)+1);

  //Node shift parameter Eq. 23
  eh->bnode=8.41*pow(OMh2,0.435);

  //Approx for the sound horizon, Eq. 26
  eh->rsound_approx=params->h*44.5*log(9.83/OMh2)/
    sqrt(1+10*pow(OBh2,0.75));

  return eh;
}

static double tkEH_0(double keq,double k,double a,double b)
{
  //////
  // Eisentstein & Hu's Tk_0
  // see astro-ph/9709112 for the relevant equations
  double q=k/(13.41*keq); //Eq 10
  double c=14.2/a+386./(1+69.9*pow(q,1.08)); //Eq 20
  double l=log(M_E+1.8*b*q); //Change of var for Eq 19
  return l/(l+c*q*q); //Returns Eq 19
}

static double tkEH_c(eh_struct *eh,double k)
{
  //////
  // Eisenstein & Hu's Tk_c
  // see astro-ph/9709112 for the relevant equations
  double f=1/(1+pow(k*eh->rsound/5.4,4)); //Eq 18
  return f*tkEH_0(eh->keq,k,1,eh->betac)+
    (1-f)*tkEH_0(eh->keq,k,eh->alphac,eh->betac); //Returns Eq 17
}

static double jbes0(double x)
{
  double jl;
  double ax2=x*x;

  if(ax2<1e-4) jl=1-ax2*(1-ax2/20.)/6.;
  else jl=sin(x)/x;

  return jl;
}

static double tkEH_b(eh_struct *eh,double k)
{
  //////
  // Eisenstein & Hu's Tk_b (Eq 21)
  // see astro-ph/9709112 for the relevant equations
  double x_bessel,part1,part2,jbes;
  double x=k*eh->rsound;

  //First term of Eq. 21
  if(k==0) x_bessel=0;
  else {
    x_bessel=x*pow(1+eh->bnode*eh->bnode*eh->bnode/(x*x*x),-1./3.);
  }

  part1=tkEH_0(eh->keq,k,1,1)/(1+pow(x/5.2,2));

  //Second term of Eq. 21
  if(k==0)
    part2=0;
  else
    part2=eh->alphab/(1+pow(eh->betab/x,3))*exp(-pow(k/eh->kSilk,1.4));

  return jbes0(x_bessel)*(part1+part2);
}

static double tsqr_EH(ccl_parameters *params,eh_struct * eh,double k,int wiggled)
{
  //////
  // Eisenstein & Hu's Tk_c
  // see astro-ph/9709112 for the relevant equations
  // Notice the last parameter in eh_power controls
  // whether to introduce wiggles (BAO) in the power spectrum.
  // We do this by default when obtaining the power spectrum.
  double tk;
  double b_frac=params->Omega_b/params->Omega_m;
  if(wiggled)
    //Case with baryons (Eq 8)
    tk=b_frac*tkEH_b(eh,k)+(1-b_frac)*tkEH_c(eh,k);
  else {
    //Zero baryon case (sec 4.2)
    double OMh2=params->Omega_m*params->h*params->h;
    // Compute Eq. 31
    double alpha_gamma=1-0.328*log(431*OMh2)*b_frac+0.38*log(22.3*OMh2)*b_frac*b_frac;
    // Compute Eq. 30
    double gamma_eff=params->Omega_m*params->h*(alpha_gamma+(1-alpha_gamma)/
						(1+pow(0.43*k*eh->rsound_approx,4)));
    // Compute Eq. 28 (assume k in h/Mpc)
    double q=k*eh->th2p7*eh->th2p7/gamma_eff;
    // Compute Eqs. 29
    double l0=log(2*M_E+1.8*q);
    double c0=14.2+731/(1+62.5*q);
    tk=l0/(l0+c0*q*q); //T_0 of Eq. 29
  }

  return tk*tk; //Return T_0^2
}

static double eh_power(ccl_parameters *params,eh_struct *eh,double k,int wiggled)
{
  //Wavenumber in units of Mpc^-1
  double kinvh=k/params->h; //Changed to h/Mpc
  return pow(k,params->n_s)*tsqr_EH(params,eh,kinvh,wiggled);
}
static void ccl_cosmology_compute_power_eh(ccl_cosmology * cosmo, int * status)
{
  //These are the limits of the splining range
  cosmo->data.k_min_lin = ccl_splines->K_MIN;
  cosmo->data.k_min_nl = ccl_splines->K_MIN;
  cosmo->data.k_max_lin = ccl_splines->K_MAX;
  cosmo->data.k_max_nl = ccl_splines->K_MAX;
  double kmin = cosmo->data.k_min_lin;
  double kmax = ccl_splines->K_MAX;

  // Compute nk from number of decades and N_K = # k per decade
  double ndecades = log10(kmax) - log10(kmin);
  int nk = (int)ceil(ndecades*ccl_splines->N_K);

  // Compute na using predefined spline spacing
  double amin = ccl_splines->A_SPLINE_MINLOG_PK;
  double amax = ccl_splines->A_SPLINE_MAX;
  int na = ccl_splines->A_SPLINE_NA_PK + ccl_splines->A_SPLINE_NLOG_PK - 1;

  // Exit if sigma8 wasn't specified
  if (isnan(cosmo->params.sigma8)) {
    *status = CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_eh(): sigma8 was not set, but is required for E&H method\n");
    return;
  }

  // New struct for EH parameters
  eh_struct *eh = eh_struct_new(&(cosmo->params));
  if (eh == NULL) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_eh(): memory allocation error\n");
    return;
  }

  // Build grids in k and a that P(k, a) will be evaluated on
  // NB: The x array is initially k, but will later be overwritten with log(k)
  double * x = ccl_log_spacing(kmin, kmax, nk);
  double * y = malloc(sizeof(double)*nk);
  double * a = ccl_linlog_spacing(amin, ccl_splines->A_SPLINE_MIN_PK,
                                  amax, ccl_splines->A_SPLINE_NLOG_PK,
                                  ccl_splines->A_SPLINE_NA_PK);
  double * y2d = malloc(nk * na * sizeof(double));
  if (a==NULL || y==NULL || x==NULL || y2d==NULL) {
    free(eh);free(x);free(y);
    free(a);free(y2d);
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_eh(): memory allocation error\n");
    return;
  }

  // Calculate P(k) on k grid. After this loop, x will contain log(k) and y
  // will contain log(pk) [which has not yet been normalized]
  // Notice the last parameter in eh_power controls
  // whether to introduce wiggles (BAO) in the power spectrum.
  // We do this by default.
  for (int i=0; i<nk; i++) {
    y[i] = log(eh_power(&cosmo->params, eh, x[i], 1));
    x[i] = log(x[i]);
  }

  // Apply growth factor, D(a), to P(k) and store in 2D (k, a) array
  double gfac, g2;
  gsl_spline2d *log_power_lin = gsl_spline2d_alloc(PLIN_SPLINE_TYPE, nk,na);
  for (int j = 0; j < na; j++) {
    gfac = ccl_growth_factor(cosmo, a[j], status);
    g2 = 2.*log(gfac);
    for (int i=0; i<nk; i++) {
      y2d[j*nk+i] = y[i]+g2;
    } // end loop over k
  } // end loop over a

  // Check that ccl_growth_factor didn't fail
  if (*status) {
    free(eh); free(x); free(y); free(a); free(y2d);
    gsl_spline2d_free(log_power_lin);
    return;
  }

  // Initialize a 2D spline over P(k, a) [which is still unnormalized by sigma8]
  int splinstatus = gsl_spline2d_init(log_power_lin, x, a, y2d,nk,na);
  if (splinstatus) {
    free(eh); free(x); free(y); free(a); free(y2d);
    gsl_spline2d_free(log_power_lin);
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_eh(): Error creating log_power_lin spline\n");
    return;
  }

  // Calculate sigma8 for the unnormalized P(k), using the standard
  // ccl_sigma8() function
  cosmo->data.p_lin = log_power_lin;
  cosmo->computed_power = true; // Temporarily set this to true
  double sigma8 = ccl_sigma8(cosmo, status);
  cosmo->computed_power = false;

  // Check that ccl_sigma8 didn't fail
  if (*status) {
    free(eh); free(x); free(y); free(a); free(y2d);
    gsl_spline2d_free(log_power_lin);
    return;
  }

  // Calculate normalization factor using computed value of sigma8, then
  // recompute P(k, a) using this normalization
  double log_normalization_factor = 2*(log(cosmo->params.sigma8) - log(sigma8));
  for (int i=0; i < nk; i++) {
    y[i] += log_normalization_factor;
  }
  for (int j = 0; j < na; j++) {
    gfac = ccl_growth_factor(cosmo, a[j], status);
    g2 = 2.*log(gfac);
    for (int i=0; i<nk; i++) {
      y2d[j*nk+i] = y[i]+g2; // Replace previous values
    } // end k loop
  } // end a loop

  splinstatus = gsl_spline2d_init(log_power_lin, x, a, y2d, nk, na);
  if (splinstatus) {
    free(eh); free(x); free(y); free(a); free(y2d);
    gsl_spline2d_free(log_power_lin);
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_eh(): Error creating log_power_lin spline\n");
    return;
  }
  else
    cosmo->data.p_lin = log_power_lin;

  // Allocate a 2D spline for the nonlinear P(k) [which is just a copy of the
  // linear one for E&H]
  gsl_spline2d * log_power_nl = gsl_spline2d_alloc(PNL_SPLINE_TYPE, nk, na);
  splinstatus = gsl_spline2d_init(log_power_nl, x, a, y2d, nk, na);

  if (splinstatus) {
    free(eh); free(x); free(y); free(a); free(y2d);
    gsl_spline2d_free(log_power_lin);
    gsl_spline2d_free(log_power_nl);
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_eh(): Error creating log_power_nl spline\n");
    return;
  }
  else
    cosmo->data.p_nl = log_power_nl;

  // Free temporary arrays
  free(eh); free(x); free(y); free(a); free(y2d);
}

/*------ ROUTINE: tsqr_BBKS -----
INPUT: ccl_parameters and k wavenumber in 1/Mpc
TASK: provide the square of the BBKS transfer function with baryonic correction
NOTE: Bardeen et al. (1986) as implemented in Sugiyama (1995)
*/

static double tsqr_BBKS(ccl_parameters * params, double k)
{
  double q = k/(params->Omega_m*params->h*params->h*exp(-params->Omega_b*(1.0+pow(2.*params->h,.5)/params->Omega_m)));
  return pow(log(1.+2.34*q)/(2.34*q),2.0)/pow(1.+3.89*q+pow(16.1*q,2.0)+pow(5.46*q,3.0)+pow(6.71*q,4.0),0.5);
}

/*------ ROUTINE: bbks_power -----
INPUT: ccl_parameters and k wavenumber in 1/Mpc
TASK: provide the BBKS power spectrum with baryonic correction at single k
*/

//Calculate Normalization see Cosmology Notes 8.105 (TODO: whose comment is this?)
static double bbks_power(ccl_parameters * params, double k)
{
  return pow(k,params->n_s)*tsqr_BBKS(params, k);
}

/*------ ROUTINE: ccl_cosmology_compute_bbks_power -----
INPUT: cosmology
TASK: provide spline for the BBKS power spectrum with baryonic correction
*/

static void ccl_cosmology_compute_power_bbks(ccl_cosmology * cosmo, int * status)
{
  //These are the limits of the splining range
  cosmo->data.k_min_lin=ccl_splines->K_MIN;
  cosmo->data.k_min_nl=ccl_splines->K_MIN;
  cosmo->data.k_max_lin=ccl_splines->K_MAX;
  cosmo->data.k_max_nl=ccl_splines->K_MAX;
  double kmin = cosmo->data.k_min_lin;
  double kmax = ccl_splines->K_MAX;
  //Compute nk from number of decades and N_K = # k per decade
  double ndecades = log10(kmax) - log10(kmin);
  int nk = (int)ceil(ndecades*ccl_splines->N_K);
  double amin = ccl_splines->A_SPLINE_MINLOG_PK;
  double amax = ccl_splines->A_SPLINE_MAX;
  int na = ccl_splines->A_SPLINE_NA_PK+ccl_splines->A_SPLINE_NLOG_PK-1;

  // Exit if sigma8 wasn't specified
  if (isnan(cosmo->params.sigma8)) {
    *status = CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_bbks(): sigma8 not set, required for BBKS\n");
    return; //Return without anything to free
  }

  // The x array is initially k, but will later
  // be overwritten with log(k)
  double * x = ccl_log_spacing(kmin, kmax, nk);
  double * y = malloc(sizeof(double)*nk);
  double * a = ccl_linlog_spacing(amin, ccl_splines->A_SPLINE_MIN_PK, amax, ccl_splines->A_SPLINE_NLOG_PK, ccl_splines->A_SPLINE_NA_PK);
  double * y2d = malloc(nk * na * sizeof(double));

  //If error, store status, we will free later
  if (a==NULL||y==NULL|| x==NULL || y2d==0) {
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_bbks(): memory allocation error\n");
  }

  //Status flags
  int splinstatus=0;

  if(!*status){

    // After this loop x will contain log(k)
    for (int i=0; i<nk; i++) {
      y[i] = log(bbks_power(&cosmo->params, x[i]));
      x[i] = log(x[i]);
    }

    for (int j = 0; j < na; j++) {
      double gfac = ccl_growth_factor(cosmo,a[j], status);
      double g2 = 2.*log(gfac);
      for (int i=0; i<nk; i++) {
	y2d[j*nk+i] = y[i]+g2;
      }
    }
  }

  if(!*status){

    // Initialize a 2D spline over P(k, a) [which is still unnormalized by sigma8]
    gsl_spline2d * log_power_lin_unnorm = gsl_spline2d_alloc(PLIN_SPLINE_TYPE, nk,na);
    splinstatus = gsl_spline2d_init(log_power_lin_unnorm, x, a, y2d,nk,na);

    //If not, proceed
    if (!splinstatus) {
      cosmo->data.p_lin=log_power_lin_unnorm;
      cosmo->computed_power=true;
    } else {
      gsl_spline2d_free(log_power_lin_unnorm);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo,"ccl_power.c: ccl_cosmology_compute_power_bbks(): Error creating log_power_lin spline\n");
    }
  }

  if(*status){
    free(x); free(y); free(a); free(y2d);
    return;
  }

  // Calculate sigma8 for the unnormalized P(k), using the standard
  // ccl_sigma8() function
  double sigma8 = ccl_sigma8(cosmo,status);
  cosmo->computed_power=false;

  // Check that ccl_sigma8 didn't fail
  if (!*status) {

    // Calculate normalization factor using computed value of sigma8, then
    // recompute P(k, a) using this normalization
    double log_normalization_factor = 2*(log(cosmo->params.sigma8) - log(sigma8));
    for (int i=0; i<nk; i++) {
      y[i] += log_normalization_factor;
    }
    for (int j = 0; j < na; j++) {
      double gfac = ccl_growth_factor(cosmo,a[j], status);
      double g2 = 2.*log(gfac);
      for (int i=0; i<nk; i++) {
	y2d[j*nk+i] = y[i]+g2;
      }
    }
  }

  //Check growth didn't fail
  if (!*status) {

    gsl_spline2d * log_power_lin = gsl_spline2d_alloc(PLIN_SPLINE_TYPE, nk,na);
    splinstatus = gsl_spline2d_init(log_power_lin, x, a, y2d,nk,na);

    if (!splinstatus) {
      cosmo->data.p_lin=log_power_lin;
    } else {
      gsl_spline2d_free(log_power_lin);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo,"ccl_power.c: ccl_cosmology_compute_power_bbks(): Error creating log_power_lin spline\n");
    }
  }

  if(!*status){
    // Allocate a 2D spline for the nonlinear P(k)
    //[which is just a copy of the linear one for BBKS]
    gsl_spline2d * log_power_nl = gsl_spline2d_alloc(PNL_SPLINE_TYPE, nk,na);
    splinstatus = gsl_spline2d_init(log_power_nl, x, a, y2d,nk,na);
    if (!splinstatus) {
      cosmo->data.p_nl = log_power_nl;
    } else {
      gsl_spline2d_free(log_power_nl);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_bbks(): Error creating log_power_nl spline\n");
    }
  }

  free(x); free(y); free(a); free(y2d);
  return;
}

/*------ ROUTINE: ccl_cosmology_compute_power_emu -----
INPUT: cosmology
TASK: provide spline for the emulated power spectrum from Cosmic EMU
*/

static void ccl_cosmology_compute_power_emu(ccl_cosmology * cosmo, int * status)
{

  struct precision pr;        // for precision parameters
  struct background ba;       // for cosmological background
  struct thermo th;           // for thermodynamics
  struct perturbs pt;         // for source functions
  struct transfers tr;        // for transfer functions
  struct primordial pm;       // for primordial spectra
  struct spectra sp;          // for output spectra
  struct nonlinear nl;        // for non-linear spectra
  struct lensing le;
  struct output op;
  struct file_content fc;

  double Omeganuh2_eq;
  double w0wacomb = -cosmo->params.w0 - cosmo->params.wa;

  ErrorMsg errmsg; // for error messages
  // generate file_content structure
  // Configuration parameters will be passed through this structure,
  // to avoid writing and reading .ini files for every call
  int parser_length = 20;
  int init_arr[7]={0,0,0,0,0,0,0};

  //Check initialization
  if (parser_init(&fc,parser_length,"none",errmsg) == _FAILURE_) {
    *status = CCL_ERROR_CLASS;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): parser init error:%s\n", errmsg);
  }

  // Check ranges to see if the cosmology is valid
  if((cosmo->params.h<0.55) || (cosmo->params.h>0.85)){
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): h is outside allowed range\n");
  }
  if(w0wacomb<8.1e-3){ //0.3^4
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): w0 and wa do not satisfy the emulator bound\n");
  }
  if(cosmo->params.Omega_n_mass*cosmo->params.h*cosmo->params.h>0.01){
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): Omega_nu does not satisfy the emulator bound\n");
  }

  // Check to see if sigma8 was defined
  if(isnan(cosmo->params.sigma8)){
    *status=CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): sigma8 is not defined; specify sigma8 instead of A_s\n");
  }

  //If one of the previous test was unsuccessful, quit:
  if(*status) return;

  // Check if the cosmology has been set up with equal neutrino masses for the emulator
  // If not, check if the user has forced redistribution of masses and if so do this.
  if(cosmo->params.N_nu_mass>0) {
    if (cosmo->config.emulator_neutrinos_method == ccl_emu_strict){
      if (cosmo->params.N_nu_mass==3){
	//double diff1 = pow((cosmo->params.mnu[0] - cosmo->params.mnu[1]) * (cosmo->params.mnu[0] - cosmo->params.mnu[1]), 0.5);
	//double diff2 = pow((cosmo->params.mnu[1] - cosmo->params.mnu[2]) * (cosmo->params.mnu[1] - cosmo->params.mnu[2]), 0.5);
	//double diff3 = pow((cosmo->params.mnu[2] - cosmo->params.mnu[0]) * (cosmo->params.mnu[2] - cosmo->params.mnu[0]), 0.5);
	//if (diff1>1e-12 || diff2>1e-12 || diff3>1e-12){
	if (cosmo->params.mnu[0] != cosmo->params.mnu[1] || cosmo->params.mnu[0] != cosmo->params.mnu[2] || cosmo->params.mnu[1] != cosmo->params.mnu[2]){
	  *status = CCL_ERROR_INCONSISTENT;
	  ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): In the default configuration, you must pass a list of 3 equal neutrino masses or pass a sum and set mnu_type = ccl_mnu_sum_equal. If you wish to over-ride this, set config->emulator_neutrinos_method = 'ccl_emu_equalize'. This will force the neutrinos to be of equal mass but will result in internal inconsistencies.\n");
	}
      } else if (cosmo->params.N_nu_mass!=3){
	*status = CCL_ERROR_INCONSISTENT;
	ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): In the default configuration, you must pass a list of 3 equal neutrino masses or pass a sum and set mnu_type = ccl_mnu_sum_equal. If you wish to over-ride this, set config->emulator_neutrinos_method = 'ccl_emu_equalize'. This will force the neutrinos to be of equal mass but will result in internal inconsistencies.\n");
      }
    } else if (cosmo->config.emulator_neutrinos_method == ccl_emu_equalize){
      // Reset the masses to equal
      double mnu_eq[3] = {cosmo->params.sum_nu_masses / 3., cosmo->params.sum_nu_masses / 3., cosmo->params.sum_nu_masses / 3.};
      Omeganuh2_eq = ccl_Omeganuh2(1.0, 3, mnu_eq, cosmo->params.T_CMB, cosmo->data.accelerator, status);
    }
  } else {
    if(fabs(cosmo->params.N_nu_rel - 3.04)>1.e-6){
      *status=CCL_ERROR_INCONSISTENT;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): Set Neff = 3.04 for cosmic emulator predictions in absence of massive neutrinos.\n");
    }
  }

  if(!*status){
    // Prepare to run CLASS for linear scales
    ccl_fill_class_parameters(cosmo,&fc,parser_length, status);
  }

  if (!*status){
    ccl_run_class(cosmo, &fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,init_arr,status);
  }

  if (*status){
    ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
    if (parser_free(&fc)== _FAILURE_) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): Error freeing CLASS parser\n");
    }
    return;
  }

  //These are the limits of the splining range
  cosmo->data.k_min_lin=2*exp(sp.ln_k[0]);
  cosmo->data.k_max_lin=ccl_splines->K_MAX_SPLINE;
  //CLASS calculations done - now allocate CCL splines
  double kmin = cosmo->data.k_min_lin;
  double kmax = ccl_splines->K_MAX_SPLINE;
  //Compute nk from number of decades and N_K = # k per decade
  double ndecades = log10(kmax) - log10(kmin);
  int nk = (int)ceil(ndecades*ccl_splines->N_K);
  double amin = ccl_splines->A_SPLINE_MINLOG_PK;
  double amax = ccl_splines->A_SPLINE_MAX;
  int na = ccl_splines->A_SPLINE_NA_PK+ccl_splines->A_SPLINE_NLOG_PK-1;

  // The x array is initially k, but will later
  // be overwritten with log(k)
  double * x = ccl_log_spacing(kmin, kmax, nk);
  double * a = ccl_linlog_spacing(amin, ccl_splines->A_SPLINE_MIN_PK, amax, ccl_splines->A_SPLINE_NLOG_PK, ccl_splines->A_SPLINE_NA_PK);
  double * y2d_lin = malloc(nk * na * sizeof(double));
  if (a==NULL|| x==NULL || y2d_lin==NULL) {
    *status = CCL_ERROR_SPLINE;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_class(): memory allocation error\n");
  }

  if(!*status){
    // After this loop x will contain log(k), y will contain log(P_nl), z will contain log(P_lin)
    // all in Mpc, not Mpc/h units!
    double psout_l,ic;
    int s=0;
    for (int i=0; i<nk; i++) {
      for (int j = 0; j < na; j++) {
	//The 2D interpolation routines access the function values y_{k_ia_j} with the following ordering:
	//y_ij = y2d[j*N_k + i]
	//with i = 0,...,N_k-1 and j = 0,...,N_a-1.
	s |= spectra_pk_at_k_and_z(&ba, &pm, &sp,x[i],1./a[j]-1., &psout_l,&ic);
	y2d_lin[j*nk+i] = log(psout_l);
      }
      x[i] = log(x[i]);
    }
    if(s) {
      *status = CCL_ERROR_CLASS;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): Error computing CLASS power spectrum\n");
    }
  }

  if(!*status){

    gsl_spline2d * log_power = gsl_spline2d_alloc(PLIN_SPLINE_TYPE, nk,na);
    int pwstatus = gsl_spline2d_init(log_power, x, a, y2d_lin,nk,na);
    if (!pwstatus) {
      cosmo->data.p_lin = log_power;
    } else {
      gsl_spline2d_free(log_power);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): Error creating log_power spline\n");
    }
  }

  if(*status){ //Linear power spectrum failed, so quit.
      free(x);
      free(a);
      free(y2d_lin);
      ccl_free_class_structs(cosmo, &ba,&th,&pt,&tr,&pm,&sp,&nl,&le,init_arr,status);
      return;
  }

  //Now start the NL computation with the emulator
  //These are the limits of the splining range
  cosmo->data.k_min_nl=K_MIN_EMU;
  cosmo->data.k_max_nl=K_MAX_EMU;
  amin = A_MIN_EMU; //limit of the emulator
  amax = ccl_splines->A_SPLINE_MAX;
  na = ccl_splines->A_SPLINE_NA_PK;
  // The x array is initially k, but will later
  // be overwritten with log(k)
  double * logx= malloc(NK_EMU*sizeof(double));
  double * y;
  double * xstar = malloc(9 * sizeof(double));
  double * aemu = ccl_linear_spacing(amin,amax, na);
  double * y2d = malloc(NK_EMU * na * sizeof(double));
  if (aemu==NULL || y2d==NULL || logx==NULL || xstar==NULL){
    *status=CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): memory allocation error\n");
  }

  if(!*status){

    //For each redshift:
    for (int j = 0; j < na; j++){

      //Turn cosmology into xstar:
      xstar[0] = (cosmo->params.Omega_c+cosmo->params.Omega_b)*cosmo->params.h*cosmo->params.h;
      xstar[1] = cosmo->params.Omega_b*cosmo->params.h*cosmo->params.h;
      xstar[2] = cosmo->params.sigma8;
      xstar[3] = cosmo->params.h;
      xstar[4] = cosmo->params.n_s;
      xstar[5] = cosmo->params.w0;
      xstar[6] = cosmo->params.wa;
      if ((cosmo->params.N_nu_mass>0) && (cosmo->config.emulator_neutrinos_method == ccl_emu_equalize)){
	xstar[7] = Omeganuh2_eq;
      } else {
        xstar[7] = cosmo->params.Omega_n_mass*cosmo->params.h*cosmo->params.h;
      }
      xstar[8] = 1./aemu[j]-1;
      //Need to have this here because otherwise overwritten by emu in each loop

      //Call emulator at this redshift
      ccl_pkemu(xstar,&y, status, cosmo);
      if (y == NULL){
	*status = CCL_ERROR_MEMORY;
	ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): Error obtaining emulator predictions\n");
      }
      for (int i=0; i<NK_EMU; i++){
	logx[i] = log(mode[i]);
	y2d[j*NK_EMU+i] = log(y[i]);
      }
    }
  }

  if(!*status){

    gsl_spline2d * log_power_nl = gsl_spline2d_alloc(PLIN_SPLINE_TYPE, NK_EMU,na);
    int splinstatus = gsl_spline2d_init(log_power_nl, logx, aemu, y2d,NK_EMU,na);
    //Note the minimum k of the spline is different from the linear one.

    if (!splinstatus){
      cosmo->data.p_nl = log_power_nl;
    } else {
      gsl_spline2d_free(log_power_nl);
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power_emu(): Error creating log_power spline\n");
    }
  }

  free(x); free(a); free(y);
  free(xstar); free(logx); free(aemu);
  free(y2d_lin); free(y2d);
}

/*------ ROUTINE: ccl_cosmology_compute_power -----
INPUT: ccl_cosmology * cosmo
TASK: compute power spectrum
*/
void ccl_cosmology_compute_power(ccl_cosmology * cosmo, int * status)
{

  if (cosmo->computed_power) return;
    switch(cosmo->config.transfer_function_method){
        case ccl_bbks:
	  ccl_cosmology_compute_power_bbks(cosmo,status);
	  break;
        case ccl_eisenstein_hu:
	  ccl_cosmology_compute_power_eh(cosmo,status);
	  break;
        case ccl_boltzmann_class:
	  ccl_cosmology_compute_power_class(cosmo,status);
	  break;
        case ccl_emulator:
	  ccl_cosmology_compute_power_emu(cosmo,status);
	  break;
        default:
	  *status = CCL_ERROR_INCONSISTENT;
	  ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_cosmology_compute_power(): Unknown or non-implemented transfer function method: %d \n", cosmo->config.transfer_function_method);
    }

    ccl_check_status(cosmo,status);
    if (*status==0){
      cosmo->computed_power = true;
    }
  return;
}

/*------ ROUTINE: ccl_power_extrapol_highk -----
INPUT: ccl_cosmology * cosmo, a, k [1/Mpc]
TASK: extrapolate power spectrum at high k
*/
static double ccl_power_extrapol_highk(ccl_cosmology * cosmo, double k, double a,
				       gsl_spline2d * powerspl, double kmax_spline, int * status)
{
  double log_p_1;
  double deltak=1e-2; //step for numerical derivative;
  double deriv_pk_kmid,deriv2_pk_kmid;
  double lkmid;
  double lpk_kmid;

  lkmid = log(kmax_spline)-2*deltak;

  int gslstatus =  gsl_spline2d_eval_e(powerspl, lkmid,a,NULL ,NULL ,&lpk_kmid);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_power_extrapol_highk():");
    *status = CCL_ERROR_SPLINE_EV;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_power_extrapol_highk(): Spline evaluation error\n");
    return NAN;
  }
  //GSL derivatives
  gslstatus = gsl_spline2d_eval_deriv_x_e (powerspl, lkmid, a, NULL,NULL,&deriv_pk_kmid);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_power_extrapol_highk():");
    *status = CCL_ERROR_SPLINE_EV;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_power_extrapol_highk(): Spline evaluation error\n");
    return NAN;
  }
  gslstatus = gsl_spline2d_eval_deriv_xx_e (powerspl, lkmid, a, NULL,NULL,&deriv2_pk_kmid);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_power_extrapol_highk():");
    *status = CCL_ERROR_SPLINE_EV;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_power_extrapol_highk(): Spline evaluation error\n");
    return NAN;
  }
  log_p_1=lpk_kmid+deriv_pk_kmid*(log(k)-lkmid)+deriv2_pk_kmid/2.*(log(k)-lkmid)*(log(k)-lkmid);

  return log_p_1;

}

/*------ ROUTINE: ccl_power_extrapol_lowk -----
INPUT: ccl_cosmology * cosmo, a, k [1/Mpc]
TASK: extrapolate power spectrum at low k
*/
static double ccl_power_extrapol_lowk(ccl_cosmology * cosmo, double k, double a,
				      gsl_spline2d * powerspl, double kmin_spline, int * status)
{
  double log_p_1;
  double deltak=1e-2; //safety step
  double lkmin=log(kmin_spline)+deltak;
  double lpk_kmin;
  int gslstatus = gsl_spline2d_eval_e(powerspl,lkmin,a,NULL,NULL,&lpk_kmin);

  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_power_extrapol_lowk():");
    *status=CCL_ERROR_SPLINE_EV;
    ccl_cosmology_set_status_message(cosmo,"ccl_power.c: ccl_power_extrapol_lowk(): Spline evaluation error\n");
    return NAN;
  }

  return lpk_kmin+cosmo->params.n_s*(log(k)-lkmin);
}

/*------ ROUTINE: ccl_linear_matter_power -----
INPUT: ccl_cosmology * cosmo, k [1/Mpc],a
TASK: compute the linear power spectrum at a given redshift
      by rescaling using the growth function
*/

double ccl_linear_matter_power(ccl_cosmology * cosmo, double k, double a, int * status)

{
  if ((cosmo->config.transfer_function_method == ccl_emulator) && (a<A_MIN_EMU)){
    *status = CCL_ERROR_INCONSISTENT;
    ccl_cosmology_set_status_message(cosmo, "ccl_power.c: the cosmic emulator cannot be used above a=%f\n",A_MIN_EMU);
    return NAN;
  }

  if (!cosmo->computed_power) ccl_cosmology_compute_power(cosmo, status);
  // Return if compilation failed
  //if (cosmo->data.p_lin == NULL) return NAN;
  if (!cosmo->computed_power) return NAN;

  double log_p_1;
  int gslstatus;

  if(a<ccl_splines->A_SPLINE_MINLOG_PK) {  //Extrapolate linearly at high redshift
    double pk0=ccl_linear_matter_power(cosmo,k,ccl_splines->A_SPLINE_MINLOG_PK,status);
    double gf=ccl_growth_factor(cosmo,a,status)/ccl_growth_factor(cosmo,ccl_splines->A_SPLINE_MINLOG_PK,status);

    return pk0*gf*gf;
  }
  if (*status!=CCL_ERROR_INCONSISTENT){
    if(k<=cosmo->data.k_min_lin) {
      log_p_1=ccl_power_extrapol_lowk(cosmo,k,a,cosmo->data.p_lin,cosmo->data.k_min_lin,status);

      return exp(log_p_1);
    }
    else if(k<cosmo->data.k_max_lin){
      gslstatus = gsl_spline2d_eval_e(cosmo->data.p_lin, log(k), a,NULL,NULL,&log_p_1);
      if(gslstatus != GSL_SUCCESS) {
        ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_linear_matter_power():");
        *status = CCL_ERROR_SPLINE_EV;
        ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_linear_matter_power(): Spline evaluation error\n");
        return NAN;
      }
      else {
        return exp(log_p_1);
      }
    }
    else { //Extrapolate using log derivative
      log_p_1 = ccl_power_extrapol_highk(cosmo,k,a,cosmo->data.p_lin,cosmo->data.k_max_lin,status);
      return exp(log_p_1);
    }
  }

  return exp(log_p_1);
}

/*------ ROUTINE: ccl_nonlin_matter_power -----
INPUT: ccl_cosmology * cosmo, a, k [1/Mpc]
TASK: compute the nonlinear power spectrum at a given redshift
*/

double ccl_nonlin_matter_power(ccl_cosmology * cosmo, double k, double a, int *status)
{
  double log_p_1, pk;

  switch(cosmo->config.matter_power_spectrum_method) {

  //If the matter PS specified was linear, then do the linear compuation
  case ccl_linear:
    return ccl_linear_matter_power(cosmo,k,a,status);

  case ccl_halofit:
    if (!cosmo->computed_power)
      ccl_cosmology_compute_power(cosmo, status);
    if (cosmo->data.p_nl == NULL) return NAN; // Return if computation failed

    if(a<ccl_splines->A_SPLINE_MINLOG_PK) { //Extrapolate linearly at high redshift
      double pk0=ccl_nonlin_matter_power(cosmo,k,ccl_splines->A_SPLINE_MINLOG_PK,status);
      double gf=ccl_growth_factor(cosmo,a,status)/ccl_growth_factor(cosmo,ccl_splines->A_SPLINE_MINLOG_PK,status);
      return pk0*gf*gf;
    }
		break;

  case ccl_emu:
    if ((cosmo->config.transfer_function_method == ccl_emulator) && (a<A_MIN_EMU)){
      *status = CCL_ERROR_EMULATOR_BOUND;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: the cosmic emulator cannot be used above z=2\n");
      return NAN;
    }

    // Compute power spectrum if needed; return if computation failed
    if (!cosmo->computed_power){
      ccl_cosmology_compute_power(cosmo,status);
    }
    if (cosmo->data.p_nl == NULL) return NAN;
		break;

  default:
    ccl_raise_warning(
			CCL_ERROR_NOT_IMPLEMENTED,
      "config.matter_power_spectrum_method = %d not yet supported "
			"continuing with linear power spectrum\n", cosmo->config.matter_power_spectrum_method);
    cosmo->config.matter_power_spectrum_method=ccl_linear;
    return ccl_linear_matter_power(cosmo,k,a,status);
  } // end switch

  // if we get here, try to evaluate the power spectrum
  // we need to account for bounds below and above
  if (k <= cosmo->data.k_min_nl) {
    // we assume no baryonic effects below k_min_nl
    log_p_1 = ccl_power_extrapol_lowk(cosmo, k, a, cosmo->data.p_nl, cosmo->data.k_min_nl, status);
    return exp(log_p_1);
  }

  if (k < cosmo->data.k_max_nl) {
    int gslstatus = gsl_spline2d_eval_e(cosmo->data.p_nl, log(k), a, NULL ,NULL, &log_p_1);
    if (gslstatus != GSL_SUCCESS) {
      ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_nonlin_matter_power():");
      *status = CCL_ERROR_SPLINE_EV;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_nonlin_matter_power(): Spline evaluation error\n");
      return NAN;
    } else {
      pk = exp(log_p_1);
    }
  } else {
    // Extrapolate NL regime using log derivative
    log_p_1 = ccl_power_extrapol_highk(cosmo, k, a, cosmo->data.p_nl, cosmo->data.k_max_nl, status);
    pk = exp(log_p_1);
  }

	// Add baryonic correction
  if (cosmo->config.baryons_power_spectrum_method == ccl_bcm) {
    int pwstatus = 0;
    double fbcm = ccl_bcm_model_fka(cosmo, k, a, &pwstatus);
    pk *= fbcm;
    if (pwstatus) {
      *status = CCL_ERROR_SPLINE_EV;
      ccl_cosmology_set_status_message(cosmo, "ccl_power.c: ccl_nonlin_matter_power(): Error in BCM correction\n");
      return NAN;
    }
  }

  return pk;
}

// Params for sigma(R) integrand
typedef struct {
  ccl_cosmology *cosmo;
  double R;
  int* status;
} SigmaR_pars;

typedef struct {
  ccl_cosmology *cosmo;
  double R;
  int* status;
} SigmaV_pars;

/* --------- ROUTINE: w_tophat ---------
INPUT: kR, ususally a wavenumber multiplied by a smoothing radius
TASK: Output W(x)=[sin(x)-x*cos(x)]*(3/x)^3
*/
static double w_tophat(double kR)
{
  double w;
  double kR2 = kR*kR;

  // This is the Maclaurin expansion of W(x)=[sin(x)-xcos(x)]*(3/x)**3 to O(x^7), with x=kR.
  // Necessary numerically because at low x W(x) relies on the fine cancellation of two terms
  if(kR<0.1) {
     w = 1. + kR2*(-0.1 +
		   kR2*(0.003561429 +
			kR2*(-6.61376e-5 +
			     kR2*(7.51563e-7))));
     /*w =1.-0.1*kR*kR+0.003571429*kR*kR*kR*kR
      -6.61376E-5*kR*kR*kR*kR*kR*kR
      +7.51563E-7*kR*kR*kR*kR*kR*kR*kR*kR;*/
  }
  else
    w = 3.*(sin(kR) - kR*cos(kR))/(kR2*kR);
  return w;
}

// Integrand for sigmaR integral
static double sigmaR_integrand(double lk,void *params)
{
  SigmaR_pars *par=(SigmaR_pars *)params;

  double k=pow(10.,lk);
  double pk=ccl_linear_matter_power(par->cosmo,k, 1.,par->status);
  double kR=k*par->R;
  double w = w_tophat(kR);

  return pk*k*k*k*w*w;
}

// Integrand for sigmaV integral
static double sigmaV_integrand(double lk,void *params)
{
  SigmaV_pars *par=(SigmaV_pars *)params;

  double k=pow(10.,lk);
  double pk=ccl_linear_matter_power(par->cosmo,k, 1.,par->status);
  double kR=k*par->R;
  double w = w_tophat(kR);

  return pk*k*w*w/3.0;
}

/* --------- ROUTINE: ccl_sigmaR ---------
INPUT: cosmology, comoving smoothing radius, scale factor
TASK: compute sigmaR, the variance in the *linear* density field
smoothed with a tophat filter of comoving size R
*/
double ccl_sigmaR(ccl_cosmology *cosmo,double R,double a,int *status)
{
  SigmaR_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  gsl_integration_cquad_workspace *workspace=gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function=&sigmaR_integrand;
  F.params=∥
  double sigma_R;
  int gslstatus = gsl_integration_cquad(&F, log10(ccl_splines->K_MIN), log10(ccl_splines->K_MAX),
				                                0.0, ccl_gsl->INTEGRATION_SIGMAR_EPSREL,
                                        workspace,&sigma_R,NULL,NULL);
  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_sigmaR():");
    *status |= gslstatus;
  }

  gsl_integration_cquad_workspace_free(workspace);

  return sqrt(sigma_R*M_LN10/(2*M_PI*M_PI))*ccl_growth_factor(cosmo, a, status);
}

/* --------- ROUTINE: ccl_sigmaV ---------
INPUT: cosmology, comoving smoothing radius, scale factor
TASK: compute sigmaV, the variance in the *linear* displacement field
smoothed with a tophat filter of comoving size R
The linear displacement field is the gradient of the linear density field
*/
double ccl_sigmaV(ccl_cosmology *cosmo,double R,double a,int *status)
{
  SigmaV_pars par;
  par.status = status;

  par.cosmo=cosmo;
  par.R=R;
  gsl_integration_cquad_workspace *workspace=gsl_integration_cquad_workspace_alloc(ccl_gsl->N_ITERATION);
  gsl_function F;
  F.function=&sigmaV_integrand;
  F.params=∥
  double sigma_V;
	int gslstatus = gsl_integration_cquad(&F, log10(ccl_splines->K_MIN), log10(ccl_splines->K_MAX),
																				0.0, ccl_gsl->INTEGRATION_SIGMAR_EPSREL,
																				workspace,&sigma_V,NULL,NULL);

  if(gslstatus != GSL_SUCCESS) {
    ccl_raise_gsl_warning(gslstatus, "ccl_power.c: ccl_sigmaV():");
    *status |= gslstatus;
  }

  gsl_integration_cquad_workspace_free(workspace);

  return sqrt(sigma_V*M_LN10/(2*M_PI*M_PI))*ccl_growth_factor(cosmo, a, status);
}

/* --------- ROUTINE: ccl_sigma8 ---------
INPUT: cosmology
TASK: compute sigma8, the variance in the *linear* density field at a=1
smoothed with a tophat filter of comoving size 8 Mpc/h
*/
double ccl_sigma8(ccl_cosmology *cosmo, int *status)
{
  return ccl_sigmaR(cosmo,8/cosmo->params.h, 1.,status);
}

max_stars_repo_name: vrastil/CCL

max_issues_repo_name: vrastil/CCL

max_forks_repo_name: vrastil/CCL


text:
/**
 * @file matrix.i
 * @brief Wrapper for GSL matrix objects.
 * @author John McDonough, Kenichi Kumatani
 */

#ifndef GSLMATRIX_H
#define GSLMATRIX_H

#include <gsl/gsl_matrix.h>

#define GSL_MATRIX_NROWS(x) ((x)->size2)
#define GSL_MATRIX_NCOLS(x) ((x)->size2)

gsl_matrix_float* gsl_matrix_float_resize(gsl_matrix_float* m, size_t size1, size_t size2);

void gsl_matrix_float_set_cosine(gsl_matrix_float* m, size_t i, size_t j, int type);

gsl_matrix_float* gsl_matrix_float_load(gsl_matrix_float* m, const char* filename, bool old = false);

gsl_vector_float* gsl_vector_float_load(gsl_vector_float* m, const char* filename, bool old = false);

#endif

ext: h


text:
// Copyright 2011-2021 GameParadiso, Inc. All Rights Reserved.

#pragma once

#include <map>
#include <unordered_map>
#include <ranges>

#include <glm/vec3.hpp>
#include <glm/vec4.hpp>
#include <glm/mat4x4.hpp>
#include <glm/gtc/quaternion.hpp>
#include <glm/ext/matrix_transform.hpp>
#include <glm/ext/matrix_clip_space.hpp>
#include <glm/ext/scalar_constants.hpp>

#include <gsl/gsl>

using namespace std::string_literals;

max_stars_repo_name: noirhero/MemoryChunk

max_issues_repo_name: noirhero/MemoryChunk

max_forks_repo_name: noirhero/MemoryChunk


text:
///
/// @file
///
/// @author Mirko Myllykoski (mirkom@cs.umu.se), Umeå University
/// @author Lars Karlsson (larsk@cs.umu.se), Umeå University
///
/// @internal LICENSE
///
/// Copyright (c) 2019-2020, Umeå Universitet
///
/// Redistribution and use in source and binary forms, with or without
/// modification, are permitted provided that the following conditions are met:
///
/// 1. Redistributions of source code must retain the above copyright notice,
///    this list of conditions and the following disclaimer.
///
/// 2. Redistributions in binary form must reproduce the above copyright notice,
///    this list of conditions and the following disclaimer in the documentation
///    and/or other materials provided with the distribution.
///
/// 3. Neither the name of the copyright holder nor the names of its
///    contributors may be used to endorse or promote products derived from this
///    software without specific prior written permission.
///
/// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
/// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
/// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
/// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
/// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
/// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
/// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
/// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
/// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
/// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
/// POSSIBILITY OF SUCH DAMAGE.
///

#include <starneig_config.h>
#include <starneig/configuration.h>
#include "node_internal.h"
#ifdef STARNEIG_ENABLE_MPI
#include "../mpi/node_internal.h"
#include "../mpi/distr_matrix_internal.h"
#endif
#include "common.h"
#include "scratch.h"
#include <starneig/node.h>
#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <hwloc.h>
#include <starpu.h>
#ifdef MKL_SET_NUM_THREADS_LOCAL_FOUND
#include <mkl.h>
#endif
#if defined(OPENBLAS_SET_NUM_THREADS_FOUND) || \
defined(GOTO_SET_NUM_THREADS_FOUND)
#include <cblas.h>
#endif

static struct {
    /// initialization flag
    bool is_init;
    // initialization flags
    starneig_flag_t flags;
    /// library mode
    enum starneig_mode mode;
    /// blas mode
    enum starneig_blas_mode blas_mode;
    // original blas thread count
    int blas_threads_original;
    // StarPU worker bind mask
    unsigned starpu_workers_bindid[STARPU_NMAXWORKERS];
    // total number of available cpu cores
    int avail_cores;
    // total number of available gpus
    int avail_gpus;
    // number of used cpu cores
    int used_cores;
    // total number of used gpus
    int used_gpus;
} state = {
    .is_init = false,
    .flags = STARNEIG_DEFAULT,
    .mode = STARNEIG_MODE_OFF,
    .blas_mode = STARNEIG_BLAS_MODE_ORIGINAL,
    .blas_threads_original = -1,
    .avail_cores = 0,
    .avail_gpus = 0,
    .used_cores = 0,
    .used_gpus = 0
};

///
/// @brief Sets the number of BLAS threads.
///
/// @param[in] threads
///         Number of BLAS threads.
///
/// @return Previous BLAS thread count (can be 0).
///
static int set_blas_threads(int threads)
{
    starneig_verbose("Setting BLAS thread count to %d.", threads);
#ifdef MKL_SET_NUM_THREADS_LOCAL_FOUND
    return mkl_set_num_threads_local(threads);
#elif defined(OPENBLAS_SET_NUM_THREADS_FOUND)
    int old = openblas_get_num_threads();
    openblas_set_num_threads(threads);
    return old;
#elif defined(GOTO_SET_NUM_THREADS_FOUND)
    goto_set_num_threads(threads);
    return 1;
#else
    return 1;
#endif
}

///
/// @brief Sets the BLAS mode.
///
/// @param[in] mode
///         BLAS mode.
///
static void set_blas_mode(enum starneig_blas_mode mode)
{
    int old = -1;
    switch (mode) {
        case STARNEIG_BLAS_MODE_PARALLEL:
            starneig_verbose("Switching to parallel BLAS.");
            old = set_blas_threads(state.used_cores);
            break;
        case STARNEIG_BLAS_MODE_SEQUENTIAL:
            starneig_verbose("Switching to sequential BLAS.");
            old = set_blas_threads(1);
            break;
        default:
            starneig_verbose("Restoring BLAS mode.");
            if (0 <= state.blas_threads_original)
                set_blas_threads(state.blas_threads_original);
            state.blas_threads_original = -1;
    }
    if (state.blas_mode == STARNEIG_BLAS_MODE_ORIGINAL && 0 <= old)
        state.blas_threads_original = old;
    state.blas_mode = mode;
}

#ifdef MKL_SET_NUM_THREADS_LOCAL_FOUND

///
/// @brief Sets the per thread BLAS thread count to 1.
///
/// @param[in] arg
///         An unused argument.
///
static void set_worker_blas_mode(void *arg)
{
    mkl_set_num_threads_local(1);
}

#endif

///
/// @brief Reconfigures the node.
///
/// @param[in] cores
///         Number of CPU cores to use.
///
/// @param[in] gpus
///         Number of GPUs to use.
///
/// @param[in] mode
///         Library mode.
///
/// @param[in] blas_mode
///         BLAS mode.
///
#define CONFIGURE(cores, gpus, mode, blas_mode) \
    node_configure(cores, gpus, mode, blas_mode, __func__)

///
/// @brief Reconfigures the node.
///
/// @param[in] cores
///         Number of CPU cores to use.
///
/// @param[in] gpus
///         Number of GPUs to use.
///
/// @param[in] mode
///         Library mode.
///
/// @param[in] blas_mode
///         BLAS mode
///
/// @param[in] func
///         Name of the calling function.
///
static void node_configure(
    int cores, int gpus, enum starneig_mode mode,
    enum starneig_blas_mode blas_mode, char const *func)
{
#ifndef STARNEIG_ENABLE_MPI
    if (mode == STARNEIG_MODE_DM)
        starneig_fatal_error("StarPU was compiled without MPI support.");
#endif

    if (cores == state.used_cores && gpus == state.used_gpus &&
    mode == state.mode && blas_mode == state.blas_mode)
        return;

    starneig_verbose("Reconfiguring the library.");

    if (cores == state.used_cores && gpus == state.used_gpus &&
    mode == state.mode) {
        set_blas_mode(blas_mode);
        return;
    }

    //
    // shutdown StarPU
    //

    if (state.mode != STARNEIG_MODE_OFF) {
        starneig_node_resume_starpu();
#ifdef STARNEIG_ENABLE_CUDA
        if (0 < state.used_gpus) {
            starneig_verbose("Shutting down cuBLAS.");
            starpu_cublas_shutdown();
        }
#endif

        starneig_verbose("Shutting down StarPU.");

        starneig_scratch_unregister();
#ifdef STARNEIG_ENABLE_MPI
        starneig_mpi_cache_clear();
    if (state.mode == STARNEIG_MODE_DM &&
    state.flags & STARNEIG_AWAKE_MPI_WORKER)
        starneig_mpi_stop_persistent_starpumpi();
#endif

        starpu_task_wait_for_all();
        starpu_shutdown();
    }

    //
    // set the number of CPU cores
    //

    if (cores == 0)
        starneig_fatal_error("At least one CPU core must be selected.");

    if (cores < 0) {
        state.used_cores = state.avail_cores;
    }
    else {
        state.used_cores = MIN(cores, state.avail_cores);
        if (state.avail_cores < cores)
            starneig_warning(
                "Failed to acquire the desired number of CPU cores. "
                "Acquired %d.", state.used_cores);
    }

    //
    // set the number of GPUs
    //

    if (gpus < 0) {
        state.used_gpus = state.avail_gpus;
    }
    else {
#ifdef STARNEIG_ENABLE_CUDA
        state.used_gpus = MIN(gpus, state.avail_gpus);

        if (state.avail_gpus < gpus)
            starneig_warning(
                "Failed to acquire the desired number of CUDA devices. "
                "Acquired %d.", state.used_gpus);
#else
        if (0 < gpus)
            starneig_warning("StarPU was compiled without CUDA support.");
#endif
    }

    //
    // set BLAS threads
    //

    set_blas_mode(blas_mode);

    //
    // set mode
    //

    state.mode = mode;
    if (state.mode == STARNEIG_MODE_OFF)
        return;

    //
    // create StarPU configuration
    //

    starneig_verbose("Configuring StarPU.");

    struct starpu_conf conf;
    starpu_conf_init(&conf);

    int cpu_workers = state.used_cores;
    if (0 < state.used_gpus)
        cpu_workers -= state.used_gpus;
    if (state.mode == STARNEIG_MODE_DM)
        cpu_workers--;

    conf.ncpus = MAX(1, cpu_workers);
    conf.ncuda = state.used_gpus;
    conf.nopencl = 0;

//#if 1 < STARPU_MAJOR_VERSION || 2 < STARPU_MINOR_VERSION
    if (getenv("STARPU_WORKERS_CPUID") == NULL)
        conf.use_explicit_workers_bindid = 1;
//#endif
    memcpy(conf.workers_bindid, state.starpu_workers_bindid,
        sizeof(state.starpu_workers_bindid));

#ifdef STARNEIG_ENABLE_CUDA
    if (0 < state.used_gpus)
        conf.sched_policy_name = "dmdas";
    else
#endif
        conf.sched_policy_name = "prio";

    //
    // setup FXT
    //

    if (state.flags & STARNEIG_FXT_DISABLE) {
        starneig_verbose("Disabling FXT traces.");
        starpu_fxt_autostart_profiling(0);
    }
    else {
        char const *starpu_fxt_trace = getenv("STARPU_FXT_TRACE");
        if (starpu_fxt_trace == NULL || atoi(starpu_fxt_trace) != 0) {
            starneig_verbose("Keeping FXT traces enabled.");
            starpu_fxt_autostart_profiling(1);
        }
    }

    //
    // initialize StarPU
    //

    starneig_verbose("Starting StarPU.");

    unsigned seed = rand();
    int ret = starpu_init(&conf);
    srand(seed);

    if (ret != 0)
        starneig_fatal_error("Failed to initialize StarPU.");

    starpu_profiling_status_set(STARPU_PROFILING_ENABLE);
    starpu_malloc_set_align(64);

    //
    // initialize persistent StarPU-MPI
    //

#ifdef STARNEIG_ENABLE_MPI
    if (state.mode == STARNEIG_MODE_DM &&
    state.flags & STARNEIG_AWAKE_MPI_WORKER)
        starneig_mpi_start_persistent_starpumpi();
#endif

    //
    // cuBLAS
    //

    if (0 < state.used_gpus) {
        starneig_verbose("Initializing cuBLAS.");
        starpu_cublas_init();
    }

    //
    // configure workers
    //

#ifdef MKL_SET_NUM_THREADS_LOCAL_FOUND
    starneig_verbose(
        "MKL detected. Setting StarPU worker BLAS thread count to 1.");
    starpu_execute_on_each_worker(
        &set_worker_blas_mode, NULL, STARPU_CPU | STARPU_CUDA);
#endif

    starneig_node_pause_starpu();
}

void starneig_node_pause_starpu()
{
    if (state.flags & STARNEIG_AWAKE_WORKERS)
        return;

    starneig_verbose("Pausing StarPU workers.");
    starpu_pause();
}

void starneig_node_resume_starpu()
{
    if (state.flags & STARNEIG_AWAKE_WORKERS)
        return;

    starneig_verbose("Waking up StarPU workers.");
    starpu_resume();
}

void starneig_node_pause_awake_starpu()
{
    if (!(state.flags & STARNEIG_AWAKE_WORKERS))
        return;

    starneig_verbose("Pausing \"awake\" StarPU workers.");
    starpu_pause();
}

void starneig_node_resume_awake_starpu()
{
    if (!(state.flags & STARNEIG_AWAKE_WORKERS))
        return;

    starneig_verbose("Waking up \"awake\" StarPU workers.");
    starpu_resume();
}

__attribute__ ((visibility ("default")))
void starneig_node_init(int cores, int gpus, starneig_flag_t flags)
{
    starneig_set_message_mode(
        !(flags & STARNEIG_NO_MESSAGES), !(flags & STARNEIG_NO_VERBOSE));

    starneig_verbose("Initializing node.");

    if (state.is_init)
        starneig_fatal_error("The node is already initialized.");

    state.flags = flags;

    //
    // set up CUDA
    //

    state.avail_gpus = 0;
    state.used_gpus = 0;

#ifdef STARNEIG_ENABLE_CUDA
    // query the number of available CUDA devices
	if (cudaGetDeviceCount(&state.avail_gpus) != cudaSuccess) {
        starneig_warning("Failed to acquire CUDA device count.");
        state.avail_gpus = 0;
    }
#endif

    // query STARPU_NCUDA environment
    char *starpu_ncuda = getenv("STARPU_NCUDA");
    int num_starpu_gpus = (starpu_ncuda ? atoi(starpu_ncuda) : -1);

    if (num_starpu_gpus != -1) {
        if (state.avail_gpus < num_starpu_gpus)
        starneig_warning(
            "A conflict between STARPU_NCUDA and cudaGetDeviceCount().");
        state.avail_gpus = MIN(state.avail_gpus, num_starpu_gpus);
    }

    //
    // set up CPU cores
    //

    state.avail_cores = 0;
    state.used_cores = 0;

    // query the SLURM environment
    char *slurm = getenv("SLURM_CPUS_PER_TASK");
    const int num_slurm_cpus = (slurm ? atoi(slurm) : -1);

    // query STARPU_NCPUS environment
    char *starpu_ncpus = getenv("STARPU_NCPUS");
    int num_starpu_cpus = -1;
    if (starpu_ncpus) {
        if (0 < state.avail_gpus)
            num_starpu_cpus = atoi(starpu_ncpus) + 1;
        else
            num_starpu_cpus = atoi(starpu_ncpus);
    }

    // query hardware topology and CPU core binding mask

    hwloc_topology_t topology;
    hwloc_topology_init(&topology);
    hwloc_topology_load(topology);

    hwloc_cpuset_t res = hwloc_bitmap_alloc();

    hwloc_cpuset_t mask = hwloc_bitmap_alloc();
    hwloc_get_cpubind(topology, mask, HWLOC_CPUBIND_THREAD);

    int depth_cores = hwloc_get_type_depth(topology, HWLOC_OBJ_CORE);
    int num_cores = hwloc_get_nbobjs_by_depth(topology, depth_cores);

    int num_hwloc_cpus = 0;

    // iterate over all COREs
    for (int i = 0; i < num_cores && num_hwloc_cpus < STARPU_NMAXWORKERS; i++) {
        hwloc_obj_t core = hwloc_get_obj_by_depth(topology, depth_cores, i);

        // if the CORE has PUs inside it, ...
        if (core->first_child && core->first_child->type == HWLOC_OBJ_PU) {
            // iterate over them
            hwloc_obj_t pu = core->first_child;
            while (pu) {
                // if the PU is in the binding mask, ...
                hwloc_bitmap_and(res, mask, pu->cpuset);
                if (!hwloc_bitmap_iszero(res)) {
                    // add it to the worker list
                    state.starpu_workers_bindid[num_hwloc_cpus++] =
                        pu->logical_index;
                    break;
                }
                pu = pu->next_sibling;
            }
        }
        else {
            // if the CORE is in the binding mask, ...
            hwloc_bitmap_and(res, mask, core->cpuset);
            if (!hwloc_bitmap_iszero(res)) {
                // add it to the worker list
                state.starpu_workers_bindid[num_hwloc_cpus++] =
                    core->logical_index;
            }
        }
    }

    hwloc_bitmap_free(mask);
    hwloc_bitmap_free(res);
    hwloc_topology_destroy(topology);

    starneig_verbose_begin("Attached CPU cores");
    for (int i = 0; i < num_hwloc_cpus; i++)
        starneig_verbose_cont(" %d", state.starpu_workers_bindid[i]);
    starneig_verbose_cont(".\n");

    // set avail_cores

    state.avail_cores = num_hwloc_cpus;
    if (0 < num_starpu_cpus) {
        if (num_hwloc_cpus < num_starpu_cpus)
            starneig_warning(
                "A conflict between STARPU_NCPUS/STARPU_NCUDA and hwloc core "
                "binding mask.");
        state.avail_cores = MIN(state.avail_cores, num_starpu_cpus);
        if (0 < num_slurm_cpus) {
            if (num_slurm_cpus < num_starpu_cpus)
                starneig_warning(
                    "A conflict between STARPU_NCPUS/STARPU_NCUDA and "
                    "SLURM_CPUS_PER_TASK.");
            state.avail_cores = MIN(state.avail_cores, num_slurm_cpus);
        }
    }
    if (0 < num_slurm_cpus) {
        if (num_hwloc_cpus < num_slurm_cpus)
            starneig_warning(
                "A conflict between SLURM_CPUS_PER_TASK and hwloc core "
                "binding mask.");
        state.avail_cores = MIN(state.avail_cores, num_slurm_cpus);
    }

    if (state.avail_cores <= 0)
        starneig_fatal_error("Something unexpected happened.");

    state.is_init   = true;

    if (state.flags & STARNEIG_HINT_DM)
        CONFIGURE(cores, gpus, STARNEIG_MODE_DM, STARNEIG_BLAS_MODE_SEQUENTIAL);
    else
        CONFIGURE(cores, gpus, STARNEIG_MODE_SM, STARNEIG_BLAS_MODE_SEQUENTIAL);
}

__attribute__ ((visibility ("default")))
int starneig_node_initialized()
{
    return state.is_init;
}

__attribute__ ((visibility ("default")))
void starneig_node_finalize(void)
{
    CHECK_INIT();

    starneig_verbose("De-initializing node.");

    CONFIGURE(-1, -1, STARNEIG_MODE_OFF, STARNEIG_BLAS_MODE_ORIGINAL);

    starneig_set_message_mode(0, 0);

    state.avail_cores = 0;
    state.avail_gpus = 0;

    state.is_init = false;
}

__attribute__ ((visibility ("default")))
int starneig_node_get_cores(void)
{
    CHECK_INIT();
    return state.used_cores;
}

__attribute__ ((visibility ("default")))
void starneig_node_set_cores(int cores)
{
    CHECK_INIT();
    CONFIGURE(cores, starneig_node_get_gpus(), state.mode, state.blas_mode);
}

__attribute__ ((visibility ("default")))
int starneig_node_get_gpus(void)
{
    CHECK_INIT();
    return state.used_gpus;
}

__attribute__ ((visibility ("default")))
void starneig_node_set_gpus(int gpus)
{
    CHECK_INIT();
    CONFIGURE(starneig_node_get_cores(), gpus, state.mode, state.blas_mode);
}

__attribute__ ((visibility ("default")))
void starneig_node_set_mode(enum starneig_mode mode)
{
    CHECK_INIT();
    CONFIGURE(starneig_node_get_cores(), starneig_node_get_gpus(),
        mode, state.blas_mode);
}

void starneig_node_set_blas_mode(enum starneig_blas_mode blas_mode)
{
    CHECK_INIT();
    CONFIGURE(starneig_node_get_cores(), starneig_node_get_gpus(),
        state.mode, blas_mode);
}

ext: c

max_stars_repo_name: NLAFET/StarNEig

max_issues_repo_name: NLAFET/StarNEig

max_forks_repo_name: NLAFET/StarNEig


text:
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_cdf.h>
#include <stdio.h>

void evaluate(double x) {
  double probability = gsl_sf_erf(x);

  printf(
    "The probability that Normal(0, 1) random variable has a value "
    "between %g and %g is: %g\n", -x, x, probability
  );

  double area = 1 - 2 * gsl_cdf_gaussian_P(-x, 1);

  printf(
    "The integral of a Normal(0, 1) distribution "
    "between %g and %g is: %g\n", -x, x, area
  );
}

int main() {
  evaluate(0.25);
  evaluate(1.96);
}

max_stars_repo_name: tardate/LittleCodingKata

max_issues_repo_name: tardate/LittleCodingKata

max_forks_repo_name: tardate/LittleCodingKata


text:
/**
 * author: Jochen K"upper
 * created: Jan 2002
 * file: pygsl/src/statistics/longmodule.c
 * $Id: longmodule.c,v 1.10 2004/03/24 08:40:45 schnizer Exp $
 *
 * optional usage of Numeric module, available at http://numpy.sourceforge.net
 *"
 */

#include <Python.h>
#include <pygsl/error_helpers.h>
#include <pygsl/block_helpers.h>
#include <gsl/gsl_statistics.h>

/* include real functions for different data-types */

#define STATMOD_APPEND_PY_TYPE(X) X ## Int
#define STATMOD_APPEND_PYC_TYPE(X) X ## LONG
#define STATMOD_FUNC_EXT(X, Y) X ## _long ## Y
#define STATMOD_PY_AS_C PyInt_AsLong
#define STATMOD_C_TYPE long int
#include "functions.c"

/* initialization */
PyGSL_STATISTICS_INIT(long, "long")

/*
 * Local Variables:
 * mode: c
 * c-file-style: "Stroustrup"
 * End:
 */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* Copyright (c) 2016 PaddlePaddle Authors. All Rights Reserve.

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License. */

#ifndef MATHFUNCTIONS_H_
#define MATHFUNCTIONS_H_

#ifdef PADDLE_USE_MKL
#include <mkl.h>
#include <mkl_lapacke.h>
#else
extern "C" {
#include <cblas.h>
}
#ifdef PADDLE_USE_ATLAS
extern "C" {
#include <clapack.h>
}
#else
#include <lapacke.h>
#endif
#endif

#include <cmath>

namespace paddle {

template <class T>
void gemm(const CBLAS_TRANSPOSE transA,
          const CBLAS_TRANSPOSE transB,
          const int M,
          const int N,
          const int K,
          const T alpha,
          const T* A,
          const int lda,
          const T* B,
          const int ldb,
          const T beta,
          T* C,
          const int ldc);

template <class T>
int getrf(const CBLAS_ORDER Order,
          const int M,
          const int N,
          T* A,
          const int lda,
          int* ipiv);

template <class T>
int getri(
    const CBLAS_ORDER Order, const int N, T* A, const int lda, const int* ipiv);

template <class T>
void axpy(const int n, const T alpha, const T* x, T* y);

template <class T>
T dotProduct(const int n, const T* x, const T* y);

template <class T>
void vExp(const int n, const T* a, T* r);

template <class T>
void vPow(const int n, const T* a, const T b, T* r);

template <class T>
void vLog(const int n, const T* a, T* r);

template <class T>
void vAdd(const int n, const T* a, const T* b, T* r);

template <class T>
void vInvSqrt(const int n, const T* a, T* r);

template <class T>
void vLog1p(const int n, const T* a, T* r);

template <class T>
void vTanh(const int n, const T* a, T* r);

}  // namespace paddle

#endif  // MATHFUNCTIONS_H_

max_stars_repo_name: emailweixu/Paddle

max_issues_repo_name: anuranrc/Paddle

max_forks_repo_name: anuranrc/Paddle


text:
/* spdgemv.c
 *
 * Copyright (C) 2012-2014 Patrick Alken
 *
 * 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.
 */

#include <config.h>
#include <stdlib.h>
#include <math.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spmatrix.h>
#include <gsl/gsl_spblas.h>
#include <gsl/gsl_blas.h>

/*
gsl_spblas_dgemv()
  Multiply a sparse matrix and a vector

Inputs: alpha - scalar factor
        A     - sparse matrix
        x     - dense vector
        beta  - scalar factor
        y     - (input/output) dense vector

Return: y = alpha*op(A)*x + beta*y
*/

int
gsl_spblas_dgemv(const CBLAS_TRANSPOSE_t TransA, const double alpha,
                 const gsl_spmatrix *A, const gsl_vector *x,
                 const double beta, gsl_vector *y)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if ((TransA == CblasNoTrans && N != x->size) ||
      (TransA == CblasTrans && M != x->size))
    {
      GSL_ERROR("invalid length of x vector", GSL_EBADLEN);
    }
  else if ((TransA == CblasNoTrans && M != y->size) ||
           (TransA == CblasTrans && N != y->size))
    {
      GSL_ERROR("invalid length of y vector", GSL_EBADLEN);
    }
  else
    {
      size_t j, p;
      size_t incX, incY;
      size_t lenX, lenY;
      double *X, *Y;
      double *Ad;
      size_t *Ap, *Ai, *Aj;

      if (TransA == CblasNoTrans)
        {
          lenX = N;
          lenY = M;
        }
      else
        {
          lenX = M;
          lenY = N;
        }

      /* form y := beta*y */

      Y = y->data;
      incY = y->stride;

      if (beta == 0.0)
        {
          size_t jy = 0;
          for (j = 0; j < lenY; ++j)
            {
              Y[jy] = 0.0;
              jy += incY;
            }
        }
      else if (beta != 1.0)
        {
          size_t jy = 0;
          for (j = 0; j < lenY; ++j)
            {
              Y[jy] *= beta;
              jy += incY;
            }
        }

      if (alpha == 0.0)
        return GSL_SUCCESS;

      /* form y := alpha*op(A)*x + y */
      Ap = A->p;
      Ad = A->data;
      X = x->data;
      incX = x->stride;

      if ((GSL_SPMATRIX_ISCCS(A) && (TransA == CblasNoTrans)) ||
          (GSL_SPMATRIX_ISCRS(A) && (TransA == CblasTrans)))
        {
          Ai = A->i;

          for (j = 0; j < lenX; ++j)
            {
              for (p = Ap[j]; p < Ap[j + 1]; ++p)
                {
                  Y[Ai[p] * incY] += alpha * Ad[p] * X[j * incX];
                }
            }
        }
      else if ((GSL_SPMATRIX_ISCCS(A) && (TransA == CblasTrans)) ||
               (GSL_SPMATRIX_ISCRS(A) && (TransA == CblasNoTrans)))
        {
          Ai = A->i;

          for (j = 0; j < lenY; ++j)
            {
              for (p = Ap[j]; p < Ap[j + 1]; ++p)
                {
                  Y[j * incY] += alpha * Ad[p] * X[Ai[p] * incX];
                }
            }
        }
      else if (GSL_SPMATRIX_ISTRIPLET(A))
        {
          if (TransA == CblasNoTrans)
            {
              Ai = A->i;
              Aj = A->p;
            }
          else
            {
              Ai = A->p;
              Aj = A->i;
            }

          for (p = 0; p < A->nz; ++p)
            {
              Y[Ai[p] * incY] += alpha * Ad[p] * X[Aj[p] * incX];
            }
        }
      else
        {
          GSL_ERROR("unsupported matrix type", GSL_EINVAL);
        }

      return GSL_SUCCESS;
    }
} /* gsl_spblas_dgemv() */

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


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

#include <config.h>

#include <stddef.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <float.h>

#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_linalg.h>

#include "enorm.c"

/* Broyden's method. It is not an efficient or modern algorithm but
   gives an example of a rank-1 update.

   C.G. Broyden, "A Class of Methods for Solving Nonlinear
   Simultaneous Equations", Mathematics of Computation, vol 19 (1965),
   p 577-593

 */

typedef struct
  {
    gsl_matrix *H;
    gsl_matrix *lu;
    gsl_permutation *permutation;
    gsl_vector *v;
    gsl_vector *w;
    gsl_vector *y;
    gsl_vector *p;
    gsl_vector *fnew;
    gsl_vector *x_trial;
    double phi;
  }
broyden_state_t;

static int broyden_alloc (void *vstate, size_t n);
static int broyden_set (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static int broyden_iterate (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
static void broyden_free (void *vstate);

static int
broyden_alloc (void *vstate, size_t n)
{
  broyden_state_t *state = (broyden_state_t *) vstate;
  gsl_vector *v, *w, *y, *fnew, *x_trial, *p;
  gsl_permutation *perm;
  gsl_matrix *m, *H;

  m = gsl_matrix_calloc (n, n);

  if (m == 0)
    {
      GSL_ERROR ("failed to allocate space for lu", GSL_ENOMEM);
    }

  state->lu = m;

  perm = gsl_permutation_calloc (n);

  if (perm == 0)
    {
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for permutation", GSL_ENOMEM);
    }

  state->permutation = perm;

  H = gsl_matrix_calloc (n, n);

  if (H == 0)
    {
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for d", GSL_ENOMEM);
    }

  state->H = H;

  v = gsl_vector_calloc (n);

  if (v == 0)
    {
      gsl_matrix_free (H);
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for v", GSL_ENOMEM);
    }

  state->v = v;

  w = gsl_vector_calloc (n);

  if (w == 0)
    {
      gsl_vector_free (v);
      gsl_matrix_free (H);
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for w", GSL_ENOMEM);
    }

  state->w = w;

  y = gsl_vector_calloc (n);

  if (y == 0)
    {
      gsl_vector_free (w);
      gsl_vector_free (v);
      gsl_matrix_free (H);
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
    }

  state->y = y;

  fnew = gsl_vector_calloc (n);

  if (fnew == 0)
    {
      gsl_vector_free (y);
      gsl_vector_free (w);
      gsl_vector_free (v);
      gsl_matrix_free (H);
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for fnew", GSL_ENOMEM);
    }

  state->fnew = fnew;

  x_trial = gsl_vector_calloc (n);

  if (x_trial == 0)
    {
      gsl_vector_free (fnew);
      gsl_vector_free (y);
      gsl_vector_free (w);
      gsl_vector_free (v);
      gsl_matrix_free (H);
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for x_trial", GSL_ENOMEM);
    }

  state->x_trial = x_trial;

  p = gsl_vector_calloc (n);

  if (p == 0)
    {
      gsl_vector_free (x_trial);
      gsl_vector_free (fnew);
      gsl_vector_free (y);
      gsl_vector_free (w);
      gsl_vector_free (v);
      gsl_matrix_free (H);
      gsl_permutation_free (perm);
      gsl_matrix_free (m);

      GSL_ERROR ("failed to allocate space for p", GSL_ENOMEM);
    }

  state->p = p;

  return GSL_SUCCESS;
}

static int
broyden_set (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
{
  broyden_state_t *state = (broyden_state_t *) vstate;
  size_t i, j, n = function->n;
  int signum = 0;

  GSL_MULTIROOT_FN_EVAL (function, x, f);

  gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->lu);
  gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
  gsl_linalg_LU_invert (state->lu, state->permutation, state->H);

  for (i = 0; i < n; i++)
    for (j = 0; j < n; j++)
      gsl_matrix_set(state->H,i,j,-gsl_matrix_get(state->H,i,j));

  for (i = 0; i < n; i++)
    {
      gsl_vector_set (dx, i, 0.0);
    }

  state->phi = enorm (f);

  return GSL_SUCCESS;
}

static int
broyden_iterate (void *vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
{
  broyden_state_t *state = (broyden_state_t *) vstate;

  double phi0, phi1, t, lambda;

  gsl_matrix *H = state->H;
  gsl_vector *p = state->p;
  gsl_vector *y = state->y;
  gsl_vector *v = state->v;
  gsl_vector *w = state->w;
  gsl_vector *fnew = state->fnew;
  gsl_vector *x_trial = state->x_trial;
  gsl_matrix *lu = state->lu;
  gsl_permutation *perm = state->permutation;

  size_t i, j, iter;

  size_t n = function->n;

  /* p = H f */

  for (i = 0; i < n; i++)
    {
      double sum = 0;

      for (j = 0; j < n; j++)
        {
          sum += gsl_matrix_get (H, i, j) * gsl_vector_get (f, j);
        }
      gsl_vector_set (p, i, sum);
    }

  t = 1;
  iter = 0;

  phi0 = state->phi;

new_step:

  for (i = 0; i < n; i++)
    {
      double pi = gsl_vector_get (p, i);
      double xi = gsl_vector_get (x, i);
      gsl_vector_set (x_trial, i, xi + t * pi);
    }

  {
    int status = GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);

    if (status != GSL_SUCCESS)
      {
        return GSL_EBADFUNC;
      }
  }

  phi1 = enorm (fnew);

  iter++ ;

  if (phi1 > phi0 && iter < 10 && t > 0.1)
    {
      /* full step goes uphill, take a reduced step instead */

      double theta = phi1 / phi0;
      t *= (sqrt (1.0 + 6.0 * theta) - 1.0) / (3.0 * theta);
      goto new_step;
    }

  if (phi1 > phi0)
    {
      /* need to recompute Jacobian */
      int signum = 0;

      gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, lu);

      for (i = 0; i < n; i++)
        for (j = 0; j < n; j++)
          gsl_matrix_set(lu,i,j,-gsl_matrix_get(lu,i,j));

      gsl_linalg_LU_decomp (lu, perm, &signum);
      gsl_linalg_LU_invert (lu, perm, H);

      gsl_linalg_LU_solve (lu, perm, f, p);

      t = 1;

      for (i = 0; i < n; i++)
        {
          double pi = gsl_vector_get (p, i);
          double xi = gsl_vector_get (x, i);
          gsl_vector_set (x_trial, i, xi + t * pi);
        }

      {
        int status = GSL_MULTIROOT_FN_EVAL (function, x_trial, fnew);

        if (status != GSL_SUCCESS)
          {
            return GSL_EBADFUNC;
          }
      }

      phi1 = enorm (fnew);
    }

  /* y = f' - f */

  for (i = 0; i < n; i++)
    {
      double yi = gsl_vector_get (fnew, i) - gsl_vector_get (f, i);
      gsl_vector_set (y, i, yi);
    }

  /* v = H y */

  for (i = 0; i < n; i++)
    {
      double sum = 0;

      for (j = 0; j < n; j++)
        {
          sum += gsl_matrix_get (H, i, j) * gsl_vector_get (y, j);
        }

      gsl_vector_set (v, i, sum);
    }

  /* lambda = p . v */

  lambda = 0;

  for (i = 0; i < n; i++)
    {
      lambda += gsl_vector_get (p, i) * gsl_vector_get (v, i);
    }

  if (lambda == 0)
    {
      GSL_ERROR ("approximation to Jacobian has collapsed", GSL_EZERODIV) ;
    }

  /* v' = v + t * p */

  for (i = 0; i < n; i++)
    {
      double vi = gsl_vector_get (v, i) + t * gsl_vector_get (p, i);
      gsl_vector_set (v, i, vi);
    }

  /* w^T = p^T H */

  for (i = 0; i < n; i++)
    {
      double sum = 0;

      for (j = 0; j < n; j++)
        {
          sum += gsl_matrix_get (H, j, i) * gsl_vector_get (p, j);
        }

      gsl_vector_set (w, i, sum);
    }

  /* Hij -> Hij - (vi wj / lambda) */

  for (i = 0; i < n; i++)
    {
      double vi = gsl_vector_get (v, i);

      for (j = 0; j < n; j++)
        {
          double wj = gsl_vector_get (w, j);
          double Hij = gsl_matrix_get (H, i, j) - vi * wj / lambda;
          gsl_matrix_set (H, i, j, Hij);
        }
    }

  /* copy fnew into f */

  gsl_vector_memcpy (f, fnew);

  /* copy x_trial into x */

  gsl_vector_memcpy (x, x_trial);

  for (i = 0; i < n; i++)
    {
      double pi = gsl_vector_get (p, i);
      gsl_vector_set (dx, i, t * pi);
    }

  state->phi = phi1;

  return GSL_SUCCESS;
}

static void
broyden_free (void *vstate)
{
  broyden_state_t *state = (broyden_state_t *) vstate;

  gsl_matrix_free (state->H);

  gsl_matrix_free (state->lu);
  gsl_permutation_free (state->permutation);

  gsl_vector_free (state->v);
  gsl_vector_free (state->w);
  gsl_vector_free (state->y);
  gsl_vector_free (state->p);

  gsl_vector_free (state->fnew);
  gsl_vector_free (state->x_trial);

}

static const gsl_multiroot_fsolver_type broyden_type =
{"broyden",                     /* name */
 sizeof (broyden_state_t),
 &broyden_alloc,
 &broyden_set,
 &broyden_iterate,
 &broyden_free};

const gsl_multiroot_fsolver_type *gsl_multiroot_fsolver_broyden = &broyden_type;

max_stars_repo_licenses: Net-SNMP Xnet

max_issues_repo_licenses: Net-SNMP Xnet

max_forks_repo_licenses: Net-SNMP Xnet