Skip to content

Commit

Permalink
Merge pull request #55 from damonge/faster_shear_b
Browse files Browse the repository at this point in the history 
Faster shear
  • Loading branch information
@damonge
damonge authored Jul 6, 2020
2 parents 7497897 + 421f0cc commit 9e357a7
Show file tree
Hide file tree
Showing 24 changed files with 2,042 additions and 515 deletions.
9 changes: 6 additions & 3 deletions Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -12,14 +12,16 @@ 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"
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"
Expand Down Expand Up @@ -120,14 +122,15 @@ 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
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
Expand Down
12 changes: 12 additions & 0 deletions examples/LSST/param.cfg
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
}
179 changes: 105 additions & 74 deletions examples/cl_test/cl_test.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Loading

0 comments on commit 9e357a7

Please sign in to comment.