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_emu.c"

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

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

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

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

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

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

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

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

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

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

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

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

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

  if (cosmology.Omega_m < 0.04 || cosmology.Omega_m > 0.7) return 0;
  if (cosmology.omb < 0.04 || cosmology.omb > 0.055) return 0;
  if (cosmology.sigma_8 < 0.4 || cosmology.sigma_8 > 1.2) return 0;
  if (cosmology.n_spec < 0.84 || cosmology.n_spec > 1.06) return 0;
  if (cosmology.w0 < -2.1 || cosmology.w0 > -0.0) return 0;
  if (cosmology.wa < -2.6 || cosmology.wa > 2.6) return 0;
  if (cosmology.h0 < 0.4 || cosmology.h0 > 0.9) return 0;

  return 1;
}

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

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

  for (i=0;i<tomo.shear_Nbin; i++){
    nuisance.sigma_zphot_shear[i]=SPS1;
    if (nuisance.sigma_zphot_shear[i]<0.000001) 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 < -1.0 || nuisance.beta_ia > 3.0) return 0;
  if (nuisance.eta_ia < -3.0 || nuisance.eta_ia> 3.0) return 0;
  if (nuisance.eta_ia_highz < -1.0 || nuisance.eta_ia_highz> 1.0) return 0;
  // if(like.IA!=0){
  //  if (check_LF()) return 0;
  // }
return 1;
}

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

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

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

return 1;
}

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

double log_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)
{
  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;
  }

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

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

  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));
    // }
  }
  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)
{

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

  if(ell==0){
    pred= create_double_vector(0, like.Ndata-1);
    ell= create_double_vector(0, like.Ncl-1);
    darg=(log(like.lmax)-log(like.lmin))/like.Ncl;
    for (l=0;l<like.Ncl;l++){
      ell[l]=exp(log(like.lmin)+(l+0.5)*darg);
    }
    ell_Cluster= create_double_vector(0, Cluster.lbin-1);
    darg=(log(Cluster.l_max)-log(Cluster.l_min))/Cluster.lbin;
    for (l=0;l<Cluster.lbin;l++){
      ell_Cluster[l]=exp(log(Cluster.l_min)+(l+0.5)*darg);
    }
  }
// for (l=0;l<like.Ncl;l++){
//   printf("%d %le\n",i,ell[l]);
// }

  set_cosmology_params(OMM,S8,NS,W0,WA,OMB,H0,MGSigma,MGmu);
  set_nuisance_shear_calib(M1,M2,M3,M4,M5,M6,M7,M8,M9,M10);
  set_nuisance_shear_photoz(SP1,SP2,SP3,SP4,SP5,SP6,SP7,SP8,SP9,SP10,SPS1);
  set_nuisance_clustering_photoz(CP1,CP2,CP3,CP4,CP5,CP6,CP7,CP8,CP9,CP10,CPS1);
  set_nuisance_ia(A_ia,beta_ia,eta_ia,eta_ia_highz,LF_alpha,LF_P,LF_Q,LF_red_alpha,LF_red_P,LF_red_Q);
  set_nuisance_gbias(B1,B2,B3,B4,B5,B6,B7,B8,B9,B10);

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

  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 ic, input_nuisance_params in)
{
  compute_data_vector(details, ic.omega_m, ic.sigma_8, ic.n_s, ic.w0, ic.wa, ic.omega_b, ic.h0, ic.MGSigma, ic.MGmu,
    in.bias[0], in.bias[1], in.bias[2], in.bias[3],in.bias[4], in.bias[5], in.bias[6], in.bias[7],in.bias[8], in.bias[9],
    in.source_z_bias[0], in.source_z_bias[1], in.source_z_bias[2], in.source_z_bias[3], in.source_z_bias[4],
    in.source_z_bias[5], in.source_z_bias[6], in.source_z_bias[7], in.source_z_bias[8], in.source_z_bias[9],
    in.source_z_s,
    in.lens_z_bias[0], in.lens_z_bias[1], in.lens_z_bias[2], in.lens_z_bias[3], in.lens_z_bias[4],
    in.lens_z_bias[5], in.lens_z_bias[6], in.lens_z_bias[7], in.lens_z_bias[8], in.lens_z_bias[9],
    in.lens_z_s,
    in.shear_m[0], in.shear_m[1], in.shear_m[2], in.shear_m[3], in.shear_m[4],
    in.shear_m[5], in.shear_m[6], in.shear_m[7], in.shear_m[8], in.shear_m[9],
    in.A_ia, in.beta_ia, in.eta_ia, in.eta_ia_highz,
    in.lf[0], in.lf[1], in.lf[2], in.lf[3], in.lf[4], in.lf[5],
    in.m_lambda[0], in.m_lambda[1], in.m_lambda[2], in.m_lambda[3],
    in.m_lambda[4], in.m_lambda[5],in.bary[0], in.bary[1], in.bary[2]);
}

