refactor: Make hypotest return CLs as 0-d tensor (#944)

* Have hypotest return the CLs as a 0-d tensor
   - Important for fully differential likelihoods
* Update the docs to reflect changes
* Update tests to use return type of 0-d tensor

README.rst

@@ -38,10 +38,12 @@ Hello World
 .. code:: python
 
    >>> import pyhf
-   >>> 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]
+   >>> model = pyhf.simplemodels.hepdata_like(signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0])
+   >>> data = [51, 48] + model.config.auxdata
+   >>> 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
 
 What does it support
 --------------------

docs/examples/notebooks/hello-world.ipynb

@@ -41,14 +41,20 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Observed: 0.05290116224852556, Expected: 0.06445521290832805\n"
+      "Observed: 0.052515541856109765, Expected: 0.06445521290832805\n"
      ]
     }
    ],
    "source": [
-    "pdf = pyhf.simplemodels.hepdata_like(signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0])\n",
-    "CLs_obs, CLs_exp = pyhf.infer.hypotest(1.0, [51, 48] + pdf.config.auxdata, pdf, return_expected=True)\n",
-    "print('Observed: {}, Expected: {}'.format(CLs_obs, CLs_exp))"
+    "model = pyhf.simplemodels.hepdata_like(\n",
+    "    signal_data=[12.0, 11.0], bkg_data=[50.0, 52.0], bkg_uncerts=[3.0, 7.0]\n",
+    ")\n",
+    "data = [51, 48] + model.config.auxdata\n",
+    "test_mu = 1.0\n",
+    "CLs_obs, CLs_exp = pyhf.infer.hypotest(\n",
+    "    test_mu, data, model, qtilde=True, return_expected=True\n",
+    ")\n",
+    "print(f\"Observed: {CLs_obs}, Expected: {CLs_exp}\")"
    ]
   },
   {
@@ -67,13 +73,15 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Observed CL_s: 0.05290116224852556, CL_sb: 0.023599998519978738, CL_b: 0.4461149342826869\n"
+      "Observed CL_s: 0.052515541856109765, CL_sb: 0.023324961200974572, CL_b: 0.44415349012077077\n"
      ]
     }
    ],
    "source": [
-    "CLs_obs, p_values = pyhf.infer.hypotest(1.0, [51, 48] + pdf.config.auxdata, pdf, return_tail_probs=True)\n",
-    "print('Observed CL_s: {}, CL_sb: {}, CL_b: {}'.format(CLs_obs, p_values[0], p_values[1]))"
+    "CLs_obs, p_values = pyhf.infer.hypotest(\n",
+    "    test_mu, data, model, qtilde=True, return_tail_probs=True\n",
+    ")\n",
+    "print(f\"Observed CL_s: {CLs_obs}, CL_sb: {p_values[0]}, CL_b: {p_values[1]}\")"
    ]
   },
   {
@@ -92,7 +100,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "assert CLs_obs == p_values[0]/p_values[1]"
+    "assert CLs_obs == p_values[0] / p_values[1]"
    ]
   },
   {
@@ -120,7 +128,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Observed CL_s: 0.05290116224852556\n",
+      "Observed CL_s: 0.052515541856109765\n",
       "\n",
       "Expected CL_s(-2 σ): 0.0026064088679947964\n",
       "Expected CL_s(-1 σ): 0.013820657528619273\n",
@@ -131,10 +139,16 @@
     }
    ],
    "source": [
-    "CLs_obs, CLs_exp_band = pyhf.infer.hypotest(1.0, [51, 48] + pdf.config.auxdata, pdf, return_expected_set=True)\n",
-    "print('Observed CL_s: {}\\n'.format(CLs_obs))\n",
-    "for p_value, n_sigma in enumerate(np.arange(-2,3)):\n",
-    "    print('Expected CL_s{}: {}'.format('      ' if n_sigma==0 else '({} σ)'.format(n_sigma),CLs_exp_band[p_value]))"
+    "CLs_obs, CLs_exp_band = pyhf.infer.hypotest(\n",
+    "    test_mu, data, model, qtilde=True, return_expected_set=True\n",
+    ")\n",
+    "print(f\"Observed CL_s: {CLs_obs}\\n\")\n",
+    "for p_value, n_sigma in enumerate(np.arange(-2, 3)):\n",
+    "    print(\n",
+    "        \"Expected CL_s{}: {}\".format(\n",
+    "            \"      \" if n_sigma == 0 else \"({} σ)\".format(n_sigma), CLs_exp_band[p_value],\n",
+    "        )\n",
+    "    )"
    ]
   }
  ],

src/pyhf/cli/infer.py

@@ -107,8 +107,8 @@ def cls(
         testpoi, ws.data(model), model, qtilde=is_qtilde, return_expected_set=True
     )
     result = {
-        'CLs_obs': tensorlib.tolist(result[0])[0],
-        'CLs_exp': tensorlib.tolist(tensorlib.reshape(result[-1], [-1])),
+        'CLs_obs': tensorlib.tolist(result[0]),
+        'CLs_exp': [tensorlib.tolist(tensor) for tensor in result[-1]],
     }
 
     if output_file is None:

