Merge pull request #364 from McKelvey-Engineering-CSE/new_response_format

Improved Full Detector Response storage and conversion

Author: israelmcmc

CHANGELOG.md

@@ -1,12 +1,59 @@
 # Changelog
 
-## Version 0.3.0
+## Version 0.4.0
 
 Under development. 
 
 Developers: please keep track of notable changes here.
 
-- 
+### New FullDetectorResponse 
+
+HDF5 response file format has changed substantially.  The new format,
+which is *not* compatible with the .h5 response format from prior
+releases, features
+  - internal compression using BitShuffle algorithm, which gives chunk
+    access times within 50% of uncompressed for smaller responses and
+	much faster for larger ones while offering same or better compression 
+	on disk.
+  - much faster conversion from .rsp.gz to .h5, with greatly reduced
+    memory usage
+  - ability to convert from .rsp.gz to .h5 and back to 
+    .rsp.gz without losing any of the original header information
+	(which facilitates creating lower-resolution versions of large
+	responses)
+  - separation of raw counts from effective area, which allows
+    FullDetectorResponse to export the latter as needed and
+    enables future performance improvements for code that
+	needs to compute and average PSRs.  The dtype of the effective 
+    area, and hence the values returned when computing PSRs etc.,
+	can now be set at load time; hence, the same response can be
+	used in float32 or float64 with no additional overhead.
+  - automatic "good chunks" chunking of .h5 responses at creation time
+  - HDF5 tweaks to avoid writing timestamps, so that (with support
+    from a forthcoming release of histpy) the MD5 signature of an
+	.h5 response does not change every time it is written.
+	
+.rsp.gz conversion has been broken out from FullDetectorResponse into
+its own class, RspConverter, which simplifies the FullDetectorResponse
+code and in particular its open() interface.
+
+The new response code deprecates and removes the following
+response-related functionality:
+  - sparse response format
+  - reading spectrum for effective area computation from a file
+     (was broken/bit-rotted in DC3 release)
+  - miniDC2 format support (unused since DC2)
+
+### New wasabi location for development files
+
+To support file format changes in the develop branch that are
+incompatible with DC3 and cosipy versions <= 0.3, we have established
+a "COSI-SMEX/develop" tree in the public wasabi bucket.  Currently,
+this tree holds the new-format .h5 detector responses and the
+corresponding .rsp.gz files.  The DC3 tutorials have been updated to
+use this new tree for files whose format has changed since DC3.
+
+
 
 ## Version 0.2.x
 

cosipy/image_deconvolution/dataIF_COSI_DC2.py

@@ -36,13 +36,10 @@ def __init__(self, name = None):
         self._image_response = None # histpy.Histogram (dense)
 
         # None if using Galactic CDS, mandotary if using local CDS
-        self._coordsys_conv_matrix = None 
-
-        # optional
-        self.is_miniDC2_format = False #should be removed in the future
+        self._coordsys_conv_matrix = None
 
     @classmethod
