IRF classes

[maxnoe]

Right now, we do not have special datastructures for IRFs, we use astropy quantities, numpy arrays.

This makes functions simpler and more flexible but an advantage of a specific datastructure would be to keep information together that belongs together (irfs with their bins).

Since gammapy already has these classes, we could think about using them, but that would again mean adding gammapy as required depedency.

An alternative would be to use some thin data class that holds the values, binnings and metadata.

[kosack]

I'd suggest a light internal data structure, since the gammapy's IRF structure may not be stable and it has some assumptions on how IRFs are dimensioned that we might not want to make initially in pyirf. For example, they assume IRFs are 2D or 3D and are already marginalized to a specific observation (which is fine for a DL3 IRF), but in pyirf might want to be able to make IRFs as a function of not only fov coordinates + energy, but also zenith angle, or other dimensions (a master IRF that can then be later marginalized to many observations).

[kosack]

could look at xarray for inspiration for the API at least. It may be overkill as an explicit dependency, though it does give you interpolation, plotting, and IO (to hdf5 netcdf) built-in (http://xarray.pydata.org/en/stable/data-structures.html).

[maxnoe]

The idea is actually pretty similar to the gammapy Maps and MapAxis classes, which are used for the IRF classes internally. We could use those with extra dimensions.

[adonath]

@maxnoe and @kosack In Gammapy we have (for "historical reasons") two nd data structures: IRF and Map, the Map class requires a sky coordinate axis first such as WCS, HEALPix or a region for spectra and has additional physical coordinate axes, the IRF class has just the physical coordinates axis. The assumption for both is that the physical coordinate axes are independent.

This year we spend some time to refactor the existing IRF classes in Gammapy and introduced a ND base class. The class takes the declaration of the axis coordinates as a list of MapAxis objects, data, unit and meta data. It offers e.g. interpolation, integration, cumsum, normalisation methods etc. all supported along specified axis names. In general the mid-term plan for Gammapy is to further unify the current IRF class and the Map class. So that other methods such as .upsample(), .downsample(), resample_axis() are supported as well for the IRF class (currently only supported for Map...) and can possibly moved to a common base class.

Using this base class the actual implementations of the IRF classes with fixed dimensions simplify to just declaring the physical axes and order and implement specific plotting functionality see e.g. EffectiveAreaTable2D. We also decided to keep the previously existing public API, that's why the base classes are not yet publicly documented.

In general I think a package like xarray is a good choice for such data structures however currently there are some limitations, which makes it non-ideal for our use-cases. Here are a few additional thoughts I have on this:

• I think currently there is no real support for coordinates "edges" or "boundaries" in xarray, but only coordinates at the bin center. When working with "integrated" or "undersampled" data it is useful to have the edge information. Together with some assumption on the spacing of the coordinate array and shape of the data inside the bin one can really improve the precision of interpolation, by using e.g. relying on log-log interpolation.
• The same goes e.g. for integration, where on can rely on the interpolation of the cumulative sum, if the bin edges are known. In the Gammapy MapAxis object we store a bit of meta data in this regard, e.g. the interp={"linear", "log", "sqrt"} option.
• Most of the IRF class code in Gammapy is simple, however there are a few non-trivial things, that are not supported by xarray, which are useful for pyirf as well: e.g. computation of PSF containment radii with broadcasting support see here or interpolation of missing values using the aforementioned scaling along a certain axis e.g. interpolate log-log in energy.
• I guess in general it make sense to share some functionality between Gammapy and pyirf. I think an obvious case might be validation of FITS data (basically what @maxnoe implemented in https://github.com/open-gamma-ray-astro/ogadf-schema/tree/master/ogadf_schema) or visualisation of IRFs. This is both required from the users as well as IRF production side. The same might be useful for the IRF data structures in general.

[kosack]

@adonath If you plan to move IRFs to use Maps, then that sounds like a good solution. It solves the dimensionality problem (assuming the spatial coordinates can be in AltAZ) and would be a nice interface. So I would suggest we just start with using gammapy.Maps. The only reason I originally hesitated to suggest that option was the uncertainty in the Science Tools selection, but now it's clear having gammapy as a dependency is not a big issue.