text:
#include <math.h>
#include <stdlib.h>
#if !defined(__APPLE__)
#include <malloc.h>
#endif
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <string.h>

#include <fftw3.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_sf_expint.h>
#include <gsl/gsl_deriv.h>

#include "../cosmolike_core/theory/basics.c"
#include "../cosmolike_core/theory/structs.c"
#include "../cosmolike_core/theory/parameters.c"
#include "../cosmolike_core/emu17/P_cb/emu.c"
#include "../cosmolike_core/theory/recompute.c"
#include "../cosmolike_core/theory/cosmo3D.c"
#include "../cosmolike_core/theory/redshift.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/HOD.c"
#include "../cosmolike_core/theory/cosmo2D_fourier.c"
#include "../cosmolike_core/theory/IA.c"
#include "../cosmolike_core/theory/cluster.c"
#include "../cosmolike_core/theory/BAO.c"
#include "../cosmolike_core/theory/external_prior.c"
#include "../cosmolike_core/theory/covariances_3D.c"
#include "../cosmolike_core/theory/covariances_fourier.c"
#include "../cosmolike_core/theory/covariances_cluster.c"
#include "init_SRD.c"

double C_shear_tomo_sys(double ell,int z1,int z2);
double C_cgl_tomo_sys(double ell_Cluster,int zl,int nN, int zs);
double C_gl_tomo_sys(double ell,int zl,int zs);
void set_data_shear(int Ncl, double *ell, double *data, int start);
void set_data_ggl(int Ncl, double *ell, double *data, int start);
void set_data_clustering(int Ncl, double *ell, double *data, int start);
void set_data_cluster_N(double *data, int start);
void set_data_cgl(double *ell_Cluster, double *data, int start);
void compute_data_vector(char *details, double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q, double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope);
double log_multi_like(double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q, double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope);
double write_vector_wrapper(char *details, input_cosmo_params ic, input_nuisance_params in);
double log_like_wrapper(input_cosmo_params ic, input_nuisance_params in);
int get_N_tomo_shear(void);
int get_N_tomo_clustering(void);
int get_N_ggl(void);
int get_N_ell(void);

int get_N_tomo_shear(void){
  return tomo.shear_Nbin;
}
int get_N_tomo_clustering(void){
  return tomo.clustering_Nbin;
}
int get_N_ggl(void){
  return tomo.ggl_Npowerspectra;
}
int get_N_ell(void){
  return like.Ncl;
}

double C_shear_tomo_sys(double ell, int z1, int z2)
{
  double C;
  // C= C_shear_tomo_nointerp(ell,z1,z2);
  // if(like.IA==1) C+=C_II_nointerp(ell,z1,z2)+C_GI_nointerp(ell,z1,z2);

  if(like.IA!=1) C= C_shear_tomo_nointerp(ell,z1,z2);
  if(like.IA==1) C= C_shear_shear_IA(ell,z1,z2);
  //if(like.IA==1) C = C_shear_tomo_nointerp(ell,z1,z2)+C_II_nointerp(ell,z1,z2)+C_GI_nointerp(ell,z1,z2);
  if(like.IA==2) C += C_II_lin_nointerp(ell,z1,z2)+C_GI_lin_nointerp(ell,z1,z2);
  if(like.shearcalib==1) C *=(1.0+nuisance.shear_calibration_m[z1])*(1.0+nuisance.shear_calibration_m[z2]);
  //printf("%le %d %d %le\n",ell,z1,z2,C_shear_tomo_nointerp(ell,z1,z2)+C_II_JB_nointerp(ell,z1,z2)+C_GI_JB_nointerp(ell,z1,z2));
return C;
}

double C_gl_tomo_sys(double ell,int zl,int zs)
{
  double C;
  // C=C_gl_tomo_nointerp(ell,zl,zs);
  // if(like.IA==1) C += C_gI_nointerp(ell,zl,zs);

  if(like.IA!=1) C=C_gl_tomo_nointerp(ell,zl,zs);
  if(like.IA==1) C = C_ggl_IA(ell,zl,zs);
  if(like.IA==2) C += C_gI_lin_nointerp(ell,zl,zs);
  if(like.shearcalib==1) C *=(1.0+nuisance.shear_calibration_m[zs]);
return C;
}

double C_cgl_tomo_sys(double ell_Cluster, int zl,int nN, int zs)
{
  double C;
  C=C_cgl_tomo_nointerp(ell_Cluster,zl,nN,zs);
  //if(like.IA!=0) C +=
  if(like.shearcalib==1) C *=(1.0+nuisance.shear_calibration_m[zs]);
return C;
}

void set_data_shear(int Ncl, double *ell, double *data, int start)
{
  int i,z1,z2,nz;
  double a;
  for (nz = 0; nz < tomo.shear_Npowerspectra; nz++){
    z1 = Z1(nz); z2 = Z2(nz);
    for (i = 0; i < Ncl; i++){
      if (ell[i]<like.lmax_shear){ data[Ncl*nz+i] = C_shear_tomo_sys(ell[i],z1,z2);}
      else {data[Ncl*nz+i] = 0.;}
    }
  }
}

void set_data_ggl(int Ncl, double *ell, double *data, int start)
{
  int i, zl,zs,nz;
  for (nz = 0; nz < tomo.ggl_Npowerspectra; nz++){
    zl = ZL(nz); zs = ZS(nz);
    for (i = 0; i < Ncl; i++){
      if (test_kmax(ell[i],zl)){
        data[start+(Ncl*nz)+i] = C_gl_tomo_sys(ell[i],zl,zs);
      }
      else{
        data[start+(Ncl*nz)+i] = 0.;
      }
    }
  }
}

void set_data_clustering(int Ncl, double *ell, double *data, int start){
  int i, nz;
  for (nz = 0; nz < tomo.clustering_Npowerspectra; nz++){
    //printf("%d %e %e\n",nz, gbias.b[nz][1],pf_photoz(gbias.b[nz][1],nz));
    for (i = 0; i < Ncl; i++){
      if (test_kmax(ell[i],nz)){data[start+(Ncl*nz)+i] = C_cl_tomo_nointerp(ell[i],nz,nz);}
      else{data[start+(Ncl*nz)+i] = 0.;}
      //printf("%d %d %le %le\n",nz,nz,ell[i],data[Ncl*(tomo.shear_Npowerspectra+tomo.ggl_Npowerspectra + nz)+i]);
    }
  }
}

void set_data_cluster_N(double *data, int start){
  int nN, nz;
  for (nz = 0; nz < tomo.cluster_Nbin; nz++){
    for (nN = 0; nN < Cluster.N200_Nbin; nN++){
      data[start+Cluster.N200_Nbin*nz+nN] = N_N200(nz, nN);
    }
  }
}

void set_data_cgl(double *ell_Cluster, double *data, int start)
{
  int zl,zs,nN,nz,i,j;
  for(nN = 0; nN < Cluster.N200_Nbin; nN++){
    for (nz = 0; nz < tomo.cgl_Npowerspectra; nz++){
      zl = ZC(nz); zs = ZSC(nz);
      for (i = 0; i < Cluster.lbin; i++){
        j = start;
        j += (nz*Cluster.N200_Nbin+nN)*Cluster.lbin +i;
        data[j] = C_cgl_tomo_sys(ell_Cluster[i],zl,nN,zs);
      }
    }
  }
}

int set_cosmology_params(double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu)
{
  cosmology.Omega_m=OMM;
  cosmology.Omega_v= 1.0-cosmology.Omega_m;
  cosmology.sigma_8=S8;
  cosmology.n_spec= NS;
  cosmology.w0=W0;
  cosmology.wa=WA;
  cosmology.omb=OMB;
  cosmology.h0=H0;
  cosmology.MGSigma=MGSigma;
  cosmology.MGmu=MGmu;

  if (cosmology.Omega_m < 0.05 || cosmology.Omega_m > 0.6) return 0;
  if (cosmology.omb < 0.04 || cosmology.omb > 0.055) return 0;
  if (cosmology.sigma_8 < 0.5 || cosmology.sigma_8 > 1.1) return 0;
  if (cosmology.n_spec < 0.84 || cosmology.n_spec > 1.06) return 0;
  if (cosmology.w0 < -2.1 || cosmology.w0 > -0.0) return 0;
  if (cosmology.wa < -2.6 || cosmology.wa > 2.6) return 0;
  if (cosmology.h0 < 0.4 || cosmology.h0 > 0.9) return 0;
  return 1;
}

void set_nuisance_shear_calib(double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10)
{
  nuisance.shear_calibration_m[0] = M1;
  nuisance.shear_calibration_m[1] = M2;
  nuisance.shear_calibration_m[2] = M3;
  nuisance.shear_calibration_m[3] = M4;
  nuisance.shear_calibration_m[4] = M5;
  nuisance.shear_calibration_m[5] = M6;
  nuisance.shear_calibration_m[6] = M7;
  nuisance.shear_calibration_m[7] = M8;
  nuisance.shear_calibration_m[8] = M9;
  nuisance.shear_calibration_m[9] = M10;
}

int set_nuisance_shear_photoz(double SP1,double SP2,double SP3,double SP4,double SP5,double SP6,double SP7,double SP8,double SP9,double SP10,double SPS1)
{
  int i;
  nuisance.bias_zphot_shear[0]=SP1;
  nuisance.bias_zphot_shear[1]=SP2;
  nuisance.bias_zphot_shear[2]=SP3;
  nuisance.bias_zphot_shear[3]=SP4;
  nuisance.bias_zphot_shear[4]=SP5;
  nuisance.bias_zphot_shear[5]=SP6;
  nuisance.bias_zphot_shear[6]=SP7;
  nuisance.bias_zphot_shear[7]=SP8;
  nuisance.bias_zphot_shear[8]=SP9;
  nuisance.bias_zphot_shear[9]=SP10;

  for (i=0;i<tomo.shear_Nbin; i++){
    nuisance.sigma_zphot_shear[i]=SPS1;
    if (nuisance.sigma_zphot_shear[i]<0.001) return 0;
  }
  return 1;
}

int set_nuisance_clustering_photoz(double CP1,double CP2,double CP3,double CP4,double CP5,double CP6,double CP7,double CP8,double CP9,double CP10,double CPS1)
{
  int i;
  nuisance.bias_zphot_clustering[0]=CP1;
  nuisance.bias_zphot_clustering[1]=CP2;
  nuisance.bias_zphot_clustering[2]=CP3;
  nuisance.bias_zphot_clustering[3]=CP4;
  nuisance.bias_zphot_clustering[4]=CP5;
  nuisance.bias_zphot_clustering[5]=CP6;
  nuisance.bias_zphot_clustering[6]=CP7;
  nuisance.bias_zphot_clustering[7]=CP8;
  nuisance.bias_zphot_clustering[8]=CP9;
  nuisance.bias_zphot_clustering[9]=CP10;

  for (i=0;i<tomo.clustering_Nbin; i++){
    nuisance.sigma_zphot_clustering[i]=CPS1;
    if (nuisance.sigma_zphot_clustering[i]<0.001) return 0;
  }
  return 1;
}

int set_nuisance_ia(double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q)
{

  nuisance.A_ia=A_ia;
  nuisance.beta_ia=beta_ia;
  nuisance.eta_ia=eta_ia;
  nuisance.eta_ia_highz=eta_ia_highz;
  nuisance.LF_alpha=LF_alpha;
  nuisance.LF_P=LF_P;
  nuisance.LF_Q=LF_Q;
  nuisance.LF_red_alpha=LF_red_alpha;
  nuisance.LF_red_P=LF_red_P;
  nuisance.LF_red_Q=LF_red_Q;

  if (nuisance.A_ia < 0.0 || nuisance.A_ia > 10.0) return 0;
  if (nuisance.beta_ia < -1.0 || nuisance.beta_ia > 3.0) return 0;
  if (nuisance.eta_ia < -3.0 || nuisance.eta_ia> 3.0) return 0;
  if (nuisance.eta_ia_highz < -1.0 || nuisance.eta_ia_highz> 1.0) return 0;
  // printf("%le %le %le %le %le %le %le %le %le %le\n",nuisance.A_ia,nuisance.beta_ia,nuisance.eta_ia,nuisance.eta_ia_highz,nuisance.LF_alpha,nuisance.LF_P, nuisance.LF_Q,nuisance.LF_red_alpha,nuisance.LF_red_P,nuisance.LF_red_Q);
  // if(like.IA!=0){
  //  if (check_LF()) return 0;
  // }
return 1;
}

