Psmap implementation in fermipy

docs/advanced/index.rst

@@ -7,6 +7,9 @@ Advanced Analysis Methods
 This page documents some of the more advanced methods and features
 available in Fermipy:
 
+* :ref:`psmap`: Generate a PS map quantifying the data/model agreement using the algorithm described in
+  Bruel P. (2021), A&A, 656, A81. (`doi:10.1051/0004-6361/202141553 <https://arxiv.org/pdf/2109.07443.pdf>`_).
+
 * :ref:`tsmap`: Generate a test statistic (TS) map for a new source
   centered at each spatial bin in the ROI.
 
@@ -52,6 +55,7 @@ available in Fermipy:
    curvature
    lightcurve
    extension
+   psmap
    tsmap
    tscube
    residmap

docs/advanced/psmap.rst

@@ -0,0 +1,137 @@
+.. _psmap:
+
+PS Map
+======
+
+:py:meth:`~fermipy.gtanalysis.GTAnalysis.psmap` quantifies the 3d data/model agreement by computing the PS
+at each spatial bin in the ROI according to the algorithm described in Bruel P. (2021), A&A, 656, A81.
+(`doi:10.1051/0004-6361/202141553 <https://arxiv.org/pdf/2109.07443.pdf>`_). For each spatial bin, the algorithm first computes
+the data and model count spectra integrated over an energy dependent region (following the PSF energy dependence)
+and then computes the p-value that the data count spectrum is compatible with the model count spectrum (using the likelihood statistics).
+The absolute value of PS is -log10(p-value) and its sign corresponds to the sign of the sum over the spectrum of the residuals in sigma.
+The algorithm also provides a PS map in sigma units.
+
+Caveat: the algorithm is currently not able to run in the context of a joint analysis with several event types.
+
+
+Examples
+--------
+
+In order to run psmap, the user must first compute the source model map using the :py:meth:`~fermipy.gtanalysis.GTAnalysis.write_model_map'
+function specifying the name of the model with the parameter ``model_name``.
+The :py:meth:`~fermipy.gtanalysis.GTAnalysis.psmap` will then be called with the same ``model_name``.
+
+.. code-block:: python
+
+   # Write the source model map (after performing the fit)
+   gta.write_model_map(model_name="model01")
+
+   # Generate the PS map
+   psmap = gta.psmap(model_name='model01', make_plots=True)
+
+
+Configuration
+-------------
+
+The default configuration of the method is controlled with the
+:ref:`config_psmap` section of the configuration file.  The default
+configuration can be overriden by passing the option as a *kwargs*
+argument to the method.
+
+If the analysis uses likelihood weights, the user can specify the likelihood weight file with the argument ``wmap``.
+
+The user can set the parameters defining the energy dependent region used in the count spectrum integration step.
+The integration radius is defined by the function sqrt(psfpar0^2*pow(E/psfpar1,-2*psfpar2) + psfpar3^2), with E in MeV.
+The default parameters (``psfpar0``=4.0,``psfpar1``=100,``psfpar2``=0.9,``psfpar3``=0.1) are optimized to look for point-source like deviations.
+When looking for extended deviations (~1deg scale), it is recommended to use (``psfpar0``=5.0,``psfpar1``=100,``psfpar2``=0.8,``psfpar3``=1.0).
+
+The user can set the energy range with the arguments ``emin`` and ``emax``. Depending on the energy binning, it can be optimal
+to rebin the count spectra thanks to the argument ``rebin``. For an analysis with 10 bins per decade, it is recommended to set ``rebin`` to 3.
+
+.. csv-table:: *psmap* Options
+   :header:    Option, Default, Description
+   :file: ../config/psmap.csv
+   :delim: tab
+   :widths: 10,10,80
+
+
+Output
+------
+
+:py:meth:`~fermipy.gtanalysis.GTAnalysis.psmap` returns a dictionary containing the following variables:
+
+============= ====================== =================================================================
+Key           Type                   Description
+============= ====================== =================================================================
+psmax         float                  maximum of the ps map
+coordname1    str                    Name of the coordinate of the ps maximum
+coordname2    str                    Name of the coordinaste of the ps maximum
+coordx        float                  Value of the X coordinate of the ps maximum
+coordy        float                  Value of the Y coordinate of the ps maximum
+ipix          int                    Value of the i pixel of the of the ps maximum
+jpix          int                    Value of the j pixel of the of the ps maximum
+wcs2d         WCS Keywords           WCS of the ps map
+psmap         np.array               PSMAP
+psmapsigma    np.array               PSMAP in sigma units
+name          str                    NAmke of the model
+ps_map        `~fermipy.skymap.Map`  WcsNDMap PSMAP
+pssigma_map   `~fermipy.skymap.Map`  WcsNDMap PSMAP in sigma units
+config        dict                   Dictionary of the input configuration
+file          str                    Name of the output file
+file_name     str                    Full path of the output file
+============= ====================== =================================================================
+
+The ``write_fits`` option can used to write the output to a FITS or numpy file. The value of the maximum of the PS map
+can be retrieved from the output dictionary:
+
+.. code-block:: python
+
+   print('PS maximum value=%.2f, at %s=%.2f, %s=%.2f' %(psmap['psmax'],
+                                                     psmap['coordname1'],float(psmap['coordx']),
+                                                     psmap['coordname2'],float(psmap['coordy'])))
+
+   PS maximum value=3.85, at GLON-AIT=86.75, GLAT-AIT=38.62
+
+Diagnostic plots can be generated by setting ``make_plots=True`` or by
+passing the output dictionary to `~fermipy.plotting.AnalysisPlotter.make_psmap_plots`:
+
+.. code-block:: python
+
+   psmap = gta.psmap(model_name='model01', make_plots=True)
+   //equivalent to:
+   gta.plotter.make_tsmap_plots(psmap, roi=gta.roi)
+
+This will generate the following plots:
+
+* ``image_psmap`` : Map of PS values.  The color map is truncated at
+  5 sigma with isocontours at 3,4,5 PS intervals indicating values
+  above this threshold.
+
+* ``image_pssigma`` : Map of PS values converted in sigma. The color map is truncated at
+  5 sigma with isocontours at 3,4,5 PS intervals indicating values
+  above this threshold.
+
+* ``image_ps_hist`` : Histogram of PS values for all points in the
+  map. Overplotted is the reference distribution for a gaussian with mean 0 and sigma=1.
+
+.. |image_psmap| image:: model01_psmap_psmap.png
+   :width: 100%
+
+.. |image_pssigma| image:: model01_psmap_pssigma.png
+   :width: 100%
+
+.. |image_ps_hist| image:: model01_psmap_ps_hist.png
+   :width: 100%
+
+.. csv-table::
+   :header: PS Map, Sigma (PS) Map, PS Histogram
+   :widths: 33, 33, 33
+
+   |image_psmap|, |image_pssigma|, |image_ps_hist|
+
+
+Reference/API
+-------------
+
+.. automethod:: fermipy.gtanalysis.GTAnalysis.psmap
+   :noindex:

docs/changelog.rst

@@ -6,8 +6,19 @@ Changelog
 This page is a changelog for releases of Fermipy.  You can also browse
 releases on `Github <https://github.com/fermiPy/fermipy/releases>`_.
 
