diff --git a/pisa/analysis/analysis.py b/pisa/analysis/analysis.py index dfa533f1e..37c66c0f7 100644 --- a/pisa/analysis/analysis.py +++ b/pisa/analysis/analysis.py @@ -7,6 +7,7 @@ from collections.abc import Sequence, Mapping from collections import OrderedDict from copy import deepcopy +from functools import partial from operator import setitem from itertools import product import re @@ -15,11 +16,15 @@ import warnings import numpy as np -import scipy.optimize as optimize +import scipy +from scipy import optimize # this is needed for the take_step option in basinhopping from scipy._lib._util import check_random_state +# to convert from dict constraint type for differential evolution +from scipy.optimize._constraints import old_constraint_to_new from iminuit import Minuit import nlopt +from pkg_resources import parse_version import pisa from pisa import EPSILON, FTYPE, ureg @@ -34,13 +39,14 @@ LLH_METRICS, CHI2_METRICS, weighted_chi2, it_got_better, is_metric_to_maximize) -__all__ = ['MINIMIZERS_USING_SYMM_GRAD', 'MINIMIZERS_USING_CONSTRAINTS', +__all__ = ['MINIMIZERS_USING_SYMM_GRAD', 'MINIMIZERS_ACCEPTING_CONSTRS', + 'scipy_constraints_to_callables', 'get_nlopt_inequality_constraint_funcs', 'set_minimizer_defaults', 'validate_minimizer_settings', 'Counter', 'Analysis', 'BasicAnalysis'] -__author__ = 'J.L. Lanfranchi, P. Eller, S. Wren, E. Bourbeau, A. Trettin' +__author__ = 'J.L. Lanfranchi, P. Eller, S. Wren, E. Bourbeau, A. Trettin, T. Ehrhardt' -__license__ = '''Copyright (c) 2014-2020, The IceCube Collaboration +__license__ = '''Copyright (c) 2014-2024, The IceCube Collaboration Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. @@ -54,14 +60,136 @@ See the License for the specific language governing permissions and limitations under the License.''' +SUPPORTED_LOCAL_SCIPY_MINIMIZERS=( + 'cobyla', 'l-bfgs-b', 'nelder-mead', 'slsqp', 'trust-constr' +) +"""Local scipy minimizers that PISA currently supports via this interface.""" MINIMIZERS_USING_SYMM_GRAD = ('l-bfgs-b', 'slsqp') """Minimizers that use symmetrical steps on either side of a point to compute gradients. See https://github.com/scipy/scipy/issues/4916""" -MINIMIZERS_USING_CONSTRAINTS = ('cobyla') -"""Minimizers that cannot use the 'bounds' argument and instead need bounds to -be formulated in terms of constraints.""" +MINIMIZERS_ACCEPTING_CONSTRS = ('cobyla', 'slsqp', 'trust-constr', 'cobyqa') +"""Minimizers that allow constraints to be passed. According to +scipy docs, cobyla and slsqp require dictionaries, whereas +trust-constr and cobyqa require constraints to be passed as `LinearConstraint` +or `NonlinearConstraint` instances. However, as of now, the conversion to the +form required by the minimizer is done internally by scipy. Hence, this global +variable merely serves documentation purposes right now. Note that cobyqa is +only added in scipy version 1.14.0.""" + +EVAL_MSG = ('Using eval() is potentially dangerous as it can execute ' + 'arbitrary code! Do not store your config file in a place ' + 'where others have writing access!') +"""Repeatedly implemented warning about evaluating expressions representing + (in)equality constraints.""" + +# TODO: Observed or known scipy minimization issues that might be fixable with scipy updates: +# * SHGO ignores various local minimizer options (https://github.com/scipy/scipy/issues/20028) +# * unreliable global scipy minimization with constraints: non-negligible constraint +# violations or stepping out of bounds + + +def scipy_constraints_to_callables(constr_dicts, hypo_maker): + """Convert constraints expressions in terms of ParamSets + into callables for scipy. Overwrites "fun" entries + in `constr_dicts`. + """ + def constr_func(x, constr_func_params): + hypo_maker._set_rescaled_free_params(x) + if hypo_maker.__class__.__name__ == "Detectors": + # updates values for ALL detectors + update_param_values_detector(hypo_maker, hypo_maker.params.free) + return constr_func_params(hypo_maker.params) + logging.warning(EVAL_MSG) + assert isinstance(constr_dicts, Sequence) + for cd in constr_dicts: + assert isinstance(cd, Mapping) + # the equality constraint is specified as a function that takes a + # ParamSet as its input + assert "fun" in cd + constr = cd["fun"] + logging.debug(f"adding scipy constraint: {constr}") + if callable(constr): + constr_func_params = constr + else: + constr_func_params = eval(constr) + t = type(constr_func_params) + if not callable(constr_func_params): + raise TypeError(f"Evaluated object not a callable, but {t}.") + # overwrite + cd["fun"] = partial(constr_func, constr_func_params=constr_func_params) + + +def make_scipy_local_minimizer_kwargs(minimizer_settings, constrs=None, bounds=None): + """Small helper function containing common logic for + creating minimizer keyword args in calls to global + routines via their scipy interface.""" + minimizer_kwargs = deepcopy(minimizer_settings) + minimizer_kwargs["method"] = minimizer_settings["method"]["value"] + minimizer_kwargs["options"] = minimizer_settings["options"]["value"] + if constrs is not None: + minimizer_kwargs["constraints"] = constrs + if bounds is not None: + minimizer_kwargs["bounds"] = bounds + return minimizer_kwargs + + +def get_nlopt_inequality_constraint_funcs(method_kwargs, hypo_maker): + """Convert constraints expressions in terms of ParamSets + into callables for nlopt. Evals expression(s) from `method_kwargs` + and returns list of callables. + """ + def ineq_func(x, grad, ineq_func_params): + if grad.size > 0: + raise RuntimeError("gradients not supported") + hypo_maker._set_rescaled_free_params(x) + if hypo_maker.__class__.__name__ == "Detectors": + # updates values for ALL detectors + update_param_values_detector(hypo_maker, hypo_maker.params.free) + # In NLOPT, the inequality function must stay negative, while in + # scipy, the inequality function must stay positive. We keep with + # the scipy convention by flipping the sign. + return -ineq_func_params(hypo_maker.params) + assert "ineq_constraints" in method_kwargs + logging.warning(EVAL_MSG) + constr_list = method_kwargs["ineq_constraints"] + if isinstance(constr_list, str): + constr_list = [constr_list] + ineq_funcs = [] + for constr in constr_list: + # the inequality function is specified as a function that takes a + # ParamSet as its input + logging.debug(f"adding nlopt constraint (must stay positive): {constr}") + if callable(constr): + ineq_func_params = constr + else: + ineq_func_params = eval(constr) + t = type(ineq_func_params) + if not callable(ineq_func_params): + raise TypeError(f"Evaluated object not a callable, but {t}.") + ineq_funcs.append(partial(ineq_func, ineq_func_params=ineq_func_params)) + return ineq_funcs + + +def make_scipy_constraint_dict(constr_type, fun, jac=None, args=None): + """Makes a constraint dictionary in the form accepted by scipy, + see e.g. https://docs.scipy.org/doc/scipy-1.13.1/reference/generated/scipy.optimize.minimize.html""" + assert constr_type in ["eq", "ineq"] + t = type(fun) + if not callable(fun): + raise TypeError("Constraint function has to be callable, not {t}.") + constr_dict = {'type': constr_type, 'fun': fun} + if jac is not None: + t = type(jac) + if not callable(jac): + raise TypeError("Jacobian has to be callable, not {t}.") + constr_dict['jac'] = jac + if args is not None: + assert isinstance(args, Sequence) + constr_dict['args'] = args + return constr_dict + def merge_mapsets_together(mapset_list=None): '''Handle merging of multiple MapSets, when they come in @@ -88,7 +216,6 @@ def merge_mapsets_together(mapset_list=None): return new_dict -# TODO: add Nelder-Mead, as it was used previously... def set_minimizer_defaults(minimizer_settings): """Fill in default values for minimizer settings. @@ -135,9 +262,19 @@ def set_minimizer_defaults(minimizer_settings): opt_defaults.update(dict( rhobeg=0.1, maxiter=1000, tol=1e-4, )) + elif method == 'cobyqa': + # just make this solver available for now + pass + elif method == 'trust-constr': + opt_defaults.update(dict( + maxiter=200, gtol=1e-4, xtol=1e-4, barrier_tol=1e-4 + )) + elif method == 'nelder-mead': + opt_defaults.update(dict( + maxfev=1000, xatol=1e-4, fatol=1e-4 + )) else: - raise ValueError('Unhandled minimizer "%s" / FTYPE=%s' - % (method, FTYPE)) + raise ValueError(f'Unhandled minimizer "{method}" / FTYPE={FTYPE}') opt_defaults.update(new_minimizer_settings['options']['value']) @@ -151,10 +288,11 @@ def set_minimizer_defaults(minimizer_settings): return new_minimizer_settings -# TODO: add Nelder-Mead, as it was used previously... def validate_minimizer_settings(minimizer_settings): """Validate minimizer settings. + Supported minimizers are the same as in `set_minimizer_defaults`. + See source for specific thresholds set. Parameters @@ -170,6 +308,7 @@ def validate_minimizer_settings(minimizer_settings): ftype_eps = np.finfo(FTYPE).eps method = minimizer_settings['method']['value'].lower() options = minimizer_settings['options']['value'] + if method == 'l-bfgs-b': must_have = ('maxcor', 'ftol', 'gtol', 'eps', 'maxfun', 'maxiter', 'maxls') @@ -181,7 +320,29 @@ def validate_minimizer_settings(minimizer_settings): 'iprint', 'disp', 'callback') elif method == 'cobyla': must_have = ('maxiter', 'rhobeg', 'tol') - may_have = must_have + ('disp', 'catol') + may_have = must_have + ('disp', 'catol', 'constraints') + elif method == 'cobyqa': + assert parse_version(scipy.__version__) >= parse_version('1.14.0') + must_have = () + may_have = must_have + ('disp', 'maxiter', 'maxfev', 'f_target', + 'feasibility_tol', 'initial_tr_radius', + 'final_tr_radius', 'scale', 'constraints') + elif method == 'trust-constr': + must_have = ('maxiter', 'gtol', 'xtol', 'barrier_tol') + may_have = must_have + ('sparse_jacobian', 'initial_tr_radius', + 'initial_constr_penalty', 'constraints', + 'initial_barrier_parameter', + 'initial_barrier_tolerance', + 'factorization_method', + 'finite_diff_rel_step', + 'verbose', 'disp') + elif method == 'nelder-mead': + must_have = ('maxfev', 'xatol', 'fatol') + may_have = must_have + ('disp', 'maxiter', 'return_all', + 'initial_simplex', 'adaptive', + 'bounds') + else: + raise ValueError(f'Unhandled minimizer "{method}" / FTYPE={FTYPE}') missing = set(must_have).difference(set(options)) excess = set(options).difference(set(may_have)) @@ -220,7 +381,7 @@ def validate_minimizer_settings(minimizer_settings): if val > warn_lim: logging.warning(eps_gt_msg, method, 'eps', val, warn_lim) - if method == 'slsqp': + elif method == 'slsqp': err_lim, warn_lim = 2, 10 val = options['ftol'] if val < err_lim * ftype_eps: @@ -243,7 +404,7 @@ def validate_minimizer_settings(minimizer_settings): if val > warn_lim: logging.warning(eps_gt_msg, method, 'eps', val, warn_lim) - if method == 'cobyla': + elif method == 'cobyla': if options['rhobeg'] > 0.5: raise ValueError('starting step-size > 0.5 will overstep boundary') if options['rhobeg'] < 1e-2: @@ -393,7 +554,7 @@ def update_param_values_detector( hypo_maker.init_params() # TODO: move this to a central location prob. in utils -class Counter(object): +class Counter(): """Simple counter object for use as a minimizer callback.""" def __init__(self, i=0): self._count = i @@ -416,7 +577,7 @@ def count(self): """int : Current count""" return self._count -class BoundedRandomDisplacement(object): +class BoundedRandomDisplacement(): """ Add a bounded random displacement of maximum size `stepsize` to each coordinate Calling this updates `x` in-place. @@ -441,7 +602,7 @@ def __call__(self, x): x = np.clip(x, *self.bounds) # bounds are automatically broadcast return x -class HypoFitResult(object): +class HypoFitResult(): """Holds all relevant information about a fit result.""" @@ -706,7 +867,7 @@ def deserialize_detailed_metric_info(info_dict): return detailed_metric_info -class BasicAnalysis(object): +class BasicAnalysis(): """A bare-bones analysis that only fits a hypothesis to data. Full analyses with functionality beyond just fitting (doing scans, for example) @@ -1060,7 +1221,7 @@ def fit_recursively( if method in ["fit_octants", "fit_ranges"]: method = method.split("_")[1] - logging.warn(f"fit method 'fit_{method}' has been re-named to '{method}'") + logging.warning(f"fit method 'fit_{method}' has been re-named to '{method}'") # If made it here, we have a fit to do... fit_function = getattr(self, f"_fit_{method}") @@ -1230,11 +1391,7 @@ def _fit_condition(self, data_dist, hypo_maker, metric, assert len(local_fit_kwargs) == 2, ("need to fit specs, first runs if True, " "second runs if false") if type(method_kwargs["condition_func"]) is str: - logging.warn( - "Using eval() is potentially dangerous as it can execute " - "arbitrary code! Do not store your config file in a place" - "where others have writing access!" - ) + logging.warning(EVAL_MSG) condition_func = eval(method_kwargs["condition_func"]) assert callable(condition_func), "evaluated object is not a valid function" elif callable(method_kwargs["condition_func"]): @@ -1408,11 +1565,7 @@ def _fit_constrained(self, data_dist, hypo_maker, metric, logging.info("entering constrained fit...") if type(method_kwargs["ineq_func"]) is str: - logging.warn( - "Using eval() is potentially dangerous as it can execute " - "arbitrary code! Do not store your config file in a place" - "where others have writing access!" - ) + logging.warning(EVAL_MSG) ineq_func = eval(method_kwargs["ineq_func"]) assert callable(ineq_func), "evaluated object is not a valid function" elif callable(method_kwargs["ineq_func"]): @@ -1665,43 +1818,62 @@ def _fit_scipy(self, data_dist, hypo_maker, metric, global_method = method_kwargs["global_method"] if local_fit_kwargs is not None and global_method not in methods_using_local_fits: - logging.warn(f"local_fit_kwargs are ignored by global method {global_method}") + logging.warning(f"local_fit_kwargs are ignored by global method {global_method}") if global_method is None: logging.info(f"entering local scipy fit using {method_kwargs['method']['value']}") else: - assert global_method in global_scipy_methods, "unsupported global fit method" + if not global_method in global_scipy_methods: + raise ValueError(f"Unsupported global fit method {global_method}") logging.info(f"entering global scipy fit using {global_method}") if not self.blindness: logging.debug("free parameters:") logging.debug(hypo_maker.params.free) if global_method in methods_using_local_fits: - minimizer_settings = set_minimizer_defaults(local_fit_kwargs) - validate_minimizer_settings(minimizer_settings) + if local_fit_kwargs is not None: + minimizer_settings = set_minimizer_defaults(local_fit_kwargs) + validate_minimizer_settings(minimizer_settings) + else: + # We put this here such that the checks farther down don't crash + minimizer_settings = {"method": {"value": "none"}, "options": {"value": {}}} elif global_method == "differential_evolution": - # Unfortunately we are not allowed to pass these, DE with polish=True always - # uses L-BFGS-B with default settings. - if ("polish" in method_kwargs["options"].keys() - and method_kwargs["options"]["polish"]): - logging.info("Differential Evolution result will be polished using L-BFGS-B") - # We need to put the method here so that the bounds will be adjusted - # below, otherwise the polishing fit can cause crashes if it hits the - # bounds. - minimizer_settings = { - "method": {"value": "L-BFGS-B"}, - "options": {"value": {"eps": 1e-8}} - } + # unfortunately we are not allowed to configure local min. + opt = method_kwargs["options"] + if local_fit_kwargs is not None and "constraints" in local_fit_kwargs["options"]["value"]: + raise RuntimeError("Pass constraints to differential evolution only via method_kwargs!") + # polish is default + if "constraints" in opt and opt["constraints"]: + # When we detect constraints (which are not empty or None), + # we assume polishing is requested (constraints useless otherwise) + opt["polish"] = True + if "polish" in opt and opt["polish"]: + if "constraints" in opt: + if opt["constraints"] is None: + opt["constraints"] = [] + if opt["constraints"]: + logging.info( + "Differential evolution result is forced to be polished " + "with trust-constr." + ) + minimizer_settings = { + "method": {"value": "trust-constr"}, + "options": {"value": {"constraints": opt.pop("constraints")}} + } + else: + logging.info( + "Differential evolution result is forced to be polished " + "with L-BFGS-B." + ) + minimizer_settings = { + "method": {"value": "l-bfgs-b"}, + "options": {"value": {"eps": 1e-8}} + } else: + # Remove any possible constraints key which is empty + opt.pop("constraints", None) # We put this here such that the checks farther down don't crash - minimizer_settings = {"method": {"value": "None"}} - elif global_method == "dual_annealing": - minimizer_settings = { - "method": {"value": "L-BFGS-B"}, - "options": {"value": {"eps": 1e-8}} - } - logging.info("Due to a scipy bug, local minimization can only use default" - "L-BFGS-B settings. The given settings are ignored.") + minimizer_settings = {"method": {"value": "none"}, "options": {"value": {}}} else: minimizer_settings = set_minimizer_defaults(method_kwargs) validate_minimizer_settings(minimizer_settings) @@ -1724,31 +1896,38 @@ def _fit_scipy(self, data_dist, hypo_maker, metric, flip_x0 = np.zeros(len(x0), dtype=bool) minimizer_method = minimizer_settings["method"]["value"].lower() - cons = () - if minimizer_method in MINIMIZERS_USING_CONSTRAINTS: + + if ("constraints" in minimizer_settings["options"]["value"] + and minimizer_settings["options"]["value"]["constraints"]): + # convert user-specified equality or inequality constraint expressions + scipy_constraints_to_callables(minimizer_settings["options"]["value"]["constraints"], hypo_maker) + constrs = minimizer_settings["options"]["value"].pop("constraints") + else: + constrs = [] + + if minimizer_method == 'cobyla' and parse_version(scipy.__version__) < parse_version('1.11.0'): logging.warning( 'Minimizer %s requires bounds to be formulated in terms of constraints.' ' Constraining functions are auto-generated now.', minimizer_method ) - cons = [] + bound_constrs = [] for idx in range(len(x0)): - l = {'type': 'ineq', - 'fun': lambda x, i=idx: x[i] - FTYPE_PREC} # lower bound at zero - u = {'type': 'ineq', - 'fun': lambda x, i=idx: 1. - x[i]} # upper bound at 1 - cons.append(l) - cons.append(u) + fl = lambda x, i=idx: x[i] - FTYPE_PREC # lower bound at zero + l = make_scipy_constraint_dict(constr_type='ineq', fun=fl) + fu = lambda x, i=idx: 1. - x[i] # upper bound at 1 + u = make_scipy_constraint_dict(constr_type='ineq', fun=fu) + bound_constrs.append(l) + bound_constrs.append(u) + constrs += bound_constrs # The minimizer begins with a step of size `rhobeg` in the positive # direction. Flipping around 0.5 ensures that this initial step will not # overstep boundaries if `rhobeg` is 0.5. flip_x0 = np.array(x0) > 0.5 # The minimizer can't handle bounds, but they still need to be passed for # the interface to be uniform even though they are not used. - bounds = [(0, 1)]*len(x0) - else: - bounds = [(0, 1)]*len(x0) + bounds = [(0, 1)]*len(x0) x0 = np.where(flip_x0, 1 - x0, x0) if global_method is None: @@ -1775,18 +1954,18 @@ def _fit_scipy(self, data_dist, hypo_maker, metric, # Before starting minimization, check if we already have a perfect match between data and template # This can happen if using pseudodata that was generated with the nominal values for parameters # (which will also be the initial values in the fit) and blah... - # If this is the case, don't both to fit and return results right away. + # If this is the case, don't both to fit and return results right away. # Grab the hypo map hypo_asimov_dist = hypo_maker.get_outputs(return_sum=True) # Check if the hypo matches data - matches = False + matches = False if isinstance(data_dist, list): matches = all( entry.allclose(hypo_asimov_dist[ie]) for ie, entry in enumerate(data_dist) ) else: matches = data_dist.allclose(hypo_asimov_dist) - + if matches: msg = 'Initial hypo matches data, no need for fit' @@ -1824,7 +2003,7 @@ def _fit_scipy(self, data_dist, hypo_maker, metric, args=(hypo_maker, data_dist, metric, counter, fit_history, flip_x0, external_priors_penalty), bounds=bounds, - constraints=cons, + constraints=constrs, method=minimizer_settings['method']['value'], options=minimizer_settings['options']['value'], callback=self._minimizer_callback @@ -1833,34 +2012,41 @@ def _fit_scipy(self, data_dist, hypo_maker, metric, optimize_result = optimize.differential_evolution( func=self._minimizer_callable, bounds=bounds, + constraints=[old_constraint_to_new(i, cd) for i, cd in enumerate(constrs)], args=(hypo_maker, data_dist, metric, counter, fit_history, flip_x0, external_priors_penalty), callback=self._minimizer_callback, **method_kwargs["options"] ) elif global_method == "basinhopping": - if "seed" in method_kwargs["options"]: - seed = method_kwargs["options"]["seed"] + opt = method_kwargs["options"] + if "seed" in opt: + seed = opt["seed"] else: seed = None rng = check_random_state(seed) - if "step_size" in method_kwargs["options"]: - step_size = method_kwargs["options"]["step_size"] + if "stepsize" in opt: + stepsize = opt["stepsize"] else: - step_size = 0.5 + stepsize = 0.5 - take_step = BoundedRandomDisplacement(step_size, bounds, rng) - minimizer_kwargs = deepcopy(local_fit_kwargs) + take_step = BoundedRandomDisplacement(stepsize, bounds, rng) + if minimizer_method != "none": + minimizer_kwargs = make_scipy_local_minimizer_kwargs( + minimizer_settings=minimizer_settings, + constrs=constrs, + bounds=bounds + ) + else: + minimizer_kwargs = {"bounds": bounds} + # we need to pass args and bounds in any case as part of minimizer_kwargs minimizer_kwargs["args"] = ( hypo_maker, data_dist, metric, counter, fit_history, flip_x0, external_priors_penalty ) if "reset_free" in minimizer_kwargs: del minimizer_kwargs["reset_free"] - minimizer_kwargs["method"] = local_fit_kwargs["method"]["value"] - minimizer_kwargs["options"] = local_fit_kwargs["options"]["value"] - minimizer_kwargs["bounds"] = bounds def basinhopping_callback(x, f, accept): self._nit += 1 optimize_result = optimize.basinhopping( @@ -1871,44 +2057,50 @@ def basinhopping_callback(x, f, accept): minimizer_kwargs=minimizer_kwargs, **method_kwargs["options"] ) - optimize_result.success = True # basinhopping doesn't set this property elif global_method == "dual_annealing": + minimizer_kwargs = make_scipy_local_minimizer_kwargs( + minimizer_settings=minimizer_settings, + constrs=constrs, + bounds=bounds + ) + if "reset_free" in minimizer_kwargs: + del minimizer_kwargs["reset_free"] def annealing_callback(x, f, context): self._nit += 1 - # TODO: Enable use of custom minimization if scipy is fixed - # The scipy implementation is buggy insofar as it doesn't apply bounds to - # the inner minimization and there is no way to pass bounds through that - # doesn't crash. This leads to evaluations outside of the bounds. optimize_result = optimize.dual_annealing( func=self._minimizer_callable, bounds=bounds, x0=x0, args=(hypo_maker, data_dist, metric, counter, fit_history, flip_x0, external_priors_penalty), + minimizer_kwargs=minimizer_kwargs if minimizer_method != "none" else {}, callback=annealing_callback, **method_kwargs["options"] ) elif global_method == "shgo": - minimizer_kwargs = deepcopy(local_fit_kwargs) - minimizer_kwargs["args"] = ( - hypo_maker, data_dist, metric, counter, fit_history, - flip_x0, external_priors_penalty + minimizer_kwargs = make_scipy_local_minimizer_kwargs( + minimizer_settings=minimizer_settings, + constrs=constrs, + bounds=bounds ) if "reset_free" in minimizer_kwargs: del minimizer_kwargs["reset_free"] - minimizer_kwargs["method"] = local_fit_kwargs["method"]["value"] - minimizer_kwargs["options"] = local_fit_kwargs["options"]["value"] + if minimizer_kwargs["method"] != "none": + logging.warning( + f"Due to a scipy bug, shgo will ignore many local minimiser " + "options. This will most likely result in unreliable " + "behaviour or even crashes. Refer to " + "https://github.com/scipy/scipy/issues/20028." + ) optimize_result = optimize.shgo( func=self._minimizer_callable, bounds=bounds, args=(hypo_maker, data_dist, metric, counter, fit_history, flip_x0, external_priors_penalty), callback=self._minimizer_callback, - minimizer_kwargs=minimizer_kwargs, + minimizer_kwargs=minimizer_kwargs if minimizer_method != "none" else {}, **method_kwargs["options"] ) - else: - raise ValueError("Unsupported global fit method") end_t = time.time() if self.pprint: @@ -1924,7 +2116,7 @@ def annealing_callback(x, f, context): msg = '' else: msg = ' ' + str(optimize_result.message) - logging.warn('Optimization failed.' + msg) + logging.warning('Optimization failed.' + msg) # Instead of crashing completely, return a fit result with an infinite # test statistic value. metadata = {"success":optimize_result.success, "message":optimize_result.message,} @@ -2036,8 +2228,8 @@ def _fit_iminuit(self, data_dist, hypo_maker, metric, logging.info("Entering local fit using Minuit") if local_fit_kwargs is not None: - logging.warn("Local fit kwargs are ignored by 'fit_minuit'." - "Use 'method_kwargs' to set Minuit options.") + logging.warning("Local fit kwargs are ignored by 'fit_minuit'." + "Use 'method_kwargs' to set Minuit options.") if method_kwargs is None: method_kwargs = {} # use all defaults @@ -2082,7 +2274,7 @@ def loss_func(x): # In rare circumstances, minuit will try setting one of the parameters # to NaN. Minuit might be able to recover when we return NaN. if np.any(~np.isfinite(x)): - logging.warn(f"Minuit tried evaluating at invalid parameters: {x}") + logging.warning(f"Minuit tried evaluating at invalid parameters: {x}") return np.nan return self._minimizer_callable(x, *args) @@ -2145,9 +2337,9 @@ def loss_func(x): ) if not m.accurate: - logging.warn("Covariance matrix invalid.") + logging.warning("Covariance matrix invalid.") if not m.valid: - logging.warn("Minimum not valid according to Minuit's criteria.") + logging.warning("Minimum not valid according to Minuit's criteria.") # Will not assume that the minimizer left the hypo maker in the # minimized state, so set the values now (also does conversion of @@ -2247,10 +2439,10 @@ def _fit_nlopt(self, data_dist, hypo_maker, metric, logging.info("Entering fit using NLOPT") if local_fit_kwargs is not None: - logging.warn("`local_fit_kwargs` are ignored by 'fit_nlopt'." - "Use `method_kwargs` to set nlopt options and use " - "`method_kwargs['local_optimizer']` to define the settings of " - "a subsidiary NLOPT optimizer.") + logging.warning("`local_fit_kwargs` are ignored by 'fit_nlopt'." + "Use `method_kwargs` to set nlopt options and use " + "`method_kwargs['local_optimizer']` to define the settings of " + "a subsidiary NLOPT optimizer.") if method_kwargs is None: raise ValueError("Need to specify at least the algorithm to run.") @@ -2291,7 +2483,7 @@ def _fit_nlopt(self, data_dist, hypo_maker, metric, def loss_func(x, grad): if np.any(~np.isfinite(x)): - logging.warn(f"NLOPT tried evaluating at invalid parameters: {x}") + logging.warning(f"NLOPT tried evaluating at invalid parameters: {x}") return np.nan if grad.size > 0: raise RuntimeError("Gradients cannot be calculated, use a gradient-free" @@ -2361,7 +2553,7 @@ def loss_func(x, grad): if self.blindness < 2: metadata["rescaled_values"] = rescaled_pvals else: - metadata["rescaled_values"] = np.full(len(m.values), np.nan) + metadata["rescaled_values"] = np.full(len(m.values), np.nan) #FIXME: m.values? # we don't get a Hessian from nlopt metadata["hess_inv"] = np.full((len(x0), len(x0)), np.nan) @@ -2432,30 +2624,8 @@ def _define_nlopt_opt(self, method_kwargs, loss_func, hypo_maker): for k, v in method_kwargs["algorithm_params"].items(): opt.set_param(k, v) if "ineq_constraints" in method_kwargs.keys(): - logging.warn( - "Using eval() is potentially dangerous as it can execute " - "arbitrary code! Do not store your config file in a place" - "where others have writing access!" - ) - constr_list = method_kwargs["ineq_constraints"] - if isinstance(constr_list, str): - constr_list = [constr_list] - for constr in constr_list: - # the inequality function is specified as a function that takes a - # ParamSet as its input - logging.info(f"adding constraint (must stay positive): {constr}") - ineq_func_params = eval(constr) - assert callable(ineq_func_params), "evaluated object is not a valid function" - def ineq_func(x, grad): - if grad.size > 0: - raise RuntimeError("gradients not supported") - hypo_maker._set_rescaled_free_params(x) - if hypo_maker.__class__.__name__ == "Detectors": - update_param_values_detector(hypo_maker, hypo_maker.params.free) #updates values for ALL detectors - # In NLOPT, the inequality function must stay negative, while in - # scipy, the inequality function must stay positive. We keep with - # the scipy convention by flipping the sign. - return -ineq_func_params(hypo_maker.params) + ineq_funcs = get_nlopt_inequality_constraint_funcs(method_kwargs, hypo_maker) + for ineq_func in ineq_funcs: opt.add_inequality_constraint(ineq_func) # Population size for stochastic search algorithms, see @@ -2670,7 +2840,7 @@ def _minimizer_callable(self, scaled_param_vals, hypo_maker, data_dist, return sign*metric_val - def _minimizer_callback(self, xk, **unused_kwargs): # pylint: disable=unused-argument + def _minimizer_callback(self, xk, *unused_args, **unused_kwargs): # pylint: disable=unused-argument """Passed as `callback` parameter to `optimize.minimize`, and is called after each iteration. Keeps track of number of iterations. @@ -3704,6 +3874,362 @@ def _fit_nonsense( logging.info('<< PASS : test_basic_analysis >>') + +def test_constrained_minimization(pprint=False): + """Test scipy solvers without or with equality and inequality + constraints. All are run with default options as set by + `set_minimizer_defaults`.""" + from pisa.core.distribution_maker import DistributionMaker + + config = 'settings/pipeline/fast_example.cfg' + dm = DistributionMaker(config) + data_dist = dm.get_outputs(return_sum=True) + + ### slsqp test with constraints ### + def slsqp_constr(): + ana = BasicAnalysis() + ana.pprint = pprint + min_sett = { + "method": {"value": "slsqp", "desc": ""}, + "options": {"value": {}, "desc": {}} + } + + min_delta_index = 5e-3 + max_aeff_scale = 0.986 + t23 = 44.2 + # constraint function can be callable or string + constrs_list = [ + {'type': 'ineq', + 'fun': lambda params: params.delta_index.m_as("dimensionless") - min_delta_index}, + {'type': 'ineq', + 'fun': 'lambda p: -p.aeff_scale.m_as("dimensionless") + %s' % max_aeff_scale}, + {'type': 'eq', + 'fun': lambda params: params.theta23.m_as("degree") - t23} + ] + min_sett["options"]["value"]["constraints"] = constrs_list + + scipy_sett = {"method": "scipy", "method_kwargs": min_sett} + + dm.params.theta23.randomize(random_state=1) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **scipy_sett + ) + + assert bf.minimizer_metadata['success'] + tol = 1e-5 + assert recursiveEquality(bf.params.theta23.m_as('degree'), t23) + assert bf.params.delta_index.m_as('dimensionless') >= min_delta_index - tol + assert bf.params.aeff_scale.m_as('dimensionless') <= max_aeff_scale + tol + + ### cobyla test with inequality constraints (doesn't support equalities) ### + def cobyla_constr(): + ana = BasicAnalysis() + ana.pprint = pprint + + min_t23 = 46. + min_aeff_scale = 1.02 + constrs_list = [ + {'type': 'ineq', + 'fun': lambda params: params.theta23.m_as("dimensionless") - min_t23}, + {'type': 'ineq', + 'fun': lambda params: params.aeff_scale.m_as("dimensionless") - min_aeff_scale} + ] + # FIXME: steps out of bounds whenever these are included and whether bounds are also implemented as constraints or not + min_sett = { + "method": {"value": "cobyla", "desc": ""}, + "options": {"value": {"constraints": constrs_list}, "desc": {}} + } + + scipy_sett = {"method": "scipy", "method_kwargs": min_sett} + + dm.params.randomize_free(random_state=1) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **scipy_sett + ) + assert bf.minimizer_metadata['success'] + tol = 1e-5 + assert bf.params.theta23.m_as('dimensionless') >= min_t23 - tol + assert bf.params.aeff_scale.m_as('dimensionless') >= min_aeff_scale - tol + + ### trust-constr test with constraints ### + def trust_constr_constr(): + ana = BasicAnalysis() + ana.pprint = pprint + + min_delta_index = 5e-3 + max_aeff_scale = 0.986 + t23 = 44.2 + constrs_list = [ + {'type': 'ineq', + 'fun': lambda params: params.delta_index.m_as("dimensionless") - min_delta_index}, + {'type': 'ineq', + 'fun': 'lambda p: -p.aeff_scale.m_as("dimensionless") + %s' % max_aeff_scale}, + {'type': 'eq', + 'fun': lambda params: params.theta23.m_as("degree") - t23} + ] + + min_sett = { + "method": {"value": "trust-constr", "desc": ""}, + "options": {"value": {"constraints": constrs_list}, "desc": {}} + } + scipy_sett = {"method": "scipy", "method_kwargs": min_sett} + + dm.params.randomize_free(random_state=5) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **scipy_sett + ) + assert bf.minimizer_metadata['success'] + tol = 1e-5 + assert recursiveEquality(bf.params.theta23.m_as('degree'), t23) + assert bf.params.delta_index.m_as('dimensionless') >= min_delta_index - tol + assert bf.params.aeff_scale.m_as('dimensionless') <= max_aeff_scale + tol + + ### Nelder-Mead test (no constraints supported) ### + def nm_unconstr(): + ana = BasicAnalysis() + ana.pprint = pprint + + min_sett = { + "method": {"value": "Nelder-Mead", "desc": ""}, + "options": {"value": {}, "desc": {}} + } + scipy_sett = {"method": "scipy", "method_kwargs": min_sett} + + dm.params.randomize_free(random_state=9) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **scipy_sett + ) + assert bf.minimizer_metadata['success'] + + ### finally an nlopt solver with constraints ### + def some_nlopt_constr(): + ana = BasicAnalysis() + ana.pprint = pprint + + min_t23 = 46. + method_kwargs = { + "algorithm": "NLOPT_LN_COBYLA", + "ftol_abs": 1e-3, + "ftol_rel": 1e-3, + "maxeval": 100, + "ineq_constraints": [ + lambda params: params.theta23.m_as("degree") - min_t23 + ] + } + nlopt_sett = {"method": "nlopt", "method_kwargs": method_kwargs} + + # fix 2/3 to make it converge faster + fix_params = ["delta_index", "aeff_scale"] + [dm.params.fix(p) for p in fix_params] # pylint: disable=expression-not-assigned + + dm.params.randomize_free(random_state=11) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **nlopt_sett + ) + assert bf.minimizer_metadata['success'] + tol = 1e-5 + assert bf.params.theta23.m_as('degree') >= min_t23 - tol + + # unfix again + [dm.params.unfix(p) for p in fix_params] # pylint: disable=expression-not-assigned + + # now run them all + slsqp_constr() + try: + cobyla_constr() + except Exception as e: + # don't fail, just document here + logging.error("Cobyla test with constraints failed: '%s'. This is under investigation..." % str(e)) + trust_constr_constr() + nm_unconstr() + some_nlopt_constr() + + logging.info('<< PASS : test_constrained_minimization >>') + + +def test_global_scipy_minimization(pprint=False): + from pisa.core.distribution_maker import DistributionMaker + config = 'settings/pipeline/fast_example.cfg' + dm = DistributionMaker(config) + dm.params.fix("theta23") # make it converge faster + data_dist = dm.get_outputs(return_sum=True) + + def run_global_with_supported_local( + global_method_kwargs, constrs_list=None, constrs_test=None + ): + ana = BasicAnalysis() + ana.pprint = pprint + + assert (constrs_list is None) == (constrs_test is None) + + glob_meth = global_method_kwargs["global_method"] + + for loc_min in set(SUPPORTED_LOCAL_SCIPY_MINIMIZERS).union(('none',)): + min_sett = { + "method": {"value": loc_min, "desc": ""}, + "options": {"value": {}, "desc": {}} + } + if loc_min in MINIMIZERS_ACCEPTING_CONSTRS and constrs_list is not None: + min_sett["options"]["value"]["constraints"] = deepcopy(constrs_list) + + glob_sett = { + "method": "scipy", + "method_kwargs": global_method_kwargs, + "local_fit_kwargs": min_sett if loc_min != 'none' else None + } + dm.reset_free() + dm.params.aeff_scale.randomize(random_state=0) + try: + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **glob_sett + ) + if not bf.minimizer_metadata['success']: + raise RuntimeError("Minimizer unsuccessful") + if loc_min in MINIMIZERS_ACCEPTING_CONSTRS and constrs_list is not None: + constrs_test(bf) + except Exception as e: + # don't fail, just document here + logging.error( + f"test of {glob_meth} + {loc_min} failed: '{e}'. This is under investigation..." + ) + + + def run_differential_evolution(constrs_list, constrs_test, polish=True): + ana = BasicAnalysis() + ana.pprint = pprint + + assert (constrs_list is None) == (constrs_test is None) + + global_method_kwargs = { + "global_method": "differential_evolution", + "options": { + "maxiter": 100, + "tol": 0.1, + "disp": False, + "seed": 0, + "polish": polish, + "constraints": deepcopy(constrs_list) + }, + "local_fit_kwargs": None + } + glob_sett = { + "method": "scipy", + "method_kwargs": global_method_kwargs, + } + dm.reset_free() + dm.params.aeff_scale.randomize(random_state=0) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **glob_sett + ) + if not bf.minimizer_metadata['success']: + raise RuntimeError("Minimizer unsuccessful") + if constrs_test is not None: + constrs_test(bf) + + + min_delta_index = 5e-3 + max_aeff_scale = 0.986 + # all solvers that accept constraints accept inequalities + constrs_list = [ + {'type': 'ineq', + 'fun': lambda params: params.delta_index.m_as("dimensionless") - min_delta_index}, + {'type': 'ineq', + 'fun': 'lambda p: -p.aeff_scale.m_as("dimensionless") + %s' % max_aeff_scale} + ] + def constrs_test(bf): + tol = 1e-5 + if (not bf.params.delta_index.m_as('dimensionless') >= min_delta_index - tol or + not bf.params.aeff_scale.m_as('dimensionless') <= max_aeff_scale + tol): + raise ValueError("inequality constraint(s) violated!") + + method_kwargs = { + "global_method": "basinhopping", + "options": { + "niter": 3, + "T": 1, + "seed": 0, + "stepsize": 0.1 + } + } + # basinhopping w/ constraints + run_global_with_supported_local(method_kwargs, constrs_list, constrs_test) + # basinhopping w/o (skip) + #run_global_with_supported_local(method_kwargs) + + method_kwargs = { + "global_method": "dual_annealing", + "options": { + "maxiter": 10, + "seed": 0 + } + } + # dual annealing w/ constraints + run_global_with_supported_local(method_kwargs, constrs_list, constrs_test) + # dual annealing w/o (skip) + #run_global_with_supported_local(method_kwargs) + + # DE with trust-constr (w/ constraints) + run_differential_evolution(constrs_list, constrs_test) + # DE with l-bfgs-b (w/o constraints) or without polishing + # entirely takes too long: + #run_differential_evolution(None, None, polish=True) + #run_differential_evolution(None, None, polish=False) + + def run_global_min_with_constraints_from_file(fit_sett): + from pisa.utils.fileio import from_file + ana = BasicAnalysis() + ana.pprint = pprint + + fit_sett = from_file(fit_sett) + + dm.params.unfix("theta23") + dm.reset_free() + dm.params.randomize_free(random_state=0) + bf = ana.fit_recursively( + data_dist=data_dist, + hypo_maker=dm, + metric="chi2", + external_priors_penalty=None, + **fit_sett + ) + if not bf.minimizer_metadata['success']: + raise RuntimeError("Minimizer unsuccessful") + + #run_global_min_with_constraints_from_file('settings/minimizer/de_2nd_octant_popsize_10_init_sobol_polish_tol1e-1_maxiter1000.json') + + logging.info('<< PASS : test_global_scipy_minimization >>') + + if __name__ == "__main__": set_verbosity(1) test_basic_analysis(pprint=True) + test_constrained_minimization(pprint=True) + test_global_scipy_minimization(pprint=True) diff --git a/pisa_examples/resources/settings/minimizer/de_2nd_octant_popsize_10_init_sobol_polish_tol1e-1_maxiter1000.json b/pisa_examples/resources/settings/minimizer/de_2nd_octant_popsize_10_init_sobol_polish_tol1e-1_maxiter1000.json new file mode 100644 index 000000000..89a2724b7 --- /dev/null +++ b/pisa_examples/resources/settings/minimizer/de_2nd_octant_popsize_10_init_sobol_polish_tol1e-1_maxiter1000.json @@ -0,0 +1,18 @@ +{ + "method": "scipy", + "method_kwargs": { + "global_method": "differential_evolution", + "options": { + "maxiter": 1000, + "tol": 0.1, + "popsize": 10, + "init": "sobol", + "disp": true, + "seed": 0, + "constraints": [ + {"type": "ineq", + "fun": "lambda p: p.theta23.m_as('degree') - 45."} + ] + } + } +}