int set_nuisance_cluster_Mobs(double cluster_Mobs_lgN0,  double cluster_Mobs_alpha, double cluster_Mobs_beta, double cluster_Mobs_sigma0, double cluster_Mobs_sigma_qm, double cluster_Mobs_sigma_qz)
{
  //  nuisance.cluster_Mobs_lgM0 = mass_obs_norm;  //fiducial : 1.72+log(1.e+14*0.7); could use e.g. sigma = 0.2 Gaussian prior
  //  nuisance.cluster_Mobs_alpha = mass_obs_slope; //fiducial: 1.08; e.g. sigma = 0.1 Gaussian prior
  //  nuisance.cluster_Mobs_beta = mass_z_slope; //fiducial: 0.0; e.g. sigma = 0.1 Gaussian prior
  //  nuisance.cluster_Mobs_sigma = mass_obs_scatter; //fiducial 0.25; e.g. sigma = 0.05 Gaussian prior

  // fiducial values and priors from Murata et al. (2018) except for redshift-related parameters
  nuisance.cluster_Mobs_lgN0 = cluster_Mobs_lgN0; //fiducial: 3.207, flat prior [0.5, 5.0]
  nuisance.cluster_Mobs_alpha = cluster_Mobs_alpha; //fiducial: 0.993, flat prior [0.0, 2.0]
  nuisance.cluster_Mobs_beta = cluster_Mobs_beta; //fiducial: 0.0, flat prior [-1.5, 1.5]
  nuisance.cluster_Mobs_sigma0 = cluster_Mobs_sigma0; //fiducial: 0.456, flat prior [0.0, 1.5]
  nuisance.cluster_Mobs_sigma_qm = cluster_Mobs_sigma_qm; //fiducial: -0.169, flat prior [-1.5, 1.5]
  nuisance.cluster_Mobs_sigma_qz = cluster_Mobs_sigma_qz; //fiducial: 0.0, flat prior [-1.5, 1.5]

  if (nuisance.cluster_Mobs_lgN0 < 0.5 || nuisance.cluster_Mobs_lgN0 > 5.0) return 0;
  if (nuisance.cluster_Mobs_alpha < 0.0 || nuisance.cluster_Mobs_alpha > 2.0) return 0;
  if (nuisance.cluster_Mobs_beta < -1.5 || nuisance.cluster_Mobs_beta > 1.5) return 0;
  if (nuisance.cluster_Mobs_sigma0 < 0.0|| nuisance.cluster_Mobs_sigma0 > 1.5) return 0;
  if (nuisance.cluster_Mobs_sigma_qm < -1.5 && nuisance.cluster_Mobs_sigma_qm > 1.5) return 0;
  if (nuisance.cluster_Mobs_sigma_qz < -1.5 && nuisance.cluster_Mobs_sigma_qz > 1.5)return 0;

return 1;
}

int set_nuisance_gbias(double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8,double B9, double B10)
{
  int i;
  gbias.b[0] = B1;
  gbias.b[1] = B2;
  gbias.b[2] = B3;
  gbias.b[3] = B4;
  gbias.b[4] = B5;
  gbias.b[5] = B6;
  gbias.b[6] = B7;
  gbias.b[7] = B8;
  gbias.b[8] = B9;
  gbias.b[9] = B10;
  if(like.bias==1){
    for (i = 0; i < 10; i++){
      if (gbias.b[i] < 0.8 || gbias.b[i] > 3.0) return 0;
    }
  }
  return 1;
}

double log_L_SRD_SN_Y1_RENEE() // computed from /Users/teifler/Dropbox/cosmolike_store/LSSTawakens/ReneeSNChains/lsst_y1_jan29_nostarts_varyM_oldomb_rand.txt
{
  double log_L = 0.;
  double param_diff[4];
  //double inv[4][4]={{3.79965832e+03,-1.67422234e+01, 1.49107053e+03, 2.03605752e+02},{ -1.67422234e+01, 4.73633438e+02,-2.97547647e+00,-2.21008677e-01},{1.49107053e+03,-2.97547647e+00, 7.06767724e+02, 7.78147349e+01},{2.03605752e+02,-2.21008677e-01, 7.78147349e+01, 1.20866157e+01}};

  // double inv[4][4]={{2.54176679e+03,2.52990975e+01,8.90935426e+02,1.36402504e+02},{2.52990975e+01,4.81470268e+02,1.01904094e+01,6.95186662e-01},{8.90935426e+02,1.01904094e+01,4.02864434e+02,4.60275519e+01},{1.36402504e+02,6.95186662e-01,4.60275519e+01,8.41721305e+00}};

  double inv[4][4]={{2.92724944e+03,-2.18583441e+01, 1.31758412e+03, 1.57579361e+02},
 {-2.18583441e+01, 4.63427974e+02,-1.59603150e+01,-8.00922379e-01},
 {1.31758412e+03,-1.59603150e+01, 7.22127857e+02, 6.67792329e+01},
 {1.57579361e+02,-8.00922379e-01, 6.67792329e+01, 9.71021170e+00}};
  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.h0-0.6727;
  param_diff[2] = cosmology.w0+1.0;
  param_diff[3] = cosmology.wa;

  log_L = -0.5*do_matrix_mult_invcov(4,inv,param_diff);
  return log_L;
 }