double log_like_wrapper(input_cosmo_params ic, input_nuisance_params in)
{
  double like = log_multi_like(ic.omega_m, ic.sigma_8, ic.n_s, ic.w0, ic.wa, ic.omega_b, ic.h0, ic.MGSigma, ic.MGmu,
    in.bias[0], in.bias[1], in.bias[2], in.bias[3],in.bias[4], in.bias[5], in.bias[6], in.bias[7],in.bias[8], in.bias[9],
    in.source_z_bias[0], in.source_z_bias[1], in.source_z_bias[2], in.source_z_bias[3], in.source_z_bias[4],
    in.source_z_bias[5], in.source_z_bias[6], in.source_z_bias[7], in.source_z_bias[8], in.source_z_bias[9],
    in.source_z_s,
    in.lens_z_bias[0], in.lens_z_bias[1], in.lens_z_bias[2], in.lens_z_bias[3], in.lens_z_bias[4],
    in.lens_z_bias[5], in.lens_z_bias[6], in.lens_z_bias[7], in.lens_z_bias[8], in.lens_z_bias[9],
    in.lens_z_s,
    in.shear_m[0], in.shear_m[1], in.shear_m[2], in.shear_m[3], in.shear_m[4],
    in.shear_m[5], in.shear_m[6], in.shear_m[7], in.shear_m[8], in.shear_m[9],
    in.A_ia, in.beta_ia, in.eta_ia, in.eta_ia_highz,
    in.lf[0], in.lf[1], in.lf[2], in.lf[3], in.lf[4], in.lf[5],
    in.m_lambda[0], in.m_lambda[1], in.m_lambda[2], in.m_lambda[3],
    in.m_lambda[4], in.m_lambda[5],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)
{
  clock_t begin, end;
  double time_spent;
  int i;
  char arg1[400],arg2[400],arg3[400];
/* here, do your time-consuming job */
  int sce=atoi(argv[1]);

  int N_scenarios=36;

  double area_table[36]={7623.22,14786.3,9931.47,8585.43,17681.8,15126.9,9747.99,8335.08,9533.42,18331.3,12867.8,17418.9,19783.1,12538.8,15260.0,16540.7,19636.8,11112.7,10385.5,16140.2,18920.1,17976.2,11352.0,9214.77,16910.7,11995.6,16199.8,14395.1,8133.86,13510.5,19122.3,15684.5,12014.8,14059.7,10919.3,13212.7};

  double nsource_table[36]={13.991,33.3975,17.069,28.4875,35.3643,10.3802,11.1105,29.5904,31.4608,15.7463,13.3066,9.07218,9.58719,16.3234,25.8327,10.7415,38.0412,32.8821,19.8893,27.0369,15.1711,14.2418,19.1851,26.926,22.0012,12.6553,18.6304,11.9787,36.8728,22.5265,17.4381,12.3424,10.0095,23.5979,20.8771,24.8956};

  double nlens_table[36]={22.9726 ,61.5948 ,28.7809 ,51.4354 ,65.7226 ,16.3777 ,17.6899 ,53.6984 ,57.562 ,26.266 ,21.703 ,14.0587 ,14.9667 ,27.36 ,46.0364 ,17.0253 ,71.39 ,60.5184 ,34.2283 ,48.4767 ,25.1812 ,23.44 ,32.8578 ,48.2513 ,38.3766 ,20.5029 ,31.783 ,19.2647 ,68.9094 ,39.4168 ,29.4874 ,19.9292 ,15.7162 ,41.5487 ,36.1616 ,44.148};

  char survey_designation[1][200]={"LSST"};
  char tomo_binning_source[1][200]={"source_std"};
  char tomo_binning_lens[1][200]={"LSST_gold"};

  char source_zfile[36][400]={"wl_redshift_model0_WLz01.880307e-01_WLalpha8.485694e-01.txt", "wl_redshift_model1_WLz01.731166e-01_WLalpha7.662434e-01.txt", "wl_redshift_model2_WLz01.846221e-01_WLalpha8.297540e-01.txt", "wl_redshift_model3_WLz01.758423e-01_WLalpha7.812893e-01.txt", "wl_redshift_model4_WLz01.721357e-01_WLalpha7.608290e-01.txt", "wl_redshift_model5_WLz01.931476e-01_WLalpha8.768147e-01.txt", "wl_redshift_model6_WLz01.919821e-01_WLalpha8.703814e-01.txt", "wl_redshift_model7_WLz01.751912e-01_WLalpha7.776952e-01.txt", "wl_redshift_model8_WLz01.741405e-01_WLalpha7.718958e-01.txt", "wl_redshift_model9_WLz01.860048e-01_WLalpha8.373865e-01.txt", "wl_redshift_model10_WLz01.888904e-01_WLalpha8.533151e-01.txt", "wl_redshift_model11_WLz01.954564e-01_WLalpha8.895592e-01.txt", "wl_redshift_model12_WLz01.945099e-01_WLalpha8.843347e-01.txt", "wl_redshift_model13_WLz01.853877e-01_WLalpha8.339802e-01.txt", "wl_redshift_model14_WLz01.775192e-01_WLalpha7.905458e-01.txt", "wl_redshift_model15_WLz01.925611e-01_WLalpha8.735775e-01.txt", "wl_redshift_model16_WLz01.708849e-01_WLalpha7.539247e-01.txt", "wl_redshift_model17_WLz01.733832e-01_WLalpha7.677150e-01.txt", "wl_redshift_model18_WLz01.820009e-01_WLalpha8.152851e-01.txt", "wl_redshift_model19_WLz01.767381e-01_WLalpha7.862345e-01.txt", "wl_redshift_model20_WLz01.866426e-01_WLalpha8.409073e-01.txt", "wl_redshift_model21_WLz01.877261e-01_WLalpha8.468882e-01.txt", "wl_redshift_model22_WLz01.826188e-01_WLalpha8.186960e-01.txt", "wl_redshift_model23_WLz01.768086e-01_WLalpha7.866234e-01.txt", "wl_redshift_model24_WLz01.802711e-01_WLalpha8.057362e-01.txt", "wl_redshift_model25_WLz01.897506e-01_WLalpha8.580632e-01.txt", "wl_redshift_model26_WLz01.831217e-01_WLalpha8.214720e-01.txt", "wl_redshift_model27_WLz01.906926e-01_WLalpha8.632630e-01.txt", "wl_redshift_model28_WLz01.714197e-01_WLalpha7.568766e-01.txt", "wl_redshift_model29_WLz01.798667e-01_WLalpha8.035040e-01.txt", "wl_redshift_model30_WLz01.842554e-01_WLalpha8.277299e-01.txt", "wl_redshift_model31_WLz01.901797e-01_WLalpha8.604322e-01.txt", "wl_redshift_model32_WLz01.937710e-01_WLalpha8.802561e-01.txt", "wl_redshift_model33_WLz01.790701e-01_WLalpha7.991072e-01.txt", "wl_redshift_model34_WLz01.811701e-01_WLalpha8.106987e-01.txt", "wl_redshift_model35_WLz01.781525e-01_WLalpha7.940419e-01.txt"};

  char lens_zfile[36][400]={"LSS_redshift_model0_LSSz02.629496e-01_LSSalpha9.285983e-01.txt","LSS_redshift_model1_LSSz02.852923e-01_LSSalpha8.985943e-01.txt","LSS_redshift_model2_LSSz02.664822e-01_LSSalpha9.186036e-01.txt","LSS_redshift_model3_LSSz02.798758e-01_LSSalpha9.014201e-01.txt","LSS_redshift_model4_LSSz02.873873e-01_LSSalpha8.978683e-01.txt","LSS_redshift_model5_LSSz02.593970e-01_LSSalpha9.470921e-01.txt","LSS_redshift_model6_LSSz02.600213e-01_LSSalpha9.425114e-01.txt","LSS_redshift_model7_LSSz02.811155e-01_LSSalpha9.006370e-01.txt","LSS_redshift_model8_LSSz02.831875e-01_LSSalpha8.995165e-01.txt","LSS_redshift_model9_LSSz02.649368e-01_LSSalpha9.224338e-01.txt","LSS_redshift_model10_LSSz02.622058e-01_LSSalpha9.314128e-01.txt","LSS_redshift_model11_LSSz02.584820e-01_LSSalpha9.568082e-01.txt","LSS_redshift_model12_LSSz02.588053e-01_LSSalpha9.527220e-01.txt","LSS_redshift_model13_LSSz02.656075e-01_LSSalpha9.206866e-01.txt","LSS_redshift_model14_LSSz02.768397e-01_LSSalpha9.037492e-01.txt","LSS_redshift_model15_LSSz02.596975e-01_LSSalpha9.447600e-01.txt","LSS_redshift_model16_LSSz02.901709e-01_LSSalpha8.971658e-01.txt","LSS_redshift_model17_LSSz02.847362e-01_LSSalpha8.988183e-01.txt","LSS_redshift_model18_LSSz02.698330e-01_LSSalpha9.121821e-01.txt","LSS_redshift_model19_LSSz02.782257e-01_LSSalpha9.026084e-01.txt","LSS_redshift_model20_LSSz02.642756e-01_LSSalpha9.243038e-01.txt","LSS_redshift_model21_LSSz02.632273e-01_LSSalpha9.276296e-01.txt","LSS_redshift_model22_LSSz02.689934e-01_LSSalpha9.135968e-01.txt","LSS_redshift_model23_LSSz02.780987e-01_LSSalpha9.027073e-01.txt","LSS_redshift_model24_LSSz02.723464e-01_LSSalpha9.085463e-01.txt","LSS_redshift_model25_LSSz02.615210e-01_LSSalpha9.343470e-01.txt","LSS_redshift_model26_LSSz02.683327e-01_LSSalpha9.147934e-01.txt","LSS_redshift_model27_LSSz02.608392e-01_LSSalpha9.376962e-01.txt","LSS_redshift_model28_LSSz02.889655e-01_LSSalpha8.974355e-01.txt","LSS_redshift_model29_LSSz02.729686e-01_LSSalpha9.077654e-01.txt","LSS_redshift_model30_LSSz02.669178e-01_LSSalpha9.176391e-01.txt","LSS_redshift_model31_LSSz02.612016e-01_LSSalpha9.358553e-01.txt","LSS_redshift_model32_LSSz02.591077e-01_LSSalpha9.496317e-01.txt","LSS_redshift_model33_LSSz02.742326e-01_LSSalpha9.063038e-01.txt","LSS_redshift_model34_LSSz02.710103e-01_LSSalpha9.103760e-01.txt","LSS_redshift_model35_LSSz02.757517e-01_LSSalpha9.047459e-01.txt"};

double shear_prior[36]={0.00891915,0.0104498 ,0.0145972 ,0.0191916 ,0.00450246 ,0.00567828 ,0.00294841 ,0.00530922 ,0.0118632 ,0.0151849 ,0.00410151 ,0.0170622 ,0.0197331 ,0.0106615 ,0.0124445 ,0.00994507 ,0.0136251 ,0.0143491 ,0.0164314 ,0.016962 ,0.0186608 ,0.00945903 ,0.0113246 ,0.0155225 ,0.00800846 ,0.00732104 ,0.00649453 ,0.00243976 ,0.0125932 ,0.0182587 ,0.00335859 ,0.00682287 ,0.0177269 ,0.0035219 ,0.00773304 ,0.0134886};

double delta_z_prior[36]={0.0032537,0.00135316,0.00168787,0.00215043,0.00406031,0.00222358,0.00334993,0.00255186,0.00266499,0.00159226,0.00183664,0.00384965,0.00427765,0.00314377,0.00456113,0.00347868,0.00487938,0.00418152,0.00469911,0.00367598,0.0028009,0.00234161,0.00194964,0.00200982,0.00122739,0.00310886,0.00275168,0.00492736,0.00437241,0.00113931,0.00104864,0.00292328,0.00452082,0.00394114,0.00150756,0.003613};

double sigma_z[36]={0.0849973 ,0.0986032 ,0.0875521 ,0.0968222 ,0.0225239 ,0.0718278 ,0.0733675 ,0.0385274 ,0.0425549 ,0.0605867 ,0.0178555 ,0.0853407 ,0.0124119 ,0.0531027 ,0.0304032 ,0.0503145 ,0.0132213 ,0.0941765 ,0.0416444 ,0.0668198 ,0.063227 ,0.0291332 ,0.0481633 ,0.0595606 ,0.0818742 ,0.0472518 ,0.0270185 ,0.0767401 ,0.0219945 ,0.0902663 ,0.0779705 ,0.0337666 ,0.0362358 ,0.0692429 ,0.0558841 ,0.0150457};

double sigma_z_prior[36]={0.00331909,0.00529541,0.00478151,0.00437497,0.00443062,0.00486333,0.00467423,0.0036723,0.00426963,0.00515357,0.0054553,0.00310132,0.00305971,0.00406327,0.00594293,0.00348709,0.00562526,0.00396025,0.00540537,0.00500447,0.00318595,0.00460592,0.00412137,0.00336418,0.00524988,0.00390092,0.00498349,0.0056667,0.0036384,0.00455861,0.00554822,0.00381061,0.0057615,0.00357705,0.00590572,0.00422393};

  init_cosmo_runmode("halofit");
  init_bary("dmo");
  init_binning_fourier(15,20.0,3000.0,3000.0,21.0,10,10);
  init_priors(shear_prior[sce],sigma_z[sce],delta_z_prior[sce],sigma_z_prior[sce],sigma_z[sce]*0.6,delta_z_prior[sce],sigma_z_prior[sce],3.0,1.2,3.8,2.0,16.0,5.0,0.8);
  init_survey(survey_designation[0],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[0],tomo_binning_lens[0]);
  init_IA("NLA_HF","GAMA");
  init_probes("3x2pt");

  //init_Pdelta("emu",0.8,0.35);
  // init_Pdelta("linear",0.8,0.35);

  // for (i =0; i< 10; i++){
  //   save_zdistr_sources(i);
  // }
  // for (i =0; i< 10; i++){
  //   save_zdistr_lenses(i);
  // }
//   double Omega;
//   char filename[300];
// for (i=0;i<200; i++){
//   Omega=0.05+i*0.002;
//   sprintf(filename,"test_fid_%d", i);
//compute_data_vector(filename,Omega,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);
// }

sprintf(arg3,"Y10_area%le",area_table[sce]);
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_z[sce],0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,sigma_z[sce]*0.6,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0,0.0,0.0,0.0);

// compute_data_vector("mu1_Sigma0",0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,1.,1.35,1.5,1.65,1.8,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,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,1.72+log(1.e+14*0.7),1.08,0.0,0.25,0.9,0.9,0.9,0.9);
  // compute_data_vector("mu1_Sigma1",0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,1.,1.,1.35,1.5,1.65,1.8,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.01,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,1.72+log(1.e+14*0.7),1.08,0.0,0.25,0.9,0.9,0.9,0.9);

  // init_data_inv("cov/WFIRST_3x2pt_clusterN_clusterWL_inv","datav/WFIRST_all_2pt_clusterN_clusterWL_fid");

  // begin = clock();
  // log_multi_like(0.3156,0.831,0.9645,-1.,0.,0.0491685,0.6727,0.,0.,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.05,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,5.92,1.1,-0.47,0.0,0.0,0.0,0.0,0.0,0.0,0.0,3.207,0.993,0.0,0.456,0.0,0.0);
  // // printf("knonlin %le\n",nonlinear_scale_computation(1.0));
  // // printf("knonlin %le\n",nonlinear_scale_computation(0.5));
  // end = clock();
  // time_spent = (double)(end - begin) / CLOCKS_PER_SEC;
  // printf("timespent %le\n",time_spent);

  return 0;
}


