From a45e4786dd6488abebcf377e5daba52b9bba2c8e Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Wed, 31 Jul 2024 12:19:02 +0200 Subject: [PATCH 1/4] fixing init package --- flip/covariance/__init__.py | 5 +++-- flip/covariance/rcrk24/__init__.py | 0 2 files changed, 3 insertions(+), 2 deletions(-) create mode 100644 flip/covariance/rcrk24/__init__.py diff --git a/flip/covariance/__init__.py b/flip/covariance/__init__.py index 6b63626..638da73 100644 --- a/flip/covariance/__init__.py +++ b/flip/covariance/__init__.py @@ -1,3 +1,4 @@ """Init file of the flip.covariance package.""" -from . import carreres23, adamsblake20, lai22, ravouxcarreres -from .covariance import CovMatrix \ No newline at end of file + +from . import adamsblake17plane, adamsblake20, carreres23, lai22, ravouxcarreres, rcrk24 +from .covariance import CovMatrix diff --git a/flip/covariance/rcrk24/__init__.py b/flip/covariance/rcrk24/__init__.py new file mode 100644 index 0000000..e69de29 From ec3ad42c7d6b44b2ac67ba8a76a7b35c01ec72b2 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Wed, 31 Jul 2024 14:27:22 +0200 Subject: [PATCH 2/4] module compatibility --- flip/covariance/generator.py | 4 ++-- flip/covariance/lai22/generator.py | 2 +- setup.cfg | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/flip/covariance/generator.py b/flip/covariance/generator.py index 482196c..da21546 100644 --- a/flip/covariance/generator.py +++ b/flip/covariance/generator.py @@ -12,8 +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.rcrk24 import flip_terms as flip_terms_rcrk24 from flip.covariance.ravouxcarreres import flip_terms as flip_terms_ravouxcarreres +from flip.covariance.rcrk24 import flip_terms as flip_terms_rcrk24 from flip.utils import create_log log = create_log() @@ -49,7 +49,7 @@ def correlation_integration(l, r, k, integrand): integrand = ( (-1) ** (l // 2) * (k**2 / (2 * np.pi**2)) * integrand * spherical_jn(l, kr).T ) - return (-1) ** (l % 2) * integrate.simps(integrand, x=k) + return (-1) ** (l % 2) * integrate.simpson(integrand, x=k) def correlation_hankel(l, r, k, integrand, hankel_overhead_coefficient=2): diff --git a/flip/covariance/lai22/generator.py b/flip/covariance/lai22/generator.py index 51fd739..a09f283 100644 --- a/flip/covariance/lai22/generator.py +++ b/flip/covariance/lai22/generator.py @@ -17,7 +17,7 @@ def compute_correlation_coefficient_simple_integration(p, q, l, r, k, pk): 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) - return integrate.simps(integrand, x=k) + return integrate.simpson(integrand, x=k) def compute_correlation_coefficient_hankel( diff --git a/setup.cfg b/setup.cfg index 3a5af9a..487c90d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -12,7 +12,7 @@ license = MIT packages = find: install_requires = pandas - numpy + numpy <= 1.26.4 scipy matplotlib importlib-metadata From e4c53e6ff4358f7d8c5cfeab89ef609a802dabe9 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Wed, 31 Jul 2024 14:33:38 +0200 Subject: [PATCH 3/4] updating notebook --- notebook/fit_velocity.ipynb | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/notebook/fit_velocity.ipynb b/notebook/fit_velocity.ipynb index 7486b19..3621529 100644 --- a/notebook/fit_velocity.ipynb +++ b/notebook/fit_velocity.ipynb @@ -235,12 +235,18 @@ "\n", "\n", "\n", - "mu, variance_mu = vectors.compute_observed_distance_modulus(\n", + "mu = vectors.compute_observed_distance_modulus(\n", " data_velocity,\n", " test_parameters[\"alpha\"],\n", " test_parameters[\"beta\"],\n", " test_parameters[\"M_0\"],\n", " )\n", + "\n", + "variance_mu = vectors.compute_observed_distance_modulus_error(\n", + " data_velocity,\n", + " test_parameters[\"alpha\"],\n", + " test_parameters[\"beta\"],\n", + " )\n", "mask = mu > 0\n", "mask &= np.sqrt(variance_mu) - 0.1* mu < 0\n", "\n", From d9d5c57d06ab3ff83f7d8bdb6cd14153e750c7fa Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Wed, 31 Jul 2024 14:33:53 +0200 Subject: [PATCH 4/4] removing tests scripts and notebooks --- notebook/fit_velocity_test_bastien.ipynb | 520 ----------------------- scripts/develop.py | 132 ------ 2 files changed, 652 deletions(-) delete mode 100644 notebook/fit_velocity_test_bastien.ipynb delete mode 100644 scripts/develop.py diff --git a/notebook/fit_velocity_test_bastien.ipynb b/notebook/fit_velocity_test_bastien.ipynb deleted file mode 100644 index 5a98161..0000000 --- a/notebook/fit_velocity_test_bastien.ipynb +++ /dev/null @@ -1,520 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 1, - "id": "7e1da453-ef47-4920-b330-a7d9289e2f01", - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "[ 000000.00 ]: 03-04 10:14 root INFO No pypower module detected, gridding with this method is unavailable\n", - "[ 000000.49 ]: 03-04 10:14 jax._src.xla_bridge INFO Unable to initialize backend 'cuda': Unable to load CUDA. Is it installed?\n", - "[ 000000.49 ]: 03-04 10:14 jax._src.xla_bridge INFO Unable to initialize backend 'rocm': NOT_FOUND: Could not find registered platform with name: \"rocm\". Available platform names are: CUDA\n", - "[ 000000.50 ]: 03-04 10:14 jax._src.xla_bridge INFO Unable to initialize backend 'tpu': INTERNAL: Failed to open libtpu.so: libtpu.so: cannot open shared object file: No such file or directory\n", - "/tmp/ipykernel_13812/968461804.py:6: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html\n", - " from pkg_resources import resource_filename\n" - ] - } - ], - "source": [ - "import os\n", - "import numpy as np\n", - "import pandas as pd\n", - "from flip import fitter, plot_utils, utils\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\")" - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "id": "f81a2bc9-1d9e-4969-ac2d-0b5e56fd189c", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "sn_data = pd.read_parquet(os.path.join(data_path, \"velocity_data.parquet\"))\n", - "\n", - "# Removing bad fits\n", - "sn_data = sn_data[np.array(sn_data[\"status\"]) != False]\n", - "sn_data = sn_data[np.array(sn_data[\"status\"]) != None]\n", - "\n", - "coordinates_velocity = np.array([sn_data[\"ra\"], sn_data[\"dec\"], sn_data[\"como_dist\"]])\n", - "\n", - "data_velocity = sn_data.to_dict(\"list\")\n", - "for key in data_velocity.keys():\n", - " data_velocity[key] = np.array(data_velocity[key])\n", - "data_velocity[\"velocity\"] = data_velocity.pop(\"vpec\")\n", - "data_velocity[\"velocity_error\"] = np.zeros_like(data_velocity[\"velocity\"])\n", - "\n", - "\n", - "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", - "sigmau_fiducial = 15\n", - "\n", - "power_spectrum_dict = {\"vv\": [[ktt, ptt * utils.Du(ktt, sigmau_fiducial) ** 2]]}\n" - ] - }, - { - "cell_type": "markdown", - "id": "31cc34e2-751d-4880-abf0-fdbe48eb9089", - "metadata": {}, - "source": [ - "# Fit true velocity with MINUIT" - ] - }, - { - "cell_type": "code", - "execution_count": 3, - "id": "a30c697f-6597-450d-ac8b-6941e9ab8219", - "metadata": { - "tags": [] - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "[ 000000.72 ]: 03-04 10:14 root INFO Covariance matrix generated from flip with carreres23 model in 6.53e-01 seconds\n" - ] - } - ], - "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": 4, - "id": "556ad2da-dbff-48b2-ad69-24695dda48f2", - "metadata": { - "tags": [] - }, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "[ 000005.13 ]: 03-04 10:14 root INFO ┌─────────────────────────────────────────────────────────────────────────┐\n", - "│ Migrad │\n", - "├──────────────────────────────────┬──────────────────────────────────────┤\n", - "│ FCN = 1.311e+04 │ Nfcn = 39 │\n", - "│ EDM = 5.59e-07 (Goal: 0.0002) │ time = 4.3 sec │\n", - "├──────────────────────────────────┼──────────────────────────────────────┤\n", - "│ Valid Minimum │ Below EDM threshold (goal x 10) │\n", - "├──────────────────────────────────┼──────────────────────────────────────┤\n", - "│ No parameters at limit │ Below call limit │\n", - "├──────────────────────────────────┼──────────────────────────────────────┤\n", - "│ Hesse ok │ Covariance accurate │\n", - "└──────────────────────────────────┴──────────────────────────────────────┘\n", - "┌───┬──────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", - "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", - "├───┼──────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", - "│ 0 │ fs8 │ 0.435 │ 0.025 │ │ │ 0 │ │ │\n", - "│ 1 │ sigv │ 172 │ 8 │ │ │ 0 │ │ │\n", - "└───┴──────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", - "┌──────┬─────────────────────┐\n", - "│ │ fs8 sigv │\n", - "├──────┼─────────────────────┤\n", - "│ fs8 │ 0.000631 -128.6e-3 │\n", - "│ sigv │ -128.6e-3 69.2 │\n", - "└──────┴─────────────────────┘\n" - ] - } - ], - "source": [ - "\n", - "likelihood_type=\"multivariate_gaussian\"\n", - "likelihood_properties = {\"inversion_method\": \"cholesky\",\n", - " \"velocity_type\": \"direct\"}\n", - "\n", - "\n", - "parameter_dict = {\"fs8\": {\"value\": 0.4,\n", - " \"limit_low\" : 0.0,\n", - " \"fixed\" : False,\n", - " },\n", - " \"sigv\": {\"value\": 200,\n", - " \"limit_low\" : 0.0,\n", - " \"fixed\" : False,\n", - " },\n", - " }\n", - "\n", - "\n", - "\n", - "minuit_fitter = fitter.FitMinuit.init_from_covariance(\n", - " covariance_fit,\n", - " data_velocity,\n", - " parameter_dict,\n", - " likelihood_type=likelihood_type,\n", - " likelihood_properties = likelihood_properties\n", - ")\n", - "\n", - "\n", - "minuit_fitter.run()\n" - ] - }, - { - "cell_type": "markdown", - "id": "bf056026", - "metadata": {}, - "source": [ - "# Fit true vel, exemple with mcmc" - ] - }, - { - "cell_type": "code", - "execution_count": 5, - "id": "4ffe40b6", - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "likelihood_type=\"multivariate_gaussian\"\n", - "likelihood_properties = {\n", - " \"inversion_method\": \"cholesky\",\n", - " \"velocity_type\": \"direct\"\n", - " }\n", - "\n", - "\n", - "parameter_dict = {\n", - " \"fs8\": {\n", - " \"randfun\": lambda size: np.random.normal(0.4, 0.001, size=size),\n", - " },\n", - " \"sigv\": {\n", - " \"randfun\": lambda size: np.random.normal(200, 10, size=size),\n", - " },\n", - " }\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": 6, - "id": "a544ca76", - "metadata": {}, - "outputs": [], - "source": [ - "covariance = covariance_fit\n", - "data = data_velocity\n", - "\n", - "mcmc_fitter = fitter.FitMCMC.init_from_covariance(\n", - " covariance=covariance,\n", - " data=data,\n", - " parameter_dict=parameter_dict,\n", - " backend_file=None\n", - ")\n", - "likelihood = mcmc_fitter.get_likelihood(\n", - " parameter_dict,\n", - " likelihood_type=likelihood_type,\n", - " likelihood_properties=likelihood_properties,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 7, - "id": "6d54a11a", - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "mcmc_fitter = fitter.FitMCMC.init_from_covariance(\n", - " covariance_fit,\n", - " data_velocity,\n", - " parameter_dict,\n", - " likelihood_type=likelihood_type,\n", - " likelihood_properties = likelihood_properties,\n", - " nwalkers=4\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": 8, - "id": "d7f9a07a", - "metadata": {}, - "outputs": [], - "source": [ - "sampler = mcmc_fitter.sampler.run_chains(10)" - ] - }, - { - "cell_type": "markdown", - "id": "a164ce4a-0e0a-48a1-af3b-ae8c7e07a1b6", - "metadata": {}, - "source": [ - "# Fit true velocity interpolation sigma_u" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "id": "9c22818b-2e8e-4bd3-a546-7f0c4e1a1b67", - "metadata": { - "tags": [] - }, - "outputs": [ - { - "ename": "AttributeError", - "evalue": "'CovMatrix' object has no attribute 'CovMatrix'", - "output_type": "error", - "traceback": [ - "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m/home/blc56/Documents/Soft/flip/notebook/fit_velocity.ipynb Cell 12\u001b[0m line \u001b[0;36m1\n\u001b[1;32m 7\u001b[0m \u001b[39mfor\u001b[39;00m sigu \u001b[39min\u001b[39;00m sigmau_list:\n\u001b[1;32m 9\u001b[0m power_spectrum_dict \u001b[39m=\u001b[39m {\u001b[39m\"\u001b[39m\u001b[39mvv\u001b[39m\u001b[39m\"\u001b[39m: [[ktt, ptt\u001b[39m*\u001b[39m utils\u001b[39m.\u001b[39mDu(ktt, sigu)\u001b[39m*\u001b[39m\u001b[39m*\u001b[39m\u001b[39m2\u001b[39m]]}\n\u001b[0;32m---> 11\u001b[0m covariance_list\u001b[39m.\u001b[39mappend(covariance\u001b[39m.\u001b[39;49mCovMatrix\u001b[39m.\u001b[39minit_from_flip(\n\u001b[1;32m 12\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mcarreres23\u001b[39m\u001b[39m\"\u001b[39m,\n\u001b[1;32m 13\u001b[0m \u001b[39m\"\u001b[39m\u001b[39mvelocity\u001b[39m\u001b[39m\"\u001b[39m,\n\u001b[1;32m 14\u001b[0m power_spectrum_dict,\n\u001b[1;32m 15\u001b[0m coordinates_velocity\u001b[39m=\u001b[39mcoordinates_velocity,\n\u001b[1;32m 16\u001b[0m size_batch\u001b[39m=\u001b[39msize_batch,\n\u001b[1;32m 17\u001b[0m number_worker\u001b[39m=\u001b[39mnumber_worker,\n\u001b[1;32m 18\u001b[0m ))\n", - "\u001b[0;31mAttributeError\u001b[0m: 'CovMatrix' object has no attribute 'CovMatrix'" - ] - } - ], - "source": [ - "sigmau_list = np.linspace(10.0, 20.0,10)\n", - "covariance_list = []\n", - "size_batch=10_000\n", - "number_worker=16\n", - "\n", - "\n", - "for sigu in sigmau_list:\n", - " \n", - " power_spectrum_dict = {\"vv\": [[ktt, ptt* utils.Du(ktt, sigu)**2]]}\n", - "\n", - " covariance_list.append(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": "5034bd62-9be1-4b41-9c12-b1e45cd7c8a7", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "likelihood_type=\"multivariate_gaussian_interp1d\"\n", - "likelihood_properties = {\"inversion_method\": \"cholesky\",\n", - " \"velocity_type\": \"direct\"}\n", - "\n", - "parameter_dict = {\"fs8\": {\"value\": 0.4,\n", - " \"limit_low\" : 0.0,\n", - " \"fixed\" : False,\n", - " },\n", - " \"sigv\": {\"value\": 200,\n", - " \"limit_low\" : 0.0,\n", - " \"fixed\" : False,\n", - " },\n", - " \"sigu\": {\"value\": 15.0,\n", - " \"limit_low\" : 13.0,\n", - " \"limit_up\" : 17.0,\n", - " \"fixed\" : False,\n", - " },\n", - " }\n", - "\n", - "\n", - "minuit_fitter = fitter.FitMinuit.init_from_covariance(\n", - " covariance_list,\n", - " data_velocity,\n", - " parameter_dict,\n", - " likelihood_type=likelihood_type,\n", - " likelihood_properties = likelihood_properties,\n", - " interpolation_value_name = \"sigu\",\n", - " interpolation_value_range = sigmau_list,\n", - ")\n", - "\n", - "\n", - "minuit_fitter.run()\n" - ] - }, - { - "cell_type": "markdown", - "id": "453720f6-4c79-494e-b83e-ba7c5f38b47b", - "metadata": {}, - "source": [ - "# Fit estimate velocity - do not converge for now" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5de61c3f-0589-4a8d-b940-86e6e6f005af", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "\n", - "likelihood_type=\"multivariate_gaussian\"\n", - "likelihood_properties = {\"inversion_method\": \"cholesky\",\n", - " \"velocity_type\": \"saltfit\",\n", - " \"velocity_estimator\": \"full\",\n", - " }\n", - "\n", - "\n", - "parameter_dict = {\"fs8\": {\"value\": 0.4,\n", - " \"limit_low\" : 0.0,\n", - " \"limit_up\" : 1.0,\n", - " \"fixed\" : False,\n", - " },\n", - " \"sigv\": {\"value\": 200,\n", - " \"limit_low\" : 0.0,\n", - " \"limit_up\" : 300,\n", - " \"fixed\" : False,\n", - " },\n", - " \"alpha\": {\"value\": 0.1,\n", - " \"limit_low\" : 0.05,\n", - " \"limit_up\" : 0.15,\n", - " \"fixed\" : False,\n", - " },\n", - " \"beta\": {\"value\": 3.0,\n", - " \"limit_low\" : 1.5,\n", - " \"limit_up\" : 4.5,\n", - " \"fixed\" : False,\n", - " },\n", - " \"M_0\": {\"value\": -19,\n", - " \"limit_low\" : -21,\n", - " \"limit_up\" : -18,\n", - " \"fixed\" : False,\n", - " },\n", - " \"sigma_M\": {\"value\": 0.1,\n", - " \"limit_low\" : 0.0,\n", - " \"limit_up\" : 1.0,\n", - " \"fixed\" : False,\n", - " }}\n", - "\n", - "minuit_fitter = fitter.FitMinuit.init_from_covariance(\n", - " covariance_fit,\n", - " data_velocity,\n", - " parameter_dict,\n", - " likelihood_type=likelihood_type,\n", - " likelihood_properties = likelihood_properties\n", - ")\n", - "\n", - "\n", - "minuit_fitter.run()\n" - ] - }, - { - "cell_type": "markdown", - "id": "41ab89c8-8772-4ff7-8c9e-422fc2801cfc", - "metadata": {}, - "source": [ - "# Some plots" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e30a13f8-2304-4649-b269-b9a24bb0501d", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "plt.style.use(os.path.join(data_path,\"style.mplstyle\"))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "62b90723-fb72-4c46-9d9e-969a83bf7f74", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "plot_utils.plot_correlation_from_likelihood(\n", - " minuit_fitter.likelihood,\n", - " parameter_dict,\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e7a38c68-4507-4494-9207-8a65e4cffb81", - "metadata": { - "tags": [] - }, - "outputs": [], - "source": [ - "\n", - "r_perpendicular = np.linspace(-200,200,200)\n", - "r_parallel = np.linspace(-200,200,200)\n", - "\n", - "r_reference_parallel = 30 \n", - "r_reference_perpendicular = 0 \n", - "\n", - "contraction_plot = contraction.Contraction.init_from_flip(\n", - " \"carreres23\",\n", - " \"velocity\",\n", - " power_spectrum_dict,\n", - " r_perpendicular,\n", - " r_parallel,\n", - " r_reference_perpendicular,\n", - " r_reference_parallel,\n", - " coordinate_type=\"rprt\",\n", - ")\n", - "\n", - "parameter_values_dict = {parameters: parameter_dict[parameters][\"value\"] for parameters in parameter_dict}\n", - "\n", - "plot_utils.plot_2d_contraction(contraction_plot, \n", - " parameter_values_dict, \n", - " rs_multiplied = True)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "930dc56e-9eb8-46b9-b256-220cd4fa999a", - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "python3.10", - "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.10.13" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/scripts/develop.py b/scripts/develop.py deleted file mode 100644 index 1de4cc6..0000000 --- a/scripts/develop.py +++ /dev/null @@ -1,132 +0,0 @@ -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]))