double log_L_SRD_SN_Y10_RENEE() // computed from /Users/teifler/Dropbox/cosmolike_store/LSSTawakens/ReneeSNChains/lsst_y10_jan29_nostarts_varyM_oldomb_rand.txt
{
  double log_L = 0.;
  double param_diff[4];

double inv[4][4]={{1.13549508e+04,-7.56474778e+01,7.30717358e+03,5.90184798e+02},
 { -7.56474778e+01,4.91930625e+02,-4.39988764e+01,-7.28574193e+00},
 {7.30717358e+03,-4.39988764e+01,5.38324084e+03,3.24899498e+02},
 {5.90184798e+02,-7.28574193e+00,3.24899498e+02,3.84274454e+01}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.h0-0.6727;
  param_diff[2] = cosmology.w0+1.0;
  param_diff[3] = cosmology.wa;

  log_L = -0.5*do_matrix_mult_invcov(4,inv, param_diff);
  return log_L;
 }

double log_L_SRD_SL_Y1_TOM() // computed from /Users/teifler/Dropbox/cosmolike_store/LSSTawakens/TomSLChains/SL_LSSTY1.txt
{
  double log_L = 0.;
  double param_diff[4];
  double inv[4][4]={{9.57953440e+01,1.90114212e+01,1.13561665e+00,2.33866782e+02},{1.90114212e+01,1.31631828e+01,-1.88534652e-02,1.45671419e+02},{1.13561665e+00,-1.88534652e-02,6.29296427e-01,-6.10433226e+00},{2.33866782e+02,1.45671419e+02,-6.10433226e+00,2.71769694e+03}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.w0+1.0;
  param_diff[2] = cosmology.wa;
  param_diff[3] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(4,inv, param_diff);
  return log_L;
 }

double log_L_SRD_SL_Y10_TOM() // computed from /Users/teifler/Dropbox/cosmolike_store/LSSTawakens/TomSLChains/SL_LSSTY10.txt
{
  double log_L = 0.;
  double param_diff[4];
  double inv[4][4]={{2.44407663e+02,1.31525702e+02,1.62494864e+01,1.03889766e+03},{1.31525702e+02,1.90748959e+02,8.03881315e+00,1.99392208e+03},{1.62494864e+01,8.03881315e+00,2.23917868e+00,-6.64709076e+00},{1.03889766e+03,1.99392208e+03,-6.64709076e+00,3.31389425e+04}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.w0+1.0;
  param_diff[2] = cosmology.wa;
  param_diff[3] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(4,inv, param_diff);
  return log_L;
 }

double log_L_SRD_LSST_Y1_clusters() // // including MOR for Y1, see fisher.py
{
  double log_L = 0.;
  double param_diff[7];
  double inv[7][7]={{3.858101e+04, 3.588867e+04, 2.231543e+03, -9.737292e+02, -2.076847e+02, -1.567836e+03, 5.165868e+02},{3.588867e+04, 4.046241e+04, 2.204910e+03, -4.306084e+02, -1.921012e+02, -1.153062e+03, 4.146891e+02},{2.231543e+03, 2.204910e+03, 4.640472e+02, -3.849041e+01, -1.164098e+01, -2.997845e+02, 8.606974e+01},{-9.737292e+02, -4.306084e+02, -3.849041e+01, 1.057829e+02, 1.234460e+01, 3.159550e+01, -1.021792e+01},{-2.076847e+02, -1.921012e+02, -1.164098e+01, 1.234460e+01, 3.965104e+00, 6.461167e+00, -2.324549e+00},{-1.567836e+03, -1.153062e+03, -2.997845e+02, 3.159550e+01, 6.461167e+00, 9.845550e+04, -2.104350e+02},{5.165868e+02, 4.146891e+02, 8.606974e+01, -1.021792e+01, -2.324549e+00, -2.104350e+02, 2.121782e+02}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.sigma_8-0.831;
  param_diff[2] = cosmology.n_spec-0.9645;
  param_diff[3] = cosmology.w0+1.0;
  param_diff[4] = cosmology.wa;
  param_diff[5] = cosmology.omb-0.0491685;
  param_diff[6] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(7,inv, param_diff);
  return log_L;
 }

double log_L_SRD_LSST_Y10_clusters() // including MOR for Y10, see fisher.py
{
  double log_L = 0.;
  double param_diff[7];
  double inv[7][7]={{8.011883e+04, 7.501415e+04, 4.972257e+03, -1.009514e+03, -2.332928e+02, -4.035026e+03, 1.280726e+03},{7.501415e+04, 9.212235e+04, 5.029562e+03, 1.256229e+03, -8.904458e+00, -3.884022e+03, 1.236633e+03},{4.972257e+03, 5.029562e+03, 8.996244e+02, -1.335280e+01, -1.271912e+01, -9.151661e+02, 2.624954e+02},{-1.009514e+03, 1.256229e+03, -1.335280e+01, 3.943470e+02, 5.659304e+01, -3.776163e+01, 7.168167e+00},{-2.332928e+02, -8.904458e+00, -1.271912e+01, 5.659304e+01, 1.167193e+01, 4.161086e+00, -1.872393e+00},{-4.035026e+03, -3.884022e+03, -9.151661e+02, -3.776163e+01, 4.161086e+00, 1.000908e+05, -6.424232e+02},{1.280726e+03, 1.236633e+03, 2.624954e+02, 7.168167e+00, -1.872393e+00, -6.424232e+02, 3.273443e+02}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.sigma_8-0.831;
  param_diff[2] = cosmology.n_spec-0.9645;
  param_diff[3] = cosmology.w0+1.0;
  param_diff[4] = cosmology.wa;
  param_diff[5] = cosmology.omb-0.0491685;
  param_diff[6] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(7,inv, param_diff);
  return log_L;
 }

double log_L_SRD_stage3_exceptw0wa()
{
  double log_L = 0.;
  double param_diff[7];
  double inv[7][7]={{8.52388933e+03,-2.38254882e-03,4.32405053e+00,0.00000000e+00,0.00000000e+00,1.22554260e+02,5.48470403e+01},{-2.38254882e-03,3.34465215e+01,3.86664237e-01,0.00000000e+00,0.00000000e+00,-8.71498515e+01,-7.09315084e+00},{4.32405053e+00,3.86664237e-01,2.47843770e+02,0.00000000e+00,0.00000000e+00,-1.21079947e+02,4.07696833e+00},{0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00},{0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00,0.00000000e+00},{1.22554260e+02,-8.71498515e+01,-1.21079947e+02,-0.00000000e+00,-0.00000000e+00,7.17126110e+04,3.07891510e+03},{5.48470403e+01,-7.09315084e+00,4.07696833e+00,0.00000000e+00,0.00000000e+00,3.07891510e+03,5.31816639e+02}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.sigma_8-0.831;
  param_diff[2] = cosmology.n_spec-0.9645;
  param_diff[3] = cosmology.w0+1.0;
  param_diff[4] = cosmology.wa;
  param_diff[5] = cosmology.omb-0.0491685;
  param_diff[6] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(7,inv, param_diff);
  return log_L;
 }

double log_L_SRD_stage3()
{
  double log_L = 0.;
  double param_diff[7];
  double inv[7][7]={{8.52889541e+03,9.34266060e-03,2.21513407e+00,1.85008170e+01,6.12185662e+00,1.17720010e+02,5.30670581e+01},{9.34266060e-03,3.34475078e+01,3.79135616e-01,-8.02192934e-02,5.88736315e-02,-8.69202865e+01,-7.08029577e+00},{2.21513407e+00,3.79135616e-01,2.48739187e+02,-7.46022981e+00,-2.69922217e+00,-1.19693897e+02,4.78085262e+00},{1.85008170e+01,-8.02192934e-02,-7.46022981e+00,8.42937127e+01,1.68856256e+01,-4.89063661e+01,-8.77194357e+00},{6.12185662e+00,5.88736315e-02,-2.69922217e+00,1.68856256e+01,9.55489400e+00,5.27704214e+00,-1.38597499e+00},{1.17720010e+02,-8.69202865e+01,-1.19693897e+02,-4.89063661e+01,5.27704214e+00,7.17777988e+04,3.08491105e+03},{5.30670581e+01,-7.08029577e+00,4.78085262e+00,-8.77194357e+00,-1.38597499e+00,3.08491105e+03,5.32751808e+02}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.sigma_8-0.831;
  param_diff[2] = cosmology.n_spec-0.9645;
  param_diff[3] = cosmology.w0+1.0;
  param_diff[4] = cosmology.wa;
  param_diff[5] = cosmology.omb-0.0491685;
  param_diff[6] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(7,inv, param_diff);
  return log_L;
 }

double log_L_SRD_LSST_Y1_3x2pt()
{
  double log_L = 0.;
  double param_diff[7];
  double inv[7][7]={{6.461019e+05, 4.105098e+05, 9.063469e+04, -2.968463e+04, -6.656297e+03, -3.266829e+05, 7.947019e+04},{4.105098e+05, 3.102396e+05, 4.061214e+04, -2.051833e+04, -5.069455e+03, -1.156574e+05, 2.911602e+04},{9.063469e+04, 4.061214e+04, 3.652491e+04, -1.033492e+03, -1.571074e+02, -9.436089e+04, 2.396154e+04},{-2.968463e+04, -2.051833e+04, -1.033492e+03, 2.081364e+03, 4.721662e+02, 7.414854e+03, -1.864708e+03},{-6.656297e+03, -5.069455e+03, -1.571074e+02, 4.721662e+02, 1.200287e+02, 1.340576e+03, -3.293614e+02},{-3.266829e+05, -1.156574e+05, -9.436089e+04, 7.414854e+03, 1.340576e+03, 5.388793e+05, -9.884641e+04},{7.947019e+04, 2.911602e+04, 2.396154e+04, -1.864708e+03, -3.293614e+02, -9.884641e+04, 2.373400e+04}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.sigma_8-0.831;
  param_diff[2] = cosmology.n_spec-0.9645;
  param_diff[3] = cosmology.w0+1.0;
  param_diff[4] = cosmology.wa;
  param_diff[5] = cosmology.omb-0.0491685;
  param_diff[6] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(7,inv, param_diff);
  return log_L;
 }

double log_L_SRD_LSST_Y10_3x2pt()
{
  double log_L = 0.;
  double param_diff[7];
  double inv[7][7]={{1.732125e+06, 1.363116e+06, 2.420885e+05, -6.384395e+04, -1.991510e+04, -7.567872e+05, 1.858361e+05},{1.363116e+06, 1.196537e+06, 1.469326e+05, -5.148233e+04, -1.796618e+04, -4.526364e+05, 1.126267e+05},{2.420885e+05, 1.469326e+05, 1.011204e+05, -1.370858e+03, -4.391526e+02, -2.001808e+05, 5.228884e+04},{-6.384395e+04, -5.148233e+04, -1.370858e+03, 4.153843e+03, 1.137490e+03, 1.286287e+04, -3.253747e+03},{-1.991510e+04, -1.796618e+04, -4.391526e+02, 1.137490e+03, 3.720224e+02, 4.198170e+03, -1.020253e+03},{-7.567872e+05, -4.526364e+05, -2.001808e+05, 1.286287e+04, 4.198170e+03, 8.507285e+05, -1.685201e+05},{1.858361e+05, 1.126267e+05, 5.228884e+04, -3.253747e+03, -1.020253e+03, -1.685201e+05, 4.064231e+04}};

  param_diff[0] = cosmology.Omega_m-0.3156;
  param_diff[1] = cosmology.sigma_8-0.831;
  param_diff[2] = cosmology.n_spec-0.9645;
  param_diff[3] = cosmology.w0+1.0;
  param_diff[4] = cosmology.wa;
  param_diff[5] = cosmology.omb-0.0491685;
  param_diff[6] = cosmology.h0-0.6727;

  log_L = -0.5*do_matrix_mult_invcov(7,inv, param_diff);
  return log_L;
 }

double log_L_DESI()
{
  double log_L = 0.;
  double param_diff[2];
  double inv[2][2]={{1411.40888852,449.46652105},{449.46652105,   291.94320759}};

  param_diff[0] = cosmology.w0+1.0;
  param_diff[1] = cosmology.wa;

  log_L = -0.5*do_matrix_mult_invcov(2,inv, param_diff);
  return log_L;
 }

double scatter(double N200, double z)
{
  double lnM, M, a;
  lnM = 1./nuisance.cluster_Mobs_alpha*(log(N200)-nuisance.cluster_Mobs_lgN0-nuisance.cluster_Mobs_beta*log(1.+z));
  M = exp(lnM)*3.e14;
  a = 1./(1.+z);
  return scatter_lgN200_model_mz(M, a);
}

double scatter_fid(double N200, double z)
{
  double cluster_Mobs_lgN0 = 3.207;
  double cluster_Mobs_alpha = 0.993;
  double cluster_Mobs_beta = 0.0;
  double cluster_Mobs_sigma0 = 0.456;
  double cluster_Mobs_sigma_qm = 0.0; // SRD suggestion from Eduardo
//  double cluster_Mobs_sigma_qm = -0.169;
  double cluster_Mobs_sigma_qz = 0.0;

  double lnM, M, a;
  lnM = 1./cluster_Mobs_alpha*(log(N200)-cluster_Mobs_lgN0-cluster_Mobs_beta*log(1.+z));
  M = exp(lnM)*3.e14;
  a = 1./(1.+z);
  return cluster_Mobs_sigma0 + log(M/(3.e+14))*cluster_Mobs_sigma_qm + log(1./a)*cluster_Mobs_sigma_qz;
}

double log_L_SRD_LSST_CL_EDUARDO() // values for scatter function are "arbitrary" but based on some experience (Hironao/Eduarod)
{
  double log_L = 0.;
  double s_100_0p2 = scatter(100., 0.2);
  double s_30_0p2 = scatter(30., 0.2);
  double s_100_0p8 = scatter(100., 0.8); // remove if redshift dependent scatter is excluded in analysis
  double f_100_0p2 = scatter_fid(100., 0.2);
  double f_30_0p2 = scatter_fid(30., 0.2);
  double f_100_0p8 = scatter_fid(100., 0.8); // remove if redshift dependent scatter is excluded in analysis
  double sig = 0.1; // 0.05 weighting the giants

  log_L = -0.5*(exp(pow((s_100_0p2-f_100_0p2)/sig, 2.0))+ exp(pow((s_30_0p2-f_30_0p2)/sig, 2.0)) + exp(pow((s_100_0p8-f_100_0p8)/sig, 2.0)));

  return log_L;
}

double log_multi_like(double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q, double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope)
{
  int i,j,k,m=0,l;
  static double *pred;
  static double *ell;
  static double *ell_Cluster;
  static double darg;
  double chisqr,a,log_L_prior=0.0, log_L_GRS=0.0;

  if(ell==0){
    pred= create_double_vector(0, like.Ndata-1);
    ell= create_double_vector(0, like.Ncl-1);
    darg=(log(like.lmax)-log(like.lmin))/like.Ncl;
    for (l=0;l<like.Ncl;l++){
      ell[l]=exp(log(like.lmin)+(l+0.5)*darg);
    }
    ell_Cluster= create_double_vector(0, Cluster.lbin-1);
    darg=(log(Cluster.l_max)-log(Cluster.l_min))/Cluster.lbin;
    for (l=0;l<Cluster.lbin;l++){
      ell_Cluster[l]=exp(log(Cluster.l_min)+(l+0.5)*darg);
    }
  }
  if (set_cosmology_params(OMM,S8,NS,W0,WA,OMB,H0,MGSigma,MGmu)==0){
    printf("Cosmology out of bounds\n");
    return -1.0e8;
  }
  set_nuisance_shear_calib(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10);
  if (set_nuisance_shear_photoz(SP1,SP2,SP3,SP4,SP5,SP6,SP7,SP8,SP9,SP10,SPS1)==0){
    printf("Shear photo-z sigma too small\n");
    return -1.0e8;
  }
  if (set_nuisance_clustering_photoz(CP1,CP2,CP3,CP4,CP5,CP6,CP7,CP8,CP9,CP10,CPS1)==0){
    printf("Clustering photo-z sigma too small\n");
    return -1.0e8;
  }
  if (set_nuisance_ia(A_ia,beta_ia,eta_ia,eta_ia_highz,LF_alpha,LF_P,LF_Q,LF_red_alpha,LF_red_P,LF_red_Q)==0){
    printf("IA parameters out of bounds\n");
    return -1.0e8;
  }
  if (set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10)==0){
    printf("Bias out of bounds\n");
    return -1.0e8;
  }
  if (set_nuisance_cluster_Mobs(mass_obs_norm, mass_obs_slope, mass_z_slope, mass_obs_scatter_norm, mass_obs_scatter_mass_slope, mass_obs_scatter_z_slope)==0){
    printf("Mobs out of bounds\n");
    return -1.0e8;
  }

  //printf("like %le %le %le %le %le %le %le %le\n",cosmology.Omega_m, cosmology.Omega_v,cosmology.sigma_8,cosmology.n_spec,cosmology.w0,cosmology.wa,cosmology.omb,cosmology.h0);
  // printf("like %le %le %le %le\n",gbias.b[0][0], gbias.b[1][0], gbias.b[2][0], gbias.b[3][0]);
  // for (i=0; i<10; i++){
  //   printf("nuisance %le %le %le\n",nuisance.shear_calibration_m[i],nuisance.bias_zphot_shear[i],nuisance.sigma_zphot_shear[i]);
  // }

  log_L_prior=0.0;
  // if(like.Aubourg_Planck_BAO_SN==1) log_L_prior+=log_L_Planck_BAO_SN();
  // if(like.SN==1) log_L_prior+=log_L_SN();
  // if(like.BAO==1) log_L_prior+=log_L_BAO();
  // if(like.Planck==1) log_L_prior+=log_L_Planck();
  // if(like.IA!=0) log_L_prior+=log_L_ia();
  // if(like.IA!=0) log_L_prior+=log_like_f_red();
  // if(like.wlphotoz!=0) log_L_prior+=log_L_wlphotoz();
  // if(like.clphotoz!=0) log_L_prior+=log_L_clphotoz();
  // if(like.shearcalib==1) log_L_prior+=log_L_shear_calib();
  // if(like.clusterMobs==1) log_L_prior+=log_L_clusterMobs();

  // if(like.SRD_SN_Y10==1) log_L_prior+=log_L_SRD_SN_Y10_RENEE();
  // if(like.SRD_SN_Y1==1) log_L_prior+=log_L_SRD_SN_Y1_RENEE();
  // if(like.SRD_SL_Y10==1) log_L_prior+=log_L_SRD_SL_Y10_TOM();
  // if(like.SRD_SL_Y1==1) log_L_prior+=log_L_SRD_SL_Y1_TOM();

  // log_L_prior+=log_L_DESI();

  if(like.SRD==1 || like.SRD==6 || like.SRD==11) log_L_prior+=log_L_SRD_stage3();
  if(like.SRD==2 || like.SRD==6) log_L_prior+=1/(1+0.7*0.7)*log_L_SRD_SN_Y10_RENEE();
  if(like.SRD==3 || like.SRD==6) log_L_prior+=1/(1+0.62*0.62)*log_L_SRD_SL_Y10_TOM();
  if(like.SRD==4 || like.SRD==6) log_L_prior+=14202.63/18000.0*1/(1+0.62*0.62)*log_L_SRD_LSST_Y10_3x2pt();
  if(like.SRD==5 || like.SRD==6) log_L_prior+=14202.63/18000.0*1/(1+0.62*0.62)*log_L_SRD_LSST_Y10_clusters();

  if(like.SRD==7 || like.SRD==11) log_L_prior+=1/(1+0.7*0.7)*log_L_SRD_SN_Y1_RENEE();
  if(like.SRD==8 || like.SRD==11) log_L_prior+=1/(1+0.62*0.62)*log_L_SRD_SL_Y1_TOM();
  if(like.SRD==9 || like.SRD==11) log_L_prior+=11808.28/18000.0*1/(1+0.62*0.62)*log_L_SRD_LSST_Y1_3x2pt();
  if(like.SRD==10 || like.SRD==11) log_L_prior+=11808.28/18000.0*1/(1+0.62*0.62)*log_L_SRD_LSST_Y1_clusters();

  // printf("%d %d %d %d\n",like.BAO,like.wlphotoz,like.clphotoz,like.shearcalib);
  // printf("logl %le %le %le %le\n",log_L_shear_calib(),log_L_wlphotoz(),log_L_clphotoz(),log_L_clusterMobs());
  // int start=0;

  // if(like.shear_shear==1) {
  //   set_data_shear(like.Ncl, ell, pred, start);
  //   start=start+like.Ncl*tomo.shear_Npowerspectra;
  // }
  // if(like.shear_pos==1){
  //   set_data_ggl(like.Ncl, ell, pred, start);
  //   start=start+like.Ncl*tomo.ggl_Npowerspectra;
  // }
  // if(like.pos_pos==1){
  //   set_data_clustering(like.Ncl,ell,pred, start);
  //   start=start+like.Ncl*tomo.clustering_Npowerspectra;
  // }
  // if(like.clusterN==1){
  //   set_data_cluster_N(pred,start);
  //   start=start+tomo.cluster_Nbin*Cluster.N200_Nbin;
  // }
  // if(like.clusterWL==1){
  //   set_data_cgl(ell_Cluster,pred, start);
  // }
  chisqr=0.0;
  // for (i=0; i<like.Ndata; i++){
  //   for (j=0; j<like.Ndata; j++){
  //     a=(pred[i]-data_read(1,i))*invcov_read(1,i,j)*(pred[j]-data_read(1,j));
  //     chisqr=chisqr+a;
  //   }
  //  //  if (fabs(data_read(1,i)/pred[i]-1.0) >1.e-4){
  //  // printf("%d %le %le %le\n",i,data_read(1,i),pred[i],data_read(1,i)/pred[i]);
  //  //  }
  // }
  // if (chisqr<0.0){
  //   printf("errror: chisqr < 0\n");
  // }
  // if (chisqr<-1.0) exit(EXIT_FAILURE);

  return -0.5*chisqr+log_L_prior;
}

void compute_data_vector(char *details, double OMM, double S8, double NS, double W0,double WA, double OMB, double H0, double MGSigma, double MGmu, double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8, double B9, double B10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double SPS1, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double CPS1, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double beta_ia, double eta_ia, double eta_ia_highz, double LF_alpha, double LF_P, double LF_Q, double LF_red_alpha, double LF_red_P, double LF_red_Q, double mass_obs_norm, double mass_obs_slope, double mass_z_slope, double mass_obs_scatter_norm, double mass_obs_scatter_mass_slope, double mass_obs_scatter_z_slope)
{

  int i,j,k,m=0,l;
  static double *pred;
  static double *ell;
  static double *ell_Cluster;
  static double darg;
  double chisqr,a,log_L_prior=0.0;

  if(ell==0){
    pred= create_double_vector(0, like.Ndata-1);
    ell= create_double_vector(0, like.Ncl-1);
    darg=(log(like.lmax)-log(like.lmin))/like.Ncl;
    for (l=0;l<like.Ncl;l++){
      ell[l]=exp(log(like.lmin)+(l+0.5)*darg);
    }
    ell_Cluster= create_double_vector(0, Cluster.lbin-1);
    darg=(log(Cluster.l_max)-log(Cluster.l_min))/Cluster.lbin;
    for (l=0;l<Cluster.lbin;l++){
      ell_Cluster[l]=exp(log(Cluster.l_min)+(l+0.5)*darg);
    }
  }
// for (l=0;l<like.Ncl;l++){
//   printf("%d %le\n",i,ell[l]);
// }
  set_cosmology_params(OMM,S8,NS,W0,WA,OMB,H0,MGSigma,MGmu);
  set_nuisance_shear_calib(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10);
  set_nuisance_shear_photoz(SP1,SP2,SP3,SP4,SP5,SP6,SP7,SP8,SP9,SP10,SPS1);
  set_nuisance_clustering_photoz(CP1,CP2,CP3,CP4,CP5,CP6,CP7,CP8,CP9,CP10,CPS1);
  set_nuisance_ia(A_ia,beta_ia,eta_ia,eta_ia_highz,LF_alpha,LF_P,LF_Q,LF_red_alpha,LF_red_P,LF_red_Q);
  set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10);
  set_nuisance_cluster_Mobs(mass_obs_norm, mass_obs_slope, mass_z_slope, mass_obs_scatter_norm, mass_obs_scatter_mass_slope, mass_obs_scatter_z_slope);

  int start=0;
  if(like.shear_shear==1) {
    set_data_shear(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.shear_Npowerspectra;
  }
  if(like.shear_pos==1){
    //printf("ggl\n");
    set_data_ggl(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.ggl_Npowerspectra;
  }
  if(like.pos_pos==1){
    //printf("clustering\n");
    set_data_clustering(like.Ncl,ell,pred, start);
    start=start+like.Ncl*tomo.clustering_Npowerspectra;
  }

  if(like.clusterN==1){
    set_data_cluster_N(pred,start);
    start= start+tomo.cluster_Nbin*Cluster.N200_Nbin;
  }
  if(like.clusterWL==1){
    set_data_cgl(ell_Cluster,pred, start);
  }
  FILE *F;
  char filename[300];
  if (strstr(details,"FM") != NULL){
    sprintf(filename,"%s",details);
  }
  else {sprintf(filename,"datav/%s_%s",like.probes,details);}
  F=fopen(filename,"w");
  for (i=0;i<like.Ndata; i++){
    fprintf(F,"%d %le\n",i,pred[i]);
    //printf("%d %le\n",i,pred[i]);
  }
  fclose(F);
}

double write_vector_wrapper(char *details, input_cosmo_params ic, input_nuisance_params in)
{
  compute_data_vector(details, ic.omega_m, ic.sigma_8, ic.n_s, ic.w0, ic.wa, ic.omega_b, ic.h0, ic.MGSigma, ic.MGmu,
    in.bias[0], in.bias[1], in.bias[2], in.bias[3],in.bias[4], in.bias[5], in.bias[6], in.bias[7],in.bias[8], in.bias[9],
    in.source_z_bias[0], in.source_z_bias[1], in.source_z_bias[2], in.source_z_bias[3], in.source_z_bias[4],
    in.source_z_bias[5], in.source_z_bias[6], in.source_z_bias[7], in.source_z_bias[8], in.source_z_bias[9],
    in.source_z_s,
    in.lens_z_bias[0], in.lens_z_bias[1], in.lens_z_bias[2], in.lens_z_bias[3], in.lens_z_bias[4],
    in.lens_z_bias[5], in.lens_z_bias[6], in.lens_z_bias[7], in.lens_z_bias[8], in.lens_z_bias[9],
    in.lens_z_s,
    in.shear_m[0], in.shear_m[1], in.shear_m[2], in.shear_m[3], in.shear_m[4],
    in.shear_m[5], in.shear_m[6], in.shear_m[7], in.shear_m[8], in.shear_m[9],
    in.A_ia, in.beta_ia, in.eta_ia, in.eta_ia_highz,
    in.lf[0], in.lf[1], in.lf[2], in.lf[3], in.lf[4], in.lf[5],
    in.m_lambda[0], in.m_lambda[1], in.m_lambda[2], in.m_lambda[3],
    in.m_lambda[4], in.m_lambda[5]);
  return 0;
}

double log_like_wrapper(input_cosmo_params ic, input_nuisance_params in)
{
  double like = log_multi_like(ic.omega_m, ic.sigma_8, ic.n_s, ic.w0, ic.wa, ic.omega_b, ic.h0, ic.MGSigma, ic.MGmu,

    in.bias[0], in.bias[1], in.bias[2], in.bias[3],in.bias[4], in.bias[5], in.bias[6], in.bias[7],in.bias[8], in.bias[9],
    in.source_z_bias[0], in.source_z_bias[1], in.source_z_bias[2], in.source_z_bias[3], in.source_z_bias[4],
    in.source_z_bias[5], in.source_z_bias[6], in.source_z_bias[7], in.source_z_bias[8], in.source_z_bias[9],
    in.source_z_s,
    in.lens_z_bias[0], in.lens_z_bias[1], in.lens_z_bias[2], in.lens_z_bias[3], in.lens_z_bias[4],
    in.lens_z_bias[5], in.lens_z_bias[6], in.lens_z_bias[7], in.lens_z_bias[8], in.lens_z_bias[9],
    in.lens_z_s,
    in.shear_m[0], in.shear_m[1], in.shear_m[2], in.shear_m[3], in.shear_m[4],
    in.shear_m[5], in.shear_m[6], in.shear_m[7], in.shear_m[8], in.shear_m[9],
    in.A_ia, in.beta_ia, in.eta_ia, in.eta_ia_highz,
    in.lf[0], in.lf[1], in.lf[2], in.lf[3], in.lf[4], in.lf[5],
    in.m_lambda[0], in.m_lambda[1], in.m_lambda[2], in.m_lambda[3],
    in.m_lambda[4], in.m_lambda[5]);
  return like;
}

 int main(int argc, char** argv)
{
  int hit=atoi(argv[1]);
  char filename[500],arg1[400],arg2[400];
  int t;
  t=hit;
  double area_table[12]={7500.0,13000.0,16000.0,10000.0,15000.0,20000.0,10000.0,15000.0,20000.0,10000.0,15000.0,20000.0};
  double nsource_table[12]={9.8,12.1,15.1,15.1,18.9,23.5,20.3,23.5,26.9,26.9,30.8,35.0};
  double nlens_table[12]={15.0,20.0,25.0,25.0,32.0,41.0,35.0,41.0,48.0,48.0,57.0,67.0};

  char survey_designation[12][200]={"LSST_Y1","LSST_Y1","LSST_Y1","LSST_Y3","LSST_Y3","LSST_Y3","LSST_Y6","LSST_Y6","LSST_Y6","LSST_Y10","LSST_Y10","LSST_Y10"};
  char source_zfile[12][400]={"WL_zdistri_model0_z0=1.940000e-01_alpha=8.830000e-01","WL_zdistri_model1_z0=1.900000e-01_alpha=8.620000e-01","WL_zdistri_model2_z0=1.860000e-01_alpha=8.410000e-01","WL_zdistri_model3_z0=1.860000e-01_alpha=8.410000e-01","WL_zdistri_model4_z0=1.830000e-01_alpha=8.210000e-01","WL_zdistri_model5_z0=1.790000e-01_alpha=8.000000e-01","WL_zdistri_model6_z0=1.810000e-01_alpha=8.140000e-01","WL_zdistri_model7_z0=1.790000e-01_alpha=8.000000e-01","WL_zdistri_model8_z0=1.760000e-01_alpha=7.860000e-01","WL_zdistri_model9_z0=1.760000e-01_alpha=7.860000e-01","WL_zdistri_model10_z0=1.740000e-01_alpha=7.720000e-01","WL_zdistri_model11_z0=1.710000e-01_alpha=7.590000e-01"};
  char lens_zfile[12][400]={"LSS_zdistri_model0_z0=2.590000e-01_alpha=9.520000e-01","LSS_zdistri_model1_z0=2.610000e-01_alpha=9.370000e-01","LSS_zdistri_model2_z0=2.640000e-01_alpha=9.250000e-01","LSS_zdistri_model3_z0=2.640000e-01_alpha=9.250000e-01","LSS_zdistri_model4_z0=2.680000e-01_alpha=9.150000e-01","LSS_zdistri_model5_z0=2.740000e-01_alpha=9.070000e-01","LSS_zdistri_model6_z0=2.700000e-01_alpha=9.120000e-01","LSS_zdistri_model7_z0=2.740000e-01_alpha=9.070000e-01","LSS_zdistri_model8_z0=2.780000e-01_alpha=9.030000e-01","LSS_zdistri_model9_z0=2.780000e-01_alpha=9.030000e-01","LSS_zdistri_model10_z0=2.830000e-01_alpha=9.000000e-01","LSS_zdistri_model11_z0=2.880000e-01_alpha=8.980000e-01"};
  int Ntomo_lens[12]={5,5,5,7,7,7,9,9,9,10,10,10};

  init_cosmo();
  init_fisher_precision();
  init_binning_fourier(20,20.0,15000.0,3000.0,21.0,5,Ntomo_lens[t]);
  init_survey(survey_designation[t]);
  sprintf(arg1,"zdistris/%s",source_zfile[t]);
  sprintf(arg2,"zdistris/%s",lens_zfile[t]);
  init_galaxies(arg1,arg2,"gaussian","gaussian","SRD");
  init_clusters();
  init_IA("NLA_HF", "GAMA");
  init_probes("3x2pt_clusterN_clusterWL");
  init_priors("none","none","none","none");

  sprintf(filename,"%s_area%le_ng%le_nl%le",survey_designation[t],area_table[t],nsource_table[t],nlens_table[t]);
  if(strstr(survey_designation[t],"LSST_Y1") != NULL) compute_data_vector(filename,0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.413566e+00,1.567919e+00,1.731037e+00,1.900583e+00,2.074809e+00,1.413566e+00,1.567919e+00,1.731037e+00,1.900583e+00,2.074809e+00,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.03,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);

  if(strstr(survey_designation[t],"LSST_Y3") != NULL) compute_data_vector(filename,0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.392398e+00,1.500535e+00,1.613747e+00,1.731037e+00,1.851600e+00,1.974761e+00,2.100003e+00,1.731037e+00,1.900583e+00,2.074809e+00,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.03,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);

  if(strstr(survey_designation[t],"LSST_Y6") != NULL) compute_data_vector(filename,0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.380752e+00,1.463865e+00,1.550281e+00,1.639495e+00,1.731037e+00,1.824565e+00,1.919738e+00,2.016299e+00,2.114025e+00,2.074809e+00,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.03,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);

  if(strstr(survey_designation[t],"LSST_Y10") != NULL) compute_data_vector(filename,0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.376695e+00,1.451179e+00,1.528404e+00,1.607983e+00,1.689579e+00,1.772899e+00,1.857700e+00,1.943754e+00,2.030887e+00,2.118943e+00,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.03,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);
  return 0;
}


text:
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_histogram.h>
#include"mpi.h"

const double PI = 3.14159;

double rz(double red);

struct galaxy {
  double x, y, z, d, w;
};

void N2divide( int * n, int N, int Ntot);

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

  FILE * DD;
  FILE * PP;

  const int Nr = 150;
  const int Nm = 100;
  const double RSTEP = 1.0;
  const double MSTEP = 0.01;

  const int maxNgal = 40000000;
  const double RMIN = 0.0;
  const double RMAX = 150.0;

  int *n;

  char DDname[400];
  char PPname[400];

  int bit;
  int bittot;

  sprintf(DDname,argv[1]);

  MPI_Init (&argc, &argv);
  int numproc;
  int totproc;
  MPI_Comm_rank (MPI_COMM_WORLD, &numproc);
  MPI_Comm_size (MPI_COMM_WORLD, &totproc);

  n = (int *) calloc(totproc + 1, sizeof(int));

  sprintf(PPname,"%s.RR",argv[2]);

  DD = fopen(DDname,"r");
  printf("%s %d\n",DDname,numproc);
  struct galaxy * gal;
  double * DDcount;
  double * DDcount_in_place;
  gal = (struct galaxy*)malloc(maxNgal*sizeof(struct galaxy));
  DDcount = (double*)calloc(Nr*Nm,sizeof(double));
  DDcount_in_place = (double*)calloc(Nr*Nm,sizeof(double));
  struct galaxy * galp;
  galp = gal;

  int N = 0;
  int Ntest = 0;
  double sumW = 0;
  double sumW_in_place = 0;
  double ra, dec, red, dist;
  double weight = 1.0;
  double weight1;
  double weight2;
  double weight3;
  double weight4;
  double weightFKP;
  int indexFKP;
  double ddummy;
  char line[2000];
  for (int i = 0; i < 0; i ++) {
  fgets (line, 2000, DD);
  printf("%s\n", line);
  }
  int counter = 0;
  while (fscanf(DD,"%le %le %le\n",&ra,&dec,&red) != EOF) {
    counter ++;
    if (counter%3 != 0) {
      continue;
    }
    ra *= PI/180.0;
    dec *= PI/180.0;
    dist = rz(red);
    galp->x = dist*cos(dec)*cos(ra);
    galp->y = dist*cos(dec)*sin(ra);
    galp->z = dist*sin(dec);
    galp->d = dist;
//    galp->w = weight1 * weight2;
    galp++;
    N++;
  }
  fclose(DD);
  n[0] = 1;
  N2divide(n, totproc, N);
  n[totproc] = N;
  n[0] = 0;
  bit = n[numproc];
  bittot = n[numproc + 1];
  for (int i = 0; i < totproc + 1; i++) {
    printf("n[%d] = %d\n", i, n[i]);
  }
  double mu;
  double d1, d2, d3;
  double x1, y1, z1, x2, y2, z2, gd1, gd2, gw1, gw2;
  double r2, rat;
  double rr;
  double xb, yb, zb, db2;
  struct galaxy * galp1, * galp2;
  galp1 = gal;
  galp1 += bit;
  for (int i = bit; i < bittot; i++) {
    x1 = galp1->x; y1 = galp1->y; z1 = galp1->z; gd1 = galp1->d; //gw1 = galp1->w;
    galp2 = gal + i + 1;
    for (int j = i + 1; j < N; j++) {
       x2 = galp2->x; y2 = galp2->y; z2 = galp2->z; gd2 = galp2->d; //gw2 = galp2->w;
       //sumW += gw1*gw2;
       sumW += 1.0;
       d1 = x1 - x2;
       d2 = y1 - y2;
       d3 = z1 - z2;
       r2 = d1*d1 + d2*d2 + d3*d3;
       if (r2 > RMAX*RMAX || r2 < RMIN*RMIN) {
	   galp2++;
         continue;
       }
       rat = gd1/gd2;
       xb = x1 + x2*rat;
       yb = y1 + y2*rat;
       zb = z1 + z2*rat;
       db2 = xb*xb + yb*yb + zb*zb;
       mu = fabs((xb*d1 + yb*d2 + zb*d3)/sqrt(r2)/sqrt(db2));
       rr = sqrt(r2);
	int binr = (int)((rr - RMIN)/RSTEP);
	int binm = (int)(mu/MSTEP);
	if (binr >= 0 && binm >= 0 && binr < Nr && binm < Nm){
	int ind = binr + Nr*binm;
//	DDcount[ind] += gw1*gw2;
        DDcount[ind] += 1.0;
	}
       galp2++;
    }
    galp1++;
  }
  free(gal);

  MPI_Reduce(DDcount, DDcount_in_place, Nr*Nm, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  MPI_Reduce(&sumW, &sumW_in_place, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
  if (numproc == 0) {
    PP = fopen(PPname,"w");
    printf("%s\n", PPname);
    fprintf(PP,"# weighted number: %lf\n",sumW_in_place);
    fprintf(PP,"# RBINS: %d\n",Nr);
    fprintf(PP,"# MBINS: %d\n",Nm);
    int k, l;
    for (k = 0; k < Nm; k++) {
      for (l = 0; l < Nr; l++) {
      fprintf(PP,"%lf ",DDcount_in_place[k*Nr + l]);
      }
      fprintf(PP,"\n");
    }
    fclose(PP);
  }

  free(DDcount);
  free(DDcount_in_place);

  free(n);
  MPI_Finalize ();

  return 0;
}

double f (double x, void * p) {
  double Om = 0.310;
  double ff = 2997.92458/sqrt(Om*(1.0 + x)*(1.0 + x)*(1.0 + x) + 1.0 - Om);
  return ff;
}

double rz (double red) {
  gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
  double result, error;
  gsl_function F;
  F.function = &f;
  gsl_integration_qags(&F, 0, red, 0, 1e-7, 1000, w, &result, &error);
  gsl_integration_workspace_free (w);
  return result;
}

void N2divide( int * n, int N, int Ntot) {
   double dummy;
   for (int i = 1; i < N; i++) {
     dummy = double(Ntot) - 0.5;
     dummy -=
sqrt((2*double(Ntot)-1)*(2*double(Ntot)-1)-4*((2*double(Ntot)-n[i-1])*(n[i-1]-1)+double(Ntot)*(double(Ntot)-1)/N))/2.0;
     n[i] = int(dummy);
   }

   return;

}

ext: c

max_stars_repo_name: echaussidon/LSS

max_issues_repo_name: echaussidon/LSS

max_forks_repo_name: echaussidon/LSS


text:
/**
 * \author Sylvain Marsat, University of Maryland - NASA GSFC
 *
 * \brief C code headers for the implementation of the Fourier domain response for LIGO-VIRGO detectors
 *
 *
 */

#ifndef _LLVFDRESPONSE_H
#define _LLVFDRESPONSE_H

#define _XOPEN_SOURCE 500

#ifdef __GNUC__
#define UNUSED __attribute__ ((unused))
#else
#define UNUSED
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <unistd.h>
#include <getopt.h>
#include <stdbool.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_complex.h>

#include "constants.h"
#include "LLVgeometry.h"
#include "struct.h"
#include "waveform.h"
#include "timeconversion.h"

/**************************************************/
/**************** Prototypes **********************/

/* Function to convert string input network string to Networktag */
Networktag ParseNetworktag(char* string);

/* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LLV response (for a given detector), for given values of the inclination, position in the sky and polarization angle */
int LLVSimFDResponse(
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,  /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
  struct tagListmodesCAmpPhaseFrequencySeries **lists,    /* Output: list of contribution of each mode in the detector signal, in Frequency-domain amplitude and phase form, for the given detector and sky position */
  const double gpstime,                                   /* GPS time (s) when the signal at coalescence reaches geocenter */
  const double ra,                                        /* Position in the sky: J2000.0 right ascension (rad) */
  const double dec,                                       /* Position in the sky: J2000.0 declination (rad) */
  const double inclination,                               /* Inclination of the source (rad) */
  const double psi,                                       /* Polarization angle (rad) */
  const Detectortag tag);                                 /* Tag identifying the detector */

  /* Core function processing a signal (in the form of a list of modes) through the Fourier-domain LLV response, for given values of the inclination, position in the sky and polarization angle */
  int LLVSimFDResponse3Det(
    struct tagListmodesCAmpPhaseFrequencySeries **listDet1,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the detector 1 */
    struct tagListmodesCAmpPhaseFrequencySeries **listDet2,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the detector 2 */
    struct tagListmodesCAmpPhaseFrequencySeries **listDet3,  /* Output: list of contribution of each mode in Frequency-domain amplitude and phase form, in the detector 3 */
    struct tagListmodesCAmpPhaseFrequencySeries **list,      /* Input: list of modes in Frequency-domain amplitude and phase form as produced by the ROM */
    const double gpstime,                                    /* GPS time (s) when the signal at coalescence reaches geocenter */
    const double ra,                                            /* First angle for the position in the sky */
    const double dec,                                           /* Second angle for the position in the sky */
    const double inclination,                                   /* Inclination of the source */
    const double psi,                                           /* Polarization angle */
    const Networktag tag);                               /* Selector for the detector network */

/* Function setting the response matrix of a given detector, in cartesian coordinates */
void SetMatrixD(
  gsl_matrix* D,                       /* Output: matrix of the detector response Dij */
  const Detectortag tag);              /* Tag identifying the detector */
/* Function setting the position of a detector, in cartesian coordinates */
void SetVectorXd(
  gsl_vector* Xd,                      /* Output: position vector of the detector */
  const Detectortag tag);              /* Tag identifying the detector */

/* Function setting the cartesian coordinates of the wave frame vectors (X,Y,Z), given the position in the sky and polarization */
void SetVectorsXYZ(
  gsl_vector* X,                       /* Output: cartesian vector of the wave frame unit vector X */
  gsl_vector* Y,                       /* Output: cartesian vector of the wave frame unit vector Y */
  gsl_vector* Z,                       /* Output: cartesian vector of the wave frame unit vector Z */
  const double theta,                  /* First angle for the position in the sky (Earth-based spherical angle) */
  const double phi,                    /* Second angle for the position in the sky (Earth-based spherical angle) */
  const double psi);                   /* Polarization angle */

#if 0
{ /* so that editors will match succeeding brace */
#elif defined(__cplusplus)
}
#endif

#endif /* _LLVFDRESPONSE_H */

ext: h

max_stars_repo_name: titodalcanton/flare

max_issues_repo_name: titodalcanton/flare

max_forks_repo_name: titodalcanton/flare


text:
#ifndef __DIRECT_KNN_H__
#define __DIRECT_KNN_H__

#include <mpi.h>
#include <utility>
#include <cmath>
#include <vector>
#include <string>
#include <blas.h>
#include <omp.h>
#include <utility>

#define KNN_MAX_BLOCK_SIZE 512
#define KNN_MAX_MATRIX_SIZE 2e7L

// MPI message tags
enum {
  TAG_R,
  TAG_Q,
  TAG_KMIN,
  TAG_ID,
  TAG_SIZE
};

enum {
  VERTICAL,
  HORIZONTAL
};

/**
 * Options struct used to specify how a rectangular-partitioned query should be performed.
 */
struct directQueryParams
{
	char queryType; ///< Perform a k-nearest ('K') or r-near ('R') search.
	int k; ///< Number of neighbors to find (if queryType is 'K').
	double r; ///< Radius for r-near search (if queryType is 'R').
	int refParts; ///< Partition reference points into this many approximately equal parts.
	int queryParts; ///< Partition query points into this many approximately equal parts.
};

struct directQueryResults
{
	std::pair<double, long> *neighbors;
	long *neighborCounts;
};

class maxheap_comp {
    public:
        bool operator() (const std::pair<double, long> &a, const std::pair<double, long> &b) {
            double diff = fabs(a.first-b.first)/a.first;
            if( std::isinf(diff) || std::isnan(diff) ) {      // 0/0 or x/0
                return a.first < b.first;
            }
            if( diff < 1.0e-8 ) {
                return a.second < b.second;
            }
            return a.first < b.first;
        }
};

namespace knn {

   /**
    * Compute the k nearest neighbors of each query point (OBSOLETE: For debugging purposes only; use directKQueryLowMem instead)..
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
    * \param n number of local reference point
    * \param m number of local query points
    * \param k number of neighbors to find for each query point.
    * \param dim dimensionality of reference/query points
    * \return An array of pairs containing the *squared* distance to each point and its
    *	index in the local array.  Note that this is the index of the /point/ rather
    *	than the index of its location in the actual array of doubles.  To access the
    *	point, use (point index)*dim to index ref.
    */
   std::pair<double, long> *directKQuery
              ( double *ref, double *query, long n, long m, long k, int dim);

   /**
    * Compute the k nearest neighbors of each query point (OBSOLETE: For debugging purposes only; use directKQueryLowMem instead)..
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
    * \param n number of local reference point
    * \param m number of local query points
    * \param k number of neighbors to find for each query point.
    * \param dim dimensionality of reference/query points
    * \param result [out] An array of m*k pairs (allocated by caller) containing the
    * *squared* distance to each point and its index in the local array.
    * Note that this is the index of the /point/ rather than the index of its location
    * in the actual array of doubles.  To access the point, use (point index)*dim to index ref.
    * Results for query point i begin at result[i*k].
    * \param dist An array used for calculating inter-point distances. If NULL or unspecified,
    * the function allocates and deallocates this array internally.  Otherwise, it must be
    * allocated size at least n*B, where B is the value returned by getBlockSize, and
    * deallocated at some point after the function returns.
    * \param sqnormr An array used for calculating inter-point distances. If NULL or unspecified,
    * the function allocates and deallocates this array internally.  Otherwise, it must be
    * allocated size at least n and deallocated at some point after the function returns.
    * \param sqnormq An array used for calculating inter-point distances. If NULL or unspecified,
    * the function allocates and deallocates this array internally.  Otherwise, it must be
    * allocated size at least B, where B is the value returned by getBlockSize, and
    * deallocated at some point after the function returns.
    */
 void directKQuery_small_a2a( double *ref, long n, int dim, int k,
                              std::pair<double, long> *result,
                              double *dist = NULL, double *sqnormr = NULL, double *sqnormq = NULL);

   /**
    * Compute the k nearest neighbors of each query point with a smaller memory footprint and usually better performance.
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
    * \param n number of local reference point
    * \param m number of local query points
    * \param k number of neighbors to find for each query point.
    * \param dim dimensionality of reference/query points
    * \param result [out] An array of m*k pairs (allocated by caller) containing the
    * *squared* distance to each point and its
    * index in the local array.  Note that this is the index of the /point/ rather
    * than the index of its location in the actual array of doubles.  To access the
    * point, use (point index)*dim to index ref. Results for query point i begin at result[i*k].
    * \param dist An array used for calculating inter-point distances. If NULL or unspecified,
    * the function allocates and deallocates this array internally.  Otherwise, it must be
    * allocated size at least n*B, where B is the value returned by getBlockSize, and
    * deallocated at some point after the function returns.
    * \param sqnormr An array used for calculating inter-point distances. If NULL or unspecified,
    * the function allocates and deallocates this array internally.  Otherwise, it must be
    * allocated size at least n and deallocated at some point after the function returns.
    * \param sqnormq An array used for calculating inter-point distances. If NULL or unspecified,
    * the function allocates and deallocates this array internally.  Otherwise, it must be
    * allocated size at least B, where B is the value returned by getBlockSize, and
    * deallocated at some point after the function returns.
    */
   void directKQueryLowMem
              ( double *ref, double *query, long n, long m, long k, int dim, std::pair<double, long> *result,
                double *dist = NULL, double* sqnormr = NULL, double* sqnormq = NULL  );

   /**
    * Find all points in ref that lie within distance sqrt(R) of each
    * point in query.
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
    * \param n number of local reference point
    * \param m number of local query points
    * \param R The *square* of the search radius.
    * \param dim dimensionality of reference/query points
    * \param glob_ids An array of length n containing the global id of each ref point.
    * \param neighbor_count [out] An array of length m containing the number of
    * neighbors found for the corresponding query point (allocated internally).
    * \param neighbors [out] An array of length sum(neighbor_count) containing
    * pairs of ids and squared distances for each neighbor found (allocated
    * internally).
    */
   void directRQuery
           ( double *ref, double *query, long n, long m, double R, int dim, long* glob_ids,
             int **neighbor_count, std::pair<double, long> **neighbors );

  /**
   * Find all points in ref that lie within distance sqrt(R[i]) of query point i.
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
   * \param n number of local reference point
   * \param m number of local query points
   * \param R Array of length m, the *square* of the search radius for each point.
   * \param dim dimensionality of reference/query points
   * \param glob_ids An array of length n containing the global id of each ref point.
   * \param neighbor_count [out] An array of length m containing the number of
   * neighbors found for the corresponding query point (allocated internally).
   * \param neighbors [out] An array of length sum(neighbor_count) containing
   * pairs of ids and squared distances for each neighbor found (allocated
   * internally).
   */
  void directRQueryIndividual
  ( double *ref, double *query, long n, long m, double *R, int dim, long* glob_ids,
   int **neighbor_count, std::pair<double, long> **neighbors );

  /**
   * Find all points in ref that lie within distance sqrt(R[i]) of query point i.
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
   * \param n number of local reference point
   * \param m number of local query points
   * \param R Array of length m, the *square* of the search radius for each point.
   * \param dim dimensionality of reference/query points
   * \param glob_ids An array of length n containing the global id of each ref point.
   * \param neighbor_count [out] An array of length m containing the number of
   * neighbors found for the corresponding query point (allocated internally).
   * \param neighbors [out] An array of length sum(neighbor_count) containing
   * pairs of ids and squared distances for each neighbor found (allocated
   * internally).
   */
void directRQueryIndividualK
( double *ref, double *query, long n, long m, int k, double *R, int dim, long *glob_ids,
 int **neighbor_count, std::pair<double, long> **neighbors  );

   /**
    * Computes the distances between all query points and all reference points.
    * \param ref an n*dim-length array of data points.
    * \param query an m*dim-length array of query points.
    * \param n number of reference point
    * \param m number of query points
    * \param dim dimensionality of reference/query points
    * \param dist [out] m*n matrix of the *squared* distances from query points to reference points.
    * \param sqnormr An array of length n used for storing the squared norm of each reference point.
    * If NULL or omitted, allocated and deallocated internally.
    * \param sqnromq An array of length m used for storing the squared norm of each reference point.
    * If NULL or omitted, allocated and deallocated internally.
    * \param useSqnormrInput If true, use precomputed squared norm values contained in sqnormr when
    * this function is called; if false, compute squared norms internally.
    */
   void compute_distances
              ( double *ref, double *query, long n, long m, int dim, double *dist,
                  double* sqnormr = NULL, double* sqnormq = NULL, bool useSqnormrInput = false );

   /**
    * Computes the square of the norm of each row of a matrix.
    * \param a The input matrix (row-major).
    * \param n The number of rows of a.
    * \param dim The number of columns of a.
    * \param b [out] The output array (length n).
    */
   void sqnorm ( double *a, long n, int dim, double *b);

   /**
    * Find the k-nearest neighbors of each query point using distributed cyclic-shift algorithm.
    * \param ref an nlocal*dim-length array of data points.
    * \param query an mlocal*dim-length array of query points.
    * \param glob_ids An array of length nlocal containing the unique ID of each reference point stored locally.
    * \param nlocal Number of local reference points.
    * \param mlocal Number of local query points.
    * \param k Number of neighbors to find for each query point.
    * \param dim Dimensionality of reference/query points.
    * \param comm The MPI communicator to use.
    * \return An array of pairs containing the distance to each neighbor and its index;
    * neighbors for point i begin at location i*k.
    * \note All input arrays MUST be allocated with malloc (or equivalent) and deallocated with free;
    * new and delete[] may result in a segmentation fault.
    */
   std::pair<double, long> * dist_directKQuery( double* &ref, double *query,
   						 long* &glob_ids,
   						 int nlocal, int mlocal,
   						 int k,
   						 int dim,
   						 MPI_Comm comm );

   /**
    * Find all reference points lying within distance sqrt(R) of each query point.
    * \param ref an nlocal*dim-length array of data points.
    * \param query an mlocal*dim-length array of query points.
    * \param n Number of global reference points
    * \param m Number of global query points
    * \param nlocal Number of local reference points
    * \param mlocal Number of local query points
    * \param R The *square* of the search radius.
    * \param dim Dimensionality of reference/query points
    * \param glob_ids An array of length nlocal containing the unique ID of each reference point stored locally.
    * \param rneighbors [out] An externally allocated array of mlocal vectors, used to store the neighbors of each
    * query point.
    * \param comm The MPI communicator to use.
    */
    void dist_directRQuery
              ( double* &ref, double *query, long n, long m, int nlocal, int mlocal, double R,
                int dim, long* &glob_ids, std::vector< std::pair<double, long> >* rneighbors,
                MPI_Comm comm  );

   /**
    * Find all reference points lying within distance sqrt(R) of each query point,
    * up to a maximum of max_neighbors (not necessarily the nearest).
    * \param ref an nlocal*dim-length array of data points.
    * \param query an mlocal*dim-length array of query points.
    * \param n Number of global reference points
    * \param m Number of global query points
    * \param nlocal Number of local reference points
    * \param mlocal Number of local query points
    * \param R The *square* of the search radius.
    * \param dim Dimensionality of reference/query points
    * \param glob_ids An array of length nlocal containing the unique ID of each reference point stored locally.
    * \param rneighbors [out] An externally allocated array of mlocal vectors, used to store the neighbors of each
    * query point.
    * \param max_neighbors The maximum number of neighbors returned for any given query point.
    * \param comm The MPI communicator to use.
    */
    void dist_directRQuery
    ( double* &ref, double *query, long n, long m, int nlocal, int mlocal, double R,
     int dim, long* &glob_ids, std::vector< std::pair<double, long> >* rneighbors, int max_neighbors,
     MPI_Comm comm  );

  /**
   * Find all reference points lying within distance sqrt(R[i]) of query point i.
   * \param ref an nlocal*dim-length array of data points.
   * \param query an mlocal*dim-length array of query points.
   * \param n Number of global reference points
   * \param m Number of global query points
   * \param nlocal Number of local reference points
   * \param mlocal Number of local query points
   * \param R Array of length m, the *square* of the search radius for each query point.
   * \param dim Dimensionality of reference/query points
   * \param glob_ids An array of length nlocal containing the unique ID of each reference point stored locally.
   * \param rneighbors [out] An externally allocated array of mlocal vectors, used to store the neighbors of each
   * query point.
   * \param comm The MPI communicator to use.
   */
  void dist_directRQueryIndividual
  ( double* &ref, double *query, long n, long m, int nlocal, int mlocal, double *R,
   int dim, long* &glob_ids, std::vector< std::pair<double, long> >* rneighbors,
   MPI_Comm comm  );

  /**
   * Find all reference points lying within distance sqrt(R[i]) of query point i.
   * \param ref an nlocal*dim-length array of data points.
   * \param query an mlocal*dim-length array of query points.
   * \param n Number of global reference points
   * \param m Number of global query points
   * \param nlocal Number of local reference points
   * \param mlocal Number of local query points
   * \param R Array of length m, the *square* of the search radius for each query point.
   * \param dim Dimensionality of reference/query points
   * \param glob_ids An array of length nlocal containing the unique ID of each reference point stored locally.
   * \param rneighbors [out] An externally allocated array of mlocal vectors, used to store the neighbors of each
   * query point.
   * \param comm The MPI communicator to use.
   */
   void dist_directRQueryIndividualK
( double* &ref, double *query, long n, long m, int nlocal, int mlocal, double *R, int k,
 int dim, long* &glob_ids, std::vector< std::pair<double, long> >* rneighbors,
 MPI_Comm comm  );

  /**
   * Find all reference points lying within distance sqrt(R[i]) of query point i,
   * up to a maximum of max_neighbors points (not necessarily the nearest).
   * \param ref an nlocal*dim-length array of data points.
   * \param query an mlocal*dim-length array of query points.
   * \param n Number of global reference points
   * \param m Number of global query points
   * \param nlocal Number of local reference points
   * \param mlocal Number of local query points
   * \param R Array of length m, the *square* of the search radius for each query point.
   * \param dim Dimensionality of reference/query points
   * \param glob_ids An array of length nlocal containing the unique ID of each reference point stored locally.
   * \param rneighbors [out] An externally allocated array of mlocal vectors, used to store the neighbors of each
   * query point.
   * \param max_neighbors The maximum number of neighbors to return for any single query point.
   * \param comm The MPI communicator to use.
   */
   void dist_directRQueryIndividual
   ( double* &ref, double *query, long n, long m, int nlocal, int mlocal, double *R,
    int dim, long* &glob_ids, std::vector< std::pair<double, long> >* rneighbors, int max_neighbors,
    MPI_Comm comm  );

   /**
    * Merges a and b and keeps only the k smallest from each of the m sets of distance/index pairs.
    * \param a one array of m*k distance+index pairs
    * \param b the other array of m*k distance+index pairs
    * \param m the number sets of size k
    * \param k k
    * \return A newly allocated array of size m*k
    */
   std::pair<double, long> *kmin_merge( std::pair<double, long> *a, std::pair<double, long> *b, long m, long k );

   /**
    * Reads reference and query points from files in parallel and partitions them for a query with the rectangular algorithm..
    * \param refFile The name of the file which contains the reference points.
    * \param queryFile The name of the file which contains the query points.
    * \param ref [out] an n*dim array of data points.
    * \param query [out] an m*dim array of query points.
    * \param n number of reference points.
    * \param m number of query points.
    * \param dim dimensionality of reference/query points.
    * \param refParts number of parts to split reference points into (refParts * queryParts = processes).
    * \param queryParts number of parts to split query points into (refParts * queryParts = processes).
    * \param comm The MPI communicator to use.
    * \return A pair object containing the number of reference points (first) and query points (second) stored at this process.
    */
   std::pair<int, int> partitionPoints(std::string refFile, std::string queryFile, double *& ref, double *& query, long n, long m, int dim, int refParts, int queryParts, MPI_Comm comm);

   /**
    * Returns the number of local elements in the partitioning used by dist_directKQuery.
    * \param rank This process's MPI rank.
    * \param size The size of the MPI communicator used.
    * \param num The global number of elements.
    * \return The number of elements this process stores locally.
    */
   inline int getNumLocal( int rank, int size, long num ) {
      if( rank < num%size)
         return std::ceil( (double)num / (double)size );
      else
         return num / size;
   }

   /**
    * Returns the offset of the locally stored elements in the global array.
    * \param rank This process's MPI rank.
    * \param size The size of the MPI communicator used.
    * \param num The global number of elements.
    * \return The global array offset.
    */
   inline long getGlobalArrayOffset( int rank, int size, long num ) {
      if( rank < num%size )
         return std::ceil( (double)num / (double)size ) * rank;
      else
         return std::ceil( (double)num / (double)size ) * (num%size) + (rank-(num%size)) * (num/size);
   }

   /**
    * Repartitions points with replication.
    * \param points The points to repartition (modified).
    * \param localPointCount The number of points the process has.
    * \param dim The dimensionality of the points.
    * \param globalIDs The global ID associated with each point (modified).
    * \param partitionCount The number of pieces to partition the points into (the number of replications is implicitly given by size / partitionCount)
    * \param comm The communicator to use.
    * \return The new number of points the process has.
   */
   int directRectRepartition(double *&points, int localPointCount, int dim, int *&globalIDs, int partitionCount, int direction, MPI_Comm comm);

   /**
    * Performs a direct kNN or rNN query using the rectangular-partitioned algorithm, automatically repartitioning points..
    * \param refPts The reference points (modified).
    * \param queryPts The query points (modified).
    * \param localRefCount The number of reference points the process has (modified).
    * \param localQueryCount The number of query points the process has (modified).
    * \param dim The dimensionality of the points.
    * \param refIDs The global IDS of the reference points the process has (modified).
    * \param queryIDs The global IDs of the query points the process has (modified).
    * \param params The parameters of the direct kNN or rNN query.
    * \param comm The communicator to use.
    * \return The result of the query.
   */
   directQueryResults directRectRepartitionAndQuery(double *&refPts, double *&queryPts, long& localRefCount, long& localQueryCount, int dim, int *&refIDs, int *&queryIDs, directQueryParams params, MPI_Comm comm);

   /**
    * Performs a rNN query on a rectangularly partitioned set of points.
    * \param refPoints The reference points.
    * \param queryPoints The query points.
    * \param localRefCount The number of reference points the process has.
    * \param localQueryCount The number of query points the process has.
    * \param dim The dimensionality of the points.
    * \param r The search radius of the query.
    * \param queryComm The communicator to use (all processes within this comm should have the same query points).
    * \return The result of the query.
   */
   directQueryResults directQueryRectR(double *refPoints, double *queryPoints, int localRefCount, int localQueryCount, int dim, double r, MPI_Comm queryComm);

   /**
    * Performs a kNN query on a rectangularly partitioned set of points.
    * \param refPoints The reference points.
    * \param queryPoints The query points.
    * \param localRefCount The number of reference points the process has.
    * \param localQueryCount The number of query points the process has.
    * \param dim The dimensionality of the points.
    * \param k The number of neighbors to find for each query point.
    * \param queryComm The communicator to use (all processes within this comm should have the same query points).
    * \return The result of the query.
   */
   directQueryResults directQueryRectK(double *refPoints, double *queryPoints, int localRefCount, int localQueryCount, int dim, int k, MPI_Comm queryComm);

   /**
    * Determine the number of query points that will be handled by each iteration of the LowMem queries.
    * \param n The number of reference points.
    * \param m The number of query points.
    * \return The block size the query will use.  Can be multiplied by n to determine size of distance matrix.
    */
   inline int getBlockSize(long n, long m) {
     int blocksize;
     if( m > KNN_MAX_BLOCK_SIZE || n > 10000L) {
       blocksize = std::min((long)KNN_MAX_BLOCK_SIZE, m); //number of query points handled in a given iteration
       if(n * (long)blocksize > KNN_MAX_MATRIX_SIZE) blocksize = std::min((long)(KNN_MAX_MATRIX_SIZE/n), (long)blocksize); //Shrink block size if n is huge.
       blocksize = std::max(blocksize, omp_get_max_threads()); //Make sure each thread has some work.
     } else {
        blocksize = m;
     }
     return blocksize;
   }

}
#endif

max_stars_repo_name: maumueller/rehashing

max_issues_repo_name: maumueller/rehashing

max_forks_repo_name: maumueller/rehashing


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

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

#include "ccl.h"

/*----- ROUTINE: dc_NakamuraSuto -----
INPUT: cosmology, scale factor
TASK: Computes the peak threshold: delta_c(z) assuming LCDM.
Cosmology dependence of the critical linear density according to the spherical-collapse model.
Fitting function from Nakamura & Suto (1997; arXiv:astro-ph/9710107).
*/
double dc_NakamuraSuto(ccl_cosmology *cosmo, double a, int *status){

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double dc0 = (3./20.)*pow(12.*M_PI,2./3.);
  double dc = dc0*(1.+0.012299*log10(Om_mz));

  return dc;

}

/*----- ROUTINE: Dv_BryanNorman -----
INPUT: cosmology, scale factor
TASK: Computes the virial collapse density contrast with respect to the matter density assuming LCDM.
Cosmology dependence of the virial collapse density according to the spherical-collapse model
Fitting function from Bryan & Norman (1998; arXiv:astro-ph/9710107)
*/
double Dv_BryanNorman(ccl_cosmology *cosmo, double a, int *status){

  double Om_mz = ccl_omega_x(cosmo, a, ccl_species_m_label, status);
  double x = Om_mz-1.;
  double Dv0 = 18.*pow(M_PI,2);
  double Dv = (Dv0+82.*x-39.*pow(x,2))/Om_mz;

  return Dv;
}

static double sigmaM_m2r(ccl_cosmology *cosmo, double halomass, int *status)
{
  double rho_m, smooth_radius;

  // Comoving matter density
  //rho_m = ccl_constants.RHO_CRITICAL*cosmo->params.Omega_m*cosmo->params.h*cosmo->params.h;
  rho_m = ccl_rho_x(cosmo, 1., ccl_species_m_label, 1, status);

  smooth_radius = pow((3.0*halomass) / (4*M_PI*rho_m), (1.0/3.0));

  return smooth_radius;
}

void ccl_cosmology_compute_sigma(ccl_cosmology *cosmo, int *status)
{
  if(cosmo->computed_sigma)
    return;

  // create linearly-spaced values of the mass.
  int nm = cosmo->spline_params.LOGM_SPLINE_NM;
  double *m = NULL;
  double *y = NULL;
  double smooth_radius;
  double na, nb;

  m = ccl_linear_spacing(cosmo->spline_params.LOGM_SPLINE_MIN, cosmo->spline_params.LOGM_SPLINE_MAX, nm);
  if (m == NULL ||
      (fabs(m[0]-cosmo->spline_params.LOGM_SPLINE_MIN)>1e-5) ||
      (fabs(m[nm-1]-cosmo->spline_params.LOGM_SPLINE_MAX)>1e-5) ||
      (m[nm-1]>10E17)) {
    *status = CCL_ERROR_MEMORY;
    ccl_cosmology_set_status_message(cosmo,"ccl_cosmology_compute_sigmas(): Error creating linear spacing in m\n");
  }

  if (*status == 0) {
    // create space for y, to be filled with sigma and dlnsigma_dlogm
    y = malloc(sizeof(double)*nm);
    if (y == NULL) {
      *status = CCL_ERROR_MEMORY;
    }
  }

  // start up of GSL pointers
  int gslstatus = 0;
  gsl_spline *logsigma = NULL;
  gsl_spline *dlnsigma_dlogm = NULL;

  // fill in sigma, if no errors have been triggered at this time.
  if (*status == 0) {
    for (int i=0; i<nm; i++) {
      smooth_radius = sigmaM_m2r(cosmo, pow(10,m[i]), status);
      y[i] = log(ccl_sigmaR(cosmo, smooth_radius, 1., status));
    }
    logsigma = gsl_spline_alloc(cosmo->spline_params.M_SPLINE_TYPE, nm);
    if (logsigma == NULL) {
      *status = CCL_ERROR_MEMORY;
    }
  }

  if (*status == 0) {
    gslstatus = gsl_spline_init(logsigma, m, y, nm);
    if (gslstatus != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error creating sigma(M) spline\n");
    }
  }

  // again, making splines assuming nothing bad has happened to this point
  if (*status == 0 ) {
    for (int i=0; i<nm; i++) {
      if(i==0) {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i], NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]+cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&nb);
        y[i] = 2.*(na-nb)*y[i] / cosmo->spline_params.LOGM_SPLINE_DELTA;
      }
      else if (i==nm-1) {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]-cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i], NULL,&nb);
        y[i] = 2.*(na-nb)*y[i] / cosmo->spline_params.LOGM_SPLINE_DELTA;
      }
      else {
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]-cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&na);
        gslstatus |= gsl_spline_eval_e(logsigma, m[i]+cosmo->spline_params.LOGM_SPLINE_DELTA/2., NULL,&nb);
        y[i] = (na-nb) / cosmo->spline_params.LOGM_SPLINE_DELTA;
      }
    }

    if(gslstatus != GSL_SUCCESS ) {
      ccl_raise_gsl_warning(
        gslstatus, "ccl_massfunc.c: ccl_cosmology_compute_sigma():");
      *status = CCL_ERROR_SPLINE;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: ccl_cosmology_compute_sigma(): "
        "Error evaluating grid points for dlnsigma/dlogM spline\n");
    }
  }

  if (*status == 0) {
    dlnsigma_dlogm = gsl_spline_alloc(cosmo->spline_params.M_SPLINE_TYPE, nm);
    if (dlnsigma_dlogm == NULL) {
      *status = CCL_ERROR_MEMORY;
    }
  }

  if (*status == 0) {
    gslstatus = gsl_spline_init(dlnsigma_dlogm, m, y, nm);
    if (gslstatus != GSL_SUCCESS) {
      *status = CCL_ERROR_SPLINE ;
      ccl_cosmology_set_status_message(
        cosmo,
        "ccl_massfunc.c: ccl_cosmology_compute_sigma(): Error creating dlnsigma/dlogM spline\n");
    }
  }

  if (*status == 0) {
    cosmo->data.logsigma = logsigma;
    cosmo->data.dlnsigma_dlogm = dlnsigma_dlogm;
    cosmo->computed_sigma = true;
  } else {
    gsl_spline_free(logsigma);
    gsl_spline_free(dlnsigma_dlogm);
  }

  free(m);
  free(y);
}

