text:
#ifndef PYGSL_STATMODULE_H
#define PYGSL_STATMODULE_H 1

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

/*
 * This api will be only exported to the various  statistics
 * modules.
 */
static void **PyGSL_STATISTICS_API = NULL;

#define PyGSL_STATISTICS_d_A_NUM    0
#define PyGSL_STATISTICS_l_A_NUM    1
#define PyGSL_STATISTICS_d_Ad_NUM   2
#define PyGSL_STATISTICS_d_AA_NUM   3
#define PyGSL_STATISTICS_d_AAd_NUM  4
#define PyGSL_STATISTICS_d_AAdd_NUM 5
#define PyGSL_STATISTICS_d_Add_NUM  6
#define PyGSL_STATISTICS_ll_A_NUM   7

#define PyGSL_STATISTICS_d_A_PROTO  (PyObject *self, PyObject *args, \
		     double (*pointer)(const void *, size_t, size_t),\
		     int array_type, int basis_type_size)

#define PyGSL_STATISTICS_l_A_PROTO   (PyObject *self, PyObject *args, \
		      size_t (*pointer)(const void *, size_t, size_t), \
		      int array_type, int basis_type_size)

#define PyGSL_STATISTICS_d_Ad_PROTO  (PyObject *self, PyObject *args, \
		       double (*pointer)(const void *, size_t, size_t, double), \
		       int array_type, int basis_type_size)

#define PyGSL_STATISTICS_d_AA_PROTO  (PyObject *self, PyObject *args, \
		      double (*pointer)(const void  *, size_t,const void *, size_t, size_t), \
		      int array_type, int basis_type_size)

#define PyGSL_STATISTICS_d_AAd_PROTO (PyObject *self, PyObject *args, \
		       double (*pointer)(const void *, size_t,const void *, size_t, size_t, double), \
		       int array_type, int basis_type_size)

#define PyGSL_STATISTICS_d_AAdd_PROTO (PyObject *self, PyObject *args, \
			double (*pointer)(const void *, size_t,const void *, size_t, size_t, double, double), \
			int array_type, int basis_type_size)

#define PyGSL_STATISTICS_d_Add_PROTO (PyObject *self, PyObject *args, \
		       double (*pointer)(const void *, size_t, size_t, double, double), \
		       int array_type, int basis_type_size)

#define PyGSL_STATISTICS_ll_A_PROTO  (PyObject *self, PyObject *args, \
		      void (*pointer)(size_t *, size_t *, const void *, size_t, size_t), \
		      int array_type, int basis_type_size)

#if defined(PyGSL_STATISTICS_IMPORT_API)
extern PyObject *  PyGSL_statistics_d_A    PyGSL_STATISTICS_d_A_PROTO;
extern PyObject *  PyGSL_statistics_l_A    PyGSL_STATISTICS_l_A_PROTO;
extern PyObject *  PyGSL_statistics_d_Ad   PyGSL_STATISTICS_d_Ad_PROTO;
extern PyObject *  PyGSL_statistics_d_AA   PyGSL_STATISTICS_d_AA_PROTO;
extern PyObject *  PyGSL_statistics_d_AAd  PyGSL_STATISTICS_d_AAd_PROTO;
extern PyObject *  PyGSL_statistics_d_AAdd PyGSL_STATISTICS_d_AAdd_PROTO;
extern PyObject *  PyGSL_statistics_d_Add  PyGSL_STATISTICS_d_Add_PROTO;
extern PyObject *  PyGSL_statistics_ll_A   PyGSL_STATISTICS_ll_A_PROTO;

#define PyGSL_statistics_d_A    (*(PyObject* (*) PyGSL_STATISTICS_d_A_PROTO)    PyGSL_STATISTICS_API[PyGSL_STATISTICS_d_A_NUM])
#define PyGSL_statistics_l_A    (*(PyObject* (*) PyGSL_STATISTICS_l_A_PROTO)    PyGSL_STATISTICS_API[PyGSL_STATISTICS_l_A_NUM])
#define PyGSL_statistics_d_Ad   (*(PyObject* (*) PyGSL_STATISTICS_d_Ad_PROTO)   PyGSL_STATISTICS_API[PyGSL_STATISTICS_d_Ad_NUM])
#define PyGSL_statistics_d_AA   (*(PyObject* (*) PyGSL_STATISTICS_d_AA_PROTO)   PyGSL_STATISTICS_API[PyGSL_STATISTICS_d_AA_NUM])
#define PyGSL_statistics_d_AAd  (*(PyObject* (*) PyGSL_STATISTICS_d_AAd_PROTO)  PyGSL_STATISTICS_API[PyGSL_STATISTICS_d_AAd_NUM])
#define PyGSL_statistics_d_AAdd (*(PyObject* (*) PyGSL_STATISTICS_d_AAdd_PROTO) PyGSL_STATISTICS_API[PyGSL_STATISTICS_d_AAdd_NUM])
#define PyGSL_statistics_d_Add  (*(PyObject* (*) PyGSL_STATISTICS_d_Add_PROTO)  PyGSL_STATISTICS_API[PyGSL_STATISTICS_d_Add_NUM])
#define PyGSL_statistics_ll_A   (*(PyObject* (*) PyGSL_STATISTICS_ll_A_PROTO)   PyGSL_STATISTICS_API[PyGSL_STATISTICS_ll_A_NUM])

#else
static PyObject *  PyGSL_statistics_d_A    PyGSL_STATISTICS_d_A_PROTO;
static PyObject *  PyGSL_statistics_l_A    PyGSL_STATISTICS_l_A_PROTO;
static PyObject *  PyGSL_statistics_d_Ad   PyGSL_STATISTICS_d_Ad_PROTO;
static PyObject *  PyGSL_statistics_d_AA   PyGSL_STATISTICS_d_AA_PROTO;
static PyObject *  PyGSL_statistics_d_AAd  PyGSL_STATISTICS_d_AAd_PROTO;
static PyObject *  PyGSL_statistics_d_AAdd PyGSL_STATISTICS_d_AAdd_PROTO;
static PyObject *  PyGSL_statistics_d_Add  PyGSL_STATISTICS_d_Add_PROTO;
static PyObject *  PyGSL_statistics_ll_A   PyGSL_STATISTICS_ll_A_PROTO;

#endif
#define import_pygsl_stats() \
{ \
   PyObject *pygsl = NULL, *c_api = NULL, *md = NULL; \
   if ( \
      (pygsl = PyImport_ImportModule("pygsl.statistics._stat"))   != NULL && \
      (md = PyModule_GetDict(pygsl))                              != NULL && \
      (c_api = PyDict_GetItemString(md, "_PYGSL_STATISTICS_API")) != NULL && \
      (PyCObject_Check(c_api))                                        \
     ) { \
	 PyGSL_STATISTICS_API = (void **)PyCObject_AsVoidPtr(c_api); \
   } else { \
        fprintf(stderr, "Could not init pygsl.statistics._stat!\n"); \
        PyGSL_STATISTICS_API = NULL; \
   } \
   DEBUG_MESS(2, "PyGSL_API points to %p and PyGSL_STATISTICS_API points to %p\n", (void *) PyGSL_API, (void *) PyGSL_STATISTICS_API);  \
}

