Implementation of NeutrinoFlux container

[Sheshuk]

Closes #214

Initial tasks (before finalizing the interface)

• Implement a flux/fluence/etc container class
• SupernovaModel.get_flux should produce Flux class object
• Implement rate calculation using this class (passing Fluence to SimpleRate or similar class)
• Implement save/load functions for this class (to do the caching, if desired)

Next tasks:

• Finalize interface
• Document the new classes
• Update SNEWPY documentation
• Add tests

---

Postponed tasks

Implement Flux.apply(flavor_xform:FlavorTransformation) method, to apply the transformation, and produce a transformed flux

Postponed till #242 is implemented, right now SupernovaModel.get_flux gets a FlavorTransformation argument to have transformed flux)

Merge RateCalculator and SimpleRate classes

I separated the common part as SnowglobesData base class, which handles the reading and accessing stuff from the snowglobes. So I think there is no need for this merge

[Sheshuk]

Applying FlavorTransformation to the flux should wait until we implement #242.
At the moment, one can get the transformed flux via method SupernovaModel.get_flux providing the flavor_xform parameter there.

[Sheshuk]

I've added a notebook describing the usage and main features of the new flux.Container class.

[Sheshuk]

The base functionality is implemented, and before I continue, I'd like to have opinions on the interface I propose.
Please check the jupyter notebook with the usage examples of the implementation, and the code as well.
Any suggestions are welcome, please let me know if I'm moving to the right direction or not.

• For the rate computation: I separated new functionality in the RateCalculator subclass of SimpleRate. It is done temporarily for clarity (to avoid having many run and compute_rates methods in one class), we'll probably want to merge these classes into one later.
• In the rate computation, I used the integration over energy, instead of using bin-center approximation, so it should be slightly more precise.
  I still implemented the old method as RateCalculator.run_simple
• A question: SimpleRate provides us with four tables of event rates: smeared,unsmeared,weighted,unweighted. Do we ever need to use unweighted variant?
  It looks like an intermediate calculation step, before multiplication by SNOwGLoBES weights (for each channel), so maybe it doesn't have any sense for the final user?

[JostMigenda]

Thanks for the great work @Sheshuk! The notebook is a nice illustration of the functionality and I think the overall direction is good.

There are a few conceptual issues I have with the Container.sum and Container.integrate functions.
Regarding Container.sum:

• Summing along the flavor axis is okay; though the use cases for that are fairly limited, I think.
• Summing along the time or energy axis of a Flux is physically meaningless: It gives results that depend on the binning in the original model file—especially so, since several model files have non-uniform binning, e.g. very narrow sub-ms time bins during the accretion phase and then wider time bins (up to ~100 ms) during the neutrino cooling phase. Simply summing those different bins without accounting for bin widths gives completely arbitrary numbers.
• Summing over the time axis of a Fluence or the energy axis of a Spectrum is meaningful and may be useful iff the integration that created that Fluence/Spectrum used >2 entries in limits.

Regarding Container.integrate:

• Integrating over the flavor axis doesn’t really make physical sense; and unfortunately, numerically it’s not the same as summing over the same flavours.
• Integrating along an axis that was already integrated should probably return an error.

I suspect all of these could be fixed by careful checks at the start of the sum and integrate functions; but I haven’t fully thought through the implementation details.

Other than that, I’ve suggested a few minor changes to the code inline.
Oh, and especially nice work on the ContainerClass helper function—that’s a very elegant way to solve this issue! 🤩

[JostMigenda]

• For the rate computation: I separated new functionality in the RateCalculator subclass of SimpleRate. It is done temporarily for clarity (to avoid having many run and compute_rates methods in one class), we'll probably want to merge these classes into one later.

For now, that’s definitely a good idea. Since this is intended to be temporary, I haven’t looked at the RateCalculator class to closely in my review above.

Before we can merge RateCalculator and SimpleRate, we’ll have to see whether we can update the implementation of snewpy.snowglobes.calculate to use this new flux class; but let’s leave that aside for now.

• In the rate computation, I used the integration over energy, instead of using bin-center approximation, so it should be slightly more precise.
  I still implemented the old method as RateCalculator.run_simple

Since integrate simply uses trapezoidal integration, does that really make a difference? I think it only matters at the edges of your integration region and only if those don’t align with bin centres, right?

• A question: SimpleRate provides us with four tables of event rates: smeared,unsmeared,weighted,unweighted. Do we ever need to use unweighted variant?
  It looks like an intermediate calculation step, before multiplication by SNOwGLoBES weights (for each channel), so maybe it doesn't have any sense for the final user?

Agree; the unweighted numbers are an intermediate step and I don’t see any use case for exposing them to the user.

[Sheshuk]

@JostMigenda thank you for your review and suggestions!

