Merge pull request #2994 from AtreyeeS/exposure-reco

Clean up estimators to use exposure in reco energy

gammapy/estimators/asmooth_map.py

@@ -7,9 +7,9 @@
 from gammapy.maps import WcsNDMap
 from gammapy.stats import CashCountsStatistic
 from gammapy.utils.array import scale_cube
-from gammapy.makers.utils import _map_spectrum_weight
 from gammapy.modeling.models import PowerLawSpectralModel
 from .core import Estimator
+from .utils import estimate_exposure_reco_energy
 
 __all__ = ["ASmoothMapEstimator"]
 
@@ -45,7 +45,14 @@ class ASmoothMapEstimator(Estimator):
 
     tag = "ASmoothMapEstimator"
 
-    def __init__(self, scales=None, kernel=Gaussian2DKernel, spectrum=None, method="lima", threshold=5):
+    def __init__(
+        self,
+        scales=None,
+        kernel=Gaussian2DKernel,
+        spectrum=None,
+        method="lima",
+        threshold=5,
+    ):
         if spectrum is None:
             spectrum = PowerLawSpectralModel()
 
@@ -136,7 +143,7 @@ def run(self, dataset):
         background = background.sum_over_axes(keepdims=False)
 
         if dataset.exposure is not None:
-            exposure = _map_spectrum_weight(dataset.exposure, self.spectrum)
+            exposure = estimate_exposure_reco_energy(dataset, self.spectrum)
             exposure = exposure.sum_over_axes(keepdims=False)
         else:
             exposure = None

gammapy/estimators/excess_map.py

@@ -8,6 +8,7 @@
 from gammapy.maps import Map
 from gammapy.stats import CashCountsStatistic, WStatCountsStatistic
 from .core import Estimator
+from .utils import estimate_exposure_reco_energy
 
 __all__ = [
     "ExcessMapEstimator",
@@ -78,6 +79,7 @@ class ExcessMapEstimator(Estimator):
         Which additional maps to estimate besides delta TS, significance and symmetric error.
         Available options are:
 
+            * "flux": estimate flux map
             * "errn-errp": estimate asymmetric errors.
             * "ul": estimate upper limits.
 
@@ -136,6 +138,7 @@ def run(self, dataset):
                 * ts : delta TS map
                 * significance : sqrt(delta TS), or Li-Ma significance map
                 * err : symmetric error map (from covariance)
+                * flux : flux map. An exposure map must be present in the dataset to compute flux map
                 * errn : negative error map
                 * errp : positive error map
                 * ul : upper limit map
@@ -170,6 +173,15 @@ def run(self, dataset):
         err = Map.from_geom(geom, data=counts_stat.error * self.n_sigma)
         result.update({"err": err})
 
+        if dataset.exposure:
+            flux = excess / estimate_exposure_reco_energy(dataset)
+            flux.quantity = flux.quantity.to("1 / (cm2 s)")
+        else:
+            flux = Map.from_geom(
+                geom=dataset.counts.geom, data=np.nan * np.ones(dataset.data_shape)
+            )
+        result.update({"flux": flux})
+
         if "errn-errp" in self.selection_optional:
             errn = Map.from_geom(geom, data=counts_stat.compute_errn(self.n_sigma))
             errp = Map.from_geom(geom, data=counts_stat.compute_errp(self.n_sigma))
@@ -180,4 +192,5 @@ def run(self, dataset):
                 geom, data=counts_stat.compute_upper_limit(self.n_sigma_ul)
             )
             result.update({"ul": ul})
+
         return result

gammapy/estimators/tests/test_excess_map.py

@@ -147,10 +147,16 @@ def test_significance_map_estimator_map_dataset_on_off(simple_dataset_on_off):
     estimator_image = ExcessMapEstimator(
         0.11 * u.deg, apply_mask_fit=True, return_image=True
     )
+    simple_dataset_on_off.exposure.data = (
+        np.ones(simple_dataset_on_off.exposure.data.shape) * 1e6
+    )
     result_image = estimator_image.run(simple_dataset_on_off)
     assert result_image["counts"].data.shape == (1, 20, 20)
     assert_allclose(result_image["significance"].data[0, 10, 10], 5.08179, atol=1e-3)
 