#endif /* PYGSL_STATMODULE_H */

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/**
 *
 * @file core_zsetvar.c
 *
 *  PLASMA core_blas kernel
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @version 2.6.0
 * @author Mark Gates
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <lapacke.h>
#include "common.h"

/***************************************************************************//**
 *
 * @ingroup CORE_PLASMA_Complex64_t
 *
 * CORE_zsetvar sets a single variable, x := alpha.
 *
 *******************************************************************************
 *
 * @param[in] alpha
 *         Scalar to set x to, passed by pointer so it can depend on runtime value.
 *
 * @param[out] x
 *         On exit, x = alpha.
 *
 ******************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_zsetvar = PCORE_zsetvar
#define CORE_zsetvar PCORE_zsetvar
#endif
void CORE_zsetvar(const PLASMA_Complex64_t *alpha, PLASMA_Complex64_t *x)
{
    *x = *alpha;
}

ext: c


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/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_spline.c"
#include "../cosmolike_core/theory/halo.c"
#include "../cosmolike_core/theory/HOD.c"
#include "../cosmolike_core/theory/pt.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/init_baryon.c"
#include "init_DESxPlanck.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/theory/CMBxLSS_fourier.c"

typedef struct input_nuisance_params_des {
    double bias[10];
    // double bias2[10];
    double b_mag[10];
    double source_z_bias[10];
    double lens_z_bias[10];
    double shear_m[10];
    double A_ia;
    double eta_ia;
    // double bary[3];
} input_nuisance_params_des;

typedef struct input_cosmo_params_des {
    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_des;

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 BMAG1, double BMAG2, double BMAG3, double BMAG4,double BMAG5, double BMAG6, double BMAG7, double BMAG8, double BMAG9, double BMAG10, double SP1, double SP2, double SP3, double SP4, double SP5,double SP6, double SP7, double SP8, double SP9, double SP10, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double eta_ia);
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 BMAG1, double BMAG2, double BMAG3, double BMAG4,double BMAG5, double BMAG6, double BMAG7, double BMAG8, double BMAG9, double BMAG10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double eta_ia);
void write_datavector_wrapper(char *details, input_cosmo_params_des ic, input_nuisance_params_des in);
double log_like_wrapper(input_cosmo_params_des ic, input_nuisance_params_des 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.IA==4){C = C_shear_shear_IA_tab(ell,z1,z2);}
  else{printf("only support IA==4!\n");exit(1);}
  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.IA==4){C = C_ggl_IA_tab(ell,zl,zs);}
  else{printf("only support IA==4!\n");exit(1);}
  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(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);}
      if (mask(Ncl*nz+i)){ 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)){
      if (mask(start+(Ncl*nz)+i)){
        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))
      if (mask(start+(Ncl*nz)+i)) {
        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)
{
   for (int nz=0; nz<tomo.clustering_Nbin; nz++){
      for (int i=0; i<like.Ncl; i++){
         // if (ell[i]<like.lmax_kappacmb && test_kmax(ell[i],nz)){
         if (mask(start+(like.Ncl*nz)+i)) {
            data[start+(like.Ncl*nz)+i] = C_gk(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) {
         if (mask(start+(like.Ncl*nz)+i)) {
            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){
      if (mask(start+i)) {
         data[start+i] = C_kk(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.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;
  //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)
{
  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)
{
  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 eta_ia)
{
  nuisance.A_ia=A_ia;
  nuisance.eta_ia=eta_ia;
  nuisance.oneplusz0_ia = 1.62;
  if (nuisance.A_ia < 0.0 || nuisance.A_ia > 10.0) return 0;
  if (nuisance.eta_ia < -10.0 || nuisance.eta_ia> 10.0) 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 < tomo.clustering_Nbin; i++){
  //    printf("in set routine %d %le\n",i,gbias.b[i]);
    if (gbias.b[i] < 0.4 || gbias.b[i] > 3.0) return 0;
  }

  return 1;
}

int set_nuisance_bmag(double B1, double B2, double B3, double B4,double B5, double B6, double B7, double B8,double B9, double B10)
{

  int i;
  gbias.b_mag[0] = B1;
  gbias.b_mag[1] = B2;
  gbias.b_mag[2] = B3;
  gbias.b_mag[3] = B4;
  gbias.b_mag[4] = B5;
  gbias.b_mag[5] = B6;
  gbias.b_mag[6] = B7;
  gbias.b_mag[7] = B8;
  gbias.b_mag[8] = B9;
  gbias.b_mag[9] = B10;
  for (i = 0; i < tomo.clustering_Nbin; i++){
  //    printf("in set routine %d %le\n",i,gbias.b[i]);
    if (gbias.b_mag[i] < -5.0 || gbias.b_mag[i] > 5.0) return 0;
  }

  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 BMAG1, double BMAG2, double BMAG3, double BMAG4,double BMAG5, double BMAG6, double BMAG7, double BMAG8, double BMAG9, double BMAG10, double SP1, double SP2, double SP3, double SP4, double SP5, double SP6, double SP7, double SP8, double SP9, double SP10, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double eta_ia)
{
  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)==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)==0){
    printf("Clustering photo-z sigma too small\n");
    return -1.0e15;
  }
  if (set_nuisance_ia(A_ia,eta_ia)==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_bmag(BMAG1,BMAG2,BMAG3,BMAG4,BMAG5,BMAG6,BMAG7,BMAG8,BMAG9,BMAG10)==0){
    printf("Bmag 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;
  // }

  // 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;

  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.gk==1) {
    set_data_gk(ell, pred, start);
    start += like.Ncl*tomo.clustering_Nbin;
  }
  if(like.ks==1) {
    set_data_ks(ell, pred, start);
    start += like.Ncl*tomo.shear_Nbin;
  }
  if(like.kk==1) {
    set_data_kk(ell, pred, start);
    start += like.Ncl;
  }

  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_mask(1,i,j)*(pred[j]-data_read(1,j));
      // if(a>10){printf("a,i,j: %le, %d, %d, %le, %le, %le\n",a,i,j,pred[i]-data_read(1,i),pred[j]-data_read(1,j), invcov_mask(1,i,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));
    // }
  }
  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 BMAG1, double BMAG2, double BMAG3, double BMAG4,double BMAG5, double BMAG6, double BMAG7, double BMAG8, double BMAG9, double BMAG10, double SP1, double SP2, double SP3, double SP4, double SP5,double SP6, double SP7, double SP8, double SP9, double SP10, double CP1, double CP2, double CP3, double CP4, double CP5, double CP6, double CP7, double CP8, double CP9, double CP10, double M1, double M2, double M3, double M4, double M5, double M6, double M7, double M8, double M9, double M10, double A_ia, double eta_ia)
{

  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);
  set_nuisance_clustering_photoz(CP1,CP2,CP3,CP4,CP5,CP6,CP7,CP8,CP9,CP10);
  set_nuisance_ia(A_ia,eta_ia);
  set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10);
  set_nuisance_bmag(BMAG1,BMAG2,BMAG3,BMAG4,BMAG5,BMAG6,BMAG7,BMAG8,BMAG9,BMAG10);
  // 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.gk==1) {
    printf("Computing data vector: gk\n");
    set_data_gk(ell, pred, start);
    start += like.Ncl * tomo.clustering_Nbin;
  }
  if(like.ks==1) {
    printf("Computing data vector: ks\n");
    set_data_ks(ell, pred, start);
    start += like.Ncl * tomo.shear_Nbin;
  }
  if (like.kk) {
    printf("Computing data vector: kk\n");
    set_data_kk(ell, pred, start);
    start += like.Ncl;
  }

  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_des ic, input_nuisance_params_des 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.b_mag[0], in.b_mag[1], in.b_mag[2], in.b_mag[3],in.b_mag[4], in.b_mag[5], in.b_mag[6], in.b_mag[7],in.b_mag[8], in.b_mag[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.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.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.eta_ia);
}

double log_like_wrapper(input_cosmo_params_des ic, input_nuisance_params_des 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.b_mag[0], in.b_mag[1], in.b_mag[2], in.b_mag[3],in.b_mag[4], in.b_mag[5], in.b_mag[6], in.b_mag[7],in.b_mag[8], in.b_mag[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.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.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.eta_ia);
  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"};
  // char tomo_binning_lens[2][200]={"LSST_gold","LSST_gold"};

  // char source_zfile[2][400]={"src_LSSTY1","src_LSSTY6"};
  // char lens_zfile[2][400]={"lens_LSSTY1", "lens_LSSTY6"};

  init_cosmo_runmode("halofit");
  // init_bary(argv[2]);
  init_binning_fourier(20,30.0,3000.0,3000.0,20.0);

  // 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,"gaussian","gaussian",tomo_binning_source[sce],tomo_binning_lens[sce]);

  double b1[10]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, b2[10] ={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},b_mag[10] ={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  b1[0] = 1.7;
  b1[1] = 1.7;
  b1[2] = 1.7;
  b1[3] = 2.0;
  b1[4] = 2.0;

  b_mag[0] = -0.19375;
  b_mag[1] = -0.6285407;
  b_mag[2] = -0.69319886;
  b_mag[3] = 1.17735723;
  b_mag[4] = 1.87509758;

  // init_source_sample_mpp("./zdistris/mcal_1101_source.nz",4);
  init_source_sample_mpp("./zdistris/source_DESY6_3bins.nz",3);
  init_lens_sample_mpp("./zdistris/mcal_1101_lens.nz",5);
  init_IA_mpp(4);
  init_probes("6x2pt");

  init_cmb("planck");
  compute_data_vector("desy6xplanck_6x2pt_3bins_2.txt",0.25,0.75,0.97,-1.,0.,0.048,0.69,0.,0.,\
    b1[0],b1[1],b1[2],b1[3],b1[4],\
    b1[5],b1[6],b1[7],b1[8],b1[9],\
    b_mag[0],b_mag[1],b_mag[2],b_mag[3],b_mag[4],\
    b_mag[5],b_mag[6],b_mag[7],b_mag[8],b_mag[9],\
    0.0,0.0,0.0,0.0,0.0,\
    0.0,0.0,0.0,0.0,0.0,\
    0.0,0.0,0.0,0.0,0.0,\
    0.0,0.0,0.0,0.0,0.0,\
    0.0,0.0,0.0,0.0,0.0,\
    0.0,0.0,0.0,0.0,0.0,\
    0.5,0.0);

  // init_data_cov_mask("cov/cov_desy3xplanck","datav/6x2pt_desy3xplanck_6x2pt.txt", "datav/xi_Y3_6x2pt.mask");
  // // init_data_cov_mask("cov/cov_desy6xplanck_3src","datav/6x2pt_desy6xplanck_6x2pt_3bins.txt", "datav/xi_Y6_6x2pt_3src.mask");

  // for(double omm=0.2;omm<0.5; omm+=0.01){
  // printf("%le, %le\n", omm, log_multi_like(omm,0.82355,0.97,-1.,0.,0.048,0.69,0.,0.,\
  //   b1[0],b1[1],b1[2],b1[3],b1[4],\
  //   b1[5],b1[6],b1[7],b1[8],b1[9],\
  //   b_mag[0],b_mag[1],b_mag[2],b_mag[3],b_mag[4],\
  //   b_mag[5],b_mag[6],b_mag[7],b_mag[8],b_mag[9],\
  //   0.0,0.0,0.0,0.0,0.0,\
  //   0.0,0.0,0.0,0.0,0.0,\
  //   0.0,0.0,0.0,0.0,0.0,\
  //   0.0,0.0,0.0,0.0,0.0,\
  //   0.0,0.0,0.0,0.0,0.0,\
  //   0.0,0.0,0.0,0.0,0.0,\
  //   0.5,0.0));
  // exit(0);
  // }
  return 0;
}

max_stars_repo_name: CosmoLike/LSSTxSO

max_issues_repo_name: CosmoLike/LSSTxSO

max_forks_repo_name: CosmoLike/LSSTxSO


text:
#include <stdio.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_multimin.h>
/* Paraboloid with a minimum at f(dp[0],dp[1]) = 30.0 */

