Interfaces and conventions for modifiers?

[lhenkelm]

I hacked together a custom modifier for a polynomial "shapefactor" using the features recently added to the pyhf master branch by @lukasheinrich (thanks a lot for implementing those!).
After a bit of trial and error, this polynomial modifier seems to work now,
but I am worried about breaking some assumptions of how a <modifier>_builder and <modifier>_combined should behave,
so I am asking about them here to avoid inadvertantly introducing silent bugs.

A lot of the interfaces will likely soon be spelled out in docstrings once #1641 is implemented, so I try to limit myself to specific points that I find hard to guess from skimming through the code:

• required parameter sets:
  • what do the "is_shared" and "paramset_type" fields mean?
  • are the "bounds" only used inside optimisation, or are they ever used to check/clip input parameters?
• Does it make sense to enforce that the same (by key) modifier should always come with the same modifier data?
• <modifier>_builders:
  • is there an expected interface of the collect(...) method? Right now it is only used inside the respective <modifier>_builder, but it seems to be present in all builders, so do you anticipate external use? Should I provide one in custom modifier builders to be futureproof?
  • What is the intent behind accessing builder_data through finalize() rather than directly?
  • What is the best way to pass channel-specfic data through builder_data? E.g. I had to pass bin centers, which logically is channel-associated information. But to avoid potential for name collisions with particular sample names, I had to make that a new entry in every affected modifier-key and sample, (even if they all share a channel, so should share a binning)
• <modifier>_combined:
  • Is the <modifier>_combined.name class attribute the type name specified in a spec? Where is it used?
    • Why is it not used in histfactory_set?
  • <modifier>_combined.apply(pars): Does the shape of the modifications tensor match the shape of the tensor being modified by it, or is there implicit broadcasting going on?

Finally I have similar questions not about modifiers, but rather some (nearly) undocumented pyhf internals:

• tensorlib shim: Why is pyhf.default_backend its own name alongside pyhf.tensorlib? When and why should I use the former rather than the latter?
• ParamViewer:
  • What is the intended meaning of the arguments to ParamViewer.init?
  • What is ParamViewer.index_selection?
• Not a question, but docstrings for ParamViewer and _TensorViewer (and its associated factory methods) would be very useful.

Sorry about the length! Answers to only one or some of the question will still be quite helpful!

[kratsg]

• what do the "is_shared" and "paramset_type" fields mean?

is_shared was supposed to keep track of when a modifier could be correlated or not... but this got lost due to time and it seems we have a bug. This will be tracked in #1908. Here's the relevant code for the meaning of is_shared:

pyhf/src/pyhf/pdf.py

Lines 50 to 63 in 121903b

	# check the shareability (e.g. for shapesys for example)
	is_shared = paramset_requirements['is_shared']
	if not (is_shared) and modifier_def['name'] in _paramsets_requirements:
	raise ValueError(
	"Trying to add unshared-paramset but other paramsets exist with the same name."
	)
	if is_shared and not (
	_paramsets_requirements.get(
	modifier_def['name'], [{'is_shared': True}]
	)[0]['is_shared']
	):
	raise ValueError(
	"Trying to add shared-paramset but other paramset of same name is indicated to be unshared."
	)

• are the "bounds" only used inside optimisation, or are they ever used to check/clip input parameters?

That partially depends on the behavior of the optimizer itself. We don't check input parameters based on the bounds in pyhf. scipy used to have a long-standing bug where it did allow the parameter to go outside of the bounds which was only fixed about 8 months ago - and this now spits out a warning of "parameter outside of bounds - clipping" or something similar. Minuit is a little bit more robust in this case. We don't use the bounds anywhere outside of optimization afaik.

• Does it make sense to enforce that the same (by key) modifier should always come with the same modifier data?

