Added code to compress empty axes, i.e., those with axis shape = (1,)

Author: avalluvan

cosipy/response/DetectorResponse.py

@@ -42,6 +42,22 @@ def __init__(self, coord, **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']
@@ -165,7 +181,7 @@ def get_interp_response(self, target: dict):
         piecewise-linear directional responses.
         """
 
-        target = dict(sorted(target.items()))       # Sort dictionary by key (XXX: assuming response matrix also sorts in the same way)
+        # 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))
@@ -202,7 +218,7 @@ def get_spectral_response(self, copy = True):
 
         # Cache the spectral response
         if self._spec is None:
-            spec = self.project(['Ei','Em'])
+            spec = self.project(['Ei', self.mapping['Em']])
             self._spec = DetectorResponse(coord = self.coord,
                                           edges = spec.axes,
                                           contents = spec.contents,
@@ -310,7 +326,7 @@ def measured_energy_axis(self):
         :py:class:`histpy.Axes`        
         """
 
-        return self.axes['Em']
+        return self.axes[self.mapping['Em']]
 
 
 

docs/tutorials/response/LMDR.ipynb

@@ -1455,18 +1455,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "# response_path = Path('/Users/penguin/Documents/Grad School/Research/COSI/COSIpy/docs/tutorials/data/SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5')\n",
-    "response_path = Path('/Users/penguin/Documents/Grad School/Research/COSI/COSIpy/docs/tutorials/data/SMEXv12.44Ti.HEALPix04.E_1150_1164keV.binnedimaging.imagingresponse.nonsparse_nside16.area.h5')\n",
-    "# response_path = Path('/Users/penguin/Documents/Grad School/Research/COSI/COSIpy/docs/tutorials/data/DC2-Continuum_N-149_17-11-2024_14-32.imagingresponse_nside8.area.h5')"
+    "DATA_DIR = Path('/Users/penguin/Documents/Grad School/Research/COSI/COSIpy/docs/tutorials/data/')\n",
+    "\n",
+    "# response_path = DATA_DIR / Path('SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5')\n",
+    "response_path = DATA_DIR / Path('SMEXv12.44Ti.HEALPix04.E_1150_1164keV.binnedimaging.imagingresponse.nonsparse_nside16.area.h5')\n",
+    "# response_path = DATA_DIR / Path('DC2-Continuum_N-149_17-11-2024_14-32.imagingresponse_nside8.area.h5')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1482,11 +1484,20 @@
     "with FullDetectorResponse.open(response_path, unbinned=True) as response:\n",
     "    print(f\"Pixel 0 centered at {response.pix2skycoord(0)}\")\n",
     "    dr = response[0]\n",
-    "    data = response._file['DRM']['CONTENTS'][0, ...]\n",
-    "    # dr = DetectorResponse(edges=response.axes[1:], contents=data, unit=response.unit, interpolated_NuLambda=True)\n",
+    "    # data = response._file['DRM']['CONTENTS'][0, ...]\n",
+    "    # dr = DetectorResponse(edges=response.axes[1:], contents=data, unit=response.unit, interpolated_NuLambda=True)     # scipy.COO requires special handling\n",
     "    dr = response.get_interp_response(SkyCoord(lon=0, lat=0, frame=SpacecraftFrame(), unit=u.deg))"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dr = dr.drop_empty_axes()"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 7,
@@ -1533,7 +1544,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1561,12 +1572,23 @@
     }
    ],
    "source": [
-    "dr.get_spectral_response().get_interp_response({'Ei':Ei0, 'Em':Em0})"
+    "dr.get_spectral_response().get_interp_response({'Ei':511*u.keV, 'Em':511*u.keV})"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dr.mapping['PsiChi'] = None\n",
+    "dr.mapping['Theta'] = 'Theta'\n",
+    "dr.mapping['Zeta'] = 'Zeta'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -1586,7 +1608,16 @@
     }
    ],
    "source": [
-    "dr[3,3,2,385]"
+    "dr[58,4,20,108,100]         # dr[3,3,2,385]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dr.contents.shape"
    ]
   },
   {
@@ -1614,7 +1645,8 @@
     "Phi0 = 12*u.deg\n",
     "PsiChi0 = 386\n",
     "\n",
-    "target = {'Ei': Ei0, 'Em': Em0, 'Phi': Phi0, 'PsiChi': PsiChi0}\n",
+    "# target = {'Ei': Ei0, 'Em': Em0, 'Phi': Phi0, 'PsiChi': PsiChi0}\n",
+    "target = {'Ei': Ei0, 'Em': Em0, 'Phi': Phi0, 'Theta': 37*u.deg, 'Zeta': 201*u.deg}\n",
     "interpolated_response_value = dr.get_interp_response(target)\n",
     "\n",
     "interpolated_response_value"
@@ -1699,8 +1731,8 @@
     "target = {'Ei': Ei0, 'Em': Em0, 'Phi': Phi0, 'PsiChi': PsiChi0}\n",
     "\n",
     "fig, ax = plt.subplots()\n",
-    "dr.project(label1, label2).draw(ax=ax)\n",
-    "ax.scatter(target[label1], target[label2], marker='*')\n",
+    "dr.project(label1, dr.mapping[label2]).draw(ax=ax)\n",
+    "# ax.scatter(target[label1], target[label2], marker='*')\n",
     "for e1 in dr.neighbors[label1]:\n",
     "    for e2 in dr.neighbors[label2]:\n",
     "        ax.scatter(e1, e2, c='r')\n",