Merge pull request #102 from ManuelHu/hist-fix

Numeric fix for histograms with small binning

src/lgdo/types/histogram.py

@@ -75,8 +75,31 @@ def __init__(
 
         @classmethod
         def from_edges(
-            cls, edges: NDArray, binedge_attrs: dict[str, Any] | None = None
+            cls,
+            edges: NDArray | Iterable[float],
+            binedge_attrs: dict[str, Any] | None = None,
+        ) -> Histogram.Axis:
+            """Create a new axis with variable binning described by ``edges``."""
+            edges = np.array(edges)
+            return cls(edges, None, None, None, True, binedge_attrs)
+
+        @classmethod
+        def from_range_edges(
+            cls,
+            edges: NDArray | Iterable[float],
+            binedge_attrs: dict[str, Any] | None = None,
         ) -> Histogram.Axis:
+            """Create a new axis from the binning described by ``edges``, but try to convert it to
+            a evenly-spaced range object first.
+
+            .. warning ::
+
+                This function might return a wrong binning, especially in the case of very small
+                magnitudes of the spacing. See the documentation of :func:`numpy.isclose` for
+                details. Use this function only with caution, if you know the binning's order of
+                magniutude.
+            """
+            edges = np.array(edges)
             edge_diff = np.diff(edges)
             if np.any(~np.isclose(edge_diff, edge_diff[0])):
                 return cls(edges, None, None, None, True, binedge_attrs)
@@ -204,7 +227,14 @@ def __init__(
                 if not isinstance(ax, (hist.axis.Regular, hist.axis.Variable)):
                     msg = "only regular or variable axes of hist.Hist can be converted"
                     raise ValueError(msg)
-                b.append(Histogram.Axis.from_edges(ax.edges, binedge_attrs))
+                if isinstance(ax, hist.axis.Regular):
+                    step = (ax.edges[-1] - ax.edges[0]) / ax.size
+                    bax = Histogram.Axis(
+                        None, ax.edges[0], ax.edges[-1], step, True, binedge_attrs
+                    )
+                    b.append(bax)
+                else:
+                    b.append(Histogram.Axis.from_edges(ax.edges, binedge_attrs))
             b = self._create_binning(b)
         else:
             if binning is None:

tests/lh5/test_lh5_write.py

@@ -434,8 +434,8 @@ def test_write_histogram(caplog, tmptestdir):
         "my_group/my_histogram", f"{tmptestdir}/write_histogram_test.lh5"
     )
     assert np.array_equal(h3.weights.nda, np.array([[10, 10], [10, 10]]))
-    assert h3.binning[0].first == 2
-    assert h3.binning[1].last == 7
+    assert h3.binning[0].edges[0] == 2
+    assert h3.binning[1].edges[-1] == 7
     assert h3.isdensity
     assert h3.binning[0].get_binedgeattrs() == {"units": "ns"}
 
@@ -475,7 +475,9 @@ def test_write_histogram_variable(caplog, tmptestdir):
         np.array([[10, 10], [10, 10]]),
         (np.array([2, 3.5, 4]), np.array([5, 6.5, 7])),
         isdensity=True,
-    )  # Same field name, different values
+    )
+
+    # Same field name, different values
     store = lh5.LH5Store()
     store.write(
         h1,

tests/types/test_histogram.py

@@ -18,6 +18,10 @@ def test_init_hist_regular():
     h = Histogram(h, None)
 
     assert len(h.binning) == 2
+    assert len(h.binning[0].edges) == 11
+    assert h.binning[0].first == 0
+    assert h.binning[0].last == 1
+    assert h.binning[0].step == 0.1
 
     h = hist.Hist(hist.axis.Integer(start=0, stop=10))
     with pytest.raises(ValueError, match="only regular or variable axes"):
@@ -92,7 +96,10 @@ def test_datatype_name():
 def test_axes():
     h = Histogram(
         np.array([[1, 1], [1, 1]]),
-        (np.array([0, 1, 2]), np.array([2.1, 2.2, 2.3])),
+        (
+            Histogram.Axis.from_range_edges([0, 1, 2]),
+            Histogram.Axis.from_range_edges([2.1, 2.2, 2.3]),
+        ),
         attrs={"units": "m"},
     )
     assert len(h.binning) == 2
@@ -112,12 +119,25 @@ def test_axes():
     assert np.isclose(ax1.step, 0.1)
     assert isinstance(ax0.nbins, int)
 
+    h = Histogram(
+        np.array([[1, 1], [1, 1]]),
+        (np.array([0, 1, 2]), np.array([2.1, 2.2, 2.3])),
+        attrs={"units": "m"},
+    )
+    assert len(h.binning) == 2
+    str(h)
+    assert not h.binning[0].is_range
+    assert not h.binning[1].is_range
+
     h = Histogram(
         np.array([[1, 1], [1, 1]]),
         [(1, 3, 1), (4, 6, 1)],
     )
     ax0 = h.binning[0]
     str(h)
+    assert ax0.first == 1
+    assert ax0.last == 3
+    assert len(ax0.edges) == 3
 
     h = Histogram(
         np.array([[1, 1], [1, 1]]),
@@ -201,7 +221,7 @@ def test_ax_attributes():
 
 
 def test_view_as_hist():
-    h = Histogram(np.array([1, 1]), (np.array([0, 1, 2]),))
+    h = Histogram(np.array([1, 1]), (Histogram.Axis.from_range_edges([0, 1, 2]),))
     hi = h.view_as("hist")
     assert isinstance(hi.axes[0], hist.axis.Regular)
 
@@ -220,9 +240,14 @@ def test_view_as_hist():
     hi = h.view_as("hist")
     assert isinstance(hi.axes[0], hist.axis.Variable)
 
+    # no auto-conversion from variable to regular.
+    h = Histogram(np.array([1, 1]), (np.array([0, 1, 2]),))
+    hi = h.view_as("hist")
+    assert isinstance(hi.axes[0], hist.axis.Variable)
+
 
 def test_view_as_np():
-    h = Histogram(np.array([1, 1]), (np.array([0, 1, 2]),))
+    h = Histogram(np.array([1, 1]), (Histogram.Axis.from_range_edges([0, 1, 2]),))
     assert h.binning[0].is_range
     assert h.binning[0].nbins == 2
     nda, axes = h.view_as("np")