Repository: gammapy/gammapy
Pull request #5532 (2024-10-30 10:26:55)

Cleanup flux and significance maps  combination functions 
Cleanup flux and significance maps functions combination such as they have the same api.


Patch:
diff --git a/gammapy/estimators/utils.py b/gammapy/estimators/utils.py
index 1f50fa863..80c59609e 100644
--- a/gammapy/estimators/utils.py
+++ b/gammapy/estimators/utils.py
@@ -25,12 +25,15 @@ from gammapy.stats.utils import ts_to_sigma
 from .map.core import FluxMaps
 
 __all__ = [
-    "get_combined_significance_maps",
+    "combine_flux_maps",
+    "combine_significance_maps",
     "estimate_exposure_reco_energy",
     "find_peaks",
     "find_peaks_in_flux_map",
     "resample_energy_edges",
     "get_rebinned_axis",
+    "get_combined_flux_maps",
+    "get_combined_significance_maps",
     "compute_lightcurve_fvar",
     "compute_lightcurve_fpp",
     "compute_lightcurve_doublingtime",
@@ -779,14 +782,13 @@ def get_rebinned_axis(fluxpoint, axis_name="energy", method=None, **kwargs):
     return axis_new
 
 
-def get_combined_significance_maps(estimator, datasets):
-    """Compute excess and significance for a set of datasets.
-
+def combine_significance_maps(maps):
+    """Computes excess and significance for a set of datasets.
     The significance computation assumes that the model contains
     one degree of freedom per valid energy bin in each dataset.
-    This method implemented here is valid under the assumption
+    The method implemented here is valid under the assumption
     that the TS in each independent bin follows a Chi2 distribution,
-    then the sum of the TS also follows a Chi2 distribution (with the sum of degree of freedom).
+    then the sum of the TS also follows a Chi2 distribution (with the sum of the degrees of freedom).
 
     See, Zhen (2014): https://www.sciencedirect.com/science/article/abs/pii/S0167947313003204,
     Lancaster (1961): https://onlinelibrary.wiley.com/doi/10.1111/j.1467-842X.1961.tb00058.x
@@ -794,10 +796,8 @@ def get_combined_significance_maps(estimator, datasets):
 
     Parameters
     ----------
-    estimator : `~gammapy.estimators.ExcessMapEstimator` or `~gammapy.estimators.TSMapEstimator`
-        Excess Map Estimator or TS Map Estimator
-    dataset : `~gammapy.datasets.Datasets`
-        Datasets containing only `~gammapy.datasets.MapDataset`.
+    maps : list of `~gammapy.estimators.FluxMaps`
+        List of maps with the same geometry.
 
     Returns
     -------
@@ -807,27 +807,21 @@ def get_combined_significance_maps(estimator, datasets):
                 * "significance" : joint significance map.
                 * "df" : degree of freedom map (one norm per valid bin).
                 * "npred_excess" : summed excess map.
-                * "estimator_results" : dictionary containing the estimator results for each dataset.
-    """
-    from .map.excess import ExcessMapEstimator
-    from .map.ts import TSMapEstimator
+                * "estimator_results" : dictionary containing the flux maps computed for each dataset.
 
-    if not isinstance(estimator, (ExcessMapEstimator, TSMapEstimator)):
-        raise TypeError(
-            f"estimator type should be ExcessMapEstimator or TSMapEstimator), got {type(estimator)} instead."
-        )
+    See also
+    --------
+    get_combined_significance_maps : same method but computing the significance maps from estimators and datasets.
+
+    """
 
-    geom = datasets[0].counts.geom.to_image()
+    geom = maps[0].ts.geom.to_image()
     ts_sum = Map.from_geom(geom)
     ts_sum_sign = Map.from_geom(geom)
     npred_excess_sum = Map.from_geom(geom)
     df = Map.from_geom(geom)
 
-    results = dict()
-    for kd, d in enumerate(datasets):
-        result = estimator.run(d)
-        results[d.name] = result
-
+    for result in maps:
         df += np.sum(result["ts"].data > 0, axis=0)  # one dof (norm) per valid bin
         ts_sum += result["ts"].reduce_over_axes()
         ts_sum_sign += (
@@ -841,10 +835,59 @@ def get_combined_significance_maps(estimator, datasets):
         significance=significance,
         df=df,
         npred_excess=npred_excess_sum,
-        estimator_results=results,
+        estimator_results=maps,
     )
 
 
+def get_combined_significance_maps(estimator, datasets):
+    """Compute excess and significance for a set of datasets.
+
+    The significance computation assumes that the model contains
+    one degree of freedom per valid energy bin in each dataset.
+    This method implemented here is valid under the assumption
+    that the TS in each independent bin follows a Chi2 distribution,
+    then the sum of the TS also follows a Chi2 distribution (with the sum of degree of freedom).
+
+    See, Zhen (2014): https://www.sciencedirect.com/science/article/abs/pii/S0167947313003204,
+    Lancaster (1961): https://onlinelibrary.wiley.com/doi/10.1111/j.1467-842X.1961.tb00058.x
+
+
+    Parameters
+    ----------
+    estimator : `~gammapy.estimators.ExcessMapEstimator` or `~gammapy.estimators.TSMapEstimator`
+        Excess Map Estimator or TS Map Estimator
+    dataset : `~gammapy.datasets.Datasets`
+        Datasets containing only `~gammapy.datasets.MapDataset`.
+
+    Returns
+    -------
+    results : dict
+        Dictionary with entries:
+
+                * "significance" : joint significance map.
+                * "df" : degree of freedom map (one norm per valid bin).
+                * "npred_excess" : summed excess map.
+                * "estimator_results" : dictionary containing the flux maps computed for each dataset.
+
+    See also
+    --------
+    combine_significance_maps : same method but using directly the significance maps from estimators
+
+    """
+    from .map.excess import ExcessMapEstimator
+    from .map.ts import TSMapEstimator
+
+    if not isinstance(estimator, (ExcessMapEstimator, TSMapEstimator)):
+        raise TypeError(
+            f"estimator type should be ExcessMapEstimator or TSMapEstimator), got {type(estimator)} instead."
+        )
+
+    results = []
+    for dataset in datasets:
+        results.append(estimator.run(dataset))
+    return combine_significance_maps(results)
+
+
 def combine_flux_maps(
     maps, method="gaussian_errors", reference_model=None, dnde_scan_axis=None
 ):
@@ -857,7 +900,7 @@ def combine_flux_maps(
     ----------
     maps : list of `~gammapy.estimators.FluxMaps`
         List of maps with the same geometry.
-    method : {"gaussian_errors"}
+    method : str
         * gaussian_errors :
             Under the gaussian error approximation the likelihood is given by the gaussian distibution.
             The product of gaussians is also a gaussian so can derive dnde, dnde_err, and ts.
@@ -886,6 +929,12 @@ def combine_flux_maps(
     -------
     flux_maps : `~gammapy.estimators.FluxMaps`
         Joint flux map.
+
+
+    See also
+    --------
+    get_combined_flux_maps : same method but using directly the flux maps from estimators
+
     """
     gtis = [map_.gti for map_ in maps if map_.gti is not None]
     if np.any(gtis):
@@ -973,6 +1022,84 @@ def combine_flux_maps(
         )
 
 
+def get_combined_flux_maps(
+    estimator,
+    datasets,
+    method="gaussian_errors",
+    reference_model=None,
+    dnde_scan_axis=None,
+):
+    """Create a `~gammapy.estimators.FluxMaps` by combining a list of flux maps with the same geometry.
+
+    This assumes the flux maps are independent measurements of the same true value.
+    The GTI is stacked in the process.
+
+    Parameters
+    ----------
+    estimator : `~gammapy.estimators.ExcessMapEstimator` or `~gammapy.estimators.TSMapEstimator`
+        Excess Map Estimator or TS Map Estimator
+    dataset : `~gammapy.datasets.Datasets` or list of `~gammapy.datasets.MapDataset`
+        Datasets containing only `~gammapy.datasets.MapDataset`.
+    method : str
+        * gaussian_errors :
+            Under the gaussian error approximation the likelihood is given by the gaussian distibution.
+            The product of gaussians is also a gaussian so can derive dnde, dnde_err, and ts.
+        * distrib :
+            Likelihood profile approximation assuming that probabilities distributions for
+            flux points correspond to asymmetric gaussians and for upper limits to complementary error functions.
+            Use available quantities among dnde, dnde_err, dnde_errp, dnde_errn, dnde_ul, and ts.
+        * profile :
+            Sum the likelihood profile maps.
+            The flux maps must contains the `stat_scan` maps.
+
+        Default is "gaussian_errors" which is the faster but least accurate solution,
+        "distrib"  will be more accurate if dnde_errp and dnde_errn are available,
+        "profile"  will be even more accurate if "stat_scan" is available.
+
+    reference_model : `~gammapy.modeling.models.SkyModel`, optional
+        Reference model to use for conversions.
+        Default is None and is will use the reference_model of the first FluxMaps in the list.
+    dnde_scan_axis : `~gammapy.maps.MapAxis`, optional
+        Map axis providing the dnde values used to compute the profile.
+        If None, it will be derived from the first FluxMaps in the list. Default is None.
+        Used only if `method` is "distrib" or "profile".
+
+    Returns
+    -------
+    results : dict
+        Dictionary with entries:
+
+                * "flux_maps" : `gammapy.estimators.FluxMaps`
+                * "estimator_results" : dictionary containing the flux maps computed for each dataset.
+
+    See also
+    --------
+    combine_flux_maps : same method but using directly the flux maps from estimators
+
+    """
+    from .map.excess import ExcessMapEstimator
+    from .map.ts import TSMapEstimator
+
+    if not isinstance(estimator, (ExcessMapEstimator, TSMapEstimator)):
+        raise TypeError(
+            f"`estimator` type should be ExcessMapEstimator or TSMapEstimator), got {type(estimator)} instead."
+        )
+
+    results = []
+    for dataset in datasets:
+        results.append(estimator.run(dataset))
+
+    output = dict()
+    output["flux_maps"] = combine_flux_maps(
+        results,
+        method=method,
+        reference_model=reference_model,
+        dnde_scan_axis=dnde_scan_axis,
+    )
+    output["estimator_results"] = results
+    return output
+
+
 def _default_scan_map(flux_map, dnde_scan_axis=None):
     if dnde_scan_axis is None:
         dnde_scan_axis = MapAxis(


Test patch:
diff --git a/gammapy/estimators/map/tests/test_excess.py b/gammapy/estimators/map/tests/test_excess.py
index be8143f1f..1630e83ca 100644
--- a/gammapy/estimators/map/tests/test_excess.py
+++ b/gammapy/estimators/map/tests/test_excess.py
@@ -418,6 +418,10 @@ def test_joint_excess_map(simple_dataset):
     simple_dataset.models = [model]
     simple_dataset.npred()
 
+    simple_dataset2 = simple_dataset.copy()
+    simple_dataset2.models = [model]
+    simple_dataset2.npred()
+
     stacked_dataset = simple_dataset.copy(name="copy")
     stacked_dataset.counts *= 2
     stacked_dataset.exposure *= 2
@@ -430,7 +434,9 @@ def test_joint_excess_map(simple_dataset):
     assert_allclose(result["npred_excess"].data.sum(), 2 * 19733.602, rtol=1e-3)
     assert_allclose(result["sqrt_ts"].data[0, 10, 10], 5.960441, rtol=1e-3)
 
-    result = get_combined_significance_maps(estimator, [simple_dataset, simple_dataset])
+    result = get_combined_significance_maps(
+        estimator, [simple_dataset, simple_dataset2]
+    )
 
     assert_allclose(result["npred_excess"].data.sum(), 2 * 19733.602, rtol=1e-3)
     assert_allclose(result["significance"].data[10, 10], 5.618187, rtol=1e-3)
diff --git a/gammapy/estimators/map/tests/test_ts.py b/gammapy/estimators/map/tests/test_ts.py
index 16235b23f..9aea86272 100644
--- a/gammapy/estimators/map/tests/test_ts.py
+++ b/gammapy/estimators/map/tests/test_ts.py
@@ -10,6 +10,7 @@ from gammapy.estimators import TSMapEstimator
 from gammapy.estimators.utils import (
     approximate_profile_map,
     combine_flux_maps,
+    get_combined_flux_maps,
     get_combined_significance_maps,
     get_flux_map_from_profile,
 )
@@ -412,6 +413,14 @@ def test_ts_map_stat_scan(fake_dataset):
     assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
     assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
 
+    combined_results = get_combined_flux_maps(
+        estimator, [dataset, dataset.copy()], method="distrib"
+    )
+    combined_map = combined_results["flux_maps"]
+    assert len(combined_results["estimator_results"]) == 2
+    assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
+    assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
+
     combined_map = combine_flux_maps([maps, maps1], method="profile")
     assert_allclose(combined_map.ts.data, 2 * ts, rtol=1e-4)
     assert_allclose(combined_map.norm.data[success], norm[success], rtol=5e-2)
@@ -558,8 +567,10 @@ def test_with_TemplateSpatialModel():
 def test_joint_ts_map(fake_dataset):
     model = fake_dataset.models["source"]
     fake_dataset = fake_dataset.copy()
+    fake_dataset2 = fake_dataset.copy()
 
     fake_dataset.models = [model]
+    fake_dataset2.models = [model]
 
     estimator = TSMapEstimator(
         model=model, selection_optional=[], sum_over_energy_groups=True
@@ -570,7 +581,7 @@ def test_joint_ts_map(fake_dataset):
     assert_allclose(result["npred_excess"].data.sum(), 902.403647, rtol=1e-3)
     assert_allclose(result["sqrt_ts"].data[0, 10, 10], 1.360219, rtol=1e-3)
 
-    result = get_combined_significance_maps(estimator, [fake_dataset, fake_dataset])
+    result = get_combined_significance_maps(estimator, [fake_dataset, fake_dataset2])
 
     assert_allclose(result["npred_excess"].data.sum(), 2 * 902.403647, rtol=1e-3)
     assert_allclose(result["significance"].data[10, 10], 1.414529, rtol=1e-3)
@@ -581,7 +592,7 @@ def test_joint_ts_map(fake_dataset):
     estimator = TSMapEstimator(
         model=model, threshold=1, selection_optional="all", sum_over_energy_groups=True
     )
-    result = estimator.run([fake_dataset, fake_dataset.copy()])
+    result = estimator.run([fake_dataset, fake_dataset2])
     assert_allclose(result["sqrt_ts"].data[0, 10, 10], 1.92364, rtol=1e-3)