double
my_f (const gsl_vector * v, void *params)
{
  double x, y;
  double *dp = (double *) params;

  x = gsl_vector_get (v, 0);
  y = gsl_vector_get (v, 1);

  return 10.0 * (x - dp[0]) * (x - dp[0]) +
    20.0 * (y - dp[1]) * (y - dp[1]) + 30.0;
}

/* The gradient of f, df = (df/dx, df/dy). */
void
my_df (const gsl_vector * v, void *params, gsl_vector * df)
{
  double x, y;
  double *dp = (double *) params;

  x = gsl_vector_get (v, 0);
  y = gsl_vector_get (v, 1);

  gsl_vector_set (df, 0, 20.0 * (x - dp[0]));
  gsl_vector_set (df, 1, 40.0 * (y - dp[1]));
}

/* Now both f and df together. */
void
my_fdf (const gsl_vector * x, void *params, double *f, gsl_vector * df)
{
  *f = my_f (x, params);
  my_df (x, params, df);
}

int
main (void)
{
  size_t iter = 0;
  int status;

  const gsl_multimin_fdfminimizer_type *T;
  gsl_multimin_fdfminimizer *s;

  /* Position of the minimum (1,2). */
  double par[2] = { 1.0, 2.0 };

  gsl_vector *x;
  gsl_multimin_function_fdf my_func;

  my_func.f = &my_f;
  my_func.df = &my_df;
  my_func.fdf = &my_fdf;
  my_func.n = 2;
  my_func.params = ∥

  /* Starting point, x = (5,7) */

  x = gsl_vector_alloc (2);
  gsl_vector_set (x, 0, 5.0);
  gsl_vector_set (x, 1, 7.0);

  T = gsl_multimin_fdfminimizer_conjugate_fr;
  s = gsl_multimin_fdfminimizer_alloc (T, 2);

  gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);

  do
    {
      iter++;
      status = gsl_multimin_fdfminimizer_iterate (s);

      if (status)
        break;

      status = gsl_multimin_test_gradient (s->gradient, 1e-3);

      if (status == GSL_SUCCESS)
        printf ("Minimum found at:\n");

      printf ("%5d %.5f %.5f %10.5f\n", iter,
              gsl_vector_get (s->x, 0), gsl_vector_get (s->x, 1), s->f);

    }
  while (status == GSL_CONTINUE && iter < 100);

  gsl_multimin_fdfminimizer_free (s);
  gsl_vector_free (x);

  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
