diff --git a/flip/fitter.py b/flip/fitter.py index 3d1ae70..d01ee35 100644 --- a/flip/fitter.py +++ b/flip/fitter.py @@ -163,7 +163,7 @@ def init_from_covariance( data, parameter_dict, likelihood_type="multivariate_gaussian", - likelihood_properties=None, + likelihood_properties={}, **kwargs, ): """ @@ -194,7 +194,7 @@ def init_from_covariance( likelihood = minuit_fitter.get_likelihood( parameter_dict, likelihood_type=likelihood_type, - likelihood_properties=likelihood_properties, + likelihood_properties={**likelihood_properties, 'nloglik': True}, **kwargs, ) minuit_fitter.likelihood = likelihood @@ -247,7 +247,7 @@ def init_from_file( data, parameter_dict, likelihood_type=likelihood_type, - likelihood_properties=likelihood_properties, + likelihood_properties={**likelihood_properties, 'nloglik': True}, ) def setup_minuit(self, parameter_dict): diff --git a/flip/likelihood.py b/flip/likelihood.py index 5cc15c2..085d647 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -11,41 +11,36 @@ def log_likelihood_gaussian_inverse(vector, covariance_sum): _, logdet = np.linalg.slogdet(covariance_sum) inverse_covariance_sum = np.linalg.inv(covariance_sum) chi2 = np.dot(vector, np.dot(inverse_covariance_sum, vector)) - return 0.5 * (vector.size * np.log(2 * np.pi) + logdet + chi2) + return -0.5 * (vector.size * np.log(2 * np.pi) + logdet + chi2) def log_likelihood_gaussian_cholesky(vector, covariance_sum): cholesky = sc.linalg.cho_factor(covariance_sum) logdet = 2 * np.sum(np.log(np.diag(cholesky[0]))) chi2 = np.dot(vector, sc.linalg.cho_solve(cholesky, vector)) - return 0.5 * (vector.size * np.log(2 * np.pi) + logdet + chi2) + return -0.5 * (vector.size * np.log(2 * np.pi) + logdet + chi2) class BaseLikelihood(object): + + _default_likelihood_properties = { + "inversion_method": "inverse", + "velocity_type": "direct", + "velocity_estimator": "full", + } + def __init__( self, covariance=None, data=None, parameter_names=None, - likelihood_properties=None, + likelihood_properties={}, ): self.covariance = covariance self.data = data self.parameter_names = parameter_names - _default_likelihood_properties = { - "inversion_method": "inverse", - "velocity_type": "direct", - "velocity_estimator": "full", - } - if likelihood_properties == None: - likelihood_properties = _default_likelihood_properties - else: - for key in _default_likelihood_properties.keys(): - if key not in likelihood_properties.keys(): - likelihood_properties[key] = _default_likelihood_properties[key] - - self.likelihood_properties = likelihood_properties + self.likelihood_properties = {**self._default_likelihood_properties, **likelihood_properties} @classmethod def init_from_covariance( @@ -53,7 +48,7 @@ def init_from_covariance( covariance, data, parameter_names, - likelihood_properties=None, + likelihood_properties={}, **kwargs, ): """ @@ -116,12 +111,18 @@ def load_data_vector( class MultivariateGaussianLikelihood(BaseLikelihood): + _default_likelihood_properties = { + 'nloglik': False, + **BaseLikelihood._default_likelihood_properties + } + def __init__( self, covariance=None, data=None, parameter_names=None, - likelihood_properties=None, + likelihood_properties={}, + negloglik=False ): super(MultivariateGaussianLikelihood, self).__init__( covariance=covariance, @@ -130,6 +131,7 @@ def __init__( likelihood_properties=likelihood_properties, ) + def verify_covariance(self): if self.covariance.full_matrix is False: self.covariance.compute_full_matrix() @@ -147,6 +149,8 @@ def __call__(self, parameter_values): likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) + if self.likelihood_properties['nloglik']: + return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum) @@ -156,7 +160,7 @@ def __init__( covariance=None, data=None, parameter_names=None, - likelihood_properties=None, + likelihood_properties={}, interpolation_value_name=None, interpolation_value_range=None, ): @@ -241,7 +245,7 @@ def __init__( covariance=None, data=None, parameter_names=None, - likelihood_properties=None, + likelihood_properties={}, interpolation_value_range_0=None, interpolation_value_range_1=None, ): diff --git a/notebook/fit_velocity_old.ipynb b/notebook/fit_velocity_old.ipynb new file mode 100644 index 0000000..5a98161 --- /dev/null +++ b/notebook/fit_velocity_old.ipynb @@ -0,0 +1,520 @@ +{ + "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 +}