+1.2.3 (03/25/2024)
+------------------
+* 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 Documentation
+
+
+1.2.2 (01/21/2024)
+------------------
+* fix the dependence of scipy due to gammapy
+
 1.2.1 (12/08/2023)
-----------------
+------------------
 * Small bug fixes.
 * pinned astropy<6
 

docs/conf.py

@@ -17,22 +17,21 @@
 #import sphinx_bootstrap_theme
 
 if sys.version_info.major == 2:
-    import mock
-    from mock import Mock as MagicMock
+    from mock import Mock
 else:
-    from unittest import mock
-    from unittest.mock import Mock as MagicMock
+    from unittest.mock import Mock
+
+try:
+    import GtApp # has fermitools, no mock needed
+    MOCK_MODULES = []##'pyLikelihood','pyIrfLoader']
+except:
+    # Mock modules needed
+    MOCK_MODULES = ['pyLikelihood','GtApp',
+                    'BinnedAnalysis','SrcModel','AnalysisBase',
+                    'LikelihoodState','pyIrfLoader',
+                    'SummedLikelihood','FluxDensity']
 
-class Mock(MagicMock):
-    @classmethod
-    def __getattr__(cls, name):
-        return Mock()
 
-#MOCK_MODULES = ['pyLikelihood','pyIrfLoader',
-#                'BinnedAnalysis','UnbinnedAnalysis','SrcModel','AnalysisBase',
-#                'SummedLikelihood','FluxDensity','LikelihoodState',
-#                'GtApp']
-MOCK_MODULES = ['pyLikelihood','pyIrfLoader']
 
 sys.modules.update((mod_name, Mock()) for mod_name in MOCK_MODULES)
 
