Resolves #148

This PR adds support for ratio plots to the .plot... API by adding .plot_ratio and refactoring .plot_pull to use .plot_ratio.

This PR is still very rough, but I thought I'd open it up as there is a minimal working example now and then revise it heavily from there.

As the commit log is huge and messy I'll rebase and squash things into more reasonable commit sections to make review more digestible.

Questions to come to consensus on

(Discussion can happen on Issue #148)

• Should the interval functions be moved to their own stat_interval module?

Answer: Yes, this is PR #176.

TODO before requesting review

• [x] Finish refactoring .plot_pull
• [x] Add kwargs and kwargs filtering support
• [x] Fix typing errors
• [x] Add tests
• [x] All tests are passing
• [x] Docstrings added
• [x] README updated
• [x] Changelog updated.

Suggested squash and merge message

* Add .plot_ratio to BastHist API
   - Adds ratio plot support for other Hists and callables
* Factor out plot_ratio and plot_pull to perform the final subplot plotting
   - Majority of logic is factored out of plot_pull into _plot_ratiolike
   - _plot_ratiolike then calls plot_ratio or plot_pull as needed
* Add plotting tests and test images for plot_ratio
* Update test image for plot_pull
* Update README and Changelog to reflect addition of plot_ratio

Co-authored-by: Henry Schreiner <[email protected]>

I was playing with Hist.plot_pull and noticed that the range of the pulls is dynamic, meaning it is calculated from max(np.abs(pulls)) / pp_num in the code. This is not ideal since (1) pull plots are basically always displayed between -5 and 5, and (2) having several pull plots side-by-side would likely provide ranges on a plot-by-plot case, which is not ideal.

This little PR fixes the issue. It also makes sure that the colour bands by default are set at 1/2/.../5 sigma.

I will post an issue with some other matters related to the function.

Next week I will be giving some lectures and will be showing Hist. If at all possible it would be great to have this in a release :-).

https://github.com/scikit-hep/hist/blob/be2c8380793807966f1ab58c4e50962fae7e45f9/src/hist/_internal/axis.py#L5-L17

Should every axis be forced to have a name? I see you let title to be an arg (cannot be omitted). Hence that name should have the same status, it is should be indispensable, too. But in the unit tests, we are supposed to pass the test where axes have no names, e.g.

https://github.com/scikit-hep/hist/blob/be2c8380793807966f1ab58c4e50962fae7e45f9/tests/test_general.py#L5

and

https://github.com/scikit-hep/hist/blob/be2c8380793807966f1ab58c4e50962fae7e45f9/tests/test_named.py#L6

P.S. The current practice, i.e., #13, is to modify the tests to enforce them to have names. An alternative plan is to make title and name both omissible, which might be subject to the issue of filling-by-name (this functionality should be thrown out for hist with anonymous axes).

Based on the discussion in https://github.com/scikit-hep/hist/issues/146, I added some features to plot_pull. The theme is making fitting more streamlined for exploration.

• Pull curve_fit initial guess (p0) from default arguments, if they exist
• Allow a string as an alternative to a lambda function (plot_pull("a+b*x"))
• Cosmetic change to the band, and embedding fit result into the legend
• Likelihood fit (plot_pull(..., likelihood=True)) (chi2 by default, as before)
• Remove uncertainties.numpy dependency and construct band by resampling covariance matrix
• Introduce iminuit (gets the covariance matrix right, unlike scipy.optimize most of the time), but the initial guesses for iminuit are seeded from scipy

Setup

import numpy as np
from hist import Hist

np.random.seed(42)
hh = Hist.new.Reg(50, -5, 5).Double().fill(np.random.normal(0,1,int(1e5)))

Before (including a bug-fix for the variances from the above issue):

from uncertainties import unumpy as unp
def func(x, constant, mean, sigma):
    exp = unp.exp if constant.dtype == np.dtype("O") else np.exp
    return constant * exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))

hh.plot_pull(func)

After:

# as before, but no need for `uncertainties.numpy` as the error band comes
# from resampling the covariance matrix
def func(x, constant, mean, sigma):
    return constant * np.exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))
