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/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/flip/covariance/rcrk24/__init__.py b/flip/covariance/rcrk24/__init__.py
new file mode 100644
index 0000000..e69de29
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",
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]))
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