From 8fedb5c8aa3a76a63dd234e9a666d61a71a554e0 Mon Sep 17 00:00:00 2001 From: Alex Kim Date: Thu, 27 Jun 2024 14:07:46 -0700 Subject: [PATCH 01/11] first implementation with redshift-dependent cofactors --- flip/covariance/agkim24/coefficients.py | 17 ++++ flip/covariance/agkim24/fisher_terms.py | 62 +++++++++++ flip/covariance/agkim24/flip_terms.py | 30 ++++++ flip/covariance/agkim24/generator.py | 130 ++++++++++++++++++++++++ flip/covariance/covariance.py | 3 + flip/covariance/generator.py | 2 + flip/fisher.py | 23 ++++- scripts/develop.py | 82 +++++++++++++++ 8 files changed, 344 insertions(+), 5 deletions(-) create mode 100644 flip/covariance/agkim24/coefficients.py create mode 100644 flip/covariance/agkim24/fisher_terms.py create mode 100644 flip/covariance/agkim24/flip_terms.py create mode 100644 flip/covariance/agkim24/generator.py create mode 100644 scripts/develop.py diff --git a/flip/covariance/agkim24/coefficients.py b/flip/covariance/agkim24/coefficients.py new file mode 100644 index 0000000..7436fab --- /dev/null +++ b/flip/covariance/agkim24/coefficients.py @@ -0,0 +1,17 @@ +import numpy as np + + +def get_coefficients( + model_type, + parameter_values_dict, + variant=None, +): + coefficients_dict = {} + coefficients_dict["vv"] = [parameter_values_dict["fs8"] ** 2] + 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/agkim24/fisher_terms.py b/flip/covariance/agkim24/fisher_terms.py new file mode 100644 index 0000000..e1ab1d9 --- /dev/null +++ b/flip/covariance/agkim24/fisher_terms.py @@ -0,0 +1,62 @@ +import numpy as np +import astropy.cosmology +from astropy import units as u +from astropy import constants as const + +class FisherTerms(object): + + def __init__(self, coordinates_velocity, cosmo=astropy.cosmology.Planck18): + self._coordinates_velocity = coordinates_velocity + self._cosmo = cosmo + + z = astropy.cosmology.z_at_value(self._cosmo.comoving_distance, self._coordinates_velocity[2]*u.Mpc) + self._cov_factors = {"v": self._cosmo.H(z)/(1+z) / self._cosmo.H0, "g": None} # need to implement "g"" + + @property + def cov_factors(self): + return self._cov_factors + + def get_partial_derivative_coefficients(self, + model_type, + parameter_values_dict, + variant=None, + ): + if variant == "growth_index": + partial_coefficients_dict = { + "Omegam": { + "vv": [ + 2 + * parameter_values_dict["Omegam"] + ** (2 * parameter_values_dict["gamma"]) + * parameter_values_dict["gamma"] + * parameter_values_dict["s8"] ** 2 + / parameter_values_dict["Omegam"], + ], + }, + "gamma": { + "vv": [ + 2 + * parameter_values_dict["Omegam"] + ** (2 * parameter_values_dict["gamma"]) + * parameter_values_dict["s8"] ** 2 + * np.log(parameter_values_dict["Omegam"]), + ], + }, + "s8": { + "vv": [ + 2 + * parameter_values_dict["Omegam"] + ** (2 * parameter_values_dict["gamma"]) + * parameter_values_dict["s8"], + ], + }, + } + else: + partial_coefficients_dict = { + "fs8": { + "vv": [ + 2 * parameter_values_dict["fs8"], + ], + }, + } + return partial_coefficients_dict diff --git a/flip/covariance/agkim24/flip_terms.py b/flip/covariance/agkim24/flip_terms.py new file mode 100644 index 0000000..f00c9d5 --- /dev/null +++ b/flip/covariance/agkim24/flip_terms.py @@ -0,0 +1,30 @@ +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) + + +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 diff --git a/flip/covariance/agkim24/generator.py b/flip/covariance/agkim24/generator.py new file mode 100644 index 0000000..997cf57 --- /dev/null +++ b/flip/covariance/agkim24/generator.py @@ -0,0 +1,130 @@ +import multiprocessing as mp +from functools import partial + +import numpy as np +from scipy.special import spherical_jn + +from flip.covariance import cov_utils + + +def angle_between(ra_0, ra_1, dec_0, dec_1): + """Compute cos of the angle between r0 and r1.""" + cos_alpha = np.cos(ra_1 - ra_0) * np.cos(dec_0) * np.cos(dec_1) + np.sin( + dec_0 + ) * np.sin(dec_1) + return cos_alpha + + +def separation(r_0, r_1, cos_alpha): + """Compute separation between r_0 and r_1.""" + return np.sqrt(r_0**2 + r_1**2 - 2 * r_0 * r_1 * cos_alpha) + + +def window(r_0, r_1, cos_alpha, sep, j0kr, j2kr): + """Note: here, the bisector angle definition is used in wide-angle""" + win = 1 / 3 * (j0kr - 2 * j2kr) * cos_alpha + win += j2kr * r_0 * r_1 / sep**2 * (1 - cos_alpha**2) + return win + + +def intp(win, k, pk): + pint = win.T * pk + return np.trapz(pint, x=k) + + +def covariance_vv( + ra_in, + dec_in, + rcomov_in, + k_in, + pk_in, + grid_window_in=None, + size_batch=100_000, + number_worker=8, +): + N = len(ra_in) + + if grid_window_in is not None: + pk = pk_in * grid_window_in**2 + else: + pk = pk_in + + n_task = int((N * (N + 1)) / 2) - N + + batches = [] + for n in range(0, n_task, size_batch): + brange = np.arange(n, np.min((n + size_batch, n_task))) + i_list, j_list = cov_utils.compute_i_j(N, brange) + r_comovi, rai, deci = rcomov_in[i_list], ra_in[i_list], dec_in[i_list] + r_comovj, raj, decj = rcomov_in[j_list], ra_in[j_list], dec_in[j_list] + batches.append([rai, raj, deci, decj, r_comovi, r_comovj]) + + with mp.Pool(number_worker) as pool: + func = partial(compute_coef, k_in, pk) + pool_results = pool.map(func, batches) + values = np.concatenate(pool_results) + + var_val = np.trapz(pk / 3, x=k_in) + cov = np.insert(values, 0, var_val) + cov = 100**2 / (2 * np.pi**2) * cov + return cov + + +def compute_coef(k, pk, coord): + cos = angle_between(coord[0], coord[1], coord[2], coord[3]) + sep = separation(coord[4], coord[5], cos) + ksep = np.outer(k, sep) + j0 = spherical_jn(0, ksep) + j2 = spherical_jn(2, ksep) + res = window(coord[4], coord[5], cos, sep, j0, j2) + res = intp(res, k, pk) + return res + + +def generate_covariance( + model_type, + power_spectrum_dict, + coordinates_density=False, + coordinates_velocity=None, + **kwargs, +): + """ + The generate_covariance function generates a covariance matrix for the velocity field. + + Args: + model_type: Specify the type of model to generate + power_spectrum_dict: Pass the power spectrum to the function + coordinates_density: Specify the coordinates of the density field + coordinates_velocity: Generate the covariance matrix + **kwargs: Pass additional parameters to the function + : Generate the covariance matrix for the velocity field + The wide angle used is the bisector. + Returns: + A dictionary with a single key "vv" + + """ + assert model_type == "velocity" + cov_utils.check_generator_need( + model_type, + coordinates_density, + coordinates_velocity, + ) + number_densities = None + number_velocities = len(coordinates_velocity[0]) + cov_vv = covariance_vv( + coordinates_velocity[0], + coordinates_velocity[1], + coordinates_velocity[2], + power_spectrum_dict["vv"][0][0], + power_spectrum_dict["vv"][0][1], + **kwargs, + ) + + los_definition = "bisector" + + return ( + {"vv": np.array([cov_vv])}, + number_densities, + number_velocities, + los_definition, + ) diff --git a/flip/covariance/covariance.py b/flip/covariance/covariance.py index 44384c2..c5283f6 100644 --- a/flip/covariance/covariance.py +++ b/flip/covariance/covariance.py @@ -21,6 +21,7 @@ def __init__( number_densities=None, number_velocities=None, variant=None, + coordinates_velocity=None, ): """ The __init__ function is called when the class is instantiated. @@ -50,6 +51,7 @@ def __init__( self.number_densities = number_densities self.number_velocities = number_velocities self.variant = variant + self.coordinates_velocity = coordinates_velocity @classmethod def init_from_flip( @@ -116,6 +118,7 @@ def init_from_flip( number_densities=number_densities, number_velocities=number_velocities, variant=variant, + coordinates_velocity=coordinates_velocity, ) @classmethod diff --git a/flip/covariance/generator.py b/flip/covariance/generator.py index 83554ba..62e587d 100644 --- a/flip/covariance/generator.py +++ b/flip/covariance/generator.py @@ -13,6 +13,7 @@ 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.ravouxcarreres import flip_terms as flip_terms_ravouxcarreres +from flip.covariance.agkim24 import flip_terms as flip_terms_agkim24 from flip.utils import create_log log = create_log() @@ -22,6 +23,7 @@ "lai22", "carreres23", "ravouxcarreres", + "agkim24", ] diff --git a/flip/fisher.py b/flip/fisher.py index 115e0fe..d81236f 100644 --- a/flip/fisher.py +++ b/flip/fisher.py @@ -97,7 +97,7 @@ def load_error_vector( def compute_covariance_derivative( self, - partial_coefficients_dict_param, + partial_coefficients_dict_param, cov_factors = None, ): if self.covariance.model_type == "density": @@ -110,9 +110,14 @@ def compute_covariance_derivative( ) elif self.covariance.model_type == "velocity": + if cov_factors is None: + cov_prefactor = 1. + else: + cov_prefactor = np.outer(cov_factors["v"],cov_factors["v"]) + covariance_derivative_sum = np.sum( [ - partial_coefficients_dict_param["vv"][i] * cov + partial_coefficients_dict_param["vv"][i] * cov_prefactor * cov for i, cov in enumerate(self.covariance.covariance_dict["vv"]) ], axis=0, @@ -168,9 +173,16 @@ def compute_fisher_matrix( variant=None, ): - coefficients = importlib.import_module( - f"flip.covariance.{self.covariance.model_name}.fisher_terms" - ) + if self.covariance.model_name == "agkim24": + FisherTerms = getattr(importlib.import_module( + f"flip.covariance.{self.covariance.model_name}.fisher_terms"),"FisherTerms") + coefficients = FisherTerms(self.covariance.coordinates_velocity) + cov_factors = coefficients.cov_factors + else: + coefficients = importlib.import_module( + f"flip.covariance.{self.covariance.model_name}.fisher_terms") + cov_factors = None + partial_coefficients_dict = coefficients.get_partial_derivative_coefficients( self.covariance.model_type, parameter_values_dict, @@ -189,6 +201,7 @@ def compute_fisher_matrix( self.inverse_covariance_sum, self.compute_covariance_derivative( partial_coefficients_dict_param, + cov_factors=cov_factors ), ) ) diff --git a/scripts/develop.py b/scripts/develop.py new file mode 100644 index 0000000..686c60a --- /dev/null +++ b/scripts/develop.py @@ -0,0 +1,82 @@ +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"]]) + + 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 + + covariance_fit = covariance.CovMatrix.init_from_flip( + "agkim24", + # 'carreres23', + "velocity", + power_spectrum_dict, + coordinates_velocity=coordinates_velocity, + size_batch=size_batch, + number_worker=number_worker, + ) + + + ### Load fitter + + fisher_properties = { + "inversion_method": "inverse", + "velocity_type": "scatter", + } + + variant = None # can be replaced by growth_index + + parameter_dict = { + "fs8": 0.4, + "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 + + +if __name__ == "__main__": + parameter_name_list, fisher_matrix = main() + print(fisher_matrix) \ No newline at end of file From 5517adece2fb893b9cc96bdbca27e5126b7538a9 Mon Sep 17 00:00:00 2001 From: Alex Date: Fri, 28 Jun 2024 09:00:39 -0700 Subject: [PATCH 02/11] whatever --- notebook/covariance.ipynb | 181 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 181 insertions(+) create mode 100644 notebook/covariance.ipynb diff --git a/notebook/covariance.ipynb b/notebook/covariance.ipynb new file mode 100644 index 0000000..0ac7b7f --- /dev/null +++ b/notebook/covariance.ipynb @@ -0,0 +1,181 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "71e15407-b85a-4152-858f-86087f5fae7d", + "metadata": {}, + "source": [ + "This is a tutorial to use the flip package: https://github.com/corentinravoux/flip \\\n", + "It is self-contained and can be used in google collab or on your environement \\\n", + "All the data used are subsampled version of a simulation. \\\n", + "The data size is small for the tutorial, do not use it for science case. \\" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ae55978f-ac85-4f46-80cf-112204329a68", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "%%capture\n", + "!pip install git+https://github.com/corentinravoux/flip" + ] + }, + { + "cell_type": "markdown", + "id": "ef379b12-2598-413c-9f05-bb7525b4e96a", + "metadata": {}, + "source": [ + "Loading of the necessary modules" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e1da453-ef47-4920-b330-a7d9289e2f01", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "import os\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import pandas as pd\n", + "from flip import fitter, plot_utils, utils, vectors\n", + "from flip.covariance import covariance, contraction\n", + "from pkg_resources import resource_filename\n", + "flip_base = resource_filename(\"flip\", \".\")\n", + "data_path = os.path.join(flip_base, \"data\")\n", + "plt.style.use(os.path.join(data_path,\"style.mplstyle\"))" + ] + }, + { + "cell_type": "markdown", + "id": "7d184b4d-4672-4bf7-8f21-c156e0e195f7", + "metadata": {}, + "source": [ + "Loading the data, located in the package itself\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f81a2bc9-1d9e-4969-ac2d-0b5e56fd189c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "ktt, ptt = np.loadtxt(os.path.join(data_path, \"power_spectrum_tt.txt\"))\n", + "kmt, pmt = np.loadtxt(os.path.join(data_path, \"power_spectrum_mt.txt\"))\n", + "kmm, pmm = np.loadtxt(os.path.join(data_path, \"power_spectrum_mm.txt\"))\n", + "\n", + "# the extra Koda et al. damping term Eq 33 from https://arxiv.org/pdf/2303.01198\n", + "sigmau_fiducial = 15\n", + "\n", + "power_spectrum_dict = {\"vv\": [[ktt, ptt * utils.Du(ktt, sigmau_fiducial) ** 2]]}\n" + ] + }, + { + "cell_type": "markdown", + "id": "a561c106-87a9-4e42-aaca-9d103e6e750d", + "metadata": {}, + "source": [ + "First part of the flip package: fast computation of the covariance matrix based on theory and SNIa coordinates" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96947c01-98ed-4793-8fc9-9a19f0830486", + "metadata": {}, + "outputs": [], + "source": [ + "filename = \"abacussample_s500_r300_b20_d0.05_i0.1_xyz000_lnone_go.hdf5\"\n", + "import h5py\n", + "import numpy\n", + "from astropy.cosmology import Planck18 as cosmo\n", + "\n", + "def read_abacus(filename):\n", + " with h5py.File(filename, \"r\") as f:\n", + " data=f['abacussample']\n", + " zlist = list(data['z']) # observed redshift\n", + " ra = list(data['ra'])\n", + " dec = list(data['dec'])\n", + " rcom_zobs = cosmo.comoving_distance(numpy.array(zlist)).value\n", + " # cosmoz = list(data['cosmoz'])\n", + " \n", + " return numpy.array([ra,dec,rcom_zobs])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8dfdb97c-8720-46b8-aa48-43fb1b6a4f83", + "metadata": {}, + "outputs": [], + "source": [ + "coordinates_velocity = read_abacus(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a30c697f-6597-450d-ac8b-6941e9ab8219", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "size_batch=10_000\n", + "number_worker=16\n", + "\n", + "\n", + "covariance_fit = covariance.CovMatrix.init_from_flip(\n", + " \"carreres23\",\n", + " \"velocity\",\n", + " power_spectrum_dict,\n", + " coordinates_velocity=coordinates_velocity,\n", + " size_batch=size_batch,\n", + " number_worker=number_worker,\n", + ")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd338dc8-69c9-40ea-ae8b-df4eb4eda84a", + "metadata": {}, + "outputs": [], + "source": [ + "covariance_fit.write(\"temp\",\"pickle\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From cf84d6392e277110b467937edbf9dcf9bd81ca41 Mon Sep 17 00:00:00 2001 From: Alex Kim Date: Mon, 22 Jul 2024 11:41:02 -0700 Subject: [PATCH 03/11] Om-gamma parameterization gives a reasonable answer. Is it correct???? --- flip/covariance/generator.py | 2 ++ flip/covariance/rcrk24/coefficients.py | 4 +++- flip/covariance/rcrk24/fisher_terms.py | 33 ++++++++++++++++++++++---- flip/covariance/rcrk24/flip_terms.py | 2 ++ scripts/develop.py | 31 ++++++++++++++++++++---- 5 files changed, 61 insertions(+), 11 deletions(-) diff --git a/flip/covariance/generator.py b/flip/covariance/generator.py index e1d0c75..f44591e 100644 --- a/flip/covariance/generator.py +++ b/flip/covariance/generator.py @@ -14,6 +14,7 @@ from flip.covariance.lai22 import flip_terms as flip_terms_lai22 from flip.covariance.ravouxcarreres import flip_terms as flip_terms_ravouxcarreres from flip.covariance.agkim24 import flip_terms as flip_terms_agkim24 +from flip.covariance.rcrk24 import flip_terms as flip_terms_rcrk24 from flip.utils import create_log log = create_log() @@ -24,6 +25,7 @@ "carreres23", "ravouxcarreres", "agkim24", + "rcrk24", ] diff --git a/flip/covariance/rcrk24/coefficients.py b/flip/covariance/rcrk24/coefficients.py index dc318e6..f7ecb60 100644 --- a/flip/covariance/rcrk24/coefficients.py +++ b/flip/covariance/rcrk24/coefficients.py @@ -13,8 +13,10 @@ def get_coefficients( redshift_velocities = redshift_dict["v"] cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) + # print(power_spectrum_amplitude_function) + # qwd 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..adcbadc 100644 --- a/flip/covariance/rcrk24/fisher_terms.py +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -13,24 +13,39 @@ def get_partial_derivative_coefficients( cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) redshift_velocities = redshift_dict["v"] + cosmoOm=np.array(cosmo.Om(redshift_velocities)) + Omega_m_partial_derivative_coefficients = ( 2 - * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"]) + * np.array(cosmo.Om(redshift_velocities)) ** (2 * parameter_values_dict["gamma"]) * parameter_values_dict["gamma"] * power_spectrum_amplitude_function(redshift_velocities) ** 2 - / cosmo.Om(redshift_velocities).value + / np.array(cosmo.Om(redshift_velocities)) ) gamma_partial_derivative_coefficients = ( 2 - * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"]) + * np.array(cosmo.Om(redshift_velocities)) ** (2 * parameter_values_dict["gamma"]) * power_spectrum_amplitude_function(redshift_velocities) ** 2 - * np.log(cosmo.Om(redshift_velocities).value) + * np.log(cosmo.Om(redshift_velocities)) ) + Omega_m_partial_derivative_coefficients = ( + cosmoOm ** (parameter_values_dict["gamma"]-1) + * parameter_values_dict["gamma"] + * power_spectrum_amplitude_function(redshift_velocities) + ) + + gamma_partial_derivative_coefficients = ( + cosmoOm ** (parameter_values_dict["gamma"]) + * power_spectrum_amplitude_function(redshift_velocities) + * np.log(cosmoOm) + ) + + s8_partial_derivative_coefficients = ( 2 - * cosmo.Om(redshift_velocities).value ** (2 * parameter_values_dict["gamma"]) + * cosmoOm ** (2 * parameter_values_dict["gamma"]) * power_spectrum_amplitude_function(redshift_velocities) ) @@ -39,6 +54,10 @@ def get_partial_derivative_coefficients( "vv": [ np.outer( Omega_m_partial_derivative_coefficients, + cosmoOm ** (parameter_values_dict["gamma"]), + ) + + np.outer( + cosmoOm ** (parameter_values_dict["gamma"]), Omega_m_partial_derivative_coefficients, ), ], @@ -47,6 +66,10 @@ def get_partial_derivative_coefficients( "vv": [ np.outer( gamma_partial_derivative_coefficients, + cosmoOm ** (parameter_values_dict["gamma"]), + ) + + np.outer( + cosmoOm ** (parameter_values_dict["gamma"]), gamma_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 index 686c60a..7b9fb90 100644 --- a/scripts/develop.py +++ b/scripts/develop.py @@ -16,8 +16,7 @@ def main(): 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"]]) + 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(): @@ -38,17 +37,20 @@ def main(): 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( - "agkim24", + "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 = { @@ -64,6 +66,13 @@ def main(): "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, @@ -78,5 +87,17 @@ def main(): if __name__ == "__main__": + parameter_dict = { + "Om0": 0.3, + "gamma": 0.55, + "sigv": 200, + "sigma_M": 0.12, + } parameter_name_list, fisher_matrix = main() - print(fisher_matrix) \ No newline at end of file + print(fisher_matrix) + cov = np.linalg.inv(fisher_matrix[0:2,0:2]) + print(cov[0:2,0:2]) + + partials = np.array([parameter_dict['gamma']*parameter_dict['Om0']**(parameter_dict['gamma']-1),np.log(parameter_dict['Om0'])*parameter_dict['Om0']**parameter_dict['gamma']]) + print(partials.T @ cov[0:2,0:2] @ partials) + print(np.sqrt(partials.T @ cov[0:2,0:2] @ partials)) \ No newline at end of file From d770b2e3f5d26666cdff1b955ce704071d206456 Mon Sep 17 00:00:00 2001 From: AlexGKim Date: Mon, 22 Jul 2024 11:55:24 -0700 Subject: [PATCH 04/11] Delete flip/covariance/agkim24 directory getting in sync with upstream --- flip/covariance/agkim24/coefficients.py | 17 ---- flip/covariance/agkim24/fisher_terms.py | 62 ----------- flip/covariance/agkim24/flip_terms.py | 30 ------ flip/covariance/agkim24/generator.py | 130 ------------------------ 4 files changed, 239 deletions(-) delete mode 100644 flip/covariance/agkim24/coefficients.py delete mode 100644 flip/covariance/agkim24/fisher_terms.py delete mode 100644 flip/covariance/agkim24/flip_terms.py delete mode 100644 flip/covariance/agkim24/generator.py diff --git a/flip/covariance/agkim24/coefficients.py b/flip/covariance/agkim24/coefficients.py deleted file mode 100644 index 7436fab..0000000 --- a/flip/covariance/agkim24/coefficients.py +++ /dev/null @@ -1,17 +0,0 @@ -import numpy as np - - -def get_coefficients( - model_type, - parameter_values_dict, - variant=None, -): - coefficients_dict = {} - coefficients_dict["vv"] = [parameter_values_dict["fs8"] ** 2] - 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/agkim24/fisher_terms.py b/flip/covariance/agkim24/fisher_terms.py deleted file mode 100644 index e1ab1d9..0000000 --- a/flip/covariance/agkim24/fisher_terms.py +++ /dev/null @@ -1,62 +0,0 @@ -import numpy as np -import astropy.cosmology -from astropy import units as u -from astropy import constants as const - -class FisherTerms(object): - - def __init__(self, coordinates_velocity, cosmo=astropy.cosmology.Planck18): - self._coordinates_velocity = coordinates_velocity - self._cosmo = cosmo - - z = astropy.cosmology.z_at_value(self._cosmo.comoving_distance, self._coordinates_velocity[2]*u.Mpc) - self._cov_factors = {"v": self._cosmo.H(z)/(1+z) / self._cosmo.H0, "g": None} # need to implement "g"" - - @property - def cov_factors(self): - return self._cov_factors - - def get_partial_derivative_coefficients(self, - model_type, - parameter_values_dict, - variant=None, - ): - if variant == "growth_index": - partial_coefficients_dict = { - "Omegam": { - "vv": [ - 2 - * parameter_values_dict["Omegam"] - ** (2 * parameter_values_dict["gamma"]) - * parameter_values_dict["gamma"] - * parameter_values_dict["s8"] ** 2 - / parameter_values_dict["Omegam"], - ], - }, - "gamma": { - "vv": [ - 2 - * parameter_values_dict["Omegam"] - ** (2 * parameter_values_dict["gamma"]) - * parameter_values_dict["s8"] ** 2 - * np.log(parameter_values_dict["Omegam"]), - ], - }, - "s8": { - "vv": [ - 2 - * parameter_values_dict["Omegam"] - ** (2 * parameter_values_dict["gamma"]) - * parameter_values_dict["s8"], - ], - }, - } - else: - partial_coefficients_dict = { - "fs8": { - "vv": [ - 2 * parameter_values_dict["fs8"], - ], - }, - } - return partial_coefficients_dict diff --git a/flip/covariance/agkim24/flip_terms.py b/flip/covariance/agkim24/flip_terms.py deleted file mode 100644 index f00c9d5..0000000 --- a/flip/covariance/agkim24/flip_terms.py +++ /dev/null @@ -1,30 +0,0 @@ -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) - - -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 diff --git a/flip/covariance/agkim24/generator.py b/flip/covariance/agkim24/generator.py deleted file mode 100644 index 997cf57..0000000 --- a/flip/covariance/agkim24/generator.py +++ /dev/null @@ -1,130 +0,0 @@ -import multiprocessing as mp -from functools import partial - -import numpy as np -from scipy.special import spherical_jn - -from flip.covariance import cov_utils - - -def angle_between(ra_0, ra_1, dec_0, dec_1): - """Compute cos of the angle between r0 and r1.""" - cos_alpha = np.cos(ra_1 - ra_0) * np.cos(dec_0) * np.cos(dec_1) + np.sin( - dec_0 - ) * np.sin(dec_1) - return cos_alpha - - -def separation(r_0, r_1, cos_alpha): - """Compute separation between r_0 and r_1.""" - return np.sqrt(r_0**2 + r_1**2 - 2 * r_0 * r_1 * cos_alpha) - - -def window(r_0, r_1, cos_alpha, sep, j0kr, j2kr): - """Note: here, the bisector angle definition is used in wide-angle""" - win = 1 / 3 * (j0kr - 2 * j2kr) * cos_alpha - win += j2kr * r_0 * r_1 / sep**2 * (1 - cos_alpha**2) - return win - - -def intp(win, k, pk): - pint = win.T * pk - return np.trapz(pint, x=k) - - -def covariance_vv( - ra_in, - dec_in, - rcomov_in, - k_in, - pk_in, - grid_window_in=None, - size_batch=100_000, - number_worker=8, -): - N = len(ra_in) - - if grid_window_in is not None: - pk = pk_in * grid_window_in**2 - else: - pk = pk_in - - n_task = int((N * (N + 1)) / 2) - N - - batches = [] - for n in range(0, n_task, size_batch): - brange = np.arange(n, np.min((n + size_batch, n_task))) - i_list, j_list = cov_utils.compute_i_j(N, brange) - r_comovi, rai, deci = rcomov_in[i_list], ra_in[i_list], dec_in[i_list] - r_comovj, raj, decj = rcomov_in[j_list], ra_in[j_list], dec_in[j_list] - batches.append([rai, raj, deci, decj, r_comovi, r_comovj]) - - with mp.Pool(number_worker) as pool: - func = partial(compute_coef, k_in, pk) - pool_results = pool.map(func, batches) - values = np.concatenate(pool_results) - - var_val = np.trapz(pk / 3, x=k_in) - cov = np.insert(values, 0, var_val) - cov = 100**2 / (2 * np.pi**2) * cov - return cov - - -def compute_coef(k, pk, coord): - cos = angle_between(coord[0], coord[1], coord[2], coord[3]) - sep = separation(coord[4], coord[5], cos) - ksep = np.outer(k, sep) - j0 = spherical_jn(0, ksep) - j2 = spherical_jn(2, ksep) - res = window(coord[4], coord[5], cos, sep, j0, j2) - res = intp(res, k, pk) - return res - - -def generate_covariance( - model_type, - power_spectrum_dict, - coordinates_density=False, - coordinates_velocity=None, - **kwargs, -): - """ - The generate_covariance function generates a covariance matrix for the velocity field. - - Args: - model_type: Specify the type of model to generate - power_spectrum_dict: Pass the power spectrum to the function - coordinates_density: Specify the coordinates of the density field - coordinates_velocity: Generate the covariance matrix - **kwargs: Pass additional parameters to the function - : Generate the covariance matrix for the velocity field - The wide angle used is the bisector. - Returns: - A dictionary with a single key "vv" - - """ - assert model_type == "velocity" - cov_utils.check_generator_need( - model_type, - coordinates_density, - coordinates_velocity, - ) - number_densities = None - number_velocities = len(coordinates_velocity[0]) - cov_vv = covariance_vv( - coordinates_velocity[0], - coordinates_velocity[1], - coordinates_velocity[2], - power_spectrum_dict["vv"][0][0], - power_spectrum_dict["vv"][0][1], - **kwargs, - ) - - los_definition = "bisector" - - return ( - {"vv": np.array([cov_vv])}, - number_densities, - number_velocities, - los_definition, - ) From 8c8cdad54cc67ecd4f94a7cf9a411da63951ac31 Mon Sep 17 00:00:00 2001 From: Alex Kim Date: Mon, 22 Jul 2024 12:00:16 -0700 Subject: [PATCH 05/11] editing to make consistent with upstream --- flip/covariance/covariance.py | 17 --- flip/covariance/rcrk24/coefficients.py | 2 - flip/fisher.py | 25 +--- notebook/covariance.ipynb | 181 ------------------------- 4 files changed, 5 insertions(+), 220 deletions(-) delete mode 100644 notebook/covariance.ipynb diff --git a/flip/covariance/covariance.py b/flip/covariance/covariance.py index d02953c..44384c2 100644 --- a/flip/covariance/covariance.py +++ b/flip/covariance/covariance.py @@ -20,10 +20,7 @@ def __init__( full_matrix=False, number_densities=None, number_velocities=None, - redshift_dict=None, - power_spectrum_amplitude_function=None, variant=None, - coordinates_velocity=None, ): """ The __init__ function is called when the class is instantiated. @@ -52,10 +49,7 @@ 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 - self.coordinates_velocity = coordinates_velocity @classmethod def init_from_flip( @@ -67,7 +61,6 @@ def init_from_flip( coordinates_velocity=None, additional_parameters_values=None, los_definition="bisector", - power_spectrum_amplitude_function=None, variant=None, **kwargs, ): @@ -100,7 +93,6 @@ def init_from_flip( covariance_dict, number_densities, number_velocities, - redshift_dict, ) = generator_flip.generate_covariance( model_name, model_type, @@ -123,10 +115,7 @@ 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, - coordinates_velocity=coordinates_velocity, ) @classmethod @@ -138,7 +127,6 @@ def init_from_generator( coordinates_velocity=None, coordinates_density=None, additional_parameters_values=None, - power_spectrum_amplitude_function=None, variant=None, **kwargs, ): @@ -174,7 +162,6 @@ def init_from_generator( number_densities, number_velocities, los_definition, - redshift_dict, ) = generator.generate_covariance( model_type, power_spectrum_dict, @@ -194,8 +181,6 @@ 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, ) @@ -314,8 +299,6 @@ 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/rcrk24/coefficients.py b/flip/covariance/rcrk24/coefficients.py index f7ecb60..e501bb6 100644 --- a/flip/covariance/rcrk24/coefficients.py +++ b/flip/covariance/rcrk24/coefficients.py @@ -13,8 +13,6 @@ def get_coefficients( redshift_velocities = redshift_dict["v"] cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) - # print(power_spectrum_amplitude_function) - # qwd coefficient_vector = ( np.array(cosmo.Om(redshift_velocities)) ** parameter_values_dict["gamma"] * cosmo.H(redshift_velocities).value diff --git a/flip/fisher.py b/flip/fisher.py index 30567dc..115e0fe 100644 --- a/flip/fisher.py +++ b/flip/fisher.py @@ -97,7 +97,7 @@ def load_error_vector( def compute_covariance_derivative( self, - partial_coefficients_dict_param, cov_factors = None, + partial_coefficients_dict_param, ): if self.covariance.model_type == "density": @@ -110,14 +110,9 @@ def compute_covariance_derivative( ) elif self.covariance.model_type == "velocity": - if cov_factors is None: - cov_prefactor = 1. - else: - cov_prefactor = np.outer(cov_factors["v"],cov_factors["v"]) - covariance_derivative_sum = np.sum( [ - partial_coefficients_dict_param["vv"][i] * cov_prefactor * cov + partial_coefficients_dict_param["vv"][i] * cov for i, cov in enumerate(self.covariance.covariance_dict["vv"]) ], axis=0, @@ -173,22 +168,13 @@ def compute_fisher_matrix( variant=None, ): - if self.covariance.model_name == "agkim24": - FisherTerms = getattr(importlib.import_module( - f"flip.covariance.{self.covariance.model_name}.fisher_terms"),"FisherTerms") - coefficients = FisherTerms(self.covariance.coordinates_velocity) - cov_factors = coefficients.cov_factors - else: - coefficients = importlib.import_module( - f"flip.covariance.{self.covariance.model_name}.fisher_terms") - cov_factors = None - + coefficients = importlib.import_module( + f"flip.covariance.{self.covariance.model_name}.fisher_terms" + ) partial_coefficients_dict = coefficients.get_partial_derivative_coefficients( 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 = [] @@ -203,7 +189,6 @@ def compute_fisher_matrix( self.inverse_covariance_sum, self.compute_covariance_derivative( partial_coefficients_dict_param, - cov_factors=cov_factors ), ) ) diff --git a/notebook/covariance.ipynb b/notebook/covariance.ipynb deleted file mode 100644 index 0ac7b7f..0000000 --- a/notebook/covariance.ipynb +++ /dev/null @@ -1,181 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "id": "71e15407-b85a-4152-858f-86087f5fae7d", - "metadata": {}, - "source": [ - "This is a tutorial to use the flip package: https://github.com/corentinravoux/flip \\\n", - "It is self-contained and can be used in google collab or on your environement \\\n", - "All the data used are subsampled version of a simulation. \\\n", - "The data size is small for the tutorial, do not use it for science case. \\" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ae55978f-ac85-4f46-80cf-112204329a68", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "%%capture\n", - "!pip install git+https://github.com/corentinravoux/flip" - ] - }, - { - "cell_type": "markdown", - "id": "ef379b12-2598-413c-9f05-bb7525b4e96a", - "metadata": {}, - "source": [ - "Loading of the necessary modules" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7e1da453-ef47-4920-b330-a7d9289e2f01", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "import os\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "import pandas as pd\n", - "from flip import fitter, plot_utils, utils, vectors\n", - "from flip.covariance import covariance, contraction\n", - "from pkg_resources import resource_filename\n", - "flip_base = resource_filename(\"flip\", \".\")\n", - "data_path = os.path.join(flip_base, \"data\")\n", - "plt.style.use(os.path.join(data_path,\"style.mplstyle\"))" - ] - }, - { - "cell_type": "markdown", - "id": "7d184b4d-4672-4bf7-8f21-c156e0e195f7", - "metadata": {}, - "source": [ - "Loading the data, located in the package itself\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f81a2bc9-1d9e-4969-ac2d-0b5e56fd189c", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "ktt, ptt = np.loadtxt(os.path.join(data_path, \"power_spectrum_tt.txt\"))\n", - "kmt, pmt = np.loadtxt(os.path.join(data_path, \"power_spectrum_mt.txt\"))\n", - "kmm, pmm = np.loadtxt(os.path.join(data_path, \"power_spectrum_mm.txt\"))\n", - "\n", - "# the extra Koda et al. damping term Eq 33 from https://arxiv.org/pdf/2303.01198\n", - "sigmau_fiducial = 15\n", - "\n", - "power_spectrum_dict = {\"vv\": [[ktt, ptt * utils.Du(ktt, sigmau_fiducial) ** 2]]}\n" - ] - }, - { - "cell_type": "markdown", - "id": "a561c106-87a9-4e42-aaca-9d103e6e750d", - "metadata": {}, - "source": [ - "First part of the flip package: fast computation of the covariance matrix based on theory and SNIa coordinates" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "96947c01-98ed-4793-8fc9-9a19f0830486", - "metadata": {}, - "outputs": [], - "source": [ - "filename = \"abacussample_s500_r300_b20_d0.05_i0.1_xyz000_lnone_go.hdf5\"\n", - "import h5py\n", - "import numpy\n", - "from astropy.cosmology import Planck18 as cosmo\n", - "\n", - "def read_abacus(filename):\n", - " with h5py.File(filename, \"r\") as f:\n", - " data=f['abacussample']\n", - " zlist = list(data['z']) # observed redshift\n", - " ra = list(data['ra'])\n", - " dec = list(data['dec'])\n", - " rcom_zobs = cosmo.comoving_distance(numpy.array(zlist)).value\n", - " # cosmoz = list(data['cosmoz'])\n", - " \n", - " return numpy.array([ra,dec,rcom_zobs])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "8dfdb97c-8720-46b8-aa48-43fb1b6a4f83", - "metadata": {}, - "outputs": [], - "source": [ - "coordinates_velocity = read_abacus(filename)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a30c697f-6597-450d-ac8b-6941e9ab8219", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "size_batch=10_000\n", - "number_worker=16\n", - "\n", - "\n", - "covariance_fit = covariance.CovMatrix.init_from_flip(\n", - " \"carreres23\",\n", - " \"velocity\",\n", - " power_spectrum_dict,\n", - " coordinates_velocity=coordinates_velocity,\n", - " size_batch=size_batch,\n", - " number_worker=number_worker,\n", - ")\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "dd338dc8-69c9-40ea-ae8b-df4eb4eda84a", - "metadata": {}, - "outputs": [], - "source": [ - "covariance_fit.write(\"temp\",\"pickle\")" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3 (ipykernel)", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.12.2" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} From aea36cc9cbfff169269e72bf73a067fb2de71f36 Mon Sep 17 00:00:00 2001 From: Alex Kim Date: Mon, 22 Jul 2024 12:20:30 -0700 Subject: [PATCH 06/11] works and should match upstream --- flip/covariance/covariance.py | 14 ++++++++++++++ flip/covariance/generator.py | 6 ++---- flip/fisher.py | 2 ++ scripts/develop.py | 2 +- 4 files changed, 19 insertions(+), 5 deletions(-) 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 f44591e..482196c 100644 --- a/flip/covariance/generator.py +++ b/flip/covariance/generator.py @@ -12,9 +12,8 @@ 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.ravouxcarreres import flip_terms as flip_terms_ravouxcarreres -from flip.covariance.agkim24 import flip_terms as flip_terms_agkim24 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 log = create_log() @@ -24,8 +23,7 @@ "lai22", "carreres23", "ravouxcarreres", - "agkim24", - "rcrk24", + "rcrk24", ] 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 index 7b9fb90..7052ae9 100644 --- a/scripts/develop.py +++ b/scripts/develop.py @@ -100,4 +100,4 @@ def main(): partials = np.array([parameter_dict['gamma']*parameter_dict['Om0']**(parameter_dict['gamma']-1),np.log(parameter_dict['Om0'])*parameter_dict['Om0']**parameter_dict['gamma']]) print(partials.T @ cov[0:2,0:2] @ partials) - print(np.sqrt(partials.T @ cov[0:2,0:2] @ partials)) \ No newline at end of file + print(np.sqrt(partials.T @ cov[0:2,0:2] @ partials)) From d7a7c138d0666af1ce9c70d38101c9b465629145 Mon Sep 17 00:00:00 2001 From: Alex Kim Date: Wed, 24 Jul 2024 13:33:11 -0700 Subject: [PATCH 07/11] implementation of Om0-gamma model. The output looks reasonable. --- flip/covariance/rcrk24/fisher_terms.py | 97 +++++++++++++++++++------- 1 file changed, 70 insertions(+), 27 deletions(-) diff --git a/flip/covariance/rcrk24/fisher_terms.py b/flip/covariance/rcrk24/fisher_terms.py index adcbadc..e130ab2 100644 --- a/flip/covariance/rcrk24/fisher_terms.py +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -1,7 +1,6 @@ import numpy as np from astropy.cosmology import FlatLambdaCDM - def get_partial_derivative_coefficients( model_type, parameter_values_dict, @@ -10,39 +9,83 @@ def get_partial_derivative_coefficients( power_spectrum_amplitude_function=None, ): - cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) - redshift_velocities = redshift_dict["v"] + # The convention is to split the power spectrum into several terms + # P_ab = A * B * P0_xy + # To first order and under certain conventions + # A is the coefficient that + # P_ab = A P_xy + # B is the cofficient that + # P_xy = P0_xy + + # P0_xy is the power spectrum for a fiducial cosmology so has no parameter dependence + # The derivative cofficients we need are + # dA/dp B + A dB/dp + # for vv + # A = (aHfs8)_1 * (aHfs8_1) + # B = s8_1 * s8_2 + + redshift_velocities = redshift_dict["v"] + cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) + # calculate one 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"] - Omega_m_partial_derivative_coefficients = ( - 2 - * np.array(cosmo.Om(redshift_velocities)) ** (2 * parameter_values_dict["gamma"]) - * parameter_values_dict["gamma"] - * power_spectrum_amplitude_function(redshift_velocities) ** 2 - / np.array(cosmo.Om(redshift_velocities)) - ) + s80 = 0.832 # from Planck for the fiducial model - gamma_partial_derivative_coefficients = ( - 2 - * np.array(cosmo.Om(redshift_velocities)) ** (2 * parameter_values_dict["gamma"]) - * power_spectrum_amplitude_function(redshift_velocities) ** 2 - * np.log(cosmo.Om(redshift_velocities)) - ) + # D normalized so D(z=0)=1 + def s8(a): + return s80* a**(f0 + f0*3*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) \ + * np.exp((1-a)*f0*3* parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) + + 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 + + # for the following it is useful and fast to expand Omega in terms of (1-a) + def ds8dOm0(a): + domegapgamma = parameter_values_dict["gamma"] *parameter_values_dict["Om0"]**(parameter_values_dict["gamma"]-1) + domegapgammap1 = (parameter_values_dict["gamma"]+1) *parameter_values_dict["Om0"]**parameter_values_dict["gamma"] + return s80* np.log(a)*(domegapgamma + 3*parameter_values_dict["gamma"]*(domegapgamma - domegapgammap1)) \ + + (1-a)*3*parameter_values_dict["gamma"]*(domegapgamma-domegapgammap1) + + def ds8dgamma(a): + omegapgamma = parameter_values_dict["Om0"]**parameter_values_dict["gamma"] + lnOm0 = np.log(parameter_values_dict["Om0"]) + return s80*np.log(a) * (lnOm0 * omegapgamma * (1 + 3*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) + omegapgamma*3) \ + + 3*(1-a)*(1-parameter_values_dict["Om0"])*(lnOm0 * omegapgamma *parameter_values_dict["gamma"] + omegapgamma) + + 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)) + + a = 1/(1+redshift_velocities) + + aHfs8s8 = aHfs8(a)*s8(a) Omega_m_partial_derivative_coefficients = ( - cosmoOm ** (parameter_values_dict["gamma"]-1) - * parameter_values_dict["gamma"] - * power_spectrum_amplitude_function(redshift_velocities) + dAdOm0(a) * s8(a) + aHfs8(a) * ds8dOm0(a) ) gamma_partial_derivative_coefficients = ( - cosmoOm ** (parameter_values_dict["gamma"]) - * power_spectrum_amplitude_function(redshift_velocities) - * np.log(cosmoOm) + dAdgamma(a) * s8(a) + aHfs8(a) * ds8dgamma(a) ) - s8_partial_derivative_coefficients = ( 2 * cosmoOm ** (2 * parameter_values_dict["gamma"]) @@ -54,10 +97,10 @@ def get_partial_derivative_coefficients( "vv": [ np.outer( Omega_m_partial_derivative_coefficients, - cosmoOm ** (parameter_values_dict["gamma"]), + aHfs8s8, ) + np.outer( - cosmoOm ** (parameter_values_dict["gamma"]), + aHfs8s8, Omega_m_partial_derivative_coefficients, ), ], @@ -66,10 +109,10 @@ def get_partial_derivative_coefficients( "vv": [ np.outer( gamma_partial_derivative_coefficients, - cosmoOm ** (parameter_values_dict["gamma"]), + aHfs8s8, ) + np.outer( - cosmoOm ** (parameter_values_dict["gamma"]), + aHfs8s8, gamma_partial_derivative_coefficients, ), ], From a830aa34d4a7b463f3ec6756857c6cb9fff9d5cd Mon Sep 17 00:00:00 2001 From: Alex Kim Date: Wed, 24 Jul 2024 16:08:04 -0700 Subject: [PATCH 08/11] calculate derivatives with wolfram. organize code --- flip/covariance/rcrk24/fisher_terms.py | 83 ++++++++++++++++---------- 1 file changed, 53 insertions(+), 30 deletions(-) diff --git a/flip/covariance/rcrk24/fisher_terms.py b/flip/covariance/rcrk24/fisher_terms.py index e130ab2..70d790c 100644 --- a/flip/covariance/rcrk24/fisher_terms.py +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -1,6 +1,24 @@ 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, @@ -9,35 +27,46 @@ def get_partial_derivative_coefficients( power_spectrum_amplitude_function=None, ): - # The convention is to split the power spectrum into several terms - # P_ab = A * B * P0_xy - # To first order and under certain conventions - # A is the coefficient that - # P_ab = A P_xy - # B is the cofficient that - # P_xy = P0_xy - - # P0_xy is the power spectrum for a fiducial cosmology so has no parameter dependence - # The derivative cofficients we need are - # dA/dp B + A dB/dp + cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) + redshift_velocities = redshift_dict["v"] - # for vv - # A = (aHfs8)_1 * (aHfs8_1) - # B = s8_1 * s8_2 + # The Om0-gamma model f=Omega(Om0)^gamma - redshift_velocities = redshift_dict["v"] - cosmo = FlatLambdaCDM(H0=100, Om0=parameter_values_dict["Om0"]) - # calculate one a few things that are used repeatedly + # 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"] - s80 = 0.832 # from Planck for the fiducial model + # In the Om0-gamma parameterization fs80 is considered to be fixed by the CMB + # and is hence not a fit parameter + + s80 = 0.832 + + # 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"]) \ + ) - # D normalized so D(z=0)=1 def s8(a): - return s80* a**(f0 + f0*3*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) \ - * np.exp((1-a)*f0*3* parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) + return s80 * np.exp(lnD(a)) def aHfs8(a): return f * s8(a) / (1+redshift_velocities) * cosmo.H(redshift_velocities)/cosmo.H0 @@ -55,18 +84,11 @@ def dfdOm0(a): def dfdgamma(a): return np.log(cosmoOm) * f - # for the following it is useful and fast to expand Omega in terms of (1-a) def ds8dOm0(a): - domegapgamma = parameter_values_dict["gamma"] *parameter_values_dict["Om0"]**(parameter_values_dict["gamma"]-1) - domegapgammap1 = (parameter_values_dict["gamma"]+1) *parameter_values_dict["Om0"]**parameter_values_dict["gamma"] - return s80* np.log(a)*(domegapgamma + 3*parameter_values_dict["gamma"]*(domegapgamma - domegapgammap1)) \ - + (1-a)*3*parameter_values_dict["gamma"]*(domegapgamma-domegapgammap1) + return s8(a)*dlnDdOm0(a) def ds8dgamma(a): - omegapgamma = parameter_values_dict["Om0"]**parameter_values_dict["gamma"] - lnOm0 = np.log(parameter_values_dict["Om0"]) - return s80*np.log(a) * (lnOm0 * omegapgamma * (1 + 3*parameter_values_dict["gamma"]*(1-parameter_values_dict["Om0"])) + omegapgamma*3) \ - + 3*(1-a)*(1-parameter_values_dict["Om0"])*(lnOm0 * omegapgamma *parameter_values_dict["gamma"] + omegapgamma) + return s8(a)*dlnDdgamma(a) def dAdOm0(a): return a*cosmo.H(a)/cosmo.H0*(dfdOm0(a)*s8(a) + f*ds8dOm0(a)) @@ -86,6 +108,7 @@ def dAdgamma(a): dAdgamma(a) * s8(a) + aHfs8(a) * ds8dgamma(a) ) + # not yet vetted s8_partial_derivative_coefficients = ( 2 * cosmoOm ** (2 * parameter_values_dict["gamma"]) From 933fbd4bc1002ae3f605a95261c5b275b041882a Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 25 Jul 2024 11:09:21 -0700 Subject: [PATCH 09/11] implement fs8 fisher --- flip/covariance/rcrk24/fisher_terms.py | 39 +++++++++++++++++--------- scripts/develop.py | 19 +++++++------ 2 files changed, 36 insertions(+), 22 deletions(-) diff --git a/flip/covariance/rcrk24/fisher_terms.py b/flip/covariance/rcrk24/fisher_terms.py index 70d790c..9a149cb 100644 --- a/flip/covariance/rcrk24/fisher_terms.py +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -27,8 +27,14 @@ 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 @@ -37,10 +43,7 @@ def get_partial_derivative_coefficients( f = cosmoOm**parameter_values_dict["gamma"] f0 = parameter_values_dict["Om0"]**parameter_values_dict["gamma"] - # In the Om0-gamma parameterization fs80 is considered to be fixed by the CMB - # and is hence not a fit parameter - s80 = 0.832 # Calculation of s8 and its derivatives requires an integral. It is useful to # expand Omega in terms of (1-a), which allows analytic solutions @@ -96,7 +99,6 @@ def dAdOm0(a): def dAdgamma(a): return a*cosmo.H(a)/cosmo.H0*(dfdgamma(a)*s8(a) + f*ds8dgamma(a)) - a = 1/(1+redshift_velocities) aHfs8s8 = aHfs8(a)*s8(a) @@ -108,13 +110,20 @@ def dAdgamma(a): dAdgamma(a) * s8(a) + aHfs8(a) * ds8dgamma(a) ) - # not yet vetted - s8_partial_derivative_coefficients = ( - 2 - * cosmoOm ** (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": [ @@ -140,11 +149,15 @@ def dAdgamma(a): ), ], }, - "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/scripts/develop.py b/scripts/develop.py index 7052ae9..a13ee03 100644 --- a/scripts/develop.py +++ b/scripts/develop.py @@ -62,16 +62,18 @@ def main(): 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, - } + # parameter_dict = { + # "Om0": 0.3, + # "gamma": 0.55, + # "sigv": 200, + # "sigma_M": 0.12, + # } Fisher = fisher.FisherMatrix.init_from_covariance( covariance_fit, @@ -94,10 +96,9 @@ def main(): "sigma_M": 0.12, } parameter_name_list, fisher_matrix = main() - print(fisher_matrix) cov = np.linalg.inv(fisher_matrix[0:2,0:2]) - print(cov[0:2,0:2]) partials = np.array([parameter_dict['gamma']*parameter_dict['Om0']**(parameter_dict['gamma']-1),np.log(parameter_dict['Om0'])*parameter_dict['Om0']**parameter_dict['gamma']]) - print(partials.T @ cov[0:2,0:2] @ partials) + print(np.sqrt(partials.T @ cov[0:2,0:2] @ partials)) + print(1/np.sqrt(fisher_matrix[2,2])) From 319ad0f6e289cb3734bdc8a0e2bb0fb8aea33c37 Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 25 Jul 2024 16:04:48 -0700 Subject: [PATCH 10/11] confirmed first order (1-a) using wolfram --- scripts/develop.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/scripts/develop.py b/scripts/develop.py index a13ee03..1de4cc6 100644 --- a/scripts/develop.py +++ b/scripts/develop.py @@ -87,6 +87,32 @@ def main(): ) 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 = { @@ -98,7 +124,9 @@ def main(): parameter_name_list, fisher_matrix = main() cov = np.linalg.inv(fisher_matrix[0:2,0:2]) - partials = np.array([parameter_dict['gamma']*parameter_dict['Om0']**(parameter_dict['gamma']-1),np.log(parameter_dict['Om0'])*parameter_dict['Om0']**parameter_dict['gamma']]) + 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])) From f6c98a33d22d0c9840cf48f2b0a54638d4385eae Mon Sep 17 00:00:00 2001 From: Alex Date: Thu, 25 Jul 2024 18:15:10 -0700 Subject: [PATCH 11/11] more debug code. Fix carreres to work in new API --- flip/covariance/carreres23/coefficients.py | 2 ++ flip/covariance/rcrk24/fisher_terms.py | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) 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/rcrk24/fisher_terms.py b/flip/covariance/rcrk24/fisher_terms.py index 9a149cb..50456c2 100644 --- a/flip/covariance/rcrk24/fisher_terms.py +++ b/flip/covariance/rcrk24/fisher_terms.py @@ -118,8 +118,8 @@ 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)) + 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)