Repository: lbl-anp/becquerel

Issue:
Rebinning uses first and last bin as overflow bins
When using `bq.rebin` or `bq.Spectrum.rebin`, the default behavior is to silently include all data outside the new binning scheme in the leftmost bin (if less than the leftmost bin edge) or in the rightmost bin (if greater than the rightmost bin edge). This behavior is unexpected and not documented. I think it should be controlled by a kwarg that defaults to disabling the behavior.

This behavior has been noted in #148 and #209, and #209 in particular notes the needed workaround: to add an extra bin at each end and then later remove them.

I just ran into a situation where I rebinned a spectrum but it still included the problematic data that I was trying to exclude, which was unhelpful and confusing.

Here is a MWE demonstrating the behavior:

```python
import numpy as np
import becquerel as bq

counts = 1000 * np.ones(10)
spec = bq.Spectrum(counts=counts)
cal = bq.Calibration("p[0] * x", [1.0])
spec.apply_calibration(cal)
print("Original spectrum:")
print("edges: ", spec.bin_edges_kev)
print("values:", spec.counts_vals)

for edges in [
    np.linspace(2, 9, num=8),
    np.linspace(2.5, 9.5, num=8),
]:
    for method in ["interpolation", "listmode"]:
        spec2 = spec.rebin(edges, method=method)
        print("")
        print(f"Rebinned spectrum ({method}):")
        print("edges: ", spec2.bin_edges_kev)
        print("values:", spec2.counts_vals)
```

The output is:

```
Original spectrum:
edges:  [ 0.  1.  2.  3.  4.  5.  6.  7.  8.  9. 10.]
values: [1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000. 1000.]

Rebinned spectrum (interpolation):
edges:  [2. 3. 4. 5. 6. 7. 8. 9.]
values: [3000. 1000. 1000. 1000. 1000. 1000. 2000.]

Rebinned spectrum (listmode):
edges:  [2. 3. 4. 5. 6. 7. 8. 9.]
values: [3000. 1000. 1000. 1000. 1000. 1000. 2000.]

Rebinned spectrum (interpolation):
edges:  [2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5]
values: [3500. 1000. 1000. 1000. 1000. 1000. 1500.]

Rebinned spectrum (listmode):
edges:  [2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5]
values: [3495. 1021. 1007.  961. 1005. 1005. 1506.]
```

Ideally all of the outputs of the rebinned spectra would be

```
[1000. 1000. 1000. 1000. 1000. 1000. 1000.]
```

or approximately those values. (Randomness will affect the final line.)

Discussion:

