Fix issue 669

docs/changelog.rst

@@ -6,6 +6,15 @@ Changelog
 This page is a changelog for releases of Fermipy.  You can also browse
 releases on `Github <https://github.com/fermiPy/fermipy/releases>`_.
 
+1.4.2 (unreleased)
+------------------
+* Fixed ROI restore consistency in `~fermipy.gtanalysis.GTAnalysis.create` and
+  `~fermipy.gtanalysis.GTAnalysis.load_roi`: when ``edisp: true`` is set and
+  ``edisp_disable`` lists diffuse sources, reloading a saved ROI now produces
+  model counts and residual maps identical to those computed before saving.
+* Added support and tests for the new FL16Y source list
+* Fixed some inconsistencies and bugs in the skymap library
+
 1.4.1 (02/25/2026)
 ------------------
 * Added tests for macos-15-intel and improved fit stability on macos
@@ -35,15 +44,16 @@ releases on `Github <https://github.com/fermiPy/fermipy/releases>`_.
 * Added PS map implementing code from Philippe Bruel in `~fermipy.gtanalysis.GTAnalysis.psmap`
 * Added PS map visualition in `~fermipy.gtanalysis.plotter.make_psmap_plots`
 * Added `PS map <https://fermipy.readthedocs.io/en/master/advanced/psmap.html>`_ to the documentation
+* Implemented Steve calibration to fix the residmap problem
 
 1.2.2 (01/21/2024)
 ------------------
-* fix the dependence of scipy due to gammapy
+* Fixed the dependence of scipy due to gammapy
 
 1.2.1 (12/08/2023)
 ------------------
 * Small bug fixes.
-* pinned astropy<6
+* Temporarily pinned `astropy<6`
 
 1.2 (09/21/2022)
 ----------------

fermipy/gtanalysis.py

@@ -7,6 +7,7 @@
 import logging
 import tempfile
 import filecmp
+import xml.etree.ElementTree as ET
 import time
 import json
 from pathlib import Path
@@ -5487,8 +5488,40 @@ def _create_binned_analysis(self, xmlfile=None, **kwargs):
         set_edisp_kwargs('gtlike', self.config, kw)
         self.logger.debug(kw)
 
-        self._like = BinnedAnalysis(binnedData=self._obs,
-                                    **utils.unicode_to_str(kw))
+        # When edisp is enabled globally and some sources have edisp disabled,
+        # loading an XML that contains apply_edisp='false' on those sources causes
+        # BinnedAnalysis to set edisp_val=0 instead of -1, leading to incorrect
+        # model counts after load_roi.  Strip the attribute before construction so
+        # the constructor always starts with edisp_val=-1; set_edisp_flag below
+        # will then correctly mark those sources as no-edisp.
+        edisp_disable = [s for s in self.config['gtlike']['edisp_disable']
+                         if self.roi.has_source(s)]
+        tmp_srcmdl = None
+        if edisp_disable and kw.get('edisp_bins', 0) != 0:
+            try:
+                tree = ET.parse(kw['srcModel'])
+                root = tree.getroot()
+                stripped = False
+                for src in root.findall('source'):
+                    if src.get('name') in edisp_disable:
+                        spectrum = src.find('spectrum')
+                        if spectrum is not None and 'apply_edisp' in spectrum.attrib:
+                            del spectrum.attrib['apply_edisp']
+                            stripped = True
+                if stripped:
+                    fd, tmp_srcmdl = tempfile.mkstemp(suffix='.xml')
+                    os.close(fd)
+                    tree.write(tmp_srcmdl, xml_declaration=True, encoding='unicode')
+                    kw['srcModel'] = tmp_srcmdl
+            except Exception:
+                pass
+
+        try:
+            self._like = BinnedAnalysis(binnedData=self._obs,
+                                        **utils.unicode_to_str(kw))
+        finally:
+            if tmp_srcmdl and os.path.exists(tmp_srcmdl):
+                os.unlink(tmp_srcmdl)
 
 #        print(self.like.logLike.use_single_fixed_map())
 #        self.like.logLike.set_use_single_fixed_map(False)

fermipy/tests/test_gtanalysis.py

@@ -202,6 +202,54 @@ def test_gtanalysis_residmap(create_diffuse_dir, create_draco_analysis):
                  make_plots=True)
 
 
+def test_gtanalysis_create_restore_residmap_consistency(create_diffuse_dir,
+                                                        create_draco_analysis):
+    gta = create_draco_analysis
+    gta.load_roi('fit1')
+    gta.optimize()
+    gta.write_roi('restore_consistency_test', make_plots=False)
+
+    ref = gta.residmap(model={}, make_plots=False,
+                       prefix='restore_consistency_ref')
+    ref_mean = np.nanmean(ref['excess'].data)
+    ref_model_means = np.array([np.nanmean(c.model_counts_map().data)
+                                for c in gta.components])
+
+    infile = os.path.join(gta.workdir, 'restore_consistency_test.npy')
+    gta_reloaded = gtanalysis.GTAnalysis.create(infile)
+    out = gta_reloaded.residmap(model={}, make_plots=False,
+                                prefix='restore_consistency_cmp')
+    out_mean = np.nanmean(out['excess'].data)
+    out_model_means = np.array([np.nanmean(c.model_counts_map().data)
+                                for c in gta_reloaded.components])
+
+    assert_allclose(out_mean, ref_mean, rtol=1E-3, atol=1E-10)
+    assert_allclose(out_model_means, ref_model_means, rtol=1E-3, atol=1E-10)
+
+
+def test_gtanalysis_load_roi_same_instance_consistency(
+        create_diffuse_dir, create_draco_analysis):
+    gta = create_draco_analysis
+    gta.load_roi('fit1')
+    gta.optimize()
+    gta.write_roi('load_roi_same_instance_test', make_plots=False)
+
+    ref = gta.residmap(model={}, make_plots=False,
+                       prefix='load_roi_same_instance_ref')
+    ref_mean = np.nanmean(ref['excess'].data)
+    ref_model_means = np.array([np.nanmean(c.model_counts_map().data)
+                                for c in gta.components])
+
+    gta.load_roi('load_roi_same_instance_test')
+    out = gta.residmap(model={}, make_plots=False,
+                       prefix='load_roi_same_instance_cmp')
+    out_mean = np.nanmean(out['excess'].data)
+    out_model_means = np.array([np.nanmean(c.model_counts_map().data)
+                                for c in gta.components])
+
+    assert_allclose(out_mean, ref_mean, rtol=1E-3, atol=1E-10)
+    assert_allclose(out_model_means, ref_model_means, rtol=1E-3, atol=1E-10)
+
 
 #@requires_git_version('02-00-00')
 def test_gtanalysis_find_sources(create_diffuse_dir, create_draco_analysis):