src/pyhf/infer/__init__.py

@@ -23,14 +23,10 @@ def hypotest(
         >>> CLs_obs, CLs_exp_band = pyhf.infer.hypotest(
         ...     test_poi, data, model, qtilde=True, return_expected_set=True
         ... )
-        >>> print(CLs_obs)
-        [0.05251554]
-        >>> print(CLs_exp_band)
-        [[0.00260641]
-         [0.01382066]
-         [0.06445521]
-         [0.23526104]
-         [0.57304182]]
+        >>> CLs_obs
+        array(0.05251554)
+        >>> CLs_exp_band
+        [array(0.00260641), array(0.01382066), array(0.06445521), array(0.23526104), array(0.57304182)]
 
     Args:
         poi_test (Number or Tensor): The value of the parameter of interest (POI)
@@ -102,10 +98,11 @@ def hypotest(
     CLsb = sig_plus_bkg_distribution.pvalue(teststat)
     CLb = b_only_distribution.pvalue(teststat)
     CLs = CLsb / CLb
+    # Ensure that all CL values are 0-d tensors
     CLsb, CLb, CLs = (
-        tensorlib.reshape(CLsb, (1,)),
-        tensorlib.reshape(CLb, (1,)),
-        tensorlib.reshape(CLs, (1,)),
+        tensorlib.astensor(CLsb),
+        tensorlib.astensor(CLb),
+        tensorlib.astensor(CLs),
     )
 
     _returns = [CLs]
@@ -120,8 +117,7 @@ def hypotest(
             CLs = sig_plus_bkg_distribution.pvalue(
                 expected_bonly_teststat
             ) / b_only_distribution.pvalue(expected_bonly_teststat)
-            CLs_exp.append(tensorlib.reshape(CLs, (1,)))
-        CLs_exp = tensorlib.astensor(CLs_exp)
+            CLs_exp.append(tensorlib.astensor(CLs))
         if kwargs.get('return_expected'):
             _returns.append(CLs_exp[2])
         _returns.append(CLs_exp)
@@ -132,7 +128,7 @@ def hypotest(
         CLs = sig_plus_bkg_distribution.pvalue(
             expected_bonly_teststat
         ) / b_only_distribution.pvalue(expected_bonly_teststat)
-        _returns.append(tensorlib.reshape(CLs, (1,)))
+        _returns.append(tensorlib.astensor(CLs))
     # Enforce a consistent return type of the observed CLs
     return tuple(_returns) if len(_returns) > 1 else _returns[0]
 

tests/test_infer.py

@@ -59,7 +59,7 @@ def test_hypotest_default(tmpdir, hypotest_args):
     kwargs = {}
     result = pyhf.infer.hypotest(*hypotest_args, **kwargs)
     # CLs_obs
-    assert len(list(result)) == 1
+    assert pyhf.tensorlib.shape(result) == ()
     assert isinstance(result, type(tb.astensor(result)))
 
 
@@ -178,7 +178,7 @@ def logpdf(self, pars, data):
         1.0, model.expected_data(model.config.suggested_init()), model
     )
 
-    assert np.isclose(cls[0], 0.7267836451638846)
+    assert np.isclose(cls, 0.7267836451638846)
 
 
 @pytest.mark.parametrize("qtilde", [True, False])

tests/test_regression.py

@@ -28,7 +28,7 @@ def calculate_CLs(bkgonly_json, signal_patch_json):
     result = pyhf.infer.hypotest(
         1.0, workspace.data(model), model, qtilde=True, return_expected_set=True
     )
-    return result[0].tolist()[0], result[-1].ravel().tolist()
+    return result[0].tolist(), result[-1]
 
 
 def test_sbottom_regionA_1300_205_60(

tests/test_validation.py

@@ -906,16 +906,16 @@ def test_shapesys_nuisparfilter_validation():
         "observations": [{"data": [100, 10], "name": "channel1"}],
         "version": "1.0.0",
     }
-    w = pyhf.Workspace(spec)
-    m = w.model(
+    ws = pyhf.Workspace(spec)
+    model = ws.model(
         modifier_settings={
             'normsys': {'interpcode': 'code4'},
             'histosys': {'interpcode': 'code4p'},
         },
     )
-    d = w.data(m)
-    obs, exp = pyhf.infer.hypotest(1.0, d, m, return_expected_set=True)
-    pyhf_results = {'CLs_obs': obs[0], 'CLs_exp': [e[0] for e in exp]}
+    data = ws.data(model)
+    obs, exp = pyhf.infer.hypotest(1.0, data, model, return_expected_set=True)
+    pyhf_results = {'CLs_obs': obs, 'CLs_exp': [e for e in exp]}
 
     assert np.allclose(
         reference_root_results['CLs_obs'], pyhf_results['CLs_obs'], atol=1e-4, rtol=1e-5