pisa.analysis package

Submodules

pisa.analysis.analysis module

Common tools for performing an analysis collected into the class Analysis that can be subclassed by specific analyses.

class pisa.analysis.analysis.Analysis[source]

Bases: object

An analysis class that performs a flexible (global, local, nested, …) fit of a hypothesis to data.

Full analyses with functionality beyond just fitting, such as scans/profiles, are outside the scope of this class and should be defined by the user (for example by subclassing this class, see section on custom fitting methods below).

Every fit is run with the fit_recursively() method, where the fit strategy is defined by the three arguments method, method_kwargs and local_fit_kwargs (see documentation of fit_recursively() below for other arguments.) The method argument determines which sub-routine should be run, method_kwargs is a dictionary with any keyword arguments of that sub-routine, and local_fit_kwargs is a dictionary (or list thereof) defining any nested sub-routines that are run within the outer sub-routine. A sub-sub-routine defined in local_fit_kwargs should again be a dictionary with the three keywords method, method_kwargs and local_fit_kwargs. In this way, sub-routines can be arbitrarily stacked to define complex fit strategies.

The fit result is returned as a HypoFitResult instance. By default, this does not include the potentially large fit history (empty) or the binned metric maps, but they need to be requested from fit_recursively explicitly, as demonstrated in the examples below.

pprint

Whether to show live updates of minimizer progress (overridden by blindness).

Type:

bool, default: True

blindness

Whether to carry out a blind analysis. If True or 1, this hides actual parameter values from display and disallows these (as well as Jacobian, Hessian, etc.) from ending up in logfiles. If given an integer > 1, the fitted parameters are also prevented from being stored in fit info dictionaries.

Type:

bool or int, default: False

Examples

A canonical standard oscillation fit fits octants in theta23 separately and then runs a scipy minimizer to optimize locally in each octant. The arguments that would produce that result when passed to fit_recursively are:

method = "octants"
method_kwargs = {
    "angle": "theta23"
    "inflection_point": 45 * ureg.deg
}
local_fit_kwargs = {
    "method": "scipy",
    "method_kwargs": minimizer_settings,
    "local_fit_kwargs": None
}

Let’s say we also have a CP violating phase deltacp24 that we want to fit separately per quadrant split at 90 degrees. We want this done within each quadrant fit for theta23, making 4 fits in total. Then we would nest the quadrant fit for deltacp24 inside the octant fit like so:

method = "octants"
method_kwargs = {
    "angle": "theta23"
    "inflection_point": 45 * ureg.deg
}
local_fit_kwargs = {
    "method": "octants",
    "method_kwargs": {
        "angle": "deltacp24",
        "inflection_point": 90 * ureg.deg,
    }
    "local_fit_kwargs": {
        "method": "scipy",
        "method_kwargs": minimizer_settings,
        "local_fit_kwargs": None
    }
}

Let’s suppose we want to apply a grid-scan global fit method to sterile mixing parameters theta24 and deltam41, but we want to marginalize over all other parameters with a usual 3-flavor fit configuration. That could be achieved as follows:

method = "grid_scan"
method_kwargs = {
    "grid": {
        "theta24": np.geomspace(1, 20, 3) * ureg.deg,
        "deltam41": np.geomspace(0.01, 0.5, 4) * ureg["eV^2"],
    },
    "fix_grid_params": False,
}
local_fit_kwargs = {
    "method": "octants",
    "method_kwargs": {
        "angle": "theta23",
        "inflection_point": 45 * ureg.deg,
    }
    "local_fit_kwargs": {
        "method": "scipy",
        "method_kwargs": minimizer_settings,
        "local_fit_kwargs": None
    }
}

Instead of scipy, we can also use iminuit and nlopt for local minimization or global searches by writing a dictionary with "method": "iminuit" or "method": "nlopt", respectively.

NLOPT Options