+    assert result_image["flux"].unit == u.Unit("cm-2s-1")
+    assert_allclose(result_image["flux"].data[0, 10, 10], 7.6e-9, rtol=1e-3)
+
 
 def test_incorrect_selection():
     with pytest.raises(ValueError):

gammapy/estimators/tests/test_ts_map.py

@@ -83,7 +83,7 @@ def fermi_dataset():
         mask_safe=mask_safe,
         psf=psfmap,
         name="fermi-3fhl-gc",
-        edisp=edisp
+        edisp=edisp,
     )
 
     return dataset.to_image()
@@ -119,11 +119,11 @@ def test_compute_ts_map_psf(fermi_dataset):
 
     assert_allclose(result["ts"].data[29, 29], 835.140605, rtol=1e-2)
     assert_allclose(result["niter"].data[29, 29], 7)
-    assert_allclose(result["flux"].data[29, 29], 1.626651e-09, rtol=1e-2)
-    assert_allclose(result["flux_err"].data[29, 29], 9.548939e-11, rtol=1e-2)
-    assert_allclose(result["flux_errp"].data[29, 29], 9.540259e-11, rtol=1e-2)
-    assert_allclose(result["flux_errn"].data[29, 29], 9.034366e-11, rtol=1e-2)
-    assert_allclose(result["flux_ul"].data[29, 29], 1.961776e-10, rtol=1e-2)
+    assert_allclose(result["flux"].data[29, 29], 1.351949e-09, rtol=1e-2)
+    assert_allclose(result["flux_err"].data[29, 29], 7.93751176e-11, rtol=1e-2)
+    assert_allclose(result["flux_errp"].data[29, 29], 7.9376134e-11, rtol=1e-2)
+    assert_allclose(result["flux_errn"].data[29, 29], 7.5180404579e-11, rtol=1e-2)
+    assert_allclose(result["flux_ul"].data[29, 29], 1.63222157e-10, rtol=1e-2)
     assert result["flux"].unit == u.Unit("cm-2s-1")
     assert result["flux_err"].unit == u.Unit("cm-2s-1")
     assert result["flux_ul"].unit == u.Unit("cm-2s-1")

gammapy/estimators/ts_map.py

@@ -9,10 +9,12 @@
 import scipy.optimize
 from astropy.coordinates import Angle
 from astropy.utils import lazyproperty
-from gammapy.datasets.map import MapEvaluator
 from gammapy.maps import Map, WcsGeom
+from gammapy.datasets.map import MapEvaluator
 from gammapy.modeling.models import (
-    PointSpatialModel, PowerLawSpectralModel, SkyModel, ConstantFluxSpatialModel
+    PointSpatialModel,
+    PowerLawSpectralModel,
+    SkyModel,
 )
 from gammapy.stats import (
     norm_bounds_cython,
@@ -21,6 +23,7 @@
 )
 from gammapy.utils.array import shape_2N, symmetric_crop_pad_width
 from .core import Estimator
+from .utils import estimate_exposure_reco_energy
 
 __all__ = ["TSMapEstimator"]
 
@@ -67,7 +70,9 @@ class TSMapEstimator(Estimator):
     Parameters
     ----------
     model : `~gammapy.modeling.model.SkyModel`
-        Source model kernel. If set to None, assume point source model, PointSpatialModel.
+        Source model kernel. If set to None,
+        assume spatail model: point source model, PointSpatialModel.
+        spectral model: PowerLawSpectral Model of index 2
     kernel_width : `~astropy.coordinates.Angle`
         Width of the kernel to use: the kernel will be truncated at this size
     n_sigma : int
@@ -126,15 +131,15 @@ def __init__(
         threshold=None,
         rtol=0.01,
         selection_optional="all",
-        n_jobs=None
+        n_jobs=None,
     ):
         self.kernel_width = Angle(kernel_width)
 
         if model is None:
             model = SkyModel(
                 spectral_model=PowerLawSpectralModel(),
                 spatial_model=PointSpatialModel(),
-                name="ts-kernel"
+                name="ts-kernel",
             )
 
         self.model = model
