From a3abc9ad60baed683dcd9a7d6ccc0fdead8163a7 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Wed, 9 Aug 2023 11:36:43 +0100 Subject: [PATCH 01/15] add autodiff machinery to jax and pytorch backends --- src/pyhf/optimize/mixins.py | 37 ++++++++++++++++++++++--- src/pyhf/tensor/jax_backend.py | 40 ++++++++++++++++++++++++++- src/pyhf/tensor/numpy_backend.py | 26 +++++++++++++++++ src/pyhf/tensor/pytorch_backend.py | 39 ++++++++++++++++++++++++++ src/pyhf/tensor/tensorflow_backend.py | 27 ++++++++++++++++++ 5 files changed, 164 insertions(+), 5 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index c0e689c512..8d10a4d9df 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -49,6 +49,18 @@ def _internal_minimize( do_grad=do_grad, par_names=par_names, ) + + # so we can check if valid for uncertainty calc later + try: + minimizer_name = minimizer.name + if minimizer_name == "minuit": + self.using_minuit = True + else: + self.using_minuit = False + except AttributeError as e: + self.using_minuit = False + + result = self._minimize( minimizer, func, @@ -66,7 +78,7 @@ def _internal_minimize( raise exceptions.FailedMinimization(result) return result - def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=False): + def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=False, uncertainties=None, hess_inv=None): """ Post-process the fit result. @@ -79,7 +91,7 @@ def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=Fal fitted_pars = stitch_pars(tensorlib.astensor(fitresult.x)) # check if uncertainties were provided (and stitch just in case) - uncertainties = getattr(fitresult, 'unc', None) + uncertainties = getattr(fitresult, 'unc', None) or uncertainties if uncertainties is not None: # extract number of fixed parameters num_fixed_pars = len(fitted_pars) - len(fitresult.x) @@ -88,7 +100,8 @@ def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=Fal # https://github.com/scikit-hep/iminuit/issues/762 # https://github.com/scikit-hep/pyhf/issues/1918 # https://github.com/scikit-hep/cabinetry/pull/346 - uncertainties = np.where(fitresult.minuit.fixed, 0.0, uncertainties) + if self.using_minuit: + uncertainties = np.where(fitresult.minuit.fixed, 0.0, uncertainties) # stitch in zero-uncertainty for fixed values uncertainties = stitch_pars( @@ -112,11 +125,19 @@ def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=Fal ] correlations = tensorlib.stack(stitched_rows, axis=1) + if hess_inv is not None: + if hasattr(fitresult, 'hess_inv'): + if fitresult.hess_inv is None: + fitresult.hess_inv = hess_inv + else: + fitresult.hess_inv = hess_inv + fitresult.x = fitted_pars fitresult.fun = tensorlib.astensor(fitresult.fun) fitresult.unc = uncertainties fitresult.corr = correlations + return fitresult def minimize( @@ -193,8 +214,16 @@ def minimize( result = self._internal_minimize( **minimizer_kwargs, options=kwargs, par_names=par_names ) + + # compute uncertainties with automatic differentiation + if not self.using_minuit and tensorlib.name in ['jax', 'pytorch']: + hess_inv = tensorlib.fisher_cov(pdf, result.x, data) + uncertainties = tensorlib.sqrt(tensorlib.diagonal(hess_inv)) + else: + hess_inv = None + uncertainties = None result = self._internal_postprocess( - result, stitch_pars, return_uncertainties=return_uncertainties + result, stitch_pars, pdf, return_uncertainties=return_uncertainties, uncertainties=uncertainties, hess_inv=hess_inv ) _returns = [result.x] diff --git a/src/pyhf/tensor/jax_backend.py b/src/pyhf/tensor/jax_backend.py index 6ae53b7bec..aad2ca2bd4 100644 --- a/src/pyhf/tensor/jax_backend.py +++ b/src/pyhf/tensor/jax_backend.py @@ -2,7 +2,7 @@ config.update('jax_enable_x64', True) -from jax import Array +from jax import Array, hessian import jax.numpy as jnp from jax.scipy.special import gammaln, xlogy from jax.scipy import special @@ -622,3 +622,41 @@ def transpose(self, tensor_in): .. versionadded:: 0.7.0 """ return tensor_in.transpose() + + + def fisher_cov(self, model, pars, data): + """Calculates the inverse of the Fisher information matrix to estimate the covariance of the maximum likelihood estimate. + See the Cramér-Rao bound for more details on the derivation of this. + + Args: + model (:obj:`pyhf.pdf.Model`): The statistical model adhering to the schema ``model.json``. + pars (:obj:`tensor`): The (mle) model parameters at which to evaluate the uncertainty. + data (:obj:`tensor`): The observed data. + + Returns: + JAX ndarray: The covariance matrix of the maximum likelihood estimate. + """ + return jnp.linalg.inv(-hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data)) + + + def diagonal(self, tensor_in): + """ + Return the diagonal elements of the tensor. + + Example: + >>> import pyhf + >>> pyhf.set_backend("jax") + >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) + >>> tensor + tensor([[1., 0.], + [0., 1.]]) + >>> pyhf.tensorlib.diagonal(tensor) + tensor([1., 1.]) + + Args: + tensor_in (:obj:`tensor`): The input tensor object. + + Returns: + JAX ndarray: The diagonal of the input tensor. + """ + return jnp.diag(tensor_in) diff --git a/src/pyhf/tensor/numpy_backend.py b/src/pyhf/tensor/numpy_backend.py index 4623ea5797..79276d04f1 100644 --- a/src/pyhf/tensor/numpy_backend.py +++ b/src/pyhf/tensor/numpy_backend.py @@ -647,3 +647,29 @@ def transpose(self, tensor_in: Tensor[T]) -> ArrayLike: .. versionadded:: 0.7.0 """ return tensor_in.transpose() + + + def fisher_cov(self, model, pars, data): + raise NotImplementedError + + + def diagonal(self, tensor_in): + """Return the diagonal elements of the tensor. + + Example: + >>> import pyhf + >>> pyhf.set_backend("numpy") + >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) + >>> tensor + tensor([[1., 0.], + [0., 1.]]) + >>> pyhf.tensorlib.diagonal(tensor) + tensor([1., 1.]) + + Args: + tensor_in (:obj:`tensor`): The input tensor object. + + Returns: + :class:`numpy.ndarray`: The diagonal of the input tensor. + """ + return np.diag(tensor_in) \ No newline at end of file diff --git a/src/pyhf/tensor/pytorch_backend.py b/src/pyhf/tensor/pytorch_backend.py index fb0cd4ad52..76dd80e850 100644 --- a/src/pyhf/tensor/pytorch_backend.py +++ b/src/pyhf/tensor/pytorch_backend.py @@ -1,6 +1,7 @@ """PyTorch Tensor Library Module.""" import torch import torch.autograd +from torch.func import hessian from torch.distributions.utils import broadcast_all import logging import math @@ -625,3 +626,41 @@ def transpose(self, tensor_in): .. versionadded:: 0.7.0 """ return tensor_in.transpose(0, 1) + + + def fisher_cov(self, model, pars, data): + """Calculates the inverse of the Fisher information matrix to estimate the covariance of the maximum likelihood estimate. + See the Cramér-Rao bound for more details on the derivation of this. + + Args: + model (:obj:`pyhf.pdf.Model`): The statistical model adhering to the schema ``model.json``. + pars (:obj:`tensor`): The (mle) model parameters at which to evaluate the uncertainty. + data (:obj:`tensor`): The observed data. + + Returns: + PyTorch FloatTensor: The covariance matrix of the maximum likelihood estimate. + """ + return torch.linalg.inv(-hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data)) + + + def diagonal(self, tensor_in): + """ + Return the diagonal elements of the tensor. + + Example: + >>> import pyhf + >>> pyhf.set_backend("pytorch") + >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) + >>> tensor + tensor([[1., 0.], + [0., 1.]]) + >>> pyhf.tensorlib.diagonal(tensor) + tensor([1., 1.]) + + Args: + tensor_in (:obj:`tensor`): The input tensor object. + + Returns: + PyTorch FloatTensor: The diagonal of the input tensor. + """ + return torch.diagonal(tensor_in) \ No newline at end of file diff --git a/src/pyhf/tensor/tensorflow_backend.py b/src/pyhf/tensor/tensorflow_backend.py index 1e51a2a885..28c0107383 100644 --- a/src/pyhf/tensor/tensorflow_backend.py +++ b/src/pyhf/tensor/tensorflow_backend.py @@ -722,3 +722,30 @@ def transpose(self, tensor_in): .. versionadded:: 0.7.0 """ return tf.transpose(tensor_in) + + + def fisher_cov(self, model, pars, data): + raise NotImplementedError + + + def diagonal(self, tensor_in): + """Return the diagonal elements of the tensor. + + Example: + >>> import pyhf + >>> pyhf.set_backend("tensorflow") + >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) + >>> tensor + tensor([[1., 0.], + [0., 1.]]) + >>> pyhf.tensorlib.diagonal(tensor) + tensor([1., 1.]) + + Args: + tensor_in (:obj:`tensor`): The input tensor object. + + Returns: + TensorFlow Tensor: The diagonal elements of the input tensor. + + """ + return tf.linalg.diag(tensor_in) \ No newline at end of file From eee1b3fc0da5c32e8f2f42160f1b68943163c847 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Wed, 9 Aug 2023 11:52:49 +0100 Subject: [PATCH 02/15] run linting --- src/pyhf/optimize/mixins.py | 20 +++++++++++++++----- src/pyhf/tensor/jax_backend.py | 6 +++--- src/pyhf/tensor/numpy_backend.py | 21 ++++++++++++++------- src/pyhf/tensor/pytorch_backend.py | 8 ++++---- src/pyhf/tensor/tensorflow_backend.py | 8 +++----- 5 files changed, 39 insertions(+), 24 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 8d10a4d9df..728d80dc60 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -57,9 +57,8 @@ def _internal_minimize( self.using_minuit = True else: self.using_minuit = False - except AttributeError as e: + except AttributeError: self.using_minuit = False - result = self._minimize( minimizer, @@ -78,7 +77,14 @@ def _internal_minimize( raise exceptions.FailedMinimization(result) return result - def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=False, uncertainties=None, hess_inv=None): + def _internal_postprocess( + self, + fitresult, + stitch_pars, + return_uncertainties=False, + uncertainties=None, + hess_inv=None, + ): """ Post-process the fit result. @@ -137,7 +143,6 @@ def _internal_postprocess(self, fitresult, stitch_pars, return_uncertainties=Fal fitresult.unc = uncertainties fitresult.corr = correlations - return fitresult def minimize( @@ -223,7 +228,12 @@ def minimize( hess_inv = None uncertainties = None result = self._internal_postprocess( - result, stitch_pars, pdf, return_uncertainties=return_uncertainties, uncertainties=uncertainties, hess_inv=hess_inv + result, + stitch_pars, + pdf, + return_uncertainties=return_uncertainties, + uncertainties=uncertainties, + hess_inv=hess_inv, ) _returns = [result.x] diff --git a/src/pyhf/tensor/jax_backend.py b/src/pyhf/tensor/jax_backend.py index aad2ca2bd4..d0dd00ca64 100644 --- a/src/pyhf/tensor/jax_backend.py +++ b/src/pyhf/tensor/jax_backend.py @@ -623,7 +623,6 @@ def transpose(self, tensor_in): """ return tensor_in.transpose() - def fisher_cov(self, model, pars, data): """Calculates the inverse of the Fisher information matrix to estimate the covariance of the maximum likelihood estimate. See the Cramér-Rao bound for more details on the derivation of this. @@ -636,8 +635,9 @@ def fisher_cov(self, model, pars, data): Returns: JAX ndarray: The covariance matrix of the maximum likelihood estimate. """ - return jnp.linalg.inv(-hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data)) - + return jnp.linalg.inv( + -hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data) + ) def diagonal(self, tensor_in): """ diff --git a/src/pyhf/tensor/numpy_backend.py b/src/pyhf/tensor/numpy_backend.py index 79276d04f1..fc46164764 100644 --- a/src/pyhf/tensor/numpy_backend.py +++ b/src/pyhf/tensor/numpy_backend.py @@ -2,7 +2,16 @@ from __future__ import annotations import logging -from typing import TYPE_CHECKING, Callable, Generic, Mapping, Sequence, TypeVar, Union +from typing import ( + TYPE_CHECKING, + Callable, + Generic, + Mapping, + Sequence, + TypeVar, + Union, + Any, +) import numpy as np @@ -648,12 +657,10 @@ def transpose(self, tensor_in: Tensor[T]) -> ArrayLike: """ return tensor_in.transpose() - - def fisher_cov(self, model, pars, data): + def fisher_cov(self, model: Any, pars: Tensor[T], data: Tensor[T]) -> ArrayLike: raise NotImplementedError - - - def diagonal(self, tensor_in): + + def diagonal(self, tensor_in: Tensor[T]) -> ArrayLike: """Return the diagonal elements of the tensor. Example: @@ -672,4 +679,4 @@ def diagonal(self, tensor_in): Returns: :class:`numpy.ndarray`: The diagonal of the input tensor. """ - return np.diag(tensor_in) \ No newline at end of file + return np.diag(tensor_in) diff --git a/src/pyhf/tensor/pytorch_backend.py b/src/pyhf/tensor/pytorch_backend.py index 76dd80e850..f597237093 100644 --- a/src/pyhf/tensor/pytorch_backend.py +++ b/src/pyhf/tensor/pytorch_backend.py @@ -627,7 +627,6 @@ def transpose(self, tensor_in): """ return tensor_in.transpose(0, 1) - def fisher_cov(self, model, pars, data): """Calculates the inverse of the Fisher information matrix to estimate the covariance of the maximum likelihood estimate. See the Cramér-Rao bound for more details on the derivation of this. @@ -640,8 +639,9 @@ def fisher_cov(self, model, pars, data): Returns: PyTorch FloatTensor: The covariance matrix of the maximum likelihood estimate. """ - return torch.linalg.inv(-hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data)) - + return torch.linalg.inv( + -hessian(lambda pars, data: model.logpdf(pars, data)[0])(pars, data) + ) def diagonal(self, tensor_in): """ @@ -663,4 +663,4 @@ def diagonal(self, tensor_in): Returns: PyTorch FloatTensor: The diagonal of the input tensor. """ - return torch.diagonal(tensor_in) \ No newline at end of file + return torch.diagonal(tensor_in) diff --git a/src/pyhf/tensor/tensorflow_backend.py b/src/pyhf/tensor/tensorflow_backend.py index 28c0107383..0a7471243b 100644 --- a/src/pyhf/tensor/tensorflow_backend.py +++ b/src/pyhf/tensor/tensorflow_backend.py @@ -723,11 +723,9 @@ def transpose(self, tensor_in): """ return tf.transpose(tensor_in) - def fisher_cov(self, model, pars, data): raise NotImplementedError - - + def diagonal(self, tensor_in): """Return the diagonal elements of the tensor. @@ -746,6 +744,6 @@ def diagonal(self, tensor_in): Returns: TensorFlow Tensor: The diagonal elements of the input tensor. - + """ - return tf.linalg.diag(tensor_in) \ No newline at end of file + return tf.linalg.diag(tensor_in) From 123f2a0dbd3d186e324869b9b98a3cbcb574e1c3 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Wed, 9 Aug 2023 12:30:15 +0100 Subject: [PATCH 03/15] fix torch by casting res to tensor --- src/pyhf/optimize/mixins.py | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 728d80dc60..645debca1a 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -50,16 +50,6 @@ def _internal_minimize( par_names=par_names, ) - # so we can check if valid for uncertainty calc later - try: - minimizer_name = minimizer.name - if minimizer_name == "minuit": - self.using_minuit = True - else: - self.using_minuit = False - except AttributeError: - self.using_minuit = False - result = self._minimize( minimizer, func, @@ -81,6 +71,7 @@ def _internal_postprocess( self, fitresult, stitch_pars, + using_minuit, return_uncertainties=False, uncertainties=None, hess_inv=None, @@ -106,7 +97,7 @@ def _internal_postprocess( # https://github.com/scikit-hep/iminuit/issues/762 # https://github.com/scikit-hep/pyhf/issues/1918 # https://github.com/scikit-hep/cabinetry/pull/346 - if self.using_minuit: + if using_minuit: uncertainties = np.where(fitresult.minuit.fixed, 0.0, uncertainties) # stitch in zero-uncertainty for fixed values @@ -189,6 +180,22 @@ def minimize( - minimum (:obj:`float`): if ``return_fitted_val`` flagged, return minimized objective value - result (:class:`scipy.optimize.OptimizeResult`): if ``return_result_obj`` flagged """ + # literally just for the minimizer name to check if we're using minuit + _minimizer = self._get_minimizer( + lambda x: x, + [0], + [0, 1], + ) + + # so we can check if valid for uncertainty calc later + if hasattr(_minimizer, "name"): + if _minimizer.name == "minuit": + using_minuit = True + else: + using_minuit = False + else: + using_minuit = False + # Configure do_grad based on backend "automagically" if not set by user tensorlib, _ = get_backend() do_grad = tensorlib.default_do_grad if do_grad is None else do_grad @@ -221,8 +228,8 @@ def minimize( ) # compute uncertainties with automatic differentiation - if not self.using_minuit and tensorlib.name in ['jax', 'pytorch']: - hess_inv = tensorlib.fisher_cov(pdf, result.x, data) + if not using_minuit and tensorlib.name in ['jax', 'pytorch']: + hess_inv = tensorlib.fisher_cov(pdf, tensorlib.astensor(result.x), data) uncertainties = tensorlib.sqrt(tensorlib.diagonal(hess_inv)) else: hess_inv = None @@ -230,7 +237,7 @@ def minimize( result = self._internal_postprocess( result, stitch_pars, - pdf, + using_minuit, return_uncertainties=return_uncertainties, uncertainties=uncertainties, hess_inv=hess_inv, From d95c4337582f83ae46a2126b9af6543dca85b9d2 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Wed, 9 Aug 2023 15:24:11 +0100 Subject: [PATCH 04/15] add tensorflow --- src/pyhf/optimize/mixins.py | 2 +- src/pyhf/tensor/tensorflow_backend.py | 20 +++++++++++++++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 645debca1a..ae519c82aa 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -228,7 +228,7 @@ def minimize( ) # compute uncertainties with automatic differentiation - if not using_minuit and tensorlib.name in ['jax', 'pytorch']: + if not using_minuit and tensorlib.name in ['tensorflow', 'jax', 'pytorch']: hess_inv = tensorlib.fisher_cov(pdf, tensorlib.astensor(result.x), data) uncertainties = tensorlib.sqrt(tensorlib.diagonal(hess_inv)) else: diff --git a/src/pyhf/tensor/tensorflow_backend.py b/src/pyhf/tensor/tensorflow_backend.py index 0a7471243b..3989b3a9d9 100644 --- a/src/pyhf/tensor/tensorflow_backend.py +++ b/src/pyhf/tensor/tensorflow_backend.py @@ -724,7 +724,25 @@ def transpose(self, tensor_in): return tf.transpose(tensor_in) def fisher_cov(self, model, pars, data): - raise NotImplementedError + """Calculates the inverse of the Fisher information matrix to estimate the covariance of the maximum likelihood estimate. + See the Cramér-Rao bound for more details on the derivation of this. + + Args: + model (:obj:`pyhf.pdf.Model`): The statistical model adhering to the schema ``model.json``. + pars (:obj:`tensor`): The (mle) model parameters at which to evaluate the uncertainty. + data (:obj:`tensor`): The observed data. + + Returns: + TensorFlow Tensor: The covariance matrix of the maximum likelihood estimate. + """ + with tf.GradientTape() as t2: + t2.watch(pars) + with tf.GradientTape() as t1: + t1.watch(pars) + lhood = model.logpdf(pars, data)[0] + g = t1.gradient(lhood, pars) + hess = t2.jacobian(g, pars) + return tf.linalg.inv(-hess) def diagonal(self, tensor_in): """Return the diagonal elements of the tensor. From 9b99a5fa27361f614f0975eb8b5c66b1bbff9417 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Fri, 11 Aug 2023 10:20:21 +0100 Subject: [PATCH 05/15] stitch in fixed pars before calculating uncerts --- src/pyhf/optimize/mixins.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index ae519c82aa..a04d1e0711 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -229,7 +229,9 @@ def minimize( # compute uncertainties with automatic differentiation if not using_minuit and tensorlib.name in ['tensorflow', 'jax', 'pytorch']: - hess_inv = tensorlib.fisher_cov(pdf, tensorlib.astensor(result.x), data) + # stitch in missing parameters (e.g. fixed parameters) + fitted_pars = stitch_pars(tensorlib.astensor(result.x)) + hess_inv = tensorlib.fisher_cov(pdf, fitted_pars, data) uncertainties = tensorlib.sqrt(tensorlib.diagonal(hess_inv)) else: hess_inv = None From e6281eeed775e12924688ead17906f4f16fef8bf Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Fri, 11 Aug 2023 15:22:26 +0100 Subject: [PATCH 06/15] patch in fisher_cov as identity for tests --- src/pyhf/optimize/mixins.py | 11 ++++++++--- tests/test_optim.py | 20 ++++++++++++++++++-- 2 files changed, 26 insertions(+), 5 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index a04d1e0711..6cfccb5b9b 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -150,6 +150,7 @@ def minimize( return_correlations=False, do_grad=None, do_stitch=False, + skip_autodiff_uncerts=False, **kwargs, ): """ @@ -228,10 +229,14 @@ def minimize( ) # compute uncertainties with automatic differentiation - if not using_minuit and tensorlib.name in ['tensorflow', 'jax', 'pytorch']: + if ( + not using_minuit + and tensorlib.name in ['tensorflow', 'jax', 'pytorch'] + and not skip_autodiff_uncerts + ): # stitch in missing parameters (e.g. fixed parameters) - fitted_pars = stitch_pars(tensorlib.astensor(result.x)) - hess_inv = tensorlib.fisher_cov(pdf, fitted_pars, data) + all_pars = stitch_pars(tensorlib.astensor(result.x)) + hess_inv = tensorlib.fisher_cov(pdf, all_pars, data) uncertainties = tensorlib.sqrt(tensorlib.diagonal(hess_inv)) else: hess_inv = None diff --git a/tests/test_optim.py b/tests/test_optim.py index b94b722cf8..54cef2fcae 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -123,8 +123,24 @@ def test_minimize_do_grad_autoconfig(mocker, backend, backend_new): mocker.patch.object(OptimizerMixin, '_internal_minimize') mocker.patch.object(OptimizerMixin, '_internal_postprocess') + # patch out uncert calculations + import_dict = dict( + jax=pyhf.tensor.jax_backend, # .jax_backend, + pytorch=pyhf.tensor.pytorch_backend, # .pytorch_backend, + tensorflow=pyhf.tensor.tensorflow_backend, # .tensorflow_backend, + numpy=pyhf.tensor.numpy_backend, + ) + + def make_backend(backend): + class MockedBackend(import_dict[backend]): + def fisher_cov(self, *args, **kwargs): + return import_dict[backend]().astensor([[1, 0], [0, 1]]) + + return MockedBackend() + # start with first backend - pyhf.set_backend(backend, 'scipy') + pyhf.set_backend(make_backend(backend), 'scipy') + m = pyhf.simplemodels.uncorrelated_background([50.0], [100.0], [10.0]) data = pyhf.tensorlib.astensor([125.0] + m.config.auxdata) @@ -135,7 +151,7 @@ def test_minimize_do_grad_autoconfig(mocker, backend, backend_new): assert shim.call_args[1]['do_grad'] != pyhf.tensorlib.default_do_grad # now switch to new backend and see what happens - pyhf.set_backend(backend_new) + pyhf.set_backend(make_backend(backend_new)) m = pyhf.simplemodels.uncorrelated_background([50.0], [100.0], [10.0]) data = pyhf.tensorlib.astensor([125.0] + m.config.auxdata) From 11f0a0f09018d315b1a1f6b55b7ef5ca3928b555 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Mon, 14 Aug 2023 10:35:28 +0100 Subject: [PATCH 07/15] lint --- src/pyhf/optimize/mixins.py | 45 ++++++++++++++------------- src/pyhf/optimize/opt_minuit.py | 2 ++ src/pyhf/tensor/tensorflow_backend.py | 2 +- tests/test_optim.py | 32 +++++++++++++++++-- 4 files changed, 56 insertions(+), 25 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 6cfccb5b9b..eb718ea40b 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -7,8 +7,11 @@ from pyhf.optimize.common import shim from pyhf.tensor.manager import get_backend + log = logging.getLogger(__name__) +__all__ = ("OptimizerMixin",) + class OptimizerMixin: """Mixin Class to build optimizers.""" @@ -75,6 +78,8 @@ def _internal_postprocess( return_uncertainties=False, uncertainties=None, hess_inv=None, + fixed_vals=None, + init_pars=None, ): """ Post-process the fit result. @@ -97,9 +102,18 @@ def _internal_postprocess( # https://github.com/scikit-hep/iminuit/issues/762 # https://github.com/scikit-hep/pyhf/issues/1918 # https://github.com/scikit-hep/cabinetry/pull/346 - if using_minuit: - uncertainties = np.where(fitresult.minuit.fixed, 0.0, uncertainties) - + if fixed_vals is not None: # check for fixed vals + if using_minuit: + uncertainties = np.where(fitresult.minuit.fixed, 0.0, uncertainties) + else: + fixed_bools = [False] * len(init_pars) + for index, _ in fixed_vals: + fixed_bools[index] = True + uncertainties = tensorlib.where( + tensorlib.astensor(fixed_bools, dtype="bool"), + tensorlib.astensor(0.0), + uncertainties, + ) # stitch in zero-uncertainty for fixed values uncertainties = stitch_pars( tensorlib.astensor(uncertainties), @@ -150,7 +164,6 @@ def minimize( return_correlations=False, do_grad=None, do_stitch=False, - skip_autodiff_uncerts=False, **kwargs, ): """ @@ -182,20 +195,8 @@ def minimize( - result (:class:`scipy.optimize.OptimizeResult`): if ``return_result_obj`` flagged """ # literally just for the minimizer name to check if we're using minuit - _minimizer = self._get_minimizer( - lambda x: x, - [0], - [0, 1], - ) - # so we can check if valid for uncertainty calc later - if hasattr(_minimizer, "name"): - if _minimizer.name == "minuit": - using_minuit = True - else: - using_minuit = False - else: - using_minuit = False + using_minuit = hasattr(self, "name") and self.name == "minuit" # Configure do_grad based on backend "automagically" if not set by user tensorlib, _ = get_backend() @@ -229,11 +230,7 @@ def minimize( ) # compute uncertainties with automatic differentiation - if ( - not using_minuit - and tensorlib.name in ['tensorflow', 'jax', 'pytorch'] - and not skip_autodiff_uncerts - ): + if not using_minuit and tensorlib.name in ['tensorflow', 'jax', 'pytorch']: # stitch in missing parameters (e.g. fixed parameters) all_pars = stitch_pars(tensorlib.astensor(result.x)) hess_inv = tensorlib.fisher_cov(pdf, all_pars, data) @@ -241,6 +238,8 @@ def minimize( else: hess_inv = None uncertainties = None + + # uncerts are set to 0 in here for fixed pars result = self._internal_postprocess( result, stitch_pars, @@ -248,6 +247,8 @@ def minimize( return_uncertainties=return_uncertainties, uncertainties=uncertainties, hess_inv=hess_inv, + fixed_vals=fixed_vals, + init_pars=init_pars, ) _returns = [result.x] diff --git a/src/pyhf/optimize/opt_minuit.py b/src/pyhf/optimize/opt_minuit.py index 01424e8ca7..e78b81a9f2 100644 --- a/src/pyhf/optimize/opt_minuit.py +++ b/src/pyhf/optimize/opt_minuit.py @@ -4,6 +4,8 @@ import scipy import iminuit +__all__ = ("minuit_optimizer",) + class minuit_optimizer(OptimizerMixin): """ diff --git a/src/pyhf/tensor/tensorflow_backend.py b/src/pyhf/tensor/tensorflow_backend.py index 3989b3a9d9..75b8510bf2 100644 --- a/src/pyhf/tensor/tensorflow_backend.py +++ b/src/pyhf/tensor/tensorflow_backend.py @@ -764,4 +764,4 @@ def diagonal(self, tensor_in): TensorFlow Tensor: The diagonal elements of the input tensor. """ - return tf.linalg.diag(tensor_in) + return tf.linalg.diag_part(tensor_in) diff --git a/tests/test_optim.py b/tests/test_optim.py index 54cef2fcae..48ae250f8b 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -379,8 +379,8 @@ def test_optim_with_value(backend, source, spec, mu): @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.only_numpy_minuit -def test_optim_uncerts(backend, source, spec, mu): - pdf = pyhf.Model(spec, poi_name="mu") +def test_optim_uncerts_minuit(backend, source, spec, mu): + pdf = pyhf.Model(spec) data = source['bindata']['data'] + pdf.config.auxdata init_pars = pdf.config.suggested_init() @@ -404,6 +404,34 @@ def test_optim_uncerts(backend, source, spec, mu): assert pytest.approx([0.26418431, 0.0]) == pyhf.tensorlib.tolist(result[:, 1]) +@pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) +@pytest.mark.skip_numpy +def test_optim_uncerts_autodiff(backend, source, spec, mu): + pdf = pyhf.Model(spec) + data = source['bindata']['data'] + pdf.config.auxdata + + init_pars = pdf.config.suggested_init() + par_bounds = pdf.config.suggested_bounds() + + optim = pyhf.optimizer + + result = optim.minimize(pyhf.infer.mle.twice_nll, data, pdf, init_pars, par_bounds) + assert pyhf.tensorlib.tolist(result) + + result = optim.minimize( + pyhf.infer.mle.twice_nll, + data, + pdf, + init_pars, + par_bounds, + fixed_vals=[(pdf.config.poi_index, mu)], + return_uncertainties=True, + ) + assert result.shape == (2, 2) + # TODO: add proper numerical test for autodiff uncerts + # assert pytest.approx([0.26418431, 0.0]) == pyhf.tensorlib.tolist(result[:, 1]) + + @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.only_numpy_minuit def test_optim_correlations(backend, source, spec, mu): From d7f3e9d83f54a878954b271f58d35d9cb4f8975b Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Mon, 14 Aug 2023 11:57:54 +0100 Subject: [PATCH 08/15] add correlation calculations, and correctly zero out covariance and correlation matricies for fixed parameters --- src/pyhf/optimize/mixins.py | 56 +++++++++++++++++++++++++++++-------- tests/test_optim.py | 37 ++++++++++++++++++++++-- 2 files changed, 78 insertions(+), 15 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index eb718ea40b..0504ae1d18 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -78,6 +78,7 @@ def _internal_postprocess( return_uncertainties=False, uncertainties=None, hess_inv=None, + calc_correlations=False, fixed_vals=None, init_pars=None, ): @@ -122,31 +123,59 @@ def _internal_postprocess( if return_uncertainties: fitted_pars = tensorlib.stack([fitted_pars, uncertainties], axis=1) - correlations = getattr(fitresult, 'corr', None) - if correlations is not None: + cov = getattr(fitresult, 'hess_inv', None) + if cov is None and hess_inv is not None: + cov = hess_inv + + # we also need to edit the covariance matrix to zero-out uncertainties! + if fixed_vals is not None: + if using_minuit: + fixed_bools = fitresult.minuit.fixed + else: + fixed_bools = [False] * len(init_pars) + # Convert fixed_bools to a numpy array and reshape to make it a column vector + fixed_mask = tensorlib.reshape( + tensorlib.astensor(fixed_bools, dtype="bool"), (-1, 1) + ) + # Create 2D masks for rows and columns + row_mask = fixed_mask + col_mask = tensorlib.transpose(fixed_mask) + + # Use logical OR to combine the masks + final_mask = row_mask | col_mask + + # Use np.where to set elements of the covariance matrix to 0 where the mask is True + cov = tensorlib.where( + final_mask, tensorlib.astensor(0.0), tensorlib.astensor(cov) + ) + + correlations_from_fit = getattr(fitresult, 'corr', None) + if correlations_from_fit is None and calc_correlations: + correlations_from_fit = cov / tensorlib.outer(uncertainties, uncertainties) + correlations_from_fit = tensorlib.where( + tensorlib.isfinite(correlations_from_fit), + correlations_from_fit, + tensorlib.astensor(0.0), + ) + + if correlations_from_fit is not None: _zeros = tensorlib.zeros(num_fixed_pars) # possibly a more elegant way to do this stitched_columns = [ stitch_pars(tensorlib.astensor(column), stitch_with=_zeros) - for column in zip(*correlations) + for column in zip(*correlations_from_fit) ] stitched_rows = [ stitch_pars(tensorlib.astensor(row), stitch_with=_zeros) for row in zip(*stitched_columns) ] - correlations = tensorlib.stack(stitched_rows, axis=1) - - if hess_inv is not None: - if hasattr(fitresult, 'hess_inv'): - if fitresult.hess_inv is None: - fitresult.hess_inv = hess_inv - else: - fitresult.hess_inv = hess_inv + correlations_from_fit = tensorlib.stack(stitched_rows, axis=1) fitresult.x = fitted_pars fitresult.fun = tensorlib.astensor(fitresult.fun) fitresult.unc = uncertainties - fitresult.corr = correlations + fitresult.hess_inv = cov + fitresult.corr = correlations_from_fit return fitresult @@ -235,9 +264,11 @@ def minimize( all_pars = stitch_pars(tensorlib.astensor(result.x)) hess_inv = tensorlib.fisher_cov(pdf, all_pars, data) uncertainties = tensorlib.sqrt(tensorlib.diagonal(hess_inv)) + calc_correlations = True else: hess_inv = None uncertainties = None + calc_correlations = False # uncerts are set to 0 in here for fixed pars result = self._internal_postprocess( @@ -247,6 +278,7 @@ def minimize( return_uncertainties=return_uncertainties, uncertainties=uncertainties, hess_inv=hess_inv, + calc_correlations=calc_correlations, fixed_vals=fixed_vals, init_pars=init_pars, ) diff --git a/tests/test_optim.py b/tests/test_optim.py index 48ae250f8b..25e0016424 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -428,14 +428,45 @@ def test_optim_uncerts_autodiff(backend, source, spec, mu): return_uncertainties=True, ) assert result.shape == (2, 2) - # TODO: add proper numerical test for autodiff uncerts + # TODO: add proper numerical test for autodiff uncerts (does not match minuit at all) # assert pytest.approx([0.26418431, 0.0]) == pyhf.tensorlib.tolist(result[:, 1]) @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.only_numpy_minuit -def test_optim_correlations(backend, source, spec, mu): - pdf = pyhf.Model(spec, poi_name="mu") +def test_optim_correlations_minuit(backend, source, spec, mu): + pdf = pyhf.Model(spec) + data = source['bindata']['data'] + pdf.config.auxdata + + init_pars = pdf.config.suggested_init() + par_bounds = pdf.config.suggested_bounds() + + optim = pyhf.optimizer + + result = optim.minimize(pyhf.infer.mle.twice_nll, data, pdf, init_pars, par_bounds) + assert pyhf.tensorlib.tolist(result) + + result, correlations = optim.minimize( + pyhf.infer.mle.twice_nll, + data, + pdf, + init_pars, + par_bounds, + [(pdf.config.poi_index, mu)], + return_correlations=True, + ) + assert result.shape == (2,) + assert correlations.shape == (2, 2) + assert pyhf.tensorlib.tolist(result) + assert pyhf.tensorlib.tolist(correlations) + + assert np.allclose([[1.0, 0.0], [0.0, 0.0]], pyhf.tensorlib.tolist(correlations)) + + +@pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) +@pytest.mark.skip_numpy +def test_optim_correlations_autodiff(backend, source, spec, mu): + pdf = pyhf.Model(spec) data = source['bindata']['data'] + pdf.config.auxdata init_pars = pdf.config.suggested_init() From 474bd1a649a35ec4cb2a60f1d45e6bad7d31884f Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Mon, 14 Aug 2023 12:58:17 +0100 Subject: [PATCH 09/15] fix docstrings --- src/pyhf/tensor/jax_backend.py | 6 +++--- src/pyhf/tensor/numpy_backend.py | 6 +++--- src/pyhf/tensor/tensorflow_backend.py | 7 ++++--- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/pyhf/tensor/jax_backend.py b/src/pyhf/tensor/jax_backend.py index d0dd00ca64..cb1c5291f9 100644 --- a/src/pyhf/tensor/jax_backend.py +++ b/src/pyhf/tensor/jax_backend.py @@ -648,10 +648,10 @@ def diagonal(self, tensor_in): >>> pyhf.set_backend("jax") >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) >>> tensor - tensor([[1., 0.], - [0., 1.]]) + Array([[1., 0.], + [0., 1.]], dtype=float64) >>> pyhf.tensorlib.diagonal(tensor) - tensor([1., 1.]) + Array([1., 1.], dtype=float64) Args: tensor_in (:obj:`tensor`): The input tensor object. diff --git a/src/pyhf/tensor/numpy_backend.py b/src/pyhf/tensor/numpy_backend.py index fc46164764..8c540b47f3 100644 --- a/src/pyhf/tensor/numpy_backend.py +++ b/src/pyhf/tensor/numpy_backend.py @@ -668,10 +668,10 @@ def diagonal(self, tensor_in: Tensor[T]) -> ArrayLike: >>> pyhf.set_backend("numpy") >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) >>> tensor - tensor([[1., 0.], - [0., 1.]]) + array([[1., 0.], + [0., 1.]]) >>> pyhf.tensorlib.diagonal(tensor) - tensor([1., 1.]) + array([1., 1.]) Args: tensor_in (:obj:`tensor`): The input tensor object. diff --git a/src/pyhf/tensor/tensorflow_backend.py b/src/pyhf/tensor/tensorflow_backend.py index 75b8510bf2..ce3030b7f9 100644 --- a/src/pyhf/tensor/tensorflow_backend.py +++ b/src/pyhf/tensor/tensorflow_backend.py @@ -752,10 +752,11 @@ def diagonal(self, tensor_in): >>> pyhf.set_backend("tensorflow") >>> tensor = pyhf.tensorlib.astensor([[1.0, 0.0], [0.0, 1.0]]) >>> tensor - tensor([[1., 0.], - [0., 1.]]) + >>> pyhf.tensorlib.diagonal(tensor) - tensor([1., 1.]) + Date: Mon, 14 Aug 2023 13:00:41 +0100 Subject: [PATCH 10/15] adjust to ignore zeroing out corr/cov for minuit, since it already handles this --- src/pyhf/optimize/mixins.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 0504ae1d18..3d78af7cce 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -128,11 +128,8 @@ def _internal_postprocess( cov = hess_inv # we also need to edit the covariance matrix to zero-out uncertainties! - if fixed_vals is not None: - if using_minuit: - fixed_bools = fitresult.minuit.fixed - else: - fixed_bools = [False] * len(init_pars) + if fixed_vals is not None and not using_minuit: # minuit already does this + fixed_bools = [False] * len(init_pars) # Convert fixed_bools to a numpy array and reshape to make it a column vector fixed_mask = tensorlib.reshape( tensorlib.astensor(fixed_bools, dtype="bool"), (-1, 1) @@ -158,7 +155,7 @@ def _internal_postprocess( tensorlib.astensor(0.0), ) - if correlations_from_fit is not None: + if correlations_from_fit is not None and not using_minuit: _zeros = tensorlib.zeros(num_fixed_pars) # possibly a more elegant way to do this stitched_columns = [ From 5537ad91ba18fecda368868db6a02657fb21b7cf Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Mon, 14 Aug 2023 13:31:04 +0100 Subject: [PATCH 11/15] quickly patch in skipping minuit for autodiff tests (silently passing) --- tests/test_optim.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/test_optim.py b/tests/test_optim.py index 25e0016424..aea935246f 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -406,6 +406,7 @@ def test_optim_uncerts_minuit(backend, source, spec, mu): @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.skip_numpy +@pytest.mark.skip_numpy_minuit def test_optim_uncerts_autodiff(backend, source, spec, mu): pdf = pyhf.Model(spec) data = source['bindata']['data'] + pdf.config.auxdata @@ -428,8 +429,10 @@ def test_optim_uncerts_autodiff(backend, source, spec, mu): return_uncertainties=True, ) assert result.shape == (2, 2) - # TODO: add proper numerical test for autodiff uncerts (does not match minuit at all) - # assert pytest.approx([0.26418431, 0.0]) == pyhf.tensorlib.tolist(result[:, 1]) + # TODO: this does not match minuit at all (0.26418431) -- is that correct here? + assert pytest.approx([0.6693815171034548, 0.0]) == pyhf.tensorlib.tolist( + result[:, 1] + ) @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @@ -465,6 +468,7 @@ def test_optim_correlations_minuit(backend, source, spec, mu): @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.skip_numpy +@pytest.mark.skip_numpy_minuit def test_optim_correlations_autodiff(backend, source, spec, mu): pdf = pyhf.Model(spec) data = source['bindata']['data'] + pdf.config.auxdata From b1c9b3d67af1bdcc6d98103c986886bdfa3af834 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Fri, 20 Oct 2023 15:24:08 +0100 Subject: [PATCH 12/15] ignore venv in gitignore --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 66b6392026..a7893c30bc 100644 --- a/.gitignore +++ b/.gitignore @@ -44,3 +44,5 @@ htmlcov # text editors .vscode/ + +venv/ From 3a227608302ba9571924e4584c662863510878c8 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Fri, 20 Oct 2023 16:05:33 +0100 Subject: [PATCH 13/15] add poi_name='mu' --- tests/test_optim.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_optim.py b/tests/test_optim.py index aea935246f..5e0ce9338b 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -380,7 +380,7 @@ def test_optim_with_value(backend, source, spec, mu): @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.only_numpy_minuit def test_optim_uncerts_minuit(backend, source, spec, mu): - pdf = pyhf.Model(spec) + pdf = pyhf.Model(spec, poi_name="mu") data = source['bindata']['data'] + pdf.config.auxdata init_pars = pdf.config.suggested_init() @@ -408,7 +408,7 @@ def test_optim_uncerts_minuit(backend, source, spec, mu): @pytest.mark.skip_numpy @pytest.mark.skip_numpy_minuit def test_optim_uncerts_autodiff(backend, source, spec, mu): - pdf = pyhf.Model(spec) + pdf = pyhf.Model(spec, poi_name="mu") data = source['bindata']['data'] + pdf.config.auxdata init_pars = pdf.config.suggested_init() @@ -438,7 +438,7 @@ def test_optim_uncerts_autodiff(backend, source, spec, mu): @pytest.mark.parametrize('mu', [1.0], ids=['mu=1']) @pytest.mark.only_numpy_minuit def test_optim_correlations_minuit(backend, source, spec, mu): - pdf = pyhf.Model(spec) + pdf = pyhf.Model(spec, poi_name="mu") data = source['bindata']['data'] + pdf.config.auxdata init_pars = pdf.config.suggested_init() @@ -470,7 +470,7 @@ def test_optim_correlations_minuit(backend, source, spec, mu): @pytest.mark.skip_numpy @pytest.mark.skip_numpy_minuit def test_optim_correlations_autodiff(backend, source, spec, mu): - pdf = pyhf.Model(spec) + pdf = pyhf.Model(spec, poi_name="mu") data = source['bindata']['data'] + pdf.config.auxdata init_pars = pdf.config.suggested_init() From 8f286baf03e7b8dd89c5de0f8e43c6d5915b44e7 Mon Sep 17 00:00:00 2001 From: Nathan Simpson Date: Fri, 20 Oct 2023 16:33:05 +0100 Subject: [PATCH 14/15] document setting things to zero, and add pointers to the minuit issue that references this behaviour being done --- src/pyhf/optimize/mixins.py | 13 ++++++++----- tests/test_optim.py | 3 ++- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 3d78af7cce..6f7adff518 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -99,14 +99,16 @@ def _internal_postprocess( # extract number of fixed parameters num_fixed_pars = len(fitted_pars) - len(fitresult.x) - # FIXME: Set uncertainties for fixed parameters to 0 manually - # https://github.com/scikit-hep/iminuit/issues/762 - # https://github.com/scikit-hep/pyhf/issues/1918 - # https://github.com/scikit-hep/cabinetry/pull/346 + # Set uncertainties for fixed parameters to 0 manually if fixed_vals is not None: # check for fixed vals if using_minuit: + # See related discussion here: + # https://github.com/scikit-hep/iminuit/issues/762 + # https://github.com/scikit-hep/pyhf/issues/1918 + # https://github.com/scikit-hep/cabinetry/pull/346 uncertainties = np.where(fitresult.minuit.fixed, 0.0, uncertainties) else: + # Not using minuit, so don't have `fitresult.minuit.fixed` here: do it manually fixed_bools = [False] * len(init_pars) for index, _ in fixed_vals: fixed_bools[index] = True @@ -128,7 +130,8 @@ def _internal_postprocess( cov = hess_inv # we also need to edit the covariance matrix to zero-out uncertainties! - if fixed_vals is not None and not using_minuit: # minuit already does this + # NOTE: minuit already does this (https://github.com/scikit-hep/iminuit/issues/762#issuecomment-1207436406) + if fixed_vals is not None and not using_minuit: fixed_bools = [False] * len(init_pars) # Convert fixed_bools to a numpy array and reshape to make it a column vector fixed_mask = tensorlib.reshape( diff --git a/tests/test_optim.py b/tests/test_optim.py index 5e0ce9338b..36032ac090 100644 --- a/tests/test_optim.py +++ b/tests/test_optim.py @@ -429,7 +429,8 @@ def test_optim_uncerts_autodiff(backend, source, spec, mu): return_uncertainties=True, ) assert result.shape == (2, 2) - # TODO: this does not match minuit at all (0.26418431) -- is that correct here? + # NOTE: this does not match minuit at all (0.26418431), and is probably an artefact of doing a fixed POI fit on this particular model? + # see discussion in https://github.com/scikit-hep/pyhf/pull/2269 for more details. differences disappear in global fit uncertainties. assert pytest.approx([0.6693815171034548, 0.0]) == pyhf.tensorlib.tolist( result[:, 1] ) From d3142fc71f3266e7a2318765ec34bd27868b2bfc Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 20 Oct 2023 15:34:57 +0000 Subject: [PATCH 15/15] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/pyhf/optimize/mixins.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pyhf/optimize/mixins.py b/src/pyhf/optimize/mixins.py index 6f7adff518..ca51fec978 100644 --- a/src/pyhf/optimize/mixins.py +++ b/src/pyhf/optimize/mixins.py @@ -131,7 +131,7 @@ def _internal_postprocess( # we also need to edit the covariance matrix to zero-out uncertainties! # NOTE: minuit already does this (https://github.com/scikit-hep/iminuit/issues/762#issuecomment-1207436406) - if fixed_vals is not None and not using_minuit: + if fixed_vals is not None and not using_minuit: fixed_bools = [False] * len(init_pars) # Convert fixed_bools to a numpy array and reshape to make it a column vector fixed_mask = tensorlib.reshape(