Repository: gammapy/gammapy
Pull request #5161 (2024-03-09 10:21:36)

Flux maps combination
Create a new `FluxMaps` by combining a list of flux maps with the same geometry.
Under the gaussian error approximation the likelihood is given by the gaussian distribution. The product of gaussians is also a gaussian so can derive dnde, dnde_err, and ts for the combined map.


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/gammapy/gammapy/pull/5161?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 75.26%. Comparing base [(`3bcae35`)](https://app.codecov.io/gh/gammapy/gammapy/commit/3bcae3569acc5452d83927c11ee05c258df7aadc?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy) to head [(`c03e3ea`)](https://app.codecov.io/gh/gammapy/gammapy/pull/5161?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).
> Report is 5 commits behind head on main.


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #5161      +/-   ##
==========================================
+ Coverage   75.22%   75.26%   +0.03%     
==========================================
  Files         234      234              
  Lines       34840    34883      +43     
==========================================
+ Hits        26210    26253      +43     
  Misses       8630     8630              
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/gammapy/gammapy/pull/5161?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=gammapy).

bkhelifi: Hello @QRemy,
Thanks a lot for this nice utility function. I have two questions:
- maybe one could add a check and a test of maps with different geometries. This might happen frequently
- I am wondering whether one could not extend this function for non-gaussian errors, when the scan profile is present. Maybe one could make that later into an other PR, and add a ToDo in this one....
QRemy: > * I am wondering whether one could not extend this function for non-gaussian errors, when the scan profile is present. Maybe one could make that later into an other PR, and add a ToDo in this one....

Yes I had that in mind but first we need to implement the stat_scan option ont the ExcessMapEstimator and TSMapEstimator.
registerrier: From the dev call, let's try to find a more explicit name for the function. It could be possibly an independent function rather than a classmethod.
registerrier: There is a merge conflict here. Can you rebase @QRemy ?
QRemy: 
> I simply propose to call the function `average_flux_maps` as this is more explicit. Joint flux map would rather be applied to flux maps extracted jointly from a group of datasets.

Maybe `combine_flux_maps` as we have `get_combined_significance_maps`. I think the term combined is often used when combining likehoods from different experiment a posteriori, especially in cosmology like in  https://arxiv.org/pdf/1807.06209 .


Patch:
diff --git a/gammapy/estimators/utils.py b/gammapy/estimators/utils.py
index 9a6a178db..bb2cd2fe6 100644
--- a/gammapy/estimators/utils.py
+++ b/gammapy/estimators/utils.py
@@ -1,6 +1,7 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 import numpy as np
 import scipy.ndimage
+from scipy import stats
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from astropy.table import Table
@@ -753,3 +754,95 @@ def get_combined_significance_maps(estimator, datasets):
         npred_excess=npred_excess_sum,
         estimator_results=results,
     )
+
+
+def combine_flux_maps(maps, method="gaussian_errors", reference_model=None):
+    """Create a 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
+    ----------
+    maps : list of `~gammapy.estimators.FluxMaps`
+        List of maps with the same geometry.
+    method : {"gaussian_errors"}
+        * 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.
+    reference_model : `~gammapy.modeling.models.SkyModel`, optional
+        Reference model to use for conversions.
+        If None, use the reference_model of the first FluxMaps in the list.
+
+    Returns
+    -------
+    flux_maps : `~gammapy.estimators.FluxMaps`
+        Joint flux map.
+    """
+    if method != "gaussian_errors":
+        raise ValueError(
+            f'Invalid method : {method}, available methods are : "gaussian_errors"'
+        )
+
+    if isinstance(maps, FluxMaps):
+        geom = maps.dnde.sum_over_axes().geom
+        means = list(maps.dnde.copy().iter_by_image(keepdims=True))
+        sigmas = list(maps.dnde_err.copy().iter_by_image(keepdims=True))
+        mean = Map.from_geom(geom, data=means[0].data, unit=means[0].unit)
+        sigma = Map.from_geom(geom, data=sigmas[0].data, unit=sigmas[0].unit)
+        gti = maps.gti
+        meta = maps.meta
+        if reference_model is None:
+            reference_model = maps.reference_model
+    else:
+        means = [map_.dnde.copy() for map_ in maps]
+        sigmas = [map_.dnde_err.copy() for map_ in maps]
+        mean = means[0]
+        sigma = sigmas[0]
+
+        gtis = [map_.gti for map_ in maps if map_.gti is not None]
+        if np.any(gtis):
+            gti = gtis[0].copy()
+            for k in range(1, len(gtis)):
+                gti.stack(gtis[k])
+        else:
+            gti = None
+
+        if reference_model is None:
+            reference_model = maps[0].reference_model
+
+        # TODO : change this once we have stackable metadata objets
+        metas = [map_.meta for map_ in maps if map_.meta is not None]
+        meta = {}
+        if np.any(metas):
+            for data in metas:
+                meta.update(data)
+
+    for k in range(1, len(means)):
+
+        mean_k = means[k].quantity.to_value(mean.unit)
+        sigma_k = sigmas[k].quantity.to_value(sigma.unit)
+
+        mask_valid = np.isfinite(mean) & np.isfinite(sigma) & (sigma.data != 0)
+        mask_valid_k = np.isfinite(mean_k) & np.isfinite(sigma_k) & (sigma_k != 0)
+        mask = mask_valid & mask_valid_k
+        mask_k = ~mask_valid & mask_valid_k
+
+        mean.data[mask] = (
+            (mean.data * sigma_k**2 + mean_k * sigma.data**2)
+            / (sigma.data**2 + sigma_k**2)
+        )[mask]
+        sigma.data[mask] = (
+            sigma.data * sigma_k / np.sqrt(sigma.data**2 + sigma_k**2)
+        )[mask]
+
+        mean.data[mask_k] = mean_k[mask_k]
+        sigma.data[mask_k] = sigma_k[mask_k]
+
+    ts = -2 * np.log(
+        stats.norm.pdf(0, loc=mean, scale=sigma)
+        / stats.norm.pdf(mean, loc=mean, scale=sigma)
+    )
+    ts = Map.from_geom(mean.geom, data=ts)
+
+    kwargs = dict(sed_type="dnde", gti=gti, reference_model=reference_model, meta=meta)
+    return FluxMaps.from_maps(dict(dnde=mean, dnde_err=sigma, ts=ts), **kwargs)


