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/generator.py b/flip/covariance/generator.py index 7b71626..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", ] diff --git a/flip/covariance/rcrk24/coefficients.py b/flip/covariance/rcrk24/coefficients.py index dc318e6..e501bb6 100644 --- a/flip/covariance/rcrk24/coefficients.py +++ b/flip/covariance/rcrk24/coefficients.py @@ -14,7 +14,7 @@ def get_coefficients( cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) coefficient_vector = ( - cosmo.Om(redshift_velocities).value ** parameter_values_dict["gamma"] + np.array(cosmo.Om(redshift_velocities)) ** parameter_values_dict["gamma"] * cosmo.H(redshift_velocities).value / cosmo.H0 * power_spectrum_amplitude_function(redshift_velocities) diff --git a/flip/covariance/rcrk24/fisher_terms.py b/flip/covariance/rcrk24/fisher_terms.py index 07be96e..50456c2 100644 --- a/flip/covariance/rcrk24/fisher_terms.py +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -1,6 +1,23 @@ 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, @@ -10,35 +27,112 @@ def get_partial_derivative_coefficients( power_spectrum_amplitude_function=None, ): - cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) + # 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 = ( - 2 - * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"]) - * parameter_values_dict["gamma"] - * power_spectrum_amplitude_function(redshift_velocities) ** 2 - / cosmo.Om(redshift_velocities).value + dAdOm0(a) * s8(a) + aHfs8(a) * ds8dOm0(a) ) gamma_partial_derivative_coefficients = ( - 2 - * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"]) - * power_spectrum_amplitude_function(redshift_velocities) ** 2 - * np.log(cosmo.Om(redshift_velocities).value) + dAdgamma(a) * s8(a) + aHfs8(a) * ds8dgamma(a) ) - s8_partial_derivative_coefficients = ( - 2 - * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"]) - * power_spectrum_amplitude_function(redshift_velocities) + # 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, ), ], @@ -47,15 +141,23 @@ def get_partial_derivative_coefficients( "vv": [ np.outer( gamma_partial_derivative_coefficients, + aHfs8s8, + ) + + np.outer( + aHfs8s8, gamma_partial_derivative_coefficients, ), ], }, - "s8": { + "fs8": { "vv": [ np.outer( - s8_partial_derivative_coefficients, - s8_partial_derivative_coefficients, + fs8_partial_derivative_coefficients, + aHfs8s8_fs8, + ) + + np.outer( + aHfs8s8_fs8, + fs8_partial_derivative_coefficients, ), ], }, diff --git a/flip/covariance/rcrk24/flip_terms.py b/flip/covariance/rcrk24/flip_terms.py index 7f44710..b2b7a8e 100644 --- a/flip/covariance/rcrk24/flip_terms.py +++ b/flip/covariance/rcrk24/flip_terms.py @@ -23,6 +23,8 @@ def func(k): 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]} 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]))