Accessing sum of events per bin per channel

[ssaha1234]

Dear pyhf community,
I am trying to access the sum of all my simulation background, per channel. This has to be done per bin, in a given channel.

I know how to get the complete array of the list of the sum of events, using the following lines of code -

json_data = json.load(file)
wspace = pyhf.Workspace(json_data)
model = wspace.model()
data = wspace.data(model)
params = model.config.suggested_init()
params[model.config.poi_index] = 0.0
sum_bkg = model.main_model.expected_data(params, False)

The sum_bkg above, prints an array of the complete number of channels, per bin. I have two questions regarding the array -

1. I would like to exclude some samples from the sum, is there a way to use the sample name and exclude it (maybe in a loop, or using list comprehension) ?

2. The other thing I would like to do is to loop over this array, per bin, per channel, to perform pruning of the samples which do not contribute much. I know how to do the pruning, but I am unable to access the sum, per bin. It would be great if you could point me to some examples to do the same.

Thank you very much !
Best,
Shreya

[kratsg]

1. You're already using the right API, so just call model.main_model.expected_data(params, True) which returns by sample (this isn't in our documentation at the moment).
2. You can make the sum per bin by identifying the samples you want to remove, so I don't understand the question here.

[ssaha1234]

Hi,
I have a quick followup question, I managed to get yields per sample, per bin, how do I access the names of the samples and channels in the array in model_yields below. I have done something like this -

params = model.config.suggested_init()
params[model.config.poi_index] = 0.0

sum_bkg = model.main_model.expected_data(params, True)

model_yields = [
        sum_bkg[:, model.config.channel_slices[ch]].tolist()
        for ch in model.config.channels
]

The other idea would be to remove the samples/channels before making model_yields ?
Best,
Shreya

[alexander-held]

You could simultaneously loop over the content per sample and the sample names, and remove pieces of that loop depending on the name. Here are a few ideas of different approaches and structures. They probably use list comprehension a bit too much and it might be more readable to break things up into pieces.

import pyhf
import numpy as np

model = pyhf.simplemodels.uncorrelated_background(
    signal=[24.0, 22.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
)

yields_combined = model.main_model.expected_data(
    model.config.suggested_init(), return_by_sample=True
)

# approach 1
# nested dictionaries indexed by channel name and sample name
model_yields = {
    ch: {
        sample_name: sample_yield
        for sample_name, sample_yield in zip(
            model.config.samples,
            yields_combined[:, model.config.channel_slices[ch]].tolist(),
        )
    }
    for ch in model.config.channels
}
print(model_yields)
# can now filter out by name
for channel_name, channel_yields in model_yields.items():
    for sample_name, sample_yields in channel_yields.items():
        if sample_name == "background":
            continue
        else:
            print(sample_yields)


# approach 2
# removing specific samples when building the list
sample_blacklist = ["background"]
channel_blacklist = [""]
model_yields = [
    [
        sample_yield
        for sample_yield, sample_name in zip(
            yields_combined[:, model.config.channel_slices[ch]].tolist(),
            model.config.samples,
        )
        if sample_name not in sample_blacklist
    ]
    for ch in model.config.channels
    if ch not in channel_blacklist
]
print(model_yields)


# approach 3
# convert to numpy array and slice
model_yields = [
    yields_combined[:, model.config.channel_slices[ch]].tolist()
    for ch in model.config.channels
]
model_yields = np.asarray(model_yields)
channel_slice = [
    idx
    for idx in np.arange(model_yields.shape[0])
    if model.config.channels[idx] not in channel_blacklist
]
sample_slice = [
    idx
    for idx in np.arange(model_yields.shape[1])
    if model.config.samples[idx] not in sample_blacklist
]
print(model_yields[(channel_slice, sample_slice)])

edit: the last approach might not be working correctly for models with multiple channels, I have not debugged it further

[ssaha1234]

Hi @alexander-held ,
Thank you very much for the piece of code, I understood the filtering of samples, and I could remove the samples I do not need, and make an array of the yields, per channel, per sample, per bin. Now I am thinking of summing up the yields in each bin, per channel. Looking at the above example, I would like to sum up signal and background events, given they are both the samples I would like to have. I do know that using sum_bkg = model.main_model.expected_data(params, False) will give me the total sum, but again the mapping to samples and channels and filtering is complicated. I would like to think that is an easy python exercise, but in case you have more insights, let me know !
Thanks again !
Best,
Shreya

[alexander-held]

Do you want to sum over bins and samples? Here is an example of how to do that:

import pyhf
import numpy as np

model = pyhf.simplemodels.uncorrelated_background(
    signal=[24.0, 22.0], bkg=[50.0, 52.0], bkg_uncertainty=[3.0, 7.0]
)

yields_combined = model.main_model.expected_data(
    model.config.suggested_init(), return_by_sample=True
)

# approach 2
# removing specific samples when building the list
sample_blacklist = ["background"]
channel_blacklist = [""]
model_yields = [
    [
        sample_yield
        for sample_yield, sample_name in zip(
            yields_combined[:, model.config.channel_slices[ch]].tolist(),
            model.config.samples,
        )
        if sample_name not in sample_blacklist
    ]
    for ch in model.config.channels
    if ch not in channel_blacklist
]
print(model_yields)

# take yields after filtering out channels / samples
print(f"per channel, sample, bin: {model_yields}")

# sum over bins (the innermost index)
# indices are now [channel, sample]
yields_per_channel_sample = np.sum(model_yields, axis=-1)
print(f"per channel and sample: {yields_per_channel_sample}")

# sum over samples (again, the innermost index as we dropped one dimension via the sum)
# remaining index is [channel]
yields_per_channel = np.sum(yields_per_channel_sample, axis=-1)
print(f"per channel: {yields_per_channel}")

output:

[[[24.0, 22.0]]]
per channel, sample, bin: [[[24.0, 22.0]]]
per channel and sample: [[46.]]
per channel: [46.]

You can skip any channels / samples you do not want to include and then sum over the dimension you want (-1 is the innermost one, see the numpy.sum documentation).

[ssaha1234]

Hi @alexander-held,
Thank you so much, this is exactly what I want. I implemented it in a slightly different manner using dictionaries, but the .sum method is much cleaner !
Best,
Shreya