From fd7bb08d53375c8d85770c4ad539704cb61749a8 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 3 Apr 2025 13:54:56 +0100 Subject: [PATCH 01/30] Add a fisher-information-matrix sampler --- bilby/core/__init__.py | 2 +- bilby/core/fisher.py | 38 ++++++---- bilby/core/sampler/fisher.py | 124 +++++++++++++++++++++++++++++++ examples/core_examples/fisher.py | 80 ++++++++++++++++++++ setup.py | 1 + 5 files changed, 230 insertions(+), 15 deletions(-) create mode 100644 bilby/core/sampler/fisher.py create mode 100644 examples/core_examples/fisher.py diff --git a/bilby/core/__init__.py b/bilby/core/__init__.py index 968f961d0..7446bd24f 100644 --- a/bilby/core/__init__.py +++ b/bilby/core/__init__.py @@ -1 +1 @@ -from . import grid, likelihood, prior, result, sampler, series, utils, fisher +from . import grid, likelihood, prior, result, sampler, series, utils diff --git a/bilby/core/fisher.py b/bilby/core/fisher.py index 11b9af386..a5cb571c8 100644 --- a/bilby/core/fisher.py +++ b/bilby/core/fisher.py @@ -3,10 +3,13 @@ import scipy.linalg from scipy.optimize import minimize +from .utils import random + class FisherMatrixPosteriorEstimator(object): - def __init__(self, likelihood, priors, parameters=None, fd_eps=1e-6, n_prior_samples=100): - """ A class to estimate posteriors using the Fisher Matrix approach + def __init__(self, likelihood, priors, parameters=None, minimization_method="Nelder-Mead", + fd_eps=1e-6, n_prior_samples=100): + """ A class to estimate posteriors using a Fisher Information Matrix approach Parameters ---------- @@ -16,6 +19,8 @@ def __init__(self, likelihood, priors, parameters=None, fd_eps=1e-6, n_prior_sam A bilby prior object parameters: list Names of parameters to sample in + minimization_method: str (Nelder-Mead) + The method to use in scipy.optimize.minimize fd_eps: float A parameter to control the size of perturbation used when finite differencing the likelihood @@ -28,6 +33,7 @@ def __init__(self, likelihood, priors, parameters=None, fd_eps=1e-6, n_prior_sam self.parameter_names = priors.non_fixed_keys else: self.parameter_names = parameters + self.minimization_method = minimization_method self.fd_eps = fd_eps self.n_prior_samples = n_prior_samples self.N = len(self.parameter_names) @@ -45,6 +51,11 @@ def __init__(self, likelihood, priors, parameters=None, fd_eps=1e-6, n_prior_sam self.prior_width_dict[key] = width def log_likelihood(self, sample): + if isinstance(sample, dict) is False: + if isinstance(sample, pd.DataFrame) and len(sample) == 1: + sample = sample.to_dict() + else: + raise ValueError() self.likelihood.parameters.update(sample) return self.likelihood.log_likelihood() @@ -60,14 +71,12 @@ def calculate_iFIM(self, sample): return iFIM def sample_array(self, sample, n=1): - from .utils.random import rng - if sample == "maxL": sample = self.get_maximum_likelihood_sample() self.mean = np.array(list(sample.values())) self.iFIM = self.calculate_iFIM(sample) - return rng.multivariate_normal(self.mean, self.iFIM, n) + return random.rng.multivariate_normal(self.mean, self.iFIM, n) def sample_dataframe(self, sample, n=1): samples = self.sample_array(sample, n) @@ -89,8 +98,8 @@ def get_second_order_derivative(self, sample, ii, jj): def get_finite_difference_xx(self, sample, ii): # Sample grid - p = self.shift_sample_x(sample, ii, 1) - m = self.shift_sample_x(sample, ii, -1) + p = self._shift_sample_x(sample, ii, 1) + m = self._shift_sample_x(sample, ii, -1) dx = .5 * (p[ii] - m[ii]) @@ -102,10 +111,10 @@ def get_finite_difference_xx(self, sample, ii): def get_finite_difference_xy(self, sample, ii, jj): # Sample grid - pp = self.shift_sample_xy(sample, ii, 1, jj, 1) - pm = self.shift_sample_xy(sample, ii, 1, jj, -1) - mp = self.shift_sample_xy(sample, ii, -1, jj, 1) - mm = self.shift_sample_xy(sample, ii, -1, jj, -1) + pp = self._shift_sample_xy(sample, ii, 1, jj, 1) + pm = self._shift_sample_xy(sample, ii, 1, jj, -1) + mp = self._shift_sample_xy(sample, ii, -1, jj, 1) + mm = self._shift_sample_xy(sample, ii, -1, jj, -1) dx = .5 * (pp[ii] - mm[ii]) dy = .5 * (pp[jj] - mm[jj]) @@ -117,7 +126,7 @@ def get_finite_difference_xy(self, sample, ii, jj): return (loglpp - loglpm - loglmp + loglmm) / (4 * dx * dy) - def shift_sample_x(self, sample, x_key, x_coef): + def _shift_sample_x(self, sample, x_key, x_coef): vx = sample[x_key] dvx = self.fd_eps * self.prior_width_dict[x_key] @@ -127,7 +136,7 @@ def shift_sample_x(self, sample, x_key, x_coef): return shift_sample - def shift_sample_xy(self, sample, x_key, x_coef, y_key, y_coef): + def _shift_sample_xy(self, sample, x_key, x_coef, y_key, y_coef): vx = sample[x_key] vy = sample[y_key] @@ -166,9 +175,10 @@ def neg_log_like(x, self, T=1): x0, args=(self, 1), bounds=self.prior_bounds, - method="L-BFGS-B", + method=self.minimization_method, ) if out.fun < minlogL: minout = out + self.minimization_metadata = minout return {key: val for key, val in zip(self.parameter_names, minout.x)} diff --git a/bilby/core/sampler/fisher.py b/bilby/core/sampler/fisher.py new file mode 100644 index 000000000..3196a9310 --- /dev/null +++ b/bilby/core/sampler/fisher.py @@ -0,0 +1,124 @@ +import datetime + +from ..fisher import FisherMatrixPosteriorEstimator +from .base_sampler import Sampler, signal_wrapper + + +class Fisher(Sampler): + """ + TBD + + Parameters + ========== + likelihood: likelihood.Likelihood + A object with a log_l method + priors: bilby.core.prior.PriorDict, dict + Priors to be used in the search. + This has attributes for each parameter to be sampled. + outdir: str, optional + Name of the output directory + label: str, optional + Naming scheme of the output files + use_ratio: bool, optional + Switch to set whether or not you want to use the log-likelihood ratio + or just the log-likelihood + + """ + + sampler_name = "fisher" + sampling_seed_key = "seed" + + def __init__( + self, + likelihood, + priors, + outdir="outdir", + label="label", + nsamples=1000, + minimization_method="Nelder-Mead", + n_prior_samples=100, + fd_eps=1e-6, + **kwargs, + ): + super(Fisher, self).__init__( + likelihood=likelihood, + priors=priors, + outdir=outdir, + label=label, + **kwargs, + ) + self.nsamples = nsamples + self.minimization_method = minimization_method + self.n_prior_samples = n_prior_samples + self.fd_eps = fd_eps + + @classmethod + def get_expected_outputs(cls, outdir=None, label=None): + """Get lists of the expected outputs directories and files. + + These are used by :code:`bilby_pipe` when transferring files via HTCondor. + + Parameters + ---------- + outdir : str + The output directory. + label : str + The label for the run. + + Returns + ------- + list + List of file names. + list + List of directory names. Will always be empty for dynesty. + """ + raise NotImplementedError() + + @property + def sampler_init(self): + raise NotImplementedError() + + @signal_wrapper + def run_sampler(self): + + self.start_time = datetime.datetime.now() + + fisher = FisherMatrixPosteriorEstimator( + likelihood=self.likelihood, + priors=self.priors, + minimization_method=self.minimization_method, + n_prior_samples=self.n_prior_samples, + fd_eps=self.fd_eps, + ) + + samples = fisher.sample_dataframe("maxL", self.nsamples) + logl = [ + fisher.log_likelihood(sample.to_dict()) for _, sample in samples.iterrows() + ] + + end_time = datetime.datetime.now() + self.sampling_time = end_time - self.start_time + + self._generate_result(samples, logl) + + return self.result + + def _generate_result(self, samples, log_likelihood_evaluations): + """ + Extract the information we need from the output. + + Parameters + ========== + out: dynesty.result.Result + The dynesty output. + """ + + self.result.samples = samples + self.result.log_likelihood_evaluations = log_likelihood_evaluations + self.result.sampling_time = self.sampling_time.total_seconds() + + self.result.meta_data["run_statistics"] = dict( + nlikelihood=None, + neffsamples=None, + sampling_time_s=self.sampling_time.total_seconds(), + ) diff --git a/examples/core_examples/fisher.py b/examples/core_examples/fisher.py new file mode 100644 index 000000000..cb2836474 --- /dev/null +++ b/examples/core_examples/fisher.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python +""" +TBD + +""" +import bilby +import numpy as np +from bilby.core.utils.random import rng, seed + +# Sets seed of bilby's generator "rng" to "123" to ensure reproducibility +seed(123) + +# A few simple setup steps +outdir = "outdir" +bilby.utils.check_directory_exists_and_if_not_mkdir(outdir) + + +# First, we define our "signal model", in this case a simple linear function +def model(time, m, c): + return time * m + c + + +# Now we define the injection parameters which we make simulated data with +injection_parameters = dict(m=0.5, c=0.2) + +# For this example, we'll use standard Gaussian noise + +# These lines of code generate the fake data. Note the ** just unpacks the +# contents of the injection_parameters when calling the model function. +sampling_frequency = 10 +time_duration = 10 +time = np.arange(0, time_duration, 1 / sampling_frequency) +N = len(time) +sigma = rng.normal(1, 0.01, N) +data = model(time, **injection_parameters) + rng.normal(0, sigma, N) + +# Now lets instantiate a version of our GaussianLikelihood, giving it +# the time, data and signal model +likelihood = bilby.likelihood.GaussianLikelihood(time, data, model, sigma) + +# From hereon, the syntax is exactly equivalent to other bilby examples +# We make a prior +priors = dict() +priors["m"] = bilby.core.prior.Uniform(0, 5, "m") +priors["c"] = bilby.core.prior.Uniform(-2, 2, "c") + +# Run the "fisher" sampler +result_fisher = bilby.run_sampler( + likelihood=likelihood, + priors=priors, + sampler="fisher", + injection_parameters=injection_parameters, + outdir=outdir, + label="example_fisher", + nsamples=5000, +) + + +# Run dynesty as a comparison +result_dynesty = bilby.run_sampler( + likelihood=likelihood, + priors=priors, + sampler="dynesty", + nlive=500, + injection_parameters=injection_parameters, + outdir=outdir, + label="example_dynesty", +) + +# Finally plot a corner plot: all outputs are stored in outdir +bilby.result.plot_multiple( + [result_fisher, result_dynesty], + filename=f"{outdir}/comparison_fisher_dynesty.png", + labels=[ + f"Fisher ({result_fisher.meta_data['run_statistics']['sampling_time_s']:0.2f})", + f"Dynesty ({result_dynesty.meta_data['run_statistics']['sampling_time_s']:0.2f})", + ], + parameters=injection_parameters, + truth_color="C3", +) diff --git a/setup.py b/setup.py index c9c00e2e6..b94123a6f 100644 --- a/setup.py +++ b/setup.py @@ -90,6 +90,7 @@ def readfile(filename): "bilby.dynesty=bilby.core.sampler.dynesty:Dynesty", "bilby.dynamic_dynesty=bilby.core.sampler.dynamic_dynesty:DynamicDynesty", "bilby.emcee=bilby.core.sampler.emcee:Emcee", + "bilby.fisher=bilby.core.sampler.fisher:Fisher", "bilby.kombine=bilby.core.sampler.kombine:Kombine", "bilby.nessai=bilby.core.sampler.nessai:Nessai", "bilby.nestle=bilby.core.sampler.nestle:Nestle", From 4dfa98f8f255c42e5cd3a01adc45266d1d553ef8 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 3 Apr 2025 13:57:29 +0100 Subject: [PATCH 02/30] Revert change --- bilby/core/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/bilby/core/__init__.py b/bilby/core/__init__.py index 7446bd24f..968f961d0 100644 --- a/bilby/core/__init__.py +++ b/bilby/core/__init__.py @@ -1 +1 @@ -from . import grid, likelihood, prior, result, sampler, series, utils +from . import grid, likelihood, prior, result, sampler, series, utils, fisher From 8cf5808edbfbc860b80ceb792e1e5a08a4c846eb Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 3 Apr 2025 13:58:28 +0100 Subject: [PATCH 03/30] Remove old example --- .../linear_regression_with_Fisher.py | 74 ------------------- 1 file changed, 74 deletions(-) delete mode 100644 examples/core_examples/linear_regression_with_Fisher.py diff --git a/examples/core_examples/linear_regression_with_Fisher.py b/examples/core_examples/linear_regression_with_Fisher.py deleted file mode 100644 index ae7f74dd7..000000000 --- a/examples/core_examples/linear_regression_with_Fisher.py +++ /dev/null @@ -1,74 +0,0 @@ -#!/usr/bin/env python -""" -An example of how to use bilby to perform parameter estimation for -non-gravitational wave data. In this case, fitting a linear function to -data with background Gaussian noise. We then compare the result to posteriors -estimated using the Fisher Information Matrix approximation. - -""" -import copy - -import bilby -from bilby.core.utils.random import rng, seed - -# sets seed of bilby's generator "rng" to "123" -seed(123) - -import numpy as np - -# A few simple setup steps -outdir = "outdir" - - -# First, we define our "signal model", in this case a simple linear function -def model(time, m, c): - return time * m + c - - -# Now we define the injection parameters which we make simulated data with -injection_parameters = dict(m=0.5, c=0.2) - -# For this example, we'll use standard Gaussian noise - -# These lines of code generate the fake data. Note the ** just unpacks the -# contents of the injection_parameters when calling the model function. -sampling_frequency = 10 -time_duration = 10 -time = np.arange(0, time_duration, 1 / sampling_frequency) -N = len(time) -sigma = rng.normal(1, 0.01, N) -data = model(time, **injection_parameters) + rng.normal(0, sigma, N) - -# Now lets instantiate a version of our GaussianLikelihood, giving it -# the time, data and signal model -likelihood = bilby.likelihood.GaussianLikelihood(time, data, model, sigma) - -# From hereon, the syntax is exactly equivalent to other bilby examples -# We make a prior -priors = dict() -priors["m"] = bilby.core.prior.Uniform(0, 5, "m") -priors["c"] = bilby.core.prior.Uniform(-2, 2, "c") -priors = bilby.core.prior.PriorDict(priors) - -# And run sampler -result = bilby.run_sampler( - likelihood=likelihood, - priors=priors, - sampler="dynesty", - nlive=1000, - injection_parameters=injection_parameters, - outdir=outdir, - label="Nested Sampling", -) - -# Finally plot a corner plot: all outputs are stored in outdir -result.plot_corner() - -fim = bilby.core.fisher.FisherMatrixPosteriorEstimator(likelihood, priors) -result_fim = copy.deepcopy(result) -result_fim.posterior = fim.sample_dataframe("maxL", 10000) -result_fim.label = "Fisher" - -bilby.core.result.plot_multiple( - [result, result_fim], parameters=injection_parameters, truth_color="k" -) From e5064cc37a93d51f6b1fc5c996b16edf680819ea Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 3 Apr 2025 14:04:10 +0100 Subject: [PATCH 04/30] Improve the example --- examples/core_examples/fisher.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/examples/core_examples/fisher.py b/examples/core_examples/fisher.py index cb2836474..5b095ab06 100644 --- a/examples/core_examples/fisher.py +++ b/examples/core_examples/fisher.py @@ -1,6 +1,7 @@ #!/usr/bin/env python """ -TBD +This example, building on the Gaussian example, demonstrates the use of the Fisher +tools within Bilby. """ import bilby @@ -67,7 +68,7 @@ def model(time, m, c): label="example_dynesty", ) -# Finally plot a corner plot: all outputs are stored in outdir +# Plot a corner plot: all outputs are stored in outdir bilby.result.plot_multiple( [result_fisher, result_dynesty], filename=f"{outdir}/comparison_fisher_dynesty.png", @@ -78,3 +79,9 @@ def model(time, m, c): parameters=injection_parameters, truth_color="C3", ) + +# Note that the `fisher` tools can also be accessed directly +fisher = bilby.core.fisher.FisherMatrixPosteriorEstimator(likelihood, priors) +samples = fisher.sample_dataframe( + "maxL", 1000 +) # Draw a set of samples as a pandas dataframe From 69470fd9e34039ef00c3b7a36d18f0ad0237fdd2 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 3 Apr 2025 14:19:24 +0100 Subject: [PATCH 05/30] Improve the docstring --- bilby/core/sampler/fisher.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/bilby/core/sampler/fisher.py b/bilby/core/sampler/fisher.py index 3196a9310..8d85eaf4b 100644 --- a/bilby/core/sampler/fisher.py +++ b/bilby/core/sampler/fisher.py @@ -6,7 +6,8 @@ class Fisher(Sampler): """ - TBD + A sampler class that estimates the maximum likelihood using scipy, then draws + posterior samples using the Fisher information matrix. Parameters ========== @@ -19,10 +20,16 @@ class Fisher(Sampler): Name of the output directory label: str, optional Naming scheme of the output files - use_ratio: bool, optional - Switch to set whether or not you want to use the log-likelihood ratio - or just the log-likelihood - + nsamples: int + The number of samples to draw in the posterior + minimization_method: str (Nelder-Mead) + The method to use in scipy.optimize.minimize + fd_eps: float + A parameter to control the size of perturbation used when finite + differencing the likelihood + n_prior_samples: int + The number of prior samples to draw and use to attempt estimatation + of the maximum likelihood sample. """ sampler_name = "fisher" From b44f532538a332451bfab65532941fdd8f630e76 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 3 Apr 2025 14:24:47 +0100 Subject: [PATCH 06/30] Improve example: add units --- examples/core_examples/fisher.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/core_examples/fisher.py b/examples/core_examples/fisher.py index 5b095ab06..179064330 100644 --- a/examples/core_examples/fisher.py +++ b/examples/core_examples/fisher.py @@ -73,8 +73,8 @@ def model(time, m, c): [result_fisher, result_dynesty], filename=f"{outdir}/comparison_fisher_dynesty.png", labels=[ - f"Fisher ({result_fisher.meta_data['run_statistics']['sampling_time_s']:0.2f})", - f"Dynesty ({result_dynesty.meta_data['run_statistics']['sampling_time_s']:0.2f})", + f"Fisher ({result_fisher.meta_data['run_statistics']['sampling_time_s']:0.2f} s)", + f"Dynesty ({result_dynesty.meta_data['run_statistics']['sampling_time_s']:0.2f} s)", ], parameters=injection_parameters, truth_color="C3", From e79639fbf44174a69eda1bcbad3dc17ba8be28d3 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Mon, 14 Apr 2025 21:37:53 +0100 Subject: [PATCH 07/30] Improvemens to FIM and the sampler - Add use of scipy.differentiate for Hessian computation - Fix bugs in sampler implementation --- bilby/bilby_mcmc/proposals.py | 2 +- bilby/core/__init__.py | 2 +- bilby/core/{fisher.py => fisher_matrix.py} | 36 ++- bilby/core/fisher_sampler.py | 260 ++++++++++++++++++ bilby/core/result.py | 9 +- bilby/core/sampler/fisher.py | 131 --------- .../{fisher.py => fisher_example.py} | 10 +- setup.py | 2 +- 8 files changed, 306 insertions(+), 146 deletions(-) rename bilby/core/{fisher.py => fisher_matrix.py} (84%) create mode 100644 bilby/core/fisher_sampler.py delete mode 100644 bilby/core/sampler/fisher.py rename examples/core_examples/{fisher.py => fisher_example.py} (92%) diff --git a/bilby/bilby_mcmc/proposals.py b/bilby/bilby_mcmc/proposals.py index 6100d75f8..90db8b92d 100644 --- a/bilby/bilby_mcmc/proposals.py +++ b/bilby/bilby_mcmc/proposals.py @@ -6,7 +6,7 @@ from scipy.spatial.distance import jensenshannon from scipy.stats import gaussian_kde -from ..core.fisher import FisherMatrixPosteriorEstimator +from ..core.fisher_matrix import FisherMatrixPosteriorEstimator from ..core.prior import PriorDict from ..core.sampler.base_sampler import SamplerError from ..core.utils import logger, random, reflect diff --git a/bilby/core/__init__.py b/bilby/core/__init__.py index 968f961d0..5d24080fd 100644 --- a/bilby/core/__init__.py +++ b/bilby/core/__init__.py @@ -1 +1 @@ -from . import grid, likelihood, prior, result, sampler, series, utils, fisher +from . import grid, likelihood, prior, result, sampler, series, utils, fisher_matrix, fisher_sampler diff --git a/bilby/core/fisher.py b/bilby/core/fisher_matrix.py similarity index 84% rename from bilby/core/fisher.py rename to bilby/core/fisher_matrix.py index a5cb571c8..7f87f315c 100644 --- a/bilby/core/fisher.py +++ b/bilby/core/fisher_matrix.py @@ -1,9 +1,12 @@ +from packaging import version + import numpy as np import pandas as pd import scipy.linalg from scipy.optimize import minimize -from .utils import random +from .utils import random, logger +from .prior import PriorDict class FisherMatrixPosteriorEstimator(object): @@ -29,6 +32,10 @@ def __init__(self, likelihood, priors, parameters=None, minimization_method="Nel of the maximum likelihood sample. """ self.likelihood = likelihood + + if isinstance(priors, PriorDict) is False: + priors = PriorDict(priors) + if parameters is None: self.parameter_names = priors.non_fixed_keys else: @@ -83,10 +90,28 @@ def sample_dataframe(self, sample, n=1): return pd.DataFrame(samples, columns=self.parameter_names) def calculate_FIM(self, sample): - FIM = np.zeros((self.N, self.N)) - for ii, ii_key in enumerate(self.parameter_names): - for jj, jj_key in enumerate(self.parameter_names): - FIM[ii, jj] = -self.get_second_order_derivative(sample, ii_key, jj_key) + if version.parse(scipy.__version__) < version.parse("1.15"): + logger.info("Scipy version < 1.15, using fallback") + FIM = np.zeros((self.N, self.N)) + for ii, ii_key in enumerate(self.parameter_names): + for jj, jj_key in enumerate(self.parameter_names): + FIM[ii, jj] = -self.get_second_order_derivative(sample, ii_key, jj_key) + return FIM + + import scipy.differentiate as sd + + def array_to_dict(x_array): + return dict(zip(self.parameter_names, x_array)) + + def wrapped_logl(x_array): + return self.log_likelihood(array_to_dict(x_array)) + + def wrapped_logl_arb(x_array): + return np.apply_along_axis(wrapped_logl, 0, x_array) + + point = np.array([sample[key] for key in self.parameter_names]) + res = sd.hessian(wrapped_logl_arb, point) + FIM = - res.ddf return FIM @@ -181,4 +206,5 @@ def neg_log_like(x, self, T=1): minout = out self.minimization_metadata = minout + logger.info(f"Maximum likelihood estimation: {minout.message}") return {key: val for key, val in zip(self.parameter_names, minout.x)} diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py new file mode 100644 index 000000000..c3bf6902f --- /dev/null +++ b/bilby/core/fisher_sampler.py @@ -0,0 +1,260 @@ +import datetime + +import numpy as np +from scipy.stats import multivariate_normal + +from .fisher_matrix import FisherMatrixPosteriorEstimator +from .sampler.base_sampler import Sampler, signal_wrapper +from .result import rejection_sample +from .utils import logger, kish_log_effective_sample_size, safe_save_figure + + +class Fisher(Sampler): + """ + A sampler class that estimates the maximum likelihood using scipy, then draws + posterior samples using the Fisher information matrix. + + Parameters + ========== + likelihood: likelihood.Likelihood + A object with a log_l method + priors: bilby.core.prior.PriorDict, dict + Priors to be used in the search. + This has attributes for each parameter to be sampled. + outdir: str, optional + Name of the output directory + label: str, optional + Naming scheme of the output files + rejection_sampling: bool + If true, utilise rejection sampling to reweight from the Fisher matrix + Gaussian posterior approximation to the true posterior. Note: if true, + the number of samples is determined by the efficiency which can be + low if the maximum likelihood is not found or the posterior is highly + non-Gaussian. + nsamples: int + The target number of samples to draw in the posterior + minimization_method: str (Nelder-Mead) + The method to use in scipy.optimize.minimize + fd_eps: float + A parameter to control the size of perturbation used when finite + differencing the likelihood + n_prior_samples: int + The number of prior samples to draw and use to attempt estimatation + of the maximum likelihood sample. + """ + + sampler_name = "fisher" + sampling_seed_key = "seed" + default_kwargs = dict( + rejection_sampling=True, + nsamples=1000, + minimization_method="Nelder-Mead", + n_prior_samples=100, + fd_eps=1e-6, + ) + + def __init__( + self, + likelihood, + priors, + outdir="outdir", + label="label", + use_ratio=False, + plot=False, + exit_code=77, + skip_import_verification=True, + **kwargs, + ): + super(Fisher, self).__init__( + likelihood=likelihood, + priors=priors, + outdir=outdir, + label=label, + use_ratio=use_ratio, + plot=plot, + skip_import_verification=skip_import_verification, + exit_code=exit_code, + **kwargs, + ) + + @classmethod + def get_expected_outputs(cls, outdir=None, label=None): + """Get lists of the expected outputs directories and files. + + These are used by :code:`bilby_pipe` when transferring files via HTCondor. + + Parameters + ---------- + outdir : str + The output directory. + label : str + The label for the run. + + Returns + ------- + list + List of file names. + list + List of directory names. Will always be empty for dynesty. + """ + raise NotImplementedError() + + @signal_wrapper + def run_sampler(self): + + self.start_time = datetime.datetime.now() + + fisher_mpe = FisherMatrixPosteriorEstimator( + likelihood=self.likelihood, + priors=self.priors, + minimization_method=self.kwargs["minimization_method"], + n_prior_samples=self.kwargs["n_prior_samples"], + fd_eps=self.kwargs["fd_eps"], + ) + + raw_samples = fisher_mpe.sample_dataframe("maxL", self.kwargs["nsamples"]) + raw_logl = np.array( + [ + fisher_mpe.log_likelihood(sample.to_dict()) + for _, sample in raw_samples.iterrows() + ] + ) + + if self.kwargs["rejection_sampling"]: + logl_norm = multivariate_normal.logpdf( + raw_samples, mean=fisher_mpe.mean, cov=fisher_mpe.iFIM + ) + ln_weights = raw_logl - logl_norm + ln_weights -= np.max(ln_weights) + weights = np.exp(ln_weights) + + samples, idxs = rejection_sample(raw_samples, weights, return_idxs=True) + logl = raw_logl[idxs] + + nsamples = len(samples) + efficiency = 100 * nsamples / len(raw_samples) + ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) + logger.info( + f"Rejection sampling Fisher posterior produced {nsamples}" + f" with an efficiency of {efficiency}% and effective sample" + f" size {ess}" + ) + if self.plot: + self.create_rejection_sample_diagnostic(samples, raw_samples) + else: + samples = raw_samples + logl = raw_logl + + end_time = datetime.datetime.now() + self.sampling_time = end_time - self.start_time + + self._generate_result(samples, logl) + + return self.result + + def create_rejection_sample_diagnostic(self, samples, raw_samples): + import corner + import matplotlib.pyplot as plt + import matplotlib.lines as mpllines + + kwargs = dict( + bins=50, + smooth=0.9, + max_n_ticks=3, + ) + + kwargs["labels"] = [ + label.replace("_", " ") for label in self.search_parameter_keys + ] + + # Create the data array to plot and pass everything to corner + xs = samples[self.search_parameter_keys].values + rxs = raw_samples[self.search_parameter_keys].values + if len(self.search_parameter_keys) > 1: + lines = [] + ls = "-" + c = "C0" + fig = corner.corner( + rxs, + color=c, + contour_kwargs={"linestyles": ls, "alpha": 0.8}, + hist_kwargs={"ls": ls, "alpha": 0.8}, + data_kwargs={"alpha": 1}, + no_fill_contours=True, + alpha=0.8, + plot_density=False, + plot_datapoints=True, + fill_contours=False, + quantiles=[0.16, 0.84], + levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), + **kwargs + ) + lines.append(mpllines.Line2D([0], [0], color=c, linestyle=ls)) + + ls = "--" + c = "C1" + fig = corner.corner( + xs, + color=c, + contour_kwargs={"linestyles": ls, "alpha": 0.8}, + contourf_kwargs={"alpha": 0.8}, + hist_kwargs={"ls": ls, "alpha": 0.8}, + no_fill_contours=True, + fig=fig, + alpha=0.1, + plot_density=True, + plot_datapoints=False, + fill_contours=False, + quantiles=[0.16, 0.84], + levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), + **kwargs, + ) + lines.append(mpllines.Line2D([0], [0], color=c, linestyle=ls)) + + axes = fig.get_axes() + ndim = int(np.sqrt(len(axes))) + axes[ndim - 1].legend(lines, ["$g(x)$", "$f(x)$"]) + + else: + fig, ax = plt.subplots() + ax.hist( + xs, + bins=kwargs["bins"], + color="C0", + histtype="step", + ls="-", + ) + ax.hist( + rxs, + bins=kwargs["bins"], + color="C1", + histtype="step", + ls="--", + ) + ax.set_xlabel(kwargs["labels"][0]) + + filename = "{}/{}_rejection_sample.png".format(self.outdir, self.label) + logger.debug("Saving rejection-sample diagnopstic plot to {}".format(filename)) + safe_save_figure(fig=fig, filename=filename, dpi=400) + plt.close(fig) + + return fig + + def _generate_result(self, samples, log_likelihood_evaluations): + """ + Extract the information we need from the output. + + Parameters + ========== + out: dynesty.result.Result + The dynesty output. + """ + + self.result.samples = samples + self.result.log_likelihood_evaluations = log_likelihood_evaluations + + self.result.meta_data["run_statistics"] = dict( + nlikelihood=None, + neffsamples=None, + sampling_time_s=self.sampling_time.total_seconds(), + ) diff --git a/bilby/core/result.py b/bilby/core/result.py index 3569bc1da..8abc42078 100644 --- a/bilby/core/result.py +++ b/bilby/core/result.py @@ -264,7 +264,7 @@ def eval_pool(this_logl): return ln_weights, new_log_likelihood_array, new_log_prior_array, old_log_likelihood_array, old_log_prior_array -def rejection_sample(posterior, weights): +def rejection_sample(posterior, weights, return_idxs=False): """ Perform rejection sampling on a posterior using weights Parameters @@ -273,6 +273,8 @@ def rejection_sample(posterior, weights): The dataframe or array containing posterior samples weights: np.ndarray An array of weights + return_idxs: bool + If true, also return the boolean indexes Returns ======= @@ -281,7 +283,10 @@ def rejection_sample(posterior, weights): """ keep = weights > random.rng.uniform(0, max(weights), weights.shape) - return posterior[keep] + if return_idxs: + return posterior[keep], keep + else: + return posterior[keep] def reweight(result, label=None, new_likelihood=None, new_prior=None, diff --git a/bilby/core/sampler/fisher.py b/bilby/core/sampler/fisher.py deleted file mode 100644 index 8d85eaf4b..000000000 --- a/bilby/core/sampler/fisher.py +++ /dev/null @@ -1,131 +0,0 @@ -import datetime - -from ..fisher import FisherMatrixPosteriorEstimator -from .base_sampler import Sampler, signal_wrapper - - -class Fisher(Sampler): - """ - A sampler class that estimates the maximum likelihood using scipy, then draws - posterior samples using the Fisher information matrix. - - Parameters - ========== - likelihood: likelihood.Likelihood - A object with a log_l method - priors: bilby.core.prior.PriorDict, dict - Priors to be used in the search. - This has attributes for each parameter to be sampled. - outdir: str, optional - Name of the output directory - label: str, optional - Naming scheme of the output files - nsamples: int - The number of samples to draw in the posterior - minimization_method: str (Nelder-Mead) - The method to use in scipy.optimize.minimize - fd_eps: float - A parameter to control the size of perturbation used when finite - differencing the likelihood - n_prior_samples: int - The number of prior samples to draw and use to attempt estimatation - of the maximum likelihood sample. - """ - - sampler_name = "fisher" - sampling_seed_key = "seed" - - def __init__( - self, - likelihood, - priors, - outdir="outdir", - label="label", - nsamples=1000, - minimization_method="Nelder-Mead", - n_prior_samples=100, - fd_eps=1e-6, - **kwargs, - ): - super(Fisher, self).__init__( - likelihood=likelihood, - priors=priors, - outdir=outdir, - label=label, - **kwargs, - ) - self.nsamples = nsamples - self.minimization_method = minimization_method - self.n_prior_samples = n_prior_samples - self.fd_eps = fd_eps - - @classmethod - def get_expected_outputs(cls, outdir=None, label=None): - """Get lists of the expected outputs directories and files. - - These are used by :code:`bilby_pipe` when transferring files via HTCondor. - - Parameters - ---------- - outdir : str - The output directory. - label : str - The label for the run. - - Returns - ------- - list - List of file names. - list - List of directory names. Will always be empty for dynesty. - """ - raise NotImplementedError() - - @property - def sampler_init(self): - raise NotImplementedError() - - @signal_wrapper - def run_sampler(self): - - self.start_time = datetime.datetime.now() - - fisher = FisherMatrixPosteriorEstimator( - likelihood=self.likelihood, - priors=self.priors, - minimization_method=self.minimization_method, - n_prior_samples=self.n_prior_samples, - fd_eps=self.fd_eps, - ) - - samples = fisher.sample_dataframe("maxL", self.nsamples) - logl = [ - fisher.log_likelihood(sample.to_dict()) for _, sample in samples.iterrows() - ] - - end_time = datetime.datetime.now() - self.sampling_time = end_time - self.start_time - - self._generate_result(samples, logl) - - return self.result - - def _generate_result(self, samples, log_likelihood_evaluations): - """ - Extract the information we need from the output. - - Parameters - ========== - out: dynesty.result.Result - The dynesty output. - """ - - self.result.samples = samples - self.result.log_likelihood_evaluations = log_likelihood_evaluations - self.result.sampling_time = self.sampling_time.total_seconds() - - self.result.meta_data["run_statistics"] = dict( - nlikelihood=None, - neffsamples=None, - sampling_time_s=self.sampling_time.total_seconds(), - ) diff --git a/examples/core_examples/fisher.py b/examples/core_examples/fisher_example.py similarity index 92% rename from examples/core_examples/fisher.py rename to examples/core_examples/fisher_example.py index 179064330..7f51e078c 100644 --- a/examples/core_examples/fisher.py +++ b/examples/core_examples/fisher_example.py @@ -53,11 +53,11 @@ def model(time, m, c): injection_parameters=injection_parameters, outdir=outdir, label="example_fisher", - nsamples=5000, + nsamples=10000, + clean=True, + plot=True, ) - -# Run dynesty as a comparison result_dynesty = bilby.run_sampler( likelihood=likelihood, priors=priors, @@ -80,8 +80,8 @@ def model(time, m, c): truth_color="C3", ) -# Note that the `fisher` tools can also be accessed directly -fisher = bilby.core.fisher.FisherMatrixPosteriorEstimator(likelihood, priors) +# Note that the `fisher` tools can also be accessed directly, e.g: +fisher = bilby.core.fisher_matrix.FisherMatrixPosteriorEstimator(likelihood, priors) samples = fisher.sample_dataframe( "maxL", 1000 ) # Draw a set of samples as a pandas dataframe diff --git a/setup.py b/setup.py index b94123a6f..85c67d24b 100644 --- a/setup.py +++ b/setup.py @@ -90,7 +90,7 @@ def readfile(filename): "bilby.dynesty=bilby.core.sampler.dynesty:Dynesty", "bilby.dynamic_dynesty=bilby.core.sampler.dynamic_dynesty:DynamicDynesty", "bilby.emcee=bilby.core.sampler.emcee:Emcee", - "bilby.fisher=bilby.core.sampler.fisher:Fisher", + "bilby.fisher=bilby.core.fisher_sampler:Fisher", "bilby.kombine=bilby.core.sampler.kombine:Kombine", "bilby.nessai=bilby.core.sampler.nessai:Nessai", "bilby.nestle=bilby.core.sampler.nestle:Nestle", From bdb59f4e9aca3d1d60bf3e7905bbb7b224ce0207 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Tue, 15 Apr 2025 21:27:29 +0100 Subject: [PATCH 08/30] Clean up of the method --- bilby/core/fisher_matrix.py | 60 +++++++++++++++++++++++++----------- bilby/core/fisher_sampler.py | 18 +++++++---- 2 files changed, 54 insertions(+), 24 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index 7f87f315c..02898e782 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -4,6 +4,7 @@ import pandas as pd import scipy.linalg from scipy.optimize import minimize +import tqdm from .utils import random, logger from .prior import PriorDict @@ -45,6 +46,7 @@ def __init__(self, likelihood, priors, parameters=None, minimization_method="Nel self.n_prior_samples = n_prior_samples self.N = len(self.parameter_names) + # Construct prior samples at initialisation so that the prior is not stored self.prior_samples = [ priors.sample_subset(self.parameter_names) for _ in range(n_prior_samples) ] @@ -174,6 +176,22 @@ def _shift_sample_xy(self, sample, x_key, x_coef, y_key, y_coef): shift_sample[y_key] = vy + y_coef * dvy return shift_sample + def _maximize_likelihood_from_initial_sample(self, initial_sample): + x0 = list(initial_sample.values()) + + def neg_log_like(x, self, T=1): + sample = {key: val for key, val in zip(self.parameter_names, x)} + return - 1 / T * self.log_likelihood(sample) + + out = minimize( + neg_log_like, + x0, + args=(self, 1), + bounds=self.prior_bounds, + method=self.minimization_method, + ) + return out + def get_maximum_likelihood_sample(self, initial_sample=None): """ A method to attempt optimization of the maximum likelihood @@ -185,25 +203,31 @@ def get_maximum_likelihood_sample(self, initial_sample=None): is large or the posterior is wide relative to the prior, the method fails to find the global maximum in high dimensional problems. """ - minlogL = np.inf - for i in range(self.n_prior_samples): - initial_sample = self.prior_samples[i] - - x0 = list(initial_sample.values()) - - def neg_log_like(x, self, T=1): - sample = {key: val for key, val in zip(self.parameter_names, x)} - return - 1 / T * self.log_likelihood(sample) - - out = minimize( - neg_log_like, - x0, - args=(self, 1), - bounds=self.prior_bounds, - method=self.minimization_method, + + if initial_sample: + out = self._maximize_likelihood_from_initial_sample(initial_sample) + else: + logger.info(f"Maximising the likelihood using {self.n_prior_samples} prior samples") + max_logL = -np.inf + logL_list = [] + successes = 0 + for sample in tqdm.tqdm(self.prior_samples): + out = self._maximize_likelihood_from_initial_sample(sample) + logL_list.append(-out.fun) + if out.success: + successes += 1 + if out.fun > max_logL: + max_logL = -out.fun + minout = out + + if np.isinf(max_logL): + raise ValueError("Maxisation of the likelihood failed") + + logger.info( + f"Finished with {100 * successes / self.n_prior_samples}% success rate " + f"Maximum log-likelihood {max_logL} " + f"Distribution range {np.min(logL_list)}:{np.max(logL_list)}" ) - if out.fun < minlogL: - minout = out self.minimization_metadata = minout logger.info(f"Maximum likelihood estimation: {minout.message}") diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index c3bf6902f..799156451 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -121,11 +121,13 @@ def run_sampler(self): ) if self.kwargs["rejection_sampling"]: + logger.info(f"Rejection sampling the posterior from {self.kwargs['nsamples']} generated samples") logl_norm = multivariate_normal.logpdf( raw_samples, mean=fisher_mpe.mean, cov=fisher_mpe.iFIM ) + # TODO: Add prior weights ln_weights = raw_logl - logl_norm - ln_weights -= np.max(ln_weights) + ln_weights -= np.mean(ln_weights) weights = np.exp(ln_weights) samples, idxs = rejection_sample(raw_samples, weights, return_idxs=True) @@ -135,12 +137,12 @@ def run_sampler(self): efficiency = 100 * nsamples / len(raw_samples) ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) logger.info( - f"Rejection sampling Fisher posterior produced {nsamples}" + f"Rejection sampling Fisher posterior produced {nsamples} samples" f" with an efficiency of {efficiency}% and effective sample" f" size {ess}" ) if self.plot: - self.create_rejection_sample_diagnostic(samples, raw_samples) + self.create_rejection_sample_diagnostic(samples, raw_samples, fisher_mpe.mean) else: samples = raw_samples logl = raw_logl @@ -152,7 +154,7 @@ def run_sampler(self): return self.result - def create_rejection_sample_diagnostic(self, samples, raw_samples): + def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL): import corner import matplotlib.pyplot as plt import matplotlib.lines as mpllines @@ -161,6 +163,8 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples): bins=50, smooth=0.9, max_n_ticks=3, + truths=maxL, + truth_color="C3", ) kwargs["labels"] = [ @@ -178,7 +182,7 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples): rxs, color=c, contour_kwargs={"linestyles": ls, "alpha": 0.8}, - hist_kwargs={"ls": ls, "alpha": 0.8}, + hist_kwargs={"density": True, "ls": ls, "alpha": 0.8}, data_kwargs={"alpha": 1}, no_fill_contours=True, alpha=0.8, @@ -198,7 +202,7 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples): color=c, contour_kwargs={"linestyles": ls, "alpha": 0.8}, contourf_kwargs={"alpha": 0.8}, - hist_kwargs={"ls": ls, "alpha": 0.8}, + hist_kwargs={"density": True, "ls": ls, "alpha": 0.8}, no_fill_contours=True, fig=fig, alpha=0.1, @@ -223,6 +227,7 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples): color="C0", histtype="step", ls="-", + density=True, ) ax.hist( rxs, @@ -230,6 +235,7 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples): color="C1", histtype="step", ls="--", + density=True, ) ax.set_xlabel(kwargs["labels"][0]) From 40ff4428b2899869da1fda012f6499ce8c67da3d Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 16 Apr 2025 15:25:10 +0100 Subject: [PATCH 09/30] Clean up implementation --- bilby/core/fisher_matrix.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index 02898e782..66228d839 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -10,6 +10,10 @@ from .prior import PriorDict +def array_to_dict(array, keys): + return dict(zip(keys, array)) + + class FisherMatrixPosteriorEstimator(object): def __init__(self, likelihood, priors, parameters=None, minimization_method="Nelder-Mead", fd_eps=1e-6, n_prior_samples=100): @@ -99,23 +103,22 @@ def calculate_FIM(self, sample): for jj, jj_key in enumerate(self.parameter_names): FIM[ii, jj] = -self.get_second_order_derivative(sample, ii_key, jj_key) return FIM + else: + import scipy.differentiate as sd + point = np.array([sample[key] for key in self.parameter_names]) + res = sd.hessian(self.log_likelihood_from_array, point) + FIM = - res.ddf - import scipy.differentiate as sd - - def array_to_dict(x_array): - return dict(zip(self.parameter_names, x_array)) + return FIM + def log_likelihood_from_array(self, x_array): def wrapped_logl(x_array): - return self.log_likelihood(array_to_dict(x_array)) + return self.log_likelihood(array_to_dict(x_array, self.parameter_names)) def wrapped_logl_arb(x_array): return np.apply_along_axis(wrapped_logl, 0, x_array) - point = np.array([sample[key] for key in self.parameter_names]) - res = sd.hessian(wrapped_logl_arb, point) - FIM = - res.ddf - - return FIM + return wrapped_logl_arb(x_array) def get_second_order_derivative(self, sample, ii, jj): if ii == jj: From d3b78ab4e3c74d8deaf01630df9b338b73fe4c89 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 16 Apr 2025 16:09:14 +0100 Subject: [PATCH 10/30] Further clean up --- bilby/core/fisher_matrix.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index 66228d839..8a1279fce 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -10,7 +10,7 @@ from .prior import PriorDict -def array_to_dict(array, keys): +def array_to_dict(keys, array): return dict(zip(keys, array)) @@ -105,6 +105,7 @@ def calculate_FIM(self, sample): return FIM else: import scipy.differentiate as sd + logger.info("Using Scipy differentiate") point = np.array([sample[key] for key in self.parameter_names]) res = sd.hessian(self.log_likelihood_from_array, point) FIM = - res.ddf @@ -113,7 +114,7 @@ def calculate_FIM(self, sample): def log_likelihood_from_array(self, x_array): def wrapped_logl(x_array): - return self.log_likelihood(array_to_dict(x_array, self.parameter_names)) + return self.log_likelihood(array_to_dict(self.parameter_names, x_array)) def wrapped_logl_arb(x_array): return np.apply_along_axis(wrapped_logl, 0, x_array) @@ -182,14 +183,12 @@ def _shift_sample_xy(self, sample, x_key, x_coef, y_key, y_coef): def _maximize_likelihood_from_initial_sample(self, initial_sample): x0 = list(initial_sample.values()) - def neg_log_like(x, self, T=1): - sample = {key: val for key, val in zip(self.parameter_names, x)} - return - 1 / T * self.log_likelihood(sample) + def neg_log_like(x): + return - self.log_likelihood_from_array(x) out = minimize( neg_log_like, x0, - args=(self, 1), bounds=self.prior_bounds, method=self.minimization_method, ) @@ -216,20 +215,21 @@ def get_maximum_likelihood_sample(self, initial_sample=None): successes = 0 for sample in tqdm.tqdm(self.prior_samples): out = self._maximize_likelihood_from_initial_sample(sample) - logL_list.append(-out.fun) + logL = -out.fun + logL_list.append(logL) if out.success: successes += 1 - if out.fun > max_logL: - max_logL = -out.fun + if logL > max_logL: + max_logL = logL minout = out if np.isinf(max_logL): raise ValueError("Maxisation of the likelihood failed") logger.info( - f"Finished with {100 * successes / self.n_prior_samples}% success rate " - f"Maximum log-likelihood {max_logL} " - f"Distribution range {np.min(logL_list)}:{np.max(logL_list)}" + f"Finished with {100 * successes / self.n_prior_samples}% success rate| " + f"Maximum log-likelihood {max_logL}| " + f"(max-mu)/sigma= {(max_logL - np.mean(logL_list)) / np.std(logL_list)} " ) self.minimization_metadata = minout From 78a49ba03d3ff2141a8d84d9b29a55e834ca7baf Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 16 Apr 2025 16:27:38 +0100 Subject: [PATCH 11/30] Include full posterior reweighting --- bilby/core/fisher_sampler.py | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 799156451..5305bee45 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -119,14 +119,19 @@ def run_sampler(self): for _, sample in raw_samples.iterrows() ] ) + raw_logpi = self.priors.ln_prob(raw_samples, axis=0) + + if self.use_ratio: + raw_logl -= self.likelihood.noise_log_likelihood() if self.kwargs["rejection_sampling"]: logger.info(f"Rejection sampling the posterior from {self.kwargs['nsamples']} generated samples") + logl_norm = multivariate_normal.logpdf( raw_samples, mean=fisher_mpe.mean, cov=fisher_mpe.iFIM ) - # TODO: Add prior weights - ln_weights = raw_logl - logl_norm + + ln_weights = raw_logl + raw_logpi - logl_norm ln_weights -= np.mean(ln_weights) weights = np.exp(ln_weights) From 49f2b5ae82d3ca6b506b6a315bdb3c1cc8f75cee Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 17 Apr 2025 10:31:23 +0100 Subject: [PATCH 12/30] Apply correction to out-of-prior samples --- bilby/core/fisher_matrix.py | 15 ++++++++++++--- bilby/core/fisher_sampler.py | 29 ++++++++++++++++++++++------- 2 files changed, 34 insertions(+), 10 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index 8a1279fce..ee08b540f 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -54,7 +54,9 @@ def __init__(self, likelihood, priors, parameters=None, minimization_method="Nel self.prior_samples = [ priors.sample_subset(self.parameter_names) for _ in range(n_prior_samples) ] - self.prior_bounds = [(priors[key].minimum, priors[key].maximum) for key in self.parameter_names] + self.prior_bounds_min = np.array([priors[key].minimum for key in self.parameter_names]) + self.prior_bounds_max = np.array([priors[key].maximum for key in self.parameter_names]) + self.prior_bounds = list(zip(self.prior_bounds_min, self.prior_bounds_max)) self.prior_width_dict = {} for key in self.parameter_names: @@ -105,15 +107,22 @@ def calculate_FIM(self, sample): return FIM else: import scipy.differentiate as sd - logger.info("Using Scipy differentiate") + logger.info("Using Scipy differentiate to estimate the Fisher information matrix (FIM)") point = np.array([sample[key] for key in self.parameter_names]) - res = sd.hessian(self.log_likelihood_from_array, point) + res = sd.hessian(self.log_likelihood_from_array, point, initial_step=0.5) FIM = - res.ddf + logger.info(f"Estimated FIM:\n{FIM}") return FIM def log_likelihood_from_array(self, x_array): def wrapped_logl(x_array): + # Map points outside the bounds to the bounds + idxs = x_array < self.prior_bounds_min + x_array[idxs] = self.prior_bounds_min[idxs] + idxs = x_array > self.prior_bounds_max + x_array[idxs] = self.prior_bounds_max[idxs] + return self.log_likelihood(array_to_dict(self.parameter_names, x_array)) def wrapped_logl_arb(x_array): diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 5305bee45..7f1c04559 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -2,6 +2,7 @@ import numpy as np from scipy.stats import multivariate_normal +import tqdm from .fisher_matrix import FisherMatrixPosteriorEstimator from .sampler.base_sampler import Sampler, signal_wrapper @@ -113,13 +114,19 @@ def run_sampler(self): ) raw_samples = fisher_mpe.sample_dataframe("maxL", self.kwargs["nsamples"]) - raw_logl = np.array( - [ - fisher_mpe.log_likelihood(sample.to_dict()) - for _, sample in raw_samples.iterrows() - ] - ) - raw_logpi = self.priors.ln_prob(raw_samples, axis=0) + + raw_logpi = [] + raw_logl = [] + logger.info("Calculating the likelihood and priors") + for _, rs in tqdm.tqdm(raw_samples.iterrows(), total=len(raw_samples)): + raw_logpi.append(self.priors.ln_prob(rs.to_dict(), axis=0)) + if np.isinf(raw_logpi[-1]): + raw_logl.append(-np.inf) + else: + raw_logl.append(fisher_mpe.log_likelihood(rs.to_dict())) + + raw_logpi = np.real(np.array(raw_logpi)) + raw_logl = np.array(raw_logl) if self.use_ratio: raw_logl -= self.likelihood.noise_log_likelihood() @@ -132,6 +139,14 @@ def run_sampler(self): ) ln_weights = raw_logl + raw_logpi - logl_norm + + # Remove impossible samples + idxs = ~np.isinf(ln_weights) + ln_weights = ln_weights[idxs] + raw_samples = raw_samples[idxs] + raw_logl = raw_logl[idxs] + raw_logpi = raw_logpi[idxs] + ln_weights -= np.mean(ln_weights) weights = np.exp(ln_weights) From 87a337da356e86d78d8bfcd30576cd863dcc6599 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Tue, 22 Apr 2025 17:01:01 +0100 Subject: [PATCH 13/30] Refactor and clean up --- bilby/core/fisher_matrix.py | 9 +- bilby/core/fisher_sampler.py | 248 +++++++++++++++++++++++------------ 2 files changed, 174 insertions(+), 83 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index ee08b540f..cb46a5531 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -76,11 +76,17 @@ def log_likelihood(self, sample): def calculate_iFIM(self, sample): FIM = self.calculate_FIM(sample) + + # Force the FIM to be symmetric by averaging off-diagonal estimates + upper_off_diagonal_average = .5 * (np.triu(FIM, 1) + np.triu(FIM.T, 1)) + FIM = np.diag(np.diag(FIM)) + upper_off_diagonal_average + upper_off_diagonal_average.T + iFIM = scipy.linalg.inv(FIM) # Ensure iFIM is positive definite min_eig = np.min(np.real(np.linalg.eigvals(iFIM))) if min_eig < 0: + logger.warning("Scaling the iFIM to ensure it is positive definite") iFIM -= 10 * min_eig * np.eye(*iFIM.shape) return iFIM @@ -111,7 +117,7 @@ def calculate_FIM(self, sample): point = np.array([sample[key] for key in self.parameter_names]) res = sd.hessian(self.log_likelihood_from_array, point, initial_step=0.5) FIM = - res.ddf - logger.info(f"Estimated FIM:\n{FIM}") + logger.debug(f"Estimated FIM:\n{FIM}") return FIM @@ -120,6 +126,7 @@ def wrapped_logl(x_array): # Map points outside the bounds to the bounds idxs = x_array < self.prior_bounds_min x_array[idxs] = self.prior_bounds_min[idxs] + idxs = x_array > self.prior_bounds_max x_array[idxs] = self.prior_bounds_max[idxs] diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 7f1c04559..de010d7b1 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -1,13 +1,14 @@ import datetime import numpy as np +import pandas as pd from scipy.stats import multivariate_normal import tqdm from .fisher_matrix import FisherMatrixPosteriorEstimator from .sampler.base_sampler import Sampler, signal_wrapper from .result import rejection_sample -from .utils import logger, kish_log_effective_sample_size, safe_save_figure +from .utils import logger, kish_log_effective_sample_size, safe_save_figure, random class Fisher(Sampler): @@ -52,6 +53,8 @@ class Fisher(Sampler): minimization_method="Nelder-Mead", n_prior_samples=100, fd_eps=1e-6, + mirror_diagnostic_plot=False, + cov_scaling=1, ) def __init__( @@ -113,59 +116,28 @@ def run_sampler(self): fd_eps=self.kwargs["fd_eps"], ) - raw_samples = fisher_mpe.sample_dataframe("maxL", self.kwargs["nsamples"]) + maxL_sample_dict = fisher_mpe.get_maximum_likelihood_sample() + maxL_sample_array = np.array(list(maxL_sample_dict.values())) + iFIM = fisher_mpe.calculate_iFIM(maxL_sample_dict) + cov = self.kwargs["cov_scaling"] * iFIM - raw_logpi = [] - raw_logl = [] - logger.info("Calculating the likelihood and priors") - for _, rs in tqdm.tqdm(raw_samples.iterrows(), total=len(raw_samples)): - raw_logpi.append(self.priors.ln_prob(rs.to_dict(), axis=0)) - if np.isinf(raw_logpi[-1]): - raw_logl.append(-np.inf) - else: - raw_logl.append(fisher_mpe.log_likelihood(rs.to_dict())) + msg = "Generation-distribution: " + "| ".join( + [f"{key}: {val:0.5f} +/- {np.sqrt(var):.5f}" for ((key, val), var) + in zip(maxL_sample_dict.items(), np.diag(cov))] + ) + logger.info(msg) - raw_logpi = np.real(np.array(raw_logpi)) - raw_logl = np.array(raw_logl) + g_samples, g_logl, g_logpi = self._draw_samples_from_generating_distribution( + maxL_sample_array, cov, fisher_mpe) if self.use_ratio: - raw_logl -= self.likelihood.noise_log_likelihood() + g_logl -= self.likelihood.noise_log_likelihood() if self.kwargs["rejection_sampling"]: - logger.info(f"Rejection sampling the posterior from {self.kwargs['nsamples']} generated samples") - - logl_norm = multivariate_normal.logpdf( - raw_samples, mean=fisher_mpe.mean, cov=fisher_mpe.iFIM - ) - - ln_weights = raw_logl + raw_logpi - logl_norm - - # Remove impossible samples - idxs = ~np.isinf(ln_weights) - ln_weights = ln_weights[idxs] - raw_samples = raw_samples[idxs] - raw_logl = raw_logl[idxs] - raw_logpi = raw_logpi[idxs] - - ln_weights -= np.mean(ln_weights) - weights = np.exp(ln_weights) - - samples, idxs = rejection_sample(raw_samples, weights, return_idxs=True) - logl = raw_logl[idxs] - - nsamples = len(samples) - efficiency = 100 * nsamples / len(raw_samples) - ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) - logger.info( - f"Rejection sampling Fisher posterior produced {nsamples} samples" - f" with an efficiency of {efficiency}% and effective sample" - f" size {ess}" - ) - if self.plot: - self.create_rejection_sample_diagnostic(samples, raw_samples, fisher_mpe.mean) + samples, logl = self._rejection_sample(g_samples, g_logl, g_logpi, maxL_sample_array, cov) else: - samples = raw_samples - logl = raw_logl + samples = g_samples + logl = g_logl end_time = datetime.datetime.now() self.sampling_time = end_time - self.start_time @@ -174,15 +146,15 @@ def run_sampler(self): return self.result - def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL): + def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL, weights): import corner import matplotlib.pyplot as plt import matplotlib.lines as mpllines kwargs = dict( bins=50, - smooth=0.9, - max_n_ticks=3, + smooth=0.7, + max_n_ticks=5, truths=maxL, truth_color="C3", ) @@ -194,67 +166,100 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL): # Create the data array to plot and pass everything to corner xs = samples[self.search_parameter_keys].values rxs = raw_samples[self.search_parameter_keys].values + + g_color = "k" + g_ls = "--" + f_color = "C0" + f_ls = "-" + if len(self.search_parameter_keys) > 1: lines = [] - ls = "-" - c = "C0" fig = corner.corner( rxs, - color=c, - contour_kwargs={"linestyles": ls, "alpha": 0.8}, - hist_kwargs={"density": True, "ls": ls, "alpha": 0.8}, + color=g_color, + contour_kwargs={"linestyles": g_ls, "alpha": 0.8}, + hist_kwargs={"density": True, "ls": g_ls, "alpha": 0.8}, data_kwargs={"alpha": 1}, no_fill_contours=True, alpha=0.8, plot_density=False, - plot_datapoints=True, - fill_contours=False, - quantiles=[0.16, 0.84], - levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), - **kwargs - ) - lines.append(mpllines.Line2D([0], [0], color=c, linestyle=ls)) - - ls = "--" - c = "C1" - fig = corner.corner( - xs, - color=c, - contour_kwargs={"linestyles": ls, "alpha": 0.8}, - contourf_kwargs={"alpha": 0.8}, - hist_kwargs={"density": True, "ls": ls, "alpha": 0.8}, - no_fill_contours=True, - fig=fig, - alpha=0.1, - plot_density=True, plot_datapoints=False, fill_contours=False, quantiles=[0.16, 0.84], levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), - **kwargs, + **kwargs ) - lines.append(mpllines.Line2D([0], [0], color=c, linestyle=ls)) + lines.append(mpllines.Line2D([0], [0], color=g_color, linestyle=g_ls)) + + if len(xs) > len(samples.keys()): + fig = corner.corner( + xs, + color=f_color, + contour_kwargs={"linestyles": f_ls, "alpha": 0.8}, + contourf_kwargs={"alpha": 0.8}, + hist_kwargs={"density": True, "ls": f_ls, "alpha": 0.8}, + no_fill_contours=True, + fig=fig, + alpha=0.1, + plot_density=True, + plot_datapoints=False, + fill_contours=False, + quantiles=[0.16, 0.84], + levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), + **kwargs, + ) + lines.append(mpllines.Line2D([0], [0], color=f_color, linestyle=f_ls)) + else: + logger.info("Too few samples to add to the diagnostic plot") axes = fig.get_axes() ndim = int(np.sqrt(len(axes))) - axes[ndim - 1].legend(lines, ["$g(x)$", "$f(x)$"]) + axes[0].legend(lines, ["$g(x)$", "$f(x)$"]) + + axes = np.array(axes).reshape((ndim, ndim)) + + base_alpha = 0.1 + alphas = base_alpha + (1 - base_alpha) * weights + for ii in range(1, ndim): + for jj in range(0, ndim - 1): + if ii <= jj: + continue + if self.kwargs["mirror_diagnostic_plot"]: + ax = axes[jj, ii] + xsamples = rxs[:, ii] + ysamples = rxs[:, jj] + else: + ax = axes[ii, jj] + xsamples = rxs[:, jj] + ysamples = rxs[:, ii] + sc = ax.scatter( + xsamples, + ysamples, + c=np.log(weights), + alpha=alphas, + edgecolor='none', + vmin=-5, + vmax=0 + ) + cbar = fig.colorbar(sc, ax=axes[0, -1]) + cbar.set_label("Log-weight") else: fig, ax = plt.subplots() ax.hist( xs, bins=kwargs["bins"], - color="C0", + color=f_color, histtype="step", - ls="-", + ls=f_ls, density=True, ) ax.hist( rxs, bins=kwargs["bins"], - color="C1", + color=g_color, histtype="step", - ls="--", + ls=g_ls, density=True, ) ax.set_xlabel(kwargs["labels"][0]) @@ -284,3 +289,82 @@ def _generate_result(self, samples, log_likelihood_evaluations): neffsamples=None, sampling_time_s=self.sampling_time.total_seconds(), ) + + def _draw_samples_from_generating_distribution(self, sample_array, cov, fisher_mpe): + samples_array = random.rng.multivariate_normal(sample_array, cov, self.kwargs["nsamples"]) + samples = pd.DataFrame(samples_array, columns=fisher_mpe.parameter_names) + + logpi = [] + logl = [] + logger.info("Calculating the likelihood and priors") + outside_prior_count = 0 + for _, rs in tqdm.tqdm(samples.iterrows(), total=len(samples)): + logpi.append(self.priors.ln_prob(rs.to_dict(), axis=0)) + if np.isinf(logpi[-1]): + outside_prior_count += 1 + logl.append(-np.inf) + else: + logl.append(fisher_mpe.log_likelihood(rs.to_dict())) + + if outside_prior_count < len(samples): + logger.info( + f"Discarding {100 * outside_prior_count / len(samples):0.3f}% of samples that" + " fall outside the prior" + ) + else: + raise ValueError("Sampling has failed: no viable samples left") + + logpi = np.real(np.array(logpi)) + logl = np.array(logl) + return samples, logl, logpi + + def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): + logger.info(f"Rejection sampling the posterior from {len(g_samples)} samples") + + g_logl_norm = multivariate_normal.logpdf( + g_samples, mean=mean, cov=cov + ) + + ln_weights = g_logl + g_logpi - g_logl_norm + + # Remove impossible samples + idxs = ~np.isinf(ln_weights) + ln_weights = ln_weights[idxs] + g_samples = g_samples[idxs] + g_logl = g_logl[idxs] + g_logpi = g_logpi[idxs] + + # Scale + ln_weights -= np.max(ln_weights) + + # Sort by weight + idxs = np.argsort(ln_weights) + g_samples = g_samples.iloc[idxs] + g_logl = g_logl[idxs] + g_logpi = g_logpi[idxs] + ln_weights = ln_weights[idxs] + + weights = np.exp(ln_weights) + + samples, idxs = rejection_sample(g_samples, weights, return_idxs=True) + logl = g_logl[idxs] + + nsamples = len(samples) + + if self.plot: + self.create_rejection_sample_diagnostic(samples, g_samples, mean, weights) + + if nsamples == 1: + raise ValueError( + "Rejection sampling has produced a single sample and therefore failed." + ) + + efficiency = 100 * nsamples / len(g_samples) + ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) + logger.info( + f"Rejection sampling Fisher posterior produced {nsamples} samples" + f" with an efficiency of {efficiency:0.3f}% and effective sample" + f" size {ess}" + ) + + return samples, logl From a0d0cecb63740aa67e4591f541c720c1e4dbf69e Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 08:59:21 +0100 Subject: [PATCH 14/30] Add option to use the injection for the maximum likelihood initialisation --- bilby/core/fisher_matrix.py | 3 ++- bilby/core/fisher_sampler.py | 7 ++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index cb46a5531..dc331a7cb 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -223,7 +223,8 @@ def get_maximum_likelihood_sample(self, initial_sample=None): """ if initial_sample: - out = self._maximize_likelihood_from_initial_sample(initial_sample) + logger.info(f"Maximising the likelihood from initial sample {initial_sample}") + minout = self._maximize_likelihood_from_initial_sample(initial_sample) else: logger.info(f"Maximising the likelihood using {self.n_prior_samples} prior samples") max_logL = -np.inf diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index de010d7b1..1d9388785 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -55,6 +55,7 @@ class Fisher(Sampler): fd_eps=1e-6, mirror_diagnostic_plot=False, cov_scaling=1, + use_injection_for_maxL=True, ) def __init__( @@ -116,7 +117,11 @@ def run_sampler(self): fd_eps=self.kwargs["fd_eps"], ) - maxL_sample_dict = fisher_mpe.get_maximum_likelihood_sample() + if self.injection_parameters and "use_injection_for_maxL" in self.kwargs: + sample = {key: self.injection_parameters[key] for key in fisher_mpe.parameter_names} + else: + sample = None + maxL_sample_dict = fisher_mpe.get_maximum_likelihood_sample(sample) maxL_sample_array = np.array(list(maxL_sample_dict.values())) iFIM = fisher_mpe.calculate_iFIM(maxL_sample_dict) cov = self.kwargs["cov_scaling"] * iFIM From f5c4d19b3ed91bb75aca520b74c3751ebe5c150e Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 09:13:50 +0100 Subject: [PATCH 15/30] Reformat the code to enable alternative resampling approaches --- bilby/core/fisher_sampler.py | 26 ++++++++++++++++++-------- 1 file changed, 18 insertions(+), 8 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 1d9388785..40b6a677e 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -48,7 +48,7 @@ class Fisher(Sampler): sampler_name = "fisher" sampling_seed_key = "seed" default_kwargs = dict( - rejection_sampling=True, + resample="importance", nsamples=1000, minimization_method="Nelder-Mead", n_prior_samples=100, @@ -138,8 +138,10 @@ def run_sampler(self): if self.use_ratio: g_logl -= self.likelihood.noise_log_likelihood() - if self.kwargs["rejection_sampling"]: + if self.kwargs["resample"] == "rejection": samples, logl = self._rejection_sample(g_samples, g_logl, g_logpi, maxL_sample_array, cov) + elif self.kwargs["resample"] == "importance": + samples, logl = self._importance_sample(g_samples, g_logl, g_logpi, maxL_sample_array, cov) else: samples = g_samples logl = g_logl @@ -323,9 +325,7 @@ def _draw_samples_from_generating_distribution(self, sample_array, cov, fisher_m logl = np.array(logl) return samples, logl, logpi - def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): - logger.info(f"Rejection sampling the posterior from {len(g_samples)} samples") - + def _calculate_weights(self, g_samples, g_logl, g_logpi, mean, cov): g_logl_norm = multivariate_normal.logpdf( g_samples, mean=mean, cov=cov ) @@ -349,8 +349,20 @@ def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): g_logpi = g_logpi[idxs] ln_weights = ln_weights[idxs] + ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) + logger.info(f"Calculated weights have an effective sample size {ess}") + weights = np.exp(ln_weights) + return weights + + def _importance_sample(self, g_samples, g_logl, g_logpi, mean, cov): + raise NotImplementedError() + + def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): + logger.info(f"Rejection sampling the posterior from {len(g_samples)} samples") + weights = self._calculate_weights(g_samples, g_logl, g_logpi, mean, cov) + samples, idxs = rejection_sample(g_samples, weights, return_idxs=True) logl = g_logl[idxs] @@ -365,11 +377,9 @@ def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): ) efficiency = 100 * nsamples / len(g_samples) - ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) logger.info( f"Rejection sampling Fisher posterior produced {nsamples} samples" - f" with an efficiency of {efficiency:0.3f}% and effective sample" - f" size {ess}" + f" with an efficiency of {efficiency:0.3f}%" ) return samples, logl From 680380e1bde1006f5010fead6add7ee7b4646a47 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 09:35:15 +0100 Subject: [PATCH 16/30] Add importance sampling methods --- bilby/core/fisher_sampler.py | 48 +++++++++++++++++++++++++----------- 1 file changed, 33 insertions(+), 15 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 40b6a677e..33c683b7e 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -53,6 +53,7 @@ class Fisher(Sampler): minimization_method="Nelder-Mead", n_prior_samples=100, fd_eps=1e-6, + plot_diagnostic=False, mirror_diagnostic_plot=False, cov_scaling=1, use_injection_for_maxL=True, @@ -153,7 +154,7 @@ def run_sampler(self): return self.result - def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL, weights): + def create_resample_diagnostic(self, samples, raw_samples, maxL, weights, method): import corner import matplotlib.pyplot as plt import matplotlib.lines as mpllines @@ -174,6 +175,11 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL, weights xs = samples[self.search_parameter_keys].values rxs = raw_samples[self.search_parameter_keys].values + # Sort by weight (only for plotting) + idxs = np.argsort(weights) + rxs = rxs[idxs] + weights = weights[idxs] + g_color = "k" g_ls = "--" f_color = "C0" @@ -271,8 +277,9 @@ def create_rejection_sample_diagnostic(self, samples, raw_samples, maxL, weights ) ax.set_xlabel(kwargs["labels"][0]) - filename = "{}/{}_rejection_sample.png".format(self.outdir, self.label) - logger.debug("Saving rejection-sample diagnopstic plot to {}".format(filename)) + fig.suptitle(f"Resampling method: {method}") + + filename = f"{self.outdir}/{self.label}_resample_{method}.png" safe_save_figure(fig=fig, filename=filename, dpi=400) plt.close(fig) @@ -342,22 +349,33 @@ def _calculate_weights(self, g_samples, g_logl, g_logpi, mean, cov): # Scale ln_weights -= np.max(ln_weights) - # Sort by weight - idxs = np.argsort(ln_weights) - g_samples = g_samples.iloc[idxs] - g_logl = g_logl[idxs] - g_logpi = g_logpi[idxs] - ln_weights = ln_weights[idxs] - - ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) - logger.info(f"Calculated weights have an effective sample size {ess}") + self.ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) + logger.info(f"Calculated weights have an effective sample size {self.ess}") weights = np.exp(ln_weights) return weights def _importance_sample(self, g_samples, g_logl, g_logpi, mean, cov): - raise NotImplementedError() + logger.info(f"Importance sampling the posterior from {len(g_samples)} samples") + weights = self._calculate_weights(g_samples, g_logl, g_logpi, mean, cov) + + normalized_weights = weights / np.sum(weights) + idxs = np.random.choice(len(g_samples), size=self.ess, p=normalized_weights) + samples = g_samples.iloc[idxs] + logl = g_logl[idxs] + + if self.kwargs["plot_diagnostic"]: + self.create_resample_diagnostic(samples, g_samples, mean, weights, method="importance") + + nsamples = len(samples) + efficiency = 100 * nsamples / len(g_samples) + logger.info( + f"Importance sampling Fisher posterior produced {nsamples} samples" + f" with an efficiency of {efficiency:0.3f}%" + ) + + return samples, logl def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): logger.info(f"Rejection sampling the posterior from {len(g_samples)} samples") @@ -368,8 +386,8 @@ def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): nsamples = len(samples) - if self.plot: - self.create_rejection_sample_diagnostic(samples, g_samples, mean, weights) + if self.kwargs["plot_diagnostic"]: + self.create_resample_diagnostic(samples, g_samples, mean, weights, method="importance") if nsamples == 1: raise ValueError( From 1088a2067563f066549dd4da8dbf2ebfc9552331 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 12:30:59 +0100 Subject: [PATCH 17/30] Add method to utilise array method if possible --- bilby/core/fisher_sampler.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 33c683b7e..7eb4388ef 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -311,14 +311,17 @@ def _draw_samples_from_generating_distribution(self, sample_array, cov, fisher_m logpi = [] logl = [] logger.info("Calculating the likelihood and priors") - outside_prior_count = 0 - for _, rs in tqdm.tqdm(samples.iterrows(), total=len(samples)): - logpi.append(self.priors.ln_prob(rs.to_dict(), axis=0)) - if np.isinf(logpi[-1]): - outside_prior_count += 1 - logl.append(-np.inf) - else: - logl.append(fisher_mpe.log_likelihood(rs.to_dict())) + + logpi = self.priors.ln_prob(samples, axis=0) + outside_prior_count = np.sum(np.isinf(logpi)) + if outside_prior_count == 0: + logl = fisher_mpe.log_likelihood_from_array(samples.values.T) + else: + for ii, rs in tqdm.tqdm(samples.iterrows(), total=len(samples)): + if np.isinf(logpi[ii]): + logl.append(-np.inf) + else: + logl.append(fisher_mpe.log_likelihood(rs.to_dict())) if outside_prior_count < len(samples): logger.info( From cb099ca6e2235c72492048022fc6c7f27c5c23ff Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 13:16:42 +0100 Subject: [PATCH 18/30] Improve error messages --- bilby/core/fisher_sampler.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 7eb4388ef..f43e024e9 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -6,7 +6,7 @@ import tqdm from .fisher_matrix import FisherMatrixPosteriorEstimator -from .sampler.base_sampler import Sampler, signal_wrapper +from .sampler.base_sampler import Sampler, signal_wrapper, SamplerError from .result import rejection_sample from .utils import logger, kish_log_effective_sample_size, safe_save_figure, random @@ -127,7 +127,7 @@ def run_sampler(self): iFIM = fisher_mpe.calculate_iFIM(maxL_sample_dict) cov = self.kwargs["cov_scaling"] * iFIM - msg = "Generation-distribution: " + "| ".join( + msg = "Generation-distribution:\n " + "\n ".join( [f"{key}: {val:0.5f} +/- {np.sqrt(var):.5f}" for ((key, val), var) in zip(maxL_sample_dict.items(), np.diag(cov))] ) @@ -226,15 +226,14 @@ def create_resample_diagnostic(self, samples, raw_samples, maxL, weights, method logger.info("Too few samples to add to the diagnostic plot") axes = fig.get_axes() - ndim = int(np.sqrt(len(axes))) axes[0].legend(lines, ["$g(x)$", "$f(x)$"]) - axes = np.array(axes).reshape((ndim, ndim)) + axes = np.array(axes).reshape((self.ndim, self.ndim)) base_alpha = 0.1 alphas = base_alpha + (1 - base_alpha) * weights - for ii in range(1, ndim): - for jj in range(0, ndim - 1): + for ii in range(1, self.ndim): + for jj in range(0, self.ndim - 1): if ii <= jj: continue if self.kwargs["mirror_diagnostic_plot"]: @@ -314,7 +313,7 @@ def _draw_samples_from_generating_distribution(self, sample_array, cov, fisher_m logpi = self.priors.ln_prob(samples, axis=0) outside_prior_count = np.sum(np.isinf(logpi)) - if outside_prior_count == 0: + if outside_prior_count == 0 and False: logl = fisher_mpe.log_likelihood_from_array(samples.values.T) else: for ii, rs in tqdm.tqdm(samples.iterrows(), total=len(samples)): @@ -329,7 +328,7 @@ def _draw_samples_from_generating_distribution(self, sample_array, cov, fisher_m " fall outside the prior" ) else: - raise ValueError("Sampling has failed: no viable samples left") + raise SamplerError("Sampling has failed: no viable samples left") logpi = np.real(np.array(logpi)) logl = np.array(logl) @@ -371,6 +370,9 @@ def _importance_sample(self, g_samples, g_logl, g_logpi, mean, cov): if self.kwargs["plot_diagnostic"]: self.create_resample_diagnostic(samples, g_samples, mean, weights, method="importance") + if self.ess < self.ndim: + raise SamplerError("Effective sample size less than ndim: sampling has failed") + nsamples = len(samples) efficiency = 100 * nsamples / len(g_samples) logger.info( @@ -392,10 +394,8 @@ def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): if self.kwargs["plot_diagnostic"]: self.create_resample_diagnostic(samples, g_samples, mean, weights, method="importance") - if nsamples == 1: - raise ValueError( - "Rejection sampling has produced a single sample and therefore failed." - ) + if nsamples < self.ndim: + raise SamplerError("Number of samples less than ndim: sampling has failed") efficiency = 100 * nsamples / len(g_samples) logger.info( From 6fa36460c55b64deeed0619584dd5e1d4bbde092 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 13:22:38 +0100 Subject: [PATCH 19/30] Refactor and rename --- bilby/core/fisher_sampler.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index f43e024e9..792ba6674 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -49,9 +49,10 @@ class Fisher(Sampler): sampling_seed_key = "seed" default_kwargs = dict( resample="importance", - nsamples=1000, + target_nsamples=10000, + batch_nsamples=10000, + prior_nsamples=100, minimization_method="Nelder-Mead", - n_prior_samples=100, fd_eps=1e-6, plot_diagnostic=False, mirror_diagnostic_plot=False, @@ -114,7 +115,7 @@ def run_sampler(self): likelihood=self.likelihood, priors=self.priors, minimization_method=self.kwargs["minimization_method"], - n_prior_samples=self.kwargs["n_prior_samples"], + n_prior_samples=self.kwargs["prior_nsamples"], fd_eps=self.kwargs["fd_eps"], ) @@ -134,7 +135,7 @@ def run_sampler(self): logger.info(msg) g_samples, g_logl, g_logpi = self._draw_samples_from_generating_distribution( - maxL_sample_array, cov, fisher_mpe) + maxL_sample_array, cov, fisher_mpe, self.kwargs["batch_nsamples"]) if self.use_ratio: g_logl -= self.likelihood.noise_log_likelihood() @@ -303,8 +304,8 @@ def _generate_result(self, samples, log_likelihood_evaluations): sampling_time_s=self.sampling_time.total_seconds(), ) - def _draw_samples_from_generating_distribution(self, sample_array, cov, fisher_mpe): - samples_array = random.rng.multivariate_normal(sample_array, cov, self.kwargs["nsamples"]) + def _draw_samples_from_generating_distribution(self, mean, cov, fisher_mpe, nsamples): + samples_array = random.rng.multivariate_normal(mean, cov, nsamples) samples = pd.DataFrame(samples_array, columns=fisher_mpe.parameter_names) logpi = [] From bf4cf8213a76404918e235d95d00a876c9c3fb7c Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 15:08:12 +0100 Subject: [PATCH 20/30] Add batch sampling --- bilby/core/fisher_sampler.py | 120 ++++++++++++++++++++--------------- 1 file changed, 69 insertions(+), 51 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 792ba6674..5fddc6534 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -124,7 +124,7 @@ def run_sampler(self): else: sample = None maxL_sample_dict = fisher_mpe.get_maximum_likelihood_sample(sample) - maxL_sample_array = np.array(list(maxL_sample_dict.values())) + mean = np.array(list(maxL_sample_dict.values())) iFIM = fisher_mpe.calculate_iFIM(maxL_sample_dict) cov = self.kwargs["cov_scaling"] * iFIM @@ -134,28 +134,67 @@ def run_sampler(self): ) logger.info(msg) - g_samples, g_logl, g_logpi = self._draw_samples_from_generating_distribution( - maxL_sample_array, cov, fisher_mpe, self.kwargs["batch_nsamples"]) + nsamples = 0 + target_nsamples = self.kwargs["target_nsamples"] + batch_nsamples = self.kwargs["batch_nsamples"] + all_g_samples = [] + all_samples = [] + all_logl = [] + all_weights = [] + resample = self.kwargs["resample"] + logger.info(f"Starting sampling in batches of {batch_nsamples} to produce {target_nsamples} samples") + pbar = tqdm.tqdm(total=target_nsamples, desc=f"{resample.capitalize()} sampling") + while nsamples < target_nsamples: + g_samples, g_logl, g_logpi = self._draw_samples_from_generating_distribution( + mean, cov, fisher_mpe, batch_nsamples + ) - if self.use_ratio: - g_logl -= self.likelihood.noise_log_likelihood() + _methods = dict(rejection=self._rejection_sample, importance=self._importance_sample) - if self.kwargs["resample"] == "rejection": - samples, logl = self._rejection_sample(g_samples, g_logl, g_logpi, maxL_sample_array, cov) - elif self.kwargs["resample"] == "importance": - samples, logl = self._importance_sample(g_samples, g_logl, g_logpi, maxL_sample_array, cov) - else: - samples = g_samples - logl = g_logl + if resample in _methods: + weights = self._calculate_weights(g_samples, g_logl, g_logpi, mean, cov) + samples, logl = _methods[resample](g_samples, g_logl, weights) + efficiency = 100 * len(samples) / len(g_samples) + else: + logger.info("No resampling applied") + samples = g_samples + logl = g_logl + weights = np.ones_like(logl) + + nsamples += len(samples) + pbar.set_postfix({ + 'eff': f"{efficiency:.3f}%", + }) + pbar.update(len(samples)) + all_g_samples.append(g_samples) + all_samples.append(samples) + all_logl.append(logl) + all_weights.append(weights) + + pbar.close() + + g_samples = pd.concat(all_g_samples) + samples = pd.concat(all_samples) + logl = np.concat(all_logl) + weights = np.concat(all_weights) + efficiency = 100 * len(samples) / len(g_samples) + + logger.info(f"Finished sampling: total efficiency is {efficiency:0.3f}%") + + if self.kwargs["plot_diagnostic"]: + self.create_resample_diagnostic(samples, g_samples, mean, weights, method=self.kwargs["resample"]) end_time = datetime.datetime.now() self.sampling_time = end_time - self.start_time - self._generate_result(samples, logl) + if self.use_ratio: + logl -= self.likelihood.noise_log_likelihood() + + self._generate_result(samples, logl, efficiency=efficiency, nlikelihood=len(g_samples)) return self.result - def create_resample_diagnostic(self, samples, raw_samples, maxL, weights, method): + def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method): import corner import matplotlib.pyplot as plt import matplotlib.lines as mpllines @@ -164,8 +203,9 @@ def create_resample_diagnostic(self, samples, raw_samples, maxL, weights, method bins=50, smooth=0.7, max_n_ticks=5, - truths=maxL, + truths=mean, truth_color="C3", + truth_width=0.5, ) kwargs["labels"] = [ @@ -285,7 +325,7 @@ def create_resample_diagnostic(self, samples, raw_samples, maxL, weights, method return fig - def _generate_result(self, samples, log_likelihood_evaluations): + def _generate_result(self, samples, log_likelihood_evaluations, **run_stats): """ Extract the information we need from the output. @@ -298,11 +338,8 @@ def _generate_result(self, samples, log_likelihood_evaluations): self.result.samples = samples self.result.log_likelihood_evaluations = log_likelihood_evaluations - self.result.meta_data["run_statistics"] = dict( - nlikelihood=None, - neffsamples=None, - sampling_time_s=self.sampling_time.total_seconds(), - ) + run_stats["sampling_time_s"] = self.sampling_time.total_seconds() + self.result.meta_data["run_statistics"] = run_stats def _draw_samples_from_generating_distribution(self, mean, cov, fisher_mpe, nsamples): samples_array = random.rng.multivariate_normal(mean, cov, nsamples) @@ -310,11 +347,11 @@ def _draw_samples_from_generating_distribution(self, mean, cov, fisher_mpe, nsam logpi = [] logl = [] - logger.info("Calculating the likelihood and priors") + logger.debug("Calculating the likelihood and priors") logpi = self.priors.ln_prob(samples, axis=0) outside_prior_count = np.sum(np.isinf(logpi)) - if outside_prior_count == 0 and False: + if outside_prior_count == 0: logl = fisher_mpe.log_likelihood_from_array(samples.values.T) else: for ii, rs in tqdm.tqdm(samples.iterrows(), total=len(samples)): @@ -323,13 +360,15 @@ def _draw_samples_from_generating_distribution(self, mean, cov, fisher_mpe, nsam else: logl.append(fisher_mpe.log_likelihood(rs.to_dict())) - if outside_prior_count < len(samples): + if outside_prior_count == len(samples): + raise SamplerError("Sampling has failed: no viable samples left") + elif outside_prior_count > 0: logger.info( f"Discarding {100 * outside_prior_count / len(samples):0.3f}% of samples that" " fall outside the prior" ) else: - raise SamplerError("Sampling has failed: no viable samples left") + logger.debug("No samples outside the prior bounds") logpi = np.real(np.array(logpi)) logl = np.array(logl) @@ -353,55 +392,34 @@ def _calculate_weights(self, g_samples, g_logl, g_logpi, mean, cov): ln_weights -= np.max(ln_weights) self.ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) - logger.info(f"Calculated weights have an effective sample size {self.ess}") + logger.debug(f"Calculated weights have an effective sample size {self.ess}") weights = np.exp(ln_weights) return weights - def _importance_sample(self, g_samples, g_logl, g_logpi, mean, cov): - logger.info(f"Importance sampling the posterior from {len(g_samples)} samples") - weights = self._calculate_weights(g_samples, g_logl, g_logpi, mean, cov) + def _importance_sample(self, g_samples, g_logl, weights): + logger.debug(f"Importance sampling the posterior from {len(g_samples)} samples") normalized_weights = weights / np.sum(weights) idxs = np.random.choice(len(g_samples), size=self.ess, p=normalized_weights) samples = g_samples.iloc[idxs] logl = g_logl[idxs] - if self.kwargs["plot_diagnostic"]: - self.create_resample_diagnostic(samples, g_samples, mean, weights, method="importance") - if self.ess < self.ndim: raise SamplerError("Effective sample size less than ndim: sampling has failed") - nsamples = len(samples) - efficiency = 100 * nsamples / len(g_samples) - logger.info( - f"Importance sampling Fisher posterior produced {nsamples} samples" - f" with an efficiency of {efficiency:0.3f}%" - ) - return samples, logl - def _rejection_sample(self, g_samples, g_logl, g_logpi, mean, cov): - logger.info(f"Rejection sampling the posterior from {len(g_samples)} samples") - weights = self._calculate_weights(g_samples, g_logl, g_logpi, mean, cov) + def _rejection_sample(self, g_samples, g_logl, weights): + logger.debug(f"Rejection sampling the posterior from {len(g_samples)} samples") samples, idxs = rejection_sample(g_samples, weights, return_idxs=True) logl = g_logl[idxs] nsamples = len(samples) - if self.kwargs["plot_diagnostic"]: - self.create_resample_diagnostic(samples, g_samples, mean, weights, method="importance") - if nsamples < self.ndim: raise SamplerError("Number of samples less than ndim: sampling has failed") - efficiency = 100 * nsamples / len(g_samples) - logger.info( - f"Rejection sampling Fisher posterior produced {nsamples} samples" - f" with an efficiency of {efficiency:0.3f}%" - ) - return samples, logl From 8214c9a33a52efbdf2b52d03fda140ad6489d451 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 15:14:38 +0100 Subject: [PATCH 21/30] Apply black to fisher files --- bilby/core/fisher_matrix.py | 60 +++++++++++++++++++++++----------- bilby/core/fisher_sampler.py | 63 ++++++++++++++++++++++++------------ 2 files changed, 84 insertions(+), 39 deletions(-) diff --git a/bilby/core/fisher_matrix.py b/bilby/core/fisher_matrix.py index dc331a7cb..7011408f2 100644 --- a/bilby/core/fisher_matrix.py +++ b/bilby/core/fisher_matrix.py @@ -15,9 +15,16 @@ def array_to_dict(keys, array): class FisherMatrixPosteriorEstimator(object): - def __init__(self, likelihood, priors, parameters=None, minimization_method="Nelder-Mead", - fd_eps=1e-6, n_prior_samples=100): - """ A class to estimate posteriors using a Fisher Information Matrix approach + def __init__( + self, + likelihood, + priors, + parameters=None, + minimization_method="Nelder-Mead", + fd_eps=1e-6, + n_prior_samples=100, + ): + """A class to estimate posteriors using a Fisher Information Matrix approach Parameters ---------- @@ -54,8 +61,12 @@ def __init__(self, likelihood, priors, parameters=None, minimization_method="Nel self.prior_samples = [ priors.sample_subset(self.parameter_names) for _ in range(n_prior_samples) ] - self.prior_bounds_min = np.array([priors[key].minimum for key in self.parameter_names]) - self.prior_bounds_max = np.array([priors[key].maximum for key in self.parameter_names]) + self.prior_bounds_min = np.array( + [priors[key].minimum for key in self.parameter_names] + ) + self.prior_bounds_max = np.array( + [priors[key].maximum for key in self.parameter_names] + ) self.prior_bounds = list(zip(self.prior_bounds_min, self.prior_bounds_max)) self.prior_width_dict = {} @@ -78,8 +89,12 @@ def calculate_iFIM(self, sample): FIM = self.calculate_FIM(sample) # Force the FIM to be symmetric by averaging off-diagonal estimates - upper_off_diagonal_average = .5 * (np.triu(FIM, 1) + np.triu(FIM.T, 1)) - FIM = np.diag(np.diag(FIM)) + upper_off_diagonal_average + upper_off_diagonal_average.T + upper_off_diagonal_average = 0.5 * (np.triu(FIM, 1) + np.triu(FIM.T, 1)) + FIM = ( + np.diag(np.diag(FIM)) + + upper_off_diagonal_average + + upper_off_diagonal_average.T + ) iFIM = scipy.linalg.inv(FIM) @@ -109,14 +124,19 @@ def calculate_FIM(self, sample): FIM = np.zeros((self.N, self.N)) for ii, ii_key in enumerate(self.parameter_names): for jj, jj_key in enumerate(self.parameter_names): - FIM[ii, jj] = -self.get_second_order_derivative(sample, ii_key, jj_key) + FIM[ii, jj] = -self.get_second_order_derivative( + sample, ii_key, jj_key + ) return FIM else: import scipy.differentiate as sd - logger.info("Using Scipy differentiate to estimate the Fisher information matrix (FIM)") + + logger.info( + "Using Scipy differentiate to estimate the Fisher information matrix (FIM)" + ) point = np.array([sample[key] for key in self.parameter_names]) res = sd.hessian(self.log_likelihood_from_array, point, initial_step=0.5) - FIM = - res.ddf + FIM = -res.ddf logger.debug(f"Estimated FIM:\n{FIM}") return FIM @@ -148,13 +168,13 @@ def get_finite_difference_xx(self, sample, ii): p = self._shift_sample_x(sample, ii, 1) m = self._shift_sample_x(sample, ii, -1) - dx = .5 * (p[ii] - m[ii]) + dx = 0.5 * (p[ii] - m[ii]) loglp = self.log_likelihood(p) logl = self.log_likelihood(sample) loglm = self.log_likelihood(m) - return (loglp - 2 * logl + loglm) / dx ** 2 + return (loglp - 2 * logl + loglm) / dx**2 def get_finite_difference_xy(self, sample, ii, jj): # Sample grid @@ -163,8 +183,8 @@ def get_finite_difference_xy(self, sample, ii, jj): mp = self._shift_sample_xy(sample, ii, -1, jj, 1) mm = self._shift_sample_xy(sample, ii, -1, jj, -1) - dx = .5 * (pp[ii] - mm[ii]) - dy = .5 * (pp[jj] - mm[jj]) + dx = 0.5 * (pp[ii] - mm[ii]) + dy = 0.5 * (pp[jj] - mm[jj]) loglpp = self.log_likelihood(pp) loglpm = self.log_likelihood(pm) @@ -200,7 +220,7 @@ def _maximize_likelihood_from_initial_sample(self, initial_sample): x0 = list(initial_sample.values()) def neg_log_like(x): - return - self.log_likelihood_from_array(x) + return -self.log_likelihood_from_array(x) out = minimize( neg_log_like, @@ -211,7 +231,7 @@ def neg_log_like(x): return out def get_maximum_likelihood_sample(self, initial_sample=None): - """ A method to attempt optimization of the maximum likelihood + """A method to attempt optimization of the maximum likelihood This uses a simple scipy optimization approach, starting from a number of draws from the prior to avoid problems with local optimization. @@ -223,10 +243,14 @@ def get_maximum_likelihood_sample(self, initial_sample=None): """ if initial_sample: - logger.info(f"Maximising the likelihood from initial sample {initial_sample}") + logger.info( + f"Maximising the likelihood from initial sample {initial_sample}" + ) minout = self._maximize_likelihood_from_initial_sample(initial_sample) else: - logger.info(f"Maximising the likelihood using {self.n_prior_samples} prior samples") + logger.info( + f"Maximising the likelihood using {self.n_prior_samples} prior samples" + ) max_logL = -np.inf logL_list = [] successes = 0 diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 5fddc6534..185d01e5d 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -120,7 +120,10 @@ def run_sampler(self): ) if self.injection_parameters and "use_injection_for_maxL" in self.kwargs: - sample = {key: self.injection_parameters[key] for key in fisher_mpe.parameter_names} + sample = { + key: self.injection_parameters[key] + for key in fisher_mpe.parameter_names + } else: sample = None maxL_sample_dict = fisher_mpe.get_maximum_likelihood_sample(sample) @@ -129,8 +132,10 @@ def run_sampler(self): cov = self.kwargs["cov_scaling"] * iFIM msg = "Generation-distribution:\n " + "\n ".join( - [f"{key}: {val:0.5f} +/- {np.sqrt(var):.5f}" for ((key, val), var) - in zip(maxL_sample_dict.items(), np.diag(cov))] + [ + f"{key}: {val:0.5f} +/- {np.sqrt(var):.5f}" + for ((key, val), var) in zip(maxL_sample_dict.items(), np.diag(cov)) + ] ) logger.info(msg) @@ -142,14 +147,22 @@ def run_sampler(self): all_logl = [] all_weights = [] resample = self.kwargs["resample"] - logger.info(f"Starting sampling in batches of {batch_nsamples} to produce {target_nsamples} samples") - pbar = tqdm.tqdm(total=target_nsamples, desc=f"{resample.capitalize()} sampling") + logger.info( + f"Starting sampling in batches of {batch_nsamples} to produce {target_nsamples} samples" + ) + pbar = tqdm.tqdm( + total=target_nsamples, desc=f"{resample.capitalize()} sampling" + ) while nsamples < target_nsamples: - g_samples, g_logl, g_logpi = self._draw_samples_from_generating_distribution( - mean, cov, fisher_mpe, batch_nsamples + g_samples, g_logl, g_logpi = ( + self._draw_samples_from_generating_distribution( + mean, cov, fisher_mpe, batch_nsamples + ) ) - _methods = dict(rejection=self._rejection_sample, importance=self._importance_sample) + _methods = dict( + rejection=self._rejection_sample, importance=self._importance_sample + ) if resample in _methods: weights = self._calculate_weights(g_samples, g_logl, g_logpi, mean, cov) @@ -162,9 +175,11 @@ def run_sampler(self): weights = np.ones_like(logl) nsamples += len(samples) - pbar.set_postfix({ - 'eff': f"{efficiency:.3f}%", - }) + pbar.set_postfix( + { + "eff": f"{efficiency:.3f}%", + } + ) pbar.update(len(samples)) all_g_samples.append(g_samples) all_samples.append(samples) @@ -182,7 +197,9 @@ def run_sampler(self): logger.info(f"Finished sampling: total efficiency is {efficiency:0.3f}%") if self.kwargs["plot_diagnostic"]: - self.create_resample_diagnostic(samples, g_samples, mean, weights, method=self.kwargs["resample"]) + self.create_resample_diagnostic( + samples, g_samples, mean, weights, method=self.kwargs["resample"] + ) end_time = datetime.datetime.now() self.sampling_time = end_time - self.start_time @@ -190,7 +207,9 @@ def run_sampler(self): if self.use_ratio: logl -= self.likelihood.noise_log_likelihood() - self._generate_result(samples, logl, efficiency=efficiency, nlikelihood=len(g_samples)) + self._generate_result( + samples, logl, efficiency=efficiency, nlikelihood=len(g_samples) + ) return self.result @@ -241,7 +260,7 @@ def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method fill_contours=False, quantiles=[0.16, 0.84], levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), - **kwargs + **kwargs, ) lines.append(mpllines.Line2D([0], [0], color=g_color, linestyle=g_ls)) @@ -290,9 +309,9 @@ def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method ysamples, c=np.log(weights), alpha=alphas, - edgecolor='none', + edgecolor="none", vmin=-5, - vmax=0 + vmax=0, ) cbar = fig.colorbar(sc, ax=axes[0, -1]) cbar.set_label("Log-weight") @@ -341,7 +360,9 @@ def _generate_result(self, samples, log_likelihood_evaluations, **run_stats): run_stats["sampling_time_s"] = self.sampling_time.total_seconds() self.result.meta_data["run_statistics"] = run_stats - def _draw_samples_from_generating_distribution(self, mean, cov, fisher_mpe, nsamples): + def _draw_samples_from_generating_distribution( + self, mean, cov, fisher_mpe, nsamples + ): samples_array = random.rng.multivariate_normal(mean, cov, nsamples) samples = pd.DataFrame(samples_array, columns=fisher_mpe.parameter_names) @@ -375,9 +396,7 @@ def _draw_samples_from_generating_distribution(self, mean, cov, fisher_mpe, nsam return samples, logl, logpi def _calculate_weights(self, g_samples, g_logl, g_logpi, mean, cov): - g_logl_norm = multivariate_normal.logpdf( - g_samples, mean=mean, cov=cov - ) + g_logl_norm = multivariate_normal.logpdf(g_samples, mean=mean, cov=cov) ln_weights = g_logl + g_logpi - g_logl_norm @@ -407,7 +426,9 @@ def _importance_sample(self, g_samples, g_logl, weights): logl = g_logl[idxs] if self.ess < self.ndim: - raise SamplerError("Effective sample size less than ndim: sampling has failed") + raise SamplerError( + "Effective sample size less than ndim: sampling has failed" + ) return samples, logl From 6386d463e7f55efecf62460c4bdb72fa0b5493e6 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Wed, 23 Apr 2025 15:58:23 +0100 Subject: [PATCH 22/30] Add option to fail on errors or continue --- bilby/core/fisher_sampler.py | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 185d01e5d..ecd0798d6 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -50,7 +50,7 @@ class Fisher(Sampler): default_kwargs = dict( resample="importance", target_nsamples=10000, - batch_nsamples=10000, + batch_nsamples=1000, prior_nsamples=100, minimization_method="Nelder-Mead", fd_eps=1e-6, @@ -58,6 +58,7 @@ class Fisher(Sampler): mirror_diagnostic_plot=False, cov_scaling=1, use_injection_for_maxL=True, + fail_on_error=False ) def __init__( @@ -382,7 +383,11 @@ def _draw_samples_from_generating_distribution( logl.append(fisher_mpe.log_likelihood(rs.to_dict())) if outside_prior_count == len(samples): - raise SamplerError("Sampling has failed: no viable samples left") + msg = "Sampling has failed: no viable samples left" + if self.kwargs["fail_on_error"]: + raise SamplerError(msg) + else: + logger.info(msg) elif outside_prior_count > 0: logger.info( f"Discarding {100 * outside_prior_count / len(samples):0.3f}% of samples that" @@ -402,15 +407,12 @@ def _calculate_weights(self, g_samples, g_logl, g_logpi, mean, cov): # Remove impossible samples idxs = ~np.isinf(ln_weights) - ln_weights = ln_weights[idxs] - g_samples = g_samples[idxs] - g_logl = g_logl[idxs] - g_logpi = g_logpi[idxs] + ln_weights_viable = ln_weights[idxs] # Scale - ln_weights -= np.max(ln_weights) + ln_weights -= np.max(ln_weights_viable) - self.ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights)))) + self.ess = int(np.floor(np.exp(kish_log_effective_sample_size(ln_weights_viable)))) logger.debug(f"Calculated weights have an effective sample size {self.ess}") weights = np.exp(ln_weights) @@ -426,9 +428,11 @@ def _importance_sample(self, g_samples, g_logl, weights): logl = g_logl[idxs] if self.ess < self.ndim: - raise SamplerError( - "Effective sample size less than ndim: sampling has failed" - ) + msg = "Effective sample size less than ndim: sampling has failed" + if self.kwargs["fail_on_error"]: + raise SamplerError(msg) + else: + logger.info(msg) return samples, logl @@ -441,6 +445,10 @@ def _rejection_sample(self, g_samples, g_logl, weights): nsamples = len(samples) if nsamples < self.ndim: - raise SamplerError("Number of samples less than ndim: sampling has failed") + msg = "Number of samples less than ndim: sampling has failed" + if self.kwargs["fail_on_error"]: + raise SamplerError(msg) + else: + logger.info(msg) return samples, logl From 872d62656c16e0eb6d00fd0f12a46d7d05f3d1a6 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 24 Apr 2025 10:09:21 +0100 Subject: [PATCH 23/30] Add the weights to the diagnostic plot --- bilby/core/fisher_sampler.py | 150 ++++++++++++++++------------------- 1 file changed, 69 insertions(+), 81 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index ecd0798d6..6cd841e00 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -223,18 +223,20 @@ def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method bins=50, smooth=0.7, max_n_ticks=5, - truths=mean, + truths=np.concat((mean, [1])), truth_color="C3", - truth_width=0.5, ) kwargs["labels"] = [ label.replace("_", " ") for label in self.search_parameter_keys ] + kwargs["labels"].append("weights") # Create the data array to plot and pass everything to corner xs = samples[self.search_parameter_keys].values + xs = np.concat((xs, np.random.uniform(0, 1, len(xs)).reshape(-1, 1)), axis=1) rxs = raw_samples[self.search_parameter_keys].values + rxs = np.concat((rxs, weights.reshape((-1, 1))), axis=1) # Sort by weight (only for plotting) idxs = np.argsort(weights) @@ -246,96 +248,82 @@ def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method f_color = "C0" f_ls = "-" - if len(self.search_parameter_keys) > 1: - lines = [] + lines = [] + fig = corner.corner( + rxs, + color=g_color, + contour_kwargs={"linestyles": g_ls, "alpha": 0.8}, + hist_kwargs={"density": True, "ls": g_ls, "alpha": 0.8}, + data_kwargs={"alpha": 1}, + no_fill_contours=True, + alpha=0.8, + plot_density=False, + plot_datapoints=False, + fill_contours=False, + levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), + **kwargs, + ) + lines.append(mpllines.Line2D([0], [0], color=g_color, linestyle=g_ls)) + + if len(xs) > len(samples.keys()): fig = corner.corner( - rxs, - color=g_color, - contour_kwargs={"linestyles": g_ls, "alpha": 0.8}, - hist_kwargs={"density": True, "ls": g_ls, "alpha": 0.8}, - data_kwargs={"alpha": 1}, + xs, + color=f_color, + contour_kwargs={"linestyles": f_ls, "alpha": 0.8}, + contourf_kwargs={"alpha": 0.8}, + hist_kwargs={"density": True, "ls": f_ls, "alpha": 0.8}, no_fill_contours=True, - alpha=0.8, - plot_density=False, + fig=fig, + alpha=0.1, + plot_density=True, plot_datapoints=False, fill_contours=False, - quantiles=[0.16, 0.84], levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), + range=[1] * self.ndim + [(0, 1)], **kwargs, ) - lines.append(mpllines.Line2D([0], [0], color=g_color, linestyle=g_ls)) - - if len(xs) > len(samples.keys()): - fig = corner.corner( - xs, - color=f_color, - contour_kwargs={"linestyles": f_ls, "alpha": 0.8}, - contourf_kwargs={"alpha": 0.8}, - hist_kwargs={"density": True, "ls": f_ls, "alpha": 0.8}, - no_fill_contours=True, - fig=fig, - alpha=0.1, - plot_density=True, - plot_datapoints=False, - fill_contours=False, - quantiles=[0.16, 0.84], - levels=(1 - np.exp(-0.5), 1 - np.exp(-2), 1 - np.exp(-9 / 2.0)), - **kwargs, - ) - lines.append(mpllines.Line2D([0], [0], color=f_color, linestyle=f_ls)) - else: - logger.info("Too few samples to add to the diagnostic plot") + # Remove the weights from the re-weighted samples axes = fig.get_axes() - axes[0].legend(lines, ["$g(x)$", "$f(x)$"]) - - axes = np.array(axes).reshape((self.ndim, self.ndim)) - - base_alpha = 0.1 - alphas = base_alpha + (1 - base_alpha) * weights - for ii in range(1, self.ndim): - for jj in range(0, self.ndim - 1): - if ii <= jj: - continue - if self.kwargs["mirror_diagnostic_plot"]: - ax = axes[jj, ii] - xsamples = rxs[:, ii] - ysamples = rxs[:, jj] - else: - ax = axes[ii, jj] - xsamples = rxs[:, jj] - ysamples = rxs[:, ii] - sc = ax.scatter( - xsamples, - ysamples, - c=np.log(weights), - alpha=alphas, - edgecolor="none", - vmin=-5, - vmax=0, - ) - cbar = fig.colorbar(sc, ax=axes[0, -1]) - cbar.set_label("Log-weight") + axes[-1].patches[-1].remove() + for ii in range(2, self.ndim + 2): + axes[-ii].collections[-1].remove() + axes[-ii].collections[-1].remove() + lines.append(mpllines.Line2D([0], [0], color=f_color, linestyle=f_ls)) else: - fig, ax = plt.subplots() - ax.hist( - xs, - bins=kwargs["bins"], - color=f_color, - histtype="step", - ls=f_ls, - density=True, - ) - ax.hist( - rxs, - bins=kwargs["bins"], - color=g_color, - histtype="step", - ls=g_ls, - density=True, - ) - ax.set_xlabel(kwargs["labels"][0]) + logger.info("Too few samples to add to the diagnostic plot") + + axes = fig.get_axes() + axes[0].legend(lines, ["$g(x)$", "$f(x)$"]) + + axes = np.array(axes).reshape((self.ndim + 1, self.ndim + 1)) + + base_alpha = 0.1 + alphas = base_alpha + (1 - base_alpha) * weights + for ii in range(1, self.ndim + 1): + for jj in range(0, self.ndim): + if ii <= jj: + continue + if self.kwargs["mirror_diagnostic_plot"]: + ax = axes[jj, ii] + xsamples = rxs[:, ii] + ysamples = rxs[:, jj] + else: + ax = axes[ii, jj] + xsamples = rxs[:, jj] + ysamples = rxs[:, ii] + sc = ax.scatter( + xsamples, + ysamples, + c=np.log(weights), + alpha=alphas, + edgecolor="none", + vmin=-5, + vmax=0, + ) + cbar = fig.colorbar(sc, ax=axes[0, -1]) + cbar.set_label("Log-weight") fig.suptitle(f"Resampling method: {method}") From c4eacc0ea7c62d9ebf7cef22fe8d50346c8f50d9 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 24 Apr 2025 15:46:56 +0100 Subject: [PATCH 24/30] Apply conversion to the antenae pattern function calculation --- bilby/gw/detector/interferometer.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/bilby/gw/detector/interferometer.py b/bilby/gw/detector/interferometer.py index ac0c14886..d6ecf5d9b 100644 --- a/bilby/gw/detector/interferometer.py +++ b/bilby/gw/detector/interferometer.py @@ -309,6 +309,9 @@ def get_detector_response(self, waveform_polarizations, parameters, frequencies= else: mask = np.ones(len(frequencies), dtype=bool) + if "psi" not in parameters: + parameters = generate_all_bbh_parameters(parameters) + signal = {} for mode in waveform_polarizations.keys(): det_response = self.antenna_response( From 778cb25d872b1956c9b4a0c9605e76594dd9a300 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 24 Apr 2025 15:47:43 +0100 Subject: [PATCH 25/30] Fix bug with sorted dataframes passed in to generate_sky_frame_parameters If a pandas dataframe with a non-monotonically increasing index (say if it has been resampled) is passed in, the new_samples data frame that is created will have a different index and when matched NaN's will be introduced. This fixes that issue by setting the column from the values --- bilby/gw/conversion.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index 2a50d7b5a..4a4c1c8b1 100644 --- a/bilby/gw/conversion.py +++ b/bilby/gw/conversion.py @@ -2540,7 +2540,8 @@ def generate_sky_frame_parameters(samples, likelihood): new_samples.append(likelihood.get_sky_frame_parameters()) new_samples = DataFrame(new_samples) for key in new_samples: - samples[key] = new_samples[key] + # Use values to avoid index-sorting (new-samples are automatically sorted properly) + samples[key] = new_samples[key].values def fill_sample(args): From d1fabc6c9aa57867d7b11a983f577b82d0f2ee3e Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 24 Apr 2025 15:49:56 +0100 Subject: [PATCH 26/30] Resolve bug in which for zero-spin aligned models the generation failed Due to the division of 0/0, this produced infs in the calculate tilt parameters. --- bilby/gw/conversion.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index 4a4c1c8b1..a1e2f011e 100644 --- a/bilby/gw/conversion.py +++ b/bilby/gw/conversion.py @@ -236,10 +236,13 @@ def convert_to_lal_binary_black_hole_parameters(parameters): converted_parameters[f"chi_{idx}"] ** 2 + converted_parameters[f"chi_{idx}_in_plane"] ** 2 ) ** 0.5 - converted_parameters[f"cos_tilt_{idx}"] = ( - converted_parameters[f"chi_{idx}"] - / converted_parameters[f"a_{idx}"] - ) + if converted_parameters[f"a_{idx}"] == 0: + converted_parameters[f"cos_tilt_{idx}"] = 0 + else: + converted_parameters[f"cos_tilt_{idx}"] = ( + converted_parameters[f"chi_{idx}"] + / converted_parameters[f"a_{idx}"] + ) elif "a_{}".format(idx) not in original_keys: converted_parameters['a_{}'.format(idx)] = abs( converted_parameters[key]) From ab25090a1b424c7b73ac45b01caec0106877d753 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 24 Apr 2025 15:51:05 +0100 Subject: [PATCH 27/30] Add a new conversion for psi mod pi/2 to remove degenerate solutions --- bilby/gw/conversion.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/bilby/gw/conversion.py b/bilby/gw/conversion.py index a1e2f011e..f5293d39f 100644 --- a/bilby/gw/conversion.py +++ b/bilby/gw/conversion.py @@ -204,7 +204,6 @@ def convert_to_lal_binary_black_hole_parameters(parameters): added_keys: list keys which are added to parameters during function call """ - converted_parameters = parameters.copy() original_keys = list(converted_parameters.keys()) if 'luminosity_distance' not in original_keys: @@ -279,6 +278,10 @@ def convert_to_lal_binary_black_hole_parameters(parameters): - np.sign(np.cos(converted_parameters["theta_jn"])) * converted_parameters["psi"], 2 * np.pi) + + if "psi_mod_pib2" in original_keys: + converted_parameters["psi"] = np.mod(converted_parameters["psi_mod_pib2"], np.pi / 2) + added_keys = [key for key in converted_parameters.keys() if key not in original_keys] From 2b866165a0eb2abc9d6ed5eed2b17465f35fc780 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 24 Apr 2025 21:16:00 +0100 Subject: [PATCH 28/30] Add empty expected outputs --- bilby/core/fisher_sampler.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 6cd841e00..005cdae21 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -101,10 +101,12 @@ def get_expected_outputs(cls, outdir=None, label=None): Returns ------- list - List of file names. + List of file names: empty (resuming not yet implemented) list - List of directory names. Will always be empty for dynesty. + List of directory names. Will always be empty for fisher. """ + return [], [] + raise NotImplementedError() @signal_wrapper From 538da6668b5b547bb23453299408b5b44e2ae717 Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Sat, 26 Apr 2025 01:25:35 -0700 Subject: [PATCH 29/30] Fix minor bus --- bilby/core/fisher_sampler.py | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 005cdae21..1c1564051 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -1,4 +1,5 @@ import datetime +import sys import numpy as np import pandas as pd @@ -122,7 +123,7 @@ def run_sampler(self): fd_eps=self.kwargs["fd_eps"], ) - if self.injection_parameters and "use_injection_for_maxL" in self.kwargs: + if self.injection_parameters and self.kwargs["use_injection_for_maxL"]: sample = { key: self.injection_parameters[key] for key in fisher_mpe.parameter_names @@ -154,7 +155,10 @@ def run_sampler(self): f"Starting sampling in batches of {batch_nsamples} to produce {target_nsamples} samples" ) pbar = tqdm.tqdm( - total=target_nsamples, desc=f"{resample.capitalize()} sampling" + total=target_nsamples, + desc=f"{resample.capitalize()} sampling", + file=sys.stdout, + initial=0, ) while nsamples < target_nsamples: g_samples, g_logl, g_logpi = ( @@ -181,20 +185,24 @@ def run_sampler(self): pbar.set_postfix( { "eff": f"{efficiency:.3f}%", - } + }, + refresh=False, ) - pbar.update(len(samples)) - all_g_samples.append(g_samples) - all_samples.append(samples) - all_logl.append(logl) - all_weights.append(weights) + if len(samples) > self.ndim: + pbar.update(len(samples)) + all_g_samples.append(g_samples) + all_samples.append(samples) + all_logl.append(logl) + all_weights.append(weights) + else: + pbar.update(0) pbar.close() g_samples = pd.concat(all_g_samples) samples = pd.concat(all_samples) - logl = np.concat(all_logl) - weights = np.concat(all_weights) + logl = np.concatenate(all_logl) + weights = np.concatenate(all_weights) efficiency = 100 * len(samples) / len(g_samples) logger.info(f"Finished sampling: total efficiency is {efficiency:0.3f}%") @@ -225,7 +233,7 @@ def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method bins=50, smooth=0.7, max_n_ticks=5, - truths=np.concat((mean, [1])), + truths=np.concatenate((mean, [1])), truth_color="C3", ) @@ -236,9 +244,9 @@ def create_resample_diagnostic(self, samples, raw_samples, mean, weights, method # Create the data array to plot and pass everything to corner xs = samples[self.search_parameter_keys].values - xs = np.concat((xs, np.random.uniform(0, 1, len(xs)).reshape(-1, 1)), axis=1) + xs = np.concatenate((xs, np.random.uniform(0, 1, len(xs)).reshape(-1, 1)), axis=1) rxs = raw_samples[self.search_parameter_keys].values - rxs = np.concat((rxs, weights.reshape((-1, 1))), axis=1) + rxs = np.concatenate((rxs, weights.reshape((-1, 1))), axis=1) # Sort by weight (only for plotting) idxs = np.argsort(weights) From 0a73a880a59aee50c1beb4dab99e0dcddc0b238b Mon Sep 17 00:00:00 2001 From: Gregory Ashton Date: Thu, 15 May 2025 02:26:36 -0700 Subject: [PATCH 30/30] Clean ups --- bilby/core/fisher_sampler.py | 41 ++++++++++++++++++------------------ 1 file changed, 21 insertions(+), 20 deletions(-) diff --git a/bilby/core/fisher_sampler.py b/bilby/core/fisher_sampler.py index 1c1564051..5bb56a979 100644 --- a/bilby/core/fisher_sampler.py +++ b/bilby/core/fisher_sampler.py @@ -114,6 +114,7 @@ def get_expected_outputs(cls, outdir=None, label=None): def run_sampler(self): self.start_time = datetime.datetime.now() + cov_scaling = self.kwargs["cov_scaling"] fisher_mpe = FisherMatrixPosteriorEstimator( likelihood=self.likelihood, @@ -133,7 +134,7 @@ def run_sampler(self): maxL_sample_dict = fisher_mpe.get_maximum_likelihood_sample(sample) mean = np.array(list(maxL_sample_dict.values())) iFIM = fisher_mpe.calculate_iFIM(maxL_sample_dict) - cov = self.kwargs["cov_scaling"] * iFIM + cov = cov_scaling * iFIM msg = "Generation-distribution:\n " + "\n ".join( [ @@ -161,7 +162,7 @@ def run_sampler(self): initial=0, ) while nsamples < target_nsamples: - g_samples, g_logl, g_logpi = ( + g_samples, g_logl, g_logpi, discard_inef = ( self._draw_samples_from_generating_distribution( mean, cov, fisher_mpe, batch_nsamples ) @@ -181,10 +182,16 @@ def run_sampler(self): logl = g_logl weights = np.ones_like(logl) + #if efficiency < 10: + # cov_scaling *= 0.9 + # cov = cov_scaling * iFIM + nsamples += len(samples) pbar.set_postfix( { "eff": f"{efficiency:.3f}%", + "de": f"{discard_inef:.1f}%", + "cs": f"{cov_scaling:.2f}", }, refresh=False, ) @@ -370,15 +377,15 @@ def _draw_samples_from_generating_distribution( logger.debug("Calculating the likelihood and priors") logpi = self.priors.ln_prob(samples, axis=0) - outside_prior_count = np.sum(np.isinf(logpi)) - if outside_prior_count == 0: - logl = fisher_mpe.log_likelihood_from_array(samples.values.T) - else: - for ii, rs in tqdm.tqdm(samples.iterrows(), total=len(samples)): - if np.isinf(logpi[ii]): - logl.append(-np.inf) - else: - logl.append(fisher_mpe.log_likelihood(rs.to_dict())) + logl = np.full_like(logpi, -np.inf) + + in_prior = ~np.isinf(logpi) + outside_prior_count = np.sum(~in_prior) + discard_inef = 100 * outside_prior_count / len(samples) + + logl[in_prior] = fisher_mpe.log_likelihood_from_array( + samples.values[in_prior].T + ) if outside_prior_count == len(samples): msg = "Sampling has failed: no viable samples left" @@ -386,17 +393,11 @@ def _draw_samples_from_generating_distribution( raise SamplerError(msg) else: logger.info(msg) - elif outside_prior_count > 0: - logger.info( - f"Discarding {100 * outside_prior_count / len(samples):0.3f}% of samples that" - " fall outside the prior" - ) else: logger.debug("No samples outside the prior bounds") logpi = np.real(np.array(logpi)) - logl = np.array(logl) - return samples, logl, logpi + return samples, logl, logpi, discard_inef def _calculate_weights(self, g_samples, g_logl, g_logpi, mean, cov): g_logl_norm = multivariate_normal.logpdf(g_samples, mean=mean, cov=cov) @@ -430,7 +431,7 @@ def _importance_sample(self, g_samples, g_logl, weights): if self.kwargs["fail_on_error"]: raise SamplerError(msg) else: - logger.info(msg) + logger.debug(msg) return samples, logl @@ -447,6 +448,6 @@ def _rejection_sample(self, g_samples, g_logl, weights): if self.kwargs["fail_on_error"]: raise SamplerError(msg) else: - logger.info(msg) + logger.debug(msg) return samples, logl