<jupyter_start><jupyter_text>XArray exampleTo run this example, an environment.yml similar to this one could be used:```yamlname: bh_xhistogramchannels:  - conda-forgedependencies:  - python==3.8  - boost-histogram  - xhistogram  - matplotlib  - netcdf4```<jupyter_code>import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import boost_histogram as bh
from xhistogram.xarray import histogram as xhistogram<jupyter_output><empty_output><jupyter_text>Let's look at using boost-histogram to imitate the xhistogram package by reading and producing xarrays. As a reminder, xarray is a sort of generalized Pandas library, supporting ND labeled and indexed data.  Simple 1D exampleWe will start with the first example from the xhistogram docs:<jupyter_code>da = xr.DataArray(np.random.randn(100, 30),
                  dims=['time', 'x'],
                  name='foo')
bins = np.linspace(-4, 4, 20)<jupyter_output><empty_output><jupyter_text>xhistogramAnd, this is what historamming and plotting looks like:<jupyter_code>h = xhistogram(da, bins=[bins])
display(h)
h.plot(); # h is an xarray<jupyter_output><empty_output><jupyter_text>boost-histogram (direct usage)Let's first just try this by hand, to see how it works. This will not return an xarray, etc.<jupyter_code>bh_bins = bh.axis.Regular(19, -4, 4)
bh_hist = bh.Histogram(bh_bins).fill(np.asarray(da).flatten())
plt.plot(bh_hist.axes[0].centers, bh_hist);<jupyter_output><empty_output><jupyter_text>boost-histogram (adapter function)Now, let's make an adaptor for boost-histogram.<jupyter_code>def bh_xhistogram(*args, bins):
    
    # Convert bins to boost-histogram axes first
    prepare_bins = (bh.axis.Variable(b) for b in bins)
    h = bh.Histogram(*prepare_bins)
    
    # We need flat NP arrays for filling
    prepare_fill = (np.asarray(a).flatten() for a in args)
    h.fill(*prepare_fill)
    
    # Now compute the xarray output.
    return xr.DataArray(h,
            name="_".join(a.name for a in args) + "_histogram",
            coords=[(f"{a.name}_bin", arr.flatten(), a.attrs)
                      for a, arr in zip(args, h.axes.centers)])
h = bh_xhistogram(da, bins=[bins])
display(h)
h.plot();<jupyter_output><empty_output><jupyter_text>More featuresLet's add a few more features to our function defined above.* Let's allow bins to be a list of axes or even a completly prepared histogram; this will allow us to take advantage of boost-histogram features later.* Let's add a weights keyword so we can do weighted histograms as well.<jupyter_code>def bh_xhistogram(*args, bins, weights=None):
    """
    bins is either a histogram, a list of axes, or a list of bins
    """
    
    if isinstance(bins, bh.Histogram):
        h = bins
    else:
        prepare_bins = (b if isinstance(b, bh.axis.Axis) else bh.axis.Variable(b) for b in bins)
        h = bh.Histogram(*prepare_bins)
    
    prepare_fill = (np.asarray(a).flatten() for a in args)
    
    if weights is None:
        h.fill(*prepare_fill)
    else:
        prepared_weights, *_ = xr.broadcast(weights, *args)
        h.fill(*prepare_fill, weight=np.asarray(prepared_weights).flatten())
    
    return xr.DataArray(h,
            name="_".join(a.name for a in args) + "_histogram",
            coords=[(f"{a.name}_bin", arr.flatten(), a.attrs)
                      for a, arr in zip(args, h.axes.centers)])<jupyter_output><empty_output><jupyter_text>2D exampleThis also comes from the xhistogram docs.<jupyter_code># Read WOA using opendap
Temp_url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/WOA/WOA13/5_deg/annual/temp'
Salt_url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/WOA/WOA13/5_deg/annual/salt'
Oxy_url = 'http://apdrc.soest.hawaii.edu:80/dods/public_data/WOA/WOA13/5_deg/annual/doxy'

ds = xr.merge([
    xr.open_dataset(Temp_url).tmn.load(),
    xr.open_dataset(Salt_url).smn.load(),
    xr.open_dataset(Oxy_url).omn.load()])
sbins = np.arange(31,38, 0.025)
tbins = np.arange(-2, 32, 0.1)<jupyter_output><empty_output><jupyter_text>xhistogram<jupyter_code>hTS = xhistogram(ds.smn, ds.tmn, bins=[sbins, tbins])
np.log10(hTS.T).plot(levels=31)<jupyter_output>/usr/local/Caskroom/miniconda/base/envs/xtest/lib/python3.8/site-packages/xarray/core/computation.py:601: RuntimeWarning: divide by zero encountered in log10
  result_data = func(*input_data)<jupyter_text>boost-histogramWe could hand in the same bin definitions, but let's use boost-histogram axes instead:<jupyter_code>sax = bh.axis.Regular(250, 31, 38)
tax = bh.axis.Regular(340, -2, 32)

hTS = bh_xhistogram(ds.smn, ds.tmn, bins=[sax, tax])
np.log10(hTS.T).plot(levels=31)<jupyter_output><empty_output><jupyter_text>Speed comparson<jupyter_code>%%timeit
hTS = xhistogram(ds.smn, ds.tmn, bins=[sbins, tbins])
%%timeit
hTS = bh_xhistogram(ds.smn, ds.tmn, bins=[sax, tax])<jupyter_output>4.7 ms ± 245 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>Weighted histogram Let's try a more complex example from the docs; the dVol weights one:<jupyter_code>dz = np.diff(ds.lev)
dz =np.insert(dz, 0, dz[0])
dz = xr.DataArray(dz, coords= {'lev':ds.lev}, dims='lev')

dVol = dz * (5*110e3) * (5*110e3*np.cos(ds.lat*np.pi/180))

hTSw = bh_xhistogram(ds.smn, ds.tmn, bins=[sax, tax], weights=dVol)
np.log10(hTSw.T).plot(levels=31, vmin=11.5, vmax=16, cmap='brg')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13135830072666294, 'fraction_numerical': 0.020122973728339856, 'mean_word_length': 4.368, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '325383', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2140, 'n_tokens_neox': 1989, 'n_words': 462, 'path': 'docs/notebooks/xarray.ipynb', 'repo_name': 'HDembinski/boost-histogram', 'size': '5367'}