diff --git a/benchmarks/data/codes/pt_bm.py b/benchmarks/data/codes/pt_bm.py index 6015a3e4d..de61c223b 100644 --- a/benchmarks/data/codes/pt_bm.py +++ b/benchmarks/data/codes/pt_bm.py @@ -39,6 +39,25 @@ 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) + + +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, :] @@ -57,7 +76,25 @@ 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, :] + +sig3nl = g4 * one_loop_dd_bias_b3nl[8][None, :] + b1 = 1.3 b2 = 1.5 bs = 1.7 @@ -66,6 +103,8 @@ c1 = 1.9 c2 = 2.1 cd = 2.3 +ct = 2.5 +ck2 = 0.1 pgg = (b1**2 * Pd1d1 + b1*b2 * Pd1d2 + @@ -81,21 +120,36 @@ 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 + - 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) + + ck2 * 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)) + + (0.5 * b3 * c1 * sig3nl) + + (0.5 * bk2 * c1 * Pd1k2)) np.savetxt("../pt_bm_z0.txt", np.transpose([ks, pgg[0], pgm[0], pgi[0], diff --git a/pyccl/nl_pt/ept.py b/pyccl/nl_pt/ept.py index 4b42b7e5a..e36e7fa69 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, @@ -176,8 +184,8 @@ def __init__(self, *, with_NC=False, with_IA=False, else: self.exp_cutoff = 1 - # Call FAST-PT 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'], @@ -189,8 +197,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 +269,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 +303,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 +402,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 +431,29 @@ 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) + + (self._g4*ct)[:,None] * gb2tij) + + 0.5*bs[:,None]*((self._g4*c1)[:,None]*s2sij + + (self._g4*cd)[:,None] * (s2dsij) + + (self._g4*c2)[:,None] * (s2sij2) + + (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 +513,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 +541,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 +568,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): @@ -682,6 +778,70 @@ def get_pk2d_template(self, kind, *, extrap_order_lok=1, # 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 + 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, diff --git a/pyccl/nl_pt/tracers.py b/pyccl/nl_pt/tracers.py index b2e5ae313..c7dfd1691 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, 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,8 +40,9 @@ 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 = ck = ct = None gz = cosmo.growth_factor(1./(1+z)) + 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 @@ -54,8 +55,12 @@ 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 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 + return c1, c1delta, c2, ck, ct class PTTracer(CCLAutoRepr): @@ -212,7 +217,7 @@ class PTIntrinsicAlignmentTracer(PTTracer): """ type = 'IA' - def __init__(self, c1, c2=None, cdelta=None): + def __init__(self, c1, c2=None, cdelta=None, ck=None, ct=None): self.biases = {} @@ -222,6 +227,10 @@ 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 ck + self.biases['ck'] = self._get_bias_function(ck) + #Initialize ct + self.biases['ct'] = self._get_bias_function(ct) @property def c1(self): @@ -240,3 +249,15 @@ def cdelta(self): """Internal overdensity bias function. """ return self.biases['cdelta'] + + @property + 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 53b0741aa..3ddab1566 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', 'ct'] 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'], ct=bdict['ct']) 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):