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 <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>

#include "../cosmolike_core_fast/theory/basics.c"
#include "../cosmolike_core_fast/theory/structs.c"
#include "../cosmolike_core_fast/theory/parameters.c"
#include "../cosmolike_core_fast/emu17/P_cb/emu.c"
#include "../cosmolike_core_fast/theory/recompute.c"
#include "../cosmolike_core_fast/theory/cosmo3D.c"
#include "../cosmolike_core_fast/theory/redshift_spline.c"
#include "../cosmolike_core_fast/theory/halo.c"
#include "../cosmolike_core_fast/theory/HOD.c"
#include "../cosmolike_core_fast/theory/pt.c"
#include "../cosmolike_core_fast/theory/cosmo2D_fourier.c"
#include "../cosmolike_core_fast/theory/IA.c"
// #include "../cosmolike_core_fast/theory/cluster.c"
#include "../cosmolike_core_fast/theory/BAO.c"
#include "../cosmolike_core_fast/theory/external_prior.c"
#include "../cosmolike_core_fast/theory/init_baryon.c"
#include "init_LSSxCMB.c"

// Naming convention:
// g = galaxy positions ("g" as in "galaxy")
// k = kappa CMB ("k" as in "kappa")
// s = kappa from source galaxies ("s" as in "shear")
// And alphabetical order

typedef double (*C_tomo_pointer)(double l, int n1, int n2);
void twopoint_via_hankel(double **xi, double *logthetamin, double *logthetamax, C_tomo_pointer C_tomo, int ni, int nj, int N_Bessel);

#include "../cosmolike_core_fast/theory/CMBxLSS_fourier.c"

typedef struct input_cosmo_params_local {
    double omega_m;
    double sigma_8;
    double n_s;
    double w0;
    double wa;
    double omega_b;
    double h0;
    double MGSigma;
    double MGmu;
} input_cosmo_params_local;

typedef struct input_nuisance_params_local {
    double bias[10];
    double source_z_bias[10];
    double source_z_s;
    double lens_z_bias[10];
    double lens_z_s;
    double shear_m[10];
    double A_ia;
    double beta_ia;
    double eta_ia;
    double eta_ia_highz;
    double lf[6];
    double m_lambda[6];
    double bary[3];
} input_nuisance_params_local;

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);
double C_ks_sys(double ell, 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_gk(double *ell, double *data, int start);
void set_data_ks(double *ell, double *data, int start);
void set_data_kk(double *ell, 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 Q1, double Q2, double Q3);
// 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 Q1, double Q2, double Q3);
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 Q1, double Q2, double Q3);
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 Q1, double Q2, double Q3);

void write_datavector_wrapper(char *details, input_cosmo_params_local ic, input_nuisance_params_local in);
double log_like_wrapper(input_cosmo_params_local ic, input_nuisance_params_local 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);
  // clock_t t1, t2; t1=clock();
  // 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==1) C = C_shear_tomo_nointerp_IA1(ell,z1,z2);
  // t2 = clock(); printf("shear %d-%d: %le\n",z1,z2, (double)(t2-t1)/CLOCKS_PER_SEC);
  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;
// }

double C_ks_sys(double ell, int zs)
{
   double C;
   C = C_ks_nointerp(ell,zs);
   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);
#ifdef ONESAMPLE
    if(zl==0){
      for (i = 0; i < Ncl; i++){data[start+(Ncl*nz)+i] = 0.;}
    }
    else{
#endif
      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.;
        }
      }
#ifdef ONESAMPLE
    }
#endif
  }
}

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_gk(double *ell, double *data, int start)
{
#ifdef ONESAMPLE
  int nz=0;
  for (int i=0; i<like.Ncl; i++) {data[start+(like.Ncl*nz)+i] = 0.;}
  for (nz=1; nz<tomo.clustering_Nbin; nz++){
#else
  for (int nz=0; nz<tomo.clustering_Nbin; nz++){
#endif
    for (int i=0; i<like.Ncl; i++){
       if (ell[i]<like.lmax_kappacmb && test_kmax(ell[i],nz)){
          data[start+(like.Ncl*nz)+i] = C_gk_nointerp(ell[i],nz);
       }
       else{
          data[start+(like.Ncl*nz)+i] = 0.;
       }
    }
  }
}

void set_data_ks(double *ell, double *data, int start)
{
   for (int nz=0; nz<tomo.shear_Nbin; nz++){
      for (int i=0; i<like.Ncl; i++){
         if (ell[i]<like.lmax_kappacmb) {
            data[start+(like.Ncl*nz)+i] = C_ks_sys(ell[i],nz);
         }
         else{
            data[start+(like.Ncl*nz)+i] = 0.;
         }
      }
   }
}

void set_data_kk(double *ell, double *data, int start)
{
   for (int i=0; i<like.Ncl; i++){
      if (ell[i]<like.lmax_kappacmb){
         data[start+i] = C_kk_nointerp(ell[i]);
      }
      else{
         data[start+i] = 0.;
      }
   }
}

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.1 || 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.85 || cosmology.n_spec > 1.05) return 0;
  if (cosmology.w0 < -2.0 || cosmology.w0 > -0.0) return 0;
  if (cosmology.wa < -2.5 || cosmology.wa > 2.5) return 0;
  if (cosmology.h0 < 0.4 || cosmology.h0 > 0.9) return 0;
  if (cosmology.MGmu < -3.0 || cosmology.MGmu > 3.0) return 0;
  if (cosmology.MGSigma < -3.0 || cosmology.MGSigma > 3.0) return 0; // DESY1 extension paper flat priors
  //CH BEGINS
  //CH: to use for running planck15_BA0_w0_wa prior alone)
  //printf("like_fourier.c from WFIRST_forecasts: cosmology bounds set for running with planck15_BA0_w0_wa prior\n");
  //if (cosmology.Omega_m < 0.05 || cosmology.Omega_m > 0.6) return 0;
  //if (cosmology.omb < 0.01 || cosmology.omb > 0.1) 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 > 1.5) return 0;
  //if (cosmology.wa < -5.0 || cosmology.wa > 2.6) return 0;
  //if (cosmology.h0 < 0.3 || cosmology.h0 > 0.9) return 0;
  //CH ENDS
  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.0001) 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.0001) 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 < -4.0 || nuisance.beta_ia > 6.0) return 0;
  if (nuisance.eta_ia < -10.0 || nuisance.eta_ia> 10.0) return 0;
  if (nuisance.eta_ia_highz < -1.0 || nuisance.eta_ia_highz> 1.0) return 0;
  // 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;
  for (i = 0; i < 10; i++){
  //    printf("in set routine %d %le\n",i,gbias.b[i]);
#ifdef ONESAMPLE
    if (gbias.b[i] < 0.4 || gbias.b[i] > 5.0) return 0;
#else
    if (gbias.b[i] < 0.4 || gbias.b[i] > 3.0) return 0;
#endif
  }

  return 1;
}