-    def load(cls, name, event_binned_data, dict_bkg_binned_data, rsp, coordsys_conv_matrix = None, is_miniDC2_format = False):
+    def load(cls, name, event_binned_data, dict_bkg_binned_data, rsp, coordsys_conv_matrix = None):
         """
         Load data
 
@@ -57,9 +54,7 @@ def load(cls, name, event_binned_data, dict_bkg_binned_data, rsp, coordsys_conv_
         rsp : :py:class:`histpy.Histogram` or :py:class:`cosipy.response.FullDetectorResponse`
             Response
         coordsys_conv_matrix : :py:class:`cosipy.image_deconvolution.CoordsysConversionMatrix`, default False
-            Coordsys conversion matrix 
-        is_miniDC2_format : bool, default False
-            Whether the file format is for mini-DC2. It will be removed in the future.
+            Coordsys conversion matrix
 
         Returns
         -------
@@ -87,24 +82,22 @@ def load(cls, name, event_binned_data, dict_bkg_binned_data, rsp, coordsys_conv_
         if new._coordsys_conv_matrix is not None and \
            new._coordsys_conv_matrix.is_sparse:
             new._coordsys_conv_matrix.contents.enable_caching()
-            
-        new.is_miniDC2_format = is_miniDC2_format
 
         if isinstance(rsp, FullDetectorResponse):
             logger.info('Loading the response matrix onto your computer memory...')
-            new._load_full_detector_response_on_memory(rsp, is_miniDC2_format)
+            new._image_response = rsp.to_dr()
             logger.info('Finished')
         elif isinstance(rsp, Histogram):
             new._image_response = rsp
-        
+
         # We modify the axes in event, bkg_models, response. This is only for DC2.
         new._modify_axes()
-        
+
         new._data_axes = new._event.axes
-        
+
         if new._coordsys_conv_matrix is None:
             axes = (new._image_response.axes['NuLambda'].copy(), new._image_response.axes['Ei']) # will mutate axes[0]
-            axes[0].label = 'lb' 
+            axes[0].label = 'lb'
             # The gamma-ray direction of pre-computed response in DC2 is in the galactic coordinate, not in the local coordinate.
             # Actually, it is labeled as 'NuLambda'. So I replace it with 'lb'.
             new._model_axes = Axes(axes, copy_axes = False)
@@ -134,39 +127,35 @@ def _modify_axes(self):
         for name in axis_name:
 
             logger.info(f"... checking the axis {name} of the event and background files...")
-            
+
             event_edges, event_unit = self._event.axes[name].edges, self._event.axes[name].unit
 
             for key in self._bkg_models:
 
                 bkg_edges, bkg_unit = self._bkg_models[key].axes[name].edges, self._bkg_models[key].axes[name].unit
 
                 if np.all(event_edges == bkg_edges):
-                    logger.info(f"    --> pass (edges)") 
+                    logger.info(f"    --> pass (edges)")
                 else:
                     logger.error(f"Warning: the edges of the axis {name} are not consistent between the event and the background model {key}!")
                     logger.error(f"         event      : {event_edges}")
                     logger.error(f"         background : {bkg_edges}")
                     raise ValueError
 
-        # check the axes of the event/response files. 
+        # check the axes of the event/response files.
         # Note that currently (2023-08-29) no unit is stored in the binned data. So only the edges are compared. This should be modified in the future.
 
         axis_name = ['Em', 'Phi', 'PsiChi']
-        
+
         for name in axis_name:
 
             logger.info(f"...checking the axis {name} of the event and response files...")
 
             event_edges, event_unit = self._event.axes[name].edges, self._event.axes[name].unit
             response_edges, response_unit = self._image_response.axes[name].edges, self._image_response.axes[name].unit
-            
-#            if type(response_edges) == u.quantity.Quantity and self.is_miniDC2_format == True:
-            if event_unit is None and response_unit is not None and self.is_miniDC2_format == True: # this is only for the old data in the miniDC2 challenge. I will remove them in the near future (or in the final dataIF).
-                response_edges = response_edges.value
 
             if np.all(event_edges == response_edges):
-                logger.info(f"    --> pass (edges)") 
+                logger.info(f"    --> pass (edges)")
             else:
                 logger.error(f"Warning: the edges of the axis {name} are not consistent between the event and background!")
                 logger.error(f"        event      : {event_edges}")
@@ -184,12 +173,12 @@ def _modify_axes(self):
                              self._image_response.axes["Phi"], \
                              self._image_response.axes["PsiChi"]),
                             copy_axes=False)
-        
+
         self._event = Histogram(axes_cds,
                                 contents = self._event.contents,
                                 unit = self._event.unit,
                                 copy_contents = False) # overwrite axes of existing Histogram
-                                
+
         for key in self._bkg_models:
             bkg_model = self._bkg_models[key]
             self._bkg_models[key] = Histogram(axes_cds,
@@ -201,36 +190,13 @@ def _modify_axes(self):
 
         return True
 
-    def _load_full_detector_response_on_memory(self, full_detector_response, is_miniDC2_format):
-        """
-        Load a response file on the computer memory.
-        """
-
-        axes_image_response = Axes((full_detector_response.axes["NuLambda"],
-                                    full_detector_response.axes["Ei"],
-                                    full_detector_response.axes["Em"],
-                                    full_detector_response.axes["Phi"],
-                                    full_detector_response.axes["PsiChi"]),
-                                   copy_axes=False)
-        
-        if is_miniDC2_format:
-            npix = axes_image_response["NuLambda"].npix
-            slices = [ np.sum(full_detector_response[ipix].to_dense(), axis = (4,5)) for ipix in tqdm(range(npix)) ] #Ei, Em, Phi, ChiPsi
-            contents = np.stack(slices)
-        else:
-            contents = full_detector_response._file['DRM']['CONTENTS']
-
-        self._image_response = Histogram(axes_image_response, contents=contents,
-                                         unit = full_detector_response.unit,
-                                         copy_contents = False)
-
     def _calc_exposure_map(self):
         """
         Calculate exposure_map, which is an intermediate matrix used in RL algorithm.
         """
 
         logger.info("Calculating an exposure map...")
-        
+
         if self._coordsys_conv_matrix is None:
             exposure_map = np.sum(self._image_response.contents, axis = (2,3,4))
         else:
@@ -275,7 +241,7 @@ def calc_expectation(self, model, dict_bkg_norm = None, almost_zero = 1e-12):
         # This is just because in DC2 the rotate response for galactic coordinate CDS does not have an axis for time/scatt binning.
         # However it is likely that it will have such an axis in the future in order to consider background variability depending on time and pointing direction etc.
         # Then, the implementation here will not work. Thus, keep in mind that we need to modify it once the response format is fixed.
-        
+
         if self._coordsys_conv_matrix is None:
             expectation = np.tensordot( model.contents, self._image_response.contents, axes = ((0,1),(0,1)))
             # ['lb', 'Ei'] x [NuLambda(lb), Ei, Em, Phi, PsiChi] -> [Em, Phi, PsiChi]
@@ -293,12 +259,12 @@ def calc_expectation(self, model, dict_bkg_norm = None, almost_zero = 1e-12):
         expectation *= model.axes['lb'].pixarea()
         expectation += almost_zero
 
-        if dict_bkg_norm is not None: 
+        if dict_bkg_norm is not None:
             for key in self.keys_bkg_models():
                 expectation += self.bkg_model(key).contents * dict_bkg_norm[key]
 
         return Histogram(self.data_axes, contents = expectation, copy_contents = False)
-        
+
     def calc_T_product(self, dataspace_histogram):
         """
         Calculate the product of the input histogram with the transonse matrix of the response function.
