Change NDDataArray default_interp_kwargs to extrapolate

[cdeil]

@joleroi - while working on CTA IRFs, @facero and @leajouvin just now had the problem that our IRFs aren't extrapolating. I think we discussed this last week and said we should change to extrapolate by default.

Currently we have this:
https://github.com/gammapy/gammapy/blame/0e3a282736308a67fdcd51cbfded123b339e4770/gammapy/utils/nddata.py#L40

default_interp_kwargs = dict(bounds_error=False)

According to the scipy docs, it should be this to extrapolate, no?

default_interp_kwargs = dict(bounds_error=False, fill_value=None)

For now, they can use this and continue with their work:

from gammapy.utils.nddata import NDDataArray
NDDataArray.default_interp_kwargs = dict(bounds_error=False, fill_value=None)

@joleroi - If you agree, could you please make the change and add one test to make sure it's extrapolating?

[joleroi]

According to the scipy docs, it should be this to extrapolate, no?

I know that you don't remember, but you wanted to change this to fill nan values outside extrapolation bounds. You even wrote it down here

For now, they can use this and continue with their work:

@jjlk @JouvinLea - you don't have to do this. All irfs classes duplicated the default_interp_kwargs attribute (e.g. https://github.com/gammapy/gammapy/blob/master/gammapy/irf/energy_dispersion.py#L598). So it should work like this

from gammapy.irf import EnergyDispersion2D
EnergyDispersion2D.default_interp_kwargs = ...

If you want to change the default behaviour for a given IRF class, feel free to make a PR.

@joleroi - If you agree, could you please make the change and add one test to make sure it's extrapolating?

I don't think we should add the default behaviour to extrapolate everywhere (see http://docs.gammapy.org/en/latest/development/howto.html#interpolation-and-extrapolation). Horrible things can happen. If you want to extrapolate and know what you're doing, it's not too hard to add the 1 liner I put above, no?

[cdeil]

@joleroi - I thought we discussed this with @adonath in my office a week or two ago, and decided to change to extrapolate, no?

Yes, I have changed my mind on this, I'm proposing we change our policy here.

There's just no other good way. E.g. CTA IRFs and also some HESS IRFs have first offset bin node at something > zero, but still it's common for people to query them for offset = 0, close to the first node. This works fine and is a reasonable thing to do, e.g. CTA does diffuse MC and has the first offset bin from 0 to 1 deg, and the first node at 0.5 deg.

---

Apart from many users having run into the same issue, the solution you propose doesn't really work:

from gammapy.irf import EnergyDispersion2D
EnergyDispersion2D.default_interp_kwargs = ...

The problem is that most users will access IRFs from pipelines or high-level code without even knowing about underlying IRF classes. It's bad coupling between low-level IRF interpolation and what end users should do if you ask them to configure the IRF interpolators.

So @joleroi - unless you are strongly opposed, I suggest we change the policy and default to extrapolate for all IRFs in Gammapy.
(the same problem with first FOV offset bin is there for all CTA IRFs for example)

Please!?

[joleroi]

E.g. CTA IRFs and also some HESS IRFs have first offset bin node at something > zero, but still it's common for people to query them for offset = 0, close to the first node

I see that the current scheme lead to a problem if the offset axis does not go to 0.

However, e.g. for the EnergyDispersion the current behaviour is to fill zeros outside the valid migra range.

https://github.com/gammapy/gammapy/blob/master/gammapy/irf/energy_dispersion.py#L45

This makes sense since the migration probability for extreme migra values is 0. I don't know what happens if you turn on interpolation here.

Other example, effective area: It's hard to see what happend if you extrapolate to low or high energies. For that case I've added the evaluate_fill_nan which puts 0s at low energies and a constant value at high energies.

Bottom line: Just because you run into trouble for one use case it might not be a good idea to change the behaviour for all of them. So I'm still 👎 on changing the behaviour of the NDDataArray.

What you could to is change the default interp_kwargs of the classes that include an offset axis (EnergyDisperson2D, EffectiveAreaTable2D, PSF classes don't use NDData I think) to extrapolate there. This would fix your problem, no?

[adonath]

Resolved by #932?

[cdeil]

My understanding is that everything now works for AEFF and EDISP.
I haven't checked for PSF and BACKGROUND.
Probably those classes aren't even rewritten to use NDDataArray, so that is just a separate / future issue.
Closing this one now.

@joleroi - Thanks!