From 124221226f2ab7cd794451460f5d9e46083ffc9d Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Wed, 21 Jun 2023 10:09:23 +0200 Subject: [PATCH 1/7] MAP aniso with simple source hardening in the noise filtering, by allowing a full covariance with a simple basic approximation --- delensalot/core/opfilt/MAP_opfilt_aniso_t.py | 29 +++++++++++++++++++- 1 file changed, 28 insertions(+), 1 deletion(-) diff --git a/delensalot/core/opfilt/MAP_opfilt_aniso_t.py b/delensalot/core/opfilt/MAP_opfilt_aniso_t.py index 2d931627..150968c4 100644 --- a/delensalot/core/opfilt/MAP_opfilt_aniso_t.py +++ b/delensalot/core/opfilt/MAP_opfilt_aniso_t.py @@ -28,7 +28,8 @@ def apply_fini(*args, **kwargs): class alm_filter_ninv_wl(opfilt_base.alm_filter_wl): def __init__(self, ninv_geom:utils_geom.Geom, ninv: np.ndarray, ffi:remapping.deflection, transf:np.ndarray, - unlalm_info:tuple, lenalm_info:tuple, sht_threads:int,verbose=False, lmin_dotop=0, tpl:tni.template_tfilt or None =None): + unlalm_info:tuple, lenalm_info:tuple, sht_threads:int,verbose=False, lmin_dotop=0, tpl:tni.template_tfilt or None =None, + extra_ninv: np.ndarray = None): r"""CMB inverse-variance and Wiener filtering instance, using unlensed E and lensing deflection Args: @@ -40,6 +41,7 @@ def __init__(self, ninv_geom:utils_geom.Geom, ninv: np.ndarray, ffi:remapping.de lenalm_info: tuple of int, lmax and mmax of lensed CMB sht_threads: number of threads for lenspyx SHTs verbose: some printout if set, defaults to False + extra_ninv: np.ndarray, extra map for inverse noise filtering to be applied to the observed map """ lmax_unl, mmax_unl = unlalm_info @@ -75,6 +77,9 @@ def __init__(self, ninv_geom:utils_geom.Geom, ninv: np.ndarray, ffi:remapping.de self.template = tpl + self.extra_ninv = extra_ninv # I do not use this one for the diagonal pre-conditioned CG + ## assume extra_inv will enter in the form: B^t extra_ninv B, then I want to take the inverse of this + def hashdict(self): return {'ninv':self._ninv_hash(), 'transf':clhash(self.b_transf_tlm), 'deflection':self.ffi.hashdict(), @@ -101,6 +106,28 @@ def apply_map(self, tmap): """ tmap *= self.n_inv + + if self.extra_ninv is not None: + """ + Here I am doing this: + N^{-1} - N^{-1}B^t \hat{S}^2 B N^{-1}, applied on X map + + This assumes that we have a small perturbation to the noise. + + Note that for now here the noise is diagonal in pixel space. This is not the case for the extra piece of here. + """ + B = tmap.copy() + B = self.geom_.map2alm(B.copy(), self.lmax_len, self.mmax_len, self.ffi.sht_tr, (-1., 1.)) + B = almxfl(B, self.b_transf_tlm, self.mmax_len, inplace = False) + B = self.geom_.alm2map(B.copy(), self.lmax_len, self.mmax_len, self.ffi.sht_tr, (-1., 1.)) + B *= self.extra_ninv + B = self.geom_.map2alm(B.copy(), self.lmax_len, self.mmax_len, self.ffi.sht_tr, (-1., 1.)) + B = almxfl(B, self.b_transf_tlm, self.mmax_len, inplace = False) + B = self.geom_.alm2map(B.copy(), self.lmax_len, self.mmax_len, self.ffi.sht_tr, (-1., 1.)) + B *= self.n_inv + tmap -= B + + if self.template is not None: ts = [self.template] # Hack, this is only meant for one template coeffs = np.concatenate(([t.dot(tmap) for t in ts])) From 9784b358e014e4b8b46d9bff4861ce583c684a6e Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Wed, 21 Jun 2023 10:17:28 +0200 Subject: [PATCH 2/7] generic sims that allows the inclusion of extra maps at the cmb level, as well as a sims_general module to try something enw --- delensalot/sims/generic.py | 23 +++++++++- delensalot/sims/sims_general.py | 75 +++++++++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 2 deletions(-) create mode 100644 delensalot/sims/sims_general.py diff --git a/delensalot/sims/generic.py b/delensalot/sims/generic.py index c57e689a..f7829a20 100644 --- a/delensalot/sims/generic.py +++ b/delensalot/sims/generic.py @@ -51,9 +51,11 @@ class sims_cmb_len(object): verbose(defaults to True): lenspyx timing info printout + extra_tlm: optional extra map to add to the CMB map, e.g. foregrounds + """ def __init__(self, lib_dir, lmax, cls_unl, lib_pha=None, offsets_plm=None, offsets_cmbunl=None, - dlmax=1024, nside_lens=4096, epsilon=1e-7, nbands=8, cache_plm=True, verbose=True): + dlmax=1024, nside_lens=4096, epsilon=1e-7, nbands=8, cache_plm=True, verbose=True, extra_tlm = None): fields = _get_fields(cls_unl) @@ -99,6 +101,8 @@ def __init__(self, lib_dir, lmax, cls_unl, lib_pha=None, offsets_plm=None, offse self.lens_module = lenspyx self.verbose=verbose + self.extra_tlm = extra_tlm + @staticmethod def offset_index(idx, block_size, offset): """Offset index by amount 'offset' cyclically within blocks of size block_size @@ -169,17 +173,32 @@ def _cache_eblm(self, idx): def get_sim_tlm(self, idx, ret=True): fname = os.path.join(self.lib_dir, 'sim_%04d_tlm.fits' % idx) + pfname = os.path.join(self.lib_dir, 'sim_%04d_plm.fits' % idx) if not os.path.exists(fname): tlm = self.unlcmbs.get_sim_tlm(self.offset_index(idx, self.offset_cmb[0], self.offset_cmb[1])) dlm = self.get_sim_plm(idx) + + hp.write_alm(pfname, dlm) + assert 'o' not in self.fields, 'not implemented' lmaxd = hp.Alm.getlmax(dlm.size) hp.almxfl(dlm, np.sqrt(np.arange(lmaxd + 1, dtype=float) * np.arange(1, lmaxd + 2)), inplace=True) Tlen = self.lens_module.alm2lenmap(np.array(tlm), [dlm, None], epsilon=self.epsilon, verbose=self.verbose) hp.write_alm(fname, hp.map2alm(Tlen, lmax=self.lmax, iter=0)) + + if (self.extra_tlm is not None): + """ + Adding an extra CMB component to the lensed CMB. + """ + extrafname = os.path.join(self.lib_dir, f'sim_{idx:04}_{self.extra_tlm.get_name()}lm.fits') + if (not os.path.exists(extrafname)): + extra_tlm = hp.map2alm(self.extra_tlm(idx), lmax=self.lmax, iter=0) + hp.write_alm(extrafname, extra_tlm) + if ret: - return hp.read_alm(fname) + total = hp.read_alm(fname) + (0 if self.extra_tlm is None else hp.read_alm(extrafname)) + return total def get_sim_elm(self, idx, ret=True): fname = os.path.join(self.lib_dir, 'sim_%04d_elm.fits' % idx) diff --git a/delensalot/sims/sims_general.py b/delensalot/sims/sims_general.py new file mode 100644 index 00000000..398af364 --- /dev/null +++ b/delensalot/sims/sims_general.py @@ -0,0 +1,75 @@ +"""Module to allow general simulations that include cmb+foregrounds+other at the signal level. +""" + +from plancklens.sims import maps + +import healpy as hp + + +class cmb_maps(maps.cmb_maps): + + """ + Class to handle multiple objects that give a field at the CMB signal level. + """ + + def __init__(self, **kwargs): + """ + Initializes the cmb_maps object. Note, for now you have to initialize it with a sims_cmb_len object. + """ + super(cmb_maps, self).__init__(**kwargs) + self.components = [] + self.lmax = self.sims_cmb_len.lmax + + + def get_sim_tmap(self,idx): + """Returns temperature healpy map for a simulation + + Args: + idx: simulation index + + Returns: + healpy map + """ + tmap = self.get_sim_tlm(idx) + hp.almxfl(tmap,self.cl_transf,inplace=True) + tmap = hp.alm2map(tmap,self.nside) + return tmap + self.get_sim_tnoise(idx) + + def get_sim_tlm(self, idx): + """ + Returns the temperature alm of the sum of the components. + + Args: + idx: simulation index + + Returns: + alm of the sum of the components. + """ + return sum([self.check(c.get_sim_tlm(idx)) for c in self.components]) + + def check(self, alms): + """ + Just a simple check to make sure the alms have the same lmax as the one initialized the object. + + Args: + alms: alms to check + + Returns: + alms + """ + assert hp.Alm.getlmax(alms.size) == self.lmax, "The alms you are trying to add have a different lmax than the one you initialized the object with!" + return alms + + def __add__(self, other): + """ + Adds a component to the cmb_maps object. + """ + return self._update(other) + + def _update(self, other): + """ + Adds a component to the cmb_maps object. + """ + assert callable(other.get_sim_tlm), "The object you are trying to add has to have a get_sim_tlm method!" + self.components.append(other) + return self \ No newline at end of file From 85365f6cd9658014ba094255887f2d06e0229644 Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Wed, 21 Jun 2023 10:19:15 +0200 Subject: [PATCH 3/7] sims extra component --- delensalot/sims/sims_extra.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 delensalot/sims/sims_extra.py diff --git a/delensalot/sims/sims_extra.py b/delensalot/sims/sims_extra.py new file mode 100644 index 00000000..e91add59 --- /dev/null +++ b/delensalot/sims/sims_extra.py @@ -0,0 +1,25 @@ +""" +Extra sims utility. +""" + +import os +import healpy as hp +import numpy as np + + +class Extra(object): + """ + Example: extra_tlm = Extra('fgs', fgnames) + """ + + def __init__(self, baseWebsky, name, fgnames): + self.name = name + self.fgnames = fgnames + self.directory = baseWebsky + + def __call__(self, idx): + return np.sum([hp.read_map(opj(self.directory, f'{fgname}.fits')) for fgname in self.fgnames], axis = 0) + + def get_name(self): + return self.name + \ No newline at end of file From 9b4bb59d75c99767ca9f23bad616a177bb372f5f Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Wed, 21 Jun 2023 10:50:07 +0200 Subject: [PATCH 4/7] generates fake simple points sources correlated with an input lensing map (convergence map) --- delensalot/sims/foregrounds/pointsources.py | 67 +++++++++++++++++++++ 1 file changed, 67 insertions(+) create mode 100644 delensalot/sims/foregrounds/pointsources.py diff --git a/delensalot/sims/foregrounds/pointsources.py b/delensalot/sims/foregrounds/pointsources.py new file mode 100644 index 00000000..7b647b75 --- /dev/null +++ b/delensalot/sims/foregrounds/pointsources.py @@ -0,0 +1,67 @@ +""" +Generates simple point sources, can be correlated with an input lensing convergence map. +""" + +import numpy as np +import healpy as hp + + +class Foreground(object): + def __init__(self): + pass + + @staticmethod + def randomizing_fg(mappa: np.ndarray): + """ + Randomizes the phase of the input map, preserving the amplitude. + """ + f = lambda z: np.abs(z) * np.exp(1j*np.random.uniform(0., 2.*np.pi, size = z.shape)) + return f(mappa) + + @staticmethod + def randomized_map(self, mappa: np.ndarray, nside: int): + """ + Randomizes the phase of the input map, preserving the amplitude. + """ + alm = hp.map2alm(mappa) + alm = self.randomizing_fg(alm) + return hp.alm2map(alm, nside) + + +class PointSourcesSimple(Foreground): + + def __init__(self, nside: int = 2048) -> None: + self.nside = nside + + @staticmethod + def phi_lm_to_kappa_lm(plm: np.ndarray) -> np.ndarray: + """ + Converts the input phi_lm to kappa_lm. + """ + lmax = hp.Alm.getlmax(plm.size) + ls = np.arange(0, lmax) + factor = (ls * (ls + 1.)) / 2. + return hp.almxfl(plm, factor) + + + def generate_ps(self, nsrc: int, amp: float = 200, seed: int = 0, plm: np.ndarray = None, factor: float = 2.) -> np.ndarray: + + rng = np.random.default_rng(seed) + + + if plm is not None: + klm = self.phi_lm_to_kappa_lm(plm) + kmap = hp.alm2map(klm, self.nside, verbose = False) + positions = np.where(rng.poisson(abs(kmap)*factor) > 0)[0] + nsources = len(positions) + mappa = np.zeros(hp.nside2npix(self.nside)) + + else: + nsources = rng.poisson(nsrc) + mappa = np.zeros(hp.nside2npix(self.nside)) + positions = np.random.randint(0, len(mappa), nsources) + + amplitudes = rng.poisson(amp, nsources) + mappa[positions] = amplitudes + + return mappa From 6fec14e4c2f8b89ff59f9c2ec9e5400cb931091a Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Wed, 21 Jun 2023 16:38:36 +0200 Subject: [PATCH 5/7] alternative way to generate positions in poisson map, simple match filter to create a mask for point sources --- delensalot/sims/foregrounds/pointsources.py | 43 +++++++++++++++++++-- 1 file changed, 39 insertions(+), 4 deletions(-) diff --git a/delensalot/sims/foregrounds/pointsources.py b/delensalot/sims/foregrounds/pointsources.py index 7b647b75..88d83360 100644 --- a/delensalot/sims/foregrounds/pointsources.py +++ b/delensalot/sims/foregrounds/pointsources.py @@ -18,7 +18,6 @@ def randomizing_fg(mappa: np.ndarray): f = lambda z: np.abs(z) * np.exp(1j*np.random.uniform(0., 2.*np.pi, size = z.shape)) return f(mappa) - @staticmethod def randomized_map(self, mappa: np.ndarray, nside: int): """ Randomizes the phase of the input map, preserving the amplitude. @@ -26,6 +25,25 @@ def randomized_map(self, mappa: np.ndarray, nside: int): alm = hp.map2alm(mappa) alm = self.randomizing_fg(alm) return hp.alm2map(alm, nside) + + + @staticmethod + def matched_filter(input_map_alm: np.ndarray, total_cl: np.ndarray, signal_cl: np.ndarray, nside: int): + """ + Returns the matched filter map. + """ + alm = hp.almxfl(input_map_alm, np.nan_to_num(1/total_cl)) + alm = hp.almxfl(alm, signal_cl) + return hp.alm2map(alm, nside) + + def mask_from_matched_filter(self, input_map_alm: np.ndarray, total_cl: np.ndarray, signal_cl: np.ndarray, nside: int, threshold: float = 0.5): + """ + Returns a mask from the matched filter map. + """ + mappa = self.matched_filter(input_map_alm, total_cl, signal_cl, nside) + SN_map = abs(mappa) / np.std(mappa) + mask = np.where(mappa > threshold, 1, 0) + return mask class PointSourcesSimple(Foreground): @@ -42,23 +60,40 @@ def phi_lm_to_kappa_lm(plm: np.ndarray) -> np.ndarray: ls = np.arange(0, lmax) factor = (ls * (ls + 1.)) / 2. return hp.almxfl(plm, factor) + + + def _get_position_from_kappa_default(self, rng, kappa: np.ndarray, factor: float = 0.5) -> np.ndarray: + """ + Returns the positions of the point sources from the input kappa map. + """ + positions = np.where(rng.poisson(1+kappa*factor) > 0)[0] + return positions + + + def _get_position_from_kappa_alternative(self, rng, kappa: np.ndarray, factor: float = 0.5) -> np.ndarray: + """ + Returns the positions of the point sources from the input kappa map. + """ + positions = np.where(rng.poisson(abs(kappa)*factor) > 0)[0] + return positions - def generate_ps(self, nsrc: int, amp: float = 200, seed: int = 0, plm: np.ndarray = None, factor: float = 2.) -> np.ndarray: + def generate_ps(self, nsrc: int, amp: float = 100, seed: int = 0, plm: np.ndarray = None, factor: float = 0.5) -> np.ndarray: rng = np.random.default_rng(seed) + mappa = np.zeros(hp.nside2npix(self.nside)) if plm is not None: + klm = self.phi_lm_to_kappa_lm(plm) kmap = hp.alm2map(klm, self.nside, verbose = False) positions = np.where(rng.poisson(abs(kmap)*factor) > 0)[0] nsources = len(positions) - mappa = np.zeros(hp.nside2npix(self.nside)) else: + nsources = rng.poisson(nsrc) - mappa = np.zeros(hp.nside2npix(self.nside)) positions = np.random.randint(0, len(mappa), nsources) amplitudes = rng.poisson(amp, nsources) From 41b128dae4d7696dd37c4de27a1a115f744f2d66 Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Wed, 21 Jun 2023 16:39:58 +0200 Subject: [PATCH 6/7] small update on the previous commit: alternative way to generate positions in poisson map, simple match filter to create a mask for point sources --- delensalot/sims/foregrounds/pointsources.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/delensalot/sims/foregrounds/pointsources.py b/delensalot/sims/foregrounds/pointsources.py index 88d83360..38fdb2eb 100644 --- a/delensalot/sims/foregrounds/pointsources.py +++ b/delensalot/sims/foregrounds/pointsources.py @@ -44,6 +44,13 @@ def mask_from_matched_filter(self, input_map_alm: np.ndarray, total_cl: np.ndarr SN_map = abs(mappa) / np.std(mappa) mask = np.where(mappa > threshold, 1, 0) return mask + + @staticmethod + def smooth_map(mappa: np.ndarray, fwhm: float, nside: int): + """ + Smooths the input map. + """ + return hp.smoothing(mappa, fwhm = np.radians(fwhm)) class PointSourcesSimple(Foreground): @@ -62,15 +69,16 @@ def phi_lm_to_kappa_lm(plm: np.ndarray) -> np.ndarray: return hp.almxfl(plm, factor) - def _get_position_from_kappa_default(self, rng, kappa: np.ndarray, factor: float = 0.5) -> np.ndarray: + @staticmethod + def _get_position_from_kappa_default(rng, kappa: np.ndarray, factor: float = 0.5) -> np.ndarray: """ Returns the positions of the point sources from the input kappa map. """ positions = np.where(rng.poisson(1+kappa*factor) > 0)[0] return positions - - def _get_position_from_kappa_alternative(self, rng, kappa: np.ndarray, factor: float = 0.5) -> np.ndarray: + @staticmethod + def _get_position_from_kappa_alternative(rng, kappa: np.ndarray, factor: float = 0.5) -> np.ndarray: """ Returns the positions of the point sources from the input kappa map. """ @@ -88,7 +96,7 @@ def generate_ps(self, nsrc: int, amp: float = 100, seed: int = 0, plm: np.ndarra klm = self.phi_lm_to_kappa_lm(plm) kmap = hp.alm2map(klm, self.nside, verbose = False) - positions = np.where(rng.poisson(abs(kmap)*factor) > 0)[0] + positions = self._get_position_from_kappa_default(rng, kmap, factor) nsources = len(positions) else: From 833407360ecf2305531cdc1c74788ed1d1e5ba2f Mon Sep 17 00:00:00 2001 From: Saladino93 Date: Sun, 9 Jul 2023 10:42:30 +0200 Subject: [PATCH 7/7] lognormal simulations utils --- delensalot/sims/lognormal/__init__.py | 0 delensalot/sims/lognormal/lognormal_utils.py | 170 +++++++++++++++++++ delensalot/sims/lognormal/sims_lognormal.py | 81 +++++++++ 3 files changed, 251 insertions(+) create mode 100644 delensalot/sims/lognormal/__init__.py create mode 100644 delensalot/sims/lognormal/lognormal_utils.py create mode 100644 delensalot/sims/lognormal/sims_lognormal.py diff --git a/delensalot/sims/lognormal/__init__.py b/delensalot/sims/lognormal/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/delensalot/sims/lognormal/lognormal_utils.py b/delensalot/sims/lognormal/lognormal_utils.py new file mode 100644 index 00000000..daaf7afe --- /dev/null +++ b/delensalot/sims/lognormal/lognormal_utils.py @@ -0,0 +1,170 @@ +"""Module with utilities for the lognormal sims module. +""" + +import numpy as np +import flt +import healpy as hp + + +def get_out_quantities_from_a_map(outmap: np.ndarray) -> tuple: + """ + Returns the mean, variance, skewness, lambda, muG and sigmaG from a given map. + """ + skewness = get_skew_from_map(outmap) + variance = get_variance_from_map(outmap) + mean = get_mean_from_map(outmap) + lamb = get_lambda_from_skew(skewness, variance, mean) + + alpha = get_alpha(mean, lamb) + sigmaG = get_sigma_gauss(alpha, variance) + muG = get_mu_gauss(alpha, variance) + return mean, variance, skewness, lamb, muG, sigmaG + + +def create_lognormal_single_map(inputcl: np.ndarray, nside: int, lmax_gen: int, mu: float = 0.0, lamb: float = 0.0): + """ + Creates a lognormal map with a given power spectrum and skewness. These are the parameters of the lognormal distribution that are specified with mu and lamb. + """ + + + alpha = get_alpha(mu, lamb) + xisinput = cl2xi(inputcl)/alpha/alpha + + xigaussian = np.log(xisinput+1) + clgaussian = xi2cl(xigaussian) + + vargauss = np.dot(np.arange(1, 2*len(clgaussian), 2), clgaussian)/(4*np.pi) + + lmax_gen = 2*nside-1 if lmax_gen is None else lmax_gen + almgaussian = hp.synalm(clgaussian, lmax = lmax_gen) #GENERATE TO HIGH LMAX + maps = hp.alm2map(almgaussian, nside = nside, pol = False) + + #vargauss = np.array([xigaussian[i, i][0] for i in range(Nfields)]) + #vargauss = np.array([np.var(m) for m in maps]) + + expmu = (mu+lamb)*np.exp(-vargauss*0.5) + maps = np.array(maps) + maps = np.exp(maps) + maps *= expmu + maps -= lamb + return maps + + +#shifted log-normal distribution +def shifted_lognormal_zero_mean(x, sigmaG, lamb): + #equation (22) of https://www.aanda.org/articles/aa/pdf/2011/12/aa17294-11.pdf + return np.exp(-(np.log(x/lamb+1)+sigmaG**2/2)**2./(2.*sigmaG**2.))/(x+lamb)/sigmaG/np.sqrt(2.*np.pi)*(x>-lamb) + + + +def cl2xi(cl: np.ndarray, closed = False): + """ + This goes from angular power spectrum to a correlation function calculated at theta points. + """ + ls = np.arange(0, len(cl)) + factorcl = (2*ls+1)/(4*np.pi) + coeffs = cl*factorcl + return flt.idlt(coeffs, closed = closed) + + +def theta(n, closed = False): + """ + Returns the theta for which the cl2xi are calculated, for a given n + """ + return flt.theta(n, closed = closed) + + +def xi2cl(xi: np.ndarray, closed = False): + """ + This goes from correlation function calculated at theta points to an angular power spectrum. + """ + ls = np.arange(0, len(xi)) + factorcl = (2*ls+1)/(4*np.pi) + return flt.dlt(xi, closed = closed)/factorcl + + + +def get_mean_from_map(mappa: np.ndarray): + """ + Mean from a map. + """ + return np.mean(mappa) + +def get_variance_from_map(mappa: np.ndarray): + """ + Variance from a map. + """ + return np.mean(mappa**2.)-np.mean(mappa)**2. + +def get_skew_from_map(mappa: np.ndarray): + """ + Skewness from a map. + """ + return np.mean((mappa-get_mean_from_map(mappa))**3.)/np.mean(mappa**2.)**1.5 + +def y_skew(skew): + """ + Formula (12) from https://arxiv.org/pdf/1602.08503.pdf + It relates the skewnees to some factor that is used to get the lambda parameter for the log-normal generation. + """ + result = 2+skew**2.+skew*np.sqrt(4+skew**2.) + result /= 2 + return np.power(result, 1/3.) + +def get_lambda_from_skew(skew, var, mu): + lmbda = np.sqrt(var)/skew*(1+y_skew(skew)+1/y_skew(skew))-mu + return lmbda + +def get_alpha(mu, lmbda): + """ + Below formula (7) from https://arxiv.org/pdf/1602.08503.pdf + """ + return mu+lmbda + +def get_mu_gauss(alpha, var): + """ + Gets the mu parameter for the Gaussian distribution for the log-normal + """ + result = np.log(alpha**2./np.sqrt(var+alpha**2.)) + return result + +def get_sigma_gauss(alpha, var): + """ + Gets the sigma parameter for the Gaussian distribution for the log-normal. + + Here the variance is the variance of the wanted log-normal field. + """ + result = np.log(1+var/alpha**2.) + result = np.sqrt(result) + return result + + +def suppress(l: np.ndarray, lsup: float = 7000, supindex: float = 10): + """ + Suppression factor at high ell. + """ + return np.exp(-1.0*np.power(l/lsup, supindex)) + +def suppress_cls(inputcl: np.ndarray, l: np.ndarray, lsup: float = 7000, supindex: float = 10): + """ + Suppresses an input angular power spectrum at high ell. + """ + return inputcl*suppress(l, lsup, supindex) + +def process_cl(inputcl: np.ndarray, function: callable = suppress_cls, **kwargs): + """ + Processes an input angular power spectrum with a given function. + + The function must take inputcl and l as arguments. + + Args: + inputcl (np.ndarray): Input angular power spectrum. + function (callable) default suppress_cls: Function to process the angular power spectrum. + **kwargs: Keyword arguments for the function. + + Returns: + np.ndarray: Processed angular power spectrum. + """ + ls = np.arange(0, len(inputcl)) + return function(inputcl = inputcl, l = ls, **kwargs) + \ No newline at end of file diff --git a/delensalot/sims/lognormal/sims_lognormal.py b/delensalot/sims/lognormal/sims_lognormal.py new file mode 100644 index 00000000..3385f1ec --- /dev/null +++ b/delensalot/sims/lognormal/sims_lognormal.py @@ -0,0 +1,81 @@ +""" +Generates a log-normal simulation with a given power spectrum and skewness. It uses methods described in e.g. https://arxiv.org/abs/1602.08503 + +NOTE: + An improvement could be done by setting some work at the unlcmb library level. Will leave this for now as an exploration. + The reason is that in the future users might want to set up their own generation of sims. There should be some simple way to ovverride + getting maps, maybe just some class inheritance somewhere. +""" + +import healpy as hp +import numpy as np +import lognormal_utils as lu +from delensalot.sims import sims_gaussian +import os + + +class sims_gaussian(sims_gaussian.sims_gaussian): + """Simulations with lognormal phi + + Args: + lib_dir: the phases of the CMB maps and the lensed CMBs will be stored there + lmax_cmb: cmb maps are generated down to this max multipole + cls_unl: dictionary of unlensed CMB spectra + dlmax, nside_lens, facres, nbands: lenspyx lensing module parameters + wcurl: include field rotation map in the lensing deflection (default to False for historical reasons) + + + This uses the cl_fid phi from sims_postborn to generate new lensing potential fields. + + """ + def __init__(self, lib_dir, lmax_cmb, cls_unl:dict, wcurl=False, + dlmax=1024, nside_lens=4096, epsilon=1e-5, cache_plm=True, mu: float = 0.0, var: float = 1.0, skew: float = 0.0, input_cl: np.ndarray = None, lmax_gen: int = 8000): + + lmax_plm = lmax_cmb + dlmax + mmax_plm = lmax_plm + + self.lmax_plm = lmax_plm + self.mmax_plm = mmax_plm + self.path = None + + self.cache_plm = cache_plm + self.wcurl = wcurl + self.epsilon = epsilon + + cmb_cls = {} + for k in cls_unl.keys(): + cmb_cls[k] = np.copy(cls_unl[k][:lmax_cmb + dlmax + 1]) + + super(sims_gaussian, self).__init__(lib_dir, lmax_cmb, cmb_cls, + dlmax=dlmax, nside_lens=nside_lens, epsilon=self.epsilon) + + if input_cl is None: + self.input_cl = cmb_cls['pp'] + + self.mu = mu + self.lamb = lu.get_lambda_from_skew(skew, var, mu) + self.lmax_gen = lmax_gen + + + @staticmethod + def kappa_lm_to_phi_lm(klm: np.ndarray) -> np.ndarray: + """ + Converts the input kappa_lm to phi_lm. + """ + lmax = hp.Alm.getlmax(klm.size) + ls = np.arange(0, lmax) + factor = np.nan_to_num(1/((ls * (ls + 1.)) / 2.)) + return hp.almxfl(klm, factor) + + def get_sim_plm(self, idx): + """ + Get a simulated lensing potential map + """ + fn = os.path.join(self.lib_dir, 'plm_in_%04d_lmax%s.fits'%(idx, self.lmax_plm)) + if not os.path.exists(fn): + klm = lu.create_lognormal_single_map(inputcl = self.input_cl, nside = self.nnside_lensside, lmax_gen = self.lmax_gen, mu = self.mu, lamb = self.lamb) + plm = self.kappa_lm_to_phi_lm(klm) + if self.cache_plm: + hp.write_alm(fn, plm) + return plm + return hp.read_alm(fn) \ No newline at end of file