hh.plot_pull(func)

# `curve_fit` `p0` extracted from defaults, if any
def func(x, constant=80, mean=0., sigma=1.):
    return constant * np.exp(-((x - mean) ** 2.0) / (2 * sigma ** 2))
hh.plot_pull(func)

# strings are allowed to allow for more compactness than a lambda
# x is assumed to be the main variable
hh.plot_pull("constant*np.exp(-(x-mean)**2. / (2*sigma**2))")

# gaussian is a common/special function, so this also works
# reasonable guesses are made for constant/mean/sigma
hh.plot_pull("gaus")

# chi2 puts `a` around 5, but likelihood puts `a` around 1e3/50 = 20
hh.plot_pull("a+b*x", likelihood=True)

Describe the problem, if any, that your feature request is related to

Don't know how to add typing hints for these two func.

https://github.com/scikit-hep/hist/blob/dc9b209dfa5d2fa934a803c9cb589b414314fc13/src/hist/named.py#L50

https://github.com/scikit-hep/hist/blob/dc9b209dfa5d2fa934a803c9cb589b414314fc13/src/hist/named.py#L66

When completed, closes #222

To implement:

• sort by passing index
• sort by passing lambda to sorted
• add tests

@henryiii could you take a look? I had to abuse some design decisions of hist to get this work, so maybe you can recommend better workarounds.

• AxisNamedTuple (and I assume axes) as well cannot be assigned to, so I had to recreate the sorted axis
• Couldn't easily copy all meta/traits from old axis to new - hence the helper with inspect (name/label don't seem to be part of traits/metadata) - also cannot create axis and edit metadata

Here is a couple of examples 2D - yax 2D - xax 1D 3D - sort on not-projected axis

Addressing #35.

Currently it seems that axes without names are broken in Hist. If I do a new pip install -e .[dev], then try, from a command line:

hist
1
2
3
2

Then press control-D, I get:

Traceback (most recent call last):
  File "/Users/henryschreiner/git/scikit-hep/hist/.env/bin/hist", line 33, in <module>
    sys.exit(load_entry_point('hist', 'console_scripts', 'hist')())
  File "/Users/henryschreiner/git/scikit-hep/hist/src/hist/classic_hist.py", line 22, in main
    h = bh.numpy.histogram(values, bins=args.buckets, histogram=hist.Hist)
  File "/Users/henryschreiner/git/scikit-hep/hist/.env/lib/python3.8/site-packages/boost_histogram/numpy.py", line 111, in histogram
    return histogramdd((a,), (bins,), (range,), normed, weights, density, **kwargs)
  File "/Users/henryschreiner/git/scikit-hep/hist/.env/lib/python3.8/site-packages/boost_histogram/numpy.py", line 74, in histogramdd
    hist = cls(*axs, storage=bh_storage).fill(*a, weight=weights, threads=threads)
  File "/Users/henryschreiner/git/scikit-hep/hist/src/hist/core.py", line 38, in __init__
    if ax.name in self.names:
  File "/Users/henryschreiner/git/scikit-hep/hist/src/hist/_internal/axis.py", line 52, in name
    return self.metadata["name"]
TypeError: 'NoneType' object is not subscriptable

@LovelyBuggies, could you take a look? This is really just running bh.numpy.histogram([1,2,3,2], histogram=hist.Hist).

Describe your questions

I have used mypy according to https://scikit-hep.org/developer/style#type-checking-new. My pre-commit profig file is like this:

repos:
- repo: https://github.com/psf/black
  rev: 19.10b0
  hooks:
  - id: black
- repo: https://github.com/pre-commit/pre-commit-hooks
  rev: v2.5.0
  hooks:
  - id: check-added-large-files
  - id: mixed-line-ending
  - id: trailing-whitespace
  - id: check-merge-conflict
  - id: check-case-conflict
  - id: check-symlinks
  - id: check-yaml
- repo: https://github.com/pre-commit/mirrors-mypy
  rev: v0.770
  hooks:
  - id: mypy
    files: all  # I tried tests and src, but none worked

