fix: Sync SciPy minimizer initialized parameter values with fixed values (#1053)

* Fix bug in the SciPy optimizer where the init_pars is not synced up with the fixed values
   - Update the initial values to the fixed value for any fixed parameter
* Add tests for constraining values bounding initialization values

Author: kratsg

README.rst

@@ -43,7 +43,7 @@ Hello World
    >>> test_mu = 1.0
    >>> CLs_obs, CLs_exp = pyhf.infer.hypotest(test_mu, data, model, qtilde=True, return_expected=True)
    >>> print(f"Observed: {CLs_obs}, Expected: {CLs_exp}")
-   Observed: 0.052515541856109765, Expected: 0.06445521290832805
+   Observed: 0.05251497423736956, Expected: 0.06445320535890459
 
 What does it support
 --------------------

src/pyhf/infer/__init__.py

@@ -24,9 +24,9 @@ def hypotest(
         ...     test_poi, data, model, qtilde=True, return_expected_set=True
         ... )
         >>> CLs_obs
-        array(0.05251554)
+        array(0.05251497)
         >>> CLs_exp_band
-        [array(0.00260641), array(0.01382066), array(0.06445521), array(0.23526104), array(0.57304182)]
+        [array(0.00260626), array(0.01382005), array(0.06445321), array(0.23525644), array(0.57303621)]
 
     Args:
         poi_test (Number or Tensor): The value of the parameter of interest (POI)

src/pyhf/optimize/opt_scipy.py

@@ -60,6 +60,9 @@ def _minimize(
         values = [v for _, v in fixed_vals]
         if fixed_vals:
             constraints = [{'type': 'eq', 'fun': lambda v: v[indices] - values}]
+            # update the initial values to the fixed value for any fixed parameter
+            for idx, fixed_val in fixed_vals:
+                x0[idx] = fixed_val
         else:
             constraints = []
 

tests/test_optim.py

@@ -434,3 +434,23 @@ def test_stitch_pars(backend):
         0,
         60,
     ]
+
+
+def test_init_pars_sync_fixed_values_scipy(mocker):
+    opt = pyhf.optimize.scipy_optimizer()
+
+    minimizer = mocker.MagicMock()
+    opt._minimize(minimizer, None, [9, 9, 9], fixed_vals=[(0, 1)])
+    assert minimizer.call_args[0] == (None, [1, 9, 9])
+
+
+def test_init_pars_sync_fixed_values_minuit(mocker):
+    opt = pyhf.optimize.minuit_optimizer()
+
+    # patch all we need
+    from pyhf.optimize import opt_minuit
+
+    minimizer = mocker.patch.object(opt_minuit, 'iminuit')
+    opt._get_minimizer(None, [9, 9, 9], [(0, 10)] * 3, fixed_vals=[(0, 1)])
+    assert minimizer.Minuit.from_array_func.call_args[1]['start'] == [1, 9, 9]
+    assert minimizer.Minuit.from_array_func.call_args[1]['fix'] == [True, False, False]

tests/test_regression.py

@@ -79,12 +79,12 @@ def test_sbottom_regionA_1400_950_60(
             np.array(CLs_exp),
             np.array(
                 [
-                    0.002644707461012826,
-                    0.013976754489151644,
-                    0.06497313811425813,
-                    0.23644505123524753,
-                    0.5744843501873754,
-                ]
+                    0.0026445531093281147,
+                    0.013976126501170727,
+                    0.06497105816950004,
+                    0.23644030478043676,
+                    0.5744785776763938,
+                ],
             ),
             rtol=1e-5,
         )