/*----- ROUTINE: ccl_sigma_M -----
INPUT: ccl_cosmology * cosmo, double halo mass in units of Msun, double scale factor
TASK: returns sigma from the sigmaM interpolation. Also computes the sigma interpolation if
necessary.
*/
double ccl_sigmaM(ccl_cosmology *cosmo, double log_halomass, double a, int *status)
{
  double sigmaM;
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    *status = CCL_ERROR_SIGMA_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_massfunc.c: ccl_sigmaM(): linear power spctrum has not been computed!");
    return NAN;
  }

  double lgsigmaM;

  int gslstatus = gsl_spline_eval_e(cosmo->data.logsigma, log_halomass, NULL, &lgsigmaM);

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

  // Interpolate to get sigma
  sigmaM = exp(lgsigmaM)*ccl_growth_factor(cosmo, a, status);
  return sigmaM;
}

/*----- ROUTINE: ccl_dlnsigM_dlogM -----
INPUT: ccl_cosmology *cosmo, double halo mass in units of Msun
TASK: returns the value of the derivative of ln(sigma^-1) with respect to log10 in halo mass.
*/
double ccl_dlnsigM_dlogM(ccl_cosmology *cosmo, double log_halomass, int *status)
{
  // Check if sigma has already been calculated
  if (!cosmo->computed_sigma) {
    *status = CCL_ERROR_SIGMA_INIT;
    ccl_cosmology_set_status_message(
      cosmo,
      "ccl_massfunc.c: ccl_sigmaM(): linear power spctrum has not been computed!");
    return NAN;
  }

  double dlsdlgm;
  int gslstatus = gsl_spline_eval_e(cosmo->data.dlnsigma_dlogm,
				    log_halomass, NULL, &dlsdlgm);
  if(gslstatus) {
    ccl_raise_gsl_warning(gslstatus, "ccl_massfunc.c: ccl_dlnsigM_dlogM():");
    *status |= gslstatus;
  }
  return dlsdlgm;
}

