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.

Sign up for GitHub

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

Jump to bottom

Change parametrisation from geom_true to energy_axis_true #2479

Merged
merged 12 commits into from Oct 24, 2019
Merged

Change parametrisation from geom_true to energy_axis_true #2479

Show file tree
Hide file tree
Changes from 4 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
2227729
Adapt MapDataset
AtreyeeS Oct 18, 2019
18447e7
Adapt MapDatasetMaker
AtreyeeS Oct 18, 2019
618e759
change BINSZ to BINSZ_DEFAULT
AtreyeeS Oct 18, 2019
acb4cc0
Adapt MapDataset
AtreyeeS Oct 18, 2019
51e7297
Adapt MapDatasetMaker
AtreyeeS Oct 18, 2019
7df9907
change BINSZ to BINSZ_DEFAULT
AtreyeeS Oct 18, 2019
f116b93
Fix high-level interface
Bultako Oct 19, 2019
d936eb9
address review comments
Bultako Oct 22, 2019
9c64a01
Merge branch 'cleanup1' of https://github.com/AtreyeeS/gammapy into c…
AtreyeeS Oct 23, 2019
4224305
Make e_true default to e_ereco
AtreyeeS Oct 23, 2019
d662561
Minor cleanup and fix of lightcurve notebook
adonath Oct 24, 2019
24c8c0d
Fix light-curve notebook
adonath Oct 24, 2019
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
83 changes: 45 additions & 38 deletions gammapy/cube/fit.py
View file
Open in desktop
Expand Up @@ -12,7 +12,7 @@
from gammapy.cube.psf_map import PSFMap
from gammapy.data import GTI
from gammapy.irf import EffectiveAreaTable, EnergyDispersion, apply_containment_fraction
from gammapy.maps import Map, MapAxis
from gammapy.maps import Map, MapAxis, WcsGeom
from gammapy.modeling import Dataset, Parameters
from gammapy.modeling.models import BackgroundModel, SkyModel, SkyModels
from gammapy.spectrum import SpectrumDataset
Expand All @@ -33,7 +33,12 @@
MIGRA_AXIS_DEFAULT = MapAxis.from_bounds(
0.2, 5, nbin=48, node_type="edges", name="migra"
)
BINSZ_IRF = 0.2
ENERGY_AXIS_DEFAULT = MapAxis.from_edges(
Copy link
Member

Choose a reason for hiding this comment

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

I think for the true energy axis we should not introduce any default, except for e_true = e_reco, as we already do in many places. The default here (looks like it was copied over from gammapy.spectrum) does not really make sense for a cube analysis. The model evaluation will be much to slow and heavily oversampled for no good reason.

Please remove the default axis here and just uses energy_axis_true = energy_axis_reco as a default, where appropriate.

Copy link
Member Author

Choose a reason for hiding this comment

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

Yes, I did indeed copy it over from gammapy.spectrum, and this is clearly a poor choice of binning (which I think we need to take from the IRFs at some point anyways). But I think it might be needed to have a default e_true binning. I remember I was getting very different indices between spectral and map based analysis once, which was solved by changing the e_true. This would be quite non-intuitive for a general user...

@registerrier @bkhelifi do you have an opinion about this?

Copy link
Member

@adonath adonath Oct 22, 2019 •
edited

Choose a reason for hiding this comment

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

Yes, the assumption of the 3D model evaluation is that we use a "fine" binning (I think we never clarified, what this actually means, probably >5 bins / decade, but we should document it somewhere). If this is not fulfilled, the user can currently get bad or unexpected results. I think in general we should change to using SpectralModel.integral() in the 3D analysis as well (just as we do for 1D), this should already improve (if not solve...) the discrepancy. Note, that for the SpectrumDatasetMaker I already removed the default for e_true. But we can certainly discuss, whether we re-introduce it in a consistent way and with a plausible choice for 1D as well as 3D.

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 that the issue is not only the fine binning but the total range.
In theory, the true energy binning is supposed to cover the 0 to + infinity range. It is of course not necessary in practice, but the range has to be larger than the e_reco one to account for spill-over.
This is far from a small approximation in many cases if users define a range close to the actual safe range.
A possibility could be to decide that by default e_true is binned starting from e_reco[0]*0.7 up to e_reco[-1]*1.3. Or even to leave as a possible parameter the leakage fraction (with default at 30% for instance).

np.logspace(-2.0, 2.5, 109), name="energy", unit=u.TeV, interp="log"
)

