From 17680f53a63eed59267a554359563ef32e9a9f1d Mon Sep 17 00:00:00 2001 From: Carter Date: Fri, 17 Nov 2023 14:06:50 -0600 Subject: [PATCH 01/20] Added s2 and deltas terms --- fastpt/FASTPT.py | 33 +++++++++++++++++++++++++++++++++ fastpt/IA_tij.py | 40 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+) create mode 100644 fastpt/IA_tij.py diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index e9401bf..2563b9b 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -49,6 +49,7 @@ from .IA_tt import IA_tt from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B +from .IA_tij import IA_tij_feG2, IA_tij_heG2 from .OV import OV from .kPol import kPol from .RSD import RSDA, RSDB @@ -172,6 +173,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.OV_do = False self.kPol_do = False self.RSD_do = False + self.IA_tij_do = False for entry in to_do: # convert to_do list to instructions for FAST-PT initialization if entry == 'one_loop_dd': @@ -212,6 +214,12 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high elif entry == 'IRres': self.dd_do = True continue + elif entry == 'tij': + self.IA_ta_do = True + self.IA_tt_do = True + self.IA_mix_do = True + self.IA_tij_do = True + continue elif entry == 'all' or entry == 'everything': self.dd_do = True self.dd_bias_do = True @@ -222,6 +230,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.kPol_do = True self.RSD_do = True self.cleft = True + self.IA_tij_do = True continue else: raise ValueError('FAST-PT does not recognize "' + entry + '" in the to_do list.') @@ -275,6 +284,14 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.X_IA_deltaE1 = tensor_stuff(p_mat_deltaE1, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_0E0E = tensor_stuff(p_mat_0E0E, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_0B0B = tensor_stuff(p_mat_0B0B, self.N, self.m, self.eta_m, self.l, self.tau_l) + + if self.IA_tij_do: + IA_tij_feG2_tab = IA_tij_feG2() + IA_tij_heG2_tab = IA_tij_heG2() + p_mat_tij_feG2 = IA_tij_feG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_tij_heG2 = IA_tij_heG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + self.X_IA_tij_feG2 = tensor_stuff(p_mat_tij_feG2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_heG2 = tensor_stuff(p_mat_tij_heG2, self.N, self.m, self.eta_m, self.l, self.tau_l) if self.OV_do: # For OV, we can use two different values for @@ -592,6 +609,22 @@ def IA_ta(self, P, P_window=None, C_window=None): def IA_der(self, P, P_window=None, C_window=None): P_der = (self.k_original**2)*P return P_der + + def IA_tij(self,P,P_window=None, C_window=None): + P_feG2, A = self.J_k_tensor(P,self.X_IA_tij_feG2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_feG2 = self.EK.PK_original(P_feG2) + P_heG2, A = self.J_k_tensor(P,self.X_IA_tij_heG2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_heG2 = self.EK.PK_original(P_heG2) + P_A00E,A,B,C = self.IA_ta(P, P_window=P_window, C_window=C_window) + P_A0E2,D,E,F = self.IA_mix(P,P_window=P_window, C_window=C_window) + P_feG2sub = np.subtract(P_feG2,P_A00E) + P_heG2sub = np.subtract(P_heG2,P_A0E2) + + return 2*P_feG2sub, 2*P_heG2sub + + def OV(self, P, P_window=None, C_window=None): P, A = self.J_k_tensor(P, self.X_OV, P_window=P_window, C_window=C_window) diff --git a/fastpt/IA_tij.py b/fastpt/IA_tij.py new file mode 100644 index 0000000..619cd61 --- /dev/null +++ b/fastpt/IA_tij.py @@ -0,0 +1,40 @@ +from __future__ import division +import numpy as np +from .J_table import J_table +import sys +from time import time +from numpy import log, sqrt, exp, pi +from scipy.signal import fftconvolve as convolve + +def IA_tij_feG2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_feG2=np.array([[0,0,0,2,0,13/21],\ + [0,0,0,2,2,8/21],\ + [1,-1,0,2,1,1/2],\ + [-1,1,0,2,1,1/2]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_feG2.shape[0]): + x=J_table(l_mat_tij_feG2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_heG2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_heG2=np.array([[0,0,0,0,0,-9/70],\ + [0,0,2,0,0,-26/63],\ + [0,0,0,0,2,-15/49],\ + [0,0,2,0,2,-16/63],\ + [0,0,1,1,1,81/70],\ + [0,0,1,1,3,12/35],\ + [0,0,0,0,4,-16/245],\ + [1,-1,0,0,1,-3/10],\ + [1,-1,2,0,1,-1/3],\ + [1,-1,1,1,0,1/2],\ + [1,-1,1,1,2,1],\ + [1,-1,0,2,1,-1/3],\ + [1,-1,0,0,3,-1/5]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_heG2.shape[0]): + x=J_table(l_mat_tij_heG2[i]) + table=np.row_stack((table,x)) + return table[1:,:] From 81e8855f1cc1dc72183d9698f45bef2cb8c1d7ce Mon Sep 17 00:00:00 2001 From: Carter Date: Tue, 5 Dec 2023 15:16:48 -0600 Subject: [PATCH 02/20] Added most Tij Terms --- fastpt/FASTPT.py | 34 ++++++++++++++++++++++++------ fastpt/FASTPT_simple.py | 12 +++++------ fastpt/IA_tij.py | 46 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 80 insertions(+), 12 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 2563b9b..b2b11f1 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -44,12 +44,12 @@ from scipy.signal import fftconvolve import scipy.integrate as integrate from .fastpt_extr import p_window, c_window, pad_left, pad_right -from .matter_power_spt import P_13_reg, Y1_reg_NL, Y2_reg_NL +from .matter_power_spt import P_13_reg, Y1_reg_NL, Y2_reg_NL, P_22 from .initialize_params import scalar_stuff, tensor_stuff from .IA_tt import IA_tt from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B -from .IA_tij import IA_tij_feG2, IA_tij_heG2 +from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2 from .OV import OV from .kPol import kPol from .RSD import RSDA, RSDB @@ -192,6 +192,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.IA_tt_do = True self.IA_ta_do = True self.IA_mix_do = True + self.IA_tij_do = True continue elif entry == 'IA_tt': self.IA_tt_do = True @@ -288,10 +289,19 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high if self.IA_tij_do: IA_tij_feG2_tab = IA_tij_feG2() IA_tij_heG2_tab = IA_tij_heG2() + IA_tij_F2F2_tab = IA_tij_F2F2() + IA_tij_G2G2_tab = IA_tij_G2G2() + IA_tij_F2G2_tab = IA_tij_F2G2() p_mat_tij_feG2 = IA_tij_feG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_heG2 = IA_tij_heG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_tij_F2F2 = IA_tij_F2F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_tij_G2G2 = IA_tij_G2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_tij_F2G2 = IA_tij_F2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] self.X_IA_tij_feG2 = tensor_stuff(p_mat_tij_feG2, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_tij_heG2 = tensor_stuff(p_mat_tij_heG2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_heG2 = tensor_stuff(p_mat_tij_heG2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_F2F2 = tensor_stuff(p_mat_tij_F2F2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_G2G2 = tensor_stuff(p_mat_tij_G2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_F2G2 = tensor_stuff(p_mat_tij_F2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) if self.OV_do: # For OV, we can use two different values for @@ -617,12 +627,24 @@ def IA_tij(self,P,P_window=None, C_window=None): P_heG2, A = self.J_k_tensor(P,self.X_IA_tij_heG2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_heG2 = self.EK.PK_original(P_heG2) + P_F2F2, A = self.J_k_tensor(P,self.X_IA_tij_F2F2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_F2F2 = self.EK.PK_original(P_F2F2) + P_G2G2, A = self.J_k_tensor(P,self.X_IA_tij_G2G2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_G2G2 = self.EK.PK_original(P_G2G2) + P_F2G2, A = self.J_k_tensor(P,self.X_IA_tij_F2G2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_F2G2 = self.EK.PK_original(P_F2G2) P_A00E,A,B,C = self.IA_ta(P, P_window=P_window, C_window=C_window) P_A0E2,D,E,F = self.IA_mix(P,P_window=P_window, C_window=C_window) - P_feG2sub = np.subtract(P_feG2,P_A00E) - P_heG2sub = np.subtract(P_heG2,P_A0E2) + P_13 = P_13_reg(self.k_original, P) + P_tijtij = P_F2F2+P_G2G2-2*P_F2G2 + P_tijsij = P_G2G2-P_F2F2-(1/2)*P_13 + P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) + P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) - return 2*P_feG2sub, 2*P_heG2sub + return 2*P_feG2sub, 2*P_heG2sub, 2*P_tijtij, 2*P_tijsij diff --git a/fastpt/FASTPT_simple.py b/fastpt/FASTPT_simple.py index 651f254..c120c6e 100644 --- a/fastpt/FASTPT_simple.py +++ b/fastpt/FASTPT_simple.py @@ -257,13 +257,13 @@ def P22(self,P,P_window=None,C_window=None): def one_loop(self,P,P_window=None,C_window=None): - Ps,P22=self.P22(P,P_window,C_window) - P13=P_13_reg(self.k_old,Ps) - if (self.extrap): - _,P=self.EK.PK_original(P22+P13) - return P + Ps,P22=self.P22(P,P_window,C_window) + P13=P_13_reg(self.k_old,Ps) + if (self.extrap): + _,P=self.EK.PK_original(P22+P13) + return P - return P22+P13 + return P22+P13 def P_bias(self,P,P_window=None,C_window=None): # Quadraric bias Legendre components diff --git a/fastpt/IA_tij.py b/fastpt/IA_tij.py index 619cd61..04be4d7 100644 --- a/fastpt/IA_tij.py +++ b/fastpt/IA_tij.py @@ -38,3 +38,49 @@ def IA_tij_heG2(): x=J_table(l_mat_tij_heG2[i]) table=np.row_stack((table,x)) return table[1:,:] + +def IA_tij_F2F2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_F2F2=np.array([[0,0,0,0,0,1219/1470],\ + [0,0,0,0,2,671/1029],\ + [0,0,0,0,4,32/1715],\ + [2,-2,0,0,0,1/6],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,62/35],\ + [1,-1,0,0,4,8/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_F2F2.shape[0]): + x=J_table(l_mat_tij_F2F2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_G2G2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_G2G2=np.array([[0,0,0,0,0,851/1470],\ + [0,0,0,0,2,871/1029],\ + [0,0,0,0,4,128/1715],\ + [2,-2,0,0,0,1/6],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,54/35],\ + [1,-1,0,0,4,16/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_G2G2.shape[0]): + x=J_table(l_mat_tij_G2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_F2G2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_F2G2=np.array([[0,0,0,0,0,1003/1470],\ + [0,0,0,0,2,803/1029],\ + [0,0,0,0,4,64/1715],\ + [2,-2,0,0,0,1/6],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,58/35],\ + [1,-1,0,0,4,12/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_F2G2.shape[0]): + x=J_table(l_mat_tij_F2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + From b4fb9f651c28627b994e97b599cfbc08d03f38a6 Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 7 Dec 2023 12:24:57 -0600 Subject: [PATCH 03/20] Adding gb2 terms --- fastpt/FASTPT.py | 48 ++++++++++++++++++++++++++++++++++++- fastpt/IA_gb2.py | 62 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 109 insertions(+), 1 deletion(-) create mode 100644 fastpt/IA_gb2.py diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index b2b11f1..b6446f7 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -50,6 +50,7 @@ from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2 +from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he, IA_gb2_base from .OV import OV from .kPol import kPol from .RSD import RSDA, RSDB @@ -174,6 +175,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.kPol_do = False self.RSD_do = False self.IA_tij_do = False + self.IA_gb2_do = False for entry in to_do: # convert to_do list to instructions for FAST-PT initialization if entry == 'one_loop_dd': @@ -193,6 +195,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.IA_ta_do = True self.IA_mix_do = True self.IA_tij_do = True + self.IA_gb2_do = True continue elif entry == 'IA_tt': self.IA_tt_do = True @@ -221,6 +224,9 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.IA_mix_do = True self.IA_tij_do = True continue + elif entry == 'gb2': + self.IA_gb2_do = True + continue elif entry == 'all' or entry == 'everything': self.dd_do = True self.dd_bias_do = True @@ -232,6 +238,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.RSD_do = True self.cleft = True self.IA_tij_do = True + self.IA_gb2_do = True continue else: raise ValueError('FAST-PT does not recognize "' + entry + '" in the to_do list.') @@ -301,7 +308,24 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.X_IA_tij_heG2 = tensor_stuff(p_mat_tij_heG2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_F2F2 = tensor_stuff(p_mat_tij_F2F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_G2G2 = tensor_stuff(p_mat_tij_G2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_tij_F2G2 = tensor_stuff(p_mat_tij_F2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_F2G2 = tensor_stuff(p_mat_tij_F2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + + if self.IA_gb2_do: + IA_gb2_fe_tab = IA_gb2_fe() + IA_gb2_he_tab = IA_gb2_he() + IA_gb2_F2_tab = IA_gb2_F2() + IA_gb2_G2_tab = IA_gb2_G2() + IA_gb2_base_tab = IA_gb2_base() + p_mat_gb2_fe = IA_gb2_fe_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_he = IA_gb2_he_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_F2 = IA_gb2_F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_G2 = IA_gb2_G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_base = IA_gb2_base_tab[:, [0,1,5,6,7,8,9]] + self.X_IA_gb2_fe = tensor_stuff(p_mat_gb2_fe, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_he = tensor_stuff(p_mat_gb2_he, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_F2 = tensor_stuff(p_mat_gb2_F2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_G2 = tensor_stuff(p_mat_gb2_G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_base = tensor_stuff(p_mat_gb2_base, self.N, self.m, self.eta_m, self.l, self.tau_l) if self.OV_do: # For OV, we can use two different values for @@ -645,6 +669,28 @@ def IA_tij(self,P,P_window=None, C_window=None): P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) return 2*P_feG2sub, 2*P_heG2sub, 2*P_tijtij, 2*P_tijsij + + def IA_gb2(self,P,P_window=None, C_window=None): + P_base, A = self.J_k_tensor(P,self.X_IA_gb2_base, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_base = self.EK.PK_original(P_base) + P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_fe = self.EK.PK_original(P_fe) + P_he, A = self.J_k_tensor(P,self.X_IA_gb2_he, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_he = self.EK.PK_original(P_he) + P_F2, A = self.J_k_tensor(P,self.X_IA_gb2_F2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_F2 = self.EK.PK_original(P_F2) + P_G2, A = self.J_k_tensor(P,self.X_IA_gb2_G2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_G2 = self.EK.PK_original(P_G2) + P_gb2sij = P_base + P_gb2sij2 = P_he + P_gb2tij = P_G2-P_F2 + P_gb2dsij = P_fe + return P_gb2sij, P_gb2dsij, P_gb2sij2, P_gb2tij diff --git a/fastpt/IA_gb2.py b/fastpt/IA_gb2.py new file mode 100644 index 0000000..689a5e8 --- /dev/null +++ b/fastpt/IA_gb2.py @@ -0,0 +1,62 @@ +from __future__ import division +import numpy as np +from .J_table import J_table +import sys +from time import time +from numpy import log, sqrt, exp, pi +from scipy.signal import fftconvolve as convolve + +def IA_gb2_fe(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_fe=np.array([[0,0,0,2,0,1]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_fe.shape[0]): + x=J_table(l_mat_gb2_fe[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_base(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_base=np.array([[0,0,0,0,0,1]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_base.shape[0]): + x=J_table(l_mat_gb2_base[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_he(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_he=np.array([[0,0,0,0,0,-1/6],\ + [0,0,0,2,0,-1/3],\ + [0,0,0,0,2,-1/3],\ + [0,0,1,1,1,3/2],\ + [0,0,0,0,2,-1/3]],dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_he.shape[0]): + x=J_table(l_mat_gb2_he[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_F2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_F2=np.array([[0,0,0,0,0,17/21],\ + [0,0,0,0,2,4/21],\ + [1,-1,0,0,1,1/2],\ + [-1,1,0,0,1,1/2]],dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_F2.shape[0]): + x=J_table(l_mat_gb2_F2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_G2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_G2=np.array([[0,0,0,0,0,13/21],\ + [0,0,0,0,2,8/21],\ + [1,-1,0,0,1,1/2],\ + [-1,1,0,0,1,1/2]],dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_G2.shape[0]): + x=J_table(l_mat_gb2_G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] \ No newline at end of file From c9c79f120b8e9eca6692184528b3c350c0b4b8c2 Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 7 Dec 2023 12:59:59 -0600 Subject: [PATCH 04/20] Corrected gb2sij correlation --- fastpt/FASTPT.py | 11 +++-------- fastpt/IA_gb2.py | 14 +++----------- 2 files changed, 6 insertions(+), 19 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index b6446f7..680f86c 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -50,7 +50,7 @@ from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2 -from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he, IA_gb2_base +from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he from .OV import OV from .kPol import kPol from .RSD import RSDA, RSDB @@ -315,17 +315,15 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high IA_gb2_he_tab = IA_gb2_he() IA_gb2_F2_tab = IA_gb2_F2() IA_gb2_G2_tab = IA_gb2_G2() - IA_gb2_base_tab = IA_gb2_base() p_mat_gb2_fe = IA_gb2_fe_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_gb2_he = IA_gb2_he_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_gb2_F2 = IA_gb2_F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_gb2_G2 = IA_gb2_G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_base = IA_gb2_base_tab[:, [0,1,5,6,7,8,9]] self.X_IA_gb2_fe = tensor_stuff(p_mat_gb2_fe, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_he = tensor_stuff(p_mat_gb2_he, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_F2 = tensor_stuff(p_mat_gb2_F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_G2 = tensor_stuff(p_mat_gb2_G2, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_gb2_base = tensor_stuff(p_mat_gb2_base, self.N, self.m, self.eta_m, self.l, self.tau_l) + if self.OV_do: # For OV, we can use two different values for @@ -671,9 +669,6 @@ def IA_tij(self,P,P_window=None, C_window=None): return 2*P_feG2sub, 2*P_heG2sub, 2*P_tijtij, 2*P_tijsij def IA_gb2(self,P,P_window=None, C_window=None): - P_base, A = self.J_k_tensor(P,self.X_IA_gb2_base, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_base = self.EK.PK_original(P_base) P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) if (self.extrap): _, P_fe = self.EK.PK_original(P_fe) @@ -686,7 +681,7 @@ def IA_gb2(self,P,P_window=None, C_window=None): P_G2, A = self.J_k_tensor(P,self.X_IA_gb2_G2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_G2 = self.EK.PK_original(P_G2) - P_gb2sij = P_base + P_gb2sij = P_F2 P_gb2sij2 = P_he P_gb2tij = P_G2-P_F2 P_gb2dsij = P_fe diff --git a/fastpt/IA_gb2.py b/fastpt/IA_gb2.py index 689a5e8..fc0f83c 100644 --- a/fastpt/IA_gb2.py +++ b/fastpt/IA_gb2.py @@ -8,29 +8,21 @@ def IA_gb2_fe(): # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_gb2_fe=np.array([[0,0,0,2,0,1]], dtype=float) + l_mat_gb2_fe=np.array([[0,0,2,0,0,1]], dtype=float) table=np.zeros(10,dtype=float) for i in range(l_mat_gb2_fe.shape[0]): x=J_table(l_mat_gb2_fe[i]) table=np.row_stack((table,x)) return table[1:,:] -def IA_gb2_base(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_gb2_base=np.array([[0,0,0,0,0,1]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_gb2_base.shape[0]): - x=J_table(l_mat_gb2_base[i]) - table=np.row_stack((table,x)) - return table[1:,:] def IA_gb2_he(): # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient l_mat_gb2_he=np.array([[0,0,0,0,0,-1/6],\ - [0,0,0,2,0,-1/3],\ + [0,0,2,0,0,-1/3],\ [0,0,0,0,2,-1/3],\ [0,0,1,1,1,3/2],\ - [0,0,0,0,2,-1/3]],dtype=float) + [0,0,2,0,0,-1/3]],dtype=float) table=np.zeros(10,dtype=float) for i in range(l_mat_gb2_he.shape[0]): x=J_table(l_mat_gb2_he[i]) From e05a1da1e5efc2b656387b939ae2fa5ba34b69ed Mon Sep 17 00:00:00 2001 From: Carter Date: Thu, 7 Dec 2023 13:44:14 -0600 Subject: [PATCH 05/20] added s2 correlation terms, ill behaved s2|sij2 --- fastpt/FASTPT.py | 34 ++++++++++++++++++++++++++++- fastpt/IA_gb2.py | 56 +++++++++++++++++++++++++++++++++++++++++++++++- 2 files changed, 88 insertions(+), 2 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 680f86c..8ec83d3 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -51,6 +51,7 @@ from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2 from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he +from .IA_gb2 import IA_gb2_S2F2, IA_gb2_S2G2, IA_gb2_S2fe, IA_gb2_S2he from .OV import OV from .kPol import kPol from .RSD import RSDA, RSDB @@ -315,14 +316,27 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high IA_gb2_he_tab = IA_gb2_he() IA_gb2_F2_tab = IA_gb2_F2() IA_gb2_G2_tab = IA_gb2_G2() + IA_gb2_S2F2_tab = IA_gb2_S2F2() + IA_gb2_S2G2_tab = IA_gb2_S2G2() + IA_gb2_S2fe_tab = IA_gb2_S2fe() + IA_gb2_S2he_tab = IA_gb2_S2he() p_mat_gb2_fe = IA_gb2_fe_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_gb2_he = IA_gb2_he_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_gb2_F2 = IA_gb2_F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_gb2_G2 = IA_gb2_G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_S2F2 = IA_gb2_S2F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_S2G2 = IA_gb2_S2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_S2fe = IA_gb2_S2fe_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_S2he = IA_gb2_S2he_tab[:, [0, 1, 5, 6, 7, 8, 9]] self.X_IA_gb2_fe = tensor_stuff(p_mat_gb2_fe, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_he = tensor_stuff(p_mat_gb2_he, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_F2 = tensor_stuff(p_mat_gb2_F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_G2 = tensor_stuff(p_mat_gb2_G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_S2F2 = tensor_stuff(p_mat_gb2_S2F2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_S2G2 = tensor_stuff(p_mat_gb2_S2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_S2fe = tensor_stuff(p_mat_gb2_S2fe, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_gb2_S2he = tensor_stuff(p_mat_gb2_S2he, self.N, self.m, self.eta_m, self.l, self.tau_l) + if self.OV_do: @@ -686,7 +700,25 @@ def IA_gb2(self,P,P_window=None, C_window=None): P_gb2tij = P_G2-P_F2 P_gb2dsij = P_fe return P_gb2sij, P_gb2dsij, P_gb2sij2, P_gb2tij - + + def IA_s2(self, P, P_window=None, C_window=None): + P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_S2F2 = self.EK.PK_original(P_S2F2) + P_S2G2, A = self.J_k_tensor(P, self.X_IA_gb2_S2G2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_S2G2 = self.EK.PK_original(P_S2G2) + P_S2fe, A = self.J_k_tensor(P, self.X_IA_gb2_S2fe, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_S2fe = self.EK.PK_original(P_S2fe) + P_S2he, A = self.J_k_tensor(P, self.X_IA_gb2_S2he, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_S2he = self.EK.PK_original(P_S2he) + P_s2sij=P_S2F2 + P_s2dsij=P_S2fe + P_s2sij2=P_S2he + P_s2tij=P_S2G2-P_S2F2 + return P_s2sij, P_s2dsij, P_s2sij2, P_s2tij def OV(self, P, P_window=None, C_window=None): diff --git a/fastpt/IA_gb2.py b/fastpt/IA_gb2.py index fc0f83c..4101b39 100644 --- a/fastpt/IA_gb2.py +++ b/fastpt/IA_gb2.py @@ -51,4 +51,58 @@ def IA_gb2_G2(): for i in range(l_mat_gb2_G2.shape[0]): x=J_table(l_mat_gb2_G2[i]) table=np.row_stack((table,x)) - return table[1:,:] \ No newline at end of file + return table[1:,:] + +def IA_gb2_S2F2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_S2F2=np.array([[0,0,0,0,0,8/315],\ + [0,0,0,0,2,254/441],\ + [0,0,0,0,4,16/245],\ + [1,-1,0,0,1,2/15],\ + [1,-1,0,0,3,1/5],\ + [-1,1,0,0,1,2/15],\ + [-1,1,0,0,1,1/5]],dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_S2F2.shape[0]): + x=J_table(l_mat_gb2_S2F2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_S2G2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_S2G2=np.array([[0,0,0,0,0,16/315],\ + [0,0,0,0,2,214/441],\ + [0,0,0,0,4,32/245],\ + [1,-1,0,0,1,2/15],\ + [1,-1,0,0,3,1/5],\ + [-1,1,0,0,1,2/15],\ + [-1,1,0,0,1,1/5]],dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_S2G2.shape[0]): + x=J_table(l_mat_gb2_S2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_S2fe(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_S2fe=np.array([[0,0,2,0,2,2/3]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_S2fe.shape[0]): + x=J_table(l_mat_gb2_S2fe[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_gb2_S2he(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_gb2_S2he=np.array([[0,0,0,0,0,-2/45],\ + [0,0,2,0,0,-11/63],\ + [0,0,2,0,2,-2/9],\ + [0,0,0,2,2,-2/9],\ + [0,0,1,1,1,2/5],\ + [0,0,1,1,3,3/5],\ + [0,0,0,0,4,-4/35]],dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_gb2_S2he.shape[0]): + x=J_table(l_mat_gb2_S2he[i]) + table=np.row_stack((table,x)) + return table[1:,:] \ No newline at end of file From ff4a46d2913f326089fabd4a15c0aaafe3958c00 Mon Sep 17 00:00:00 2001 From: Carter Date: Fri, 8 Dec 2023 11:32:33 -0600 Subject: [PATCH 06/20] Fixed s2sij2 term --- fastpt/FASTPT.py | 1 + fastpt/IA_gb2.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 8ec83d3..3f1e12f 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -695,6 +695,7 @@ def IA_gb2(self,P,P_window=None, C_window=None): P_G2, A = self.J_k_tensor(P,self.X_IA_gb2_G2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_G2 = self.EK.PK_original(P_G2) + sig4 = np.trapz(self.k_original ** 3 * P ** 2, x=np.log(self.k_original)) / (2. * pi ** 2) P_gb2sij = P_F2 P_gb2sij2 = P_he P_gb2tij = P_G2-P_F2 diff --git a/fastpt/IA_gb2.py b/fastpt/IA_gb2.py index 4101b39..87c30e7 100644 --- a/fastpt/IA_gb2.py +++ b/fastpt/IA_gb2.py @@ -95,7 +95,7 @@ def IA_gb2_S2fe(): def IA_gb2_S2he(): # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient l_mat_gb2_S2he=np.array([[0,0,0,0,0,-2/45],\ - [0,0,2,0,0,-11/63],\ + [0,0,0,0,2,-11/63],\ [0,0,2,0,2,-2/9],\ [0,0,0,2,2,-2/9],\ [0,0,1,1,1,2/5],\ From 53da4ee6ebb78e50fe95fc5315ba437d54f2dc96 Mon Sep 17 00:00:00 2001 From: Carter Date: Fri, 15 Dec 2023 16:57:07 -0600 Subject: [PATCH 07/20] Added s2 terms, fixing tij terms --- fastpt/FASTPT.py | 9 ++++---- fastpt/IA_tij.py | 58 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 4 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 3f1e12f..38370a0 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -49,7 +49,7 @@ from .IA_tt import IA_tt from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B -from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2 +from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he from .IA_gb2 import IA_gb2_S2F2, IA_gb2_S2G2, IA_gb2_S2fe, IA_gb2_S2he from .OV import OV @@ -674,13 +674,14 @@ def IA_tij(self,P,P_window=None, C_window=None): _, P_F2G2 = self.EK.PK_original(P_F2G2) P_A00E,A,B,C = self.IA_ta(P, P_window=P_window, C_window=C_window) P_A0E2,D,E,F = self.IA_mix(P,P_window=P_window, C_window=C_window) - P_13 = P_13_reg(self.k_original, P) + P_13F = P_IA_13F(self.k_original, P) + P_13G = P_IA_13G(self.k_original,P) P_tijtij = P_F2F2+P_G2G2-2*P_F2G2 - P_tijsij = P_G2G2-P_F2F2-(1/2)*P_13 + P_tijsij = P_G2G2+P_13G-P_F2F2-P_13F P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) - return 2*P_feG2sub, 2*P_heG2sub, 2*P_tijtij, 2*P_tijsij + return 2*P_feG2sub, 2*P_heG2sub, 2*P_13F, 2*P_13G def IA_gb2(self,P,P_window=None, C_window=None): P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) diff --git a/fastpt/IA_tij.py b/fastpt/IA_tij.py index 04be4d7..90f459f 100644 --- a/fastpt/IA_tij.py +++ b/fastpt/IA_tij.py @@ -84,3 +84,61 @@ def IA_tij_F2G2(): table=np.row_stack((table,x)) return table[1:,:] +def P_IA_13G(k,P): + N=k.size + n = np.arange(-N+1,N ) + dL=log(k[1])-log(k[0]) + s=n*dL + cut=7 + 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)] + + # For Zbar + Z1=lambda r : -82. + 12/r**2 + 4*r**2 -6*r**4 +log((r+1.)/np.absolute(r-1.))*(-6/r**3+15/r-9*r-3*r**3+3*r**5) + Z1_high=lambda r : -108 + (72/5)*r**2+ (12/7)*r**4-(352/105)*r**6-(164/105)*r**8 + (58/35)*r**10-(4/21)*r**12-(2/3)*r**14 + Z1_low=lambda r: -316/5+4/(35*r**8)+26/(21*r**6)+608/(105*r**4)-408/(35*r**2)+4/(3*r**12)-34/(21*r**10) + + f_mid_low=Z1(exp(-mid_low_s)) + f_mid_high=Z1(exp(-mid_high_s)) + f_high = Z1_high(exp(-high_s)) + f_low = Z1_low(exp(-low_s)) + + f=np.hstack((f_low,f_mid_low,-72.,f_mid_high,f_high)) + # print(f) + + g= convolve(P, f) * dL + g_k=g[N-1:2*N-1] + deltaE2= k**3/(336.*pi**2) * P*g_k + return deltaE2 + +def P_IA_13F(k,P): + N=k.size + n= np.arange(-N+1,N ) + dL=log(k[1])-log(k[0]) + s=n*dL + + cut=7 + 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 : (12./r**2 -158. + 100.*r**2-42.*r**4 \ + + 3./r**3*(r**2-1.)**3*(7*r**2+2.)*log((r+1.)/np.absolute(r-1.)) ) *r + Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8) *r + Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8) *r + + 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= convolve(P, f) * dL + g_k=g[N-1:2*N-1] + P_bar= 1./252.* k**3/(2*pi)**2*P*g_k + + return P_bar From 811e1d46cbf459582d753aea8e7a0edac19521a6 Mon Sep 17 00:00:00 2001 From: Carter Date: Fri, 15 Dec 2023 17:04:30 -0600 Subject: [PATCH 08/20] Multiplied gb2 and s2 terms by 2 --- fastpt/FASTPT.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 38370a0..776c7d7 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -701,7 +701,7 @@ def IA_gb2(self,P,P_window=None, C_window=None): P_gb2sij2 = P_he P_gb2tij = P_G2-P_F2 P_gb2dsij = P_fe - return P_gb2sij, P_gb2dsij, P_gb2sij2, P_gb2tij + return 2*P_gb2sij, 2*P_gb2dsij, 2*P_gb2sij2, 2*P_gb2tij def IA_s2(self, P, P_window=None, C_window=None): P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) @@ -720,7 +720,7 @@ def IA_s2(self, P, P_window=None, C_window=None): P_s2dsij=P_S2fe P_s2sij2=P_S2he P_s2tij=P_S2G2-P_S2F2 - return P_s2sij, P_s2dsij, P_s2sij2, P_s2tij + return 2*P_s2sij, 2*P_s2dsij, 2*P_s2sij2, 2*P_s2tij def OV(self, P, P_window=None, C_window=None): From 32dd2de224b5a01f806c61d5645a6b4017d21f9b Mon Sep 17 00:00:00 2001 From: CarterDW Date: Fri, 1 Mar 2024 16:00:08 -0500 Subject: [PATCH 09/20] Fixing TijSij Cross Correlation --- fastpt/FASTPT.py | 54 ++++++++++++++++++++++++++++++++++++--- fastpt/IA_tij.py | 66 +++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 107 insertions(+), 13 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 776c7d7..bf7b1b7 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -49,9 +49,10 @@ from .IA_tt import IA_tt from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B -from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F +from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, P_22F_reg, P_22G_reg from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he from .IA_gb2 import IA_gb2_S2F2, IA_gb2_S2G2, IA_gb2_S2fe, IA_gb2_S2he +from .J_k import J_k from .OV import OV from .kPol import kPol from .RSD import RSDA, RSDB @@ -424,6 +425,44 @@ def one_loop_dd(self, P, P_window=None, C_window=None): _, P_1loop = self.EK.PK_original(P_1loop) return P_1loop, Ps + + def tester_function(self,P,P_window=None,C_window=None): + one_loop_coef = np.array( + [2 * 1219 / 1470., 2 * 671 / 1029., 2 * 32 / 1715., 2 * 1 / 3., 2 * 62 / 35., 2 * 8 / 35., 1 / 3.]) + + # get the roundtrip Fourier power spectrum, i.e. P=IFFT[FFT[P]] + # get the matrix for each J_k component + nu = -2 + Ps, mat = self.J_k_scalar(P, self.X_spt, nu, P_window=P_window, C_window=C_window) + + P22_mat = np.multiply(one_loop_coef, np.transpose(mat)) + P22 = np.sum(P22_mat, 1) + P13 = P_13_reg(self.k_old, Ps) + #All new terms + P_13F = P_IA_13F(self.k_original, P) + P_13G = P_IA_13G(self.k_original,P) + #P_F2F2, A = self.J_k_tensor(P,self.X_IA_tij_F2F2, P_window=P_window, C_window=C_window) + #if (self.extrap): + # _, P_F2F2 = self.EK.PK_original(P_F2F2) + #P_G2G2, A = self.J_k_tensor(P,self.X_IA_tij_G2G2, P_window=P_window, C_window=C_window) + #if (self.extrap): + # _, P_G2G2 = self.EK.PK_original(P_G2G2) + #param_matrix=np.array([[2,-2,0,1]]) + #Power, mat=J_k(self.k_original,P,param_matrix,P_window=P_window, C_window=C_window, n_pad=n_pad) + #regF = 1/3*mat[0,:] + #regG = 1/9*mat[0,:] + #P_22F=P_F2F2+regF + #P_22G=P_G2G2+regG + P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) + P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) + if (self.extrap): + _, P22 = self.EK.PK_original(P22) + _, P13 = self.EK.PK_original(P13) + return P22, P13, P_22F, P_13F, P_22G, P_13G + + + + def one_loop_dd_bias(self, P, P_window=None, C_window=None): nu = -2 @@ -675,13 +714,20 @@ def IA_tij(self,P,P_window=None, C_window=None): P_A00E,A,B,C = self.IA_ta(P, P_window=P_window, C_window=C_window) P_A0E2,D,E,F = self.IA_mix(P,P_window=P_window, C_window=C_window) P_13F = P_IA_13F(self.k_original, P) - P_13G = P_IA_13G(self.k_original,P) + P_13G = P_IA_13G(self.k_original,P,) + P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) + P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) + + #param_matrix=np.array([[2,-2,0,1]]) + #Power, mat=J_k(self.k_original,P,param_matrix,P_window=P_window, C_window=C_window, n_pad=n_pad) + #regF = 1/3*mat[0,:] + #regG = 1/9*mat[0,:] P_tijtij = P_F2F2+P_G2G2-2*P_F2G2 - P_tijsij = P_G2G2+P_13G-P_F2F2-P_13F + P_tijsij = P_22G-P_22F+P_13G-P_13F P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) - return 2*P_feG2sub, 2*P_heG2sub, 2*P_13F, 2*P_13G + return P_tijsij,P_feG2sub,P_heG2sub,P_tijtij def IA_gb2(self,P,P_window=None, C_window=None): P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) diff --git a/fastpt/IA_tij.py b/fastpt/IA_tij.py index 90f459f..1585a4d 100644 --- a/fastpt/IA_tij.py +++ b/fastpt/IA_tij.py @@ -1,6 +1,7 @@ from __future__ import division import numpy as np -from .J_table import J_table +from .J_table import J_table +from .J_k import J_k import sys from time import time from numpy import log, sqrt, exp, pi @@ -96,23 +97,70 @@ def P_IA_13G(k,P): mid_low_s=s[ (s >= -cut) & (s < 0)] # For Zbar - Z1=lambda r : -82. + 12/r**2 + 4*r**2 -6*r**4 +log((r+1.)/np.absolute(r-1.))*(-6/r**3+15/r-9*r-3*r**3+3*r**5) - Z1_high=lambda r : -108 + (72/5)*r**2+ (12/7)*r**4-(352/105)*r**6-(164/105)*r**8 + (58/35)*r**10-(4/21)*r**12-(2/3)*r**14 - Z1_low=lambda r: -316/5+4/(35*r**8)+26/(21*r**6)+608/(105*r**4)-408/(35*r**2)+4/(3*r**12)-34/(21*r**10) + Z1=lambda r : (12/r**2-26+4*r**2-6*r**4+(3/r**3)*(r**2-1)**3*(r**2+2)*log((r+1.)/np.absolute(r-1.))) *r + Z1_low=lambda r : (-224/5+1248/(35*r**2)-608/(105*r**4)-26/(21*r**6)-4/(35*r**8)+34/(21*r**10) -4/(3*r**12)) *r + Z1_high=lambda r: ((-32*r**2)/5-(96*r**4)/7+(352*r**6)/105+(164*r**8)/105-(58*r**10)/35+(4*r**12)/21+(2*r**14)/3) *r f_mid_low=Z1(exp(-mid_low_s)) f_mid_high=Z1(exp(-mid_high_s)) f_high = Z1_high(exp(-high_s)) f_low = Z1_low(exp(-low_s)) - f=np.hstack((f_low,f_mid_low,-72.,f_mid_high,f_high)) + f=np.hstack((f_low,f_mid_low,-16.,f_mid_high,f_high)) # print(f) g= convolve(P, f) * dL g_k=g[N-1:2*N-1] - deltaE2= k**3/(336.*pi**2) * P*g_k + deltaE2= 1/84 * k**3/(2*pi)**2 * P*g_k return deltaE2 +def P_22F_reg(k,P, P_window, C_window, n_pad): + # P_22 Legendre components + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. + + param_matrix=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] ]) + + + Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) + A=1219/1470.*mat[0,:] + B=671/1029.*mat[1,:] + C=32/1715.*mat[2,:] + D=1/3.*mat[3,:] + E=62/35.*mat[4,:] + F=8/35.*mat[5,:] + reg=1/3.*mat[6,:] + + return 2*(A+B+C+D+E+F)+ reg + +def P_22G_reg(k,P, P_window, C_window, n_pad): + # P_22 Legendre components + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. + + param_matrix=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], [2,-2,0,0]]) + + + Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) + A=1003/1470.*mat[0,:] + B=803/1029.*mat[1,:] + C=64/1715.*mat[2,:] + D=1/3.*mat[3,:] + E=58/35.*mat[4,:] + F=12/35.*mat[5,:] + reg=1/3*mat[6,:] + #G=1/6.*mat[7,:] + + return 2*(A+B+C+D+E+F)+ reg + def P_IA_13F(k,P): N=k.size n= np.arange(-N+1,N ) @@ -125,10 +173,10 @@ def P_IA_13F(k,P): mid_high_s=s[ (s <= cut) & (s > 0)] mid_low_s=s[ (s >= -cut) & (s < 0)] - Z=lambda r : (12./r**2 -158. + 100.*r**2-42.*r**4 \ + Z=lambda r : (12./r**2 +10. + 100.*r**2-42.*r**4 \ + 3./r**3*(r**2-1.)**3*(7*r**2+2.)*log((r+1.)/np.absolute(r-1.)) ) *r - Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8) *r - Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8) *r + Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8-50./21./r**10-4./3./r**12) *r + Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8+74./35.*r**10-20./3.*r**12+14./3.*r**14) *r f_mid_low=Z(exp(-mid_low_s)) f_mid_high=Z(exp(-mid_high_s)) From c9bd07afd9724bc12d47e6f3d0f5f1fb03a524f8 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Thu, 14 Mar 2024 15:23:07 -0400 Subject: [PATCH 10/20] Still fixing tijsij --- fastpt/FASTPT.py | 21 +++++-- fastpt/IA_tij.py | 143 +++++++++++++++++++++++++++-------------------- 2 files changed, 97 insertions(+), 67 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index bf7b1b7..2c19983 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -49,7 +49,7 @@ from .IA_tt import IA_tt from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B -from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, P_22F_reg, P_22G_reg +from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, P_22F_reg, P_22G_reg, IA_tij_F2G2reg from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he from .IA_gb2 import IA_gb2_S2F2, IA_gb2_S2G2, IA_gb2_S2fe, IA_gb2_S2he from .J_k import J_k @@ -301,16 +301,19 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high IA_tij_F2F2_tab = IA_tij_F2F2() IA_tij_G2G2_tab = IA_tij_G2G2() IA_tij_F2G2_tab = IA_tij_F2G2() + IA_tij_F2G2reg_tab =IA_tij_F2G2reg() p_mat_tij_feG2 = IA_tij_feG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_heG2 = IA_tij_heG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_F2F2 = IA_tij_F2F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_G2G2 = IA_tij_G2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_F2G2 = IA_tij_F2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_tij_F2G2reg_tab = IA_tij_F2G2reg_tab[:, [0, 1, 5, 6, 7, 8, 9]] self.X_IA_tij_feG2 = tensor_stuff(p_mat_tij_feG2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_heG2 = tensor_stuff(p_mat_tij_heG2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_F2F2 = tensor_stuff(p_mat_tij_F2F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_G2G2 = tensor_stuff(p_mat_tij_G2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_tij_F2G2 = tensor_stuff(p_mat_tij_F2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_F2G2 = tensor_stuff(p_mat_tij_F2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_IA_tij_F2G2reg = tensor_stuff(p_mat_tij_F2G2reg_tab, self.N, self.m, self.eta_m, self.l, self.tau_l) if self.IA_gb2_do: IA_gb2_fe_tab = IA_gb2_fe() @@ -455,9 +458,19 @@ def tester_function(self,P,P_window=None,C_window=None): #P_22G=P_G2G2+regG P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) + #P_22Gnonreg, A = self.J_k_tensor(P,self.X_IA_tij_F2G2reg, P_window=P_window, C_window=C_window) + #param_matrix=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]]) + + + #Power, mat=J_k(self.k_original,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=self.n_pad) + #reg = 1/3*mat[6,:] + if (self.extrap): _, P22 = self.EK.PK_original(P22) _, P13 = self.EK.PK_original(P13) + #_, P_22Gnonreg = self.EK.PK_original(P_22Gnonreg) + #P_22G=P_22Gnonreg+reg return P22, P13, P_22F, P_13F, P_22G, P_13G @@ -718,10 +731,6 @@ def IA_tij(self,P,P_window=None, C_window=None): P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) - #param_matrix=np.array([[2,-2,0,1]]) - #Power, mat=J_k(self.k_original,P,param_matrix,P_window=P_window, C_window=C_window, n_pad=n_pad) - #regF = 1/3*mat[0,:] - #regG = 1/9*mat[0,:] P_tijtij = P_F2F2+P_G2G2-2*P_F2G2 P_tijsij = P_22G-P_22F+P_13G-P_13F P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) diff --git a/fastpt/IA_tij.py b/fastpt/IA_tij.py index 1585a4d..1de9e45 100644 --- a/fastpt/IA_tij.py +++ b/fastpt/IA_tij.py @@ -116,77 +116,98 @@ def P_IA_13G(k,P): def P_22F_reg(k,P, P_window, C_window, n_pad): # P_22 Legendre components - # We calculate a regularized version of P_22 - # by omitting the J_{2,-2,0} term so that the - # integral converges. In the final power spectrum - # we add the asymptotic portions of P_22 and P_13 so - # that we get a convergent integral. See section of XXX. + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. - param_matrix=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] ]) + param_matrix=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] ]) - Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) - A=1219/1470.*mat[0,:] - B=671/1029.*mat[1,:] - C=32/1715.*mat[2,:] - D=1/3.*mat[3,:] - E=62/35.*mat[4,:] - F=8/35.*mat[5,:] - reg=1/3.*mat[6,:] + Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) + A=1219/1470.*mat[0,:] + B=671/1029.*mat[1,:] + C=32/1715.*mat[2,:] + D=1/3.*mat[3,:] + E=62/35.*mat[4,:] + F=8/35.*mat[5,:] + reg=1/3.*mat[6,:] - return 2*(A+B+C+D+E+F)+ reg + return 2*(A+B+C+D+E+F)+ reg def P_22G_reg(k,P, P_window, C_window, n_pad): # P_22 Legendre components - # We calculate a regularized version of P_22 - # by omitting the J_{2,-2,0} term so that the - # integral converges. In the final power spectrum - # we add the asymptotic portions of P_22 and P_13 so - # that we get a convergent integral. See section of XXX. + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. - param_matrix=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], [2,-2,0,0]]) + param_matrix=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]]) - Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) - A=1003/1470.*mat[0,:] - B=803/1029.*mat[1,:] - C=64/1715.*mat[2,:] - D=1/3.*mat[3,:] - E=58/35.*mat[4,:] - F=12/35.*mat[5,:] - reg=1/3*mat[6,:] - #G=1/6.*mat[7,:] + Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) + A=1003/1470.*mat[0,:] + B=803/1029.*mat[1,:] + C=64/1715.*mat[2,:] + D=1/3.*mat[3,:] + E=58/35.*mat[4,:] + F=12/35.*mat[5,:] + reg=1/3.*mat[6,:] + + + return 2*(A+B+C+D+E+F)+ reg + +def IA_tij_F2G2reg(): + # P_22 Legendre components + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. + + l_mat_tij_F2G2=np.array([[0,0,0,0,0,1003/1470],\ + [0,0,0,0,2,803/1029],\ + [0,0,0,0,4,64/1715],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,58/35],\ + [1,-1,0,0,4,12/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_F2G2.shape[0]): + x=J_table(l_mat_tij_F2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] - return 2*(A+B+C+D+E+F)+ reg def P_IA_13F(k,P): - N=k.size - n= np.arange(-N+1,N ) - dL=log(k[1])-log(k[0]) - s=n*dL - - cut=7 - 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 : (12./r**2 +10. + 100.*r**2-42.*r**4 \ - + 3./r**3*(r**2-1.)**3*(7*r**2+2.)*log((r+1.)/np.absolute(r-1.)) ) *r - Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8-50./21./r**10-4./3./r**12) *r - Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8+74./35.*r**10-20./3.*r**12+14./3.*r**14) *r - - 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= convolve(P, f) * dL - g_k=g[N-1:2*N-1] - P_bar= 1./252.* k**3/(2*pi)**2*P*g_k - - return P_bar + N=k.size + n= np.arange(-N+1,N ) + dL=log(k[1])-log(k[0]) + s=n*dL + + cut=7 + 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 : (12./r**2 +10. + 100.*r**2-42.*r**4 \ + + 3./r**3*(r**2-1.)**3*(7*r**2+2.)*log((r+1.)/np.absolute(r-1.)) ) *r + Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8-50./21./r**10-4./3./r**12) *r + Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8+74./35.*r**10-20./3.*r**12+14./3.*r**14) *r + + 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= convolve(P, f) * dL + g_k=g[N-1:2*N-1] + P_bar= 1./252.* k**3/(2*pi)**2*P*g_k + + return P_bar From ad6861a6bd0a2f93d6ecd5b2fd4640f04fbff5f9 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 30 Apr 2024 14:30:38 -0400 Subject: [PATCH 11/20] Updated factors of 2 to reflect appropriate values --- fastpt/FASTPT.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 2c19983..632efb7 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -736,7 +736,7 @@ def IA_tij(self,P,P_window=None, C_window=None): P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) - return P_tijsij,P_feG2sub,P_heG2sub,P_tijtij + return 2*P_tijsij,2*P_feG2sub,2*P_heG2sub,2*P_tijtij def IA_gb2(self,P,P_window=None, C_window=None): P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) From dc0010ca059fbbdb196b5468fa25004b772da87a Mon Sep 17 00:00:00 2001 From: CarterDW Date: Mon, 6 May 2024 13:40:57 -0400 Subject: [PATCH 12/20] Changed structure of tijsij calculation --- fastpt/FASTPT.py | 32 +++++++++++++++++++++++++++++--- 1 file changed, 29 insertions(+), 3 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 632efb7..6a3cd38 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -221,6 +221,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.dd_do = True continue elif entry == 'tij': + self.IA_dd_do = True self.IA_ta_do = True self.IA_tt_do = True self.IA_mix_do = True @@ -257,6 +258,7 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.X_spt = scalar_stuff(p_mat, nu, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_lpt = scalar_stuff(p_mat_lpt, nu, self.N, self.m, self.eta_m, self.l, self.tau_l) + self.X_sptG = scalar_stuff(p_mat, nu, self.N, self.m, self.eta_m, self.l, self.tau_l) if self.cleft: nu = -2 @@ -391,6 +393,8 @@ def one_loop_dd(self, P, P_window=None, C_window=None): P13 = P_13_reg(self.k_old, Ps) P_1loop = P22 + P13 + + if (self.dd_bias_do): # if dd_bias is in to_do, this function acts like one_loop_dd_bias @@ -441,6 +445,7 @@ def tester_function(self,P,P_window=None,C_window=None): P22_mat = np.multiply(one_loop_coef, np.transpose(mat)) P22 = np.sum(P22_mat, 1) P13 = P_13_reg(self.k_old, Ps) + #All new terms P_13F = P_IA_13F(self.k_original, P) P_13G = P_IA_13G(self.k_original,P) @@ -458,6 +463,14 @@ def tester_function(self,P,P_window=None,C_window=None): #P_22G=P_G2G2+regG P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) + + #Test of different P_22G + one_loop_coefG= np.array( + [2*1003/1470, 2*803/1029, 2*64/1715, 2*1/3, 2*58/35, 2*12/35, 1/3]) + + PsG, matG = self.J_k_scalar(P, self.X_sptG, nu, P_window=P_window, C_window=C_window) + P22G_mat = np.multiply(one_loop_coefG, np.transpose(matG)) + P_22G2 = np.sum(P22G_mat, 1) #P_22Gnonreg, A = self.J_k_tensor(P,self.X_IA_tij_F2G2reg, P_window=P_window, C_window=C_window) #param_matrix=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]]) @@ -469,9 +482,10 @@ def tester_function(self,P,P_window=None,C_window=None): if (self.extrap): _, P22 = self.EK.PK_original(P22) _, P13 = self.EK.PK_original(P13) + _, P_22G2 = self.EK.PK_original(P_22G2) #_, P_22Gnonreg = self.EK.PK_original(P_22Gnonreg) #P_22G=P_22Gnonreg+reg - return P22, P13, P_22F, P_13F, P_22G, P_13G + return P22, P13, P_22F, P_13F, P_22G, P_13G, P_22G2 @@ -728,9 +742,21 @@ def IA_tij(self,P,P_window=None, C_window=None): P_A0E2,D,E,F = self.IA_mix(P,P_window=P_window, C_window=C_window) P_13F = P_IA_13F(self.k_original, P) P_13G = P_IA_13G(self.k_original,P,) - P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) - P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) + nu=-2 + Ps, mat = self.J_k_scalar(P, self.X_spt, nu, P_window=P_window, C_window=C_window) + one_loop_coef = np.array( + [2 * 1219 / 1470., 2 * 671 / 1029., 2 * 32 / 1715., 2 * 1 / 3., 2 * 62 / 35., 2 * 8 / 35., 1 / 3.]) + P22_mat = np.multiply(one_loop_coef, np.transpose(mat)) + P_22F = np.sum(P22_mat, 1) + one_loop_coefG= np.array( + [2*1003/1470, 2*803/1029, 2*64/1715, 2*1/3, 2*58/35, 2*12/35, 1/3]) + PsG, matG = self.J_k_scalar(P, self.X_sptG, nu, P_window=P_window, C_window=C_window) + P22G_mat = np.multiply(one_loop_coefG, np.transpose(matG)) + P_22G = np.sum(P22G_mat, 1) + if (self.extrap): + _, P_22F=self.EK.PK_original(P_22F) + _, P_22G=self.EK.PK_original(P_22G) P_tijtij = P_F2F2+P_G2G2-2*P_F2G2 P_tijsij = P_22G-P_22F+P_13G-P_13F P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) From 0ba548fc491143f8a46a60289effb71fdb57a17b Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 6 Aug 2024 11:20:06 -0400 Subject: [PATCH 13/20] Fixed S2F2 term --- fastpt/FASTPT.py | 7 +++++-- fastpt/IA_gb2.py | 35 ++++++++++++++++++++++++++++++++++- 2 files changed, 39 insertions(+), 3 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 6a3cd38..d0d5e76 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -50,7 +50,7 @@ from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, P_22F_reg, P_22G_reg, IA_tij_F2G2reg -from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he +from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he, P_IA_13S2F2 from .IA_gb2 import IA_gb2_S2F2, IA_gb2_S2G2, IA_gb2_S2fe, IA_gb2_S2he from .J_k import J_k from .OV import OV @@ -788,6 +788,9 @@ def IA_s2(self, P, P_window=None, C_window=None): P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_S2F2 = self.EK.PK_original(P_S2F2) + + P_13S2F2 = P_IA_13S2F2(self.k_original, P) + P_S2G2, A = self.J_k_tensor(P, self.X_IA_gb2_S2G2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_S2G2 = self.EK.PK_original(P_S2G2) @@ -797,7 +800,7 @@ def IA_s2(self, P, P_window=None, C_window=None): P_S2he, A = self.J_k_tensor(P, self.X_IA_gb2_S2he, P_window=P_window, C_window=C_window) if (self.extrap): _, P_S2he = self.EK.PK_original(P_S2he) - P_s2sij=P_S2F2 + P_s2sij=P_S2F2+2*P_13S2F2 P_s2dsij=P_S2fe P_s2sij2=P_S2he P_s2tij=P_S2G2-P_S2F2 diff --git a/fastpt/IA_gb2.py b/fastpt/IA_gb2.py index 87c30e7..18872a3 100644 --- a/fastpt/IA_gb2.py +++ b/fastpt/IA_gb2.py @@ -6,6 +6,39 @@ from numpy import log, sqrt, exp, pi from scipy.signal import fftconvolve as convolve +def P_IA_13S2F2(k,P): + + N=k.size + n= np.arange(-N+1,N ) + dL=log(k[1])-log(k[0]) + s=n*dL + cut=7 + 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)] + + + + Z1=lambda r : (((4.*r*(45.-165.*r**2+379.*r**4+45.*r**6)+45.*(-1.+r**2)**4*log((-1.+r)**2)-90.*(-1.+r**2)**4*log(np.absolute(1.+r)))/(2016.*r**3))-68./63*r**2)/2. + Z1_high=lambda r : ((-16*r**2)/21. + (16*r**4)/49. - (16*r**6)/441. - (16*r**8)/4851. - 16*r**10/21021.- 16*r**12/63063.)/2. + Z1_low=lambda r: (-16/21.+ 16/(49.*r**2) - 16/(441.*r**4) - 16/(4851.*r**6) - 16/(21021.*r**8) - 16/(63063.*r**10)- 16/(153153.*r**12) )/2. + + + f_mid_low=Z1(exp(-mid_low_s))*exp(-mid_low_s) + f_mid_high=Z1(exp(-mid_high_s))*exp(-mid_high_s) + f_high = Z1_high(exp(-high_s))*exp(-high_s) + f_low = Z1_low(exp(-low_s))*exp(-low_s) + + f=np.hstack((f_low,f_mid_low,-0.2381002916036672,f_mid_high,f_high)) + + + g= convolve(P, f) * dL + g_k=g[N-1:2*N-1] + S2F2= k**3/(2.*pi**2) * P*g_k + return S2F2 + + def IA_gb2_fe(): # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient l_mat_gb2_fe=np.array([[0,0,2,0,0,1]], dtype=float) @@ -61,7 +94,7 @@ def IA_gb2_S2F2(): [1,-1,0,0,1,2/15],\ [1,-1,0,0,3,1/5],\ [-1,1,0,0,1,2/15],\ - [-1,1,0,0,1,1/5]],dtype=float) + [-1,1,0,0,3,1/5]],dtype=float) table=np.zeros(10,dtype=float) for i in range(l_mat_gb2_S2F2.shape[0]): x=J_table(l_mat_gb2_S2F2[i]) From 5cd7a93753374e5a89fc758b213e4e9866e31dce Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 6 Aug 2024 12:59:19 -0400 Subject: [PATCH 14/20] Moved galaxy bias x tij terms into a seperate function --- fastpt/FASTPT.py | 35 +++++++++++++++++++++++++---------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index d0d5e76..1ec0bc7 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -764,6 +764,29 @@ def IA_tij(self,P,P_window=None, C_window=None): return 2*P_tijsij,2*P_feG2sub,2*P_heG2sub,2*P_tijtij + + def IA_gb2tij(self,P,P_window=None, C_window=None): + P_F2, A = self.J_k_tensor(P,self.X_IA_gb2_F2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_F2 = self.EK.PK_original(P_F2) + P_G2, A = self.J_k_tensor(P,self.X_IA_gb2_G2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_G2 = self.EK.PK_original(P_G2) + P_gb2tij = P_G2-P_F2 + P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_S2F2 = self.EK.PK_original(P_S2F2) + + P_13S2F2 = P_IA_13S2F2(self.k_original, P) + + P_S2G2, A = self.J_k_tensor(P, self.X_IA_gb2_S2G2, P_window=P_window, C_window=C_window) + if (self.extrap): + _, P_S2G2 = self.EK.PK_original(P_S2G2) + P_s2tij=P_S2G2-P_S2F2 + + return 2*P_gb2tij,2*P_s2tij + + def IA_gb2(self,P,P_window=None, C_window=None): P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) if (self.extrap): @@ -774,15 +797,11 @@ def IA_gb2(self,P,P_window=None, C_window=None): P_F2, A = self.J_k_tensor(P,self.X_IA_gb2_F2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_F2 = self.EK.PK_original(P_F2) - P_G2, A = self.J_k_tensor(P,self.X_IA_gb2_G2, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_G2 = self.EK.PK_original(P_G2) sig4 = np.trapz(self.k_original ** 3 * P ** 2, x=np.log(self.k_original)) / (2. * pi ** 2) P_gb2sij = P_F2 P_gb2sij2 = P_he - P_gb2tij = P_G2-P_F2 P_gb2dsij = P_fe - return 2*P_gb2sij, 2*P_gb2dsij, 2*P_gb2sij2, 2*P_gb2tij + return 2*P_gb2sij, 2*P_gb2dsij, 2*P_gb2sij2 def IA_s2(self, P, P_window=None, C_window=None): P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) @@ -791,9 +810,6 @@ def IA_s2(self, P, P_window=None, C_window=None): P_13S2F2 = P_IA_13S2F2(self.k_original, P) - P_S2G2, A = self.J_k_tensor(P, self.X_IA_gb2_S2G2, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_S2G2 = self.EK.PK_original(P_S2G2) P_S2fe, A = self.J_k_tensor(P, self.X_IA_gb2_S2fe, P_window=P_window, C_window=C_window) if (self.extrap): _, P_S2fe = self.EK.PK_original(P_S2fe) @@ -803,8 +819,7 @@ def IA_s2(self, P, P_window=None, C_window=None): P_s2sij=P_S2F2+2*P_13S2F2 P_s2dsij=P_S2fe P_s2sij2=P_S2he - P_s2tij=P_S2G2-P_S2F2 - return 2*P_s2sij, 2*P_s2dsij, 2*P_s2sij2, 2*P_s2tij + return 2*P_s2sij, 2*P_s2dsij, 2*P_s2sij2 def OV(self, P, P_window=None, C_window=None): From ee52d2b2b3dd0a952b74a1a9480546046c8982d9 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 29 Oct 2024 01:53:12 -0400 Subject: [PATCH 15/20] Updated naming convention for tij terms, seperated out tij and gb2 terms --- fastpt/FASTPT.py | 42 ++++-------------------------------------- 1 file changed, 4 insertions(+), 38 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 1ec0bc7..aa01205 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -722,7 +722,7 @@ def IA_der(self, P, P_window=None, C_window=None): P_der = (self.k_original**2)*P return P_der - def IA_tij(self,P,P_window=None, C_window=None): + def IA_ct(self,P,P_window=None, C_window=None): P_feG2, A = self.J_k_tensor(P,self.X_IA_tij_feG2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_feG2 = self.EK.PK_original(P_feG2) @@ -762,10 +762,10 @@ def IA_tij(self,P,P_window=None, C_window=None): P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) - return 2*P_tijsij,2*P_feG2sub,2*P_heG2sub,2*P_tijtij + return 2*P_0TE,2*P_0ETE,2*P_E2TE,2*P_TETE - def IA_gb2tij(self,P,P_window=None, C_window=None): + def IA_ct_mix(self,P,P_window=None, C_window=None): P_F2, A = self.J_k_tensor(P,self.X_IA_gb2_F2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_F2 = self.EK.PK_original(P_F2) @@ -784,43 +784,9 @@ def IA_gb2tij(self,P,P_window=None, C_window=None): _, P_S2G2 = self.EK.PK_original(P_S2G2) P_s2tij=P_S2G2-P_S2F2 - return 2*P_gb2tij,2*P_s2tij + return 2*P_d2TE,2*P_s2TE - def IA_gb2(self,P,P_window=None, C_window=None): - P_fe, A = self.J_k_tensor(P,self.X_IA_gb2_fe, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_fe = self.EK.PK_original(P_fe) - P_he, A = self.J_k_tensor(P,self.X_IA_gb2_he, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_he = self.EK.PK_original(P_he) - P_F2, A = self.J_k_tensor(P,self.X_IA_gb2_F2, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_F2 = self.EK.PK_original(P_F2) - sig4 = np.trapz(self.k_original ** 3 * P ** 2, x=np.log(self.k_original)) / (2. * pi ** 2) - P_gb2sij = P_F2 - P_gb2sij2 = P_he - P_gb2dsij = P_fe - return 2*P_gb2sij, 2*P_gb2dsij, 2*P_gb2sij2 - - def IA_s2(self, P, P_window=None, C_window=None): - P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_S2F2 = self.EK.PK_original(P_S2F2) - - P_13S2F2 = P_IA_13S2F2(self.k_original, P) - - P_S2fe, A = self.J_k_tensor(P, self.X_IA_gb2_S2fe, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_S2fe = self.EK.PK_original(P_S2fe) - P_S2he, A = self.J_k_tensor(P, self.X_IA_gb2_S2he, P_window=P_window, C_window=C_window) - if (self.extrap): - _, P_S2he = self.EK.PK_original(P_S2he) - P_s2sij=P_S2F2+2*P_13S2F2 - P_s2dsij=P_S2fe - P_s2sij2=P_S2he - return 2*P_s2sij, 2*P_s2dsij, 2*P_s2sij2 - def OV(self, P, P_window=None, C_window=None): P, A = self.J_k_tensor(P, self.X_OV, P_window=P_window, C_window=C_window) From 2da24e8fa09426a0f117a5df987e1731bf64587c Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 12 Nov 2024 00:52:09 -0500 Subject: [PATCH 16/20] Continuing seperation of tij and galaxy bias --- fastpt/FASTPT.py | 91 +++++------------------------------------------- 1 file changed, 9 insertions(+), 82 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index aa01205..070462d 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -227,9 +227,6 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.IA_mix_do = True self.IA_tij_do = True continue - elif entry == 'gb2': - self.IA_gb2_do = True - continue elif entry == 'all' or entry == 'everything': self.dd_do = True self.dd_bias_do = True @@ -241,7 +238,6 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.RSD_do = True self.cleft = True self.IA_tij_do = True - self.IA_gb2_do = True continue else: raise ValueError('FAST-PT does not recognize "' + entry + '" in the to_do list.') @@ -304,45 +300,30 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high IA_tij_G2G2_tab = IA_tij_G2G2() IA_tij_F2G2_tab = IA_tij_F2G2() IA_tij_F2G2reg_tab =IA_tij_F2G2reg() + IA_gb2_F2_tab = IA_gb2_F2() + IA_gb2_G2_tab = IA_gb2_G2() + IA_gb2_S2F2_tab = IA_gb2_S2F2() + IA_gb2_S2G2_tab = IA_gb2_S2G2() p_mat_tij_feG2 = IA_tij_feG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_heG2 = IA_tij_heG2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_F2F2 = IA_tij_F2F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_G2G2 = IA_tij_G2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_F2G2 = IA_tij_F2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] p_mat_tij_F2G2reg_tab = IA_tij_F2G2reg_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_F2 = IA_gb2_F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_G2 = IA_gb2_G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_S2F2 = IA_gb2_S2F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] + p_mat_gb2_S2G2 = IA_gb2_S2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] self.X_IA_tij_feG2 = tensor_stuff(p_mat_tij_feG2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_heG2 = tensor_stuff(p_mat_tij_heG2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_F2F2 = tensor_stuff(p_mat_tij_F2F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_G2G2 = tensor_stuff(p_mat_tij_G2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_F2G2 = tensor_stuff(p_mat_tij_F2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_tij_F2G2reg = tensor_stuff(p_mat_tij_F2G2reg_tab, self.N, self.m, self.eta_m, self.l, self.tau_l) - - if self.IA_gb2_do: - IA_gb2_fe_tab = IA_gb2_fe() - IA_gb2_he_tab = IA_gb2_he() - IA_gb2_F2_tab = IA_gb2_F2() - IA_gb2_G2_tab = IA_gb2_G2() - IA_gb2_S2F2_tab = IA_gb2_S2F2() - IA_gb2_S2G2_tab = IA_gb2_S2G2() - IA_gb2_S2fe_tab = IA_gb2_S2fe() - IA_gb2_S2he_tab = IA_gb2_S2he() - p_mat_gb2_fe = IA_gb2_fe_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_he = IA_gb2_he_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_F2 = IA_gb2_F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_G2 = IA_gb2_G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_S2F2 = IA_gb2_S2F2_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_S2G2 = IA_gb2_S2G2_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_S2fe = IA_gb2_S2fe_tab[:, [0, 1, 5, 6, 7, 8, 9]] - p_mat_gb2_S2he = IA_gb2_S2he_tab[:, [0, 1, 5, 6, 7, 8, 9]] - self.X_IA_gb2_fe = tensor_stuff(p_mat_gb2_fe, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_gb2_he = tensor_stuff(p_mat_gb2_he, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_F2 = tensor_stuff(p_mat_gb2_F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_G2 = tensor_stuff(p_mat_gb2_G2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_S2F2 = tensor_stuff(p_mat_gb2_S2F2, self.N, self.m, self.eta_m, self.l, self.tau_l) self.X_IA_gb2_S2G2 = tensor_stuff(p_mat_gb2_S2G2, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_gb2_S2fe = tensor_stuff(p_mat_gb2_S2fe, self.N, self.m, self.eta_m, self.l, self.tau_l) - self.X_IA_gb2_S2he = tensor_stuff(p_mat_gb2_S2he, self.N, self.m, self.eta_m, self.l, self.tau_l) - if self.OV_do: @@ -432,60 +413,6 @@ def one_loop_dd(self, P, P_window=None, C_window=None): _, P_1loop = self.EK.PK_original(P_1loop) return P_1loop, Ps - - def tester_function(self,P,P_window=None,C_window=None): - one_loop_coef = np.array( - [2 * 1219 / 1470., 2 * 671 / 1029., 2 * 32 / 1715., 2 * 1 / 3., 2 * 62 / 35., 2 * 8 / 35., 1 / 3.]) - - # get the roundtrip Fourier power spectrum, i.e. P=IFFT[FFT[P]] - # get the matrix for each J_k component - nu = -2 - Ps, mat = self.J_k_scalar(P, self.X_spt, nu, P_window=P_window, C_window=C_window) - - P22_mat = np.multiply(one_loop_coef, np.transpose(mat)) - P22 = np.sum(P22_mat, 1) - P13 = P_13_reg(self.k_old, Ps) - - #All new terms - P_13F = P_IA_13F(self.k_original, P) - P_13G = P_IA_13G(self.k_original,P) - #P_F2F2, A = self.J_k_tensor(P,self.X_IA_tij_F2F2, P_window=P_window, C_window=C_window) - #if (self.extrap): - # _, P_F2F2 = self.EK.PK_original(P_F2F2) - #P_G2G2, A = self.J_k_tensor(P,self.X_IA_tij_G2G2, P_window=P_window, C_window=C_window) - #if (self.extrap): - # _, P_G2G2 = self.EK.PK_original(P_G2G2) - #param_matrix=np.array([[2,-2,0,1]]) - #Power, mat=J_k(self.k_original,P,param_matrix,P_window=P_window, C_window=C_window, n_pad=n_pad) - #regF = 1/3*mat[0,:] - #regG = 1/9*mat[0,:] - #P_22F=P_F2F2+regF - #P_22G=P_G2G2+regG - P_22F = P_22F_reg(self.k_original,P, P_window=P_window, C_window=C_window, n_pad=self.n_pad) - P_22G = P_22G_reg(self.k_original,P,P_window=P_window,C_window=C_window, n_pad=self.n_pad) - - #Test of different P_22G - one_loop_coefG= np.array( - [2*1003/1470, 2*803/1029, 2*64/1715, 2*1/3, 2*58/35, 2*12/35, 1/3]) - - PsG, matG = self.J_k_scalar(P, self.X_sptG, nu, P_window=P_window, C_window=C_window) - P22G_mat = np.multiply(one_loop_coefG, np.transpose(matG)) - P_22G2 = np.sum(P22G_mat, 1) - #P_22Gnonreg, A = self.J_k_tensor(P,self.X_IA_tij_F2G2reg, P_window=P_window, C_window=C_window) - #param_matrix=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]]) - - - #Power, mat=J_k(self.k_original,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=self.n_pad) - #reg = 1/3*mat[6,:] - - if (self.extrap): - _, P22 = self.EK.PK_original(P22) - _, P13 = self.EK.PK_original(P13) - _, P_22G2 = self.EK.PK_original(P_22G2) - #_, P_22Gnonreg = self.EK.PK_original(P_22Gnonreg) - #P_22G=P_22Gnonreg+reg - return P22, P13, P_22F, P_13F, P_22G, P_13G, P_22G2 @@ -765,7 +692,7 @@ def IA_ct(self,P,P_window=None, C_window=None): return 2*P_0TE,2*P_0ETE,2*P_E2TE,2*P_TETE - def IA_ct_mix(self,P,P_window=None, C_window=None): + def IA_ctbias(self,P,P_window=None, C_window=None): P_F2, A = self.J_k_tensor(P,self.X_IA_gb2_F2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_F2 = self.EK.PK_original(P_F2) From a56cf4758d9429dc7f9780ec1196acf5d4ebebab Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 12 Nov 2024 01:00:29 -0500 Subject: [PATCH 17/20] Renamed tij files to align with function, fully seperated galaxy bias and tij terms --- fastpt/FASTPT.py | 5 +- fastpt/IA_ct.py | 213 +++++++++++++++++++++++++++++ fastpt/{IA_gb2.py => IA_ctbias.py} | 0 3 files changed, 215 insertions(+), 3 deletions(-) create mode 100644 fastpt/IA_ct.py rename fastpt/{IA_gb2.py => IA_ctbias.py} (100%) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 070462d..0126f2e 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -49,9 +49,8 @@ from .IA_tt import IA_tt from .IA_ABD import IA_A, IA_DEE, IA_DBB, P_IA_B from .IA_ta import IA_deltaE1, P_IA_deltaE2, IA_0E0E, IA_0B0B -from .IA_tij import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, P_22F_reg, P_22G_reg, IA_tij_F2G2reg -from .IA_gb2 import IA_gb2_F2, IA_gb2_fe, IA_gb2_G2, IA_gb2_he, P_IA_13S2F2 -from .IA_gb2 import IA_gb2_S2F2, IA_gb2_S2G2, IA_gb2_S2fe, IA_gb2_S2he +from .IA_ct import IA_tij_feG2, IA_tij_heG2, IA_tij_F2F2, IA_tij_G2G2, IA_tij_F2G2, P_IA_13G, P_IA_13F, P_22F_reg, P_22G_reg, IA_tij_F2G2reg +from .IA_ctbias import IA_gb2_F2, IA_gb2_G2, IA_gb2_S2F2, IA_gb2_S2G2 from .J_k import J_k from .OV import OV from .kPol import kPol diff --git a/fastpt/IA_ct.py b/fastpt/IA_ct.py new file mode 100644 index 0000000..1de9e45 --- /dev/null +++ b/fastpt/IA_ct.py @@ -0,0 +1,213 @@ +from __future__ import division +import numpy as np +from .J_table import J_table +from .J_k import J_k +import sys +from time import time +from numpy import log, sqrt, exp, pi +from scipy.signal import fftconvolve as convolve + +def IA_tij_feG2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_feG2=np.array([[0,0,0,2,0,13/21],\ + [0,0,0,2,2,8/21],\ + [1,-1,0,2,1,1/2],\ + [-1,1,0,2,1,1/2]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_feG2.shape[0]): + x=J_table(l_mat_tij_feG2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_heG2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_heG2=np.array([[0,0,0,0,0,-9/70],\ + [0,0,2,0,0,-26/63],\ + [0,0,0,0,2,-15/49],\ + [0,0,2,0,2,-16/63],\ + [0,0,1,1,1,81/70],\ + [0,0,1,1,3,12/35],\ + [0,0,0,0,4,-16/245],\ + [1,-1,0,0,1,-3/10],\ + [1,-1,2,0,1,-1/3],\ + [1,-1,1,1,0,1/2],\ + [1,-1,1,1,2,1],\ + [1,-1,0,2,1,-1/3],\ + [1,-1,0,0,3,-1/5]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_heG2.shape[0]): + x=J_table(l_mat_tij_heG2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_F2F2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_F2F2=np.array([[0,0,0,0,0,1219/1470],\ + [0,0,0,0,2,671/1029],\ + [0,0,0,0,4,32/1715],\ + [2,-2,0,0,0,1/6],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,62/35],\ + [1,-1,0,0,4,8/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_F2F2.shape[0]): + x=J_table(l_mat_tij_F2F2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_G2G2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_G2G2=np.array([[0,0,0,0,0,851/1470],\ + [0,0,0,0,2,871/1029],\ + [0,0,0,0,4,128/1715],\ + [2,-2,0,0,0,1/6],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,54/35],\ + [1,-1,0,0,4,16/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_G2G2.shape[0]): + x=J_table(l_mat_tij_G2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def IA_tij_F2G2(): + # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient + l_mat_tij_F2G2=np.array([[0,0,0,0,0,1003/1470],\ + [0,0,0,0,2,803/1029],\ + [0,0,0,0,4,64/1715],\ + [2,-2,0,0,0,1/6],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,58/35],\ + [1,-1,0,0,4,12/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_F2G2.shape[0]): + x=J_table(l_mat_tij_F2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + +def P_IA_13G(k,P): + N=k.size + n = np.arange(-N+1,N ) + dL=log(k[1])-log(k[0]) + s=n*dL + cut=7 + 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)] + + # For Zbar + Z1=lambda r : (12/r**2-26+4*r**2-6*r**4+(3/r**3)*(r**2-1)**3*(r**2+2)*log((r+1.)/np.absolute(r-1.))) *r + Z1_low=lambda r : (-224/5+1248/(35*r**2)-608/(105*r**4)-26/(21*r**6)-4/(35*r**8)+34/(21*r**10) -4/(3*r**12)) *r + Z1_high=lambda r: ((-32*r**2)/5-(96*r**4)/7+(352*r**6)/105+(164*r**8)/105-(58*r**10)/35+(4*r**12)/21+(2*r**14)/3) *r + + f_mid_low=Z1(exp(-mid_low_s)) + f_mid_high=Z1(exp(-mid_high_s)) + f_high = Z1_high(exp(-high_s)) + f_low = Z1_low(exp(-low_s)) + + f=np.hstack((f_low,f_mid_low,-16.,f_mid_high,f_high)) + # print(f) + + g= convolve(P, f) * dL + g_k=g[N-1:2*N-1] + deltaE2= 1/84 * k**3/(2*pi)**2 * P*g_k + return deltaE2 + +def P_22F_reg(k,P, P_window, C_window, n_pad): + # P_22 Legendre components + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. + + param_matrix=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] ]) + + + Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) + A=1219/1470.*mat[0,:] + B=671/1029.*mat[1,:] + C=32/1715.*mat[2,:] + D=1/3.*mat[3,:] + E=62/35.*mat[4,:] + F=8/35.*mat[5,:] + reg=1/3.*mat[6,:] + + return 2*(A+B+C+D+E+F)+ reg + +def P_22G_reg(k,P, P_window, C_window, n_pad): + # P_22 Legendre components + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. + + param_matrix=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]]) + + + Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) + A=1003/1470.*mat[0,:] + B=803/1029.*mat[1,:] + C=64/1715.*mat[2,:] + D=1/3.*mat[3,:] + E=58/35.*mat[4,:] + F=12/35.*mat[5,:] + reg=1/3.*mat[6,:] + + + return 2*(A+B+C+D+E+F)+ reg + +def IA_tij_F2G2reg(): + # P_22 Legendre components + # We calculate a regularized version of P_22 + # by omitting the J_{2,-2,0} term so that the + # integral converges. In the final power spectrum + # we add the asymptotic portions of P_22 and P_13 so + # that we get a convergent integral. See section of XXX. + + l_mat_tij_F2G2=np.array([[0,0,0,0,0,1003/1470],\ + [0,0,0,0,2,803/1029],\ + [0,0,0,0,4,64/1715],\ + [2,-2,0,0,2,1/3],\ + [1,-1,0,0,1,58/35],\ + [1,-1,0,0,4,12/35]], dtype=float) + table=np.zeros(10,dtype=float) + for i in range(l_mat_tij_F2G2.shape[0]): + x=J_table(l_mat_tij_F2G2[i]) + table=np.row_stack((table,x)) + return table[1:,:] + + +def P_IA_13F(k,P): + N=k.size + n= np.arange(-N+1,N ) + dL=log(k[1])-log(k[0]) + s=n*dL + + cut=7 + 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 : (12./r**2 +10. + 100.*r**2-42.*r**4 \ + + 3./r**3*(r**2-1.)**3*(7*r**2+2.)*log((r+1.)/np.absolute(r-1.)) ) *r + Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8-50./21./r**10-4./3./r**12) *r + Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8+74./35.*r**10-20./3.*r**12+14./3.*r**14) *r + + 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= convolve(P, f) * dL + g_k=g[N-1:2*N-1] + P_bar= 1./252.* k**3/(2*pi)**2*P*g_k + + return P_bar diff --git a/fastpt/IA_gb2.py b/fastpt/IA_ctbias.py similarity index 100% rename from fastpt/IA_gb2.py rename to fastpt/IA_ctbias.py From 915d0bf2dae1b212eb7b80a8e1a279e674d56d53 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Tue, 12 Nov 2024 01:00:54 -0500 Subject: [PATCH 18/20] Still updating ctbias --- fastpt/IA_ctbias.py | 81 +---------------- fastpt/IA_tij.py | 213 -------------------------------------------- 2 files changed, 1 insertion(+), 293 deletions(-) delete mode 100644 fastpt/IA_tij.py diff --git a/fastpt/IA_ctbias.py b/fastpt/IA_ctbias.py index 18872a3..7a189cd 100644 --- a/fastpt/IA_ctbias.py +++ b/fastpt/IA_ctbias.py @@ -6,62 +6,6 @@ from numpy import log, sqrt, exp, pi from scipy.signal import fftconvolve as convolve -def P_IA_13S2F2(k,P): - - N=k.size - n= np.arange(-N+1,N ) - dL=log(k[1])-log(k[0]) - s=n*dL - cut=7 - 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)] - - - - Z1=lambda r : (((4.*r*(45.-165.*r**2+379.*r**4+45.*r**6)+45.*(-1.+r**2)**4*log((-1.+r)**2)-90.*(-1.+r**2)**4*log(np.absolute(1.+r)))/(2016.*r**3))-68./63*r**2)/2. - Z1_high=lambda r : ((-16*r**2)/21. + (16*r**4)/49. - (16*r**6)/441. - (16*r**8)/4851. - 16*r**10/21021.- 16*r**12/63063.)/2. - Z1_low=lambda r: (-16/21.+ 16/(49.*r**2) - 16/(441.*r**4) - 16/(4851.*r**6) - 16/(21021.*r**8) - 16/(63063.*r**10)- 16/(153153.*r**12) )/2. - - - f_mid_low=Z1(exp(-mid_low_s))*exp(-mid_low_s) - f_mid_high=Z1(exp(-mid_high_s))*exp(-mid_high_s) - f_high = Z1_high(exp(-high_s))*exp(-high_s) - f_low = Z1_low(exp(-low_s))*exp(-low_s) - - f=np.hstack((f_low,f_mid_low,-0.2381002916036672,f_mid_high,f_high)) - - - g= convolve(P, f) * dL - g_k=g[N-1:2*N-1] - S2F2= k**3/(2.*pi**2) * P*g_k - return S2F2 - - -def IA_gb2_fe(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_gb2_fe=np.array([[0,0,2,0,0,1]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_gb2_fe.shape[0]): - x=J_table(l_mat_gb2_fe[i]) - table=np.row_stack((table,x)) - return table[1:,:] - - -def IA_gb2_he(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_gb2_he=np.array([[0,0,0,0,0,-1/6],\ - [0,0,2,0,0,-1/3],\ - [0,0,0,0,2,-1/3],\ - [0,0,1,1,1,3/2],\ - [0,0,2,0,0,-1/3]],dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_gb2_he.shape[0]): - x=J_table(l_mat_gb2_he[i]) - table=np.row_stack((table,x)) - return table[1:,:] - def IA_gb2_F2(): # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient l_mat_gb2_F2=np.array([[0,0,0,0,0,17/21],\ @@ -94,6 +38,7 @@ def IA_gb2_S2F2(): [1,-1,0,0,1,2/15],\ [1,-1,0,0,3,1/5],\ [-1,1,0,0,1,2/15],\ + [-1,1,0,0,1,1/5]],dtype=float) [-1,1,0,0,3,1/5]],dtype=float) table=np.zeros(10,dtype=float) for i in range(l_mat_gb2_S2F2.shape[0]): @@ -114,28 +59,4 @@ def IA_gb2_S2G2(): for i in range(l_mat_gb2_S2G2.shape[0]): x=J_table(l_mat_gb2_S2G2[i]) table=np.row_stack((table,x)) - return table[1:,:] - -def IA_gb2_S2fe(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_gb2_S2fe=np.array([[0,0,2,0,2,2/3]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_gb2_S2fe.shape[0]): - x=J_table(l_mat_gb2_S2fe[i]) - table=np.row_stack((table,x)) - return table[1:,:] - -def IA_gb2_S2he(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_gb2_S2he=np.array([[0,0,0,0,0,-2/45],\ - [0,0,0,0,2,-11/63],\ - [0,0,2,0,2,-2/9],\ - [0,0,0,2,2,-2/9],\ - [0,0,1,1,1,2/5],\ - [0,0,1,1,3,3/5],\ - [0,0,0,0,4,-4/35]],dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_gb2_S2he.shape[0]): - x=J_table(l_mat_gb2_S2he[i]) - table=np.row_stack((table,x)) return table[1:,:] \ No newline at end of file diff --git a/fastpt/IA_tij.py b/fastpt/IA_tij.py deleted file mode 100644 index 1de9e45..0000000 --- a/fastpt/IA_tij.py +++ /dev/null @@ -1,213 +0,0 @@ -from __future__ import division -import numpy as np -from .J_table import J_table -from .J_k import J_k -import sys -from time import time -from numpy import log, sqrt, exp, pi -from scipy.signal import fftconvolve as convolve - -def IA_tij_feG2(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_tij_feG2=np.array([[0,0,0,2,0,13/21],\ - [0,0,0,2,2,8/21],\ - [1,-1,0,2,1,1/2],\ - [-1,1,0,2,1,1/2]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_tij_feG2.shape[0]): - x=J_table(l_mat_tij_feG2[i]) - table=np.row_stack((table,x)) - return table[1:,:] - -def IA_tij_heG2(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_tij_heG2=np.array([[0,0,0,0,0,-9/70],\ - [0,0,2,0,0,-26/63],\ - [0,0,0,0,2,-15/49],\ - [0,0,2,0,2,-16/63],\ - [0,0,1,1,1,81/70],\ - [0,0,1,1,3,12/35],\ - [0,0,0,0,4,-16/245],\ - [1,-1,0,0,1,-3/10],\ - [1,-1,2,0,1,-1/3],\ - [1,-1,1,1,0,1/2],\ - [1,-1,1,1,2,1],\ - [1,-1,0,2,1,-1/3],\ - [1,-1,0,0,3,-1/5]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_tij_heG2.shape[0]): - x=J_table(l_mat_tij_heG2[i]) - table=np.row_stack((table,x)) - return table[1:,:] - -def IA_tij_F2F2(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_tij_F2F2=np.array([[0,0,0,0,0,1219/1470],\ - [0,0,0,0,2,671/1029],\ - [0,0,0,0,4,32/1715],\ - [2,-2,0,0,0,1/6],\ - [2,-2,0,0,2,1/3],\ - [1,-1,0,0,1,62/35],\ - [1,-1,0,0,4,8/35]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_tij_F2F2.shape[0]): - x=J_table(l_mat_tij_F2F2[i]) - table=np.row_stack((table,x)) - return table[1:,:] - -def IA_tij_G2G2(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_tij_G2G2=np.array([[0,0,0,0,0,851/1470],\ - [0,0,0,0,2,871/1029],\ - [0,0,0,0,4,128/1715],\ - [2,-2,0,0,0,1/6],\ - [2,-2,0,0,2,1/3],\ - [1,-1,0,0,1,54/35],\ - [1,-1,0,0,4,16/35]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_tij_G2G2.shape[0]): - x=J_table(l_mat_tij_G2G2[i]) - table=np.row_stack((table,x)) - return table[1:,:] - -def IA_tij_F2G2(): - # Ordering is \alpha, \beta, l_1, l_2, l, A coeficient - l_mat_tij_F2G2=np.array([[0,0,0,0,0,1003/1470],\ - [0,0,0,0,2,803/1029],\ - [0,0,0,0,4,64/1715],\ - [2,-2,0,0,0,1/6],\ - [2,-2,0,0,2,1/3],\ - [1,-1,0,0,1,58/35],\ - [1,-1,0,0,4,12/35]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_tij_F2G2.shape[0]): - x=J_table(l_mat_tij_F2G2[i]) - table=np.row_stack((table,x)) - return table[1:,:] - -def P_IA_13G(k,P): - N=k.size - n = np.arange(-N+1,N ) - dL=log(k[1])-log(k[0]) - s=n*dL - cut=7 - 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)] - - # For Zbar - Z1=lambda r : (12/r**2-26+4*r**2-6*r**4+(3/r**3)*(r**2-1)**3*(r**2+2)*log((r+1.)/np.absolute(r-1.))) *r - Z1_low=lambda r : (-224/5+1248/(35*r**2)-608/(105*r**4)-26/(21*r**6)-4/(35*r**8)+34/(21*r**10) -4/(3*r**12)) *r - Z1_high=lambda r: ((-32*r**2)/5-(96*r**4)/7+(352*r**6)/105+(164*r**8)/105-(58*r**10)/35+(4*r**12)/21+(2*r**14)/3) *r - - f_mid_low=Z1(exp(-mid_low_s)) - f_mid_high=Z1(exp(-mid_high_s)) - f_high = Z1_high(exp(-high_s)) - f_low = Z1_low(exp(-low_s)) - - f=np.hstack((f_low,f_mid_low,-16.,f_mid_high,f_high)) - # print(f) - - g= convolve(P, f) * dL - g_k=g[N-1:2*N-1] - deltaE2= 1/84 * k**3/(2*pi)**2 * P*g_k - return deltaE2 - -def P_22F_reg(k,P, P_window, C_window, n_pad): - # P_22 Legendre components - # We calculate a regularized version of P_22 - # by omitting the J_{2,-2,0} term so that the - # integral converges. In the final power spectrum - # we add the asymptotic portions of P_22 and P_13 so - # that we get a convergent integral. See section of XXX. - - param_matrix=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] ]) - - - Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) - A=1219/1470.*mat[0,:] - B=671/1029.*mat[1,:] - C=32/1715.*mat[2,:] - D=1/3.*mat[3,:] - E=62/35.*mat[4,:] - F=8/35.*mat[5,:] - reg=1/3.*mat[6,:] - - return 2*(A+B+C+D+E+F)+ reg - -def P_22G_reg(k,P, P_window, C_window, n_pad): - # P_22 Legendre components - # We calculate a regularized version of P_22 - # by omitting the J_{2,-2,0} term so that the - # integral converges. In the final power spectrum - # we add the asymptotic portions of P_22 and P_13 so - # that we get a convergent integral. See section of XXX. - - param_matrix=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]]) - - - Power, mat=J_k(k,P,param_matrix, P_window=P_window, C_window=C_window, n_pad=n_pad) - A=1003/1470.*mat[0,:] - B=803/1029.*mat[1,:] - C=64/1715.*mat[2,:] - D=1/3.*mat[3,:] - E=58/35.*mat[4,:] - F=12/35.*mat[5,:] - reg=1/3.*mat[6,:] - - - return 2*(A+B+C+D+E+F)+ reg - -def IA_tij_F2G2reg(): - # P_22 Legendre components - # We calculate a regularized version of P_22 - # by omitting the J_{2,-2,0} term so that the - # integral converges. In the final power spectrum - # we add the asymptotic portions of P_22 and P_13 so - # that we get a convergent integral. See section of XXX. - - l_mat_tij_F2G2=np.array([[0,0,0,0,0,1003/1470],\ - [0,0,0,0,2,803/1029],\ - [0,0,0,0,4,64/1715],\ - [2,-2,0,0,2,1/3],\ - [1,-1,0,0,1,58/35],\ - [1,-1,0,0,4,12/35]], dtype=float) - table=np.zeros(10,dtype=float) - for i in range(l_mat_tij_F2G2.shape[0]): - x=J_table(l_mat_tij_F2G2[i]) - table=np.row_stack((table,x)) - return table[1:,:] - - -def P_IA_13F(k,P): - N=k.size - n= np.arange(-N+1,N ) - dL=log(k[1])-log(k[0]) - s=n*dL - - cut=7 - 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 : (12./r**2 +10. + 100.*r**2-42.*r**4 \ - + 3./r**3*(r**2-1.)**3*(7*r**2+2.)*log((r+1.)/np.absolute(r-1.)) ) *r - Z_low=lambda r : (352./5.+96./.5/r**2 -160./21./r**4 - 526./105./r**6 +236./35./r**8-50./21./r**10-4./3./r**12) *r - Z_high=lambda r: (928./5.*r**2 - 4512./35.*r**4 +416./21.*r**6 +356./105.*r**8+74./35.*r**10-20./3.*r**12+14./3.*r**14) *r - - 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= convolve(P, f) * dL - g_k=g[N-1:2*N-1] - P_bar= 1./252.* k**3/(2*pi)**2*P*g_k - - return P_bar From 34efd06ce8bb8244d8ed8c6bc7fb1785cb9852b1 Mon Sep 17 00:00:00 2001 From: CarterDW Date: Thu, 21 Nov 2024 02:18:59 -0500 Subject: [PATCH 19/20] Edited naming convention for nlxia tij terms --- fastpt/FASTPT.py | 20 +++++++++----------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/fastpt/FASTPT.py b/fastpt/FASTPT.py index 0126f2e..ae92313 100644 --- a/fastpt/FASTPT.py +++ b/fastpt/FASTPT.py @@ -176,7 +176,6 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.kPol_do = False self.RSD_do = False self.IA_tij_do = False - self.IA_gb2_do = False for entry in to_do: # convert to_do list to instructions for FAST-PT initialization if entry == 'one_loop_dd': @@ -196,7 +195,6 @@ def __init__(self, k, nu=None, to_do=None, param_mat=None, low_extrap=None, high self.IA_ta_do = True self.IA_mix_do = True self.IA_tij_do = True - self.IA_gb2_do = True continue elif entry == 'IA_tt': self.IA_tt_do = True @@ -683,12 +681,12 @@ def IA_ct(self,P,P_window=None, C_window=None): if (self.extrap): _, P_22F=self.EK.PK_original(P_22F) _, P_22G=self.EK.PK_original(P_22G) - P_tijtij = P_F2F2+P_G2G2-2*P_F2G2 - P_tijsij = P_22G-P_22F+P_13G-P_13F - P_feG2sub = np.subtract(P_feG2,(1/2)*P_A00E) - P_heG2sub = np.subtract(P_heG2,(1/2)*P_A0E2) + P_tEtE = P_F2F2+P_G2G2-2*P_F2G2 + P_0tE = P_22G-P_22F+P_13G-P_13F + P_0EtE = np.subtract(P_feG2,(1/2)*P_A00E) + P_E2tE = np.subtract(P_heG2,(1/2)*P_A0E2) - return 2*P_0TE,2*P_0ETE,2*P_E2TE,2*P_TETE + return 2*P_0tE,2*P_0EtE,2*P_E2tE,2*P_tEtE def IA_ctbias(self,P,P_window=None, C_window=None): @@ -698,19 +696,19 @@ def IA_ctbias(self,P,P_window=None, C_window=None): P_G2, A = self.J_k_tensor(P,self.X_IA_gb2_G2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_G2 = self.EK.PK_original(P_G2) - P_gb2tij = P_G2-P_F2 + P_d2tE = P_G2-P_F2 P_S2F2, A = self.J_k_tensor(P, self.X_IA_gb2_S2F2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_S2F2 = self.EK.PK_original(P_S2F2) - P_13S2F2 = P_IA_13S2F2(self.k_original, P) + #P_13S2F2 = P_IA_13S2F2(self.k_original, P) P_S2G2, A = self.J_k_tensor(P, self.X_IA_gb2_S2G2, P_window=P_window, C_window=C_window) if (self.extrap): _, P_S2G2 = self.EK.PK_original(P_S2G2) - P_s2tij=P_S2G2-P_S2F2 + P_s2tE=P_S2G2-P_S2F2 - return 2*P_d2TE,2*P_s2TE + return 2*P_d2tE,2*P_s2tE From b388cf0a395df7f5d886bd8d0524820b3f2ded3f Mon Sep 17 00:00:00 2001 From: CarterDW Date: Wed, 19 Feb 2025 10:11:39 -0500 Subject: [PATCH 20/20] removed testing file which should not of been present --- fastpt/fastpttest.py | 87 -------------------------------------------- 1 file changed, 87 deletions(-) delete mode 100644 fastpt/fastpttest.py diff --git a/fastpt/fastpttest.py b/fastpt/fastpttest.py deleted file mode 100644 index b137b31..0000000 --- a/fastpt/fastpttest.py +++ /dev/null @@ -1,87 +0,0 @@ -import fastpt as pt -from time import time -import matplotlib.pyplot as plt -import numpy as np - - # Version check - - # load the data file - -d = np.loadtxt('/home/carterw/FAST-PT/examples/Pk_test.dat') - # declare k and the power spectrum -k = d[:, 0]; -P = d[:, 1] - - # set the parameters for the power spectrum window and - # Fourier coefficient window - # P_window=np.array([.2,.2]) -C_window = .75 - # document this better in the user manual - - # padding length -n_pad = int(0.5 * len(k)) - # to_do=['one_loop_dd','IA_tt'] -to_do = ['one_loop_dd'] - # to_do=['dd_bias','IA_all'] - # to_do=['all'] - - # initialize the FASTPT class - # including extrapolation to higher and lower k -t1 = time() -fpt = pt(k, to_do=to_do, low_extrap=-5, high_extrap=3, n_pad=n_pad) - -t2 = time() - # calculate 1loop SPT (and time the operation) - # P_spt=fastpt.one_loop_dd(P,C_window=C_window) -P_lpt = fpt.one_loop_dd_bias_lpt(P, C_window=C_window) -P_der = fpt.IA_der(P,C_window=C_window) - - # for M = 10**14 M_sun/h -b1L = 1.02817 -b2L = -0.0426292 -b3L = -2.55751 -b1E = 1 + b1L - - # for M = 10**14 M_sun/h - # b1_lag = 1.1631 - # b2_lag = 0.1162 - - # [Ps, Pnb, Pb1L, Pb1L_2, Pb1L_b2L, Pb2L, Pb2L_2, Pb3L, Pb1L_b3L] = [P_lpt[0],P_lpt[1],P_lpt[2],P_lpt[3],P_lpt[4],P_lpt[5],P_lpt[6],P_lpt[7],P_lpt[8]] -[Ps, Pnb, Pb1L, Pb1L_2, Pb1L_b2L, Pb2L, Pb2L_2] = [P_lpt[0], P_lpt[1], P_lpt[2], P_lpt[3], P_lpt[4], P_lpt[5], - P_lpt[6]] - - # Pgg_lpt = (b1E**2)*P + Pnb + (b1L)*(Pb1L) + (b1L**2)*Pb1L_2 + (b1L*b2L)*Pb1L_b2L + (b2L)*(Pb2L) + (b2L**2)*Pb2L_2 + (b3L)*(Pb3L) + (b1L*b3L)*Pb1L_b3L -Pgg_lpt = (b1E ** 2) * P + Pnb + (b1L) * (Pb1L) + (b1L ** 2) * Pb1L_2 + (b1L * b2L) * Pb1L_b2L + (b2L) * (Pb2L) + ( - b2L ** 2) * Pb2L_2 - - # print([pnb,pb1L,pb1L_2,pb2L,Pb1L_b2L]) - -t3 = time() -print('initialization time for', to_do, "%10.3f" % (t2 - t1), 's') -print('one_loop_dd recurring time', "%10.3f" % (t3 - t2), 's') - - # calculate tidal torque EE and BB P(k) - # P_IA_tt=fastpt.IA_tt(P,C_window=C_window) - # P_IA_ta=fastpt.IA_ta(P,C_window=C_window) - # P_IA_mix=fastpt.IA_mix(P,C_window=C_window) - # P_RSD=fastpt.RSD_components(P,1.0,C_window=C_window) - # P_kPol=fastpt.kPol(P,C_window=C_window) - # P_OV=fastpt.OV(P,C_window=C_window) - # P_IRres=fastpt.IRres(P,C_window=C_window) - # make a plot of 1loop SPT results - - -ax = plt.subplot(111) -ax.set_xscale('log') -ax.set_yscale('log') -ax.set_ylabel(r'$P(k)$', size=30) -ax.set_xlabel(r'$k$', size=30) - -ax.plot(k, P, label='linear') -# ax.plot(k,P_spt[0], label=r'$P_{22}(k) + P_{13}(k)$' ) -ax.plot(k, Pgg_lpt, label='P_lpt') -ax.plot(k,P_der, label='P_der') - -plt.legend(loc=3) -plt.grid() -plt.show() \ No newline at end of file