#ifndef SCALED_VARIANCE_H
#define SCALED_VARIANCE_H

#include "smash/decaymodes.h"
#include "smash/particletype.h"

#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_vector.h>

bool load_particle_types();

class ScaledVarianceCalculator {
 public:
  ScaledVarianceCalculator(const smash::ParticleTypePtrList& types_in_the_box,
                           double T, double mub, double mus, double muq,
                           bool energy_conservation,
                           bool B_conservation,
                           bool S_conservation,
                           bool Q_conservation,
                           bool quantum_statistics) :
    T_ (T),
    mub_ (mub),
    mus_ (mus),
    muq_ (muq),
    energy_conservation_ (energy_conservation),
    B_conservation_ (B_conservation),
    S_conservation_ (S_conservation),
    Q_conservation_ (Q_conservation),
    quantum_statistics_ (quantum_statistics),
    quantum_series_max_terms_(100),
    quantum_series_rel_precision_(1e-12) {
    all_types_in_the_box_.clear();
    for (const smash::ParticleTypePtr t : types_in_the_box) {
      all_types_in_the_box_.push_back(t);
    }
  }
  void set_T (double T) { T_ = T; }
  void set_mub (double mub) { mub_ = mub; }
  void set_mus (double mus) { mus_ = mus; }
  void set_muq (double muq) { muq_ = muq; }
  void set_B_conservation (bool B_cons) { B_conservation_ = B_cons; }
  void set_S_conservation (bool S_cons) { S_conservation_ = S_cons; }
  void set_Q_conservation (bool Q_cons) { Q_conservation_ = Q_cons; }
  void set_quantum_statistics (bool qs) { quantum_statistics_ = qs; }

  struct solver_params {
    smash::ParticleTypePtrList* types;
    double e;   // energy density
    double nb;  // baryon number density
    double ns;  // strangeness density
    double nq;  // charge density
  };
  static int set_eos_solver_equations(const gsl_vector* x, void* params,
                                      gsl_vector* f);
  std::string print_solver_state(size_t iter,
                                 gsl_multiroot_fsolver* solver) const;
  void setTmu_from_conserved(double Etot, double V,
                             double B, double S, double Q);

  /**
   * Assumes that there is a gas of particles species defined
   * by all_types_in_the_box_. This function computes fluctuations
   * of the total number of particles belonging to a subset of species
   * defined by type_of_interest function. Fluctuations depend on
   * the ensemble and which conservation laws it respects. This is
   * determined by variables energy_conservation_, B_conservation_,
   * S_conservation_, and Q_conservation_.
   *
   * For reference see arxiv.org/pdf/0706.3290.pdf, Eq. (29-40).
   * The approximations of this paper are only valid for large
   * enough volume.
   *
   * \param[in] type_of_interest Function true for particle species, for
   *            which the mean and the scaled variance are computed
   * \param[out] mean density and scaled variance of type_of_interest species
   */
  std::pair<double, double> scaled_variance(
    std::function<bool(const smash::ParticleTypePtr)> type_of_interest);

  void prepare_decays();

  friend std::ostream& operator<< (std::ostream&,
                                   const ScaledVarianceCalculator&);

 private:
  /// List of species included in the gas
  smash::ParticleTypePtrList all_types_in_the_box_;
  /**
   *  Structure for holding all decay final state in a form convenient
   *  for statistical caulculations:
   *  resonance ->
   *       vector of decays
   *       each decay = (branching ratio and vector (product -> count))
   *  Sum of branching ratios should be 1 for all particles.
   *  By convention stable particles decay into themselves with brancing ratio 1.
   */
  std::map<smash::ParticleTypePtr,
           std::vector<std::pair<double,
                               std::map<smash::ParticleTypePtr, int>>>>
       all_decay_final_states_;
  /// Temperature of the gas [GeV]
  double T_;
  /// Baryo-chemical potential of the gas [GeV]
  double mub_;
  /// Strangeness chemical potential of the gas [GeV]
  double mus_;
  /// Charge chemical potential of the gas [GeV]
  double muq_;
  /// If ensemble includes energy conservation
  bool energy_conservation_;
  /// If ensemble includes baryon number conservation
  bool B_conservation_;
  /// If ensemble includes strangeness conservation
  bool S_conservation_;
  /// If ensemble includes electric charge conservation
  bool Q_conservation_;
  /// If quantum statistics is included
  bool quantum_statistics_;
  /// Maximal number of terms in the series for quantum formulas
  const unsigned int quantum_series_max_terms_;
  /// Relative precision, at which quantum series summation stops
  const double quantum_series_rel_precision_;
};

#endif // SCALED_VARIANCE_H

ext: h


text:
#ifndef __LISAUTILS_H__
#define __LISAUTILS_H__ 1

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

#include "constants.h"
#include "struct.h"
#include "EOBNRv2HMROMstruct.h"
#include "EOBNRv2HMROM.h"
#include "waveform.h"
#include "wip.h"
#include "splinecoeffs.h"
#include "fresnel.h"
#include "likelihood.h"
#include "LISAFDresponse.h"
#include "LISAnoise.h"

#if defined(__cplusplus)
extern "C" {
#define complex _Complex
#elif 0
} /* so that editors will match preceding brace */
#endif

/***************** Enumerator to choose what masses/time set to sample for *****************/

typedef enum SampleMassParamstag {
  m1m2,
  Mchirpeta
} SampleMassParamstag;
/* Superseded by sampleLframe */
// typedef enum SampleTimeParamtag {
//   tSSB,
//   tL
// } SampleTimeParamtag;

/* Function to convert string input SampleMassParams to tag */
SampleMassParamstag ParseSampleMassParamstag(char* string);
/* Function to convert string input SampleTimeParams to tag */
//SampleTimeParamtag ParseSampleTimeParamtag(char* string);

/***************** Structure definitions *****************/

/* Parameters for the generation of a LISA waveform (in the form of a list of modes) */
typedef struct tagLISAParams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double m1;                 /* mass of companion 1 (solar masses, default 2e6) */
  double m2;                 /* mass of companion 2 (solar masses, default 1e6) */
  double distance;           /* distance of source (Mpc, default 1e3) */
  double lambda;             /* first angle for the position in the sky (rad, default 0) */
  double beta;               /* second angle for the position in the sky (rad, default 0) */
  double inclination;        /* inclination of L relative to line of sight (rad, default PI/3) */
  double polarization;       /* polarization angle (rad, default 0) */
  int nbmode;                /* number of modes to generate (starting with 22) - defaults to 5 (all modes) */
} LISAParams;