Patch:
diff --git a/becquerel/core/rebin.py b/becquerel/core/rebin.py
index ea39c99..2622c80 100644
--- a/becquerel/core/rebin.py
+++ b/becquerel/core/rebin.py
@@ -164,7 +164,7 @@ def _rebin_interpolation(
       in_spectrum: iterable of input spectrum counts in each bin
       in_edges: iterable of input bin edges (len = len(in_spectrum) + 1)
       out_edges_no_rightmost: iterable of output bin edges
-                              (sans the rightmost bin)
+                              (sans the rightmost bin edge)
       slopes: the slopes of each histogram bin, with the lines drawn between
               each bin edge. (len = len(in_spectrum))
       out_spectrum: for nb.guvectorize; This is the return array, do not
@@ -278,6 +278,7 @@ def rebin(
     method="interpolation",
     slopes=None,
     zero_pad_warnings=True,
+    include_overflows=False,
 ):
     """
     Spectra rebinning via deterministic or stochastic methods.
@@ -307,6 +308,9 @@ def rebin(
             [num_spectra, num_bins_in + 1] or [num_bins_in + 1]
         zero_pad_warnings (boolean): warn when edge overlap results in
             appending empty bins
+        include_overflows (boolean): whether to include counts below
+            out_edges[0] and above out_edges[-1] in the first and last
+            of the "out" bins, respectively, or to discard them (default).
 
     Raises:
         AssertionError: for bad input arguments
@@ -362,12 +366,35 @@ def rebin(
     _check_any_overlap(in_edges, out_edges)
     if zero_pad_warnings:
         _check_partial_overlap(in_edges, out_edges)
+
+    # add an extra bin at each end if we don't want overflow counts in
+    # the first and last bins
+    if not include_overflows:
+        if not np.isfinite(out_edges[0]):
+            raise RebinError(
+                "Lowest output edge must be finite if not including overflows"
+            )
+        if not np.isfinite(out_edges[-1]):
+            raise RebinError(
+                "Highest output edge must be finite if not including overflows"
+            )
+        out_edges_temp = np.concatenate(
+            (np.array([-np.inf]), out_edges, np.array([np.inf]))
+        )
+    else:
+        out_edges_temp = out_edges
+
     # Remove highest output edge, necessary for guvectorize
-    out_edges = out_edges[..., :-1]
+    out_edges = out_edges_temp[..., :-1]
 
     # Specific calls to (wrapped) nb.guvectorize'd rebinning methods
     if method == "interpolation":
-        return _rebin_interpolation(in_spectra, in_edges, out_edges, slopes)
+        hist = _rebin_interpolation(in_spectra, in_edges, out_edges, slopes)
     elif method == "listmode":
-        return _rebin_listmode(in_spectra, in_edges, out_edges)
-    raise ValueError(f"{method} is not a valid rebinning method")
+        hist = _rebin_listmode(in_spectra, in_edges, out_edges)
+    else:
+        raise ValueError(f"{method} is not a valid rebinning method")
+
+    if not include_overflows:
+        hist = hist[..., 1:-1]
+    return hist
diff --git a/becquerel/core/spectrum.py b/becquerel/core/spectrum.py
index 92178ff..76d0a88 100644
--- a/becquerel/core/spectrum.py
+++ b/becquerel/core/spectrum.py
@@ -1252,7 +1252,12 @@ class Spectrum:
         return obj
 
     def rebin(
-        self, out_edges, method="interpolation", slopes=None, zero_pad_warnings=True
+        self,
+        out_edges,
+        method="interpolation",
+        slopes=None,
+        zero_pad_warnings=True,
+        include_overflows=False,
     ):
         """
         Spectra rebinning via deterministic or stochastic methods.
@@ -1269,8 +1274,11 @@ class Spectrum:
                 for quadratic interpolation
                 (only applies for "interpolation" method)
                 [len(spectrum) + 1]
-        zero_pad_warnings (boolean): warn when edge overlap results in
-            appending empty bins
+            zero_pad_warnings (boolean): warn when edge overlap results in
+                appending empty bins
+            include_overflows (boolean): whether to include counts below
+                out_edges[0] and above out_edges[-1] in the first and last
+                of the "out" bins, respectively, or to discard them (default).
 
         Raises:
             SpectrumError: for bad input arguments
@@ -1297,6 +1305,7 @@ class Spectrum:
             method=method,
             slopes=slopes,
             zero_pad_warnings=zero_pad_warnings,
+            include_overflows=include_overflows,
         )
         return Spectrum(
             counts=out_spec,


Test patch:
diff --git a/tests/rebin_test.py b/tests/rebin_test.py
index 0f6ea16..67a1237 100644
--- a/tests/rebin_test.py
+++ b/tests/rebin_test.py
@@ -33,18 +33,35 @@ def old_edges(request):
         np.linspace(-1, 30, 31),
         np.linspace(0, 30, 7),
         np.linspace(0, 30, 45),
+        np.linspace(0, 25, 26),
+        np.linspace(5, 30, 26),
+        np.linspace(5, 25, 21),
     ],
     ids=[
         "1 keV bins",
         "1 keV bins with negative starting edge",
         "5 keV bins",
         ".66 keV bins",
+        "1 keV bins, left edge == 0, right edge < 30",
+        "1 keV bins, left edge > 0, right edge == 30",
+        "1 keV bins, left edge > 0, right edge < 30",
     ],
 )
 def new_edges(request):
     return request.param.astype(float)
 
 
+@pytest.fixture(
+    params=[True, False],
+    ids=[
+        "include overflow events outside of the new binning in the first and last bins",
+        "discard any events outside of the new binning",
+    ],
+)
+def include_overflows(request):
+    return request.param
+
+
 @pytest.fixture(params=[1, 5, 25], ids=["sparse counts", "~5 counts", "~25 counts"])
 def lam(request):
     """
@@ -72,13 +89,17 @@ def make_fake_spec_array(lam, size, dtype=float):
 class TestRebin:
     """Tests for core.rebin()"""
 
-    def test_total_counts(self, lam, old_edges, new_edges):
+    def test_total_counts(self, lam, old_edges, new_edges, include_overflows):
         """Check total counts in spectrum data before and after rebin
 
         interpolation and listmode are separate as the latter requires ints
         """
-        self._test_total_counts(lam, old_edges, new_edges, "interpolation", float)
-        self._test_total_counts(lam, old_edges, new_edges, "listmode", int)
+        self._test_total_counts(
+            lam, old_edges, new_edges, "interpolation", float, include_overflows
+        )
+        self._test_total_counts(
+            lam, old_edges, new_edges, "listmode", int, include_overflows
+        )
 
     def test_total_counts_listmode_float(self, lam, old_edges, new_edges):
         """Check that listmode rebinning raises a warning if counts are floats
