Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add Background 2D class #1312

Merged
merged 5 commits into from Feb 16, 2018
Merged
Show file tree
Hide file tree
Changes from 3 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
129 changes: 127 additions & 2 deletions gammapy/irf/background.py
Expand Up @@ -7,10 +7,11 @@
from ..utils.nddata import NDDataArray, BinnedDataAxis
from ..utils.scripts import make_path
from ..utils.fits import fits_table_to_table, table_to_fits_table

import numpy as np

__all__ = [
'Background3D',
'Background2D',
]


Expand Down Expand Up @@ -132,4 +133,128 @@ def to_table(self):

def to_fits(self, name='BACKGROUND'):
"""Convert to `~astropy.io.fits.BinTable`."""
return table_to_fits_table(self.to_table(), name)
return table_to_fits_table(self.to_table(), name)


class Background2D(object):
"""Background 2D.

Data format specification: :ref:`gadf:bkg_3d`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

3d -> 2d


Parameters
-----------
energy_lo, energy_hi : `~astropy.units.Quantity`
Energy binning
offset_lo, offset_hi : `~astropy.units.Quantity`
FOV coordinate offset-axis binning
data : `~astropy.units.Quantity`
Background rate (usually: ``s^-1 MeV^-1 sr^-1``)
"""
default_interp_kwargs = dict(bounds_error=False, fill_value=None)
"""Default Interpolation kwargs for `~NDDataArray`. Extrapolate."""

def __init__(self, energy_lo, energy_hi,
offset_lo, offset_hi,
data, meta=None, interp_kwargs=None):

if interp_kwargs is None:
interp_kwargs = self.default_interp_kwargs
axes = [
BinnedDataAxis(
energy_lo, energy_hi,
interpolation_mode='log', name='energy'),
BinnedDataAxis(
offset_lo, offset_hi,
interpolation_mode='linear', name='offset'),
]
self.data = NDDataArray(axes=axes, data=data,
interp_kwargs=interp_kwargs)
self.meta = OrderedDict(meta) if meta else OrderedDict()

def __str__(self):
ss = self.__class__.__name__
ss += '\n{}'.format(self.data)
return ss

@classmethod
def from_table(cls, table):
"""Read from `~astropy.table.Table`."""
# Spec says key should be "BKG", but there are files around
# (e.g. CTA 1DC) that use "BGD". For now we support both
if 'BKG' in table.colnames:
bkg_name = 'BKG'
elif 'BGD' in table.colnames:
bkg_name = 'BGD'
else:
raise ValueError('Invalid column names. Need "BKG" or "BGD".')

# Currently some files (e.g. CTA 1DC) contain unit in the FITS file
# '1/s/MeV/sr', which is invalid ( try: astropy.unit.Unit('1/s/MeV/sr')
# This should be corrected.
# For now, we hard-code the unit here:
data_unit = u.Unit('s-1 MeV-1 sr-1')
return cls(
energy_lo=table['ENERG_LO'].quantity[0],
energy_hi=table['ENERG_HI'].quantity[0],
offset_lo=table['THETA_LO'].quantity[0],
offset_hi=table['THETA_HI'].quantity[0],
data=table[bkg_name].data[0] * data_unit,
meta=table.meta,
)

@classmethod
def from_hdulist(cls, hdulist, hdu='BACKGROUND'):
"""Create from `~astropy.io.fits.HDUList`."""
fits_table = hdulist[hdu]
table = fits_table_to_table(fits_table)
return cls.from_table(table)

@classmethod
def read(cls, filename, hdu='BACKGROUND'):
"""Read from file."""
filename = make_path(filename)
hdulist = fits.open(str(filename))
return cls.from_hdulist(hdulist, hdu=hdu)

def to_table(self):
"""Convert to `~astropy.table.Table`."""
meta = self.meta.copy()
table = Table(meta=meta)

table['THETA_LO'] = self.data.axis('offset').lo[np.newaxis]
table['THETA_HI'] = self.data.axis('offset').hi[np.newaxis]
table['ENERG_LO'] = self.data.axis('energy').lo[np.newaxis]
table['ENERG_HI'] = self.data.axis('energy').hi[np.newaxis]
table['BKG'] = self.data.data[np.newaxis]
return table

def to_fits(self, name='BACKGROUND'):
"""Convert to `~astropy.io.fits.BinTable`."""
return table_to_fits_table(self.to_table(), name)

def evaluate(self, fov_offset, fov_phi=None, energy_reco=None, **kwargs):
"""
Evaluate the `Background2D` at a given offset and energy.

Parameters
----------
fov_offset : `~astropy.coordinates.Angle`
offset in the FOV
fov_phi: `~astropy.coordinates.Angle`
azimuth angle in the FOV. Not used for this class since the background model is radially symmetric
energy_reco : `~astropy.units.Quantity`
Vector of energy (1D) on which the model is evaluated
kwargs : dict
option for interpolation for `~scipy.interpolate.RegularGridInterpolator`

Returns
-------
array : `~astropy.units.Quantity`
Interpolated values, axis order is the same as for the NDData array

"""
if energy_reco is None:
energy_reco = self.data.axis('energy').nodes

