This notebook shows, for a given redshift distribution (provided by Joe Zuntz), how we (Patricia Larsen and Nan Li) calculate the effective convergence and shear maps by using Argonne Simulations, which is designed for cosmoDC2. We also attached some codes for the calculation of the power spectrums of weak lensing maps at the bottom of this notebook. Some cross-checks will be implemented later, for example, comparing our results to the theoretical predictions created by pyccl. Should you have any questions or are interested, please feel free to slack us (@linan7788626, @plarsen, @joezuntz). The maps are located on NERSC, and if you want to take a look at the data, please let me know.
%matplotlib inline
import numpy as np
import pylab as pl
import healpy as hp
from scipy import interpolate
from scipy.integrate import quad
# path to the simulation data on Cooley, which is an Argonne cluster, please contact Patricia for more information
direc_shears = </path/to/the/simulations/data/>
def loadinWeakLensingMaps(step):
'''
Load in the data in the directory above.
The data files contain integrated (0 to zs) weak lensing maps for a given source plane.
The input parameter <step> stands for the snapshotID corresponding to the source plane.
'''
try:
a00 = hp.read_map(direc_shears+'rt_'+step+'_A1_2.fits',0)
a01 = hp.read_map(direc_shears+'rt_'+step+'_A1_2.fits',1)
a10 = hp.read_map(direc_shears+'rt_'+step+'_A2_2.fits',0)
a11 = hp.read_map(direc_shears+'rt_'+step+'_A2_2.fits',1)
except:
print("Can not find SnapShot %s"%step)
return None
kappa_m = 1.0 - (a00 + a11)*0.5 # kappa map
shear1_m = (a11-a00)*0.5 # shear1 map
shear2_m = -(a01+a10)*0.5 # shear2 map
return kappa_m, shear1_m, shear2_m
def weights_nz_dz(zs_arr,steps):
'''
Calculate the weights of the effective weak lensing maps for
a given redshift distribution of sources, i.e., n(z).
There are <nstep> redshift bins for sources,
- the zeroth bin is [0, (zs[0]+zs[1])/2); zbin = zs[0].
- the last bin is [(zs[-2]+zs[-1])/2), 1.0); zbin = zs[-1]
- the ith bin is [(zs[i-1]+z[i])/2, (zs[i]+z[i+1])/2); zbin = zs[i]
By integrating the normalized n(z) in the bins,
we can obtain the weights for the source planes.
'''
zs_in, ng_in = np.loadtxt("./source_1.txt", usecols=(0, 1), unpack=True) # n(z) from Joe Zuntz
dzs = zs_in[1:]-zs_in[:-1]
hng = (ng_in[1:]+ng_in[:-1])*0.5
area_tot = np.sum(hng*dzs)
ng_rescale = ng_in/area_tot # normalization
f = interpolate.interp1d(zs_in, ng_rescale) # the function created by interpolation
fnz_norm = f(zs_arr)
zs_min = 0.0
zs_max = 3.0
wt_arr = zs_arr*0.0
zinit = 200.
nsteps = 500
a = np.linspace(1./(zinit+1.),1.,nsteps)
z = 1./a - 1.
for i in range(len(zs_arr)):
if i==0:
zs_1 = zs_min
zs_m = zs_arr[0]
zs_2 = z[int(steps[i])]
wt_arr[i] = quad(f, zs_1,zs_2)[0]
elif i==(len(zs_arr)-1):
zs_1 = z[int(steps[i-1])]
zs_m= zs_arr[-1]
zs_2 = zs_max
wt_arr[i] = quad(f, zs_1,zs_2)[0]
else:
zs_1 = z[int(steps[i-1])]
zs_m= zs_arr[i]
zs_2 = z[int(steps[i])]
wt_arr[i] = quad(f, zs_1,zs_2)[0]
return wt_arr
def cal_eff_wl_maps():
'''
calculate the total effective weak lensing maps by summing the list of maps with the corresponding
weights.
'''
# Initial the redshifts of source planes
data_list = {}
data_list['steps'] = np.array(['487', '475', '464', '453', '442', '432', '421', '411',
'401', '392', '382', '373', '365', '355', '347', '338',
'331', '323', '315', '307', '300', '293', '286', '279',
'272', '266', '259', '253', '247', '241', '235', '230',
'224', '219', '213', '208', '203', '198', '194', '189',
'184', '180', '176', '171', '167', '163', '159', '155',
'151', '148', '144', '141', '137', '134', '131', '127',
'124', '121'])
data_list['zs'] = np.array([0.019587576, 0.04249394, 0.06348896, 0.08982432, 0.11601174,
0.14167726, 0.17023492, 0.19876349, 0.22816956, 0.25715542,
0.28774858, 0.31976724, 0.3496554, 0.38333392, 0.4181242,
0.45316005, 0.4872712, 0.5210507, 0.5587207, 0.5984776,
0.63728297, 0.67555356, 0.7154497, 0.7575482, 0.8018985,
0.8445264, 0.8898808, 0.9368595, 0.9829167, 1.030879,
1.0813811, 1.1297112, 1.1807017, 1.2338057, 1.2901752,
1.3487234, 1.4050395, 1.4638758, 1.5192988, 1.5775385,
1.6453991, 1.7094843, 1.7690942, 1.83975, 1.9136083,
1.9828148, 2.0554118, 2.131457, 2.2114959, 2.2847886,
2.3618498, 2.4422584, 2.5272007, 2.6155977, 2.6953585,
2.7935054, 2.895901, 2.9888024])
nsteps = len(data_list['zs'])
weights = weights_nz_dz(data_list['zs'], data_list['steps'])
kappa_map_last, shear1_map_last, shear2_map_last = loadinWeakLensingMaps(data_list['steps'][-1])
kappa_map = kappa_map_last*weights[-1]
shear1_map = shear1_map_last*weights[-1]
shear2_map = shear2_map_last*weights[-1]
for i in range(nsteps-1):
res_tmp = loadinWeakLensingMaps(data_list['steps'][i])
if res_tmp==None:
continue
else:
kappa_m_tmp, \
shear1_m_tmp, \
shear2_m_tmp = res_tmp
kappa_map = kappa_map + kappa_m_tmp*weights[i]
shear1_map = shear1_map + shear1_m_tmp*weights[i]
shear2_map = shear2_map + shear2_m_tmp*weights[i]
return kappa_map, shear1_map, shear2_map
def vis_kappa_eff(kappa_map):
'''
A simple function to visualize an octant convergence map with healpix.
'''
nside = 4096
#-----------------------------------
# Create masks according to the octant properties
#
x,y,z = hp.pix2vec(nside, np.arange(hp.nside2npix(nside)))
mask_octant = (x>0)&(y>0)&(z<0)
#-----------------------------------
# setup the colormap
#
cmap = pl.cm.jet
cmap.set_over(cmap(1.0))
cmap.set_under('w')
cmap.set_bad('gray')
#-----------------------------------
# make a plot in healpix map
#
kappa_map = hp.ma(kappa_map)
kappa_map.mask = np.logical_not(mask_octant)
hp.mollview(kappa_map, cmap=cmap)
return 0
'''
The main function to calculate kappa and shear maps. It takes about 8 minutes.
'''
ka_map, s1_map, s2_map = cal_eff_wl_maps()
'''
Write the maps to a fits file with healpy
'''
hp.write_map('WL_MAPs_3.fits',m=(ka_map,s1_map,s2_map),overwrite=True)
'''
Functions for the calculation power spectrum, please slack Patricia (@plarsen) for more details.
'''
def upgrade_pixels(pix_list):
pix_list_new=np.array([],dtype=np.int64)
for i in range(4):
pix_list_new = np.concatenate((pix_list_new,pix_list*4+i))
return pix_list_new
def compute_masks():
nside=4096
x,y,z = hp.pix2vec(nside, np.arange(hp.nside2npix(nside)))
mask_octant = (x>0)&(y>0)&(z<0)
pix_list_image = [8786, 8787, 8788, 8789, 8790, 8791, 8792, 8793, 8794, 8913, 8914, 8915, 8916, 8917, 8918, 8919, 8920, 8921, 9042, 9043, 9044, 9045, 9046, 9047, 9048, 9049,
9050, 9169, 9170, 9171, 9172, 9173, 9174, 9175, 9176, 9177, 9178, 9298, 9299, 9300, 9301, 9302, 9303, 9304, 9305, 9306, 9425, 9426, 9427, 9428, 9429, 9430,
9431, 9432, 9433, 9434, 9554, 9555, 9556, 9557, 9558, 9559, 9560, 9561, 9562, 9681, 9682, 9683, 9684, 9685, 9686, 9687, 9688, 9689, 9690, 9810, 9811, 9812,
9813, 9814, 9815, 9816, 9817, 9818, 9937, 9938, 9939, 9940, 9941, 9942, 9943, 9944, 9945, 9946, 10066, 10067, 10068, 10069, 10070, 10071, 10072, 10073, 10074, 10193,
10194, 10195, 10196, 10197, 10198, 10199, 10200, 10201, 10202, 10321, 10322, 10323, 10324, 10325, 10326, 10327, 10328, 10329, 10444, 10445, 10446, 10447, 10448, 10449, 10450,
10451, 10452]
pix_list_nest = hp.ring2nest(32,pix_list_image)
for i in range(int(np.log2(nside/32))):
pix_list_nest = upgrade_pixels(pix_list_nest)
pix_list_image = hp.nest2ring(nside,pix_list_nest)
mask_image = np.zeros(hp.nside2npix(nside))
mask_image[pix_list_image]=1.0
return mask_octant,mask_image
def compute_power(map_array,mask,lmax=1500):
'''compute power spectra for a given array of maps with a given mask (anafast with fsky mask correction)'''
fsky = np.sum(mask)/(len(mask)+0.0)
alms_wl = hp.map2alm(map_array,lmax=lmax)
cls_wl = hp.alm2cl(alms_wl)/fsky/(hp.pixwin(4096)[:lmax+1]**2) # assuming octant map
return cls_wl
'''
compute and plot the power spectrums of kappa and shear maps.
'''
# hp.smoothing(kappa,lmax=lmax,fwhm= (arcmin_smoothing/60.)*(np.pi/180.)) #smoothing if needed
# read maps from the fits file
kappa0 = hp.read_map('WL_MAPs_3.fits',0)
shear1 = hp.read_map('WL_MAPs_3.fits',1)
shear2 = hp.read_map('WL_MAPs_3.fits',2)
# create masks in healpix space
mask_oct, mask_image = compute_masks()
# compute the power spectrums of the maps with masks.
lmax=1500
cls = compute_power((kappa0*mask_oct,shear1*mask_oct,shear2*mask_oct),mask_oct,lmax=lmax)
# plot the power spectrums
pl.figure()
l = np.arange(lmax+1)
pl.loglog(l,cls[0],label='convergence power')
pl.loglog(l,cls[1],label='shear power')