max_stars_repo_name: timothydmorton/CCL

max_issues_repo_name: timothydmorton/CCL

max_forks_repo_name: timothydmorton/CCL


text:
#pragma once
#include <string>
#include <type_traits>
#include <tuple>
#include <gsl.h>
#include "roster.h"
#include "timestamp.h"
#include <boost/multi_index/sequenced_index.hpp>

class Grade
{
public:
	class Key
	{
	private:
		Roster const& roster_;
		std::string grader_email_;
		std::string student_name_;

	private:
		friend bool operator==(const Key lhs, const Key rhs)
		{
			const auto result{ lhs.grader_email_ == rhs.grader_email_ && lhs.student_name_ == rhs.student_name_ };
			return result;
		}

	public:
		explicit Key(Roster const& roster);

		std::string getGraderEmail() const
		{
			Expects(!grader_email_.empty() && roster_.is_email_good(grader_email_));
			return grader_email_;
		}

		void setGraderEmail(const std::string grader_email)
		{
			Expects(!grader_email.empty() && roster_.is_email_good(grader_email));
			grader_email_ = grader_email;
		}

		std::string getStudentName() const
		{
			Expects(!student_name_.empty());
			return student_name_;
		}

		void setStudentName(const std::string student_name)
		{
			Expects(!student_name.empty());
			student_name_ = student_name;
		}