// 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 Q1, double Q2, double Q3)
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 Q1, double Q2, double Q3)
{
  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=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.0e15;
  }
  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.0e15;
  }
  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.0e15;
  }
  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.0e15;
  }
  if (set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10)==0){
    printf("Bias out of bounds\n");
    return -1.0e15;
  }
  // 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.0e15;
  // }

  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.Planck15_BAO_w0wa==1) log_L_prior+=log_L_Planck15_BAO_w0wa();//CH
  //if(like.Planck15_BAO_H070p6_JLA_w0wa==1) log_L_prior+=log_L_Planck15_BAO_H070p6_JLA_w0wa();//CH
  // 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.IA!=0) {
    log_L = 0.0;
    log_L -= pow((nuisance.A_ia - prior.A_ia[0])/prior.A_ia[1],2.0);
    log_L -= pow((nuisance.beta_ia - prior.beta_ia[0])/prior.beta_ia[1],2.0);
    log_L -= pow((nuisance.eta_ia - prior.eta_ia[0])/prior.eta_ia[1],2.0);
    log_L -= pow((nuisance.eta_ia_highz - prior.eta_ia_highz[0])/prior.eta_ia_highz[1],2.0);
    log_L_prior+=0.5*log_L;
  }
  if(like.baryons==1){;
    log_L = 0.0;
    log_L -= pow((Q1 - prior.bary_Q1[0])/prior.bary_Q1[1],2.0);
    log_L -= pow((Q2 - prior.bary_Q2[0])/prior.bary_Q2[1],2.0);
    log_L -= pow((Q3 - prior.bary_Q3[0])/prior.bary_Q3[1],2.0);
    log_L_prior+=0.5*log_L;
  }

  // if(like.clusterMobs==1) log_L_prior+=log_L_clusterMobs();

  // 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;
  // clock_t t1, t2;
  // t1 = clock();
  if(like.shear_shear==1) {
    set_data_shear(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.shear_Npowerspectra;
  }
  // t2 = clock(); printf("shear: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.shear_pos==1){
    set_data_ggl(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.ggl_Npowerspectra;
  }
  // t2 = clock(); printf("ggl: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.pos_pos==1){
    set_data_clustering(like.Ncl,ell,pred, start);
    start=start+like.Ncl*tomo.clustering_Npowerspectra;
  }
  // t2 = clock(); printf("gg: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.gk==1) {
    set_data_gk(ell, pred, start);
    start += like.Ncl*tomo.clustering_Nbin;
  }
  // t2 = clock(); printf("gk: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.ks==1) {
    set_data_ks(ell, pred, start);
    start += like.Ncl*tomo.shear_Nbin;
  }
  // t2 = clock(); printf("ks: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.kk==1) {
    set_data_kk(ell, pred, start);
    start += like.Ncl;
  }
  // t2 = clock(); printf("kk: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  chisqr=0.0;
  for (i=0; i<like.Ndata; i++){
    for (j=0; j<like.Ndata; j++){
      a=(pred[i]-data_read(1,i)+Q1*bary_read(1,0,i)+Q2*bary_read(1,1,i)+Q3*bary_read(1,2,i))*invcov_read(1,i,j)*(pred[j]-data_read(1,j)+Q1*bary_read(1,0,j)+Q2*bary_read(1,1,j)+Q3*bary_read(1,2,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)) < 1.e-25){
    //    printf("%d %le %le %le\n",i,data_read(1,i),pred[i],invcov_read(1,i,i));
    // }
  }
  // t2 = clock(); printf("chisqr: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if (chisqr<0.0){
    printf("error: chisqr = %le\n",chisqr);
    //exit(EXIT_FAILURE);
  }
  printf("%le\n",chisqr);
  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, double Q1, double Q2, double Q3)
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 Q1, double Q2, double Q3)
{

  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]);
// }
  // clock_t t1, t2;

  // t1 = clock();
  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);

  // t2 = clock(); printf("setting: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  int start=0;
  if(like.shear_shear==1) {
    set_data_shear(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.shear_Npowerspectra;
  }
  // t2 = clock(); printf("shear: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.shear_pos==1){
    //printf("ggl\n");
    set_data_ggl(like.Ncl, ell, pred, start);
    start=start+like.Ncl*tomo.ggl_Npowerspectra;
  }
  // t2 = clock(); printf("ggl: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.pos_pos==1){
    //printf("clustering\n");
    set_data_clustering(like.Ncl,ell,pred, start);
    start=start+like.Ncl*tomo.clustering_Npowerspectra;
  }
  // t2 = clock(); printf("gg: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;

  if(like.gk==1) {
    printf("Computing data vector: gk\n");
    set_data_gk(ell, pred, start);
    start += like.Ncl * tomo.clustering_Nbin;
  }
  // t2 = clock(); printf("gk: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if(like.ks==1) {
    printf("Computing data vector: ks\n");
    set_data_ks(ell, pred, start);
    start += like.Ncl * tomo.shear_Nbin;
  }
  // t2 = clock(); printf("ks: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;
  if (like.kk) {
    printf("Computing data vector: kk\n");
    set_data_kk(ell, pred, start);
    start += like.Ncl;
  }
  // t2 = clock(); printf("kk: %le\n", (double)(t2-t1)/CLOCKS_PER_SEC); t1 = t2;

  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);
  // printf("&gbias.b1_function %p\n",&gbias.b1_function);
  // printf("gbias.b1_function  %p\n",gbias.b1_function);
  // printf("bgal_z   %p\n",bgal_z);
  // printf("&bgal_z  %p\n",&bgal_z);
  // printf("b1_per_bin   %p\n",b1_per_bin);
  // printf("&b1_per_bin  %p\n",&b1_per_bin);

}

void write_datavector_wrapper(char *details, input_cosmo_params_local ic, input_nuisance_params_local 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],in.bary[0], in.bary[1], in.bary[2]);
  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.bary[0], in.bary[1], in.bary[2]);
}

double log_like_wrapper(input_cosmo_params_local ic, input_nuisance_params_local 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],in.bary[0], in.bary[1], in.bary[2]);
  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.bary[0], in.bary[1], in.bary[2]);
  return like;
}

double log_like_wrapper_1sample(input_cosmo_params_local ic, input_nuisance_params_local in) // Only used in sampling, not in datav
{
  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.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.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.bary[0], in.bary[1], in.bary[2]);
  return like;
}

void save_zdistr_sources(int zs){
  double z,dz =(redshift.shear_zdistrpar_zmax-redshift.shear_zdistrpar_zmin)/300.0;
  printf("Printing redshift distribution n(z) for source redshift bin %d\n",zs);

   FILE *F1;
   char filename[300];
   sprintf(filename,"zdistris/zdist_sources_bin%d.txt",zs);
   F1 = fopen(filename,"w");
   for (z =redshift.shear_zdistrpar_zmin; z< redshift.shear_zdistrpar_zmax; z+= dz){
      fprintf(F1,"%e %e\n", z, zdistr_photoz(z,zs));
   }
}

