diff --git a/Makefile b/Makefile index 45bb8cd..39a785f 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) @@ -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" @@ -19,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" @@ -120,6 +122,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 +130,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/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.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/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/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..fd3e34f --- /dev/null +++ b/examples/shear_test/param.cfg @@ -0,0 +1,61 @@ +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= true + store_skewers= false +} + +kappa: +{ + z_out= [2.076379E-01] + nside= 128 +} + +shear: +{ + nside= 256 + # 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" + write=false +} 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/beaming.c b/src/beaming.c index 0c9a3fa..334549b 100644 --- a/src/beaming.c +++ b/src/beaming.c @@ -298,14 +298,18 @@ 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); +#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) + srcs_beams_preproc(par); + if(par->do_imap) + imap_beams_preproc(par); if(NodeThis==0) timer(0); @@ -329,27 +333,35 @@ 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); +#ifdef _USE_NEW_LENSING + if(par->do_shear) + shear_get_beam_properties(par); +#endif //_USE_NEW_LENSING 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); +#ifdef _USE_NEW_LENSING + if(par->do_shear) + shear_beams_postproc(par); +#endif //_USE_NEW_LENSING 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 5bc858b..c427254 100644 --- a/src/common.c +++ b/src/common.c @@ -330,141 +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_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_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; @@ -498,13 +363,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)); @@ -546,7 +413,96 @@ void catalog_free(Catalog *cat) free(cat); } -HealpixShells *hp_shell_alloc(int nside,int nside_base,int nr) +void hp_shell_adaptive_free(HealpixShellsAdaptive *shell) +{ + int ib,ir; + free(shell->r); + free(shell->nside); + free(shell->num_pix_per_beam); + for(ib=0;ibnbeams;ib++) { + for(ir=0;irnr;ir++) + free(shell->data[ib][ir]); + free(shell->data[ib]); + free(shell->pos[ib]); + } + free(shell->data); + free(shell->pos); + 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; irNSIDE_MAX_HPX) + report_error(1,"Can't go beyond nside=%d\n",NSIDE_MAX_HPX); + if(nside_maxnbeams=0; + for(ib=NodeThis;ibnbeams++; + + shell->nq=nq; + shell->nr=nr; + 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_arr[ir]; + nside_ratio=shell->nside[ir]/nside_base; + shell->num_pix_per_beam[ir]=nside_ratio*nside_ratio; + } + + 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; +} + +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 +521,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 +545,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; @@ -610,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; +} diff --git a/src/common.h b/src/common.h index 3c9a408..f917c5b 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 @@ -175,6 +179,7 @@ typedef struct { double rmax; double dr; double idr; + int has_shear; int has_skw; int skw_gauss; float *d_skw; @@ -183,6 +188,18 @@ typedef struct { } Catalog; typedef struct { + int nq; + int nr; //Number of spherical shells + flouble *r; //r in this shell + int *nside; //Resolution parameter + long *num_pix_per_beam; + int nbeams; //Number of different bases + double **pos; //3D positions of the pixels in the furthest shell + flouble ***data; +} HealpixShellsAdaptive; + +typedef struct { + int nq; //Number of quantities to map int nside; //Resolution parameter long num_pix; long *listpix; @@ -274,7 +291,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) @@ -315,6 +332,13 @@ typedef struct { flouble **fl_mean_extra_kappa; flouble **cl_extra_kappa; #endif //_ADD_EXTRA_KAPPA + // Shear + 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; + HealpixShellsAdaptive *smap; //Shear maps at each redshift // ISW int do_isw; //Do you want to output isw maps? int n_isw; //How many maps? @@ -359,11 +383,14 @@ 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); +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_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) @@ -390,6 +417,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 @@ -399,6 +427,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); @@ -462,6 +491,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..3384d46 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); @@ -769,3 +788,22 @@ void compute_tracer_cosmo(ParamCoLoRe *par) csm_params_free(pars); } + +flouble *compute_shear_spacing(ParamCoLoRe *par) +{ + flouble *rarr; + int ir,nr=par->n_shear; + rarr=my_malloc(nr*sizeof(flouble)); + if(par->shear_spacing_type==SPACING_R) { + flouble dr=par->r_max/nr; + for(ir=0;irshear_spacing_type==SPACING_LOGZ) { + flouble dlogz=log(1+par->z_max)/nr; + for(ir=0;irn_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 c748aa1..4f561a7 100644 --- a/src/io.c +++ b/src/io.c @@ -88,22 +88,28 @@ 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 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->shear_spacing_type=SPACING_R; 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; + par->write_shear=0; for(ii=0;iifnameBzSrcs[ii],"default"); sprintf(par->fnameNzSrcs[ii],"default"); @@ -132,6 +138,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 @@ -315,7 +322,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; @@ -354,11 +361,33 @@ 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; 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)); } +#ifdef _USE_NEW_LENSING + //Shear maps + cset=config_lookup(conf,"shear"); + if(cset!=NULL) { + par->do_shear=1; + par->do_srcs_shear=1; + 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)); + conf_read_bool(conf,"shear","write",&(par->write_shear)); + } + // 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) { par->do_isw=1; @@ -398,12 +427,13 @@ 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_srcs_shear+par->do_kappa+par->do_shear+par->do_isw+par->do_skewers; init_fftw(par); + cosmo_set(par); + get_max_memory(par,test_memory+par->just_do_pred); - cosmo_set(par); print_info("\n"); double dk=2*M_PI/par->l_box; @@ -426,10 +456,12 @@ 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) - 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"); @@ -463,19 +495,31 @@ 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); + +#ifdef _USE_NEW_LENSING + if(par->do_shear) { + flouble *r_arr=compute_shear_spacing(par); + 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.); + free(r_arr); + } +#endif //_USE_NEW_LENSING 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); + free(conf); + return par; } @@ -741,6 +785,71 @@ void write_kappa(ParamCoLoRe *par) print_info("\n"); } +void write_shear(ParamCoLoRe *par) +{ + int ir; + char fname[256]; + 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)); + if(NodeThis==0) timer(0); + print_info("*** Writing shear source maps\n"); + for(ir=0;irsmap->nr;ir++) { + long ip, ib; + + //Write local pixels to dummy map + for(ip=0;ipprefixOut,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++) { + 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]; + } + } + + //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); +#endif //_HAVE_MPI + + //Write dummy map + if(NodeThis==0) + he_write_healpix_map(map_write,2,par->smap->nside[ir],fname,1); + } + free(map_write[0]); + free(map_write[1]); + free(map_write); + + 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\n",ir, + par->smap->r[ir], + get_bg(par, par->smap->r[ir], BG_Z, 0)); + } + fclose(f); + } + + if(NodeThis==0) timer(2); + print_info("\n"); +} + void write_isw(ParamCoLoRe *par) { int ir; @@ -868,7 +977,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->cats[ipop]->has_shear) nfields=7; print_info(" %d-th population (FITS)\n",ipop); @@ -901,7 +1010,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->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; } @@ -911,7 +1020,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->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); } @@ -994,7 +1103,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->cats[ipop]->has_shear) fprintf(fil,"#[6]e1, [7]e2\n"); else fprintf(fil,"\n"); @@ -1002,7 +1111,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->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"); @@ -1082,6 +1191,13 @@ void param_colore_free(ParamCoLoRe *par) #endif //_ADD_EXTRA_KAPPA } +#ifdef _USE_NEW_LENSING + if(par->do_shear) { + if(par->smap!=NULL) + hp_shell_adaptive_free(par->smap); + } +#endif //_USE_NEW_LENSING + if(par->do_isw) { if(par->pd_map!=NULL) hp_shell_free(par->pd_map); @@ -1102,4 +1218,5 @@ void param_colore_free(ParamCoLoRe *par) #ifdef _DEBUG fclose(par->f_dbg); #endif //_DEBUG + free(par); } diff --git a/src/main.c b/src/main.c index 6ea08dd..f5cff83 100644 --- a/src/main.c +++ b/src/main.c @@ -70,51 +70,67 @@ int main(int argc,char **argv) //Compute normalization of density field for biasing 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); +#ifdef _USE_NEW_LENSING + if(par->do_shear) + shear_set_cartesian(par); +#endif //_USE_NEW_LENSING 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); +#ifdef _USE_NEW_LENSING + if(par->do_shear) + shear_distribute(par); +#endif //_USE_NEW_LENSING 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); +#ifdef _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) + 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) 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); +#ifdef _USE_NEW_LENSING + if(par->do_shear && par->write_shear) + write_shear(par); +#endif ///_USE_NEW_LENSING if(par->do_isw) write_isw(par); + if(par->do_srcs) + write_srcs(par); + if(par->do_imap) + write_imap(par); } print_info("\n"); 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; diff --git a/src/shear.c b/src/shear.c new file mode 100644 index 0000000..2ab66f8 --- /dev/null +++ b/src/shear.c @@ -0,0 +1,218 @@ +/////////////////////////////////////////////////////////////////////// +// // +// 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 *r_i=my_malloc(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->r); + for(ir=0;irsmap->nr;ir++) { + par->smap->r[ir]=r_i[i_sorted[ir]]; + } + free(r_i); + free(i_sorted); + + //Zero all data and clear +#ifdef _HAVE_OMP +#pragma omp parallel default(none) \ + shared(par) +#endif //_HAVE_OMP + { + long ir, ib, ipp; + +#ifdef _HAVE_OMP +#pragma omp for +#endif //_HAVE_OMP + 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; + } + } //end omp for + } //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); + free(fac_r_2); + 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; +} + +void shear_beams_postproc(ParamCoLoRe *par) +{ + //Set rf to end of cell + return; +} diff --git a/src/srcs.c b/src/srcs.c index 307c25a..ce2ea14 100644 --- a/src/srcs.c +++ b/src/srcs.c @@ -21,6 +21,50 @@ /////////////////////////////////////////////////////////////////////// #include "common.h" +#ifdef _USE_NEW_LENSING +static int get_r_index_shear(HealpixShellsAdaptive *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(rr[1]) + gotit=1; + else + ir0++; + } + else if(ir0==sh->nr-1) { + if(r>=sh->r[sh->nr-2]) { + ir0=sh->nr-2; + gotit=1; + } + else + ir0--; + } + else { + if(rr[ir0]) + ir0--; + else { + if(r>=sh->r[ir0+1]) + ir0++; + else + gotit=1; + } + } + } + + return ir0; +} +#endif //_USE_NEW_LENSING + + static inline void cart2sph(double x,double y,double z,double *r,double *cth,double *phi) { *r=sqrt(x*x+y*y+z*z); @@ -208,17 +252,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); @@ -340,7 +384,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 @@ -413,10 +460,11 @@ 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 - if(par->do_lensing) { +#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)); @@ -427,6 +475,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 @@ -459,12 +508,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) { @@ -478,8 +527,9 @@ static void srcs_get_beam_properties_single(ParamCoLoRe *par,int ipop) } } +#ifndef _USE_NEW_LENSING //Compute lensing shear - if(par->do_lensing) { + 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; @@ -531,11 +581,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 - if(par->do_lensing) { +#ifndef _USE_NEW_LENSING + if(cat->has_shear) { free(fac_r_1); free(fac_r_2); } +#endif //_USE_NEW_LENSING }//end omp parallel } @@ -549,14 +602,21 @@ 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); +#ifdef _USE_NEW_LENSING + int ir_s=0; + HealpixShellsAdaptive *smap=NULL; + if(cat->has_shear) + smap = par->smap; +#endif //_USE_NEW_LENSING #ifdef _HAVE_OMP #pragma omp for @@ -570,7 +630,53 @@ 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_up,ipix_lo,ibase,ibase_here; + double u[3]; + double g1_lo,g1_up,g2_lo,g2_up; + //Find lower shell index + ir_s=get_r_index_shear(smap, r, ir_s); + if(ir_s>=smap->nr-1) + report_error(1, "SHIT\n"); + + //Find intervals + 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 + 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; + + //Find pixel in lower edge + vec2pix_nest(smap->nside[ir_s],u,&ipix_lo); + //Offset to pixel indices in this node + 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_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_lo*(1-h)+g1_up*h; + cat->srcs[ii].e2=g2_lo*(1-h)+g2_up*h; +#endif //_USE_NEW_LENSING + } //Skewers if(cat->has_skw) {