        bool is_good() const;
	};

private:
    Roster const& roster_;
	Key key_;
	Timestamp timestamp_{};
	unsigned score_{ 0 };

public:
	bool is_good() const;

	Key getKey() const
	{
		return key_;
	}

	Timestamp getTimestamp() const
	{
		return timestamp_;
	}

	unsigned getScore() const
	{
		return score_;
	}

private:
	void setTimestamp(const Timestamp timestamp)
	{
		Expects(timestamp.is_good());
		timestamp_ = timestamp;
	}

	friend std::istream& operator>>(std::istream& is, Grade& grade)
	{
		Timestamp timestamp;

		is >> timestamp;
		Ensures(timestamp.is_good());
		grade.setTimestamp(timestamp);

		char grader_email_cstr[256]{};
		constexpr auto max_grader_len{ std::extent<decltype(grader_email_cstr)>::value };
		is.getline(grader_email_cstr, max_grader_len, ',');
		auto const grader_len{ strlen(grader_email_cstr) };
		Ensures(grader_len > 0 && grader_len < max_grader_len);

		grade.key_.setGraderEmail(grader_email_cstr);

		is.ignore(std::numeric_limits<std::streamsize>::max(), ',');

		char student_name_cstr[256]{};
		constexpr auto max_student_name_len{ std::extent<decltype(student_name_cstr)>::value };
		is.getline(student_name_cstr, max_student_name_len, ',');
		auto const student_len{ strlen(student_name_cstr) };
		Ensures(student_len > 0 && student_len < max_student_name_len);

		grade.key_.setStudentName(student_name_cstr);

		Ensures(grade.is_good());

		is.seekg(-2, std::ios_base::end);

		auto const is_comma{ is.peek()==',' };

		if (is_comma)
		{
			is.ignore();
		}

		auto grade_value{ std::numeric_limits<int>::max() };

		is >> grade_value;

		Ensures(grade_value >= 0 && grade_value<(10+1) && is.eof());

		grade.score_ = 10 * grade_value;

		Ensures(grade.score_ >= 0 && grade.score_ < (100+1));

		return is;
	}

public:
	explicit Grade(Roster const& roster);
};

