Merge branch 'develop' into MAP_RichardsonLucy

.github/workflows/unit_tests.yml

@@ -27,6 +27,8 @@ jobs:
           pip install pytest pytest-cov
           pytest tests --junitxml=junit/test-results.xml --cov=cosipy --cov-report=xml --cov-report=html
       - name: Codecov
-        uses: codecov/[email protected]
-        env:
-          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}            
+        uses: codecov/codecov-action@v4
+        with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: true
+

cosipy/image_deconvolution/allskyimage.py

@@ -39,7 +39,7 @@ def __init__(self,
                  coordsys = 'galactic',
                  label_image = 'lb',
                  label_energy = 'Ei',
-                 unit = 1 / u.s / u.cm**2 / u.sr
+                 unit = '1/(s*cm*cm*sr)'
                  ):
 
         if energy_edges.unit != u.keV:
@@ -175,18 +175,21 @@ def set_values_from_extendedmodel(self, extendedmodel):
 
         self[:] = np.tensordot(normalized_map, integrated_flux.contents, axes = 0) 
 
-    def smoothing(self, fwhm = 0.0 * u.deg, sigma = None):
+    def smoothing(self, fwhm = None, sigma = None):
         """
         Smooth a map with a Gaussian filter
 
         Parameters
         ----------
         fwhm: :py:class:`astropy.units.quantity.Quantity`
-            The FWHM of the Gaussian (with a unit of deg or rad).
+            The FWHM of the Gaussian (with a unit of deg or rad). Default: 0 deg
         sigma: :py:class:`astropy.units.quantity.Quantity`
             The sigma of the Gaussian (with a unit of deg or rad). Override fwhm.
         """
 
+        if fwhm is None:
+            fwhm = 0.0 * u.deg
+
         if sigma is not None:
             fwhm = 2.354820 * sigma
 

cosipy/image_deconvolution/dataIF_COSI_DC2.py

@@ -5,6 +5,8 @@
 import logging
 logger = logging.getLogger(__name__)
 
+import warnings
+
 from histpy import Histogram, Axes
 
 from cosipy.response import FullDetectorResponse
@@ -100,7 +102,7 @@ def _modify_axes(self):
         Modify the axes of data. This method will be removed in the future.
         """
 
-        logger.warning("Note that _modify_axes() in DataIF_COSI_DC2 was implemented for a temporary use. It will be removed in the future.")
+        warnings.warn("Note that _modify_axes() in DataIF_COSI_DC2 was implemented for a temporary use. It will be removed in the future.", DeprecationWarning)
 
         if self._coordsys_conv_matrix is None:
             axis_name = ['Em', 'Phi', 'PsiChi']

cosipy/image_deconvolution/deconvolution_algorithm_base.py

@@ -59,7 +59,11 @@ def __init__(self, initial_model, dataset, mask, parameter):
         logger.debug(f'dict_dataset_indexlist_for_bkg_models: {self.dict_dataset_indexlist_for_bkg_models}')
 
         # minimum flux
-        self.minimum_flux = parameter.get('minimum_flux:value', 0.0) * u.Unit(parameter.get('minimum_flux:unit', initial_model.unit))
+        self.minimum_flux = parameter.get('minimum_flux:value', 0.0)
+
+        minimum_flux_unit = parameter.get('minimum_flux:unit', initial_model.unit)
+        if minimum_flux_unit is not None:
+            self.minimum_flux = self.minimum_flux*u.Unit(minimum_flux_unit)
 
         # parameters of the iteration
         self.iteration_count = 0

cosipy/image_deconvolution/model_base.py

@@ -62,7 +62,7 @@ def mask_pixels(self, mask, fill_value = 0):
         fill_value: float or :py:class:`astropy.units.quantity.Quantity`
         """
 
-        if not isinstance(fill_value, u.quantity.Quantity):
+        if not isinstance(fill_value, u.quantity.Quantity) and self.unit is not None:
             fill_value *= self.contents.unit
 
         model_new = copy.deepcopy(self)

cosipy/response/FullDetectorResponse.py

@@ -11,6 +11,9 @@
 import numpy as np
 import mhealpy as hp
 from mhealpy import HealpixBase, HealpixMap
+
+from scipy.special import erf
+
 from yayc import Configurator
 
 from scoords import SpacecraftFrame, Attitude
@@ -200,9 +203,12 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
                     if norm =="Linear" :
                         emin = int(line[2])
                         emax = int(line[3])
-
-
-
+
+                    if norm == "Gaussian" :
+                        Gauss_mean = float(line[2])   
+                        Gauss_sig = float(line[3])
+                        Gauss_cutoff = float(line[4])          
+
                 elif key == "MS":
                     if line[1] == "true" :
                         sparse = True
@@ -243,7 +249,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
                     break
 
         #check if the type of spectrum is known
-        assert norm=="powerlaw" or norm=="Mono" or norm=="Linear","unknown normalisation !" 
+        assert norm=="powerlaw" or norm=="Mono" or norm=="Linear" or norm=="Gaussian",f"unknown normalisation ! {norm}" 
 
         #check the number of simulated events is not 0
         assert nevents_sim != 0,"number of simulated events is 0 !" 
@@ -256,7 +262,6 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
         else :
             print("Sparse matrice ? {0}".format(sparse))
         edges = ()
-        #print(axes_edges)
 
         for axis_edges, axis_type in zip(axes_edges, axes_types):
 
@@ -412,6 +417,16 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
         #print(ewidth)
         #print(ecenters)
 
+        #if we have one single bin, treat the gaussian norm like the mono one
+        #also check that the gaussian spectrum is fully contained in that bin 
+        if len(ewidth) == 1 and norm == "Gaussian":
+            edges = dr.axes['Ei'].edges
+            gauss_int = 0.5 * (1 + erf( (edges[0]-Gauss_mean)/(4*np.sqrt(2)) ) ) + 0.5 * (1 + erf( (edges[1]-Gauss_mean)/(4*np.sqrt(2)) ) )
+
+            assert gauss_int == 1, "The gaussian spectrum is not fully contained in this single bin !"
+            print("Only one bin so we will use the Mono normalisation")
+            norm ="Mono"
+
         if Spectrumfile is not None and norm=="file":
             print("normalisation : spectrum file")
             # From spectrum file
@@ -458,7 +473,11 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
             print("normalisation : mono")
 
             nperchannel_norm = np.array([1.])
-
+
+        elif norm == "Gaussian" :
+            raise NotImplementedError("Gausssian norm for multiple bins not yet implemented")
+
+
         nperchannel = nperchannel_norm * nevents_sim
         # Full-sky?
         if not single_pixel:
@@ -487,9 +506,7 @@ def _open_rsp(cls, filename, Spectrumfile=None,norm="Linear" ,single_pixel = Fal
             pass
 
         # create a .h5 file with the good structure
-        filename = filename.replace(
-        ".rsp.gz", "_nside{0}.area.h5".format(nside))
-
+        filename = filename.replace(".rsp.gz","_nside{0}.area.h5".format(nside))
         f = h5.File(filename, mode='w')
 
         drm = f.create_group('DRM')