Implementation of NeutrinoFlux container

Sheshuk commented 1 year ago

Closes #214

Initial tasks (before finalizing the interface)

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

Next tasks:

[x] Finalize interface
[x] Document the new classes
[x] Update SNEWPY documentation
[x] 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 commented 1 year ago

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 commented 1 year ago

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

Sheshuk commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

@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 commented 1 year ago

I have:

Implemented the checks which axes can be summed, and which - integrated.
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.
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.

JostMigenda commented 1 year ago

Sorry about the delay in reviewing this!

Thanks for the updates on the summation/integration checks! And your other points seem fine to me, too. This looks almost ready to merge, there’s just two small points remaining:

First of all, the integration tests are currently failing due to an off-by-one error in the array dimensions introduced in the last few commits; full traceback copied below. @Sheshuk, can you figure out what the issue is or should I take a look?

Full traceback

``` ERROR: test_simplerate (python.snewpy.test.simplerate_integrationtest.TestSimpleRate.test_simplerate) Integration test based on SNEWS2.0_rate_table_singleexample.py ---------------------------------------------------------------------- Traceback (most recent call last): File "/home/runner/work/snewpy/snewpy/python/snewpy/test/simplerate_integrationtest.py", line 32, in test_simplerate snowglobes.simulate(SNOwGLoBES_path, tarredfile, detector_input=detector, detector_effects=True) File "/opt/hostedtoolcache/Python/3.11.3/x64/lib/python3.11/site-packages/snewpy/snowglobes.py", line 351, in simulate res=sng.run(flux_files, det) ^^^^^^^^^^^^^^^^^^^^^^^^ File "/opt/hostedtoolcache/Python/3.11.3/x64/lib/python3.11/site-packages/snewpy/snowglobes_interface.py", line 307, in run result = self.compute_rates(detector, material, fluxes, energies) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "/opt/hostedtoolcache/Python/3.11.3/x64/lib/python3.11/site-packages/snewpy/snowglobes_interface.py", line 252, in compute_rates rates = np.dot(smear,rates) * effic ^^^^^^^^^^^^^^^^^^^ File "<__array_function__ internals>", line 200, in dot ValueError: shapes (200,200) and (199,) not aligned: 200 (dim 1) != 199 (dim 0) ```

And second, as discussed previously (see point 6), we should first use and carefully test these container classes ourselves before we make them part of the public API. Some of that will be in #251; and additionally in the pre-SN implementation. Until then, I would either omit these container classes from the public documentation completely (or at least add a warning box at the top of flux.rst that this API is still unstable).

Sheshuk commented 1 year ago

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 ee69beb13e64362e30f49c770b0d62011bb1e57d - 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 commented 1 year ago

@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 https://github.com/SNEWS2/snewpy/pull/236/commits/a42ac5438827d86756473bf0b4d707eaaf2b7e8a

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 commented 1 year ago

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
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
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 commented 1 year ago

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

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

JostMigenda commented 1 year ago

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

SNEWS2 / snewpy