Deal with undefined behavior around spectral norm in RspConverter#463
Deal with undefined behavior around spectral norm in RspConverter#463jdbuhler wants to merge 9 commits intocositools:developfrom
Conversation
any Ei bin when computing effective area correction for response. * Don't allow inf or nan in effective area in the above case -- set the correction to zero for these bins. * Actually honor emin and emax in Linear normalization, as it is honored for the powerlaw case * Consistently strip units when computing effective area correction * Try to correct (currently untestable) validation code for gaussian spectral normalization
Codecov Report❌ Patch coverage is
🚀 New features to boost your workflow:
|
|
Hi @jdbuhler , could you try to test it with this : If it is working, add this file to the |
|
This is a sparse response, which is not supported by RspConverter. (It should cleanly report this fact based on the header, but it is crashing instead, which I'll need to fix.) I believe we decided when the RspConverter patch was first accepted that sparse responses had been deprecated. |
|
Also, a follow-up question: this file has Mono normalization, but it specifies a parameter (511) to the norm. The RspConverter code (which follows what was in cospy 0.3.x) believes that Mono normalization does not take any parameters, and my parameter parsing code consequently throws an error. What does the parameter to Mono mean? I'm guessing it means a mono energetic beam (in this case, 511 keV), in which case the existing normalization code is totally broken for this case unless the Ei axis has only a single bin that includes the target energy. |
it's a valid file * slightly more robust stripping of extensions from input filename (but we still don't support uncompressed .rsp, though perhaps we should)
both reading and writing. * Add test cases for uncompressed read/write * FDR error message should not mention .rsp, as it only opens .h5 files
|
@jdbuhler Ah right I forgot about the sparse thing sorry. And yes also this is not a gaussian but a monoenergetic beam my bad. Yes the Mono beam has a parameter telling which energy is used but we don't need it for any normalization since we use it for only one single Ei bin (for 511, 1809,.. response) . So all the events ends up in this bin. I will try to simulate a small response with a gaussian beam. Note that we never used it actually |
mono-energetic energy. If not specified, normalization is ill-defined if Ei axis has multiple bins. If specified, normalization is zero for bins not containing the energy (which is weird and throws a warning, but at least is well-defined and should not crash). * test case for Mono in fact used an Ei axis with multiple bins but did not specify an energy! Fix to provide an energy.
without a suitable response file
|
Thanks, Savitri! I would definitely like a small .rsp file with only a single Ei bin. That should be enough to exercise the Gaussian code (and also the Mono code that I can't otherwise test right now). |
|
I'm not sure we actually need a simulation done with a Gaussian beam for testing purposes. |
|
Ok then let just have a test with the mono energetic one. This one is not sparse : https://seafile.rlp.net/seafhttp/f/54540b5547d9411b92db/?op=view |
|
OK, update: I took @GallegoSav 's latest monoenergetic response, reduced it to a workable size for a test case, and removed the norm from the header so I could try whatever I want. It does indeed barf as expected with the Gaussian code that is there now, because of the infinite support issue. However, I note that Gaussian normalization is supposed to take 3 parameters (the third is a "cutoff"), but the test only looks at two. I'm not sure what the cutoff is, but I'm guessing it might be how many std deviations to go out from the mean, so that the support is finite. So the question now is, how do we want to address this? Reasonable options I can think of include: (1) Declare that nobody actually cares about Gaussian normalization anymore, and remove it entirely. (2) Clarify the meaning of the cutoff parameter, which would tell me what the range of energies is for which the normalization is nonzero (the equivalent of emin and emax for linear or powerlaw normalization), and use that to set the bin weighting for effective area. This is well-defined even if the input contains multiple Ei bins. (3) Replace the cutoff parameter with an explicit emin and emax, and again use that to set the bin weighting. I note that for a response with a single Ei bin, the normalization for this bin isn't necessarily going to be 1 under schemes 2 or 3, which is intentional -- some of the incident photons are expected to fall outside the range of energies represented on the Ei axis. Do you have views on how we should proceed? @israelmcmc , do you want to weigh in as well? |
* Gaussian normalization takes 3 arguments
|
@jdbuhler Looking at the Gaussian function definition in cosima documentation, it show that indeed the third parameter is the cutoff in number of sigma . It also shows the energy range as function of the parameters. |
|
Ok, that is sufficiently well-defined that I can implement it and remove the special-casing for Gaussian. I'll work on that. |
|
It took me a while for get to this one, but it seems everything has been figured out. Thanks @jdbuhler and @GallegoSav! |
|
Should be ready to go in a couple of hours at most (I am adding a test case for @GallegoSav 's mono response) |
…energy * make sure we use an empty norm parameter list, not [None], to describe a Mono norm without energy internally, so that we don't emit 'None' in the saved response headers
|
OK, I think this is done. @GallegoSav, can you please review? |
3 participants
The RspConverter must compute an effective area correction when converting a response from MegaLib to HDF5 format. This correction depends on the energy spectrum from which incident photons were drawn when performing the simulations that produced the response matrix. That spectrum (the spectral norm) is specified either in the .rsp file or as an argument to the RspConverter's init function.
It is possible to specify a norm that implies that some Ei bins receive zero incident photons. In such a case, the effective area correction is undefined for those bins. This patch logs a warning when this case happens and sets eff_area for the affected bins to zero, rather than to infinity or some other invalid result.
The behavior of some of the norm options could be unexpected or undefined. The new code enforces the emin and emax parameters of the spectral norm in the Linear case (the powerlaw case already enforced these) and acts sensibly for the Mono norm whether or not an energy is provided as an argument. It also adds support for computing an effective area correction for the Gaussian norm, wehther or not the input contains multiple Ei bins.
Notably, many of the test cases for the RspConverter were using spectral norms that did not cover the entire range of Ei, which before this patch resulted in divide-by-zero warnings from Numpy while running the test suite. Moreover, the Mono test case did not specify an energy yet was applied to a response with multiple Ei bins. These cases have now been fixed. Thanks to Savitri Gallego for providing a test response with a single Ei bin, which has been added to test_data.
In responding to a bug report, I got annoyed that the RspConverter only reads compressed .rsp.gz files. I added support for reading and writing either compressed or uncompressed .rsp files.