From 3afc6e9d100bef16d940e77137f84e15a5bfaa0b Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Mon, 20 Jul 2020 12:25:54 -0400 Subject: [PATCH 01/10] rewrite opt_scipy to use tensor stitching instead --- src/pyhf/__init__.py | 12 ++++++--- src/pyhf/optimize/opt_scipy.py | 48 ++++++++++++++++++++++++++-------- 2 files changed, 45 insertions(+), 15 deletions(-) diff --git a/src/pyhf/__init__.py b/src/pyhf/__init__.py index 8ff379360b..f25093b6f0 100644 --- a/src/pyhf/__init__.py +++ b/src/pyhf/__init__.py @@ -4,10 +4,8 @@ from .exceptions import InvalidBackend from . import events -tensorlib = tensor.numpy_backend() -default_backend = tensorlib -optimizer = optimize.scipy_optimizer() -default_optimizer = optimizer +tensorlib = None +optimizer = None def get_backend(): @@ -27,6 +25,12 @@ def get_backend(): return tensorlib, optimizer +tensorlib = tensor.numpy_backend() +default_backend = tensorlib +optimizer = optimize.scipy_optimizer() +default_optimizer = optimizer + + @events.register('change_backend') def set_backend(backend, custom_optimizer=None): """ diff --git a/src/pyhf/optimize/opt_scipy.py b/src/pyhf/optimize/opt_scipy.py index 408754f6d2..25fb8d4c45 100644 --- a/src/pyhf/optimize/opt_scipy.py +++ b/src/pyhf/optimize/opt_scipy.py @@ -1,4 +1,6 @@ """scipy.optimize-based Optimizer using finite differences.""" +from .. import get_backend, default_backend +from ..tensor.common import _TensorViewer from scipy.optimize import minimize import logging @@ -28,19 +30,37 @@ def minimize( Returns: bestfit parameters - + """ + tensorlib, _ = get_backend() + + all_idx = default_backend.astensor(range(pdf.config.npars), dtype='int') + all_init = default_backend.astensor(init_pars) + fixed_vals = fixed_vals or [] - indices = [i for i, _ in fixed_vals] - values = [v for _, v in fixed_vals] - constraints = [{'type': 'eq', 'fun': lambda v: v[indices] - values}] + fixed_values = [x[1] for x in fixed_vals] + fixed_idx = [x[0] for x in fixed_vals] + + variable_idx = [x for x in all_idx if x not in fixed_idx] + variable_init = all_init[variable_idx] + variable_bounds = [par_bounds[i] for i in variable_idx] + + tv = _TensorViewer([fixed_idx, variable_idx]) + + data = tensorlib.astensor(data) + fixed_values_tensor = tensorlib.astensor(fixed_values, dtype='float') + + def func(pars): + pars = tensorlib.astensor(pars) + constrained_pars = tv.stitch([fixed_values_tensor, pars]) + return objective(constrained_pars, data, pdf) + result = minimize( - objective, - init_pars, - constraints=constraints, + func, + variable_init, method='SLSQP', - args=(data, pdf), - bounds=par_bounds, + jac=False, + bounds=variable_bounds, options=dict(maxiter=self.maxiter), ) try: @@ -48,6 +68,12 @@ def minimize( except AssertionError: log.error(result) raise + + nonfixed_vals = result.x + fitted_val = result.fun + fitted_pars = tv.stitch( + [fixed_values_tensor, tensorlib.astensor(nonfixed_vals)] + ) if return_fitted_val: - return result.x, result.fun - return result.x + return fitted_pars, fitted_val + return fitted_pars From 288576473511c1661b7a80a58374401d08419a6d Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Mon, 20 Jul 2020 12:32:54 -0400 Subject: [PATCH 02/10] fix up all doctests associated with this change --- README.rst | 2 +- src/pyhf/infer/__init__.py | 12 ++++++------ src/pyhf/infer/mle.py | 2 +- src/pyhf/infer/test_statistics.py | 2 +- 4 files changed, 9 insertions(+), 9 deletions(-) diff --git a/README.rst b/README.rst index a59ac85467..634c0b5389 100644 --- a/README.rst +++ b/README.rst @@ -41,7 +41,7 @@ Hello World >>> pdf = pyhf.simplemodels.hepdata_like(signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0]) >>> CLs_obs, CLs_exp = pyhf.infer.hypotest(1.0, [51, 48] + pdf.config.auxdata, pdf, return_expected=True) >>> print('Observed: {}, Expected: {}'.format(CLs_obs, CLs_exp)) - Observed: [0.05290116], Expected: [0.06445521] + Observed: [0.05290051], Expected: [0.0644532] What does it support -------------------- diff --git a/src/pyhf/infer/__init__.py b/src/pyhf/infer/__init__.py index a63b48484b..6e282caa4c 100644 --- a/src/pyhf/infer/__init__.py +++ b/src/pyhf/infer/__init__.py @@ -24,13 +24,13 @@ def hypotest( ... test_poi, data, model, qtilde=True, return_expected_set=True ... ) >>> print(CLs_obs) - [0.05251554] + [0.05251497] >>> print(CLs_exp_band) - [[0.00260641] - [0.01382066] - [0.06445521] - [0.23526104] - [0.57304182]] + [[0.00260626] + [0.01382005] + [0.0644532 ] + [0.23525643] + [0.57303619]] Args: poi_test (Number or Tensor): The value of the parameter of interest (POI) diff --git a/src/pyhf/infer/mle.py b/src/pyhf/infer/mle.py index 7f1cb445af..c3903ef6a2 100644 --- a/src/pyhf/infer/mle.py +++ b/src/pyhf/infer/mle.py @@ -84,7 +84,7 @@ def fixed_poi_fit(poi_val, data, pdf, init_pars=None, par_bounds=None, **kwargs) ... test_poi, data, model, return_fitted_val=True ... ) >>> bestfit_pars - array([1. , 0.97224597, 0.87553894]) + array([1. , 0.97224594, 0.87553889]) >>> twice_nll array([28.92218013]) >>> -2 * model.logpdf(bestfit_pars, data) == twice_nll diff --git a/src/pyhf/infer/test_statistics.py b/src/pyhf/infer/test_statistics.py index 8b767dbe51..926b6777c6 100644 --- a/src/pyhf/infer/test_statistics.py +++ b/src/pyhf/infer/test_statistics.py @@ -31,7 +31,7 @@ def qmu(mu, data, pdf, init_pars, par_bounds): >>> init_pars = model.config.suggested_init() >>> par_bounds = model.config.suggested_bounds() >>> pyhf.infer.test_statistics.qmu(test_mu, data, model, init_pars, par_bounds) - 3.938244920380498 + 3.9382449203593524 Args: mu (Number or Tensor): The signal strength parameter From 4986a5dc192aaf02321b356b360f50dd00e6439c Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Mon, 20 Jul 2020 20:14:55 -0400 Subject: [PATCH 03/10] fix up validation tests --- tests/test_validation.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/test_validation.py b/tests/test_validation.py index 3a59823e57..a4f4a095f0 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -709,13 +709,13 @@ def validate_hypotest(pdf, data, mu_test, expected_result, tolerance=1e-6): @pytest.mark.parametrize( 'setup_and_tolerance', [ - (setup_1bin_shapesys(), 1e-6), - (setup_1bin_lumi(), 4e-6), - (setup_1bin_normsys(), 2e-9), - (setup_2bin_histosys(), 8e-5), - (setup_2bin_2channel(), 1e-6), - (setup_2bin_2channel_couplednorm(), 1e-6), - (setup_2bin_2channel_coupledhistosys(), 1e-6), + (setup_1bin_shapesys(), 4e-5), + (setup_1bin_lumi(), 6e-5), + (setup_1bin_normsys(), 2e-4), + (setup_2bin_histosys(), 9e-5), + (setup_2bin_2channel(), 3e-5), + (setup_2bin_2channel_couplednorm(), 3e-5), + (setup_2bin_2channel_coupledhistosys(), 4e-5), (setup_2bin_2channel_coupledshapefactor(), 2.5e-6), ], ids=[ From 707a6d4b7d9ba91b37400dead952da64ebe6ef91 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Fri, 17 Jul 2020 14:50:13 -0400 Subject: [PATCH 04/10] fix regression tests --- tests/test_regression.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_regression.py b/tests/test_regression.py index 7cca1822c4..e5f5785519 100644 --- a/tests/test_regression.py +++ b/tests/test_regression.py @@ -43,7 +43,7 @@ def test_sbottom_regionA_1300_205_60( CLs_obs, CLs_exp = calculate_CLs( sbottom_regionA_bkgonly_json, sbottom_regionA_1300_205_60_patch_json ) - assert CLs_obs == pytest.approx(0.24443627759085326, rel=1e-5) + assert CLs_obs == pytest.approx(0.24443627759085326, rel=1.2e-5) assert np.all( np.isclose( np.array(CLs_exp), @@ -56,7 +56,7 @@ def test_sbottom_regionA_1300_205_60( 0.8910420971601081, ] ), - rtol=1e-5, + rtol=1.2e-5, ) ) @@ -73,7 +73,7 @@ def test_sbottom_regionA_1400_950_60( CLs_obs, CLs_exp = calculate_CLs( sbottom_regionA_bkgonly_json, sbottom_regionA_1400_950_60_patch_json ) - assert CLs_obs == pytest.approx(0.021373283911064852, rel=1e-5) + assert CLs_obs == pytest.approx(0.021373283911064852, rel=1.1e-5) assert np.all( np.isclose( np.array(CLs_exp), @@ -86,7 +86,7 @@ def test_sbottom_regionA_1400_950_60( 0.5744843501873754, ] ), - rtol=1e-5, + rtol=7e-5, ) ) @@ -103,7 +103,7 @@ def test_sbottom_regionA_1500_850_60( CLs_obs, CLs_exp = calculate_CLs( sbottom_regionA_bkgonly_json, sbottom_regionA_1500_850_60_patch_json ) - assert CLs_obs == pytest.approx(0.04536774062150508, rel=1e-5) + assert CLs_obs == pytest.approx(0.04536774062150508, rel=1.3e-5) assert np.all( np.isclose( np.array(CLs_exp), @@ -116,7 +116,7 @@ def test_sbottom_regionA_1500_850_60( 0.6553686728646031, ] ), - rtol=1e-5, + rtol=7e-5, ) ) From 97cd430f59b3d9916f1d049d129c9f0dacb9d442 Mon Sep 17 00:00:00 2001 From: Matthew Feickert Date: Mon, 20 Jul 2020 21:42:20 -0500 Subject: [PATCH 05/10] Rename this doens't change the technical code, so don't bother the PR author with the request --- src/pyhf/optimize/opt_scipy.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pyhf/optimize/opt_scipy.py b/src/pyhf/optimize/opt_scipy.py index 25fb8d4c45..204f9ff73b 100644 --- a/src/pyhf/optimize/opt_scipy.py +++ b/src/pyhf/optimize/opt_scipy.py @@ -38,12 +38,12 @@ def minimize( all_init = default_backend.astensor(init_pars) fixed_vals = fixed_vals or [] - fixed_values = [x[1] for x in fixed_vals] fixed_idx = [x[0] for x in fixed_vals] + fixed_values = [x[1] for x in fixed_vals] - variable_idx = [x for x in all_idx if x not in fixed_idx] + variable_idx = [idx for idx in all_idx if idx not in fixed_idx] variable_init = all_init[variable_idx] - variable_bounds = [par_bounds[i] for i in variable_idx] + variable_bounds = [par_bounds[idx] for idx in variable_idx] tv = _TensorViewer([fixed_idx, variable_idx]) From cda70639826b20c894c4f1d745e1e0374de5ee3a Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Tue, 21 Jul 2020 09:25:51 -0400 Subject: [PATCH 06/10] add test against vanilla constraints for each optimizer --- tests/test_validation.py | 84 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 84 insertions(+) diff --git a/tests/test_validation.py b/tests/test_validation.py index a4f4a095f0..68bd77cd5b 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -5,6 +5,8 @@ import pytest from pathlib import Path import numpy as np +import scipy +import iminuit @pytest.fixture(scope='module') @@ -701,6 +703,9 @@ def validate_hypotest(pdf, data, mu_test, expected_result, tolerance=1e-6): return_expected_set=True, qtilde=False, ) + + breakpoint() + assert abs(CLs_obs - expected_result['obs']) / expected_result['obs'] < tolerance for result, expected in zip(CLs_exp_set, expected_result['exp']): assert abs(result - expected) / expected < tolerance @@ -923,3 +928,82 @@ def test_shapesys_nuisparfilter_validation(): assert np.allclose( reference_root_results['CLs_exp'], pyhf_results['CLs_exp'], atol=1e-4, rtol=1e-5 ) + + +def test_scipy_constraints(): + poi_val = 1.0 + pdf = pyhf.simplemodels.hepdata_like([5.0], [10.0], [3.5]) + data = [10.0] + pdf.config.auxdata + + result = pyhf.infer.mle.fixed_poi_fit(poi_val, data, pdf) + + init = pdf.config.suggested_init() + bounds = pdf.config.suggested_bounds() + fixed_vals = [(pdf.config.poi_index, poi_val)] + + indices = [i for i, _ in fixed_vals] + values = [v for _, v in fixed_vals] + constraints = [{'type': 'eq', 'fun': lambda v: v[indices] - values}] + + def objective(pars): + return pyhf.infer.mle.twice_nll(pars, data, pdf)[0] + + expected_result = scipy.optimize.minimize( + objective, + init, + constraints=constraints, + method='SLSQP', + args=(data, pdf), + bounds=bounds, + ) + + assert expected_result.success + assert np.allclose(result, expected_result.x) + + +def test_minuit_constraints(): + pyhf.set_backend('numpy', pyhf.optimize.minuit_optimizer()) + poi_val = 1.0 + pdf = pyhf.simplemodels.hepdata_like([5.0], [10.0], [3.5]) + data = [10.0] + pdf.config.auxdata + + result = pyhf.infer.mle.fixed_poi_fit(poi_val, data, pdf) + + init = pdf.config.suggested_init() + bounds = pdf.config.suggested_bounds() + fixed_vals = [(pdf.config.poi_index, poi_val)] + + indices = [i for i, _ in fixed_vals] + values = [v for _, v in fixed_vals] + constraints = [{'type': 'eq', 'fun': lambda v: v[indices] - values}] + + parnames = ['p{}'.format(i) for i in range(len(init))] + kw = {'limit_p{}'.format(i): b for i, b in enumerate(bounds)} + initvals = {'p{}'.format(i): v for i, v in enumerate(init)} + step_sizes = { + 'error_p{}'.format(i): (b[1] - b[0]) / float(1000) for i, b in enumerate(bounds) + } + + constraints = {} + for index, value in fixed_vals: + constraints['fix_p{}'.format(index)] = True + initvals['p{}'.format(index)] = value + + def objective(pars): + return pyhf.infer.mle.twice_nll(pars, data, pdf)[0] + + mm = iminuit.Minuit( + objective, + errordef=1, + use_array_call=True, + name=parnames, + **kw, + **constraints, + **initvals, + **step_sizes, + ) + + expected_result = mm.migrad() + + assert expected_result + assert np.allclose(result, mm.np_values()) From 4b77a759893e378a9402a2ece9962886885a293b Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Tue, 21 Jul 2020 09:37:00 -0400 Subject: [PATCH 07/10] remove breakpoint --- tests/test_validation.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_validation.py b/tests/test_validation.py index 68bd77cd5b..dd2e1f0107 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -704,8 +704,6 @@ def validate_hypotest(pdf, data, mu_test, expected_result, tolerance=1e-6): qtilde=False, ) - breakpoint() - assert abs(CLs_obs - expected_result['obs']) / expected_result['obs'] < tolerance for result, expected in zip(CLs_exp_set, expected_result['exp']): assert abs(result - expected) / expected < tolerance From 7995d0838ead80ce0ae933189ccf449beebccb5b Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Tue, 21 Jul 2020 09:48:21 -0400 Subject: [PATCH 08/10] fix it up --- tests/test_validation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_validation.py b/tests/test_validation.py index dd2e1f0107..7f27c80386 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -951,7 +951,6 @@ def objective(pars): init, constraints=constraints, method='SLSQP', - args=(data, pdf), bounds=bounds, ) From 63c82080c65d2059f20a68ae542a96921cef0519 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Tue, 21 Jul 2020 10:02:22 -0400 Subject: [PATCH 09/10] what the hell gha From a0db844360f65b00849ae941979e252522b41be6 Mon Sep 17 00:00:00 2001 From: Giordon Stark Date: Tue, 21 Jul 2020 10:04:16 -0400 Subject: [PATCH 10/10] black it up --- tests/test_validation.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/tests/test_validation.py b/tests/test_validation.py index 7f27c80386..86f8f5c512 100644 --- a/tests/test_validation.py +++ b/tests/test_validation.py @@ -947,11 +947,7 @@ def objective(pars): return pyhf.infer.mle.twice_nll(pars, data, pdf)[0] expected_result = scipy.optimize.minimize( - objective, - init, - constraints=constraints, - method='SLSQP', - bounds=bounds, + objective, init, constraints=constraints, method='SLSQP', bounds=bounds, ) assert expected_result.success