@@ -97,7 +96,7 @@ def __getattr__(cls, name):
 # General information about the project.
 project = u'Fermipy'
 author = u'Fermipy Developers'
-copyright = u'2016-2022, ' + author
+copyright = u'2016-2024, ' + author
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
@@ -113,7 +112,7 @@ def __getattr__(cls, name):
 #
 # This is also used if you do content translation via gettext catalogs.
 # Usually you set "language" from the command line for these cases.
-language = None
+language = 'en'
 
 # There are two options for replacing |today|: either, you set today to some
 # non-false value, then it is used:
@@ -156,10 +155,10 @@ def __getattr__(cls, name):
 
 on_rtd = os.environ.get('READTHEDOCS', None) == 'True'
 
-if not on_rtd:  # only import and set the theme if we're building docs locally
-    import sphinx_rtd_theme
-    html_theme = 'sphinx_rtd_theme'
-    html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
+#if not on_rtd:  # only import and set the theme if we're building docs locally
+import sphinx_rtd_theme
+html_theme = 'sphinx_rtd_theme'
+html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # The theme to use for HTML and HTML Help pages.  See the documentation for
 # a list of builtin themes.

docs/config.rst

@@ -433,6 +433,21 @@ about using this method see the :ref:`findsources` page.
    :delim: tab
    :widths: 10,10,80
 
+.. _config_psmap:
+
+psmap
+-----
+
+The options in *psmap* control the default behavior of the
+`~fermipy.gtanalysis.GTAnalysis.psmap` method.  For more information
+about using this method see the :ref:`psmap` page.
+
+.. csv-table:: *psmap* Options
+   :header:    Option, Default, Description
+   :file: config/psmap.csv
+   :delim: tab
+   :widths: 10,10,80
+
 .. _config_tsmap:
 
 tsmap

docs/config/psmap.csv

@@ -0,0 +1,20 @@
+``chatter``	1	output verbosity
+``emax``	1000000000.0	maximum energy/MeV
+``emin``	1.0	minimum energy/MeV
+``fixedradius``	-1.0	Fixed radius
+``ipix``	-1	number of pixel i axis
+``jpix``	-1	number of pixel j axis
+``make_plots``	False	Generate diagnostic plots.
+``maxpoissoncount``	100	Maximum number of poisson counts
+``model_name``	None	Model name
+``nbinpdf``	50	Number of bin of the PSF
+``outfile``	psmap.fits	outfile name
+``prob_epsilon``	1e-07	precision parameter
+``psfpar0``	4.0	PSF parameter 0
+``psfpar1``	100	PSF parameter 1
+``psfpar2``	0.9	PSF parameter 2
+``psfpar3``	0.1	PSF parameter 3
+``rebin``	1	Rebin
+``scaleaxis``	20	SCale axis
+``wmap``		weight 3d map
+``write_fits``	True	Write the output to a FITS file.

