diff --git a/fastpt/IA/IA_EFT.py b/fastpt/IA/IA_EFT.py new file mode 100644 index 0000000..11fd969 --- /dev/null +++ b/fastpt/IA/IA_EFT.py @@ -0,0 +1,44 @@ +import numpy as np +import sys + +nu = -2 +p_mat_11 = np.array([[0, 0, 0, 0], [0, 0, 2, 0], [0, 0, 4, 0], [2, -2, 2, 0], + [1, -1, 1, 0], [1, -1, 3, 0], [2, -2, 0, 1]]) +p_mat_12 = np.array([[0, 0, 0, 0], [0, 0, 2, 0], [0, 0, 4, 0], [1, -1, 1, 0], + [1, -1, 3, 0], [-1, 1, 1, 0], [-1, 1, 3, 0] +]) +p_mat_13 = np.array([[0, 0, 0, 0], [0, 0, 2, 0], [1, -1, 1, 0], [-1, 1, 1, 0]]) +p_mat_22 = np.array([[0, 0, 0, 0], [0, 0, 2, 0], [0, 0, 4, 0]]) +p_mat_23 = np.array([[0, 0, 0, 0], [0, 0, 2, 0]]) +p_mat_24 = np.array([[2, 0, 0, 0], [2, 0, 2, 0], [2, 0, 4, 0], [1, 1, 1, 0], + [1, 1, 3, 0], [1, 1, 5, 0], [0, 2, 0, 0], [0, 2, 2, 0], + [0, 2, 4, 0]]) +p_mat_33 = np.array([[0, 0, 0, 0]]) +p_mat_34 = np.array([[2, 0, 0, 0], [2, 0, 2, 0], [1, 1, 1, 0,], [1, 1, 3, 0], + [0, 2, 0, 0], [0, 2, 2, 0]]) +p_mat_44 = np.array([[4, 0, 0, 0,], [4, 0, 2, 0], [4, 0, 4, 0], [3, 1, 1, 0], + [3, 1, 3, 0], [3, 1, 5, 0], [2, 2, 0, 0], [2, 2, 2, 0], + [2, 2, 4, 0], [2, 2, 6, 0], [1, 3, 1, 0], [1, 3, 3, 0], + [1, 3, 5, 0], [0, 4, 0, 0], [0, 4, 2, 0], [0, 4, 4, 0]]) +p_mat_55 = np.array([[4, 0, 0, 0], [4, 0, 2, 0], [4, 0, 4, 0], [2, 2, 0, 0], + [2, 2, 2, 0], [2, 2, 4, 0], [0, 4, 0, 0], [0, 4, 2, 0], + [0, 4, 4, 0]]) + +coef_11 = [1219/1470., 671/1029., 32/1715., 1/3., 62/35., 8/35., 1/6.] +coef_12 = [31/105., 94/147., 16/245., 3/10., 1/5., 3/10., 1/5.] +coef_13 = [17/21., 4/21., 1/2., 1/2.] +coef_22 = [1/5., 4/7., 8/35.] +coef_23 = [1/3., 2/3.] +coef_24 = np.sqrt(2/3) * np.array([1/5., 4/7., 8/35., 3*13/35., 37/45, + 4/63., 1/5., 4/7., 8/35.]) +coef_33 = 1. +coef_34 = np.sqrt(2/3) * np.array([1/3., 2/3., 3*3/5., 1/5., 1/3., 2/3.]) +coef_44 = [2/15., 8/21., 16/105., 52/35., 148/135., 16/189., 104/105., 152/63., 676/1155, + 8/693., 52/35., 148/135., 16/189., 2/15., 8/21., 16/105.] +coef_55 = [1/30., 1/42., -2/35., -1/15., -1/21., 4/35., 1/30., 1/42., -2/35.] + +def IA_EFT_mat(): + return p_mat_11, p_mat_12, p_mat_13, p_mat_22, p_mat_23, p_mat_24, p_mat_33, p_mat_34, p_mat_44, p_mat_55 + +def IA_EFT_coef(): + return coef_11, coef_12, coef_13, coef_22, coef_23, coef_24, coef_33, coef_34, coef_44, coef_55 \ No newline at end of file diff --git a/fastpt/core/FASTPT.py b/fastpt/core/FASTPT.py index 952ad07..3cbd933 100644 --- a/fastpt/core/FASTPT.py +++ b/fastpt/core/FASTPT.py @@ -40,7 +40,7 @@ import numpy as np from numpy import exp, log, cos, sin, pi from ..utils.fastpt_extr import p_window, c_window -from ..utils.matter_power_spt import P_13_reg, Y1_reg_NL, Y2_reg_NL +from ..utils.matter_power_spt import P_13_reg, Y1_reg_NL, Y2_reg_NL, J2_integral, J3_integral from ..utils.initialize_params import scalar_stuff, tensor_stuff from ..IA.IA_tt import IA_tt from ..IA.IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B @@ -49,6 +49,7 @@ from ..IA.IA_gb2 import IA_gb2_S2F2, IA_gb2_S2fe, IA_gb2_S2he from ..IA.IA_ct import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, IA_tij_F2G2reg from ..IA.IA_ctbias import IA_gb2_F2, IA_gb2_G2, IA_gb2_S2F2, IA_gb2_S2G2 +from ..IA.IA_EFT import IA_EFT_mat, IA_EFT_coef from ..utils.OV import OV from ..utils.kPol import kPol from ..rsd.RSD import RSDA, RSDB @@ -244,7 +245,8 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high 'IA_mix': False, 'OV': False, 'kPol': False, 'RSD': False, 'IRres': False, 'tij': False, 'gb2': False, - 'all': False, 'everything': False + 'all': False, 'everything': False, + 'EFT': False } if to_do: @@ -257,6 +259,8 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high elif entry in {'IA_all', 'IA'}: for key in ['IA_tt', 'IA_ta', 'IA_mix', 'gb2', 'tij']: self.todo_dict[key] = True + elif entry == 'EFT': + self.todo_dict['EFT'] = True elif entry == 'dd_bias': self.todo_dict['one_loop_dd'] = True self.todo_dict['dd_bias'] = True @@ -318,6 +322,21 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high if self.todo_dict['RSD']: self.X_RSDA self.X_RSDB + + if self.todo_dict['EFT']: + # TODO: I don't like that this is so big, maybe move these to the bottom and have one call here. + nu = -2 + self.EFT_matrices = IA_EFT_mat() + self.Jabl_I11 = scalar_stuff(self.EFT_matrices[0], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I12 = scalar_stuff(self.EFT_matrices[1], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I13 = scalar_stuff(self.EFT_matrices[2], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I22 = scalar_stuff(self.EFT_matrices[3], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I23 = scalar_stuff(self.EFT_matrices[4], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I24 = scalar_stuff(self.EFT_matrices[5], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I33 = scalar_stuff(self.EFT_matrices[6], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I34 = scalar_stuff(self.EFT_matrices[7], nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I44 = scalar_stuff(self.EFT_matrices[8], -1.6, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.Jabl_I55 = scalar_stuff(self.EFT_matrices[9], -1.6, self.N, self.m, self.eta_m, self.l, self.tau_l) @property def k_original(self): @@ -1013,7 +1032,76 @@ def IA_ta(self, P, P_window=None, C_window=None): P_0E0E = self.compute_term("P_0E0E", self.X_IA_0E0E, P=P, P_window=P_window, C_window=C_window) P_0B0B = self.compute_term("P_0B0B", self.X_IA_0B0B, P=P, P_window=P_window, C_window=C_window) return P_deltaE1, P_deltaE2, P_0E0E, P_0B0B - + + def eft_integrals(self, P, P_window=None, C_window=None, remove_lowk=False): + # Returns the I_nm, J_n integrals based on arXiv:2303.15565. See Eqs. (2.39), (2.40) and + # corresponding kernels in Eqs. (A.1)-(A.3). Power spectra can be obtained using (2.41) and (2.63). + + # Coefficients for the (22)-type integrals: + IA_coef = IA_EFT_coef() + + def get_Inm(P, coef, Jabl, nu, P_window=P_window, C_window=C_window): + # Function to compute the final (22)-integrals: + Ps, mat = self.J_k_scalar(P, Jabl, nu, P_window, C_window) + P_mat = np.multiply(coef, np.transpose(mat)) + Inm = np.sum(P_mat, 1) + return Inm, Ps + + # Compute the (22)-integrals + I11, Ps = get_Inm(P, IA_coef[0], self.Jabl_I11, -2, C_window) + I12, _ = get_Inm(P, IA_coef[1], self.Jabl_I12, -2, C_window) + I13, _ = get_Inm(P, IA_coef[2], self.Jabl_I13, -2, C_window) + I22, _ = get_Inm(P, IA_coef[3], self.Jabl_I22, -2, C_window) + I23, _ = get_Inm(P, IA_coef[4], self.Jabl_I23, -2, C_window) + I24, _ = get_Inm(P, IA_coef[5], self.Jabl_I24, -2, C_window) + I33, _ = get_Inm(P, IA_coef[6], self.Jabl_I33, -2, C_window) + I34, _ = get_Inm(P, IA_coef[7], self.Jabl_I34, -2, C_window) + I44, _ = get_Inm(P, IA_coef[8], self.Jabl_I44, -1.6, C_window) + I55, _ = get_Inm(P, IA_coef[9], self.Jabl_I55, -1.6, C_window) + + I24 /= self.k_extrap ** 2 + I34 /= self.k_extrap ** 2 + I44 /= self.k_extrap ** 4 + I55 /= self.k_extrap ** 4 + + # subtract the low-k limit of the integral: + if remove_lowk: + # FIXME: Change this to the analytical expressions for the low-k limit once available. + print( + "Warning: Removing low-k from predictions based on initial k-value at this point.") + I22 -= I22[0] + I23 -= I23[0] + I24 -= I24[0] + I33 -= I33[0] + I44 -= I44[0] + I55 -= I55[0] + + # Compute the (13)-integrals: + J1 = P_13_reg(self.k_extrap, Ps) / 2 + J2 = J2_integral(self.k_extrap, Ps) + J3 = J3_integral(self.k_extrap, Ps) + + _, I11 = self.EK.PK_original(I11) + _, I12 = self.EK.PK_original(I12) + _, I13 = self.EK.PK_original(I13) + _, I22 = self.EK.PK_original(I22) + _, I23 = self.EK.PK_original(I23) + _, I24 = self.EK.PK_original(I24) + _, I33 = self.EK.PK_original(I33) + _, I34 = self.EK.PK_original(I34) + _, I44 = self.EK.PK_original(I44) + _, I55 = self.EK.PK_original(I55) + _, J1 = self.EK.PK_original(J1) + _, J2 = self.EK.PK_original(J2) + _, J3 = self.EK.PK_original(J3) + + I14 = ((28 * I12 - I22 + I23) / 2 / np.sqrt(6) - 5 * I24 + 5 * I34) / 7 + I66 = I22 / 9 - np.sqrt(6) / 9 * I24 + I44 / 6 + I67 = I22 / 36 + I23 / 12 - 5 * np.sqrt(6) / 72 * I24 - np.sqrt(6) / 24 * I34 + I44 / 6 + I77 = I22 / 144 + I23 / 24 + I33 / 16 - np.sqrt(6) * I24 / 36 - np.sqrt(6) * I34 / 12 + I44 / 6 + + return I11, I12, I13, I14, I22, I23, I24, I33, I34, I44, I55, I66, I67, I77, J1, J2, J3 + def _get_P_deltaE2(self, P): hash_key, P_hash = self._create_hash_key("P_deltaE2", None, P, None, None) result = self.cache.get("P_deltaE2", hash_key) diff --git a/fastpt/utils/matter_power_spt.py b/fastpt/utils/matter_power_spt.py index 6e603d1..005f4ab 100644 --- a/fastpt/utils/matter_power_spt.py +++ b/fastpt/utils/matter_power_spt.py @@ -154,3 +154,70 @@ def one_loop(k,P,P_window=None,C_window=None,n_pad=None): return P1,P22_reg,P13_reg +def J2_integral(k, P): + # calculates the J_2 integral in the EFT of IA + # via a discrete convolution integral + + N = k.size + n = np.arange(-N+1, N) + dL = log(k[1])-log(k[0]) + s = n*dL + + cut = 4 + high_s = s[s > cut] + low_s = s[s < -cut] + mid_high_s = s[(s <= cut) & (s > 0)] + mid_low_s = s[(s >= -cut) & (s < 0)] + + Z = lambda r: (15/r-55*r-55*r**3+15*r**5 + + (60-15/r**2-90*r**2+60*r**4-15*r**6)*log((r+1)/np.absolute(r-1))/2) + Z_low = lambda r: -128*r+384/7/r-128/21/r**3-128/231/r**5-128/1001/r**7 + Z_high = lambda r: -128*r**3+384/7*r**5-128/21*r**7 + + f_mid_low = Z(exp(-mid_low_s)) + f_mid_high = Z(exp(-mid_high_s)) + f_high = Z_high(exp(-high_s)) + f_low = Z_low(exp(-low_s)) + + f = np.hstack((f_low, f_mid_low, -80, f_mid_high, f_high)) + + g = fftconvolve(P, f) * dL + g_k = g[N-1:2*N-1] + P_bar = 1/42*k**3/(2*pi)**2*P*g_k + + return P_bar + + +def J3_integral(k, P): + # calculates the J_3 integral in the EFT of IA + # via a discrete convolution integral + + N = k.size + n = np.arange(-N+1, N) + dL = log(k[1])-log(k[0]) + s = n*dL + + cut = 4 + high_s = s[s > cut] + low_s = s[s < -cut] + mid_high_s = s[(s <= cut) & (s > 0)] + mid_low_s = s[(s >= -cut) & (s < 0)] + + Z = lambda r: (15/r-10*r+164*r**3-150*r**5+45*r**7+ + (15-15/r**2+90*r**2-210*r**4+165*r**6- + 45*r**8)*log((r+1)/np.absolute(r-1))/2) + Z_low = lambda r: 256/7*r+256/7/r-256/33/r**3-256/273/r**5-256/1001/r**7 + Z_high = lambda r: 256*r**3-2304/7*r**5+3328/21*r**7 + + f_mid_low = Z(exp(-mid_low_s)) + f_mid_high = Z(exp(-mid_high_s)) + f_high = Z_high(exp(-high_s)) + f_low = Z_low(exp(-low_s)) + + f = np.hstack((f_low, f_mid_low, 64, f_mid_high, f_high)) + + g = fftconvolve(P, f) * dL + g_k = g[N-1:2*N-1] + P_bar = 1/168*k**3/(2*pi)**2*P*g_k + + return P_bar