From b4855c1893141aaccb3cb8e8eb46cb70d9cb4e03 Mon Sep 17 00:00:00 2001 From: Fabio Pintore Date: Wed, 27 Nov 2019 17:08:15 +0100 Subject: [PATCH 1/6] Add sample_background in MapDatasetEventSampler --- gammapy/cube/simulate.py | 63 ++++++++++++++++++++++++++++- gammapy/cube/tests/test_simulate.py | 60 ++++++++++++++++++++++++++- 2 files changed, 120 insertions(+), 3 deletions(-) diff --git a/gammapy/cube/simulate.py b/gammapy/cube/simulate.py index bc16b069dc..030f199f43 100644 --- a/gammapy/cube/simulate.py +++ b/gammapy/cube/simulate.py @@ -1,5 +1,7 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst """Simulate observations""" +import numpy as np +from astropy.table import Table import astropy.units as u from gammapy.cube import ( MapDataset, @@ -7,11 +9,18 @@ make_map_background_irf, make_map_exposure_true_energy, ) +from gammapy.data import EventList from gammapy.maps import WcsNDMap from gammapy.modeling.models import BackgroundModel +from gammapy.modeling.models import ( + ConstantTemporalModel, + LightCurveTemplateTemporalModel, + PhaseCurveTemplateTemporalModel, +) + from gammapy.utils.random import get_random_state -__all__ = ["simulate_dataset"] +__all__ = ["simulate_dataset", "MapDatasetEventSampler"] def simulate_dataset( @@ -89,3 +98,55 @@ def simulate_dataset( dataset.counts = WcsNDMap(geom, counts) return dataset + + +class MapDatasetEventSampler: + """Sample events from a map dataset + + Parameters + ---------- + random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`} + Defines random number generator initialisation. + Passed to `~gammapy.utils.random.get_random_state`. + """ + + def __init__(self, random_state="random-seed"): + self.random_state = get_random_state(random_state) + + def sample_background(self, dataset): + """Sample background + Parameters + ---------- + dataset : `MapDataset` + Map dataset. + + Returns + ------- + events : `EventList` + Background events + """ + table = Table() + + background = dataset.background_model.evaluate() + n_events = self.random_state.poisson(np.sum(background.data)) + + # sample position + coords = background.sample_coord(n_events, self.random_state) + table["ENERGY"] = coords["energy"] + table["RA"] = coords["lon"] + table["DEC"] = coords["lat"] + table["MC_ID"] = 0 + + # sample time + t_start, t_stop, t_ref = ( + dataset.gti.time_start, + dataset.gti.time_stop, + dataset.gti.time_ref, + ) + model = ConstantTemporalModel() + time = model.sample_time(n_events, t_start, t_stop, self.random_state) + table["TIME"] = u.Quantity( + ((time.mjd - dataset.gti.time_ref.mjd) * u.day).to(u.s) + ).value + + return EventList(table) diff --git a/gammapy/cube/tests/test_simulate.py b/gammapy/cube/tests/test_simulate.py index 4342497bc6..8e6cb36563 100644 --- a/gammapy/cube/tests/test_simulate.py +++ b/gammapy/cube/tests/test_simulate.py @@ -1,11 +1,14 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import numpy as np from numpy.testing import assert_allclose +from astropy.time import Time import astropy.units as u from astropy.coordinates import SkyCoord -from gammapy.cube import MapDataset, simulate_dataset +from gammapy.cube import MapDataset, simulate_dataset, MapDatasetEventSampler +from gammapy.cube.tests.test_fit import get_map_dataset +from gammapy.data import GTI from gammapy.irf import load_cta_irfs -from gammapy.maps import MapAxis, WcsGeom +from gammapy.maps import Map, MapAxis, WcsGeom from gammapy.modeling.models import ( GaussianSpatialModel, PowerLawSpectralModel, @@ -58,3 +61,56 @@ def test_simulate(): ) assert_allclose(dataset.psf.data[5, 32, 32], 0.04203219) assert_allclose(dataset.edisp.data.data[10, 10], 0.85944298, rtol=1e-5) + + +def dataset_maker(): + position = SkyCoord(0.0, 0.0, frame="galactic", unit="deg") + energy_axis = MapAxis.from_bounds( + 1, 100, nbin=30, unit="TeV", name="energy", interp="log" + ) + + exposure = Map.create( + binsz=0.02, + map_type="wcs", + skydir=position, + width="2 deg", + axes=[energy_axis], + coordsys="GAL", + unit="cm2 s", + ) + + spatial_model = GaussianSpatialModel( + lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg", frame="galactic" + ) + + spectral_model = PowerLawSpectralModel(amplitude="1e-11 cm-2 s-1 TeV-1") + skymodel = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model) + + geom = WcsGeom.create( + skydir=position, binsz=0.02, width="5 deg", coordsys="GAL", axes=[energy_axis] + ) + + t_ref = Time("2010-01-01T00:00:00") + t_min = 0 * u.s + t_max = 30000 * u.s + + gti = GTI.create(start=t_min, stop=t_max) + + dataset = get_map_dataset( + sky_model=skymodel, geom=geom, geom_etrue=geom, edisp=True + ) + dataset.gti = gti + + return dataset + + +def test_MDE_sample_background(): + dataset = dataset_maker() + sampler = MapDatasetEventSampler(random_state=0) + bkg_evt = sampler.sample_background(dataset=dataset) + + assert len(bkg_evt.table["ENERGY"]) == 375084 + assert_allclose(bkg_evt.table["ENERGY"][0], 2.1613281656472028, rtol=1e-5) + assert_allclose(bkg_evt.table["RA"][0], 0.7167667097717688, rtol=1e-5) + assert_allclose(bkg_evt.table["DEC"][0], 1.1422945173936792, rtol=1e-5) + assert_allclose(bkg_evt.table["MC_ID"][0], 0, rtol=1e-5) From 35a175a187e528cf265424ece1951e4c14b7239d Mon Sep 17 00:00:00 2001 From: Fabio Pintore Date: Wed, 27 Nov 2019 17:12:25 +0100 Subject: [PATCH 2/6] Add sample_background in MapDatasetEventSampler --- gammapy/cube/tests/test_simulate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/gammapy/cube/tests/test_simulate.py b/gammapy/cube/tests/test_simulate.py index 8e6cb36563..699d6b9e92 100644 --- a/gammapy/cube/tests/test_simulate.py +++ b/gammapy/cube/tests/test_simulate.py @@ -111,6 +111,7 @@ def test_MDE_sample_background(): assert len(bkg_evt.table["ENERGY"]) == 375084 assert_allclose(bkg_evt.table["ENERGY"][0], 2.1613281656472028, rtol=1e-5) + assert_allclose(bkg_evt.table["ENERGY"][1], 1.7671637269378804, rtol=1e-5) assert_allclose(bkg_evt.table["RA"][0], 0.7167667097717688, rtol=1e-5) assert_allclose(bkg_evt.table["DEC"][0], 1.1422945173936792, rtol=1e-5) assert_allclose(bkg_evt.table["MC_ID"][0], 0, rtol=1e-5) From 71d32e21f5aac37ed18f48f73692ade3035b4111 Mon Sep 17 00:00:00 2001 From: Fabio Pintore Date: Wed, 27 Nov 2019 17:17:50 +0100 Subject: [PATCH 3/6] Add sample_background in MapDatasetEventSampler --- gammapy/cube/tests/test_simulate.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gammapy/cube/tests/test_simulate.py b/gammapy/cube/tests/test_simulate.py index 699d6b9e92..8e6cb36563 100644 --- a/gammapy/cube/tests/test_simulate.py +++ b/gammapy/cube/tests/test_simulate.py @@ -111,7 +111,6 @@ def test_MDE_sample_background(): assert len(bkg_evt.table["ENERGY"]) == 375084 assert_allclose(bkg_evt.table["ENERGY"][0], 2.1613281656472028, rtol=1e-5) - assert_allclose(bkg_evt.table["ENERGY"][1], 1.7671637269378804, rtol=1e-5) assert_allclose(bkg_evt.table["RA"][0], 0.7167667097717688, rtol=1e-5) assert_allclose(bkg_evt.table["DEC"][0], 1.1422945173936792, rtol=1e-5) assert_allclose(bkg_evt.table["MC_ID"][0], 0, rtol=1e-5) From 5a2a5f981cf19edd329afe1cc289bfdf049cbef8 Mon Sep 17 00:00:00 2001 From: Fabio Pintore Date: Thu, 28 Nov 2019 11:12:22 +0100 Subject: [PATCH 4/6] Implementing Axel comments --- gammapy/cube/simulate.py | 21 ++++++++------------- gammapy/cube/tests/test_simulate.py | 16 +++------------- 2 files changed, 11 insertions(+), 26 deletions(-) diff --git a/gammapy/cube/simulate.py b/gammapy/cube/simulate.py index 030f199f43..e96161d7bd 100644 --- a/gammapy/cube/simulate.py +++ b/gammapy/cube/simulate.py @@ -12,11 +12,7 @@ from gammapy.data import EventList from gammapy.maps import WcsNDMap from gammapy.modeling.models import BackgroundModel -from gammapy.modeling.models import ( - ConstantTemporalModel, - LightCurveTemplateTemporalModel, - PhaseCurveTemplateTemporalModel, -) +from gammapy.modeling.models import ConstantTemporalModel from gammapy.utils.random import get_random_state @@ -115,10 +111,11 @@ def __init__(self, random_state="random-seed"): def sample_background(self, dataset): """Sample background + Parameters ---------- dataset : `MapDataset` - Map dataset. + Map dataset. Returns ------- @@ -133,20 +130,18 @@ def sample_background(self, dataset): # sample position coords = background.sample_coord(n_events, self.random_state) table["ENERGY"] = coords["energy"] - table["RA"] = coords["lon"] - table["DEC"] = coords["lat"] + table["RA"] = coords.skycoord.icrs.ra.deg + table["DEC"] = coords.skycoord.icrs.dec.deg table["MC_ID"] = 0 # sample time - t_start, t_stop, t_ref = ( + time_start, time_stop, time_ref = ( dataset.gti.time_start, dataset.gti.time_stop, dataset.gti.time_ref, ) model = ConstantTemporalModel() - time = model.sample_time(n_events, t_start, t_stop, self.random_state) - table["TIME"] = u.Quantity( - ((time.mjd - dataset.gti.time_ref.mjd) * u.day).to(u.s) - ).value + time = model.sample_time(n_events, time_start, time_stop, self.random_state) + table["TIME"] = u.Quantity(((time.mjd - time_ref.mjd) * u.day).to(u.s)).value return EventList(table) diff --git a/gammapy/cube/tests/test_simulate.py b/gammapy/cube/tests/test_simulate.py index 8e6cb36563..6092ad014e 100644 --- a/gammapy/cube/tests/test_simulate.py +++ b/gammapy/cube/tests/test_simulate.py @@ -63,22 +63,13 @@ def test_simulate(): assert_allclose(dataset.edisp.data.data[10, 10], 0.85944298, rtol=1e-5) +@requires_data() def dataset_maker(): position = SkyCoord(0.0, 0.0, frame="galactic", unit="deg") energy_axis = MapAxis.from_bounds( 1, 100, nbin=30, unit="TeV", name="energy", interp="log" ) - exposure = Map.create( - binsz=0.02, - map_type="wcs", - skydir=position, - width="2 deg", - axes=[energy_axis], - coordsys="GAL", - unit="cm2 s", - ) - spatial_model = GaussianSpatialModel( lon_0="0 deg", lat_0="0 deg", sigma="0.2 deg", frame="galactic" ) @@ -90,7 +81,6 @@ def dataset_maker(): skydir=position, binsz=0.02, width="5 deg", coordsys="GAL", axes=[energy_axis] ) - t_ref = Time("2010-01-01T00:00:00") t_min = 0 * u.s t_max = 30000 * u.s @@ -111,6 +101,6 @@ def test_MDE_sample_background(): assert len(bkg_evt.table["ENERGY"]) == 375084 assert_allclose(bkg_evt.table["ENERGY"][0], 2.1613281656472028, rtol=1e-5) - assert_allclose(bkg_evt.table["RA"][0], 0.7167667097717688, rtol=1e-5) - assert_allclose(bkg_evt.table["DEC"][0], 1.1422945173936792, rtol=1e-5) + assert_allclose(bkg_evt.table["RA"][0], 265.7253792887848, rtol=1e-5) + assert_allclose(bkg_evt.table["DEC"][0], -27.727581635186304, rtol=1e-5) assert_allclose(bkg_evt.table["MC_ID"][0], 0, rtol=1e-5) From 0e416adf74ed1fe87655a17204d6c061f7c065b4 Mon Sep 17 00:00:00 2001 From: Fabio Pintore Date: Thu, 28 Nov 2019 11:17:02 +0100 Subject: [PATCH 5/6] Implementing Axel comments --- gammapy/cube/tests/test_simulate.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gammapy/cube/tests/test_simulate.py b/gammapy/cube/tests/test_simulate.py index 6092ad014e..9483b1f4b2 100644 --- a/gammapy/cube/tests/test_simulate.py +++ b/gammapy/cube/tests/test_simulate.py @@ -1,14 +1,13 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import numpy as np from numpy.testing import assert_allclose -from astropy.time import Time import astropy.units as u from astropy.coordinates import SkyCoord from gammapy.cube import MapDataset, simulate_dataset, MapDatasetEventSampler from gammapy.cube.tests.test_fit import get_map_dataset from gammapy.data import GTI from gammapy.irf import load_cta_irfs -from gammapy.maps import Map, MapAxis, WcsGeom +from gammapy.maps import MapAxis, WcsGeom from gammapy.modeling.models import ( GaussianSpatialModel, PowerLawSpectralModel, From 04fecb27273570bef97f7e49fe3afd92b0abcf11 Mon Sep 17 00:00:00 2001 From: Axel Donath Date: Thu, 28 Nov 2019 14:13:55 +0100 Subject: [PATCH 6/6] Move requires_data decorator --- gammapy/cube/tests/test_simulate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gammapy/cube/tests/test_simulate.py b/gammapy/cube/tests/test_simulate.py index 9483b1f4b2..061e579d5a 100644 --- a/gammapy/cube/tests/test_simulate.py +++ b/gammapy/cube/tests/test_simulate.py @@ -62,7 +62,6 @@ def test_simulate(): assert_allclose(dataset.edisp.data.data[10, 10], 0.85944298, rtol=1e-5) -@requires_data() def dataset_maker(): position = SkyCoord(0.0, 0.0, frame="galactic", unit="deg") energy_axis = MapAxis.from_bounds( @@ -93,6 +92,7 @@ def dataset_maker(): return dataset +@requires_data() def test_MDE_sample_background(): dataset = dataset_maker() sampler = MapDatasetEventSampler(random_state=0)