From 641a4b23c1f47abea8206257c9b27b5ddb58cbdb Mon Sep 17 00:00:00 2001 From: Bastien Carreres Date: Wed, 13 Mar 2024 16:18:19 -0400 Subject: [PATCH 1/8] Add a negloglik feature --- flip/fitter.py | 6 +- flip/likelihood.py | 44 +-- notebook/fit_velocity_old.ipynb | 520 ++++++++++++++++++++++++++++++++ 3 files changed, 547 insertions(+), 23 deletions(-) create mode 100644 notebook/fit_velocity_old.ipynb 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 +} From 31fbf8e9eb83af3e7e7e1c1a4761d0f1ba7dde39 Mon Sep 17 00:00:00 2001 From: Bastien Carreres Date: Mon, 18 Mar 2024 16:57:51 -0400 Subject: [PATCH 2/8] add scalar vmean --- flip/likelihood.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/flip/likelihood.py b/flip/likelihood.py index 085d647..40a8068 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -122,7 +122,7 @@ def __init__( data=None, parameter_names=None, likelihood_properties={}, - negloglik=False + negloglik=False, ): super(MultivariateGaussianLikelihood, self).__init__( covariance=covariance, @@ -130,8 +130,7 @@ def __init__( parameter_names=parameter_names, likelihood_properties=likelihood_properties, ) - - + def verify_covariance(self): if self.covariance.full_matrix is False: self.covariance.compute_full_matrix() @@ -150,8 +149,8 @@ def __call__(self, parameter_values): 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) + return -likelihood_function(vector - parameter_values_dict["vmean"], covariance_sum) + return likelihood_function(vector - parameter_values_dict["vmean"], covariance_sum) class MultivariateGaussianLikelihoodInterpolate1D(BaseLikelihood): From 654949e6fcf28bd393983f693ce4ca1e12b84362 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 8 Apr 2024 13:23:55 +0200 Subject: [PATCH 3/8] formating --- flip/likelihood.py | 31 +++++++++++++++++++------------ 1 file changed, 19 insertions(+), 12 deletions(-) diff --git a/flip/likelihood.py b/flip/likelihood.py index 40a8068..4976832 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -24,10 +24,10 @@ def log_likelihood_gaussian_cholesky(vector, covariance_sum): class BaseLikelihood(object): _default_likelihood_properties = { - "inversion_method": "inverse", - "velocity_type": "direct", - "velocity_estimator": "full", - } + "inversion_method": "inverse", + "velocity_type": "direct", + "velocity_estimator": "full", + } def __init__( self, @@ -40,7 +40,10 @@ def __init__( self.data = data self.parameter_names = parameter_names - self.likelihood_properties = {**self._default_likelihood_properties, **likelihood_properties} + self.likelihood_properties = { + **self._default_likelihood_properties, + **likelihood_properties, + } @classmethod def init_from_covariance( @@ -112,9 +115,9 @@ def load_data_vector( class MultivariateGaussianLikelihood(BaseLikelihood): _default_likelihood_properties = { - 'nloglik': False, - **BaseLikelihood._default_likelihood_properties - } + "nloglik": False, + **BaseLikelihood._default_likelihood_properties, + } def __init__( self, @@ -130,7 +133,7 @@ def __init__( parameter_names=parameter_names, likelihood_properties=likelihood_properties, ) - + def verify_covariance(self): if self.covariance.full_matrix is False: self.covariance.compute_full_matrix() @@ -148,9 +151,13 @@ 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 - parameter_values_dict["vmean"], covariance_sum) - return likelihood_function(vector - parameter_values_dict["vmean"], covariance_sum) + if self.likelihood_properties["nloglik"]: + return -likelihood_function( + vector - parameter_values_dict["vmean"], covariance_sum + ) + return likelihood_function( + vector - parameter_values_dict["vmean"], covariance_sum + ) class MultivariateGaussianLikelihoodInterpolate1D(BaseLikelihood): From 218beed98148026c5281f94f1d239ebf3e34cc4a Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 8 Apr 2024 13:24:05 +0200 Subject: [PATCH 4/8] renaming --- .../{fit_velocity_old.ipynb => fit_velocity_test_bastien.ipynb} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename notebook/{fit_velocity_old.ipynb => fit_velocity_test_bastien.ipynb} (100%) diff --git a/notebook/fit_velocity_old.ipynb b/notebook/fit_velocity_test_bastien.ipynb similarity index 100% rename from notebook/fit_velocity_old.ipynb rename to notebook/fit_velocity_test_bastien.ipynb From 9be897b6650b361bbd66b3f7c2331e6a792feed9 Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 8 Apr 2024 13:36:33 +0200 Subject: [PATCH 5/8] switching vmean to the computation of vectors --- flip/likelihood.py | 8 ++------ flip/vectors.py | 3 +++ 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/flip/likelihood.py b/flip/likelihood.py index 4976832..d016d00 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -152,12 +152,8 @@ def __call__(self, parameter_values): f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) if self.likelihood_properties["nloglik"]: - return -likelihood_function( - vector - parameter_values_dict["vmean"], covariance_sum - ) - return likelihood_function( - vector - parameter_values_dict["vmean"], covariance_sum - ) + return -likelihood_function(vector, covariance_sum) + return likelihood_function(vector, covariance_sum) class MultivariateGaussianLikelihoodInterpolate1D(BaseLikelihood): diff --git a/flip/vectors.py b/flip/vectors.py index 67d2ca5..6a919e6 100644 --- a/flip/vectors.py +++ b/flip/vectors.py @@ -77,6 +77,9 @@ def load_velocity_vector( f"""Please choose a velocity_type among {_avail_velocity_type}""" ) + if "vmean" in parameter_values_dict: + velocity = velocity - parameter_values_dict["vmean"] + return velocity, velocity_error From 515b7649707cf249520d8366d6260f58dc0a727e Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 8 Apr 2024 13:36:43 +0200 Subject: [PATCH 6/8] format --- flip/likelihood.py | 11 +++++++++++ flip/vectors.py | 4 +--- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/flip/likelihood.py b/flip/likelihood.py index d016d00..eb7c2ba 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -195,6 +195,17 @@ def __init__( self.interpolation_value_range = interpolation_value_range def verify_covariance(self): + """ + The verify_covariance function is used to ensure that the covariance matrix of each + parameter in the model has been computed. If it has not, then this function will compute + it and store it as a full matrix. + + Args: + self: Bind the method to the object + + Returns: + Nothing + """ for i in range(len(self.covariance)): if self.covariance[i].full_matrix is False: self.covariance[i].compute_full_matrix() diff --git a/flip/vectors.py b/flip/vectors.py index 6a919e6..23ef53e 100644 --- a/flip/vectors.py +++ b/flip/vectors.py @@ -128,9 +128,7 @@ def compute_observed_distance_modulus( mu = data["mb"] + alpha * data["x1"] - beta * data["c"] - M0 variance_mu = ( - data["e_mb"] ** 2 - + alpha**2 * data["e_x1"] ** 2 - + beta**2 * data["e_c"] ** 2 + data["e_mb"] ** 2 + alpha**2 * data["e_x1"] ** 2 + beta**2 * data["e_c"] ** 2 ) variance_mu += ( 2 * alpha * data["cov_mb_x1"] From 8249855dfe9ee310801327f948d7b5b495d12a5d Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 8 Apr 2024 13:41:27 +0200 Subject: [PATCH 7/8] switch negloglik to True by default and put it in basis likelihood --- flip/likelihood.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/flip/likelihood.py b/flip/likelihood.py index eb7c2ba..826b92d 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -27,6 +27,7 @@ class BaseLikelihood(object): "inversion_method": "inverse", "velocity_type": "direct", "velocity_estimator": "full", + "negloglik": True, } def __init__( @@ -114,18 +115,12 @@ 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={}, - negloglik=False, ): super(MultivariateGaussianLikelihood, self).__init__( covariance=covariance, @@ -151,7 +146,7 @@ def __call__(self, parameter_values): likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) - if self.likelihood_properties["nloglik"]: + if self.likelihood_properties["negloglik"]: return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum) @@ -249,6 +244,8 @@ def __call__(self, parameter_values): likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) + if self.likelihood_properties["negloglik"]: + return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum) @@ -348,4 +345,6 @@ def __call__( likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) + if self.likelihood_properties["negloglik"]: + return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum) From 362517adb44f0d2d0029aecd65ab7c519b1a602e Mon Sep 17 00:00:00 2001 From: corentinravoux Date: Mon, 8 Apr 2024 13:42:45 +0200 Subject: [PATCH 8/8] renaming for PEP8 standard --- flip/fitter.py | 4 ++-- flip/likelihood.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/flip/fitter.py b/flip/fitter.py index d01ee35..7976fcc 100644 --- a/flip/fitter.py +++ b/flip/fitter.py @@ -194,7 +194,7 @@ def init_from_covariance( likelihood = minuit_fitter.get_likelihood( parameter_dict, likelihood_type=likelihood_type, - likelihood_properties={**likelihood_properties, 'nloglik': True}, + likelihood_properties={**likelihood_properties, 'negative_log_likelihood': 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, 'nloglik': True}, + likelihood_properties={**likelihood_properties, 'negative_log_likelihood': True}, ) def setup_minuit(self, parameter_dict): diff --git a/flip/likelihood.py b/flip/likelihood.py index 826b92d..5fdc2fa 100644 --- a/flip/likelihood.py +++ b/flip/likelihood.py @@ -27,7 +27,7 @@ class BaseLikelihood(object): "inversion_method": "inverse", "velocity_type": "direct", "velocity_estimator": "full", - "negloglik": True, + "negative_log_likelihood": True, } def __init__( @@ -146,7 +146,7 @@ def __call__(self, parameter_values): likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) - if self.likelihood_properties["negloglik"]: + if self.likelihood_properties["negative_log_likelihood"]: return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum) @@ -244,7 +244,7 @@ def __call__(self, parameter_values): likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) - if self.likelihood_properties["negloglik"]: + if self.likelihood_properties["negative_log_likelihood"]: return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum) @@ -345,6 +345,6 @@ def __call__( likelihood_function = eval( f"log_likelihood_gaussian_{self.likelihood_properties['inversion_method']}" ) - if self.likelihood_properties["negloglik"]: + if self.likelihood_properties["negative_log_likelihood"]: return -likelihood_function(vector, covariance_sum) return likelihood_function(vector, covariance_sum)