diff --git a/setup.cfg b/setup.cfg index 4583a94a..a663ffeb 100644 --- a/setup.cfg +++ b/setup.cfg @@ -56,7 +56,7 @@ filterwarnings = ignore:no code associated:UserWarning:typeguard: [flake8] -max-complexity = 16 +max-complexity = 18 max-line-length = 88 exclude = docs/conf.py count = True diff --git a/src/cabinetry/fit/__init__.py b/src/cabinetry/fit/__init__.py index 8b36acbf..6a02d5ba 100644 --- a/src/cabinetry/fit/__init__.py +++ b/src/cabinetry/fit/__init__.py @@ -425,6 +425,7 @@ def ranking( data: List[float], *, fit_results: Optional[FitResults] = None, + poi_name: Optional[str] = None, init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, @@ -442,6 +443,8 @@ def ranking( data (List[float]): data (including auxdata) the model is fit to fit_results (Optional[FitResults], optional): nominal fit results to use for ranking, if not specified will repeat nominal fit, defaults to None + poi_name (Optional[str], optional): impact is calculated with respect to this + parameter, defaults to None (use POI specified in workspace) 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 @@ -451,6 +454,9 @@ def ranking( 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 + Returns: RankingResults: fit results for parameters, and pre- and post-fit impacts """ @@ -466,7 +472,13 @@ def ranking( labels = model.config.par_names() prefit_unc = model_utils.prefit_uncertainties(model) - nominal_poi = fit_results.bestfit[model.config.poi_index] + + # 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: + raise ValueError("no POI specified, cannot calculate ranking") + + nominal_poi = fit_results.bestfit[poi_index] # need to get values for parameter settings, as they will be partially changed # during the ranking (init/fix changes) @@ -476,9 +488,9 @@ def ranking( all_impacts = [] for i_par, label in enumerate(labels): - if label == model.config.poi_name: + if i_par == poi_index: continue # do not calculate impact of POI on itself - log.info(f"calculating impact of {label} on {labels[model.config.poi_index]}") + log.info(f"calculating impact of {label} on {labels[poi_index]}") # hold current parameter constant fix_pars_ranking = fix_pars.copy() @@ -509,7 +521,7 @@ def ranking( par_bounds=par_bounds, custom_fit=custom_fit, ) - poi_val = fit_results_ranking.bestfit[model.config.poi_index] + poi_val = fit_results_ranking.bestfit[poi_index] parameter_impact = poi_val - nominal_poi log.debug( f"POI is {poi_val:.6f}, difference to nominal is " @@ -526,9 +538,9 @@ def ranking( # remove parameter of interest from bestfit / uncertainty / labels # such that their entries match the entries of the impacts - bestfit = np.delete(fit_results.bestfit, model.config.poi_index) - uncertainty = np.delete(fit_results.uncertainty, model.config.poi_index) - labels = np.delete(fit_results.labels, model.config.poi_index).tolist() + bestfit = np.delete(fit_results.bestfit, poi_index) + uncertainty = np.delete(fit_results.uncertainty, poi_index) + labels = np.delete(fit_results.labels, poi_index).tolist() ranking_results = RankingResults( bestfit, uncertainty, labels, prefit_up, prefit_down, postfit_up, postfit_down @@ -582,7 +594,7 @@ def scan( # get index of parameter with name par_name par_index = model_utils._parameter_index(par_name, labels) - if par_index == -1: + if par_index is None: raise ValueError(f"parameter {par_name} not found in model") # run a fit with the parameter not held constant, to find the best-fit point @@ -644,6 +656,7 @@ def limit( tolerance: float = 0.01, maxiter: int = 100, confidence_level: float = 0.95, + poi_name: Optional[str] = None, init_pars: Optional[List[float]] = None, fix_pars: Optional[List[bool]] = None, par_bounds: Optional[List[Tuple[float, float]]] = None, @@ -652,7 +665,9 @@ def limit( Limits are calculated for the parameter of interest (POI) defined in the model. Brent's algorithm is used to automatically determine POI values to be tested. The - desired confidence level can be configured, and defaults to 95%. + desired confidence level can be configured, and defaults to 95%. In order to support + setting the POI directly without model recompilation, this temporarily changes the + POI index and name in the model configuration. Args: model (pyhf.pdf.Model): model to use in fits @@ -668,6 +683,8 @@ def limit( 100 confidence_level (float, optional): confidence level for calculation, defaults to 0.95 (95%) + poi_name (Optional[str], optional): limit is calculated for this parameter, + defaults to None (use POI specified in workspace) 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 @@ -676,13 +693,27 @@ def limit( parameter bounds for fit, defaults to None (use ``pyhf`` suggested bounds) Raises: + ValueError: if no POI is found ValueError: if lower and upper bracket value are the same + ValueError: if starting brackets do not enclose the limit Returns: LimitResults: observed and expected limits, CLs values, and scanned points """ pyhf.set_backend(pyhf.tensorlib, pyhf.optimize.minuit_optimizer(verbose=1)) + # 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: + raise ValueError("no POI specified, cannot calculate limit") + + # set POI index / name in model config to desired value, hypotest will pick this up + # save original values to reset model later + original_model_poi_index = model.config.poi_index + original_model_poi_name = model.config.poi_name + model.config.poi_index = poi_index + model.config.poi_name = model.config.par_names()[poi_index] + # show two decimals only if confidence level in percent is not an integer cl_label = ( f"{confidence_level:.{0 if (confidence_level * 100).is_integer() else 2}%}" @@ -708,12 +739,15 @@ def limit( if bracket is None: bracket = (bracket_left_default, bracket_right_default) elif bracket[0] == bracket[1]: + # set POI index / name in model back to their original values + model.config.poi_index = original_model_poi_index + model.config.poi_name = original_model_poi_name raise ValueError(f"the two bracket values must not be the same: {bracket}") cache_CLs: Dict[float, tuple] = {} # cache storing all relevant results def _cls_minus_threshold( - poi: float, + poi_val: float, model: pyhf.pdf.Model, data: List[float], cls_target: float, @@ -726,7 +760,7 @@ def _cls_minus_threshold( cache to avoid re-fitting known POI values and to store all relevant values. Args: - poi (float): value for parameter of interest + poi_val (float): value for parameter of interest model (pyhf.pdf.Model): model to use in fits data (List[float]): data (including auxdata) the model is fit to cls_target (float): target CLs value to find by varying POI @@ -738,19 +772,20 @@ def _cls_minus_threshold( Returns: float: absolute value of difference to CLs=``cls_target`` """ - if poi <= 0: + if poi_val <= 0: # no fit needed for negative POI value, return a default value log.debug( - f"skipping fit for {model.config.poi_name} = {poi:.4f}, setting CLs = 1" + f"skipping fit for {model.config.poi_name} = {poi_val:.4f}, setting " + "CLs = 1" ) return 1 - cls_target # distance of CLs = 1 to target CLs - cache = cache_CLs.get(poi) + cache = cache_CLs.get(poi_val) if cache: observed, expected = cache # use result from cache else: # calculate CLs results = pyhf.infer.hypotest( - poi, + poi_val, data, model, init_pars=init_pars, @@ -761,10 +796,10 @@ def _cls_minus_threshold( ) observed = float(results[0]) # 1 value per scan point expected = np.asarray(results[1]) # 5 per point (with 1 and 2 sigma bands) - cache_CLs.update({poi: (observed, expected)}) + cache_CLs.update({poi_val: (observed, expected)}) current_CLs = np.hstack((observed, expected))[which_limit] log.debug( - f"{model.config.poi_name} = {poi:.4f}, {limit_label} CLs = " + f"{model.config.poi_name} = {poi_val:.4f}, {limit_label} CLs = " f"{current_CLs:.4f}{' (cached)' if cache else ''}" ) return current_CLs - cls_target @@ -801,6 +836,9 @@ def _cls_minus_threshold( f"CLs values at {bracket[0]:.4f} and {bracket[1]:.4f} do not bracket " f"CLs={cls_target:.4f}, try a different starting bracket" ) + # set POI index / name in model back to their original values + model.config.poi_index = original_model_poi_index + model.config.poi_name = original_model_poi_name raise if not res.converged: @@ -842,6 +880,10 @@ def _cls_minus_threshold( bracket = (bracket_left, bracket_right) + # set POI index / name in model back to their original values + model.config.poi_index = original_model_poi_index + model.config.poi_name = original_model_poi_name + # report all results log.info(f"total of {steps_total} steps to calculate all limits") if not all_converged: diff --git a/src/cabinetry/model_utils.py b/src/cabinetry/model_utils.py index fb40bf3a..1ad15c3d 100644 --- a/src/cabinetry/model_utils.py +++ b/src/cabinetry/model_utils.py @@ -424,12 +424,14 @@ def unconstrained_parameter_count(model: pyhf.pdf.Model) -> int: return n_pars -def _parameter_index(par_name: str, labels: Union[List[str], Tuple[str, ...]]) -> int: +def _parameter_index( + par_name: str, labels: Union[List[str], Tuple[str, ...]] +) -> Optional[int]: """Returns the position of a parameter with a given name in the list of parameters. Useful together with ``pyhf.pdf._ModelConfig.par_names`` to find the position of a parameter when the name is known. If the parameter is not found, logs an error and - returns a default value of -1. + returns a default value of None. Args: par_name (str): name of parameter to find in list @@ -437,14 +439,48 @@ def _parameter_index(par_name: str, labels: Union[List[str], Tuple[str, ...]]) - names in the model Returns: - int: index of parameter + Optional[int]: index of parameter, or None if parameter was not found """ - par_index = next((i for i, label in enumerate(labels) if label == par_name), -1) - if par_index == -1: + par_index = next((i for i, label in enumerate(labels) if label == par_name), None) + if par_index is None: log.error(f"parameter {par_name} not found in model") return par_index +def _poi_index( + model: pyhf.pdf.Model, *, poi_name: Optional[str] = None +) -> Optional[int]: + """Returns the index of the POI specified in the argument or the model default. + + If a string is given as argument, this takes priority. Otherwise the POI from the + model is used. If no POI is specified there either, logs an error and returns a + default value of None. + + Args: + model (pyhf.pdf.Model): model for which to find the POI index + poi_name (Optional[str], optional): name of the POI, defaults to None + + Raises: + ValueError: if the specified POI name cannot be found in the model + + Returns: + Optional[int]: POI index, or None if no POI could be found + """ + if poi_name is not None: + # use POI given by kwarg if specified + poi_index = _parameter_index(poi_name, model.config.par_names()) + if poi_index is None: + raise ValueError(f"parameter {poi_name} not found in model") + elif model.config.poi_index is not None: + # use POI specified in model + poi_index = model.config.poi_index + else: + log.error("could not find POI for model") + poi_index = None + + return poi_index + + def _strip_auxdata(model: pyhf.pdf.Model, data: List[float]) -> List[float]: """Always returns observed yields, no matter whether data includes auxdata. diff --git a/src/cabinetry/schemas/config.json b/src/cabinetry/schemas/config.json index 67b7adf5..d841da9e 100644 --- a/src/cabinetry/schemas/config.json +++ b/src/cabinetry/schemas/config.json @@ -62,14 +62,14 @@ "$$target": "#/definitions/general", "description": "general settings", "type": "object", - "required": ["Measurement", "POI", "InputPath", "HistogramFolder"], + "required": ["Measurement", "InputPath", "HistogramFolder"], "properties": { "Measurement": { "description": "name of measurement", "type": "string" }, "POI": { - "description": "name of parameter of interest", + "description": "name of parameter of interest, defaults to empty string", "type": "string" }, "InputPath": { diff --git a/src/cabinetry/workspace.py b/src/cabinetry/workspace.py index 60f8d522..4f3d0a12 100644 --- a/src/cabinetry/workspace.py +++ b/src/cabinetry/workspace.py @@ -358,7 +358,8 @@ def measurements(self) -> List[Dict[str, Any]]: parameters = {"parameters": parameters_list} config_dict.update(parameters) - config_dict.update({"poi": self.config["General"]["POI"]}) + # POI defaults to "" (interpreted as "no POI" by pyhf) if not specified + config_dict.update({"poi": self.config["General"].get("POI", "")}) measurement.update({"config": config_dict}) measurements.append(measurement) return measurements diff --git a/tests/fit/test_fit.py b/tests/fit/test_fit.py index 8877fa99..8f2eaa2b 100644 --- a/tests/fit/test_fit.py +++ b/tests/fit/test_fit.py @@ -416,7 +416,7 @@ def test_ranking(mock_fit, example_spec): example_spec["measurements"][0]["config"]["parameters"][0]["fixed"] = False bestfit = np.asarray([1.0, 0.9]) uncertainty = np.asarray([0.1, 0.02]) - labels = ["mu", "staterror"] + labels = ["Signal strength", "staterror"] fit_results = fit.FitResults(bestfit, uncertainty, labels, np.empty(0), 0.0) model, data = model_utils.model_and_data(example_spec) ranking_results = fit.ranking(model, data, fit_results=fit_results) @@ -445,10 +445,17 @@ def test_ranking(mock_fit, example_spec): assert np.allclose(ranking_results.postfit_up, [0.2]) assert np.allclose(ranking_results.postfit_down, [-0.2]) - # fixed parameter in ranking, custom fit + # fixed parameter in ranking, custom fit, POI via kwarg example_spec["measurements"][0]["config"]["parameters"][0]["fixed"] = True + example_spec["measurements"][0]["config"]["poi"] = "" model, data = model_utils.model_and_data(example_spec) - ranking_results = fit.ranking(model, data, fit_results=fit_results, custom_fit=True) + ranking_results = fit.ranking( + model, + data, + fit_results=fit_results, + poi_name="Signal strength", + custom_fit=True, + ) # expect two calls in this ranking (and had 4 before, so 6 total): pre-fit # uncertainty is 0 since parameter is fixed, mock post-fit uncertainty is not 0 assert mock_fit.call_count == 6 @@ -462,6 +469,7 @@ def test_ranking(mock_fit, example_spec): ranking_results = fit.ranking( model, data, + poi_name="Signal strength", init_pars=[1.5, 1.0], fix_pars=[False, False], par_bounds=[(0, 5), (0.1, 10)], @@ -495,6 +503,10 @@ def test_ranking(mock_fit, example_spec): assert np.allclose(ranking_results.postfit_up, [0.3]) assert np.allclose(ranking_results.postfit_down, [-0.3]) + # no POI specified anywhere + with pytest.raises(ValueError, match="no POI specified, cannot calculate ranking"): + fit.ranking(model, data, fit_results=fit_results) + @mock.patch( "cabinetry.fit._fit_model", @@ -621,21 +633,32 @@ def test_limit(example_spec_with_background, caplog): ] caplog.clear() - # Asimov dataset with nominal signal strength of 0 + # Asimov dataset with nominal signal strength of 0, POI via kwarg example_spec_with_background["measurements"][0]["config"]["parameters"][0][ "inits" ] = [0.0] + example_spec_with_background["measurements"][0]["config"]["poi"] = "" model, data = model_utils.model_and_data(example_spec_with_background, asimov=True) - limit_results = fit.limit(model, data) + assert model.config.poi_index is None # no POI set before calculation + assert model.config.poi_name is None + limit_results = fit.limit(model, data, poi_name="Signal strength") assert np.allclose(limit_results.observed_limit, 0.586, rtol=2e-2) assert np.allclose(limit_results.expected_limit, expected_limit, rtol=2e-2) + assert model.config.poi_index is None # model config is preserved + assert model.config.poi_name is None caplog.clear() # bracket does not contain root, custom confidence level with pytest.raises( ValueError, match=re.escape("f(a) and f(b) must have different signs") ): - fit.limit(model, data, bracket=(1.0, 2.0), confidence_level=0.9) + fit.limit( + model, + data, + bracket=(1.0, 2.0), + confidence_level=0.9, + poi_name="Signal strength", + ) assert ( "CLs values at 1.0000 and 2.0000 do not bracket CLs=0.1000, try a different " "starting bracket" in [rec.message for rec in caplog.records] @@ -644,7 +667,7 @@ def test_limit(example_spec_with_background, caplog): # bracket with identical values with pytest.raises(ValueError, match="the two bracket values must not be the same"): - fit.limit(model, data, bracket=(3.0, 3.0)) + fit.limit(model, data, bracket=(3.0, 3.0), poi_name="Signal strength") caplog.clear() # init/fixed pars, par bounds, lower POI bound below 0 @@ -655,6 +678,7 @@ def test_limit(example_spec_with_background, caplog): fit.limit( model, data, + poi_name="Signal strength", init_pars=[0.9, 1.0], fix_pars=[False, True], par_bounds=[(-1, 5), (0.1, 10.0)], @@ -676,6 +700,13 @@ def test_limit(example_spec_with_background, caplog): ] caplog.clear() + # new model, error raised in previous step caused changes in poi_index / poi_name + model, data = model_utils.model_and_data(example_spec_with_background) + + # no POI specified anywhere + with pytest.raises(ValueError, match="no POI specified, cannot calculate limit"): + fit.limit(model, data) + def test_significance(example_spec_with_background): # increase observed data for smaller observed p-value diff --git a/tests/test_configuration.py b/tests/test_configuration.py index c5deb23d..98cb88a0 100644 --- a/tests/test_configuration.py +++ b/tests/test_configuration.py @@ -18,7 +18,6 @@ def test_validate(): config_valid = { "General": { "Measurement": "", - "POI": "", "HistogramFolder": "", "InputPath": "", }, @@ -32,7 +31,6 @@ def test_validate(): config_multiple_data_samples = { "General": { "Measurement": "", - "POI": "", "HistogramFolder": "", "InputPath": "", }, diff --git a/tests/test_model_utils.py b/tests/test_model_utils.py index 4fc1477f..219df065 100644 --- a/tests/test_model_utils.py +++ b/tests/test_model_utils.py @@ -4,6 +4,7 @@ import numpy as np import pyhf +import pytest from cabinetry import model_utils from cabinetry.fit.results_containers import FitResults @@ -332,11 +333,34 @@ def test__parameter_index(caplog): assert model_utils._parameter_index(par_name, tuple(labels)) == 1 caplog.clear() - assert model_utils._parameter_index("x", labels) == -1 + assert model_utils._parameter_index("x", labels) is None assert "parameter x not found in model" in [rec.message for rec in caplog.records] caplog.clear() +def test__poi_index(example_spec, caplog): + caplog.set_level(logging.DEBUG) + model = pyhf.Workspace(example_spec).model() + + # POI given by name + assert model_utils._poi_index(model, poi_name="staterror_Signal-Region[0]") == 1 + + # parameter not found in model + with pytest.raises(ValueError, match="parameter x not found in model"): + model_utils._poi_index(model, poi_name="x") + + # POI specified in model config + assert model_utils._poi_index(model) == 0 + caplog.clear() + + # no POI + example_spec["measurements"][0]["config"]["poi"] = "" + model = pyhf.Workspace(example_spec).model() + assert model_utils._poi_index(model) is None + assert "could not find POI for model" in [rec.message for rec in caplog.records] + caplog.clear() + + def test__strip_auxdata(example_spec): model = pyhf.Workspace(example_spec).model() data_with_aux = list(model.expected_data([1.0, 1.0], include_auxdata=True)) diff --git a/tests/test_workspace.py b/tests/test_workspace.py index a30c4ca1..27c5cce6 100644 --- a/tests/test_workspace.py +++ b/tests/test_workspace.py @@ -315,9 +315,9 @@ def test_WorkspaceBuilder_channels(mock_contains, mock_histogram): assert mock_histogram.call_count == 1 -def test_WorkspaceBuilder_get_measurement(): +def test_WorkspaceBuilder_measurements(): example_config = { - "General": {"Measurement": "fit", "POI": "mu", "HistogramFolder": "path"}, + "General": {"Measurement": "fit", "HistogramFolder": "path"}, "NormFactors": [ {"Name": "NF", "Nominal": 1.0, "Bounds": [0.0, 5.0], "Fixed": False} ], @@ -326,7 +326,7 @@ def test_WorkspaceBuilder_get_measurement(): { "name": "fit", "config": { - "poi": "mu", + "poi": "", "parameters": [{"name": "NF", "inits": [1.0], "bounds": [[0.0, 5.0]]}], }, } @@ -334,7 +334,7 @@ def test_WorkspaceBuilder_get_measurement(): ws_builder = workspace.WorkspaceBuilder(example_config) assert ws_builder.measurements() == expected_measurement - # no norm factor settings + # no norm factor settings, POI specified example_config_no_NF_settings = { "General": {"Measurement": "fit", "POI": "mu", "HistogramFolder": "path"}, "NormFactors": [{"Name": "NF"}], @@ -372,7 +372,6 @@ def test_WorkspaceBuilder_get_measurement(): example_config_const_sys = { "General": { "Measurement": "fit", - "POI": "mu", "Fixed": [{"Name": "par_A", "Value": 1.2}], "HistogramFolder": "path", }, @@ -382,7 +381,7 @@ def test_WorkspaceBuilder_get_measurement(): { "name": "fit", "config": { - "poi": "mu", + "poi": "", "parameters": [{"name": "par_A", "fixed": True, "inits": [1.2]}], }, } @@ -397,11 +396,11 @@ def test_WorkspaceBuilder_get_measurement(): return_value=None, ) as mock_find_const: example_config_sys = { - "General": {"Measurement": "fit", "POI": "mu", "HistogramFolder": "path"}, + "General": {"Measurement": "fit", "HistogramFolder": "path"}, "Systematics": [{"Name": "par_A"}], } expected_measurement_sys = [ - {"name": "fit", "config": {"poi": "mu", "parameters": []}} + {"name": "fit", "config": {"poi": "", "parameters": []}} ] ws_builder = workspace.WorkspaceBuilder(example_config_sys) assert ws_builder.measurements() == expected_measurement_sys diff --git a/tests/visualize/test_visualize_plot_model.py b/tests/visualize/test_visualize_plot_model.py index 3a8b778d..881f8090 100644 --- a/tests/visualize/test_visualize_plot_model.py +++ b/tests/visualize/test_visualize_plot_model.py @@ -108,7 +108,7 @@ def test_data_mc(tmp_path, caplog): histo_dict_list[1]["yields"] = np.asarray([0, 5]) caplog.clear() with mock.patch("cabinetry.visualize.utils._save_and_close") as mock_close_safe: - with pytest.warns(None) as warn_record: + with pytest.warns() as warn_record: fig = fig = plot_model.data_mc( histo_dict_list, total_model_unc_log,