@@ -94,22 +115,33 @@ class TestRebin:
                 zero_pad_warnings=False,
             )
 
-    def _test_total_counts(self, lam, old_edges, new_edges, method, dtype):
+    def _test_total_counts(
+        self, lam, old_edges, new_edges, method, dtype, include_overflows
+    ):
         """Check total counts in spectrum data before and after rebin"""
         old_counts, new_counts = self._gen_old_new_counts(
-            lam, old_edges, new_edges, method, dtype
+            lam, old_edges, new_edges, method, dtype, include_overflows
         )
-        assert np.isclose(old_counts.sum(), new_counts.sum())
+        if include_overflows:
+            assert np.isclose(old_counts.sum(), new_counts.sum())
+        else:
+            assert int(old_counts.sum()) >= int(new_counts.sum())
 
-    def test_rebin2d_counts(self, lam, old_edges, new_edges):
+    def test_rebin2d_counts(self, lam, old_edges, new_edges, include_overflows):
         """Check total counts in spectrum data before and after rebin
 
         interpolation and listmode are separate as the latter requires ints
         """
-        self._test_total_counts_2d(lam, old_edges, new_edges, "interpolation", "float")
-        self._test_total_counts_2d(lam, old_edges, new_edges, "listmode", "int")
+        self._test_total_counts_2d(
+            lam, old_edges, new_edges, "interpolation", "float", include_overflows
+        )
+        self._test_total_counts_2d(
+            lam, old_edges, new_edges, "listmode", "int", include_overflows
+        )
 
-    def _test_total_counts_2d(self, lam, old_edges, new_edges, method, dtype):
+    def _test_total_counts_2d(
+        self, lam, old_edges, new_edges, method, dtype, include_overflows
+    ):
         """Check total counts in spectra data before and after rebin"""
 
         nspectra = 20
@@ -123,38 +155,63 @@ class TestRebin:
             new_edges,
             method=method,
             zero_pad_warnings=False,
+            include_overflows=include_overflows,
         )
-        assert np.allclose(old_counts_2d.sum(axis=1), new_counts_2d.sum(axis=1))
+        if include_overflows:
+            assert np.allclose(old_counts_2d.sum(axis=1), new_counts_2d.sum(axis=1))
+        else:
+            assert np.all(
+                old_counts_2d.sum(axis=1).astype(int)
+                >= new_counts_2d.sum(axis=1).astype(int)
+            )
 
-    def test_unchanged_binning(self, lam, old_edges):
+    def test_unchanged_binning(self, lam, old_edges, include_overflows):
         """Check counts in spectrum is the same after a identity/null rebin
         (i.e. same bin edges for old and new)
 
         interpolation and listmode are separate as the latter requires ints
         """
-        self._test_unchanged_binning(lam, old_edges, "interpolation", "float")
-        self._test_unchanged_binning(lam, old_edges, "listmode", "int")
+        self._test_unchanged_binning(
+            lam, old_edges, "interpolation", "float", include_overflows
+        )
+        self._test_unchanged_binning(
+            lam, old_edges, "listmode", "int", include_overflows
+        )
 
-    def _test_unchanged_binning(self, lam, edges, method, dtype):
+    def _test_unchanged_binning(self, lam, edges, method, dtype, include_overflows):
         old_counts, new_counts = self._gen_old_new_counts(
-            lam, edges, edges, method, dtype
+            lam, edges, edges, method, dtype, include_overflows
         )
         assert np.allclose(old_counts, new_counts)
 
-    def _gen_old_new_counts(self, lam, old_edges, new_edges, method, dtype):
+    def _gen_old_new_counts(
+        self, lam, old_edges, new_edges, method, dtype, include_overflows
+    ):
         """generate old and new counts (1D)"""
         old_counts = make_fake_spec_array(lam, len(old_edges) - 1, dtype=dtype)
         new_counts = bq.core.rebin.rebin(
-            old_counts, old_edges, new_edges, method=method, zero_pad_warnings=False
+            old_counts,
+            old_edges,
+            new_edges,
+            method=method,
+            zero_pad_warnings=False,
+            include_overflows=include_overflows,
         )
         return (old_counts, new_counts)
 
