From dc400d5620c31401c560fd81cdaed79abad7e503 Mon Sep 17 00:00:00 2001 From: L. Tibaldo Date: Tue, 19 Dec 2017 16:32:20 +0100 Subject: [PATCH] ported residual computation to obsutils --- cscripts/csresmap.py | 66 +++++++----------------------------------------------------------- cscripts/obsutils.py | 89 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 59 deletions(-) diff --git a/cscripts/csresmap.py b/cscripts/csresmap.py index a5aab3b..bc3f397 100755 --- a/cscripts/csresmap.py +++ b/cscripts/csresmap.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(): diff --git a/cscripts/obsutils.py b/cscripts/obsutils.py index 10af62c..5e4eb95 100755 --- a/cscripts/obsutils.py +++ b/cscripts/obsutils.py @@ -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 -- libgit2 0.21.4