Skip to content

EFT implementation #9

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
chrgeorgiou wants to merge 3 commits into
base: master
Choose a base branch
from
Open

EFT implementation #9

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1b759e0
Compute EFT (22) and (13)+(31) integrals.
chrgeorgiou Aug 20, 2024
3821a19
Merge branch 'master' into EFT_implementation
chrgeorgiou Mar 26, 2026
f3206ae
Added EFT implementation changes to v4.
chrgeorgiou Mar 26, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
44 changes: 44 additions & 0 deletions fastpt/IA/IA_EFT.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
94 changes: 91 additions & 3 deletions fastpt/core/FASTPT.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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)
Expand Down
67 changes: 67 additions & 0 deletions fastpt/utils/matter_power_spt.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Loading