From afc8abf78f714d1232d15b9e645330bd9a0288ec Mon Sep 17 00:00:00 2001 From: data_vadar <46063942+sayansaha1996@users.noreply.github.com> Date: Fri, 24 Jun 2022 18:36:21 +0530 Subject: [PATCH 01/23] Updated profile.py Updated to have the xmax option for analitic_kappa_ft(). --- lensit/clusterlens/profile.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py index 9be6ad9..890a154 100644 --- a/lensit/clusterlens/profile.py +++ b/lensit/clusterlens/profile.py @@ -166,10 +166,16 @@ def sigma_int(self, M200, z, R, xmax = 100, npoints=1000): sigma[i] = np.trapz( 2 * r_arr * rho_arr / np.sqrt(r_arr**2 - iR**2), r_arr) return sigma - def analitic_kappa_ft(self, M200, z, ell, const_c=None): + def analitic_kappa_ft(self, M200, z, ell, const_c=None, xmax=None): """Analytic Fourier transform of the convergence fiels for a NFW profile from Oguri&Takada 2010, Eq.28""" - c = self.get_concentration(M200, z, const_c) + """ If xmax not given it calculates the function for a NFW profile truncated + at R200 = c*rs. If xmax is given, then it calculates the function + for a NFW profile truncated at xmax*rs""" + if xmax is None: + c = self.get_concentration(M200, z, const_c) + else: + c=xmax mu_nfw = np.log(1. + c) - c / (1. + c) rs = self.get_rs(M200, z) chi = self.cosmo.comoving_radial_distance(z) From 186a6019cade7d1e84caaf3543a56ccfce4eb96b Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Tue, 1 Nov 2022 14:58:54 +0100 Subject: [PATCH 02/23] Update lensingmap.py --- lensit/clusterlens/lensingmap.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index b3283ce..7170ae3 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -17,7 +17,7 @@ def get_cluster_libdir(cambinifile, profilename, npix, lpix_amin, ellmax_sky, M2 class cluster_maps(object): - def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, cmb_exp='5muKamin_1amin', cache_maps=False): + def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, ellmin_data=10, cmb_exp='5muKamin_1amin', cache_maps=False): """Library for flat-sky CMB simulations lensed by a galaxy cluster. Args: @@ -75,10 +75,12 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd # Set the lmax of the data maps self.ellmax_data = ellmax_data + self.ellmin_data = ellmin_data self.cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=self.ellmax_data) # The resolution of the data map could be lower than the resolution of the sky map (=the true map) - self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= self.ellmax_data, num_threads=self.num_threads) + #self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= self.ellmax_data, num_threads=self.num_threads) + self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: np.logical_and(ell <= self.ellmax_data , ell>=self.ellmin_data), num_threads=self.num_threads) vcell_amin2 = np.prod(self.lib_datalm.ell_mat.lsides) / np.prod(self.lib_datalm.ell_mat.shape) * (180 * 60. / np.pi) ** 2 nTpix = sN_uKamin / np.sqrt(vcell_amin2) From 77c934686f4a576c73dcaa90eda6089b58bb4113 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Tue, 1 Nov 2022 15:13:35 +0100 Subject: [PATCH 03/23] Update profile.py --- lensit/clusterlens/profile.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py index 3ed5a00..6053b11 100644 --- a/lensit/clusterlens/profile.py +++ b/lensit/clusterlens/profile.py @@ -166,16 +166,10 @@ def sigma_int(self, M200, z, R, xmax = 100, npoints=1000): sigma[i] = np.trapz( 2 * r_arr * rho_arr / np.sqrt(r_arr**2 - iR**2), r_arr) return sigma - def analitic_kappa_ft(self, M200, z, ell, const_c=None, xmax=None): + def analitic_kappa_ft(self, M200, z, ell, const_c=None): """Analytic Fourier transform of the convergence fiels for a NFW profile from Oguri&Takada 2010, Eq.28""" - """ If xmax not given it calculates the function for a NFW profile truncated - at R200 = c*rs. If xmax is given, then it calculates the function - for a NFW profile truncated at xmax*rs""" - if xmax is None: - c = self.get_concentration(M200, z, const_c) - else: - c=xmax + c = self.get_concentration(M200, z, const_c) mu_nfw = np.log(1. + c) - c / (1. + c) rs = self.get_rs(M200, z) chi = self.cosmo.comoving_radial_distance(z) From 58d0a7b5fb000f59c67b620d17a14eacf3616de4 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Sat, 7 Jan 2023 13:07:12 +0100 Subject: [PATCH 04/23] Update lensingmap.p Update lensing map to add ellmindata, and class sim_cmb_len_lss_plus_cluster --- lensit/clusterlens/lensingmap.py | 176 +++++++++++++++++++++++++++++-- 1 file changed, 168 insertions(+), 8 deletions(-) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index 7170ae3..4ded4b5 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -1,3 +1,5 @@ +from __future__ import print_function + from lensit.clusterlens.profile import profile import numpy as np from lensit import ell_mat, pbs @@ -9,19 +11,30 @@ import os from camb import CAMBdata +import os +import pickle as pk + +import numpy as np + +from lensit.sims.ffs_cmbs import sim_cmb_unl, sims_cmb_len +from lensit.pbs import pbs +from lensit.misc.misc_utils import npy_hash +from lensit.sims import ffs_phas, sims_generic +from lensit.ffs_deflect import ffs_deflect -def get_cluster_libdir(cambinifile, profilename, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp): - return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename, +def get_cluster_libdir(cambinifile, profilename, key, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp): + return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename,"lensed_by_%s" % key, 'npix%s_lpix_%samin_lmaxsky%s' % (npix, lpix_amin, ellmax_sky), 'M200_%.6E_z%s' % (M200, z), '%s_sims' % nsims) class cluster_maps(object): - def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, ellmin_data=10, cmb_exp='5muKamin_1amin', cache_maps=False): + def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, ellmin_data=10, cmb_exp='5muKamin_1amin', cache_maps=False): """Library for flat-sky CMB simulations lensed by a galaxy cluster. Args: lib_dir: various things will be cached there + key: a keyword among "lss", "cluster" and "lss_plus_cluster", which is gonna decide how the patch has been lensed npix: number of pixels on one side of the square box lpix_amin: physical size (in arcmin) of the pixels nsims: number of CMB maps to simulate @@ -30,7 +43,9 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd profilename: string defining the density profile of the cluster (e.g. nfw) ellmax_sky: maximum multipole of CMB spectra used to generate the CMB maps """ + assert key in ['lss', 'cluster', 'lss_plus_cluster'], key self.libdir = libdir + self.key = key self.cosmo = cosmo if profilename == 'nfw': self.M200 = profparams['M200c'] @@ -50,10 +65,20 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd self.num_threads=int(os.environ.get('OMP_NUM_THREADS', 1)) self.lib_skyalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, num_threads=self.num_threads, filt_func=lambda ell: ell <= self.ellmax_sky) - nfields = 3 + #nfields = 4 + if self.key == 'cluster': + nfields = 3 + else: + nfields = 4 camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T - self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]} + cpp_fid = cosmo.get_lens_potential_cls(lmax=self.ellmax_sky, raw_cl=True).T[0] + + if self.key == 'cluster': + self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]} + else: + self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} + for cl in self.cls_unl.values(): cl[6000:] = 0 @@ -71,8 +96,19 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd # Get the noise and beam of the experiment sN_uKamin, sN_uKaminP, Beam_FWHM_amin, ellmin, ellmax = li.get_config(self.cmb_exp) - self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) - + #self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + #self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + #self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha) + #self.len_cmbs = sim_cmb_unl(self.cls_unl, lib_pha=skypha) + if self.key == 'lss': + self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha) + elif self.key == 'cluster': + self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + elif self.key == 'lss_plus_cluster': + self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + else: + assert 0 + # Set the lmax of the data maps self.ellmax_data = ellmax_data self.ellmin_data = ellmin_data @@ -179,13 +215,25 @@ def get_kappa_map(self, M200, z, xmax=None): """ return self.haloprofile.kappa_map(M200, z, self.lib_skyalm.shape, self.lib_skyalm.lsides, xmax) +verbose = False + +def get_fields(cls): + fields = ['p', 't', 'e', 'b', 'o'] + ret = ['p', 't', 'e', 'b', 'o'] + for f in fields: + if not ((f + f) in cls.keys()): ret.remove(f) + for k in cls.keys(): + for f in k: + if f not in ret: ret.append(f) + return ret class sim_cmb_len_cluster(ffs_cmbs.sims_cmb_len): def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): super(sim_cmb_len_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) - self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) + self.conc_par = haloprofile.get_concentration(M200,z) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) self.defl_map = haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." @@ -193,3 +241,115 @@ def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:f def _get_f(self, idx=None): """We refine the displacement field using the convergence map of the cluster""" return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + +class sim_cmb_len_lss_plus_cluster: + def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): + #super(sim_cmb_len_lss_plus_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) + + self.M200 = M200 + self.z = z + self.haloprofile = haloprofile + self.conc_par = self.haloprofile.get_concentration(M200,z) + self.kappa_map = self.haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) + self.defl_map = self.haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) + + if not os.path.exists(lib_dir) and pbs.rank == 0: + os.makedirs(lib_dir) + pbs.barrier() + self.lib_skyalm = lib_skyalm + fields = get_fields(cls_unl) + if lib_pha is None and pbs.rank == 0: + lib_pha = ffs_phas.ffs_lib_phas(os.path.join(lib_dir, 'phas'), len(fields), lib_skyalm) + else: # Check that the lib_alms are compatible : + assert lib_pha.lib_alm == lib_skyalm + pbs.barrier() + + self.unlcmbs = sim_cmb_unl(cls_unl, lib_pha) + self.Pool = use_Pool + self.cache_lens = cache_lens + fn_hash = os.path.join(lib_dir, 'sim_hash.pk') + if not os.path.exists(fn_hash) and pbs.rank == 0: + pk.dump(self.hashdict(), open(fn_hash, 'wb'), protocol=2) + pbs.barrier() + sims_generic.hash_check(self.hashdict(), pk.load(open(fn_hash, 'rb'))) + self.lib_dir = lib_dir + self.fields = fields + + #assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." + + def _get_f_cluster(self, idx=None): + """We refine the displacement field using the convergence map of the cluster""" + return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + + def hashdict(self): + return {'unl_cmbs': self.unlcmbs.hashdict(), + 'M200': self.M200, + 'z': self.z}#, + #'haloprofile': self.haloprofile} + + def is_full(self): + return self.unlcmbs.lib_pha.is_full() + + def get_sim_plm(self, idx): + return self.unlcmbs.get_sim_plm(idx) + + def get_sim_olm(self, idx): + return self.unlcmbs.get_sim_olm(idx) + + def _get_f(self, idx): + if 'p' in self.unlcmbs.fields and 'o' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_frompolm(self.lib_skyalm, plm, olm, verbose=False) + elif 'p' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + return ffs_deflect.displacement_fromplm(self.lib_skyalm, plm) + elif 'o' in self.unlcmbs.fields: + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False) + else: + assert 0 + + def get_sim_alm(self, idx, field): + if field == 't': + return self.get_sim_tlm(idx) + elif field == 'p': + return self.get_sim_plm(idx) + elif field == 'o': + return self.get_sim_olm(idx) + elif field == 'q': + return self.get_sim_qulm(idx)[0] + elif field == 'u': + return self.get_sim_qulm(idx)[1] + elif field == 'e': + return self.lib_skyalm.QUlms2EBalms(self.get_sim_qulm(idx))[0] + elif field == 'b': + return self.lib_skyalm.QUlms2EBalms(self.get_sim_qulm(idx))[1] + else: + assert 0, (field, self.fields) + + def get_sim_tlm(self, idx): + fname = os.path.join(self.lib_dir, 'sim_%04d_tlm.npy'%idx) + if not os.path.exists(fname): + Tlm = self._get_f(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool) + Tlm = self._get_f_cluster(idx).lens_alm(self.lib_skyalm, Tlm, use_Pool=self.Pool) + if not self.cache_lens: return Tlm + np.save(fname, Tlm) + return np.load(fname) + + def get_sim_qulm(self, idx): + fname = os.path.join(self.lib_dir, 'sim_%04d_qulm.npy'%idx) + if not os.path.exists(fname): + Qlm, Ulm = self.lib_skyalm.EBlms2QUalms( + np.array([self.unlcmbs.get_sim_elm(idx), self.unlcmbs.get_sim_blm(idx)])) + f = self._get_f(idx) + Qlm = f.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool) + Ulm = f.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool) + f_cluster = self._get_f_cluster(idx) + Qlm = f_cluster.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool) + Ulm = f_cluster.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool) + if not self.cache_lens: return np.array([Qlm, Ulm]) + np.save(fname, np.array([Qlm, Ulm])) + return np.load(fname) + + From 691d66d18e4aa76f238a24ed462f1c63675f9334 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Fri, 12 May 2023 20:01:24 +0200 Subject: [PATCH 05/23] Update lensingmap.py Final update --- lensit/clusterlens/lensingmap.py | 40 +++++++++++++++++++++++++------- 1 file changed, 32 insertions(+), 8 deletions(-) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index 4ded4b5..c3d7ed4 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -67,17 +67,22 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co #nfields = 4 if self.key == 'cluster': - nfields = 3 + nfields = 4 else: nfields = 4 camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T + #camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T/100 + + #camb_cls = cosmo.get_lensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T + cpp_fid = cosmo.get_lens_potential_cls(lmax=self.ellmax_sky, raw_cl=True).T[0] - if self.key == 'cluster': - self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]} - else: - self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} + #if self.key == 'cluster': + # self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]} + #else: + # self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} + self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} for cl in self.cls_unl.values(): cl[6000:] = 0 @@ -233,10 +238,12 @@ def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:f super(sim_cmb_len_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) self.conc_par = haloprofile.get_concentration(M200,z) - self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= 3 * self.conc_par) + #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) + #self.kappa_map = 1e2*haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) self.defl_map = haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) - assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." + #assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." def _get_f(self, idx=None): """We refine the displacement field using the convergence map of the cluster""" @@ -250,7 +257,8 @@ def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:f self.z = z self.haloprofile = haloprofile self.conc_par = self.haloprofile.get_concentration(M200,z) - self.kappa_map = self.haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= 3 * self.conc_par) + #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) self.defl_map = self.haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) if not os.path.exists(lib_dir) and pbs.rank == 0: @@ -309,6 +317,20 @@ def _get_f(self, idx): return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False) else: assert 0 + + def _get_f_lss_plus_cluster(self, idx): + if 'p' in self.unlcmbs.fields and 'o' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_frompolm(self.lib_skyalm, plm, olm, verbose=False) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + elif 'p' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + return ffs_deflect.displacement_fromplm(self.lib_skyalm, plm) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + elif 'o' in self.unlcmbs.fields: + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + else: + assert 0 def get_sim_alm(self, idx, field): if field == 't': @@ -333,6 +355,7 @@ def get_sim_tlm(self, idx): if not os.path.exists(fname): Tlm = self._get_f(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool) Tlm = self._get_f_cluster(idx).lens_alm(self.lib_skyalm, Tlm, use_Pool=self.Pool) + #Tlm = self._get_f_lss_plus_cluster(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool) if not self.cache_lens: return Tlm np.save(fname, Tlm) return np.load(fname) @@ -343,6 +366,7 @@ def get_sim_qulm(self, idx): Qlm, Ulm = self.lib_skyalm.EBlms2QUalms( np.array([self.unlcmbs.get_sim_elm(idx), self.unlcmbs.get_sim_blm(idx)])) f = self._get_f(idx) + #f = self._get_f_lss_plus_cluster(idx) Qlm = f.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool) Ulm = f.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool) f_cluster = self._get_f_cluster(idx) From 7e0979a641a92bb742d1b42fc3a712a526e5aa19 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Fri, 12 May 2023 20:09:12 +0200 Subject: [PATCH 06/23] Update profile.py corrected the centre problem in numpy.fft, analytical fourier transform expression, and added cap for fx and gx --- lensit/clusterlens/profile.py | 55 +++++++++++++++++++++++++++-------- 1 file changed, 43 insertions(+), 12 deletions(-) diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py index 6053b11..8cf326a 100644 --- a/lensit/clusterlens/profile.py +++ b/lensit/clusterlens/profile.py @@ -41,7 +41,7 @@ def get_rho_s(self, M200, z, const_c=None): def get_rs(self, M200, z, const_c=None): return self.get_r200(M200, z) / self.get_concentration(M200, z, const_c) - def get_thetas_amin(M200,z, const_c=None): + def get_thetas_amin(self, M200,z, const_c=None): return self.r_to_theta(z, self.get_rs(M200, z, const_c)) def get_r200(self, M200, z): @@ -101,26 +101,36 @@ def sigma_nfw(self, M200, z, R, xmax=None): def fx(self, x, tol=6): """This integral of the NFW profile along the line of sight is integrated up to infinity See Bartelmann 1996 or Wright et al 1999""" + ## Need to fix for M200, z, now it is for M200, z = 2*1e14, 0.7 f = np.zeros_like(x) xp = np.where(x>1) xo = np.where(np.abs(x-1.) < 10**(-tol)) - xm = np.where(x<1) + xm = np.logical_and(x<1, x>=self.theta_amin_to_x(2*1e14, 0.7, 0.12)) + #xm = np.where(x<1) + xd = np.where(xxmax, i.e. R>rs*xmax See equation 27 of Takada and Jain 2003""" + ## Need to fix for M200, z, now it is for M200, z = 2*1e14, 0.7 g = np.zeros_like(x) - xp = np.where(np.logical_and(x>1, x1, x<=xmax) xo = np.where(np.abs(x-1.) < 10**(-tol)) - xm = np.where(x<1) + xm = np.logical_and(x<1, x>=self.theta_amin_to_x(2*1e14, 0.7, 0.12)) + #xm = np.where(x<1) + xd = np.where(x Date: Sat, 20 May 2023 18:19:27 +0200 Subject: [PATCH 07/23] Update __init__.py --- lensit/__init__.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/lensit/__init__.py b/lensit/__init__.py index 859f1b6..4f5cc7d 100644 --- a/lensit/__init__.py +++ b/lensit/__init__.py @@ -249,6 +249,16 @@ def get_config(exp): Beam_FWHM_amin = 1. ellmin=10 ellmax=3000 + elif exp == 'S4_sayan': + sN_uKamin = 1. + Beam_FWHM_amin = 1. + ellmin=100 + ellmax=4000 + elif exp == 'CMB_HD': + sN_uKamin = 0.5 + Beam_FWHM_amin = 0.5 + ellmin=100 + ellmax=4000 elif exp == 'Planck45': sN_uKamin = 45. Beam_FWHM_amin = 5 From f077c0fa1be2e46805fe2bd5a4444b74e3d6094c Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Sat, 27 May 2023 18:25:59 +0200 Subject: [PATCH 08/23] Update ffs_deflect.py --- lensit/ffs_deflect/ffs_deflect.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/lensit/ffs_deflect/ffs_deflect.py b/lensit/ffs_deflect/ffs_deflect.py index fe437e7..222b16f 100644 --- a/lensit/ffs_deflect/ffs_deflect.py +++ b/lensit/ffs_deflect/ffs_deflect.py @@ -275,6 +275,7 @@ def lens_alm_adjoint(self, lib_alm: ffs_alm, alm: np.ndarray, lib_alm_out: ffs_a plan.setpts(ys, xs) self._bwdnufftplan = plan m = lib_alm.alm2map(alm).reshape(sy * sx) + m = m.astype(np.complex128, copy=False) ret = self._bwdnufftplan.execute(m) return lib_alm_out.rfftmap2alm(ret[:ry, :rx]) @@ -716,4 +717,4 @@ def displacement_fromolm(lib_plm, olm, **kwargs): def displacement_frompolm(lib_plm, plm, olm, **kwargs): return ffs_displacement(lib_plm.alm2map(plm * lib_plm.get_ikx() - olm * lib_plm.get_iky()), lib_plm.alm2map(plm * lib_plm.get_iky() + olm * lib_plm.get_ikx()), - lib_plm.ell_mat.lsides, **kwargs) \ No newline at end of file + lib_plm.ell_mat.lsides, **kwargs) From 493b5119ef84e26947a581112a2bbba4aae72fc4 Mon Sep 17 00:00:00 2001 From: data_vadar Date: Tue, 6 Jun 2023 14:56:12 +0200 Subject: [PATCH 09/23] added a new line in README --- README.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 814d783..ba5f730 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,7 @@ # Lensit +This is a new line added for this branch clusterlens + [![PyPI version](https://badge.fury.io/py/lensit.svg)](https://badge.fury.io/py/lensit)[![alt text](https://readthedocs.org/projects/lensit/badge/?version=latest)](https://lensit.readthedocs.io/en/latest)[![Build Status](https://travis-ci.com/carronj/lensit.svg?branch=master)](https://travis-ci.com/carronj/lensit) This is a set of python tools dedicated to CMB lensing and CMB delensing, by Julien Carron. @@ -38,4 +40,4 @@ The notebook 'demo_lensit.ipynb' shows an example of iterative lensing map recon Other example and tests scripts might follow, or you may just write to me. ![alt text](https://erc.europa.eu/sites/default/files/content/erc_banner-vertical.jpg) -![SNSF logo](./docs/SNF_logo_standard_web_color_neg_e.svg) \ No newline at end of file +![SNSF logo](./docs/SNF_logo_standard_web_color_neg_e.svg) From 7dbfd9e6a5780f742afaf238b180d58e584b24d9 Mon Sep 17 00:00:00 2001 From: data_vadar Date: Tue, 6 Jun 2023 15:15:05 +0200 Subject: [PATCH 10/23] updated everything for clusterlens --- lensit/clusterlens/lensingmap.py | 206 ++++++++++++++++++++++++++++-- lensit/clusterlens/profile.py | 55 ++++++-- lensit/ffs_deflect/ffs_deflect.py | 3 +- 3 files changed, 241 insertions(+), 23 deletions(-) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index b3283ce..c3d7ed4 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -1,3 +1,5 @@ +from __future__ import print_function + from lensit.clusterlens.profile import profile import numpy as np from lensit import ell_mat, pbs @@ -9,19 +11,30 @@ import os from camb import CAMBdata +import os +import pickle as pk + +import numpy as np + +from lensit.sims.ffs_cmbs import sim_cmb_unl, sims_cmb_len +from lensit.pbs import pbs +from lensit.misc.misc_utils import npy_hash +from lensit.sims import ffs_phas, sims_generic +from lensit.ffs_deflect import ffs_deflect -def get_cluster_libdir(cambinifile, profilename, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp): - return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename, +def get_cluster_libdir(cambinifile, profilename, key, npix, lpix_amin, ellmax_sky, M200, z, nsims, cmbexp): + return opj(li._get_lensitdir()[0], 'temp', 'clustermaps', 'camb_%s' % cambinifile, 'cmbexp_%s' % cmbexp, '%s_profile' % profilename,"lensed_by_%s" % key, 'npix%s_lpix_%samin_lmaxsky%s' % (npix, lpix_amin, ellmax_sky), 'M200_%.6E_z%s' % (M200, z), '%s_sims' % nsims) class cluster_maps(object): - def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, cmb_exp='5muKamin_1amin', cache_maps=False): + def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBdata, profparams:dict, profilename='nfw', ellmax_sky = 6000, ellmax_data = 3000, ellmin_data=10, cmb_exp='5muKamin_1amin', cache_maps=False): """Library for flat-sky CMB simulations lensed by a galaxy cluster. Args: lib_dir: various things will be cached there + key: a keyword among "lss", "cluster" and "lss_plus_cluster", which is gonna decide how the patch has been lensed npix: number of pixels on one side of the square box lpix_amin: physical size (in arcmin) of the pixels nsims: number of CMB maps to simulate @@ -30,7 +43,9 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd profilename: string defining the density profile of the cluster (e.g. nfw) ellmax_sky: maximum multipole of CMB spectra used to generate the CMB maps """ + assert key in ['lss', 'cluster', 'lss_plus_cluster'], key self.libdir = libdir + self.key = key self.cosmo = cosmo if profilename == 'nfw': self.M200 = profparams['M200c'] @@ -50,10 +65,25 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd self.num_threads=int(os.environ.get('OMP_NUM_THREADS', 1)) self.lib_skyalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, num_threads=self.num_threads, filt_func=lambda ell: ell <= self.ellmax_sky) - nfields = 3 + #nfields = 4 + if self.key == 'cluster': + nfields = 4 + else: + nfields = 4 camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T - self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]} + #camb_cls = cosmo.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T/100 + + #camb_cls = cosmo.get_lensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=self.ellmax_sky).T + + cpp_fid = cosmo.get_lens_potential_cls(lmax=self.ellmax_sky, raw_cl=True).T[0] + + #if self.key == 'cluster': + # self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3]} + #else: + # self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} + self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} + for cl in self.cls_unl.values(): cl[6000:] = 0 @@ -71,14 +101,27 @@ def __init__(self, libdir:str, npix:int, lpix_amin:float, nsims:int, cosmo:CAMBd # Get the noise and beam of the experiment sN_uKamin, sN_uKaminP, Beam_FWHM_amin, ellmin, ellmax = li.get_config(self.cmb_exp) - self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) - + #self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + #self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + #self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha) + #self.len_cmbs = sim_cmb_unl(self.cls_unl, lib_pha=skypha) + if self.key == 'lss': + self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha) + elif self.key == 'cluster': + self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + elif self.key == 'lss_plus_cluster': + self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + else: + assert 0 + # Set the lmax of the data maps self.ellmax_data = ellmax_data + self.ellmin_data = ellmin_data self.cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=self.ellmax_data) # The resolution of the data map could be lower than the resolution of the sky map (=the true map) - self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= self.ellmax_data, num_threads=self.num_threads) + #self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: ell <= self.ellmax_data, num_threads=self.num_threads) + self.lib_datalm = ell_mat.ffs_alm_pyFFTW(self.ellmat, filt_func=lambda ell: np.logical_and(ell <= self.ellmax_data , ell>=self.ellmin_data), num_threads=self.num_threads) vcell_amin2 = np.prod(self.lib_datalm.ell_mat.lsides) / np.prod(self.lib_datalm.ell_mat.shape) * (180 * 60. / np.pi) ** 2 nTpix = sN_uKamin / np.sqrt(vcell_amin2) @@ -177,17 +220,160 @@ def get_kappa_map(self, M200, z, xmax=None): """ return self.haloprofile.kappa_map(M200, z, self.lib_skyalm.shape, self.lib_skyalm.lsides, xmax) +verbose = False + +def get_fields(cls): + fields = ['p', 't', 'e', 'b', 'o'] + ret = ['p', 't', 'e', 'b', 'o'] + for f in fields: + if not ((f + f) in cls.keys()): ret.remove(f) + for k in cls.keys(): + for f in k: + if f not in ret: ret.append(f) + return ret class sim_cmb_len_cluster(ffs_cmbs.sims_cmb_len): def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): super(sim_cmb_len_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) - self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) + self.conc_par = haloprofile.get_concentration(M200,z) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= 3 * self.conc_par) + #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) + #self.kappa_map = 1e2*haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) self.defl_map = haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) - assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." + #assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." def _get_f(self, idx=None): """We refine the displacement field using the convergence map of the cluster""" return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + +class sim_cmb_len_lss_plus_cluster: + def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): + #super(sim_cmb_len_lss_plus_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) + + self.M200 = M200 + self.z = z + self.haloprofile = haloprofile + self.conc_par = self.haloprofile.get_concentration(M200,z) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= 3 * self.conc_par) + #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) + self.defl_map = self.haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) + + if not os.path.exists(lib_dir) and pbs.rank == 0: + os.makedirs(lib_dir) + pbs.barrier() + self.lib_skyalm = lib_skyalm + fields = get_fields(cls_unl) + if lib_pha is None and pbs.rank == 0: + lib_pha = ffs_phas.ffs_lib_phas(os.path.join(lib_dir, 'phas'), len(fields), lib_skyalm) + else: # Check that the lib_alms are compatible : + assert lib_pha.lib_alm == lib_skyalm + pbs.barrier() + + self.unlcmbs = sim_cmb_unl(cls_unl, lib_pha) + self.Pool = use_Pool + self.cache_lens = cache_lens + fn_hash = os.path.join(lib_dir, 'sim_hash.pk') + if not os.path.exists(fn_hash) and pbs.rank == 0: + pk.dump(self.hashdict(), open(fn_hash, 'wb'), protocol=2) + pbs.barrier() + sims_generic.hash_check(self.hashdict(), pk.load(open(fn_hash, 'rb'))) + self.lib_dir = lib_dir + self.fields = fields + + #assert 'p' not in self.fields, "Remove the lensing potential power spectrum from the input cls_unl dictionnary to avoid errors." + + def _get_f_cluster(self, idx=None): + """We refine the displacement field using the convergence map of the cluster""" + return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + + def hashdict(self): + return {'unl_cmbs': self.unlcmbs.hashdict(), + 'M200': self.M200, + 'z': self.z}#, + #'haloprofile': self.haloprofile} + + def is_full(self): + return self.unlcmbs.lib_pha.is_full() + + def get_sim_plm(self, idx): + return self.unlcmbs.get_sim_plm(idx) + + def get_sim_olm(self, idx): + return self.unlcmbs.get_sim_olm(idx) + + def _get_f(self, idx): + if 'p' in self.unlcmbs.fields and 'o' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_frompolm(self.lib_skyalm, plm, olm, verbose=False) + elif 'p' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + return ffs_deflect.displacement_fromplm(self.lib_skyalm, plm) + elif 'o' in self.unlcmbs.fields: + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False) + else: + assert 0 + + def _get_f_lss_plus_cluster(self, idx): + if 'p' in self.unlcmbs.fields and 'o' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_frompolm(self.lib_skyalm, plm, olm, verbose=False) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + elif 'p' in self.unlcmbs.fields: + plm = self.get_sim_plm(idx) + return ffs_deflect.displacement_fromplm(self.lib_skyalm, plm) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + elif 'o' in self.unlcmbs.fields: + olm = self.get_sim_olm(idx) + return ffs_deflect.displacement_fromolm(self.lib_skyalm, olm, verbose=False) + ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) + else: + assert 0 + + def get_sim_alm(self, idx, field): + if field == 't': + return self.get_sim_tlm(idx) + elif field == 'p': + return self.get_sim_plm(idx) + elif field == 'o': + return self.get_sim_olm(idx) + elif field == 'q': + return self.get_sim_qulm(idx)[0] + elif field == 'u': + return self.get_sim_qulm(idx)[1] + elif field == 'e': + return self.lib_skyalm.QUlms2EBalms(self.get_sim_qulm(idx))[0] + elif field == 'b': + return self.lib_skyalm.QUlms2EBalms(self.get_sim_qulm(idx))[1] + else: + assert 0, (field, self.fields) + + def get_sim_tlm(self, idx): + fname = os.path.join(self.lib_dir, 'sim_%04d_tlm.npy'%idx) + if not os.path.exists(fname): + Tlm = self._get_f(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool) + Tlm = self._get_f_cluster(idx).lens_alm(self.lib_skyalm, Tlm, use_Pool=self.Pool) + #Tlm = self._get_f_lss_plus_cluster(idx).lens_alm(self.lib_skyalm, self.unlcmbs.get_sim_tlm(idx), use_Pool=self.Pool) + if not self.cache_lens: return Tlm + np.save(fname, Tlm) + return np.load(fname) + + def get_sim_qulm(self, idx): + fname = os.path.join(self.lib_dir, 'sim_%04d_qulm.npy'%idx) + if not os.path.exists(fname): + Qlm, Ulm = self.lib_skyalm.EBlms2QUalms( + np.array([self.unlcmbs.get_sim_elm(idx), self.unlcmbs.get_sim_blm(idx)])) + f = self._get_f(idx) + #f = self._get_f_lss_plus_cluster(idx) + Qlm = f.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool) + Ulm = f.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool) + f_cluster = self._get_f_cluster(idx) + Qlm = f_cluster.lens_alm(self.lib_skyalm, Qlm, use_Pool=self.Pool) + Ulm = f_cluster.lens_alm(self.lib_skyalm, Ulm, use_Pool=self.Pool) + if not self.cache_lens: return np.array([Qlm, Ulm]) + np.save(fname, np.array([Qlm, Ulm])) + return np.load(fname) + + diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py index 6053b11..8cf326a 100644 --- a/lensit/clusterlens/profile.py +++ b/lensit/clusterlens/profile.py @@ -41,7 +41,7 @@ def get_rho_s(self, M200, z, const_c=None): def get_rs(self, M200, z, const_c=None): return self.get_r200(M200, z) / self.get_concentration(M200, z, const_c) - def get_thetas_amin(M200,z, const_c=None): + def get_thetas_amin(self, M200,z, const_c=None): return self.r_to_theta(z, self.get_rs(M200, z, const_c)) def get_r200(self, M200, z): @@ -101,26 +101,36 @@ def sigma_nfw(self, M200, z, R, xmax=None): def fx(self, x, tol=6): """This integral of the NFW profile along the line of sight is integrated up to infinity See Bartelmann 1996 or Wright et al 1999""" + ## Need to fix for M200, z, now it is for M200, z = 2*1e14, 0.7 f = np.zeros_like(x) xp = np.where(x>1) xo = np.where(np.abs(x-1.) < 10**(-tol)) - xm = np.where(x<1) + xm = np.logical_and(x<1, x>=self.theta_amin_to_x(2*1e14, 0.7, 0.12)) + #xm = np.where(x<1) + xd = np.where(xxmax, i.e. R>rs*xmax See equation 27 of Takada and Jain 2003""" + ## Need to fix for M200, z, now it is for M200, z = 2*1e14, 0.7 g = np.zeros_like(x) - xp = np.where(np.logical_and(x>1, x1, x<=xmax) xo = np.where(np.abs(x-1.) < 10**(-tol)) - xm = np.where(x<1) + xm = np.logical_and(x<1, x>=self.theta_amin_to_x(2*1e14, 0.7, 0.12)) + #xm = np.where(x<1) + xd = np.where(x Date: Tue, 6 Jun 2023 12:28:05 -0300 Subject: [PATCH 11/23] filtering grad leg --- lensit/ffs_covs/ffs_cov.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/lensit/ffs_covs/ffs_cov.py b/lensit/ffs_covs/ffs_cov.py index f11d592..82b2eb9 100644 --- a/lensit/ffs_covs/ffs_cov.py +++ b/lensit/ffs_covs/ffs_cov.py @@ -632,7 +632,7 @@ def apply_conddiagcl(self, typ, alms, use_cls_len=True, use_Pool=0): def apply_condpseudiagcl(self, typ, alms, use_Pool=0): return self.apply_conddiagcl(typ, alms, use_Pool=use_Pool) - def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, **kwargs): + def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, ellmax_gradleg=None, **kwargs): r"""Unormalized quadratic estimates (potential and curl). Note: @@ -646,7 +646,7 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, **kwargs): iblms: inverse-variance filtered CMB alm arrays lib_qlm: *ffs_alm* instance describing the lensing alm arrays use_cls_len: use lensed or unlensed cls in QE weights (numerator), defaults to lensed cls - + ellmax_gradleg: cut scales above this value on the gradient leg (used for unbiased cluster lensing reconstruction) """ assert iblms.shape == self._skyalms_shape(typ), (iblms.shape, self._skyalms_shape(typ)) @@ -660,6 +660,12 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, **kwargs): for _j in range(len(typ)): clms[_i] += pmat.get_unlPmat_ij(typ, self.lib_skyalm, weights_cls, _i, _j) * iblms[_j] + if ellmax_gradleg is not None: + fl = np.ones(self.lib_skyalm.ellmax+1) + fl[ellmax_gradleg+1:] = 0 + for _i in range(len(typ)): + lib_qlm.almxfl(clms[_i], fl, inplace=True) + t.checkpoint(" get_qlms::mult with %s Pmat" % ({True: 'len', False: 'unl'}[use_cls_len])) _map = lambda alm: self.lib_skyalm.alm2map(alm) From 1f9f967265f3d7667e23cf36fc8ad4f750c21edf Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Wed, 7 Jun 2023 14:39:32 +0200 Subject: [PATCH 12/23] Update README.md --- README.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/README.md b/README.md index ba5f730..e473049 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,5 @@ # Lensit -This is a new line added for this branch clusterlens - [![PyPI version](https://badge.fury.io/py/lensit.svg)](https://badge.fury.io/py/lensit)[![alt text](https://readthedocs.org/projects/lensit/badge/?version=latest)](https://lensit.readthedocs.io/en/latest)[![Build Status](https://travis-ci.com/carronj/lensit.svg?branch=master)](https://travis-ci.com/carronj/lensit) This is a set of python tools dedicated to CMB lensing and CMB delensing, by Julien Carron. From 661e66f469392e9fdfc66b58267b2ff9cba4dcb3 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Wed, 7 Jun 2023 14:41:07 +0200 Subject: [PATCH 13/23] Update __init__.py --- lensit/__init__.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/lensit/__init__.py b/lensit/__init__.py index 4f5cc7d..859f1b6 100644 --- a/lensit/__init__.py +++ b/lensit/__init__.py @@ -249,16 +249,6 @@ def get_config(exp): Beam_FWHM_amin = 1. ellmin=10 ellmax=3000 - elif exp == 'S4_sayan': - sN_uKamin = 1. - Beam_FWHM_amin = 1. - ellmin=100 - ellmax=4000 - elif exp == 'CMB_HD': - sN_uKamin = 0.5 - Beam_FWHM_amin = 0.5 - ellmin=100 - ellmax=4000 elif exp == 'Planck45': sN_uKamin = 45. Beam_FWHM_amin = 5 From 6fd1b3183392a7d8f9853577b987966f06f4f8f1 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Wed, 7 Jun 2023 14:44:19 +0200 Subject: [PATCH 14/23] Update README.md From c56dbd49f94de3a449912f52c49349438e0e978a Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Wed, 7 Jun 2023 14:47:52 +0200 Subject: [PATCH 15/23] Update README.md From 9a34d139c4f9d183b648064d9932ef9cd13a088f Mon Sep 17 00:00:00 2001 From: data_vadar Date: Wed, 7 Jun 2023 15:06:32 +0200 Subject: [PATCH 16/23] Update ffs_maps.py Added unlensed observed map in simulations --- lensit/sims/ffs_maps.py | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) diff --git a/lensit/sims/ffs_maps.py b/lensit/sims/ffs_maps.py index 1fc14e8..e9e43c3 100644 --- a/lensit/sims/ffs_maps.py +++ b/lensit/sims/ffs_maps.py @@ -231,6 +231,19 @@ def _loadQnoise(self): def _loadUnoise(self): return np.load(self.nUpix) + def _build_sim_tmap_unl(self, idx): + tmap = self.lib_skyalm.almxfl(self.lencmbs.unlcmbs.get_sim_tlm(idx), self.cl_transf) + tmap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, tmap)) + return tmap + self.get_noise_sim_tmap(idx) + + def _build_sim_qumap_unl(self, idx): + qmap, umap = self.lencmbs.unlcmbs.get_sim_qulm(idx) + self.lib_skyalm.almxfl(qmap, self.cl_transf, inplace=True) + self.lib_skyalm.almxfl(umap, self.cl_transf, inplace=True) + qmap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, qmap)) + umap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, umap)) + return np.array([qmap + self.get_noise_sim_qmap(idx), umap + self.get_noise_sim_umap(idx)]) + def _build_sim_tmap(self, idx): tmap = self.lib_skyalm.almxfl(self.lencmbs.get_sim_tlm(idx), self.cl_transf) tmap = self.lib_datalm.alm2map(self.lib_datalm.udgrade(self.lencmbs.lib_skyalm, tmap)) @@ -253,6 +266,28 @@ def get_noise_sim_qmap(self, idx): def get_noise_sim_umap(self, idx): return self._loadUnoise() * self.pix_pha.get_sim(idx, 2) + def get_sim_tmap_unl(self, idx): + fn = os.path.join(self.lib_dir, 'sim_tmap_unl_%05d.npy'%idx) + if self.cache_sims and not os.path.exists(fn): + if pbs.rank == 0: + np.save(fn, self._build_sim_tmap_unl(idx)) + pbs.barrier() + if self.cache_sims: + return np.load(fn) + else: + return self._build_sim_tmap_unl(idx) + + def get_sim_qumap_unl(self, idx): + fn = os.path.join(self.lib_dir, 'sim_qumap_unl_%05d.npy'%idx) + if self.cache_sims and not os.path.exists(fn): + if pbs.rank == 0: + np.save(fn, np.array(self._build_sim_qumap_unl(idx))) + pbs.barrier() + if self.cache_sims: + return np.load(fn) + else: + return self._build_sim_qumap_unl(idx) + def get_sim_tmap(self, idx): fn = os.path.join(self.lib_dir, 'sim_tmap_%05d.npy'%idx) if self.cache_sims and not os.path.exists(fn): From 2bd9dcec2a7b1c1d78ae3211b58d4b272fd1181c Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Sat, 17 Jun 2023 20:37:27 +0200 Subject: [PATCH 17/23] Update requirements.txt --- requirements.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 86eb142..f5474be 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,3 +1,4 @@ numpy scipy -pyfftw \ No newline at end of file +pyfftw +finufftpy From 726fd65bac081b26b113759a95dd283b1950fae7 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Sat, 17 Jun 2023 20:46:33 +0200 Subject: [PATCH 18/23] Update requirements.txt --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index f5474be..9d09f37 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ numpy scipy pyfftw -finufftpy +finufft From 0cd7ce0615f2d9112b5510e78f9a09f0c9d364fc Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Fri, 23 Jun 2023 16:38:02 +0200 Subject: [PATCH 19/23] Update lensingmap.py --- lensit/clusterlens/lensingmap.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index c3d7ed4..ef51c94 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -84,8 +84,11 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co # self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} + # Set the high ell to zero escept for the lensing potential cpp_fid for cl in self.cls_unl.values(): cl[6000:] = 0 + self.cls_unl['pp'][6000:] = cpp_fid[6000:] + # Generate the CMB random phases skypha_libdir = opj(self.libdir, 'len_alms', 'skypha') From ff454848a9c0cfeca2be0ba5943128b45ef2176e Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Sun, 25 Jun 2023 13:12:52 +0200 Subject: [PATCH 20/23] fixed cpp_fid in lensingmap.py --- lensit/clusterlens/lensingmap.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index ef51c94..7f86581 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -85,9 +85,9 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} # Set the high ell to zero escept for the lensing potential cpp_fid - for cl in self.cls_unl.values(): - cl[6000:] = 0 - self.cls_unl['pp'][6000:] = cpp_fid[6000:] + for key, cl in zip(self.cls_unl.keys(), self.cls_unl.values()): + if key != 'pp': + cl[6000:] = 0 # Generate the CMB random phases From b3d3890d496ad53980fe02ea0c48213cff3171e8 Mon Sep 17 00:00:00 2001 From: Sayan Saha <46063942+s-Sayan@users.noreply.github.com> Date: Sat, 8 Jul 2023 12:07:45 +0200 Subject: [PATCH 21/23] Update lensingmap.py --- lensit/clusterlens/lensingmap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index 7f86581..c83fe5a 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -84,7 +84,7 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co # self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} self.cls_unl = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid} - # Set the high ell to zero escept for the lensing potential cpp_fid + # Set the high ell to zero except for the lensing potential cpp_fid for key, cl in zip(self.cls_unl.keys(), self.cls_unl.values()): if key != 'pp': cl[6000:] = 0 From a5d47543cbba526952f778d69ff1c324b4c96970 Mon Sep 17 00:00:00 2001 From: data_vadar Date: Tue, 25 Jul 2023 17:05:26 +0200 Subject: [PATCH 22/23] final updates --- lensit/__init__.py | 17 ++++++++++++++++- lensit/clusterlens/lensingmap.py | 6 +++--- lensit/clusterlens/profile.py | 4 ++-- lensit/ffs_covs/ffs_cov.py | 5 +++-- lensit/ffs_deflect/ffs_deflect.py | 2 +- lensit/ffs_iterators/ffs_iterator.py | 18 ++++++++++++------ lensit/ffs_iterators/ffs_iterator_nufft.py | 6 +++--- 7 files changed, 40 insertions(+), 18 deletions(-) diff --git a/lensit/__init__.py b/lensit/__init__.py index 859f1b6..a4b2e86 100644 --- a/lensit/__init__.py +++ b/lensit/__init__.py @@ -127,7 +127,7 @@ def get_maps_lib(exp, LDres, HDres=14, wrotation=False, cache_lenalms=True, cac len_cmbs = get_lencmbs_lib(res=HDres, cache_sims=cache_lenalms, nsims=nsims, wrotation=wrotation) lmax_sky = len_cmbs.lib_skyalm.ellmax cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=lmax_sky) - lib_datalm = ffs_covs.ell_mat.ffs_alm_pyFFTW(get_ellmat(LDres, HDres), filt_func=lambda ell: ell <= lmax_sky, + lib_datalm = ell_mat.ffs_alm_pyFFTW(get_ellmat(LDres, HDres), filt_func=lambda ell: ell <= lmax_sky, num_threads=num_threads) fsky = int(np.round(np.prod(len_cmbs.lib_skyalm.ell_mat.lsides) / 4. / np.pi * 1000.)) vcell_amin2 = np.prod(lib_datalm.ell_mat.lsides) / np.prod(lib_datalm.ell_mat.shape) * (180 * 60. / np.pi) ** 2 @@ -249,6 +249,21 @@ def get_config(exp): Beam_FWHM_amin = 1. ellmin=10 ellmax=3000 + elif exp == 'S4_sayan': + sN_uKamin = 1. + Beam_FWHM_amin = 1. + ellmin=100 + ellmax=4000 + elif exp == 'S4_no_noise': + sN_uKamin = 0. + Beam_FWHM_amin = 1. + ellmin=100 + ellmax=4000 + elif exp == 'CMB_HD': + sN_uKamin = 0.5 + Beam_FWHM_amin = 0.5 + ellmin=100 + ellmax=4000 elif exp == 'Planck45': sN_uKamin = 45. Beam_FWHM_amin = 5 diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index c83fe5a..4e21e94 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -88,8 +88,8 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co for key, cl in zip(self.cls_unl.keys(), self.cls_unl.values()): if key != 'pp': cl[6000:] = 0 - - + + # Generate the CMB random phases skypha_libdir = opj(self.libdir, 'len_alms', 'skypha') skypha = ffs_phas.ffs_lib_phas(skypha_libdir, nfields, self.lib_skyalm, nsims_max=nsims) @@ -379,4 +379,4 @@ def get_sim_qulm(self, idx): np.save(fname, np.array([Qlm, Ulm])) return np.load(fname) - + \ No newline at end of file diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py index 8cf326a..15d375c 100644 --- a/lensit/clusterlens/profile.py +++ b/lensit/clusterlens/profile.py @@ -286,8 +286,8 @@ def kappa_map(self, M200, z, shape, lsides, xmax=None): # Center of the cluster #x0 = shape[0]/2 #y0 = shape[1]/2 - x0 = 0 - y0 = 0 + x0 = 0. + y0 = 0. #X, Y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1])) if shape[0] % 2 == 0 and shape[1] % 2 == 0: diff --git a/lensit/ffs_covs/ffs_cov.py b/lensit/ffs_covs/ffs_cov.py index 82b2eb9..7f570fe 100644 --- a/lensit/ffs_covs/ffs_cov.py +++ b/lensit/ffs_covs/ffs_cov.py @@ -663,8 +663,9 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, ellmax_gradleg=None, * if ellmax_gradleg is not None: fl = np.ones(self.lib_skyalm.ellmax+1) fl[ellmax_gradleg+1:] = 0 - for _i in range(len(typ)): - lib_qlm.almxfl(clms[_i], fl, inplace=True) + print("Using grad_cut for T") + if 'T' in typ: + lib_qlm.almxfl(clms[0], fl, inplace=True) t.checkpoint(" get_qlms::mult with %s Pmat" % ({True: 'len', False: 'unl'}[use_cls_len])) diff --git a/lensit/ffs_deflect/ffs_deflect.py b/lensit/ffs_deflect/ffs_deflect.py index 222b16f..46e02b5 100644 --- a/lensit/ffs_deflect/ffs_deflect.py +++ b/lensit/ffs_deflect/ffs_deflect.py @@ -717,4 +717,4 @@ def displacement_fromolm(lib_plm, olm, **kwargs): def displacement_frompolm(lib_plm, plm, olm, **kwargs): return ffs_displacement(lib_plm.alm2map(plm * lib_plm.get_ikx() - olm * lib_plm.get_iky()), lib_plm.alm2map(plm * lib_plm.get_iky() + olm * lib_plm.get_ikx()), - lib_plm.ell_mat.lsides, **kwargs) + lib_plm.ell_mat.lsides, **kwargs) \ No newline at end of file diff --git a/lensit/ffs_iterators/ffs_iterator.py b/lensit/ffs_iterators/ffs_iterator.py index 8e92020..c547174 100644 --- a/lensit/ffs_iterators/ffs_iterator.py +++ b/lensit/ffs_iterators/ffs_iterator.py @@ -88,10 +88,10 @@ def __init__(self, lib_dir, typ, filt, dat_maps, lib_qlm, Plm0, H0, cpp_prior, def newton_step_length(it, norm_incr): # FIXME # Just trying if half the step is better for S4 QU - if filt.Nlev_uKamin('t') > 2.1: return 1.0 + if filt.Nlev_uKamin('t') > 2.1: return 0.5 if filt.Nlev_uKamin('t') <= 2.1 and norm_incr >= 0.5: - return 0.5 - return 0.5 + return 0.1 + return 0.1 self.newton_step_length = newton_step_length self.soltn0 = soltn0 @@ -662,9 +662,15 @@ def __init__(self, lib_dir, typ, filt, dat_maps, lib_qlm, Plm0, H0, cpp_prior, # 'u': (filt.Nlev_uKamin('u') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_sky_ivf + 1)} lmax_ivf = filt.lib_datalm.ellmax iso_libdat = filt.lib_datalm - cls_noise = {'t': (filt.Nlev_uKamin('t') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1), - 'q': (filt.Nlev_uKamin('q') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1), - 'u': (filt.Nlev_uKamin('u') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1)} + cl_t = (filt.Nlev_uKamin('t') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1) + cl_q = (filt.Nlev_uKamin('q') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1) + cl_u = (filt.Nlev_uKamin('u') / 60. / 180. * np.pi) ** 2 * np.ones(lmax_ivf + 1) + #cl_t[3000:] = 1e9 + #cl_q[3000:] = 1e9 + #cl_u[3000:] = 1e9 + cls_noise = {'t': cl_t, + 'q': cl_q, + 'u': cl_u} self.isocov = ffs_cov.ffs_diagcov_alm(os.path.join(lib_dir, 'isocov'), iso_libdat, filt.cls, filt.cls, filt.cl_transf, cls_noise, lib_skyalm=filt.lib_skyalm, init_rank=init_rank, diff --git a/lensit/ffs_iterators/ffs_iterator_nufft.py b/lensit/ffs_iterators/ffs_iterator_nufft.py index 0dc4e25..cd17298 100644 --- a/lensit/ffs_iterators/ffs_iterator_nufft.py +++ b/lensit/ffs_iterators/ffs_iterator_nufft.py @@ -91,10 +91,10 @@ def __init__(self, lib_dir, typ, filt: filtr, dat_maps, lib_qlm, Plm0, H0, cpp_p def newton_step_length(it, norm_incr): # FIXME # Just trying if half the step is better for S4 QU - if filt.Nlev_uKamin('t') > 2.1: return 1.0 + if filt.Nlev_uKamin('t') > 2.1: return 0.5 if filt.Nlev_uKamin('t') <= 2.1 and norm_incr >= 0.5: - return 0.5 - return 0.5 + return 0.1 + return 0.1 self.newton_step_length = newton_step_length self.soltn0 = soltn0 From f9482a4a135f22687dbddd5582f9114723bd2157 Mon Sep 17 00:00:00 2001 From: data_vadar Date: Wed, 26 Jul 2023 15:17:53 +0200 Subject: [PATCH 23/23] added kappa_0 function in lensingmap.py --- lensit/clusterlens/lensingmap.py | 45 +- lensit/clusterlens/profile.py | 7 +- lensit/misc/misc_utils.py | 12 + notebooks/demo_kappa0.ipynb | 1025 ++++++++++++++++++++++++++++++ 4 files changed, 1077 insertions(+), 12 deletions(-) create mode 100644 notebooks/demo_kappa0.ipynb diff --git a/lensit/clusterlens/lensingmap.py b/lensit/clusterlens/lensingmap.py index 4e21e94..e3d2d70 100644 --- a/lensit/clusterlens/lensingmap.py +++ b/lensit/clusterlens/lensingmap.py @@ -4,7 +4,7 @@ import numpy as np from lensit import ell_mat, pbs from lensit.sims import ffs_phas, ffs_cmbs, ffs_maps -from lensit.misc.misc_utils import enumerate_progress, gauss_beam +from lensit.misc.misc_utils import enumerate_progress, gauss_beam, pp_to_kk, p_to_k from lensit.ffs_deflect import ffs_deflect import lensit as li from os.path import join as opj @@ -39,7 +39,7 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co lpix_amin: physical size (in arcmin) of the pixels nsims: number of CMB maps to simulate cosmo: Instantiated camb.results.CAMBdata object - profparams: dict containing the parameters defining the profile + profparams: dict containing the parameters defining the profile M200, z, xmaxn profilename: string defining the density profile of the cluster (e.g. nfw) ellmax_sky: maximum multipole of CMB spectra used to generate the CMB maps """ @@ -50,8 +50,12 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co if profilename == 'nfw': self.M200 = profparams['M200c'] self.z = profparams['z'] + self.xmaxn = profparams['xmaxn'] self.npix = npix self.lpix_amin = lpix_amin + lbox_amin = npix*lpix_amin + lbox_rad = (lbox_amin/60)*(np.pi/180) + self.lbox_rad = lbox_rad shape = (self.npix, self.npix) lpix_rad = self.lpix_amin * np.pi / 180 / 60 lsides = (lpix_rad*self.npix, lpix_rad*self.npix) @@ -99,6 +103,8 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co pbs.barrier() self.haloprofile = profile(self.cosmo, profilename) + self.xmax = self.xmaxn * self.haloprofile.get_concentration(self.M200, self.z) + self.kappa0 = self.haloprofile.get_kappa0(self.M200, self.z, self.xmax) lencmbs_libdir = opj(self.libdir, 'len_alms') @@ -111,9 +117,9 @@ def __init__(self, libdir:str, key:str, npix:int, lpix_amin:float, nsims:int, co if self.key == 'lss': self.len_cmbs = sims_cmb_len(lencmbs_libdir, self.lib_skyalm, self.cls_unl, lib_pha=skypha) elif self.key == 'cluster': - self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + self.len_cmbs = sim_cmb_len_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.xmaxn, self.haloprofile, lib_pha=skypha) elif self.key == 'lss_plus_cluster': - self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.haloprofile, lib_pha=skypha) + self.len_cmbs = sim_cmb_len_lss_plus_cluster(lencmbs_libdir, self.lib_skyalm, self.cls_unl, self.M200, self.z, self.xmaxn, self.haloprofile, lib_pha=skypha) else: assert 0 @@ -222,6 +228,29 @@ def get_kappa_map(self, M200, z, xmax=None): convergence map: numpy array of shape self.lib_skyalm.shape """ return self.haloprofile.kappa_map(M200, z, self.lib_skyalm.shape, self.lib_skyalm.lsides, xmax) + + def get_kappa0_from_sim(self, lmin, lmax, phi_obs_lm, NL): + """ Get the kappa_0 estimate from one simulation + Args: + lmin, lmax: the lrange for the sum in the numerator and denominator of the estimator + phi_obs_lm: phi (lensing potential) estimate of a single simulation/observation + NL: reconstruction niose of the phi estimate + Returns: + kappa_0 estiamte for the single phi simulation/observation + + P.S. the arguments is taken as the phi (not kappa) of the observation and accordingly the noise + """ + ell = np.where(self.lib_skyalm.get_Nell()[:lmax+1] > 1)[0] + ell = ell[ell >= lmin] + #kappa_l = kappa_ell_lensit(ell)*lbox_rad + kappa_l = self.haloprofile.analitic_kappa_ft(self.M200, self.z, self.xmax, ell) + num_l = self.lib_skyalm.get_Nell()[ell] + phi_obs_el = self.lib_skyalm.bin_realpart_inell(phi_obs_lm) #the average + denom_int = (num_l)*(kappa_l/self.lbox_rad/self.kappa0)**2 / (NL[ell]*pp_to_kk(ell)) + num_int = (num_l)*(kappa_l/self.lbox_rad/self.kappa0) * phi_obs_el[ell]*p_to_k(ell) / (NL[ell]*pp_to_kk(ell)) ## since the observed quantity is phi and we need kappa + denom1 = np.sum(denom_int) + numer1 = np.sum(num_int) + return numer1/denom1 verbose = False @@ -237,11 +266,11 @@ def get_fields(cls): class sim_cmb_len_cluster(ffs_cmbs.sims_cmb_len): - def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): + def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, xmaxn:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): super(sim_cmb_len_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) self.conc_par = haloprofile.get_concentration(M200,z) - self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= 3 * self.conc_par) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= xmaxn * self.conc_par) #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) #self.kappa_map = 1e2*haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax=self.conc_par) self.defl_map = haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) @@ -253,14 +282,14 @@ def _get_f(self, idx=None): return ffs_deflect.ffs_displacement(self.defl_map[0], self.defl_map[1], self.lib_skyalm.lsides) class sim_cmb_len_lss_plus_cluster: - def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): + def __init__(self, lib_dir:str, lib_skyalm:ell_mat.ffs_alm, cls_unl:dict, M200:float, z:float, xmaxn:float, haloprofile:profile, lib_pha=None, use_Pool=0, cache_lens=False): #super(sim_cmb_len_lss_plus_cluster, self).__init__(lib_dir, lib_skyalm, cls_unl, lib_pha=lib_pha, use_Pool=0, cache_lens=False) self.M200 = M200 self.z = z self.haloprofile = haloprofile self.conc_par = self.haloprofile.get_concentration(M200,z) - self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= 3 * self.conc_par) + self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides, xmax= xmaxn * self.conc_par) #self.kappa_map = haloprofile.kappa_map(M200, z, lib_skyalm.shape, lib_skyalm.lsides) self.defl_map = self.haloprofile.kmap2deflmap(self.kappa_map, lib_skyalm.shape, lib_skyalm.lsides) diff --git a/lensit/clusterlens/profile.py b/lensit/clusterlens/profile.py index 15d375c..cbeb360 100644 --- a/lensit/clusterlens/profile.py +++ b/lensit/clusterlens/profile.py @@ -70,11 +70,10 @@ def sigma_crit(self, z): Dang_LS = self.cosmo.angular_diameter_distance2(z, self.zcmb) return const.c_Mpcs**2 /(4*np.pi*const.G_Mpc3_pMsol_ps2) * Dang_OS/Dang_OL/Dang_LS - def get_kappa0(self, M200, z): + def get_kappa0(self, M200, z, xmax=None): """Return the value of the lensing convergence profile for x = 1, i.e. R = rs""" - rs = self.get_rs(M200, z) - rhos = self.get_rho_s(M200, z) - return 2*rhos*rs/3 / self.sigma_crit(z) + thsc = self.get_thetas_amin(M200,z) + return self.kappa_theta(M200, z, thsc, xmax=xmax)[0] def sigma_nfw(self, M200, z, R, xmax=None): """Analytic expression for the surface mass desinty of a NFW profile diff --git a/lensit/misc/misc_utils.py b/lensit/misc/misc_utils.py index f4be233..f63a6bf 100644 --- a/lensit/misc/misc_utils.py +++ b/lensit/misc/misc_utils.py @@ -92,11 +92,23 @@ def pp_to_kk(ls): """ return ls ** 2 * (ls+1) ** 2 * 0.25 +def p_to_k(ls): + """ Converts lensing potential el into lensing convergence el. + :math:`C_\ell^{\kappa \kappa} = (\ell (\ell +1))^2 /4 C_\ell^{\phi \phi}` + + """ + return ls * (ls+1.) * 0.5 + def kk_to_pp(ls): """Converts lensing convergence power spectrum into lensing potential power spectrum. """ return cl_inverse(pp_to_kk(ls)) +def k_to_p(ls): + """Converts lensing convergence el into lensing potential el. + """ + return cl_inverse(p_to_k(ls)) + class timer: def __init__(self, verbose, prefix='', suffix=''): self.t0 = time.time() diff --git a/notebooks/demo_kappa0.ipynb b/notebooks/demo_kappa0.ipynb new file mode 100644 index 0000000..9bc1848 --- /dev/null +++ b/notebooks/demo_kappa0.ipynb @@ -0,0 +1,1025 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "import numpy as np \n", + "import camb \n", + "from numpy.fft import fftshift\n", + "from tqdm import tqdm\n", + "\n", + "import lensit as li\n", + "from lensit.clusterlens import lensingmap, profile \n", + "from lensit.misc.misc_utils import gauss_beam, pp_to_kk, p_to_k, cl_inverse\n", + "from lensit.ffs_covs import ffs_cov, ell_mat\n", + "from plancklens.wigners import wigners\n", + "from plancklens import n0s, nhl\n", + "from plancklens.n1 import n1\n", + "\n", + "import os\n", + "import os.path as op\n", + "import matplotlib as mpl\n", + "from matplotlib import pyplot as plt\n", + "\n", + "from scipy.interpolate import UnivariateSpline as spline" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# We start by intiating CAMB which will give us the relevant cosmology \n", + "cambinifile = 'planck_2018_acc'\n", + "\n", + "pars = camb.read_ini(op.join(op.dirname('/Users/sayan/CMB_WORK/CAMB-1.1.3/inifiles'), 'inifiles', cambinifile + '.ini'))\n", + "results = camb.get_results(pars)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The size of the data patch is 307.2 X 307.2 arcmin central box\n", + "/Users/sayan/Project_Geneva/modules/lensit_sims/temp/clustermaps/camb_planck_2018_acc/cmbexp_S4_sayan/nfw_profile/lensed_by_cluster/npix1024_lpix_0.3amin_lmaxsky6002/M200_2.000000E+14_z0.7/1000_sims\n" + ] + } + ], + "source": [ + "# We define here the parameters for the profile of the cluster\n", + "M200, z, xmaxn = 2 * 1e14, 0.7, 3\n", + "profname = 'nfw'\n", + "key = \"cluster\" # \"lss\"/\"cluster\"/\"lss_plus_cluster\"\n", + "profparams={'M200c':M200, 'z':z, 'xmaxn': xmaxn}\n", + "hprofile = profile.profile(results, profname)\n", + "xmax = xmaxn*hprofile.get_concentration(M200, z)\n", + "# Define here the map square patches\n", + "npix = 1024 # Number of pixels\n", + "lpix_amin = 0.3 # Physical size of a pixel in arcmin (There is bug when <0.2 amin, due to low precision in Cl_TE at )\n", + "\n", + "print(\"The size of the data patch is %0.1f X %0.1f arcmin central box\"%(npix*lpix_amin, npix*lpix_amin))\n", + "\n", + "# Maximum multipole used to generate the CMB maps from the CMB power spectra\n", + "# ellmaxsky = 6000 # (bug when ellmax>6300 because of low precision in Cl_TE of CAMB )\n", + "ellmaxsky = 6002 \n", + "\n", + "# Set the maximum ell observed in the CMB data maps\n", + "ellmaxdat = 4000\n", + "ellmindat = 100\n", + "\n", + "# Number of simulated maps \n", + "nsims = 1000\n", + "\n", + "# Set CMB experiment for noise level and beam\n", + "cmb_exp='S4_sayan'\n", + "\n", + "# We will cache things in this directory \n", + "libdir = lensingmap.get_cluster_libdir(cambinifile, profname, key, npix, lpix_amin, ellmaxsky, M200, z, nsims, cmb_exp)\n", + "#libdir = op.join(libdir,\"trunc\")\n", + "print(libdir)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "lmax = ellmaxsky\n", + "cpp_fid = results.get_lens_potential_cls(lmax=lmax, raw_cl=True).T[0]\n", + "\n", + "camb_cls = results.get_unlensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=lmax).T\n", + "cls_unl_fid = {'tt':camb_cls[0], 'ee':camb_cls[1], 'bb':camb_cls[2], 'te':camb_cls[3], 'pp':cpp_fid}\n", + "\n", + "camb_cls_len = results.get_lensed_scalar_cls(CMB_unit='muK', raw_cl=True, lmax=lmax).T\n", + "cls_len_fid = {'tt':camb_cls_len[0], 'ee':camb_cls_len[1], 'bb':camb_cls_len[2], 'te':camb_cls_len[3], 'pp':cpp_fid}" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/sayan/Project_Geneva/git_repos/LensIt/lensit/clusterlens/profile.py:320: RuntimeWarning: invalid value encountered in true_divide\n", + " dx_lm = 2 * rfft_kappa * 1.j * KX / (KX**2+KY**2)\n", + "/Users/sayan/Project_Geneva/git_repos/LensIt/lensit/clusterlens/profile.py:321: RuntimeWarning: invalid value encountered in true_divide\n", + " dy_lm = 2 * rfft_kappa * 1.j * KY / (KX**2+KY**2)\n" + ] + } + ], + "source": [ + "np.random.seed(seed=20)\n", + "clustermaps = lensingmap.cluster_maps(libdir, key, npix, lpix_amin, nsims, results, profparams, profilename=profname, ellmax_sky = ellmaxsky, ellmax_data=ellmaxdat, ellmin_data=ellmindat, cmb_exp=cmb_exp)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "ellmax_sky = clustermaps.ellmax_sky\n", + "sN_uKamin, sN_uKaminP, Beam_FWHM_amin, ellmin, ellmax = li.get_config(clustermaps.cmb_exp)\n", + "\n", + "cls_noise = {'t': (sN_uKamin * np.pi / 180. / 60.) ** 2 * np.ones(clustermaps.ellmax_sky + 1),\n", + " 'q':(sN_uKaminP * np.pi / 180. / 60.) ** 2 * np.ones(clustermaps.ellmax_sky + 1),\n", + " 'u':(sN_uKaminP * np.pi / 180. / 60.) ** 2 * np.ones(clustermaps.ellmax_sky + 1)} # simple flat noise Cls\n", + "# cl_transf = gauss_beam(Beam_FWHM_amin / 60. * np.pi / 180., lmax=ellmax_sky)\n", + "# lib_alm = ell_mat.ffs_alm_pyFFTW(get_ellmat(LD_res, HD_res=HD_res),\n", + " # filt_func=lambda ell: (ell >= ellmin) & (ell <= ellmax), num_threads=pyFFTWthreads)\n", + "# lib_skyalm = ell_mat.ffs_alm_pyFFTW(clustermaps.ellmat,\n", + " # filt_func=lambda ell: (ell <= ellmax_sky), num_threads=clustermaps.num_threads)\n", + "\n", + "cl_transf = clustermaps.cl_transf\n", + "lib_skyalm = clustermaps.lib_skyalm\n", + "\n", + "typ = 'T'\n", + "\n", + "lib_dir = op.join(clustermaps.dat_libdir, typ)\n", + "# isocov = ffs_cov.ffs_diagcov_alm(lib_dir, clustermaps.lib_datalm, clustermaps.cls_unl, cls_len, cl_transf, cls_noise, lib_skyalm=lib_skyalm)\n", + "isocov = ffs_cov.ffs_diagcov_alm(lib_dir, clustermaps.lib_datalm, cls_unl_fid, cls_unl_fid, cl_transf, cls_noise, lib_skyalm=lib_skyalm)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " lensit.ffs_covs.ffs_cov [00:00:00] (total [00:00:00]) inverse len Pmats curvpOlm \n", + " lensit.ffs_covs.ffs_cov [00:00:01] (total [00:00:01]) Fxx , part 1 curvpOlm \n", + " lensit.ffs_covs.ffs_cov [00:00:00] (total [00:00:01]) Fyy , part 1 curvpOlm \n", + " lensit.ffs_covs.ffs_cov [00:00:00] (total [00:00:01]) Fxy , part 1 curvpOlm \n", + " lensit.ffs_covs.ffs_cov [00:00:00] (total [00:00:01]) Fxx , part 2 curvpOlm \n", + " lensit.ffs_covs.ffs_cov [00:00:00] (total [00:00:01]) Fyy , part 2 curvpOlm \n", + " lensit.ffs_covs.ffs_cov [00:00:00] (total [00:00:01]) Fxy , part 2 curvpOlm \n" + ] + } + ], + "source": [ + "H0len = isocov.get_response(typ, lib_skyalm, use_cls_len=True)[0]\n", + "def get_starting_point(idx, typ, clustermaps, keyword=None): \n", + " \"\"\"\n", + " This returns initial data for simulation index 'idx' from a CMB-S4 simulation library.\n", + " On first call the simulation library will generate all simulations phases, hence might take a little while.\n", + " \"\"\" \n", + "\n", + " print(\" I will be using data from ell=%s to ell=%s only\"%(isocov.lib_datalm.ellmin, isocov.lib_datalm.ellmax))\n", + " print(\" The sky band-limit is ell=%s\"%(isocov.lib_skyalm.ellmax))\n", + " \n", + " lib_qlm = lib_skyalm #: This means we will reconstruct the lensing potential for all unlensed sky modes.\n", + " ellmax_sky = lib_skyalm.ellmax\n", + " ell = np.arange(ellmax_sky+1)\n", + "\n", + " if typ=='QU':\n", + " datalms = np.array([isocov.lib_datalm.map2alm(m) for m in clustermaps.maps_lib.get_sim_qumap(idx)]) \n", + " elif typ =='T':\n", + " datalms = np.array([isocov.lib_datalm.map2alm(clustermaps.maps_lib.get_sim_tmap(idx))]) \n", + " elif typ =='TQU':\n", + " datalms = np.array([isocov.lib_datalm.map2alm(m) for m in np.array([clustermaps.maps_lib.get_sim_tmap(idx), clustermaps.maps_lib.get_sim_qumap(idx)[0], clustermaps.maps_lib.get_sim_qumap(idx)[1]])]) \n", + " \n", + " use_cls_len = True\n", + " \n", + " plm1 = 0.5 * isocov.get_qlms(typ, isocov.get_iblms(typ, datalms, use_cls_len=use_cls_len)[0], lib_qlm, \n", + " use_cls_len=use_cls_len)[0]\n", + " \n", + " plmqe = lib_skyalm.almxfl(plm1, cl_inverse(H0len), inplace=False)\n", + "\n", + " return plmqe" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "1\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "2\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "3\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "4\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "5\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "6\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "7\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "8\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "9\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "10\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "11\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "12\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "13\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "14\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "15\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "16\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "17\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "18\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "19\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "20\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "21\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "22\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "23\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "24\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "25\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "26\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "27\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "28\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "29\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "30\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "31\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "32\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "33\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "34\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "35\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "36\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "37\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "38\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "39\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "40\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "41\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "42\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "43\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "44\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "45\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "46\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "47\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "48\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "49\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "50\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "51\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "52\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "53\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "54\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "55\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "56\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "57\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "58\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "59\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "60\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "61\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "62\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "63\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "64\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "65\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "66\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "67\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "68\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "69\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "70\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "71\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "72\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "73\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "74\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "75\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "76\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "77\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "78\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "79\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "80\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "81\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "82\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "83\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "84\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "85\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "86\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "87\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "88\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "89\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "90\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "91\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "92\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "93\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "94\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "95\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "96\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "97\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "98\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n", + "99\n", + " I will be using data from ell=140 to ell=3999 only\n", + " The sky band-limit is ell=6002\n", + "Building nufft fwd plan\n", + " [00:00:00] (total [00:00:00]) get_qlms::mult with len Pmat \n", + " [00:00:00] (total [00:00:00]) get_qlms::cartesian gradients \n", + " [00:00:00] (total [00:00:00]) get_qlms::rotation to phi Omega \n" + ] + } + ], + "source": [ + "nmaps = 100\n", + "plmqes = [None]*nmaps\n", + "\n", + "if nsims >1:\n", + " for idx in range(nmaps):\n", + " print(idx)\n", + " plmqes[idx] = get_starting_point(idx, typ, clustermaps)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 100/100 [00:00<00:00, 133.07it/s]\n" + ] + } + ], + "source": [ + "auto_corr_1 = nmaps*[None]\n", + "for idx in tqdm(range(nmaps)):\n", + " auto_corr_1[idx] = lib_skyalm.alm2cl(plmqes[idx])\n", + "\n", + "phi_var_1 = np.mean(auto_corr_1, axis=0)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "100%|██████████| 100/100 [00:01<00:00, 88.48it/s]\n" + ] + } + ], + "source": [ + "kappa_0_sims_QE_T_1 = [None]*nmaps\n", + "lmax, lmin= 5000, ellmindat\n", + "ell = np.arange(ellmaxsky+1)\n", + "for idx in tqdm(range(nmaps)):\n", + " kappa_0_sims_QE_T_1[idx] = clustermaps.get_kappa0_from_sim(lmin=lmin, lmax=lmax, phi_obs_lm=plmqes[idx], NL=phi_var_1)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "========================QE Results========================\n", + "input kappa_0 0.128488\n", + "mean kappa 0 from 100 simulations is 0.120446\n", + "error from 100 simulations 0.019439\n" + ] + } + ], + "source": [ + "print(\"========================QE Results========================\")\n", + "print(\"input kappa_0 %f\"%(hprofile.get_kappa0(M200, z, xmax)))\n", + "print(\"mean kappa 0 from %i simulations is %f\"%(nmaps,np.mean(kappa_0_sims_QE_T_1)))\n", + "print(\"error from %i simulations %f\"%(nmaps,np.std(kappa_0_sims_QE_T_1)/np.sqrt(nmaps)))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.9.13 ('test')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.13" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "2a0a0d92d61eb81c1bfacf4848f33179e421b61b01830eeaa00bb11be584f4b6" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}