[FEATURE] "group by" for cat axes / index based slicing

[andrzejnovak]

I don't think this is currently implemented, but would be super useful, allowing to merge samples that were processed separately.

I am imagining syntax like:
h[{'category: {'merged': ['sampleA', 'sampleB'], ...}}]

Also I thought

h[...,[0, 2], ...]

would work, but it doesn't seem to be possibly currently.

[andrzejnovak]

@henryiii is there currently no way to do this even manually by assignment due to https://github.com/scikit-hep/boost-histogram/blob/6f4813f1e4b326ca14074b739f2214383f46bec6/src/boost_histogram/_internal/hist.py#L819 ?

[henryiii]

h2[...] = h1.view(flow=True) allows assignment. Histogram histogram assignment would probably need some checking on the axes which is not implemented yet.

[henryiii]

h[...,[0, 2], ...] would probably need (or be best with) support in Boost.Histogram, see scikit-hep/boost-histogram#296. @HDembinski, is this something that can be supported upstream? If we have to, we can implement it in boost-histogram via a workaround.

PS: Assuming this is for unordered axes only.

[henryiii]

PS: The issue that got opened and fixed in Boost.Histogram was for slicing on categorical axes, which enabled h[..., 0:2, ...] to work, but not selecting a subset of a categorical axes as is requested here and mentioned on the original boost-histogram issue.

[andrzejnovak]

Ok thanks. In case anyone stumbles here. This seems to be the workaround, thanks @henryiii

def groupby(h, groupmap, axis='dataset'):
    new = Hist(*[ax for ax in h.axes if ax.name != axis], 
                hist.axis.StrCategory(groupmap.keys(), name=axis, growth=True), 
                hist.storage.Weight()
          )

    for name, cats in groupmap.items():
        grouped = sum([h[{axis: name}] for name in cats])
        new[{axis: name}] = grouped.view(flow=True)
    return new

[HDembinski]

Seems like a nice feature for boost-histogram.

[andrzejnovak]

Related issue. This new[{axis: name}] = grouped.view(flow=True) syntax fails when growth axis dimensions don't match.

[henryiii]

Seems like a nice feature for boost-histogram.

boost-histogram doesn't have named axes, so it wouldn't be as pretty, and would need another layer of wrapping in Hist anyway, just like fill, project, ... (not against it, but probably best to implement it here first)

This new[{axis: name}] = grouped.view(flow=True) syntax fails when growth axis dimensions don't match.

How would it know what entries to add?

h[...,[0, 2], ...]

This is almost implementable on top of scikit-hep/boost-histogram#576, save for the caveats mentioned there.

[andrzejnovak]

How would it know what entries to add?

Admittedly I didn't think about it too deeply, but it could just pad zeros to the dimensions along the missing categorical entries? Should be equivalent to adding two histograms where the growth/cat axis has different entries?

[nsmith-]

Since we know the new axis elements already (the dictionary keys) I think we could have a workaround without growth as follows:

import hist

def group(h: hist.Hist, oldname: str, newname: str, grouping: dict[str, list[str]]):
    hnew = hist.Hist(
        hist.axis.StrCategory(grouping, name=newname),
        *(ax for ax in h.axes if ax.name != oldname),
        storage=h._storage_type,
    )
    for i, indices in enumerate(grouping.values()):
        hnew.view(flow=True)[i] = h[{oldname: indices}][{oldname: sum}].view(flow=True)

    return hnew

Note that the new axis is put at the beginning (for convenience in implementation). I couldn't find a public accessor for the storage type though.

An example

h = (
    hist.Hist.new
    .StrCat("abcde", name="letter")
    .Reg(10, 0, 1, name="number")
    .Double()
)

grouping = {
    "vowel": ["a", "e"],
    "consonant": ["b", "c", "d"],
}

print(group(h, "letter", "type", grouping))

returning

Hist(
  StrCategory(['vowel', 'consonant'], name='type', label='type'),
  Regular(10, 0, 1, name='number', label='number'),
  storage=Double())

[nsmith-]

A small update to my previous comment: the workaround now needs h._storage_type() due to a warning about passing the type and not an instance. Perhaps we can have a public accessor for the storage type that is stable?