max_stars_repo_name: labermt/SRSGrade

max_issues_repo_name: labermt/SRSGrade

max_forks_repo_name: labermt/SRSGrade


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "parmt_utils.h"
#ifdef PARMT_USE_INTEL
#include <mkl_lapacke.h>
#include <mkl_cblas.h>
#else
#include <lapacke.h>
#include <cblas.h>
#endif
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"
#include "iscl/random/random.h"
#include "iscl/signal/convolve.h"

int parmt_utils_getEffectiveRank(const double relError,
                                 const struct parmtNoiseBasis_struct basis)
{
    double frac;
    int i, rank;
    rank = 1;
    for (i=1; i<basis.nvals; i++)
    {
        // Modification of Demmel 3.9.3 to deal with sqrt(eig)
        if (100.0*pow(basis.sqrtEvals[i], 2)/pow(basis.sqrtEvals[0], 2) < relError)
        {
            rank = i;
            break;
        }
    }
    return rank;
}

int parmt_utils_makeKLNoise(const int rank,
                            const int npts,
                            const struct parmtNoiseBasis_struct basis,
                            double *__restrict__ xn)
{
    double *r, xscal;
    int ir, ierr, ldz, rankUse;
    if (rank < 1 || npts < 1 || basis.npts != npts || xn == NULL)
    {
        if (rank < 1){fprintf(stderr, "%s: rank is too small\n", __func__);}
        if (npts < 1){fprintf(stderr, "%s: no points\n", __func__);}
        if (basis.npts != npts)
        {
            fprintf(stderr, "%s: npts != basis.npts\n", __func__);
        }
        if (xn == NULL){fprintf(stderr, "%s: error x is NULL\n", __func__);}
        return -1;
    }
    rankUse = MAX(1, MIN(basis.nvals, rank));
    ldz = basis.lde;
    r = random_rand64f(rankUse, &ierr);
    array_zeros64f_work(npts, xn);
    for (ir=0; ir<rank; ir++)
    {
        xscal = r[ir]*basis.sqrtEvals[ir];
printf("%f %f\n", r[ir], basis.sqrtEvals[ir]);
        cblas_daxpy(npts, xscal, &basis.evects[ir*ldz], 1, xn, 1);
    }
    memory_free64f(&r);
    return 0;
}
//============================================================================//
/*!
 * @brief Releases memory on the noise basis struture
 *
 * @param[out] basis   on exit all memory has been released and all variables
 *                     set to NULL or 0.
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_utils_freeNoiseBasis(struct parmtNoiseBasis_struct *basis)
{
    memory_free64f(&basis->sqrtEvals);
    memory_free64f(&basis->evects);
    memset(basis, 0, sizeof(struct parmtNoiseBasis_struct));
    return 0;
}
//============================================================================//
/*!
 * @brief Computes the eigenvectors and eigenvalues of the noise
 *        autocorrelation matrix which will be further used in the
 *        Karhunen-Loeve expansion.
 *
 * @param[in] npts     number of data points
 * @param[in] data     noise signal from which to compute noise basis.
 *
 * @param[out] basis   on successful exit contains the eigendecomposition
 *                     of the autocorrelation noise matrix for use in the KL
 *                     expansion.
 *
 * @result 0 indicates success
 *
 * @author Ben Baker
 *
 * @copyright ISTI distributed under Apache 2
 *
 */