/* Global parameters for the waveform generation and overlap computation */
typedef struct tagLISAGlobalParams {
  double fRef;               /* reference frequency (Hz, default 0 which is interpreted as Mf=0.14) */
  double deltatobs;          /* max duration of observation (years, default 2) - the start of the signals might be cut in time instead of cut in frequency */
  double minf;               /* Minimal frequency (Hz, default=0) - when set to 0, use the lowest frequency where the detector noise model is trusted __LISASimFD_Noise_fLow (set somewhat arbitrarily)*/
  double maxf;               /* Maximal frequency (Hz, default=0) - when set to 0, use the highest frequency where the detector noise model is trusted __LISASimFD_Noise_fHigh (set somewhat arbitrarily)*/
  int tagextpn;              /* Tag to allow PN extension of the waveform at low frequencies */
  int tagtRefatLISA;         /* Tag to signal time to be referenced to arrival at LISA rather than at SSB */
  double Mfmatch;            /* When PN extension allowed, geometric matching frequency: will use ROM above this value. If <=0, use ROM down to the lowest covered frequency */
  int setphiRefatfRef;       /* Flag for adjusting the FD phase at phiRef at the given fRef, which depends also on tRef - if false, treat phiRef simply as an orbital phase shift (minus an observer phase shift) (default=1) */
  int nbmodeinj;             /* number of modes to include in the injection (starting with 22) - defaults to 5 (all modes) */
  int nbmodetemp;            /* number of modes to include in the templates (starting with 22) - defaults to 5 (all modes) */
  int tagint;                /* Tag choosing the integrator: 0 for wip (default), 1 for linear integration */
  TDItag tagtdi;             /* Tag choosing the TDI variables to use */
  int nbptsoverlap;          /* Number of points to use in loglinear overlaps (default 32768) */
  LISAconstellation *variant;  /* A structure defining the LISA constellation features */
  int zerolikelihood;        /* Tag to zero out the likelihood, to sample from the prior for testing purposes (default 0) */
  int frozenLISA;            /* Freeze the orbital configuration to the time of peak of the injection (default 0) */
  ResponseApproxtag responseapprox;    /* Approximation in the GAB and orb response - choices are full (full response, default), lowfL (keep orbital delay frequency-dependence but simplify constellation response) and lowf (simplify constellation and orbital response) - WARNING : at the moment noises are not consistent, and TDI combinations from the GAB are unchanged */
  int tagsimplelikelihood22; /* Tag to use simplified, frozen-LISA and lowf likelihood where mode overlaps are precomputed - 22-mode only - can only be used when the masses and time (tL) are pinned to injection values (Note: when using --snr, distance adjustment done using responseapprox, not the simple response) */
  int tagsimplelikelihoodHM; /* Tag to use simplified, frozen-LISA and lowf likelihood where mode overlaps are precomputed - set of modes - can only be used when the masses and time (tL) are pinned to injection values (Note: when using --snr, distance adjustment done using responseapprox, not the simple response) */
} LISAGlobalParams;

typedef struct tagLISASignalCAmpPhase
{
  struct tagListmodesCAmpPhaseFrequencySeries* TDI1Signal;   /* Signal in the TDI channel 1, in the form of a list of the contribution of each mode */
  struct tagListmodesCAmpPhaseFrequencySeries* TDI2Signal;   /* Signal in the TDI channel 2, in the form of a list of the contribution of each mode */
  struct tagListmodesCAmpPhaseFrequencySeries* TDI3Signal;   /* Signal in the TDI channel 3, in the form of a list of the contribution of each mode */
  double TDI123hh;                                           /* Combined Inner product (h|h) for TDI channels 123 */
} LISASignalCAmpPhase;

typedef struct tagLISAInjectionCAmpPhase
{
  struct tagListmodesCAmpPhaseSpline* TDI1Splines;   /* Signal in the TDI channel 1, in the form of a list of splines for the contribution of each mode */
  struct tagListmodesCAmpPhaseSpline* TDI2Splines;   /* Signal in the TDI channel 2, in the form of a list of splines for the contribution of each mode */
  struct tagListmodesCAmpPhaseSpline* TDI3Splines;   /* Signal in the TDI channel 3, in the form of a list of splines for the contribution of each mode */
  double TDI123ss;                                   /* Combined Inner product (s|s) for TDI channels 123 */
} LISAInjectionCAmpPhase;

typedef struct tagLISASignalReIm /* We don't store the SNRs here, as we will use -1/2(h-s|h-s) for the likelihood */
{
  struct tagReImFrequencySeries* TDI1Signal;   /* Signal in the TDI channel 1, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* TDI2Signal;   /* Signal in the TDI channel 2, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* TDI3Signal;   /* Signal in the TDI channel 3, in the form of a Re/Im frequency series where the modes have been summed */
} LISASignalReIm;

typedef struct tagLISAInjectionReIm /* Storing the vectors of frequencies and noise values - We don't store the SNRs here, as we will use -1/2(h-s|h-s) for the likelihood */
{
  struct tagReImFrequencySeries* TDI1Signal;   /* Signal in the TDI channel 1, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* TDI2Signal;   /* Signal in the TDI channel 2, in the form of a Re/Im frequency series where the modes have been summed */
  struct tagReImFrequencySeries* TDI3Signal;   /* Signal in the TDI channel 3, in the form of a Re/Im frequency series where the modes have been summed */
  gsl_vector* freq;                            /* Vector of frequencies of the injection (assumed to be the same for A,E,T) */
  gsl_vector* noisevalues1;                    /* Vector of noise values on freq for TDI channel 1 */
  gsl_vector* noisevalues2;                    /* Vector of noise values on freq for TDI channel 2 */
  gsl_vector* noisevalues3;                    /* Vector of noise values on freq for TDI channel 3 */
} LISAInjectionReIm;

typedef struct tagLISAPrior {
  SampleMassParamstag samplemassparams;   /* Choose the set of mass params to sample from - options are m1m2 and Mchirpeta (default m1m2) */
  //SampleTimeParamtag sampletimeparam;     /* Choose the time param to sample from - options are tSSB and tL (default tSSB) */
  int sampleLframe;          /* flag to sample L-frame params tL, lambdaL, betaL, psiL instead of SSB-frame params -- priors are interpreted for those L-frame params - no phase transformation */
  double deltaT;             /* width of time prior centered on injected value (s) (default 1e5) */
  double comp_min;           /* minimum component mass (solar masses) (default 1e4) */
  double comp_max;           /* maximum component mass (solar masses) (default 1e8) */
  double mtot_min;           /* minimum total mass (solar masses) (default 5*1e4) */
  double mtot_max;           /* maximum total mass (solar masses) (default 1e8) */
  double qmax;               /* maximum asymmetric mass ratio (>=1) (default 11.98) */
  double Mchirp_min;         /* Minimum chirp mass in Solar masses - when sampling Mchirpeta (default=2e4) */
  double Mchirp_max;         /* Maximum chirp mass in Solar masses - when sampling Mchirpeta (default=4e7) */
  double eta_min;            /* Minimum symmetric mass ratio eta - when sampling Mchirpeta (default=0.072) */
  double eta_max;            /* Maximum symmetric mass ratio eta - when sampling Mchirpeta (default=0.25) */
  double dist_min;           /* minimum distance of source (Mpc) (default 100) */
  double dist_max;           /* maximum distance of source (Mpc) (default 40*1e3) */
  double lambda_min;         /* minimum lambda (rad, default 0) - for testing */
  double lambda_max;         /* maximum lambda (rad, default 2pi) - for testing */
  double beta_min;           /* minimum beta (rad, default 0) - for testing */
  double beta_max;           /* maximum beta (rad, default pi) - for testing */
  double phase_min;          /* minimum phase (rad, default 0) - for testing */
  double phase_max;          /* maximum phase (rad, default 2pi) - for testing */
  double pol_min;            /* minimum polarization (rad, default 0) - for testing */
  double pol_max;            /* maximum polarization (rad, default 2pi) - for testing */
  double inc_min;            /* minimum inclination (rad, default 0) - for testing */
  double inc_max;            /* maximum inclination (rad, default pi) - for testing */
  double fix_m1;
  double fix_m2;
  double fix_Mchirp;
  double fix_eta;
  double fix_time;
  double fix_lambda;
  double fix_beta;
  double fix_phase;
  double fix_pol;
  double fix_dist;
  double fix_inc;
  int pin_m1;
  int pin_m2;
  int pin_Mchirp;
  int pin_eta;
  int pin_time;
  int pin_lambda;
  int pin_beta;
  int pin_phase;
  int pin_pol;
  int pin_dist;
  int pin_inc;
  double snr_target;
  int rescale_distprior;
  int flat_distprior;
  int logflat_massprior;
} LISAPrior;

