Clusterlens

README.md

@@ -38,4 +38,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)
+![SNSF logo](./docs/SNF_logo_standard_web_color_neg_e.svg)

lensit/__init__.py

@@ -358,6 +358,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

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):
@@ -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 
@@ -101,26 +100,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(x<self.theta_amin_to_x(2*1e14, 0.7, 0.12))
+        idx = self.theta_amin_to_x(2*1e14, 0.7, 0.12)
         f[xp] = (1 - 2/np.sqrt(x[xp]**2 - 1) * np.arctan(np.sqrt((x[xp] - 1)/(x[xp] + 1))))/(x[xp]**2-1)
         f[xm] = (1 - 2/np.sqrt(1 - x[xm]**2) * np.arctanh(np.sqrt((1 - x[xm])/(1 + x[xm]))))/(x[xm]**2-1)
         f[xo] = 1/3
+        f[xd] = (1 - 2/np.sqrt(1 - idx**2) * np.arctanh(np.sqrt((1 - idx)/(1 + idx))))/(idx**2-1)
         return f
 
 
     def gx(self, x, xmax, tol=5):
         """We apply a cutoff in the halo profile for x>xmax, 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, x<xmax))
+        xp = np.logical_and(x>1, 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<self.theta_amin_to_x(2*1e14, 0.7, 0.12))
+        idx = self.theta_amin_to_x(2*1e14, 0.7, 0.12)
         g[xp] = (np.sqrt(xmax**2 - x[xp]**2)/(1+xmax) - 1/np.sqrt(x[xp]**2 - 1) * np.arccos((x[xp]**2 + xmax) / (x[xp] * (1+xmax))) )/(x[xp]**2-1)
         g[xm] = (np.sqrt(xmax**2 - x[xm]**2)/(1+xmax) - 1/np.sqrt(1-x[xm]**2) * np.arccosh((x[xm]**2 + xmax) / (x[xm] * (1+xmax))) )/(x[xm]**2-1)
         g[xo] = np.sqrt(xmax**2 - 1)/(3*(1+xmax)) * (1+ 1/ (1+xmax))
+        g[xd] = (np.sqrt(xmax**2 - idx**2)/(1+xmax) - 1/np.sqrt(1-idx**2) * np.arccosh((idx**2 + xmax) / (idx * (1+xmax))) )/(idx**2-1)
         return g
 
 
@@ -166,20 +175,20 @@ 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, xmax, ell, const_c=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)
-        mu_nfw = np.log(1. + c) - c / (1. + c)
+        c1 = self.get_concentration(M200, z, const_c)
+        mu_nfw = np.log(1. + c1) - c1 / (1. + c1)
         rs = self.get_rs(M200, z)
         chi = self.cosmo.comoving_radial_distance(z)
         k = ell / chi
         x = ((1. + z) * k * rs)
         Six, Cix = sici(x)
-        Sixpc, Cixpc = sici(x * (1. + c))
+        Sixpc, Cixpc = sici(x * (1. + xmax))
         Sidiff = Sixpc - Six
         Cidiff = Cixpc - Cix
-        u0 = np.sin(x) * Sidiff + np.cos(x) * Cidiff - np.sin(x * c) / (x * (1. + c))
+        u0 = np.sin(x) * Sidiff + np.cos(x) * Cidiff - np.sin(x * xmax) / (x * (1. + xmax))
         ufft = 1. / mu_nfw * u0
 
         kappaft = M200 * ufft * (1+z)**2 /chi**2 / self.sigma_crit(z)
@@ -224,6 +233,18 @@ def theta_amin_to_r(self, z, theta_amin):
         R = Dang * theta_amin * np.pi/180/60
         return R
 
+    def theta_amin_to_x(self, M200, z, th_amin, const_c=None):
+        """Angle substended at chararcteric scale x = r/rs
+
+            Args:
+                M200: cluster M200 mass in solar masses (?)
+                z: cluster redshift
+                x:  dimensionless R / Rs
+
+        """
+        return th_amin/self.get_thetas_amin(M200, z, const_c=const_c)
+
+
 
     def x_to_theta_amin(self,M200, z, x, const_c=None):
         """Angle substended at chararcteric scale x = r/rs
@@ -262,9 +283,18 @@ def kappa_map(self, M200, z, shape, lsides, xmax=None):
         dtheta_x = lsides[0]/shape[0] * 180/np.pi*60
         dtheta_y = lsides[1]/shape[1] * 180/np.pi*60
         # Center of the cluster
-        x0 = (shape[0]+1.)/2. 
-        y0 = (shape[1]+1.)/2
-        X, Y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]))
+        #x0 = shape[0]/2
+        #y0 = shape[1]/2       
+        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:
+            X, Y =  np.meshgrid(np.concatenate((np.arange(0,shape[0]//2), np.arange(-shape[0]//2,0))), np.concatenate((np.arange(0,shape[1]//2), np.arange(-shape[1]//2,0))))
+        else:
+            assert 0, "pixel size is not right"
+
+        #X, Y = np.meshgrid(np.arange(shape[0]), np.arange(shape[1]))
         theta_amin = self.pix_to_theta(X, Y, (dtheta_x,dtheta_y),  (x0, y0))
 
         return self.kappa_theta(M200, z, theta_amin, xmax=xmax)

lensit/ffs_covs/ffs_cov.py

@@ -672,7 +672,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, use_cls_grad=False, iblms2=None, **kwargs):
+    def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, iblms2=None,  ellmax_gradleg=None, **kwargs):
         r"""Unormalized quadratic estimates (potential and curl).
 
         Note:
@@ -688,6 +688,7 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, ib
             use_cls_len: use lensed or unlensed cls in QE weights (numerator), defaults to lensed cls
             use_cls_grad: use grad cls in QE weights (numerator), defaults to use_cls_len
             iblms2: second leg inverse variance filtered CMB maps (useful for RDN0 estimates)
+            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))
@@ -714,6 +715,13 @@ def get_qlms(self, typ, iblms, lib_qlm, use_cls_len=True, use_cls_grad=False, ib
                 for _j in range(len(typ)):
                     clms2[_i] += pmat.get_unlPmat_ij(typ, self.lib_skyalm, weights_cls, _i, _j) * iblms2[_j]
 
+        if ellmax_gradleg is not None:
+            fl = np.ones(self.lib_skyalm.ellmax+1)
+            fl[ellmax_gradleg+1:] = 0
+            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]))
 
         _map = lambda alm: self.lib_skyalm.alm2map(alm)

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,

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

lensit/misc/misc_utils.py

@@ -119,11 +119,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()

lensit/sims/ffs_maps.py

@@ -233,6 +233,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))
@@ -255,6 +268,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):