void save_zdistr_lenses(int zl){
   double z,dz =(redshift.clustering_zdistrpar_zmax-redshift.clustering_zdistrpar_zmin)/300.0;
  printf("Printing redshift distribution n(z) and bias b(z) for lens redshift bin %d\n",zl);

   FILE *F1;
   char filename[300];
   sprintf(filename,"zdistris/zdist_lenses_bin%d.txt", zl);
   F1 = fopen(filename,"w");
   for (z =redshift.clustering_zdistrpar_zmin; z< redshift.clustering_zdistrpar_zmax; z+= dz){
      fprintf(F1,"%e %e\n", z, pf_photoz(z,zl));
   }
}

 int main(int argc, char** argv)
{
  int i;
  char arg1[400],arg2[400],arg3[400];
/* here, do your time-consuming job */
  int sce=atoi(argv[1]);

  int N_scenarios=2;

  double sigma_zphot_shear[3]={0.05,0.05};
  double sigma_zphot_clustering[3]={0.03,0.03};

  double area_table[2]={12300.0,16500.0}; // Y1 corresponds to DESC SRD Y1, Y6 corresponds to assuming that we cover the full SO area=0.4*fsky and at a depth of 26.1 which is in a range of reasonable scenarios (see https://github.com/LSSTDESC/ObsStrat/tree/static/static )
  double nsource_table[2]={11.0,23.0};
  double nlens_table[2]={18.0,41.0};

  char survey_designation[2][200]={"LSSTxSO_Y1","LSSTxSO_Y6"};
  char tomo_binning_source[2][200]={"source_std","source_std"};
  // even for lens=src, this setting is valid, because the input lens/src zfiles are generated with the same z-cut:(0.2,1.2)

  char source_zfile[2][400]={"src_LSSTY1ext_2-island","src_LSSTY6ext_2-island"};
#ifdef ONESAMPLE
  char lens_zfile[2][400]={"src_LSSTY1ext_2-island","src_LSSTY6ext_2-island"};
  char tomo_binning_lens[2][200]={"lens=src","lens=src"};
  nlens_table[0] = nsource_table[0];
  nlens_table[1] = nsource_table[1];
  sigma_zphot_clustering[0] = sigma_zphot_shear[0];
  sigma_zphot_clustering[1] = sigma_zphot_shear[1];
#else
  char lens_zfile[2][400]={"lens_LSSTY1ext_2-island", "lens_LSSTY6ext_2-island"};
  char tomo_binning_lens[2][200]={"LSST_gold","LSST_gold"};
#endif
  char cmb_yr[2][100]={"so_Y1","so_Y5"};

  init_cosmo_runmode("halofit");
  init_bary(argv[2]);
  init_binning_fourier(15,20.0,3000.0,3000.0,21.0,10,10);

  init_priors(0.002,sigma_zphot_shear[sce],0.001,0.001,sigma_zphot_clustering[sce],0.001,0.001,3.0,1.2,3.8,2.0,16.0,5.0,0.8);
  init_survey(survey_designation[sce],nsource_table[sce],nlens_table[sce],area_table[sce]);
  sprintf(arg1,"zdistris/%s",source_zfile[sce]);
  sprintf(arg2,"zdistris/%s",lens_zfile[sce]);
  init_galaxies(arg1,arg2,"outlier_sim","outlier_sim",tomo_binning_source[sce],tomo_binning_lens[sce]);
  init_IA("NLA_HF","GAMA");
  init_probes("6x2pt");
  init_cmb(cmb_yr[sce]);

#ifdef ONESAMPLE
  sprintf(arg3,"%s_1sample_%s_outext_2-island",survey_designation[sce],argv[2]);
#else
  sprintf(arg3,"%s_%s_outext_2-island",survey_designation[sce],argv[2]);
#endif
  compute_data_vector(arg3,0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,\
    gbias.b[0],gbias.b[1],gbias.b[2],gbias.b[3],gbias.b[4],\
    gbias.b[5],gbias.b[6],gbias.b[7],gbias.b[8],gbias.b[9],\
    0.0,0.0,0.0,0.0,0.0,\
    0.0,0.0,0.0,0.0,0.0,\
    sigma_zphot_shear[sce],0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,sigma_zphot_clustering[sce],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,0.0,0.0,0.0);

  return 0;
}

max_stars_repo_name: CosmoLike/LSSTxSO

max_issues_repo_name: CosmoLike/LSSTxSO

max_forks_repo_name: CosmoLike/LSSTxSO


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

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

#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_filter.h>

static int filter_impulse(const double scale, const double epsilon, const double t, const gsl_vector * x, const gsl_vector * xmedian,
                          gsl_vector * y, gsl_vector * xsigma, size_t * noutlier, gsl_vector_int * ioutlier);

/*
gsl_filter_impulse_alloc()
  Allocate a workspace for impulse detection filtering.

Inputs: K - number of samples in window; if even, it is rounded up to
            the next odd, to have a symmetric window

Return: pointer to workspace
*/

gsl_filter_impulse_workspace *
gsl_filter_impulse_alloc(const size_t K)
{
  gsl_filter_impulse_workspace *w;

  w = calloc(1, sizeof(gsl_filter_impulse_workspace));
  if (w == 0)
    {
      GSL_ERROR_NULL ("failed to allocate space for workspace", GSL_ENOMEM);
    }

  w->movstat_workspace_p = gsl_movstat_alloc(K);
  if (w->movstat_workspace_p == 0)
    {
      gsl_filter_impulse_free(w);
      return NULL;
    }

  return w;
}

void
gsl_filter_impulse_free(gsl_filter_impulse_workspace * w)
{
  if (w->movstat_workspace_p)
    gsl_movstat_free(w->movstat_workspace_p);

  free(w);
}

/*
gsl_filter_impulse()
  Apply an impulse detection filter to an input vector. The filter output is

y_i = { x_i, |x_i - m_i| <= t * S_i
      { m_i, |x_i - m_i| > t * S_i

where m_i is the median of the window W_i^H and S_i is the scale estimate (MAD, IQR, S_n, Q_n)

Inputs: endtype    - how to handle signal end points
        scale_type - which statistic to use for scale estimate (MAD, IQR, etc)
        t          - number of standard deviations required to identity outliers (>= 0)
        x          - input vector, size n
        y          - (output) filtered vector, size n
        xmedian    - (output) vector of median values of x, size n
                     xmedian_i = median of window centered on x_i
        xsigma     - (output) vector of estimated local standard deviations of x, size n
                     xsigma_i = sigma for i-th window: scale*MAD
        noutlier   - (output) number of outliers detected
        ioutlier   - (output) int array indicating outliers identified, size n; may be NULL
                     ioutlier_i = 1 if outlier detected, 0 if not
        w          - workspace

Notes:
*/

