refactor: Use TensorViewer in SciPy optimizer, instead of explicit constraints

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
 --------------------

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):
     """

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)

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

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

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,26 +30,50 @@ 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_idx = [x[0] for x in fixed_vals]
+        fixed_values = [x[1] for x in fixed_vals]
+
+        variable_idx = [idx for idx in all_idx if idx not in fixed_idx]
+        variable_init = all_init[variable_idx]
+        variable_bounds = [par_bounds[idx] for idx 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:
             assert result.success
         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

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,
         )
     )
 

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,7 @@ def validate_hypotest(pdf, data, mu_test, expected_result, tolerance=1e-6):
         return_expected_set=True,
         qtilde=False,
     )
+
     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
@@ -709,13 +712,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=[
@@ -923,3 +926,77 @@ 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', 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())