docs/config/sed.csv

@@ -5,6 +5,7 @@
 ``free_radius``	None	Free normalizations of background sources within this angular distance in degrees from the source of interest.  If None then no sources will be freed.
 ``make_plots``	False	Generate diagnostic plots.
 ``ul_confidence``	0.95	Confidence level for flux upper limit.
+``ul_ts_threshold``	4	Minimum threshold of TS for displaying flux point below this vlaue an Upper Limit is calculated.
 ``use_local_index``	False	Use a power-law approximation to the shape of the global spectrum in each bin.  If this is false then a constant index set to `bin_index` will be used.
 ``write_fits``	True	Write the output to a FITS file.
 ``write_npy``	True	Write the output dictionary to a numpy file.

docs/developer.rst

@@ -77,6 +77,7 @@ The following are steps for creating a new release:
 .. code-block:: bash
 
    $ python setup.py sdist upload -r pypi
+   $ twine upload dist/fermipy-XX.YY.ZZ.tar.gz
 
 6. Create a new release on conda-forge by opening a PR on the
    `fermipy-feedstock

docs/index.rst

@@ -27,7 +27,7 @@ performing common analysis tasks:
 
 * Extracting a spectral energy distribution (SED) of a source.
 
-* Generating TS and residual maps for a region of interest.
+* Generating PS/TS and residual maps for a region of interest.
 
 * Finding new source candidates.
 

docs/requirements.txt

@@ -5,3 +5,4 @@ astropy
 pyyaml
 healpy
 gammapy
+sphinx_rtd_theme

fermipy/defaults.py

@@ -332,6 +332,31 @@ def make_attrs_class(typename, d):
     'init_lambda': (0, 'Initial value of damping parameter for newton step size calculation.   A value of zero disables damping.', float),
 }
 
+# PS map
+psmap = {
+    'model_name':(None, 'Model name',str),
+    'wmap':('', 'weight 3d map',str),
+    'outfile':('psmap.fits', 'outfile name',str),
+    'fixedradius':(-1.0, 'Fixed radius',float),
+    'psfpar0':(4.0, "PSF parameter 0", float),
+    'psfpar1':(100,"PSF parameter 1", float),
+    'psfpar2':(0.9,"PSF parameter 2", float),
+    'psfpar3':(0.1,"PSF parameter 3", float),
+    'maxpoissoncount':(100, "Maximum number of poisson counts", float),
+    'prob_epsilon':(1e-7, 'precision parameter', float),
+    'nbinpdf':(50, "Number of bin of the PSF", int),
+    'scaleaxis':(20,"SCale axis", float),
+    'emin':(1.0, "minimum energy/MeV", float),
+    'emax':(1e9, "maximum energy/MeV", float),
+    'chatter':(1, "output verbosity", int),
+    'ipix':(-1, "number of pixel i axis", int),
+    'jpix':(-1, "number of pixel j axis", int),
+    'rebin':(1, "Rebin", int),
+    'make_plots': common['make_plots'],
+    'write_fits': common['write_fits'],
+}
+
+
 # Options for Source Finder
 sourcefind = {
     'model': common['model'],
@@ -412,6 +437,9 @@ def make_attrs_class(typename, d):
     'free_radius': common['free_radius'],
     'free_pars': (None, 'Set the parameters of the source of interest that will be freed when performing '
                   'the global fit.  By default all parameters will be freed.', list),
+    'ul_ts_threshold': (4, 'Minimum threshold of TS for displaying flux point '
+                        'below this vlaue an Upper Limit is calculated.', float),
+
     'ul_confidence': (0.95, 'Confidence level for flux upper limit.',
                       float),
     'cov_scale': (3.0, 'Scale factor that sets the strength of the prior on nuisance '