int
gsl_filter_impulse(const gsl_filter_end_t endtype, const gsl_filter_scale_t scale_type, const double t,
                   const gsl_vector * x, gsl_vector * y, gsl_vector * xmedian, gsl_vector * xsigma, size_t * noutlier,
                   gsl_vector_int * ioutlier, gsl_filter_impulse_workspace * w)
{
  const size_t n = x->size;

  if (n != y->size)
    {
      GSL_ERROR("input and output vectors must have same length", GSL_EBADLEN);
    }
  else if (xmedian->size != n)
    {
      GSL_ERROR("xmedian vector must match input size", GSL_EBADLEN);
    }
  else if (xsigma->size != n)
    {
      GSL_ERROR("xsigma vector must match input size", GSL_EBADLEN);
    }
  else if ((ioutlier != NULL) && (ioutlier->size != n))
    {
      GSL_ERROR("ioutlier vector must match input size", GSL_EBADLEN);
    }
  else if (t < 0.0)
    {
      GSL_ERROR("t must be non-negative", GSL_EDOM);
    }
  else
    {
      int status;
      double scale = 1.0;

      switch (scale_type)
        {
          case GSL_FILTER_SCALE_MAD:
            {
              /* compute window medians and MADs */
              gsl_movstat_mad(endtype, x, xmedian, xsigma, w->movstat_workspace_p);

              break;
            }

          case GSL_FILTER_SCALE_IQR:
            {
              /* multiplication factor for IQR to estimate stddev for Gaussian signal */
              scale = 0.741301109252801;

              /* calculate the window medians */
              gsl_movstat_median(endtype, x, xmedian, w->movstat_workspace_p);

              /* calculate window IQRs */
              gsl_movstat_qqr(endtype, x, 0.25, xsigma, w->movstat_workspace_p);

              break;
            }

          case GSL_FILTER_SCALE_SN:
            {
              /* calculate the window medians */
              gsl_movstat_median(endtype, x, xmedian, w->movstat_workspace_p);

              /* calculate window S_n values */
              gsl_movstat_Sn(endtype, x, xsigma, w->movstat_workspace_p);

              break;
            }

          case GSL_FILTER_SCALE_QN:
            {
              /* calculate the window medians */
              gsl_movstat_median(endtype, x, xmedian, w->movstat_workspace_p);

              /* calculate window Q_n values */
              gsl_movstat_Qn(endtype, x, xsigma, w->movstat_workspace_p);

              break;
            }

          default:
            GSL_ERROR("unknown scale type", GSL_EDOM);
            break;
        }

      /* apply impulse detecting filter using previously computed scale estimate */
      status = filter_impulse(scale, 0.0, t, x, xmedian, y, xsigma, noutlier, ioutlier);

      return status;
    }
}

/*
filter_impulse()
  Apply an impulse detection filter to an input vector. The filter output is

y_i = { x_i, |x_i - m_i| <= t * S_i OR S_i < epsilon
      { m_i, |x_i - m_i| > t * S_i

where m_i is the median of the window W_i^H and S_i is the scale estimate (MAD, IQR, etc)

Inputs: scale    - scale factor to multiply xsigma to get unbiased estimate of stddev for Gaussian data
        epsilon  - minimum allowed scale estimate for identifying outliers
        t        - number of standard deviations required to identity outliers (>= 0)
        x        - input vector, size n
        xmedian  - vector of median values of x, size n
                   xmedian_i = median of window centered on x_i
        y        - (output) filtered vector, size n
        xsigma   - (output) vector of estimated local standard deviations of x, size n
                   xsigma_i = S_n for i-th window
        noutlier - (output) number of outliers detected
        ioutlier - (output) int array indicating outliers identified, size n; may be NULL
                   ioutlier_i = 1 if outlier detected, 0 if not

Notes:
1) If S_i = 0 or is very small for a particular sample, then the filter may erroneously flag the
sample as an outlier, since it will act as a standard median filter. To avoid this scenario, the
parameter epsilon specifies the minimum value of S_i which can be used in the filter test. Any
samples for which S_i < epsilon are passed through unchanged.
*/

static int
filter_impulse(const double scale, const double epsilon, const double t, const gsl_vector * x, const gsl_vector * xmedian,
               gsl_vector * y, gsl_vector * xsigma, size_t * noutlier, gsl_vector_int * ioutlier)
{
  const size_t n = x->size;

  if (n != y->size)
    {
      GSL_ERROR("input and output vectors must have same length", GSL_EBADLEN);
    }
  else if (xmedian->size != n)
    {
      GSL_ERROR("xmedian vector must match input size", GSL_EBADLEN);
    }
  else if (xsigma->size != n)
    {
      GSL_ERROR("xsigma vector must match input size", GSL_EBADLEN);
    }
  else if ((ioutlier != NULL) && (ioutlier->size != n))
    {
      GSL_ERROR("ioutlier vector must match input size", GSL_EBADLEN);
    }
  else if (t < 0.0)
    {
      GSL_ERROR("t must be non-negative", GSL_EDOM);
    }
  else
    {
      size_t i;

      *noutlier = 0;

      /* build output vector */
      for (i = 0; i < n; ++i)
        {
          double xi = gsl_vector_get(x, i);
          double xmedi = gsl_vector_get(xmedian, i);
          double absdevi = fabs(xi - xmedi); /* absolute deviation for this sample */
          double *xsigmai = gsl_vector_ptr(xsigma, i);

          /* multiply by scale factor to get estimate of standard deviation */
          *xsigmai *= scale;

          /*
           * If the absolute deviation for this sample is more than t stddevs
           * for this window (and S_i is sufficiently large to avoid scale implosion),
           * set the output value to the window median, otherwise use the original sample
           */
          if ((*xsigmai >= epsilon) && (absdevi > t * (*xsigmai)))
            {
              gsl_vector_set(y, i, xmedi);
              ++(*noutlier);
              if (ioutlier)
                gsl_vector_int_set(ioutlier, i, 1);
            }
          else
            {
              gsl_vector_set(y, i, xi);
              if (ioutlier)
                gsl_vector_int_set(ioutlier, i, 0);
            }
        }

      return GSL_SUCCESS;
    }
}

ext: c

max_stars_repo_name: karanbirsandhu/nu-sense

max_issues_repo_name: karanbirsandhu/nu-sense

max_forks_repo_name: karanbirsandhu/nu-sense


text:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <unistd.h>
#include <math.h>
#include <float.h>
#include <iniparser.h>
#include "tdsearch_gridsearch.h"
#ifdef TDSEARCH_USE_INTEL
#include <mkl_cblas.h>
#else
#include <cblas.h>
#endif
#include "tdsearch_greens.h"
#include "sacio.h"
#include "parmt_utils.h"
#include "iscl/array/array.h"
#include "iscl/memory/memory.h"
#include "iscl/os/os.h"
#include "iscl/signal/convolve.h"

#define NTSTAR 12   /*!< Default number of t*'s in the grid search. */
#define TSTAR0 0.1  /*!< Default first t* in grid search. */
#define TSTAR1 1.2  /*!< Default last t* in grid search. */
#define NDEPTH 12   /*!< Default number of depths in grid search. */
#define DEPTH0 0.0  /*!< First depth in grid search (km). */
#define DEPTH1 3.0  /*!< Last depth in grid search (km). */
#define DOUBLE_PAGE_SIZE sysconf(_SC_PAGE_SIZE)/sizeof(double)
#define PAGE_SIZE sysconf(_SC_PAGE_SIZE)
#define LDG 6