typedef struct tagLISARunParams {
  double eff;                /* target efficiency (default 0.1) */
  double tol;                /* logZ tolerance (default 0.5) */
  int    consteff;           /* constant efficiency mode (default 0) */
  int    nlive;              /* number of live points (default 1000) */
  int    writeparams;        /* Write params - if 1, write run parameters to file (default 1) */
  char   outroot[200];       /* output root (default "chains/LISAinference_") */
  int    bambi;              /* run BAMBI? (default 0) */
  int    resume;             /* resume form previous run? (default 0) */
  int    maxiter;            /* max number of iterations (default 0 - ignore) */
  char   netfile[200];       /* NN settings file (default "LISAinference.inp") */
  int    mmodal;             /* use multimodal decomposition ? */
  int    maxcls;             /* max number of modes in multimodal decomposition */
  int    nclspar;            /* number of parameters to use for multimodal decomposition - in the order of the cube */
  double ztol;               /* in multimodal decomposition, modes with lnZ lower than Ztol are ignored */
  int    seed;               /* seed the inference by setting one of the live points to the injection ? */
} LISARunParams;

/* Parameters for the generation of a LISA waveform (in the form of a list of modes) */
typedef struct tagLISAAddParams {
  double tRef;               /* reference time (s) - GPS time at the frequency representing coalescence */
  double phiRef;             /* reference phase (rad) - phase at the frequency representing coalescence (or at fRef if specified) */
  double m1;                 /* mass of companion 1 (solar masses, default 2e6) */
  double m2;                 /* mass of companion 2 (solar masses, default 1e6) */
  double distance;           /* distance of source (Mpc, default 1e3) */
  double lambda;             /* first angle for the position in the sky (rad, default 0) */
  double beta;               /* second angle for the position in the sky (rad, default 0) */
  double inclination;        /* inclination of L relative to line of sight (rad, default PI/3) */
  double polarization;       /* polarization angle (rad, default 0) */
  int loadparamsfile;        /* Option to load physical parameters from file for LISAlikelihood and to output resulting likelihoods to file (default 0) */
  int nlinesparams;          /* Number of lines in params file for LISAlikelihood */
  char indir[256];           /* Input directory for LISAlikelihood */
  char infile[256];          /* Input file for LISAlikelihood */
  char outdir[256];          /* Output directory for LISAlikelihood */
  char outfile[256];         /* Output file for LISAlikelihood */
} LISAAddParams;

/* Saved precomputed values for the injection , when using the simplified frozenLISA lowf likelihood - here for now assumes 22 mode only, old convention */
/* NOTE: difference in convention, new convention has a factor 3 in the integrand defining <lm|l'm'> */
typedef struct tagSimpleLikelihoodPrecomputedValues22 {
  double normalization;
  double complex sa;
  double complex se;
} SimpleLikelihoodPrecomputedValues22;
/* Saved precomputed values for the injection , when using the simplified frozenLISA lowf likelihood - here for set of modes, new convention */
/* NOTE: difference in convention, new convention has a factor 3 in the integrand defining <lm|l'm'> */
typedef struct tagSimpleLikelihoodPrecomputedValuesHM {
  gsl_matrix* Lambda_lm_lpmp; /* Matrix <lm|l'm'> */
  gsl_vector* sa_lm_real;         /* Vector s_a^{lm}, real part*/
  gsl_vector* sa_lm_imag;         /* Vector s_a^{lm}, imaginary part */
  gsl_vector* se_lm_real;         /* Vector s_e^{lm}, real part*/
  gsl_vector* se_lm_imag;         /* Vector s_e^{lm}, imaginary part */
} SimpleLikelihoodPrecomputedValuesHM;

/************ Functions for LISA parameters, injection, likelihood, prior ************/

/* Parse command line to initialize LISAParams, LISAGlobalParams, LISAPrior, and LISARunParams objects */
/* Masses are input in solar masses and distances in Mpc - converted in SI for the internals */
void parse_args_LISA(ssize_t argc, char **argv,
  LISAParams* params,
  LISAGlobalParams* globalparams,
  LISAPrior* prior,
  LISARunParams* run,
  LISAAddParams* addparams);

/* Functions to print the parameters of the run in files for reference */
int print_parameters_to_file_LISA(
  LISAParams* params,
  LISAGlobalParams* globalparams,
  LISAPrior* prior,
  LISARunParams* run);
int print_rescaleddist_to_file_LISA(
  LISAParams* params,
  LISAGlobalParams* globalparams,
  LISAPrior* prior,
  LISARunParams* run);
int print_snrlogZ_to_file_LISA(LISARunParams* run, double SNR, double logZ);
/* Function printing injection/signal parameters to stdout */
void report_LISAParams(LISAParams* params);

/* Initialization and clean-up for LISASignal structures */
void LISASignalCAmpPhase_Cleanup(LISASignalCAmpPhase* signal);
void LISASignalCAmpPhase_Init(LISASignalCAmpPhase** signal);
void LISAInjectionCAmpPhase_Cleanup(LISAInjectionCAmpPhase* signal);
void LISAInjectionCAmpPhase_Init(LISAInjectionCAmpPhase** signal);
void LISASignalReIm_Cleanup(LISASignalReIm* signal);
void LISASignalReIm_Init(LISASignalReIm** signal);
void LISAInjectionReIm_Cleanup(LISAInjectionReIm* signal);
void LISAInjectionReIm_Init(LISAInjectionReIm** signal);

//Function to restrict range of the signal/injection to within desired limits.
int listmodesCAmpPhaseTrim(ListmodesCAmpPhaseFrequencySeries* listSeries);

/* Function generating a LISA signal as a list of modes in CAmp/Phase form, from LISA parameters */
int LISAGenerateSignalCAmpPhase(
  struct tagLISAParams* params,                 /* Input: set of LISA parameters of the signal */
  struct tagLISASignalCAmpPhase* signal);  /* Output: structure for the generated signal */
/* Function generating a LISA injection as a list of modes, given as preinterpolated splines, from LISA parameters */
int LISAGenerateInjectionCAmpPhase(
  struct tagLISAParams* injectedparams,    /* Input: set of LISA parameters of the signal */
  struct tagLISAInjectionCAmpPhase* signal);  /* Output: structure for the generated signal */
