ported residual computation to obsutils

L. Tibaldo
1 parent e55c60c7
Showing 2 changed files with 96 additions and 59 deletions
cscripts/csresmap.py
cscripts/obsutils.py
@@ -19,9 +19,9 @@
 #
 # ==========================================================================
 import sys
-import math
 import gammalib
 import ctools
+from cscripts import obsutils
  
  
 # ============== #
@@ -209,67 +209,15 @@ class csresmap(ctools.csobservation):
             # Get model map into GSkyMap object
             modelmap = model.cube().counts().copy()
  
-        # Store counts map as residual map. Note that we need a special
+        # Calculate residual maps
+        # Note that we need a special
         # construct here to avoid memory leaks. This seems to be a SWIG feature
         # as SWIG creates a new object when calling binning.cube()
-        self._resmap = countmap.copy()
-        self._resmap.stack_maps()
+        countmap1 = countmap.copy()
+        countmap1.stack_maps()
+        print("TEST",countmap1.classname())
         modelmap.stack_maps()
-
-        # Get residual map algorithm type
-        algorithm = self['algorithm'].string()
-
-        # Subtract maps
-        if algorithm == 'SUB':
-            self._resmap -= modelmap
-
-        # Subtract and divide by model map
-        elif algorithm == 'SUBDIV':
-            self._resmap -= modelmap
-            self._resmap /= modelmap
-
-        # Subtract and divide by sqrt of model map
-        elif algorithm == 'SUBDIVSQRT':
-            self._resmap -= modelmap
-            self._resmap /= modelmap.sqrt()
-
-        # Calculate significance from Li&Ma derivation
-        elif algorithm == 'SIGNIFICANCE':
-
-            # Compute sign map
-            signmap = (self._resmap - modelmap).sign()
-
-            # Loop over every bin in the map
-            for i in range(self._resmap.npix()):
-
-                # If the model value > 0.0 do the computation as normal ...
-                model_val = modelmap[i]
-                if model_val > 0.0:
-
-                    # If the data value is also > 0 then compute the
-                    # significance^2 and save it to the map ...
-                    data_val = self._resmap[i]
-                    if data_val > 0.0:
-                        log_val         = math.log(data_val/model_val)
-                        self._resmap[i] = (data_val*log_val) + model_val  - data_val
-
-                    # ... otherwise compute the reduced value of the above
-                    # expression. This is necessary to avoid computing log(0).
-                    else:
-                        self._resmap[i] = model_val
-
-                # ... otherwise hard-code the significance to 0
-                else:
-                    self._resmap[i] = 0.0
-
-            # Compute significance map
-            self._resmap *= 2.0
-            self._resmap  = self._resmap.sqrt()
-            self._resmap  = self._resmap * signmap
-
-        # Raise exception if algorithm is unknown
-        else:
-            raise TypeError('Algorithm "'+algorithm+'" not known')
+        self._resmap = obsutils.residuals(self,countmap1,modelmap)
  
         # Optionally publish map
         if self['publish'].boolean():
@@ -828,3 +828,92 @@ def get_onoff_obs(cls, obs):
  
     # Return On/Off oberservation container
     return onoff_obs
+
+# ========================================= #
+# Calculate residuals from counts and model #
+# ========================================= #
+def residuals(cls,counts,model):
+    """
+    Calculate residuals given counts and models,
+    according to algorithm specified by user.
+    Can handle GSkyMap or GNdarray objects
+    :param cls: `~ctools.cscript`
+        cscript class
+    :param counts: `~gammalib.GSkyMap/~gammalib.GNdarray' data counts
+    :param model: `~gammalib.GSkyMap/~gammalib.GNdarray' model counts
+    :return residuals: `~gammalib.GSkyMap/~gammalib.GNdarray' residuals
+    """
+    # Find type of objects we are manipulating
+    # and set size to iterate later
+
+    # If GNdarray
+    if counts.classname() == 'GNdarray':
+        nelem = counts.size()
+    elif counts.classname() == 'GSkyMap':
+        nelem = counts.npix()
+    else:
+        msg = 'obsutils.residuals only handles ' \
+              + 'GNdarray or GSkyMap objects.\n'
+        raise RuntimeError(msg)
+
+    # Copy counts to initialise residuals
+    residuals = counts.copy()
+
+    # Get residual map algorithm type
+    algorithm = cls['algorithm'].string()
+
+    # Subtract
+    if algorithm == 'SUB':
+        residuals -= model
+
+    # Subtract and divide by model
+    elif algorithm == 'SUBDIV':
+        residuals -= model
+        residuals /= model
+
+    # Subtract and divide by sqrt of model
+    elif algorithm == 'SUBDIVSQRT':
+        residuals -= model
+        residuals /= model.sqrt()
+
+    # Calculate significance from Li&Ma
+    elif algorithm == 'SIGNIFICANCE':
+
+        # Compute sign
+        sign = (residuals - model).sign()
+
+        # Loop over every bin
+        for i in range(nelem):
+
+            # If the model value > 0.0 do the computation as normal ...
+            model_val = model[i]
+            if model_val > 0.0:
+
+                # If the data value is also > 0 then compute the
+                # significance^2 and save it ...
+                data_val = residuals[i]
+                if data_val > 0.0:
+                    log_val = math.log(data_val / model_val)
+                    residuals[i] = (
+                                      data_val * log_val) + model_val - data_val
+
+                # ... otherwise compute the reduced value of the above
+                # expression. This is necessary to avoid computing log(0).
+                else:
+                    residuals[i] = model_val
+
+            # ... otherwise hard-code the significance to 0
+            else:
+                residuals[i] = 0.0
+
+        # Compute significance map
+        residuals *= 2.0
+        residuals = residuals.sqrt()
+        residuals = residuals * sign
+
+    # Raise exception if algorithm is unknown
+    else:
+        raise TypeError('Algorithm "' + algorithm + '" not known')
+
+    # Return
+    return residuals
 \ No newline at end of file