Interpolation API

[joleroi]

@TarekHC

I added an interpolator to the ND data scratch class. Let's us this to discuss/improve it. At the moment you have to add an interpolator manually to the array by calling

nddata.add_linear_interpolator()

this creates a scipy RegularGridInterpolator

http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.RegularGridInterpolator.html

then you can evaluate it using the evaluate method, which has a keyword for the interpolation method. At the moment 'linear' and 'nearest' (these 2 are covered by the scipy RegularGridInterpolator).
The idea was to add more methods to this evaluate function, that call into other interpolators (if we want other methods at all). So you would for example have

hist.add_spline_interpolator()
hist.evaluate(..., method=spline)

An example how the interpolation works is also in https://github.com/joleroi/nddata/blob/master/NDData_demo.ipynb

Do you think that's a good idea?
Do you think we should add the regular grid interpolator by default?

cc @cdeil

[TarekHC]

Hi @joleroi ,

Yes, after going through the link you provide, this RegularGridInterpolator seems like a good option (indeed, you could set it by default). I'm not very involved with the analysis tools right now, but the interpolator needs to be quite fast. Maybe it would be a good time to test this implementation with the benchmarks @cfeil suggested.

The interpolation scheme keyword will be associated with each IRF element (effective area, PSF, etc..) or with each axis? Would it be possible to have the keyword define the interpolation method applied to each axis of the multidimensional array? Not sure if it would be possible to implement...

Unfortunately, until the end of the CTA Kashiwa meeting (20/5) I won't have much time to contribute... Sorry about that! But I promise as soon as I have some time I will try to get more involved.

[cdeil]

@joleroi - Thanks for continuing with this.

I think it's probably impossible to to get this right now (even if we have thought about this a lot and prototyped), because we don't know what we want / need in the final product in terms of interpolation.

Here's just some of the questions I have:

• How will this class be used by concrete IRFs, like the AEFF, EDISP, PSF and BKG (and maybe other non-IRF classes like SkyCube) that we have?
• Do we have interpolators where __init__ is compute-intensive? Or can interpolator construction (and thus configuration) always happen in evaluate?
• How flexible do we want to be? Give the user a fixed set of pre-implemented ways how to do interpolation (with some configuration options) or design the system to be extensible (e.g. user can add their own interpolator to NDData)?
• Do we support different interpolation methods along different axes? (E.g. nearest-neighbor in energy, spline in offset and migra for EDISP)?
• Do the interpolators need to support axes that integrate exactly to 1 (e.g. for the PSF or EDISP)? Where is this interpolation implemented so that it's efficient and precise and knows about special axis scales (e.g. use log(flux)-log(energy) trapezoidal integration)? I.e. is there some axis-specific interpolation info / configuration or is it global on nddata?

So my suggestions would be to use a development method that's more focused on concrete use cases and incremental improvements / refactoring:

• Don't try to get interpolation right now.
• Instead start with what you have, move it into gammapy.irf very soon and try to use it AEFF, a simple IRF where no integration is needed. This will already be a big job (a week?), because you need to figure out how this class and the AEFF class are related (inheritance? mix-in? configuration?) and how FITS I/O should work and implement FITS I/O and add more tests.
• Then try to use this for an IRF class with axes that have to integrate to 1, e.g. tabulated EDISP or tabulated PSF.
• Then gammapy.irf is ready for 1.0.
• Then we can start thinking about how to make it more flexible (again to support concrete use cases, like e.g. a cubic spline axis along MIGRA.

That said ... your particular API here of making the user remember to call nddata.add_linear_interpolator is not very nice. Since RegularGridInterpolator.__init__ doesn't compute anything (or at least it's very fast compared to evaluate), you can just instantiate it in the nddata.evaluate call or in nddata.__init__. If you had interpolators where __init__ takes time (which you don't at the moment) you could delay their creation to first use, by making them properties nddata.spline_interpolator that call the interpolator constructor on first access and then cache the object. But it might be non-trivial to configure the interpolator __init__ in this case, i.e. expose this is a simple-to use, yet flexible way.

[joleroi]

@cdeil thanks for your comments. Since interpolation and FITS I/O are actually implemented and tested I agree that the next step is to start with a concrete IRF example. I will make a PR to gammapy next week. Let's discuss detailed questions then.