-    def test_subset_bin_edges(self, lam, old_edges):
+    def test_subset_bin_edges(self, lam, old_edges, include_overflows):
         old_counts = make_fake_spec_array(lam, len(old_edges) - 1, dtype=float)
         new_counts = bq.core.rebin.rebin(
-            old_counts, old_edges, old_edges[1:-1], method="interpolation"
+            old_counts,
+            old_edges,
+            old_edges[1:-1],
+            method="interpolation",
+            include_overflows=include_overflows,
         )
-        assert np.isclose(np.sum(old_counts), np.sum(new_counts))
+        if include_overflows:
+            assert np.isclose(np.sum(old_counts), np.sum(new_counts))
+        else:
+            assert int(np.sum(old_counts)) >= int(np.sum(new_counts))
         assert np.allclose(old_counts[2:-2], new_counts[1:-1])
 
     def test_overlap_warnings(self, old_spec_float):
@@ -217,3 +274,56 @@ class TestRebin:
             label="rebinned",
         )
         plt.show()
+
+
+@pytest.mark.parametrize(
+    "edges2, counts2, method, include_overflows",
+    [
+        (
+            np.linspace(2, 9, num=8),
+            np.array([3000, 1000, 1000, 1000, 1000, 1000, 2000]),
+            "interpolation",
+            True,
+        ),
+        (
+            np.linspace(2, 9, num=8),
+            np.array([1000, 1000, 1000, 1000, 1000, 1000, 1000]),
+            "interpolation",
+            False,
+        ),
+        (
+            np.linspace(2, 9, num=8),
+            np.array([3000, 1000, 1000, 1000, 1000, 1000, 2000]),
+            "listmode",
+            True,
+        ),
+        (
+            np.linspace(2, 9, num=8),
+            np.array([1000, 1000, 1000, 1000, 1000, 1000, 1000]),
+            "listmode",
+            False,
+        ),
+        (
+            np.linspace(2.5, 9.5, num=8),
+            np.array([3500, 1000, 1000, 1000, 1000, 1000, 1500]),
+            "interpolation",
+            True,
+        ),
+        (
+            np.linspace(2.5, 9.5, num=8),
+            np.array([1000, 1000, 1000, 1000, 1000, 1000, 1000]),
+            "interpolation",
+            False,
+        ),
+    ],
+)
+def test_rebin_include_overflows(edges2, counts2, method, include_overflows):
+    """Perform specific numerical rebinning tests."""
+
+    counts = 1000 * np.ones(10)
+    spec = bq.Spectrum(counts=counts)
+    cal = bq.Calibration("p[0] * x", [1.0])
+    spec.apply_calibration(cal)
+
+    spec2 = spec.rebin(edges2, method=method, include_overflows=include_overflows)
+    assert np.allclose(counts2, spec2.counts_vals)
diff --git a/tests/spectrum_test.py b/tests/spectrum_test.py
index c4258d5..b725c83 100644
--- a/tests/spectrum_test.py
+++ b/tests/spectrum_test.py
@@ -1122,16 +1122,39 @@ def rebin_spectrum_success(request):
     return make_spec(request.param[0], lt=request.param[1])
 
 
-def test_spectrum_rebin_success(rebin_spectrum_success, rebin_new_edges, rebin_method):
+@pytest.fixture(
+    params=[True, False],
+    ids=[
+        "include overflow events outside of the new binning in the first and last bins",
+        "discard any events outside of the new binning",
+    ],
+)
+def include_overflows(request):
+    return request.param
+
+
+def test_spectrum_rebin_success(
+    rebin_spectrum_success, rebin_new_edges, rebin_method, include_overflows
+):
     kwargs = dict(
-        out_edges=rebin_new_edges, method=rebin_method, zero_pad_warnings=False
+        out_edges=rebin_new_edges,
+        method=rebin_method,
+        zero_pad_warnings=False,
+        include_overflows=include_overflows,
     )
     if (rebin_spectrum_success._counts is None) and (rebin_method == "listmode"):
         with pytest.warns(bq.SpectrumWarning):
             spec = rebin_spectrum_success.rebin(**kwargs)
     else:
         spec = rebin_spectrum_success.rebin(**kwargs)
-    assert np.isclose(rebin_spectrum_success.counts_vals.sum(), spec.counts_vals.sum())
+    if include_overflows:
+        assert np.isclose(
+            rebin_spectrum_success.counts_vals.sum(), spec.counts_vals.sum()
+        )
+    else:
+        assert int(rebin_spectrum_success.counts_vals.sum()) >= int(
+            spec.counts_vals.sum()
+        )
     if rebin_spectrum_success.livetime is None:
         assert spec.livetime is None
     else: