vide_public/c_tools/hod/buh.c

2289 lines
52 KiB
C

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>
#ifdef PARALLEL
#include <mpi.h>
#endif
#include "header.h"
#define MAGSOL 5.33
/* Data structures to hold multiple correlation functions
* and their corresponding covariance matrices.
*/
struct WP_ARRAY {
double *x;
double *r;
double *e;
double *eigen;
double **covar;
double mag1;
double mag2;
double mean_mag;
double ngal;
int np;
double M_min;
double sigma_logM;
double M1;
double M_low;
double M_cut;
double mmin_hi;
} wp_array[15];
struct gao_parameters {
double f0;
double s0;
} GAO;
double **covar_array[15];
/* Local globals for qromo
*/
double *mass_g9,
*lum_g9,
*sig_g9,
*y_lum_g9,
*y_sig_g9,
curlum_g9,
OMEGA_TEMP;
int n_g9,
VARY_GAMMA=0,
VARY_CVIR=0,
VARY_ALPHA=0;
/* External function for density dependence.
*/
void dd_hod_functions(float *hmass, float *hdensity, int nhalo);
/* Local functions.
*/
void ml_input_data(int);
void switch_wp_data(int n);
double schechter(double mag, double mstar, double phi0,double alpha);
double func_integrate_schechter(double m);
double func_integrate_lum_schechter(double m);
double ml_initialize(int nwp, int nhod, double *a, double **cov1, double *avg1, double *start_dev);
double chi2_ml_function(int count);
void mass2light_functions(int n_wp, double *a);
double total_halo_light(double mass, int n_wp, double *a, double *ss);
double func_ml_lum(double m);
double func_ml_norm(double m);
double chi2_CNOC_data(int count);
double chi2_ml_wrapper(double *a);
void ml_function_test(int n_wp, double *a);
void ml_function_test_monte_carlo(int n_wp, double *a);
int poisson_deviate(double nave);
double poisson_prob(int n, double nave);
void dispersion_test(int n_wp, double *a, double mass);
double random_luminosity(double m1, double m2);
void calc_multiplicity_function(int n_wp, int n_start);
double chi2_multiplicity_function(int count);
void analyze_chain(void);
void output_wp_fits(int n_wp, double *a);
/* Local Functions for simulation population.
*/
double NFW_density(double r, double rs, double ps);
double NFW_velocity(double mass, double v[]);
void NFW_position(double mass, double x[]);
void ml_mock_catalog(int n_wp, double *a);
double random_2df_magnitude_error(void);
void ml_mock_catalog_gao(int n_wp, double *a);
double N_cen_gao(double m, int ii, double *mmin);
double N_cen_hiden(double m);
/* External functions.
*/
double chi2_wp_wrapper(double *a);
void ml_minimization()
{
int i,j,k,n_wp=9,n_hod=0,n1,n2,n3;
double stepfac=1;
double error=1,tolerance=0,**cov1,**tmp,*a,*a1,*avg1,chi2,chi2prev,*start_dev,
**evect,*eval,*aprev,*atemp,**tmp1,*opar,x1,fsat,chi2array[15],**chain,x,x2,x3;
int n,nrot,niter=0,count=0,imax_chain=30000,NSTEP = 50;
long IDUM=-555;
int *pcheck,pcnt,ptot=10,nancheck,astart[10];
fprintf(stderr,"\n\nCHI2 MINIMIZATION OF M/L DATA..........\n");
fprintf(stderr, "--------------------------------------------\n\n");
ml_input_data(n_wp);
Work.imodel=2;
Work.chi2=1;
HOD.pdfc = 9;
/* In Zehavi et al, this is 40 Mpc/h,
* but Norberg et al use 50.
*/
wp.pi_max=50.0;
srand48(32498793);
astart[0] = 1;
for(i=1;i<=5;++i)
astart[i] = astart[i-1] + 2;
for(i=6;i<=8;++i)
astart[i] = astart[i-1] + 3;
n = astart[8]+2;
if(HOD.free[3]) {
VARY_ALPHA=1;
HOD.free[3] = 0;
n++;
}
if(HOD.free[6]) {
VARY_CVIR=1;
HOD.free[6] = 0;
n++;
}
if(HOD.free[10]==2) {
VARY_GAMMA=1;
HOD.free[10]=0;
}
for(i=8;i<100;++i)
if(HOD.free[i]) { n++; MCMC++; }
OMEGA_TEMP = 1.651693e-01;
printf("mcmc_min> %d free parameters\n",n);
pcheck=calloc(ptot,sizeof(int));
a=dvector(1,n);
start_dev=dvector(1,n);
a1=dvector(1,n_hod);
aprev=dvector(1,n);
atemp=dvector(1,n);
cov1=dmatrix(1,n,1,n);
avg1=dvector(1,n);
tmp=dmatrix(1,n,1,n);
tmp1=dmatrix(1,n,1,1);
evect=dmatrix(1,n,1,n);
eval=dvector(1,n);
chain=dmatrix(1,imax_chain,1,n);
chi2prev=ml_initialize(n_wp,n_hod,a,cov1,avg1,start_dev);
if(IDUM_MCMC==-1)
mass2light_functions(n_wp,a);
if(IDUM_MCMC==-2)
ml_mock_catalog(n_wp,a);
if(IDUM_MCMC==-3)
{
calc_multiplicity_function(n_wp,2);
calc_multiplicity_function(n_wp,5);
exit(0);
}
if(IDUM_MCMC==-4)
{
output_wp_fits(n_wp,a);
exit(0);
}
if(IDUM_MCMC==-5)
ml_mock_catalog_gao(n_wp,a);
if(IDUM_MCMC==-6)
ml_function_test(n_wp,a);
if(IDUM_MCMC==-7)
ml_function_test_monte_carlo(n_wp,a);
if(IDUM_MCMC==-8)
analyze_chain();
niter++;
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
IDUM=IDUM_MCMC;
pcnt=0;
pcheck[pcnt]=1;
stepfac=1;
while(niter<NSTEP)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
stepfac = stepfac*pow(0.9,5-j);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
if(VARY_GAMMA)
stepfac = 0.7;
for(i=1;i<=n;++i)
a[i] = (1+gasdev(&IDUM)*start_dev[i]*stepfac)*aprev[i];
if(MCMC>1)
{
RESET_COSMOLOGY++;
i = 7;
j = astart[8]+2+VARY_ALPHA+VARY_CVIR;
if(HOD.free[++i])OMEGA_TEMP = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])SIGV = a[++j];
if(HOD.free[++i])SPECTRAL_INDX = a[++j];
if(VARY_GAMMA)
GAMMA = gasdev(&IDUM)*0.02 + 0.15;
printf("COSMO %d %f %f %f %d\n",count+1,OMEGA_TEMP,SIGMA_8,GAMMA,j);
}
chi2=0;
/* Restrict CVIR_FAC
*/
i = astart[8]+3+VARY_ALPHA;
if(VARY_CVIR)
if(a[i]<0.3 || a[i]>1.2)continue;
if(HOD.free[10])
if(GAMMA<0.11 && GAMMA>0.19)continue;
++count;
for(i=n_wp-1;i>=0;--i)
{
n_hod = 3;
HOD.free[5] = 1;
if(i<6){ n_hod = 2; HOD.free[5] = 0; HOD.sigma_logM = 0.2; }
for(k=0,j=astart[i];j<=astart[i]+n_hod-1;++j)
{
a1[++k] = a[j];
}
if(VARY_ALPHA)
HOD.alpha = a[astart[8]+3];
if(VARY_CVIR)
CVIR_FAC = a[astart[8]+4];
wp.ncf = n_hod;
switch_wp_data(i);
chi2+=chi2array[i]=chi2_wp_wrapper(a1);
wp_array[i].M_min = HOD.M_min;
wp_array[i].sigma_logM = HOD.sigma_logM;
wp_array[i].M1 = HOD.M1;
wp_array[i].M_low = HOD.M_low;
wp_array[i].M_cut = HOD.M_cut;
if(!ThisTask)
{
printf("TRY%d %d ",i,count);
for(j=1;j<=n_hod;++j)
printf("%.4e ",a1[j]);
printf("%e\n",chi2array[i]);
fflush(stdout);
}
}
chi2+=chi2_ml_function(count);
chi2+=chi2_multiplicity_function(count);
if(!ThisTask)
{
printf("TRY_ALL %d %e %e\n",count,chi2,chi2prev);
fflush(stdout);
}
pcheck[pcnt]=0;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
continue;
pcheck[pcnt]=1;
niter++;
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask){
printf("ACC_ALL %d %d %e\n",niter,count,chi2);
for(j=0;j<n_wp;++j)
{
printf("ACC%d %d %d ",j,niter,count);
n_hod = 3;
if(j<6){ n_hod = 2; }
for(k=0,i=astart[j];i<=astart[j]+n_hod-1;++i)
printf("%e ",a[i]);
for(i=astart[8]+3;i<=n;++i)
printf("%e ",a[i]);
printf("%e %e\n",chi2array[j],chi2);fflush(stdout);
}
}
}
stepfac=1.0;
pcnt=-1;
while(niter<imax_chain)
{
pcnt++;
if(pcnt==ptot)
{
for(j=i=0;i<ptot;++i)j+=pcheck[i];
stepfac = stepfac*pow(0.9,5-j);
if(!ThisTask)printf("STEPFAC %f %d %d\n",stepfac,j,count);
pcnt=0;
}
if(VARY_GAMMA)
stepfac = 0.7;
for(j=1;j<=n;++j)
{
avg1[j]=0;
for(k=1;k<=n;++k)
cov1[j][k]=0;
}
for(i=niter-NSTEP+1;i<=niter;++i)
{
for(j=1;j<=n;++j)
{
avg1[j]+=chain[i][j];
for(k=1;k<=n;++k)
cov1[j][k]+=chain[i][j]*chain[i][k];
}
}
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
tmp[i][j] = cov1[i][j]/NSTEP - avg1[i]*avg1[j]/NSTEP/NSTEP;
jacobi(tmp,n,eval,evect,&nrot);
gaussj(evect,n,tmp1,1);
/* Check to make sure the eigenvalues are all there.
*/
nancheck=0;
for(i=1;i<=n;++i)
if(isnan(eval[i])|| eval[i]<=0)nancheck=1;
for(i=1;i<=n;++i)
printf("%d %e\n",i,eval[i]);
if(nancheck)
{
printf("NANCHECK %d\n",count);
for(i=1;i<=n;++i)
a[i] = gasdev(&IDUM)*sqrt(tmp[i][i])*stepfac;
}
else
{
for(i=1;i<=n;++i)
atemp[i] = gasdev(&IDUM)*sqrt(eval[i])*stepfac;
for(i=1;i<=n;++i)
for(a[i]=0,j=1;j<=n;++j)
a[i] += atemp[j]*evect[j][i];
}
for(i=1;i<=n;++i)
a[i] += aprev[i];
if(MCMC>1)
{
RESET_COSMOLOGY++;
i=7;
j = astart[8]+2+VARY_ALPHA+VARY_CVIR;
if(HOD.free[++i])OMEGA_TEMP = a[++j];
if(HOD.free[++i])SIGMA_8 = a[++j];
if(HOD.free[++i])GAMMA = a[++j];
if(HOD.free[++i])VBIAS = a[++j];
if(HOD.free[++i])VBIAS_C = a[++j];
if(HOD.free[++i])SIGV = a[++j];
if(VARY_GAMMA)
GAMMA = gasdev(&IDUM)*0.02 + 0.15;
printf("COSMO %f %f %f\n",OMEGA_TEMP,SIGMA_8,GAMMA);
}
chi2=0;
/* Restrict CVIR_FAC and GAMMA
*/
i = astart[8]+3+VARY_ALPHA;
if(VARY_CVIR)
if(a[i]<0.3 || a[i]>1.2)continue;
if(HOD.free[10])
if(GAMMA<0.11 && GAMMA>0.19)continue;
++count;
for(i=n_wp-1;i>=0;--i)
{
n_hod = 3;
HOD.free[5] = 1;
if(i<6){ n_hod = 2; HOD.free[5] = 0; HOD.sigma_logM = 0.2; }
for(k=0,j=astart[i];j<=astart[i]+n_hod-1;++j)
{
a1[++k] = a[j];
}
if(VARY_ALPHA)
HOD.alpha = a[astart[8]+3];
if(VARY_CVIR)
CVIR_FAC = a[astart[8]+4];
wp.ncf = n_hod;
switch_wp_data(i);
chi2+=chi2array[i]=chi2_wp_wrapper(a1);
wp_array[i].M_min = HOD.M_min;
wp_array[i].sigma_logM = HOD.sigma_logM;
wp_array[i].M1 = HOD.M1;
wp_array[i].M_low = HOD.M_low;
wp_array[i].M_cut = HOD.M_cut;
if(!ThisTask)
{
printf("TRY%d %d ",i,count);
for(j=1;j<=n_hod;++j)
printf("%.4e ",a1[j]);
printf("%e\n",chi2array[i]);
fflush(stdout);
}
}
chi2+=chi2_ml_function(count);
chi2+=chi2_multiplicity_function(count);
if(!ThisTask)
{
printf("TRY_ALL %d %e %e\n",count,chi2,chi2prev);
fflush(stdout);
}
pcheck[pcnt]=0;
if(!(chi2<chi2prev || drand48() <= exp(-(chi2-chi2prev)/2)))
continue;
pcheck[pcnt]=1;
niter++;
for(i=1;i<=n;++i)
chain[niter][i]=a[i];
for(i=1;i<=n;++i)
avg1[i] += a[i];
for(i=1;i<=n;++i)
aprev[i] = a[i];
for(i=1;i<=n;++i)
for(j=1;j<=n;++j)
cov1[i][j] += a[i]*a[j];
chi2prev=chi2;
if(!ThisTask){
printf("ACC_ALL %d %d %e\n",niter,count,chi2);
for(j=0;j<n_wp;++j)
{
printf("ACC%d %d %d ",j,niter,count);
n_hod = 3;
if(j<6){ n_hod = 2; }
for(k=0,i=astart[j];i<=astart[j]+n_hod-1;++i)
printf("%e ",a[i]);
for(i=astart[8]+3;i<=n;++i)
printf("%e ",a[i]);
printf("%e %e\n",chi2array[j],chi2);fflush(stdout);
}
}
}
exit(0);
}
/* This routine loops over the standard wp_input function that inputs wp data,
* and inserts each set of data into the wp_array structure.
* Currently I have hard-wired the magnitude limits and the filename.
*/
void ml_input_data(int n_wp)
{
int i,j,k;
FILE *fp;
double mstar=-19.66,phi0=0.0161,alpha=-1.21,low_magnitude=-17.0;
/*fp=openfile("galaxy_density_array.dat");*/
for(i=0;i<n_wp;++i)
{
wp_array[i].mag1 = low_magnitude-0.5*i;
wp_array[i].mag2 = low_magnitude-0.5*(i+1);
wp_array[i].mean_mag = low_magnitude-0.5*(i+.5); /* TEMP */
/* Old data -> no PCA, plus with innermost data point
*
sprintf(wp.fname_wp,"/gal/sr2/tinker/ml_ratio/data/xi_june04.%.02f.%.02f.keML.txt",
wp_array[i].mag2,wp_array[i].mag1);
*/
sprintf(wp.fname_wp,"/gal/sr2/tinker/voids/data/half_mag_samples/cov.10001.%.02f.%.02f.keML.com.xi2.2.log2",
wp_array[i].mag2,wp_array[i].mag1);
if(i==n_wp-1)
sprintf(wp.fname_wp,"/gal/sr2/tinker/voids/data/threshold_samples/cov.10001.-22.00.-21.00.keML.com.xi2.2.log2");
sprintf(wp.fname_covar,"xi_covar.%d",i);
wp_input();
wp_array[i].x=dvector(1,wp.np);
wp_array[i].r=dvector(1,wp.np);
wp_array[i].e=dvector(1,wp.np);
wp_array[i].eigen=dvector(1,wp.np);
wp_array[i].covar=dmatrix(1,wp.np,1,wp.np);
for(j=1;j<=wp.np;++j)
{
wp_array[i].x[j]=wp.x[j];
wp_array[i].r[j]=wp.r[j];
wp_array[i].e[j]=wp.e[j];
wp_array[i].eigen[j]=wp.eigen[j];
for(k=1;k<=wp.np;++k)
wp_array[i].covar[j][k] = wp.covar[j][k];
}
wp_array[i].np=wp.np;
free_dvector(wp.x,1,wp.np);
free_dvector(wp.r,1,wp.np);
free_dvector(wp.e,1,wp.np);
/* This is TEMP as well.
*/
if(i==8)wp_array[i].mag2 = -23;
wp_array[i].ngal = qromo(func_integrate_schechter,
wp_array[i].mag2,wp_array[i].mag1,midpnt);
if(i==8)wp_array[i].ngal = 1.25e-4;
fprintf(stderr,"%f %e\n",wp_array[i].mean_mag,wp_array[i].ngal);
/*fscanf(fp,"%lf",&wp_array[i].ngal);*/
}
/*fclose(fp);*/
wp.x=dvector(1,wp.np);
wp.r=dvector(1,wp.np);
wp.e=dvector(1,wp.np);
}
/* This is just the value of the Schechter Function fit to
* 2F at magnitude m--> for getting total number density of galaxies.
*/
double func_integrate_schechter(double m)
{
double mstar=-19.66,phi0=0.0161,alpha=-1.21;
return(schechter(m,mstar,phi0,alpha));
}
/* This is the luminosity-weighted Schechter function value
* at magnitdue m--> to get the total amount of light.
*/
double func_integrate_lum_schechter(double m)
{
double mstar=-19.66,phi0=0.0161,alpha=-1.21,lum;
lum = pow(10.0,-0.4*(m-MAGSOL));
return(schechter(m,mstar,phi0,alpha)*lum);
}
/* Function to return the value of the Schechter function
* given the parameters.
*/
double schechter(double mag, double mstar, double phi0,double alpha)
{
return(0.4*log(10.0)*phi0*pow(10.0,-0.4*(mag-mstar)*(alpha+1))*
exp(-pow(10.0,-0.4*(mag-mstar))));
}
/* Routine to switch the global wp data variables with
* the current dataset in the chi^2 loop.
*/
void switch_wp_data(int n)
{
int i,j;
i=n;
for(j=1;j<=wp_array[i].np;++j)
{
wp.x[j]=wp_array[i].x[j];
wp.r[j]=wp_array[i].r[j];
wp.e[j]=wp_array[i].e[j];
}
wp.np=wp_array[i].np;
GALAXY_DENSITY=wp_array[i].ngal;
HOD.i_wp = n;
}
double ml_initialize(int nwp, int nhod, double *a, double **cov1, double *avg1, double *start_dev)
{
int i,j=0,k,astart[10],np;
double x1,x2,m1[50],alpha[50],mcut[50],sig[50],*a1,chi2,chi2array[15];
long IDUM = -556;
FILE *fp;
astart[0] = 1;
for(i=1;i<=5;++i)
astart[i] = astart[i-1] + 2;
for(i=6;i<=8;++i)
astart[i] = astart[i-1] + 3;
a1=dvector(1,nhod);
fp=openfile("initial_values.dat");
for(i=1;i<=nwp;++i)
fscanf(fp,"%lf %lf %lf %lf",&m1[i],&alpha[i],&mcut[i],&sig[i]);
fclose(fp);
j=0;
for(k=1;k<=nwp;++k)
{
i=0;
if(HOD.free[++i])
{
/*
* This one will eventually be set set to one.
*
*/
}
if(HOD.free[++i])
{
a[++j]=log10(m1[k]);
start_dev[j]=0.001;
}
/*
* Reserve alpha to get the same for all wp
*/
if(HOD.free[++i])
{
/*
a[++j]=alpha[k];
start_dev[j]=0.01;
*/
}
if(HOD.free[++i])
{
a[++j]=log10(mcut[k]);
start_dev[j]=0.001;
}
/* Only leave sigma_logM free for i_wp>=6
*/
if(HOD.free[++i] && k>5)
{
a[++j]=log10(sig[k]);
start_dev[j]=0.01;
}
}
if(VARY_ALPHA)
{
a[++j]=HOD.alpha;
start_dev[j]=0.01;
}
if(VARY_CVIR)
{
a[++j]=CVIR_FAC;
start_dev[j]=0.01;
}
np = j;
/* If using Powell's method or amoeba, then stop here
*/
for(i=0;i<nwp;++i)
{
printf("IV %d ",i);
nhod = 3;
if(i<6){ nhod = 2; }
for(j=astart[i];j<=astart[i]+nhod-1;++j)
printf("%f ",a[j]);
printf("%f ",HOD.alpha);
printf("\n");
fflush(stdout);
}
if(!MCMC)return 0;
if(MCMC>1)
{
i=7;
j=np;
if(HOD.free[++i]){ a[++j]=OMEGA_TEMP; start_dev[j]=0.005; np++; }
if(HOD.free[++i]){ a[++j]=SIGMA_8; start_dev[j]=0.01; np++; }
if(HOD.free[++i]){ a[++j]=GAMMA; start_dev[j]=0.01; np++; }
if(HOD.free[++i]){ a[++j]=VBIAS; start_dev[j]=0.01; np++; }
if(HOD.free[++i]){ a[++j]=VBIAS_C; start_dev[j]=0.01; np++; }
if(HOD.free[++i]){ a[++j]=SIGV; start_dev[j]=0.01; np++; }
}
for(i=1;i<=np;++i)
{
avg1[i]=a[i];
printf("BEGIN %d %f\n",i,a[i]);
for(j=1;j<=nhod*nwp;++j)
cov1[i][j]=a[i]*a[j];
}
chi2=0;
for(i=nwp-1;i>=0;--i)
{
nhod = 3;
HOD.free[5] = 1;
if(i<6){ nhod = 2; HOD.free[5] = 0; HOD.sigma_logM = 0.2; }
wp.ncf=nhod;
for(k=0,j=astart[i];j<=astart[i]+nhod-1;++j)
{
a1[++k] = a[j];
}
if(VARY_ALPHA)
HOD.alpha = a[astart[8]+3];
if(VARY_CVIR)
CVIR_FAC = a[astart[8]+4];
printf("ALPHA = %f CVIR = %f\n",HOD.alpha,CVIR_FAC);
switch_wp_data(i);
chi2+=chi2array[i]=chi2_wp_wrapper(a1);
printf("CHI %d %e\n",i,chi2array[i]);
wp_array[i].M_min = HOD.M_min;
wp_array[i].sigma_logM = HOD.sigma_logM;
wp_array[i].M1 = HOD.M1;
wp_array[i].M_low = HOD.M_low;
wp_array[i].M_cut = HOD.M_cut;
}
chi2+=chi2_ml_function(0) + chi2_multiplicity_function(0);
printf("INITIAL CHI2 %e\n",chi2);
return(chi2);
}
/* This function takes the 2PIGG multiplicity function and returns a chi^2 value
* based on the multplicity function implies by the current point in the chain.
*/
double chi2_multiplicity_function(int count)
{
int i,j,k,nmass=200,nmulti=145,n_start,n_wp=9;
double dlogm,mass,mlo,mhi,dndM,nsat,ncen,chi2a=0,chi2b=0,x,nlo,nhi;
static double *multi,*dy,*ngal;
FILE *fp;
static int flag=0,NP;
static float *N_2PIGG,*M1_2PIGG,*E1_2PIGG,*M2_2PIGG,*E2_2PIGG;
if(!flag)
{
flag = 1;
fp = openfile("/gal/sr2/tinker/ml_ratio/2PIGG_MULTI.dat");
NP = filesize(fp);
N_2PIGG = vector(1,NP);
E1_2PIGG = vector(1,NP);
M1_2PIGG = vector(1,NP);
E2_2PIGG = vector(1,NP);
M2_2PIGG = vector(1,NP);
for(i=1;i<=NP;++i)
{
fscanf(fp,"%f %f %f %f %f",&N_2PIGG[i],&M2_2PIGG[i],&E2_2PIGG[i],&M1_2PIGG[i],&E1_2PIGG[i]);
if(M1_2PIGG[i]>-99) {
M1_2PIGG[i] = pow(10.0,M1_2PIGG[i]);
E1_2PIGG[i] = log(10.0)*M1_2PIGG[i]*E1_2PIGG[i];
}
if(M2_2PIGG[i]>-99) {
M2_2PIGG[i] = pow(10.0,M2_2PIGG[i]);
E2_2PIGG[i] = log(10.0)*M2_2PIGG[i]*E2_2PIGG[i];
}
}
fclose(fp);
multi = dvector(1,nmulti);
dy = dvector(1,nmulti);
ngal = dvector(1,nmulti);
for(i=1;i<=nmulti;++i)
ngal[i] = i;
}
/* First do the multiplicity function for galaxies M_b<-18
*/
n_start = 2;
for(i=1;i<=nmulti;++i)
multi[i] = 0;
mlo = log(wp_array[n_start].M_low);
mhi = log(HOD.M_max);
dlogm = (mhi - mlo)/(nmass-1);
for(j=1;j<=nmass;++j)
{
mass = exp((j-1)*dlogm + mlo);
dndM = dndM_interp(mass);
nsat = ncen = 0;
for(i=n_start;i<n_wp;++i)
{
HOD.i_wp = i;
HOD.M1 = wp_array[i].M1;
HOD.M_min = wp_array[i].M_min;
HOD.M_low = wp_array[i].M_low;
HOD.M_cut = wp_array[i].M_cut;
HOD.sigma_logM = wp_array[i].sigma_logM;
nsat += N_sat(mass);
ncen += N_cen(mass);
}
for(k=1;k<=nmulti;++k)
multi[k] += poisson_prob(k-1,nsat)*dndM*dlogm*mass*ncen;
}
for(i=2;i<=NP;++i)
{
if(M1_2PIGG[i]<0)continue;
nlo = pow(10.0,N_2PIGG[i]-0.14137/2.0);
nhi = pow(10.0,N_2PIGG[i]+0.14137/2.0);
x=0;
for(j=(int)nlo-1;j<=(int)nhi+1;++j)
if(j>0 && j>=nlo && j<nhi)
x+=multi[j];
x/=0.14137;
chi2a += (x-M1_2PIGG[i])*(x-M1_2PIGG[i])/E1_2PIGG[i]/E1_2PIGG[i];
printf("MULTI%d %d %e %e %e %e\n",n_start,count,N_2PIGG[i],x,M1_2PIGG[i],E1_2PIGG[i]);
}
/* First do the multiplicity function for galaxies M_b<-19.5
*/
n_start = 5;
for(i=1;i<=nmulti;++i)
multi[i] = 0;
mlo = log(wp_array[n_start].M_low);
mhi = log(HOD.M_max);
dlogm = (mhi - mlo)/(nmass-1);
for(j=1;j<=nmass;++j)
{
mass = exp((j-1)*dlogm + mlo);
dndM = dndM_interp(mass);
nsat = ncen = 0;
for(i=n_start;i<n_wp;++i)
{
HOD.i_wp = i;
HOD.M1 = wp_array[i].M1;
HOD.M_min = wp_array[i].M_min;
HOD.M_low = wp_array[i].M_low;
HOD.M_cut = wp_array[i].M_cut;
HOD.sigma_logM = wp_array[i].sigma_logM;
nsat += N_sat(mass);
ncen += N_cen(mass);
}
for(k=1;k<=nmulti;++k)
multi[k] += poisson_prob(k-1,nsat)*dndM*dlogm*mass*ncen;
}
for(i=2;i<=NP;++i)
{
if(M2_2PIGG[i]<0)continue;
nlo = pow(10.0,N_2PIGG[i]-0.14137/2.0);
nhi = pow(10.0,N_2PIGG[i]+0.14137/2.0);
x=0;
for(j=(int)nlo-1;j<=(int)nhi+1;++j)
if(j>0 && j>=nlo && j<nhi)
x+=multi[j];
x/=0.14137;
chi2b += (x-M2_2PIGG[i])*(x-M2_2PIGG[i])/E2_2PIGG[i]/E2_2PIGG[i];
printf("MULTI%d %d %e %e %e %e\n",n_start,count,N_2PIGG[i],x,M2_2PIGG[i],E2_2PIGG[i]);
}
/* Sum up the chi2 values and return.
*/
printf("CHI2_MULTI %d %e %e %e\n",count,chi2a+chi2b,chi2a,chi2b);
return(chi2a+chi2b);
}
/* This function takes in the 2PIGG (corrected) data and returns a chi^2 value
* of the M/L-L curve
*/
double chi2_ml_function(int count)
{
FILE *fp;
int i,j,k,n=100,n_wp=9;
double *lum,*sig,mlo,mhi,dlogm,*mass,*ml_lum,*y_lum, *y_sig, *a,chi2,ml;
double t0,t1;
static int flag=0,N_2PIGG,N_LENSING;
static float *ML_2PIGG, *L_2PIGG,*E1_2PIGG,*E2_2PIGG,*ML_LENSING,*L_LENSING,*E_LENSING;
/* return(0); */
if(!flag)
{
flag = 1;
fp = openfile("/gal/sr2/tinker/ml_ratio/2PIGG_ML_linear.dat");
N_2PIGG = filesize(fp);
muh(N_2PIGG);
L_2PIGG = vector(1,N_2PIGG);
E1_2PIGG = vector(1,N_2PIGG);
E2_2PIGG = vector(1,N_2PIGG);
ML_2PIGG = vector(1,N_2PIGG);
for(i=1;i<=N_2PIGG;++i)
fscanf(fp,"%f %f %f %f",&L_2PIGG[i],&ML_2PIGG[i],&E1_2PIGG[i],&E2_2PIGG[i]);
fclose(fp);
fp = openfile("/gal/sr2/tinker/ml_ratio/LENSING_ML.dat");
N_LENSING = filesize(fp);
muh(N_LENSING);
L_LENSING = vector(1,N_LENSING);
E_LENSING = vector(1,N_LENSING);
ML_LENSING = vector(1,N_LENSING);
for(i=1;i<=N_LENSING;++i)
fscanf(fp,"%f %f %f",&L_LENSING[i],&ML_LENSING[i],&E_LENSING[i]);
fclose(fp);
}
n_g9=n;
mlo=log(1.0e11);
mhi=log(5.0e15);
dlogm=(mhi-mlo)/(n-1);
mass_g9=dvector(1,n);
lum_g9=dvector(1,n);
sig_g9=dvector(1,n);
ml_lum=dvector(1,n);
y_lum_g9=dvector(1,n);
y_sig_g9=dvector(1,n);
/* Tabulate <M/L>_M and sig_L(M) curves
*/
for(i=1;i<=n;++i)
{
curlum_g9=1;
if(i==1)curlum_g9=-1;
mass_g9[i] = exp((i-1)*dlogm+mlo);
lum_g9[i] = total_halo_light(mass_g9[i],n_wp,a,&sig_g9[i]);
if(mass_g9[i]<wp_array[0].M1/2)sig_g9[i]/=2;
}
spline(mass_g9,lum_g9,n,1.0E+30,1.0E+30,y_lum_g9);
spline(mass_g9,sig_g9,n,1.0E+30,1.0E+30,y_sig_g9);
/* Calculate the (M/L)_L curve at the 2PIGG locations.
* Here is where we use OMEGA_TEMP: scale the M/L values by
* the new value of OMEGA.
*/
chi2=0;
for(i=1;i<=N_2PIGG;++i)
{
curlum_g9 = L_2PIGG[i];
ml = qromo(func_ml_lum,mlo,mhi,midpnt)/qromo(func_ml_norm,mlo,mhi,midpnt);
ml = ml*OMEGA_TEMP/OMEGA_M;
if(ml>ML_2PIGG[i])
chi2 += (ml-ML_2PIGG[i])*(ml-ML_2PIGG[i])/E1_2PIGG[i]/E1_2PIGG[i];
else
chi2 += (ml-ML_2PIGG[i])*(ml-ML_2PIGG[i])/E2_2PIGG[i]/E2_2PIGG[i];
printf("MLC %d %e %e %e\n",count,curlum_g9,ML_2PIGG[i],ml);
}
for(i=1;i<=N_LENSING;++i)
{
curlum_g9 = L_LENSING[i];
ml = qromo(func_ml_lum,mlo,mhi,midpnt)/qromo(func_ml_norm,mlo,mhi,midpnt);
ml = ml*OMEGA_TEMP/OMEGA_M;
chi2 += (ml-ML_LENSING[i])*(ml-ML_LENSING[i])/E_LENSING[i]/E_LENSING[i];
printf("MLC %d %e %e %e\n",count,curlum_g9,ML_LENSING[i],ml);
}
printf("ML_CHI2 %d %e\n",count,chi2);
chi2+=chi2_CNOC_data(count);
return(chi2);
}
double chi2_CNOC_data(int count)
{
static int flag=0,n;
static float *mass,x1,x2,lum_avg=0;
double chi2,dx1,*ax1,lum,ml_ave=0;
FILE *fp;
int i;
if(!flag)
{
flag = 1;
fp = openfile("/gal/sr2/tinker/ml_ratio/CNOC_ML.dat");
n = filesize(fp);
mass = vector(1,n);
for(i=1;i<=n;++i)
{
fscanf(fp,"%f %f %f",&mass[i],&x1,&x2);
lum_avg += mass[i]/x1;
mass[i] /= (OMEGA_TEMP/OMEGA_M);
}
lum_avg=lum_avg/n/1.07*0.887;
fclose(fp);
}
for(i=1;i<=n;++i)
{
lum = total_halo_light(mass[i],9,ax1,&dx1);
ml_ave += mass[i]/lum/(OMEGA_M/OMEGA_TEMP);
}
ml_ave/=n;
chi2 = (ml_ave - 382)*(ml_ave - 382)/34./34.;
printf("MLC %d %e %e %e %f\n",count,lum_avg,ml_ave,382.0,OMEGA_TEMP/OMEGA_M);
return(chi2);
}
/* This function converts the (M/L)_M curve to the (M/L)_L curve.
* It assumes that the dispersion about <M/L> at fixed M is a Gaussian
* with variance as the sum of the individual variances for each mag bin.
*
* NB- Each sat function is Poisson, and each cen function is Bernoulli,
* but with nine of each we're reaching the central limit theorem. Works great
* for M>10^14 M_sol, but more of an approximation for 10^13 M_sol.
*/
void mass2light_functions(int n_wp, double *a)
{
int i,j,k,n=100;
double *lum,*sig,mlo,mhi,dlogm,*mass,*ml_lum,*y_lum, *y_sig;
double t0,t1;
n_g9=n;
mlo=log(1.0e11);
mhi=log(5.0e15);
dlogm=(mhi-mlo)/(n-1);
mass_g9=dvector(1,n);
lum_g9=dvector(1,n);
sig_g9=dvector(1,n);
ml_lum=dvector(1,n);
y_lum_g9=dvector(1,n);
y_sig_g9=dvector(1,n);
t0 = clock();
/* Tabulate <M/L>_M and sig_L(M) curves
*/
for(i=1;i<=n;++i)
{
mass_g9[i] = exp((i-1)*dlogm+mlo);
lum_g9[i] = total_halo_light(mass_g9[i],n_wp,a,&sig_g9[i]);
if(mass_g9[i]<wp_array[0].M1/2)sig_g9[i]/=2;
/* printf("ML %e %e %e\n",log10(mass_g9[i]),log10(lum_g9[i]),mass_g9[i]/lum_g9[i]); */
}
/* exit(0); */
fprintf(stderr,"TOTAL LUM DENSITY: %e\n",qromo(func_integrate_lum_schechter,-24.0,-13.0,midpnt));
spline(mass_g9,lum_g9,n,1.0E+30,1.0E+30,y_lum_g9);
spline(mass_g9,sig_g9,n,1.0E+30,1.0E+30,y_sig_g9);
/* Calculate the (M/L)_L curve
*/
for(i=1;i<=n;++i)
{
curlum_g9 = lum_g9[i];
ml_lum[i] = qromo(func_ml_lum,mlo,mhi,midpnt)/qromo(func_ml_norm,mlo,mhi,midpnt);
printf("ML %e %e %e %e %e\n",log10(mass_g9[i]),log10(lum_g9[i]),
mass_g9[i]/lum_g9[i],ml_lum[i],log10(sig_g9[i]));
}
t1 = clock();
fprintf(stderr,"%.2f\n", difftime(t1,t0)/CLOCKS_PER_SEC);
exit(0);
}
/* Integrand of (M/L)_L function.
*/
double func_ml_lum(double m)
{
double s,l,x;
m=exp(m);
splint(mass_g9,lum_g9,y_lum_g9,n_g9,m,&l);
splint(mass_g9,sig_g9,y_sig_g9,n_g9,m,&s);
x=(curlum_g9-l)/s;
x=exp(-x*x/2)/(RT2PI*s)*m/l*dndM_interp(m);
return(x*m);
}
/* Integrand of (M/L)_L normalization.
*/
double func_ml_norm(double m)
{
double s,l,x;
m=exp(m);
splint(mass_g9,lum_g9,y_lum_g9,n_g9,m,&l);
splint(mass_g9,sig_g9,y_sig_g9,n_g9,m,&s);
x=(curlum_g9-l)/s;
x=exp(-x*x/2)/(RT2PI*s)*dndM_interp(m);
return(x*m);
}
/* Sums up the total amount of light in each halo, also tabulates dispersion.
*/
double total_halo_light(double mass, int n_wp, double *a, double *ss)
{
int i,j,k;
double ncen,nsat,mean=0,err=0,norm,corr_mean,p=1,mean_cen=0;
static double mmin[15];
if(OUTPUT==2)
printf("OCC %e ",log10(mass));
mean=err=0;
for(j=n_wp-1;j>=0;--j)
{
HOD.i_wp=j;
HOD.M1 = wp_array[j].M1;
HOD.M_cut = wp_array[j].M_cut;
HOD.M_min = wp_array[j].M_min;
HOD.sigma_logM = wp_array[j].sigma_logM;
HOD.M_low = wp_array[j].M_low;
ncen = N_cen(mass);
nsat = N_sat(mass);
/* Sum up probability of no satellite galaxies in this halos.
*/
p*=exp(-nsat);
mean_cen+=ncen*pow(10.0,-0.4*(wp_array[j].mean_mag-MAGSOL));
/* Sum the squares of the errors
*/
err+=(nsat+ncen*(1-ncen))*pow(10.0,-0.8*(wp_array[j].mean_mag-MAGSOL));
mean+=(ncen+nsat)*pow(10.0,-0.4*(wp_array[j].mean_mag-MAGSOL));
if(OUTPUT==2)
printf("%e %e ",ncen,nsat);
}
norm = nsat/func_integrate_schechter(-17.25);
if(OUTPUT==2){
printf("\n");
fflush(stdout);
}
*ss = sqrt(err);
/* Use the value of N_sat at -17>M_b>-17.5 to normalize the schechter function
* for extrapolating to lower luminosities.
*/
norm = nsat/func_integrate_schechter(-17.25);
corr_mean = mean + qromo(func_integrate_lum_schechter,-17.0,-13.0,midpnt)*norm;
mean = corr_mean;
return(mean);
/* Now remove systems with no satellite galaxies from the means.
*/
if(p==1)p=0;
corr_mean = (mean - mean_cen*p)/(1-p);
mean = corr_mean;
fflush(stdout);
return(mean);
}
/********************************************************************************
* Code for getting the central occupation funciton based on the magnitude
* bin in question, and the next highest magnitude bin. (The lower limit of the
* brighter mag bin will define the upper limit of the fainter bin.)
*
* The general form of each will be pdfc==2.
*
* The routine assumes that the values in the global HOD structure
* represent the HOD parameters to be used for the ii node in the array
* (which is relevent since M_min must be set after the other params are chosen).
*
* If ii==8, then we are at the largest luminosity bin, which we'll treat
* currently as a threshold sample.
*/
double N_cen_i(double m, int ii)
{
int i;
double n1,n2,n,logm,n2max=0;
logm=log10(m);
n1=0.5*(1+erf((logm - log10(HOD.M_min))/HOD.sigma_logM));
if(ii==8)return(n1);
for(n2=0,i=ii+1;i<9;++i)
if(m>wp_array[i].M_low)
{
n2=0.5*(1+erf((logm - log10(wp_array[i].M_min))/wp_array[i].sigma_logM));
if(n2>n2max)n2max=n2;
}
n = n1-n2max;
/*
if(n2>1)n2=1;
if(n1+n2>1) n = 1-n2;
else n = n1;
if(ii==6 && (logm<13.01 && logm>12.99))
{
printf("NCEN %e %e %e %e %e %e\n",n1,n2,HOD.M_min,
HOD.sigma_logM,wp_array[ii+1].M_min,wp_array[ii+1].sigma_logM);
}
*/
if(n<0)return(0);
return(n);
}
/*********************************************************************************/
/*********************************************************************************
*
* The following functions are to calculate the multiplicity function
* for a given magnitude threshold.
*/
/*********************************************************************************/
/*********************************************************************************/
void calc_multiplicity_function(int n_wp, int n_start)
{
int i,j,k,nmass=2000,nmulti=200;
double dlogm,mass,*multi,mlo,mhi,dndM,nsat,*ncumu,ncen;
multi = dvector(1,nmulti);
ncumu = dvector(1,nmulti);
for(i=1;i<=nmulti;++i)
multi[i] = 0;
mlo = log(wp_array[n_start].M_low);
mhi = log(HOD.M_max);
dlogm = (mhi - mlo)/(nmass-1);
for(j=1;j<=nmass;++j)
{
mass = exp((j-1)*dlogm + mlo);
dndM = dndM_interp(mass);
nsat = ncen = 0;
for(i=n_start;i<n_wp;++i)
{
HOD.i_wp = i;
HOD.M1 = wp_array[i].M1;
HOD.M_min = wp_array[i].M_min;
HOD.M_low = wp_array[i].M_low;
HOD.M_cut = wp_array[i].M_cut;
HOD.sigma_logM = wp_array[i].sigma_logM;
nsat += N_sat(mass);
ncen += N_cen(mass);
}
for(k=1;k<=nmulti;++k)
{
multi[k] += poisson_prob(k-1,nsat)*dndM*dlogm*mass*ncen;
}
}
ncumu[nmulti]=multi[nmulti];
for(i=nmulti-1;i>0;--i)
ncumu[i] = ncumu[i+1]+multi[i];
for(i=1;i<=nmulti;++i)
printf("MULTI%d %d %e %e\n",n_start,i,multi[i],ncumu[i]);
}
/*********************************************************************************/
/*********************************************************************************
*
* The following functions are for using powell/amoeba to minimze the set of
* correlation functions.
*
*/
/*********************************************************************************/
/*********************************************************************************/
void ml_powell()
{
int n,niter,i,j,n_wp=9,n_hod,n_mag,k,astart[10];
double *a,**pp,*yy,FTOL=1.0E-3,chi2min,s1,dlogm,m,d[100],*avg1,**cov1,*start_dev,*a1;
FILE *fp;
char aa[1000];
fprintf(stderr,"\n\nCHI2 MINIMIZATION OF BINNED W_P(R_P) DATA..........\n");
fprintf(stderr, "---------------------------------------------------\n\n");
if(POWELL)
FTOL=1.0E-3;
else
FTOL=1.0E-4;
/* In Zehavi et al, this is 40 Mpc/h,
* but Norberg et al use 50.
*/
wp.pi_max=50.0;
ml_input_data(n_wp);
Work.imodel=2;
Work.chi2=1;
HOD.pdfc = 9;
MCMC=0;
srand48(32498793);
/* Find the number of free parameters in the minimization
* for the real-space correlation function.
*/
n_hod=0;
for(n=0,i=1;i<100;++i)
{
if(i<=7)
n+=HOD.free[i]*n_wp;
else
n+=HOD.free[i];
if(i<=7)
n_hod+=HOD.free[i];
if(OUTPUT)
printf("mcmc_min> free[%i] = %d\n",i,HOD.free[i]);
}
wp.ncf=n_hod;
wp.ncf_tot=n;
a=dvector(1,n);
start_dev=dvector(1,n);
ml_initialize(n_wp,n_hod,a,cov1,avg1,start_dev);
astart[0] = 1;
for(i=1;i<=5;++i)
astart[i] = astart[i-1] + 2;
for(i=6;i<=8;++i)
astart[i] = astart[i-1] + 3;
for(i=1;i<=n;++i)
a[i]=pow(10.0,a[i]);
for(i=1;i<=n;++i)
printf("a[%d] = %e\n",i,a[i]);
n_mag=n_wp-1;
/* Loop through all magnitudes.
*/
MAGNITUDE_LOOP:
n = n_hod = 3;
HOD.free[5] = 1;
if(n_mag<6){ n = n_hod = 2; HOD.free[5] = 0; HOD.sigma_logM = 0.2;}
n = n_hod;
a1=dvector(1,n);
if(POWELL)
pp=dmatrix(1,n,1,n);
else
pp=dmatrix(1,n+1,1,n);
yy=dvector(1,n+1);
for(k=0,j=astart[n_mag];j<=astart[n_mag]+n_hod-1;++j)
a1[++k] = a[j];
switch_wp_data(n_mag);
/* Make the starting stepsize 10% of the initial values.
*/
for(i=1;i<=n;++i)
d[i]=a1[i]*0.1;
if(POWELL)
{
for(i=1;i<=n;++i)
{
for(j=1;j<=n;++j)
{
pp[i][j]=0;
if(i==j)pp[i][j]+=d[j];
}
}
}
else
{
for(j=1;j<=n;++j)
pp[1][j]=a1[j];
yy[1]=chi2_ml_wrapper(a1);
for(i=1;i<=n;++i)
{
a1[i]+=d[i];
if(i>1)a1[i-1]-=d[i-1];
yy[i+1]=chi2_ml_wrapper(a1);
for(j=1;j<=n;++j)
pp[i+1][j]=a1[j];
}
a1[wp.ncf_tot]-=d[wp.ncf_tot];
}
if(POWELL)
{
if(OUTPUT)printf("wp_min> starting powell.\n");
powell(a1,pp,n,FTOL,&niter,&chi2min,chi2_wp);
chi2min = chi2_wp(a1);
}
else
{
if(OUTPUT)printf("wp_min> starting amoeba.\n");
amoeba(pp,yy,n,FTOL,chi2_wp,&niter);
for(i=1;i<=n;++i)a1[i]=pp[1][i];
chi2min = chi2_wp(a1);
}
printf("POWELL%d %e %e ",n_mag,chi2min,HOD.M_min);
for(i=1;i<=n;++i)
printf("%e ",a1[i]);
printf("%e\n",GALAXY_BIAS);
wp_array[n_mag].M_min = HOD.M_min;
wp_array[n_mag].M_low = HOD.M_low;
wp_array[n_mag].M_cut = HOD.M_cut;
wp_array[n_mag].M1 = HOD.M1;
wp_array[n_mag].sigma_logM = HOD.sigma_logM;
n_mag--;
if(n_mag<0)
exit(0);
free_dvector(a1,1,n);
free_dvector(yy,1,n+1);
if(POWELL)
free_dmatrix(pp,1,n,1,n);
else
free_dmatrix(pp,1,n+1,1,n);
goto MAGNITUDE_LOOP;
}
/* This routine loops through all 9 magnitude bins
* to get the total chi^2 for a set of parameters.
*/
double chi2_ml_wrapper(double *a)
{
static double *a1;
static int flag=1,niter=0;
double chi2,chi2i,t0,t1;
int i,j,nwp,nhod;
t0=clock();
nwp=9;
nhod=wp.ncf;
if(flag)
{
flag=0;
a1=dvector(1,nhod);
}
chi2=0;
++niter;
for(i=nwp-1;i>=0;--i)
{
for(j=1;j<=nhod;++j)
a1[j] = a[i+(j-1)*nwp+1];
switch_wp_data(i);
chi2+=chi2i=chi2_wp(a1);
printf("TRY%d %d ",i,niter);
for(j=1;j<=nhod;++j)
printf("%.4e ",a1[j]);
printf("%e\n",chi2i);
fflush(stdout);
wp_array[i].M_min = HOD.M_min;
wp_array[i].sigma_logM = HOD.sigma_logM;
}
t1=clock();
printf("ITERALL %d %e %.2f\n",niter,chi2,difftime(t1,t0)/CLOCKS_PER_SEC);
return(chi2);
}
/*********************************************************************************
*
* This is for Monte Carlo testing of some of the routines.
*
*/
/*********************************************************************************/
/*********************************************************************************/
/*********************************************************************************/
void ml_function_test(int n_wp, double *a)
{
FILE *fp;
int i,j,n,imass;
double mass;
char aa[1000];
/* fp=openfile("/gal/sr2/tinker/voids/LANL/halo.LANL.400"); */
fp=openfile("/gal/sr2/tinker/voids/LANL/halo.400_z.5.dat");
while(!feof(fp))
{
fscanf(fp,"%d %d",&i,&imass);
fgets(aa,1000,fp);
mass=imass*RHO_CRIT*OMEGA_M*pow(0.3125,3.0);
dispersion_test(n_wp,a,mass);
if(feof(fp))break;
}
exit(0);
}
/* does the same as the above function, but instead of using a mass
* function for a simulation, it uses a mass funcation created by
* monte carlo sampling of the analytic function.
*/
void ml_function_test_monte_carlo(int n_wp, double *a)
{
FILE *fp;
int i,j,n,imass;
double mass,pmax,p;
char aa[1000];
long IDUMx=-442123;
n=1E6;
i=0;
pmax = dndM_interp(1.0E11);
while(i<n)
{
mass = (1.0e16)*ran2(&IDUMx) + 1.0e11;
p = dndM_interp(mass)/pmax;
if(ran2(&IDUMx)>p)continue;
i++;
dispersion_test(n_wp,a,mass);
}
}
void ml_mock_catalog(int n_wp, double *a)
{
FILE *fp,*fpa[9],*fp2,*fpb[9];
int i,j,k,n,imass,n1,j_start=0,i1,galcnt[9][1000],halocnt[1000],imag;
double mass,xg[3],vg[3],nsat,nc[10],ncen,mlo,mag,err1,err2;
char aa[1000];
float x1,xh[3],vh[3];
long IDUM3 = -445;
fp=openfile("/gal/sr2/tinker/voids/LANL/halo.LANL.400");
for(i=j_start;i<n_wp;++i)
{
sprintf(aa,"%s.mock.%d",Task.root_filename,i);
fpa[i] = fopen(aa,"w");
sprintf(aa,"%s.mag.%d",Task.root_filename,i);
fpb[i] = fopen(aa,"w");
}
for(i=0;i<1000;++i)
halocnt[i]=0;
for(j=0;j<n_wp;++j)
for(i=0;i<1000;++i)
galcnt[j][i]=0;
j=j_start;
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M1 = pow(10.0,a[j+0*n_wp+1]);
HOD.M_cut = pow(10.0,a[j+1*n_wp+1]);
HOD.sigma_logM = pow(10.0,a[j+2*n_wp+1]);
mlo = set_low_mass();
printf("MLO %e %e %f\n",mlo,HOD.M_min,HOD.sigma_logM);
fflush(stdout);
while(!feof(fp))
{
fscanf(fp,"%d %d %e %e %e %e %e %e %e",
&i,&imass,&xh[0],&xh[1],&xh[2],&x1,&vh[0],&vh[1],&vh[2]);
mass=imass*RHO_CRIT*OMEGA_M*pow(0.3125,3.0);
if(mass<mlo)continue;
i1 = (int)(log10(mass)/0.1);
for(j=j_start;j<n_wp;++j)
{
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M1 = pow(10.0,a[j+0*n_wp+1]);
HOD.M_cut = pow(10.0,a[j+1*n_wp+1]);
HOD.sigma_logM = pow(10.0,a[j+2*n_wp+1]);
HOD.alpha = 1.0;
nsat = N_sat(mass);
if(j>j_start)
nc[j] = nc[j-1]+N_cen(mass);
else
nc[j] = N_cen(mass);
n1 = poisson_deviate(nsat);
for(i=1;i<=n1;++i)
{
NFW_position(mass,xg);
NFW_velocity(mass,vg);
for(k=0;k<3;++k)
{
xg[k]+=xh[k];
if(xg[k]<0)xg[k]+=BOX_SIZE;
if(xg[k]>BOX_SIZE)xg[k]-=BOX_SIZE;
vg[k]+=vh[k];
}
fprintf(fpa[j],"%e %e %e %e %e %e\n",xg[0],xg[1],xg[2],vg[0],vg[1],vg[2]);
mag = random_luminosity(wp_array[j].mag1,wp_array[j].mag2);
fprintf(fpb[j],"%f\n",mag);
/* random_2df_magnitude_error(); */
imag = (-mag-17)*2;
if(imag>=n_wp)imag = n_wp;
/* Bin up the galaxies by halo mass to check the HOD
*/
galcnt[imag][i1]++;
}
}
ncen=drand48();
for(j=j_start;j<n_wp;++j)
if(ncen<nc[j])break;
if(j<n_wp) {
fprintf(fpa[j],"%e %e %e %e %e %e %d\n",xh[0],xh[1],xh[2],vh[0],vh[1],vh[2],j);
mag = random_luminosity(wp_array[j].mag1,wp_array[j].mag2);
fprintf(fpb[j],"%f\n",mag);
imag = (-mag-17)*2;
if(imag>=n_wp)imag = n_wp;
galcnt[imag][i1]++;
}
halocnt[i1]++;
if(feof(fp))break;
}
/* output the binned HOD
*/
for(j=0;j<n_wp;++j)
{
sprintf(aa,"binned_HOD_err.%d",j);
fp2=fopen(aa,"w");
for(i=0;i<1000;++i)
if(galcnt[j][i]>0)
fprintf(fp2,"%d %f %f %d %d\n",
i,(i+0.5)*0.1,(float)galcnt[j][i]/halocnt[i],galcnt[j][i],halocnt[i]);
fclose(fp2);
}
exit(0);
}
double random_2df_magnitude_error()
{
double p,m;
p=0;
while(3*drand48()>p)
{
m = drand48()*2-1;
p = 1/RT2PI/.14*exp(-m*m/2/.14/.14)*0.7 +
1/RT2PI/.235*exp(-pow(log(1+m),2.0)/2/.235/.235)/(1+m)*0.3;
/*
printf("BUH %f %f %f %f\n",m,p,
1/RT2PI/.14*exp(-m*m/2/.14/14),1/RT2PI/.235*exp(-pow(log(1+m),2.0)/2/.235/.235)/(1+m)*0.3 );
*/
}
return(m);
}
void NFW_position(double mass, double x[])
{
double r,pr,max_p,costheta,sintheta,phi1,signs,rvir,rs,cvir;
cvir=halo_concentration(mass)*CVIR_FAC;
rvir=pow(3*mass/(4*DELTA_HALO*PI*RHO_CRIT*OMEGA_M),1.0/3.0);
rs=rvir/cvir;
max_p=NFW_density(rs,rs,1.0)*rs*rs*4.0*PI;
for(;;) {
r=drand48()*rvir;
pr=NFW_density(r,rs,1.0)*r*r*4.0*PI/max_p;
if(drand48()<=pr)
{
costheta=2.*(drand48()-.5);
sintheta=sqrt(1.-costheta*costheta);
signs=2.*(drand48()-.5);
costheta=signs*costheta/fabs(signs);
phi1=2.0*PI*drand48();
x[0]=r*sintheta*cos(phi1);
x[1]=r*sintheta*sin(phi1);
x[2]=r*costheta;
return;
}
}
}
/* This is the NFW density profile
*/
double NFW_density(double r, double rs, double ps)
{
return(ps*rs/(r*(1+r/rs)*(1+r/rs)));
}
/* This sets the velocity to be isotropic Gaussian.
*/
double NFW_velocity(double mass, double v[])
{
static long IDUM2=-455;
static double fac = -1;
double sigv;
int i;
if(fac<0)
fac=sqrt(4.499E-48)*pow(4*DELTA_HALO*PI*OMEGA_M*RHO_CRIT/3,1.0/6.0)*3.09E19;
sigv=fac*pow(mass,1.0/3.0)/sqrt(2.0);
for(i=0;i<3;++i)
v[i]=gasdev(&IDUM2)*sigv*VBIAS;
}
void dispersion_test(int n_wp, double *a, double mass)
{
int i,j,k,n=1,n1,n1t;
double ml_ratio,x1,nsat,ncen,err=0,mean=0,nc[10],mag;
for(i=1;i<=n;++i)
{
n1t =mean=err=x1=0;
for(j=0;j<n_wp;++j)
{
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M_low=wp_array[j].M_low;
HOD.M1 = wp_array[j].M1;
HOD.M_cut = wp_array[j].M_cut;
HOD.sigma_logM = wp_array[j].sigma_logM;
nsat = N_sat(mass);
if(j>0)
{
if(mass>HOD.M_low)
nc[j] = nc[j-1]+N_cen(mass);
else
nc[j] = nc[j-1]+0;
}
else
{
if(mass>HOD.M_low)
nc[j] = N_cen(mass);
else
nc[j] = 0;
}
n1t += n1 = poisson_deviate(nsat);
mag = random_luminosity(wp_array[j].mag1,wp_array[j].mag2);
x1 += (n1)*pow(10.0,-0.4*(mag-MAGSOL));
}
ncen=drand48();
for(j=0;j<n_wp;++j)
nc[j]/=nc[8];
for(j=0;j<n_wp;++j)
if(ncen<nc[j])break;
if(j>=n_wp)j=n_wp-1;
mag = random_luminosity(wp_array[j].mag1,wp_array[j].mag2);
x1+=pow(10.0,-0.4*(mag-MAGSOL));
if(x1>0) {
HOD.i_wp = j;
printf("LUM %f %f %f %d %f %e %d %e\n",
log10(x1),log10(mass),mass/x1,n1t,ncen,nc[j],j,N_cen(mass));
fflush(stdout);
}
}
}
double random_luminosity(double m1, double m2)
{
double p1,p2,pn,m;
pn = func_integrate_schechter(m1);
for(;;)
{
m = drand48()*(m2-m1)+m1;
p1 = func_integrate_schechter(m)/pn;
p2 = drand48();
if(p2<p1)break;
}
return(m);
}
int poisson_deviate(double nave)
{
static int flag=0;
double p,pp;
int n;
p=0;
pp=1;
if(!flag){
srand48(42398);
flag++;
}
while(p<pp)
{
if(nave<1)
n=(int)(drand48()*20);
else
n=(int)(drand48()*30*nave);
p=poisson_prob(n,nave);
pp=drand48();
}
return(n);
}
double poisson_prob(int n, double nave)
{
int i;
double fac=1;
if(n>0)
for(i=1;i<=n;++i)
fac*=nave/i;
return((float)(fac*exp(-nave)));
}
/*********************************************************************************/
/*********************************************************************************/
/*
* For outputting the fits and the data all in one file.
*
*/
/*********************************************************************************/
/*********************************************************************************/
void output_wp_fits(int n_wp, double *a)
{
int i,j,nr=50;
double *x[10],rlo=0.1,rhi=50,dlogr,r;
FILE *fp;
char fname[100];
dlogr = (log(rhi)-log(rlo))/(nr-1);
for(i=0;i<n_wp;++i)
x[i] = dvector(1,nr);
for(i=0;i<n_wp;++i)
{
HOD.i_wp=i;
HOD.M_low = wp_array[i].M_low;
HOD.M_min = wp_array[i].M_min;
HOD.M1 = wp_array[i].M1;
HOD.M_cut = wp_array[i].M_cut;
HOD.sigma_logM = wp_array[i].sigma_logM;
GALAXY_DENSITY=wp_array[i].ngal;
/*
printf("%d %e %e %e %e %e %e\n",
i,HOD.M_low,HOD.M_min,HOD.M1,HOD.alpha,HOD.M_cut,HOD.sigma_logM);
*/
RESET_FLAG_1H = RESET_FLAG_2H = 1;
for(j=1;j<=nr;++j)
{
r = exp(dlogr*(j-1))*rlo;
x[i][j] = projected_xi(r);
}
}
sprintf(fname,"%s.wp_fits",Task.root_filename);
fp=fopen(fname,"w");
for(j=1;j<=nr;++j)
{
fprintf(fp,"%e ",exp(dlogr*(j-1))*rlo);
for(i=0;i<n_wp;++i)
fprintf(fp,"%e ",x[i][j]);
fprintf(fp,"\n");
}
fclose(fp);
sprintf(fname,"%s.wp_data",Task.root_filename);
fp=fopen(fname,"w");
for(j=1;j<=wp_array[0].np;++j)
{
fprintf(fp,"%e ",wp_array[0].r[j]);
for(i=0;i<n_wp;++i)
fprintf(fp,"%e %e ",wp_array[i].x[j],wp_array[i].e[j]);
fprintf(fp,"\n");
}
fclose(fp);
}
/***********************************************************************************
*
* This is to make mocks based on the Gao effect. This requires that the
* local halo densities are read in from a file, and that M_min is
* recalculated such that the number density is help constant.
*
***********************************************************************************/
void ml_mock_catalog_gao(int n_wp, double *a)
{
FILE *fp,*fpa[9],*fp2,*fpb[9];
int i,j,k,n,imass,n1,j_start=0,i1,galcnt[9][1000],halocnt[1000],imag,nh,density_flag=1,ih;
double mass,xg[3],vg[3],nsat,nc[10],ncen,mlo,mag,err1,err2,
mmin_hi[10],mmin_lo[10],mmin_fac[10],*mmin_tmp;
float *hmass,*hden;
char aa[1000];
float x1,xh[3],vh[3];
long IDUM3 = -445;
srand48(4323432);
GAO_EFFECT = 1;
DENSITY_THRESHOLD = 2.0;
GAO.f0 = 2.0;
GAO.s0 = 2.0;
fprintf(stderr,"BEGINNING ML_MOCK GAO...\n");
fp=openfile("/gal/sr2/tinker/voids/LANL/halo.LANL.400");
fp2=openfile("/gal/sr2/tinker/voids/LANL/den.r5.LANL.400");
for(i=j_start;i<n_wp;++i)
{
sprintf(aa,"%s.mock.%d",Task.root_filename,i);
fpa[i] = fopen(aa,"w");
sprintf(aa,"%s.mag.%d",Task.root_filename,i);
fpb[i] = fopen(aa,"w");
}
j=2;
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M1 = pow(10.0,a[j+0*n_wp+1]);
HOD.M_cut = pow(10.0,a[j+1*n_wp+1]);
HOD.sigma_logM = pow(10.0,a[j+2*n_wp+1]);
mlo = set_low_mass();
printf("NCEN %e\n",N_avg(5.0e11));
/* Read in the masses and the densities.
*/
nh = filesize(fp);
hmass = vector(1,nh);
hden = vector(1,nh);
for(i=1;i<=nh;++i)
{
fscanf(fp,"%d %d %e %e %e %e %e %e %e",
&i,&imass,&xh[0],&xh[1],&xh[2],&x1,&vh[0],&vh[1],&vh[2]);
fscanf(fp2,"%e",&hden[i]);
mass=imass*RHO_CRIT*OMEGA_M*pow(0.3125,3.0);
hmass[i] = mass;
}
rewind(fp);
fclose(fp2);
for(i=0;i<1000;++i)
halocnt[i]=0;
for(j=0;j<n_wp;++j)
for(i=0;i<1000;++i)
galcnt[j][i]=0;
j=j_start;
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M1 = pow(10.0,a[j+0*n_wp+1]);
HOD.M_cut = pow(10.0,a[j+1*n_wp+1]);
HOD.sigma_logM = pow(10.0,a[j+2*n_wp+1]);
mlo = set_low_mass();
printf("MLO %e %e %f\n",mlo,HOD.M_min,HOD.sigma_logM);
fflush(stdout);
for(j=n_wp-1;j>=j_start;--j)
{
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M1 = pow(10.0,a[j+0*n_wp+1]);
HOD.M_cut = pow(10.0,a[j+1*n_wp+1]);
HOD.sigma_logM = pow(10.0,a[j+2*n_wp+1]);
GALAXY_DENSITY = wp_array[j].ngal;
dd_hod_functions(hmass,hden,nh);
mmin_lo[j] = HOD.M_min_loden;
mmin_hi[j] = HOD.M_min_loden*(1+GAO.f0*exp(-j*j/2.0/GAO.s0/GAO.s0));
wp_array[j].M_min = HOD.M_min_loden;
wp_array[j].M_low = HOD.M_low;
fprintf(stdout,"MMIN %d %e %e %e\n",j,mmin_lo[j],mmin_hi[j],HOD.M_low);
}
fflush(stdout);
for(ih=1;ih<=nh;++ih)
{
fscanf(fp,"%d %d %e %e %e %e %e %e %e",
&i,&imass,&xh[0],&xh[1],&xh[2],&x1,&vh[0],&vh[1],&vh[2]);
mass=imass*RHO_CRIT*OMEGA_M*pow(0.3125,3.0);
if(mass<mlo)continue;
i1 = (int)(log10(mass)/0.1);
/*
mass = 5e11;
hden[1] = 10;
hden[2] = 0;
*/
for(j=j_start;j<n_wp;++j)
{
HOD.i_wp=j;
HOD.M_min=wp_array[j].M_min;
HOD.M_low=wp_array[j].M_low;
HOD.M1 = pow(10.0,a[j+0*n_wp+1]);
HOD.M_cut = pow(10.0,a[j+1*n_wp+1]);
HOD.sigma_logM = pow(10.0,a[j+2*n_wp+1]);
if(hden[ih]>DENSITY_THRESHOLD)
{
if(j>j_start)
nc[j] = nc[j-1]+N_cen_hiden(mass);
else
nc[j] = N_cen_hiden(mass);
/* printf("NCEN %d %d %e %e\n",i,ih,N_cen_hiden(mass),nc[j]);*/
}
else
{
if(j>j_start)
nc[j] = nc[j-1]+N_cen(mass);
else
nc[j] = N_cen(mass);
/* printf("NCEN %d %d %e %e\n",i,ih,N_cen_hiden(mass),nc[j]);*/
}
nsat = N_sat(mass);
n1 = poisson_deviate(nsat);
for(i=1;i<=n1;++i)
{
NFW_position(mass,xg);
NFW_velocity(mass,vg);
for(k=0;k<3;++k)
{
xg[k]+=xh[k];
if(xg[k]<0)xg[k]+=BOX_SIZE;
if(xg[k]>BOX_SIZE)xg[k]-=BOX_SIZE;
vg[k]+=vh[k];
}
fprintf(fpa[j],"%e %e %e %e %e %e\n",xg[0],xg[1],xg[2],vg[0],vg[1],vg[2]);
mag = random_luminosity(wp_array[j].mag1,wp_array[j].mag2);
fprintf(fpb[j],"%f\n",mag);
imag = (-mag-17)*2;
if(imag>=n_wp)imag = n_wp;
/* Bin up the galaxies by halo mass to check the HOD
*/
galcnt[imag][i1]++;
}
}
density_flag=0;
ncen=drand48();
for(j=j_start;j<n_wp;++j)
if(ncen<nc[j])break;
if(j<n_wp) {
/*
if(j==4) {
for(k=j_start;k<n_wp;++k)
printf("BOO %d %d %e %e\n",i,k,nc[k],ncen);
}
*/
fprintf(fpa[j],"%e %e %e %e %e %e\n",xh[0],xh[1],xh[2],vh[0],vh[1],vh[2]);
mag = random_luminosity(wp_array[j].mag1,wp_array[j].mag2);
fprintf(fpb[j],"%f\n",mag);
imag = (-mag-17)*2;
if(imag>=n_wp)imag = n_wp;
galcnt[imag][i1]++;
}
halocnt[i1]++;
if(feof(fp))break;
}
/* output the binned HOD
*/
for(j=0;j<n_wp;++j)
{
sprintf(aa,"binned_HOD_err.%d",j);
fp2=fopen(aa,"w");
for(i=0;i<1000;++i)
if(galcnt[j][i]>0)
fprintf(fp2,"%d %f %f %d %d\n",
i,(i+0.5)*0.1,(float)galcnt[j][i]/halocnt[i],galcnt[j][i],halocnt[i]);
fclose(fp2);
}
exit(0);
}
double N_cen_gao(double m, int ii, double *mmin)
{
int i;
double n1,n2,n,logm,n2max=0;
logm=log10(m);
n1=0.5*(1+erf((logm - log10(HOD.M_min))/HOD.sigma_logM));
if(ii==8)return(n1);
for(n2=0,i=ii+1;i<9;++i)
if(m>wp_array[i].M_low)
{
n2=0.5*(1+erf((logm - log10(mmin[i]))/wp_array[i].sigma_logM));
if(n2>n2max)n2max=n2;
}
n = n1-n2max;
if(n<0)return(0);
return(n);
}
double N_cen_hiden(double m)
{
int i,ii;
double n1,n2,n,logm,n2max=0,f;
if(m<HOD.M_low)return(0);
ii = HOD.i_wp;
f=1 + GAO.f0*exp(-ii*ii/2.0/GAO.s0/GAO.s0);
logm=log10(m);
n1=0.5*(1+erf((logm - log10(HOD.M_min*f))/HOD.sigma_logM));
if(ii==8)return(n1);
for(n2=0,i=ii+1;i<9;++i)
if(m>wp_array[i].M_low)
{
f=1 + GAO.f0*exp(-i*i/2.0/GAO.s0/GAO.s0);
n2=0.5*(1+erf((logm - log10(wp_array[i].M_min*f))/wp_array[i].sigma_logM));
if(n2>n2max)n2max=n2;
}
n = n1-n2max;
if(n<0)return(0);
return(n);
}
/*************************************************************************************88
* This is to read in a chain and do some other calculations with it.
*/
void analyze_chain()
{
int i,j,k,n,n_wp=9,i1,i2;
FILE *fp,*fpa[9];
char aa[1000];
float x1,x2,x3,x4,x5,x6,x7;
/*
for(i=0;i<n_wp;++i)
{
sprintf(aa,"grep ACC%d out.560 > acc.%d",i,i);
fprintf(stderr,"%s\n",aa);
system(aa);
}
*/
for(i=0;i<n_wp;++i)
{
sprintf(aa,"acc.%d",i);
fpa[i] = openfile(aa);
}
n = filesize(fpa[0]);
for(i=1;i<=n;++i)
{
for(j=n_wp-1;j>=0;--j)
{
if(j==n_wp-1)RESET_COSMOLOGY++;
if(j>=6)
{
fscanf(fpa[j],"%4s %d %d %e %e %e %e %e %e %e",aa,&i1,&i2,&x1,&x2,&x3,&x4,&x5,&x6,&x7);
}
else
{
fscanf(fpa[j],"%4s %d %d %e %e %e %e %e %e",aa,&i1,&i2,&x1,&x2,&x4,&x5,&x6,&x7);
x3 = log10(0.2);
}
printf("%d %f %f %f %f %f %f %f\n",i,x1,x2,x3,x4,x5,x6,x7);
fgets(aa,1000,fpa[j]);
HOD.M1 = pow(10.0,x1);
HOD.M_cut = pow(10.0,x2);
HOD.sigma_logM = pow(10.0,x3);
HOD.alpha = x4;
CVIR_FAC = x5;
OMEGA_TEMP = x6;
SIGMA_8 = x7;
set_HOD_params();
wp_array[j].M_min = HOD.M_min;
wp_array[j].sigma_logM = HOD.sigma_logM;
wp_array[j].M1 = HOD.M1;
wp_array[j].M_low = HOD.M_low;
wp_array[j].M_cut = HOD.M_cut;
}
printf("SIGMA %f %f %f\n",SIGMA_8,OMEGA_TEMP,OMEGA_M);
fflush(stdout);
chi2_multiplicity_function(i);
printf("CHI2_CNOC %e\n",chi2_CNOC_data(i));
}
exit(0);
}