Test patch:
diff --git a/gammapy/estimators/map/tests/test_core.py b/gammapy/estimators/map/tests/test_core.py
index 0311b8735..57ea51745 100644
--- a/gammapy/estimators/map/tests/test_core.py
+++ b/gammapy/estimators/map/tests/test_core.py
@@ -2,12 +2,14 @@
 import pytest
 import numpy as np
 from numpy.testing import assert_allclose
+from scipy import stats
 import astropy.units as u
 from astropy.coordinates import SkyCoord
 from astropy.io import fits
 from astropy.time import Time
 from gammapy.data import GTI
 from gammapy.estimators import FluxMaps
+from gammapy.estimators.utils import combine_flux_maps
 from gammapy.maps import MapAxis, Maps, RegionGeom, TimeMapAxis, WcsNDMap
 from gammapy.modeling.models import (
     LogParabolaSpectralModel,
@@ -210,6 +212,56 @@ def test_map_properties(map_flux_estimate):
     assert_allclose(fe.eflux_ul.quantity.value[:, 2, 2], [4.60517e-10, 4.60517e-10])
 
 
+def test_combine_flux_maps(map_flux_estimate, wcs_flux_map, reference_model):
+    start = u.Quantity([1, 2], "min")
+    stop = u.Quantity([1.5, 2.5], "min")
+    gti = GTI.create(start, stop)
+
+    model = SkyModel(PowerLawSpectralModel(amplitude="1e-10 cm-2s-1TeV-1", index=2))
+    fe = FluxMaps(data=map_flux_estimate, reference_model=model, gti=gti)
+    iE = 0
+    energy = fe.geom.axes[0].center[iE]
+
+    fe_new = combine_flux_maps([fe, fe])
+    assert_allclose(fe_new.dnde.quantity, fe.dnde.quantity)
+    assert_allclose(
+        fe_new.dnde_err.quantity.value, fe.dnde_err.quantity.value / np.sqrt(2)
+    )
+
+    fe_new = combine_flux_maps(
+        [fe, fe], reference_model=FluxMaps.reference_model_default.spectral_model
+    )
+    assert_allclose(fe_new.dnde.quantity, fe.dnde.quantity)
+    assert_allclose(
+        fe_new.dnde_err.quantity.value, fe.dnde_err.quantity.value / np.sqrt(2)
+    )
+
+    ratio = model.spectral_model(
+        energy
+    ) / FluxMaps.reference_model_default.spectral_model(energy)
+
+    assert_allclose(
+        fe_new.norm.quantity.value[iE, :, :] / fe.norm.quantity.value[iE, :, :], ratio
+    )
+    assert fe_new.meta == fe.meta
+
+    fe = FluxMaps(wcs_flux_map, reference_model)
+    fe_new = combine_flux_maps([fe, fe], reference_model=reference_model)
+    assert_allclose(fe_new.dnde.quantity, fe.dnde.quantity)
+    assert_allclose(
+        fe_new.dnde_err.quantity.value, fe.dnde_err.quantity.value / np.sqrt(2)
+    )
+    assert_allclose(fe_new.norm.quantity, fe.norm.quantity)
+
+    ts = -2 * np.log(
+        stats.norm.pdf(0, loc=fe.dnde.data, scale=fe.dnde_err.data)
+        / stats.norm.pdf(fe.dnde.data, loc=fe.dnde.data, scale=fe.dnde_err.data)
+    )
+    assert_allclose(fe_new.ts, ts * 2)
+
+    fe_new = combine_flux_maps(fe)
+
+
 def test_flux_map_properties(wcs_flux_map, reference_model):
     fluxmap = FluxMaps(wcs_flux_map, reference_model)