diff --git a/flip/covariance/adamsblake17plane/coefficients.py b/flip/covariance/adamsblake17plane/coefficients.py index 3ad548d..1cdee6a 100644 --- a/flip/covariance/adamsblake17plane/coefficients.py +++ b/flip/covariance/adamsblake17plane/coefficients.py @@ -2,6 +2,8 @@ def get_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): coefficients_dict = {} if model_type in ["density", "full", "density_velocity"]: diff --git a/flip/covariance/adamsblake17plane/fisher_terms.py b/flip/covariance/adamsblake17plane/fisher_terms.py index 2265c14..3a8cb58 100644 --- a/flip/covariance/adamsblake17plane/fisher_terms.py +++ b/flip/covariance/adamsblake17plane/fisher_terms.py @@ -5,6 +5,8 @@ def get_partial_derivative_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): if model_type == "density": return get_partial_derivative_coefficients_density( diff --git a/flip/covariance/adamsblake17plane/flip_terms.py b/flip/covariance/adamsblake17plane/flip_terms.py index 513c5a8..e0d88fc 100644 --- a/flip/covariance/adamsblake17plane/flip_terms.py +++ b/flip/covariance/adamsblake17plane/flip_terms.py @@ -57,3 +57,4 @@ def N_vv_0_2_0(theta, phi): "vv_0_2": 1, } multi_index_model = False +redshift_dependent_model = False diff --git a/flip/covariance/adamsblake17plane/generator.py b/flip/covariance/adamsblake17plane/generator.py index c4c2d51..86529f9 100644 --- a/flip/covariance/adamsblake17plane/generator.py +++ b/flip/covariance/adamsblake17plane/generator.py @@ -253,5 +253,11 @@ def generate_covariance( ) los_definition = "bisector" - - return covariance_dict, number_densities, number_velocities, los_definition + redshift_dict = None + return ( + covariance_dict, + number_densities, + number_velocities, + los_definition, + redshift_dict, + ) diff --git a/flip/covariance/adamsblake20/coefficients.py b/flip/covariance/adamsblake20/coefficients.py index dfe84ce..cfa8395 100644 --- a/flip/covariance/adamsblake20/coefficients.py +++ b/flip/covariance/adamsblake20/coefficients.py @@ -2,6 +2,8 @@ def get_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): coefficients_dict = {} if model_type in ["density", "full", "density_velocity"]: diff --git a/flip/covariance/adamsblake20/fisher_terms.py b/flip/covariance/adamsblake20/fisher_terms.py index b185a0b..110e4c7 100644 --- a/flip/covariance/adamsblake20/fisher_terms.py +++ b/flip/covariance/adamsblake20/fisher_terms.py @@ -5,6 +5,8 @@ def get_partial_derivative_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): if model_type == "density": return get_partial_derivative_coefficients_density( diff --git a/flip/covariance/adamsblake20/flip_terms.py b/flip/covariance/adamsblake20/flip_terms.py index 601770f..453ad45 100644 --- a/flip/covariance/adamsblake20/flip_terms.py +++ b/flip/covariance/adamsblake20/flip_terms.py @@ -72,11 +72,7 @@ def M_gg_1_2_0(sig_g): def func(k): return ( -5 * np.exp(-(k**2) * sig_g**2) / (k**2 * sig_g**2) - - 5 - / 4 - * np.sqrt(np.pi) - * scipy.special.erf(k * sig_g) - / (k**3 * sig_g**3) + - 5 / 4 * np.sqrt(np.pi) * scipy.special.erf(k * sig_g) / (k**3 * sig_g**3) - 45 / 4 * np.exp(-(k**2) * sig_g**2) / (k**4 * sig_g**4) + (45 / 8) * np.sqrt(np.pi) @@ -359,3 +355,4 @@ def N_vv_0_2_0(theta, phi): "vv_0_2": 1, } multi_index_model = False +redshift_dependent_model = False diff --git a/flip/covariance/carreres23/coefficients.py b/flip/covariance/carreres23/coefficients.py index 7436fab..36b5dc7 100644 --- a/flip/covariance/carreres23/coefficients.py +++ b/flip/covariance/carreres23/coefficients.py @@ -5,6 +5,8 @@ def get_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): coefficients_dict = {} coefficients_dict["vv"] = [parameter_values_dict["fs8"] ** 2] diff --git a/flip/covariance/carreres23/fisher_terms.py b/flip/covariance/carreres23/fisher_terms.py index ed2e119..6c16d23 100644 --- a/flip/covariance/carreres23/fisher_terms.py +++ b/flip/covariance/carreres23/fisher_terms.py @@ -5,6 +5,8 @@ def get_partial_derivative_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): if variant == "growth_index": partial_coefficients_dict = { diff --git a/flip/covariance/carreres23/flip_terms.py b/flip/covariance/carreres23/flip_terms.py index f00c9d5..c4c23d2 100644 --- a/flip/covariance/carreres23/flip_terms.py +++ b/flip/covariance/carreres23/flip_terms.py @@ -28,3 +28,4 @@ def N_vv_0_2_0(theta, phi): dictionary_lmax = {"vv": [2]} dictionary_subterms = {"vv_0_0": 1, "vv_0_1": 0, "vv_0_2": 1} multi_index_model = False +redshift_dependent_model = False \ No newline at end of file diff --git a/flip/covariance/carreres23/generator.py b/flip/covariance/carreres23/generator.py index 997cf57..ccdc07d 100644 --- a/flip/covariance/carreres23/generator.py +++ b/flip/covariance/carreres23/generator.py @@ -121,10 +121,12 @@ def generate_covariance( ) los_definition = "bisector" + redshift_dict = None return ( {"vv": np.array([cov_vv])}, number_densities, number_velocities, los_definition, + redshift_dict, ) diff --git a/flip/covariance/contraction.py b/flip/covariance/contraction.py index 7a7082d..89f6fac 100644 --- a/flip/covariance/contraction.py +++ b/flip/covariance/contraction.py @@ -17,6 +17,8 @@ def __init__( coordinates_dict=None, basis_definition=None, endpoint_los_definition=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, variant=None, ): self.model_name = model_name @@ -25,6 +27,8 @@ def __init__( self.coordinates_dict = coordinates_dict self.basis_definition = basis_definition self.endpoint_los_definition = endpoint_los_definition + self.redshift_dict = redshift_dict + self.power_spectrum_amplitude_function = power_spectrum_amplitude_function self.variant = variant @classmethod @@ -41,12 +45,15 @@ def init_from_flip( additional_parameters_values=None, basis_definition="bisector", endpoint_los_definition="bisector", + redshift=None, + power_spectrum_amplitude_function=None, variant=None, **kwargs, ): ( contraction_dict, coordinates_dict, + redshift_dict, ) = contract_covariance( model_name, model_type, @@ -59,6 +66,7 @@ def init_from_flip( additional_parameters_values=additional_parameters_values, basis_definition=basis_definition, endpoint_los_definition=endpoint_los_definition, + redshift=redshift, **kwargs, ) @@ -69,6 +77,8 @@ def init_from_flip( coordinates_dict=coordinates_dict, basis_definition=basis_definition, endpoint_los_definition=endpoint_los_definition, + redshift_dict=redshift_dict, + power_spectrum_amplitude_function=power_spectrum_amplitude_function, variant=variant, ) @@ -125,7 +135,9 @@ def compute_contraction_sum( coefficients_dict = coefficients.get_coefficients( self.model_type, parameter_values_dict, - self.variant, + variant=self.variant, + redshift_dict=self.redshift_dict, + power_spectrum_amplitude_function=self.power_spectrum_amplitude_function, ) contraction_covariance_sum_dict = {} if self.model_type == "density": @@ -293,6 +305,7 @@ def contract_covariance( additional_parameters_values=None, basis_definition="bisector", endpoint_los_definition="bisector", + redshift=None, number_worker=8, hankel=True, ): @@ -339,5 +352,10 @@ def contract_covariance( number_worker=number_worker, hankel=hankel, )[:, 1:].reshape(-1, len(coord_1), len(coord_2)) - - return contraction_dict, coordinates_dict + redshift_dict = generator_flip.generate_redshift_dict( + model_name, + model_type, + redshift_velocity=redshift, + redshift_density=redshift, + ) + return contraction_dict, coordinates_dict, redshift_dict diff --git a/flip/covariance/covariance.py b/flip/covariance/covariance.py index 44384c2..e37bc8a 100644 --- a/flip/covariance/covariance.py +++ b/flip/covariance/covariance.py @@ -20,6 +20,8 @@ def __init__( full_matrix=False, number_densities=None, number_velocities=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, variant=None, ): """ @@ -49,6 +51,8 @@ def __init__( self.full_matrix = full_matrix self.number_densities = number_densities self.number_velocities = number_velocities + self.redshift_dict = redshift_dict + self.power_spectrum_amplitude_function = power_spectrum_amplitude_function self.variant = variant @classmethod @@ -61,6 +65,7 @@ def init_from_flip( coordinates_velocity=None, additional_parameters_values=None, los_definition="bisector", + power_spectrum_amplitude_function=None, variant=None, **kwargs, ): @@ -93,6 +98,7 @@ def init_from_flip( covariance_dict, number_densities, number_velocities, + redshift_dict, ) = generator_flip.generate_covariance( model_name, model_type, @@ -115,6 +121,8 @@ def init_from_flip( full_matrix=False, number_densities=number_densities, number_velocities=number_velocities, + redshift_dict=redshift_dict, + power_spectrum_amplitude_function=power_spectrum_amplitude_function, variant=variant, ) @@ -127,6 +135,7 @@ def init_from_generator( coordinates_velocity=None, coordinates_density=None, additional_parameters_values=None, + power_spectrum_amplitude_function=None, variant=None, **kwargs, ): @@ -162,6 +171,7 @@ def init_from_generator( number_densities, number_velocities, los_definition, + redshift_dict, ) = generator.generate_covariance( model_type, power_spectrum_dict, @@ -181,6 +191,8 @@ def init_from_generator( full_matrix=False, number_densities=number_densities, number_velocities=number_velocities, + redshift_dict=redshift_dict, + power_spectrum_amplitude_function=power_spectrum_amplitude_function, variant=variant, ) @@ -299,6 +311,8 @@ def compute_covariance_sum( self.model_type, parameter_values_dict, variant=self.variant, + redshift_dict=self.redshift_dict, + power_spectrum_amplitude_function=self.power_spectrum_amplitude_function, ) coefficients_dict_diagonal = coefficients.get_diagonal_coefficients( self.model_type, diff --git a/flip/covariance/generator.py b/flip/covariance/generator.py index 83554ba..482196c 100644 --- a/flip/covariance/generator.py +++ b/flip/covariance/generator.py @@ -12,6 +12,7 @@ from flip.covariance.adamsblake20 import flip_terms as flip_terms_adamsblake20 from flip.covariance.carreres23 import flip_terms as flip_terms_carreres23 from flip.covariance.lai22 import flip_terms as flip_terms_lai22 +from flip.covariance.rcrk24 import flip_terms as flip_terms_rcrk24 from flip.covariance.ravouxcarreres import flip_terms as flip_terms_ravouxcarreres from flip.utils import create_log @@ -22,6 +23,7 @@ "lai22", "carreres23", "ravouxcarreres", + "rcrk24", ] @@ -45,10 +47,7 @@ def correlation_integration(l, r, k, integrand): """ kr = np.outer(k, r) integrand = ( - (-1) ** (l // 2) - * (k**2 / (2 * np.pi**2)) - * integrand - * spherical_jn(l, kr).T + (-1) ** (l // 2) * (k**2 / (2 * np.pi**2)) * integrand * spherical_jn(l, kr).T ) return (-1) ** (l % 2) * integrate.simps(integrand, x=k) @@ -409,6 +408,49 @@ def compute_cov( return covariance +def generate_redshift_dict( + model_name, + model_type, + redshift_velocity=None, + redshift_density=None, + coordinates_velocity=None, + coordinates_density=None, +): + redshift_dependent_model = eval(f"flip_terms_{model_name}.redshift_dependent_model") + if redshift_dependent_model: + redshift_dict = {} + else: + return None + + if model_type in ["density", "full", "density_velocity"]: + if redshift_dependent_model: + if redshift_density is not None: + redshift_dict["g"] = redshift_density + else: + if len(coordinates_density) < 4: + raise ValueError( + "You are using a model which is redshift dependent." + "Please provide redshifts as the fourth field" + "of the coordinates_density value" + ) + else: + redshift_dict["g"] = coordinates_density[3] + if model_type in ["velocity", "full", "density_velocity"]: + if redshift_dependent_model: + if redshift_velocity is not None: + redshift_dict["v"] = redshift_velocity + else: + if len(coordinates_velocity) < 4: + raise ValueError( + "You are using a model which is redshift dependent." + "Please provide redshifts as the fourth field" + "of the coordinates_velocity value" + ) + else: + redshift_dict["v"] = coordinates_velocity[3] + return redshift_dict + + def generate_covariance( model_name, model_type, @@ -446,6 +488,7 @@ def generate_covariance( coordinates_velocity, ) covariance_dict = {} + if model_type in ["density", "full", "density_velocity"]: covariance_dict["gg"] = compute_cov( model_name, @@ -493,4 +536,11 @@ def generate_covariance( hankel=hankel, los_definition=los_definition, ) - return covariance_dict, number_densities, number_velocities + + redshift_dict = generate_redshift_dict( + model_name, + model_type, + coordinates_velocity=coordinates_velocity, + coordinates_density=coordinates_density, + ) + return covariance_dict, number_densities, number_velocities, redshift_dict diff --git a/flip/covariance/lai22/coefficients.py b/flip/covariance/lai22/coefficients.py index 4f857fb..01d5134 100644 --- a/flip/covariance/lai22/coefficients.py +++ b/flip/covariance/lai22/coefficients.py @@ -7,6 +7,8 @@ def get_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): coefficients_dict = {} if model_type in ["density", "full", "density_velocity"]: diff --git a/flip/covariance/lai22/fisher_terms.py b/flip/covariance/lai22/fisher_terms.py index b57b975..4910fbe 100644 --- a/flip/covariance/lai22/fisher_terms.py +++ b/flip/covariance/lai22/fisher_terms.py @@ -5,6 +5,8 @@ def get_partial_derivative_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): if model_type == "density": return get_partial_derivative_coefficients_density( diff --git a/flip/covariance/lai22/flip_terms.py b/flip/covariance/lai22/flip_terms.py index a6d7bb7..5271062 100644 --- a/flip/covariance/lai22/flip_terms.py +++ b/flip/covariance/lai22/flip_terms.py @@ -15984,3 +15984,4 @@ def N_vv_0_0_2_0(theta, phi): "vv_0_0_2": 1, } multi_index_model = True +redshift_dependent_model = False \ No newline at end of file diff --git a/flip/covariance/lai22/generator.py b/flip/covariance/lai22/generator.py index dba81fa..51fd739 100644 --- a/flip/covariance/lai22/generator.py +++ b/flip/covariance/lai22/generator.py @@ -16,9 +16,7 @@ def compute_correlation_coefficient_simple_integration(p, q, l, r, k, pk): """ " Here the sigma_u is added to pk later. The (2*np.pi**2) is added here in the Lai et al. formalism.""" kr = np.outer(k, r) - integrand = ( - spherical_jn(l, kr).T * k**2 * k ** (2 * (p + q)) * pk / (2 * np.pi**2) - ) + integrand = spherical_jn(l, kr).T * k**2 * k ** (2 * (p + q)) * pk / (2 * np.pi**2) return integrate.simps(integrand, x=k) @@ -834,4 +832,12 @@ def generate_covariance( los_definition=los_definition, **kwargs, ) - return covariance_dict, number_densities, number_velocities, los_definition + + redshift_dict = None + return ( + covariance_dict, + number_densities, + number_velocities, + los_definition, + redshift_dict, + ) diff --git a/flip/covariance/ravouxcarreres/coefficients.py b/flip/covariance/ravouxcarreres/coefficients.py index dfe84ce..cfa8395 100644 --- a/flip/covariance/ravouxcarreres/coefficients.py +++ b/flip/covariance/ravouxcarreres/coefficients.py @@ -2,6 +2,8 @@ def get_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): coefficients_dict = {} if model_type in ["density", "full", "density_velocity"]: diff --git a/flip/covariance/ravouxcarreres/fisher_terms.py b/flip/covariance/ravouxcarreres/fisher_terms.py index b185a0b..110e4c7 100644 --- a/flip/covariance/ravouxcarreres/fisher_terms.py +++ b/flip/covariance/ravouxcarreres/fisher_terms.py @@ -5,6 +5,8 @@ def get_partial_derivative_coefficients( model_type, parameter_values_dict, variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, ): if model_type == "density": return get_partial_derivative_coefficients_density( diff --git a/flip/covariance/ravouxcarreres/flip_terms.py b/flip/covariance/ravouxcarreres/flip_terms.py index c3f5b6c..99418cf 100644 --- a/flip/covariance/ravouxcarreres/flip_terms.py +++ b/flip/covariance/ravouxcarreres/flip_terms.py @@ -883,3 +883,4 @@ def N_vv_0_2_0(theta, phi): "vv_0_2": 1, } multi_index_model = False +redshift_dependent_model = False diff --git a/flip/covariance/rcrk24/coefficients.py b/flip/covariance/rcrk24/coefficients.py new file mode 100644 index 0000000..e501bb6 --- /dev/null +++ b/flip/covariance/rcrk24/coefficients.py @@ -0,0 +1,31 @@ +import numpy as np +from astropy.cosmology import FlatLambdaCDM + + +def get_coefficients( + model_type, + parameter_values_dict, + variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, +): + coefficients_dict = {} + redshift_velocities = redshift_dict["v"] + cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) + + coefficient_vector = ( + np.array(cosmo.Om(redshift_velocities)) ** parameter_values_dict["gamma"] + * cosmo.H(redshift_velocities).value + / cosmo.H0 + * power_spectrum_amplitude_function(redshift_velocities) + / (1 + redshift_velocities) + ) + + coefficients_dict["vv"] = [np.outer(coefficient_vector, coefficient_vector)] + return coefficients_dict + + +def get_diagonal_coefficients(model_type, parameter_values_dict): + coefficients_dict = {} + coefficients_dict["vv"] = parameter_values_dict["sigv"] ** 2 + return coefficients_dict diff --git a/flip/covariance/rcrk24/fisher_terms.py b/flip/covariance/rcrk24/fisher_terms.py new file mode 100644 index 0000000..b41f52a --- /dev/null +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -0,0 +1,209 @@ +import numpy as np +from astropy.cosmology import FlatLambdaCDM + +# The flip convention is to split the power spectrum into several terms +# where linearity assumptions are made +# P_ab = A * B * P0_xy +# +# A is the coefficients where +# P_ab = A P_xy +# B is the cofficient where +# P_xy = B P0_xy +# and +# P0_xy is the power spectrum for a fiducial cosmology at z=0 + +# The derivative cofficients we need are +# dA/dp B + A dB/dp + +# for vv +# A = (aHfs8)_1 * (aHfs8_1) +# B = s8_1 * s8_2 + + +def get_partial_derivative_coefficients( + model_type, + parameter_values_dict, + variant=None, + redshift_dict=None, + power_spectrum_amplitude_function=None, +): + + # s80 is considered to be fixed by the CMB and is hence not a fit parameter + + s80 = 0.832 + redshift_velocities = redshift_dict["v"] + a = 1 / (1 + redshift_velocities) + + cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) + + # The Om0-gamma model f=Omega(Om0)^gamma + + # pre-calculate a few things that are used repeatedly + cosmoOm = np.array(cosmo.Om(redshift_velocities)) + f = cosmoOm ** parameter_values_dict["gamma"] + f0 = parameter_values_dict["Om0"] ** parameter_values_dict["gamma"] + + # Calculation of s8 and its derivatives requires an integral. It is useful to + # expand Omega in terms of (1-a), which allows analytic solutions + + # approximation + def lnD(a): + return np.log(a) * ( + f0 + + f0 + * 3 + * parameter_values_dict["gamma"] + * (1 - parameter_values_dict["Om0"]) + ) + (1 - a) * f0 * 3 * parameter_values_dict["gamma"] * ( + 1 - parameter_values_dict["Om0"] + ) + + def dlnDdOm0(a): + return ( + parameter_values_dict["gamma"] + * parameter_values_dict["Om0"] ** (parameter_values_dict["gamma"] - 1) + * ( + 3 + * (a - 1) + * ( + parameter_values_dict["gamma"] * (parameter_values_dict["Om0"] - 1) + + parameter_values_dict["Om0"] + ) + + np.log(a) + * ( + -3 + * parameter_values_dict["gamma"] + * (parameter_values_dict["Om0"] - 1) + - 3 * parameter_values_dict["Om0"] + + 1 + ) + ) + ) + + def dlnDdgamma(a): + return ( + 3 * (1 - a) * (1 - parameter_values_dict["Om0"]) * f0 + + 3 + * (1 - a) + * parameter_values_dict["gamma"] + * (1 - parameter_values_dict["Om0"]) + * f0 + * np.log(parameter_values_dict["Om0"]) + + np.log(a) + * ( + 3 * (1 - parameter_values_dict["Om0"]) * f0 + + 3 + * parameter_values_dict["gamma"] + * (1 - parameter_values_dict["Om0"]) + * f0 + * np.log(parameter_values_dict["Om0"]) + + f0 * np.log(parameter_values_dict["Om0"]) + ) + ) + + def s8(a): + return s80 * np.exp(lnD(a)) + + def aHfs8(a): + return ( + f + * s8(a) + / (1 + redshift_velocities) + * cosmo.H(redshift_velocities) + / cosmo.H0 + ) + + # now for the partials + + def dOmdOm0(a): + numerator = parameter_values_dict["Om0"] * a ** (-3) + denominator = numerator + 1 - parameter_values_dict["Om0"] + return a ** (-3) / denominator - numerator / denominator**2 * (a ** (-3) - 1) + + def dfdOm0(a): + return parameter_values_dict["gamma"] * f / cosmoOm * dOmdOm0(a) + + def dfdgamma(a): + return np.log(cosmoOm) * f + + def ds8dOm0(a): + return s8(a) * dlnDdOm0(a) + + def ds8dgamma(a): + return s8(a) * dlnDdgamma(a) + + def dAdOm0(a): + return a * cosmo.H(a) / cosmo.H0 * (dfdOm0(a) * s8(a) + f * ds8dOm0(a)) + + def dAdgamma(a): + return a * cosmo.H(a) / cosmo.H0 * (dfdgamma(a) * s8(a) + f * ds8dgamma(a)) + + aHfs8s8 = aHfs8(a) * s8(a) + + Omega_m_partial_derivative_coefficients = dAdOm0(a) * s8(a) + aHfs8(a) * ds8dOm0(a) + + gamma_partial_derivative_coefficients = dAdgamma(a) * s8(a) + aHfs8(a) * ds8dgamma( + a + ) + + # in the fs8 case + def s8_fs8(a): + return s80 + parameter_values_dict["fs8"] * np.log(a) + + def ds8dfs8(a): + return np.log(a) + + fs8_partial_derivative_coefficients = ( + a + * cosmo.H(redshift_velocities) + / cosmo.H0 + * (s8_fs8(a) + parameter_values_dict["fs8"] * ds8dfs8(a)) + ) + + aHfs8s8_fs8 = ( + a + * cosmo.H(redshift_velocities) + / cosmo.H0 + * parameter_values_dict["fs8"] + * s8_fs8(a) + ) + partial_coefficients_dict = { + "Omegam": { + "vv": [ + np.outer( + Omega_m_partial_derivative_coefficients, + aHfs8s8, + ) + + np.outer( + aHfs8s8, + Omega_m_partial_derivative_coefficients, + ), + ], + }, + "gamma": { + "vv": [ + np.outer( + gamma_partial_derivative_coefficients, + aHfs8s8, + ) + + np.outer( + aHfs8s8, + gamma_partial_derivative_coefficients, + ), + ], + }, + "fs8": { + "vv": [ + np.outer( + fs8_partial_derivative_coefficients, + aHfs8s8_fs8, + ) + + np.outer( + aHfs8s8_fs8, + fs8_partial_derivative_coefficients, + ), + ], + }, + } + + return partial_coefficients_dict diff --git a/flip/covariance/rcrk24/flip_terms.py b/flip/covariance/rcrk24/flip_terms.py new file mode 100644 index 0000000..b2b7a8e --- /dev/null +++ b/flip/covariance/rcrk24/flip_terms.py @@ -0,0 +1,33 @@ +import numpy as np +import scipy + + +def M_vv_0_0_0(): + def func(k): + return (10000 / 9) / k**2 + + return func + + +def N_vv_0_0_0(theta, phi): + return 3 * np.cos(theta) + + +def M_vv_0_2_0(): + def func(k): + return (10000 / 9) / k**2 + + return func + + +def N_vv_0_2_0(theta, phi): + return (9 / 2) * np.cos(2 * phi) + (3 / 2) * np.cos(theta) + +def power_spectrum_amplitude_function(r): + return 1. + +dictionary_terms = {"vv": ["0"]} +dictionary_lmax = {"vv": [2]} +dictionary_subterms = {"vv_0_0": 1, "vv_0_1": 0, "vv_0_2": 1} +multi_index_model = False +redshift_dependent_model = True diff --git a/flip/covariance/symbolic.py b/flip/covariance/symbolic.py index 74daada..b539164 100644 --- a/flip/covariance/symbolic.py +++ b/flip/covariance/symbolic.py @@ -160,6 +160,7 @@ def write_output( l1max_list=None, l2max_list=None, multi_index_model=False, + redshift_dependent_model=False, ): """ The write_output function takes the following arguments: @@ -259,6 +260,8 @@ def write_output( f.write(f"multi_index_model = {multi_index_model}") f.write("\n") + f.write(f"redshift_dependent_model = {redshift_dependent_model}") + f.write("\n") f.close() @@ -274,6 +277,7 @@ def write_M_N_functions( l1max_list=None, l2max_list=None, multi_index_model=False, + redshift_dependent_model=False, ): """ The write_M_N_functions function is used to generate the M_N functions for a given model. @@ -350,6 +354,7 @@ def write_M_N_functions( l1max_list=l1max_list, l2max_list=l2max_list, multi_index_model=multi_index_model, + redshift_dependent_model=redshift_dependent_model, ) @@ -637,6 +642,43 @@ def generate_generalized_ravouxcarreres_functions( ) +def generate_generalized_rcrk24_functions( + filename="./rcrk24/flip_terms.py", number_worker=8 +): + """ + The generate_generalized_rcrk24_functions function generates the flip_terms.py file in the carreres23 directory, which contains functions that calculate M and N terms for a generalized version of Carreres' (2012) model 2 and 3. + + Args: + filename: Specify the name of the file that will be generated + number_worker: Determine the number of processes to use for multiprocessing + + Returns: + A list of functions, + + """ + mu1, mu2 = sy.symbols("mu1 mu2") + k = sy.symbols("k", positive=True, finite=True, real=True) + type_list = ["vv"] + term_index_list = [["0"]] + lmax_list = [[2]] + l1max_list = [[1]] + l2max_list = [[1]] + dict_B = {"B_vv_0": 100**2 * mu1 * mu2 / k**2} + + write_M_N_functions( + filename, + type_list, + term_index_list, + lmax_list, + dict_B, + number_worker=number_worker, + wide_angle=True, + l1max_list=l1max_list, + l2max_list=l2max_list, + redshift_dependent_model=True, + ) + + def compute_partial_derivative_dictionnary( name_models, components, @@ -710,7 +752,7 @@ def write_partial_derivatives( ) f.write( - "def get_partial_derivative_coefficients(model_type,parameter_values_dict,variant=None,):\n" + "def get_partial_derivative_coefficients(model_type,parameter_values_dict,variant=None,redshift_dict=None,power_spectrum_amplitude_function=None,):\n" ) write_one_function( f, @@ -1104,6 +1146,7 @@ def generate_files(): generate_generalized_adamsblake20_functions() generate_generalized_lai22_functions() generate_generalized_ravouxcarreres_functions() + generate_generalized_rcrk24_functions() def generate_fisher_files(): diff --git a/flip/fisher.py b/flip/fisher.py index 115e0fe..dafd2e3 100644 --- a/flip/fisher.py +++ b/flip/fisher.py @@ -175,6 +175,8 @@ def compute_fisher_matrix( self.covariance.model_type, parameter_values_dict, variant=variant, + redshift_dict=self.covariance.redshift_dict, + power_spectrum_amplitude_function=self.covariance.power_spectrum_amplitude_function, ) parameter_name_list = [] covariance_derivative_sum_list = [] diff --git a/scripts/develop.py b/scripts/develop.py new file mode 100644 index 0000000..1de4cc6 --- /dev/null +++ b/scripts/develop.py @@ -0,0 +1,132 @@ +import os + +import numpy as np +import pandas as pd +from pkg_resources import resource_filename + +from flip import fisher, utils +from flip.covariance import covariance + +def main(): + flip_base = resource_filename("flip", ".") + data_path = os.path.join(flip_base, "data") + + ### Load data + sn_data = pd.read_parquet(os.path.join(data_path, "velocity_data.parquet")) + + sn_data = sn_data[np.array(sn_data["status"]) != False] + sn_data = sn_data[np.array(sn_data["status"]) != None] + coordinates_velocity = np.array([sn_data["ra"], sn_data["dec"], sn_data["rcom_zobs"], sn_data["zobs"]]) + + data_velocity = sn_data.to_dict("list") + for key in data_velocity.keys(): + data_velocity[key] = np.array(data_velocity[key]) + data_velocity["velocity"] = data_velocity.pop("vpec") + data_velocity["velocity_error"] = np.zeros_like(data_velocity["velocity"]) + + + ktt, ptt = np.loadtxt(os.path.join(data_path, "power_spectrum_tt.txt")) + kmt, pmt = np.loadtxt(os.path.join(data_path, "power_spectrum_mt.txt")) + kmm, pmm = np.loadtxt(os.path.join(data_path, "power_spectrum_mm.txt")) + + sigmau_fiducial = 15 + + power_spectrum_dict = {"vv": [[ktt, ptt * utils.Du(ktt, sigmau_fiducial) ** 2]]} + + ### Compute covariance + size_batch = 10_000 + number_worker = 16 + + + from flip.covariance.rcrk24.flip_terms import power_spectrum_amplitude_function + covariance_fit = covariance.CovMatrix.init_from_flip( + "rcrk24", + # "agk24" + # 'carreres23', + "velocity", + power_spectrum_dict, + coordinates_velocity=coordinates_velocity, + size_batch=size_batch, + number_worker=number_worker, + power_spectrum_amplitude_function=power_spectrum_amplitude_function, + ) + + ### Load fitter + + fisher_properties = { + "inversion_method": "inverse", + "velocity_type": "scatter", + } + + variant = None # can be replaced by growth_index + + parameter_dict = { + "fs8": 0.4, + "Om0": 0.3, + "gamma": 0.55, + "sigv": 200, + "sigma_M": 0.12, + } + + # parameter_dict = { + # "Om0": 0.3, + # "gamma": 0.55, + # "sigv": 200, + # "sigma_M": 0.12, + # } + + Fisher = fisher.FisherMatrix.init_from_covariance( + covariance_fit, + data_velocity, + parameter_dict, + fisher_properties=fisher_properties, + ) + + parameter_name_list, fisher_matrix = Fisher.compute_fisher_matrix( + parameter_dict, variant=variant + ) + return parameter_name_list, fisher_matrix + +def dlnDdOm0(parameter_values_dict): + a=(1/(1+0.05)) + lna=np.log(a) + return ( + parameter_values_dict["gamma"] * parameter_values_dict["Om0"]**(parameter_values_dict["gamma"]-1) * + ( + 3 * parameter_values_dict["gamma"] * (parameter_values_dict["Om0"]-1) * (a - lna -1) + + 3 * (a-1) * parameter_values_dict["Om0"] - + 3 * np.log(a) * parameter_values_dict["Om0"] + lna + ) + ) + +def dlnDdgamma(parameter_values_dict): + a=(1/(1+0.05)) + lna=np.log(a) + f0 = parameter_values_dict["Om0"]**parameter_values_dict["gamma"] + return ( + f0 * + ( + np.log(parameter_values_dict["Om0"]) * + ( + 3 * parameter_values_dict["gamma"] * (parameter_values_dict["Om0"]-1) * (a - lna -1) + lna + ) + + 3 * (parameter_values_dict["Om0"]-1) * (a - lna -1) + ) + ) + +if __name__ == "__main__": + parameter_dict = { + "Om0": 0.3, + "gamma": 0.55, + "sigv": 200, + "sigma_M": 0.12, + } + parameter_name_list, fisher_matrix = main() + cov = np.linalg.inv(fisher_matrix[0:2,0:2]) + + s80 = 0.832 + partials = s80*np.array([parameter_dict['gamma']*parameter_dict['Om0']**(parameter_dict['gamma']-1),np.log(parameter_dict['Om0'])*parameter_dict['Om0']**parameter_dict['gamma']]) + partials = partials + parameter_dict['Om0']**parameter_dict['gamma'] *s80 * np.array([dlnDdOm0(parameter_dict), dlnDdgamma(parameter_dict)]) + + print(np.sqrt(partials.T @ cov[0:2,0:2] @ partials)) + print(1/np.sqrt(fisher_matrix[2,2])) diff --git a/scripts/flip_compute_power_spectra.py b/scripts/flip_compute_power_spectra.py index cb206d6..1404026 100644 --- a/scripts/flip_compute_power_spectra.py +++ b/scripts/flip_compute_power_spectra.py @@ -15,11 +15,16 @@ maximal_wavenumber = 1.000 number_points = 1000 logspace = True -redshift = 0.11 -normalize_power_spectrum = True # Normalize by sigma_8 -power_spectrum_non_linear_model = ( - "halofit" # If you want to add non linearity, what model to use -) +redshift = 0.0 + +# Can be changed by "growth_rate" (fsigma_8^fid normalization) or "no_normalization" +# Here, it normalize by sigma_8^fid +normalization_power_spectrum = "growth_amplitude" + +# If you want to add non linearity, what model to use +power_spectrum_non_linear_model = "halofit" + +# Model used to compute matter and velocity divergence power spectra power_spectrum_model = "linearbel" save_path = "./" # If not None, will save all calculated power spectra in this folder @@ -37,7 +42,7 @@ maximal_wavenumber, number_points, logspace=logspace, - normalize_power_spectrum=normalize_power_spectrum, + normalization_power_spectrum=normalization_power_spectrum, power_spectrum_non_linear_model=power_spectrum_non_linear_model, power_spectrum_model=power_spectrum_model, save_path=save_path,