@@ -375,4 +341,4 @@ def calc_log_likelihood(self, expectation):
         """
         log_likelihood = np.sum( self.event * np.log(expectation) ) - np.sum(expectation)
 
-        return log_likelihood 
+        return log_likelihood

cosipy/response/DetectorResponse.py

@@ -69,7 +69,7 @@ def get_spectral_response(self, copy = True):
         ----------
         copy : bool
             If true, a copy of the cached spectral response will be returned.
-        
+
         Returns
         -------
         :py:class:`DetectorResponse`
@@ -80,13 +80,14 @@ def get_spectral_response(self, copy = True):
             spec = self.project(['Ei','Em'])
             self._spec = DetectorResponse(spec.axes,
                                           contents = spec.contents,
-                                          unit = spec.unit)
+                                          unit = spec.unit,
+                                          copy_contents = False)
 
         if copy:
             return self._spec.copy()
         else:
             return self._spec
-        
+
     def get_effective_area(self, energy = None, copy = True):
         """
         Compute the effective area at a given energy. If no energy is specified, the
@@ -98,14 +99,14 @@ def get_effective_area(self, energy = None, copy = True):
             Energy/energies at which to interpolate the linearly effective area
         copy : bool
             If true, a copy of the cached effective will be returned.
-        
+
         Returns
         -------
         :py:class:`astropy.units.Quantity` or :py:class:`histpy.Histogram`
         """
-        
+
         if self._aeff is None:
-            self._aeff = self.get_spectral_response(copy = False).project('Ei').to_dense()
+            self._aeff = self.get_spectral_response(copy = False).project('Ei')
 
         if energy is None:
             if copy:
@@ -125,7 +126,7 @@ def get_dispersion_matrix(self):
         -------
         :py:class:`histpy.Histogram`
         """
-        
+
         # Get spectral response and effective area normalization
         spec = self.get_spectral_response(copy = False)
         norm = self.get_effective_area().full_contents
@@ -139,10 +140,10 @@ def get_dispersion_matrix(self):
             norm[norm == 0] = 1*norm.unit
 
             logger.warn("Null effective area, cannot properly compute dispersion matrix.")
-        
+
         # "Broadcast" such that it has the compatible dimensions with the 2D matrix
         norm = spec.expand_dims(norm, 'Ei')
-        
+
         # Normalize column-by-column
         return (spec / norm)
 
@@ -155,22 +156,18 @@ def photon_energy_axis(self):
         -------
         :py:class:`histpy.Axes`
         """
-        
+
         return self.axes['Ei']
 
-    
+
     @property
     def measured_energy_axis(self):
         """
         Measured energy bins (``Em``).
 
         Returns
         -------
-        :py:class:`histpy.Axes`        
+        :py:class:`histpy.Axes`
         """
-        
-        return self.axes['Em']
-        
 
-        
-    
+        return self.axes['Em']