Add n-dim data base class for gammapy.irf

[joleroi]

This PR introduces gammapy.utils.nddata
This module contains a base class for an NDData array. It is supposed to be subclassed by all classes in `gammapy.irf`` and maybe in other places.

An scratch EffectiveArea2D class is also part of this PR for the sake of illustration.

[cdeil]

Use this:

from astropy.tests.helper import pytest

[cdeil]

I've done a superficial review and left some comments inline.

Do you want to apply it to AEFF or some other IRF in this PR, or add some simple dummy IRF in the test or as an example? (I'm not sure what is best ... either way this will be a large PR, bit it'll still be OK to review.)

[joleroi]

Thanks for you comments. I added an minimal aeff example file. Next thing is to test if the NDDataArray FITS I/O works for the test datastore

cc @JouvinLea Since you're also working on IRF stuff maybe you want to have a look, too? If it comes to the worst you will have to use this class at some point 😉

[JouvinLea]

Sorry I really didn't follow this PR... What is the purpose exactly?

[joleroi]

I'll answer via email

[joleroi]

@cdeil I added an example EffectiveArea2D class, that uses the NDData base class. Please review. I will continue with moving functionality from the old EffectiveAreaTable2D class in the next days.

[joleroi]

The axis names are currently just the column names in the Fits file, we might want to map this to energy and offset?

[cdeil]

Yes, probably a mechanism to map FITS axis names to something else would be nice, for cases where the FITS names are cryptic.

[cdeil]

There's a problem here and a solution here. 😄

[cdeil]

I looked at this for 10 minutes, but I would need at least an hour to read through your implementation and make good suggestions. I can do this later if you want, but not this week.

Naively I would have expected some scheme where a given IRF is a subclass and has class-level attributes to declare which axes are present and maybe some other things.
This is how e.g. Django models work or Astropy models or Astropy coordinate frames.
See https://docs.djangoproject.com/en/1.9/topics/db/models/

The axes are class-specific and not instance-specific in your implementation, right?
Then probably the axis configuration declaration for a given IRF should probably be class-level attributes like in the Django model example I linked to, no?

(The alternative would be that there is a factory function to create specific IRFs or a single class that is configurable on construction.)

Does this comment make any sense? Or should I give a code example how I would have expected the AEFF class implementation to look like?

[cdeil]

This is an example of a single class that is configurable on construction:
http://xarray.pydata.org/en/stable/data-structures.html#creating-a-dataarray

For our IRFs I don't know if one configurable class or a sub-class zoo would work better.

[joleroi]

@cdeil thanks. you already helped a lot. Did I understand you correctly that you would have expected something like

class EffArea
    __init__(self, energy, offset, ...):
        self.energy = DataAxis(energy)
        self.offset = DataAxis(offset)

This is how I had it in mind originally. I changed to the 'weird' scheme because I have a from_table constructor in the base class (NDDataArray) that interferes with that (I already adds axes, so the subclasses cannot call this constructor). But as you say, the from_table methods should be more specific I guess (read a fixes set of columns, instead of trying to guess). I put it to the base class in order to be able to read/write any NDDataArray and not only the subclasses.

I will change this, and then add the plotting methods etc. to make the review easier, ok?

UPDATE: The NDDataArray has an axes attribute, that is use to set up the interpolator etc.This would still be present in the EffectiveArea2d. So the constructor would be more like

class EffArea
    __init__(self, energy, offset, ...):
         super.__init__()
        self.energy = self.axes[0]
        self.offset = self.axes[1]

[cdeil]

@joleroi - To be honest I don't know what I want / expect here.

The goal should be that the class is nice to use that that each individual IRF class is small, basically declarative which axes and binning are present, the boilerplate code should be in the base class.

Something like this:

In [3]: Disk2D.__init__??
Signature: Disk2D.__init__(self, amplitude=1, x_0=0, y_0=0, R_0=1, **kwargs)
Source:
            def __init__(self, *params, **kwargs):
                return super(cls, self).__init__(*params, **kwargs)

File:      ~/Library/Python/3.5/lib/python/site-packages/astropy-1.2.dev15393-py3.5-macosx-10.11-x86_64.egg/astropy/modeling/core.py
Type:      function

In [4]: Disk2D
Out[4]: 
<class 'astropy.modeling.functional_models.Disk2D'>
Name: Disk2D
Inputs: ('x', 'y')
Outputs: ('z',)
Fittable parameters: ('amplitude', 'x_0', 'y_0', 'R_0')

In [5]: from astropy.modeling.models import Disk2D

There the parameters (axes in your case) are class-level attributed (see here) and then some dark magic is generating __init__.

I now see that in Sherpa the parameters are constructed in __init__ (like you're constructing the axes in __init__ now:
https://github.com/sherpa/sherpa/blob/3dc0c33146ea5c93b0dcff07be804e25769d7642/sherpa/astro/models/__init__.py#L582

I'll try to think about this and read up on class-level attributes in Python over the weekend.

But yes ... adding tests or examples that show how your current class is used and that it works would be very helpful for review.

[joleroi]

@cdeil
I start to understand the 2 different schemes you propose. Right now I don't know how I should make the base class work with class level attributes. All the interpolation etc. work with calls to self.axes and loops. How could this be changed to use class level attributes like self.offset? I will take 10 min. to try to understand this 😄

Otherwise I will continue with the 'construction on initialization' scheme. If we encounter severe problems later, we have to find a way to switch to the other scheme.

[joleroi]

I worked in this with @adonath and we think the subclassing API is sufficiently elegant for now (not as fancy as astropy.models but it does the job, if someone has time it can always be improved). Axes are class level attributes defined in each subclass now.
https://github.com/gammapy/gammapy/pull/527/files#diff-cfb4e8b67c434385b17b6cd254d4c5b0R42
The order has to be defined in a separate attribute axis_names (astropy magically gets the order from the order the class level attributes are written in)

The drawback is that NDData cannot be instanciated anymore (since it does not have any axes defined).

[joleroi]

@cdeil I will merge this now. The example class EffectiveArea2d covers all the functionality from the old class EffectiveAreaTable2d but is not used anywhere. Next I will add a 1D effective area class based on NDDataArray. You can have a look next week and then I will remove the old irf.effective_area_table module and replace it with the new classes, if nobody raises any concerns.