/* Function generating a LISA signal as a frequency series in Re/Im form where the modes have been summed, from LISA parameters - takes as argument the frequencies on which to evaluate */
int LISAGenerateSignalReIm(
  struct tagLISAParams* params,       /* Input: set of LISA parameters of the template */
  gsl_vector* freq,                   /* Input: frequencies on which evaluating the waveform (from the injection) */
  struct tagLISASignalReIm* signal);  /* Output: structure for the generated signal */
/* Function generating a LISA injection as a frequency series in Re/Im form where the modes have been summed, from LISA parameters - frequencies on which to evaluate are to be determined internally */
int LISAGenerateInjectionReIm(
  struct tagLISAParams* injectedparams,       /* Input: set of LISA parameters of the injection */
  double fLow,                                /* Input: starting frequency */
  int nbpts,                                  /* Input: number of frequency samples */
  int tagsampling,                            /* Input: tag for using linear (0) or logarithmic (1) sampling */
  struct tagLISAInjectionReIm* signal);       /* Output: structure for the generated signal */

/*Wrapper for waveform generation with possibly a combination of EOBNRv2HMROM and TaylorF2*/
/* Note: GenerateWaveform accepts masses and distances in SI units, whereas LISA params is in solar masses and Mpc */
int GenerateWaveform(
  struct tagListmodesCAmpPhaseFrequencySeries **listhlm,  /* Output: list of modes in Frequency-domain amplitude and phase form */
  int nbmode,                                    /* Number of modes to generate (starting with the 22) */
  double f_match,                                /* Minimum frequency using EOBNRv2HMROM */
  double f_min,                                  /* Minimum frequency required */
  double deltatRef,                              /* Time shift so that the peak of the 22 mode occurs at deltatRef */
  double phiRef,                                 /* Phase at reference frequency */
  double fRef,                                   /* Reference frequency (Hz); 0 defaults to fLow */
  double m1SI,                                   /* Mass of companion 1 (kg) */
  double m2SI,                                   /* Mass of companion 2 (kg) */
  double distance                               /* Distance of source (m) */
		     );

/* Non-spinning merger TaylorF2 waveform, copied and condensed from LAL */
void TaylorF2nonspin(
		double *amp,                            /**< FD waveform amplitude (modulus)*/
		double *phase,                          /**< FD waveform phase */
		const double *freqs,                    /**< frequency points at which to evaluate the waveform (Hz) */
		const int size,                         /** number of freq samples */
		const double m1_SI,                     /**< mass of companion 1 (kg) */
		const double m2_SI,                     /**< mass of companion 2 (kg) */
		const double distance,                  /** distance (m) */
		const double imatch                     /**< index at which to match phase;
							   assumes arrays are preloaded at imatch and imatch+1
							   with the required result */
		     );

/* checks prior boundaires */
int PriorBoundaryCheckm1m2(LISAPrior *prior, double *Cube);
int PriorBoundaryCheckMchirpeta(LISAPrior *prior, double *Cube);

/* Prior functions from Cube to physical parameters
   x1 is min, x2 is max when specified
   r is Cube value */
double CubeToFlatPrior(double r, double x1, double x2);
double CubeToLogFlatPrior(double r, double x1, double x2);
double CubeToPowerPrior(double p, double r, double x1, double x2);
double CubeToGaussianPrior(double r, double mean, double sigma);
double CubeToSinPrior(double r, double x1, double x2);
double CubeToCosPrior(double r, double x1, double x2);
/* Prior functions from physical parameters to Cube
   x1 is min, x2 is max when specified
   y is physical value */
double FlatPriorToCube(double y, double x1, double x2);
double LogFlatPriorToCube(double y, double x1, double x2);
double PowerPriorToCube(double p, double y, double x1, double x2);
double SinPriorToCube(double y, double x1, double x2);
double CosPriorToCube(double y, double x1, double x2);

/* log-Likelihood functions */
double CalculateLogLCAmpPhase(LISAParams *params, LISAInjectionCAmpPhase* injection);
double CalculateLogLReIm(LISAParams *params, LISAInjectionReIm* injection);

/* Functions for simplified likelihood using precomputing relevant values */
int LISAComputeSimpleLikelihoodPrecomputedValues22(SimpleLikelihoodPrecomputedValues22* simplelikelihoodvals22, LISAParams* params);
int LISAComputeSimpleLikelihoodPrecomputedValuesHM(SimpleLikelihoodPrecomputedValuesHM* simplelikelihoodvalsHM, LISAParams* params);
double CalculateLogLSimpleLikelihood22(SimpleLikelihoodPrecomputedValues22* simplelikelihoodvals22, LISAParams* params);
double CalculateLogLSimpleLikelihoodHM(SimpleLikelihoodPrecomputedValuesHM* simplelikelihoodvalsHM, LISAParams* params);

/************ Global Parameters ************/

extern LISAParams* injectedparams;
extern LISAGlobalParams* globalparams;
extern LISAPrior* priorParams;
extern LISAAddParams* addparams;
extern double logZdata; /* TODO: not used */
extern SimpleLikelihoodPrecomputedValues22* simplelikelihoodinjvals22;
extern SimpleLikelihoodPrecomputedValuesHM* simplelikelihoodinjvalsHM;

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

#endif

max_stars_repo_name: JohnGBaker/flare

max_issues_repo_name: JohnGBaker/flare

max_forks_repo_name: JohnGBaker/flare


text:
#pragma once

#define NOMINMAX
#include <windows.h>

#include <unknwn.h>
#include <restrictederrorinfo.h>
#include <hstring.h>
#include <MemoryBuffer.h>

// DirectX
#include <windows.ui.xaml.media.dxinterop.h>
#include <d3d11.h>
#include <dxgi.h>
#include <dxgi1_2.h>
#include <dxgi1_3.h>
#include <DirectXMath.h>

// C++/WinRT
#include "winrt/Microsoft.UI.Xaml.Automation.Peers.h"
#include "winrt/Microsoft.UI.Xaml.Controls.Primitives.h"
#include "winrt/Microsoft.UI.Xaml.Media.h"
#include "winrt/Microsoft.UI.Xaml.XamlTypeInfo.h"
#include "winrt/Windows.ApplicationModel.h"
#include "winrt/Windows.ApplicationModel.Activation.h"
#include "winrt/Windows.ApplicationModel.Resources.h"
#include "winrt/Windows.Foundation.h"
#include "winrt/Windows.Foundation.Collections.h"
#include "winrt/Windows.Globalization.NumberFormatting.h"
#include "winrt/Windows.Media.h"
#include "winrt/Windows.Media.Audio.h"
#include "winrt/Windows.Media.MediaProperties.h"
#include "winrt/Windows.Media.Render.h"
#include "winrt/Windows.Storage.h"
#include "winrt/Windows.Storage.Pickers.h"
#include "winrt/Windows.Storage.Streams.h"
#include "winrt/Windows.System.Threading.h"
#include "winrt/Windows.UI.Core.h"
#include "winrt/Windows.UI.Xaml.h"
#include "winrt/Windows.UI.Xaml.Controls.h"
#include "winrt/Windows.UI.Xaml.Controls.Primitives.h"
#include "winrt/Windows.UI.Xaml.Data.h"
#include "winrt/Windows.UI.Xaml.Documents.h"
#include "winrt/Windows.UI.Xaml.Interop.h"
#include "winrt/Windows.UI.Xaml.Input.h"
#include "winrt/Windows.UI.Xaml.Markup.h"
#include "winrt/Windows.UI.Xaml.Navigation.h"