text:
/**
\ * External definitions to parametrize spectrum traces
 * in aXe grism exposures
 */

#ifndef _SPC_TRACE_FUNCTIONS_H

#define _SPC_TRACE_FUNCTIONS_H

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

#include <gsl/gsl_vector.h>
#include <gsl/gsl_integration.h>

#include "aXe_errors.h"

/* public */

/** A function for parametrizing spectrum traces; deriv and/or path_len
   (i.e. its derivative or the path length from parameter 0) may be
   NULL if they are not known analytically.
   func should return the the y offset for a given pixel column x of
   the spectrum trace
*/
typedef struct
{
  int type;			/* A unique type ID for the trace function */
  double (*func) (const double, const void *const pars);
  double (*deriv) (const double, const void *const pars);
  double (*path_len) (const double, const void *const pars);
  void *data;		/* private data for trace function, should
			   be passed as a second parameter */
}
trace_func;

/* The create_xxx functions return a pointer to an allocated trace_func
	 or NULL if the allocation failed.  */
extern trace_func *
create_poly2 (const double a0, const double a1, const double a2);

extern trace_func *
create_polyN (gsl_vector *v);

extern trace_func *
vector_to_trace_poly2 (gsl_vector *v);

extern trace_func *
vector_to_trace_polyN(gsl_vector *v);

extern double
polyN_ds(double x, void *pars);

extern void
free_poly2(trace_func * func);

extern void
free_polyN(trace_func * func);

#endif /* ! _SPC_TRACE_FUNCTIONS_H */

ext: h

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/****************************************************************************
 *                                                                          *
 *  Author : lukasz.iwaszkiewicz@gmail.com                                  *
 *  ~~~~~~~~                                                                *
 *  License : see COPYING file for details.                                 *
 *  ~~~~~~~~~                                                               *
 ****************************************************************************/

#pragma once
#include <etl/cstring.h>
// #include <fmt/core.h>
// #include <fmt/format.h>
#include <gsl/gsl>

namespace logging {

/// Line size including tick number and task name.
static constexpr size_t MAX_LINE_SIZE = 256;

/// Have to be initialized.
void init ();

/**
 * Add a log. Do not call from an ISR.
 */
bool log (gsl::czstring<> str);

// template <typename Format, typename... Args> bool log2 (Format &&format, Args &&... args)
// {
//         etl::string<MAX_LINE_SIZE> buf;
//         fmt::format_to (std::back_inserter (buf), std::forward<Format> (format), std::forward<Args> (args)...);
//         log (buf.c_str ());
//         return true;
// }

} // namespace logging

max_stars_repo_name: iwasz/freertos-test

max_issues_repo_name: iwasz/freertos-test

max_forks_repo_name: iwasz/freertos-test


text:
#pragma once

#include "shapes3D.h"

#include <gsl/gsl>

namespace ad {

template <int N_baseEdges>
std::vector<vertex3D::Data> make_pyramid(float aBaseRadius, float aHeight)
{
    std::vector<vertex3D::Data> result;
    // There are twice the number of edges in the pyramid, each edge being two vertices
    constexpr std::size_t ElementsCount{N_baseEdges*2*2};
    result.reserve(ElementsCount);

    Vec<3> topPosition = {0.0f, 0.0f, aHeight};

    auto getBasePosition = [](int aId, float aRadius) -> Vec<3>
    {
        return {
            std::cos(2*pi<>*aId/N_baseEdges) * aRadius,
            std::sin(2*pi<>*aId/N_baseEdges) * aRadius,
            0.0f,
        };
    };

    Vec<3> position_prev = getBasePosition(-1, aBaseRadius);
    for(int i=0; i<N_baseEdges; ++i)
    {
        Vec<3> position_a = getBasePosition(i, aBaseRadius);
        Vec<3> position_b = getBasePosition(i+1, aBaseRadius);

        // The base segment
        {
            // The normal of BA vector (not AB)
            Vec<3> ba = (position_a - position_b).normalize();
            Vec<3> base_normal = {-ba.y(), ba.x(), 0.0f};

            result.emplace_back(vertex3D::Data{{position_a}, {base_normal}});
            result.emplace_back(vertex3D::Data{{position_b}, {base_normal}});
        }

        // The segment connecting base vertex "i" to the top
        {
            float teta = std::atan(aBaseRadius/aHeight); // angle between Z and a side edge
            Vec<3> b_prev = (position_prev - position_b).normalize();
            Vec<3> side_normal = {
                -b_prev.y() * std::cos(teta),
                b_prev.x() * std::cos(teta),
                std::sin(teta),
            };
            result.emplace_back(vertex3D::Data{{position_a}, {side_normal}});
            result.emplace_back(vertex3D::Data{{topPosition}, {side_normal}});
        }

        position_prev = position_a;
    }

    Expects(result.size() == ElementsCount);
    return result;
}

} // namespace ad

max_stars_repo_name: FranzPoize/shmurp

max_issues_repo_name: FranzPoize/shmurp

max_forks_repo_name: FranzPoize/shmurp


text:
#pragma once

#ifndef _NOGSL
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <vector>

using namespace std;

typedef double (* GD_F_TYPE)(const gsl_vector*, void*);
typedef void (* GD_DF_TYPE)(const gsl_vector*, void*, gsl_vector*);
typedef void (* GD_FDF_TYPE)(const gsl_vector*, void*, double*, gsl_vector*);

class GradientDescent {
	int N; // Total number of components
	gsl_vector *x;  // The current state
	gsl_vector *ds; // The step size

	gsl_multimin_fdfminimizer* minidf;  // The GSL minimizer for when derivatives needed
	const gsl_multimin_fdfminimizer_type *Tdf ;
	gsl_multimin_function_fdf myfundf;

public:
	double INIT_STEP_SIZE = 2.0;
	double TOLERANCE = 0.1;
	double PRECISION = 1e-3;
	double GRAD_PRECISION = 1e-3;
	double ITERATIONS = 4000;

	GradientDescent(int N_p): N(N_p) {
		Tdf = gsl_multimin_fdfminimizer_vector_bfgs2;
		minidf = gsl_multimin_fdfminimizer_alloc(Tdf, N);
		myfundf.n = N;
		x = gsl_vector_alloc(N);
	}

	~GradientDescent() {
		gsl_multimin_fdfminimizer_free(minidf);
		gsl_vector_free(x);
	}

	gsl_vector* getResults() {
		return minidf->x;
	}
	void init(GD_F_TYPE f, GD_DF_TYPE df, GD_FDF_TYPE fdf, void* p);
	double optimize(gsl_vector* initX, bool stopAtZero = true); // optimizes until value is 0
};

#endif

max_stars_repo_name: natebragg/sketch-backend

max_issues_repo_name: natebragg/sketch-backend

max_forks_repo_name: natebragg/sketch-backend


text:
/**
 *
 * @file example_zpotrf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example of Cholesky factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>

int check_factorization(int, PLASMA_Complex64_t*, PLASMA_Complex64_t*, int, int);

int IONE=1;
int ISEED[4] = {0,0,0,1};   /* initial seed for zlarnv() */

