From d36791da74c588f8d4d8c261264db2e8c88cc735 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 20 Jan 2020 16:52:20 +0000 Subject: [PATCH 01/29] nq parameter --- src/common.c | 5 +++-- src/common.h | 3 ++- src/io.c | 6 +++--- 3 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/common.c b/src/common.c index a9543ea..1c4314e 100644 --- a/src/common.c +++ b/src/common.c @@ -546,7 +546,7 @@ void catalog_free(Catalog *cat) free(cat); } -HealpixShells *hp_shell_alloc(int nside,int nside_base,int nr) +HealpixShells *hp_shell_alloc(int nq, int nside,int nside_base,int nr) { if(nside>NSIDE_MAX_HPX) report_error(1,"Can't go beyond nside=%d\n",NSIDE_MAX_HPX); @@ -565,6 +565,7 @@ HealpixShells *hp_shell_alloc(int nside,int nside_base,int nr) for(ib=NodeThis;ibnq=nq; shell->nside=nside; shell->num_pix=nside_ratio*nside_ratio*nbeams_here; shell->listpix=my_malloc(shell->num_pix*sizeof(long)); @@ -588,7 +589,7 @@ HealpixShells *hp_shell_alloc(int nside,int nside_base,int nr) shell->rf=my_malloc(shell->nr*sizeof(flouble)); //Zero all data and clear - shell->data=my_calloc(shell->nr*shell->num_pix,sizeof(flouble)); + shell->data=my_calloc(shell->nq*shell->nr*shell->num_pix,sizeof(flouble)); shell->nadd=my_calloc(shell->nr*shell->num_pix,sizeof(int)); return shell; diff --git a/src/common.h b/src/common.h index 3c9a408..328de18 100644 --- a/src/common.h +++ b/src/common.h @@ -183,6 +183,7 @@ typedef struct { } Catalog; typedef struct { + int nq; //Number of quantities to map int nside; //Resolution parameter long num_pix; long *listpix; @@ -359,7 +360,7 @@ void end_rng(gsl_rng *rng); unsigned long long get_max_memory(ParamCoLoRe *par,int just_test); void get_radial_params(double rmax,int ngrid,int *nr,double *dr); //void get_random_angles(gsl_rng *rng,int ipix_nest,int ipix0,int nside,double *th,double *phi); -HealpixShells *hp_shell_alloc(int nside,int nside_base,int nr); +HealpixShells *hp_shell_alloc(int nq,int nside,int nside_base,int nr); void hp_shell_free(HealpixShells *shell); CatalogCartesian *catalog_cartesian_alloc(int nsrcs); void catalog_cartesian_free(CatalogCartesian *cat); diff --git a/src/io.c b/src/io.c index c748aa1..24e7ad3 100644 --- a/src/io.c +++ b/src/io.c @@ -463,16 +463,16 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) for(ii=0;iin_imap;ii++) { FILE *fnu=fopen(par->fnameNuImap[ii],"r"); if(fnu==NULL) error_open_file(par->fnameNuImap[ii]); - par->imap[ii]=hp_shell_alloc(par->nside_imap[ii],par->nside_base,linecount(fnu)); + par->imap[ii]=hp_shell_alloc(1,par->nside_imap[ii],par->nside_base,linecount(fnu)); fclose(fnu); } } if(par->do_kappa) - par->kmap=hp_shell_alloc(par->nside_kappa,par->nside_base,par->n_kappa); + par->kmap=hp_shell_alloc(1,par->nside_kappa,par->nside_base,par->n_kappa); if(par->do_isw) - par->pd_map=hp_shell_alloc(par->nside_isw,par->nside_base,par->n_isw); + par->pd_map=hp_shell_alloc(1,par->nside_isw,par->nside_base,par->n_isw); compute_tracer_cosmo(par); From 691c09bfba5915874b24ccdb12b119aa754fe935 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 20 Jan 2020 17:43:11 +0000 Subject: [PATCH 02/29] shear like kappa --- Makefile | 3 +- src/beaming.c | 30 ++++++++++------- src/common.c | 23 ++++++++++++- src/common.h | 17 ++++++++++ src/cosmo.c | 38 +++++++++++++++++++++- src/io.c | 90 ++++++++++++++++++++++++++++++++++++++++++++++++++- src/main.c | 40 ++++++++++++++--------- 7 files changed, 209 insertions(+), 32 deletions(-) diff --git a/Makefile b/Makefile index 45bb8cd..cba6053 100644 --- a/Makefile +++ b/Makefile @@ -120,6 +120,7 @@ DENSO = src/density.o SRCSO = src/srcs.o IMAPO = src/imap.o KAPPAO = src/kappa.o +SHEARO = src/shear.o ISWO = src/isw.o IOO = src/io.o HPIXO = src/healpix_extra.o @@ -127,7 +128,7 @@ BEAMO = src/beaming.o PREDICTO = src/predictions.o FFTLOGO = src/fftlog.o MAIN = src/main.c -OFILES = $(COMMONO) $(COSMOMADO) $(COSMOO) $(FOURIERO) $(DENSO) $(BEAMO) $(IOO) $(HPIXO) $(SRCSO) $(IMAPO) $(KAPPAO) $(ISWO) $(FFTLOGO) $(PREDICTO) +OFILES = $(COMMONO) $(COSMOMADO) $(COSMOO) $(FOURIERO) $(DENSO) $(BEAMO) $(IOO) $(HPIXO) $(SRCSO) $(IMAPO) $(KAPPAO) $(SHEARO) $(ISWO) $(FFTLOGO) $(PREDICTO) #OFILES = $(COMMONO) $(COSMOMADO) $(COSMOO) $(FOURIERO) $(DENSO) $(LCO) $(IOO) $(HPIXO) $(PIXO) $(PREDICTO) $(FFTLOGO) EXEC = CoLoRe diff --git a/src/beaming.c b/src/beaming.c index fd4e6f5..c172b7a 100644 --- a/src/beaming.c +++ b/src/beaming.c @@ -298,14 +298,16 @@ void get_beam_properties(ParamCoLoRe *par) return; } - if(par->do_srcs) - srcs_beams_preproc(par); - if(par->do_imap) - imap_beams_preproc(par); if(par->do_kappa) kappa_beams_preproc(par); + if(par->do_shear) + shear_beams_preproc(par); if(par->do_isw) isw_beams_preproc(par); + if(par->do_srcs) + srcs_beams_preproc(par); + if(par->do_imap) + imap_beams_preproc(par); if(NodeThis==0) timer(0); @@ -329,27 +331,31 @@ void get_beam_properties(ParamCoLoRe *par) par->nz_here=par->nz_all[node_i_am_now]; par->iz0_here=par->iz0_all[node_i_am_now]; - if(par->do_srcs) - srcs_get_beam_properties(par); - if(par->do_imap) - imap_get_beam_properties(par); if(par->do_kappa) kappa_get_beam_properties(par); + if(par->do_shear) + shear_get_beam_properties(par); if(par->do_isw) isw_get_beam_properties(par); + if(par->do_srcs) + srcs_get_beam_properties(par); + if(par->do_imap) + imap_get_beam_properties(par); } #ifdef _HAVE_MPI free(buffer_sr); #endif //_HAVE_MPI - if(par->do_srcs) - srcs_beams_postproc(par); - if(par->do_imap) - imap_beams_postproc(par); if(par->do_kappa) kappa_beams_postproc(par); + if(par->do_shear) + shear_beams_postproc(par); if(par->do_isw) isw_beams_postproc(par); + if(par->do_srcs) + srcs_beams_postproc(par); + if(par->do_imap) + imap_beams_postproc(par); if(NodeThis==0) timer(2); print_info("\n"); diff --git a/src/common.c b/src/common.c index 1c4314e..8cfe787 100644 --- a/src/common.c +++ b/src/common.c @@ -419,6 +419,19 @@ unsigned long long get_max_memory(ParamCoLoRe *par,int just_test) total_GB_kappa+=size_map*nr*(sizeof(flouble)+sizeof(int)); } + unsigned long long total_GB_shear=0; + if(par->do_shear) { + int ib; + int nr=par->n_shear; + int nbases=he_nside2npix(par->nside_base); + int nside_ratio=par->nside_shear/par->nside_base; + int npix_perbeam=nside_ratio*nside_ratio; + unsigned long long size_map=0; + for(ib=NodeThis;ibdo_isw) { int ib; @@ -431,7 +444,13 @@ unsigned long long get_max_memory(ParamCoLoRe *par,int just_test) size_map+=npix_perbeam; total_GB_isw+=size_map*nr*(sizeof(flouble)+sizeof(int)); } - total_GB=total_GB_gau+total_GB_lpt+total_GB_srcs+total_GB_imap+total_GB_kappa+total_GB_isw; + total_GB=total_GB_gau+ + total_GB_lpt+ + total_GB_srcs+ + total_GB_imap+ + total_GB_kappa+ + total_GB_shear+ + total_GB_isw; #ifdef _DEBUG int jj; @@ -447,6 +466,8 @@ unsigned long long get_max_memory(ParamCoLoRe *par,int just_test) printf(", %.3lf GB (imap)",(double)(total_GB_imap/pow(1024.,3))); if(par->do_kappa) printf(", %.3lf MB (kappa)",(double)(total_GB_kappa/pow(1024.,2))); + if(par->do_shear) + printf(", %.3lf MB (shear)",(double)(total_GB_shear/pow(1024.,2))); if(par->do_isw) printf(", %.3lf MB (isw)",(double)(total_GB_isw/pow(1024.,2))); printf("]\n"); diff --git a/src/common.h b/src/common.h index 328de18..b47978b 100644 --- a/src/common.h +++ b/src/common.h @@ -316,6 +316,12 @@ typedef struct { flouble **fl_mean_extra_kappa; flouble **cl_extra_kappa; #endif //_ADD_EXTRA_KAPPA + // Shear + int do_shear; //Do you want to output shear maps? + int n_shear; //How many maps? + double z_shear_out[NPLANES_MAX]; //Array of source plane redshifts + int nside_shear; + HealpixShells *smap; //Shear maps at each redshift // ISW int do_isw; //Do you want to output isw maps? int n_isw; //How many maps? @@ -400,6 +406,7 @@ void write_lpt(ParamCoLoRe *par,unsigned long long npart,flouble *x,flouble *y,f void write_srcs(ParamCoLoRe *par); void write_imap(ParamCoLoRe *par); void write_kappa(ParamCoLoRe *par); +void write_shear(ParamCoLoRe *par); void write_isw(ParamCoLoRe *par); void param_colore_free(ParamCoLoRe *par); @@ -463,6 +470,16 @@ void kappa_get_beam_properties(ParamCoLoRe *par); void kappa_beams_postproc(ParamCoLoRe *par); +////// +// Functions defined in shear.c +void shear_set_cartesian(ParamCoLoRe *par); +void shear_distribute(ParamCoLoRe *par); +void shear_get_local_properties(ParamCoLoRe *par); +void shear_beams_preproc(ParamCoLoRe *par); +void shear_get_beam_properties(ParamCoLoRe *par); +void shear_beams_postproc(ParamCoLoRe *par); + + ////// // Functions defined in isw.c void isw_set_cartesian(ParamCoLoRe *par); diff --git a/src/cosmo.c b/src/cosmo.c index f0f01e5..badbfa0 100644 --- a/src/cosmo.c +++ b/src/cosmo.c @@ -24,7 +24,8 @@ #include #define CL_TYPE_KAPPA 0 -#define CL_TYPE_ISW 1 +#define CL_TYPE_SHEAR 1 +#define CL_TYPE_ISW 2 static double f_of_r_linear(ParamCoLoRe *par,double r,double *f,double f0,double ff) { @@ -127,6 +128,20 @@ static double window_kappa_limber(ParamCoLoRe *par,int l,double k,double chi_0,d } } +static double window_shear_limber(ParamCoLoRe *par,int l,double k,double chi_0,double chi_f,double chi_s) +{ + double chi_l=(l+0.5)/k; + + if((chi_l<=0) || (chi_lchi_f)) + return 0; + else { + double gf=get_bg(par,chi_l,BG_D1,0); + double aa=1/(1+get_bg(par,chi_l,BG_Z,0)); + + return par->prefac_lensing*sqrt((l-1)*l*(l+1.)*(l+2.))*gf*(chi_s-chi_l)/(k*k*chi_s*chi_l*aa); + } +} + static double window_isw_limber(ParamCoLoRe *par,int l,double k,double chi_0,double chi_f,double chi_s) { double chi_l=(l+0.5)/k; @@ -147,6 +162,10 @@ static double cl_integrand(double lk,void *params) wa=window_kappa_limber(p->par,p->l,k,p->chi_0_a,p->chi_f_a,p->chi_s); wb=window_kappa_limber(p->par,p->l,k,p->chi_0_b,p->chi_f_b,p->chi_s); } + else if(p->cl_type==CL_TYPE_SHEAR) { + wa=window_shear_limber(p->par,p->l,k,p->chi_0_a,p->chi_f_a,p->chi_s); + wb=window_shear_limber(p->par,p->l,k,p->chi_0_b,p->chi_f_b,p->chi_s); + } else if(p->cl_type==CL_TYPE_ISW) { wa=window_isw_limber(p->par,p->l,k,p->chi_0_a,p->chi_f_a,p->chi_s); wb=window_isw_limber(p->par,p->l,k,p->chi_0_b,p->chi_f_b,p->chi_s); @@ -693,6 +712,15 @@ void cosmo_set(ParamCoLoRe *par) } } + if(par->do_shear) { + for(ii=0;iin_shear;ii++) { + double z=par->z_shear_out[ii]; + if(z>par->z_max) { + report_error(1,"Source plane %d outside redshift range\n",ii+1); + } + } + } + if(par->do_isw) { for(ii=0;iin_isw;ii++) { double z=par->z_isw_out[ii]; @@ -758,6 +786,14 @@ void compute_tracer_cosmo(ParamCoLoRe *par) par->kmap->rf[ii]=csm_radial_comoving_distance(pars,1./(1+z)); } } + if(par->do_shear) { + int ii; + for(ii=0;iin_shear;ii++) { + double z=par->z_shear_out[ii]; + par->kmap->r0[ii]=csm_radial_comoving_distance(pars,1./(1+z)); + par->kmap->rf[ii]=csm_radial_comoving_distance(pars,1./(1+z)); + } + } if(par->do_isw) { int ii; for(ii=0;iin_isw;ii++) { diff --git a/src/io.c b/src/io.c index 24e7ad3..3b4973b 100644 --- a/src/io.c +++ b/src/io.c @@ -97,12 +97,15 @@ static ParamCoLoRe *param_colore_new(void) par->do_isw=0; par->do_imap=0; par->do_kappa=0; + par->do_shear=0; par->do_isw=0; par->n_srcs=-1; par->n_imap=-1; par->n_kappa=-1; + par->n_shear=-1; par->n_isw=-1; par->nside_kappa=-1; + par->nside_shear=-1; par->nside_isw=-1; for(ii=0;iifnameBzSrcs[ii],"default"); @@ -124,6 +127,7 @@ static ParamCoLoRe *param_colore_new(void) } for(ii=0;iiz_kappa_out[ii]=-1; + par->z_shear_out[ii]=-1; par->z_isw_out[ii]=-1; } par->cats_c=NULL; @@ -132,6 +136,7 @@ static ParamCoLoRe *param_colore_new(void) par->nsources_this=NULL; par->imap=NULL; par->kmap=NULL; + par->smap=NULL; par->pd_map=NULL; //Beam distribution @@ -359,6 +364,15 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) conf_read_int(conf,"kappa","nside",&(par->nside_kappa)); } + //Shear maps + cset=config_lookup(conf,"shear"); + if(cset!=NULL) { + par->do_shear=1; + par->do_lensing=1; + conf_read_double_array(conf,"shear","z_out",par->z_shear_out,&(par->n_shear),NPLANES_MAX); + conf_read_int(conf,"shear","nside",&(par->nside_shear)); + } + cset=config_lookup(conf,"isw"); if(cset!=NULL) { par->do_isw=1; @@ -398,7 +412,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) else par->do_smoothing=0; - par->need_beaming=par->do_lensing+par->do_kappa+par->do_isw+par->do_skewers; + par->need_beaming=par->do_lensing+par->do_kappa+par->do_shear+par->do_isw+par->do_skewers; init_fftw(par); get_max_memory(par,test_memory+par->just_do_pred); @@ -426,6 +440,8 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) print_info(" %d intensity mapping species\n",par->n_imap); if(par->do_kappa) print_info(" %d lensing source planes\n",par->n_kappa); + if(par->do_shear) + print_info(" %d lensing source planes\n",par->n_shear); if(par->do_isw) print_info(" %d ISW source planes\n",par->n_isw); if(par->do_lensing) @@ -471,6 +487,9 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) if(par->do_kappa) par->kmap=hp_shell_alloc(1,par->nside_kappa,par->nside_base,par->n_kappa); + if(par->do_shear) + par->smap=hp_shell_alloc(2,par->nside_shear,par->nside_base,par->n_shear); + if(par->do_isw) par->pd_map=hp_shell_alloc(1,par->nside_isw,par->nside_base,par->n_isw); @@ -741,6 +760,70 @@ void write_kappa(ParamCoLoRe *par) print_info("\n"); } +void write_shear(ParamCoLoRe *par) +{ + int ir; + char fname[256]; + long npx=he_nside2npix(par->nside_shear); + flouble **map_write=my_malloc(2*sizeof(flouble *)); + int *map_nadd=my_malloc(npx*sizeof(int)); + map_write[0]=my_malloc(npx*sizeof(flouble)); + map_write[1]=my_malloc(npx*sizeof(flouble)); + if(NodeThis==0) timer(0); + print_info("*** Writing shear source maps\n"); + for(ir=0;irsmap->nr;ir++) { + long ip; + long ir_t=ir*par->smap->num_pix; + + //Write local pixels to dummy map + for(ip=0;ipprefixOut,ir); + for(ip=0;ipsmap->num_pix;ip++) { + int id_pix=par->smap->listpix[ip]; + map_write[0][id_pix]+=par->smap->data[0 + 2*(ip + ir_t)]; + map_write[1][id_pix]+=par->smap->data[1 + 2*(ip + ir_t)]; + map_nadd[id_pix]+=par->smap->nadd[ir_t+ip]; + } + + //Collect all dummy maps +#ifdef _HAVE_MPI + if(NodeThis==0) + MPI_Reduce(MPI_IN_PLACE,map_write[0],npx,FLOUBLE_MPI,MPI_SUM,0,MPI_COMM_WORLD); + else + MPI_Reduce(map_write[0],NULL ,npx,FLOUBLE_MPI,MPI_SUM,0,MPI_COMM_WORLD); + if(NodeThis==0) + MPI_Reduce(MPI_IN_PLACE,map_write[1],npx,FLOUBLE_MPI,MPI_SUM,0,MPI_COMM_WORLD); + else + MPI_Reduce(map_write[1],NULL ,npx,FLOUBLE_MPI,MPI_SUM,0,MPI_COMM_WORLD); + if(NodeThis==0) + MPI_Reduce(MPI_IN_PLACE,map_nadd,npx,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD); + else + MPI_Reduce(map_nadd ,NULL ,npx,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD); +#endif //_HAVE_MPI + + for(ip=0;ip0) { + map_write[0][ip]/=map_nadd[ip]; + map_write[1][ip]/=map_nadd[ip]; + } + } + + //Write dummy map + if(NodeThis==0) + he_write_healpix_map(map_write,2,par->nside_shear,fname,1); + } + free(map_write[0]); + free(map_write[1]); + free(map_write); + free(map_nadd); + if(NodeThis==0) timer(2); + print_info("\n"); +} + void write_isw(ParamCoLoRe *par) { int ir; @@ -1082,6 +1165,11 @@ void param_colore_free(ParamCoLoRe *par) #endif //_ADD_EXTRA_KAPPA } + if(par->do_shear) { + if(par->smap!=NULL) + hp_shell_free(par->smap); + } + if(par->do_isw) { if(par->pd_map!=NULL) hp_shell_free(par->pd_map); diff --git a/src/main.c b/src/main.c index 6ea08dd..59ce82f 100644 --- a/src/main.c +++ b/src/main.c @@ -72,34 +72,40 @@ int main(int argc,char **argv) compute_density_normalization(par); //Get information from slabs - if(par->do_srcs) - srcs_set_cartesian(par); - if(par->do_imap) - imap_set_cartesian(par); if(par->do_kappa) kappa_set_cartesian(par); + if(par->do_shear) + shear_set_cartesian(par); if(par->do_isw) isw_set_cartesian(par); - - //Distribute information across if(par->do_srcs) - srcs_distribute(par); + srcs_set_cartesian(par); if(par->do_imap) - imap_distribute(par); + imap_set_cartesian(par); + + //Distribute information across if(par->do_kappa) kappa_distribute(par); + if(par->do_shear) + shear_distribute(par); if(par->do_isw) isw_distribute(par); - - //Postprocess after if(par->do_srcs) - srcs_get_local_properties(par); + srcs_distribute(par); if(par->do_imap) - imap_get_local_properties(par); + imap_distribute(par); + + //Postprocess after if(par->do_kappa) kappa_get_local_properties(par); + if(par->do_shear) + shear_get_local_properties(par); if(par->do_isw) isw_get_local_properties(par); + if(par->do_srcs) + srcs_get_local_properties(par); + if(par->do_imap) + imap_get_local_properties(par); //All-to-all communication of density field //and computation of all required quantities @@ -107,14 +113,16 @@ int main(int argc,char **argv) get_beam_properties(par); //Write output - if(par->do_srcs) - write_srcs(par); - if(par->do_imap) - write_imap(par); if(par->do_kappa) write_kappa(par); + if(par->do_shear) + write_shear(par); if(par->do_isw) write_isw(par); + if(par->do_srcs) + write_srcs(par); + if(par->do_imap) + write_imap(par); } print_info("\n"); From 95faa5e78cb5eb727a78257aa801c9fc4fef4f83 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 20 Jan 2020 18:38:35 +0000 Subject: [PATCH 03/29] fixed bug --- src/cosmo.c | 4 +- src/shear.c | 194 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 196 insertions(+), 2 deletions(-) create mode 100644 src/shear.c diff --git a/src/cosmo.c b/src/cosmo.c index badbfa0..1e41ea8 100644 --- a/src/cosmo.c +++ b/src/cosmo.c @@ -790,8 +790,8 @@ void compute_tracer_cosmo(ParamCoLoRe *par) int ii; for(ii=0;iin_shear;ii++) { double z=par->z_shear_out[ii]; - par->kmap->r0[ii]=csm_radial_comoving_distance(pars,1./(1+z)); - par->kmap->rf[ii]=csm_radial_comoving_distance(pars,1./(1+z)); + par->smap->r0[ii]=csm_radial_comoving_distance(pars,1./(1+z)); + par->smap->rf[ii]=csm_radial_comoving_distance(pars,1./(1+z)); } } if(par->do_isw) { diff --git a/src/shear.c b/src/shear.c new file mode 100644 index 0000000..7e009f6 --- /dev/null +++ b/src/shear.c @@ -0,0 +1,194 @@ +/////////////////////////////////////////////////////////////////////// +// // +// Copyright 2012 David Alonso // +// // +// // +// This file is part of CoLoRe. // +// // +// CoLoRe 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. // +// // +// CoLoRe 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 CoLoRe. If not, see . // +// // +/////////////////////////////////////////////////////////////////////// +#include "common.h" + +void shear_set_cartesian(ParamCoLoRe *par) +{ + return; +} + +void shear_get_local_properties(ParamCoLoRe *par) +{ + return; +} + +void shear_distribute(ParamCoLoRe *par) +{ + return; +} + +void shear_beams_preproc(ParamCoLoRe *par) +{ + //Sort radii in ascending order + int ir; + flouble *r0_i=my_malloc(par->smap->nr*sizeof(flouble)); + flouble *rf_i=my_malloc(par->smap->nr*sizeof(flouble)); + + memcpy(r0_i,par->smap->r0,par->smap->nr*sizeof(flouble)); + memcpy(rf_i,par->smap->rf,par->smap->nr*sizeof(flouble)); + + int *i_sorted=ind_sort(par->smap->nr,par->smap->r0); + for(ir=0;irsmap->nr;ir++) { + par->smap->r0[ir]=r0_i[i_sorted[ir]]; + par->smap->rf[ir]=rf_i[i_sorted[ir]]; + } + free(r0_i); + free(rf_i); + free(i_sorted); + + //Zero all data and clear +#ifdef _HAVE_OMP +#pragma omp parallel default(none) \ + shared(par) +#endif //_HAVE_OMP + { + long ipp; + +#ifdef _HAVE_OMP +#pragma omp for +#endif //_HAVE_OMP + for(ipp=0;ippsmap->num_pix*par->smap->nr;ipp++) { + par->smap->data[ipp]=0; + par->smap->nadd[ipp]=1; + } //end omp for + } //end omp parallel + + return; +} + +void shear_get_beam_properties(ParamCoLoRe *par) +{ + HealpixShells *smap=par->smap; + +#ifdef _HAVE_OMP +#pragma omp parallel default(none) \ + shared(par,smap) +#endif //_HAVE_OMP + { + double idx=par->n_grid/par->l_box; + + //Setup radial decomposition + int i_r,nr; + double idr,dr; + get_radial_params(par->r_max,par->n_grid,&nr,&dr); + idr=1./dr; + + //Compute index of each source plane + int *i_r_max_arr=my_malloc(smap->nr*sizeof(int)); + int *i_r_min_arr=my_malloc(smap->nr*sizeof(int)); + double *inv_r_max=my_malloc(smap->nr*sizeof(double)); + for(i_r=0;i_rnr;i_r++) { + int i_r_here=(int)(smap->rf[i_r]*idr+0.5); + inv_r_max[i_r]=1./(i_r_here*dr); + i_r_max_arr[i_r]=MIN(i_r_here,nr-1); + } + i_r_min_arr[0]=0; + for(i_r=1;i_rnr;i_r++) + i_r_min_arr[i_r]=i_r_max_arr[i_r-1]+1; + + //Fill up integral kernels + double *fac_r_1=my_malloc(nr*sizeof(double)); + double *fac_r_2=my_malloc(nr*sizeof(double)); + for(i_r=0;i_rnum_pix;ip++) { + int ax,added; + flouble t[6]; + double r1[6],r2[6],xn[3]; + double shear1_1=0,shear1_2=0; + double shear2_1=0,shear2_2=0; + double *u=&(smap->pos[3*ip]); + double prefac=idx*idx; + double cth_h=1,sth_h=0,cph_h=1,sph_h=0; + + cth_h=u[2]; + if(cth_h>=1) cth_h=1; + if(cth_h<=-1) cth_h=-1; + sth_h=sqrt((1-cth_h)*(1+cth_h)); + if(sth_h!=0) { + cph_h=u[0]/sth_h; + sph_h=u[1]/sth_h; + } + + r1[0]=(cth_h*cth_h*cph_h*cph_h-sph_h*sph_h)*prefac; + r1[1]=(2*cph_h*sph_h*(cth_h*cth_h+1))*prefac; + r1[2]=(-2*cth_h*sth_h*cph_h)*prefac; + r1[3]=(cth_h*cth_h*sph_h*sph_h-cph_h*cph_h)*prefac; + r1[4]=(-2*cth_h*sth_h*sph_h)*prefac; + r1[5]=(sth_h*sth_h)*prefac; + + r2[0]=(-2*cth_h*cph_h*sph_h)*prefac; + r2[1]=(2*cth_h*(cph_h*cph_h-sph_h*sph_h))*prefac; + r2[2]=(2*sth_h*sph_h)*prefac; + r2[3]=(2*cth_h*sph_h*cph_h)*prefac; + r2[4]=(-2*sth_h*cph_h)*prefac; + r2[5]=0; + for(i_r=0;i_rnr;i_r++) { + int irr; + int irmin=i_r_min_arr[i_r]; + int irmax=i_r_max_arr[i_r]; + for(irr=irmin;irr<=irmax;irr++) { + double rm=(irr+0.5)*dr; + for(ax=0;ax<3;ax++) + xn[ax]=(rm*u[ax]+par->pos_obs[ax])*idx; + added=interpolate_from_grid(par,xn,NULL,NULL,t,NULL,NULL,RETURN_TID,INTERP_TYPE_SHEAR); + if(added) { + double dotp1=0,dotp2=0; + for(ax=0;ax<6;ax++) { + dotp1+=r1[ax]*t[ax]; + dotp2+=r2[ax]*t[ax]; + } + shear1_1+=dotp1*fac_r_1[irr]; + shear1_2+=dotp1*fac_r_2[irr]; + shear2_1+=dotp2*fac_r_1[irr]; + shear2_2+=dotp2*fac_r_2[irr]; + } + } + smap->data[2*(i_r*smap->num_pix+ip)+0]+=(shear1_1-inv_r_max[i_r]*shear1_2); + smap->data[2*(i_r*smap->num_pix+ip)+1]+=(shear2_1-inv_r_max[i_r]*shear2_2); + } + } //end omp for + + free(fac_r_1); + free(fac_r_2); + free(i_r_max_arr); + free(i_r_min_arr); + free(inv_r_max); + } //end omp parallel + + return; +} + +void shear_beams_postproc(ParamCoLoRe *par) +{ + //Set rf to end of cell + return; +} From c8711179a0c5f5678164849729940b4a94141f2a Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 20 Jan 2020 19:34:00 +0000 Subject: [PATCH 04/29] dr_shear --- src/common.h | 6 ++++-- src/cosmo.c | 32 +++++++++++++++----------------- src/io.c | 46 +++++++++++++++++++++++++++++++--------------- src/srcs.c | 6 +++--- 4 files changed, 53 insertions(+), 37 deletions(-) diff --git a/src/common.h b/src/common.h index b47978b..2e82439 100644 --- a/src/common.h +++ b/src/common.h @@ -275,7 +275,7 @@ typedef struct { // Sources int do_srcs; //Do we include sources? int do_skewers; //Do we include skewer information? - int do_lensing; //Do we need to compute the lensing potential? + int do_srcs_shear; //Do we need to compute the lensing potential? int n_srcs; //Number of source types char fnameBzSrcs[NPOP_MAX][256]; //Files containing b(z) for each source type char fnameNzSrcs[NPOP_MAX][256]; //Files containing dN/dzdOmega (in deg^-2) @@ -318,8 +318,9 @@ typedef struct { #endif //_ADD_EXTRA_KAPPA // Shear int do_shear; //Do you want to output shear maps? + double dr_shear; //Separation between shear maps. + double idr_shear; //1/dr_shear int n_shear; //How many maps? - double z_shear_out[NPLANES_MAX]; //Array of source plane redshifts int nside_shear; HealpixShells *smap; //Shear maps at each redshift // ISW @@ -397,6 +398,7 @@ void cosmo_set(ParamCoLoRe *par); double r_of_z(ParamCoLoRe *par,double z); double get_bg(ParamCoLoRe *par,double r,int tag,int ipop); void compute_tracer_cosmo(ParamCoLoRe *par); +flouble *compute_shear_spacing(ParamCoLoRe *par); ////// // Functions defined in io.c diff --git a/src/cosmo.c b/src/cosmo.c index 1e41ea8..b6d4411 100644 --- a/src/cosmo.c +++ b/src/cosmo.c @@ -712,15 +712,6 @@ void cosmo_set(ParamCoLoRe *par) } } - if(par->do_shear) { - for(ii=0;iin_shear;ii++) { - double z=par->z_shear_out[ii]; - if(z>par->z_max) { - report_error(1,"Source plane %d outside redshift range\n",ii+1); - } - } - } - if(par->do_isw) { for(ii=0;iin_isw;ii++) { double z=par->z_isw_out[ii]; @@ -786,14 +777,6 @@ void compute_tracer_cosmo(ParamCoLoRe *par) par->kmap->rf[ii]=csm_radial_comoving_distance(pars,1./(1+z)); } } - if(par->do_shear) { - int ii; - for(ii=0;iin_shear;ii++) { - double z=par->z_shear_out[ii]; - par->smap->r0[ii]=csm_radial_comoving_distance(pars,1./(1+z)); - par->smap->rf[ii]=csm_radial_comoving_distance(pars,1./(1+z)); - } - } if(par->do_isw) { int ii; for(ii=0;iin_isw;ii++) { @@ -805,3 +788,18 @@ void compute_tracer_cosmo(ParamCoLoRe *par) csm_params_free(pars); } + +flouble *compute_shear_spacing(ParamCoLoRe *par) +{ + //Figure out appropriate radial sampling + flouble *rarr; + int ir,nr=(int)(par->r_max/par->dr_shear+1.); + par->dr_shear=par->r_max/nr; + par->idr_shear=1./par->dr_shear; + rarr=my_malloc(nr*sizeof(flouble)); + for(ir=0;irdr_shear; + + par->n_shear=nr; + return rarr; +} diff --git a/src/io.c b/src/io.c index 3b4973b..13c0fe3 100644 --- a/src/io.c +++ b/src/io.c @@ -93,11 +93,13 @@ static ParamCoLoRe *param_colore_new(void) //Tracers par->do_srcs=0; par->do_skewers=0; - par->do_lensing=0; + par->do_srcs_shear=0; par->do_isw=0; par->do_imap=0; par->do_kappa=0; par->do_shear=0; + par->dr_shear=100000.; + par->idr_shear=1E-5; par->do_isw=0; par->n_srcs=-1; par->n_imap=-1; @@ -127,7 +129,6 @@ static ParamCoLoRe *param_colore_new(void) } for(ii=0;iiz_kappa_out[ii]=-1; - par->z_shear_out[ii]=-1; par->z_isw_out[ii]=-1; } par->cats_c=NULL; @@ -320,7 +321,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) conf_read_string(conf,c_dum,"bias_filename",par->fnameBzSrcs[ii]); conf_read_bool(conf,c_dum,"include_shear",&(par->shear_srcs[ii])); if(par->shear_srcs[ii]) - par->do_lensing=1; + par->do_srcs_shear=1; conf_read_bool(conf,c_dum,"store_skewers",&(par->skw_srcs[ii])); if(par->skw_srcs[ii]) { par->do_skewers=1; @@ -359,7 +360,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) cset=config_lookup(conf,"kappa"); if(cset!=NULL) { par->do_kappa=1; - par->do_lensing=1; + par->do_srcs_shear=1; conf_read_double_array(conf,"kappa","z_out",par->z_kappa_out,&(par->n_kappa),NPLANES_MAX); conf_read_int(conf,"kappa","nside",&(par->nside_kappa)); } @@ -368,10 +369,15 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) cset=config_lookup(conf,"shear"); if(cset!=NULL) { par->do_shear=1; - par->do_lensing=1; - conf_read_double_array(conf,"shear","z_out",par->z_shear_out,&(par->n_shear),NPLANES_MAX); + par->do_srcs_shear=1; + conf_read_double(conf,"shear","dr_shear",&(par->dr_shear)); conf_read_int(conf,"shear","nside",&(par->nside_shear)); + if(par->dr_shear<=par->l_box/par->n_grid) + report_error(1,"Spacing between shear source planes is smaller than box resolution\n"); } + // Check shear exists if requested with catalog + if(par->do_srcs_shear && !(par->do_shear)) + report_error(1,"Include a \"shear\" section if you want shear with your galaxies\n"); cset=config_lookup(conf,"isw"); if(cset!=NULL) { @@ -412,7 +418,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) else par->do_smoothing=0; - par->need_beaming=par->do_lensing+par->do_kappa+par->do_shear+par->do_isw+par->do_skewers; + par->need_beaming=par->do_srcs_shear+par->do_kappa+par->do_shear+par->do_isw+par->do_skewers; init_fftw(par); get_max_memory(par,test_memory+par->just_do_pred); @@ -444,8 +450,8 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) print_info(" %d lensing source planes\n",par->n_shear); if(par->do_isw) print_info(" %d ISW source planes\n",par->n_isw); - if(par->do_lensing) - print_info(" Will include lensing shear\n"); + if(par->do_srcs_shear) + print_info(" Will include lensing shear in source catalog\n"); if(!par->need_beaming) print_info(" Will NOT need to all-to-all communicate fields\n"); print_info("\n"); @@ -487,8 +493,18 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) if(par->do_kappa) par->kmap=hp_shell_alloc(1,par->nside_kappa,par->nside_base,par->n_kappa); - if(par->do_shear) + if(par->do_shear) { + flouble *r_arr=compute_shear_spacing(par); par->smap=hp_shell_alloc(2,par->nside_shear,par->nside_base,par->n_shear); + int i_r; + for(i_r=0;i_rn_shear;i_r++) { + par->smap->r0[i_r]=r_arr[i_r]; + par->smap->rf[i_r]=r_arr[i_r]; + printf("%d %lE %lE\n",i_r,r_arr[i_r], + get_bg(par,r_arr[i_r],BG_Z,-1)); + } + free(r_arr); + } if(par->do_isw) par->pd_map=hp_shell_alloc(1,par->nside_isw,par->nside_base,par->n_isw); @@ -951,7 +967,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) char *ttype[]={"TYPE","RA" ,"DEC","Z_COSMO","DZ_RSD","E1","E2"}; char *tform[]={"1J" ,"1E" ,"1E" ,"1E" ,"1E" ,"1E","1E"}; char *tunit[]={"NA" ,"DEG","DEG","NA" ,"NA" ,"NA","NA"}; - if(par->do_lensing) + if(par->do_srcs_shear) nfields=7; print_info(" %d-th population (FITS)\n",ipop); @@ -984,7 +1000,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) dec_arr[ii]=par->cats[ipop]->srcs[row_here+ii].dec; z0_arr[ii]=par->cats[ipop]->srcs[row_here+ii].z0; rsd_arr[ii]=par->cats[ipop]->srcs[row_here+ii].dz_rsd; - if(par->do_lensing) { + if(par->do_srcs_shear) { e1_arr[ii]=par->cats[ipop]->srcs[row_here+ii].e1; e2_arr[ii]=par->cats[ipop]->srcs[row_here+ii].e2; } @@ -994,7 +1010,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) fits_write_col(fptr,TFLOAT,3,row_here+1,1,nrw_here,dec_arr,&status); fits_write_col(fptr,TFLOAT,4,row_here+1,1,nrw_here,z0_arr,&status); fits_write_col(fptr,TFLOAT,5,row_here+1,1,nrw_here,rsd_arr,&status); - if(par->do_lensing) { + if(par->do_srcs_shear) { fits_write_col(fptr,TFLOAT,6,row_here+1,1,nrw_here,e1_arr,&status); fits_write_col(fptr,TFLOAT,7,row_here+1,1,nrw_here,e2_arr,&status); } @@ -1077,7 +1093,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) FILE *fil=fopen(fname,"w"); if(fil==NULL) error_open_file(fname); fprintf(fil,"#[1]type [2]RA, [3]dec, [4]z0, [5]dz_RSD "); - if(par->do_lensing) + if(par->do_srcs_shear) fprintf(fil,"#[6]e1, [7]e2\n"); else fprintf(fil,"\n"); @@ -1085,7 +1101,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) fprintf(fil,"%d %E %E %E %E ", ipop,par->cats[ipop]->srcs[jj].ra,par->cats[ipop]->srcs[jj].dec, par->cats[ipop]->srcs[jj].z0,par->cats[ipop]->srcs[jj].dz_rsd); - if(par->do_lensing) + if(par->do_srcs_shear) fprintf(fil,"%E %E \n",par->cats[ipop]->srcs[jj].e1,par->cats[ipop]->srcs[jj].e2); else fprintf(fil,"\n"); diff --git a/src/srcs.c b/src/srcs.c index f2d1d10..3373b04 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -416,7 +416,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) double *fac_r_1=NULL,*fac_r_2=NULL; //Kernels for the LOS integrals - if(par->do_lensing) { + if(par->do_srcs_shear) { fac_r_1=my_malloc(cat->nr*sizeof(double)); fac_r_2=my_malloc(cat->nr*sizeof(double)); @@ -479,7 +479,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) } //Compute lensing shear - if(par->do_lensing) { + if(par->do_srcs_shear) { //Compute linear transformations needed for shear double r1[6],r2[6]; double cth_h=1,sth_h=0,cph_h=1,sph_h=0; @@ -532,7 +532,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) cat->srcs[ip].e2+=e2; } }//end omp for - if(par->do_lensing) { + if(par->do_srcs_shear) { free(fac_r_1); free(fac_r_2); } From b1a9c25b73a0a50348b720f38e9c50298a11bd1c Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 20 Jan 2020 19:45:28 +0000 Subject: [PATCH 05/29] recomputed r_arr --- src/common.h | 1 - src/cosmo.c | 1 - src/io.c | 1 - src/shear.c | 9 +++++++++ 4 files changed, 9 insertions(+), 3 deletions(-) diff --git a/src/common.h b/src/common.h index 2e82439..c1958b2 100644 --- a/src/common.h +++ b/src/common.h @@ -319,7 +319,6 @@ typedef struct { // Shear int do_shear; //Do you want to output shear maps? double dr_shear; //Separation between shear maps. - double idr_shear; //1/dr_shear int n_shear; //How many maps? int nside_shear; HealpixShells *smap; //Shear maps at each redshift diff --git a/src/cosmo.c b/src/cosmo.c index b6d4411..6cf7cec 100644 --- a/src/cosmo.c +++ b/src/cosmo.c @@ -795,7 +795,6 @@ flouble *compute_shear_spacing(ParamCoLoRe *par) flouble *rarr; int ir,nr=(int)(par->r_max/par->dr_shear+1.); par->dr_shear=par->r_max/nr; - par->idr_shear=1./par->dr_shear; rarr=my_malloc(nr*sizeof(flouble)); for(ir=0;irdr_shear; diff --git a/src/io.c b/src/io.c index 13c0fe3..056d893 100644 --- a/src/io.c +++ b/src/io.c @@ -99,7 +99,6 @@ static ParamCoLoRe *param_colore_new(void) par->do_kappa=0; par->do_shear=0; par->dr_shear=100000.; - par->idr_shear=1E-5; par->do_isw=0; par->n_srcs=-1; par->n_imap=-1; diff --git a/src/shear.c b/src/shear.c index 7e009f6..12ce826 100644 --- a/src/shear.c +++ b/src/shear.c @@ -101,6 +101,7 @@ void shear_get_beam_properties(ParamCoLoRe *par) inv_r_max[i_r]=1./(i_r_here*dr); i_r_max_arr[i_r]=MIN(i_r_here,nr-1); } + i_r_min_arr[0]=0; for(i_r=1;i_rnr;i_r++) i_r_min_arr[i_r]=i_r_max_arr[i_r-1]+1; @@ -177,6 +178,14 @@ void shear_get_beam_properties(ParamCoLoRe *par) } } //end omp for +#pragma omp single + { + for(i_r=0;i_rnr;i_r++) { + smap->r0[i_r]=1./inv_r_max[i_r]; + smap->rf[i_r]=1./inv_r_max[i_r]; + } + } + free(fac_r_1); free(fac_r_2); free(i_r_max_arr); From 2794d2b0d0581dc2c23f088fe936f11529fb049f Mon Sep 17 00:00:00 2001 From: David Alonso Date: Tue, 21 Jan 2020 09:58:05 +0000 Subject: [PATCH 06/29] made compatible with other faster_shear branch --- examples/cl_test/cl_test.py | 179 ++++++++++-------- examples/read_colore.py | 326 +++++++++++++++++++------------- examples/read_grid.py | 66 ++++--- examples/shear_test/Bz_test.txt | 256 +++++++++++++++++++++++++ examples/shear_test/Nz_test.txt | 256 +++++++++++++++++++++++++ examples/shear_test/param.cfg | 56 ++++++ examples/simple/mk_files.py | 78 ++++---- examples/simple/param.cfg | 2 +- examples/simple/read_skewers.py | 52 ++--- src/common.c | 4 +- src/common.h | 3 +- src/io.c | 10 +- src/srcs.c | 11 +- 13 files changed, 987 insertions(+), 312 deletions(-) create mode 100644 examples/shear_test/Bz_test.txt create mode 100644 examples/shear_test/Nz_test.txt create mode 100644 examples/shear_test/param.cfg diff --git a/examples/cl_test/cl_test.py b/examples/cl_test/cl_test.py index 796ba7e..5b1b899 100644 --- a/examples/cl_test/cl_test.py +++ b/examples/cl_test/cl_test.py @@ -4,90 +4,121 @@ from astropy.io import fits import os.path -nside=128 -npix=hp.nside2npix(nside) +nside = 128 +npix = hp.nside2npix(nside) -#Analyze sources -ifile=0 -nmap=np.zeros(hp.nside2npix(nside)) -e1map=np.zeros(hp.nside2npix(nside)) -e2map=np.zeros(hp.nside2npix(nside)) -while os.path.isfile('examples/cl_test/out_srcs_s1_%d.fits'%ifile) : - hdulist = fits.open('examples/cl_test/out_srcs_s1_%d.fits'%ifile) +# Analyze sources +ifile = 0 +nmap = np.zeros(hp.nside2npix(nside)) +e1map = np.zeros(hp.nside2npix(nside)) +e2map = np.zeros(hp.nside2npix(nside)) +while os.path.isfile('examples/cl_test/out_srcs_s1_%d.fits' % ifile): + hdulist = fits.open('examples/cl_test/out_srcs_s1_%d.fits' % ifile) tbdata = hdulist[1].data - - pix=hp.ang2pix(nside,(90-tbdata['DEC'])*np.pi/180,tbdata['RA']*np.pi/180) - n=np.bincount(pix,minlength=npix) - e1=np.bincount(pix,minlength=npix,weights=tbdata['E1']) - e2=np.bincount(pix,minlength=npix,weights=tbdata['E2']) - nmap+=n - e1map+=e1 - e2map+=e2 - ifile+=1 -ndens=(np.sum(nmap)+0.0)/(4*np.pi) -mp_e1=e1map/nmap; mp_e1[nmap<=0]=0 -mp_e2=e2map/nmap; mp_e2[nmap<=0]=0 -mp_d=(nmap+0.0)/np.mean(nmap+0.0)-1 -mp_db,mp_E,mp_B=hp.alm2map(hp.map2alm(np.array([mp_d,mp_e1,mp_e2]),pol=True),pol=False,nside=nside) + pix = hp.ang2pix(nside, + np.radians(90-tbdata['DEC']), + np.radians(tbdata['RA'])) + n = np.bincount(pix, minlength=npix) + e1 = np.bincount(pix, minlength=npix, weights=tbdata['E1']) + e2 = np.bincount(pix, minlength=npix, weights=tbdata['E2']) + nmap += n + e1map += e1 + e2map += e2 + ifile += 1 -lt,cls_dd=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_dd.txt',unpack=True); -lt,clt_dl=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_d1l2.txt',unpack=True); -lt,clt_ll=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_ll.txt',unpack=True); -lt,clt_kd=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_dc.txt',unpack=True); -lt,clt_kk=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_cc.txt',unpack=True); -lt,clt_id=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_di.txt',unpack=True); -lt,clt_ii=np.loadtxt('examples/cl_test/pred_lj/outlj_cl_ii.txt',unpack=True); -cln_dd=np.ones_like(lt)/ndens -clt_dd=cls_dd+cln_dd -cld_dd,cld_ee,cld_bb,cld_de,cld_eb,cld_db=hp.anafast(np.array([mp_d,mp_e1,mp_e2]),pol=True); -ld=np.arange(len(cld_dd)); +ndens = (np.sum(nmap)+0.0)/(4*np.pi) +mp_e1 = e1map/nmap +mp_e1[nmap <= 0] = 0 +mp_e2 = e2map / nmap +mp_e2[nmap <= 0] = 0 +mp_d = (nmap + 0.0) / np.mean(nmap + 0.0) - 1 +mp_db, mp_E, mp_B = hp.alm2map(hp.map2alm(np.array([mp_d, mp_e1, mp_e2]), + pol=True), + pol=False, + nside=nside) -#Analyze kappa -mp_k=hp.read_map("examples/cl_test/out_kappa_z000.fits") -cld_kk=hp.anafast(mp_k); ld=np.arange(len(cld_kk)) -cld_kd=hp.anafast(mp_k,map2=mp_d) +lt, cls_dd = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_dd.txt', + unpack=True) +lt, clt_dl = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_d1l2.txt', + unpack=True) +lt, clt_ll = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_ll.txt', + unpack=True) +lt, clt_kd = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_dc.txt', + unpack=True) +lt, clt_kk = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_cc.txt', + unpack=True) +lt, clt_id = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_di.txt', + unpack=True) +lt, clt_ii = np.loadtxt('examples/cl_test/pred_lj/outlj_cl_ii.txt', + unpack=True) +cln_dd = np.ones_like(lt) / ndens +clt_dd = cls_dd + cln_dd +d = hp.anafast(np.array([mp_d, mp_e1, mp_e2]), pol=True) +cld_dd, cld_ee, cld_bb, cld_de, cld_eb, cld_db = d +ld = np.arange(len(cld_dd)) -#Analyze ISW -mp_i=hp.read_map("examples/cl_test/out_isw_z000.fits") -cld_ii=hp.anafast(mp_i); ld=np.arange(len(cld_ii)) -cld_id=hp.anafast(mp_i,map2=mp_d) +# Analyze kappa +mp_k = hp.read_map("examples/cl_test/out_kappa_z000.fits") +cld_kk = hp.anafast(mp_k) +ld = np.arange(len(cld_kk)) +cld_kd = hp.anafast(mp_k, map2=mp_d) -#Plots -hp.mollview(mp_d,title='$\\delta_g$'); -hp.mollview(mp_E,title='$\\gamma^E_g$'); -hp.mollview(mp_B,title='$\\gamma^B_g$'); -hp.mollview(mp_e1,title='$e_1$'); -hp.mollview(mp_e2,title='$e_2$'); -hp.mollview(mp_k,title='$\\kappa$'); -hp.mollview(mp_i,title='$\\dot{\\phi}$'); +# Analyze ISW +mp_i = hp.read_map("examples/cl_test/out_isw_z000.fits") +cld_ii = hp.anafast(mp_i) +ld = np.arange(len(cld_ii)) +cld_id = hp.anafast(mp_i, map2=mp_d) + +# Plots +hp.mollview(mp_d, title='$\\delta_g$') +hp.mollview(mp_E, title='$\\gamma^E_g$') +hp.mollview(mp_B, title='$\\gamma^B_g$') +hp.mollview(mp_e1, title='$e_1$') +hp.mollview(mp_e2, title='$e_2$') +hp.mollview(mp_k, title='$\\kappa$') +hp.mollview(mp_i, title='$\\dot{\\phi}$') plt.figure() -plt.hist(mp_e1,bins=100,histtype='step'); -plt.hist(mp_e2,bins=100,histtype='step'); -plt.xlabel('$e_1,\\,e_2$',fontsize=16) +plt.hist(mp_e1, bins=100, histtype='step') +plt.hist(mp_e2, bins=100, histtype='step') +plt.xlabel('$e_1, \\, e_2$', fontsize=16) plt.figure() -plt.plot(ld,cld_dd,'r-',label='$\\delta_g\\times\\delta_g$',lw=2) -plt.plot(ld,cld_ee,'y-',label='$\\gamma^E_g\\times\\gamma^E_g$',lw=2) -plt.plot(ld,cld_de,'c-',label='$\\gamma^E_g\\times\\delta_g$',lw=2) -plt.plot(ld,cld_bb,'m',linestyle='solid' ,label='$\\gamma^B_g\\times\\gamma^B_g$',lw=2) -plt.plot(ld,cld_eb,'m',linestyle='dashed',label='$\\gamma^E_g\\times\\gamma^B_g$',lw=2) -plt.plot(ld,cld_db,'m',linestyle='dotted',label='$\\gamma^B_g\\times\\delta_g$',lw=2) -plt.plot(ld,cld_kd,'g-',label='$\\kappa-\\delta_g$',lw=2) -plt.plot(ld,cld_kk,'b-',label='$\\kappa-\\kappa$',lw=2) -plt.plot(ld,cld_id,'r-',label='$\\dot{\\phi}-\\delta_g$',lw=2) -plt.plot(ld,cld_ii,'y-',label='$\\dot{\\phi}-\\dot{\\phi}$',lw=2) -plt.plot(lt,clt_dd,'k-',lw=2) -plt.plot(lt,2*clt_dl,'k-',lw=2) -plt.plot(lt,4*clt_ll,'k-',lw=2) -plt.plot(lt,clt_kd,'k-',lw=2) -plt.plot(lt,clt_kk,'k-',lw=2) -plt.plot(lt,clt_id,'k--',lw=2) -plt.plot(lt,clt_ii,'k--',lw=2) +plt.plot(ld, cld_dd, 'r-', + label='$\\delta_g\\times\\delta_g$', lw=2) +plt.plot(ld, cld_ee, 'y-', + label='$\\gamma^E_g\\times\\gamma^E_g$', lw=2) +plt.plot(ld, cld_de, 'c-', + label='$\\gamma^E_g\\times\\delta_g$', lw=2) +plt.plot(ld, cld_bb, 'm', + linestyle='solid', + label='$\\gamma^B_g\\times\\gamma^B_g$', lw=2) +plt.plot(ld, cld_eb, 'm', + linestyle='dashed', + label='$\\gamma^E_g\\times\\gamma^B_g$', lw=2) +plt.plot(ld, cld_db, 'm', + linestyle='dotted', + label='$\\gamma^B_g\\times\\delta_g$', lw=2) +plt.plot(ld, cld_kd, 'g-', + label='$\\kappa-\\delta_g$', lw=2) +plt.plot(ld, cld_kk, 'b-', + label='$\\kappa-\\kappa$', lw=2) +plt.plot(ld, cld_id, 'r-', + label='$\\dot{\\phi}-\\delta_g$', lw=2) +plt.plot(ld, cld_ii, 'y-', + label='$\\dot{\\phi}-\\dot{\\phi}$', lw=2) +plt.plot(lt, clt_dd, 'k-', lw=2) +plt.plot(lt, 2*clt_dl, 'k-', lw=2) +plt.plot(lt, 4*clt_ll, 'k-', lw=2) +plt.plot(lt, clt_kd, 'k-', lw=2) +plt.plot(lt, clt_kk, 'k-', lw=2) +plt.plot(lt, clt_id, 'k--', lw=2) +plt.plot(lt, clt_ii, 'k--', lw=2) plt.loglog() -plt.xlabel('$\\ell$',fontsize=16) -plt.ylabel('$C_\\ell$',fontsize=16) -plt.xlim([2,192]) -plt.legend(loc='lower left',frameon=False,labelspacing=0.1,ncol=2) +plt.xlabel('$\\ell$', fontsize=16) +plt.ylabel('$C_\\ell$', fontsize=16) +plt.xlim([2, 192]) +plt.legend(loc='lower left', frameon=False, + labelspacing=0.1, ncol=2) plt.show() diff --git a/examples/read_colore.py b/examples/read_colore.py index d9c7e37..8e933f4 100644 --- a/examples/read_colore.py +++ b/examples/read_colore.py @@ -6,141 +6,195 @@ import sys import os -def read_ascii(fname) : - ifile=0 - - type_arr,ra_arr,dec_arr,z0_arr,rsd_arr,e1_arr,e2_arr=np.loadtxt(fname+"_%d.txt"%ifile,unpack=True) - - ifile=1 - while os.path.isfile(fname+"_%d.h5"%ifile) : - typ,ra,dec,z0,rsd,e1,e2=np.loadtxt(fname+"_%d.txt"%ifile,unpack=True) - type_arr=np.concatenate((type_arr,typ_arr)) - ra_arr=np.concatenate((ra_arr,ra_arr)) - dec_arr=np.concatenate((dec_arr,dec_arr)) - z0_arr=np.concatenate((z0_arr,z0_arr)) - rsd_arr=np.concatenate((rsd_arr,rsd_arr)) - e1_arr=np.concatenate((e1_arr,e1_arr)) - e2_arr=np.concatenate((e2_arr,e2_arr)) - ifile+=1 - - return ra_arr,dec_arr,z0_arr,rsd_arr,type_arr,e1_arr,e2_arr,None,None,None - -def read_hdf5(fname,ipop) : - ifile=0 - - ff=h5.File(fname+"_%d.h5"%ifile,"r") - tabname='/sources%d'%ipop - - ra_arr=ff[tabname]['RA'] - dec_arr=ff[tabname]['DEC'] - z0_arr=ff[tabname]['Z_COSMO'] - rsd_arr=ff[tabname]['DZ_RSD'] - e1_arr=ff[tabname]['E1'] - e2_arr=ff[tabname]['E2'] - - ifile=1 - while os.path.isfile(fname+"_%d.h5"%ifile) : - ff=h5.File(fname+"_%d.h5"%ifile,"r") - tabname='/sources%d'%ipop - - ra_arr=np.concatenate((ra_arr,ff[tabname]['RA'])) - dec_arr=np.concatenate((dec_arr,ff[tabname]['DEC'])) - z0_arr=np.concatenate((z0_arr,ff[tabname]['Z_COSMO'])) - rsd_arr=np.concatenate((rsd_arr,ff[tabname]['DZ_RSD'])) - e1_arr=np.concatenate((e1_arr,ff[tabname]['E1'])) - e2_arr=np.concatenate((e2_arr,ff[tabname]['E2'])) - ifile+=1 - - type_arr=ipop*np.ones(len(ra_arr)) - - return ra_arr,dec_arr,z0_arr,rsd_arr,type_arr,e1_arr,e2_arr,None,None,None - -def read_fits(fname) : - ifile=0 - - fi=fits.open(fname+"_%d.fits"%ifile) - data=(fi[1]).data - if len(fi)>2 : with_skewers=True - else : with_skewers=False - - if with_skewers : - dens_skewers=(fi[2]).data - vrad_skewers=(fi[3]).data - data_skewers=(fi[4]).data - else : - dens_skewers=None - vrad_skewers=None - data_skewers=None - - ra_arr=data['RA'] - dec_arr=data['DEC'] - z0_arr=data['Z_COSMO'] - rsd_arr=data['DZ_RSD'] - e1_arr=data['E1'] - e2_arr=data['E2'] - type_arr=data['TYPE'] - - ifile=1 - while os.path.isfile(fname+"_%d.fits"%ifile) : - fi=fits.open(fname+"_%d.fits"%ifile) - data=(fi[1]).data - if with_skewers : - dens_skw=(fi[2]).data - vrad_skw=(fi[3]).data - - ra_arr=np.concatenate((ra_arr,data['RA'])) - dec_arr=np.concatenate((dec_arr,data['DEC'])) - z0_arr=np.concatenate((z0_arr,data['Z_COSMO'])) - rsd_arr=np.concatenate((rsd_arr,data['DZ_RSD'])) - e1_arr=np.concatenate((e1_arr,data['E1'])) - e2_arr=np.concatenate((e2_arr,data['E2'])) - type_arr=np.concatenate((type_arr,data['TYPE'])) - if with_skewers : - dens_skewers=np.concatenate((dens_skewers,dens_skw)) - vrad_skewers=np.concatenate((vrad_skewers,vrad_skw)) - ifile+=1 - - return ra_arr,dec_arr,z0_arr,rsd_arr,type_arr,e1_arr,e2_arr,dens_skewers,vrad_skewers,data_skewers - -if len(sys.argv)!= 3 : - print("Usage: read_colore.py file_name file_format ('ASCII', 'FITS' or 'HDF5')") + +def read_ascii(fname): + ifile = 0 + + d = np.loadtxt(fname + "_%d.txt" % ifile, unpack=True) + has_shear = False + if len(d) == 7: + has_shear = True + type_arr, ra_arr, dec_arr, z0_arr, rsd_arr, e1_arr, e2_arr = d + else: + type_arr, ra_arr, dec_arr, z0_arr, rsd_arr = d + e1_arr = None + e2_arr = None + + ifile = 1 + while os.path.isfile(fname + "_%d.h5" % ifile): + d = np.loadtxt(fname + "_%d.txt" % ifile, unpack=True) + if has_shear: + typ, ra, dec, z0, rsd, e1, e2 = d + else: + typ, ra, dec, z0, rsd = d + type_arr = np.concatenate((type_arr, type_arr)) + ra_arr = np.concatenate((ra_arr, ra_arr)) + dec_arr = np.concatenate((dec_arr, dec_arr)) + z0_arr = np.concatenate((z0_arr, z0_arr)) + rsd_arr = np.concatenate((rsd_arr, rsd_arr)) + if has_shear: + e1_arr = np.concatenate((e1_arr, e1_arr)) + e2_arr = np.concatenate((e2_arr, e2_arr)) + ifile += 1 + + return ra_arr, dec_arr, z0_arr, rsd_arr, \ + type_arr, e1_arr, e2_arr, None, None, None + + +def read_hdf5(fname, ipop): + ifile = 0 + + ff = h5.File(fname + "_%d.h5" % ifile, "r") + tabname = '/sources%d' % ipop + + ra_arr = ff[tabname]['RA'] + dec_arr = ff[tabname]['DEC'] + z0_arr = ff[tabname]['Z_COSMO'] + rsd_arr = ff[tabname]['DZ_RSD'] + + has_shear = 'E1' in ff[tabname].keys() + if has_shear: + e1_arr = ff[tabname]['E1'] + e2_arr = ff[tabname]['E2'] + else: + e1_arr = None + e2_arr = None + + ifile = 1 + while os.path.isfile(fname + "_%d.h5" % ifile): + ff = h5.File(fname + "_%d.h5" % ifile, "r") + tabname = '/sources%d' % ipop + + ra_arr = np.concatenate((ra_arr, ff[tabname]['RA'])) + dec_arr = np.concatenate((dec_arr, ff[tabname]['DEC'])) + z0_arr = np.concatenate((z0_arr, ff[tabname]['Z_COSMO'])) + rsd_arr = np.concatenate((rsd_arr, ff[tabname]['DZ_RSD'])) + if has_shear: + e1_arr = np.concatenate((e1_arr, ff[tabname]['E1'])) + e2_arr = np.concatenate((e2_arr, ff[tabname]['E2'])) + ifile += 1 + + type_arr = ipop*np.ones(len(ra_arr)) + + return ra_arr, dec_arr, z0_arr, rsd_arr, \ + type_arr, e1_arr, e2_arr, None, None, None + + +def read_fits(fname): + ifile = 0 + + fi = fits.open(fname + "_%d.fits" % ifile) + data = (fi[1]).data + with_skewers = len(fi) > 2 + + if with_skewers: + dens_skewers = (fi[2]).data + vrad_skewers = (fi[3]).data + data_skewers = (fi[4]).data + else: + dens_skewers = None + vrad_skewers = None + data_skewers = None + + ra_arr = data['RA'] + dec_arr = data['DEC'] + z0_arr = data['Z_COSMO'] + rsd_arr = data['DZ_RSD'] + has_shear = 'E1' in data.names + if has_shear: + e1_arr = data['E1'] + e2_arr = data['E2'] + else: + e1_arr = None + e2_arr = None + type_arr = data['TYPE'] + + ifile = 1 + while os.path.isfile(fname + "_%d.fits" % ifile): + fi = fits.open(fname + "_%d.fits" % ifile) + data = (fi[1]).data + if with_skewers: + dens_skw = (fi[2]).data + vrad_skw = (fi[3]).data + + ra_arr = np.concatenate((ra_arr, data['RA'])) + dec_arr = np.concatenate((dec_arr, data['DEC'])) + z0_arr = np.concatenate((z0_arr, data['Z_COSMO'])) + rsd_arr = np.concatenate((rsd_arr, data['DZ_RSD'])) + if has_shear: + e1_arr = np.concatenate((e1_arr, data['E1'])) + e2_arr = np.concatenate((e2_arr, data['E2'])) + type_arr = np.concatenate((type_arr, data['TYPE'])) + if with_skewers: + dens_skewers = np.concatenate((dens_skewers, dens_skw)) + vrad_skewers = np.concatenate((vrad_skewers, vrad_skw)) + ifile += 1 + + return ra_arr, dec_arr, z0_arr, rsd_arr, type_arr, \ + e1_arr, e2_arr, dens_skewers, vrad_skewers, data_skewers + + +if len(sys.argv) != 3: + print("Usage: read_colore.py file_name file_format " + "('ASCII', 'FITS' or 'HDF5')") exit(1) -fname=sys.argv[1] -fmt=sys.argv[2] - -if fmt=='ASCII' : - ra_arr,dec_arr,z_arr,rsd_arr,type_arr,e1_arr,e2_arr,dskw_arr,vskw_arr,data_skw=read_ascii(fname) -elif fmt=='FITS' : - ra_arr,dec_arr,z_arr,rsd_arr,type_arr,e1_arr,e2_arr,dskw_arr,vskw_arr,data_skw=read_fits(fname) -elif fmt=='HDF5' : - ra_arr,dec_arr,z_arr,rsd_arr,type_arr,e1_arr,e2_arr,dskw_arr,vskw_arr,data_skw=read_hdf5(fname,1) - -nside=64 -npix=hp.nside2npix(nside) -mp=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix])[0] -mpr=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix],weights=rsd_arr)[0] -mpe1=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix],weights=e1_arr)[0] -mpe2=np.histogram(hp.ang2pix(nside,np.pi*(90-dec_arr)/180,np.pi*ra_arr/180),bins=npix,range=[0,npix],weights=e2_arr)[0] -plt.figure(); -plt.hist(z_arr,bins=100,histtype='step'); -plt.xlabel('$z$',fontsize=16); plt.ylabel('$N(z)$',fontsize=16); -hp.mollview(mp,title='$N_g(\\hat{\\bf n})$') -hp.mollview(mpr/mp,title='$v_r$') -hp.mollview(mpe1/mp,title='$e_1$') -hp.mollview(mpe2/mp,title='$e_2$') - -if dskw_arr is not None : - #Plot skewers - nside_skw=64 - ipix_skw=1000 - id_in_pix=np.where(hp.ang2pix(nside_skw,(90-dec_arr)*np.pi/180,ra_arr*np.pi/180)==ipix_skw)[0] - cols=['r','g','b','y','k'] - plt.figure(); plt.title('Skewers'); plt.xlabel('$z$'); plt.ylabel('$\\delta$'); - for i in id_in_pix : - plt.plot(data_skw['Z'],dskw_arr[i,:],cols[i%5]+'-') - plt.plot([z_arr[i],z_arr[i]],[1.1*np.amin(dskw_arr[i,:]),1.1*np.amax(dskw_arr[i,:])],cols[i%5]+'--') - plt.figure(); plt.title('Skewers'); plt.xlabel('$z$'); plt.ylabel('$v_r$'); - for i in id_in_pix : - plt.plot(data_skw['Z'],vskw_arr[i,:],cols[i%5]+'-') - plt.plot([z_arr[i],z_arr[i]],[1.1*np.amin(vskw_arr[i,:]),1.1*np.amax(vskw_arr[i,:])],cols[i%5]+'--') +fname = sys.argv[1] +fmt = sys.argv[2] + +if fmt == 'ASCII': + ra_arr, dec_arr, z_arr, rsd_arr, type_arr, \ + e1_arr, e2_arr, dskw_arr, vskw_arr, data_skw = read_ascii(fname) +elif fmt == 'FITS': + ra_arr, dec_arr, z_arr, rsd_arr, type_arr, \ + e1_arr, e2_arr, dskw_arr, vskw_arr, data_skw = read_fits(fname) +elif fmt == 'HDF5': + ra_arr, dec_arr, z_arr, rsd_arr, type_arr, \ + e1_arr, e2_arr, dskw_arr, vskw_arr, data_skw = read_hdf5(fname, 1) + +nside = 64 +npix = hp.nside2npix(nside) +ipixs = hp.ang2pix(nside, np.radians(90-dec_arr), np.radians(ra_arr)) +mp = np.histogram(ipixs, bins=npix, range=[0, npix])[0] +hp.mollview(mp, title='$N_g(\\hat{\\bf n})$') +mpr = np.histogram(ipixs, bins=npix, range=[0, npix], weights=rsd_arr)[0] +hp.mollview(mpr/mp, title='$v_r$') +if e1_arr is not None: + mpe1 = np.histogram(ipixs, bins=npix, range=[0, npix], weights=e1_arr)[0] + mpe2 = np.histogram(ipixs, bins=npix, range=[0, npix], weights=e2_arr)[0] + hp.mollview(mpe1/mp, title='$e_1$') + hp.mollview(mpe2/mp, title='$e_2$') +plt.figure() +plt.hist(z_arr, bins=100, histtype='step') +plt.xlabel('$z$', fontsize=16) +plt.ylabel('$N(z)$', fontsize=16) + +if dskw_arr is not None: + # Plot skewers + nside_skw = 64 + ipix_skw = 1000 + id_in_pix = np.where(hp.ang2pix(nside_skw, + np.radians(90 - dec_arr), + np.radians(ra_arr)) + == ipix_skw)[0] + cols = ['r', 'g', 'b', 'y', 'k'] + plt.figure() + plt.title('Skewers') + plt.xlabel('$z$') + plt.ylabel('$\\delta$') + for i in id_in_pix: + plt.plot(data_skw['Z'], dskw_arr[i, :], cols[i % 5]+'-') + plt.plot([z_arr[i], z_arr[i]], + [1.1 * np.amin(dskw_arr[i, :]), + 1.1 * np.amax(dskw_arr[i, :])], + cols[i % 5]+'--') + plt.figure() + plt.title('Skewers') + plt.xlabel('$z$') + plt.ylabel('$v_r$') + for i in id_in_pix: + plt.plot(data_skw['Z'], vskw_arr[i, :], cols[i % 5]+'-') + plt.plot([z_arr[i], z_arr[i]], + [1.1 * np.amin(vskw_arr[i, :]), + 1.1 * np.amax(vskw_arr[i, :])], + cols[i % 5]+'--') plt.show() diff --git a/examples/read_grid.py b/examples/read_grid.py index 38f65f7..407519f 100644 --- a/examples/read_grid.py +++ b/examples/read_grid.py @@ -1,45 +1,51 @@ import numpy as np import sys as sys import matplotlib.pyplot as plt -from matplotlib import cm -def read_grid(prefix) : - f=open(prefix+"_0.dat","rb") - nfiles,size_float=np.fromfile(f,dtype=np.int32,count=2) - lbox=np.fromfile(f,dtype=np.float64,count=1)[0] - ngrid=np.fromfile(f,dtype=np.int32,count=1)[0] + +def read_grid(prefix): + f = open(prefix+"_0.dat", "rb") + nfiles, size_float = np.fromfile(f, dtype=np.int32, count=2) + lbox = np.fromfile(f, dtype=np.float64, count=1)[0] + ngrid = np.fromfile(f, dtype=np.int32, count=1)[0] f.close() - if size_float==4 : - f_type=np.float32 - else : - f_type=np.float64 - - grid_out=np.zeros([ngrid,ngrid,ngrid]) - for ifil in np.arange(nfiles) : - f=open(prefix+"_%d.dat"%ifil,"rb") - nf,sz=np.fromfile(f,dtype=np.int32,count=2) - lb=np.fromfile(f,dtype=np.float64,count=1) - ng,nz_here,iz0_here=np.fromfile(f,dtype=np.int32,count=3) - for iz in np.arange(nz_here) : - grid_out[iz0_here+iz,:,:]=np.fromfile(f,dtype=f_type,count=ng*ng).reshape([ng,ng]) + if size_float == 4: + f_type = np.float32 + else: + f_type = np.float64 + + grid_out = np.zeros([ngrid, ngrid, ngrid]) + for ifil in np.arange(nfiles): + f = open(prefix+"_%d.dat" % ifil, "rb") + nf, sz = np.fromfile(f, dtype=np.int32, count=2) + _ = np.fromfile(f, dtype=np.float64, count=1) + ng, nz_here, iz0_here = np.fromfile(f, dtype=np.int32, count=3) + for iz in np.arange(nz_here): + d = np.fromfile(f, dtype=f_type, count=ng*ng).reshape([ng, ng]) + grid_out[iz0_here+iz, :, :] = d f.close() - return ngrid,lbox,np.array(grid_out) + return ngrid, lbox, np.array(grid_out) -if len(sys.argv)!= 2 : + +if len(sys.argv) != 2: print("Usage: read_grid.py prefix") exit(1) -prefix=sys.argv[1] +prefix = sys.argv[1] + +ng, lb, dens = read_grid(prefix) +print("Ngrid = %d" % ng) +print("Lbox = %.3lf Mpc/h" % lb) -ng,lb,dens=read_grid(prefix) -print("Ngrid=%d"%ng) -print("Lbox=%.3lf Mpc/h"%lb) -def plot_slice(slic) : + +def plot_slice(slic): plt.figure() - plt.imshow(slic,origin='lower',interpolation='nearest'); + plt.imshow(slic, origin='lower', interpolation='nearest') plt.colorbar() -plot_slice(dens[ng//2,:,:]) -plot_slice(dens[:,ng//2,:]) -plot_slice(dens[:,:,ng//2]) + + +plot_slice(dens[ng//2, :, :]) +plot_slice(dens[:, ng//2, :]) +plot_slice(dens[:, :, ng//2]) plt.show() diff --git a/examples/shear_test/Bz_test.txt b/examples/shear_test/Bz_test.txt new file mode 100644 index 0000000..f6b048d --- /dev/null +++ b/examples/shear_test/Bz_test.txt @@ -0,0 +1,256 @@ +0.000000000000000000e+00 1.000000000000000000e+00 +3.921568627450980338e-03 1.000000000000000000e+00 +7.843137254901960675e-03 1.000000000000000000e+00 +1.176470588235294101e-02 1.000000000000000000e+00 +1.568627450980392135e-02 1.000000000000000000e+00 +1.960784313725490169e-02 1.000000000000000000e+00 +2.352941176470588203e-02 1.000000000000000000e+00 +2.745098039215686236e-02 1.000000000000000000e+00 +3.137254901960784270e-02 1.000000000000000000e+00 +3.529411764705882304e-02 1.000000000000000000e+00 +3.921568627450980338e-02 1.000000000000000000e+00 +4.313725490196078372e-02 1.000000000000000000e+00 +4.705882352941176405e-02 1.000000000000000000e+00 +5.098039215686274439e-02 1.000000000000000000e+00 +5.490196078431372473e-02 1.000000000000000000e+00 +5.882352941176470507e-02 1.000000000000000000e+00 +6.274509803921568540e-02 1.000000000000000000e+00 +6.666666666666666574e-02 1.000000000000000000e+00 +7.058823529411764608e-02 1.000000000000000000e+00 +7.450980392156862642e-02 1.000000000000000000e+00 +7.843137254901960675e-02 1.000000000000000000e+00 +8.235294117647058709e-02 1.000000000000000000e+00 +8.627450980392156743e-02 1.000000000000000000e+00 +9.019607843137254777e-02 1.000000000000000000e+00 +9.411764705882352811e-02 1.000000000000000000e+00 +9.803921568627450844e-02 1.000000000000000000e+00 +1.019607843137254888e-01 1.000000000000000000e+00 +1.058823529411764691e-01 1.000000000000000000e+00 +1.098039215686274495e-01 1.000000000000000000e+00 +1.137254901960784298e-01 1.000000000000000000e+00 +1.176470588235294101e-01 1.000000000000000000e+00 +1.215686274509803905e-01 1.000000000000000000e+00 +1.254901960784313708e-01 1.000000000000000000e+00 +1.294117647058823650e-01 1.000000000000000000e+00 +1.333333333333333315e-01 1.000000000000000000e+00 +1.372549019607843257e-01 1.000000000000000000e+00 +1.411764705882352922e-01 1.000000000000000000e+00 +1.450980392156862864e-01 1.000000000000000000e+00 +1.490196078431372528e-01 1.000000000000000000e+00 +1.529411764705882470e-01 1.000000000000000000e+00 +1.568627450980392135e-01 1.000000000000000000e+00 +1.607843137254902077e-01 1.000000000000000000e+00 +1.647058823529411742e-01 1.000000000000000000e+00 +1.686274509803921684e-01 1.000000000000000000e+00 +1.725490196078431349e-01 1.000000000000000000e+00 +1.764705882352941291e-01 1.000000000000000000e+00 +1.803921568627450955e-01 1.000000000000000000e+00 +1.843137254901960898e-01 1.000000000000000000e+00 +1.882352941176470562e-01 1.000000000000000000e+00 +1.921568627450980504e-01 1.000000000000000000e+00 +1.960784313725490169e-01 1.000000000000000000e+00 +2.000000000000000111e-01 1.000000000000000000e+00 +2.039215686274509776e-01 1.000000000000000000e+00 +2.078431372549019718e-01 1.000000000000000000e+00 +2.117647058823529382e-01 1.000000000000000000e+00 +2.156862745098039325e-01 1.000000000000000000e+00 +2.196078431372548989e-01 1.000000000000000000e+00 +2.235294117647058931e-01 1.000000000000000000e+00 +2.274509803921568596e-01 1.000000000000000000e+00 +2.313725490196078538e-01 1.000000000000000000e+00 +2.352941176470588203e-01 1.000000000000000000e+00 +2.392156862745098145e-01 1.000000000000000000e+00 +2.431372549019607809e-01 1.000000000000000000e+00 +2.470588235294117752e-01 1.000000000000000000e+00 +2.509803921568627416e-01 1.000000000000000000e+00 +2.549019607843137081e-01 1.000000000000000000e+00 +2.588235294117647300e-01 1.000000000000000000e+00 +2.627450980392156965e-01 1.000000000000000000e+00 +2.666666666666666630e-01 1.000000000000000000e+00 +2.705882352941176294e-01 1.000000000000000000e+00 +2.745098039215686514e-01 1.000000000000000000e+00 +2.784313725490196179e-01 1.000000000000000000e+00 +2.823529411764705843e-01 1.000000000000000000e+00 +2.862745098039215508e-01 1.000000000000000000e+00 +2.901960784313725727e-01 1.000000000000000000e+00 +2.941176470588235392e-01 1.000000000000000000e+00 +2.980392156862745057e-01 1.000000000000000000e+00 +3.019607843137254721e-01 1.000000000000000000e+00 +3.058823529411764941e-01 1.000000000000000000e+00 +3.098039215686274606e-01 1.000000000000000000e+00 +3.137254901960784270e-01 1.000000000000000000e+00 +3.176470588235293935e-01 1.000000000000000000e+00 +3.215686274509804154e-01 1.000000000000000000e+00 +3.254901960784313819e-01 1.000000000000000000e+00 +3.294117647058823484e-01 1.000000000000000000e+00 +3.333333333333333148e-01 1.000000000000000000e+00 +3.372549019607843368e-01 1.000000000000000000e+00 +3.411764705882353033e-01 1.000000000000000000e+00 +3.450980392156862697e-01 1.000000000000000000e+00 +3.490196078431372362e-01 1.000000000000000000e+00 +3.529411764705882582e-01 1.000000000000000000e+00 +3.568627450980392246e-01 1.000000000000000000e+00 +3.607843137254901911e-01 1.000000000000000000e+00 +3.647058823529411575e-01 1.000000000000000000e+00 +3.686274509803921795e-01 1.000000000000000000e+00 +3.725490196078431460e-01 1.000000000000000000e+00 +3.764705882352941124e-01 1.000000000000000000e+00 +3.803921568627450789e-01 1.000000000000000000e+00 +3.843137254901961009e-01 1.000000000000000000e+00 +3.882352941176470673e-01 1.000000000000000000e+00 +3.921568627450980338e-01 1.000000000000000000e+00 +3.960784313725490002e-01 1.000000000000000000e+00 +4.000000000000000222e-01 1.000000000000000000e+00 +4.039215686274509887e-01 1.000000000000000000e+00 +4.078431372549019551e-01 1.000000000000000000e+00 +4.117647058823529216e-01 1.000000000000000000e+00 +4.156862745098039436e-01 1.000000000000000000e+00 +4.196078431372549100e-01 1.000000000000000000e+00 +4.235294117647058765e-01 1.000000000000000000e+00 +4.274509803921568429e-01 1.000000000000000000e+00 +4.313725490196078649e-01 1.000000000000000000e+00 +4.352941176470588314e-01 1.000000000000000000e+00 +4.392156862745097978e-01 1.000000000000000000e+00 +4.431372549019607643e-01 1.000000000000000000e+00 +4.470588235294117863e-01 1.000000000000000000e+00 +4.509803921568627527e-01 1.000000000000000000e+00 +4.549019607843137192e-01 1.000000000000000000e+00 +4.588235294117646856e-01 1.000000000000000000e+00 +4.627450980392157076e-01 1.000000000000000000e+00 +4.666666666666666741e-01 1.000000000000000000e+00 +4.705882352941176405e-01 1.000000000000000000e+00 +4.745098039215686070e-01 1.000000000000000000e+00 +4.784313725490196290e-01 1.000000000000000000e+00 +4.823529411764705954e-01 1.000000000000000000e+00 +4.862745098039215619e-01 1.000000000000000000e+00 +4.901960784313725283e-01 1.000000000000000000e+00 +4.941176470588235503e-01 1.000000000000000000e+00 +4.980392156862745168e-01 1.000000000000000000e+00 +5.019607843137254832e-01 1.000000000000000000e+00 +5.058823529411764497e-01 1.000000000000000000e+00 +5.098039215686274161e-01 1.000000000000000000e+00 +5.137254901960783826e-01 1.000000000000000000e+00 +5.176470588235294601e-01 1.000000000000000000e+00 +5.215686274509804266e-01 1.000000000000000000e+00 +5.254901960784313930e-01 1.000000000000000000e+00 +5.294117647058823595e-01 1.000000000000000000e+00 +5.333333333333333259e-01 1.000000000000000000e+00 +5.372549019607842924e-01 1.000000000000000000e+00 +5.411764705882352589e-01 1.000000000000000000e+00 +5.450980392156862253e-01 1.000000000000000000e+00 +5.490196078431373028e-01 1.000000000000000000e+00 +5.529411764705882693e-01 1.000000000000000000e+00 +5.568627450980392357e-01 1.000000000000000000e+00 +5.607843137254902022e-01 1.000000000000000000e+00 +5.647058823529411686e-01 1.000000000000000000e+00 +5.686274509803921351e-01 1.000000000000000000e+00 +5.725490196078431016e-01 1.000000000000000000e+00 +5.764705882352940680e-01 1.000000000000000000e+00 +5.803921568627451455e-01 1.000000000000000000e+00 +5.843137254901961120e-01 1.000000000000000000e+00 +5.882352941176470784e-01 1.000000000000000000e+00 +5.921568627450980449e-01 1.000000000000000000e+00 +5.960784313725490113e-01 1.000000000000000000e+00 +5.999999999999999778e-01 1.000000000000000000e+00 +6.039215686274509443e-01 1.000000000000000000e+00 +6.078431372549019107e-01 1.000000000000000000e+00 +6.117647058823529882e-01 1.000000000000000000e+00 +6.156862745098039547e-01 1.000000000000000000e+00 +6.196078431372549211e-01 1.000000000000000000e+00 +6.235294117647058876e-01 1.000000000000000000e+00 +6.274509803921568540e-01 1.000000000000000000e+00 +6.313725490196078205e-01 1.000000000000000000e+00 +6.352941176470587870e-01 1.000000000000000000e+00 +6.392156862745097534e-01 1.000000000000000000e+00 +6.431372549019608309e-01 1.000000000000000000e+00 +6.470588235294117974e-01 1.000000000000000000e+00 +6.509803921568627638e-01 1.000000000000000000e+00 +6.549019607843137303e-01 1.000000000000000000e+00 +6.588235294117646967e-01 1.000000000000000000e+00 +6.627450980392156632e-01 1.000000000000000000e+00 +6.666666666666666297e-01 1.000000000000000000e+00 +6.705882352941175961e-01 1.000000000000000000e+00 +6.745098039215686736e-01 1.000000000000000000e+00 +6.784313725490196401e-01 1.000000000000000000e+00 +6.823529411764706065e-01 1.000000000000000000e+00 +6.862745098039215730e-01 1.000000000000000000e+00 +6.901960784313725394e-01 1.000000000000000000e+00 +6.941176470588235059e-01 1.000000000000000000e+00 +6.980392156862744724e-01 1.000000000000000000e+00 +7.019607843137254388e-01 1.000000000000000000e+00 +7.058823529411765163e-01 1.000000000000000000e+00 +7.098039215686274828e-01 1.000000000000000000e+00 +7.137254901960784492e-01 1.000000000000000000e+00 +7.176470588235294157e-01 1.000000000000000000e+00 +7.215686274509803821e-01 1.000000000000000000e+00 +7.254901960784313486e-01 1.000000000000000000e+00 +7.294117647058823151e-01 1.000000000000000000e+00 +7.333333333333332815e-01 1.000000000000000000e+00 +7.372549019607843590e-01 1.000000000000000000e+00 +7.411764705882353255e-01 1.000000000000000000e+00 +7.450980392156862919e-01 1.000000000000000000e+00 +7.490196078431372584e-01 1.000000000000000000e+00 +7.529411764705882248e-01 1.000000000000000000e+00 +7.568627450980391913e-01 1.000000000000000000e+00 +7.607843137254901578e-01 1.000000000000000000e+00 +7.647058823529411242e-01 1.000000000000000000e+00 +7.686274509803922017e-01 1.000000000000000000e+00 +7.725490196078431682e-01 1.000000000000000000e+00 +7.764705882352941346e-01 1.000000000000000000e+00 +7.803921568627451011e-01 1.000000000000000000e+00 +7.843137254901960675e-01 1.000000000000000000e+00 +7.882352941176470340e-01 1.000000000000000000e+00 +7.921568627450980005e-01 1.000000000000000000e+00 +7.960784313725489669e-01 1.000000000000000000e+00 +8.000000000000000444e-01 1.000000000000000000e+00 +8.039215686274510109e-01 1.000000000000000000e+00 +8.078431372549019773e-01 1.000000000000000000e+00 +8.117647058823529438e-01 1.000000000000000000e+00 +8.156862745098039102e-01 1.000000000000000000e+00 +8.196078431372548767e-01 1.000000000000000000e+00 +8.235294117647058432e-01 1.000000000000000000e+00 +8.274509803921568096e-01 1.000000000000000000e+00 +8.313725490196078871e-01 1.000000000000000000e+00 +8.352941176470588536e-01 1.000000000000000000e+00 +8.392156862745098200e-01 1.000000000000000000e+00 +8.431372549019607865e-01 1.000000000000000000e+00 +8.470588235294117530e-01 1.000000000000000000e+00 +8.509803921568627194e-01 1.000000000000000000e+00 +8.549019607843136859e-01 1.000000000000000000e+00 +8.588235294117646523e-01 1.000000000000000000e+00 +8.627450980392157298e-01 1.000000000000000000e+00 +8.666666666666666963e-01 1.000000000000000000e+00 +8.705882352941176627e-01 1.000000000000000000e+00 +8.745098039215686292e-01 1.000000000000000000e+00 +8.784313725490195957e-01 1.000000000000000000e+00 +8.823529411764705621e-01 1.000000000000000000e+00 +8.862745098039215286e-01 1.000000000000000000e+00 +8.901960784313724950e-01 1.000000000000000000e+00 +8.941176470588235725e-01 1.000000000000000000e+00 +8.980392156862745390e-01 1.000000000000000000e+00 +9.019607843137255054e-01 1.000000000000000000e+00 +9.058823529411764719e-01 1.000000000000000000e+00 +9.098039215686274384e-01 1.000000000000000000e+00 +9.137254901960784048e-01 1.000000000000000000e+00 +9.176470588235293713e-01 1.000000000000000000e+00 +9.215686274509803377e-01 1.000000000000000000e+00 +9.254901960784314152e-01 1.000000000000000000e+00 +9.294117647058823817e-01 1.000000000000000000e+00 +9.333333333333333481e-01 1.000000000000000000e+00 +9.372549019607843146e-01 1.000000000000000000e+00 +9.411764705882352811e-01 1.000000000000000000e+00 +9.450980392156862475e-01 1.000000000000000000e+00 +9.490196078431372140e-01 1.000000000000000000e+00 +9.529411764705881804e-01 1.000000000000000000e+00 +9.568627450980392579e-01 1.000000000000000000e+00 +9.607843137254902244e-01 1.000000000000000000e+00 +9.647058823529411908e-01 1.000000000000000000e+00 +9.686274509803921573e-01 1.000000000000000000e+00 +9.725490196078431238e-01 1.000000000000000000e+00 +9.764705882352940902e-01 1.000000000000000000e+00 +9.803921568627450567e-01 1.000000000000000000e+00 +9.843137254901960231e-01 1.000000000000000000e+00 +9.882352941176471006e-01 1.000000000000000000e+00 +9.921568627450980671e-01 1.000000000000000000e+00 +9.960784313725490335e-01 1.000000000000000000e+00 +1.000000000000000000e+00 1.000000000000000000e+00 diff --git a/examples/shear_test/Nz_test.txt b/examples/shear_test/Nz_test.txt new file mode 100644 index 0000000..442ed84 --- /dev/null +++ b/examples/shear_test/Nz_test.txt @@ -0,0 +1,256 @@ +0.000000000000000000e+00 0.000000000000000000e+00 +3.921568627450980338e-03 1.068331278683803731e+00 +7.843137254901960675e-03 4.202845534734766630e+00 +1.176470588235294101e-02 9.254930211039754440e+00 +1.568627450980392135e-02 1.603891875026000946e+01 +1.960784313725490169e-02 2.434626750571996823e+01 +2.352941176470588203e-02 3.395501257952240337e+01 +2.745098039215686236e-02 4.463698985186998414e+01 +3.137254901960784270e-02 5.616361167977535018e+01 +3.529411764705882304e-02 6.831056002116484649e+01 +3.921568627450980338e-02 8.086159528748243019e+01 +4.313725490196078372e-02 9.361160969752239680e+01 +4.705882352941176405e-02 1.063690194580922395e+02 +5.098039215686274439e-02 1.189575719353623526e+02 +5.490196078431372473e-02 1.312176337047463619e+02 +5.882352941176470507e-02 1.430070190462103028e+02 +6.274509803921568540e-02 1.542014141184155278e+02 +6.666666666666666574e-02 1.646944486711946354e+02 +7.058823529411764608e-02 1.743974641945524695e+02 +7.450980392156862642e-02 1.832390245656867194e+02 +7.843137254901960675e-02 1.911642124002464982e+02 +8.235294117647058709e-02 1.981337513798920611e+02 +8.627450980392156743e-02 2.041229918026451742e+02 +9.019607843137254777e-02 2.091207935000871032e+02 +9.411764705882352811e-02 2.131283371153007806e+02 +9.803921568627450844e-02 2.161578915734747852e+02 +1.019607843137254888e-01 2.182315624411413069e+02 +1.058823529411764691e-01 2.193800427967983921e+02 +1.098039215686274495e-01 2.196413852588563884e+02 +1.137254901960784298e-01 2.190598109660013790e+02 +1.176470588235294101e-01 2.176845686050466213e+02 +1.215686274509803905e-01 2.155688540521034042e+02 +1.254901960784313708e-01 2.127687988495256093e+02 +1.294117647058823650e-01 2.093425335940052321e+02 +1.333333333333333315e-01 2.053493303665117651e+02 +1.372549019607843257e-01 2.008488265946616025e+02 +1.411764705882352922e-01 1.959003312012475533e+02 +1.450980392156862864e-01 1.905622125547342876e+02 +1.490196078431372528e-01 1.848913665916337550e+02 +1.529411764705882470e-01 1.789427625178360586e+02 +1.568627450980392135e-01 1.727690627055255845e+02 +1.607843137254902077e-01 1.664203127723483249e+02 +1.647058823529411742e-01 1.599436973472175509e+02 +1.686274509803921684e-01 1.533833566791917349e+02 +1.725490196078431349e-01 1.467802590186515488e+02 +1.764705882352941291e-01 1.401721235799316787e+02 +1.803921568627450955e-01 1.335933888683011617e+02 +1.843137254901960898e-01 1.270752212087877098e+02 +1.882352941176470562e-01 1.206455584374708963e+02 +1.921568627450980504e-01 1.143291838958573123e+02 +1.960784313725490169e-01 1.081478260949468222e+02 +2.000000000000000111e-01 1.021202796775702950e+02 +2.039215686274509776e-01 9.626254359641313840e+01 +2.078431372549019718e-01 9.058797273261032501e+01 +2.117647058823529382e-01 8.510743949860280111e+01 +2.156862745098039325e-01 7.982950229266842257e+01 +2.196078431372548989e-01 7.476057799562924799e+01 +2.235294117647058931e-01 6.990511601797362573e+01 +2.274509803921568596e-01 6.526577171406843547e+01 +2.313725490196078538e-01 6.084357727605765831e+01 +2.352941176470588203e-01 5.663810850090050053e+01 +2.392156862745098145e-01 5.264764608785791467e+01 +2.431372549019607809e-01 4.886933036919739237e+01 +2.470588235294117752e-01 4.529930860305245233e+01 +2.509803921568627416e-01 4.193287416377489762e+01 +2.549019607843137081e-01 3.876459715164177311e+01 +2.588235294117647300e-01 3.578844611065179038e+01 +2.627450980392156965e-01 3.299790069083670829e+01 +2.666666666666666630e-01 3.038605522071579301e+01 +2.705882352941176294e-01 2.794571326711452386e+01 +2.745098039215686514e-01 2.566947335453546231e+01 +2.784313725490196179e-01 2.354980609572345074e+01 +2.823529411764705843e-01 2.157912305015505083e+01 +2.862745098039215508e-01 1.974983767911588828e+01 +2.901960784313725727e-01 1.805441880600607973e+01 +2.941176470588235392e-01 1.648543701974709919e+01 +2.980392156862745057e-01 1.503560447882364848e+01 +3.019607843137254721e-01 1.369780858472158513e+01 +3.058823529411764941e-01 1.246513999740432865e+01 +3.098039215686274606e-01 1.133091546303356800e+01 +3.137254901960784270e-01 1.028869591634609826e+01 +3.176470588235293935e-01 9.332300307841872922e+00 +3.215686274509804154e-01 8.455815590035038909e+00 +3.254901960784313819e-01 7.653603278218902517e+00 +3.294117647058823484e-01 6.920302980168123064e+00 +3.333333333333333148e-01 6.250833266550138489e+00 +3.372549019607843368e-01 5.640390230072614308e+00 +3.411764705882353033e-01 5.084444057022494690e+00 +3.450980392156862697e-01 4.578733910255725625e+00 +3.490196078431372362e-01 4.119261398221039983e+00 +3.529411764705882582e-01 3.702282880533814158e+00 +3.568627450980392246e-01 3.324300837198350234e+00 +3.607843137254901911e-01 2.982054506001525773e+00 +3.647058823529411575e-01 2.672509971020624597e+00 +3.686274509803921795e-01 2.392849864719143405e+00 +3.725490196078431460e-01 2.140462826828922971e+00 +3.764705882352941124e-01 1.912932845190227216e+00 +3.803921568627450789e-01 1.708028586972206053e+00 +3.843137254901961009e-01 1.523692813232671028e+00 +3.882352941176470673e-01 1.358031955588851591e+00 +3.921568627450980338e-01 1.209305920834413017e+00 +3.960784313725490002e-01 1.075918177615324156e+00 +4.000000000000000222e-01 9.564061687214104879e-01 +4.039215686274509887e-01 8.494320831066718425e-01 +4.078431372549019551e-01 7.537740133600058234e-01 +4.117647058823529216e-01 6.683175169441297747e-01 +4.156862745098039436e-01 5.920475930385402652e-01 +4.196078431372549100e-01 5.240410811941168534e-01 +4.235294117647058765e-01 4.634594831653214286e-01 +4.274509803921568429e-01 4.095422051641463557e-01 +4.313725490196078649e-01 3.616002143129065649e-01 +4.352941176470588314e-01 3.190101001981088813e-01 +4.392156862745097978e-01 2.812085300848371494e-01 +4.431372549019607643e-01 2.476870844834671015e-01 +4.470588235294117863e-01 2.179874583144098543e-01 +4.509803921568627527e-01 1.916970118418715729e-01 +4.549019607843137192e-01 1.684446547976964359e-01 +4.588235294117646856e-01 1.478970466482674717e-01 +4.627450980392157076e-01 1.297550957316845610e-01 +4.666666666666666741e-01 1.137507399731469832e-01 +4.705882352941176405e-01 9.964399204092909412e-02 +4.745098039215686070e-01 8.722023210434050322e-02 +4.784313725490196290e-01 7.628773177213844792e-02 +4.823529411764705954e-01 6.667539330178526291e-02 +4.862745098039215619e-01 5.823068875588902765e-02 +4.901960784313725283e-01 5.081778442405252699e-02 +4.941176470588235503e-01 4.431583651030910481e-02 +4.980392156862745168e-01 3.861744479463392193e-02 +5.019607843137254832e-01 3.362725169997348579e-02 +5.058823529411764497e-01 2.926067492381023358e-02 +5.098039215686274161e-01 2.544276251712017031e-02 +5.137254901960783826e-01 2.210716000672708487e-02 +5.176470588235294601e-01 1.919517985382571817e-02 +5.215686274509804266e-01 1.665496421729987728e-02 +5.254901960784313930e-01 1.444073264183820191e-02 +5.294117647058823595e-01 1.251210691507246599e-02 +5.333333333333333259e-01 1.083350593309243340e-02 +5.372549019607842924e-01 9.373603978450788590e-03 +5.411764705882352589e-01 8.104846348429190633e-03 +5.450980392156862253e-01 7.003016773632588475e-03 +5.490196078431373028e-01 6.046851538039785257e-03 +5.529411764705882693e-01 5.217695651920266527e-03 +5.568627450980392357e-01 4.499196839238384094e-03 +5.607843137254902022e-01 3.877033482232177582e-03 +5.647058823529411686e-01 3.338673018851970667e-03 +5.686274509803921351e-01 2.873157614882566619e-03 +5.725490196078431016e-01 2.470914233135632597e-03 +5.764705882352940680e-01 2.123586498432896742e-03 +5.803921568627451455e-01 1.823886010579858861e-03 +5.843137254901961120e-01 1.565460989537466362e-03 +5.882352941176470784e-01 1.342780348910457833e-03 +5.921568627450980449e-01 1.151031487039061701e-03 +5.960784313725490113e-01 9.860302607219010359e-04 +5.999999999999999778e-01 8.441417661875919628e-04 +6.039215686274509443e-01 7.222106965933727590e-04 +6.078431372549019107e-01 6.175001762303809760e-04 +6.117647058823529882e-01 5.276380898647893081e-04 +6.156862745098039547e-01 4.505700322908598367e-04 +6.196078431372549211e-01 3.845180992011551914e-04 +6.235294117647058876e-01 3.279448268143635463e-04 +6.274509803921568540e-01 2.795216652042206685e-04 +6.313725490196078205e-01 2.381014397453117546e-04 +6.352941176470587870e-01 2.026943172763924223e-04 +6.392156862745097534e-01 1.724468491658003207e-04 +6.431372549019608309e-01 1.466237130795824376e-04 +6.470588235294117974e-01 1.245918194827885263e-04 +6.509803921568627638e-01 1.058064882811769216e-04 +6.549019607843137303e-01 8.979943602210142535e-05 +6.588235294117646967e-01 7.616834516440906485e-05 +6.627450980392156632e-01 6.456781450297631626e-05 +6.666666666666666297e-01 5.470151426188767574e-05 +6.705882352941175961e-01 4.631539098492538426e-05 +6.745098039215686736e-01 3.919178645477776422e-05 +6.784313725490196401e-01 3.314435173585888355e-05 +6.823529411764706065e-01 2.801365230556775355e-05 +6.862745098039215730e-01 2.366337333644052165e-05 +6.901960784313725394e-01 1.997704571556435543e-05 +6.941176470588235059e-01 1.685522351611885245e-05 +6.980392156862744724e-01 1.421305252851508929e-05 +7.019607843137254388e-01 1.197817725800544497e-05 +7.058823529411765163e-01 1.008894062922757655e-05 +7.098039215686274828e-01 8.492836619446896473e-06 +7.137254901960784492e-01 7.145181272354444791e-06 +7.176470588235294157e-01 6.007972113007999240e-06 +7.215686274509803821e-01 5.048909971613468607e-06 +7.254901960784313486e-01 4.240560700044185322e-06 +7.294117647058823151e-01 3.559637292874652092e-06 +7.333333333333332815e-01 2.986385559613448928e-06 +7.372549019607843590e-01 2.504058785635889422e-06 +7.411764705882353255e-01 2.098468809150214666e-06 +7.450980392156862919e-01 1.757602668363236797e-06 +7.490196078431372584e-01 1.471295470264241541e-06 +7.529411764705882248e-01 1.230951429487096504e-06 +7.568627450980391913e-01 1.029306148332517077e-06 +7.607843137254901578e-01 8.602241798444023556e-07 +7.647058823529411242e-01 7.185267546450238944e-07 +7.686274509803922017e-01 5.998452763758125087e-07 +7.725490196078431682e-01 5.004968152132526169e-07 +7.764705882352941346e-01 4.173783672596259759e-07 +7.803921568627451011e-01 3.478771111857133678e-07 +7.843137254901960675e-01 2.897942923791173966e-07 +7.882352941176470340e-01 2.412807077806807514e-07 +7.921568627450980005e-01 2.007820591819846381e-07 +7.960784313725489669e-01 1.669926956217001312e-07 +8.000000000000000444e-01 1.388164823959390845e-07 +8.039215686274510109e-01 1.153337200636120377e-07 +8.078431372549019773e-01 9.577319599287327950e-08 +8.117647058823529438e-01 7.948858718965636986e-08 +8.156862745098039102e-01 6.593854959994755781e-08 +8.196078431372548767e-01 5.466992857062798071e-08 +8.235294117647058432e-01 4.530361009480943141e-08 +8.274509803921568096e-01 3.752260493086025512e-08 +8.313725490196078871e-01 3.106201945687400082e-08 +8.352941176470588536e-01 2.570061974272315322e-08 +8.392156862745098200e-01 2.125374011326026221e-08 +8.431372549019607865e-01 1.756732557583797718e-08 +8.470588235294117530e-01 1.451292986954370597e-08 +8.509803921568627194e-01 1.198351840006535751e-08 +8.549019607843136859e-01 9.889948671011918400e-09 +8.588235294117646523e-01 8.158020624846170579e-09 +8.627450980392157298e-01 6.726006090731369703e-09 +8.666666666666666963e-01 5.542580752661560259e-09 +8.705882352941176627e-01 4.565094083835654263e-09 +8.745098039215686292e-01 3.758122870723129987e-09 +8.784313725490195957e-01 3.092262552773033892e-09 +8.823529411764705621e-01 2.543117870044440305e-09 +8.862745098039215286e-01 2.090460444577560984e-09 +8.901960784313724950e-01 1.717526095212250243e-09 +8.941176470588235725e-01 1.410429047029721564e-09 +8.980392156862745390e-01 1.157673870777880511e-09 +9.019607843137255054e-01 9.497490808050514523e-10 +9.058823529411764719e-01 7.787889223463260382e-10 +9.098039215686274384e-01 6.382920669460394523e-10 +9.137254901960784048e-01 5.228877731521405033e-10 +9.176470588235293713e-01 4.281416132457059891e-10 +9.215686274509803377e-01 3.503941620920148764e-10 +9.254901960784314152e-01 2.866271304675949711e-10 +9.294117647058823817e-01 2.343523355936235116e-10 +9.333333333333333481e-01 1.915196640904271920e-10 +9.372549019607843146e-01 1.564408207850728153e-10 +9.411764705882352811e-01 1.277261906741403528e-10 +9.450980392156862475e-01 1.042325876499478982e-10 +9.490196078431372140e-01 8.502003647272707900e-11 +9.529411764705881804e-01 6.931604580390590614e-11 +9.568627450980392579e-01 5.648608990538612681e-11 +9.607843137254902244e-01 4.600923325774715359e-11 +9.647058823529411908e-01 3.745801291169126001e-11 +9.686274509803921573e-01 3.048184378203772682e-11 +9.725490196078431238e-01 2.479333728695467987e-11 +9.764705882352940902e-01 2.015702788779434400e-11 +9.803921568627450567e-01 1.638008868160185728e-11 +9.843137254901960231e-01 1.330468915579979407e-11 +9.882352941176471006e-01 1.080170797232554461e-11 +9.921568627450980671e-01 8.765563245858261760e-12 +9.960784313725490335e-01 7.109963921283296938e-12 +1.000000000000000000e+00 5.764419961355172104e-12 diff --git a/examples/shear_test/param.cfg b/examples/shear_test/param.cfg new file mode 100644 index 0000000..6029a6e --- /dev/null +++ b/examples/shear_test/param.cfg @@ -0,0 +1,56 @@ +global: +{ + prefix_out= "examples/shear_test/out"; + output_format= "FITS"; + output_density= false + pk_filename= "examples/shear_test/Pk_CAMB_test.dat" + z_min= 0.001 + z_max= 0.400 + seed= 1001 + write_pred=false + just_write_pred=false + pred_dz=0.1 +} + +field_par: +{ + r_smooth= 5. + smooth_potential= true + n_grid= 512 + dens_type= 0 + lpt_buffer_fraction= 0.6 + lpt_interp_type= 1 + output_lpt= 0 +} + +cosmo_par: +{ + omega_M= 0.3 + omega_L= 0.7 + omega_B= 0.05 + h= 0.7 + w= -1.0 + ns= 0.96 + sigma_8= 0.8 +} + +srcs1: +{ + nz_filename= "examples/shear_test/Nz_test.txt" + bias_filename= "examples/shear_test/Bz_test.txt" + include_shear= false + dr_shear= 100. + store_skewers= false +} + +kappa: +{ + z_out= [2.076379E-01] + nside= 128 +} + +shear: +{ + nside= 128 + dr_shear= 100. +} diff --git a/examples/simple/mk_files.py b/examples/simple/mk_files.py index eaeca98..30b2fa2 100755 --- a/examples/simple/mk_files.py +++ b/examples/simple/mk_files.py @@ -2,44 +2,48 @@ import numpy as np import matplotlib.pyplot as plt -nz=256 -z0=0 -zf=1.0 -zx=0.09 -pow1=1.5 -ngal=1.4E6 -zx2=0.60 -ngal2=2e6 -pow2=4.0 +nz = 256 +z0 = 0 +zf = 1.0 +zx = 0.09 +pow1 = 1.5 +ngal = 1.4E6 +zx2 = 0.60 +ngal2 = 2e6 +pow2 = 4.0 -nux=1420. -n_nu=15 -z0nu=0.1 -zfnu=0.5 -nu0=nux/(1+zfnu) -nuf=nux/(1+z0nu) -dnu=(nuf-nu0)/n_nu -nu0_arr=nu0+dnu*np.arange(n_nu) -nuf_arr=nu0+dnu*(1+np.arange(n_nu)) +nux = 1420. +n_nu = 15 +z0nu = 0.1 +zfnu = 0.5 +nu0 = nux/(1+zfnu) +nuf = nux/(1+z0nu) +dnu = (nuf-nu0)/n_nu +nu0_arr = nu0+dnu*np.arange(n_nu) +nuf_arr = nu0+dnu*(1+np.arange(n_nu)) -zarr=z0+(zf-z0)*np.arange(nz)/(nz-1.) -nzarr=(zarr/zx)**2*np.exp(-(zarr/zx)**pow1); -nzarr2=(zarr/zx2)**2*np.exp(-(zarr/zx2)**pow2); -bzarr=np.ones_like(zarr) -bzarr2=2*np.ones_like(zarr) -tzarr=5.5919E-2+2.3242E-1*zarr-2.4136E-2*zarr**2. -norm=ngal/(4*np.pi*(180/np.pi)**2*np.sum(nzarr)*(zf-z0)/nz) -nzarr*=norm -norm2=ngal2/(4*np.pi*(180/np.pi)**2*np.sum(nzarr2)*(zf-z0)/nz) -nzarr2*=norm2 +zarr = z0+(zf-z0)*np.arange(nz)/(nz-1.) +nzarr = (zarr/zx)**2*np.exp(-(zarr/zx)**pow1) +nzarr2 = (zarr/zx2)**2*np.exp(-(zarr/zx2)**pow2) +bzarr = np.ones_like(zarr) +bzarr2 = 2*np.ones_like(zarr) +tzarr = 5.5919E-2+2.3242E-1*zarr-2.4136E-2*zarr**2. +norm = ngal/(4*np.pi*(180/np.pi)**2*np.sum(nzarr)*(zf-z0)/nz) +nzarr *= norm +norm2 = ngal2/(4*np.pi*(180/np.pi)**2*np.sum(nzarr2)*(zf-z0)/nz) +nzarr2 *= norm2 -np.savetxt("Nz_test.txt",np.transpose([zarr,nzarr])) -np.savetxt("Nz2_test.txt",np.transpose([zarr,nzarr2])) -np.savetxt("Tz_test.txt",np.transpose([zarr,tzarr])) -np.savetxt("Bz_test.txt",np.transpose([zarr,bzarr])) -np.savetxt("Bz2_test.txt",np.transpose([zarr,bzarr2])) -np.savetxt("nuTable.txt",np.transpose([nu0_arr,nuf_arr])) +np.savetxt("Nz_test.txt", np.transpose([zarr, nzarr])) +np.savetxt("Nz2_test.txt", np.transpose([zarr, nzarr2])) +np.savetxt("Tz_test.txt", np.transpose([zarr, tzarr])) +np.savetxt("Bz_test.txt", np.transpose([zarr, bzarr])) +np.savetxt("Bz2_test.txt", np.transpose([zarr, bzarr2])) +np.savetxt("nuTable.txt", np.transpose([nu0_arr, nuf_arr])) -plt.plot(zarr,nzarr); plt.plot(zarr,nzarr2); plt.show() -plt.plot(zarr,bzarr); plt.show() -plt.plot(zarr,tzarr); plt.show() +plt.plot(zarr, nzarr) +plt.plot(zarr, nzarr2) +plt.show() +plt.plot(zarr, bzarr) +plt.show() +plt.plot(zarr, tzarr) +plt.show() diff --git a/examples/simple/param.cfg b/examples/simple/param.cfg index a731d9d..c044723 100644 --- a/examples/simple/param.cfg +++ b/examples/simple/param.cfg @@ -83,7 +83,7 @@ srcs1: # 1-> z, 2-> b(z) bias_filename= "examples/simple/Bz_test.txt" #Do you want to include shear ellipticities? - include_shear= false + include_shear= true #Do you want to store line-of-sight skewers for each object? store_skewers= true #You can also store the Gaussian density field (instead of the non-Gaussian one) diff --git a/examples/simple/read_skewers.py b/examples/simple/read_skewers.py index 961c82d..c920a7a 100644 --- a/examples/simple/read_skewers.py +++ b/examples/simple/read_skewers.py @@ -4,40 +4,46 @@ hdulist = fits.open('out_srcs_s1_0.fits') -#First HDU contains the source catalog +# First HDU contains the source catalog print(hdulist[1].header.keys) -plt.figure(); plt.hist(hdulist[1].data['Z_COSMO'],bins=100) +plt.figure() +plt.hist(hdulist[1].data['Z_COSMO'], bins=100) print(" ") -#Second HDU contains the density skewers as a FITS image -#The skewers have the same ordering as the sources in the catalog -#(i.e. skewer hdulist[2].data[i,:] corresponds to source hdulist[1].data[i]) -id=np.argmax(hdulist[1].data['Z_COSMO']) +# Second HDU contains the density skewers as a FITS image +# The skewers have the same ordering as the sources in the catalog +# (i.e. skewer hdulist[2].data[i, :] corresponds to source hdulist[1].data[i]) +id = np.argmax(hdulist[1].data['Z_COSMO']) print(hdulist[2].header.keys) -plt.figure(); plt.plot(hdulist[4].data['R'],hdulist[2].data[id]); -plt.xlabel('$r\\,\\,[{\\rm Mpc}/h]$',fontsize=18); -plt.ylabel('$\\delta$',fontsize=18) +plt.figure() +plt.plot(hdulist[4].data['R'], hdulist[2].data[id]) +plt.xlabel('$r\\, \\, [{\\rm Mpc}/h]$', fontsize=18) +plt.ylabel('$\\delta$', fontsize=18) print(" ") -#Third HDU contains the velocity skewers. The units of the velocity are -#such that the skewers contain the redshift distortion associated with -#the peculiar velocity field +# Third HDU contains the velocity skewers. The units of the velocity are +# such that the skewers contain the redshift distortion associated with +# the peculiar velocity field print(hdulist[3].header.keys) -plt.figure(); plt.plot(hdulist[4].data['R'],hdulist[3].data[id]); -plt.xlabel('$r\\,\\,[{\\rm Mpc}/h]$',fontsize=18); -plt.ylabel('$\\delta z_{\\rm RSD}$',fontsize=18) +plt.figure() +plt.plot(hdulist[4].data['R'], hdulist[3].data[id]) +plt.xlabel('$r\\, \\, [{\\rm Mpc}/h]$', fontsize=18) +plt.ylabel('$\\delta z_{\\rm RSD}$', fontsize=18) print(" ") -#Fourth HDU is a table containing background cosmological quantities at -#the distances where the skewers are sampled (see the use of -#hdulist[4].data['R'] in the previous examples +# Fourth HDU is a table containing background cosmological quantities at +# the distances where the skewers are sampled (see the use of +# hdulist[4].data['R'] in the previous examples print(hdulist[4].header.keys) -plt.figure(); -plt.plot(hdulist[4].data['Z'],hdulist[4].data['R']*0.001,label='$r(z)\\,[{\\rm Gpc}/h]$') -plt.plot(hdulist[4].data['Z'],hdulist[4].data['D'],label='$D_\\delta(z)$') -plt.plot(hdulist[4].data['Z'],hdulist[4].data['V'],label='$D_v(z)$') +plt.figure() +plt.plot(hdulist[4].data['Z'], hdulist[4].data['R']*0.001, + label='$r(z)\\, [{\\rm Gpc}/h]$') +plt.plot(hdulist[4].data['Z'], hdulist[4].data['D'], + label='$D_\\delta(z)$') +plt.plot(hdulist[4].data['Z'], hdulist[4].data['V'], + label='$D_v(z)$') plt.legend(loc='lower right') -plt.xlabel('$z$',fontsize=18) +plt.xlabel('$z$', fontsize=18) print(" ") plt.show() diff --git a/src/common.c b/src/common.c index 8cfe787..0b58d10 100644 --- a/src/common.c +++ b/src/common.c @@ -519,13 +519,15 @@ void catalog_cartesian_free(CatalogCartesian *cat) free(cat); } -Catalog *catalog_alloc(int nsrcs,int has_skw,int skw_gauss,double rmax,int ng) +Catalog *catalog_alloc(int nsrcs,int has_shear,int has_skw, + int skw_gauss,double rmax,int ng) { Catalog *cat=my_malloc(sizeof(Catalog)); if(nsrcs>0) { cat->nsrc=nsrcs; cat->srcs=my_malloc(nsrcs*sizeof(Src)); + cat->has_shear=has_shear; cat->has_skw=has_skw; cat->skw_gauss=skw_gauss; get_radial_params(rmax,ng,&(cat->nr),&(cat->dr)); diff --git a/src/common.h b/src/common.h index c1958b2..a4d7f75 100644 --- a/src/common.h +++ b/src/common.h @@ -175,6 +175,7 @@ typedef struct { double rmax; double dr; double idr; + int has_shear; int has_skw; int skw_gauss; float *d_skw; @@ -370,7 +371,7 @@ HealpixShells *hp_shell_alloc(int nq,int nside,int nside_base,int nr); void hp_shell_free(HealpixShells *shell); CatalogCartesian *catalog_cartesian_alloc(int nsrcs); void catalog_cartesian_free(CatalogCartesian *cat); -Catalog *catalog_alloc(int nsrcs,int has_skw,int skw_gauss,double rmax,int ng); +Catalog *catalog_alloc(int nsrcs,int has_shear,int has_skw,int skw_gauss,double rmax,int ng); void catalog_free(Catalog *cat); static inline double bias_model(double d,double b) diff --git a/src/io.c b/src/io.c index 056d893..de4013e 100644 --- a/src/io.c +++ b/src/io.c @@ -966,7 +966,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) char *ttype[]={"TYPE","RA" ,"DEC","Z_COSMO","DZ_RSD","E1","E2"}; char *tform[]={"1J" ,"1E" ,"1E" ,"1E" ,"1E" ,"1E","1E"}; char *tunit[]={"NA" ,"DEG","DEG","NA" ,"NA" ,"NA","NA"}; - if(par->do_srcs_shear) + if(par->cats[ipop]->has_shear) nfields=7; print_info(" %d-th population (FITS)\n",ipop); @@ -999,7 +999,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) dec_arr[ii]=par->cats[ipop]->srcs[row_here+ii].dec; z0_arr[ii]=par->cats[ipop]->srcs[row_here+ii].z0; rsd_arr[ii]=par->cats[ipop]->srcs[row_here+ii].dz_rsd; - if(par->do_srcs_shear) { + if(par->cats[ipop]->has_shear) { e1_arr[ii]=par->cats[ipop]->srcs[row_here+ii].e1; e2_arr[ii]=par->cats[ipop]->srcs[row_here+ii].e2; } @@ -1009,7 +1009,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) fits_write_col(fptr,TFLOAT,3,row_here+1,1,nrw_here,dec_arr,&status); fits_write_col(fptr,TFLOAT,4,row_here+1,1,nrw_here,z0_arr,&status); fits_write_col(fptr,TFLOAT,5,row_here+1,1,nrw_here,rsd_arr,&status); - if(par->do_srcs_shear) { + if(par->cats[ipop]->has_shear) { fits_write_col(fptr,TFLOAT,6,row_here+1,1,nrw_here,e1_arr,&status); fits_write_col(fptr,TFLOAT,7,row_here+1,1,nrw_here,e2_arr,&status); } @@ -1092,7 +1092,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) FILE *fil=fopen(fname,"w"); if(fil==NULL) error_open_file(fname); fprintf(fil,"#[1]type [2]RA, [3]dec, [4]z0, [5]dz_RSD "); - if(par->do_srcs_shear) + if(par->cats[ipop]->has_shear) fprintf(fil,"#[6]e1, [7]e2\n"); else fprintf(fil,"\n"); @@ -1100,7 +1100,7 @@ static void write_catalog(ParamCoLoRe *par,int ipop) fprintf(fil,"%d %E %E %E %E ", ipop,par->cats[ipop]->srcs[jj].ra,par->cats[ipop]->srcs[jj].dec, par->cats[ipop]->srcs[jj].z0,par->cats[ipop]->srcs[jj].dz_rsd); - if(par->do_srcs_shear) + if(par->cats[ipop]->has_shear) fprintf(fil,"%E %E \n",par->cats[ipop]->srcs[jj].e1,par->cats[ipop]->srcs[jj].e2); else fprintf(fil,"\n"); diff --git a/src/srcs.c b/src/srcs.c index 3373b04..1af985e 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -340,7 +340,10 @@ void srcs_distribute(ParamCoLoRe *par) static void srcs_get_local_properties_single(ParamCoLoRe *par,int ipop) { - par->cats[ipop]=catalog_alloc(par->cats_c[ipop]->nsrc,par->skw_srcs[ipop],par->skw_gauss[ipop], + par->cats[ipop]=catalog_alloc(par->cats_c[ipop]->nsrc, + par->shear_srcs[ipop], + par->skw_srcs[ipop], + par->skw_gauss[ipop], par->r_max,par->n_grid); #ifdef _HAVE_OMP @@ -416,7 +419,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) double *fac_r_1=NULL,*fac_r_2=NULL; //Kernels for the LOS integrals - if(par->do_srcs_shear) { + if(cat->has_shear) { fac_r_1=my_malloc(cat->nr*sizeof(double)); fac_r_2=my_malloc(cat->nr*sizeof(double)); @@ -479,7 +482,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) } //Compute lensing shear - if(par->do_srcs_shear) { + if(cat->has_shear) { //Compute linear transformations needed for shear double r1[6],r2[6]; double cth_h=1,sth_h=0,cph_h=1,sph_h=0; @@ -532,7 +535,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) cat->srcs[ip].e2+=e2; } }//end omp for - if(par->do_srcs_shear) { + if(cat->has_shear) { free(fac_r_1); free(fac_r_2); } From bbd3b874fab134362414f02e5826b8ae5e3f9f15 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Tue, 21 Jan 2020 10:58:55 +0000 Subject: [PATCH 07/29] implemented --- Makefile | 2 + examples/shear_test/param.cfg | 5 +- src/io.c | 2 - src/shear.c | 2 +- src/srcs.c | 99 +++++++++++++++++++++++++++++++++-- 5 files changed, 101 insertions(+), 9 deletions(-) diff --git a/Makefile b/Makefile index cba6053..59221e2 100644 --- a/Makefile +++ b/Makefile @@ -12,6 +12,8 @@ DEFINEFLAGS += -D_LONGIDS #DEFINEFLAGS += -D_BIAS_MODEL_2 #Use linear bias model #DEFINEFLAGS += -D_BIAS_MODEL_3 +#Use new shear method +DEFINEFLAGS += -D_USE_NEW_LENSING #Generate debug help. Only useful for development DEFINEFLAGS += -D_DEBUG #Use double precision floating point? Set to "yes" or "no" diff --git a/examples/shear_test/param.cfg b/examples/shear_test/param.cfg index 6029a6e..89db182 100644 --- a/examples/shear_test/param.cfg +++ b/examples/shear_test/param.cfg @@ -38,8 +38,7 @@ srcs1: { nz_filename= "examples/shear_test/Nz_test.txt" bias_filename= "examples/shear_test/Bz_test.txt" - include_shear= false - dr_shear= 100. + include_shear= true store_skewers= false } @@ -51,6 +50,6 @@ kappa: shear: { - nside= 128 + nside= 256 dr_shear= 100. } diff --git a/src/io.c b/src/io.c index de4013e..64310c0 100644 --- a/src/io.c +++ b/src/io.c @@ -499,8 +499,6 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) for(i_r=0;i_rn_shear;i_r++) { par->smap->r0[i_r]=r_arr[i_r]; par->smap->rf[i_r]=r_arr[i_r]; - printf("%d %lE %lE\n",i_r,r_arr[i_r], - get_bg(par,r_arr[i_r],BG_Z,-1)); } free(r_arr); } diff --git a/src/shear.c b/src/shear.c index 12ce826..335c37f 100644 --- a/src/shear.c +++ b/src/shear.c @@ -66,7 +66,7 @@ void shear_beams_preproc(ParamCoLoRe *par) #ifdef _HAVE_OMP #pragma omp for #endif //_HAVE_OMP - for(ipp=0;ippsmap->num_pix*par->smap->nr;ipp++) { + for(ipp=0;ipp<2*par->smap->num_pix*par->smap->nr;ipp++) { par->smap->data[ipp]=0; par->smap->nadd[ipp]=1; } //end omp for diff --git a/src/srcs.c b/src/srcs.c index 1af985e..d3a8c13 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -21,6 +21,48 @@ /////////////////////////////////////////////////////////////////////// #include "common.h" +static int get_r_index_smap(HealpixShells *sh,double r,int ir_start) +{ + int gotit=0; + int ir0; + if(ir_start<0) + ir0=0; + else if(ir_start>=sh->nr) + ir0=sh->nr-1; + else + ir0=ir_start; + + while(!gotit) { + if(ir0==0) { + if(rr0[1]) + gotit=1; + else + ir0++; + } + else if(ir0==sh->nr-1) { + if(r>=sh->rf[sh->nr-2]) { + ir0=sh->nr-2; + gotit=1; + } + else + ir0--; + } + else { + if(rr0[ir0]) + ir0--; + else { + if(r>=sh->rf[ir0+1]) + ir0++; + else + gotit=1; + } + } + } + + return ir0; +} + + static inline void cart2sph(double x,double y,double z,double *r,double *cth,double *phi) { *r=sqrt(x*x+y*y+z*z); @@ -419,6 +461,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) double *fac_r_1=NULL,*fac_r_2=NULL; //Kernels for the LOS integrals +#ifndef _USE_NEW_LENSING if(cat->has_shear) { fac_r_1=my_malloc(cat->nr*sizeof(double)); fac_r_2=my_malloc(cat->nr*sizeof(double)); @@ -430,6 +473,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) fac_r_2[ip]=rm*rm*pg*cat->dr; } } +#endif //_USE_NEW_LENSING #ifdef _HAVE_OMP #pragma omp for @@ -481,6 +525,7 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) } } +#ifndef _USE_NEW_LENSING //Compute lensing shear if(cat->has_shear) { //Compute linear transformations needed for shear @@ -534,11 +579,14 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) cat->srcs[ip].e1+=e1; cat->srcs[ip].e2+=e2; } +#endif //_USE_NEW_LENSING }//end omp for +#ifndef _USE_NEW_LENSING if(cat->has_shear) { free(fac_r_1); free(fac_r_2); } +#endif //_USE_NEW_LENSING }//end omp parallel } @@ -552,14 +600,24 @@ void srcs_get_beam_properties(ParamCoLoRe *par) static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) { Catalog *cat=par->cats[ipop]; + CatalogCartesian *catc=par->cats_c[ipop]; #ifdef _HAVE_OMP -#pragma omp parallel default(none) \ - shared(par,cat) +#pragma omp parallel default(none) \ + shared(par,cat,catc,NNodes,NodeThis) #endif //_HAVE_OMP { int ii; double factor_vel=-par->fgrowth_0/(1.5*par->hubble_0*par->OmegaM); + int ir_s=0; +#ifdef _USE_NEW_LENSING + int nside_ratio=1; + HealpixShells *smap=NULL; + if(cat->has_shear) { + smap = par->smap; + nside_ratio=smap->nside/par->nside_base; + } +#endif //_USE_NEW_LENSING #ifdef _HAVE_OMP #pragma omp for @@ -573,7 +631,42 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) cat->srcs[ii].dz_rsd*=vg*factor_vel; //Shear - //Nothing else to do here + if(cat->has_shear) { +#ifdef _USE_NEW_LENSING + int ax; + long ipix,ibase,ibase_here; + double u[3]; + double g1_0,g1_f,g2_0,g2_f; + //Find lower shell index + ir_s=get_r_index_smap(smap, r, ir_s); + + //Find intervals + double h = (r - smap->r0[ir_s]) / (smap->r0[ir_s] - smap->r0[ir_s+1]); + + //Find pixel index + //Find base this galaxy belongs to + for(ax=0;ax<3;ax++) + u[ax]=catc->pos[NPOS_CC*ii+ax]; + vec2pix_nest(par->nside_base,u,&ibase); + if(ibase%NNodes!=NodeThis) + report_error(1,"Bad base!!\n"); + //Find base index within this node + ibase_here=(ibase-NodeThis)/NNodes; + vec2pix_nest(smap->nside,u,&ipix); + //Offset to pixel indices in this node + ipix-=(ibase-ibase_here)*nside_ratio*nside_ratio; + + //Find shear at edges + g1_0=smap->data[2*(ir_s*smap->num_pix+ipix)+0]; + g2_0=smap->data[2*(ir_s*smap->num_pix+ipix)+1]; + g1_f=smap->data[2*((ir_s+1)*smap->num_pix+ipix)+0]; + g2_f=smap->data[2*((ir_s+1)*smap->num_pix+ipix)+1]; + + //Interpolate + cat->srcs[ii].e1=g1_0*(1-h)+g1_f*h; + cat->srcs[ii].e2=g2_0*(1-h)+g2_f*h; +#endif //_USE_NEW_LENSING + } //Skewers if(cat->has_skw) { From 02b2416a5bc8812122256ec5917b91800c310cf3 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 22 Jan 2020 19:56:42 +0000 Subject: [PATCH 08/29] fixed a bunch of bugs --- src/density.c | 2 ++ src/healpix_extra.c | 1 + src/io.c | 12 +++++++++++- src/shear.c | 8 ++++++++ 4 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/density.c b/src/density.c index ba15876..1ce16a1 100644 --- a/src/density.c +++ b/src/density.c @@ -1331,10 +1331,12 @@ void compute_density_normalization(ParamCoLoRe *par) free(norm_srcs_arr[ipop]); gsl_spline_free(spline_norm_srcs[ipop]); } + free(norm_srcs_arr); gsl_interp_accel_free(intacc_srcs); for(ipop=0;ipopn_imap;ipop++) { free(norm_imap_arr[ipop]); gsl_spline_free(spline_norm_imap[ipop]); } + free(norm_imap_arr); gsl_interp_accel_free(intacc_imap); } diff --git a/src/healpix_extra.c b/src/healpix_extra.c index f292a03..f17fd97 100644 --- a/src/healpix_extra.c +++ b/src/healpix_extra.c @@ -53,6 +53,7 @@ void he_write_healpix_map(flouble **tmap,int nfields,long nside,char *fname,int free(ttype); free(tform); free(tunit); + free(map_dum); } flouble *he_read_healpix_map(char *fname,long *nside,int nfield) diff --git a/src/io.c b/src/io.c index 64310c0..5e9f39e 100644 --- a/src/io.c +++ b/src/io.c @@ -88,6 +88,7 @@ static ParamCoLoRe *param_colore_new(void) sprintf(par->prefixOut,"default"); par->output_format=0; par->do_pred=0; + par->just_do_pred=0; par->do_pred=1; //Tracers @@ -420,9 +421,15 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) par->need_beaming=par->do_srcs_shear+par->do_kappa+par->do_shear+par->do_isw+par->do_skewers; init_fftw(par); + cosmo_set(par); + + // We need to compute the number of maps here + if(par->do_shear) { + flouble *r_arr=compute_shear_spacing(par); + free(r_arr); + } get_max_memory(par,test_memory+par->just_do_pred); - cosmo_set(par); print_info("\n"); double dk=2*M_PI/par->l_box; @@ -508,6 +515,8 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) compute_tracer_cosmo(par); + free(conf); + return par; } @@ -1203,4 +1212,5 @@ void param_colore_free(ParamCoLoRe *par) #ifdef _DEBUG fclose(par->f_dbg); #endif //_DEBUG + free(par); } diff --git a/src/shear.c b/src/shear.c index 335c37f..7bdac30 100644 --- a/src/shear.c +++ b/src/shear.c @@ -68,6 +68,12 @@ void shear_beams_preproc(ParamCoLoRe *par) #endif //_HAVE_OMP for(ipp=0;ipp<2*par->smap->num_pix*par->smap->nr;ipp++) { par->smap->data[ipp]=0; + } //end omp for + +#ifdef _HAVE_OMP +#pragma omp for +#endif //_HAVE_OMP + for(ipp=0;ippsmap->num_pix*par->smap->nr;ipp++) { par->smap->nadd[ipp]=1; } //end omp for } //end omp parallel @@ -178,7 +184,9 @@ void shear_get_beam_properties(ParamCoLoRe *par) } } //end omp for +#ifdef _HAVE_OMP #pragma omp single +#endif //_HAVE_OMP { for(i_r=0;i_rnr;i_r++) { smap->r0[i_r]=1./inv_r_max[i_r]; From 46861c20274590667d33ae669cf086d49e2ce91e Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 12 Mar 2020 17:34:54 +0000 Subject: [PATCH 09/29] factor 2 --- src/shear.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shear.c b/src/shear.c index 7bdac30..254131a 100644 --- a/src/shear.c +++ b/src/shear.c @@ -133,7 +133,7 @@ void shear_get_beam_properties(ParamCoLoRe *par) double shear1_1=0,shear1_2=0; double shear2_1=0,shear2_2=0; double *u=&(smap->pos[3*ip]); - double prefac=idx*idx; + double prefac=2*idx*idx; //2/Dx^2 double cth_h=1,sth_h=0,cph_h=1,sph_h=0; cth_h=u[2]; From e47a5d0313b06ad757ac0692e2e009401c927d9f Mon Sep 17 00:00:00 2001 From: David Alonso Date: Sat, 14 Mar 2020 18:15:23 +0000 Subject: [PATCH 10/29] new spacing scheme --- examples/shear_test/param.cfg | 7 ++++++- src/common.h | 6 +++++- src/cosmo.c | 17 +++++++++++------ src/io.c | 14 ++++++++++---- 4 files changed, 32 insertions(+), 12 deletions(-) diff --git a/examples/shear_test/param.cfg b/examples/shear_test/param.cfg index 89db182..fef20fa 100644 --- a/examples/shear_test/param.cfg +++ b/examples/shear_test/param.cfg @@ -51,5 +51,10 @@ kappa: shear: { nside= 256 - dr_shear= 100. + # Select number of slices + # For these parameters, 11 corresponds roughly + # to dr_shear = 100 Mpc/h. + n_shear= 11 + # Select spacing type ("r" or "log(1+z)") + spacing_type= "r" } diff --git a/src/common.h b/src/common.h index a4d7f75..fb08619 100644 --- a/src/common.h +++ b/src/common.h @@ -73,6 +73,10 @@ #define RETURN_PDOT 8 #define RETURN_GAUSS 16 +// Radial spacing +#define SPACING_R 0 +#define SPACING_LOGZ 1 + //Interpolation type #ifndef INTERP_TYPE_SKW #define INTERP_TYPE_SKW INTERP_CIC @@ -319,7 +323,7 @@ typedef struct { #endif //_ADD_EXTRA_KAPPA // Shear int do_shear; //Do you want to output shear maps? - double dr_shear; //Separation between shear maps. + int shear_spacing_type; //log(1+z)? r? int n_shear; //How many maps? int nside_shear; HealpixShells *smap; //Shear maps at each redshift diff --git a/src/cosmo.c b/src/cosmo.c index 6cf7cec..d909f26 100644 --- a/src/cosmo.c +++ b/src/cosmo.c @@ -791,14 +791,19 @@ void compute_tracer_cosmo(ParamCoLoRe *par) flouble *compute_shear_spacing(ParamCoLoRe *par) { - //Figure out appropriate radial sampling flouble *rarr; - int ir,nr=(int)(par->r_max/par->dr_shear+1.); - par->dr_shear=par->r_max/nr; + int ir,nr=par->n_shear; rarr=my_malloc(nr*sizeof(flouble)); - for(ir=0;irdr_shear; + if(par->shear_spacing_type==SPACING_R) { + flouble dr=par->r_max/nr; + for(ir=0;irz_max)/nr; + for(ir=0;irn_shear=nr; return rarr; } diff --git a/src/io.c b/src/io.c index 5e9f39e..016b1f4 100644 --- a/src/io.c +++ b/src/io.c @@ -99,7 +99,7 @@ static ParamCoLoRe *param_colore_new(void) par->do_imap=0; par->do_kappa=0; par->do_shear=0; - par->dr_shear=100000.; + par->shear_spacing_type=SPACING_R; par->do_isw=0; par->n_srcs=-1; par->n_imap=-1; @@ -370,10 +370,16 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) if(cset!=NULL) { par->do_shear=1; par->do_srcs_shear=1; - conf_read_double(conf,"shear","dr_shear",&(par->dr_shear)); + conf_read_int(conf,"shear","n_shear",&(par->n_shear)); + char spacing_string[256]="default"; + conf_read_string(conf,"shear","spacing_type",spacing_string); + if(!strcmp(spacing_string,"r")) + par->shear_spacing_type=SPACING_R; + else if(!strcmp(spacing_string,"log(1+z)")) + par->shear_spacing_type=SPACING_LOGZ; + else + report_error(1,"Unknown spacing type %s\n",spacing_string); conf_read_int(conf,"shear","nside",&(par->nside_shear)); - if(par->dr_shear<=par->l_box/par->n_grid) - report_error(1,"Spacing between shear source planes is smaller than box resolution\n"); } // Check shear exists if requested with catalog if(par->do_srcs_shear && !(par->do_shear)) From ac4fdbee31e52353cc784f888c254588b6af553b Mon Sep 17 00:00:00 2001 From: David Alonso Date: Sat, 14 Mar 2020 18:31:26 +0000 Subject: [PATCH 11/29] small bug --- src/cosmo.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/cosmo.c b/src/cosmo.c index d909f26..3384d46 100644 --- a/src/cosmo.c +++ b/src/cosmo.c @@ -799,10 +799,10 @@ flouble *compute_shear_spacing(ParamCoLoRe *par) for(ir=0;irshear_spacing_type==SPACING_LOGZ) { flouble dlogz=log(1+par->z_max)/nr; for(ir=0;ir Date: Thu, 19 Mar 2020 19:08:59 +0000 Subject: [PATCH 12/29] stupid bug --- src/io.c | 15 +++++++++++++++ src/srcs.c | 2 +- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/io.c b/src/io.c index 016b1f4..d2b9969 100644 --- a/src/io.c +++ b/src/io.c @@ -848,6 +848,21 @@ void write_shear(ParamCoLoRe *par) free(map_write[1]); free(map_write); free(map_nadd); + + if(NodeThis==0) { + sprintf(fname,"%s_shear_r.txt", par->prefixOut); + FILE *f=fopen(fname, "w"); + if(f==NULL) + report_error(1,"Couldn't open file %s\n", fname); + for(ir=0;irsmap->nr;ir++) { + fprintf(f,"%d %lE %lE %lE %lE\n",ir, + par->smap->r0[ir], par->smap->rf[ir], + get_bg(par, par->smap->r0[ir], BG_Z, 0), + get_bg(par, par->smap->rf[ir], BG_Z, 0)); + } + fclose(f); + } + if(NodeThis==0) timer(2); print_info("\n"); } diff --git a/src/srcs.c b/src/srcs.c index d3a8c13..2e5c558 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -641,7 +641,7 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) ir_s=get_r_index_smap(smap, r, ir_s); //Find intervals - double h = (r - smap->r0[ir_s]) / (smap->r0[ir_s] - smap->r0[ir_s+1]); + double h = (r - smap->r0[ir_s]) / (smap->r0[ir_s+1] - smap->r0[ir_s]); //Find pixel index //Find base this galaxy belongs to From 4e22129a1be0a6b70a7164765e8ad532e9c59c51 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 30 Mar 2020 18:03:52 +0100 Subject: [PATCH 13/29] deactivate shear if not needed --- src/beaming.c | 10 ++++++++-- src/io.c | 8 ++++++++ src/main.c | 8 ++++++++ 3 files changed, 24 insertions(+), 2 deletions(-) diff --git a/src/beaming.c b/src/beaming.c index c172b7a..a208c75 100644 --- a/src/beaming.c +++ b/src/beaming.c @@ -300,8 +300,10 @@ void get_beam_properties(ParamCoLoRe *par) if(par->do_kappa) kappa_beams_preproc(par); +#ifdef _USE_NEW_LENSING if(par->do_shear) shear_beams_preproc(par); +#endif //_USE_NEW_LENSING if(par->do_isw) isw_beams_preproc(par); if(par->do_srcs) @@ -333,8 +335,10 @@ void get_beam_properties(ParamCoLoRe *par) if(par->do_kappa) kappa_get_beam_properties(par); +#ifdef _USE_NEW_LENSING if(par->do_shear) - shear_get_beam_properties(par); + shear_get_beam_properties(par); +#endif //_USE_NEW_LENSING if(par->do_isw) isw_get_beam_properties(par); if(par->do_srcs) @@ -348,8 +352,10 @@ void get_beam_properties(ParamCoLoRe *par) if(par->do_kappa) kappa_beams_postproc(par); +#ifdef _USE_NEW_LENSING if(par->do_shear) - shear_beams_postproc(par); + shear_beams_postproc(par); +#endif //_USE_NEW_LENSING if(par->do_isw) isw_beams_postproc(par); if(par->do_srcs) diff --git a/src/io.c b/src/io.c index d2b9969..925aed2 100644 --- a/src/io.c +++ b/src/io.c @@ -365,6 +365,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) conf_read_int(conf,"kappa","nside",&(par->nside_kappa)); } +#ifdef _USE_NEW_LENSING //Shear maps cset=config_lookup(conf,"shear"); if(cset!=NULL) { @@ -384,6 +385,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) // Check shear exists if requested with catalog if(par->do_srcs_shear && !(par->do_shear)) report_error(1,"Include a \"shear\" section if you want shear with your galaxies\n"); +#endif //_USE_NEW_LENSING cset=config_lookup(conf,"isw"); if(cset!=NULL) { @@ -430,10 +432,12 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) cosmo_set(par); // We need to compute the number of maps here +#ifdef _USE_NEW_LENSING if(par->do_shear) { flouble *r_arr=compute_shear_spacing(par); free(r_arr); } +#endif //_USE_NEW_LENSING get_max_memory(par,test_memory+par->just_do_pred); print_info("\n"); @@ -505,6 +509,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) if(par->do_kappa) par->kmap=hp_shell_alloc(1,par->nside_kappa,par->nside_base,par->n_kappa); +#ifdef _USE_NEW_LENSING if(par->do_shear) { flouble *r_arr=compute_shear_spacing(par); par->smap=hp_shell_alloc(2,par->nside_shear,par->nside_base,par->n_shear); @@ -515,6 +520,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) } free(r_arr); } +#endif //_USE_NEW_LENSING if(par->do_isw) par->pd_map=hp_shell_alloc(1,par->nside_isw,par->nside_base,par->n_isw); @@ -1208,10 +1214,12 @@ void param_colore_free(ParamCoLoRe *par) #endif //_ADD_EXTRA_KAPPA } +#ifdef _USE_NEW_LENSING if(par->do_shear) { if(par->smap!=NULL) hp_shell_free(par->smap); } +#endif //_USE_NEW_LENSING if(par->do_isw) { if(par->pd_map!=NULL) diff --git a/src/main.c b/src/main.c index 59ce82f..8a01800 100644 --- a/src/main.c +++ b/src/main.c @@ -74,8 +74,10 @@ int main(int argc,char **argv) //Get information from slabs if(par->do_kappa) kappa_set_cartesian(par); +#ifndef _USE_NEW_LENSING if(par->do_shear) shear_set_cartesian(par); +#endif //_USE_NEW_LENSING if(par->do_isw) isw_set_cartesian(par); if(par->do_srcs) @@ -86,8 +88,10 @@ int main(int argc,char **argv) //Distribute information across if(par->do_kappa) kappa_distribute(par); +#ifndef _USE_NEW_LENSING if(par->do_shear) shear_distribute(par); +#endif //_USE_NEW_LENSING if(par->do_isw) isw_distribute(par); if(par->do_srcs) @@ -98,8 +102,10 @@ int main(int argc,char **argv) //Postprocess after if(par->do_kappa) kappa_get_local_properties(par); +#ifndef _USE_NEW_LENSING if(par->do_shear) shear_get_local_properties(par); +#endif //_USE_NEW_LENSING if(par->do_isw) isw_get_local_properties(par); if(par->do_srcs) @@ -115,8 +121,10 @@ int main(int argc,char **argv) //Write output if(par->do_kappa) write_kappa(par); +#ifndef _USE_NEW_LENSING if(par->do_shear) write_shear(par); +#endif ///_USE_NEW_LENSING if(par->do_isw) write_isw(par); if(par->do_srcs) From 672a0729ae85df1415a233e74203befa1d4b2538 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Sat, 4 Apr 2020 08:56:04 +0100 Subject: [PATCH 14/29] bugs, write_shear --- examples/shear_test/param.cfg | 1 + src/common.h | 3 ++- src/io.c | 3 ++- src/main.c | 10 +++++----- 4 files changed, 10 insertions(+), 7 deletions(-) diff --git a/examples/shear_test/param.cfg b/examples/shear_test/param.cfg index fef20fa..fd3e34f 100644 --- a/examples/shear_test/param.cfg +++ b/examples/shear_test/param.cfg @@ -57,4 +57,5 @@ shear: n_shear= 11 # Select spacing type ("r" or "log(1+z)") spacing_type= "r" + write=false } diff --git a/src/common.h b/src/common.h index fb08619..c66fb50 100644 --- a/src/common.h +++ b/src/common.h @@ -322,7 +322,8 @@ typedef struct { flouble **cl_extra_kappa; #endif //_ADD_EXTRA_KAPPA // Shear - int do_shear; //Do you want to output shear maps? + int do_shear; //Do you want to create shear maps? + int write_shear; //Do you want to output shear maps? int shear_spacing_type; //log(1+z)? r? int n_shear; //How many maps? int nside_shear; diff --git a/src/io.c b/src/io.c index 925aed2..3e9d756 100644 --- a/src/io.c +++ b/src/io.c @@ -109,6 +109,7 @@ static ParamCoLoRe *param_colore_new(void) par->nside_kappa=-1; par->nside_shear=-1; par->nside_isw=-1; + par->write_shear=0; for(ii=0;iifnameBzSrcs[ii],"default"); sprintf(par->fnameNzSrcs[ii],"default"); @@ -360,7 +361,6 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) cset=config_lookup(conf,"kappa"); if(cset!=NULL) { par->do_kappa=1; - par->do_srcs_shear=1; conf_read_double_array(conf,"kappa","z_out",par->z_kappa_out,&(par->n_kappa),NPLANES_MAX); conf_read_int(conf,"kappa","nside",&(par->nside_kappa)); } @@ -381,6 +381,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) else report_error(1,"Unknown spacing type %s\n",spacing_string); conf_read_int(conf,"shear","nside",&(par->nside_shear)); + conf_read_bool(conf,"shear","write",&(par->write_shear)); } // Check shear exists if requested with catalog if(par->do_srcs_shear && !(par->do_shear)) diff --git a/src/main.c b/src/main.c index 8a01800..9c6d805 100644 --- a/src/main.c +++ b/src/main.c @@ -74,7 +74,7 @@ int main(int argc,char **argv) //Get information from slabs if(par->do_kappa) kappa_set_cartesian(par); -#ifndef _USE_NEW_LENSING +#ifdef _USE_NEW_LENSING if(par->do_shear) shear_set_cartesian(par); #endif //_USE_NEW_LENSING @@ -88,7 +88,7 @@ int main(int argc,char **argv) //Distribute information across if(par->do_kappa) kappa_distribute(par); -#ifndef _USE_NEW_LENSING +#ifdef _USE_NEW_LENSING if(par->do_shear) shear_distribute(par); #endif //_USE_NEW_LENSING @@ -102,7 +102,7 @@ int main(int argc,char **argv) //Postprocess after if(par->do_kappa) kappa_get_local_properties(par); -#ifndef _USE_NEW_LENSING +#ifdef _USE_NEW_LENSING if(par->do_shear) shear_get_local_properties(par); #endif //_USE_NEW_LENSING @@ -121,8 +121,8 @@ int main(int argc,char **argv) //Write output if(par->do_kappa) write_kappa(par); -#ifndef _USE_NEW_LENSING - if(par->do_shear) +#ifdef _USE_NEW_LENSING + if(par->do_shear && par->write_shear) write_shear(par); #endif ///_USE_NEW_LENSING if(par->do_isw) From 1cd17098bcbea749bce605c4030ace7ebf9ad61d Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 20 May 2020 10:26:11 +0100 Subject: [PATCH 15/29] debugging --- Makefile | 2 +- examples/LSST/param.cfg | 16 ++++++++++++---- src/main.c | 2 +- 3 files changed, 14 insertions(+), 6 deletions(-) diff --git a/Makefile b/Makefile index 59221e2..90159c6 100644 --- a/Makefile +++ b/Makefile @@ -21,7 +21,7 @@ USE_SINGLE_PRECISION = yes #Add random perturbations to kappa from redshifts outside the box ADD_EXTRA_KAPPA = yes #Compile with HDF5 capability? Set to "yes" or "no" -USE_HDF5 = yes +USE_HDF5 = no #Use OMP parallelization? Set to "yes" or "no" USE_OMP = yes #Use MPI parallelization? Set to "yes" or "no" diff --git a/examples/LSST/param.cfg b/examples/LSST/param.cfg index 278acc0..cdf5b22 100644 --- a/examples/LSST/param.cfg +++ b/examples/LSST/param.cfg @@ -16,7 +16,7 @@ field_par: { r_smooth= 2.0 smooth_potential= true - n_grid= 4096 + n_grid= 512 dens_type= 0 lpt_buffer_fraction= 0.6 lpt_interp_type= 1 @@ -36,7 +36,7 @@ cosmo_par: srcs1: { - nz_filename= "examples/LSST/NzRed.txt" + nz_filename= "examples/LSST/NzRed_use.txt" bias_filename= "examples/LSST/BzRed.txt" include_shear= false store_skewers= false @@ -44,7 +44,7 @@ srcs1: srcs2: { - nz_filename= "examples/LSST/NzBlue.txt" + nz_filename= "examples/LSST/NzBlue_use.txt" bias_filename= "examples/LSST/BzBlue.txt" include_shear= true store_skewers= false @@ -52,6 +52,14 @@ srcs2: kappa: { - z_out= [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6] + z_out= [0.1] nside= 512 } + +shear: +{ + nside= 512 + n_shear= 11 + spacing_type= "r" + write=false +} \ No newline at end of file diff --git a/src/main.c b/src/main.c index 9c6d805..d4400f3 100644 --- a/src/main.c +++ b/src/main.c @@ -117,7 +117,7 @@ int main(int argc,char **argv) //and computation of all required quantities if(par->need_beaming) get_beam_properties(par); - + //Write output if(par->do_kappa) write_kappa(par); From b04e219d4f3164835915f24b1cf494288e5f0cb6 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 20 May 2020 12:47:35 +0100 Subject: [PATCH 16/29] mk_new_nz --- examples/LSST/mk_new_nz.py | 8 ++++++++ 1 file changed, 8 insertions(+) create mode 100644 examples/LSST/mk_new_nz.py diff --git a/examples/LSST/mk_new_nz.py b/examples/LSST/mk_new_nz.py new file mode 100644 index 0000000..84321f0 --- /dev/null +++ b/examples/LSST/mk_new_nz.py @@ -0,0 +1,8 @@ +import numpy as np +import sys + +factor = float(sys.argv[1]) +z, nz = np.loadtxt("NzBlue.txt", unpack=True) +np.savetxt("NzBlue_use.txt", np.transpose([z, factor * nz])) +z, nz = np.loadtxt("NzRed.txt", unpack=True) +np.savetxt("NzRed_use.txt", np.transpose([z, factor * nz])) From 83e50d116e49666066e7fd58bc970a3b8f973559 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 20 May 2020 19:40:12 +0100 Subject: [PATCH 17/29] fixed floating point bug --- Makefile | 6 +++--- src/srcs.c | 26 ++++++++++++++------------ 2 files changed, 17 insertions(+), 15 deletions(-) diff --git a/Makefile b/Makefile index 90159c6..38f93ba 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ ###Compiler and compilation options COMP_SER = gcc COMP_MPI = mpicc -OPTIONS = -Wall -O3 -std=c99 +OPTIONS = -Wall -Wno-format-overflow -O3 -std=c99 # ### Behavioural flags #Use double precision integer (enable in general) @@ -33,8 +33,8 @@ USE_MPI = yes #GSL #GSL_INC = -I/add/path #GSL_LIB = -L/add/path -GSL_INC = -I/home/alonso/include -GSL_LIB = -L/home/alonso/lib +GSL_INC = -I/users/damonge/include +GSL_LIB = -L/users/damonge/lib #FFTW FFTW_INC = FFTW_LIB = diff --git a/src/srcs.c b/src/srcs.c index 2e5c558..42eed46 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -250,17 +250,17 @@ static void srcs_set_cartesian_single(ParamCoLoRe *par,int ipop) double rvel=factor_vel*get_rvel(par,ix,iy,iz,x0,y0,z0,rr); double dz_rsd=rvel*get_bg(par,rr,BG_V1,0); for(ip=0;ipcats_c[ipop]->pos[NPOS_CC*pid+0]=pos[0]; - par->cats_c[ipop]->pos[NPOS_CC*pid+1]=pos[1]; - par->cats_c[ipop]->pos[NPOS_CC*pid+2]=pos[2]; + par->cats_c[ipop]->pos[NPOS_CC*pid+0]=x0+dx*(rng_01(rng_thr)-0.5); + par->cats_c[ipop]->pos[NPOS_CC*pid+1]=y0+dx*(rng_01(rng_thr)-0.5); + par->cats_c[ipop]->pos[NPOS_CC*pid+2]=z0+dx*(rng_01(rng_thr)-0.5); par->cats_c[ipop]->pos[NPOS_CC*pid+3]=dz_rsd; + for(ax=0;ax<3;ax++) + pos[ax]=par->cats_c[ipop]->pos[NPOS_CC*pid+ax]; vec2pix_ring(par->nside_base,pos,&pix_id_ring); ring2nest(par->nside_base,pix_id_ring,&pix_id_nest); @@ -506,12 +506,12 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) double rm=(i_r+0.5)*cat->dr; for(ax=0;ax<3;ax++) xn[ax]=(rm*u[ax]+par->pos_obs[ax])*idx; - if(cat->skw_gauss) { - added=interpolate_from_grid(par,xn,&dens,v,NULL,NULL,&gauss,RETURN_GAUSS | RETURN_VEL,INTERP_TYPE_SKW); - } - else { - added=interpolate_from_grid(par,xn,&dens,v,NULL,NULL,&gauss,RETURN_DENS | RETURN_VEL,INTERP_TYPE_SKW); - } + if(cat->skw_gauss) { + added=interpolate_from_grid(par,xn,&dens,v,NULL,NULL,&gauss,RETURN_GAUSS | RETURN_VEL,INTERP_TYPE_SKW); + } + else { + added=interpolate_from_grid(par,xn,&dens,v,NULL,NULL,&gauss,RETURN_DENS | RETURN_VEL,INTERP_TYPE_SKW); + } if(added) { vr=0.5*idx*(v[0]*u[0]+v[1]*u[1]+v[2]*u[2]); if(cat->skw_gauss) { @@ -655,6 +655,8 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) vec2pix_nest(smap->nside,u,&ipix); //Offset to pixel indices in this node ipix-=(ibase-ibase_here)*nside_ratio*nside_ratio; + if((ipix<0) || (ipix>=smap->num_pix)) + report_error(1,"Bad pixel!!\n"); //Find shear at edges g1_0=smap->data[2*(ir_s*smap->num_pix+ipix)+0]; From 7b2618c0a31b881a98f59640e3f52967686cdc8b Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 20 May 2020 19:42:51 +0100 Subject: [PATCH 18/29] spurious changes --- examples/LSST/mk_new_nz.py | 8 -------- examples/LSST/param.cfg | 16 ++++------------ 2 files changed, 4 insertions(+), 20 deletions(-) delete mode 100644 examples/LSST/mk_new_nz.py diff --git a/examples/LSST/mk_new_nz.py b/examples/LSST/mk_new_nz.py deleted file mode 100644 index 84321f0..0000000 --- a/examples/LSST/mk_new_nz.py +++ /dev/null @@ -1,8 +0,0 @@ -import numpy as np -import sys - -factor = float(sys.argv[1]) -z, nz = np.loadtxt("NzBlue.txt", unpack=True) -np.savetxt("NzBlue_use.txt", np.transpose([z, factor * nz])) -z, nz = np.loadtxt("NzRed.txt", unpack=True) -np.savetxt("NzRed_use.txt", np.transpose([z, factor * nz])) diff --git a/examples/LSST/param.cfg b/examples/LSST/param.cfg index cdf5b22..278acc0 100644 --- a/examples/LSST/param.cfg +++ b/examples/LSST/param.cfg @@ -16,7 +16,7 @@ field_par: { r_smooth= 2.0 smooth_potential= true - n_grid= 512 + n_grid= 4096 dens_type= 0 lpt_buffer_fraction= 0.6 lpt_interp_type= 1 @@ -36,7 +36,7 @@ cosmo_par: srcs1: { - nz_filename= "examples/LSST/NzRed_use.txt" + nz_filename= "examples/LSST/NzRed.txt" bias_filename= "examples/LSST/BzRed.txt" include_shear= false store_skewers= false @@ -44,7 +44,7 @@ srcs1: srcs2: { - nz_filename= "examples/LSST/NzBlue_use.txt" + nz_filename= "examples/LSST/NzBlue.txt" bias_filename= "examples/LSST/BzBlue.txt" include_shear= true store_skewers= false @@ -52,14 +52,6 @@ srcs2: kappa: { - z_out= [0.1] + z_out= [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6] nside= 512 } - -shear: -{ - nside= 512 - n_shear= 11 - spacing_type= "r" - write=false -} \ No newline at end of file From 21ca9033422b882faf9b2a3c652ca74dc2b9565d Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 17 Jun 2020 18:26:00 +0100 Subject: [PATCH 19/29] adaptive creator written --- Makefile | 4 ++-- src/common.c | 62 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/common.h | 14 ++++++++++++ src/io.c | 3 +++ 4 files changed, 81 insertions(+), 2 deletions(-) diff --git a/Makefile b/Makefile index 38f93ba..39a785f 100644 --- a/Makefile +++ b/Makefile @@ -33,8 +33,8 @@ USE_MPI = yes #GSL #GSL_INC = -I/add/path #GSL_LIB = -L/add/path -GSL_INC = -I/users/damonge/include -GSL_LIB = -L/users/damonge/lib +GSL_INC = -I/home/alonso/include +GSL_LIB = -L/home/alonso/lib #FFTW FFTW_INC = FFTW_LIB = diff --git a/src/common.c b/src/common.c index 0b58d10..ea36fb2 100644 --- a/src/common.c +++ b/src/common.c @@ -569,6 +569,68 @@ void catalog_free(Catalog *cat) free(cat); } +HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_base,int nr, flouble *r_arr, flouble dx, flouble dx_fraction) +{ + if(nside_max>NSIDE_MAX_HPX) + report_error(1,"Can't go beyond nside=%d\n",NSIDE_MAX_HPX); + if(nside_maxnbeams=0; + for(ib=NodeThis;ibnbeams++; + + shell->nr=nr; + shell->nside=my_malloc(nr*sizeof(int)); + shell->nside_ratio=my_malloc(nr*sizeof(int)); + shell->num_pix_per_beam=my_malloc(nr*sizeof(long)); + shell->r=my_malloc(nr*sizeof(flouble)); + shell->ipix_0=my_malloc(nr*sizeof(long *)); + shell->beam_id=my_malloc(shell->nbeams*sizeof(long)); + for(ib=0;ibnbeams;ib++) + shell->beam_id[ib]=NodeThis+NNodes*ib; + for(ir=0; irr[ir]=r; + shell->nside[ir]=nside_here; + nside_ratio=nside_here/nside_base; + shell->nside_ratio[ir]=nside_ratio; + shell->num_pix_per_beam[ir]=nside_ratio*nside_ratio; + printf("%d %d %d %ld %lE %lE %lE",ir, shell->nside[ir], + shell->nside_ratio[ir], shell->num_pix_per_beam[ir], + shell->r[ir], + sqrt(4*M_PI/he_nside2npix(nside_here))*r,dx); + printf(" ["); + shell->ipix_0[ir]=my_malloc(shell->nbeams*sizeof(long)); + for(ib=0;ibnbeams;ib++) { + shell->ipix_0[ir][ib]=shell->beam_id[ib]*nside_ratio*nside_ratio; + printf(" %ld", shell->ipix_0[ir][ib]); + } + printf("]\n"); + } + + shell->data=my_malloc(nr*sizeof(flouble **)); + for(ir=0; irdata[ir]=my_malloc(shell->nbeams*sizeof(flouble *)); + for(ib=0;ibnbeams;ib++) + shell->data[ir][ib]=my_malloc(shell->nq*shell->num_pix_per_beam[ir]*sizeof(flouble)); + } + exit(1); +} + HealpixShells *hp_shell_alloc(int nq, int nside,int nside_base,int nr) { if(nside>NSIDE_MAX_HPX) diff --git a/src/common.h b/src/common.h index c66fb50..ad71728 100644 --- a/src/common.h +++ b/src/common.h @@ -187,6 +187,19 @@ typedef struct { float *v_skw; } Catalog; +typedef struct { + int nq; + int nr; //Number of spherical shells + flouble *r; //r in this shell + int *nside; //Resolution parameter + int *nside_ratio; + long *num_pix_per_beam; + int nbeams; //Number of different bases + long *beam_id; + long **ipix_0; + flouble ***data; +} HealpixShellsAdaptive; + typedef struct { int nq; //Number of quantities to map int nside; //Resolution parameter @@ -373,6 +386,7 @@ unsigned long long get_max_memory(ParamCoLoRe *par,int just_test); void get_radial_params(double rmax,int ngrid,int *nr,double *dr); //void get_random_angles(gsl_rng *rng,int ipix_nest,int ipix0,int nside,double *th,double *phi); HealpixShells *hp_shell_alloc(int nq,int nside,int nside_base,int nr); +HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_base,int nr, flouble *r_arr, flouble dx, flouble dx_fraction); void hp_shell_free(HealpixShells *shell); CatalogCartesian *catalog_cartesian_alloc(int nsrcs); void catalog_cartesian_free(CatalogCartesian *cat); diff --git a/src/io.c b/src/io.c index 3e9d756..be29902 100644 --- a/src/io.c +++ b/src/io.c @@ -513,6 +513,9 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) #ifdef _USE_NEW_LENSING if(par->do_shear) { flouble *r_arr=compute_shear_spacing(par); + HealpixShellsAdaptive *ss=hp_shell_adaptive_alloc(2, 4096, par->nside_base, + par->n_shear, r_arr, par->l_box/par->n_grid, + 1.); par->smap=hp_shell_alloc(2,par->nside_shear,par->nside_base,par->n_shear); int i_r; for(i_r=0;i_rn_shear;i_r++) { From cd091c91b41f1256a62ba7a155ef418777b92a88 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 17 Jun 2020 18:27:32 +0100 Subject: [PATCH 20/29] adaptive creator written --- src/common.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/common.c b/src/common.c index ea36fb2..7ed886b 100644 --- a/src/common.c +++ b/src/common.c @@ -628,7 +628,8 @@ HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_ for(ib=0;ibnbeams;ib++) shell->data[ir][ib]=my_malloc(shell->nq*shell->num_pix_per_beam[ir]*sizeof(flouble)); } - exit(1); + + return shell; } HealpixShells *hp_shell_alloc(int nq, int nside,int nside_base,int nr) From 476ab0d0424a5f93225d3bae7ff6787dddd0b49c Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 17 Jun 2020 18:41:03 +0100 Subject: [PATCH 21/29] small bug in creator and destructor --- src/common.c | 21 +++++++++++++++++++++ src/common.h | 1 + src/io.c | 2 ++ 3 files changed, 24 insertions(+) diff --git a/src/common.c b/src/common.c index 7ed886b..8ced4d3 100644 --- a/src/common.c +++ b/src/common.c @@ -569,6 +569,26 @@ void catalog_free(Catalog *cat) free(cat); } +void hp_shell_adaptive_free(HealpixShellsAdaptive *shell) +{ + int ir; + free(shell->r); + free(shell->nside); + free(shell->nside_ratio); + free(shell->num_pix_per_beam); + free(shell->beam_id); + for(ir=0;irnr;ir++) { + int ib; + free(shell->ipix_0[ir]); + for(ib=0;ibnbeams;ib++) + free(shell->data[ir][ib]); + free(shell->data[ir]); + } + free(shell->ipix_0); + free(shell->data); + free(shell); +} + HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_base,int nr, flouble *r_arr, flouble dx, flouble dx_fraction) { if(nside_max>NSIDE_MAX_HPX) @@ -584,6 +604,7 @@ HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_ for(ib=NodeThis;ibnbeams++; + shell->nq=nq; shell->nr=nr; shell->nside=my_malloc(nr*sizeof(int)); shell->nside_ratio=my_malloc(nr*sizeof(int)); diff --git a/src/common.h b/src/common.h index ad71728..ffef6e0 100644 --- a/src/common.h +++ b/src/common.h @@ -388,6 +388,7 @@ void get_radial_params(double rmax,int ngrid,int *nr,double *dr); HealpixShells *hp_shell_alloc(int nq,int nside,int nside_base,int nr); HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_base,int nr, flouble *r_arr, flouble dx, flouble dx_fraction); void hp_shell_free(HealpixShells *shell); +void hp_shell_adaptive_free(HealpixShellsAdaptive *shell); CatalogCartesian *catalog_cartesian_alloc(int nsrcs); void catalog_cartesian_free(CatalogCartesian *cat); Catalog *catalog_alloc(int nsrcs,int has_shear,int has_skw,int skw_gauss,double rmax,int ng); diff --git a/src/io.c b/src/io.c index be29902..8862238 100644 --- a/src/io.c +++ b/src/io.c @@ -516,6 +516,8 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) HealpixShellsAdaptive *ss=hp_shell_adaptive_alloc(2, 4096, par->nside_base, par->n_shear, r_arr, par->l_box/par->n_grid, 1.); + hp_shell_adaptive_free(ss); + exit(1); par->smap=hp_shell_alloc(2,par->nside_shear,par->nside_base,par->n_shear); int i_r; for(i_r=0;i_rn_shear;i_r++) { From 4716d733c83607efa35032a33181dc5f62c115d9 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 17 Jun 2020 18:42:26 +0100 Subject: [PATCH 22/29] small bug in creator and destructor --- src/common.h | 2 +- src/io.c | 14 ++++---------- 2 files changed, 5 insertions(+), 11 deletions(-) diff --git a/src/common.h b/src/common.h index ffef6e0..83ad6a8 100644 --- a/src/common.h +++ b/src/common.h @@ -340,7 +340,7 @@ typedef struct { int shear_spacing_type; //log(1+z)? r? int n_shear; //How many maps? int nside_shear; - HealpixShells *smap; //Shear maps at each redshift + HealpixShellsAdaptive *smap; //Shear maps at each redshift // ISW int do_isw; //Do you want to output isw maps? int n_isw; //How many maps? diff --git a/src/io.c b/src/io.c index 8862238..fbd7e57 100644 --- a/src/io.c +++ b/src/io.c @@ -513,17 +513,11 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) #ifdef _USE_NEW_LENSING if(par->do_shear) { flouble *r_arr=compute_shear_spacing(par); - HealpixShellsAdaptive *ss=hp_shell_adaptive_alloc(2, 4096, par->nside_base, - par->n_shear, r_arr, par->l_box/par->n_grid, - 1.); - hp_shell_adaptive_free(ss); + par->smap = hp_shell_adaptive_alloc(2, 4096, par->nside_base, + par->n_shear, r_arr, par->l_box/par->n_grid, + 1.); + hp_shell_adaptive_free(par->smap); exit(1); - par->smap=hp_shell_alloc(2,par->nside_shear,par->nside_base,par->n_shear); - int i_r; - for(i_r=0;i_rn_shear;i_r++) { - par->smap->r0[i_r]=r_arr[i_r]; - par->smap->rf[i_r]=r_arr[i_r]; - } free(r_arr); } #endif //_USE_NEW_LENSING From b7486730ec42309d57f0c3ebc4d7e842fc32a426 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Wed, 17 Jun 2020 19:11:10 +0100 Subject: [PATCH 23/29] writer --- src/io.c | 39 ++++++++++++--------------------------- 1 file changed, 12 insertions(+), 27 deletions(-) diff --git a/src/io.c b/src/io.c index fbd7e57..8285017 100644 --- a/src/io.c +++ b/src/io.c @@ -513,7 +513,7 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) #ifdef _USE_NEW_LENSING if(par->do_shear) { flouble *r_arr=compute_shear_spacing(par); - par->smap = hp_shell_adaptive_alloc(2, 4096, par->nside_base, + par->smap = hp_shell_adaptive_alloc(2, par->nside_shear, par->nside_base, par->n_shear, r_arr, par->l_box/par->n_grid, 1.); hp_shell_adaptive_free(par->smap); @@ -800,27 +800,25 @@ void write_shear(ParamCoLoRe *par) char fname[256]; long npx=he_nside2npix(par->nside_shear); flouble **map_write=my_malloc(2*sizeof(flouble *)); - int *map_nadd=my_malloc(npx*sizeof(int)); map_write[0]=my_malloc(npx*sizeof(flouble)); map_write[1]=my_malloc(npx*sizeof(flouble)); if(NodeThis==0) timer(0); print_info("*** Writing shear source maps\n"); for(ir=0;irsmap->nr;ir++) { - long ip; - long ir_t=ir*par->smap->num_pix; + long ip, ib; //Write local pixels to dummy map for(ip=0;ipprefixOut,ir); - for(ip=0;ipsmap->num_pix;ip++) { - int id_pix=par->smap->listpix[ip]; - map_write[0][id_pix]+=par->smap->data[0 + 2*(ip + ir_t)]; - map_write[1][id_pix]+=par->smap->data[1 + 2*(ip + ir_t)]; - map_nadd[id_pix]+=par->smap->nadd[ir_t+ip]; + for(ib=0;ibsmap->nbeams;ib++) { + for(ip=0;ipsmap->num_pix_per_beam[ir];ip++) { + int id_pix=par->smap->ipix_0[ir][ib]+ip; + map_write[0][id_pix]+=par->smap->data[ir][ib][2*ip+0]; + map_write[1][id_pix]+=par->smap->data[ir][ib][2*ip+1]; + } } //Collect all dummy maps @@ -833,19 +831,8 @@ void write_shear(ParamCoLoRe *par) MPI_Reduce(MPI_IN_PLACE,map_write[1],npx,FLOUBLE_MPI,MPI_SUM,0,MPI_COMM_WORLD); else MPI_Reduce(map_write[1],NULL ,npx,FLOUBLE_MPI,MPI_SUM,0,MPI_COMM_WORLD); - if(NodeThis==0) - MPI_Reduce(MPI_IN_PLACE,map_nadd,npx,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD); - else - MPI_Reduce(map_nadd ,NULL ,npx,MPI_INT,MPI_SUM,0,MPI_COMM_WORLD); #endif //_HAVE_MPI - for(ip=0;ip0) { - map_write[0][ip]/=map_nadd[ip]; - map_write[1][ip]/=map_nadd[ip]; - } - } - //Write dummy map if(NodeThis==0) he_write_healpix_map(map_write,2,par->nside_shear,fname,1); @@ -853,7 +840,6 @@ void write_shear(ParamCoLoRe *par) free(map_write[0]); free(map_write[1]); free(map_write); - free(map_nadd); if(NodeThis==0) { sprintf(fname,"%s_shear_r.txt", par->prefixOut); @@ -861,10 +847,9 @@ void write_shear(ParamCoLoRe *par) if(f==NULL) report_error(1,"Couldn't open file %s\n", fname); for(ir=0;irsmap->nr;ir++) { - fprintf(f,"%d %lE %lE %lE %lE\n",ir, - par->smap->r0[ir], par->smap->rf[ir], - get_bg(par, par->smap->r0[ir], BG_Z, 0), - get_bg(par, par->smap->rf[ir], BG_Z, 0)); + fprintf(f,"%d %lE %lE\n",ir, + par->smap->r[ir], + get_bg(par, par->smap->r[ir], BG_Z, 0)); } fclose(f); } @@ -1217,7 +1202,7 @@ void param_colore_free(ParamCoLoRe *par) #ifdef _USE_NEW_LENSING if(par->do_shear) { if(par->smap!=NULL) - hp_shell_free(par->smap); + hp_shell_adaptive_free(par->smap); } #endif //_USE_NEW_LENSING From 012a51b5b51101022c45b15832a5b46b69060b9d Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 18 Jun 2020 12:48:05 +0100 Subject: [PATCH 24/29] beamer --- src/beaming.c | 5 +- src/common.c | 48 ++++++++------ src/common.h | 1 + src/io.c | 2 - src/main.c | 4 +- src/shear.c | 180 +++++++++++++++++++++++++++++++++++++++++++------- src/srcs.c | 10 +-- 7 files changed, 193 insertions(+), 57 deletions(-) diff --git a/src/beaming.c b/src/beaming.c index a208c75..887fff6 100644 --- a/src/beaming.c +++ b/src/beaming.c @@ -337,7 +337,7 @@ void get_beam_properties(ParamCoLoRe *par) kappa_get_beam_properties(par); #ifdef _USE_NEW_LENSING if(par->do_shear) - shear_get_beam_properties(par); + shear_get_beam_properties(par); #endif //_USE_NEW_LENSING if(par->do_isw) isw_get_beam_properties(par); @@ -349,7 +349,8 @@ void get_beam_properties(ParamCoLoRe *par) #ifdef _HAVE_MPI free(buffer_sr); #endif //_HAVE_MPI - + exit(1); + if(par->do_kappa) kappa_beams_postproc(par); #ifdef _USE_NEW_LENSING diff --git a/src/common.c b/src/common.c index 8ced4d3..c3fee1d 100644 --- a/src/common.c +++ b/src/common.c @@ -571,21 +571,23 @@ void catalog_free(Catalog *cat) void hp_shell_adaptive_free(HealpixShellsAdaptive *shell) { - int ir; + int ib,ir; free(shell->r); free(shell->nside); free(shell->nside_ratio); free(shell->num_pix_per_beam); free(shell->beam_id); - for(ir=0;irnr;ir++) { - int ib; - free(shell->ipix_0[ir]); - for(ib=0;ibnbeams;ib++) - free(shell->data[ir][ib]); - free(shell->data[ir]); + for(ib=0;ibnbeams;ib++) { + for(ir=0;irnr;ir++) + free(shell->data[ib][ir]); + free(shell->data[ib]); + free(shell->pos[ib]); } + for(ir=0;irnr;ir++) + free(shell->ipix_0[ir]); free(shell->ipix_0); free(shell->data); + free(shell->pos); free(shell); } @@ -630,26 +632,30 @@ HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_ nside_ratio=nside_here/nside_base; shell->nside_ratio[ir]=nside_ratio; shell->num_pix_per_beam[ir]=nside_ratio*nside_ratio; - printf("%d %d %d %ld %lE %lE %lE",ir, shell->nside[ir], - shell->nside_ratio[ir], shell->num_pix_per_beam[ir], - shell->r[ir], - sqrt(4*M_PI/he_nside2npix(nside_here))*r,dx); - printf(" ["); shell->ipix_0[ir]=my_malloc(shell->nbeams*sizeof(long)); - for(ib=0;ibnbeams;ib++) { + for(ib=0;ibnbeams;ib++) shell->ipix_0[ir][ib]=shell->beam_id[ib]*nside_ratio*nside_ratio; - printf(" %ld", shell->ipix_0[ir][ib]); - } - printf("]\n"); } - shell->data=my_malloc(nr*sizeof(flouble **)); - for(ir=0; irdata[ir]=my_malloc(shell->nbeams*sizeof(flouble *)); - for(ib=0;ibnbeams;ib++) - shell->data[ir][ib]=my_malloc(shell->nq*shell->num_pix_per_beam[ir]*sizeof(flouble)); + long npix_hi=shell->num_pix_per_beam[nr-1]; + shell->data=my_malloc(shell->nbeams*sizeof(flouble **)); + shell->pos=my_malloc(shell->nbeams*sizeof(double **)); + for(ib=0;ibnbeams;ib++) { + shell->data[ib]=my_malloc(nr*sizeof(flouble *)); + shell->pos[ib]=my_malloc(npix_hi*3*sizeof(double)); + for(ir=0; irdata[ib][ir]=my_malloc(shell->nq*shell->num_pix_per_beam[ir]*sizeof(flouble)); } + for(ib=0;ibnbeams;ib++) { + long ip, ip0=(ib*NNodes+NodeThis)*npix_hi; + double *u=shell->pos[ib]; + for(ip=0;ipnside[nr-1],id_nest,u); + u+=3; + } + } return shell; } diff --git a/src/common.h b/src/common.h index 83ad6a8..12cd4b8 100644 --- a/src/common.h +++ b/src/common.h @@ -197,6 +197,7 @@ typedef struct { int nbeams; //Number of different bases long *beam_id; long **ipix_0; + double **pos; flouble ***data; } HealpixShellsAdaptive; diff --git a/src/io.c b/src/io.c index 8285017..6e5fead 100644 --- a/src/io.c +++ b/src/io.c @@ -516,8 +516,6 @@ ParamCoLoRe *read_run_params(char *fname,int test_memory) par->smap = hp_shell_adaptive_alloc(2, par->nside_shear, par->nside_base, par->n_shear, r_arr, par->l_box/par->n_grid, 1.); - hp_shell_adaptive_free(par->smap); - exit(1); free(r_arr); } #endif //_USE_NEW_LENSING diff --git a/src/main.c b/src/main.c index d4400f3..f5cff83 100644 --- a/src/main.c +++ b/src/main.c @@ -70,7 +70,7 @@ int main(int argc,char **argv) //Compute normalization of density field for biasing compute_density_normalization(par); - + //Get information from slabs if(par->do_kappa) kappa_set_cartesian(par); @@ -112,7 +112,7 @@ int main(int argc,char **argv) srcs_get_local_properties(par); if(par->do_imap) imap_get_local_properties(par); - + //All-to-all communication of density field //and computation of all required quantities if(par->need_beaming) diff --git a/src/shear.c b/src/shear.c index 254131a..78c93eb 100644 --- a/src/shear.c +++ b/src/shear.c @@ -40,19 +40,15 @@ void shear_beams_preproc(ParamCoLoRe *par) { //Sort radii in ascending order int ir; - flouble *r0_i=my_malloc(par->smap->nr*sizeof(flouble)); - flouble *rf_i=my_malloc(par->smap->nr*sizeof(flouble)); + flouble *r_i=my_malloc(par->smap->nr*sizeof(flouble)); - memcpy(r0_i,par->smap->r0,par->smap->nr*sizeof(flouble)); - memcpy(rf_i,par->smap->rf,par->smap->nr*sizeof(flouble)); + memcpy(r_i,par->smap->r,par->smap->nr*sizeof(flouble)); - int *i_sorted=ind_sort(par->smap->nr,par->smap->r0); + int *i_sorted=ind_sort(par->smap->nr,par->smap->r); for(ir=0;irsmap->nr;ir++) { - par->smap->r0[ir]=r0_i[i_sorted[ir]]; - par->smap->rf[ir]=rf_i[i_sorted[ir]]; + par->smap->r[ir]=r_i[i_sorted[ir]]; } - free(r0_i); - free(rf_i); + free(r_i); free(i_sorted); //Zero all data and clear @@ -61,29 +57,26 @@ void shear_beams_preproc(ParamCoLoRe *par) shared(par) #endif //_HAVE_OMP { - long ipp; + long ir, ib, ipp; #ifdef _HAVE_OMP #pragma omp for #endif //_HAVE_OMP - for(ipp=0;ipp<2*par->smap->num_pix*par->smap->nr;ipp++) { - par->smap->data[ipp]=0; - } //end omp for - -#ifdef _HAVE_OMP -#pragma omp for -#endif //_HAVE_OMP - for(ipp=0;ippsmap->num_pix*par->smap->nr;ipp++) { - par->smap->nadd[ipp]=1; + for(ib=0;ibsmap->nbeams;ib++) { + for(ir=0;irsmap->nr;ir++) { + for(ipp=0;ipp<2*par->smap->num_pix_per_beam[ir];ipp++) + par->smap->data[ib][ir][ipp]=0; + } } //end omp for } //end omp parallel return; } +/* void shear_get_beam_properties(ParamCoLoRe *par) { - HealpixShells *smap=par->smap; + HealpixShellsAdaptive *smap=par->smap; #ifdef _HAVE_OMP #pragma omp parallel default(none) \ @@ -103,7 +96,7 @@ void shear_get_beam_properties(ParamCoLoRe *par) int *i_r_min_arr=my_malloc(smap->nr*sizeof(int)); double *inv_r_max=my_malloc(smap->nr*sizeof(double)); for(i_r=0;i_rnr;i_r++) { - int i_r_here=(int)(smap->rf[i_r]*idr+0.5); + int i_r_here=(int)(smap->r[i_r]*idr+0.5); inv_r_max[i_r]=1./(i_r_here*dr); i_r_max_arr[i_r]=MIN(i_r_here,nr-1); } @@ -188,10 +181,145 @@ void shear_get_beam_properties(ParamCoLoRe *par) #pragma omp single #endif //_HAVE_OMP { - for(i_r=0;i_rnr;i_r++) { - smap->r0[i_r]=1./inv_r_max[i_r]; - smap->rf[i_r]=1./inv_r_max[i_r]; - } + for(i_r=0;i_rnr;i_r++) + smap->r[i_r]=1./inv_r_max[i_r]; + } + + free(fac_r_1); + free(fac_r_2); + free(i_r_max_arr); + free(i_r_min_arr); + free(inv_r_max); + } //end omp parallel + + return; +} +*/ + +void shear_get_beam_properties(ParamCoLoRe *par) +{ + HealpixShellsAdaptive *smap=par->smap; + +#ifdef _HAVE_OMP +#pragma omp parallel default(none) \ + shared(par,smap) +#endif //_HAVE_OMP + { + double idx=par->n_grid/par->l_box; + long npix_hi=smap->num_pix_per_beam[smap->nr-1]; + + //Setup radial decomposition + int i_r,nr; + double idr,dr; + get_radial_params(par->r_max,par->n_grid,&nr,&dr); + idr=1./dr; + + //Compute index of each source plane + int *i_r_max_arr=my_malloc(smap->nr*sizeof(int)); + int *i_r_min_arr=my_malloc(smap->nr*sizeof(int)); + long *npix_ratio=my_malloc(smap->nr*sizeof(long)); + double *inv_npix_ratio=my_malloc(smap->nr*sizeof(double)); + double *inv_r_max=my_malloc(smap->nr*sizeof(double)); + for(i_r=0;i_rnr;i_r++) { + int i_r_here=(int)(smap->r[i_r]*idr+0.5); + inv_r_max[i_r]=1./(i_r_here*dr); + i_r_max_arr[i_r]=MIN(i_r_here,nr-1); + npix_ratio[i_r]=npix_hi/smap->num_pix_per_beam[i_r]; + inv_npix_ratio[i_r]=1./((double)(npix_ratio[i_r])); + } + + i_r_min_arr[0]=0; + for(i_r=1;i_rnr;i_r++) + i_r_min_arr[i_r]=i_r_max_arr[i_r-1]+1; + + //Fill up integral kernels + double *fac_r_1=my_malloc(nr*sizeof(double)); + double *fac_r_2=my_malloc(nr*sizeof(double)); + for(i_r=0;i_rnbeams;ib++) { + //We loop over all pixels in the last plane (since this has the + //highest resolution). + long ip; +#ifdef _HAVE_OMP +#pragma omp for +#endif //_HAVE_OMP + for(ip=0;ippos[ib][3*ip]); + double prefac=2*idx*idx; //2/Dx^2 + double cth_h=1,sth_h=0,cph_h=1,sph_h=0; + + cth_h=u[2]; + if(cth_h>=1) cth_h=1; + if(cth_h<=-1) cth_h=-1; + sth_h=sqrt((1-cth_h)*(1+cth_h)); + if(sth_h!=0) { + cph_h=u[0]/sth_h; + sph_h=u[1]/sth_h; + } + + r1[0]=(cth_h*cth_h*cph_h*cph_h-sph_h*sph_h)*prefac; + r1[1]=(2*cph_h*sph_h*(cth_h*cth_h+1))*prefac; + r1[2]=(-2*cth_h*sth_h*cph_h)*prefac; + r1[3]=(cth_h*cth_h*sph_h*sph_h-cph_h*cph_h)*prefac; + r1[4]=(-2*cth_h*sth_h*sph_h)*prefac; + r1[5]=(sth_h*sth_h)*prefac; + + r2[0]=(-2*cth_h*cph_h*sph_h)*prefac; + r2[1]=(2*cth_h*(cph_h*cph_h-sph_h*sph_h))*prefac; + r2[2]=(2*sth_h*sph_h)*prefac; + r2[3]=(2*cth_h*sph_h*cph_h)*prefac; + r2[4]=(-2*sth_h*cph_h)*prefac; + r2[5]=0; + for(i_r=0;i_rnr;i_r++) { + int irr; + int irmin=i_r_min_arr[i_r]; + int irmax=i_r_max_arr[i_r]; + long ipix_this=ip/npix_ratio[i_r]; + for(irr=irmin;irr<=irmax;irr++) { + double rm=(irr+0.5)*dr; + for(ax=0;ax<3;ax++) + xn[ax]=(rm*u[ax]+par->pos_obs[ax])*idx; + added=interpolate_from_grid(par,xn,NULL,NULL,t,NULL,NULL,RETURN_TID,INTERP_TYPE_SHEAR); + if(added) { + double dotp1=0,dotp2=0; + for(ax=0;ax<6;ax++) { + dotp1+=r1[ax]*t[ax]; + dotp2+=r2[ax]*t[ax]; + } + shear1_1+=dotp1*fac_r_1[irr]; + shear1_2+=dotp1*fac_r_2[irr]; + shear2_1+=dotp2*fac_r_1[irr]; + shear2_2+=dotp2*fac_r_2[irr]; + } + } + if(ipix_this>=smap->num_pix_per_beam[i_r]) { + printf("SHIT\n"); + exit(1); + } + smap->data[ib][i_r][2*ipix_this+0]+=(shear1_1-inv_r_max[i_r]*shear1_2)*inv_npix_ratio[i_r]; + smap->data[ib][i_r][2*ipix_this+1]+=(shear2_1-inv_r_max[i_r]*shear2_2)*inv_npix_ratio[i_r]; + } + } //end omp for + } + +#ifdef _HAVE_OMP +#pragma omp single +#endif //_HAVE_OMP + { + for(i_r=0;i_rnr;i_r++) + smap->r[i_r]=1./inv_r_max[i_r]; } free(fac_r_1); @@ -199,6 +327,8 @@ void shear_get_beam_properties(ParamCoLoRe *par) free(i_r_max_arr); free(i_r_min_arr); free(inv_r_max); + free(npix_ratio); + free(inv_npix_ratio); } //end omp parallel return; diff --git a/src/srcs.c b/src/srcs.c index 42eed46..bacc7cf 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -21,7 +21,7 @@ /////////////////////////////////////////////////////////////////////// #include "common.h" -static int get_r_index_smap(HealpixShells *sh,double r,int ir_start) +static int get_r_index_smap(HealpixShellsAdaptive *sh,double r,int ir_start) { int gotit=0; int ir0; @@ -34,13 +34,13 @@ static int get_r_index_smap(HealpixShells *sh,double r,int ir_start) while(!gotit) { if(ir0==0) { - if(rr0[1]) + if(rr[1]) gotit=1; else ir0++; } else if(ir0==sh->nr-1) { - if(r>=sh->rf[sh->nr-2]) { + if(r>=sh->r[sh->nr-2]) { ir0=sh->nr-2; gotit=1; } @@ -48,10 +48,10 @@ static int get_r_index_smap(HealpixShells *sh,double r,int ir_start) ir0--; } else { - if(rr0[ir0]) + if(rr[ir0]) ir0--; else { - if(r>=sh->rf[ir0+1]) + if(r>=sh->r[ir0+1]) ir0++; else gotit=1; From 492e60ef31fdd87f95a3546e932e089132303f21 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 18 Jun 2020 15:04:20 +0100 Subject: [PATCH 25/29] first working version --- src/beaming.c | 1 - src/common.c | 18 +++--------------- src/common.h | 3 --- src/io.c | 11 ++++++----- src/srcs.c | 48 ++++++++++++++++++++++++++++-------------------- 5 files changed, 37 insertions(+), 44 deletions(-) diff --git a/src/beaming.c b/src/beaming.c index 887fff6..6dd0dfc 100644 --- a/src/beaming.c +++ b/src/beaming.c @@ -349,7 +349,6 @@ void get_beam_properties(ParamCoLoRe *par) #ifdef _HAVE_MPI free(buffer_sr); #endif //_HAVE_MPI - exit(1); if(par->do_kappa) kappa_beams_postproc(par); diff --git a/src/common.c b/src/common.c index c3fee1d..311bee6 100644 --- a/src/common.c +++ b/src/common.c @@ -574,24 +574,21 @@ void hp_shell_adaptive_free(HealpixShellsAdaptive *shell) int ib,ir; free(shell->r); free(shell->nside); - free(shell->nside_ratio); free(shell->num_pix_per_beam); - free(shell->beam_id); for(ib=0;ibnbeams;ib++) { for(ir=0;irnr;ir++) free(shell->data[ib][ir]); free(shell->data[ib]); free(shell->pos[ib]); } - for(ir=0;irnr;ir++) - free(shell->ipix_0[ir]); - free(shell->ipix_0); free(shell->data); free(shell->pos); free(shell); } -HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_base,int nr, flouble *r_arr, flouble dx, flouble dx_fraction) +HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_base, + int nr, flouble *r_arr, flouble dx, + flouble dx_fraction) { if(nside_max>NSIDE_MAX_HPX) report_error(1,"Can't go beyond nside=%d\n",NSIDE_MAX_HPX); @@ -609,13 +606,8 @@ HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_ shell->nq=nq; shell->nr=nr; shell->nside=my_malloc(nr*sizeof(int)); - shell->nside_ratio=my_malloc(nr*sizeof(int)); shell->num_pix_per_beam=my_malloc(nr*sizeof(long)); shell->r=my_malloc(nr*sizeof(flouble)); - shell->ipix_0=my_malloc(nr*sizeof(long *)); - shell->beam_id=my_malloc(shell->nbeams*sizeof(long)); - for(ib=0;ibnbeams;ib++) - shell->beam_id[ib]=NodeThis+NNodes*ib; for(ir=0; irr[ir]=r; shell->nside[ir]=nside_here; nside_ratio=nside_here/nside_base; - shell->nside_ratio[ir]=nside_ratio; shell->num_pix_per_beam[ir]=nside_ratio*nside_ratio; - shell->ipix_0[ir]=my_malloc(shell->nbeams*sizeof(long)); - for(ib=0;ibnbeams;ib++) - shell->ipix_0[ir][ib]=shell->beam_id[ib]*nside_ratio*nside_ratio; } long npix_hi=shell->num_pix_per_beam[nr-1]; diff --git a/src/common.h b/src/common.h index 12cd4b8..a92a0e8 100644 --- a/src/common.h +++ b/src/common.h @@ -192,11 +192,8 @@ typedef struct { int nr; //Number of spherical shells flouble *r; //r in this shell int *nside; //Resolution parameter - int *nside_ratio; long *num_pix_per_beam; int nbeams; //Number of different bases - long *beam_id; - long **ipix_0; double **pos; flouble ***data; } HealpixShellsAdaptive; diff --git a/src/io.c b/src/io.c index 6e5fead..9ccf0d8 100644 --- a/src/io.c +++ b/src/io.c @@ -796,7 +796,7 @@ void write_shear(ParamCoLoRe *par) { int ir; char fname[256]; - long npx=he_nside2npix(par->nside_shear); + long npx=he_nside2npix(par->smap->nside[par->smap->nr-1]); flouble **map_write=my_malloc(2*sizeof(flouble *)); map_write[0]=my_malloc(npx*sizeof(flouble)); map_write[1]=my_malloc(npx*sizeof(flouble)); @@ -812,10 +812,11 @@ void write_shear(ParamCoLoRe *par) } sprintf(fname,"!%s_shear_z%03d.fits",par->prefixOut,ir); for(ib=0;ibsmap->nbeams;ib++) { + long ipix_0=(NodeThis+NNodes*ib)*par->smap->num_pix_per_beam[ir]; for(ip=0;ipsmap->num_pix_per_beam[ir];ip++) { - int id_pix=par->smap->ipix_0[ir][ib]+ip; - map_write[0][id_pix]+=par->smap->data[ir][ib][2*ip+0]; - map_write[1][id_pix]+=par->smap->data[ir][ib][2*ip+1]; + long id_pix=ipix_0+ip; + map_write[0][id_pix]+=par->smap->data[ib][ir][2*ip+0]; + map_write[1][id_pix]+=par->smap->data[ib][ir][2*ip+1]; } } @@ -833,7 +834,7 @@ void write_shear(ParamCoLoRe *par) //Write dummy map if(NodeThis==0) - he_write_healpix_map(map_write,2,par->nside_shear,fname,1); + he_write_healpix_map(map_write,2,par->smap->nside[ir],fname,1); } free(map_write[0]); free(map_write[1]); diff --git a/src/srcs.c b/src/srcs.c index bacc7cf..6a98d43 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -21,6 +21,7 @@ /////////////////////////////////////////////////////////////////////// #include "common.h" +#ifdef _USE_NEW_LENSING static int get_r_index_smap(HealpixShellsAdaptive *sh,double r,int ir_start) { int gotit=0; @@ -61,6 +62,7 @@ static int get_r_index_smap(HealpixShellsAdaptive *sh,double r,int ir_start) return ir0; } +#endif //_USE_NEW_LENSING static inline void cart2sph(double x,double y,double z,double *r,double *cth,double *phi) @@ -458,10 +460,10 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) { long ip; double idx=par->n_grid/par->l_box; - double *fac_r_1=NULL,*fac_r_2=NULL; //Kernels for the LOS integrals #ifndef _USE_NEW_LENSING + double *fac_r_1=NULL,*fac_r_2=NULL; if(cat->has_shear) { fac_r_1=my_malloc(cat->nr*sizeof(double)); fac_r_2=my_malloc(cat->nr*sizeof(double)); @@ -609,14 +611,11 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) { int ii; double factor_vel=-par->fgrowth_0/(1.5*par->hubble_0*par->OmegaM); - int ir_s=0; #ifdef _USE_NEW_LENSING - int nside_ratio=1; - HealpixShells *smap=NULL; - if(cat->has_shear) { + int ir_s=0; + HealpixShellsAdaptive *smap=NULL; + if(cat->has_shear) smap = par->smap; - nside_ratio=smap->nside/par->nside_base; - } #endif //_USE_NEW_LENSING #ifdef _HAVE_OMP @@ -634,14 +633,16 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) if(cat->has_shear) { #ifdef _USE_NEW_LENSING int ax; - long ipix,ibase,ibase_here; + long ipix_up,ipix_lo,ibase,ibase_here; double u[3]; - double g1_0,g1_f,g2_0,g2_f; + double g1_lo,g1_up,g2_lo,g2_up; //Find lower shell index ir_s=get_r_index_smap(smap, r, ir_s); + if(ir_s>=smap->nr-1) + report_error(1, "SHIT\n"); //Find intervals - double h = (r - smap->r0[ir_s]) / (smap->r0[ir_s+1] - smap->r0[ir_s]); + double h = (r - smap->r[ir_s]) / (smap->r[ir_s+1] - smap->r[ir_s]); //Find pixel index //Find base this galaxy belongs to @@ -652,21 +653,28 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) report_error(1,"Bad base!!\n"); //Find base index within this node ibase_here=(ibase-NodeThis)/NNodes; - vec2pix_nest(smap->nside,u,&ipix); + + //Find pixel in lower edge + vec2pix_nest(smap->nside[ir_s],u,&ipix_lo); //Offset to pixel indices in this node - ipix-=(ibase-ibase_here)*nside_ratio*nside_ratio; - if((ipix<0) || (ipix>=smap->num_pix)) + ipix_lo-=ibase*smap->num_pix_per_beam[ir_s]; + if((ipix_lo<0) || (ipix_lo>=smap->num_pix_per_beam[ir_s])) report_error(1,"Bad pixel!!\n"); - //Find shear at edges - g1_0=smap->data[2*(ir_s*smap->num_pix+ipix)+0]; - g2_0=smap->data[2*(ir_s*smap->num_pix+ipix)+1]; - g1_f=smap->data[2*((ir_s+1)*smap->num_pix+ipix)+0]; - g2_f=smap->data[2*((ir_s+1)*smap->num_pix+ipix)+1]; + g1_lo=smap->data[ibase_here][ir_s][2*ipix_lo+0]; + g2_lo=smap->data[ibase_here][ir_s][2*ipix_lo+1]; + + //Same in upper edge + vec2pix_nest(smap->nside[ir_s+1],u,&ipix_up); + ipix_up-=ibase*smap->num_pix_per_beam[ir_s+1]; + if((ipix_up<0) || (ipix_up>=smap->num_pix_per_beam[ir_s+1])) + report_error(1,"Bad pixel!!\n"); + g1_up=smap->data[ibase_here][ir_s+1][2*ipix_up+0]; + g2_up=smap->data[ibase_here][ir_s+1][2*ipix_up+1]; //Interpolate - cat->srcs[ii].e1=g1_0*(1-h)+g1_f*h; - cat->srcs[ii].e2=g2_0*(1-h)+g2_f*h; + cat->srcs[ii].e1=g1_lo*(1-h)+g1_up*h; + cat->srcs[ii].e2=g2_lo*(1-h)+g2_up*h; #endif //_USE_NEW_LENSING } From 66ea6560c7509dd9656f8bbbfe0fc3785e0a0412 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 18 Jun 2020 18:23:13 +0100 Subject: [PATCH 26/29] final --- src/common.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/common.c b/src/common.c index 311bee6..4d5f76b 100644 --- a/src/common.c +++ b/src/common.c @@ -612,7 +612,7 @@ HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_ int nside_ratio; flouble r=r_arr[ir]; int nside_here=nside_base; - while(nside_here<=nside_max) { + while(nside_here Date: Thu, 18 Jun 2020 18:46:58 +0100 Subject: [PATCH 27/29] cleanup --- src/beaming.c | 2 +- src/common.c | 2 +- src/common.h | 5 +- src/io.c | 7 --- src/shear.c | 127 +------------------------------------------------- src/srcs.c | 4 +- 6 files changed, 9 insertions(+), 138 deletions(-) diff --git a/src/beaming.c b/src/beaming.c index 6dd0dfc..9041e55 100644 --- a/src/beaming.c +++ b/src/beaming.c @@ -349,7 +349,7 @@ void get_beam_properties(ParamCoLoRe *par) #ifdef _HAVE_MPI free(buffer_sr); #endif //_HAVE_MPI - + if(par->do_kappa) kappa_beams_postproc(par); #ifdef _USE_NEW_LENSING diff --git a/src/common.c b/src/common.c index 4d5f76b..059e468 100644 --- a/src/common.c +++ b/src/common.c @@ -615,7 +615,7 @@ HealpixShellsAdaptive *hp_shell_adaptive_alloc(int nq, int nside_max, int nside_ while(nside_heredo_shear) { - flouble *r_arr=compute_shear_spacing(par); - free(r_arr); - } -#endif //_USE_NEW_LENSING get_max_memory(par,test_memory+par->just_do_pred); print_info("\n"); diff --git a/src/shear.c b/src/shear.c index 78c93eb..2ab66f8 100644 --- a/src/shear.c +++ b/src/shear.c @@ -62,8 +62,8 @@ void shear_beams_preproc(ParamCoLoRe *par) #ifdef _HAVE_OMP #pragma omp for #endif //_HAVE_OMP - for(ib=0;ibsmap->nbeams;ib++) { - for(ir=0;irsmap->nr;ir++) { + for(ir=0;irsmap->nr;ir++) { + for(ib=0;ibsmap->nbeams;ib++) { for(ipp=0;ipp<2*par->smap->num_pix_per_beam[ir];ipp++) par->smap->data[ib][ir][ipp]=0; } @@ -73,129 +73,6 @@ void shear_beams_preproc(ParamCoLoRe *par) return; } -/* -void shear_get_beam_properties(ParamCoLoRe *par) -{ - HealpixShellsAdaptive *smap=par->smap; - -#ifdef _HAVE_OMP -#pragma omp parallel default(none) \ - shared(par,smap) -#endif //_HAVE_OMP - { - double idx=par->n_grid/par->l_box; - - //Setup radial decomposition - int i_r,nr; - double idr,dr; - get_radial_params(par->r_max,par->n_grid,&nr,&dr); - idr=1./dr; - - //Compute index of each source plane - int *i_r_max_arr=my_malloc(smap->nr*sizeof(int)); - int *i_r_min_arr=my_malloc(smap->nr*sizeof(int)); - double *inv_r_max=my_malloc(smap->nr*sizeof(double)); - for(i_r=0;i_rnr;i_r++) { - int i_r_here=(int)(smap->r[i_r]*idr+0.5); - inv_r_max[i_r]=1./(i_r_here*dr); - i_r_max_arr[i_r]=MIN(i_r_here,nr-1); - } - - i_r_min_arr[0]=0; - for(i_r=1;i_rnr;i_r++) - i_r_min_arr[i_r]=i_r_max_arr[i_r-1]+1; - - //Fill up integral kernels - double *fac_r_1=my_malloc(nr*sizeof(double)); - double *fac_r_2=my_malloc(nr*sizeof(double)); - for(i_r=0;i_rnum_pix;ip++) { - int ax,added; - flouble t[6]; - double r1[6],r2[6],xn[3]; - double shear1_1=0,shear1_2=0; - double shear2_1=0,shear2_2=0; - double *u=&(smap->pos[3*ip]); - double prefac=2*idx*idx; //2/Dx^2 - double cth_h=1,sth_h=0,cph_h=1,sph_h=0; - - cth_h=u[2]; - if(cth_h>=1) cth_h=1; - if(cth_h<=-1) cth_h=-1; - sth_h=sqrt((1-cth_h)*(1+cth_h)); - if(sth_h!=0) { - cph_h=u[0]/sth_h; - sph_h=u[1]/sth_h; - } - - r1[0]=(cth_h*cth_h*cph_h*cph_h-sph_h*sph_h)*prefac; - r1[1]=(2*cph_h*sph_h*(cth_h*cth_h+1))*prefac; - r1[2]=(-2*cth_h*sth_h*cph_h)*prefac; - r1[3]=(cth_h*cth_h*sph_h*sph_h-cph_h*cph_h)*prefac; - r1[4]=(-2*cth_h*sth_h*sph_h)*prefac; - r1[5]=(sth_h*sth_h)*prefac; - - r2[0]=(-2*cth_h*cph_h*sph_h)*prefac; - r2[1]=(2*cth_h*(cph_h*cph_h-sph_h*sph_h))*prefac; - r2[2]=(2*sth_h*sph_h)*prefac; - r2[3]=(2*cth_h*sph_h*cph_h)*prefac; - r2[4]=(-2*sth_h*cph_h)*prefac; - r2[5]=0; - for(i_r=0;i_rnr;i_r++) { - int irr; - int irmin=i_r_min_arr[i_r]; - int irmax=i_r_max_arr[i_r]; - for(irr=irmin;irr<=irmax;irr++) { - double rm=(irr+0.5)*dr; - for(ax=0;ax<3;ax++) - xn[ax]=(rm*u[ax]+par->pos_obs[ax])*idx; - added=interpolate_from_grid(par,xn,NULL,NULL,t,NULL,NULL,RETURN_TID,INTERP_TYPE_SHEAR); - if(added) { - double dotp1=0,dotp2=0; - for(ax=0;ax<6;ax++) { - dotp1+=r1[ax]*t[ax]; - dotp2+=r2[ax]*t[ax]; - } - shear1_1+=dotp1*fac_r_1[irr]; - shear1_2+=dotp1*fac_r_2[irr]; - shear2_1+=dotp2*fac_r_1[irr]; - shear2_2+=dotp2*fac_r_2[irr]; - } - } - smap->data[2*(i_r*smap->num_pix+ip)+0]+=(shear1_1-inv_r_max[i_r]*shear1_2); - smap->data[2*(i_r*smap->num_pix+ip)+1]+=(shear2_1-inv_r_max[i_r]*shear2_2); - } - } //end omp for - -#ifdef _HAVE_OMP -#pragma omp single -#endif //_HAVE_OMP - { - for(i_r=0;i_rnr;i_r++) - smap->r[i_r]=1./inv_r_max[i_r]; - } - - free(fac_r_1); - free(fac_r_2); - free(i_r_max_arr); - free(i_r_min_arr); - free(inv_r_max); - } //end omp parallel - - return; -} -*/ - void shear_get_beam_properties(ParamCoLoRe *par) { HealpixShellsAdaptive *smap=par->smap; diff --git a/src/srcs.c b/src/srcs.c index 6a98d43..af3b116 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -22,7 +22,7 @@ #include "common.h" #ifdef _USE_NEW_LENSING -static int get_r_index_smap(HealpixShellsAdaptive *sh,double r,int ir_start) +static int get_r_index_shear(HealpixShellsAdaptive *sh,double r,int ir_start) { int gotit=0; int ir0; @@ -637,7 +637,7 @@ static void srcs_beams_postproc_single(ParamCoLoRe *par,int ipop) double u[3]; double g1_lo,g1_up,g2_lo,g2_up; //Find lower shell index - ir_s=get_r_index_smap(smap, r, ir_s); + ir_s=get_r_index_shear(smap, r, ir_s); if(ir_s>=smap->nr-1) report_error(1, "SHIT\n"); From 42a30b70145793c68fb0e30a24e3e2dd593b9706 Mon Sep 17 00:00:00 2001 From: David Alonso Date: Thu, 18 Jun 2020 19:17:22 +0100 Subject: [PATCH 28/29] memory calculation --- src/common.c | 362 +++++++++++++++++++++++++++------------------------ 1 file changed, 193 insertions(+), 169 deletions(-) diff --git a/src/common.c b/src/common.c index 059e468..5e30f27 100644 --- a/src/common.c +++ b/src/common.c @@ -330,162 +330,6 @@ void get_random_angles(gsl_rng *rng,int ipix_nest,int ipix0,int nside,double *th } */ -unsigned long long get_max_memory(ParamCoLoRe *par,int just_test) -{ - unsigned long long total_GB=0; - unsigned long long total_GB_gau=0; - unsigned long long total_GB_lpt=0; - int fac_gau=2; - if(par->need_beaming) fac_gau=3; - - total_GB_gau=(fac_gau*(par->nz_here+1)*((long)(par->n_grid*(par->n_grid/2+1))))*sizeof(dftw_complex); - - if(par->dens_type==DENS_TYPE_1LPT) { - total_GB_lpt=(unsigned long long)(3*(1+par->lpt_buffer_fraction)*par->nz_here* - ((long)((par->n_grid/2+1)*par->n_grid))*sizeof(dftw_complex)); - } - else if(par->dens_type==DENS_TYPE_2LPT) { - total_GB_lpt=0; - total_GB_lpt=(unsigned long long)(8*(1+par->lpt_buffer_fraction)*par->nz_here* - ((long)((par->n_grid/2+1)*par->n_grid))*sizeof(dftw_complex)); - } - - unsigned long long total_GB_srcs=0; - if(par->do_srcs) { - int ipop; - for(ipop=0;ipopn_srcs;ipop++) { - int nz,ii; - long nsrc=0; - double nztot=0; - FILE *fi=fopen(par->fnameNzSrcs[ipop],"r"); - - double *zarr,*nzarr; - if(fi==NULL) error_open_file(par->fnameNzSrcs[ipop]); - nz=linecount(fi); rewind(fi); - zarr=my_malloc(nz*sizeof(double)); - nzarr=my_malloc(nz*sizeof(double)); - for(ii=0;iifnameNzSrcs[ipop],ii+1); - nzarr[ii]*=RTOD*RTOD; - } - for(ii=0;ii=par->z_min) && (zarr[ii]<=par->z_max)) - nztot+=nzarr[ii]*(zarr[ii+1]-zarr[ii]); - } - if((zarr[ii]>=par->z_min) && (zarr[ii]<=par->z_max)) - nztot+=nzarr[ii]*(zarr[ii+1]-zarr[ii]); - nztot*=4*M_PI/NNodes; - nsrc+=(long)(nztot); - if(just_test) - print_info(" Expect %ld type-%d sources\n",(long)(nztot*NNodes),ipop); - free(zarr); - free(nzarr); - fclose(fi); - - long size_source=sizeof(Src)+NPOS_CC*sizeof(float)+sizeof(int); - if(par->skw_srcs[ipop]) { - int nr=NSAMP_RAD*par->n_grid/2; - size_source+=2*nr*sizeof(float); - } - total_GB_srcs+=size_source*nsrc; - } - } - - unsigned long long total_GB_imap=0; - if(par->do_imap) { - int ipop; - for(ipop=0;ipopn_imap;ipop++) { - int nr; - unsigned long long size_map=he_nside2npix(par->nside_imap[ipop]); - FILE *fnu=fopen(par->fnameNuImap[ipop],"r"); - if(fnu==NULL) error_open_file(par->fnameNuImap[ipop]); - nr=linecount(fnu); - fclose(fnu); - total_GB_imap+=size_map*nr*(sizeof(flouble)+sizeof(int)); - } - } - - unsigned long long total_GB_kappa=0; - if(par->do_kappa) { - int ib; - int nr=par->n_kappa; - int nbases=he_nside2npix(par->nside_base); - int nside_ratio=par->nside_kappa/par->nside_base; - int npix_perbeam=nside_ratio*nside_ratio; - unsigned long long size_map=0; - for(ib=NodeThis;ibdo_shear) { - int ib; - int nr=par->n_shear; - int nbases=he_nside2npix(par->nside_base); - int nside_ratio=par->nside_shear/par->nside_base; - int npix_perbeam=nside_ratio*nside_ratio; - unsigned long long size_map=0; - for(ib=NodeThis;ibdo_isw) { - int ib; - int nr=par->n_isw; - int nbases=he_nside2npix(par->nside_base); - int nside_ratio=par->nside_isw/par->nside_base; - int npix_perbeam=nside_ratio*nside_ratio; - unsigned long long size_map=0; - for(ib=NodeThis;ibdens_type==DENS_TYPE_1LPT) || (par->dens_type==DENS_TYPE_2LPT)) - printf(", %.3lf GB (%dLPT)",(double)(total_GB_lpt/pow(1024.,3)),par->dens_type); - if(par->do_srcs) - printf(", %.3lf GB (srcs)",(double)(total_GB_srcs/pow(1024.,3))); - if(par->do_imap) - printf(", %.3lf GB (imap)",(double)(total_GB_imap/pow(1024.,3))); - if(par->do_kappa) - printf(", %.3lf MB (kappa)",(double)(total_GB_kappa/pow(1024.,2))); - if(par->do_shear) - printf(", %.3lf MB (shear)",(double)(total_GB_shear/pow(1024.,2))); - if(par->do_isw) - printf(", %.3lf MB (isw)",(double)(total_GB_isw/pow(1024.,2))); - printf("]\n"); - } -#ifdef _HAVE_MPI - MPI_Barrier(MPI_COMM_WORLD); -#endif //_HAVE_MPI - } -#endif //_DEBUG - - if(just_test==0) { - void *ptest=my_malloc(total_GB); - free(ptest); - } - - return total_GB; -} - void get_radial_params(double rmax,int ngrid,int *nr,double *dr) { *nr=NSAMP_RAD*ngrid/2; @@ -586,6 +430,26 @@ void hp_shell_adaptive_free(HealpixShellsAdaptive *shell) free(shell); } +static int *get_adaptive_nside(int nside_max, int nside_base, + int nr, flouble *r_arr, flouble dx, flouble dx_fraction) +{ + int ir; + int *nsides=my_malloc(nr*sizeof(int)); + for(ir=0; irnq=nq; shell->nr=nr; - shell->nside=my_malloc(nr*sizeof(int)); + shell->nside=get_adaptive_nside(nside_max, nside_base, nr, + r_arr, dx, dx_fraction); shell->num_pix_per_beam=my_malloc(nr*sizeof(long)); shell->r=my_malloc(nr*sizeof(flouble)); for(ir=0; irr[ir]=r; - shell->nside[ir]=nside_here; - nside_ratio=nside_here/nside_base; + shell->r[ir]=r_arr[ir]; + nside_ratio=shell->nside[ir]/nside_base; shell->num_pix_per_beam[ir]=nside_ratio*nside_ratio; } @@ -712,3 +567,172 @@ void hp_shell_free(HealpixShells *shell) free(shell->nadd); free(shell); } + +unsigned long long get_max_memory(ParamCoLoRe *par,int just_test) +{ + unsigned long long total_GB=0; + unsigned long long total_GB_gau=0; + unsigned long long total_GB_lpt=0; + int fac_gau=2; + if(par->need_beaming) fac_gau=3; + + total_GB_gau=(fac_gau*(par->nz_here+1)*((long)(par->n_grid*(par->n_grid/2+1))))*sizeof(dftw_complex); + + if(par->dens_type==DENS_TYPE_1LPT) { + total_GB_lpt=(unsigned long long)(3*(1+par->lpt_buffer_fraction)*par->nz_here* + ((long)((par->n_grid/2+1)*par->n_grid))*sizeof(dftw_complex)); + } + else if(par->dens_type==DENS_TYPE_2LPT) { + total_GB_lpt=0; + total_GB_lpt=(unsigned long long)(8*(1+par->lpt_buffer_fraction)*par->nz_here* + ((long)((par->n_grid/2+1)*par->n_grid))*sizeof(dftw_complex)); + } + + unsigned long long total_GB_srcs=0; + if(par->do_srcs) { + int ipop; + for(ipop=0;ipopn_srcs;ipop++) { + int nz,ii; + long nsrc=0; + double nztot=0; + FILE *fi=fopen(par->fnameNzSrcs[ipop],"r"); + + double *zarr,*nzarr; + if(fi==NULL) error_open_file(par->fnameNzSrcs[ipop]); + nz=linecount(fi); rewind(fi); + zarr=my_malloc(nz*sizeof(double)); + nzarr=my_malloc(nz*sizeof(double)); + for(ii=0;iifnameNzSrcs[ipop],ii+1); + nzarr[ii]*=RTOD*RTOD; + } + for(ii=0;ii=par->z_min) && (zarr[ii]<=par->z_max)) + nztot+=nzarr[ii]*(zarr[ii+1]-zarr[ii]); + } + if((zarr[ii]>=par->z_min) && (zarr[ii]<=par->z_max)) + nztot+=nzarr[ii]*(zarr[ii+1]-zarr[ii]); + nztot*=4*M_PI/NNodes; + nsrc+=(long)(nztot); + if(just_test) + print_info(" Expect %ld type-%d sources\n",(long)(nztot*NNodes),ipop); + free(zarr); + free(nzarr); + fclose(fi); + + long size_source=sizeof(Src)+NPOS_CC*sizeof(float)+sizeof(int); + if(par->skw_srcs[ipop]) { + int nr=NSAMP_RAD*par->n_grid/2; + size_source+=2*nr*sizeof(float); + } + total_GB_srcs+=size_source*nsrc; + } + } + + unsigned long long total_GB_imap=0; + if(par->do_imap) { + int ipop; + for(ipop=0;ipopn_imap;ipop++) { + int nr; + unsigned long long size_map=he_nside2npix(par->nside_imap[ipop]); + FILE *fnu=fopen(par->fnameNuImap[ipop],"r"); + if(fnu==NULL) error_open_file(par->fnameNuImap[ipop]); + nr=linecount(fnu); + fclose(fnu); + total_GB_imap+=size_map*nr*(sizeof(flouble)+sizeof(int)); + } + } + + unsigned long long total_GB_kappa=0; + if(par->do_kappa) { + int ib; + int nr=par->n_kappa; + int nbases=he_nside2npix(par->nside_base); + int nside_ratio=par->nside_kappa/par->nside_base; + int npix_perbeam=nside_ratio*nside_ratio; + unsigned long long size_map=0; + for(ib=NodeThis;ibdo_shear) { + int ib,ir; + int nr=par->n_shear; + int nbases=he_nside2npix(par->nside_base); + flouble *rs=compute_shear_spacing(par); + int *nsides=get_adaptive_nside(par->nside_shear, + par->nside_base, nr, rs, + par->l_box/par->n_grid, 1.); + int nbeams_here=0; + for(ib=NodeThis;ibnside_base; + npix_total+=nside_ratio*nside_ratio; + } + npix_total*=nbeams_here; + int nside_ratio_hi=nsides[nr-1]/par->nside_base; + unsigned long long npix_hi=nside_ratio_hi*nside_ratio_hi*nbeams_here; + total_GB_shear+=npix_total*2*sizeof(flouble)+npix_hi*3*sizeof(double); + free(rs); + free(nsides); + } + + unsigned long long total_GB_isw=0; + if(par->do_isw) { + int ib; + int nr=par->n_isw; + int nbases=he_nside2npix(par->nside_base); + int nside_ratio=par->nside_isw/par->nside_base; + int npix_perbeam=nside_ratio*nside_ratio; + unsigned long long size_map=0; + for(ib=NodeThis;ibdens_type==DENS_TYPE_1LPT) || (par->dens_type==DENS_TYPE_2LPT)) + printf(", %.3lf GB (%dLPT)",(double)(total_GB_lpt/pow(1024.,3)),par->dens_type); + if(par->do_srcs) + printf(", %.3lf GB (srcs)",(double)(total_GB_srcs/pow(1024.,3))); + if(par->do_imap) + printf(", %.3lf GB (imap)",(double)(total_GB_imap/pow(1024.,3))); + if(par->do_kappa) + printf(", %.3lf MB (kappa)",(double)(total_GB_kappa/pow(1024.,2))); + if(par->do_shear) + printf(", %.3lf MB (shear)",(double)(total_GB_shear/pow(1024.,2))); + if(par->do_isw) + printf(", %.3lf MB (isw)",(double)(total_GB_isw/pow(1024.,2))); + printf("]\n"); + } +#ifdef _HAVE_MPI + MPI_Barrier(MPI_COMM_WORLD); +#endif //_HAVE_MPI + } +#endif //_DEBUG + + if(just_test==0) { + void *ptest=my_malloc(total_GB); + free(ptest); + } + + return total_GB; +} From 7ee73e42ee90e4c3576558a415776b1f5fde507c Mon Sep 17 00:00:00 2001 From: David Alonso Date: Mon, 6 Jul 2020 14:08:57 +0100 Subject: [PATCH 29/29] CCL test files --- examples/LSST/param.cfg | 12 +++ examples/cl_test/cl_test_ccl.py | 154 ++++++++++++++++++++++++++++++++ examples/cl_test/param.cfg | 11 ++- src/predictions.c | 4 +- 4 files changed, 178 insertions(+), 3 deletions(-) create mode 100644 examples/cl_test/cl_test_ccl.py diff --git a/examples/LSST/param.cfg b/examples/LSST/param.cfg index 278acc0..ffb4e6a 100644 --- a/examples/LSST/param.cfg +++ b/examples/LSST/param.cfg @@ -55,3 +55,15 @@ kappa: z_out= [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6] nside= 512 } + +shear: +{ + nside= 2048 + # Select number of slices + # For these parameters, 11 corresponds roughly + # to dr_shear = 100 Mpc/h. + n_shear= 50 + # Select spacing type ("r" or "log(1+z)") + spacing_type= "r" + write=false +} diff --git a/examples/cl_test/cl_test_ccl.py b/examples/cl_test/cl_test_ccl.py new file mode 100644 index 0000000..2e8b476 --- /dev/null +++ b/examples/cl_test/cl_test_ccl.py @@ -0,0 +1,154 @@ +import numpy as np +import healpy as hp +import matplotlib.pyplot as plt +import pyccl as ccl +import os +from astropy.io import fits + +# Hubble constant +h = 0.7 + +nside = 128 +npix = hp.nside2npix(nside) + +# We will make two bins (z_photo < 0.15 and z_photo > 0.15) +nbins = 2 +nz_h = 50 +zsplit = 0.15 + +# This will contain the N(z) of the different bins +nz_tot = np.zeros([nbins, nz_h]) +sigz = 0.03 + +# These will be the density and ellipticity maps +dmap = np.zeros([nbins, npix]) +e1map = np.zeros([nbins, npix]) +e2map = np.zeros([nbins, npix]) + +# Now we loop over all source files +ifile = 0 +while os.path.isfile('out_srcs_s1_%d.fits' % ifile): + print("HI") + hdulist = fits.open('out_srcs_s1_%d.fits' % ifile) + d = hdulist[1].data + n_g = len(d) + + # Generate random photo-z + z_photo = d['Z_COSMO'] + sigz*(1+d['Z_COSMO'])*np.random.randn(n_g) + + # Split into 2 + msk1 = z_photo <= zsplit + msk2 = z_photo > zsplit + + # For each bin, add to the number and ellipticity maps + for ibin, msk in enumerate([msk1, msk2]): + dd = d[msk] + + pix = hp.ang2pix(nside, + np.radians(90-dd['DEC']), + np.radians(dd['RA'])) + n = np.bincount(pix, minlength=npix) + e1 = np.bincount(pix, minlength=npix, weights=dd['E1']) + e2 = np.bincount(pix, minlength=npix, weights=dd['E2']) + dmap[ibin, :] += n + e1map[ibin, :] += e1 + e2map[ibin, :] += e2 + + # Add also to N(z) + nz, z_edges = np.histogram(dd['Z_COSMO'], bins=nz_h, + range=[0., 0.5]) + nz_tot[ibin, :] += nz + ifile += 1 + +# Midpoint of N(z) histogram +z_nz = 0.5*(z_edges[1:] + z_edges[:-1]) + +# Compute map and overdensity map +# Compute also shot noise level +shotnoise = np.zeros(nbins) +for ib in range(nbins): + ndens = (np.sum(dmap[ib])+0.0)/(4*np.pi) + shotnoise[ib] = 1./ndens + e1map[ib, :] = e1map[ib]/dmap[ib] + e1map[ib, dmap[ib] <= 0] = 0 + e2map[ib, :] = e2map[ib]/dmap[ib] + e2map[ib, dmap[ib] <= 0] = 0 + dmap[ib, :] = (dmap[ib, :] + 0.0) / np.mean(dmap[ib] + 0.0) - 1 + +# Read P(k) theory prediction +zs = [] +pks_dd = [] +pks_dm = [] +pks_mm = [] +z_pk = 0. +while os.path.isfile("out_pk_srcs_pop0_z%.3lf.txt" % z_pk): + ks, pdd, pdm, pmm = np.loadtxt("out_pk_srcs_pop0_z%.3lf.txt" % z_pk, unpack=True) + # The delta-delta prediction involves some Fourier transforms that make it unstable + # at high-k, so we just set it to zero. + pdd[pmm<1E-30] = 0 + pks_dd.append(pdd) + pks_dm.append(pdm) + pks_mm.append(pmm) + zs.append(z_pk) + z_pk += 0.1 +# Reverse order (because CCL needs increasing scale factor - not redshift). +# Also, CCL uses non-h units. +zs = np.array(zs)[::-1] +ks = np.array(ks) * h +pks_dd = np.array(pks_dd)[::-1, :] / h**3 +pks_dm = np.array(pks_dm)[::-1, :] / h**3 +pks_mm = np.array(pks_mm)[::-1, :] / h**3 + +# Create CCL P(k) structures +cosmo = ccl.Cosmology(Omega_c=0.25, Omega_b=0.05, h=h, n_s=0.96, sigma8=0.8) +pk2d_dd = ccl.Pk2D(a_arr=1./(1+zs), lk_arr=np.log(ks), pk_arr=pks_dd, + is_logp=False, extrap_order_hik=0, extrap_order_lok=0) +pk2d_dm = ccl.Pk2D(a_arr=1./(1+zs), lk_arr=np.log(ks), pk_arr=pks_dm, + is_logp=False, extrap_order_hik=0, extrap_order_lok=0) +pk2d_mm = ccl.Pk2D(a_arr=1./(1+zs), lk_arr=np.log(ks), pk_arr=pks_mm, + is_logp=False, extrap_order_hik=0, extrap_order_lok=0) +# These can be evaluated as follows: +print(pk2d_dd.eval(k=0.1, a=1/(1+0.2), cosmo=cosmo)) + +# Create a number counts and a weak lensing tracer for each of the bins +tr_d = [ccl.NumberCountsTracer(cosmo, False, (z_nz, nz_tot[i]), bias=(z_nz, np.ones_like(z_nz))) + for i in range(nbins)] +tr_l = [ccl.WeakLensingTracer(cosmo, (z_nz, nz_tot[i])) for i in range(nbins)] + +# Compute power spectra. I'm only doing delta-delta here. +pairs=[(0,0), (0,1), (1,1)] +cl_dd_d = np.array([hp.anafast(dmap[p1], dmap[p2]) for p1, p2 in pairs]) +larr = np.arange(3*nside) +cl_dd_t = np.array([ccl.angular_cl(cosmo, tr_d[p1], tr_d[p2], larr, p_of_k_a=pk2d_mm) + for p1, p2 in pairs]) + +# Compare C_ells with theory +for i, (p1, p2) in enumerate(pairs): + plt.figure() + plt.title('d-%d d-%d' % (p1, p2)) + if p1 == p2: + nl = shotnoise[p1] + plt.plot(larr, nl*np.ones_like(larr), 'k--', label='Shot noise') + else: + nl = 0 + + # Rebin to reduce noise + cld = np.mean((cl_dd_d[i]-nl).reshape([-1, 6]), axis=1) + clt = np.mean(cl_dd_t[i].reshape([-1, 6]), axis=1) + l = np.mean(larr.reshape([-1, 6]), axis=1) + msk = l<2*nside + plt.plot(l[msk], cld[msk], 'k.', label='Sim') + plt.plot(l[msk], clt[msk], 'r-', label='Prediction') + plt.loglog() + plt.xlabel(r'$\ell$', fontsize=15) + plt.ylabel(r'$C_\ell$', fontsize=15) + plt.legend(loc='lower left') + +# Plot N(z)s +plt.figure() +for n in nz_tot: + plt.plot(z_nz, n) +plt.plot(z_nz, np.sum(nz_tot, axis=0), 'k-') +plt.xlabel(r'$z$', fontsize=15) +plt.ylabel(r'$N(z)$', fontsize=15) +plt.show() diff --git a/examples/cl_test/param.cfg b/examples/cl_test/param.cfg index 480404a..78e7322 100644 --- a/examples/cl_test/param.cfg +++ b/examples/cl_test/param.cfg @@ -7,7 +7,7 @@ global: z_min= 0.001 z_max= 0.400 seed= 1001 - write_pred=false + write_pred=true just_write_pred=false pred_dz=0.1 } @@ -53,3 +53,12 @@ isw: z_out= [0.38] nside= 128 } + +shear: +{ + nside= 128 + n_shear= 11 + spacing_type= "r" + write= true + } + \ No newline at end of file diff --git a/src/predictions.c b/src/predictions.c index a7f8b09..0a2b25c 100644 --- a/src/predictions.c +++ b/src/predictions.c @@ -31,8 +31,8 @@ void write_predictions(ParamCoLoRe *par) // first generate k array, sufficiently finely spaced // note that we need to sufficiently pad on both ends const int Nk=10000; - const double kmin=1e-3; - const double kmax=50;; + const double kmin=1e-4; + const double kmax=100; const double kminout=kmin; const double kmaxout=kmax; const double rminout=0.5;