Source code for threeML.minimizer.minuit_minimizer

from __future__ import print_function

import collections
from builtins import range

import numpy as np
from iminuit import Minuit

from import setup_logger
from threeML.minimizer.minimization import (CannotComputeCovariance,
                                            CannotComputeErrors, FitFailed,

log = setup_logger(__name__)

[docs]class MINOSFailed(Exception): pass
# This is a function to add a method to a class # We will need it in the MinuitMinimizer
[docs]def add_method(self, method, name=None): if name is None: name = method.__name__ setattr(self.__class__, name, method)
[docs]class MinuitMinimizer(LocalMinimizer): valid_setup_keys = ("ftol",) # @TODO: Is this still relevant? # NOTE: this class is built to be able to work both with iMinuit and with a boost interface to SEAL # minuit, i.e., it does not rely on functionality that iMinuit provides which is not of the original # minuit. This makes the implementation a little bit more cumbersome, but more adaptable if we want # to switch back to the bare bone SEAL minuit def __init__(self, function, parameters, verbosity=0, setup_dict=None): # This will contain the results of the last call to Migrad self._last_migrad_results = None super(MinuitMinimizer, self).__init__( function, parameters, verbosity, setup_dict ) def _setup(self, user_setup_dict): # Prepare the dictionaries for the parameters which will be used by iminuit iminuit_init_parameters = collections.OrderedDict() iminuit_errors = collections.OrderedDict() iminuit_limits = collections.OrderedDict() iminuit_fixed_parameters = collections.OrderedDict() # List of variable names that will be used for iminuit variable_names_for_iminuit = [] # NOTE: we use the internal_ versions of value, min_value and max_value because they don't have # units, and they are transformed to make the fit easier (for example in log scale) for ( parameter_path, (value, delta, minimum, maximum), ) in self._internal_parameters.items(): current_name = self._parameter_name_to_minuit_name(parameter_path) variable_names_for_iminuit.append(current_name) # Initial value iminuit_init_parameters["%s" % current_name] = value # Initial delta iminuit_errors["%s" % current_name] = delta # Limits iminuit_limits["%s" % current_name] = (minimum, maximum) # This is useless, since all parameters here are free, # but do it anyway for clarity iminuit_fixed_parameters["%s" % current_name] = False # Tell imnuit what parameter names are iminuit_init_parameters["name"] = variable_names_for_iminuit # Finally we can instance the Minuit class self.minuit = Minuit(self.function, **iminuit_init_parameters) for param, value in iminuit_errors.items(): self.minuit.errors[param] = value for param, value in iminuit_limits.items(): self.minuit.limits[param] = value for param, value in iminuit_fixed_parameters.items(): self.minuit.fixed[param] = value # This is to tell Minuit that we are dealing with likelihoods, # not chi square self.minuit.errordef = Minuit.LIKELIHOOD self.minuit.print_level = self.verbosity if user_setup_dict is not None: if "ftol" in user_setup_dict: self.minuit.tol = user_setup_dict["ftol"] else: # Do nothing and leave the default in iminuit pass self._best_fit_parameters = None self._function_minimum_value = None @staticmethod def _parameter_name_to_minuit_name(parameter): """ Translate the name of the parameter to the format accepted by Minuit :param parameter: the parameter name, of the form source.component.shape.parname :return: a minuit-friendly name for the parameter, such as source_component_shape_parname """ return parameter.replace(".", "_") # Override this because minuit uses different names
[docs] def restore_best_fit(self): """ Set the parameters back to their best fit value :return: none """ # Reset the internal value of all parameters super(MinuitMinimizer, self).restore_best_fit() # Update also the internal iminuit dictionary for k, par in self.parameters.items(): minuit_name = self._parameter_name_to_minuit_name(k) self.minuit.values[minuit_name] = par._get_internal_value()
def _print_current_status(self): """ To be used to print info before raising an exception :return: """ log.error("Last status:") for line in str(self._last_migrad_results.fmin).splitlines(): log.error(line) # Print params to get some info about the failure for line in str(self.minuit.params).splitlines(): log.error(line) def _minimize(self): """ Minimize the function using MIGRAD :param compute_covar: whether to compute the covariance (and error estimates) or not :return: best_fit: a dictionary containing the parameters at their best fit values function_minimum : the value for the function at the minimum NOTE: if the minimization fails, the dictionary will be empty and the function_minimum will be set to minimization.FIT_FAILED """ # Try a maximum of 10 times and break as soon as the fit is ok self.minuit.reset() self._last_migrad_results = self.minuit.migrad() for i in range(9): if self.minuit.valid: break else: # Try again self._last_migrad_results = self.minuit.migrad() if not self.minuit.valid: self._print_current_status() raise FitFailed( "MIGRAD call failed. This is usually due to unconstrained parameters." ) else: # Gather the optimized values for all parameters from the internal # iminuit dictionary best_fit_values = [] for k, par in self.parameters.items(): minuit_name = self._parameter_name_to_minuit_name(k) best_fit_values.append(self.minuit.values[minuit_name]) return best_fit_values, self.minuit.fval # Override the default _compute_covariance_matrix def _compute_covariance_matrix(self, best_fit_values): self.minuit.hesse() try: covariance = np.array(self.minuit.covariance) except RuntimeError: # Covariance computation has failed # Print current status self._print_current_status() log.error( "HESSE failed. Most probably some of your parameters are unconstrained.") raise CannotComputeCovariance( ) return covariance
[docs] def get_errors(self): """ Compute asymmetric errors using MINOS (slow, but accurate) and print them. NOTE: this should be called immediately after the minimize() method :return: a dictionary containing the asymmetric errors for each parameter. """ self.restore_best_fit() if not self.minuit.valid: raise CannotComputeErrors( "MIGRAD results not valid, cannot compute errors." ) try: self.minuit.minos() except: self._print_current_status() raise MINOSFailed( "MINOS has failed. This is not necessarily a problem if:\n\n" "* There are unconstrained parameters (the error is undefined). This is usually signaled " "by an approximated error, printed after the fit, larger than the best fit value\n\n" "* The fit is very difficult, because of high correlation between parameters. This is " "signaled by values close to 1.0 or -1.0 in the correlation matrix printed after the " "fit step.\n\n" "In this cases you can check the contour plots with get_contours(). If you are using a " "user-defined model, you can also try to reformulate your model with less correlated " "parameters." ) # Make a list for the results errors = collections.OrderedDict() for k, par in self.parameters.items(): minuit_name = self._parameter_name_to_minuit_name(k) minus_error = self.minuit.merrors[minuit_name].lower plus_error = self.minuit.merrors[minuit_name].upper if par.has_transformation(): # Need to transform in the external reference best_fit_value_internal = self._fit_results.loc[par.path, "value"] _, minus_error_external = par.internal_to_external_delta( best_fit_value_internal, minus_error ) _, plus_error_external = par.internal_to_external_delta( best_fit_value_internal, plus_error ) else: minus_error_external = minus_error plus_error_external = plus_error errors[k] = (minus_error_external, plus_error_external) return errors