int main ()
{

    int cores = 2;
    int N     = 10 ;
    int LDA   = 10 ;
    int info_factorization;

    PLASMA_Complex64_t *A1   = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *A2   = (PLASMA_Complex64_t *)malloc(LDA*N*sizeof(PLASMA_Complex64_t));

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /* Plasma Initialize */
    PLASMA_Init(cores);
    printf("-- PLASMA is initialized to run on %d cores. \n",cores);

    /* Initialize A1 and A2 for Symmetric Positive Matrix */
    PLASMA_zplghe( (double)N, N, A1, LDA, 51 );
    PLASMA_zlacpy( PlasmaUpperLower, N, N, A1, LDA, A2, LDA );

    /* Plasma routines */
    PLASMA_zpotrf(PlasmaUpper, N, A2, LDA);

    /* Check the factorization */
    info_factorization = check_factorization( N, A1, A2, LDA, PlasmaUpper);

    if ( info_factorization != 0 )
       printf("-- Error in ZPOTRF example ! \n");
    else
       printf("-- Run of ZPOTRF example successful ! \n");

    free(A1); free(A2);

    PLASMA_Finalize();

    return EXIT_SUCCESS;
}

/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */

int check_factorization(int N, PLASMA_Complex64_t *A1, PLASMA_Complex64_t *A2, int LDA, int uplo)
{
    double Anorm, Rnorm;
    PLASMA_Complex64_t alpha;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    PLASMA_Complex64_t *Residual = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L1       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    PLASMA_Complex64_t *L2       = (PLASMA_Complex64_t *)malloc(N*N*sizeof(PLASMA_Complex64_t));
    double *work              = (double *)malloc(N*sizeof(double));

    memset((void*)L1, 0, N*N*sizeof(PLASMA_Complex64_t));
    memset((void*)L2, 0, N*N*sizeof(PLASMA_Complex64_t));

    alpha= 1.0;

    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
    if (uplo == PlasmaUpper){
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Residual, N, work);
    Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A1, LDA, work);

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}

ext: c


text:
/**
 *
 * Copyright 2018, Planet Labs, Inc.
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include <stdio.h>

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_multifit.h>
#include <math.h>

/*  A wrapper around the GSL multi-variate robust regression routine.

    This routine should only ever be called from the Python wrapper function, so if the
    parameters need to be changed, only these two files have to be changed.

    The array variables are given as pointers to doubles, and are referenced as
    1-d arrays. However, the indexes are calculated as though they were multi-dimensional
    arrays, and the dimensions are described below.

    Input Variables
    ***************

    numParams is the number of parameters to be fitted, also equal
        to the number of independent variables
    numImages is the number of images in the stack, and corresponds to the number
        of points in each fit (before removing nulls)
    numRows is the number of rows in the image
    numCols is the number of columns in the image

    The x parameter contains the values of all the independent variables.
    It comes in one of two flavours, depending on the value of the
    perPixelX parameter. If perPixelX is 1, then the x parameter is a 4-dimensional
    array of doubles, with shape
        (numParams, numImages, numRows, numCols)
    otherwise it is a 2-dimensional array with shape
        (numParams, numImages)
    and it is assumed that the independent variables are constant over all pixels.
    This latter is the most likely case.

    The Y variable is notionally a 3-dimensional array of doubles.
    The shape should be
        (numImages, numRows, numCols)
    It contains the values of the dependent variable.

    The nullVal parameter is a scalar double value. Any occurrence of this value in the
    y array will exclude that point from the fit.

    Output Variables
    ****************

    The following variables are calculated within this routine, and passed back via the
    pointers in the parameter list.

    The c array is a 3-dimensional array of doubles. It corresponds to an image
    stack of the fitted coefficients, and its shape is
        (numParams, numRows, numCols)

    The adj_Rsqrd array stores the adjusted R^2 coefficient of determination
        statistic.

*/
void wrap_gsl_multifit_robust(double *x, double *y, double *c, double *adj_Rsqrd,
        int *numIter, double *rmse, int method, int perPixelX,
        int numRows, int numCols, int numImages, int numParams,
        int numRowsX, int numColsX, double nullVal) {
    int row, col, img, param, n, xNdx, yNdx, pixNdx;
    gsl_matrix *gslX, *gslCov;
    gsl_vector *gslY, *gslC;
    gsl_multifit_robust_workspace *workspace;
    gsl_multifit_robust_type *regressionType;
    gsl_multifit_robust_stats stats;
    int gslErrorCode;

    /* Turn off the default error handler, which aborts at the first error. */
    gsl_set_error_handler_off();

    /* Translate the integer type given into the corresponding GSL pointer */
    switch (method) {
        case 1: regressionType = gsl_multifit_robust_bisquare; break;
        case 2: regressionType = gsl_multifit_robust_cauchy; break;
        case 3: regressionType = gsl_multifit_robust_fair; break;
        case 4: regressionType = gsl_multifit_robust_huber; break;
        case 5: regressionType = gsl_multifit_robust_ols; break;
        case 6: regressionType = gsl_multifit_robust_welsch; break;
    }

    /* These structures can be allocated outside the pixel loop */
    gslCov = gsl_matrix_calloc(numParams, numParams);
    gslC = gsl_vector_calloc(numParams);

    /* Loop over all pixels */
    for (row=0; row<numRows; row++) {
        for (col=0; col<numCols; col++) {
            /* Count how many non-null y values we have. */
            n = 0;
            for (img=0; img<numImages; img++) {
                if (y[img*numRows*numCols+row*numCols+col] != nullVal) n++;
            }

            /* Allocate various structures, now we know how many non-nulls */
            if (n >= numParams) {
                workspace = gsl_multifit_robust_alloc(regressionType, n, numParams);
                gslX = gsl_matrix_calloc(n, numParams);
                gslY = gsl_vector_calloc(n);

                /* Copy the data from this pixel into the relevant GSL structures. Note
                   that we skip over null values, based on nulls in the y variable.
                */
                n = 0;
                for (img=0; img<numImages; img++) {
                    yNdx = img*numRows*numCols+row*numCols+col;
                    if (y[yNdx] != nullVal) {
                        gslY->data[n*gslY->stride] = y[yNdx];
                        for (param=0; param<numParams; param++) {
                            if (perPixelX == 0) {
                                xNdx = param * numImages + img;
                            } else {
                                xNdx = param * numRowsX * numColsX * numImages +
                                    img * numRowsX * numColsX + row * numColsX + col;
                            }
                            gslX->data[n*gslX->tda + param] = x[xNdx];
                        }
                        n++;
                    }
                }

                /* Do the regression fit */
                gslErrorCode = gsl_multifit_robust(gslX, gslY, gslC, gslCov, workspace);

                if (gslErrorCode == 0) {
                    /* Copy the coefficients back into the image stack of coefficients */
                    for (param=0; param<numParams; param++) {
                        c[param*numRows*numCols + row*numCols + col] = gslC->data[param*gslC->stride];
                    }

                    /* Copy some useful statistics into their arrays */
                    stats = gsl_multifit_robust_statistics(workspace);
                    pixNdx = row*numCols + col;
                    adj_Rsqrd[pixNdx] = stats.adj_Rsq;
                    numIter[pixNdx] = stats.numit;
                    rmse[pixNdx] = stats.rmse;
                }

                /* Free per-pixel structures */
                gsl_matrix_free(gslX);
                gsl_vector_free(gslY);
                gsl_multifit_robust_free(workspace);
            }
        }
    }

    /* Free the other structures */
    gsl_vector_free(gslC);
    gsl_matrix_free(gslCov);
}

max_stars_repo_name: planetlabs/planet-tmask

max_issues_repo_name: planetlabs/planet-tmask

max_forks_repo_name: planetlabs/planet-tmask


text:
/*
 * optimize_tnc.h
 *
 *  Created on: Feb 9, 2010
 *      Author: smitty
 */

#ifndef _OPTIMIZE_STATE_RECONSTRUCTOR_NLOPT_H_
#define _OPTIMIZE_STATE_RECONSTRUCTOR_NLOPT_H_

#include <nlopt.h>

#include "state_reconstructor.h"
#include "rate_model.h"

#include <armadillo>
using namespace arma;

void optimize_sr_nlopt(RateModel * _rm,StateReconstructor * _sr, mat * _free_mask, int _nfree);

#endif /* _OPTIMIZE_STATE_RECONSTRUCTOR_NLOPT_H_ */


text:
#include <gsl/gsl_roots.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_errno.h>
#include "eos_rho.h"
#include "root_func.h"
#include "global.h"

// "reverse" EOS. input: P, T. output: rho
double
eos_rho (double pres, void *params)
{
  struct param *eos_rho_params = (struct param *) params;
  int status;
  double rho_low, rho_hi, root;
  // root-finding stuff for inverting EOS
  const gsl_root_fsolver_type *T;
  gsl_root_fsolver *s = NULL;
  gsl_function F;
  // set up root-finding machinery for inverting EOS
  F.function = &root_func;
  F.params = params;
  T = gsl_root_fsolver_brent;
  s = gsl_root_fsolver_alloc (T);
  status = gsl_root_fsolver_set (s, &F, rho_root_min, rho_root_max);

  tmp_pres = pres;
  do
    {
      status = gsl_root_fsolver_iterate (s);
      root = gsl_root_fsolver_root (s);
      rho_low = gsl_root_fsolver_x_lower (s);
      rho_hi = gsl_root_fsolver_x_upper (s);
      status = gsl_root_test_interval (rho_low, rho_hi, 0, 0.001);
    }
  while (status == GSL_CONTINUE);

  gsl_root_fsolver_free (s);
  return (root);
}


text:
/*
 * File:   Random.h
 * Author: nguyentran
 *
 * Created on May 27, 2013, 10:46 AM
 */

#ifndef RANDOM_H
#define    RANDOM_H

#include <gsl/gsl_rng.h>
#include "PropertyMacro.h"
#include "Strategies/AdaptiveCyclingStrategy.h"

class Model;