I deliberately put a type error in a test file, but pre-commit did not find it. How to make mypy CI work for unit tests?

Closes #156, adds a .profile method. Based heavily on @aminnj's example.

TODO:

• [x] Needs tests
• This probably could be a WeightedMean, probably give a choice? Just got the basics working for now. (Can revisit later)
• [x] Should hide or handle warnings (true for some of our tests, too)

Followup:

• histoprint should handle showing kind=MEAN histograms (it doesn't like the NaN's, I think) CC @ast0815 (edit: done!)
• mplhep should show a better plot style for kind=MEAN, and should also ensure it handles kind=MEAN's requirements (we could also dispatch differently, but the best fix would be in mplhep) CC @andrzejnovak
• Boost-histogram should describe accumulators/storages a bit more in the docs, especially setting with a stack, including the correct order. CC Me

If kind=MEAN, you are only supposed to plot variances if counts() > 1, and values if counts() >= 1.

• Does named indexing work on all Hist's? I don’t see the code for it in BaseHist, only in NamedHist. Something like this: h[{“x”:2}] should be converted into h[{0:2}] assuming that axes 0 has name "x".
• NamedHist should simply verify that the index item is a dict, and that it has only string keys. Then it should just call BaseHist's indexing.
• NamedHist's fill should simply be: def fill(self, **kwargs): return super().fill(**kwargs) (maybe with a few more keyword arguments listed to be nice to the user's inspection tools). BaseHist should be able to fill via kwargs or via position.

Overall, NamedHist should be pretty short, it is just taking away functionally (position-based access) from BaseHist. I'm not sure how strict we want to be - we can also have an OnlyNamedAxesTuple too.

Also we should make sure you can cast between Hist and NamedHist. Hist(NamedHist(...)), etc. I can do that one.

Describe the bug

The error bars produced in the ratio of an efficiency style plot (hist.plot_ratio(..., rp_uncertainty_type = "poisson-ratio", ...)) can extend above unity. This is unexpected and not meaningful if the numerator is a true subset of the denominator.

Steps to reproduce

This is observed using hist==2.4.0.

Following the example from the docs,

hist0 = hist.Hist(hist.axis.Regular(50, -5, 5), underflow = False, overflow = False).fill(np.random.normal(size=1700))
hist1 = hist0.copy() * 0.98

fig = plt.figure(figsize = (10,8))
_ = hist1.plot_ratio(
    hist0,
    rp_ylabel="Efficiency",
    rp_num_label="hist1",
    rp_denom_label="hist0",
    rp_uncert_draw_type="line",
    rp_uncertainty_type="poisson-ratio",
    rp_ylim=[0.9,1.1],
)

This produces something like the following,

where we see the error bars going above unity.

If I take the same data used in hist0 and hist1 above, and produce a TEfficiency object in ROOT, I get something like the following, which is something more like what I would expect since I believe hist is using a Clopper-Pearson interval (default in ROOT).

Describe the problem, if any, that your feature request is related to

Right now the https://github.com/dask-contrib/dask-histogram package exposes an interface that is very nearly 1:1 with boost histogram, which I think is the correct scope for this package. However, many users of columnar analysis tools in HEP prefer hist since its UI is more intuitive and comfortable. It would be best to have a sub-module within hist that exposed a dask-histogram backed interface and suite of tools.

Describe the feature you'd like

I would like users to be able to do something like:

import hist.dask as hist

or

import hist

h = hist.Hist(**axes).as_dask()

and have access to the hist interface with (nearly) all of its features but backed by computation in dask.

This is very close to practice since dask_histogram already very closely follows (and mostly directly uses) the boost-histogram interface, which means that most of what hist does should be easy to adapt.

This will require a few pieces of interface to be smoothened out in dask_histogram since it's not 1:1 with boost-histogram in some cases. (e.g. what is mentioned in the bottom of the comment here: https://github.com/dask-contrib/dask-histogram/issues/35#issuecomment-1368108690)

Describe alternatives, if any, you've considered

