From 1c81292f43a66069ccf39fd9d454cb062b34c454 Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 12:37:55 +0200 Subject: [PATCH 1/7] add option to use scipy or minuit modifier. Removed the custom modifier --- src/cabinetry/fit/__init__.py | 316 +++++++++------------------------- 1 file changed, 77 insertions(+), 239 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index ebb62d52..56a7e73f 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -41,7 +41,7 @@ def print_results(fit_results: FitResults) -> None: ) -def _fit_model_pyhf( +def _fit_model( model: pyhf.pdf.Model, data: List[float], *, @@ -49,7 +49,9 @@ def _fit_model_pyhf( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, + minimizer: Optional[Literal['minuit','scipy']] = 'minuit', strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> FitResults: @@ -65,32 +67,43 @@ def _fit_model_pyhf( data (List[float]): the data to fit the model to minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS algorithm for all parameters specified, defaults to None (does not run - MINOS) + MINOS). Requires ``minuit`` as the minimizer. init_pars (Optional[List[float]], optional): list of initial parameter settings, defaults to None (use ``pyhf`` suggested inits) fix_pars (Optional[List[bool]], optional): list of booleans specifying which parameters are held constant, defaults to None (use ``pyhf`` suggestion) par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``minuit``. strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details - see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to - None (use ``iminuit`` default of 0.1) - + see either ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance) if using + minuit, or ``scipy.optimize.minimize``. defaults to None (use ``iminuit`` + default of 0.1, or `scipy`` method specific defaults) Returns: FitResults: object storing relevant fit results """ _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings - pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(verbose=1)) - - # strategy=None is currently not supported in pyhf - # https://github.com/scikit-hep/pyhf/issues/1785 - strategy_kwarg = {"strategy": strategy} if strategy is not None else {} + # set up the optimiser + if minimizer == 'minuit': + pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(verbose=1)) + # strategy=None is currently not supported in pyhf + # https://github.com/scikit-hep/pyhf/issues/1785 + optimizer_kwarg = {"strategy": strategy} if strategy is not None else {} + elif minimizer == 'scipy': + pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.scipy_optimizer(verbose=1)) + optimizer_kwarg = {"solver_options": solver_options} if solver_options is not None else {} + else: + raise ValueError(f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}") + result, corr_mat, best_twice_nll, result_obj = pyhf.infer.mle.fit( data, model, @@ -103,22 +116,30 @@ def _fit_model_pyhf( return_result_obj=True, maxiter=maxiter, tolerance=tolerance, - **strategy_kwarg, + **optimizer_kwarg, ) - log.info(f"Migrad status:\n{result_obj.minuit.fmin}") - bestfit = pyhf.tensorlib.to_numpy(result[:, 0]) - # set errors for fixed parameters to 0 (see iminuit#762) - uncertainty = np.where( - result_obj.minuit.fixed, 0.0, pyhf.tensorlib.to_numpy(result[:, 1]) - ) labels = model.config.par_names - corr_mat = pyhf.tensorlib.to_numpy(corr_mat) best_twice_nll = float(best_twice_nll) # convert 0-dim np.ndarray to float - minos_results = ( - _run_minos(result_obj.minuit, minos, labels) if minos is not None else {} - ) + if minimizer == 'minuit': + log.info(f"Migrad status:\n{result_obj.minuit.fmin}") + bestfit = pyhf.tensorlib.to_numpy(result[:, 0]) + # set errors for fixed parameters to 0 (see iminuit#762) + uncertainty = np.where( + result_obj.minuit.fixed, 0.0, pyhf.tensorlib.to_numpy(result[:, 1]) + ) + corr_mat = pyhf.tensorlib.to_numpy(corr_mat) + minos_results = ( + _run_minos(result_obj.minuit, minos, labels) if minos is not None else {} + ) + + elif minimizer == 'scipy': + bestfit = pyhf.tensorlib.to_numpy(result) + # scipy does not return uncertainty or correlation results + uncertainty = np.zeros(bestfit.shape) + corr_mat = np.diag(no.ones(bestfit.shape)) + minos_results = None fit_results = FitResults( bestfit, @@ -128,198 +149,13 @@ def _fit_model_pyhf( best_twice_nll, minos_uncertainty=minos_results, ) - pyhf.set_backend(pyhf.tensorlib, initial_optimizer) # restore optimizer settings - return fit_results - -def _fit_model_custom( - model: pyhf.pdf.Model, - data: List[float], - *, - minos: Optional[Union[List[str], Tuple[str, ...]]] = None, - init_pars: Optional[List[float]] = None, - fix_pars: Optional[List[bool]] = None, - par_bounds: Optional[List[Tuple[float, float]]] = None, - strategy: Optional[Literal[0, 1, 2]] = None, - maxiter: Optional[int] = None, - tolerance: Optional[float] = None, -) -> FitResults: - """Uses ``iminuit`` directly to perform a maximum likelihood fit. - - Parameters set to be fixed in the model are held constant. The ``init_pars`` - argument allows to override the ``pyhf`` default initial parameter settings, the - ``fix_pars`` argument overrides which parameters are held constant, ``par_bounds`` - sets parameter bounds. - - Args: - model (pyhf.pdf.Model): the model to use in the fit - data (List[float]): the data to fit the model to - minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS - algorithm for all parameters specified, defaults to None (does not run - MINOS) - init_pars (Optional[List[float]], optional): list of initial parameter settings, - defaults to None (use ``pyhf`` suggested inits) - fix_pars (Optional[List[bool]], optional): list of booleans specifying which - parameters are held constant, defaults to None (use ``pyhf`` suggestion) - par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with - parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) - strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by - Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior - of strategy 0 with user-provided gradients and 1 otherwise) - maxiter (Optional[int], optional): allowed number of calls for minimization, - defaults to None (use ``pyhf`` default of 100,000) - tolerance (Optional[float]), optional): tolerance for convergence, for details - see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to - None (use ``iminuit`` default of 0.1) - - Raises: - ValueError: if minimization fails - - Returns: - FitResults: object storing relevant fit results - """ - _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings - pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(verbose=1)) - - # use parameter settings provided in function arguments if they exist, else defaults - init_pars = init_pars or model.config.suggested_init() - fix_pars = fix_pars or model.config.suggested_fixed() - par_bounds = par_bounds or model.config.suggested_bounds() - - labels = model.config.par_names - - def twice_nll_func(pars: np.ndarray) -> Any: - """The objective for minimization: twice the negative log-likelihood. - - The return value is float-like, but not always a float. The actual type depends - on the active ``pyhf`` backend. - - Args: - pars (np.ndarray): parameter values at which the NLL is evaluated - - Returns: - Any: twice the negative log-likelihood - """ - twice_nll = -2 * model.logpdf(pars, data) - return twice_nll[0] - - m = iminuit.Minuit(twice_nll_func, init_pars, name=labels) - m.fixed = fix_pars - m.limits = par_bounds - m.errordef = 1 - m.print_level = 1 - - if strategy is not None: - m.strategy = strategy - else: - # pick strategy like pyhf: 0 if backend provides autodiff gradients, otherwise 1 - m.strategy = 0 if pyhf.tensorlib.default_do_grad else 1 - - maxiter = maxiter or 100_000 - m.tol = tolerance or 0.1 # goal: EDM < 0.002*tol*errordef - - m.migrad(ncall=maxiter) - m.hesse() # use default call limit (consistent with pyhf) - - log.info(f"MINUIT status:\n{m.fmin}") - if not m.valid: - raise ValueError("Minimization failed, minimum is invalid.") - - bestfit = np.asarray(m.values) - # set errors for fixed parameters to 0 (see iminuit#762) - uncertainty = np.where(m.fixed, 0.0, m.errors) - corr_mat = m.covariance.correlation() # iminuit.util.Matrix, subclass of np.ndarray - best_twice_nll = m.fval - - minos_results = _run_minos(m, minos, labels) if minos is not None else {} - - fit_results = FitResults( - bestfit, - uncertainty, - labels, - corr_mat, - best_twice_nll, - minos_uncertainty=minos_results, - ) + log.debug(f"-2 log(L) = {fit_results.best_twice_nll:.6f} at best-fit point") + pyhf.set_backend(pyhf.tensorlib, initial_optimizer) # restore optimizer settings return fit_results -def _fit_model( - model: pyhf.pdf.Model, - data: List[float], - *, - minos: Optional[Union[List[str], Tuple[str, ...]]] = None, - init_pars: Optional[List[float]] = None, - fix_pars: Optional[List[bool]] = None, - par_bounds: Optional[List[Tuple[float, float]]] = None, - strategy: Optional[Literal[0, 1, 2]] = None, - maxiter: Optional[int] = None, - tolerance: Optional[float] = None, - custom_fit: bool = False, -) -> FitResults: - """Interface for maximum likelihood fits through ``pyhf.infer`` API or ``iminuit``. - - Parameters set to be fixed in the model are held constant. The ``init_pars`` - argument allows to override the ``pyhf`` default initial parameter settings, the - ``fix_pars`` argument overrides which parameters are held constant, ``par_bounds`` - sets parameter bounds. - - Args: - model (pyhf.pdf.Model): the model to use in the fit - data (List[float]): the data to fit the model to - minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS - algorithm for all parameters specified, defaults to None (does not run - MINOS) - init_pars (Optional[List[float]], optional): list of initial parameter settings, - defaults to None (use ``pyhf`` suggested inits) - fix_pars (Optional[List[bool]], optional): list of booleans specifying which - parameters are held constant, defaults to None (use ``pyhf`` suggestion) - par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with - parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) - strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by - Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior - of strategy 0 with user-provided gradients and 1 otherwise) - maxiter (Optional[int], optional): allowed number of calls for minimization, - defaults to None (use ``pyhf`` default of 100,000) - tolerance (Optional[float]), optional): tolerance for convergence, for details - see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to - None (use ``iminuit`` default of 0.1) - custom_fit (bool, optional): whether to use the ``pyhf.infer`` API or - ``iminuit``, defaults to False (using ``pyhf.infer``) - - Returns: - FitResults: object storing relevant fit results - """ - if not custom_fit: - # use pyhf infer API - fit_results = _fit_model_pyhf( - model, - data, - minos=minos, - init_pars=init_pars, - fix_pars=fix_pars, - par_bounds=par_bounds, - strategy=strategy, - maxiter=maxiter, - tolerance=tolerance, - ) - else: - # use iminuit directly - fit_results = _fit_model_custom( - model, - data, - minos=minos, - init_pars=init_pars, - fix_pars=fix_pars, - par_bounds=par_bounds, - strategy=strategy, - maxiter=maxiter, - tolerance=tolerance, - ) - log.debug(f"-2 log(L) = {fit_results.best_twice_nll:.6f} at best-fit point") - return fit_results - def _run_minos( minuit_obj: iminuit.Minuit, @@ -432,46 +268,43 @@ def fit( model: pyhf.pdf.Model, data: List[float], *, - minos: Optional[Union[str, List[str], Tuple[str, ...]]] = None, - goodness_of_fit: bool = False, + minos: Optional[Union[List[str], Tuple[str, ...]]] = None, init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, + minimizer: Optional[Literal['minuit','scipy']] = 'minuit', strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, - custom_fit: bool = False, -) -> FitResults: + ) -> FitResults: """Performs a maximum likelihood fit, reports and returns the results. - Depending on the ``custom_fit`` keyword argument, this uses either the - ``pyhf.infer`` API or ``iminuit`` directly. - - Args: - model (pyhf.pdf.Model): model to use in fit - data (List[float]): data (including auxdata) the model is fit to - minos (Optional[Union[str, List[str], Tuple[str, ...]]], optional): runs the - MINOS algorithm for all parameters specified, defaults to None (does not run - MINOS) - goodness_of_fit (bool, optional): calculate goodness of fit with a saturated - model (perfectly fits data with shapefactors in all bins), defaults to False + Args: + model (pyhf.pdf.Model): the model to use in the fit + data (List[float]): the data to fit the model to + minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS + algorithm for all parameters specified, defaults to None (does not run + MINOS). Requires ``minuit`` as the minimizer. init_pars (Optional[List[float]], optional): list of initial parameter settings, defaults to None (use ``pyhf`` suggested inits) fix_pars (Optional[List[bool]], optional): list of booleans specifying which parameters are held constant, defaults to None (use ``pyhf`` suggestion) par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``minuit``. strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to scipy + optimizer. See ``scipy.optimize.show_options()`` for all available options. maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details - see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to - None (use ``iminuit`` default of 0.1) - custom_fit (bool, optional): whether to use the ``pyhf.infer`` API or - ``iminuit``, defaults to False (using ``pyhf.infer``) - + see either ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance) if using + minuit, or ``scipy.optimize.minimize``. defaults to None (use ``iminuit`` + default of 0.1, or `scipy`` method specific defaults) Returns: FitResults: object storing relevant fit results """ @@ -489,10 +322,11 @@ def fit( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, + minimizer=minimizer, strategy=strategy, + solver_options=solver_options, maxiter=maxiter, tolerance=tolerance, - custom_fit=custom_fit, ) print_results(fit_results) @@ -515,9 +349,9 @@ def ranking( fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, - custom_fit: bool = False, ) -> RankingResults: """Calculates the impact of nuisance parameters on the parameter of interest (POI). @@ -542,13 +376,13 @@ def ranking( strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to None (use ``iminuit`` default of 0.1) - custom_fit (bool, optional): whether to use the ``pyhf.infer`` API or - ``iminuit``, defaults to False (using ``pyhf.infer``) Raises: ValueError: if no POI is found @@ -563,10 +397,11 @@ def ranking( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, + minimizer='minuit', strategy=strategy, + solver_options=solver_options maxiter=maxiter, tolerance=tolerance, - custom_fit=custom_fit, ) labels = model.config.par_names @@ -618,10 +453,11 @@ def ranking( init_pars=init_pars_ranking, fix_pars=fix_pars_ranking, par_bounds=par_bounds, + minimizer='minuit', strategy=strategy, + solver_options=solver_options, maxiter=maxiter, tolerance=tolerance, - custom_fit=custom_fit, ) poi_val = fit_results_ranking.bestfit[poi_index] parameter_impact = poi_val - nominal_poi @@ -661,9 +497,9 @@ def scan( fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, - custom_fit: bool = False, ) -> ScanResults: """Performs a likelihood scan over the specified parameter. @@ -688,13 +524,13 @@ def scan( strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to None (use ``iminuit`` default of 0.1) - custom_fit (bool, optional): whether to use the ``pyhf.infer`` API or - ``iminuit``, defaults to False (using ``pyhf.infer``) Raises: ValueError: if parameter is not found in model @@ -717,10 +553,11 @@ def scan( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, + minimizer='minuit', strategy=strategy, + solver_options=solver_options, maxiter=maxiter, tolerance=tolerance, - custom_fit=custom_fit, ) nominal_twice_nll = fit_results.best_twice_nll par_mle = fit_results.bestfit[par_index] @@ -755,10 +592,11 @@ def scan( init_pars=init_pars_scan, fix_pars=fix_pars, par_bounds=par_bounds, + minimizer='scipy', strategy=strategy, + solver_options=solver_options, maxiter=maxiter, tolerance=tolerance, - custom_fit=custom_fit, ) # subtract best-fit delta_nlls[i_par] = scan_fit_results.best_twice_nll - nominal_twice_nll From 5147623ddec533874369b7fed2896c32293632ed Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 14:08:35 +0200 Subject: [PATCH 2/7] add minimizer option to limit and significance --- src/cabinetry/fit/__init__.py | 52 +++++++++++++++++++++++++---------- 1 file changed, 38 insertions(+), 14 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index 56a7e73f..8d651e1f 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -399,7 +399,7 @@ def ranking( par_bounds=par_bounds, minimizer='minuit', strategy=strategy, - solver_options=solver_options + solver_options=solver_options, maxiter=maxiter, tolerance=tolerance, ) @@ -617,7 +617,9 @@ def limit( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, + minimizer: Optional[Literal['minuit','scipy']] = 'scipy', strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> LimitResults: @@ -651,9 +653,13 @@ def limit( parameters are held constant, defaults to None (use ``pyhf`` suggestion) par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``scipy``. strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details @@ -669,12 +675,18 @@ def limit( LimitResults: observed and expected limits, CLs values, and scanned points """ _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings - pyhf.set_backend( - pyhf.tensorlib, - pyhf.optimize.minuit_optimizer( - verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance - ), - ) + + if minimizer == 'minuit': + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance) + + elif minimizer == 'scipy': + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, solver_options=solver_options, maxiter=maxiter, tolerance=tolerance) + + else: + raise ValueError(f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}") + pyhf.set_backend(pyhf.tensorlib, optimizer) # use POI given by kwarg, fall back to POI specified in model poi_index = model_utils._poi_index(model, poi_name=poi_name) @@ -891,7 +903,9 @@ def significance( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, + minimizer: Optional[Literal['minuit','scipy']] = 'scipy', strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> SignificanceResults: @@ -910,9 +924,14 @@ def significance( parameters are held constant, defaults to None (use ``pyhf`` suggestion) par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``scipy``. strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details @@ -923,12 +942,17 @@ def significance( SignificanceResults: observed and expected p-values and significances """ _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings - pyhf.set_backend( - pyhf.tensorlib, - pyhf.optimize.minuit_optimizer( - verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance - ), - ) + if minimizer == 'minuit': + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance) + + elif minimizer == 'scipy': + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, solver_options=solver_options, maxiter=maxiter, tolerance=tolerance) + else: + raise ValueError(f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}") + pyhf.set_backend(pyhf.tensorlib, optimizer) + # use POI given by kwarg, fall back to POI specified in model poi_index = model_utils._poi_index(model, poi_name=poi_name) @@ -974,4 +998,4 @@ def significance( obs_p_val, obs_significance, exp_p_val, exp_significance ) pyhf.set_backend(pyhf.tensorlib, initial_optimizer) # restore optimizer settings - return significance_results + return significance_results \ No newline at end of file From 6c8cd81e66822646f62d21c33b4882b4c9a34050 Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 14:13:11 +0200 Subject: [PATCH 3/7] add option to override default minuit scipy combo for ranking and limit setting --- src/cabinetry/fit/__init__.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index 8d651e1f..60a52c59 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -348,6 +348,7 @@ def ranking( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, + minimizer: Optional[Literal['minuit','scipy']] = None, strategy: Optional[Literal[0, 1, 2]] = None, solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, @@ -373,6 +374,10 @@ def ranking( parameters are held constant, defaults to None (use ``pyhf`` suggestion) par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to None (use minuit for the + initial fit to find the post fit uncertainties, then scipy for all fits used + to determine the impact). strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) @@ -397,7 +402,7 @@ def ranking( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, - minimizer='minuit', + minimizer='minuit' if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -453,7 +458,7 @@ def ranking( init_pars=init_pars_ranking, fix_pars=fix_pars_ranking, par_bounds=par_bounds, - minimizer='minuit', + minimizer='scipy' if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -496,6 +501,7 @@ def scan( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, + minimizer: Optional[Literal['minuit','scipy']] = None, strategy: Optional[Literal[0, 1, 2]] = None, solver_options: Optional[Dict] = None, maxiter: Optional[int] = None, @@ -521,6 +527,10 @@ def scan( parameters are held constant, defaults to None (use ``pyhf`` suggestion) par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to None (use minuit for the + initial fit to find the post fit uncertainties, then scipy for all fits used + for the likelihood scan). strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior of strategy 0 with user-provided gradients and 1 otherwise) @@ -553,7 +563,7 @@ def scan( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, - minimizer='minuit', + minimizer='minuit' if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -592,7 +602,7 @@ def scan( init_pars=init_pars_scan, fix_pars=fix_pars, par_bounds=par_bounds, - minimizer='scipy', + minimizer='scipy' if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, From 11596277d95316e055cd605b6c5f9343a01ccd6e Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 14:48:56 +0200 Subject: [PATCH 4/7] run precommit --- src/cabinetry/fit/__init__.py | 213 +++++++++++++++--------- src/cabinetry/fit/results_containers.py | 4 +- 2 files changed, 133 insertions(+), 84 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index 60a52c59..8229f0c5 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -49,9 +49,9 @@ def _fit_model( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, - minimizer: Optional[Literal['minuit','scipy']] = 'minuit', + minimizer: Optional[Literal["minuit", "scipy"]] = "minuit", strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict] = None, + solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> FitResults: @@ -84,8 +84,8 @@ def _fit_model( maxiter (Optional[int], optional): allowed number of calls for minimization, defaults to None (use ``pyhf`` default of 100,000) tolerance (Optional[float]), optional): tolerance for convergence, for details - see either ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance) if using - minuit, or ``scipy.optimize.minimize``. defaults to None (use ``iminuit`` + see either ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance) if using + minuit, or ``scipy.optimize.minimize``. defaults to None (use ``iminuit`` default of 0.1, or `scipy`` method specific defaults) Returns: FitResults: object storing relevant fit results @@ -93,17 +93,22 @@ def _fit_model( _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings # set up the optimiser - if minimizer == 'minuit': + optimizer_kwarg: Optional[Dict[str, Any]] = None + if minimizer == "minuit": pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(verbose=1)) # strategy=None is currently not supported in pyhf # https://github.com/scikit-hep/pyhf/issues/1785 optimizer_kwarg = {"strategy": strategy} if strategy is not None else {} - elif minimizer == 'scipy': + elif minimizer == "scipy": pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.scipy_optimizer(verbose=1)) - optimizer_kwarg = {"solver_options": solver_options} if solver_options is not None else {} + optimizer_kwarg = ( + {"solver_options": solver_options} if solver_options is not None else {} + ) else: - raise ValueError(f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}") - + raise ValueError( + f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}" + ) + result, corr_mat, best_twice_nll, result_obj = pyhf.infer.mle.fit( data, model, @@ -122,7 +127,7 @@ def _fit_model( labels = model.config.par_names best_twice_nll = float(best_twice_nll) # convert 0-dim np.ndarray to float - if minimizer == 'minuit': + if minimizer == "minuit": log.info(f"Migrad status:\n{result_obj.minuit.fmin}") bestfit = pyhf.tensorlib.to_numpy(result[:, 0]) # set errors for fixed parameters to 0 (see iminuit#762) @@ -133,12 +138,12 @@ def _fit_model( minos_results = ( _run_minos(result_obj.minuit, minos, labels) if minos is not None else {} ) - - elif minimizer == 'scipy': + + elif minimizer == "scipy": bestfit = pyhf.tensorlib.to_numpy(result) # scipy does not return uncertainty or correlation results uncertainty = np.zeros(bestfit.shape) - corr_mat = np.diag(no.ones(bestfit.shape)) + corr_mat = np.diag(np.ones(bestfit.shape)) minos_results = None fit_results = FitResults( @@ -150,13 +155,12 @@ def _fit_model( minos_uncertainty=minos_results, ) - log.debug(f"-2 log(L) = {fit_results.best_twice_nll:.6f} at best-fit point") - + log.debug(f"-2 log(L) = {fit_results.best_twice_nll:.6f} at best-fit point") + pyhf.set_backend(pyhf.tensorlib, initial_optimizer) # restore optimizer settings return fit_results - def _run_minos( minuit_obj: iminuit.Minuit, minos: Union[List[str], Tuple[str, ...]], @@ -268,45 +272,46 @@ def fit( model: pyhf.pdf.Model, data: List[float], *, - minos: Optional[Union[List[str], Tuple[str, ...]]] = None, + minos: Optional[Union[str, List[str], Tuple[str, ...]]] = None, + goodness_of_fit: bool = False, init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, - minimizer: Optional[Literal['minuit','scipy']] = 'minuit', + minimizer: Optional[Literal["minuit", "scipy"]] = "minuit", strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict] = None, + solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, - ) -> FitResults: +) -> FitResults: """Performs a maximum likelihood fit, reports and returns the results. - Args: - model (pyhf.pdf.Model): the model to use in the fit - data (List[float]): the data to fit the model to - minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS - algorithm for all parameters specified, defaults to None (does not run - MINOS). Requires ``minuit`` as the minimizer. - init_pars (Optional[List[float]], optional): list of initial parameter settings, - defaults to None (use ``pyhf`` suggested inits) - fix_pars (Optional[List[bool]], optional): list of booleans specifying which - parameters are held constant, defaults to None (use ``pyhf`` suggestion) - par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with - parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) - minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for - the fit. Can be 'minuit' or 'scipy', defaults to ``minuit``. - strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by - Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior - of strategy 0 with user-provided gradients and 1 otherwise) - solver_options (Optional[Dict], optional): solver options passed to scipy - optimizer. See ``scipy.optimize.show_options()`` for all available options. - maxiter (Optional[int], optional): allowed number of calls for minimization, - defaults to None (use ``pyhf`` default of 100,000) - tolerance (Optional[float]), optional): tolerance for convergence, for details - see either ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance) if using - minuit, or ``scipy.optimize.minimize``. defaults to None (use ``iminuit`` - default of 0.1, or `scipy`` method specific defaults) - Returns: - FitResults: object storing relevant fit results + Args: + model (pyhf.pdf.Model): the model to use in the fit + data (List[float]): the data to fit the model to + minos (Optional[Union[List[str], Tuple[str, ...]]], optional): runs the MINOS + algorithm for all parameters specified, defaults to None (does not run + MINOS). Requires ``minuit`` as the minimizer. + init_pars (Optional[List[float]], optional): list of initial parameter + settings, defaults to None (use ``pyhf`` suggested inits) + fix_pars (Optional[List[bool]], optional): list of booleans specifying which + parameters are held constant, defaults to None (use ``pyhf`` suggestion) + par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with + parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``minuit``. + strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by + Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior + of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to scipy + optimizer. See ``scipy.optimize.show_options()`` for all available options. + maxiter (Optional[int], optional): allowed number of calls for minimization, + defaults to None (use ``pyhf`` default of 100,000) + tolerance (Optional[float]), optional): tolerance for convergence, for details + see either ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance) if using + minuit, or ``scipy.optimize.minimize``. defaults to None (use ``iminuit`` + default of 0.1, or `scipy`` method specific defaults) + Returns: + FitResults: object storing relevant fit results """ log.info("performing maximum likelihood fit") @@ -348,9 +353,9 @@ def ranking( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, - minimizer: Optional[Literal['minuit','scipy']] = None, + minimizer: Optional[Literal["minuit", "scipy"]] = None, strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict] = None, + solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> RankingResults: @@ -375,7 +380,7 @@ def ranking( par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for - the fit. Can be 'minuit' or 'scipy', defaults to None (use minuit for the + the fit. Can be 'minuit' or 'scipy', defaults to None (use minuit for the initial fit to find the post fit uncertainties, then scipy for all fits used to determine the impact). strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by @@ -402,7 +407,7 @@ def ranking( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, - minimizer='minuit' if minimizer is None else minimizer, + minimizer="minuit" if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -458,7 +463,7 @@ def ranking( init_pars=init_pars_ranking, fix_pars=fix_pars_ranking, par_bounds=par_bounds, - minimizer='scipy' if minimizer is None else minimizer, + minimizer="scipy" if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -501,9 +506,9 @@ def scan( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, - minimizer: Optional[Literal['minuit','scipy']] = None, + minimizer: Optional[Literal["minuit", "scipy"]] = None, strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict] = None, + solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> ScanResults: @@ -528,7 +533,7 @@ def scan( par_bounds (Optional[List[Tuple[float, float]]], optional): list of tuples with parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for - the fit. Can be 'minuit' or 'scipy', defaults to None (use minuit for the + the fit. Can be 'minuit' or 'scipy', defaults to None (use minuit for the initial fit to find the post fit uncertainties, then scipy for all fits used for the likelihood scan). strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by @@ -563,7 +568,7 @@ def scan( init_pars=init_pars, fix_pars=fix_pars, par_bounds=par_bounds, - minimizer='minuit' if minimizer is None else minimizer, + minimizer="minuit" if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -602,7 +607,7 @@ def scan( init_pars=init_pars_scan, fix_pars=fix_pars, par_bounds=par_bounds, - minimizer='scipy' if minimizer is None else minimizer, + minimizer="scipy" if minimizer is None else minimizer, strategy=strategy, solver_options=solver_options, maxiter=maxiter, @@ -615,6 +620,55 @@ def scan( return scan_results +def _get_optimizer( + minimizer: Optional[Literal["minuit", "scipy"]] = "scipy", + strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict[str, Any]] = None, + maxiter: Optional[int] = None, + tolerance: Optional[float] = None, +) -> pyhf.optimize.mixins.OptimizerMixin: + """Creates an optimizer instance for either `scipy` or `minuit`. + + Args: + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``scipy``. + strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by + Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior + of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. + maxiter (Optional[int], optional): allowed number of calls for minimization, + defaults to None (use ``pyhf`` default of 100,000) + tolerance (Optional[float]), optional): tolerance for convergence, for details + see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to + None (use ``iminuit`` default of 0.1) + + Raises: + ValueError: if minimizer is not 'minuit' or 'scipy' + + Returns: + pyhf.optimize.mixins.OptimizerMixin: pyhf scipy optimizer object. + + """ + if minimizer == "minuit": + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance + ) + return optimizer + + if minimizer == "scipy": + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, + solver_options=solver_options, + maxiter=maxiter, + tolerance=tolerance, + ) + return optimizer + raise ValueError( + f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}" + ) + + def limit( model: pyhf.pdf.Model, data: List[float], @@ -627,9 +681,9 @@ def limit( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, - minimizer: Optional[Literal['minuit','scipy']] = 'scipy', + minimizer: Optional[Literal["minuit", "scipy"]] = "scipy", strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict] = None, + solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> LimitResults: @@ -685,17 +739,15 @@ def limit( LimitResults: observed and expected limits, CLs values, and scanned points """ _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings - - if minimizer == 'minuit': - optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance) - - elif minimizer == 'scipy': - optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, solver_options=solver_options, maxiter=maxiter, tolerance=tolerance) - else: - raise ValueError(f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}") + optimizer = _get_optimizer( + minimizer=minimizer, + strategy=strategy, + solver_options=solver_options, + maxiter=maxiter, + tolerance=tolerance, + ) + pyhf.set_backend(pyhf.tensorlib, optimizer) # use POI given by kwarg, fall back to POI specified in model @@ -913,9 +965,9 @@ def significance( init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, - minimizer: Optional[Literal['minuit','scipy']] = 'scipy', + minimizer: Optional[Literal["minuit", "scipy"]] = "scipy", strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict] = None, + solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, ) -> SignificanceResults: @@ -952,18 +1004,15 @@ def significance( SignificanceResults: observed and expected p-values and significances """ _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings - if minimizer == 'minuit': - optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance) - - elif minimizer == 'scipy': - optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, solver_options=solver_options, maxiter=maxiter, tolerance=tolerance) - else: - raise ValueError(f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}") + optimizer = _get_optimizer( + minimizer=minimizer, + strategy=strategy, + solver_options=solver_options, + maxiter=maxiter, + tolerance=tolerance, + ) pyhf.set_backend(pyhf.tensorlib, optimizer) - # use POI given by kwarg, fall back to POI specified in model poi_index = model_utils._poi_index(model, poi_name=poi_name) if poi_index is None: @@ -1008,4 +1057,4 @@ def significance( obs_p_val, obs_significance, exp_p_val, exp_significance ) pyhf.set_backend(pyhf.tensorlib, initial_optimizer) # restore optimizer settings - return significance_results \ No newline at end of file + return significance_results diff --git a/src/cabinetry/fit/results_containers.py b/src/cabinetry/fit/results_containers.py index f35f0a86..e7ccc61c 100644 --- a/src/cabinetry/fit/results_containers.py +++ b/src/cabinetry/fit/results_containers.py @@ -1,6 +1,6 @@ """Provides containers for inference results.""" -from typing import Dict, List, NamedTuple, Tuple +from typing import Dict, List, NamedTuple, Optional, Tuple import numpy as np @@ -26,7 +26,7 @@ class FitResults(NamedTuple): corr_mat: np.ndarray best_twice_nll: float goodness_of_fit: float = -1 - minos_uncertainty: Dict[str, Tuple[float, float]] = {} + minos_uncertainty: Optional[Dict[str, Tuple[float, float]]] = {} class RankingResults(NamedTuple): From 3657c886d5bb36c6d3505782bfb162d4f8e3d617 Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 14:52:21 +0200 Subject: [PATCH 5/7] move to _get_optimizer everywhere --- src/cabinetry/fit/__init__.py | 125 ++++++++++++++++------------------ 1 file changed, 58 insertions(+), 67 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index 8229f0c5..471e9d40 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -41,6 +41,55 @@ def print_results(fit_results: FitResults) -> None: ) +def _get_optimizer( + minimizer: Optional[Literal["minuit", "scipy"]] = "scipy", + strategy: Optional[Literal[0, 1, 2]] = None, + solver_options: Optional[Dict[str, Any]] = None, + maxiter: Optional[int] = None, + tolerance: Optional[float] = None, +) -> pyhf.optimize.mixins.OptimizerMixin: + """Creates an optimizer instance for either `scipy` or `minuit`. + + Args: + minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for + the fit. Can be 'minuit' or 'scipy', defaults to ``scipy``. + strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by + Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior + of strategy 0 with user-provided gradients and 1 otherwise) + solver_options (Optional[Dict], optional): solver options passed to ``scipy`` + optimizer. See ``scipy.optimize.show_options()`` for all available options. + maxiter (Optional[int], optional): allowed number of calls for minimization, + defaults to None (use ``pyhf`` default of 100,000) + tolerance (Optional[float]), optional): tolerance for convergence, for details + see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to + None (use ``iminuit`` default of 0.1) + + Raises: + ValueError: if minimizer is not 'minuit' or 'scipy' + + Returns: + pyhf.optimize.mixins.OptimizerMixin: pyhf scipy optimizer object. + + """ + if minimizer == "minuit": + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance + ) + return optimizer + + if minimizer == "scipy": + optimizer = pyhf.optimize.minuit_optimizer( + verbose=1, + solver_options=solver_options, + maxiter=maxiter, + tolerance=tolerance, + ) + return optimizer + raise ValueError( + f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}" + ) + + def _fit_model( model: pyhf.pdf.Model, data: List[float], @@ -93,21 +142,15 @@ def _fit_model( _, initial_optimizer = pyhf.get_backend() # store initial optimizer settings # set up the optimiser - optimizer_kwarg: Optional[Dict[str, Any]] = None - if minimizer == "minuit": - pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(verbose=1)) - # strategy=None is currently not supported in pyhf - # https://github.com/scikit-hep/pyhf/issues/1785 - optimizer_kwarg = {"strategy": strategy} if strategy is not None else {} - elif minimizer == "scipy": - pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.scipy_optimizer(verbose=1)) - optimizer_kwarg = ( - {"solver_options": solver_options} if solver_options is not None else {} - ) - else: - raise ValueError( - f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}" - ) + optimizer = _get_optimizer( + minimizer=minimizer, + strategy=strategy, + solver_options=solver_options, + maxiter=maxiter, + tolerance=tolerance, + ) + + pyhf.set_backend(pyhf.tensorlib, optimizer) result, corr_mat, best_twice_nll, result_obj = pyhf.infer.mle.fit( data, @@ -119,9 +162,6 @@ def _fit_model( return_correlations=True, return_fitted_val=True, return_result_obj=True, - maxiter=maxiter, - tolerance=tolerance, - **optimizer_kwarg, ) labels = model.config.par_names @@ -620,55 +660,6 @@ def scan( return scan_results -def _get_optimizer( - minimizer: Optional[Literal["minuit", "scipy"]] = "scipy", - strategy: Optional[Literal[0, 1, 2]] = None, - solver_options: Optional[Dict[str, Any]] = None, - maxiter: Optional[int] = None, - tolerance: Optional[float] = None, -) -> pyhf.optimize.mixins.OptimizerMixin: - """Creates an optimizer instance for either `scipy` or `minuit`. - - Args: - minimizer (Optional[Literal['minuit','scipy']]), optional): minimizer used for - the fit. Can be 'minuit' or 'scipy', defaults to ``scipy``. - strategy (Optional[Literal[0, 1, 2]], optional): minimization strategy used by - Minuit, can be 0/1/2, defaults to None (then uses ``pyhf`` default behavior - of strategy 0 with user-provided gradients and 1 otherwise) - solver_options (Optional[Dict], optional): solver options passed to ``scipy`` - optimizer. See ``scipy.optimize.show_options()`` for all available options. - maxiter (Optional[int], optional): allowed number of calls for minimization, - defaults to None (use ``pyhf`` default of 100,000) - tolerance (Optional[float]), optional): tolerance for convergence, for details - see ``iminuit.Minuit.tol`` (uses EDM < 0.002*tolerance), defaults to - None (use ``iminuit`` default of 0.1) - - Raises: - ValueError: if minimizer is not 'minuit' or 'scipy' - - Returns: - pyhf.optimize.mixins.OptimizerMixin: pyhf scipy optimizer object. - - """ - if minimizer == "minuit": - optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance - ) - return optimizer - - if minimizer == "scipy": - optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, - solver_options=solver_options, - maxiter=maxiter, - tolerance=tolerance, - ) - return optimizer - raise ValueError( - f"minimizer supports only 'minuit' and 'scipy', received : {minimizer}" - ) - - def limit( model: pyhf.pdf.Model, data: List[float], From a744c527011735881ba8bb6b043dbd7b8ef7148c Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 15:15:56 +0200 Subject: [PATCH 6/7] fix typing issue --- src/cabinetry/fit/__init__.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index 471e9d40..fabad371 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -47,7 +47,7 @@ def _get_optimizer( solver_options: Optional[Dict[str, Any]] = None, maxiter: Optional[int] = None, tolerance: Optional[float] = None, -) -> pyhf.optimize.mixins.OptimizerMixin: +) -> Union[pyhf.optimize.scipy_optimizer, pyhf.optimize.minuit_optimizer]: """Creates an optimizer instance for either `scipy` or `minuit`. Args: @@ -68,7 +68,8 @@ def _get_optimizer( ValueError: if minimizer is not 'minuit' or 'scipy' Returns: - pyhf.optimize.mixins.OptimizerMixin: pyhf scipy optimizer object. + Union[pyhf.optimize.scipy_optimizer, pyhf.optimize.minuit_optimizer]: pyhf + scipy or minuit optimizer object. """ if minimizer == "minuit": From 120ea045ca284d860ddec66b16b25192846b6c4f Mon Sep 17 00:00:00 2001 From: Marcel Date: Sat, 24 Jun 2023 19:28:06 +0200 Subject: [PATCH 7/7] fix typing issue --- src/cabinetry/fit/__init__.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index fabad371..1214a12b 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -74,15 +74,18 @@ def _get_optimizer( """ if minimizer == "minuit": optimizer = pyhf.optimize.minuit_optimizer( - verbose=1, strategy=strategy, maxiter=maxiter, tolerance=tolerance + verbose=1, + strategy={"strategy": strategy} if strategy is not None else {}, + maxiter=maxiter if maxiter is not None else 100000, + tolerance=tolerance, ) return optimizer if minimizer == "scipy": - optimizer = pyhf.optimize.minuit_optimizer( + optimizer = pyhf.optimize.scipy_optimizer( verbose=1, - solver_options=solver_options, - maxiter=maxiter, + solver_options=solver_options if solver_options is not None else {}, + maxiter=maxiter if maxiter is not None else 100000, tolerance=tolerance, ) return optimizer