Repository: scikit-hep/hist

Issue:
[FEATURE] profile histograms
I personally find myself constructing profiles (`TH2::ProfileX()/TH2::ProfileY()`) during exploratory analysis. Is this something best left to the user, or would `hist` benefit from having the ability to do something like `h.profile(0).plot()`? To consider a specific case of 2 numerical axes,

```python
import numpy as np

def calculate_profile(self, axis=0):
    # matches the default error option for `TH2::ProfileX/Y()`
    values = self.values()
    variances = self.variances()
    centers = self.axes[axis].centers
    if axis == 0:
        num = np.matmul(values, centers)
        num_err = np.sqrt(np.matmul(variances, centers ** 2))
    else:
        num = np.matmul(values.T, centers)
        num_err = np.sqrt(np.matmul(variances.T, centers ** 2))
    den = np.sum(values.T, axis=axis)
    den_err = np.sqrt(np.sum(variances.T, axis=axis))
    new_values = num / den
    new_variances = ((num_err / den) ** 2 - (den_err * num / den ** 2.0) ** 2.0)
    return new_values, new_variances
```
(apologies for the convoluted numpy operations. I'm sure it can be simplified with some more thought)
and
```python
import hist

np.random.seed(42)
xy = np.random.multivariate_normal([0, 1], cov=[[1.,0.8],[0.8,1]], size=1000)
h = hist.Hist(hist.axis.Regular(50, -3, 3),hist.axis.Regular(50, -3, 3)).fill(*xy.T)

fig, ax = plt.subplots()
h.plot()

profile_axis = 0
values, variances = calculate_profile(h, axis=profile_axis)
ax.errorbar(h.axes[profile_axis].centers, values, yerr=variances**0.5, marker=".", color="r", alpha=0.5)
```
would yield
![image](https://user-images.githubusercontent.com/5760027/111580825-72178700-8775-11eb-9b2f-92620bf3e38d.png)




Discussion:
aminnj: After fixing the sign above, I get consistency between the function I have in https://github.com/scikit-hep/hist/issues/156 as well as ROOT (`TH2F::ProfileX()`), but the error bars from `h.profile(1).plot(label="hist")` look to be too small:
![image](https://user-images.githubusercontent.com/5760027/111738247-3a284680-883e-11eb-832b-8c52db43b7a8.png).
I'm not sure what the issue could be. I verified that the calculations up to, and including, `new_variances` in
 https://github.com/scikit-hep/hist/pull/160/files/3a290a9f8fe2a7023686a74c80e46c8fd1d9929e#diff-dd1c19dfb377c8c01bb8892201933a2101eef8d01c705cb55ab3c93530cfec78R313-R317 
are correct, but the return histogram from `h.profile(1)` has variances that are too small.

aminnj: Here's a (currently failing) test case.
```python
import hist
import numpy as np

xy = np.array([
    [-2, 1.5], [-2, -3.5], [-2, 1.5], # x = -2
    [0., -2.], [0., -2.], [0., 0.], [0., 2.], [0., 4.], # x = 0
    [2, 1.5], # x = +2
])
h = hist.Hist(hist.axis.Regular(5, -5, 5),hist.axis.Regular(5, -5, 5)).fill(*xy.T)

# Profile out the y-axis
hp = h.profile(1)

# Exclude edge bins since no values from above will fall into them
# When there are values there, ROOT does something funky in those bins,
# despite these bins not being in the axis that is profiled out, and
# despite there being no overflow... to be understood.
print("hist values:   ", hp.values()[1:-1])
print("hist variances:", hp.variances()[1:-1])
# hist values:    [-1.48029737e-16  4.00000000e-01  2.00000000e+00]
# hist variances: [0.88888889 0.2176            nan]


print("root values:   ", np.array([0.0, 0.4, 2.0]))
print("root variances:", np.array([2.66666667, 1.088, 0.0]))
# root values:    [0.  0.4 2. ]
# root variances: [2.66666667 1.088      0.        ]
```

The ROOT part is
```python
import ROOT
h = ROOT.TH2F("h","", 5,-5,5, 5,-5,5)
for x,y in xy:
    h.Fill(x,y)
hp = h.ProfileX()
for i in [2,3,4]:
    print(hp.GetBinContent(i), hp.GetBinError(i)**2)
```
henryiii: I added a property-based test, which kept failing, and finally highlighted the mistake - my mistake, will fix. (I'll probably put the pbt in boost-histogram, rather than here).
henryiii: Now that should be closer?
aminnj: Indeed, errors look perfect now.
![image](https://user-images.githubusercontent.com/5760027/111854530-e456ab80-88dc-11eb-930c-ee3ed95751d5.png)

For the test case above, I get
```
hist variances: [2.66666667 1.088             nan]
root variances: [2.66666667 1.088      0.        ]
```
I guess the `0` vs `nan` is a convention. There's only 1 value from which to compute a std in that slice, so if nan suppresses the error bar when plotting, then I don't see a practical difference between the two.
ast0815: I have just released v2.0.1 of Histoprint. It should now work with `nan` in the values. They are treated as `0`.
sourcery-ai[bot]: # Sourcery Code Quality Report


✅ &nbsp;Merging this PR will **increase** code quality in the affected files by **0.14%**.

| Quality metrics | Before | After | Change |
| --- | :---: | :---: | :---: |
| [Complexity](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#complexity) | 9.80 🙂 | 9.62 🙂 | -0.18 👍 |
| [Method Length](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#method-length) | 164.33 😞 | 164.16 😞 | -0.17 👍 |
| [Working memory](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#working-memory) | 15.99 ⛔ | 15.94 ⛔ | -0.05 👍 |
| **[Quality](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#quality-score)** | **45.09%** 😞| **45.23%** 😞| **0.14%** 👍 |



| Other metrics | Before | After | Change |
| --- | :---: | :---: | :---: |
| Lines | 1992 | 2032 | 40 |



| Changed files | Quality Before | Quality After | Quality Change |
| --- | :---: | :---: | :---: |
| src/hist/basehist.py | 63.32% 🙂 | 61.43% 🙂 | -1.89% 👎|
| tests/test_general.py | 45.08% 😞 | 45.09% 😞 | 0.01% 👍|
| tests/test_named.py | 38.04% 😞 | 38.05% 😞 | 0.01% 👍|



Here are some functions in these files that still need a tune-up:

| File | Function | Complexity | Length | Working Memory | Quality | Recommendation |
| --- | --- | :---: | :---: | :---: | :---: | --- |
| tests/test_named.py | test\_named\_fill | 48 ⛔ | 1405 ⛔ | 19 ⛔ | 8.33% ⛔ | Refactor to reduce nesting. Try splitting into smaller methods. Extract out complex expressions |
| tests/test_general.py | test\_general\_fill | 48 ⛔ | 987 ⛔ | 19 ⛔ | 8.33% ⛔ | Refactor to reduce nesting. Try splitting into smaller methods. Extract out complex expressions |
| tests/test_named.py | test\_named\_project.test\_named\_index\_access | 0 ⭐ | 401 ⛔ | 27 ⛔ | 36.05% 😞 | Try splitting into smaller methods. Extract out complex expressions |
| tests/test_named.py | test\_named\_init | 6 ⭐ | 701 ⛔ | 14 😞 | 37.25% 😞 | Try splitting into smaller methods. Extract out complex expressions |
| src/hist/basehist.py | BaseHist.\_\_init\_\_ | 18 🙂 | 176 😞 | 14 😞 | 37.34% 😞 | Try splitting into smaller methods. Extract out complex expressions |


## Legend and Explanation

The emojis denote the absolute quality of the code:

 - ⭐ excellent
 - 🙂 good
 - 😞 poor
 - ⛔ very poor


The 👍 and 👎 indicate whether the quality has improved or gotten worse with this pull request.

***

Please see our documentation [here](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report) for details on how these metrics are calculated.

We are actively working on this report - lots more documentation and extra metrics to come!

Let us know what you think of it by mentioning @sourcery-ai in a comment.
sourcery-ai[bot]: # Sourcery Code Quality Report


✅ &nbsp;Merging this PR will **increase** code quality in the affected files by **0.14%**.

| Quality metrics | Before | After | Change |
| --- | :---: | :---: | :---: |
| [Complexity](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#complexity) | 9.37 🙂 | 9.20 🙂 | -0.17 👍 |
| [Method Length](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#method-length) | 163.71 😞 | 163.55 😞 | -0.16 👍 |
| [Working memory](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#working-memory) | 15.94 ⛔ | 15.89 ⛔ | -0.05 👍 |
| **[Quality](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report#quality-score)** | **45.58%** 😞| **45.72%** 😞| **0.14%** 👍 |



| Other metrics | Before | After | Change |
| --- | :---: | :---: | :---: |
| Lines | 1968 | 2008 | 40 |



| Changed files | Quality Before | Quality After | Quality Change |
| --- | :---: | :---: | :---: |
| src/hist/basehist.py | 62.38% 🙂 | 60.74% 🙂 | -1.64% 👎|
| tests/test_general.py | 45.52% 😞 | 45.54% 😞 | 0.02% 👍|
| tests/test_named.py | 38.44% 😞 | 38.45% 😞 | 0.01% 👍|



Here are some functions in these files that still need a tune-up:

| File | Function | Complexity | Length | Working Memory | Quality | Recommendation |
| --- | --- | :---: | :---: | :---: | :---: | --- |
| tests/test_named.py | test\_named\_fill | 48 ⛔ | 1329 ⛔ | 19 ⛔ | 8.33% ⛔ | Refactor to reduce nesting. Try splitting into smaller methods. Extract out complex expressions |
| tests/test_general.py | test\_general\_fill | 48 ⛔ | 943 ⛔ | 19 ⛔ | 8.33% ⛔ | Refactor to reduce nesting. Try splitting into smaller methods. Extract out complex expressions |
| tests/test_named.py | test\_named\_project.test\_named\_index\_access | 0 ⭐ | 401 ⛔ | 27 ⛔ | 36.05% 😞 | Try splitting into smaller methods. Extract out complex expressions |
| tests/test_named.py | test\_named\_init | 6 ⭐ | 671 ⛔ | 14 😞 | 37.27% 😞 | Try splitting into smaller methods. Extract out complex expressions |
| src/hist/basehist.py | BaseHist.\_\_init\_\_ | 18 🙂 | 176 😞 | 14 😞 | 37.34% 😞 | Try splitting into smaller methods. Extract out complex expressions |


## Legend and Explanation

The emojis denote the absolute quality of the code:

 - ⭐ excellent
 - 🙂 good
 - 😞 poor
 - ⛔ very poor


The 👍 and 👎 indicate whether the quality has improved or gotten worse with this pull request.

***

Please see our documentation [here](https://github.com/sourcery-ai/sourcery/wiki/Quality-Report) for details on how these metrics are calculated.

We are actively working on this report - lots more documentation and extra metrics to come!

Let us know what you think of it by mentioning @sourcery-ai in a comment.

Patch:
diff --git a/src/hist/basehist.py b/src/hist/basehist.py
index 626e3f0..179921f 100644
--- a/src/hist/basehist.py
+++ b/src/hist/basehist.py
@@ -282,6 +282,41 @@ class BaseHist(bh.Histogram, metaclass=MetaConstructor, family=hist):
 
         return super().__setitem__(self._index_transform(index), value)
 
+    def profile(self: T, axis: Union[int, str]) -> T:
+        """
+        Returns a profile (Mean/WeightedMean) histogram from a normal histogram
+        with N-1 axes. The axis given is profiled over and removed from the
+        final histogram.
+        """
+
+        if self.kind != bh.Kind.COUNT:
+            raise TypeError("Profile requires a COUNT histogram")
+
+        axes = list(self.axes)
+        iaxis = axis if isinstance(axis, int) else self._name_to_index(axis)
+        axes.pop(iaxis)
+
+        values = self.values()
+        tmp_variances = self.variances()
+        variances = tmp_variances if tmp_variances is not None else values
+        centers = self.axes[iaxis].centers
+
+        count = np.sum(values, axis=iaxis)
+
+        num = np.tensordot(values, centers, ([iaxis], [0]))
+        num_err = np.sqrt(np.tensordot(variances, centers ** 2, ([iaxis], [0])))
+
+        den = np.sum(values, axis=iaxis)
+        den_err = np.sqrt(np.sum(variances, axis=iaxis))
+
+        with np.errstate(invalid="ignore"):
+            new_values = num / den
+            new_variances = (num_err / den) ** 2 - (den_err * num / den ** 2) ** 2
+
+        retval = self.__class__(*axes, storage=hist.storage.Mean())
+        retval[...] = np.stack([count, new_values, count * new_variances], axis=-1)
+        return retval
+
     def density(self) -> np.ndarray:
         """
         Density NumPy array.


Test patch:
diff --git a/tests/test_general.py b/tests/test_general.py
index 1de998c..78ebf48 100644
--- a/tests/test_general.py
+++ b/tests/test_general.py
@@ -597,6 +597,7 @@ def test_general_transform_proxy():
         Hist.new.Sqrt(3, -4, 25)
 
     h1 = Hist.new.Log(4, 1, 10_000).Log(3, 1 / 1_000, 1).Double()
+
     h1.fill([2, 11, 101, 1_001], [1 / 999, 1 / 99, 1 / 9, 1 / 9])
     assert h1[0, 0] == 1
     assert h1[1, 1] == 1
@@ -604,11 +605,11 @@ def test_general_transform_proxy():
     assert h1[3, 2] == 1
 
     # Missing arguments
-    with pytest.raises(Exception):
-        Hist.new.Regular(4, 1, 10_000).Log()
+    with pytest.raises(TypeError):
+        Hist.new.Reg(4, 1, 10_000).Log()
 
-    # wrong value
-    with pytest.raises(Exception):
+    # Wrong value
+    with pytest.raises(ValueError):
         Hist.new.Log(3, -1, 10_000)
 
     h2 = Hist.new.Pow(24, 1, 5, power=2).Pow(124, 1, 5, power=3).Int64()
@@ -636,25 +637,28 @@ def test_general_transform_proxy():
             4, 1, 5, forward=ftype(np.log), inverse=ftype(np.exp)
         )
     ).Int64()
+
     h3.fill([1, 2, 3, 4], [1, 2, 3, 4])
     assert h3[0, 0] == 1
     assert h3[1, 1] == 1
     assert h3[2, 2] == 1
     assert h3[3, 3] == 1
 
-    # wrong value
-    with pytest.raises(Exception):
-        Hist().Regular(24, 1, 5).Func(ftype(math.log), ftype(math.exp))
+    # wrong args
+    with pytest.raises(TypeError):
+        Hist.new.Reg(24, 1, 5).Func(ftype(math.log), ftype(math.exp))
 
-    # wrong value
-    assert Hist.new.Func(
-        4, -1, 5, forward=ftype(math.log), inverse=ftype(math.log)
-    ).Double()
-    with pytest.raises(Exception):
-        Hist.new.Func(4, -1, 5, forward=ftype(np.log), inverse=ftype(np.log))
+    # wrong value raises uncatchable warning
+    # Hist.new.Func(
+    #     4, -1, 5, forward=ftype(math.log), inverse=ftype(math.log)
+    # ).Double()
+
+    with pytest.raises(ValueError):
+        with pytest.warns(RuntimeWarning):
+            Hist.new.Func(4, -1, 5, forward=ftype(np.log), inverse=ftype(np.log))
 
     # lack args
-    with pytest.raises(Exception):
+    with pytest.raises(TypeError):
         Hist.new.Func(4, 1, 5)
 
 
diff --git a/tests/test_named.py b/tests/test_named.py
index e0076b2..b401dcf 100644
--- a/tests/test_named.py
+++ b/tests/test_named.py
@@ -649,12 +649,12 @@ def test_named_transform_proxy():
     assert h1[2, 2] == 1
     assert h1[3, 2] == 1
 
-    # wrong value
-    with pytest.raises(Exception):
-        NamedHist.new.Regular(4, 1, 10_000, name="x").Log()
+    # wrong arguments
+    with pytest.raises(TypeError):
+        NamedHist.new.Reg(4, 1, 10_000, name="x").Log()
 
     # wrong value
-    with pytest.raises(Exception):
+    with pytest.raises(ValueError):
         NamedHist.new.Log(3, -1, 10_000, name="x")
 
     h2 = (
@@ -669,15 +669,15 @@ def test_named_transform_proxy():
     assert h2[15, 63] == 1
 
     # based on existing axis
-    with pytest.raises(Exception):
-        NamedHist.new.Regular(24, 1, 5, name="x").Pow(2)
+    with pytest.raises(TypeError):
+        NamedHist.new.Reg(24, 1, 5, name="x").Pow(2)
 
     # wrong value
-    with pytest.raises(Exception):
+    with pytest.raises(ValueError):
         NamedHist.new.Pow(24, -1, 5, power=1 / 2, name="x")
 
     # lack args
-    with pytest.raises(Exception):
+    with pytest.raises(TypeError):
         NamedHist.new.Pow(24, 1, 5, name="x")
 
     ftype = ctypes.CFUNCTYPE(ctypes.c_double, ctypes.c_double)
@@ -695,20 +695,22 @@ def test_named_transform_proxy():
     assert h3[3, 3] == 1
 
     # based on existing axis
-    with pytest.raises(Exception):
-        NamedHist.new.Regular(24, 1, 5, name="x").Func(ftype(math.log), ftype(math.exp))
+    with pytest.raises(TypeError):
+        NamedHist.new.Reg(24, 1, 5, name="x").Func(ftype(math.log), ftype(math.exp))
 
-    # wrong value
-    assert NamedHist.new.Func(
-        4, -1, 5, name="x", forward=ftype(math.log), inverse=ftype(math.log)
-    )
-    with pytest.raises(Exception):
-        NamedHist.new.Func(
-            4, -1, 5, name="x", forward=ftype(np.log), inverse=ftype(np.log)
-        )
+    # Uncatchable warning
+    # assert NamedHist.new.Func(
+    #     4, -1, 5, name="x", forward=ftype(math.log), inverse=ftype(math.log)
+    # )
+
+    with pytest.raises(ValueError):
+        with pytest.warns(RuntimeWarning):
+            NamedHist.new.Func(
+                4, -1, 5, name="x", forward=ftype(np.log), inverse=ftype(np.log)
+            )
 
     # lack args
-    with pytest.raises(Exception):
+    with pytest.raises(TypeError):
         NamedHist.new.Func(4, 1, 5, name="x")
 
 
@@ -720,7 +722,7 @@ def test_named_hist_proxy():
     h = NamedHist.new.Reg(10, 0, 1, name="x").Double().fill(x=[0.5, 0.5])
     assert h[0.5j] == 2
 
-    assert type(h) == NamedHist
+    assert type(h) is NamedHist
 
     with pytest.raises(AttributeError):
         NamedHist().new
diff --git a/tests/test_profile.py b/tests/test_profile.py
new file mode 100644
index 0000000..0d10913
--- /dev/null
+++ b/tests/test_profile.py
@@ -0,0 +1,35 @@
+import numpy as np
+from pytest import approx
+
+import hist
+
+
+def test_example():
+    xy = np.array(
+        [
+            [-2, 1.5],
+            [-2, -3.5],
+            [-2, 1.5],  # x = -2
+            [0.0, -2.0],
+            [0.0, -2.0],
+            [0.0, 0.0],
+            [0.0, 2.0],
+            [0.0, 4.0],  # x = 0
+            [2, 1.5],  # x = +2
+        ]
+    )
+    h = hist.Hist(
+        hist.axis.Regular(5, -5, 5, name="x"), hist.axis.Regular(5, -5, 5, name="y")
+    ).fill(*xy.T)
+
+    # Profile out the y-axis
+    hp = h.profile("y")
+
+    # Exclude edge bins since no values from above will fall into them
+    # When there are values there, ROOT does something funky in those bins,
+    # despite these bins not being in the axis that is profiled out, and
+    # despite there being no overflow... to be understood.
+    assert hp.values()[1:-1] == approx(np.array([0.0, 0.4, 2.0]))
+    assert hp.variances()[1:-1] == approx(
+        np.array([2.66666667, 1.088, float("nan")]), nan_ok=True
+    )