/*!
 * @brief Sets the default grid search parameters.
 *
 * @param[out] tds   On exit contains the default grid search settings such
 *                   as the t* range, depth range, and toggles off the max
 *                   time shift.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_setDefaults(struct tdSearch_struct *tds)
{
    int ierr;
    memset(tds, 0, sizeof(struct tdSearch_struct));
    tds->maxShiftTime = DBL_MAX;
    ierr = tdsearch_gridSearch_defineTstarGrid(NTSTAR, TSTAR0, TSTAR1, tds);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to set t* grid\n", __func__);
        return -1;
    }
    ierr = tdsearch_gridSearch_defineDepthGrid(NDEPTH, DEPTH0, DEPTH1, tds);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to set depth grid\n", __func__);
        return -1;
    }
    return ierr;
}
//============================================================================//
/*!
 * @brief Frees memory on the grid-search structure.
 *
 * @param[out] tds    On exit all memory has been freed and parameters set to
 *                    0 or NULL.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_free(struct tdSearch_struct *tds)
{
    memory_free64f(&tds->G);
    memory_free64f(&tds->Gxc);
    memory_free64f(&tds->tstar);
    memory_free64f(&tds->depths);
    memory_free64f(&tds->synthetic);
    memory_free64f(&tds->xcorr);
    memory_free32i(&tds->lags);
    memory_free32i(&tds->grnsMatrixPtr);
    memory_free32i(&tds->grnsXCMatrixPtr);
    memset(tds, 0, sizeof(struct tdSearch_struct));
    return 0;
}
//============================================================================//
/*!
 * @brief Defines the t* grid in grid search.
 *
 * @param[in] ntstars    Number of t*'s in grid search (> 0).
 * @param[in] tstar0     First t* in grid search (>= 0).
 * @param[in] tstar1     Final t* in grid search (>= tstar0).
 *
 * @param[out] tds       Grid search structure with t* grid.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_defineTstarGrid(const int ntstars, const double tstar0,
                                        const double tstar1,
                                        struct tdSearch_struct *tds)
{
    int ierr;
    tds->ntstar = 0;
    if (ntstars < 1 || tstar0 < 0.0 || tstar0 > tstar1)
    {
        if (ntstars < 1)
        {
            fprintf(stderr, "%s: Error no t*'s %d in search\n", __func__, ntstars);
        }
        if (tstar0 < 0.0)
        {
            fprintf(stderr, "%s: Error tstar0 %f must be non-negative\n",
                    __func__, tstar0);
        }
        if (tstar0 > tstar1)
        {
            fprintf(stderr, "%s: Error tstar0 %f > tstar1 %f\n",
                    __func__, tstar0, tstar1);
        }
        return -1;
    }
    tds->ntstar = ntstars;
    memory_free64f(&tds->tstar);
    tds->tstar = array_linspace64f(tstar0, tstar1, ntstars, &ierr);
    return ierr;
}
//============================================================================//
/*!
 * @brief Defines the depth grid in grid search.
 *
 * @param[in] ndepths    Number of depths in grid search (> 0).
 * @param[in] depth0     First depth (km) in grid search (>= 0).
 * @param[in] depth1     Final depth (km) in grid search (>= depth0).
 *
 * @param[out] tds       Grid search structure with depth grid.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_defineDepthGrid(const int ndepths, const double depth0,
                                        const double depth1,
                                        struct tdSearch_struct *tds)
{
    int ierr;
    tds->ndepth = 0;
    if (ndepths < 1 || depth0 < 0.0 || depth0 > depth1)
    {
        if (ndepths < 1)
        {
            fprintf(stderr, "%s: Error no depthss %d in search\n",
                    __func__, ndepths);
        }
        if (depth0 < 0.0)
        {
            fprintf(stderr, "%s: Error depth0 %f must be non-negative\n",
                    __func__, depth0);
        }
        if (depth0 > depth1)
        {
            fprintf(stderr, "%s: Error depth0 %f > depth1 %f\n",
                    __func__, depth0, depth1);
        }
        return -1;
    }
    tds->ndepth = ndepths;
    memory_free64f(&tds->depths);
    tds->depths = array_linspace64f(depth0, depth1, ndepths, &ierr);
    return ierr;
}
//============================================================================//
/*!
 * @brief Reads the grid search parameters from the ini file.
 *
 * @param[in] iniFile    Name of ini file with parameters.
 *
 * @param[out] tds       On successful exit contains the grid search parameters.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_initializeFromFile(const char *iniFile,
                                           struct tdSearch_struct *tds)
{
    double depth0, depth1, tstar0, tstar1;
    dictionary *ini;
    int ierr, ndepth, ntstar;
    // Set the defaults
    ierr = 0;
    tdsearch_gridSearch_setDefaults(tds);
    // Verify the file exists
    if (!os_path_isfile(iniFile))
    {
        fprintf(stderr, "%s: Error ini file %s doesn't exist\n", __func__, iniFile);
        return -1;
    }
    ini = iniparser_load(iniFile);
    // Unpack the parameter file
    tds->maxShiftTime
      = iniparser_getdouble(ini, "tdSearch:gridSearch:maxShiftTime\0", DBL_MAX);
    if (tds->maxShiftTime < DBL_MAX)
    {
        if (tds->maxShiftTime <= 0.0)
        {
            fprintf(stderr, "%s: Max shift time must be positive\n", __func__);
            ierr = 1;
            goto ERROR;
        }
    }
    // Number of t*'s
    ntstar = iniparser_getint(ini, "tdSearch:gridSearch:ntstar\0", NTSTAR);
    if (ntstar < 1)
    {
        fprintf(stderr, "%s: number of t*'s must be positive\n", __func__);
        ierr = 1;
        goto ERROR;
    }
    // Number of depths
    ndepth = iniparser_getint(ini, "tdSearch:gridSearch:ndepths\0", NDEPTH);
    if (ndepth < 1)
    {
        fprintf(stderr, "%s: number of depths must be positive\n", __func__);
        ierr = 1;
        goto ERROR;
    }
    // Depth grid search range
    depth0 = iniparser_getdouble(ini, "tdSearch:gridSearch:depth0\0", DEPTH0);
    depth1 = iniparser_getdouble(ini, "tdSearch:gridSearch:depth1\0", DEPTH1);
    if (depth0 >= depth1)
    {
        fprintf(stderr, "%s: depth0 >= depth1 is invalid\n", __func__);
        ierr = 1;
        goto ERROR;
    }
    // t* grid search range
    tstar0 = iniparser_getdouble(ini, "tdSearch:gridSearch:tstar0\0", TSTAR0);
    tstar1 = iniparser_getdouble(ini, "tdSearch:gridSearch:tstar1\0", TSTAR1);
    if (tstar0 >= tstar1)
    {
        fprintf(stderr, "%s: tstar0 >= tstar1 is invalid\n", __func__);
        ierr = 1;
        goto ERROR;
    }
    // Create the grid
    ierr = tdsearch_gridSearch_defineTstarGrid(ntstar,
                                               tstar0, tstar1,
                                               tds);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error defining t* grid\n", __func__);
        goto ERROR;
    }
    ierr = tdsearch_gridSearch_defineDepthGrid(ndepth,
                                               depth0, depth1,
                                               tds);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error defining depth grid\n", __func__);
        goto ERROR;
    }
ERROR:;
    iniparser_freedict(ini);
    return ierr;
}
//============================================================================//
/*!
 * @brief Returns the number of depths in the grid search.
 *
 * @param[in] tds     Holds the number of t*'s.
 *
 * @result The number of t*'s in the grid search.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_getNumberOfTstars(const struct tdSearch_struct tds)
{
    return tds.ntstar;
}
//============================================================================//
/*!
 * @brief Returns the number of depths in the grid search.
 *
 * @param[in] tds     Holds the number of depths.
 *
 * @result The number of depths in the grid search.
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_getNumberOfDepths(const struct tdSearch_struct tds)
{
    return tds.ndepth;
}
//============================================================================//
/*!
 * @brief Convenience utility for computing the index corresponding to the
 *        2D row major depth x t* grid.
 *
 * @param[in] it      The t* index [0,ntstar-1].
 * @param[in] id      The depth indx [0,ndepth-1].
 * @param[in] ntstar  The number of t*'s in the grid search.
 *
 * @result The grid search index (= id*ntstar + it).
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_gridToIndex2(const int it, const int id,
                                     const int ntstar)
{
    int idt;
    idt = id*ntstar + it;
    return idt;
}
//============================================================================//
/*!
 * @brief Convenience utility for computing the index corresponding to the
 *        2D row major depth x t* grid.
 *
 * @param[in] it      The t* index [0,ntstar-1].
 * @param[in] id      The depth indx [0,ndepth-1].
 * @param[in] tds     Grid search structure holding the number of ntstar's
 *                    in the grid search.
 *
 * @result The grid search index (= id*tds.ntstar + it).
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_gridToIndex(const int it, const int id,
                                    const struct tdSearch_struct tds)
{
    int idt, ntstar;
    ntstar = tdsearch_gridSearch_getNumberOfTstars(tds);
    idt = tdsearch_gridSearch_gridToIndex2(it, id, ntstar);
    return idt;
}
//============================================================================//
/*!
 * @brief Computes the forward modeling matrix, G, so that estimates, e, can
 *        can be efficiently generated via \f$ G m = e \f$ for any moment
 *        tensor, m.  The columns of G represent cross-correlations of the
 *        Green's functions corresponding to the \f$ m_{ij} \f$ moment tensor
 *        component so that the estimate, e, is a cross-correlogram.
 *
 * @param[in] iobs     Observation index in the range [0,nobs-1].
 * @param[in] data     Structure containing the data.
 * @param[in] grns     Structure containing the Green's functions.
 *
 * @param[in,out] tds  On input contains the initialized tds structure.
 * @param[in,out] tds  On exit the forward modeling matrix G has been set.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdsearch_gridSearch_setForwardModelingMatrices(
    const int iobs,
    const struct tdSearchData_struct data,
    const struct tdSearchGreens_struct grns,
    struct tdSearch_struct *tds)
{
    double *G, *Gxc, *work;
    int *grnsMatrixPtr, *grnsXCMatrixPtr, indices[6],
        i, idep, idt, ierr, indx, it, jndx, kndx,
        lxc, lxcMax, lxc0, ntstar, ndepth, npgrns, npts;
    // Reset pointers
    memory_free64f(&tds->G);
    memory_free64f(&tds->Gxc);
    memory_free32i(&tds->grnsMatrixPtr);
    memory_free32i(&tds->grnsXCMatrixPtr);
    // Check sizes
    ndepth = grns.ndepth;
    ntstar = grns.ntstar;
    if (ndepth != tds->ndepth || ntstar != tds->ntstar)
    {
        if (ndepth != tds->ndepth)
        {
            fprintf(stderr, "%s: Error ndepth %d != tds->ndepth %d\n",
                    __func__, ndepth, tds->ndepth);
        }
        if (ntstar != tds->ntstar)
        {
            fprintf(stderr, "%s: Error ntstar %d != tds->ntstar %d\n",
                    __func__, ntstar, tds->ntstar);
        }
        return -1;
    }
    npts = data.obs[iobs].npts;
    if (npts < 1)
    {
        fprintf(stderr, "%s: Error data has no points\n", __func__);
        return -1;
    }
    tds->dnorm = cblas_dnrm2(npts, data.obs[iobs].data, 1);
    if (fabs(tds->dnorm) < DBL_EPSILON)
    {
        fprintf(stderr, "%s: Error - data is identically zero\n", __func__);
        return -1;
    }
    tds->npts = npts;
    sacio_getFloatHeader(SAC_FLOAT_DELTA, data.obs[iobs].header, &tds->dt);
    // Set the Green's functions forward modeling matrix size
    lxcMax = 0;
    grnsMatrixPtr = memory_calloc32i(ndepth*ntstar+1);
    grnsXCMatrixPtr = memory_calloc32i(ndepth*ntstar+1);
    grnsMatrixPtr[0] = 0;
    grnsXCMatrixPtr[0] = 0;
    for (idep=0; idep<ndepth; idep++)
    {
        for (it=0; it<ntstar; it++)
        {
            // Get the Green's functions indices
            ierr = tdsearch_greens_getGreensFunctionsIndices(iobs, it,
                                                             idep, grns,
                                                             indices);
            if (ierr != 0)
            {
                fprintf(stderr, "%s: Failed to get indices\n", __func__);
                return -1;
            }
            // Compute the cross-correlation length at (dep, t*)
            idt = tdsearch_gridSearch_gridToIndex2(it, idep, ntstar);
            lxc0 = 0;
            for (i=0; i<6; i++)
            {
                indx = indices[i];
                npgrns = grns.grns[indx].npts;
                if (npgrns < 1)
                {
                    fprintf(stderr, "%s: Error - no points in grns fn\n", __func__);
                    return -1;
                }
                lxc = npgrns + npts - 1;
                lxcMax = MAX(lxcMax, lxc);
                if (i == 0)
                {
                    lxc0 = lxc;
                    grnsMatrixPtr[idt+1] = grnsMatrixPtr[idt] + npgrns;
                    grnsXCMatrixPtr[idt+1] = grnsXCMatrixPtr[idt] + lxc;
                }
                if (lxc != lxc0)
                {
                    fprintf(stderr, "%s: Inconsistent column sizes\n", __func__);
                    return -1;
                }
            }
        }
    }
    if (lxcMax < 1)
    {
        fprintf(stderr, "%s: Error - cross-correlations are empty\n", __func__);
        return -1;
    }
    // Set the forward modeling matrices
    tds->nptsGrns = grnsMatrixPtr[ndepth*ntstar];
    G = memory_calloc64f(LDG*grnsMatrixPtr[ndepth*ntstar]);
    Gxc = memory_calloc64f(LDG*grnsXCMatrixPtr[ndepth*ntstar]);
    work = memory_calloc64f(lxcMax);
    for (idep=0; idep<ndepth; idep++)
    {
        for (it=0; it<ntstar; it++)
        {
            idt = tdsearch_gridSearch_gridToIndex2(it, idep, ntstar);
            // Get the Green's functions indices
            ierr = tdsearch_greens_getGreensFunctionsIndices(iobs, it,
                                                             idep, grns,
                                                             indices);
            // Compute the cross-correlation length at (dep, t*)
            idt = tdsearch_gridSearch_gridToIndex2(it, idep, ntstar);
            jndx = LDG*grnsMatrixPtr[idt];
            kndx = LDG*grnsXCMatrixPtr[idt];
            for (i=0; i<6; i++)
            {
                indx = indices[i];
                npgrns = grns.grns[indx].npts;
                signal_convolve_correlate64f_work(npts, data.obs[iobs].data,
                                                  npgrns, grns.grns[indx].data,
                                                  CONVCOR_FULL,
                                                  lxc, work);
                //printf("%e %e\n",cblas_dnrm2(npgrns, grns.grns[indx].data, 1),
                //               cblas_dnrm2(lxc, work, 1));
                cblas_dcopy(npgrns, grns.grns[indx].data, 1, &G[jndx+i],   LDG);
                cblas_dcopy(lxc,    work,                 1, &Gxc[kndx+i], LDG);
            }
            //printf("%d\n", grnsMatrixPtr[idt+1]);
        }
    }
    tds->G = G;
    tds->Gxc = Gxc;
    tds->grnsMatrixPtr = grnsMatrixPtr;
    tds->grnsXCMatrixPtr = grnsXCMatrixPtr;
    memory_free64f(&work);
    return 0;
}
//============================================================================//
/*!
 * @brief Sets the moment tensor on the tdSearch structure.
 *
 * @param[in] m6    moment tensor in NED system packed:
 *                   \f$ \{ m_{xx}, m_{yy}, m_{zz},
 *                         m_{xy}, m_{xz}, m_{yz} \} \f$.
 *
 * @param[out] tds  On successful exit contains the moment tensor.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 * @author Ben Baker, ISTI
 *
 */
