Replacing superfluous edits with files from cositools/develop

cosipy/response/DetectorResponse.py

@@ -1,12 +1,12 @@
 import logging
 logger = logging.getLogger(__name__)
 
-from pathlib import Path
-import itertools
+import numpy as np
+
 from copy import deepcopy
 
-from histpy import Histogram
-import numpy as np
+from histpy import Histogram, Axes, Axis
+
 import astropy.units as u
 
 class DetectorResponse(Histogram):
@@ -31,124 +31,12 @@ class DetectorResponse(Histogram):
         Physical area units, if not specified as part of ``contents``
     """
 
-    def __init__(self, interpolated_NuLambda=False, **kwargs):
+    def __init__(self, *args, **kwargs):
 
-        super().__init__(**kwargs)
-
-        self.interpolated_NuLambda = interpolated_NuLambda
-        self._set_mapping()
+        super().__init__(*args, **kwargs)
 
         self._spec = None
         self._aeff = None
-
-    def drop_empty_axes(self):
-
-        axis_names = self.axes.labels
-        axes_projection = []
-
-        for name in axis_names:
-            if (name in ['Ei', 'Em', 'eps', 'Eps']) | (self.axes[name].edges.size > 2):
-                axes_projection.append(name)
-
-        spec = self.project(axes_projection)
-
-        return DetectorResponse(edges = spec.axes,
-                                contents = spec.contents,
-                                unit = spec.unit,
-                                interpolated_NuLambda = self.interpolated_NuLambda)
-
-    def _set_mapping(self):
-        self.mapping = {}
-        target_names = ['Ei', 'Em', 'Phi', 'PsiChi', 'SigmaTau', 'Dist']
-        numlabels = len(self.axes.labels)
-
-        for key, label in zip(target_names[:numlabels], self.axes.labels):
-            self.mapping[key] = label           # Format: key_target : label
-                                                # Example: {'Ei': 'Ei', 'Em': 'eps', 'Phi': 'Phi', 'PsiChi': 'PsiChi'}
-
-    def _get_all_interp_weights(self, target: dict):
-
-        indices = []
-        weights = []
-
-        for key in target:
-            label = self.mapping[key]
-            axis = self.axes[label]
-            axis_scale = axis._scale
-            axis_type = str(type(axis)).split('.')[-1].strip("'>")      # XXX: Could probably be simplified using `isinstance()`
-
-            # Scale
-            if axis_scale in ['linear', 'log']:   # To ensure nonlinear binning parametrizations are converted to a linear scale.
-
-                # Axis Type
-                if axis_type in ['Axis', 'HealpixAxis']:
-                    if key == label:    # If key and label are the same, then there was no reparametrization along this axis
-                        idx, w = axis.interp_weights(target[key])
-                    else:
-                        centers = self.transform_eps_to_Em(axis.centers, target['Ei'])      # Transform coordinates to more physical units      # TODO: Generalize this
-                        absdiff = np.abs(centers - target[key])                             # Calculate absolute difference to given target
-                        idx = np.argpartition(absdiff, (1,2))[:2]                               # Find indices corresponding to two smallest absdiff
-                        w = 1 - np.partition(absdiff, (1,2))[:2] / (centers[1] - centers[0])    # Calculate weights corresponding to two smallest absdiff
-                else:
-                    raise ValueError(f'Axis type: {axis_type} is not supported')
-
-            # elif axis_scale == 'nonlinear':   
-            #     pass
-
-            else:
-                raise ValueError(f'{axis_scale} binning / scale scheme is not supported')
-
-            indices.append(idx)
-            weights.append(w)
-
-        return (indices, weights)
-
-    def transform_eps_to_Em(self, eps, Ei0):
-        return (eps + 1) * Ei0
-
-    def transform_Em_to_eps(self, Em, Ei0):
-        return Em/Ei0 - 1
-
-    def get_nearest_neighbors(self, target: dict, indices=None):
-        if indices is not None:
-            neighbors = {}
-            for idx, key in zip(indices, target):
-                label = self.mapping[key]
-                neighbors[label] = self.axes[label].centers[idx]
-            return neighbors
-
-        else:
-            target = dict(sorted(target.items()))
-            indices, _ = self._get_all_interp_weights(target)
-            return self.get_nearest_neighbors(target, indices)
-
-    def get_interp_response(self, target: dict):
-        """
-        Currently only supports nonlinear spectral responses (
-        and for a particular parametrization)
-        TODO: In the future, this will also support nonlinear / 
-        piecewise-linear directional responses.
-        """
-
-        # target = dict(sorted(target.items()))       # Sort dictionary by key (XXX: assuming response matrix also sorts in the same way)
-        indices, weights = self._get_all_interp_weights(target)
-        perm_indices = list(itertools.product(*indices))
-        perm_weights = list(itertools.product(*weights))
-
-        if len(target) == len(self.axes):
-            interpolated_response_value = 0
-            for idx, w in zip(perm_indices, perm_weights):
-                interpolated_response_value += np.prod(w) * self.contents[idx]
-
-        else:
-            interpolated_response_value = np.zeros(len(self.axes['Ei']) - 1) * self.contents.unit     # XXX: Assuming all measured variables require interpolation
-            for idx, w in zip(perm_indices, perm_weights):
-                i = (Ellipsis,) + idx   # XXX: Assuming 'Ei' is the first index
-                interpolated_response_value += np.prod(w) * self.contents[i]
-
-        self.neighbors = self.get_nearest_neighbors(target, indices)
-
-        return interpolated_response_value
 
     def get_spectral_response(self, copy = True):
         """