NLOPT can be dropped in place of scipy and iminuit by writing a dictionary with "method": "nlopt" and choosing the algorithm by its name of the form NLOPT_{G,L}{N}_XXXX. PISA supports all of the derivative-free global (https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#global-optimization) and local (https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#local-derivative-free-optimization) algorithms. Algorithms requiring gradients such as BFGS are not supported. To use the Nelder-Mead algorithm, for example, the following settings could be used:

nlopt_settings = {
    "method": "nlopt",
    "method_kwargs": {
        "algorithm": "NLOPT_LN_NELDERMEAD",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
        # other options that can be set here:
        # xtol_abs, xtol_rel, stopval, maxeval, maxtime
        # after maxtime seconds, stop and return best result so far
        "maxtime": 60
    },
    "local_fit_kwargs": None  # no further nesting available
}

and then run a "chi2" fit with

best_fit_info = ana.fit_recursively(
    data_dist,
    dm,
    "chi2",
    None,
    store_fit_history=False,
    include_metric_maps=True,
    **nlopt_settings
)

, where no fit history will be kept in memory or be part of best_fit_info. Of course, you can also nest the nlopt_settings dictionary in any of the octants, ranges and so on by passing it as local_fit_kwargs.

Adding constraints

Adding inequality constraints to algorithms that support it is possible by writing a lambda function in a string that expects to get the current parameters as a ParamSet and returns a float. The result will ensure that the passed function stays positive (to be consistent with SciPy, cf. below). The string will be passed to eval() to build the callable function. For example, a silly way to bound delta_index > 0.1 would be:

"method_kwargs": {
    "algorithm": "NLOPT_LN_COBYLA",
    "ftol_abs": 1e-5,
    "ftol_rel": 1e-5,
    "maxtime": 30,
    "ineq_constraints": [
        # be sure to convert parameters to their magnitude
        "lambda params: params.delta_index.m - 0.1"
    ]
}

See also

Code example in test_constrained_minimization()

COBYLA fit using example pipeline with inequality constraint on theta23

Adding inequality constraints to algorithms that don’t support it can be done by either nesting the local fit in the constrained_fit method or to use NLOPT’s AUGLAG method that adds a penalty for constraint violations internally. For example, we could do this to fulfill the same constraint with the PRAXIS algorithm:

"method_kwargs": {
    "algorithm": "NLOPT_AUGLAG",
    "ineq_constraints":[
        "lambda params: params.delta_index.m - 0.1"
    ],
    "local_optimizer": {
        # supports all the same options as above
        "algorithm": "NLOPT_LN_PRAXIS",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
    }
}

Using global searches with local subsidiary minimizers

Some global searches, like evolutionary strategies, use local subsidiary minimizers. These can be defined just as above by passing a dictionary with the settings to the local_optimizer keyword. Note that, again, only gradient-free methods are supported. Here is an example for the “Multi-Level single linkage” (MLSL) algorithm, using PRAXIS as the local optimizer:

"method_kwargs": {
    "algorithm": "NLOPT_G_MLSL_LDS",
    "local_optimizer": {
        "algorithm": "NLOPT_LN_PRAXIS",
        "ftol_abs": 1e-5,
        "ftol_rel": 1e-5,
    }
}

For some evolutionary strategies such as ISRES, the population option can also be set.

"method_kwargs": {
    "algorithm": "NLOPT_GN_ISRES",
    "population": 100,
}

SciPy Options

PISA supports the local algorithms in SUPPORTED_LOCAL_SCIPY_MINIMIZERS. It also supports the global algorithms “differential_evolution”, “basinhopping”, “dual_annealing”, and “shgo”. Of these, all but “differential_evolution” allow configuring local subsidiary algorithms. Both inequality and equality constraints can be added to algorithms that support them.

Step by step, let us define a strategy consisting of Dual Annealing together with local SLSQP minimization, with both equality and inequality constraints for some fit parameters. First, set up the “global” algorithm only:

scipy_settings = {
    "method": "scipy",
    "method_kwargs": {
        "global_method": "dual_annealing",
        "options": {
            "maxiter": 10,
            "seed": 0,
            # other options that can be set here:
            # initial_temp, restart_temp_ratio, visit, accept, maxfun, no_local_search, rng
        },
    }
}

Second, load an SLSQP configuration from a file (included in PISA):

local_settings = from_file("settings/minimizer/slsqp_ftol1e-6_eps1e-4_maxiter1000.json")

Third, require theta23 to remain > 42 degrees via an inequality constraint and delta_index to stay fixed at -0.1 via an equality constraint, by inserting a list of constraint dicts into the SLSQP settings:

# constraint function can be callable or string
constrs_list = [
    {'type': 'eq',
     'fun': f'lambda params: params.delta_index.m_as("dimensionless") + 0.1'},
    {'type': 'ineq',
     'fun': lambda params: params.theta23.m_as("degree") - 42}
]
local_settings["options"]["value"]["constraints"] = constrs_list

Fourth, insert the SLSQP configuration with constraints into the global settings:

scipy_settings["local_fit_kwargs"] = local_settings

Finally, run a "chi2" fit with

best_fit_info = ana.fit_recursively(
    data_dist,
    dm,
    "chi2",
    None,
    store_fit_history=False,
    include_metric_maps=True,
    **scipy_settings
)

, where no fit history will be kept in memory or be part of best_fit_info.

Caution

Non-negligible constraint violations or stepping out of bounds may occur when global SciPy optimization is used in combination with constraints! This will hopefully get fixed soon.

See also

Code examples in test_constrained_minimization() and global_scipy_minimization()

Fits using example pipeline and different scipy solvers and constraint types

Custom fitting methods

Custom fitting methods are added by subclassing the analysis. The fit function name has to follow the scheme _fit_{method} where method is the name of the fit method. For instance, the function for scipy is called _fit_scipy and can be called by setting “method”: “scipy” in the fit strategy dict.

The function has to accept the parameters data_dist, hypo_maker, metric, external_priors_penalty, method_kwargs, and local_fit_kwargs. See fit_recursively() for descriptions of these arguments. The return value of the function must be a HypoFitResult. As an example, the following sub-class of Analysis has a custom fit method that, nonsensically, always sets 42 degrees as the starting value for theta23:

class SubclassedAnalysis(Analysis):

    def _fit_nonsense(
        self, data_dist, hypo_maker, metric,
        external_priors_penalty, method_kwargs, local_fit_kwargs,
        store_fit_history, include_metric_maps
    ):
        logging.info("Starting nonsense fit (setting theta23 to 42 deg)...")

        for pipeline in hypo_maker:
            if "theta23" in pipeline.params.free.names:
                pipeline.params.theta23.value = 42 * ureg["deg"]

        best_fit_info = self.fit_recursively(
            data_dist, hypo_maker, metric, external_priors_penalty,
            local_fit_kwargs["method"], local_fit_kwargs["method_kwargs"],
            local_fit_kwargs["local_fit_kwargs"],
            store_fit_history=store_fit_history,
            include_metric_maps=include_metric_maps
        )

        return best_fit_info

Now, the nonsense fit method can be used and nested with any other fit method like so:

ana = SubclassedAnalysis()
local_minuit = OrderedDict(
    method="iminuit",
    method_kwargs={
        "tol": 10,
    },
    local_fit_kwargs=None
)

local_nonsense_minuit = OrderedDict(
    method="nonsense",
    method_kwargs=None,
    local_fit_kwargs=local_minuit
)

fit_result = ana.fit_recursively(
    data_dist,
    distribution_maker,
    "chi2",
    None,
    store_fit_history=True,
    include_metric_maps=True,
    **local_nonsense_minuit
)
fit_recursively(data_dist, hypo_maker, metric, external_priors_penalty, method, method_kwargs=None, local_fit_kwargs=None, store_fit_history=False, include_metric_maps=False)[source]

Recursively apply global search strategies with local sub-fits.

Parameters:

  • data_dist (Sequence of MapSets or MapSet) – Data distribution to be fit. Can be an actual-, Asimov-, or pseudo-data distribution (where the latter two are derived from simulation and so aren’t technically “data”).

  • hypo_maker (Detectors or DistributionMaker) – Creates the per-bin expectation values per map based on its param values. Free params in the hypo_maker are modified by the minimizer to achieve a “best” fit.

  • metric (string or iterable of strings) – Metric by which to evaluate the fit. See documentation of Map.

  • external_priors_penalty (func) – User defined prior penalty function, which takes hypo_maker and metric as arguments and returns numerical value of penalty to the metric value. It is expected sign of the penalty is correctly specified inside the external_priors_penalty (e.g. negative for llh or positive for chi2).

  • method (str) – Name of the sub-routine to be run. Currently, the options are scipy, octants, best_of, grid_scan, constrained, ranges, condition, iminuit, and nlopt.

  • method_kwargs (dict) – Any keyword arguments taken by the sub-routine. May be None if the sub-routine takes no additional arguments.

  • local_fit_kwargs (dict or list thereof) – A dictionary defining subsidiary sub-routines with the keywords method, method_kwargs and local_fit_kwargs. May be None if the sub-routine is itself a local or global fit that runs no further subsidiary fits.

  • store_fit_history (bool (default: False)) – Whether to keep track of the fit history (sampled metric and parameter values) when minimizing, to include it in the result.

  • include_metric_maps (bool (default: False)) – Whether to include the binned metric contributions at the best fit in the result.

pisa.analysis.analysis.BasicAnalysis

Simple alias of Analysis for backwards compatibility. Note: both __name__ and __qualname__ of BasicAnalysis evaluate to “Analysis”.

class pisa.analysis.analysis.Counter(i=0)[source]

Bases: object

Simple counter object for use as a minimizer callback.

property count

Current count

Type:

int

reset()[source]

Reset counter

class pisa.analysis.analysis.HypoFitResult(metric=None, metric_val=None, data_dist=None, hypo_maker=None, minimizer_time=None, num_distributions_generated=None, minimizer_metadata=None, fit_history=None, other_metrics=None, include_detailed_metric_info=False, include_maps_binned=False)[source]

Bases: object

Holds all relevant information about a fit result.

__getitem__(i)[source]
static deserialize_detailed_metric_info(info_dict)[source]

Re-instantiate all PISA objects that used to be in the dictionary.

property detailed_metric_info
classmethod from_state(state)[source]

Initialize a HypoFitResult from a dictionary. Requires all entries in _state_attrs to be present as keys.

static get_detailed_metric_info(data_dist, hypo_maker, hypo_asimov_dist, params, metric, generalized_poisson_hypo=None, other_metrics=None, detector_name=None, include_maps_binned=False)[source]

Get detailed fit information, including e.g. maps that yielded the metric.

Parameters:

  • data_dist

  • hypo_maker

  • hypo_asimov_dist

  • params

  • metric

  • generalized_poisson_hypo

  • other_metrics

  • detector_name

  • include_maps_binned

Returns:

detailed_metric_info

Return type:

OrderedDict

property hypo_asimov_dist
property params
property serializable_state
property state
pisa.analysis.analysis.MINIMIZERS_ACCEPTING_CONSTRS = ('cobyla', 'slsqp', 'trust-constr', 'cobyqa')

Minimizers that allow constraints to be passed. All of these accept inequalities, and ‘slsqp’ and ‘trust-constr’ in addition accept equalities. As of SciPy 1.16.0, slsqp requires dictionaries, whereas trust-constr, cobyla 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.

pisa.analysis.analysis.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

pisa.analysis.analysis.SUPPORTED_LOCAL_SCIPY_MINIMIZERS = ('cobyla', 'l-bfgs-b', 'nelder-mead', 'slsqp', 'trust-constr')

Local SciPy minimizers that PISA currently supports via this interface.

pisa.analysis.analysis.global_scipy_minimization(pprint=False)[source]

Test global SciPy solvers with constraints.

pisa.analysis.analysis.test_analysis(pprint=False)[source]

Test recursive fit strategies on an Analysis sub-class.

pisa.analysis.analysis.test_constrained_minimization(pprint=False)[source]

Test SciPy solvers without or with equality and inequality constraints. All are run with default options as set by set_minimizer_defaults().

pisa.analysis.bayesian_analysis module

Common tools for performing Bayesian analysis.

pisa.analysis.bayesian_analysis.MCMC_sampling(data_dist, hypo_maker, *, metric, nwalkers, burnin, nsteps, pprint=True, return_burn_in=False, random_state=None, sampling_algorithm=None)[source]

Performs MCMC sampling. Only supports serial (single CPU) execution at the moment. See issue #830.

Parameters:

  • data_dist (Sequence of MapSets or MapSet) – Data distribution to be fit. Can be an actual-, Asimov-, or pseudo-data distribution (where the latter two are derived from simulation and so aren’t technically “data”).

  • hypo_maker (Detectors or DistributionMaker) – Creates the per-bin expectation values per map based on its param values. Free params in the hypo_maker are modified by the minimizer to achieve a “best” fit.

  • metric (string or iterable of strings) – Metric by which to evaluate the fit. See documentation of Map.

  • nwalkers (int) – Number of walkers

  • burnin (int) – Number of steps in burn in phase

  • nSteps (int) – Number of steps after burn in

  • pprint (bool (default: True)) – Whether to show live updates of the sampling progress.

  • return_burn_in (bool (default: False)) – Also return the steps of the burn in phase.

  • random_state (None or type accepted by utils.random_numbers.get_random_state (default: None)) – Random state of the walker starting points.

  • sampling_algorithm (None or emcee.moves object (default: None)) – Sampling algorithm used by the emcee sampler. None means to use the default which is a Goodman & Weare “stretch move” with parallelization. See https://emcee.readthedocs.io/en/stable/user/moves/#moves-user to learn more about the emcee sampling algorithms.

Returns:

  • scaled_chain (numpy array) – Array containing all points in the parameter space visited by each walker. It is sorted by steps, so all the first steps of all walkers come first. To for example get all values of the Nth parameter and the ith walker, use scaled_chain[i::nwalkers, N].

  • scaled_chain_burnin (numpy array (optional)) – Same as scaled_chain, but for the burn in phase.

pisa.analysis.configure_nlopt_minimization module

Utilities for validating and configuring nlopt minimization.

pisa.analysis.configure_nlopt_minimization.get_nlopt_inequality_constraint_funcs(method_kwargs, hypo_maker)[source]

Convert constraints expressions in terms of ParamSets into callables for nlopt. Evals expression(s) from method_kwargs and returns list of callables.

pisa.analysis.configure_scipy_minimization module

Utilities for validating and configuring scipy minimization settings.

pisa.analysis.configure_scipy_minimization.make_scipy_constraint_dict(constr_type, fun, jac=None, args=None)[source]

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

pisa.analysis.configure_scipy_minimization.make_scipy_local_minimizer_kwargs(minimizer_settings, constrs=None, bounds=None)[source]

Small helper function containing common logic for creating minimizer keyword args in calls to global routines via their scipy interface.

pisa.analysis.configure_scipy_minimization.scipy_constraints_to_callables(constr_dicts, hypo_maker)[source]

Convert constraints expressions in terms of ParamSets into callables for scipy. Overwrites “fun” entries in constr_dicts.

pisa.analysis.configure_scipy_minimization.set_minimizer_defaults(minimizer_settings)[source]

Fill in default values for minimizer settings.

Parameters:

minimizer_settings (dict)

Returns:

new_minimizer_settings

Return type:

dict

pisa.analysis.configure_scipy_minimization.validate_minimizer_settings(minimizer_settings)[source]

Validate minimizer settings.

Supported minimizers are the same as in set_minimizer_defaults.

See source for specific thresholds set.

Parameters:

minimizer_settings (dict)

Raises:

ValueError – If any minimizer settings are deemed to be invalid.

pisa.analysis.manipulate_params module

Utilities for manipulating parameters when fitting.

pisa.analysis.manipulate_params.get_separate_octant_params(hypo_maker, angle_name, inflection_point, tolerance=None)[source]

This function creates versions of the angle param that are confined to a single octant. It does this for both octant cases. This is used to allow fits to be done where only one of the octants is allowed. The fit can then be done for the two octant cases and compared to find the best fit.

Parameters:

  • hypo_maker (DistributionMaker or Detectors) – The hypothesis maker being used by the fitter

  • angle_name (string) – Name of the angle for which to create separate octant params.

  • inflection_point (quantity) – Point distinguishing between the two octants, e.g. 45 degrees

  • tolerance (quantity) – If the starting value is closer to the inflection point than the value of the tolerance, it is offset away from the inflection point by this amount.

Returns:

  • angle_orig (Param) – angle param as it was before applying the octant separation

  • angle_case1 (Param) – angle param confined to first octant

  • angle_case2 (Param) – angle param confined to second octant

pisa.analysis.manipulate_params.update_param_values(hypo_maker, params, update_nominal_values=False, update_range=False, update_is_fixed=False)[source]

Update just the values of parameters of a DistributionMaker without replacing the memory references inside.

This should be used in place of hypo_maker.update_params(params) unless one explicitly wants to replace the memory references to which the parameters in the DistributionMaker are pointing.

pisa.analysis.manipulate_params.update_param_values_detector(hypo_maker, params, update_nominal_values=False, update_range=False, update_is_fixed=False)[source]

Modification of the update_param_values function to use with the Detectors class.

Module contents