The alternative is using the dask histogram interface directly, which while effective, is significantly less widely used and less pleasant than hist. This would require significant readoption and generate attrition risk in our user base.

Bumps pypa/gh-action-pypi-publish from 1.5.1 to 1.6.4.

Dependabot will resolve any conflicts with this PR as long as you don't alter it yourself. You can also trigger a rebase manually by commenting @dependabot rebase.

In:

import matplotlib.pyplot as plt
import numpy
from   hist import Hist

#------------------------------------------------------
def get_hist():
    h = (
      Hist.new
      .Reg(10, -2 ,2, name='x', label='x-axis')
      .Int64()
    )
    data=numpy.random.normal(size=10000)
    h.fill(data)

    return h
#------------------------------------------------------
def main():
    h_1 = get_hist()
    h_2 = get_hist()
    h_3 = get_hist()

    fig, (ax, rax) = plt.subplots(nrows=2, gridspec_kw={"height_ratios": (3, 1)})
    axs = {"main_ax": ax, "ratio_ax": rax}

    h_2.plot_ratio(h_1,
        rp_ylabel          ='Ratio',
        rp_num_label       ='data',
        rp_denom_label     ='sim 1',
        rp_uncert_draw_type='bar',
        ax_dict = axs,
    )

    h_3.plot_ratio(h_1,
        rp_ylabel          ='Ratio',
        rp_num_label       ='data',
        rp_denom_label     ='sim 2',
        rp_uncert_draw_type='bar',
        ax_dict = axs,
    )

    plt.show()
#------------------------------------------------------
if __name__ == '__main__':
    main()

The code is unable to:

1. Plot the ratio of data to both simulations and there is no easy way to make this work. I see repeated histograms.
2. The axes are not aligned by default
3. I do not see documentation (or at least it is not easy to find) fr the rp_* arguments.

Describe the feature you'd like

The user should be able to do:

    h_1 = get_hist()
    h_2 = get_hist()
    h_3 = get_hist()

   h_1.plot_ratio([h_2, h_3])

   plt.show()

and we should get two figures, upper one with h_* overlaid. Lower with two ratios, to h_2 and h_3. The axes should be aligned, the labels and legends should be taken from the histograms themselves and we should not have to do any manipulation of artists.

Describe alternatives, if any, you've considered

The way the code is implemented is bad, it's too complicated, and I do not have time to make it work the way I need to, so I am moving back to pure matplotlib. The plots I need do not need to be perfect and matplotlib is good enough for me. It would be nice if hist can do quickly what we need though.

Cheers.

updates:

• github.com/psf/black: 22.10.0 → 22.12.0
• github.com/pre-commit/pre-commit-hooks: v4.3.0 → v4.4.0
• github.com/PyCQA/isort: 5.10.1 → 5.11.4
• github.com/asottile/pyupgrade: v3.2.0 → v3.3.1
• github.com/nbQA-dev/nbQA: 1.5.3 → 1.6.0
• github.com/pycqa/flake8: 5.0.4 → 6.0.0
• github.com/pre-commit/mirrors-mypy: v0.982 → v0.991
• github.com/mgedmin/check-manifest: 0.48 → 0.49
• github.com/shellcheck-py/shellcheck-py: v0.8.0.4 → v0.9.0.2

Describe the bug

Adding two hists of different storage types just results in having an empty hist without raising an error:

Steps to reproduce

import numpy as np
from hist import Hist
import hist
a_hist = hist.Hist(hist.axis.Regular(3, -1, 1), storage=hist.storage.Int64())
a_hist.sum()

0.0

b_hist = hist.Hist(hist.axis.Regular(3, -1, 1), storage=hist.storage.Double())
b_hist.fill(np.random.normal(size=1000))
b_hist.sum()

681.0

a_hist += b_hist
a_hist.sum()

0.0

Describe the problem, if any, that your feature request is related to I believe the current plot_ratio() method does not take into account weights of the histograms when calculating the errors.

Describe the feature you'd like

Propagate errors correctly into the ratio uncertainty, taking into account the weights

Describe alternatives, if any, you've considered

Using coffea.hist.plotratio instead