int tdSearch_gridSearch_setMomentTensor(const double *__restrict__ m6,
                                        struct tdSearch_struct *tds)
{
    tds->lhaveMT = false;
    array_zeros64f_work(8, tds->mt); // NOTE - this is padded
    if (m6 == NULL)
    {
        fprintf(stderr, "%s: Error moment tensor cannot be NULL\n", __func__);
        return -1;
    }
    array_copy64f_work(6, m6, tds->mt);
    tds->lhaveMT = true;
    return 0;
}
//============================================================================//
/*!
 * @brief Sets the moment tensor elements on the grid search structure.
 *
 * @param[in] m11      m11 moment tensor in Newton meters.
 * @param[in] m22      m22 moment tensor in Newton meters.
 * @param[in] m33      m33 moment tensor in Newton meters.
 * @param[in] m12      m12 moment tensor in Newton meters.
 * @param[in] m13      m13 moment tensor in Newton meters.
 * @param[in] m23      m23 moment tensor in Newton meters.
 * @param[in] basis    The appropriate moment tensor basis.
 * @param[in] basis    CE_USE indicates the moment tensor is in up, south
 *                     east coordinates like what one would obtain from
 *                     the GCMT catalog.
 * @param[in] basis    CE_NED indicates the moment tensor is in north, east,
 *                     down coordinates and is consistent with Programs in
 *                     Seismology.
 *
 * @param[in,out] tds  On input contains the initialized grid search
 *                     structure.
 * @param[in,out] tds  On exit contains the input moment tensor.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 */