class Random {
 DISALLOW_COPY_AND_ASSIGN(Random)

 DISALLOW_MOVE(Random)

 VIRTUAL_PROPERTY(unsigned long, seed)

 public:
  gsl_rng *G_RNG;

  explicit Random(gsl_rng *g_rng = nullptr);

  virtual ~Random();

  void initialize(const unsigned long &seed = 0);

  void release() const;

  virtual int random_poisson(const double &poisson_mean);

  virtual unsigned long random_uniform(unsigned long range);

  virtual unsigned long random_uniform_int(const unsigned long &from, const unsigned long &to);

  virtual double random_uniform_double(const double &from, const double &to);

  /*
   * This function will return a random number in [0,1)
   */
  virtual double random_uniform();

  virtual double random_normal(const double &mean, const double &sd);

  virtual double random_normal_truncated(const double &mean, const double &sd);

  virtual int random_normal(const int &mean, const int &sd);

  virtual int random_normal_truncated(const int &mean, const int &sd);

  virtual double random_beta(const double &alpha, const double &beta);

  virtual double random_gamma(const double &shape, const double &scale);

  virtual double cdf_gamma_distribution(const double &x, const double &alpha, const double &beta);

  virtual double cdf_gamma_distribution_inverse(const double &p, const double &alpha, const double &beta);

  virtual double random_flat(const double &from, const double &to);

  virtual void random_multinomial(const size_t &K, const unsigned &N, double p[], unsigned n[]);

  virtual void random_shuffle(void *base, size_t base_length, size_t size_of_type);

  virtual double cdf_standard_normal_distribution(const double &p);

  virtual int random_binomial(const double &p, const unsigned int &n);

  void shuffle(void *base, const size_t &n, const size_t &size);
};

#endif    /* RANDOM_H */


text:
/* vector/gsl_vector_complex_float.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 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.
 */

#ifndef __GSL_VECTOR_COMPLEX_FLOAT_H__
#define __GSL_VECTOR_COMPLEX_FLOAT_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_float.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_float.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  float *data;
  gsl_block_complex_float *block;
  int owner;
} gsl_vector_complex_float;

typedef struct
{
  gsl_vector_complex_float vector;
} _gsl_vector_complex_float_view;

typedef _gsl_vector_complex_float_view gsl_vector_complex_float_view;

typedef struct
{
  gsl_vector_complex_float vector;
} _gsl_vector_complex_float_const_view;

typedef const _gsl_vector_complex_float_const_view gsl_vector_complex_float_const_view;

/* Allocation */

GSL_FUN gsl_vector_complex_float *gsl_vector_complex_float_alloc (const size_t n);
GSL_FUN gsl_vector_complex_float *gsl_vector_complex_float_calloc (const size_t n);