int parmt_utils_getNoiseBasis64f(const int npts,
                                 const double *__restrict__ data,
                                 struct parmtNoiseBasis_struct *basis)
{
    double *C, *CtC, *s, xnorm;
    int i, ierr, j, ldc, ldctc, m, n, nrows;
    const enum corrMatrixType_enum type = CORRMTX_AUTOCORRELATION;
    //------------------------------------------------------------------------//
    //
    // error checking
    ierr = 0;
    memset(basis, 0, sizeof(struct parmtNoiseBasis_struct));
    if (npts < 1 || data == NULL)
    {
        if (npts < 1)
        {
            fprintf(stderr, "%s: No points in noise signal\n", __func__);
        }
        if (data == NULL)
        {
            fprintf(stderr, "%s; Noise signal is NULL\n", __func__);
        }
        return -1;
    }
    C = NULL;
    CtC = NULL;
    // Normalize the noise energy in the noise signal
    xnorm = cblas_dnrm2(npts, data, 1);
    if (fabs(xnorm) < 1.e-15)
    {
        fprintf(stderr, "%s: Error division by zero\n", __func__);
        return -1;
    }
    xnorm = 1.0/xnorm;
    s = array_copy64f(npts, data, &ierr);
    cblas_dscal(npts, xnorm, s, 1);
    // Figure out the size of the noise correlation matrix
    n = npts;
    m = npts - 1;
    ldc = n;
    if (type == CORRMTX_AUTOCORRELATION)
    {
        nrows = n + m;
        C = memory_calloc64f(ldc*nrows);
    }
    else
    {
        fprintf(stderr, "%s: Invalid correlation matrix type\n", __func__);
        return -1;
    }
    ldctc = n;
    CtC = memory_calloc64f(n*n);
    // Compute the correlation matrix and C'C
    ierr = convolve_corrmtx64f_work(npts, s,
                                    m, ldc, C,
                                    type, true, ldctc, CtC);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to compute correlation matrix\n", __func__);
        goto ERROR;
    }
    ierr = parmt_utils_getNoiseBasisFromCtC64f(npts, ldctc, CtC, basis);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to compute noise basis\n", __func__);
    }
ERROR:;
    memory_free64f(&C);
    memory_free64f(&CtC);
    return ierr;
}

int parmt_utils_getNoiseBasisFromCtC64f(const int npts, const int ldctc,
                                        const double *__restrict__ CtC,
                                        struct parmtNoiseBasis_struct *basis)
{
    double *CtCw, *w;
    int ierr, i, info, j;
    ierr = 0;
    memset(basis, 0, sizeof(struct parmtNoiseBasis_struct));
    CtCw = memory_calloc64f(npts*ldctc);
    for (i=0; i<npts; i++)
    {
        cblas_dcopy(npts, &CtC[i*ldctc], 1, &CtCw[i*ldctc], 1);
    }
    // Compute the eigendecompsition of the symmetric matrix
    w = memory_calloc64f(npts);
    info = LAPACKE_dsyev(LAPACK_ROW_MAJOR, 'V', 'U', npts, CtCw, ldctc, w);
    if (info != 0)
    {
        ierr = 1;
        if (info > 0)
        {
            fprintf(stderr, "%s: Failure to converge on %d'th element\n",
                    __func__, info);
        }
        else
        {
            fprintf(stderr, "%s: %d'th parameter is invalid\n", __func__, info);
        }
        goto ERROR;
    }
    // Set the basis
    basis->lde = npts;
    basis->nvals = npts;
    basis->npts = npts;
    basis->sqrtEvals = memory_calloc64f(basis->nvals);
    basis->evects = memory_calloc64f(basis->lde*basis->nvals);
    // Put into descending order
    for (i=0; i<npts; i++)
    {
        basis->sqrtEvals[i] = sqrt(w[npts-1-i]);
    }
    // Copy corresponding eigenvectors in descending order
    for (j=0; j<npts; j++)
    {
        for (i=0; i<npts; i++)
        {
            basis->evects[j*basis->lde+i] = CtCw[i*ldctc+(npts-1-j)];
        }
    }
ERROR:;
    memory_free64f(&w);
    memory_free64f(&CtCw);
    return ierr;
}