BINSZ_IRF_DEFAULT = 0.2
MARGIN_IRF_DEFAULT = 0.5
# TODO: Choose optimal binnings depending on IRFs


Expand Down Expand Up @@ -363,9 +368,11 @@ def from_geoms(
def create(
cls,
geom,
geom_irf=None,
energy_axis_true=None,
migra_axis=None,
rad_axis=None,
binsz_irf=None,
margin_irf=None,
reference_time="2000-01-01",
name="",
**kwargs,
Expand All @@ -376,54 +383,54 @@ def create(
----------
geom: `~gammapy.maps.WcsGeom`
Reference target geometry in reco energy, used for counts and background maps
geom_irf: `~gammapy.maps.WcsGeom`
Reference image geometry in true energy, used for IRF maps.
energy_axis_true: `~gammapy.maps.MapAxis`
True energy axis used for IRF maps
migra_axis: `~gammapy.maps.MapAxis`
Migration axis for the energy dispersion map
rad_axis: `~gammapy.maps.MapAxis`
Rad axis for the psf map
binsz_irf_maps: float
IRF Map pixel size in degrees.
margin_irf: float
IRF map margin size in degrees
reference_time: `~astropy.time.Time`
the reference time to use in GTI definition
name : str
Name of the dataset.

Returns
--------
empty_maps: `MapDataset`
A MapDataset containing zero filled maps
"""
geom_irf = geom_irf or geom.to_binsz(BINSZ_IRF)

migra_axis = migra_axis or MIGRA_AXIS_DEFAULT
rad_axis = rad_axis or RAD_AXIS_DEFAULT
energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT
binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT
margin_irf = margin_irf or MARGIN_IRF_DEFAULT
margin_irf = margin_irf * u.deg

wcs = geom.to_image()
geom_exposure = wcs.to_cube([energy_axis_true])

wcs_irf = WcsGeom.create(
binsz=binsz_irf,
width=wcs.width + margin_irf,
skydir=wcs.center_skydir,
proj=wcs.projection,
coordsys=wcs.coordsys,
)

counts = Map.from_geom(geom, unit="")

background = Map.from_geom(geom, unit="")
background_model = BackgroundModel(background)

energy_axis = geom_irf.get_axis_by_name("ENERGY")

exposure_geom = geom.to_image().to_cube([energy_axis])
exposure = Map.from_geom(exposure_geom, unit="m2 s")
exposure_irf = Map.from_geom(geom_irf, unit="m2 s")

mask_safe = np.zeros(geom.data_shape, dtype=bool)

gti = GTI.create([] * u.s, [] * u.s, reference_time=reference_time)

geom_migra = geom_irf.to_image().to_cube([migra_axis, energy_axis])
edisp_map = Map.from_geom(geom_migra, unit="")
loc = migra_axis.edges.searchsorted(1.0)
edisp_map.data[:, loc, :, :] = 1.0
edisp = EDispMap(edisp_map, exposure_irf)

geom_rad = geom_irf.to_image().to_cube([rad_axis, energy_axis])
psf_map = Map.from_geom(geom_rad, unit="sr-1")
psf = PSFMap(psf_map, exposure_irf)
geom_psf = wcs_irf.to_cube([rad_axis, energy_axis_true])
geom_edisp = wcs_irf.to_cube([migra_axis, energy_axis_true])

return cls(
counts=counts,
exposure=exposure,
psf=psf,
edisp=edisp,
background_model=background_model,
gti=gti,
mask_safe=mask_safe,
return cls.from_geoms(
geom,
geom_exposure,
geom_psf,
geom_edisp,
reference_time=reference_time,
name=name,
**kwargs,
)
Expand Down
65 changes: 49 additions & 16 deletions gammapy/cube/make.py
View file
Open in desktop
Expand Up @@ -2,16 +2,24 @@
import logging
from functools import lru_cache
import numpy as np
import astropy.units as u
from astropy.coordinates import Angle
from astropy.nddata.utils import NoOverlapError
from astropy.utils import lazyproperty
from gammapy.irf import EnergyDependentMultiGaussPSF
from gammapy.maps import Map
from gammapy.maps import Map, WcsGeom
from gammapy.modeling.models import BackgroundModel
from .background import make_map_background_irf
from .edisp_map import make_edisp_map
from .exposure import _map_spectrum_weight, make_map_exposure_true_energy
from .fit import BINSZ_IRF, MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, MapDataset
from .fit import (
BINSZ_IRF_DEFAULT,
MARGIN_IRF_DEFAULT,
MIGRA_AXIS_DEFAULT,
RAD_AXIS_DEFAULT,
ENERGY_AXIS_DEFAULT,
MapDataset,
)
from .psf_map import make_psf_map

__all__ = ["MapDatasetMaker", "MapMakerRing"]
Expand All @@ -28,14 +36,16 @@ class MapDatasetMaker:
Reference image geometry in reco energy, used for counts and background maps
offset_max : `~astropy.coordinates.Angle`
Maximum offset angle
geom_true : `~gammapy.maps.WcsGeom`
Reference image geometry in true energy, used for IRF maps. It can have a coarser
spatial bins than the counts geom.
If none, the same as geom is assumed
energy_axis_true: `~gammapy.maps.MapAxis`
True energy axis used for IRF maps
migra_axis : `~gammapy.maps.MapAxis`
Migration axis for edisp map
rad_axis : `~gammapy.maps.MapAxis`
Radial axis for psf map.
binsz_irf: float
IRF Map pixel size in degrees.
margin_irf: float
IRF map margin size in degrees
cutout : bool
Whether to cutout the observation.
cutout_mode : {'trim', 'partial', 'strict'}
Expand All @@ -47,19 +57,27 @@ def __init__(
self,
geom,
offset_max,
geom_true=None,
background_oversampling=None,
energy_axis_true=None,
migra_axis=None,
rad_axis=None,
binsz_irf=None,
margin_irf=None,
cutout_mode="trim",
cutout=True,
):

self.geom = geom
self.geom_true = geom_true if geom_true else geom.to_binsz(BINSZ_IRF)
self.offset_max = Angle(offset_max)
self.background_oversampling = background_oversampling
self.migra_axis = migra_axis if migra_axis else MIGRA_AXIS_DEFAULT
self.rad_axis = rad_axis if rad_axis else RAD_AXIS_DEFAULT
self.energy_axis_true = energy_axis_true or ENERGY_AXIS_DEFAULT
self.binsz_irf = binsz_irf or BINSZ_IRF_DEFAULT

self.margin_irf = margin_irf or MARGIN_IRF_DEFAULT
self.margin_irf = self.margin_irf * u.deg

self.cutout_mode = cutout_mode
self.cutout_width = 2 * self.offset_max
self.cutout = cutout
Expand All @@ -75,25 +93,39 @@ def _cutout_geom(self, geom, observation):
else:
return geom

@lazyproperty
def wcs_irf(self):
Copy link
Member

Choose a reason for hiding this comment

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

Please rename to geom_image_irf...

"""Spatial geometry of IRF Maps (`Geom`)"""
wcs = self.geom.to_image()
return WcsGeom.create(
binsz=self.binsz_irf,
width=wcs.width + self.margin_irf,
skydir=wcs.center_skydir,
proj=wcs.projection,
coordsys=wcs.coordsys,
)

@lazyproperty
def geom_irf(self):
Copy link
Member

Choose a reason for hiding this comment

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

Please rename to geom_exposure_irf

"""Geom of Exposure map associated with IRFs (`Geom`)"""
return self.wcs_irf.to_cube([self.energy_axis_true])

@lazyproperty
def geom_exposure(self):
"""Exposure map geom (`Geom`)"""
energy_axis = self.geom_true.get_axis_by_name("energy")
geom_exposure = self.geom.to_image().to_cube([energy_axis])
geom_exposure = self.geom.to_image().to_cube([self.energy_axis_true])
return geom_exposure

@lazyproperty
def geom_psf(self):
"""PSFMap geom (`Geom`)"""
energy_axis = self.geom_true.get_axis_by_name("ENERGY")
geom_psf = self.geom_true.to_image().to_cube([self.rad_axis, energy_axis])
geom_psf = self.wcs_irf.to_cube([self.rad_axis, self.energy_axis_true])
return geom_psf

@lazyproperty
def geom_edisp(self):
"""EdispMap geom (`Geom`)"""
energy_axis = self.geom_true.get_axis_by_name("ENERGY")
geom_edisp = self.geom_true.to_image().to_cube([self.migra_axis, energy_axis])
geom_edisp = self.wcs_irf.to_cube([self.migra_axis, self.energy_axis_true])
return geom_edisp

def make_counts(self, observation):
Expand Down Expand Up @@ -150,7 +182,8 @@ def make_exposure_irf(self, observation):
exposure : `Map`
Exposure map.
"""
geom = self._cutout_geom(self.geom_true, observation)

geom = self._cutout_geom(self.geom_irf, observation)
exposure = make_map_exposure_true_energy(
pointing=observation.pointing_radec,
livetime=observation.observation_live_time_duration,
Expand Down Expand Up @@ -281,7 +314,7 @@ def make_mask_safe_irf(self, observation):
mask : `~numpy.ndarray`
Mask
"""
geom = self._cutout_geom(self.geom_true, observation)
geom = self._cutout_geom(self.geom_irf, observation)
offset = geom.separation(observation.pointing_radec)
return offset >= self.offset_max

Expand Down
28 changes: 17 additions & 11 deletions gammapy/cube/tests/test_fit.py
View file
Open in desktop
Expand Up @@ -360,22 +360,29 @@ def test_map_fit_one_energy_bin(sky_model, geom_image):

def test_create(geom, geom_etrue):
# tests empty datasets created

migra_axis = MapAxis(nodes=np.linspace(0.0, 3.0, 51), unit="", name="migra")
rad_axis = MapAxis(nodes=np.linspace(0.0, 1.0, 51), unit="deg", name="theta")
e_reco = MapAxis.from_edges(
np.logspace(-1.0, 1.0, 3), name="energy", unit=u.TeV, interp="log"
)
e_true = MapAxis.from_edges(
np.logspace(-1.0, 1.0, 4), name="energy", unit=u.TeV, interp="log"
)
geom = WcsGeom.create(binsz=0.02, width=(2, 2), axes=[e_reco])
empty_dataset = MapDataset.create(
geom=geom, energy_axis_true=e_true, migra_axis=migra_axis, rad_axis=rad_axis
)

empty_dataset = MapDataset.create(geom, geom_etrue, migra_axis, rad_axis)

assert_allclose(empty_dataset.counts.data.sum(), 0.0)
assert_allclose(empty_dataset.background_model.map.data.sum(), 0.0)
assert empty_dataset.counts.data.shape == (2, 100, 100)

assert empty_dataset.psf.psf_map.data.shape == (3, 50, 100, 100)
assert empty_dataset.psf.exposure_map.data.shape == (3, 1, 100, 100)
assert empty_dataset.exposure.data.shape == (3, 100, 100)

assert empty_dataset.edisp.edisp_map.data.shape == (3, 50, 100, 100)
assert empty_dataset.edisp.exposure_map.data.shape == (3, 1, 100, 100)
assert empty_dataset.psf.psf_map.data.shape == (3, 50, 12, 12)
assert empty_dataset.psf.exposure_map.data.shape == (3, 1, 12, 12)

assert_allclose(empty_dataset.edisp.edisp_map.data.sum(), 30000)
assert empty_dataset.edisp.edisp_map.data.shape == (3, 50, 12, 12)
assert empty_dataset.edisp.exposure_map.data.shape == (3, 1, 12, 12)
assert_allclose(empty_dataset.edisp.edisp_map.data.sum(), 432)

assert_allclose(empty_dataset.gti.time_delta, 0.0 * u.s)

Expand All @@ -394,7 +401,6 @@ def test_from_geoms():
wcs_irf = WcsGeom.create(binsz=0.1, width=(2.5, 2.5))
geom = wcs.to_cube([e_reco])
geom_exposure = wcs.to_cube([e_true])
geom_irf = wcs_irf.to_cube([e_true])
geom_psf = wcs_irf.to_cube([rad_axis, e_true])
geom_edisp = wcs_irf.to_cube([migra_axis, e_true])

Expand Down