@@ -103,17 +103,17 @@ 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.045367205665400624, rel=1e-5)
     assert np.all(
         np.isclose(
             np.array(CLs_exp),
             np.array(
                 [
-                    0.0059847029077065295,
-                    0.026103516126601122,
-                    0.10093985752614597,
-                    0.3101988586187604,
-                    0.6553686728646031,
+                    0.00598431785676406,
+                    0.026102240062850574,
+                    0.10093641492218848,
+                    0.31019245951964736,
+                    0.6553623337518385,
                 ]
             ),
             rtol=1e-5,

tests/test_validation.py

@@ -50,13 +50,13 @@ def expected_result_1bin_shapesys(mu=1.0):
     if mu == 1:
         expected_result = {
             "exp": [
-                0.06371799398864626,
-                0.15096503398048894,
-                0.3279606950533305,
-                0.6046087303039118,
-                0.8662627605298466,
+                0.06372011644331387,
+                0.1509686618126131,
+                0.3279657430196915,
+                0.604613569829645,
+                0.8662652332047568,
             ],
-            "obs": 0.4541865416107029,
+            "obs": 0.45418892944576333,
         }
     return expected_result
 
@@ -121,8 +121,14 @@ def spec_1bin_lumi():
 def expected_result_1bin_lumi(mu=1.0):
     if mu == 1:
         expected_result = {
-            "exp": [0.01060338, 0.04022273, 0.13614217, 0.37078321, 0.71104119],
-            "obs": 0.01047275,
+            "exp": [
+                0.01060400765567206,
+                0.04022451457730529,
+                0.13614632580079802,
+                0.37078985531427255,
+                0.7110468540175344,
+            ],
+            "obs": 0.010473144401519705,
         }
     return expected_result
 
@@ -189,13 +195,13 @@ def expected_result_1bin_normsys(mu=1.0):
     if mu == 1:
         expected_result = {
             "exp": [
-                7.47169462e-10,
-                5.7411551509088054e-08,
-                3.6898088062731205e-06,
-                0.00016965731538267896,
-                0.004392708998555453,
+                7.472581399417304e-10,
+                5.741738272450336e-08,
+                3.690120950161796e-06,
+                0.00016966882793076826,
+                0.004392935288879465,
             ],
-            "obs": 0.0006735317023683173,
+            "obs": 0.0006735336290569807,
         }
     return expected_result
 
@@ -262,7 +268,7 @@ def expected_result_2bin_histosys(mu=1):
     if mu == 1:
         expected_result = {
             "exp": [
-                7.134513306138892e-06,
+                7.133904244038431e-06,
                 0.00012547100627138575,
                 0.001880010666437615,
                 0.02078964907605385,
@@ -352,13 +358,13 @@ def expected_result_2bin_2channel(mu=1.0):
     if mu == 1:
         expected_result = {
             "exp": [
-                0.00043491354821983556,
-                0.0034223000502860606,
-                0.02337423265831151,
-                0.1218654225510158,
-                0.40382074249477845,
+                0.0004349234603527283,
+                0.003422361539161119,
+                0.02337454317608372,
+                0.12186650297311125,
+                0.40382274594391104,
             ],
-            "obs": 0.056332621064982304,
+            "obs": 0.0563327694384318,
         }
     return expected_result
 
@@ -447,16 +453,20 @@ def spec_2bin_2channel_couplednorm(source=source_2bin_2channel_couplednorm()):
 
 @pytest.fixture(scope='module')
 def expected_result_2bin_2channel_couplednorm(mu=1.0):
+    # NB: mac/linux differ for exp[0]
+    # Mac:   0.055222676184648795
+    # Linux: 0.05522273289103311
+    # Fill with midpoint of both values
     if mu == 1:
         expected_result = {
             "exp": [
-                0.055223914655538435,
-                0.13613239925395315,
-                0.3068720101493323,
-                0.5839470093910164,
-                0.8554725461337025,
+                0.05522270453784095,
+                0.1361301880753241,
+                0.30686879632329855,
+                0.5839437910061168,
+                0.8554708284963864,
             ],
-            "obs": 0.5906228034705155,
+            "obs": 0.5906216823766879,
         }
     return expected_result
 
@@ -569,11 +579,11 @@ def expected_result_2bin_2channel_coupledhistosys(mu=1.0):
     if mu == 1:
         expected_result = {
             "exp": [
-                1.7653746536962154e-05,
-                0.00026265644807799805,
-                0.00334003612780065,
-                0.031522353024659715,
-                0.17907742915143962,
+                1.7654378902209275e-05,
+                0.00026266409358853543,
+                0.0033401113778672156,
+                0.03152286332324451,
+                0.17907927340107824,
             ],
             "obs": 0.07967400132261188,
         }
@@ -703,7 +713,7 @@ def validate_hypotest(pdf, data, mu_test, expected_result, tolerance=1e-6):
     )
     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
+        assert abs(result - expected) / expected < tolerance, result
 
 
 @pytest.mark.parametrize(