#include <algorithm>
#include <array>
#include <chrono>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <mutex>
#include <optional>
#include <string>
#include <string_view>
#include <stdexcept>
#include <thread>
#include <sstream>

#include <gsl/gsl>


text:
/*
 * Copyright 2012-2016 M. Andersen and L. Vandenberghe.
 * Copyright 2010-2011 L. Vandenberghe.
 * Copyright 2004-2009 J. Dahl and L. Vandenberghe.
 *
 * This file is part of CVXOPT.
 *
 * CVXOPT 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.
 *
 * CVXOPT 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, see <http://www.gnu.org/licenses/>.
 */

#include "cvxopt.h"

#include <complex.h>

#include "misc.h"

#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#include <time.h>

PyDoc_STRVAR(gsl__doc__,"Random Module.");

static unsigned long seed = 0;
static const gsl_rng_type *rng_type;
static gsl_rng *rng;

static char doc_getseed[] =
  "Returns the seed value for the random number generator.\n\n"
  "getseed()";

static PyObject * getseed(PyObject *self)
{
  return Py_BuildValue("l",seed);
}

static char doc_setseed[] =
  "Sets the seed value for the random number generator.\n\n"
  "setseed(value = 0)\n\n"
  "ARGUMENTS\n"
  "value     integer seed. If the value is 0, then the system clock\n"
  "          measured in seconds is used instead";

static PyObject * setseed(PyObject *self, PyObject *args)
{
  unsigned long seed_ = 0;
  time_t seconds;

  if (!PyArg_ParseTuple(args, "|l", &seed_))
    return NULL;

  if (!seed_) {
    time(&seconds);
    seed = (unsigned long)seconds;
  }
  else seed = seed_;

  return Py_BuildValue("");
}

static char doc_normal[] =
  "Randomly generates a matrix with normally distributed entries.\n\n"
  "normal(nrows, ncols=1, mean=0, std=1)\n\n"
  "PURPOSE\n"
  "Returns a matrix with typecode 'd' and size nrows by ncols, with\n"
  "its entries randomly generated from a normal distribution with mean\n"
  "m and standard deviation std.\n\n"
  "ARGUMENTS\n"
  "nrows     number of rows\n\n"
  "ncols     number of columns\n\n"
  "mean      approximate mean of the distribution\n\n"
  "std       standard deviation of the distribution";
static PyObject *
normal(PyObject *self, PyObject *args, PyObject *kwrds)
{
  matrix *obj;
  int i, nrows, ncols = 1;
  double m = 0, s = 1;
  char *kwlist[] = {"nrows", "ncols", "mean", "std",  NULL};

  if (!PyArg_ParseTupleAndKeywords(args, kwrds, "i|idd", kwlist,
	  &nrows, &ncols, &m, &s)) return NULL;

  if (s < 0.0) PY_ERR(PyExc_ValueError, "std must be non-negative");

  if ((nrows<0) || (ncols<0)) {
    PyErr_SetString(PyExc_TypeError, "dimensions must be non-negative");
    return NULL;
  }

  if (!(obj = Matrix_New(nrows, ncols, DOUBLE)))
    return PyErr_NoMemory();

  gsl_rng_env_setup();
  rng_type = gsl_rng_default;
  rng = gsl_rng_alloc (rng_type);
  gsl_rng_set(rng, seed);

  for (i = 0; i < nrows*ncols; i++)
    MAT_BUFD(obj)[i] = gsl_ran_gaussian (rng, s) + m;

  seed = gsl_rng_get (rng);
  gsl_rng_free(rng);

  return (PyObject *)obj;
}

static char doc_uniform[] =
  "Randomly generates a matrix with uniformly distributed entries.\n\n"
  "uniform(nrows, ncols=1, a=0, b=1)\n\n"
  "PURPOSE\n"
  "Returns a matrix with typecode 'd' and size nrows by ncols, with\n"
  "its entries randomly generated from a uniform distribution on the\n"
  "interval (a,b).\n\n"
  "ARGUMENTS\n"
  "nrows     number of rows\n\n"
  "ncols     number of columns\n\n"
  "a         lower bound\n\n"
  "b         upper bound";

static PyObject *
uniform(PyObject *self, PyObject *args, PyObject *kwrds)
{
  matrix *obj;
  int i, nrows, ncols = 1;
  double a = 0, b = 1;

  char *kwlist[] = {"nrows", "ncols", "a", "b", NULL};

  if (!PyArg_ParseTupleAndKeywords(args, kwrds, "i|idd", kwlist,
	  &nrows, &ncols, &a, &b)) return NULL;

  if (a>b) PY_ERR(PyExc_ValueError, "a must be less than b");

  if ((nrows<0) || (ncols<0))
    PY_ERR_TYPE("dimensions must be non-negative");

  if (!(obj = (matrix *)Matrix_New(nrows, ncols, DOUBLE)))
    return PyErr_NoMemory();

  gsl_rng_env_setup();
  rng_type = gsl_rng_default;
  rng = gsl_rng_alloc (rng_type);
  gsl_rng_set(rng, seed);

  for (i= 0; i < nrows*ncols; i++)
    MAT_BUFD(obj)[i] = gsl_ran_flat (rng, a, b);

  seed = gsl_rng_get (rng);
  gsl_rng_free(rng);

  return (PyObject *)obj;
}

static PyMethodDef gsl_functions[] = {
{"getseed", (PyCFunction)getseed, METH_VARARGS|METH_KEYWORDS, doc_getseed},
{"setseed", (PyCFunction)setseed, METH_VARARGS|METH_KEYWORDS, doc_setseed},
{"normal", (PyCFunction)normal, METH_VARARGS|METH_KEYWORDS, doc_normal},
{"uniform", (PyCFunction)uniform, METH_VARARGS|METH_KEYWORDS, doc_uniform},
{NULL}  /* Sentinel */
};

#if PY_MAJOR_VERSION >= 3

static PyModuleDef gsl_module = {
    PyModuleDef_HEAD_INIT,
    "gsl",
    gsl__doc__,
    -1,
    gsl_functions,
    NULL, NULL, NULL, NULL
};

PyMODINIT_FUNC PyInit_gsl(void)
{
  PyObject *m;
  if (!(m = PyModule_Create(&gsl_module))) return NULL;
  if (import_cvxopt() < 0) return NULL;
  return m;
}

#else

PyMODINIT_FUNC initgsl(void)
{
  PyObject *m;
  m = Py_InitModule3("cvxopt.gsl", gsl_functions, gsl__doc__);
  if (import_cvxopt() < 0) return;
}

#endif

max_stars_repo_name: Hap-Hugh/quicksel

max_issues_repo_name: Hap-Hugh/quicksel

max_forks_repo_name: Hap-Hugh/quicksel