int tdSearch_gridSearch_setMomentTensorFromElements(
    const double m11, const double m22, const double m33,
    const double m12, const double m13, const double m23,
    const enum compearthCoordSystem_enum basis,
    struct tdSearch_struct *tds)
{
    double mt[6], mtNED[6];
    int ierr;
    tds->lhaveMT = false;
    mt[0] = m11;
    mt[1] = m22;
    mt[2] = m33;
    mt[3] = m12;
    mt[4] = m13;
    mt[5] = m23;
    ierr = compearth_convertMT(1, basis, CE_NED, mt, mtNED);
    if (ierr != 0)
    {
        printf("%s: Failed to set moment tensor\n", __func__);
        return -1;
    }
    ierr = tdSearch_gridSearch_setMomentTensor(mtNED, tds);
    if (ierr != 0)
    {
        printf("%s: Failed to set moment tensor\n", __func__);
        return -1;
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Sets the moment tensor information from an event structure.
 *
 * @param[in] event    Contains the event on which the moment tensor is defined.
 * @param[in,out] tds  On input this is the initialized grid search structure.
 * @param[in,out] tds  On exit contains the given moment tensor.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 */
int tdSearch_gridSearch_setMomentTensorFromEvent(
    const struct tdSearchEventParms_struct event,
    struct tdSearch_struct *tds)
{
    int ierr;
    ierr = tdSearch_gridSearch_setMomentTensorFromElements(
                    event.m11, event.m22, event.m33,
                    event.m12, event.m13, event.m23,
                    event.basis, tds);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Error setting moment tensor\n", __func__);
        return -1;
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Computes the cross-correlations between the data and synthetics
 *        at each (t*,depth) in the grid search.
 *
 * @param[in,out] tds   On input contains the initialized grid search
 *                      structure with the Green's functions and source
 *                      properties.
 * @param[in,out] tds   On exit contains the result of the cross-correlations
 *                      and lags at each (t*,depth) in the grid search.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 */
int tdSearch_gridSearch_performGridSearch(struct tdSearch_struct *tds)
{
    double *u, *ud, dnorm, unorm;
    int id, idt, indx, iopt, it, jndx, l1, l2, lxc, maxlags,
        ndepth, ndt, nlag, npgrns, nrows, ntstar;
    const int ncols = 6;
    enum isclError_enum ierr;
    tds->itopt =-1;
    tds->idopt =-1;
    ndepth = tds->ndepth;
    ntstar = tds->ntstar;
    ndt = ntstar*ndepth;
    dnorm = tds->dnorm;
    memory_free64f(&tds->xcorr);
    memory_free64f(&tds->synthetic);
    memory_free32i(&tds->lags);
    // Some easy checks
    if (!tds->lhaveMT)
    {
        fprintf(stderr, "%s: Error moment tensor not set\n", __func__);
        return -1;
    }
    if (tds->G == NULL || tds->Gxc == NULL ||
        tds->grnsMatrixPtr == NULL || tds->grnsXCMatrixPtr == NULL)
    {
        fprintf(stderr, "%s: You need to set the forward modeling matrices\n", __func__);
        return -1;
    }
    maxlags =-1;
    if (tds->maxShiftTime < DBL_MAX)
    {
        maxlags = (int) (tds->maxShiftTime/tds->dt + 0.5);
    }
    // Set output space
    tds->xcorr = memory_calloc64f(ndt);
    tds->lags = memory_calloc32i(ndt);
    // Compute the cross-correlation
    nrows = tds->grnsXCMatrixPtr[ndt];
    ud = memory_calloc64f(nrows);
    cblas_dgemv(CblasRowMajor, CblasNoTrans,
                nrows, ncols, 1.0, tds->Gxc, LDG,
                tds->mt, 1, 0.0, ud, 1);
    // Compute the synthetics
    nrows = tds->grnsMatrixPtr[ndt];
    u = memory_calloc64f(nrows);
    cblas_dgemv(CblasRowMajor, CblasNoTrans,
                nrows, ncols, 1.0, tds->G, LDG,
                tds->mt, 1, 0.0, u, 1);
    // Now do the grid search
    for (id=0; id<ndepth; id++)
    {
        for (it=0; it<ntstar; it++)
        {
            idt = tdsearch_gridSearch_gridToIndex(it, id, *tds);
            indx = tds->grnsMatrixPtr[idt];
            jndx = tds->grnsXCMatrixPtr[idt];
            lxc = tds->grnsXCMatrixPtr[idt+1]
                - tds->grnsXCMatrixPtr[idt];
            npgrns = tds->grnsMatrixPtr[idt+1]
                   - tds->grnsMatrixPtr[idt];
            unorm = cblas_dnrm2(npgrns, &u[indx], 1);
            if (maxlags < 0)
            {
                iopt = array_argmax64f(lxc, &ud[jndx], &ierr);
            }
            else
            {
                l1 = MAX(0,  npgrns - maxlags);
                l2 = MIN(lxc-1, npgrns + maxlags);
                nlag = l2 - l1 + 1;
                iopt = l1 + array_argmax64f(nlag, &ud[jndx+l1], &ierr);
 /*
for (int k=0; k<lxc; k++)
{
printf("%d %e\n",k, ud[jndx+k]);
}
getchar();
*/
            }
            tds->xcorr[idt] = ud[jndx+iopt]/(dnorm*unorm);
            // Compute the lag where the times range from [-npts+1:npgrns-1]
            tds->lags[idt] =-npgrns + 1 + iopt;
            //printf("%f %f %f %d\n", tds->depths[id], tds->tstar[it], xc, idelay);
        }
        //printf("\n");
    }
    tds->synthetic = u;
    memory_free64f(&ud);
    // Compute the optimum indices
    iopt = array_argmax64f(ndt, tds->xcorr, &ierr);
    tds->idopt = iopt/tds->ntstar;
    tds->itopt = iopt - tds->idopt*tds->ntstar;
    if (tdsearch_gridSearch_gridToIndex(tds->itopt, tds->idopt, *tds) != iopt)
    {
        fprintf(stdout, "%s: Switching to slow optimum computation\n", __func__);
        for (id=0; id<ndepth; id++)
        {
            for (it=0; it<ntstar; it++)
            {
                idt = tdsearch_gridSearch_gridToIndex(it, id, *tds);
                if (iopt == idt)
                {
                    tds->idopt = id;
                    tds->itopt = it;
                }
            }
        }
    }
    return 0;
}
//============================================================================//
/*!
 * @brief Writes a map of the max of the cross-correlograms as a function
 *        of t* and depth in a format useful for Gnpulot as well as a
 *        corresponding Gnuplot plotting script.
 *
 * @param[in] dirnm   Directory where the heat map is to be written.
 * @param[in] froot   The root name of the heat map.
 * @param[in] tds     Grid search structure with the tabulated
 *                    cross-correlogram maximums.
 *
 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 *
 */
int tdsearch_gridSearch_writeHeatMap(const char *dirnm, const char *froot,
                                     const struct tdSearch_struct tds)
{
    FILE *ofl;
    char fname[PATH_MAX], fgnu[PATH_MAX], cmd[256];
    int id, idt, it;
    memset(fname, 0, PATH_MAX*sizeof(char));
    if (!tds.xcorr || !tds.lags)
    {
        fprintf(stderr, "%s: Error gridsearch likely not performed\n", __func__);
        return -1;
    }
    if (!tds.tstar || !tds.depths)
    {
        fprintf(stderr, "%s: Error tds likely not initialized\n", __func__);
        return -1;
    }
    if (dirnm != NULL)
    {
        sprintf(fname, "%s/%s.surf_txt", dirnm, froot);
    }
    else
    {
        sprintf(fname, "./%s.surf_txt", froot);
    }
    ofl = fopen(fname, "w");
    for (id=0; id<tds.ndepth; id++)
    {
        for (it=0; it<tds.ntstar; it++)
        {
            idt = tdsearch_gridSearch_gridToIndex(it, id, tds);
            fprintf(ofl, "%f %f %e %d\n",
                    tds.tstar[it], tds.depths[id],
                    tds.xcorr[idt], tds.lags[idt]);
        }
       fprintf(ofl, "\n");
    }
    fclose(ofl);
    // and the gnuplot file
    memset(fgnu, 0, PATH_MAX*sizeof(char));
    if (dirnm != NULL)
    {
        sprintf(fgnu, "%s/%s.gnu", dirnm, froot);
    }
    else
    {
        sprintf(fgnu, "./%s.gnu", froot);
    }
    ofl = fopen(fgnu, "w");
    fprintf(ofl, "#!/usr/bin/gnuplot -persist\n");
    fprintf(ofl, "set pm3d map\n");
    fprintf(ofl, "set xlabel 't*'\n");
    fprintf(ofl, "set ylabel 'Depth (km)'\n");
    fprintf(ofl, "set title 'Cross-Correlation'\n");
    fprintf(ofl, "set grid\n");
    fprintf(ofl, "set xrange [%f:%f]\n", tds.tstar[0],
                                         tds.tstar[MAX(0, tds.ntstar-1)]);
    fprintf(ofl, "set yrange [%f:%f]\n", tds.depths[MAX(0, tds.ndepth-1)],
                                         tds.depths[0]);
    fprintf(ofl, "splot '%s.surf_txt' u 1:2:3\n", froot);
    fclose(ofl);
    // change permissions
    memset(cmd, 0, 256*sizeof(char));
    sprintf(cmd, "chmod 0755 %s", fgnu);
    system(cmd);
    return 0;
}
//============================================================================//
/*!
 * @brief Makes a shifted synthetic seismogram corresponding to the
 *        given observation, t*, and depth.
 *
 * @param[in] iobs    Observation index in the range [0,nobs-1].
 * @param[in] it      t* index in the range [0,tds.ntstar-1].
 * @param[in] id      Depth index in the range [0,tds.ndepth].
 * @param[in] data    Contains the observed seismograms.
 * @param[in] grns    Contains the Green's functions used to make
 *                    the synthetic seisogram.
 * @param[in] tds     Contains the time lags computed in the grid search.
 *
 * @param[out] synth  Structure containing the synthetic seismogram.

 * @result 0 indicates success.
 *
 * @ingroup tdsearch_gridsearch
 */
int tdSearch_gridSearch_makeSACSynthetic(
    const int iobs, const int it, const int id,
    const struct tdSearchData_struct data,
    const struct tdSearchGreens_struct grns,
    const struct tdSearch_struct tds,
    struct sacData_struct *synth)
{
    double epoch;
    int idt, ierr, igrns, npGrns;
    memset(synth, 0, sizeof(struct sacData_struct));
    // Pick off the header
    if (it < 0 || it >= tds.ntstar || id < 0 || id >= tds.ndepth)
    {
        fprintf(stderr, "%s: Invalid (t*,depth) index (%d,%d)\n", __func__, it, id);
        return - 1;
    }
    if (tds.synthetic == NULL || tds.lags == NULL)
    {
        if (tds.synthetic == NULL)
        {
            fprintf(stderr, "%s: Error tds.synethetic is NULL\n", __func__);
        }
        if (tds.lags == NULL)
        {
            fprintf(stderr, "%s: Error tds.lags is NULL\n", __func__);
        }
        return -1;
    }
    idt = tdsearch_gridSearch_gridToIndex(it, id, tds);
    igrns = tdsearch_greens_getGreensFunctionIndex(G11_GRNS, iobs,
                                                   it, id, grns);
    if (idt < 0 || igrns < 0)
    {
        fprintf(stderr, "%s: Failed getting an index\n", __func__);
        return -1;
    }
    // Copy header
    sacio_copyHeader(grns.grns[igrns].header, &synth->header);
    // Copy corresponding synthetic
    synth->npts = grns.grns[igrns].npts;
    synth->data = sacio_malloc64f(synth->npts);
    npGrns = tds.grnsMatrixPtr[idt+1] - tds.grnsMatrixPtr[idt];
    if (npGrns != synth->npts)
    {
        fprintf(stderr, "%s: Size mismatch\n", __func__);
        sacio_free(synth);
        return -1;
    }
    array_copy64f_work(npGrns, &tds.synthetic[tds.grnsMatrixPtr[idt]],
                       synth->data);
    // Get the component inline with the estimate
    sacio_setCharacterHeader(SAC_CHAR_KCMPNM, data.obs[iobs].header.kcmpnm,
                             &synth->header);
    // Fix the timing
    ierr = sacio_getEpochalStartTime(synth->header, &epoch);
    if (ierr != 0)
    {
        fprintf(stderr, "%s: Failed to get start time\n", __func__);
        return -1;
    }
    epoch = epoch + (double) tds.lags[idt]*tds.dt;
    sacio_setEpochalStartTime(epoch, &synth->header);
    return 0;
}
//============================================================================//