array = self.data.evaluate(offset=fov_offset, energy=energy_reco, **kwargs)
return array
62 changes: 61 additions & 1 deletion gammapy/irf/tests/test_background.py
Expand Up @@ -2,10 +2,15 @@
from __future__ import absolute_import, division, print_function, unicode_literals
import pytest
import astropy.units as u
import numpy as np
from numpy.testing import assert_allclose, assert_equal
from astropy.tests.helper import assert_quantity_allclose
from astropy.coordinates import Angle
from astropy.io import fits
from ...utils.testing import requires_dependency, requires_data
from ..background import Background3D
from ..background import Background3D, Background2D
from ...utils.fits import table_to_fits_table
from ...utils.energy import EnergyBounds


@pytest.fixture(scope='session')
Expand Down Expand Up @@ -49,3 +54,58 @@ def test_background_3d_write(bkg_3d):
hdu = table_to_fits_table(bkg_3d.to_table())
assert_equal(hdu.data['DETX_LO'][0], bkg_3d.data.axis('detx').lo.value)
assert hdu.header['TUNIT1'] == bkg_3d.data.axis('detx').lo.unit


def make_test_array():
# Create a dummy `Background2D`
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it would be better to have an even simpler test object. (even less code, easier to think about in the test methods below)

How about this?

energy = [1, 3, 5, 10] * u.TeV
offset = [0, 1, 5] * u.deg
data = energy.value * offset.value * u.Unit('s-1 MeV-1 sr-1')
return Background2D(
    energy_lo=energy[:-1], energy_hi=energy[1:],
    offset_lo=offset[:-1], offset_hi=offset[1:],
    data=data,
)

This is small, i.e. fast. And it's less code. And it's simpler to see what the expected value is in evaluate calls, if you choose a few points carefully (e.g. at a node and in the middle between two nodes, and throw in a third case that leads to extrapolation (e.g. offset=0) to establish behaviour).

I think we should do all our IRF tests like this: have super-simple test cases and fill them with data values that is easy and generic to write when making the test case, and also to know what to expect in evaluate.

energy = [1, 10, 100] * u.TeV
offset = [0, 1, 2, 3] * u.deg
data = np.zeros((len(energy)-1,len(offset)-1)) * u.Unit('s-1 MeV-1 sr-1')
data.value[1,0]=2
data.value[1,1]=4
return Background2D(
energy_lo=energy[:-1], energy_hi=energy[1:],
offset_lo=offset[:-1], offset_hi=offset[1:],
data=data,
)


def test_background2d_read_write(tmpdir):
bkg_2d_1 = make_test_array()
filename = str(tmpdir / "bkg2d.fits")
prim_hdu = fits.PrimaryHDU()
hdu_bkg = bkg_2d_1.to_fits()
hdulist = fits.HDUList([prim_hdu, hdu_bkg])
hdulist.writeto(filename)
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cdeil
Still get the same issue that if I am not writing the file, with only using from_hdulist(), I lost the dimensions quantity... I checked for example for the effective area, they are writting the file and using the function read then.

bkg_2d_2 = Background2D.read(filename)

axis = bkg_2d_2.data.axis('energy')
assert axis.nbins == 2
assert axis.unit == 'TeV'

axis = bkg_2d_2.data.axis('offset')
assert axis.nbins == 3
assert axis.unit == 'deg'

data = bkg_2d_2.data.data
assert data.shape == (2, 3)
assert data.unit == u.Unit('s-1 MeV-1 sr-1')


@requires_dependency('scipy')
def test_background2d_evaluate():
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This needs a requires_dependency('scipy').
Generally I would suggest you follow the coding pattern from the 3D case, i.e. return the test case as a pytest fixture and inject it into tests.

bkg_2d = make_test_array()
data_unit = u.Unit('s-1 MeV-1 sr-1')

# Interoplate at the energy and offset bin of the bgk_2d axes
off = bkg_2d.data.axis('offset').nodes
e_reco = bkg_2d.data.axis('energy').nodes
res = bkg_2d.evaluate(fov_offset=off, energy_reco=e_reco)
assert_quantity_allclose(res, bkg_2d.data.data)

#Test linear interpolation for offset
Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@cdeil
This is easy to to a easy test example regarding the offset. But since the nodes will be the log center for the energy then I can not use just integer anymore

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I would manually compute sqrt(e_1 * e_2), i.e. log center, and then use that in the evaluate call, see if it's giving the mean as value between the IRF value at the two nodes in one case.

off_tab=Angle(np.array([1]),"deg")
e_reco_tab=bkg_2d.data.axis('energy').nodes[1]
res = bkg_2d.evaluate(fov_offset=off_tab, energy_reco=e_reco_tab)
assert_quantity_allclose(res, 3 * data_unit)