From bbb8be6a42716e7365d972fa793aa961531a9397 Mon Sep 17 00:00:00 2001 From: Carter Date: Tue, 17 Oct 2023 15:20:57 -0500 Subject: [PATCH 01/18] Added a cder term (not in cross correlations) --- pyccl/nl_pt/ept.py | 2 ++ pyccl/nl_pt/tracers.py | 19 +++++++++++++++---- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index a50acc6a0..15fdba86b 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -260,6 +260,8 @@ def reshape_fastpt(tupl): reshape_fastpt(self.ia_tt) self.ia_mix = self.pt.IA_mix(**kw) reshape_fastpt(self.ia_mix) + self.ia_der = self.pt.IA_der(**kw) + reshape_fastpt(self.ia_der) # b1/bk power spectrum pks = {} diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index b2e5ae313..703bab2ad 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -10,7 +10,7 @@ from ..pyutils import _check_array_params -def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, +def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ader=None, Om_m2_for_c2=False, Om_m_fid=0.3): """ Function to convert from :math:`A_{ia}` values to :math:`c_{ia}` values, @@ -40,8 +40,10 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, Om_m = cosmo['Omega_m'] rho_crit = physical_constants.RHO_CRITICAL - c1 = c1delta = c2 = None + c1 = c1delta = c2 = cder = None gz = cosmo.growth_factor(1./(1+z)) + knorm = 1 #Units of Mpc/h, normalizes units out of derivative term + if a1 is not None: c1 = -1*a1*5e-14*rho_crit*Om_m/gz @@ -54,8 +56,10 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, c2 = a2*5*5e-14*rho_crit*Om_m**2/(Om_m_fid*gz**2) else: # DES convention c2 = a2*5*5e-14*rho_crit*Om_m/(gz**2) + if ader is not None: + cder= ader*knorm**2*5e-14*rho_crit*Om_m/gz - return c1, c1delta, c2 + return c1, c1delta, c2, cder class PTTracer(CCLAutoRepr): @@ -212,7 +216,7 @@ class PTIntrinsicAlignmentTracer(PTTracer): """ type = 'IA' - def __init__(self, c1, c2=None, cdelta=None): + def __init__(self, c1, c2=None, cdelta=None, cder=None): self.biases = {} @@ -222,6 +226,8 @@ def __init__(self, c1, c2=None, cdelta=None): self.biases['c2'] = self._get_bias_function(c2) # Initialize cdelta self.biases['cdelta'] = self._get_bias_function(cdelta) + # Initialize cder + self.biases['cder'] = self._get_bias_function(cder) @property def c1(self): @@ -240,3 +246,8 @@ def cdelta(self): """Internal overdensity bias function. """ return self.biases['cdelta'] + + @property + def cder(self): + """Internal derivative bias function + """ From 67b0985ff3cf0d5e4c8b1a9e70838fe6add23ee6 Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 19 Oct 2023 15:32:01 -0500 Subject: [PATCH 02/18] Added derivative contribution to tracers --- pyccl/nl_pt/ept.py | 40 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 15fdba86b..f7ca80c8b 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -141,6 +141,7 @@ def __init__(self, *, with_NC=False, with_IA=False, log10k_min=-4, log10k_max=2, nk_per_decade=20, a_arr=None, k_cutoff=None, n_exp_cutoff=4, b1_pk_kind='nonlinear', bk2_pk_kind='nonlinear', + ak2_pk_kind='nonlinear', pad_factor=1.0, low_extrap=-5.0, high_extrap=3.0, P_window=None, C_window=0.75, sub_lowk=False): self.with_matter_1loop = with_matter_1loop @@ -191,9 +192,12 @@ def __init__(self, *, with_NC=False, with_IA=False, raise ValueError(f"Unknown P(k) prescription {b1_pk_kind}") if bk2_pk_kind not in ['linear', 'nonlinear', 'pt']: raise ValueError(f"Unknown P(k) prescription {bk2_pk_kind}") + if ak2_pk_kind not in ['linear', 'nonlinear', 'pt']: + raise ValueError(f"Unknown P(k) prescription in {ak2_pk_kind}") self.b1_pk_kind = b1_pk_kind self.bk2_pk_kind = bk2_pk_kind - if (self.b1_pk_kind == 'pt') or (self.bk2_pk_kind == 'pt'): + self.ak2_pk_kind = bk2_pk_kind + if (self.b1_pk_kind == 'pt') or (self.bk2_pk_kind == 'pt') or (self.ak2_pk_kind == 'pt'): self.with_matter_1loop = True # Initialize all expensive arrays to ``None``. @@ -283,6 +287,24 @@ def reshape_fastpt(tupl): self.pk_b1 = pks[self.b1_pk_kind] self.pk_bk = pks[self.bk2_pk_kind] + # a1/ak power spectrum + pksa = {} + if 'nonlinear' in [self.ak2_pk_kind]: + pksa['nonlinear'] = np.array([cosmo.nonlin_matter_power(self.k_s, a) + for a in self.a_s]) + if 'linear' in [self.ak2_pk_kind]: + pksa['linear'] = np.array([cosmo.linear_matter_power(self.k_s, a) + for a in self.a_s]) + if 'pt' in [self.ak2_pk_kind]: + if 'linear' in pksa: + pka = pks['linear'] + else: + pka = np.array([cosmo.linear_matter_power(self.k_s, a) for a in self.a_s]) + pk += self._g4T*self.one_loop_dd[0] + pks['pt'] = pk + self.pk_ak=pks[self.ak2_pk_kind] + + # Reset template power spectra self._pk2d_temp = {} self._cosmo = cosmo @@ -366,6 +388,7 @@ def _get_pgi(self, trg, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases b1 = trg.b1(self.z_s) @@ -381,10 +404,12 @@ def _get_pgi(self, trg, tri): c1 = tri.c1(self.z_s) c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) + ck = tri.cder(self.z_s) pgi = b1[:, None] * (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + - (self._g4*c2)[:, None] * (a0e2 + b0e2)) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + + ck[:, None]*(Pak2)) return pgi*self.exp_cutoff def _get_pgm(self, trg): @@ -444,14 +469,17 @@ def _get_pii(self, tr1, tr2, return_bb=False): a00e, c00e, a0e0e, a0b0b = self.ia_ta ae2e2, ab2b2 = self.ia_tt a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases c11 = tr1.c1(self.z_s) c21 = tr1.c2(self.z_s) cd1 = tr1.cdelta(self.z_s) + ck1 = tr1.cder(self.z_s) c12 = tr2.c1(self.z_s) c22 = tr2.c2(self.z_s) cd2 = tr2.cdelta(self.z_s) + ck2 = tr2.cder(self.z_s) if return_bb: pii = ((cd1*cd2*self._g4)[:, None]*a0b0b + @@ -463,7 +491,8 @@ def _get_pii(self, tr1, tr2, return_bb=False): (cd1*cd2*self._g4)[:, None]*a0e0e + (c21*c22*self._g4)[:, None]*ae2e2 + ((c11*c22+c21*c12)*self._g4)[:, None]*(a0e2+b0e2) + - ((cd1*c22+cd2*c21)*self._g4)[:, None]*d0ee2) + ((cd1*c22+cd2*c21)*self._g4)[:, None]*d0ee2 + + (ck1*c12 + ck2*c11)[:,None]* (Pak2)) return pii*self.exp_cutoff @@ -485,15 +514,18 @@ def _get_pim(self, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases c1 = tri.c1(self.z_s) c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) + ck = tri.cder(self.z_s) pim = (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + - (self._g4*c2)[:, None] * (a0e2 + b0e2)) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + + (ck)[:,None]*Pak2) return pim*self.exp_cutoff def _get_pmm(self): From 0a1222ea4485aec4ca2a0d14575bb5cc7b2a783d Mon Sep 17 00:00:00 2001 From: Carter Date: Mon, 23 Oct 2023 11:20:23 -0500 Subject: [PATCH 03/18] Changed cder to ck --- pyccl/nl_pt/ept.py | 20 ++++++++++---------- pyccl/nl_pt/tracers.py | 15 ++++++++------- 2 files changed, 18 insertions(+), 17 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index f7ca80c8b..3ffc404f4 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -196,8 +196,8 @@ def __init__(self, *, with_NC=False, with_IA=False, raise ValueError(f"Unknown P(k) prescription in {ak2_pk_kind}") self.b1_pk_kind = b1_pk_kind self.bk2_pk_kind = bk2_pk_kind - self.ak2_pk_kind = bk2_pk_kind - if (self.b1_pk_kind == 'pt') or (self.bk2_pk_kind == 'pt') or (self.ak2_pk_kind == 'pt'): + self.ak2_pk_kind = ak2_pk_kind + if (self.b1_pk_kind == 'pt') or (self.bk2_pk_kind == 'pt'): self.with_matter_1loop = True # Initialize all expensive arrays to ``None``. @@ -297,12 +297,12 @@ def reshape_fastpt(tupl): for a in self.a_s]) if 'pt' in [self.ak2_pk_kind]: if 'linear' in pksa: - pka = pks['linear'] + pka = pksa['linear'] else: pka = np.array([cosmo.linear_matter_power(self.k_s, a) for a in self.a_s]) - pk += self._g4T*self.one_loop_dd[0] - pks['pt'] = pk - self.pk_ak=pks[self.ak2_pk_kind] + pka += self._g4T*self.one_loop_dd[0] + pksa['pt'] = pka + self.pk_ak=pksa[self.ak2_pk_kind] # Reset template power spectra @@ -404,7 +404,7 @@ def _get_pgi(self, trg, tri): c1 = tri.c1(self.z_s) c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) - ck = tri.cder(self.z_s) + ck = tri.ck(self.z_s) pgi = b1[:, None] * (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + @@ -475,11 +475,11 @@ def _get_pii(self, tr1, tr2, return_bb=False): c11 = tr1.c1(self.z_s) c21 = tr1.c2(self.z_s) cd1 = tr1.cdelta(self.z_s) - ck1 = tr1.cder(self.z_s) + ck1 = tr1.ck(self.z_s) c12 = tr2.c1(self.z_s) c22 = tr2.c2(self.z_s) cd2 = tr2.cdelta(self.z_s) - ck2 = tr2.cder(self.z_s) + ck2 = tr2.ck(self.z_s) if return_bb: pii = ((cd1*cd2*self._g4)[:, None]*a0b0b + @@ -520,7 +520,7 @@ def _get_pim(self, tri): c1 = tri.c1(self.z_s) c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) - ck = tri.cder(self.z_s) + ck = tri.ck(self.z_s) pim = (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index 703bab2ad..1973c789c 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -10,7 +10,7 @@ from ..pyutils import _check_array_params -def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ader=None, +def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, Om_m2_for_c2=False, Om_m_fid=0.3): """ Function to convert from :math:`A_{ia}` values to :math:`c_{ia}` values, @@ -56,10 +56,10 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ader=None, c2 = a2*5*5e-14*rho_crit*Om_m**2/(Om_m_fid*gz**2) else: # DES convention c2 = a2*5*5e-14*rho_crit*Om_m/(gz**2) - if ader is not None: - cder= ader*knorm**2*5e-14*rho_crit*Om_m/gz + if ak is not None: + ck= ak*knorm**2*5e-14*rho_crit*Om_m/gz - return c1, c1delta, c2, cder + return c1, c1delta, c2, ck class PTTracer(CCLAutoRepr): @@ -216,7 +216,7 @@ class PTIntrinsicAlignmentTracer(PTTracer): """ type = 'IA' - def __init__(self, c1, c2=None, cdelta=None, cder=None): + def __init__(self, c1, c2=None, cdelta=None, ck=None): self.biases = {} @@ -227,7 +227,7 @@ def __init__(self, c1, c2=None, cdelta=None, cder=None): # Initialize cdelta self.biases['cdelta'] = self._get_bias_function(cdelta) # Initialize cder - self.biases['cder'] = self._get_bias_function(cder) + self.biases['ck'] = self._get_bias_function(ck) @property def c1(self): @@ -248,6 +248,7 @@ def cdelta(self): return self.biases['cdelta'] @property - def cder(self): + def ck(self): """Internal derivative bias function """ + return self.biases['ck'] From b3f3d918e27f3b90e4716afdf93fd1d8b5b699ca Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 9 Nov 2023 12:56:24 -0600 Subject: [PATCH 04/18] Added unit test functionality for der term --- pyccl/nl_pt/ept.py | 68 ++++++++++++++++++++++++----------- pyccl/nl_pt/tracers.py | 2 +- pyccl/tests/test_ept_power.py | 21 ++++++----- 3 files changed, 61 insertions(+), 30 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 3ffc404f4..25949b7ed 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -13,19 +13,19 @@ 'm:m': 'm:m', 'm:b1': 'm:m', 'm:b2': 'm:b2', 'm:b3nl': 'm:b3nl', 'm:bs': 'm:bs', 'm:bk2': 'm:bk2', 'm:c1': 'm:m', 'm:c2': 'm:c2', 'm:cdelta': 'm:cdelta', - 'b1:b1': 'm:m', 'b1:b2': 'm:b2', 'b1:b3nl': 'm:b3nl', - 'b1:bs': 'm:bs', 'b1:bk2': 'm:bk2', 'b1:c1': 'm:m', - 'b1:c2': 'm:c2', 'b1:cdelta': 'm:cdelta', 'b2:b2': 'b2:b2', - 'b2:b3nl': 'zero', 'b2:bs': 'b2:bs', 'b2:bk2': 'zero', - 'b2:c1': 'zero', 'b2:c2': 'zero', 'b2:cdelta': 'zero', - 'b3nl:b3nl': 'zero', 'b3nl:bs': 'zero', + 'm:ck' : 'm:ck', 'b1:b1': 'm:m', 'b1:b2': 'm:b2', + 'b1:b3nl': 'm:b3nl', 'b1:bs': 'm:bs', 'b1:bk2': 'm:bk2', + 'b1:c1': 'm:m', 'b1:c2': 'm:c2', 'b1:cdelta': 'm:cdelta', 'b1:ck' : 'm:ck', + 'b2:b2': 'b2:b2', 'b2:b3nl': 'zero', 'b2:bs': 'b2:bs', + 'b2:bk2': 'zero','b2:c1': 'zero', 'b2:c2': 'zero', + 'b2:cdelta': 'zero', 'b3nl:b3nl': 'zero', 'b3nl:bs': 'zero', 'b3nl:bk2': 'zero', 'b3nl:c1': 'zero', 'b3nl:c2': 'zero', 'b3nl:cdelta': 'zero', 'bs:bs': 'bs:bs', 'bs:bk2': 'zero', 'bs:c1': 'zero', 'bs:c2': 'zero', 'bs:cdelta': 'zero', 'bk2:bk2': 'zero', 'bk2:c1': 'zero', 'bk2:c2': 'zero', 'bk2:cdelta': 'zero', 'c1:c1': 'm:m', - 'c1:c2': 'm:c2', 'c1:cdelta': 'm:cdelta', 'c2:c2': 'c2:c2', - 'c2:cdelta': 'c2:cdelta', 'cdelta:cdelta': 'cdelta:cdelta'} + 'c1:c2': 'm:c2', 'c1:cdelta': 'm:cdelta', 'c1:ck' : 'm:ck', 'c2:c2': 'c2:c2', + 'c2:cdelta': 'c2:cdelta', 'cdelta:cdelta': 'cdelta:cdelta', 'ck:ck' : 'zero'} class EulerianPTCalculator(CCLAutoRepr): @@ -47,7 +47,7 @@ class EulerianPTCalculator(CCLAutoRepr): .. math:: s^I_{ij}=c_1\\,s_{ij}+c_2(s_{ik}s_{jk}-s^2\\delta_{ik}/3) - +c_\\delta\\,\\delta\\,s_{ij} + +c_\\delta\\,\\delta\\,s_{ij + c_k\\k^2\\,s_{ij}} (note that the higher-order terms are not divided by 2!). @@ -115,6 +115,9 @@ class EulerianPTCalculator(CCLAutoRepr): bk2_pk_kind (:obj:`str`): power spectrum to use for the non-local bias terms in the expansion. Same options and default as ``b1_pk_kind``. + ak2_pk_kind (:obj:`str`): power spectrum to use for the derivative + term of the IA expansion. Same options and default as + ``b1_pk_kind``. pad_factor (:obj:`float`): fraction of the :math:`\\log_{10}(k)` interval you to add as padding for FFTLog calculations. low_extrap (:obj:`float`): decimal logaritm of the minimum Fourier @@ -131,6 +134,8 @@ class EulerianPTCalculator(CCLAutoRepr): documentation for more details. sub_lowk (:obj:`bool`): if ``True``, the small-scale white noise contribution to some of the terms will be subtracted. + usefptk (::obj:`bool`):) if `True``, will use the FASTPT IA k2 + term instead of the CCL IA k2 term """ __repr_attrs__ = __eq_attrs__ = ('with_NC', 'with_IA', 'with_matter_1loop', 'k_s', 'a_s', 'exp_cutoff', 'b1_pk_kind', @@ -143,11 +148,11 @@ def __init__(self, *, with_NC=False, with_IA=False, b1_pk_kind='nonlinear', bk2_pk_kind='nonlinear', ak2_pk_kind='nonlinear', pad_factor=1.0, low_extrap=-5.0, high_extrap=3.0, - P_window=None, C_window=0.75, sub_lowk=False): + P_window=None, C_window=0.75, sub_lowk=False, usefptk=False): self.with_matter_1loop = with_matter_1loop self.with_NC = with_NC self.with_IA = with_IA - + self.ufpt=usefptk # Set FAST-PT parameters self.fastpt_par = {'pad_factor': pad_factor, 'low_extrap': low_extrap, @@ -180,7 +185,16 @@ def __init__(self, *, with_NC=False, with_IA=False, self.exp_cutoff = 1 # Call FAST-PT - import fastpt as fpt + try: + import fastpt as fpt + except: + raise ImportError("Your attempted import of FAST-PT has failed. You either dont have fast-pt installed, or have the wrong version. Try running pip install fast-pt or conda install fast-pt, then try again") + + if( not hasattr(fpt, "IA_ta")): + raise ValueError(f"Your FAST-PT version lacks a required attribute. You may have the wrong fast-pt install. Try running pip install fast-pt or conda install fast-pt, then try again") + + if( not hasattr(fpt, "IA_der") and self.ufpt): + raise ValueError(f"You need a newer version of FAST-PT to use fpt for k2 term") n_pad = int(self.fastpt_par['pad_factor'] * len(self.k_s)) self.pt = fpt.FASTPT(self.k_s, to_do=to_do, low_extrap=self.fastpt_par['low_extrap'], @@ -264,8 +278,9 @@ def reshape_fastpt(tupl): reshape_fastpt(self.ia_tt) self.ia_mix = self.pt.IA_mix(**kw) reshape_fastpt(self.ia_mix) - self.ia_der = self.pt.IA_der(**kw) - reshape_fastpt(self.ia_der) + if(self.ufpt): + self.ia_der = self.pt.IA_der(**kw) + reshape_fastpt(self.ia_der) # b1/bk power spectrum pks = {} @@ -287,7 +302,7 @@ def reshape_fastpt(tupl): self.pk_b1 = pks[self.b1_pk_kind] self.pk_bk = pks[self.bk2_pk_kind] - # a1/ak power spectrum + # ak power spectrum pksa = {} if 'nonlinear' in [self.ak2_pk_kind]: pksa['nonlinear'] = np.array([cosmo.nonlin_matter_power(self.k_s, a) @@ -388,7 +403,10 @@ def _get_pgi(self, trg, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix - Pak2 = self.pk_ak*(self.k_s**2)[None, :] + if(self.ufpt): + Pak2 = self.ia_der + else: + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases b1 = trg.b1(self.z_s) @@ -409,7 +427,7 @@ def _get_pgi(self, trg, tri): pgi = b1[:, None] * (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:, None]*(Pak2)) + ck[:, None] * Pak2) return pgi*self.exp_cutoff def _get_pgm(self, trg): @@ -469,7 +487,10 @@ def _get_pii(self, tr1, tr2, return_bb=False): a00e, c00e, a0e0e, a0b0b = self.ia_ta ae2e2, ab2b2 = self.ia_tt a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix - Pak2 = self.pk_ak*(self.k_s**2)[None, :] + if(self.ufpt): + Pak2 = self.ia_der + else: + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases c11 = tr1.c1(self.z_s) @@ -492,7 +513,7 @@ def _get_pii(self, tr1, tr2, return_bb=False): (c21*c22*self._g4)[:, None]*ae2e2 + ((c11*c22+c21*c12)*self._g4)[:, None]*(a0e2+b0e2) + ((cd1*c22+cd2*c21)*self._g4)[:, None]*d0ee2 + - (ck1*c12 + ck2*c11)[:,None]* (Pak2)) + (ck1*c12 + ck2*c11)[:,None] * (Pak2)) return pii*self.exp_cutoff @@ -514,7 +535,10 @@ def _get_pim(self, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix - Pak2 = self.pk_ak*(self.k_s**2)[None, :] + if(self.ufpt): + Pak2 = self.ia_der + else: + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases c1 = tri.c1(self.z_s) @@ -525,7 +549,7 @@ def _get_pim(self, tri): pim = (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + - (ck)[:,None]*Pak2) + ck[:,None] * Pak2) return pim*self.exp_cutoff def _get_pmm(self): @@ -690,6 +714,8 @@ def get_pk2d_template(self, kind, *, extrap_order_lok=1, pk = self._g4T * (self.ia_mix[0]+self.ia_mix[1]) elif pk_name == 'm:cdelta': pk = self._g4T * (self.ia_ta[0]+self.ia_ta[1]) + elif pk_name == 'm:ck': + pk = self.pk_ak*(self.k_s**2) elif pk_name == 'b2:b2': if self.fastpt_par['sub_lowk']: s4 = self.dd_bias[7] diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index 1973c789c..48608b4a0 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -42,7 +42,7 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, rho_crit = physical_constants.RHO_CRITICAL c1 = c1delta = c2 = cder = None gz = cosmo.growth_factor(1./(1+z)) - knorm = 1 #Units of Mpc/h, normalizes units out of derivative term + knorm = 1 #Units of Mpc/h, normalizes units out of ck term if a1 is not None: diff --git a/pyccl/tests/test_ept_power.py b/pyccl/tests/test_ept_power.py index 9f4a95224..b9cf7703d 100644 --- a/pyccl/tests/test_ept_power.py +++ b/pyccl/tests/test_ept_power.py @@ -3,13 +3,14 @@ import pytest from contextlib import nullcontext + COSMO = ccl.Cosmology( Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96, transfer_function='bbks', matter_power_spectrum='linear') TRS = {'TG': ccl.nl_pt.PTNumberCountsTracer(b1=2.0, b2=2.0, bs=2.0), 'TI': ccl.nl_pt.PTIntrinsicAlignmentTracer(c1=2.0, c2=2.0, - cdelta=2.0), + cdelta=2.0, ck=2.0), 'TM': ccl.nl_pt.PTMatterTracer()} PTC = ccl.nl_pt.EulerianPTCalculator(with_NC=True, with_IA=True, with_matter_1loop=True, @@ -58,7 +59,7 @@ def test_ept_pk2d_bb_smoke(): def test_ept_get_pk2d_nl(nl): ptc = ccl.nl_pt.EulerianPTCalculator( with_NC=True, with_IA=True, with_matter_1loop=True, - b1_pk_kind=nl, bk2_pk_kind=nl, cosmo=COSMO) + b1_pk_kind=nl, bk2_pk_kind=nl, ak2_pk_kind=nl,cosmo=COSMO) pk = ptc.get_biased_pk2d(TRS['TG']) assert isinstance(pk, ccl.Pk2D) @@ -72,7 +73,7 @@ def test_ept_k2pk_types(typ_nlin, typ_nloc): tm = ccl.nl_pt.PTNumberCountsTracer(1., 0., 0.) ptc1 = ccl.nl_pt.EulerianPTCalculator( with_NC=True, with_IA=True, with_matter_1loop=True, - b1_pk_kind=typ_nlin, bk2_pk_kind=typ_nloc, + b1_pk_kind=typ_nlin, bk2_pk_kind=typ_nloc, ak2_pk_kind=typ_nloc, cosmo=COSMO) ptc2 = ccl.nl_pt.EulerianPTCalculator( with_NC=True, with_IA=True, with_matter_1loop=True, @@ -92,7 +93,7 @@ def test_ept_deconstruction(kind): with_matter_1loop=True, cosmo=COSMO, sub_lowk=True) b_nc = ['b1', 'b2', 'b3nl', 'bs', 'bk2'] - b_ia = ['c1', 'c2', 'cdelta'] + b_ia = ['c1', 'c2', 'cdelta', 'ck'] pk1 = ptc.get_pk2d_template(kind) def get_tr(tn): @@ -110,14 +111,14 @@ def get_tr(tn): bdict[tn] = 1.0 return ccl.nl_pt.PTIntrinsicAlignmentTracer( c1=bdict['c1'], c2=bdict['c2'], - cdelta=bdict['cdelta']) + cdelta=bdict['cdelta'], ck=bdict['ck']) tn1, tn2 = kind.split(':') t1 = get_tr(tn1) t2 = get_tr(tn2) is_nl = tn1 in ["b2", "bs", "bk2", "b3nl"] - is_g = tn2 in ["c1", "c2", "cdelta"] + is_g = tn2 in ["c1", "c2", "cdelta", 'ck'] with pytest.warns(ccl.CCLWarning) if is_nl and is_g else nullcontext(): pk2 = ptc.get_biased_pk2d(t1, tracer2=t2) if pk1 is None: @@ -135,7 +136,7 @@ def get_tr(tn): @pytest.mark.parametrize('kind', ['c2:c2', 'c2:cdelta', 'cdelta:cdelta']) def test_ept_deconstruction_bb(kind): - b_ia = ['c1', 'c2', 'cdelta'] + b_ia = ['c1', 'c2', 'cdelta', 'ck'] pk1 = PTC.get_pk2d_template(kind, return_ia_bb=True) def get_tr(tn): @@ -143,7 +144,7 @@ def get_tr(tn): bdict[tn] = 1.0 return ccl.nl_pt.PTIntrinsicAlignmentTracer( c1=bdict['c1'], c2=bdict['c2'], - cdelta=bdict['cdelta']) + cdelta=bdict['cdelta'], ck=bdict['ck']) tn1, tn2 = kind.split(':') t1 = get_tr(tn1) @@ -213,6 +214,10 @@ def test_ept_calculator_raises(): # Wrong type of b1 and bk2 power spectra with pytest.raises(ValueError): ccl.nl_pt.EulerianPTCalculator(bk2_pk_kind='non-linear') + + #wrong type of ak2 power spectra + with pytest.raises(ValueError): + ccl.nl_pt.EulerianPTCalculator(ak2_pk_kind="non-linear") # Uninitialized templates with pytest.raises(ccl.CCLError): From 159f3330d3efbe3d65406338656f74e086f2ab2a Mon Sep 17 00:00:00 2001 From: Carter Date: Tue, 14 Nov 2023 10:17:13 -0600 Subject: [PATCH 05/18] changed cder to ck in tracers.py --- pyccl/nl_pt/tracers.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index 48608b4a0..73c6daa8b 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -40,10 +40,9 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, Om_m = cosmo['Omega_m'] rho_crit = physical_constants.RHO_CRITICAL - c1 = c1delta = c2 = cder = None + c1 = c1delta = c2 = ck = None gz = cosmo.growth_factor(1./(1+z)) - knorm = 1 #Units of Mpc/h, normalizes units out of ck term - + knorm = 1 # Units of Mpc/h, normalizes units out of ck term if a1 is not None: c1 = -1*a1*5e-14*rho_crit*Om_m/gz @@ -57,7 +56,7 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, else: # DES convention c2 = a2*5*5e-14*rho_crit*Om_m/(gz**2) if ak is not None: - ck= ak*knorm**2*5e-14*rho_crit*Om_m/gz + ck = ak*knorm**2*5e-14*rho_crit*Om_m/gz return c1, c1delta, c2, ck @@ -246,7 +245,7 @@ def cdelta(self): """Internal overdensity bias function. """ return self.biases['cdelta'] - + @property def ck(self): """Internal derivative bias function From 385d68638707a7d4c57f0b851fe013d1bdb0c21b Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 7 Dec 2023 12:25:59 -0600 Subject: [PATCH 06/18] Adding tij to CCL IA terms --- pyccl/nl_pt/ept.py | 45 ++++++++++++++++++++++++++++------- pyccl/nl_pt/tracers.py | 20 ++++++++++++---- pyccl/tests/test_ept_power.py | 6 ++--- 3 files changed, 55 insertions(+), 16 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 25949b7ed..44fd065f9 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -13,19 +13,22 @@ 'm:m': 'm:m', 'm:b1': 'm:m', 'm:b2': 'm:b2', 'm:b3nl': 'm:b3nl', 'm:bs': 'm:bs', 'm:bk2': 'm:bk2', 'm:c1': 'm:m', 'm:c2': 'm:c2', 'm:cdelta': 'm:cdelta', - 'm:ck' : 'm:ck', 'b1:b1': 'm:m', 'b1:b2': 'm:b2', + 'm:ck': 'm:ck', 'b1:b1': 'm:m', 'b1:b2': 'm:b2', 'b1:b3nl': 'm:b3nl', 'b1:bs': 'm:bs', 'b1:bk2': 'm:bk2', - 'b1:c1': 'm:m', 'b1:c2': 'm:c2', 'b1:cdelta': 'm:cdelta', 'b1:ck' : 'm:ck', + 'b1:c1': 'm:m', 'b1:c2': 'm:c2', 'b1:cdelta': 'm:cdelta', 'b1:ck': 'm:ck', 'b2:b2': 'b2:b2', 'b2:b3nl': 'zero', 'b2:bs': 'b2:bs', - 'b2:bk2': 'zero','b2:c1': 'zero', 'b2:c2': 'zero', + 'b2:bk2': 'zero', 'b2:c1': 'zero', 'b2:c2': 'zero', 'b2:cdelta': 'zero', 'b3nl:b3nl': 'zero', 'b3nl:bs': 'zero', 'b3nl:bk2': 'zero', 'b3nl:c1': 'zero', 'b3nl:c2': 'zero', 'b3nl:cdelta': 'zero', 'bs:bs': 'bs:bs', 'bs:bk2': 'zero', 'bs:c1': 'zero', 'bs:c2': 'zero', 'bs:cdelta': 'zero', 'bk2:bk2': 'zero', 'bk2:c1': 'zero', 'bk2:c2': 'zero', 'bk2:cdelta': 'zero', 'c1:c1': 'm:m', - 'c1:c2': 'm:c2', 'c1:cdelta': 'm:cdelta', 'c1:ck' : 'm:ck', 'c2:c2': 'c2:c2', - 'c2:cdelta': 'c2:cdelta', 'cdelta:cdelta': 'cdelta:cdelta', 'ck:ck' : 'zero'} + 'c1:c2': 'm:c2', 'c1:cdelta': 'm:cdelta', 'c1:ck': 'm:ck', + 'c2:c2': 'c2:c2', 'c2:cdelta': 'c2:cdelta', + 'cdelta:cdelta': 'cdelta:cdelta', 'ck:ck': 'zero', 'm:ct': 'm:ct', + 'b1:ct': 'm:ct', 'c1:ct': 'm:ct','c2:ct': 'c2:ct', + 'cdelta:ct': 'cdelta:ct', 'ct:ct': 'ct:ct', 'ck:ct': 'm:ct'} class EulerianPTCalculator(CCLAutoRepr): @@ -195,6 +198,8 @@ def __init__(self, *, with_NC=False, with_IA=False, if( not hasattr(fpt, "IA_der") and self.ufpt): raise ValueError(f"You need a newer version of FAST-PT to use fpt for k2 term") + if (not hasattr(fpt, "IA_tij")): + raise ValueError(f"You are using an older version of FAST-PT, please run pip install fast-pt or conda install fast-pt, then try again") n_pad = int(self.fastpt_par['pad_factor'] * len(self.k_s)) self.pt = fpt.FASTPT(self.k_s, to_do=to_do, low_extrap=self.fastpt_par['low_extrap'], @@ -281,6 +286,7 @@ def reshape_fastpt(tupl): if(self.ufpt): self.ia_der = self.pt.IA_der(**kw) reshape_fastpt(self.ia_der) + self.ia_tij = self.pt.IA_tij(**kw) # b1/bk power spectrum pks = {} @@ -403,6 +409,8 @@ def _get_pgi(self, trg, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + tijdsij, tij2sij, tijtij, tijsij = self.ia_ta + if(self.ufpt): Pak2 = self.ia_der else: @@ -423,11 +431,13 @@ def _get_pgi(self, trg, tri): c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) ck = tri.ck(self.z_s) + ct = tri.ct(self.z_s) pgi = b1[:, None] * (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:, None] * Pak2) + ck[:, None] * Pak2 + + ct[:, None] * tijsij) return pgi*self.exp_cutoff def _get_pgm(self, trg): @@ -487,6 +497,7 @@ def _get_pii(self, tr1, tr2, return_bb=False): a00e, c00e, a0e0e, a0b0b = self.ia_ta ae2e2, ab2b2 = self.ia_tt a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + tijdsij, tij2sij, tijtij, tijsij = self.ia_ta if(self.ufpt): Pak2 = self.ia_der else: @@ -497,10 +508,12 @@ def _get_pii(self, tr1, tr2, return_bb=False): c21 = tr1.c2(self.z_s) cd1 = tr1.cdelta(self.z_s) ck1 = tr1.ck(self.z_s) + ct1 = tr1.ct(self.z_s) c12 = tr2.c1(self.z_s) c22 = tr2.c2(self.z_s) cd2 = tr2.cdelta(self.z_s) ck2 = tr2.ck(self.z_s) + ct2 = tr2.ct(self.z_s) if return_bb: pii = ((cd1*cd2*self._g4)[:, None]*a0b0b + @@ -513,7 +526,12 @@ def _get_pii(self, tr1, tr2, return_bb=False): (c21*c22*self._g4)[:, None]*ae2e2 + ((c11*c22+c21*c12)*self._g4)[:, None]*(a0e2+b0e2) + ((cd1*c22+cd2*c21)*self._g4)[:, None]*d0ee2 + - (ck1*c12 + ck2*c11)[:,None] * (Pak2)) + (ck1*c12 + ck2*c11)[:,None] * (Pak2) + + (ct1*c12 + ct2*c11)[:,None]*(tijsij) + + (ct1*c22 + ct2*c21)[:,None] * (tij2sij) + + (ct1*cd2 + ct2*cd1)[:,None] * (tijdsij) + + (ct1*ct2)[:,None] * (tijtij) + + (ct1*ck2 + ct2*ck1)[:,None] * (tijsij)) return pii*self.exp_cutoff @@ -535,6 +553,7 @@ def _get_pim(self, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + tijdsij, tij2sij, tijtij, tijsij = self.ia_ta if(self.ufpt): Pak2 = self.ia_der else: @@ -545,11 +564,13 @@ def _get_pim(self, tri): c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) ck = tri.ck(self.z_s) + ct = tri.ct(self.z_s) pim = (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:,None] * Pak2) + ck[:,None] * Pak2 + + ct[:,None] * tijsij) return pim*self.exp_cutoff def _get_pmm(self): @@ -740,6 +761,14 @@ def get_pk2d_template(self, kind, *, extrap_order_lok=1, pk = self._g4T * self.ia_ta[2] elif pk_name == 'cdelta:cdelta_bb': pk = self._g4T * self.ia_ta[3] + elif pk_name == 'm:ct': + pk = self.ia_tij[3] + elif pk_name == 'c2:ct': + pk = self._g4T * self.ia_tij[1] + elif pk_name == 'cdelta:ct': + pk = self._g4T * self.ia_tij[0] + elif pk_name == 'ct:ct': + pk = self.ia_tij[2] elif pk_name == 'zero': # If zero, store None and return self._pk2d_temp[pk_name] = None diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index 73c6daa8b..394ae6851 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -10,7 +10,7 @@ from ..pyutils import _check_array_params -def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, +def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, at=None, Om_m2_for_c2=False, Om_m_fid=0.3): """ Function to convert from :math:`A_{ia}` values to :math:`c_{ia}` values, @@ -40,7 +40,7 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, Om_m = cosmo['Omega_m'] rho_crit = physical_constants.RHO_CRITICAL - c1 = c1delta = c2 = ck = None + c1 = c1delta = c2 = ck = ct = None gz = cosmo.growth_factor(1./(1+z)) knorm = 1 # Units of Mpc/h, normalizes units out of ck term @@ -57,8 +57,10 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, c2 = a2*5*5e-14*rho_crit*Om_m/(gz**2) if ak is not None: ck = ak*knorm**2*5e-14*rho_crit*Om_m/gz + if at is not None: + ct=at*5e-14*rho_crit(Om_m)/gz - return c1, c1delta, c2, ck + return c1, c1delta, c2, ck, ct class PTTracer(CCLAutoRepr): @@ -215,7 +217,7 @@ class PTIntrinsicAlignmentTracer(PTTracer): """ type = 'IA' - def __init__(self, c1, c2=None, cdelta=None, ck=None): + def __init__(self, c1, c2=None, cdelta=None, ck=None, ct=None): self.biases = {} @@ -225,8 +227,10 @@ def __init__(self, c1, c2=None, cdelta=None, ck=None): self.biases['c2'] = self._get_bias_function(c2) # Initialize cdelta self.biases['cdelta'] = self._get_bias_function(cdelta) - # Initialize cder + # Initialize ck self.biases['ck'] = self._get_bias_function(ck) + #Initialize ct + self.biases['ct'] = self._get_bias_function(ct) @property def c1(self): @@ -251,3 +255,9 @@ def ck(self): """Internal derivative bias function """ return self.biases['ck'] + + @property + def ct(self): + """Internal velocity bias function + """ + return self.biases['ct'] diff --git a/pyccl/tests/test_ept_power.py b/pyccl/tests/test_ept_power.py index b9cf7703d..ae13ab010 100644 --- a/pyccl/tests/test_ept_power.py +++ b/pyccl/tests/test_ept_power.py @@ -93,7 +93,7 @@ def test_ept_deconstruction(kind): with_matter_1loop=True, cosmo=COSMO, sub_lowk=True) b_nc = ['b1', 'b2', 'b3nl', 'bs', 'bk2'] - b_ia = ['c1', 'c2', 'cdelta', 'ck'] + b_ia = ['c1', 'c2', 'cdelta', 'ck', 'ct'] pk1 = ptc.get_pk2d_template(kind) def get_tr(tn): @@ -111,14 +111,14 @@ def get_tr(tn): bdict[tn] = 1.0 return ccl.nl_pt.PTIntrinsicAlignmentTracer( c1=bdict['c1'], c2=bdict['c2'], - cdelta=bdict['cdelta'], ck=bdict['ck']) + cdelta=bdict['cdelta'], ck=bdict['ck'], ct=bdict['ct']) tn1, tn2 = kind.split(':') t1 = get_tr(tn1) t2 = get_tr(tn2) is_nl = tn1 in ["b2", "bs", "bk2", "b3nl"] - is_g = tn2 in ["c1", "c2", "cdelta", 'ck'] + is_g = tn2 in ["c1", "c2", "cdelta", 'ck', 'ct'] with pytest.warns(ccl.CCLWarning) if is_nl and is_g else nullcontext(): pk2 = ptc.get_biased_pk2d(t1, tracer2=t2) if pk1 is None: From 7fc8039c1a2bfd4a019cb0d216542b613def11c1 Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 7 Dec 2023 14:06:09 -0600 Subject: [PATCH 07/18] Added first galaxy bias terms to pgi tracer --- pyccl/nl_pt/ept.py | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 44fd065f9..f62c6c300 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -287,6 +287,11 @@ def reshape_fastpt(tupl): self.ia_der = self.pt.IA_der(**kw) reshape_fastpt(self.ia_der) self.ia_tij = self.pt.IA_tij(**kw) + reshape_fastpt(self.ia_tij) + self.ia_gb2 = self.pt.IA_gb2(**kw) + reshape_fastpt(self.ia_gb2) + self.ia_s2 = self.pt.IA_s2(**kw) + reshape_fastpt(self.ia_s2) # b1/bk power spectrum pks = {} @@ -407,9 +412,12 @@ def _get_pgi(self, trg, tri): self._check_init() # Get Pk templates Pd1d1 = self.pk_b1 + a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix tijdsij, tij2sij, tijtij, tijsij = self.ia_ta + gb2sij, gb2dsij, gb2sij2, gb2tij = self.ia_gb2 + s2sij, s2dsij, s2sij2, s2tij = self.ia_s2 if(self.ufpt): Pak2 = self.ia_der @@ -438,6 +446,22 @@ def _get_pgi(self, trg, tri): (self._g4*c2)[:, None] * (a0e2 + b0e2) + ck[:, None] * Pak2 + ct[:, None] * tijsij) + pgi = (b1[:,None]*(c1[:, None] * Pd1d1 + + (self._g4*cd)[:, None] * (a00e + c00e) + + (self._g4*c2)[:, None] * (a0e2 + b0e2) + + ck[:, None] * Pak2 + ct[:, None] * tijsij) + + b2[:,None]*(c1[:,None]*gb2sij + + (self._g4*cd)[:,None] * (gb2dsij) + + (self._g4*c2)[:, None] * (gb2sij2) + + ck[:, None] * gb2sij*self.k_s**2 + + ct[:,None] * gb2tij) + + bs[:,None]*(c1[:,None]*s2sij + + (self._g4*cd)[:,None] * (s2dsij) + + (self._g4*c2)[:,None] * (s2sij2) + + ck[:, None] * s2sij*self.k_s**2 + + ct[:, None] *s2tij)) + + return pgi*self.exp_cutoff def _get_pgm(self, trg): From fb6d0a3c7af7c9e8f50e62441d47c44028f788cf Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 30 Apr 2024 14:30:05 -0400 Subject: [PATCH 08/18] Adeed 3nl component to pgi tracer, updated pii tracer to one loop order --- pyccl/nl_pt/ept.py | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index f62c6c300..e3166fe73 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -292,6 +292,8 @@ def reshape_fastpt(tupl): reshape_fastpt(self.ia_gb2) self.ia_s2 = self.pt.IA_s2(**kw) reshape_fastpt(self.ia_s2) + self.ia_one_loop_dd_bias_b3nl = self.pt.one_loop_dd_bias_b3nl(**kw) + reshape_fastpt(self.ia_one_loop_dd_bias_b3nl) # b1/bk power spectrum pks = {} @@ -415,9 +417,10 @@ def _get_pgi(self, trg, tri): a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix - tijdsij, tij2sij, tijtij, tijsij = self.ia_ta + tijdsij, tij2sij, tijtij, tijsij = self.ia_tij gb2sij, gb2dsij, gb2sij2, gb2tij = self.ia_gb2 s2sij, s2dsij, s2sij2, s2tij = self.ia_s2 + d1,d2,d3,d4,d5,d6,d7,d8,sig3nl = self.ia_one_loop_dd_bias_b3nl if(self.ufpt): Pak2 = self.ia_der @@ -441,25 +444,31 @@ def _get_pgi(self, trg, tri): ck = tri.ck(self.z_s) ct = tri.ct(self.z_s) - pgi = b1[:, None] * (c1[:, None] * Pd1d1 + - (self._g4*cd)[:, None] * (a00e + c00e) + - (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:, None] * Pak2 + - ct[:, None] * tijsij) + #pgi = b1[:, None] * (c1[:, None] * Pd1d1 + + # (self._g4*cd)[:, None] * (a00e + c00e) + + # (self._g4*c2)[:, None] * (a0e2 + b0e2) + + # ck[:, None] * Pak2 + + # self._g4*ct[:, None] * tijsij) pgi = (b1[:,None]*(c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:, None] * Pak2 + ct[:, None] * tijsij) + - b2[:,None]*(c1[:,None]*gb2sij + + ck[:, None] * Pak2 + self._g4*ct[:, None] * tijsij) + + b2[:,None]*(self._g4*c1[:,None]*gb2sij + (self._g4*cd)[:,None] * (gb2dsij) + (self._g4*c2)[:, None] * (gb2sij2) + ck[:, None] * gb2sij*self.k_s**2 + - ct[:,None] * gb2tij) + - bs[:,None]*(c1[:,None]*s2sij + + self._g4*ct[:,None] * gb2tij) + + bs[:,None]*(self._g4*c1[:,None]*s2sij + (self._g4*cd)[:,None] * (s2dsij) + (self._g4*c2)[:,None] * (s2sij2) + ck[:, None] * s2sij*self.k_s**2 + - ct[:, None] *s2tij)) + self._g4*ct[:, None] *s2tij) + + bk2[:,None]*self.k_s**2*(c1[:, None] * Pd1d1 + + (self._g4*cd)[:, None] * (a00e + c00e) + + (self._g4*c2)[:, None] * (a0e2 + b0e2) + + ck[:, None] * Pak2 + + self._g4*ct[:, None] * tijsij) + + b3nl[:,None]*(self._g4*c1[:,None]*sig3nl)) return pgi*self.exp_cutoff @@ -555,7 +564,7 @@ def _get_pii(self, tr1, tr2, return_bb=False): (ct1*c22 + ct2*c21)[:,None] * (tij2sij) + (ct1*cd2 + ct2*cd1)[:,None] * (tijdsij) + (ct1*ct2)[:,None] * (tijtij) + - (ct1*ck2 + ct2*ck1)[:,None] * (tijsij)) + (ct1*ck2 + ct2*ck1)[:,None] * (tijsij)*(self.k_s**2)[None, :]) return pii*self.exp_cutoff From c045c0086b7e14438f76bf230069a2db42186ae9 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 7 May 2024 13:30:57 -0400 Subject: [PATCH 09/18] Updated tracers to properly take in prefactors --- pyccl/nl_pt/ept.py | 18 +++++++++--------- pyccl/nl_pt/tracers.py | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index e3166fe73..91103cfd3 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -452,23 +452,23 @@ def _get_pgi(self, trg, tri): pgi = (b1[:,None]*(c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:, None] * Pak2 + self._g4*ct[:, None] * tijsij) + - b2[:,None]*(self._g4*c1[:,None]*gb2sij + + ck[:, None] * Pak2 + (self._g4*ct)[:, None] * tijsij) + + b2[:,None]*((self._g4*c1)[:,None]*gb2sij + (self._g4*cd)[:,None] * (gb2dsij) + (self._g4*c2)[:, None] * (gb2sij2) + - ck[:, None] * gb2sij*self.k_s**2 + - self._g4*ct[:,None] * gb2tij) + - bs[:,None]*(self._g4*c1[:,None]*s2sij + + ck[:, None] * (gb2sij*self.k_s**2) + + (self._g4*ct)[:,None] * gb2tij) + + bs[:,None]*((self._g4*c1)[:,None]*s2sij + (self._g4*cd)[:,None] * (s2dsij) + (self._g4*c2)[:,None] * (s2sij2) + - ck[:, None] * s2sij*self.k_s**2 + - self._g4*ct[:, None] *s2tij) + + ck[:, None] * (s2sij*self.k_s**2) + + (self._g4*ct)[:, None] *s2tij) + bk2[:,None]*self.k_s**2*(c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + ck[:, None] * Pak2 + - self._g4*ct[:, None] * tijsij) + - b3nl[:,None]*(self._g4*c1[:,None]*sig3nl)) + (self._g4*ct)[:, None] * tijsij) + + b3nl[:,None]*((self._g4*c1)[:,None]*sig3nl)) return pgi*self.exp_cutoff diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index 394ae6851..c7dfd1691 100644 --- a/pyccl/nl_pt/tracers.py +++ b/pyccl/nl_pt/tracers.py @@ -56,9 +56,9 @@ def translate_IA_norm(cosmo, *, z, a1=1.0, a1delta=None, a2=None, ak=None, at=No else: # DES convention c2 = a2*5*5e-14*rho_crit*Om_m/(gz**2) if ak is not None: - ck = ak*knorm**2*5e-14*rho_crit*Om_m/gz + ck = ak*(knorm**2)*5e-14*rho_crit*Om_m/gz if at is not None: - ct=at*5e-14*rho_crit(Om_m)/gz + ct=at*5e-14*rho_crit*Om_m/gz return c1, c1delta, c2, ck, ct From 0456c5051ca85595a534c67cd2fd2cdf3eeb43b6 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 19 Nov 2024 02:29:39 -0500 Subject: [PATCH 10/18] Fixing tracers coefficients, removing non one loop terms --- pyccl/nl_pt/ept.py | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 91103cfd3..4f5089d62 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -195,9 +195,6 @@ def __init__(self, *, with_NC=False, with_IA=False, if( not hasattr(fpt, "IA_ta")): raise ValueError(f"Your FAST-PT version lacks a required attribute. You may have the wrong fast-pt install. Try running pip install fast-pt or conda install fast-pt, then try again") - - if( not hasattr(fpt, "IA_der") and self.ufpt): - raise ValueError(f"You need a newer version of FAST-PT to use fpt for k2 term") if (not hasattr(fpt, "IA_tij")): raise ValueError(f"You are using an older version of FAST-PT, please run pip install fast-pt or conda install fast-pt, then try again") n_pad = int(self.fastpt_par['pad_factor'] * len(self.k_s)) @@ -453,23 +450,17 @@ def _get_pgi(self, trg, tri): (self._g4*cd)[:, None] * (a00e + c00e) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + ck[:, None] * Pak2 + (self._g4*ct)[:, None] * tijsij) + - b2[:,None]*((self._g4*c1)[:,None]*gb2sij + + 0.5*b2[:,None]*((self._g4*c1)[:,None]*gb2sij + (self._g4*cd)[:,None] * (gb2dsij) + (self._g4*c2)[:, None] * (gb2sij2) + ck[:, None] * (gb2sij*self.k_s**2) + (self._g4*ct)[:,None] * gb2tij) + - bs[:,None]*((self._g4*c1)[:,None]*s2sij + + 0.5*bs[:,None]*((self._g4*c1)[:,None]*s2sij + (self._g4*cd)[:,None] * (s2dsij) + (self._g4*c2)[:,None] * (s2sij2) + ck[:, None] * (s2sij*self.k_s**2) + (self._g4*ct)[:, None] *s2tij) + - bk2[:,None]*self.k_s**2*(c1[:, None] * Pd1d1 + - (self._g4*cd)[:, None] * (a00e + c00e) + - (self._g4*c2)[:, None] * (a0e2 + b0e2) + - ck[:, None] * Pak2 + - (self._g4*ct)[:, None] * tijsij) + - b3nl[:,None]*((self._g4*c1)[:,None]*sig3nl)) - + 0.5*b3nl[:,None]*((self._g4*c1)[:,None]*sig3nl)) return pgi*self.exp_cutoff @@ -563,8 +554,7 @@ def _get_pii(self, tr1, tr2, return_bb=False): (ct1*c12 + ct2*c11)[:,None]*(tijsij) + (ct1*c22 + ct2*c21)[:,None] * (tij2sij) + (ct1*cd2 + ct2*cd1)[:,None] * (tijdsij) + - (ct1*ct2)[:,None] * (tijtij) + - (ct1*ck2 + ct2*ck1)[:,None] * (tijsij)*(self.k_s**2)[None, :]) + (ct1*ct2)[:,None] * (tijtij)) return pii*self.exp_cutoff From 5e9fa0cbc3140dc0f1d22a6711d229f18dc03c0f Mon Sep 17 00:00:00 2001 From: Vincent Schacknies <136936124+vschac@users.noreply.github.com> Date: Thu, 23 Jan 2025 15:15:21 -0500 Subject: [PATCH 11/18] new PK_ALIAS terms --- pyccl/nl_pt/ept.py | 253 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 190 insertions(+), 63 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 4b42b7e5a..ce190d1cd 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -11,19 +11,21 @@ 'm:m': 'm:m', 'm:b1': 'm:m', 'm:b2': 'm:b2', 'm:b3nl': 'm:b3nl', 'm:bs': 'm:bs', 'm:bk2': 'm:bk2', 'm:c1': 'm:m', 'm:c2': 'm:c2', 'm:cdelta': 'm:cdelta', - 'b1:b1': 'm:m', 'b1:b2': 'm:b2', 'b1:b3nl': 'm:b3nl', - 'b1:bs': 'm:bs', 'b1:bk2': 'm:bk2', 'b1:c1': 'm:m', - 'b1:c2': 'm:c2', 'b1:cdelta': 'm:cdelta', 'b2:b2': 'b2:b2', - 'b2:b3nl': 'zero', 'b2:bs': 'b2:bs', 'b2:bk2': 'zero', - 'b2:c1': 'zero', 'b2:c2': 'zero', 'b2:cdelta': 'zero', - 'b3nl:b3nl': 'zero', 'b3nl:bs': 'zero', - 'b3nl:bk2': 'zero', 'b3nl:c1': 'zero', 'b3nl:c2': - 'zero', 'b3nl:cdelta': 'zero', 'bs:bs': 'bs:bs', - 'bs:bk2': 'zero', 'bs:c1': 'zero', 'bs:c2': 'zero', - 'bs:cdelta': 'zero', 'bk2:bk2': 'zero', 'bk2:c1': 'zero', + 'm:ck': 'm:ck', 'b1:b1': 'm:m', 'b1:b2': 'm:b2', + 'b1:b3nl': 'm:b3nl', 'b1:bs': 'm:bs', 'b1:bk2': 'm:bk2', + 'b1:c1': 'm:m', 'b1:c2': 'm:c2', 'b1:cdelta': 'm:cdelta', 'b1:ck': 'm:ck', + 'b2:b2': 'b2:b2', 'b2:bs': 'b2:bs', + 'b2:bk2': 'zero', 'b2:c1': 'b2:m', 'b2:c2': 'b2:c2', + 'b2:cdelta': 'b2:cdelta','b3nl:c1': 'b3nl:m','bs:bs': 'bs:bs', + 'bs:bk2': 'zero', 'c1:bs': 'm:bs', 'bs:c2': 'bs:c2', + 'bs:cdelta': 'bs:cdelta', 'bk2:bk2': 'zero', 'c1:bk2': 'm:bk2', 'bk2:c2': 'zero', 'bk2:cdelta': 'zero', 'c1:c1': 'm:m', - 'c1:c2': 'm:c2', 'c1:cdelta': 'm:cdelta', 'c2:c2': 'c2:c2', - 'c2:cdelta': 'c2:cdelta', 'cdelta:cdelta': 'cdelta:cdelta'} + 'c1:c2': 'm:c2', 'c1:cdelta': 'm:cdelta', 'c1:ck': 'm:ck', + 'c2:c2': 'c2:c2', 'c2:cdelta': 'c2:cdelta', + 'cdelta:cdelta': 'cdelta:cdelta', 'ck:ck': 'zero', 'm:ct': 'm:ct', + 'b1:ct': 'm:ct', 'c1:ct': 'm:ct','c2:ct': 'c2:ct', + 'cdelta:ct': 'cdelta:ct', 'ct:ct': 'ct:ct', 'ck:ct': 'zero'} + class EulerianPTCalculator(CCLAutoRepr): @@ -45,7 +47,7 @@ class EulerianPTCalculator(CCLAutoRepr): .. math:: s^I_{ij}=c_1\\,s_{ij}+c_2(s_{ik}s_{jk}-s^2\\delta_{ik}/3) - +c_\\delta\\,\\delta\\,s_{ij} + +c_\\delta\\,\\delta\\,s_{ij + c_k\\k^2\\,s_{ij}} (note that the higher-order terms are not divided by 2!). @@ -113,6 +115,9 @@ class EulerianPTCalculator(CCLAutoRepr): bk2_pk_kind (:obj:`str`): power spectrum to use for the non-local bias terms in the expansion. Same options and default as ``b1_pk_kind``. + ak2_pk_kind (:obj:`str`): power spectrum to use for the derivative + term of the IA expansion. Same options and default as + ``b1_pk_kind``. pad_factor (:obj:`float`): fraction of the :math:`\\log_{10}(k)` interval you to add as padding for FFTLog calculations. low_extrap (:obj:`float`): decimal logaritm of the minimum Fourier @@ -129,6 +134,8 @@ class EulerianPTCalculator(CCLAutoRepr): documentation for more details. sub_lowk (:obj:`bool`): if ``True``, the small-scale white noise contribution to some of the terms will be subtracted. + usefptk (::obj:`bool`):) if `True``, will use the FASTPT IA k2 + term instead of the CCL IA k2 term """ __repr_attrs__ = __eq_attrs__ = ('with_NC', 'with_IA', 'with_matter_1loop', 'k_s', 'a_s', 'exp_cutoff', 'b1_pk_kind', @@ -139,12 +146,13 @@ def __init__(self, *, with_NC=False, with_IA=False, log10k_min=-4, log10k_max=2, nk_per_decade=20, a_arr=None, k_cutoff=None, n_exp_cutoff=4, b1_pk_kind='nonlinear', bk2_pk_kind='nonlinear', + ak2_pk_kind='nonlinear', pad_factor=1.0, low_extrap=-5.0, high_extrap=3.0, - P_window=None, C_window=0.75, sub_lowk=False): + P_window=None, C_window=0.75, sub_lowk=False, usefptk=False): self.with_matter_1loop = with_matter_1loop self.with_NC = with_NC self.with_IA = with_IA - + self.ufpt=usefptk # Set FAST-PT parameters self.fastpt_par = {'pad_factor': pad_factor, 'low_extrap': low_extrap, @@ -177,7 +185,16 @@ def __init__(self, *, with_NC=False, with_IA=False, self.exp_cutoff = 1 # Call FAST-PT - import fastpt as fpt + try: + import fastpt as fpt + except: + raise ImportError("Your attempted import of FAST-PT has failed. You either dont have fast-pt installed, or have the wrong version. Try running pip install fast-pt or conda install fast-pt, then try again") + + if( not hasattr(fpt, "IA_ta")): + raise ValueError(f"Your FAST-PT version lacks a required attribute. You may have the wrong fast-pt install. Try running pip install fast-pt or conda install fast-pt, then try again") + if (not hasattr(fpt, "IA_tij")): + pass + #raise ValueError(f"You are using an older version of FAST-PT, please run pip install fast-pt or conda install fast-pt, then try again") n_pad = int(self.fastpt_par['pad_factor'] * len(self.k_s)) self.pt = fpt.FASTPT(self.k_s, to_do=to_do, low_extrap=self.fastpt_par['low_extrap'], @@ -189,8 +206,11 @@ def __init__(self, *, with_NC=False, with_IA=False, raise ValueError(f"Unknown P(k) prescription {b1_pk_kind}") if bk2_pk_kind not in ['linear', 'nonlinear', 'pt']: raise ValueError(f"Unknown P(k) prescription {bk2_pk_kind}") + if ak2_pk_kind not in ['linear', 'nonlinear', 'pt']: + raise ValueError(f"Unknown P(k) prescription in {ak2_pk_kind}") self.b1_pk_kind = b1_pk_kind self.bk2_pk_kind = bk2_pk_kind + self.ak2_pk_kind = ak2_pk_kind if (self.b1_pk_kind == 'pt') or (self.bk2_pk_kind == 'pt'): self.with_matter_1loop = True @@ -258,6 +278,19 @@ def reshape_fastpt(tupl): reshape_fastpt(self.ia_tt) self.ia_mix = self.pt.IA_mix(**kw) reshape_fastpt(self.ia_mix) + if(self.ufpt): + self.ia_der = self.pt.IA_der(**kw) + reshape_fastpt(self.ia_der) + self.ia_ct = self.pt.IA_ct(**kw) + reshape_fastpt(self.ia_ct) + self.ia_ctbias = self.pt.IA_ctbias(**kw) + reshape_fastpt(self.ia_ctbias) + self.ia_d2 = self.pt.IA_d2(**kw) + reshape_fastpt(self.ia_d2) + self.ia_s2 = self.pt.IA_s2(**kw) + reshape_fastpt(self.ia_s2) + self.ia_one_loop_dd_bias_b3nl = self.pt.one_loop_dd_bias_b3nl(**kw) + reshape_fastpt(self.ia_one_loop_dd_bias_b3nl) # b1/bk power spectrum pks = {} @@ -279,6 +312,24 @@ def reshape_fastpt(tupl): self.pk_b1 = pks[self.b1_pk_kind] self.pk_bk = pks[self.bk2_pk_kind] + # ak power spectrum + pksa = {} + if 'nonlinear' in [self.ak2_pk_kind]: + pksa['nonlinear'] = np.array([cosmo.nonlin_matter_power(self.k_s, a) + for a in self.a_s]) + if 'linear' in [self.ak2_pk_kind]: + pksa['linear'] = np.array([cosmo.linear_matter_power(self.k_s, a) + for a in self.a_s]) + if 'pt' in [self.ak2_pk_kind]: + if 'linear' in pksa: + pka = pksa['linear'] + else: + pka = np.array([cosmo.linear_matter_power(self.k_s, a) for a in self.a_s]) + pka += self._g4T*self.one_loop_dd[0] + pksa['pt'] = pka + self.pk_ak=pksa[self.ak2_pk_kind] + + # Reset template power spectra self._pk2d_temp = {} self._cosmo = cosmo @@ -360,8 +411,20 @@ def _get_pgi(self, trg, tri): self._check_init() # Get Pk templates Pd1d1 = self.pk_b1 + a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + tijsij, tijdsij, tij2sij, tijtij = self.ia_ct + gb2tij, s2tij = self.ia_ctbias + gb2sij, gb2dsij, gb2sij2 = self.ia_d2 + s2sij, s2dsij, s2sij2= self.ia_s2 + d1,d2,d3,d4,d5,d6,d7,d8,sig3nl = self.ia_one_loop_dd_bias_b3nl + Pd1s2 = self._g4T * self.dd_bias[4] + + if(self.ufpt): + Pak2 = self.ia_der + else: + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases b1 = trg.b1(self.z_s) @@ -377,10 +440,31 @@ def _get_pgi(self, trg, tri): c1 = tri.c1(self.z_s) c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) - - pgi = b1[:, None] * (c1[:, None] * Pd1d1 + + ck = tri.ck(self.z_s) + ct = tri.ct(self.z_s) + + #pgi = b1[:, None] * (c1[:, None] * Pd1d1 + + # (self._g4*cd)[:, None] * (a00e + c00e) + + # (self._g4*c2)[:, None] * (a0e2 + b0e2) + + # ck[:, None] * Pak2 + + # self._g4*ct[:, None] * tijsij) + pgi = (b1[:,None]*(c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + - (self._g4*c2)[:, None] * (a0e2 + b0e2)) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + + ck[:, None] * Pak2 + (self._g4*ct)[:, None] * tijsij) + + 0.5*b2[:,None]*((self._g4*c1)[:,None]*gb2sij + + (self._g4*cd)[:,None] * (gb2dsij) + + (self._g4*c2)[:, None] * (gb2sij2) + + ck[:, None] * (gb2sij*self.k_s**2) + + (self._g4*ct)[:,None] * gb2tij) + + 0.5*bs[:,None]*((self._g4*c1)[:,None]*Pd1s2 + + (self._g4*cd)[:,None] * (s2dsij) + + (self._g4*c2)[:,None] * (s2sij2) + + ck[:, None] * (s2sij*self.k_s**2) + + (self._g4*ct)[:, None] *s2tij) + + 0.5*b3nl[:,None]*((self._g4*c1)[:,None]*sig3nl) + + 0.5*bk2[:,None]*((self._g4*c1)[:,None]*(Pd1d1*self.k_s**2))) + return pgi*self.exp_cutoff def _get_pgm(self, trg): @@ -440,14 +524,23 @@ def _get_pii(self, tr1, tr2, return_bb=False): a00e, c00e, a0e0e, a0b0b = self.ia_ta ae2e2, ab2b2 = self.ia_tt a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + tijsij, tijdsij, tij2sij, tijtij = self.ia_ct + if(self.ufpt): + Pak2 = self.ia_der + else: + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases c11 = tr1.c1(self.z_s) c21 = tr1.c2(self.z_s) cd1 = tr1.cdelta(self.z_s) + ck1 = tr1.ck(self.z_s) + ct1 = tr1.ct(self.z_s) c12 = tr2.c1(self.z_s) c22 = tr2.c2(self.z_s) cd2 = tr2.cdelta(self.z_s) + ck2 = tr2.ck(self.z_s) + ct2 = tr2.ct(self.z_s) if return_bb: pii = ((cd1*cd2*self._g4)[:, None]*a0b0b + @@ -459,7 +552,12 @@ def _get_pii(self, tr1, tr2, return_bb=False): (cd1*cd2*self._g4)[:, None]*a0e0e + (c21*c22*self._g4)[:, None]*ae2e2 + ((c11*c22+c21*c12)*self._g4)[:, None]*(a0e2+b0e2) + - ((cd1*c22+cd2*c21)*self._g4)[:, None]*d0ee2) + ((cd1*c22+cd2*c21)*self._g4)[:, None]*d0ee2 + + (ck1*c12 + ck2*c11)[:,None] * (Pak2) + + (ct1*c12 + ct2*c11)[:,None]*(tijsij) + + (ct1*c22 + ct2*c21)[:,None] * (tij2sij) + + (ct1*cd2 + ct2*cd1)[:,None] * (tijdsij) + + (ct1*ct2)[:,None] * (tijtij)) return pii*self.exp_cutoff @@ -481,15 +579,24 @@ def _get_pim(self, tri): Pd1d1 = self.pk_b1 a00e, c00e, a0e0e, a0b0b = self.ia_ta a0e2, b0e2, d0ee2, d0bb2 = self.ia_mix + tijsij, tijdsij, tij2sij, tijtij = self.ia_ct + if(self.ufpt): + Pak2 = self.ia_der + else: + Pak2 = self.pk_ak*(self.k_s**2)[None, :] # Get biases c1 = tri.c1(self.z_s) c2 = tri.c2(self.z_s) cd = tri.cdelta(self.z_s) + ck = tri.ck(self.z_s) + ct = tri.ct(self.z_s) pim = (c1[:, None] * Pd1d1 + (self._g4*cd)[:, None] * (a00e + c00e) + - (self._g4*c2)[:, None] * (a0e2 + b0e2)) + (self._g4*c2)[:, None] * (a0e2 + b0e2) + + ck[:,None] * Pak2 + + ct[:,None] * tijsij) return pim*self.exp_cutoff def _get_pmm(self): @@ -640,48 +747,68 @@ def get_pk2d_template(self, kind, *, extrap_order_lok=1, # Construct power spectrum array s4 = 0. - if pk_name == 'm:m': - pk = self.pk_b1 - elif pk_name == 'm:b2': - pk = 0.5*self._g4T*self.dd_bias[2] - elif pk_name == 'm:b3nl': - pk = 0.5*self._g4T*self.dd_bias[8] - elif pk_name == 'm:bs': - pk = 0.5*self._g4T*self.dd_bias[4] - elif pk_name == 'm:bk2': - pk = 0.5*self.pk_bk*(self.k_s**2) - elif pk_name == 'm:c2': - pk = self._g4T * (self.ia_mix[0]+self.ia_mix[1]) - elif pk_name == 'm:cdelta': - pk = self._g4T * (self.ia_ta[0]+self.ia_ta[1]) - elif pk_name == 'b2:b2': - if self.fastpt_par['sub_lowk']: - s4 = self.dd_bias[7] - pk = 0.25*self._g4T*(self.dd_bias[3] - 2*s4) - elif pk_name == 'b2:bs': - if self.fastpt_par['sub_lowk']: - s4 = self.dd_bias[7] - pk = 0.25*self._g4T*(self.dd_bias[5] - 4*s4/3) - elif pk_name == 'bs:bs': - if self.fastpt_par['sub_lowk']: - s4 = self.dd_bias[7] - pk = 0.25*self._g4T*(self.dd_bias[6] - 8*s4/9) - elif pk_name == 'c2:c2': - pk = self._g4T * self.ia_tt[0] - elif pk_name == 'c2:c2_bb': - pk = self._g4T * self.ia_tt[1] - elif pk_name == 'c2:cdelta': - pk = self._g4T * self.ia_mix[2] - elif pk_name == 'c2:cdelta_bb': - pk = self._g4T * self.ia_mix[3] - elif pk_name == 'cdelta:cdelta': - pk = self._g4T * self.ia_ta[2] - elif pk_name == 'cdelta:cdelta_bb': - pk = self._g4T * self.ia_ta[3] - elif pk_name == 'zero': - # If zero, store None and return - self._pk2d_temp[pk_name] = None - return None + match pk_name: + case 'm:m': + pk = self.pk_b1 + case 'm:b2'| 'b2:m': + pk = 0.5 * self._g4T * self.dd_bias[2] + case 'm:b3nl'| 'b3nl:m': + pk = 0.5 * self._g4T * self.dd_bias[8] + case 'm:bs'| 'bs:m': + pk = 0.5 * self._g4T * self.dd_bias[4] + case 'm:bk2'| 'bk2:m': + pk = 0.5 * self.pk_bk * (self.k_s**2) + case 'm:c2': + pk = self._g4T * (self.ia_mix[0] + self.ia_mix[1]) + case 'm:cdelta': + pk = self._g4T * (self.ia_ta[0] + self.ia_ta[1]) + case 'm:ck': + pk = self.pk_ak * (self.k_s**2) + case 'b2:b2': + if self.fastpt_par['sub_lowk']: + s4 = self.dd_bias[7] + pk = 0.25 * self._g4T * (self.dd_bias[3] - 2 * s4) + case 'b2:bs': + if self.fastpt_par['sub_lowk']: + s4 = self.dd_bias[7] + pk = 0.25 * self._g4T * (self.dd_bias[5] - 4 * s4 / 3) + case 'bs:bs': + if self.fastpt_par['sub_lowk']: + s4 = self.dd_bias[7] + pk = 0.25 * self._g4T * (self.dd_bias[6] - 8 * s4 / 9) + case 'c2:c2': + pk = self._g4T * self.ia_tt[0] + case 'c2:c2_bb': + pk = self._g4T * self.ia_tt[1] + case 'c2:cdelta': + pk = self._g4T * self.ia_mix[2] + case 'c2:cdelta_bb': + pk = self._g4T * self.ia_mix[3] + case 'cdelta:cdelta': + pk = self._g4T * self.ia_ta[2] + case 'cdelta:cdelta_bb': + pk = self._g4T * self.ia_ta[3] + case 'm:ct': + pk = self._g4T * self.ia_ct[0] + case 'c2:ct': + pk = self._g4T * self.ia_ct[2] + case 'cdelta:ct': + pk = self._g4T * self.ia_ct[1] + case 'ct:ct': + pk = self._g4T * self.ia_ct[3] + case 'b2:c2': + pk = 0.5*self._g4T * self.ia_d2[2] + case 'b2:cdelta': + pk = 0.5*self._g4T * self.ia_d2[1] + case 'bs:cdelta': + pk = 0.5*self._g4T * self.ia_s2[1] + case 'bs:c2': + pk = 0.5*self._g4T * self.ia_s2[2] + case 'zero': + # If zero, store None and return + self._pk2d_temp[pk_name] = None + return None + # Build interpolator pk2d = Pk2D(a_arr=self.a_s, From e5850278eec700c932a5b6efb321d0b15785e960 Mon Sep 17 00:00:00 2001 From: Vincent Schacknies <136936124+vschac@users.noreply.github.com> Date: Thu, 23 Jan 2025 15:16:07 -0500 Subject: [PATCH 12/18] fixed warning verbosity --- pyccl/tests/test_ept_power.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/pyccl/tests/test_ept_power.py b/pyccl/tests/test_ept_power.py index 53b0741aa..42f2257ee 100644 --- a/pyccl/tests/test_ept_power.py +++ b/pyccl/tests/test_ept_power.py @@ -3,13 +3,14 @@ import pytest from contextlib import nullcontext + COSMO = ccl.Cosmology( Omega_c=0.27, Omega_b=0.045, h=0.67, sigma8=0.8, n_s=0.96, transfer_function='bbks', matter_power_spectrum='linear') TRS = {'TG': ccl.nl_pt.PTNumberCountsTracer(b1=2.0, b2=2.0, bs=2.0), 'TI': ccl.nl_pt.PTIntrinsicAlignmentTracer(c1=2.0, c2=2.0, - cdelta=2.0), + cdelta=2.0, ck=2.0), 'TM': ccl.nl_pt.PTMatterTracer()} PTC = ccl.nl_pt.EulerianPTCalculator(with_NC=True, with_IA=True, with_matter_1loop=True, @@ -60,7 +61,7 @@ def test_ept_pk2d_bb_smoke(): def test_ept_get_pk2d_nl(nl): ptc = ccl.nl_pt.EulerianPTCalculator( with_NC=True, with_IA=True, with_matter_1loop=True, - b1_pk_kind=nl, bk2_pk_kind=nl, cosmo=COSMO) + b1_pk_kind=nl, bk2_pk_kind=nl, ak2_pk_kind=nl,cosmo=COSMO) pk = ptc.get_biased_pk2d(TRS['TG']) assert isinstance(pk, ccl.Pk2D) @@ -74,7 +75,7 @@ def test_ept_k2pk_types(typ_nlin, typ_nloc): tm = ccl.nl_pt.PTNumberCountsTracer(1., 0., 0.) ptc1 = ccl.nl_pt.EulerianPTCalculator( with_NC=True, with_IA=True, with_matter_1loop=True, - b1_pk_kind=typ_nlin, bk2_pk_kind=typ_nloc, + b1_pk_kind=typ_nlin, bk2_pk_kind=typ_nloc, ak2_pk_kind=typ_nloc, cosmo=COSMO) ptc2 = ccl.nl_pt.EulerianPTCalculator( with_NC=True, with_IA=True, with_matter_1loop=True, @@ -94,7 +95,7 @@ def test_ept_deconstruction(kind): with_matter_1loop=True, cosmo=COSMO, sub_lowk=True) b_nc = ['b1', 'b2', 'b3nl', 'bs', 'bk2'] - b_ia = ['c1', 'c2', 'cdelta'] + b_ia = ['c1', 'c2', 'cdelta', 'ck', 'ct'] pk1 = ptc.get_pk2d_template(kind) def get_tr(tn): @@ -112,23 +113,31 @@ def get_tr(tn): bdict[tn] = 1.0 return ccl.nl_pt.PTIntrinsicAlignmentTracer( c1=bdict['c1'], c2=bdict['c2'], - cdelta=bdict['cdelta']) + cdelta=bdict['cdelta'], ck=bdict['ck'], ct=bdict['ct']) tn1, tn2 = kind.split(':') t1 = get_tr(tn1) t2 = get_tr(tn2) - + is_nl = tn1 in ["b2", "bs", "bk2", "b3nl"] - is_g = tn2 in ["c1", "c2", "cdelta"] + is_g = tn2 in ["c1", "c2", "cdelta", 'ck', 'ct'] ccl.update_warning_verbosity('high') with pytest.warns(ccl.CCLWarning) if is_nl and is_g else nullcontext(): pk2 = ptc.get_biased_pk2d(t1, tracer2=t2) ccl.update_warning_verbosity('low') + if pk1 is None: assert pk2(0.5, 1.0, cosmo=COSMO) == 0.0 else: v1 = pk1(0.5, 1.0, cosmo=COSMO) v2 = pk2(0.5, 1.0, cosmo=COSMO) + if not np.allclose(v1, v2, atol=0, rtol=1e-6): + pytest.fail(f""" + Failed comparison for {kind} + v1 (template): {v1} + v2 (biased): {v2} + relative diff: {abs(v1-v2)/abs(v1)} + """) assert np.allclose(v1, v2, atol=0, rtol=1e-6) # Check cached pk3 = ptc._pk2d_temp[ccl.nl_pt.ept._PK_ALIAS[kind]] @@ -139,7 +148,7 @@ def get_tr(tn): @pytest.mark.parametrize('kind', ['c2:c2', 'c2:cdelta', 'cdelta:cdelta']) def test_ept_deconstruction_bb(kind): - b_ia = ['c1', 'c2', 'cdelta'] + b_ia = ['c1', 'c2', 'cdelta', 'ck'] pk1 = PTC.get_pk2d_template(kind, return_ia_bb=True) def get_tr(tn): @@ -147,7 +156,7 @@ def get_tr(tn): bdict[tn] = 1.0 return ccl.nl_pt.PTIntrinsicAlignmentTracer( c1=bdict['c1'], c2=bdict['c2'], - cdelta=bdict['cdelta']) + cdelta=bdict['cdelta'], ck=bdict['ck']) tn1, tn2 = kind.split(':') t1 = get_tr(tn1) @@ -217,6 +226,10 @@ def test_ept_calculator_raises(): # Wrong type of b1 and bk2 power spectra with pytest.raises(ValueError): ccl.nl_pt.EulerianPTCalculator(bk2_pk_kind='non-linear') + + #wrong type of ak2 power spectra + with pytest.raises(ValueError): + ccl.nl_pt.EulerianPTCalculator(ak2_pk_kind="non-linear") # Uninitialized templates with pytest.raises(ccl.CCLError): From 49ad959a3a7850e6450ad987e568261028b8ffbb Mon Sep 17 00:00:00 2001 From: Vincent Schacknies Date: Mon, 27 Jan 2025 14:39:07 -0500 Subject: [PATCH 13/18] updated benchmarking --- benchmarks/data/codes/pt_bm.py | 54 ++++++++++++++++++++++++++++++---- 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/benchmarks/data/codes/pt_bm.py b/benchmarks/data/codes/pt_bm.py index 6015a3e4d..e12f35867 100644 --- a/benchmarks/data/codes/pt_bm.py +++ b/benchmarks/data/codes/pt_bm.py @@ -39,6 +39,20 @@ P_window=P_window, C_window=C_window) +#New IA terms +ia_ct = pt_ob.IA_ct(pk, + P_window=P_window, + C_window=C_window) +ia_ctbias = pt_ob.IA_ctbias(pk, + P_window=P_window, + C_window=C_window) +ia_d2 = pt_ob.IA_d2(pk, + P_window=P_window, + C_window=C_window) +ia_s2 = pt_ob.IA_s2(pk, + P_window=P_window, + C_window=C_window) + g4 = g4[:, None] Pd1d2 = g4 * dd_bias[2][None, :] Pd2d2 = g4 * dd_bias[3][None, :] @@ -57,7 +71,23 @@ b0e2 = g4 * ia_mix[1][None, :] d0ee2 = g4 * ia_mix[2][None, :] d0bb2 = g4 * ia_mix[3][None, :] - +#New terms +tijsij = g4 * ia_ct[0][None, :] +tijdsij = g4 * ia_ct[1][None, :] +tij2sij = g4 * ia_ct[2][None, :] +tijtij = g4 * ia_ct[3][None, :] + +gb2tij = g4 * ia_ctbias[0][None, :] +s2tij = g4 * ia_ctbias[1][None, :] + +gb2sij = g4 * ia_d2[0][None, :] +gb2dsij = g4 * ia_d2[1][None, :] +gb2sij2 = g4 * ia_d2[2][None, :] + +s2sij = g4 * ia_s2[0][None, :] +s2dsij = g4 * ia_s2[1][None, :] +s2sij2 = g4 * ia_s2[2][None, :] + b1 = 1.3 b2 = 1.5 bs = 1.7 @@ -66,6 +96,7 @@ c1 = 1.9 c2 = 2.1 cd = 2.3 +ct = 2.5 pgg = (b1**2 * Pd1d1 + b1*b2 * Pd1d2 + @@ -81,9 +112,6 @@ 0.5 * bs * Pd1s2 + 0.5 * b3 * Pd1d3 + 0.5 * bk2 * Pd1k2) -pgi = b1 * (c1 * Pd1d1 + - cd * (a00e + c00e) + - c2 * (a0e2 + b0e2)) pii = (c1**2 * Pd1d1 + 2 * c1 * cd * (a00e + c00e) + cd**2 * a0e0e + @@ -96,12 +124,26 @@ pim = (c1 * Pd1d1 + cd * (a00e + c00e) + c2 * (a0e2 + b0e2)) +#Updated pgi +pgi = b1 * (c1 * Pd1d1 + + g4 * cd * (a00e + c00e) + + g4 * c2 * (a0e2 + b0e2) + + bk2 * Pd1k2 + + ct * tijsij + + 0.5*b2*((c1*gb2sij) + + (cd*gb2dsij) + + (c2*gb2sij2) + + (ct*gb2tij)) + + 0.5*bs*((c1*s2sij) + + (cd*s2dsij) + + (c2*s2sij2) + + (ct*s2tij))) -np.savetxt("../pt_bm_z0.txt", +np.savetxt("../pt_bm_test.txt", np.transpose([ks, pgg[0], pgm[0], pgi[0], pii[0], pii_bb[0], pim[0]]), header='[0]-k [1]-GG [2]-GM [3]-GI [4]-II [5]-II_BB [6]-IM') -np.savetxt("../pt_bm_z1.txt", +np.savetxt("../pt_bm_test2.txt", np.transpose([ks, pgg[1], pgm[1], pgi[1], pii[1], pii_bb[1], pim[1]]), header='[0]-k [1]-GG [2]-GM [3]-GI [4]-II [5]-II_BB [6]-IM') From 393543daa2ce8839cb0c818ee461a148fe6b5955 Mon Sep 17 00:00:00 2001 From: Vincent Schacknies Date: Tue, 28 Jan 2025 14:39:24 -0500 Subject: [PATCH 14/18] updated benchmarking --- benchmarks/data/codes/pt_bm.py | 18 ++++++++++++------ pyccl/nl_pt/ept.py | 4 +--- 2 files changed, 13 insertions(+), 9 deletions(-) diff --git a/benchmarks/data/codes/pt_bm.py b/benchmarks/data/codes/pt_bm.py index e12f35867..3dbb82219 100644 --- a/benchmarks/data/codes/pt_bm.py +++ b/benchmarks/data/codes/pt_bm.py @@ -97,6 +97,7 @@ c2 = 2.1 cd = 2.3 ct = 2.5 +ck2 = 0.1 pgg = (b1**2 * Pd1d1 + b1*b2 * Pd1d2 + @@ -114,21 +115,23 @@ 0.5 * bk2 * Pd1k2) pii = (c1**2 * Pd1d1 + 2 * c1 * cd * (a00e + c00e) + - cd**2 * a0e0e + - c2**2 * ae2e2 + + 2 * cd**2 * a0e0e + #maybe multiple by 2? + 2 * c2**2 * ae2e2 + #maybe multiple by 2? 2 * c1 * c2 * (a0e2 + b0e2) + - 2 * cd * c2 * d0ee2) + (2 * cd * c2 * d0ee2) + + 2 * c1 * ck2 * Pd1k2) pii_bb = (cd**2 * a0b0b + c2**2 * ab2b2 + 2 * cd * c2 * d0bb2) pim = (c1 * Pd1d1 + cd * (a00e + c00e) + - c2 * (a0e2 + b0e2)) + c2 * (a0e2 + b0e2) + + ck2 * Pd1k2) #Updated pgi pgi = b1 * (c1 * Pd1d1 + g4 * cd * (a00e + c00e) + g4 * c2 * (a0e2 + b0e2) + - bk2 * Pd1k2 + + ck2 * Pd1k2 + ct * tijsij + 0.5*b2*((c1*gb2sij) + (cd*gb2dsij) + @@ -137,8 +140,11 @@ 0.5*bs*((c1*s2sij) + (cd*s2dsij) + (c2*s2sij2) + - (ct*s2tij))) + (ct*s2tij)) + + (0.5 * b3nl * c1 * sig3nl) + + (0.5 * bk2 * c1 * Pd1k2)) +#Need to change these names and override old benchmark text file np.savetxt("../pt_bm_test.txt", np.transpose([ks, pgg[0], pgm[0], pgi[0], pii[0], pii_bb[0], pim[0]]), diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 97b04766a..ef8803f3f 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -453,13 +453,11 @@ def _get_pgi(self, trg, tri): ck[:, None] * Pak2 + (self._g4*ct)[:, None] * tijsij) + 0.5*b2[:,None]*((self._g4*c1)[:,None]*gb2sij + (self._g4*cd)[:,None] * (gb2dsij) + - (self._g4*c2)[:, None] * (gb2sij2) + - ck[:, None] * (gb2sij*self.k_s**2) + + (self._g4*c2)[:, None] * (gb2sij2) + (self._g4*ct)[:,None] * gb2tij) + 0.5*bs[:,None]*((self._g4*c1)[:,None]*Pd1s2 + (self._g4*cd)[:,None] * (s2dsij) + (self._g4*c2)[:,None] * (s2sij2) + - ck[:, None] * (s2sij*self.k_s**2) + (self._g4*ct)[:, None] *s2tij) + 0.5*b3nl[:,None]*((self._g4*c1)[:,None]*sig3nl) + 0.5*bk2[:,None]*((self._g4*c1)[:,None]*(Pd1d1*self.k_s**2))) From bebf0f7200a8d99667b621c143c978cb389508c6 Mon Sep 17 00:00:00 2001 From: Vincent Schacknies Date: Tue, 28 Jan 2025 16:03:58 -0500 Subject: [PATCH 15/18] change in get_pgi --- pyccl/nl_pt/ept.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index ef8803f3f..4e1721a41 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -193,7 +193,8 @@ def __init__(self, *, with_NC=False, with_IA=False, if( not hasattr(fpt, "IA_ta")): raise ValueError(f"Your FAST-PT version lacks a required attribute. You may have the wrong fast-pt install. Try running pip install fast-pt or conda install fast-pt, then try again") if (not hasattr(fpt, "IA_tij")): - raise ValueError(f"You are using an older version of FAST-PT, please run pip install fast-pt or conda install fast-pt, then try again") + pass + #raise ValueError(f"You are using an older version of FAST-PT, please run pip install fast-pt or conda install fast-pt, then try again") n_pad = int(self.fastpt_par['pad_factor'] * len(self.k_s)) self.pt = fpt.FASTPT(self.k_s, to_do=to_do, low_extrap=self.fastpt_par['low_extrap'], @@ -455,7 +456,7 @@ def _get_pgi(self, trg, tri): (self._g4*cd)[:,None] * (gb2dsij) + (self._g4*c2)[:, None] * (gb2sij2) + (self._g4*ct)[:,None] * gb2tij) + - 0.5*bs[:,None]*((self._g4*c1)[:,None]*Pd1s2 + + 0.5*bs[:,None]*((self._g4*c1)[:,None]*Pd1s2+ #other term we tried: s2sij + (self._g4*cd)[:,None] * (s2dsij) + (self._g4*c2)[:,None] * (s2sij2) + (self._g4*ct)[:, None] *s2tij) + From a53fbd9a0276960e55f8e5f31e1c2e8992967d86 Mon Sep 17 00:00:00 2001 From: Vincent Schacknies Date: Fri, 31 Jan 2025 12:01:40 -0500 Subject: [PATCH 16/18] fixed benchmarking --- benchmarks/data/codes/pt_bm.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/benchmarks/data/codes/pt_bm.py b/benchmarks/data/codes/pt_bm.py index 3dbb82219..de61c223b 100644 --- a/benchmarks/data/codes/pt_bm.py +++ b/benchmarks/data/codes/pt_bm.py @@ -53,6 +53,11 @@ P_window=P_window, C_window=C_window) + +one_loop_dd_bias_b3nl = pt_ob.one_loop_dd_bias_b3nl(pk, + P_window=P_window, + C_window=C_window) + g4 = g4[:, None] Pd1d2 = g4 * dd_bias[2][None, :] Pd2d2 = g4 * dd_bias[3][None, :] @@ -88,6 +93,8 @@ s2dsij = g4 * ia_s2[1][None, :] s2sij2 = g4 * ia_s2[2][None, :] +sig3nl = g4 * one_loop_dd_bias_b3nl[8][None, :] + b1 = 1.3 b2 = 1.5 bs = 1.7 @@ -141,15 +148,14 @@ (cd*s2dsij) + (c2*s2sij2) + (ct*s2tij)) + - (0.5 * b3nl * c1 * sig3nl) + + (0.5 * b3 * c1 * sig3nl) + (0.5 * bk2 * c1 * Pd1k2)) -#Need to change these names and override old benchmark text file -np.savetxt("../pt_bm_test.txt", +np.savetxt("../pt_bm_z0.txt", np.transpose([ks, pgg[0], pgm[0], pgi[0], pii[0], pii_bb[0], pim[0]]), header='[0]-k [1]-GG [2]-GM [3]-GI [4]-II [5]-II_BB [6]-IM') -np.savetxt("../pt_bm_test2.txt", +np.savetxt("../pt_bm_z1.txt", np.transpose([ks, pgg[1], pgm[1], pgi[1], pii[1], pii_bb[1], pim[1]]), header='[0]-k [1]-GG [2]-GM [3]-GI [4]-II [5]-II_BB [6]-IM') From a2886677adb89b88d12f57eab9f844d94b0e3120 Mon Sep 17 00:00:00 2001 From: Vincent Schacknies Date: Fri, 31 Jan 2025 12:21:27 -0500 Subject: [PATCH 17/18] add ct in test_ept_decon --- pyccl/nl_pt/ept.py | 15 +++------------ pyccl/tests/test_ept_power.py | 4 ++-- 2 files changed, 5 insertions(+), 14 deletions(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 4e1721a41..c1b159c57 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -184,17 +184,8 @@ def __init__(self, *, with_NC=False, with_IA=False, else: self.exp_cutoff = 1 - # Call FAST-PT - try: - import fastpt as fpt - except: - raise ImportError("Your attempted import of FAST-PT has failed. You either dont have fast-pt installed, or have the wrong version. Try running pip install fast-pt or conda install fast-pt, then try again") - - if( not hasattr(fpt, "IA_ta")): - raise ValueError(f"Your FAST-PT version lacks a required attribute. You may have the wrong fast-pt install. Try running pip install fast-pt or conda install fast-pt, then try again") - if (not hasattr(fpt, "IA_tij")): - pass - #raise ValueError(f"You are using an older version of FAST-PT, please run pip install fast-pt or conda install fast-pt, then try again") + import fastpt as fpt + n_pad = int(self.fastpt_par['pad_factor'] * len(self.k_s)) self.pt = fpt.FASTPT(self.k_s, to_do=to_do, low_extrap=self.fastpt_par['low_extrap'], @@ -456,7 +447,7 @@ def _get_pgi(self, trg, tri): (self._g4*cd)[:,None] * (gb2dsij) + (self._g4*c2)[:, None] * (gb2sij2) + (self._g4*ct)[:,None] * gb2tij) + - 0.5*bs[:,None]*((self._g4*c1)[:,None]*Pd1s2+ #other term we tried: s2sij + + 0.5*bs[:,None]*((self._g4*c1)[:,None]*s2sij + (self._g4*cd)[:,None] * (s2dsij) + (self._g4*c2)[:,None] * (s2sij2) + (self._g4*ct)[:, None] *s2tij) + diff --git a/pyccl/tests/test_ept_power.py b/pyccl/tests/test_ept_power.py index 42f2257ee..3ddab1566 100644 --- a/pyccl/tests/test_ept_power.py +++ b/pyccl/tests/test_ept_power.py @@ -148,7 +148,7 @@ def get_tr(tn): @pytest.mark.parametrize('kind', ['c2:c2', 'c2:cdelta', 'cdelta:cdelta']) def test_ept_deconstruction_bb(kind): - b_ia = ['c1', 'c2', 'cdelta', 'ck'] + b_ia = ['c1', 'c2', 'cdelta', 'ck', 'ct'] pk1 = PTC.get_pk2d_template(kind, return_ia_bb=True) def get_tr(tn): @@ -156,7 +156,7 @@ def get_tr(tn): bdict[tn] = 1.0 return ccl.nl_pt.PTIntrinsicAlignmentTracer( c1=bdict['c1'], c2=bdict['c2'], - cdelta=bdict['cdelta'], ck=bdict['ck']) + cdelta=bdict['cdelta'], ck=bdict['ck'], ct=bdict['ct']) tn1, tn2 = kind.split(':') t1 = get_tr(tn1) From eba1c8d826febe6df855f4675f9ee50054fe502c Mon Sep 17 00:00:00 2001 From: Vincent Schacknies Date: Fri, 31 Jan 2025 12:53:10 -0500 Subject: [PATCH 18/18] revert match-case back --- pyccl/nl_pt/ept.py | 46 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index c1b159c57..e36e7fa69 100644 --- a/pyccl/nl_pt/ept.py +++ b/pyccl/nl_pt/ept.py @@ -736,6 +736,50 @@ def get_pk2d_template(self, kind, *, extrap_order_lok=1, # Construct power spectrum array s4 = 0. + if pk_name == 'm:m': + pk = self.pk_b1 + elif pk_name == 'm:b2': + pk = 0.5*self._g4T*self.dd_bias[2] + elif pk_name == 'm:b3nl': + pk = 0.5*self._g4T*self.dd_bias[8] + elif pk_name == 'm:bs': + pk = 0.5*self._g4T*self.dd_bias[4] + elif pk_name == 'm:bk2': + pk = 0.5*self.pk_bk*(self.k_s**2) + elif pk_name == 'm:c2': + pk = self._g4T * (self.ia_mix[0]+self.ia_mix[1]) + elif pk_name == 'm:cdelta': + pk = self._g4T * (self.ia_ta[0]+self.ia_ta[1]) + elif pk_name == 'b2:b2': + if self.fastpt_par['sub_lowk']: + s4 = self.dd_bias[7] + pk = 0.25*self._g4T*(self.dd_bias[3] - 2*s4) + elif pk_name == 'b2:bs': + if self.fastpt_par['sub_lowk']: + s4 = self.dd_bias[7] + pk = 0.25*self._g4T*(self.dd_bias[5] - 4*s4/3) + elif pk_name == 'bs:bs': + if self.fastpt_par['sub_lowk']: + s4 = self.dd_bias[7] + pk = 0.25*self._g4T*(self.dd_bias[6] - 8*s4/9) + elif pk_name == 'c2:c2': + pk = self._g4T * self.ia_tt[0] + elif pk_name == 'c2:c2_bb': + pk = self._g4T * self.ia_tt[1] + elif pk_name == 'c2:cdelta': + pk = self._g4T * self.ia_mix[2] + elif pk_name == 'c2:cdelta_bb': + pk = self._g4T * self.ia_mix[3] + elif pk_name == 'cdelta:cdelta': + pk = self._g4T * self.ia_ta[2] + elif pk_name == 'cdelta:cdelta_bb': + pk = self._g4T * self.ia_ta[3] + elif pk_name == 'zero': + # If zero, store None and return + self._pk2d_temp[pk_name] = None + return None + ''' + Match-case requires python 3.10, functionally the same as if-elif-else though slightly faster match pk_name: case 'm:m': pk = self.pk_b1 @@ -797,7 +841,7 @@ def get_pk2d_template(self, kind, *, extrap_order_lok=1, # If zero, store None and return self._pk2d_temp[pk_name] = None return None - + ''' # Build interpolator pk2d = Pk2D(a_arr=self.a_s,