@@ -151,7 +156,7 @@ def __init__(
             n_sigma=self.n_sigma,
             n_sigma_ul=self.n_sigma_ul,
             selection_optional=selection_optional,
-            ts_threshold=threshold
+            ts_threshold=threshold,
         )
 
     def estimate_kernel(self, dataset):
@@ -204,44 +209,6 @@ def estimate_kernel(self, dataset):
             )
         return kernel
 
-    def estimate_exposure(self, dataset):
-        """Estimate exposure map in reco energy
-
-
-        Parameters
-        ----------
-        dataset : `~gammapy.datasets.MapDataset`
-            Input dataset.
-
-        Returns
-        -------
-        exposure : `Map`
-            Exposure map
-
-        """
-        # TODO: clean this up a bit...
-        models = dataset.models
-
-        model = SkyModel(
-            spectral_model=self.model.spectral_model,
-            spatial_model=ConstantFluxSpatialModel(),
-        )
-        model.apply_irf["psf"] = False
-
-        energy_axis = dataset.exposure.geom.get_axis_by_name("energy_true")
-        energy = energy_axis.edges
-
-        flux = model.spectral_model.integral(
-            emin=energy.min(), emax=energy.max()
-        )
-
-        self._flux_estimator.flux_ref = flux.to_value("cm-2 s-1")
-        dataset.models = [model]
-        npred = dataset.npred()
-        dataset.models = models
-        data = (npred.data / flux).to("cm2 s")
-        return npred.copy(data=data.value, unit=data.unit)
-
     def estimate_flux_default(self, dataset, kernel, exposure=None):
         """Estimate default flux map using a given kernel.
 
@@ -258,10 +225,11 @@ def estimate_flux_default(self, dataset, kernel, exposure=None):
             Approximate flux map.
         """
         if exposure is None:
-            exposure = self.estimate_exposure(dataset)
+            exposure = estimate_exposure_reco_energy(dataset, self.model.spectral_model)
 
         kernel = kernel / np.sum(kernel ** 2)
         flux = (dataset.counts - dataset.npred()) / exposure
+        flux.quantity = flux.quantity.to("1 / (cm2 s)")
         flux = flux.convolve(kernel)
         return flux.sum_over_axes()
 
@@ -289,9 +257,7 @@ def estimate_mask_default(dataset, kernel):
         slice_y = slice(kernel.shape[1] // 2, -kernel.shape[1] // 2 + 1)
         mask[slice_y, slice_x] = 1
 
-        mask &= dataset.mask_safe.reduce_over_axes(
-                func=np.logical_or, keepdims=False
-            )
+        mask &= dataset.mask_safe.reduce_over_axes(func=np.logical_or, keepdims=False)
 
         # in some image there are pixels, which have exposure, but zero
         # background, which doesn't make sense and causes the TS computation
@@ -362,7 +328,7 @@ def run(self, dataset):
         counts = dataset.counts
         background = dataset.npred()
 
-        exposure = self.estimate_exposure(dataset)
+        exposure = estimate_exposure_reco_energy(dataset, self.model.spectral_model)
 
         kernel = self.estimate_kernel(dataset)
 
@@ -377,7 +343,7 @@ def run(self, dataset):
             background=background.data.astype(float),
             kernel=kernel.data,
             flux=flux.data,
-            flux_estimator=self._flux_estimator
+            flux_estimator=self._flux_estimator,
         )
 
         x, y = np.where(np.squeeze(mask.data))
@@ -412,6 +378,7 @@ def run(self, dataset):
             m.data[j, i] = [_[name.replace("flux", "norm")] for _ in results]
             if "flux" in name:
                 m.data *= self._flux_estimator.flux_ref
+                m.quantity = m.quantity.to("1 / (cm2 s)")
             result[name] = m
 
         result["sqrt_ts"] = self.estimate_sqrt_ts(result["ts"])