@@ -167,8 +55,8 @@ def get_spectral_response(self, copy = True):
 
         # Cache the spectral response
         if self._spec is None:
-            spec = self.project(['Ei', self.mapping['Em']])
-            self._spec = DetectorResponse(edges=spec.axes,
+            spec = self.project(['Ei','Em'])
+            self._spec = DetectorResponse(spec.axes,
                                           contents = spec.contents,
                                           unit = spec.unit)
 
@@ -235,21 +123,6 @@ def get_dispersion_matrix(self):
 
         # Normalize column-by-column
         return (spec / norm)
-
-        # NOTE: When histpy is updated to do away without overflow and
-        # underflow bins, the following lines of code can replace this
-        # function body.
-
-        # # Get spectral response and effective area normalization
-        # spec = self.get_spectral_response(copy = False)
-        # norm = self.get_effective_area().contents
-
-        # # "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)    # XXX: Runtime warning needs to be dealt with in histpy. 
-        #                         # Overflow and underflow bins functionality should be removed.
 
     @property
     def photon_energy_axis(self):
@@ -274,7 +147,7 @@ def measured_energy_axis(self):
         :py:class:`histpy.Axes`        
         """
 
-        return self.axes[self.mapping['Em']]
+        return self.axes['Em']
 
 
 

cosipy/response/FullDetectorResponse.py

@@ -1,6 +1,5 @@
 from .PointSourceResponse import PointSourceResponse
 from .DetectorResponse import DetectorResponse
-from .ListModePSR import ListModePSR
 from astromodels.core.model_parser import ModelParser
 import matplotlib.pyplot as plt
 from astropy.time import Time
@@ -50,8 +49,7 @@ def __init__(self, *args, **kwargs):
         pass
 
     @classmethod
-    def open(cls, filename, unbinned=False, 
-             Spectrumfile=None,norm="Linear" ,single_pixel = False,alpha=0,emin=90,emax=10000):
+    def open(cls, filename,Spectrumfile=None,norm="Linear" ,single_pixel = False,alpha=0,emin=90,emax=10000):
         """
         Open a detector response file.
 
@@ -80,7 +78,7 @@ def open(cls, filename, unbinned=False,
         """
 
         filename = Path(filename)
-        cls.unbinned = unbinned
+
 
         if filename.suffix == ".h5":
             return cls._open_h5(filename)
@@ -757,15 +755,23 @@ def __getitem__(self, pix):
         if not isinstance(pix, (int, np.integer)) or pix < 0 or not pix < self.npix:
             raise IndexError("Pixel number out of range, or not an integer")
 
+        #check if we have sparse matrice or not
+
         if self._sparse:
-            coords = np.reshape(self._file['DRM']['BIN_NUMBERS'][pix], (self.ndim-1, -1))
+            coords = np.reshape(
+                self._file['DRM']['BIN_NUMBERS'][pix], (self.ndim-1, -1))
             data = np.array(self._file['DRM']['CONTENTS'][pix])
-            contents = COO(coords=coords, data=data, shape=tuple(self.axes.nbins[1:]))
-        else:
-            data = self._file['DRM']['CONTENTS'][pix]
-            contents = data
 
-        return DetectorResponse(edges=self.axes[1:], contents=contents, unit=self.unit, interpolated_NuLambda=self.unbinned)
+            return DetectorResponse(self.axes[1:],
+                                contents=COO(coords=coords,
+                                             data=data,
+                                             shape=tuple(self.axes.nbins[1:])),
+                                unit=self.unit)
+
+        else :
+            data = self._file['DRM']['CONTENTS'][pix]
+            return DetectorResponse(self.axes[1:],
+                                contents=data, unit=self.unit)                 
 
     def close(self):
         """
@@ -816,10 +822,12 @@ def get_interp_response(self, coord):
 
         pixels, weights = self.get_interp_weights(coord)
 
-        dr = DetectorResponse(edges=self.axes[1:],
-                            sparse=self._sparse,
-                            unit=self.unit,
-                            interpolated_NuLambda=self.unbinned)
+
+
+        dr = DetectorResponse(self.axes[1:],
+                              sparse=self._sparse,
+                              unit=self.unit)
+
 
         for p, w in zip(pixels, weights):
 
@@ -836,7 +844,7 @@ def get_point_source_response(self,
         Convolve the all-sky detector response with exposure for a source at a given
         sky location.
 
-        Provide either a exposure map (aka dwell time map) or a combination of a 
+        Provide either a exposure map (aka dweel time map) or a combination of a 
         sky coordinate and a spacecraft attitude map.
 
         Parameters
@@ -871,8 +879,8 @@ def get_point_source_response(self,
                     "Exposure map has a different grid than the detector response")
 
             psr = PointSourceResponse(self.axes[1:],
-                                    sparse=self._sparse,
-                                    unit=u.cm*u.cm*u.s)
+                                      sparse=self._sparse,
+                                      unit=u.cm*u.cm*u.s)
 
             for p in range(self.npix):
 
@@ -894,60 +902,42 @@ def get_point_source_response(self,
             if self.is_sparse:
                 raise ValueError("Coord +  scatt_map currently only supported for dense responses")
 
-            axis_label = "PsiChi"
+            axis = "PsiChi"
 
-            coords_axis = Axis(np.arange(coord.size+1), label = 'coords')   # Create axis of length number of input coords + 1
+            coords_axis = Axis(np.arange(coord.size+1), label = 'coords')
 
-            psr = Histogram([coords_axis] + list(deepcopy(self.axes[1:])),  # Create new "NuLambda" axis
+            psr = Histogram([coords_axis] + list(deepcopy(self.axes[1:])), 
                             unit = self.unit * scatt_map.unit)
 
-            psr.axes[axis_label].coordsys = coord.frame     # Set coordinate system of PsiChi axis to input coordinate frame. Axis coordsys was set when response file was opened and initialized using HealpixBase.__init__
+            psr.axes[axis].coordsys = coord.frame
 
-            for pixels, exposure in \
-                            zip(scatt_map.contents.coords.transpose(),        # Arrays of [[x's], [y's]] --> [[x,y]'s]
-                            scatt_map.contents.data):
+            for i,(pixels, exposure) in \
+                enumerate(zip(scatt_map.contents.coords.transpose(),
+                              scatt_map.contents.data)):
 
                 #gc.collect() # HDF5 cache issues
 
-                # Calculate attitude of coordinates in scatt_map
                 att = Attitude.from_axes(x = scatt_map.axes['x'].pix2skycoord(pixels[0]),
-                                        y = scatt_map.axes['y'].pix2skycoord(pixels[1]))
-
-                coord.attitude = att    # Set attitude to given input coordinate. It is used to transform Spacecraft coordinates (scoords) from ICRS to mhealpy pixel
-
-                if self.unbinned:
-                    dr_list = []
-
-                    if coord.size == 1:
-                        dr = self.get_interp_response(coord)
-                        dr_list.append(dr)
+                                         y = scatt_map.axes['y'].pix2skycoord(pixels[1]))
 
-                    else:
-                        for c in coord:
-                            dr = self.get_interp_response(c)
-                            dr_list.append(dr)                    
-
-                    dr_pix_2 = Histogram.concatenate(coords_axis, dr_list)
-                    dr_pix = deepcopy(dr_pix_2)     # See XXX note below
+                coord.attitude = att
 
-                else:
-                    loc_nulambda_pixels = np.array(self.axes['NuLambda'].find_bin(coord), ndmin = 1)  # Pixel corresponding to coord
-                    dr_pix = Histogram.concatenate(coords_axis, [self[pixel] for pixel in loc_nulambda_pixels])     # Concatenate Axis object with DetectorResponse Histogram (Histograms if multiple target_coords are provided)
+                #TODO: Change this to interpolation
+                loc_nulambda_pixels = np.array(self.axes['NuLambda'].find_bin(coord),
+                                               ndmin = 1)
+
+                dr_pix = Histogram.concatenate(coords_axis, [self[i] for i in loc_nulambda_pixels])
 
-                dr_pix.axes['PsiChi'].coordsys = SpacecraftFrame(attitude = att)    # Update the attitude of the PsiChi axis from `None` to `att`
-                                                                                    # XXX: this line is messing up the coordsys of response. I was 
-                                                                                    # forced to include a `deepcopy` statement to strictly divide 
-                                                                                    # their relationship. Alternate implementations are welcome. 
-                                                                                    # Even `copy` is not sufficient.
+                dr_pix.axes['PsiChi'].coordsys = SpacecraftFrame(attitude = att)
 
-                self._sum_rot_hist(dr_pix, psr, exposure)   # Rotate PsiChi from local SC to galactic coordinates. (Function only affects Healpix Axis of psr Histogram)
+                self._sum_rot_hist(dr_pix, psr, exposure)
 
             # Convert to PSR
             psr = tuple([PointSourceResponse(psr.axes[1:],
-                                            contents = data,
-                                            sparse = psr.is_sparse,
-                                            unit = psr.unit)
-                        for data in psr[:]])
+                                             contents = data,
+                                             sparse = psr.is_sparse,
+                                             unit = psr.unit)
+                         for data in psr[:]])
 
             if coord.size == 1:
                 return psr[0]

cosipy/response/PointSourceResponse.py

@@ -61,7 +61,6 @@ def get_expectation(self, spectrum):
 
         energy_axis = self.photon_energy_axis
 
-
         flux = get_integrated_spectral_model(spectrum, energy_axis)
 
         expectation = np.tensordot(self.contents, flux.contents, axes = ([0], [0]))

cosipy/spacecraftfile/SpacecraftFile.py

@@ -134,7 +134,6 @@ def parse_from_file(cls, file):
         orientation_file = np.loadtxt(file, usecols=(1, 2, 3, 4, 5, 6, 7, 8),delimiter=' ', skiprows=1, comments=("#", "EN"))
         time_stamps = orientation_file[:, 0]
         axis_1 = orientation_file[:, [2, 1]]
-
         axis_2 = orientation_file[:, [4, 3]]
         axis_3 = orientation_file[:, [7, 6]]
         altitude = np.array(orientation_file[:, 5]) 
@@ -1115,7 +1114,3 @@ def plot_rmf(self, file_name = None, save_name = None, dpi = 300):
         #plt.show()
 
         return
-
-    # len() functionality for SpacecraftFile / ori objects
-    def __len__(self):
-        return len(self._time)