About the integration/summation - I totally agree, in some cases it doesn't have physical sense. I thought that it's up to a user to choose physically meaningful operations, but having checks would be nice. I need to think about how to implement it. Probably keeping a set of integrable axes in the class?

Since integrate simply uses trapezoidal integration, does that really make a difference? I think it only matters at the edges of your integration region and only if those don’t align with bin centres, right?

You're correct - if the integration bins are the same as the initial energy sampling points. But the sampling points can be much more detailed - this is especially the case I want for preSN models - while the integration bins are defined by the SNOwGLoBES smearing matrix bins.
In this case sampling cross-sections in the same energy points as flux, and integrating these in the needed energy bins give a more precise result.
I don't know if I describe it clearly.

[Sheshuk]

I have:

1. Implemented the checks which axes can be summed, and which - integrated.
2. Separated the reading of SNOwGLoBES data files from SimpleRate to SnowglobesData class.
   Both RateCalculator and SimpleRate classes inherit from the SnowglobesData, and implement just the rate calculation in their own ways.
   So I would prefer not to merge RateCalculator and SimpleRate, and just keep the calculator in snewpy.rate_calculator.RateCalculator, and snewpy.snowglobes_interface.SnowglobesData separately, because the latter does not need to be exposed to user.
3. Postponed implementing the Container.apply(transformation).

Before we can merge RateCalculator and SimpleRate, we’ll have to see whether we can update the implementation of snewpy.snowglobes.calculate to use this new flux class; but let’s leave that aside for now.

So now most other tasks were finished or postponed, I think I'm ready to work on this.

[Sheshuk]

Thanks for looking into this!

the integration tests are currently failing due to an off-by-one error in the array dimensions

Right, it was introduced in ee69beb - it looks like we were treating the SNOWGLOBES binning incorrectly all this time, but the tests are still for the old version.
I will try to revert it in this branch, so the tests pass, and will probably have to reimplement it again in #251

[JostMigenda]

@Sheshuk and I have talked about this PR and #251 earlier today and fixed some remaining bugs. Since both PRs are closely related, it makes sense to first merge #251 (where these SNOwGLoBES binning issues & related tests are both fixed) into this one, and then merge this PR here into the main branch.

For the record, the final commit on this PR, before the first merge, was a42ac54

And aside from the functional reasons for this PR, which we’ve discussed in detail, I’d just like to note that there are some nice performance benefits as well: In my testing using the example in SNOwGLoBES_usage.ipynb (Zha_2017_s17 progenitor, AdiabaticMSW transformation, 60 time bins between 0.742 and 0.762 s), the total time for generate_fluence, simulate and collate went from ~35 s to ~10 s. ⚡

[JostMigenda]

Some more issues with the integration tests. The first one was because of a tiny change in the keys for the results dictionary; that’s fixed now.

The current issue requires a bit more thought, because this PR changes the binning and that leads to slightly different event counts:

Unsmeared:

Channel	Before	After	Change
ibd	12498.485	12506.335	0.6 ‰
nue_O16	110.328	110.435	1.0 ‰
nuebar_O16	137.608	137.741	1.0 ‰
nc	275.077	275.397	1.2 ‰
e	1006.425	1006.915	0.5 ‰
Total in Super-K¹	4488.94	4491.78	0.6‰

The unsmeared numbers typically change by less than 1 per mille; I think that’s nothing to worry about.

Smeared:

Channel	Before	After	Change
ibd	11916.581	11969.904	4.5 ‰
nue_O16	109.625	109.772	1.3 ‰
nuebar_O16	135.467	135.737	2.0 ‰
nc	28.176	28.208	1.1 ‰
e	455.939	461.574	12.2 ‰
Total in Super-K¹	4046.65	4065.66	4.7 ‰

After smearing², the changes are much bigger, mainly in the electron scattering and IBD channels. In water Cherenkov detectors, those are the two channels which have the most events near the detector threshold (~5 MeV or so, depending on the detector); so any changes to the binning that cause events to slip above/below the detector threshold would predominantly affect those two channels.

So I think these changes are within what we’d expect; but maybe we should wait a little to let people comment? I’ll also try to double-check the binning code tomorrow, just to feel confident I didn’t miss anything.

¹ This is the sum of the rows above, multiplied by 0.32 (because the SNOwGLoBES detector config has 100 kt volume, while SK’s inner detector has 32 kt).
² smearing and detector efficiency, to be precise

[evanoconnor]

@JostMigenda Are the percentages off by a factor of 10?

edit: I see, per mil, a weird unit …

[JostMigenda]

@evanoconnor Changes are all in per mille, not per cent.

[JostMigenda]

I did go through the rate calculation code once again; thanks to @Sheshuk for pointing me to pp. 100&105 of the GLoBES manual for more on the energy binning and smearing files.
It looks good to me and I would be happy to merge this in a moment; I’ll just need to submit a quick commit to fix the number of events required in the integration tests first.