@@ -441,6 +408,7 @@ class SimpleMapDataset:
         Kernel array
 
     """
+
     def __init__(self, model, counts, background, x_guess):
         self.model = model
         self.counts = counts
@@ -460,9 +428,7 @@ def npred(self, norm):
 
     def stat_sum(self, norm):
         """Stat sum"""
-        return cash_sum_cython(
-            self.counts.ravel(), self.npred(norm).ravel()
-        )
+        return cash_sum_cython(self.counts.ravel(), self.npred(norm).ravel())
 
     def stat_derivative(self, norm):
         """Stat derivative"""
@@ -473,7 +439,11 @@ def stat_derivative(self, norm):
     def stat_2nd_derivative(self, norm):
         """Stat 2nd derivative"""
         with np.errstate(invalid="ignore", divide="ignore"):
-            return (self.model ** 2 * self.counts / (self.background + norm * self.model) ** 2).sum()
+            return (
+                self.model ** 2
+                * self.counts
+                / (self.background + norm * self.model) ** 2
+            ).sum()
 
     @classmethod
     def from_arrays(cls, counts, background, exposure, flux, position, kernel):
@@ -486,17 +456,26 @@ def from_arrays(cls, counts, background, exposure, flux, position, kernel):
             counts=counts_cutout,
             background=background_cutout,
             model=kernel * exposure_cutout,
-            x_guess=x_guess
+            x_guess=x_guess,
         )
 
 
 # TODO: merge with `FluxEstimator`?
 class BrentqFluxEstimator(Estimator):
     """Single parameter flux estimator"""
+
     _available_selection_optional = ["errn-errp", "ul"]
     tag = "BrentqFluxEstimator"
 
-    def __init__(self, rtol, n_sigma, n_sigma_ul, selection_optional=None, max_niter=20, ts_threshold=None):
+    def __init__(
+        self,
+        rtol,
+        n_sigma,
+        n_sigma_ul,
+        selection_optional=None,
+        max_niter=20,
+        ts_threshold=None,
+    ):
         self.rtol = rtol
         self.n_sigma = n_sigma
         self.n_sigma_ul = n_sigma_ul
@@ -537,7 +516,9 @@ def estimate_best_fit(self, dataset):
         result["ts"] = (stat_null - stat) * np.sign(norm)
         result["norm"] = norm
         result["niter"] = niter
-        result["norm_err"] = np.sqrt(1 / dataset.stat_2nd_derivative(norm)) * self.n_sigma
+        result["norm_err"] = (
+            np.sqrt(1 / dataset.stat_2nd_derivative(norm)) * self.n_sigma
+        )
         result["stat"] = stat
         return result
 
@@ -548,7 +529,7 @@ def nan_result(self):
             "stat": np.nan,
             "norm_err": np.nan,
             "ts": np.nan,
-            "niter": 0
+            "niter": 0,
         }
 
         if "errn-errp" in self.selection_optional:
@@ -581,11 +562,7 @@ def ts_diff(x):
             warnings.simplefilter("ignore")
             try:
                 result_fit = scipy.optimize.brentq(
-                    ts_diff,
-                    min_norm,
-                    max_norm,
-                    maxiter=self.max_niter,
-                    rtol=self.rtol,
+                    ts_diff, min_norm, max_norm, maxiter=self.max_niter, rtol=self.rtol,
                 )
                 return (result_fit - norm) * factor
             except (RuntimeError, ValueError):
@@ -635,15 +612,7 @@ def run(self, dataset):
         return result
 
 
-def _ts_value(
-    position,
-    counts,
-    exposure,
-    background,
-    kernel,
-    flux,
-    flux_estimator
-):
+def _ts_value(position, counts, exposure, background, kernel, flux, flux_estimator):
     """Compute TS value at a given pixel position.
 
     Uses approach described in Stewart (2009).
@@ -675,6 +644,6 @@ def _ts_value(
         exposure=exposure,
         kernel=kernel * flux_estimator.flux_ref,
         position=position,
-        flux=flux
+        flux=flux,
     )
     return flux_estimator.run(dataset)