GSL_FUN gsl_vector_complex_float *
gsl_vector_complex_float_alloc_from_block (gsl_block_complex_float * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

GSL_FUN gsl_vector_complex_float *
gsl_vector_complex_float_alloc_from_vector (gsl_vector_complex_float * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

GSL_FUN void gsl_vector_complex_float_free (gsl_vector_complex_float * v);

/* Views */

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_view_array (float *base,
                                     size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_view_array_with_stride (float *base,
                                                 size_t stride,
                                                 size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_view_array (const float *base,
                                           size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_view_array_with_stride (const float *base,
                                                       size_t stride,
                                                       size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_subvector (gsl_vector_complex_float *base,
                                         size_t i,
                                         size_t n);

GSL_FUN _gsl_vector_complex_float_view
gsl_vector_complex_float_subvector_with_stride (gsl_vector_complex_float *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_subvector (const gsl_vector_complex_float *base,
                                               size_t i,
                                               size_t n);

GSL_FUN _gsl_vector_complex_float_const_view
gsl_vector_complex_float_const_subvector_with_stride (const gsl_vector_complex_float *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

GSL_FUN _gsl_vector_float_view
gsl_vector_complex_float_real (gsl_vector_complex_float *v);

GSL_FUN _gsl_vector_float_view
gsl_vector_complex_float_imag (gsl_vector_complex_float *v);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_complex_float_const_real (const gsl_vector_complex_float *v);

GSL_FUN _gsl_vector_float_const_view
gsl_vector_complex_float_const_imag (const gsl_vector_complex_float *v);

/* Operations */

GSL_FUN void gsl_vector_complex_float_set_zero (gsl_vector_complex_float * v);
GSL_FUN void gsl_vector_complex_float_set_all (gsl_vector_complex_float * v,
                                       gsl_complex_float z);
GSL_FUN int gsl_vector_complex_float_set_basis (gsl_vector_complex_float * v, size_t i);

GSL_FUN int gsl_vector_complex_float_fread (FILE * stream,
                                    gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_fwrite (FILE * stream,
                                     const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_fscanf (FILE * stream,
                                     gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_fprintf (FILE * stream,
                                      const gsl_vector_complex_float * v,
                                      const char *format);

GSL_FUN int gsl_vector_complex_float_memcpy (gsl_vector_complex_float * dest, const gsl_vector_complex_float * src);

GSL_FUN int gsl_vector_complex_float_reverse (gsl_vector_complex_float * v);

GSL_FUN int gsl_vector_complex_float_swap (gsl_vector_complex_float * v, gsl_vector_complex_float * w);
GSL_FUN int gsl_vector_complex_float_swap_elements (gsl_vector_complex_float * v, const size_t i, const size_t j);

GSL_FUN int gsl_vector_complex_float_equal (const gsl_vector_complex_float * u,
                                    const gsl_vector_complex_float * v);

GSL_FUN int gsl_vector_complex_float_isnull (const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_ispos (const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_isneg (const gsl_vector_complex_float * v);
GSL_FUN int gsl_vector_complex_float_isnonneg (const gsl_vector_complex_float * v);

GSL_FUN int gsl_vector_complex_float_add (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_sub (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_mul (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_div (gsl_vector_complex_float * a, const gsl_vector_complex_float * b);
GSL_FUN int gsl_vector_complex_float_scale (gsl_vector_complex_float * a, const gsl_complex_float x);
GSL_FUN int gsl_vector_complex_float_add_constant (gsl_vector_complex_float * a, const gsl_complex_float x);
GSL_FUN int gsl_vector_complex_float_axpby (const gsl_complex_float alpha, const gsl_vector_complex_float * x, const gsl_complex_float beta, gsl_vector_complex_float * y);

GSL_FUN INLINE_DECL gsl_complex_float gsl_vector_complex_float_get (const gsl_vector_complex_float * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_complex_float_set (gsl_vector_complex_float * v, const size_t i, gsl_complex_float z);
GSL_FUN INLINE_DECL gsl_complex_float *gsl_vector_complex_float_ptr (gsl_vector_complex_float * v, const size_t i);
GSL_FUN INLINE_DECL const gsl_complex_float *gsl_vector_complex_float_const_ptr (const gsl_vector_complex_float * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex_float
gsl_vector_complex_float_get (const gsl_vector_complex_float * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      gsl_complex_float zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_FLOAT_AT (v, i);
}

INLINE_FUN
void
gsl_vector_complex_float_set (gsl_vector_complex_float * v,
                              const size_t i, gsl_complex_float z)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_FLOAT_AT (v, i) = z;
}

INLINE_FUN
gsl_complex_float *
gsl_vector_complex_float_ptr (gsl_vector_complex_float * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_FLOAT_AT (v, i);
}

INLINE_FUN
const gsl_complex_float *
gsl_vector_complex_float_const_ptr (const gsl_vector_complex_float * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_FLOAT_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_FLOAT_H__ */

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:
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>

#define NRAND 1000000

#define HIST_SIZE 25
#define HIST_MIN 0
#define HIST_MAX (HIST_SIZE - 1)

int hist[HIST_SIZE];

int main (void)
{
	const gsl_rng_type * T;
	gsl_rng * r;
	double mu = 2.0;

	/* create a generator chosen by the
	   environment variable GSL_RNG_TYPE */
	gsl_rng_env_setup();
	T = gsl_rng_default;
	r = gsl_rng_alloc (T);

	printf("Distribution with mean of %f\n", mu);

	/* print N random variates chosen from
	   the poisson distribution with mean
	   parameter mu */
	for (int i = 0; i < NRAND; i++) {
		/* unsigned int k = gsl_ran_poisson(r, mu) + 10; */
		unsigned int k = gsl_ran_gaussian(r, mu) + 10;
		if (HIST_MIN <= k && k <= HIST_MAX) {
			hist[k]++;
		} else {
			printf ("outside history range: %u\n", k);
		}
	}

	for (int i = 0; i < HIST_SIZE; i++) {
		printf("[%d] ", i);
		for (int j = 0; j < hist[i] / (NRAND / 100); j++) {
			putchar('*');
		}
		printf("\n");
	}

	gsl_rng_free (r);
	return 0;
}

max_stars_repo_name: Rugang/scraps

max_issues_repo_name: Rugang/scraps

max_forks_repo_name: Rugang/scraps


text:
/* linalg/condest.c
 *
 * Copyright (C) 2016 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 <string.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>

/*
 * This module contains routines for estimating the condition number
 * of matrices in the 1-norm. The algorithm is based on the paper,
 *
 * [1] N. J. Higham, "FORTRAN codes for estimating the one-norm of
 * a real or complex matrix, with applications to condition estimation",
 * ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
 */

static double condest_tri_norm1(CBLAS_UPLO_t Uplo, const gsl_matrix * A);
static int condest_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A,
                             double * rcond, gsl_vector * work);
static int condest_same_sign(const gsl_vector * x, const gsl_vector * y);
static int condest_invtriu(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);
static int condest_invtril(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params);

/*
gsl_linalg_tri_upper_rcond()
  Estimate reciprocal condition number of upper triangular matrix

Inputs: A     - upper triangular matrix, N-by-N
        rcond - (output) reciprocal condition number estimate
        work  - workspace, length 3*N

Return: success/error
*/

int
gsl_linalg_tri_upper_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work)
{
  int status = condest_tri_rcond(CblasUpper, A, rcond, work);
  return status;
}

/*
gsl_linalg_tri_lower_rcond()
  Estimate reciprocal condition number of lower triangular matrix

Inputs: A     - lower triangular matrix, N-by-N
        rcond - (output) reciprocal condition number estimate
        work  - workspace, length 3*N

Return: success/error
*/

int
gsl_linalg_tri_lower_rcond(const gsl_matrix * A, double * rcond, gsl_vector * work)
{
  int status = condest_tri_rcond(CblasLower, A, rcond, work);
  return status;
}

/*
gsl_linalg_invnorm1()
  Estimate the 1-norm of ||A^{-1}||, where A is a square
N-by-N matrix

Inputs: N        - size of matrix
        Ainvx    - pointer to function which calculates:
                   x := A^{-1} x or x := A^{-t} x
        params   - parameters to pass to Ainvx
        Ainvnorm - (output) estimate of ||A^{-1}||_1
        work     - workspace, length 3*N
*/

int
gsl_linalg_invnorm1(const size_t N,
                    int (* Ainvx)(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params),
                    void * params, double * Ainvnorm, gsl_vector * work)
{
  if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      const size_t maxit = 5;
      gsl_vector_view x = gsl_vector_subvector(work, 0, N);
      gsl_vector_view v = gsl_vector_subvector(work, N, N);
      gsl_vector_view xi = gsl_vector_subvector(work, 2*N, N);
      double gamma, gamma_old, temp;
      size_t i, k;

      for (i = 0; i < N; ++i)
        gsl_vector_set(&x.vector, i, 1.0 / (double) N);

      /* compute v = A^{-1} x */
      gsl_vector_memcpy(&v.vector, &x.vector);
      (*Ainvx)(CblasNoTrans, &v.vector, params);

      /* gamma = ||v||_1 */
      gamma = gsl_blas_dasum(&v.vector);

      /* xi = sign(v) */
      for (i = 0; i < N; ++i)
        {
          double vi = gsl_vector_get(&v.vector, i);
          gsl_vector_set(&xi.vector, i, GSL_SIGN(vi));
        }

      /* x = A^{-t} xi */
      gsl_vector_memcpy(&x.vector, &xi.vector);
      (*Ainvx)(CblasTrans, &x.vector, params);

      for (k = 0; k < maxit; ++k)
        {
          size_t j = (size_t) gsl_blas_idamax(&x.vector);

          /* v := A^{-1} e_j */
          gsl_vector_set_zero(&v.vector);
          gsl_vector_set(&v.vector, j, 1.0);
          (*Ainvx)(CblasNoTrans, &v.vector, params);

          gamma_old = gamma;
          gamma = gsl_blas_dasum(&v.vector);

          /* check for repeated sign vector (algorithm has converged) */
          if (condest_same_sign(&v.vector, &xi.vector) || (gamma < gamma_old))
            break;

          /* xi = sign(v) */
          for (i = 0; i < N; ++i)
            {
              double vi = gsl_vector_get(&v.vector, i);
              gsl_vector_set(&xi.vector, i, GSL_SIGN(vi));
            }

          /* x = A^{-t} sign(v) */
          gsl_vector_memcpy(&x.vector, &xi.vector);
          (*Ainvx)(CblasTrans, &x.vector, params);
        }

      temp = 1.0; /* (-1)^i */
      for (i = 0; i < N; ++i)
        {
          double term = 1.0 + (double) i / (N - 1.0);
          gsl_vector_set(&x.vector, i, temp * term);
          temp = -temp;
        }

      /* x := A^{-1} x */
      (*Ainvx)(CblasNoTrans, &x.vector, params);

      temp = 2.0 * gsl_blas_dasum(&x.vector) / (3.0 * N);
      if (temp > gamma)
        {
          gsl_vector_memcpy(&v.vector, &x.vector);
          gamma = temp;
        }

      *Ainvnorm = gamma;

      return GSL_SUCCESS;
    }
}

static int
condest_tri_rcond(CBLAS_UPLO_t Uplo, const gsl_matrix * A, double * rcond, gsl_vector * work)
{
  const size_t M = A->size1;
  const size_t N = A->size2;

  if (M != N)
    {
      GSL_ERROR ("matrix must be square", GSL_ENOTSQR);
    }
  else if (work->size != 3 * N)
    {
      GSL_ERROR ("work vector must have length 3*N", GSL_EBADLEN);
    }
  else
    {
      int status;
      double Anorm = condest_tri_norm1(Uplo, A); /* ||A||_1 */
      double Ainvnorm;                           /* ||A^{-1}||_1 */

      *rcond = 0.0;

      /* don't continue if matrix is singular */
      if (Anorm == 0.0)
        return GSL_SUCCESS;

      /* estimate ||A^{-1}||_1 */
      if (Uplo == CblasUpper)
        status = gsl_linalg_invnorm1(N, condest_invtriu, (void *) A, &Ainvnorm, work);
      else
        status = gsl_linalg_invnorm1(N, condest_invtril, (void *) A, &Ainvnorm, work);

      if (status)
        return status;

      if (Ainvnorm != 0.0)
        *rcond = (1.0 / Anorm) / Ainvnorm;

      return GSL_SUCCESS;
    }
}

/* calculate 1 norm of triangular matrix */
static double
condest_tri_norm1(CBLAS_UPLO_t Uplo, const gsl_matrix * A)
{
  const size_t N = A->size2;
  double max = 0.0;
  size_t i, j;

  if (Uplo == CblasUpper)
    {
      for (j = 0; j < N; ++j)
        {
          double sum = 0.0;
          for (i = 0; i <= j; ++i)
            {
              double Aij = gsl_matrix_get(A, i, j);
              sum += fabs(Aij);
            }

          max = GSL_MAX(max, sum);
        }
    }
  else
    {
      for (j = 0; j < N; ++j)
        {
          double sum = 0.0;
          for (i = j; i < N; ++i)
            {
              double Aij = gsl_matrix_get(A, i, j);
              sum += fabs(Aij);
            }

          max = GSL_MAX(max, sum);
        }
    }

  return max;
}

/* return 1 if sign(x) = sign(y), 0 otherwise */
static int
condest_same_sign(const gsl_vector * x, const gsl_vector * y)
{
  const size_t n = x->size;
  size_t i;

  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(y, i);
      if (GSL_SIGN(xi) != GSL_SIGN(yi))
        return 0;
    }

  return 1;
}

/* x := A^{-1} x, A upper triangular */
static int
condest_invtriu(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
{
  gsl_matrix * A = (gsl_matrix *) params;
  return gsl_blas_dtrsv(CblasUpper, TransA, CblasNonUnit, A, x);
}

/* x := A^{-1} x, A lower triangular */
static int
condest_invtril(CBLAS_TRANSPOSE_t TransA, gsl_vector * x, void * params)
{
  gsl_matrix * A = (gsl_matrix *) params;
  return gsl_blas_dtrsv(CblasLower, TransA, CblasNonUnit, A, x);
}

max_issues_repo_name: Brian-ning/HMNE

max_forks_repo_name: Brian-ning/HMNE


text:
#define MAXLOCI 1000000
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

/* data definition- holds sample size of each species, segsites in each species, D, length of locus */
struct hkaData{
  int na, nb, length, lArray[30], sArray[30];
  double d, sa, sb;
  char name[81];
};

struct hkaParams{
  int *lengths;
};

double harmSum(int n);
double harmSumSquare(int n);
void getData(int argc, char *argv[]);
void usage();
double calculateHKA(gsl_vector *params,gsl_vector *chiSquares, gsl_vector *expSA, gsl_vector *expD);
int chopByWhite(char *in, char *outArray[], int outSize);
int setHKAFunction(const gsl_vector *x, void *p, gsl_vector *f);
int solveSystem(gsl_vector *params);

max_stars_repo_name: kern-lab/hka

max_issues_repo_name: kern-lab/hka

max_forks_repo_name: kern-lab/hka


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

#include <config.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_ieee_utils.h>

/* A table of sign characters, 0=positive, 1=negative. We print a space
   instead of a unary + sign for compatibility with bc */

static char signs[2]={' ','-'} ;

void
gsl_ieee_fprintf_float (FILE * stream, const float * x) {
  gsl_ieee_float_rep r ;
  gsl_ieee_float_to_rep(x, &r) ;

  switch (r.type)
    {
    case GSL_IEEE_TYPE_NAN:
      fprintf(stream, "NaN") ;
      break ;
    case GSL_IEEE_TYPE_INF:
      fprintf(stream, "%cInf", signs[r.sign]) ;
      break ;
    case GSL_IEEE_TYPE_NORMAL:
      fprintf(stream, "%c1.%s*2^%d", signs[r.sign], r.mantissa, r.exponent) ;
      break ;
    case GSL_IEEE_TYPE_DENORMAL:
      fprintf(stream, "%c0.%s*2^%d", signs[r.sign], r.mantissa, r.exponent + 1) ;
      break ;
    case GSL_IEEE_TYPE_ZERO:
      fprintf(stream, "%c0", signs[r.sign]) ;
      break ;
    default:
      fprintf(stream, "[non-standard IEEE float]") ;
    }
}

void
gsl_ieee_printf_float (const float * x)
{
  gsl_ieee_fprintf_float (stdout,x);
}

void
gsl_ieee_fprintf_double (FILE * stream, const double * x) {
  gsl_ieee_double_rep r ;
  gsl_ieee_double_to_rep (x, &r) ;

  switch (r.type)
    {
    case GSL_IEEE_TYPE_NAN:
      fprintf(stream, "NaN") ;
      break ;
    case GSL_IEEE_TYPE_INF:
      fprintf(stream, "%cInf", signs[r.sign]) ;
      break ;
    case GSL_IEEE_TYPE_NORMAL:
      fprintf(stream, "%c1.%s*2^%d", signs[r.sign], r.mantissa, r.exponent) ;
      break ;
    case GSL_IEEE_TYPE_DENORMAL:
      fprintf(stream, "%c0.%s*2^%d", signs[r.sign], r.mantissa, r.exponent + 1) ;
      break ;
    case GSL_IEEE_TYPE_ZERO:
      fprintf(stream, "%c0", signs[r.sign]) ;
      break ;
    default:
      fprintf(stream, "[non-standard IEEE double]") ;
    }
}

void
gsl_ieee_printf_double (const double * x)
{
  gsl_ieee_fprintf_double (stdout,x);
}

max_stars_repo_name: ielomariala/Hex-Game


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

#include <gsl/gsl_math.h>
#include <gsl/gsl_movstat.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>

int
main(void)
{
  const size_t N = 1000;                                 /* length of time series */
  const double sigma[] = { 1.0, 5.0, 1.0, 3.0, 5.0 };    /* variances */
  const size_t N_sigma[] = { 200, 450, 600, 850, 1000 }; /* samples where variance changes */
  const size_t K = 41;                                   /* window size */
  gsl_vector *x = gsl_vector_alloc(N);
  gsl_vector *xmedian = gsl_vector_alloc(N);
  gsl_vector *xmad = gsl_vector_alloc(N);
  gsl_vector *xiqr = gsl_vector_alloc(N);
  gsl_vector *xSn = gsl_vector_alloc(N);
  gsl_vector *xQn = gsl_vector_alloc(N);
  gsl_vector *xsd = gsl_vector_alloc(N);
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);
  gsl_movstat_workspace * w = gsl_movstat_alloc(K);
  size_t idx = 0;
  size_t i;

  for (i = 0; i < N; ++i)
    {
      double gi = gsl_ran_gaussian(r, sigma[idx]);
      double u = gsl_rng_uniform(r);
      double outlier = (u < 0.01) ? 15.0*GSL_SIGN(gi) : 0.0;
      double xi = gi + outlier;

      gsl_vector_set(x, i, xi);

      if (i == N_sigma[idx] - 1)
        ++idx;
    }

  /* compute moving statistics */
  gsl_movstat_mad(GSL_MOVSTAT_END_TRUNCATE, x, xmedian, xmad, w);
  gsl_movstat_qqr(GSL_MOVSTAT_END_TRUNCATE, x, 0.25, xiqr, w);
  gsl_movstat_Sn(GSL_MOVSTAT_END_TRUNCATE, x, xSn, w);
  gsl_movstat_Qn(GSL_MOVSTAT_END_TRUNCATE, x, xQn, w);
  gsl_movstat_sd(GSL_MOVSTAT_END_TRUNCATE, x, xsd, w);

  /* scale IQR by factor to approximate standard deviation */
  gsl_vector_scale(xiqr, 0.7413);

  /* print results */
  idx = 0;
  for (i = 0; i < N; ++i)
    {
      printf("%zu %f %f %f %f %f %f %f\n",
             i,
             gsl_vector_get(x, i),
             sigma[idx],
             gsl_vector_get(xmad, i),
             gsl_vector_get(xiqr, i),
             gsl_vector_get(xSn, i),
             gsl_vector_get(xQn, i),
             gsl_vector_get(xsd, i));

      if (i == N_sigma[idx] - 1)
        ++idx;
    }

  gsl_vector_free(x);
  gsl_vector_free(xmedian);
  gsl_vector_free(xmad);
  gsl_vector_free(xiqr);
  gsl_vector_free(xSn);
  gsl_vector_free(xQn);
  gsl_vector_free(xsd);
  gsl_rng_free(r);
  gsl_movstat_free(w);

  return 0;
}

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/**
 * copyright  (C) 2004
 * the icecube collaboration
 * $Id: I3GSLRandomService.h 161127 2018-02-20 14:27:18Z kjmeagher $
 *
 * @brief An implementation of the I3RandomService interface.
 *
 * Uses the gsl library for the random numbers
 *
 * @version $Revision: 161127 $
 * @date $Date: 2018-02-20 07:27:18 -0700 (Tue, 20 Feb 2018) $
 * @author pretz
 */

#ifndef I3GSLRANDOMSERVICE_H
#define I3GSLRANDOMSERVICE_H

#include "phys-services/I3RandomService.h"

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

/**
 * This is (the state for) a shim which allows us to count the calls to the
 * RNG, but otherwise hands all work off to the real GSL implmentations.
 */
typedef struct
{
  unsigned long int seed;
  uint64_t icalls;
  uint64_t dcalls;
  gsl_rng* rng;
} gsl_rng_wrapper_state;

extern const gsl_rng_type gsl_rng_counting_wrapper;

class I3GSLRandomService : public I3RandomService{
 public:
  /**
   * default constructor
   */
  I3GSLRandomService();

  /**
   * constructor
   */
  explicit I3GSLRandomService(unsigned long int seed, bool track_state=true);

  /**
   * destructor
   */
  virtual ~I3GSLRandomService();

  /**
   * a number drawn from a binomial distribution
   */
  virtual int Binomial(int ntot, double prob);

  /**
   * A number from an Exponential distribution
   */
  virtual double Exp(double tau);

  /**
   * An integer drawn uniformly from [0,imax)
   */
  virtual unsigned int Integer(unsigned int imax);

  /**
   * An integer drawn from a Poisson distribution
   */
  virtual int Poisson(double mean);

  /**
   * A number drawn from a Poisson distribution, as a double
   */
  virtual double PoissonD(double mean);

  /**
   * a double drawn from a uniform distribution (0,x1)
   */
  virtual double Uniform(double x1 = 1);

  /**
   * a double drawn from a uniform distribution (x1,x2)
   */
  virtual double Uniform(double x1, double x2);

  /**
   * a double drawn from a Gaussian distribution with given
   * mean and standard deviation
   */
  virtual double Gaus(double mean, double stddev);

  /**
   * get all information necessary to restore the internal
   * state of the generator
   */
  virtual I3FrameObjectPtr GetState() const;

  /**
   * restore the internal state of the generator
   */
  virtual void RestoreState(I3FrameObjectConstPtr state);

 private:

  // private copy constructors and assignment
  I3GSLRandomService(const I3GSLRandomService& );
  I3GSLRandomService operator=(const I3GSLRandomService& );

  /**
   * Helper function which constructs our preferred GSL RNG.
   */
  static void construct(gsl_rng*& r);
  /**
   * Helper function which constructs a GSL RNG wrapped with the counting shim
   */
  static void construct_counted(gsl_rng*& r);

  gsl_rng* r;
  bool track_state;

  SET_LOGGER("I3GSLRandomService");

  //let gsl_wrapper_set use construct()
  friend void gsl_wrapper_set(void* vstate, unsigned long int s);
};

I3_POINTER_TYPEDEFS(I3GSLRandomService);

#endif //I3GSLRANDOMSERVICE_H


text:
/* specfunc/legendre.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman */

/* Declare private but non-local support functions
 * used in various Legendre function evaluations.
 */

#include <gsl/gsl_sf_result.h>

/* Large negative mu asymptotic
 * P^{-mu}_{-1/2 + I tau}, mu -> Inf
 * |x| < 1
 */
int
gsl_sf_conicalP_xlt1_large_neg_mu_e(double mu, double tau, double x,
                                       gsl_sf_result * result, double * ln_multiplier);

/* Large tau uniform asymptotics
 * P^{-mu}_{-1/2 + I tau}, tau -> Inf
 * 1 < x
 */
int
gsl_sf_conicalP_xgt1_neg_mu_largetau_e(const double mu, const double tau,
                                          const double x, double acosh_x,
                                          gsl_sf_result * result, double * ln_multiplier);

/* Large tau uniform asymptotics
 * P^{-mu}_{-1/2 + I tau}, tau -> Inf
 * -1 < x < 1
 */
int
gsl_sf_conicalP_xlt1_neg_mu_largetau_e(const double mu, const double tau,
                                          const double x, const double acos_x,
                                          gsl_sf_result * result, double * ln_multiplier);

/* P^{mu}_{-1/2 + I tau}
 * x->Inf
 *
 *  * This is effective to precision EPS for
 *
 *    (mu^2 + tau^2)/((1 + tau^2)^(1/2) x^2) < EPS^{1/3}
 *
 * since it goes only to a fixed order, based on the
 * representation in terms of hypegeometric functions
 * of argument 1/x^2.
 * [Zhurina+Karmazina, (3.8)]
 */
int
gsl_sf_conicalP_large_x_e(const double mu, const double tau, const double x,
                             gsl_sf_result * result, double * ln_multiplier);


text:
#ifndef __COMPLEARN_REALCOMPRESSOR_H
#define __COMPLEARN_REALCOMPRESSOR_H

#define COMPLEARN_REAL_COMPRESSOR_TYPE        (real_compressor_get_type ())
#define COMPLEARN_TYPE_REAL_COMPRESSOR        (real_compressor_get_type ())
#define COMPLEARN_REAL_COMPRESSOR(obj)        (G_TYPE_CHECK_INSTANCE_CAST ((obj), COMPLEARN_REAL_COMPRESSOR_TYPE, CompLearnRealCompressor))
#define IS_COMPLEARN_REAL_COMPRESSOR(obj)        (G_TYPE_CHECK_INSTANCE_TYPE ((obj), COMPLEARN_REAL_COMPRESSOR_TYPE))
#define COMPLEARN_REAL_COMPRESSOR_GET_INTERFACE(inst) (G_TYPE_INSTANCE_GET_INTERFACE((inst), COMPLEARN_REAL_COMPRESSOR_TYPE, CompLearnRealCompressorIface))

#include <glib.h>
#include <glib-object.h>
#include <glib/garray.h>
#include <gsl/gsl_matrix.h>

#define COMPLEARN_ERROR 1
#define COMPLEARN_ERROR_NO_COMPRESSOR_SET 100

typedef struct _CompLearnRealCompressor CompLearnRealCompressor;
typedef struct _CompLearnRealCompressorIface CompLearnRealCompressorIface;

struct _CompLearnRealCompressorIface {
  GTypeInterface parent;

  GString *(*compress)(CompLearnRealCompressor *self, const GString *input);
  GString *(*decompress)(CompLearnRealCompressor *self, const GString *input);
  GString *(*blurb)(CompLearnRealCompressor *self);
  GString *(*canonical_extension)(CompLearnRealCompressor *self);
  GString *(*name)(CompLearnRealCompressor *self);
  GString *(*compressor_version)(CompLearnRealCompressor *self);
  GString *(*binding_version)(CompLearnRealCompressor *self);
  gboolean (*is_threadsafe)(CompLearnRealCompressor *self);
  gboolean (*is_compressible)(CompLearnRealCompressor *self, const GString *input);
  gboolean (*is_decompressible)(CompLearnRealCompressor *self, const GString *input);
  gboolean (*is_just_size)(CompLearnRealCompressor *self);
  gboolean (*is_hash_function)(CompLearnRealCompressor *self);
  GString *(*hash)(CompLearnRealCompressor *self, const GString *input);
  gboolean (*is_operational)(CompLearnRealCompressor *self);
  gboolean (*is_private_property)(CompLearnRealCompressor *self, const char *propname);
  gdouble (*compressed_size)(CompLearnRealCompressor *self, const GString *input);
  guint64 (*window_size)(CompLearnRealCompressor *self);
  CompLearnRealCompressor *(*clone)(CompLearnRealCompressor *self);
};

GType real_compressor_get_type(void);

GString *real_compressor_compress(CompLearnRealCompressor *self,const GString *input);
GString *real_compressor_hash(CompLearnRealCompressor *self,const GString *input);
GString *real_compressor_decompress(CompLearnRealCompressor *self,const GString *input);
gdouble real_compressor_compressed_size(CompLearnRealCompressor *self,const GString *input);
GString *real_compressor_blurb(CompLearnRealCompressor *self);
GString *real_compressor_name(CompLearnRealCompressor *self);
GString *real_compressor_compressor_version(CompLearnRealCompressor *self);
GString *real_compressor_binding_version(CompLearnRealCompressor *self);
gboolean real_compressor_is_compressible(CompLearnRealCompressor *self, const GString *input);
gboolean real_compressor_is_decompressible(CompLearnRealCompressor *self, const GString *input);
gboolean real_compressor_is_private_property(CompLearnRealCompressor *self, const char *input);
guint64 real_compressor_window_size(CompLearnRealCompressor *self);
gboolean real_compressor_is_threadsafe(CompLearnRealCompressor *self);
gboolean real_compressor_is_just_size(CompLearnRealCompressor *self);
gboolean real_compressor_is_hash_function(CompLearnRealCompressor *self);
GString *real_compressor_canonical_extension(CompLearnRealCompressor *rc);
gboolean real_compressor_is_operational(CompLearnRealCompressor *self);
CompLearnRealCompressor *real_compressor_clone(CompLearnRealCompressor *self);

#define SET_DEFAULT_PROPS(groupname, clt, mobj) \
  do { \
  GParamSpec **gps, **cur; \
  g_assert(mobj != NULL); \
    if (complearn_environment_get_nameable(groupname) == NULL) { \
      complearn_environment_register_nameable(groupname, G_OBJECT(mobj)); \
      break; \
    } \
  gps =g_object_class_list_properties(G_OBJECT_CLASS(clt(mobj)), NULL); \
  for (cur = gps; *cur; cur += 1) { \
    GValue v = {0,}; \
    g_value_init(&v, (*cur)->value_type); \
    g_param_value_set_default(*cur, &v); \
    g_object_set_property(G_OBJECT(mobj), (*cur)->name, &v); \
    complearn_environment_register_property(G_OBJECT(mobj), (*cur), &v); \
  } } while(0)

#define SET_DEFAULT_COMPRESSOR_PROPS(groupname, clt, mobj) \
  do { \
      SET_DEFAULT_PROPS(groupname, clt, mobj); \
    if (complearn_environment_get_nameable(groupname) == NULL) { \
      complearn_environment_register_compressor(mobj); \
    } \
  } while (0);

#endif

#define G_LOG_LEVEL_NOTICE G_LOG_LEVEL_USER_SHIFT
#define g_notice(...) g_log(G_LOG_DOMAIN, G_LOG_LEVEL_NOTICE, __VA_ARGS__)

void real_compressor_interface_init (gpointer g_iface, gpointer iface_data);

max_stars_repo_name: rudi-cilibrasi/classic-complearn

max_issues_repo_name: rudi-cilibrasi/classic-complearn

max_forks_repo_name: rudi-cilibrasi/classic-complearn


text:
/*

Excited States software: KGS
Contributors: See CONTRIBUTORS.txt
Contact: kgs-contact@simtk.org

Copyright (C) 2009-2017 Stanford University

Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies
of the Software, and to permit persons to whom the Software is furnished to do
so, subject to the following conditions:

This entire text, including the above copyright notice and this permission notice
shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
IN THE SOFTWARE.

*/

#ifndef KGS_LSNRELATIVEDIRECTION_H
#define KGS_LSNRELATIVEDIRECTION_H

#include <gsl/gsl_vector.h>

#include "Selection.h"
#include "Direction.h"
#include "core/Configuration.h"

class LSNrelativeDirection: public Direction  {
 public:
  LSNrelativeDirection(Selection& atomsMoving,
                       const std::vector< std::tuple<Atom*, Atom*, double> >& goal_distances);

  const std::vector< std::tuple<Atom*, Atom*, double> > getGoalDistances(){return m_goalDistances;}
  void setGoalDistances(std::vector< std::tuple<Atom*, Atom*, double> > goal_distances){m_goalDistances=goal_distances;}

 protected:
  void computeGradient(Configuration* conf, Configuration* conf2, gsl_vector* ret) override;
  void computeGradientExternalN(gsl_matrix* fullN, Configuration* conf, Configuration* conf2, gsl_vector* ret);

 private:

  void fillmatrices(Configuration* current_q, gsl_matrix* targetJacobian, gsl_matrix* targetDirection);
  // void clashFreeGradient(gsl_vector* gradient, gsl_vector* admissible_gradient, Molecule* protein);
  gsl_matrix* determineBestMove(gsl_matrix* N, gsl_matrix* targetJacobian, gsl_matrix* TargetPosition);
  Selection& m_atomsMovingSelection;
  std::vector< std::tuple<Atom*, Atom*, double> > m_goalDistances;

  friend class LSNclashAvoidingMove; //requires the computeGradientExternalN function
};

#endif //KGS_LSNRELATIVEDIRECTION_H


text:
#include <math.h>
#include <matrix.h>
#include <mex.h>
#include <string.h>
#include <blas.h>

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double * A = mxGetPr(prhs[0]);
    double * B = mxGetPr(prhs[1]);
    double * rows = mxGetPr(prhs[2]);
    double * cols = mxGetPr(prhs[3]);

    mwSize n_inds = mxGetM(prhs[2]);
    mwSize m_len_A = mxGetM(prhs[0]);
    mwSize m_len_B = mxGetM(prhs[1]);

    plhs[0] = mxCreateDoubleMatrix(n_inds, 1, mxREAL);
    double *pr = mxGetPr(plhs[0]);

    mwSize cur_row;
    mwSize cur_col;

    mwSignedIndex one = 1;

    mwSignedIndex dot_length = (mwSignedIndex)m_len_B;

    /* For each element we wish to calculate */
    for (int i = 0; i < n_inds; i++)
    {
        cur_row = rows[i];
        cur_col = cols[i];

        pr[i] = ddot(&dot_length, &A[m_len_A * (cur_row-1)], &one, &B[m_len_A * (cur_col-1)], &one);
    }

}

max_stars_repo_name: sjtrny/kssc

max_issues_repo_name: sjtrny/kssc

max_forks_repo_name: sjtrny/kssc