gammapy.modeling.fit.FitStepResult.method() - Code Metrics - Inspection of "Prepare changelog for v0.19" - gammapy/gammapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — master (#3318)

unknown

created 2021-11-22 15:35 UTC

gammapy.modeling.fit.FitStepResult.method() A

↳ Parent: gammapy.modeling.fit

Complexity

Conditions

Size

Total Lines	4
Code Lines	3

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	1
eloc	3
nop	1
dl	0
loc	4
rs	10
c	0
b	0
f	0

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import itertools
import logging
import numpy as np
from gammapy.utils.pbar import progress_bar
from gammapy.utils.table import table_from_row_data
from .covariance import Covariance
from .iminuit import (
    confidence_iminuit,
    contour_iminuit,
    covariance_iminuit,
    optimize_iminuit,
)
from .scipy import confidence_scipy, optimize_scipy
from .sherpa import optimize_sherpa

__all__ = ["Fit"]

log = logging.getLogger(__name__)


class Registry:
    """Registry of available backends for given tasks.

    Gives users the power to extend from their scripts.
    Used by `Fit` below.

    Not sure if we should call it "backend" or "method" or something else.
    Probably we will code up some methods, e.g. for profile analysis ourselves,
    using scipy or even just Python / Numpy?
    """

    register = {
        "optimize": {
            "minuit": optimize_iminuit,
            "sherpa": optimize_sherpa,
            "scipy": optimize_scipy,
        },
        "covariance": {
            "minuit": covariance_iminuit,
            # "sherpa": covariance_sherpa,
            # "scipy": covariance_scipy,
        },
        "confidence": {
            "minuit": confidence_iminuit,
            # "sherpa": confidence_sherpa,
            "scipy": confidence_scipy,
        },
    }

    @classmethod
    def get(cls, task, backend):
        if task not in cls.register:
            raise ValueError(f"Unknown task {task!r}")

        backend_options = cls.register[task]

        if backend not in backend_options:
            raise ValueError(f"Unknown backend {backend!r} for task {task!r}")

        return backend_options[backend]


registry = Registry()


class Fit:
    """Fit class.

    The fit class provides a uniform interface to multiple fitting backends.
    Currently available: "minuit", "sherpa" and "scipy"

    Parameters
    ----------
    backend : {"minuit", "scipy" "sherpa"}
        Global backend used for fitting, default : minuit
    optimize_opts : dict
        Keyword arguments passed to the optimizer. For the `"minuit"` backend
        see https://iminuit.readthedocs.io/en/latest/api.html#iminuit.Minuit
        for a detailed description of the available options. If there is an entry
        'migrad_opts', those options will be passed to `iminuit.Minuit.migrad()`.

        For the `"sherpa"` backend you can from the options `method = {"simplex",  "levmar", "moncar", "gridsearch"}`
        Those methods are described and compared in detail on
        http://cxc.cfa.harvard.edu/sherpa/methods/index.html. The available
        options of the optimization methods are described on the following
        pages in detail:

            * http://cxc.cfa.harvard.edu/sherpa/ahelp/neldermead.html
            * http://cxc.cfa.harvard.edu/sherpa/ahelp/montecarlo.html
            * http://cxc.cfa.harvard.edu/sherpa/ahelp/gridsearch.html
            * http://cxc.cfa.harvard.edu/sherpa/ahelp/levmar.html

        For the `"scipy"` backend the available options are described in detail here:
        https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html

    covariance_opts : dict
        Covariance options passed to the given backend.
    confidence_opts : dict
        Extra arguments passed to the backend. E.g. `iminuit.Minuit.minos` supports
        a ``maxcall`` option. For the scipy backend ``confidence_opts`` are forwarded
        to `~scipy.optimize.brentq`. If the confidence estimation fails, the bracketing
        interval can be adapted by modifying the the upper bound of the interval (``b``) value.
    store_trace : bool
        Whether to store the trace of the fit
    """

    def __init__(
        self,
        backend="minuit",
        optimize_opts=None,
        covariance_opts=None,
        confidence_opts=None,
        store_trace=False,
    ):
        self.store_trace = store_trace
        self.backend = backend

        if optimize_opts is None:
            optimize_opts = {"backend": backend}

        if covariance_opts is None:
            covariance_opts = {"backend": backend}

        if confidence_opts is None:
            confidence_opts = {"backend": backend}

        self.optimize_opts = optimize_opts
        self.covariance_opts = covariance_opts
        self.confidence_opts = confidence_opts
        self._minuit = None

    @property
    def minuit(self):
        """Iminuit object"""
        return self._minuit

    @staticmethod
    def _parse_datasets(datasets):
        from gammapy.datasets import Datasets

        datasets = Datasets(datasets)
        return datasets, datasets.parameters

    def run(self, datasets):
        """Run all fitting steps.

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.

        Returns
        -------
        fit_result : `FitResult`
            Fit result
        """
        optimize_result = self.optimize(datasets=datasets)

        if self.backend not in registry.register["covariance"]:
            log.warning("No covariance estimate - not supported by this backend.")
            return optimize_result

        covariance_result = self.covariance(datasets=datasets)

        return FitResult(
            optimize_result=optimize_result,
            covariance_result=covariance_result,
        )

    def optimize(self, datasets):
        """Run the optimization.

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.

        Returns
        -------
        optimize_result : `OptimizeResult`
            Optimization result
        """
        datasets, parameters = self._parse_datasets(datasets=datasets)
        datasets.parameters.check_limits()

        parameters.autoscale()

        kwargs = self.optimize_opts.copy()
        backend = kwargs.pop("backend", self.backend)

        compute = registry.get("optimize", backend)
        # TODO: change this calling interface!
        # probably should pass a fit statistic, which has a model, which has parameters
        # and return something simpler, not a tuple of three things
        factors, info, optimizer = compute(
            parameters=parameters,
            function=datasets.stat_sum,
            store_trace=self.store_trace,
            **kwargs,
        )

        if backend == "minuit":
            self._minuit = optimizer
            kwargs["method"] = "migrad"

        trace = table_from_row_data(info.pop("trace"))

        if self.store_trace:
            idx = [
                parameters.index(par)
                for par in parameters.unique_parameters.free_parameters
            ]
            unique_names = np.array(datasets.models.parameters_unique_names)[idx]
            trace.rename_columns(trace.colnames[1:], list(unique_names))

        # Copy final results into the parameters object
        parameters.set_parameter_factors(factors)
        parameters.check_limits()
        return OptimizeResult(
            parameters=parameters,
            total_stat=datasets.stat_sum(),
            backend=backend,
            method=kwargs.get("method", backend),
            trace=trace,
            **info,
        )

    def covariance(self, datasets):
        """Estimate the covariance matrix.

        Assumes that the model parameters are already optimised.

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.

        Returns
        -------
        result : `CovarianceResult`
            Results
        """
        datasets, parameters = self._parse_datasets(datasets=datasets)

        kwargs = self.covariance_opts.copy()
        kwargs["minuit"] = self.minuit
        backend = kwargs.pop("backend", self.backend)
        compute = registry.get("covariance", backend)

        with parameters.restore_status():
            if self.backend == "minuit":
                method = "hesse"
            else:
                method = ""

            factor_matrix, info = compute(
                parameters=parameters, function=datasets.stat_sum, **kwargs
            )

            datasets.models.covariance = Covariance.from_factor_matrix(
                parameters=parameters, matrix=factor_matrix
            )

        # TODO: decide what to return, and fill the info correctly!
        return CovarianceResult(
            backend=backend,
            method=method,
            success=info["success"],
            message=info["message"],
        )

    def confidence(self, datasets, parameter, sigma=1, reoptimize=True):
        """Estimate confidence interval.

        Extra ``kwargs`` are passed to the backend.
        E.g. `iminuit.Minuit.minos` supports a ``maxcall`` option.

        For the scipy backend ``kwargs`` are forwarded to `~scipy.optimize.brentq`. If the
        confidence estimation fails, the bracketing interval can be adapted by modifying the
        the upper bound of the interval (``b``) value.

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.
        parameter : `~gammapy.modeling.Parameter`
            Parameter of interest
        sigma : float
            Number of standard deviations for the confidence level
        reoptimize : bool
            Re-optimize other parameters, when computing the confidence region.

        Returns
        -------
        result : dict
            Dictionary with keys "errp", 'errn", "success" and "nfev".
        """
        datasets, parameters = self._parse_datasets(datasets=datasets)

        kwargs = self.confidence_opts.copy()
        backend = kwargs.pop("backend", self.backend)

        compute = registry.get("confidence", backend)
        parameter = parameters[parameter]

        with parameters.restore_status():
            result = compute(
                parameters=parameters,
                parameter=parameter,
                function=datasets.stat_sum,
                sigma=sigma,
                reoptimize=reoptimize,
                **kwargs,
            )

        result["errp"] *= parameter.scale
        result["errn"] *= parameter.scale
        return result

    def stat_profile(self, datasets, parameter, reoptimize=False):
        """Compute fit statistic profile.

        The method used is to vary one parameter, keeping all others fixed.
        So this is taking a "slice" or "scan" of the fit statistic.

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.
        parameter : `~gammapy.modeling.Parameter`
            Parameter of interest. The specification for the scan, such as bounds
            and number of values is taken from the parameter object.
        reoptimize : bool
            Re-optimize other parameters, when computing the confidence region.

        Returns
        -------
        results : dict
            Dictionary with keys "values", "stat" and "fit_results". The latter contains an
            empty list, if `reoptimize` is set to False
        """
        datasets, parameters = self._parse_datasets(datasets=datasets)
        parameter = parameters[parameter]
        values = parameter.scan_values

        stats = []
        fit_results = []
        with parameters.restore_status():
            for value in progress_bar(values, desc="Scan values"):
                parameter.value = value
                if reoptimize:
                    parameter.frozen = True
                    result = self.optimize(datasets=datasets)
                    stat = result.total_stat
                    fit_results.append(result)
                else:
                    stat = datasets.stat_sum()
                stats.append(stat)

        return {
            f"{parameter.name}_scan": values,
            "stat_scan": np.array(stats),
            "fit_results": fit_results,
        }

    def stat_surface(self, datasets, x, y, reoptimize=False):
        """Compute fit statistic surface.

        The method used is to vary two parameters, keeping all others fixed.
        So this is taking a "slice" or "scan" of the fit statistic.

        Caveat: This method can be very computationally intensive and slow

        See also: `Fit.stat_contour`

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.
        x, y : `~gammapy.modeling.Parameter`
            Parameters of interest
        reoptimize : bool
            Re-optimize other parameters, when computing the confidence region.

        Returns
        -------
        results : dict
            Dictionary with keys "x_values", "y_values", "stat" and "fit_results". The latter contains an
            empty list, if `reoptimize` is set to False
        """
        datasets, parameters = self._parse_datasets(datasets=datasets)

        x, y = parameters[x], parameters[y]

        stats = []
        fit_results = []

        with parameters.restore_status():
            for x_value, y_value in progress_bar(
                itertools.product(x.scan_values, y.scan_values), desc="Trial values"
            ):
                x.value, y.value = x_value, y_value

                if reoptimize:
                    x.frozen, y.frozen = True, True
                    result = self.optimize(datasets=datasets)
                    stat = result.total_stat
                    fit_results.append(result)
                else:
                    stat = datasets.stat_sum()

                stats.append(stat)

        shape = (len(x.scan_values), len(y.scan_values))
        stats = np.array(stats).reshape(shape)

        if reoptimize:
            fit_results = np.array(fit_results).reshape(shape)

        return {
            f"{x.name}_scan": x.scan_values,
            f"{y.name}_scan": y.scan_values,
            "stat_scan": stats,
            "fit_results": fit_results,
        }

    def stat_contour(self, datasets, x, y, numpoints=10, sigma=1):
        """Compute stat contour.

        Calls ``iminuit.Minuit.mncontour``.

        This is a contouring algorithm for a 2D function
        which is not simply the fit statistic function.
        That 2D function is given at each point ``(par_1, par_2)``
        by re-optimising all other free parameters,
        and taking the fit statistic at that point.

        Very compute-intensive and slow.

        Parameters
        ----------
        datasets : `Datasets` or list of `Dataset`
            Datasets to optimize.
        x, y : `~gammapy.modeling.Parameter`
            Parameters of interest
        numpoints : int
            Number of contour points
        sigma : float
            Number of standard deviations for the confidence level

        Returns
        -------
        result : dict
            Dictionary containing the parameter values defining the contour, with the
            boolean flag "success" and the info objects from ``mncontour``.
        """
        datasets, parameters = self._parse_datasets(datasets=datasets)

        x = parameters[x]
        y = parameters[y]

        with parameters.restore_status():
            result = contour_iminuit(
                parameters=parameters,
                function=datasets.stat_sum,
                x=x,
                y=y,
                numpoints=numpoints,
                sigma=sigma,
            )

        x_name = x.name
        y_name = y.name
        x = result["x"] * x.scale
        y = result["y"] * y.scale

        return {
            x_name: x,
            y_name: y,
            "success": result["success"],
        }


class FitStepResult:
    """Fit result base class"""

    def __init__(self, backend, method, success, message):
        self._success = success
        self._message = message
        self._backend = backend
        self._method = method

    @property
    def backend(self):
        """Optimizer backend used for the fit."""
        return self._backend

    @property
    def method(self):
        """Optimizer method used for the fit."""
        return self._method

    @property
    def success(self):
        """Fit success status flag."""
        return self._success

    @property
    def message(self):
        """Optimizer status message."""
        return self._message

    def __repr__(self):
        return (
            f"{self.__class__.__name__}\n\n"
            f"\tbackend    : {self.backend}\n"
            f"\tmethod     : {self.method}\n"
            f"\tsuccess    : {self.success}\n"
            f"\tmessage    : {self.message}\n"
        )


class CovarianceResult(FitStepResult):
    """Covariance result object."""

    pass


class OptimizeResult(FitStepResult):
    """Optimize result object."""
    def __init__(self, parameters, nfev, total_stat, trace, **kwargs):
        self._parameters = parameters
        self._nfev = nfev
        self._total_stat = total_stat
        self._trace = trace
        super().__init__(**kwargs)

    @property
    def parameters(self):
        """Best fit parameters"""
        return self._parameters

    @property
    def trace(self):
        """Parameter trace from the optimisation"""
        return self._trace

    @property
    def nfev(self):
        """Number of function evaluations."""
        return self._nfev

    @property
    def total_stat(self):
        """Value of the fit statistic at minimum."""
        return self._total_stat

    def __repr__(self):
        str_ = super().__repr__()
        str_ += f"\tnfev       : {self.nfev}\n"
        str_ += f"\ttotal stat : {self.total_stat:.2f}\n\n"
        return str_


class FitResult:
    """Fit result class

    Parameters
    ----------
    optimize_result : `OptimizeResult`
        Result of the optimization step.
    covariance_result : `CovarianceResult`
        Result of the covariance step.
    """
    def __init__(self, optimize_result=None, covariance_result=None):
        self._optimize_result = optimize_result
        self._covariance_result = covariance_result

    # TODO: is the convenience access needed?
    @property
    def parameters(self):
        """Best fit parameters of the optimization step"""
        return self.optimize_result.parameters

    # TODO: is the convenience access needed?
    @property
    def total_stat(self):
        """Total stat of the optimization step"""
        return self.optimize_result.total_stat

    # TODO: is the convenience access needed?
    @property
    def trace(self):
        """Parameter trace of the optimisation step"""
        return self.optimize_result.trace

    # TODO: is the convenience access needed?
    @property
    def nfev(self):
        """Number of function evaluations of the optimisation step"""
        return self.optimize_result.nfev

    # TODO: is the convenience access needed?
    @property
    def backend(self):
        """Optimizer backend used for the fit."""
        return self.optimize_result.backend

    # TODO: is the convenience access needed?
    @property
    def method(self):
        """Optimizer method used for the fit."""
        return self.optimize_result.method

    # TODO: is the convenience access needed?
    @property
    def message(self):
        """Optimizer status message."""
        return self.optimize_result.message

    @property
    def success(self):
        """Total success flag"""
        success = self.optimize_result.success and self.covariance_result.success
        return success

    @property
    def optimize_result(self):
        """Optimize result"""
        return self._optimize_result

    @property
    def covariance_result(self):
        """Optimize result"""
        return self._optimize_result

    def __repr__(self):
        str_ = ""
        if self.optimize_result:
            str_ += str(self.optimize_result)

        if self.covariance_result:
            str_ += str(self.covariance_result)

        return str_



1			# Licensed under a 3-clause BSD style license - see LICENSE.rst
2			import itertools
3			import logging
4			import numpy as np
5			from gammapy.utils.pbar import progress_bar
6			from gammapy.utils.table import table_from_row_data
7			from .covariance import Covariance
8			from .iminuit import (
9			confidence_iminuit,
10			contour_iminuit,
11			covariance_iminuit,
12			optimize_iminuit,
13			)
14			from .scipy import confidence_scipy, optimize_scipy
15			from .sherpa import optimize_sherpa
16
17			__all__ = ["Fit"]
18
19			log = logging.getLogger(__name__)
20
21
22			class Registry:
23			"""Registry of available backends for given tasks.
24
25			Gives users the power to extend from their scripts.
26			Used by `Fit` below.
27
28			Not sure if we should call it "backend" or "method" or something else.
29			Probably we will code up some methods, e.g. for profile analysis ourselves,
30			using scipy or even just Python / Numpy?
31			"""
32
33			register = {
34			"optimize": {
35			"minuit": optimize_iminuit,
36			"sherpa": optimize_sherpa,
37			"scipy": optimize_scipy,
38			},
39			"covariance": {
40			"minuit": covariance_iminuit,
41			# "sherpa": covariance_sherpa,
42			# "scipy": covariance_scipy,
43			},
44			"confidence": {
45			"minuit": confidence_iminuit,
46			# "sherpa": confidence_sherpa,
47			"scipy": confidence_scipy,
48			},
49			}
50
51			@classmethod
52			def get(cls, task, backend):
53			if task not in cls.register:
54			raise ValueError(f"Unknown task {task!r}")
55
56			backend_options = cls.register[task]
57
58			if backend not in backend_options:
59			raise ValueError(f"Unknown backend {backend!r} for task {task!r}")
60
61			return backend_options[backend]
62
63
64			registry = Registry()
65
66
67			class Fit:
68			"""Fit class.
69
70			The fit class provides a uniform interface to multiple fitting backends.
71			Currently available: "minuit", "sherpa" and "scipy"
72
73			Parameters
74			----------
75			backend : {"minuit", "scipy" "sherpa"}
76			Global backend used for fitting, default : minuit
77			optimize_opts : dict
78			Keyword arguments passed to the optimizer. For the `"minuit"` backend
79			see https://iminuit.readthedocs.io/en/latest/api.html#iminuit.Minuit
80			for a detailed description of the available options. If there is an entry
81			'migrad_opts', those options will be passed to `iminuit.Minuit.migrad()`.
82
83			For the `"sherpa"` backend you can from the options `method = {"simplex", "levmar", "moncar", "gridsearch"}`
84			Those methods are described and compared in detail on
85			http://cxc.cfa.harvard.edu/sherpa/methods/index.html. The available
86			options of the optimization methods are described on the following
87			pages in detail:
88
89			* http://cxc.cfa.harvard.edu/sherpa/ahelp/neldermead.html
90			* http://cxc.cfa.harvard.edu/sherpa/ahelp/montecarlo.html
91			* http://cxc.cfa.harvard.edu/sherpa/ahelp/gridsearch.html
92			* http://cxc.cfa.harvard.edu/sherpa/ahelp/levmar.html
93
94			For the `"scipy"` backend the available options are described in detail here:
95			https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
96
97			covariance_opts : dict
98			Covariance options passed to the given backend.
99			confidence_opts : dict
100			Extra arguments passed to the backend. E.g. `iminuit.Minuit.minos` supports
101			a ``maxcall`` option. For the scipy backend ``confidence_opts`` are forwarded
102			to `~scipy.optimize.brentq`. If the confidence estimation fails, the bracketing
103			interval can be adapted by modifying the the upper bound of the interval (``b``) value.
104			store_trace : bool
105			Whether to store the trace of the fit
106			"""
107
108			def __init__(
109			self,
110			backend="minuit",
111			optimize_opts=None,
112			covariance_opts=None,
113			confidence_opts=None,
114			store_trace=False,
115			):
116			self.store_trace = store_trace
117			self.backend = backend
118
119			if optimize_opts is None:
120			optimize_opts = {"backend": backend}
121
122			if covariance_opts is None:
123			covariance_opts = {"backend": backend}
124
125			if confidence_opts is None:
126			confidence_opts = {"backend": backend}
127
128			self.optimize_opts = optimize_opts
129			self.covariance_opts = covariance_opts
130			self.confidence_opts = confidence_opts
131			self._minuit = None
132
133			@property
134			def minuit(self):
135			"""Iminuit object"""
136			return self._minuit
137
138			@staticmethod
139			def _parse_datasets(datasets):
140			from gammapy.datasets import Datasets
141
142			datasets = Datasets(datasets)
143			return datasets, datasets.parameters
144
145			def run(self, datasets):
146			"""Run all fitting steps.
147
148			Parameters
149			----------
150			datasets : `Datasets` or list of `Dataset`
151			Datasets to optimize.
152
153			Returns
154			-------
155			fit_result : `FitResult`
156			Fit result
157			"""
158			optimize_result = self.optimize(datasets=datasets)
159
160			if self.backend not in registry.register["covariance"]:
161			log.warning("No covariance estimate - not supported by this backend.")
162			return optimize_result
163
164			covariance_result = self.covariance(datasets=datasets)
165
166			return FitResult(
167			optimize_result=optimize_result,
168			covariance_result=covariance_result,
169			)
170
171			def optimize(self, datasets):
172			"""Run the optimization.
173
174			Parameters
175			----------
176			datasets : `Datasets` or list of `Dataset`
177			Datasets to optimize.
178
179			Returns
180			-------
181			optimize_result : `OptimizeResult`
182			Optimization result
183			"""
184			datasets, parameters = self._parse_datasets(datasets=datasets)
185			datasets.parameters.check_limits()
186
187			parameters.autoscale()
188
189			kwargs = self.optimize_opts.copy()
190			backend = kwargs.pop("backend", self.backend)
191
192			compute = registry.get("optimize", backend)
193			# TODO: change this calling interface!
194			# probably should pass a fit statistic, which has a model, which has parameters
195			# and return something simpler, not a tuple of three things
196			factors, info, optimizer = compute(
197			parameters=parameters,
198			function=datasets.stat_sum,
199			store_trace=self.store_trace,
200			**kwargs,
201			)
202
203			if backend == "minuit":
204			self._minuit = optimizer
205			kwargs["method"] = "migrad"
206
207			trace = table_from_row_data(info.pop("trace"))
208
209			if self.store_trace:
210			idx = [
211			parameters.index(par)
212			for par in parameters.unique_parameters.free_parameters
213			]
214			unique_names = np.array(datasets.models.parameters_unique_names)[idx]
215			trace.rename_columns(trace.colnames[1:], list(unique_names))
216
217			# Copy final results into the parameters object
218			parameters.set_parameter_factors(factors)
219			parameters.check_limits()
220			return OptimizeResult(
221			parameters=parameters,
222			total_stat=datasets.stat_sum(),
223			backend=backend,
224			method=kwargs.get("method", backend),
225			trace=trace,
226			**info,
227			)
228
229			def covariance(self, datasets):
230			"""Estimate the covariance matrix.
231
232			Assumes that the model parameters are already optimised.
233
234			Parameters
235			----------
236			datasets : `Datasets` or list of `Dataset`
237			Datasets to optimize.
238
239			Returns
240			-------
241			result : `CovarianceResult`
242			Results
243			"""
244			datasets, parameters = self._parse_datasets(datasets=datasets)
245
246			kwargs = self.covariance_opts.copy()
247			kwargs["minuit"] = self.minuit
248			backend = kwargs.pop("backend", self.backend)
249			compute = registry.get("covariance", backend)
250
251			with parameters.restore_status():
252			if self.backend == "minuit":
253			method = "hesse"
254			else:
255			method = ""
256
257			factor_matrix, info = compute(
258			parameters=parameters, function=datasets.stat_sum, **kwargs
259			)
260
261			datasets.models.covariance = Covariance.from_factor_matrix(
262			parameters=parameters, matrix=factor_matrix
263			)
264
265			# TODO: decide what to return, and fill the info correctly!
266			return CovarianceResult(
267			backend=backend,
268			method=method,
269			success=info["success"],
270			message=info["message"],
271			)
272
273			def confidence(self, datasets, parameter, sigma=1, reoptimize=True):
274			"""Estimate confidence interval.
275
276			Extra ``kwargs`` are passed to the backend.
277			E.g. `iminuit.Minuit.minos` supports a ``maxcall`` option.
278
279			For the scipy backend ``kwargs`` are forwarded to `~scipy.optimize.brentq`. If the
280			confidence estimation fails, the bracketing interval can be adapted by modifying the
281			the upper bound of the interval (``b``) value.
282
283			Parameters
284			----------
285			datasets : `Datasets` or list of `Dataset`
286			Datasets to optimize.
287			parameter : `~gammapy.modeling.Parameter`
288			Parameter of interest
289			sigma : float
290			Number of standard deviations for the confidence level
291			reoptimize : bool
292			Re-optimize other parameters, when computing the confidence region.
293
294			Returns
295			-------
296			result : dict
297			Dictionary with keys "errp", 'errn", "success" and "nfev".
298			"""
299			datasets, parameters = self._parse_datasets(datasets=datasets)
300
301			kwargs = self.confidence_opts.copy()
302			backend = kwargs.pop("backend", self.backend)
303
304			compute = registry.get("confidence", backend)
305			parameter = parameters[parameter]
306
307			with parameters.restore_status():
308			result = compute(
309			parameters=parameters,
310			parameter=parameter,
311			function=datasets.stat_sum,
312			sigma=sigma,
313			reoptimize=reoptimize,
314			**kwargs,
315			)
316
317			result["errp"] *= parameter.scale
318			result["errn"] *= parameter.scale
319			return result
320
321			def stat_profile(self, datasets, parameter, reoptimize=False):
322			"""Compute fit statistic profile.
323
324			The method used is to vary one parameter, keeping all others fixed.
325			So this is taking a "slice" or "scan" of the fit statistic.
326
327			Parameters
328			----------
329			datasets : `Datasets` or list of `Dataset`
330			Datasets to optimize.
331			parameter : `~gammapy.modeling.Parameter`
332			Parameter of interest. The specification for the scan, such as bounds
333			and number of values is taken from the parameter object.
334			reoptimize : bool
335			Re-optimize other parameters, when computing the confidence region.
336
337			Returns
338			-------
339			results : dict
340			Dictionary with keys "values", "stat" and "fit_results". The latter contains an
341			empty list, if `reoptimize` is set to False
342			"""
343			datasets, parameters = self._parse_datasets(datasets=datasets)
344			parameter = parameters[parameter]
345			values = parameter.scan_values
346
347			stats = []
348			fit_results = []
349			with parameters.restore_status():
350			for value in progress_bar(values, desc="Scan values"):
351			parameter.value = value
352			if reoptimize:
353			parameter.frozen = True
354			result = self.optimize(datasets=datasets)
355			stat = result.total_stat
356			fit_results.append(result)
357			else:
358			stat = datasets.stat_sum()
359			stats.append(stat)
360
361			return {
362			f"{parameter.name}_scan": values,
363			"stat_scan": np.array(stats),
364			"fit_results": fit_results,
365			}
366
367			def stat_surface(self, datasets, x, y, reoptimize=False):
368			"""Compute fit statistic surface.
369
370			The method used is to vary two parameters, keeping all others fixed.
371			So this is taking a "slice" or "scan" of the fit statistic.
372
373			Caveat: This method can be very computationally intensive and slow
374
375			See also: `Fit.stat_contour`
376
377			Parameters
378			----------
379			datasets : `Datasets` or list of `Dataset`
380			Datasets to optimize.
381			x, y : `~gammapy.modeling.Parameter`
382			Parameters of interest
383			reoptimize : bool
384			Re-optimize other parameters, when computing the confidence region.
385
386			Returns
387			-------
388			results : dict
389			Dictionary with keys "x_values", "y_values", "stat" and "fit_results". The latter contains an
390			empty list, if `reoptimize` is set to False
391			"""
392			datasets, parameters = self._parse_datasets(datasets=datasets)
393
394			x, y = parameters[x], parameters[y]
395
396			stats = []
397			fit_results = []
398
399			with parameters.restore_status():
400			for x_value, y_value in progress_bar(
401			itertools.product(x.scan_values, y.scan_values), desc="Trial values"
402			):
403			x.value, y.value = x_value, y_value
404
405			if reoptimize:
406			x.frozen, y.frozen = True, True
407			result = self.optimize(datasets=datasets)
408			stat = result.total_stat
409			fit_results.append(result)
410			else:
411			stat = datasets.stat_sum()
412
413			stats.append(stat)
414
415			shape = (len(x.scan_values), len(y.scan_values))
416			stats = np.array(stats).reshape(shape)
417
418			if reoptimize:
419			fit_results = np.array(fit_results).reshape(shape)
420
421			return {
422			f"{x.name}_scan": x.scan_values,
423			f"{y.name}_scan": y.scan_values,
424			"stat_scan": stats,
425			"fit_results": fit_results,
426			}
427
428			def stat_contour(self, datasets, x, y, numpoints=10, sigma=1):
429			"""Compute stat contour.
430
431			Calls ``iminuit.Minuit.mncontour``.
432
433			This is a contouring algorithm for a 2D function
434			which is not simply the fit statistic function.
435			That 2D function is given at each point ``(par_1, par_2)``
436			by re-optimising all other free parameters,
437			and taking the fit statistic at that point.
438
439			Very compute-intensive and slow.
440
441			Parameters
442			----------
443			datasets : `Datasets` or list of `Dataset`
444			Datasets to optimize.
445			x, y : `~gammapy.modeling.Parameter`
446			Parameters of interest
447			numpoints : int
448			Number of contour points
449			sigma : float
450			Number of standard deviations for the confidence level
451
452			Returns
453			-------
454			result : dict
455			Dictionary containing the parameter values defining the contour, with the
456			boolean flag "success" and the info objects from ``mncontour``.
457			"""
458			datasets, parameters = self._parse_datasets(datasets=datasets)
459
460			x = parameters[x]
461			y = parameters[y]
462
463			with parameters.restore_status():
464			result = contour_iminuit(
465			parameters=parameters,
466			function=datasets.stat_sum,
467			x=x,
468			y=y,
469			numpoints=numpoints,
470			sigma=sigma,
471			)
472
473			x_name = x.name
474			y_name = y.name
475			x = result["x"] * x.scale
476			y = result["y"] * y.scale
477
478			return {
479			x_name: x,
480			y_name: y,
481			"success": result["success"],
482			}
483
484
485			class FitStepResult:
486			"""Fit result base class"""
487
488			def __init__(self, backend, method, success, message):
489			self._success = success
490			self._message = message
491			self._backend = backend
492			self._method = method
493
494			@property
495			def backend(self):
496			"""Optimizer backend used for the fit."""
497			return self._backend
498
499			@property
500			def method(self):
501			"""Optimizer method used for the fit."""
502			return self._method
503
504			@property
505			def success(self):
506			"""Fit success status flag."""
507			return self._success
508
509			@property
510			def message(self):
511			"""Optimizer status message."""
512			return self._message
513
514			def __repr__(self):
515			return (
516			f"{self.__class__.__name__}\n\n"
517			f"\tbackend : {self.backend}\n"
518			f"\tmethod : {self.method}\n"
519			f"\tsuccess : {self.success}\n"
520			f"\tmessage : {self.message}\n"
521			)
522
523
524			class CovarianceResult(FitStepResult):
525			"""Covariance result object."""
526
527			pass
528
529
530			class OptimizeResult(FitStepResult):
531			"""Optimize result object."""
532			def __init__(self, parameters, nfev, total_stat, trace, **kwargs):
533			self._parameters = parameters
534			self._nfev = nfev
535			self._total_stat = total_stat
536			self._trace = trace
537			super().__init__(**kwargs)
538
539			@property
540			def parameters(self):
541			"""Best fit parameters"""
542			return self._parameters
543
544			@property
545			def trace(self):
546			"""Parameter trace from the optimisation"""
547			return self._trace
548
549			@property
550			def nfev(self):
551			"""Number of function evaluations."""
552			return self._nfev
553
554			@property
555			def total_stat(self):
556			"""Value of the fit statistic at minimum."""
557			return self._total_stat
558
559			def __repr__(self):
560			str_ = super().__repr__()
561			str_ += f"\tnfev : {self.nfev}\n"
562			str_ += f"\ttotal stat : {self.total_stat:.2f}\n\n"
563			return str_
564
565
566			class FitResult:
567			"""Fit result class
568
569			Parameters
570			----------
571			optimize_result : `OptimizeResult`
572			Result of the optimization step.
573			covariance_result : `CovarianceResult`
574			Result of the covariance step.
575			"""
576			def __init__(self, optimize_result=None, covariance_result=None):
577			self._optimize_result = optimize_result
578			self._covariance_result = covariance_result
579
580			# TODO: is the convenience access needed?
581			@property
582			def parameters(self):
583			"""Best fit parameters of the optimization step"""
584			return self.optimize_result.parameters
585
586			# TODO: is the convenience access needed?
587			@property
588			def total_stat(self):
589			"""Total stat of the optimization step"""
590			return self.optimize_result.total_stat
591
592			# TODO: is the convenience access needed?
593			@property
594			def trace(self):
595			"""Parameter trace of the optimisation step"""
596			return self.optimize_result.trace
597
598			# TODO: is the convenience access needed?
599			@property
600			def nfev(self):
601			"""Number of function evaluations of the optimisation step"""
602			return self.optimize_result.nfev
603
604			# TODO: is the convenience access needed?
605			@property
606			def backend(self):
607			"""Optimizer backend used for the fit."""
608			return self.optimize_result.backend
609
610			# TODO: is the convenience access needed?
611			@property
612			def method(self):
613			"""Optimizer method used for the fit."""
614			return self.optimize_result.method
615
616			# TODO: is the convenience access needed?
617			@property
618			def message(self):
619			"""Optimizer status message."""
620			return self.optimize_result.message
621
622			@property
623			def success(self):
624			"""Total success flag"""
625			success = self.optimize_result.success and self.covariance_result.success
626			return success
627
628			@property
629			def optimize_result(self):
630			"""Optimize result"""
631			return self._optimize_result
632
633			@property
634			def covariance_result(self):
635			"""Optimize result"""
636			return self._optimize_result
637
638			def __repr__(self):
639			str_ = ""
640			if self.optimize_result:
641			str_ += str(self.optimize_result)
642
643			if self.covariance_result:
644			str_ += str(self.covariance_result)
645
646			return str_
647
648

gammapy / gammapy

Pull Request — master (#3318)

gammapy.modeling.fit.FitStepResult.method() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like