Not sure what you mean "by key" here. The fundamental point is that two modifiers of the same type can exist on a sample if they have different names. If they have the same name, we take the last defined modifier type/name (#1899) which is a bug and needs to be fixed. If two modifiers of the same type have different names on the same sample, that's also fine. But there is a relatively fundamental assumption that you can't have multiple modifiers of the same type and name on the same sample.

• is there an expected interface of the collect(...) method? Right now it is only used inside the respective <modifier>_builder, but it seems to be present in all builders, so do you anticipate external use? Should I provide one in custom modifier builders to be futureproof?

Yes please. Try to be as consistent in the API as possible. We'll end up trying to add typing to all of this, and likely we'll need to add collect as an expected interface, even if it's pass-thru.

• What is the intent behind accessing builder_data through finalize() rather than directly?

Directly meaning what? Lukas can comment better, but we need to sometimes modify or concatenate after the fact after we build up combined modifiers. Until you call finalize() the builder_data is not ready. We could make finalize() pass-through and return None by default, and then you can always ask for the combined_mod.builder_data instead. This was a choice, but it's probably useful to refactor to be cleaner as part of #1641.

• What is the best way to pass channel-specfic data through builder_data? E.g. I had to pass bin centers, which logically is channel-associated information. But to avoid potential for name collisions with particular sample names, I had to make that a new entry in every affected modifier-key and sample, (even if they all share a channel, so should share a binning)

Why do you need bin-centers? HistFactory doesn't care about the actual binning/bin edges, just that you have n-bins. More generally, not sure what the best way is. We'll have to see what you want to do and approach it on a case-by-case basis... definitely avoid globals though.

• Is the <modifier>_combined.name class attribute the type name specified in a spec? Where is it used?

Yes, it's a bit convoluted. It's needed here

pyhf/src/pyhf/pdf.py

Lines 543 to 547 in acde7f4

	for modifier_applier in self.modifiers_appliers.values():
	if modifier_applier.op_code == "addition":
	self._delta_mods.append(modifier_applier.name)
	elif modifier_applier.op_code == "multiplication":
	self._factor_mods.append(modifier_applier.name)

and is defined in

pyhf/src/pyhf/modifiers/__init__.py

Lines 39 to 47 in acde7f4

	histfactory_set = {
	"histosys": (histosys_builder, histosys_combined),
	"lumi": (lumi_builder, lumi_combined),
	"normfactor": (normfactor_builder, normfactor_combined),
	"normsys": (normsys_builder, normsys_combined),
	"shapefactor": (shapefactor_builder, shapefactor_combined),
	"shapesys": (shapesys_builder, shapesys_combined),
	"staterror": (staterror_builder, staterror_combined),
	}

for the "spec". The idea is that we pick up by "name" and the combined is just a unique name since there should only ever be one combined modifier for that name.

* Why is it not used in [`histfactory_set`](https://github.com/scikit-hep/pyhf/blob/b798ef23c77db6b1d6d606f4b360a5f47fd63dd2/src/pyhf/modifiers/__init__.py#L39)?

It should. Needs to be fixed in #1641.

• <modifier>_combined.apply(pars): Does the shape of the modifications tensor match the shape of the tensor being modified by it, or is there implicit broadcasting going on?

There can be broadcasting depending on the batch size as well. pars is more like (nbatch,npars) in terms of the shape, which is typically (1,npars) by default if there's no batching.

• tensorlib shim: Why is pyhf.default_backend its own name alongside pyhf.tensorlib? When and why should I use the former rather than the latter?

pyhf.tensorlib will always refer to the current tensorlib being used at the time of import. This is useful for users, terrible for developers.

From a user perspective:

>>> import pyhf
>>> from pyhf import tensorlib
>>> pyhf.set_backend('jax')
>>> tensorlib
<pyhf.tensor.numpy_backend.numpy_backend object at 0x10f96f300>
>>> pyhf.tensorlib
<pyhf.tensor.jax_backend.jax_backend object at 0x10f9633c0>

Note that users won't explicitly import tensorlib but if they do, it's imported as a cached object. The problem for us as developers is that typically, we do something like from . import tensorlib which keeps it cached at import time, but it might have changed by the time it gets used in the code because the user would've changed the backend. This is why you commonly see a lot of tensorlib, _ = get_backend() because that forces us to make sure we're not using a cache'd version.

Now, default_backend is not always the same, and will be the default backend (typically numpy). Sometimes, e.g. in the case for auto-diff, the default backend needs to switch to jax because you want to allow for end-to-end auto-diff-ability, and that is another complication. The idea behind default backend is that it's typically something that you can load/translate into the currently-defined backend "easily" or trivially. In the case of numpy, this is a guarantee. In the case of jax or something else as default backend, then the user is more expert.

• What is the intended meaning of the arguments to ParamViewer.init?

This is when you may or may not have batched input (see my answer above) so you initialize a ParamViewer to allow you to extract the parameters in a potentially-batched way without caring if you're batched or not.

• What is ParamViewer.index_selection?

The indices in the potentially batched pars that are relevant for the initialized ParamViewer object (think of it like an actual view of an array, but not quite). It's a bunch of bookkeeping.

• Not a question, but docstrings for ParamViewer and _TensorViewer (and its associated factory methods) would be very useful.

Yeah. Can you file an issue for that?

[lhenkelm]

Thanks a bunch for the explanations! Some follow-up questions below

• is there an expected interface of the collect(...) method? Right now it is only used inside the respective _builder, but it seems to be present in all builders, so do you anticipate external use? Should I provide one in custom modifier builders to be futureproof?

Yes please. Try to be as consistent in the API as possible. We'll end up trying to add typing to all of this, and likely we'll need to add collect as an expected interface, even if it's pass-thru.

Thanks for the guidance here. Could you specify (or is it maybe already spelled out somewhere?) what the expected behavior of a general <modifier>_builder.collect method is? In the one built-in class that I think I understand (the shapefactor_builder) its just a means of preparing the contents of the masking tensor, but the way it is implemented (taking lots more data as inputs than needed for the mask, returning a dict with the mask-list rather than the list itself) makes me think it is just a special case of a more general pattern. It would be very helpful to have that general behavior stated explicitly.

• What is the best way to pass channel-specfic data through builder_data? E.g. I had to pass bin centers, which logically is channel-associated information. But to avoid potential for name collisions with particular sample names, I had to make that a new entry in every affected modifier-key and sample, (even if they all share a channel, so should share a binning)

Why do you need bin-centers? HistFactory doesn't care about the actual binning/bin edges, just that you have n-bins. More generally, not sure what the best way is. We'll have to see what you want to do and approach it on a case-by-case basis... definitely avoid globals though.

E.g. if the custom modifier expresses a multiplication of the nominal sample yield by some parametric function of whatever observable the bins in the channel are binned in.

In my particular example, I am modifying the yield of a sample binned in the mass of some resonance candidate jet by a truncated polynomial series in the mass. The polynomial evaliuation takes the polynomial series coefficients (which are part of the model parameters) and the points at which to evaluate the polynomial as inputs. Different bins have different widths, so I need the bin centers to get sensible, reproducible results that don't change in surprising ways when the binning is adjusted.

The way I am handling this right now is by having every JSON spec for a custom polynomial modifier pass the bin centers along as part of it's data attribute,
which the <...>_builder then checks for consistency if two modifiers (of different names and/or samples) share a channel.
It is a bit of an awkward duplication, but it keeps the disruption to the schema minimal (it is contained in the new modifier/its "data" field in the JSON spec. )

Now, default_backend is not always the same, and will be the default backend (typically numpy). Sometimes, e.g. in the case for auto-diff, the default backend needs to switch to jax because you want to allow for end-to-end auto-diff-ability, and that is another complication. The idea behind default backend is that it's typically something that you can load/translate into the currently-defined backend "easily" or trivially.

So, in practical terms I am not doing anything wrong if I ignore the "default_tensorlib_changed" event/ the default_backend, and just consistently refer to tb, _ = get_backend() and subscribe to the "tensorlib_changed" event in a custom modifier as needed, right?

I am not sure how relevant this is (it seems right now there is no code in master that is actually subscribing to the "default_tensorlib_changed" and "default_optimizer_changed" events, I suppose that will change?),
but why would someone want to switch the backend yet leave the default_backend at some old value?
I am probably missing something very obvious, but also for your example of auto-diff, would it not be best to have every tensor consistently switch to the format expected/supported by the auto-diff implementation? Or is it valuable for the auto-diff-ing user to (in a programmable way) find out what the default backend used to be?

[kratsg]

Thanks for the guidance here. Could you specify (or is it maybe already spelled out somewhere?) what the expected behavior of a general <modifier>_builder.collect method is? In the one built-in class that I think I understand (the shapefactor_builder) its just a means of preparing the contents of the masking tensor, but the way it is implemented (taking lots more data as inputs than needed for the mask, returning a dict with the mask-list rather than the list itself) makes me think it is just a special case of a more general pattern. It would be very helpful to have that general behavior stated explicitly.

We haven't figured out the general pattern yet. This is likely something that will end up inspiring or being refactored for the general "open-world" implementation (basically a pythonic version of RooFit/RooStats).

The way you're doing it is fine as long as pyhf is able to handle it via the custom_mods approach. We're only opinionated about certain things within reason.

So, in practical terms I am not doing anything wrong if I ignore the "default_tensorlib_changed" event/ the default_backend, and just consistently refer to tb, _ = get_backend() and subscribe to the "tensorlib_changed" event in a custom modifier as needed, right?

Yes. No reason to subscribe to it if you don't need to listen to it.

I am not sure how relevant this is (it seems right now there is no code in master that is actually subscribing to the "default_tensorlib_changed" and "default_optimizer_changed" events, I suppose that will change?),
but why would someone want to switch the backend yet leave the default_backend at some old value?
I am probably missing something very obvious, but also for your example of auto-diff, would it not be best to have every tensor consistently switch to the format expected/supported by the auto-diff implementation? Or is it valuable for the auto-diff-ing user to (in a programmable way) find out what the default backend used to be?

Some implementations on top of pyhf might need it. The default backend is used for calculating a cache set of values, so that when one switches the backend, only a subset of the calculations need to be re-done on the fly. It's a niche case that mostly comes up related to autodifferentiability.