Merge branch '1867-csresmap-add-log-map-GSkyMap' into integration

Jürgen Knödlseder
2 parents a02887fb da87d300
Showing 4 changed files with 217 additions and 12 deletions
include/GSkyMap.hpp
pyext/GSkyMap.i
src/sky/GSkyMap.cpp
test/test_GSky.py
 /***************************************************************************
  *                       GSkyMap.hpp - Sky map class                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -85,6 +85,10 @@ class GVector;
 class GSkyMap : public GBase {
  
 	friend GSkyMap sqrt(const GSkyMap& map);
+	friend GSkyMap log(const GSkyMap& map);
+	friend GSkyMap log10(const GSkyMap& map);
+	friend GSkyMap abs(const GSkyMap& map);
+	friend GSkyMap sign(const GSkyMap& map);
  
 public:
     // Constructors and destructors
@@ -245,6 +245,18 @@ public:
     GSkyMap sqrt() {
         return sqrt(*self);
     }
+    GSkyMap log() {
+        return log(*self);
+    }
+    GSkyMap log10() {
+        return log10(*self);
+    }
+    GSkyMap abs() {
+        return abs(*self);
+    }
+    GSkyMap sign() {
+        return sign(*self);
+    }
     // Python 2.x operator/=
     GSkyMap __idiv__(const GSkyMap& map) {
         self->operator/=(map);
@@ -76,6 +76,9 @@
 #define G_READ_HEALPIX                   "GSkyMap::read_healpix(GFitsTable*)"
 #define G_READ_WCS                           "GSkyMap::read_wcs(GFitsImage*)"
 #define G_ALLOC_WCS                         "GSkyMap::alloc_wcs(GFitsImage*)"
+#define G_SQRT                                       "GSkyMap sqrt(GSkyMap&)"
+#define G_LOG                                         "GSkyMap log(GSkyMap&)"
+#define G_LOG10                                     "GSkyMap log10(GSkyMap&)"
  
 /* __ Macros _____________________________________________________________ */
  
@@ -3019,10 +3022,10 @@ bool GSkyMap::is_wcs(const GFitsHDU&amp; hdu) const
  ==========================================================================*/
  
 /***********************************************************************//**
- * @brief Computes square root of skymap elements
+ * @brief Computes square root of sky map elements
  *
- * @param[in] map Skymap.
- * @return Skymap containing the square root of every element.
+ * @param[in] map Sky map.
+ * @return Sky map containing the square root of every element.
  ***************************************************************************/
 GSkyMap sqrt(const GSkyMap& map)
 {
@@ -3038,20 +3041,177 @@ GSkyMap sqrt(const GSkyMap&amp; map)
             // Get the content from the bin
             double content = map(j,i);
  
-            // Check if content is not negative
+            // Throw an exception if content is negative
             if (content < 0.0) {
-                std::string msg = "Negative value envountered."
-                                  " Cannot take the sqrt from a negative value";
-                throw GException::invalid_value(G_OP_UNARY_DIV, msg);
+                std::string msg = "Negative value encountered. "
+                                  "Cannot take the sqrt from a negative "
+                                  "value";
+                throw GException::invalid_value(G_SQRT, msg);
             }
  
             // Set content of the result map
             result(j,i) = std::sqrt(content);
  
-        } // endfor: Loop over all bins
+        } // endfor: looped over all bins
+
+    } // endfor: looped over all maps
+
+    // Return sky map
+    return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes the natural logarithm of sky map elements
+ *
+ * @param[in] map Sky map.
+ * @return Sky map containing the natural logarithm of every element.
+ ***************************************************************************/
+GSkyMap log(const GSkyMap& map)
+{
+    // Initialise result vector
+    GSkyMap result(map);
+
+    // Loop over all maps
+    for (int i = 0; i < map.nmaps(); ++i) {
+
+        // Loop over all bins
+        for (int j = 0; j < map.npix(); ++j) {
+
+            // Get the content from the bin
+            double content = map(j,i);
+
+            // Throw an exception if content is negative
+            if (content <= 0.0) {
+                std::string msg = "Negative or zero value encountered. "
+                                  "Cannot calculate the logarithm of a "
+                                  "negative or zero value";
+                throw GException::invalid_value(G_LOG, msg);
+            }
+
+            // Set content of the result map
+            result(j,i) = std::log(content);
+
+        } // endfor: looped over all bins
+
+    } // endfor: looped over all maps
+
+    // Return sky map
+    return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes the base 10 logarithm of sky map elements
+ *
+ * @param[in] map Sky map.
+ * @return Sky map containing the base 10 logarithm of every element.
+ ***************************************************************************/
+GSkyMap log10(const GSkyMap& map)
+{
+    // Initialise result vector
+    GSkyMap result(map);
+
+    // Loop over all maps
+    for (int i = 0; i < map.nmaps(); ++i) {
+
+        // Loop over all bins
+        for (int j = 0; j < map.npix(); ++j) {
+
+            // Get the content from the bin
+            double content = map(j,i);
+
+            // Throw an exception if content is negative
+            if (content <= 0.0) {
+                std::string msg = "Negative or zero value encountered. "
+                                  "Cannot calculate the logarithm of a "
+                                  "negative or zero value";
+                throw GException::invalid_value(G_LOG10, msg);
+            }
+
+            // Set content of the result map
+            result(j,i) = std::log10(content);
+
+        } // endfor: looped over all bins
+
+    } // endfor: looped over all maps
+
+    // Return sky map
+    return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes the absolute value of sky map elements
+ *
+ * @param[in] map Sky map.
+ * @return Sky map containing the absolute value of every element.
+ ***************************************************************************/
+GSkyMap abs(const GSkyMap& map)
+{
+    // Initialise result vector
+    GSkyMap result(map);
+
+    // Loop over all maps
+    for (int i = 0; i < map.nmaps(); ++i) {
+
+        // Loop over all bins
+        for (int j = 0; j < map.npix(); ++j) {
+
+            // Get the content from the bin
+            double content = map(j,i);
+
+            // Set content of the result map
+            result(j,i) = std::abs(content);
+
+        } // endfor: looped over all bins
+
+    } // endfor: looped over all maps
+
+    // Return sky map
+    return result;
+}
+
+
+/***********************************************************************//**
+ * @brief Computes the sign value of sky map elements
+ *
+ * @param[in] map Sky map.
+ * @return Sky map containing the sign value of every pixel.
+ *
+ * This method returns a sky map filled with a value of 1 if the pixel
+ * is positive, a value of -1 if the pixel is negative or a value of 0
+ * if the pixel is 0.
+ ***************************************************************************/
+GSkyMap sign(const GSkyMap& map)
+{
+    // Initialise result vector
+    GSkyMap result(map);
+
+    // Loop over all maps
+    for (int i = 0; i < map.nmaps(); ++i) {
+
+        // Loop over all bins
+        for (int j = 0; j < map.npix(); ++j) {
+
+            // Get the content from the bin
+            double content = map(j,i);
+
+            // Handle the 3 cases
+            if (content < 0.0) {
+                result(j,i) = -1.0;
+            }
+            else if (content > 0.0) {
+                result(j,i) = +1.0;
+            }
+            else {
+                result(j,i) = 0.0;
+            }
+
+        } // endfor: looped over all bins
  
-    } // endfor: Loop over all maps
+    } // endfor: looed over all maps
  
-    // Return vector
+    // Return sky map
     return result;
 }
@@ -72,6 +72,7 @@ class Test(gammalib.GPythonTestSuite):
  
         # Append tests
         self.append(self.test_methods, "Test sky map methods")
+        self.append(self.test_friend_methods, "Test sky map friend methods")
         self.append(self.test_operators, "Test sky map operators")
         self.append(self.test_skymap_healpix, "Test HEALPix map")
         self.append(self.test_skymap_ait, "Test AIT projection map")
@@ -123,7 +124,35 @@ class Test(gammalib.GPythonTestSuite):
  
         # Return
         return
-    
+
+    # Test sky map friend methods
+    def test_friend_methods(self):
+        """
+        Test the sky map friend methods
+        """
+        # Setup sky maps
+        m = gammalib.GSkyMap('CAR', 'CEL', 83.63, 22.01, -3.7, 2.6, 10, 8, 12)
+
+        # Fill all values of skymaps with numbers > 0
+        values = range(1,m.npix()+1)
+        for j in range(m.nmaps()):
+            for i in range(m.npix()):
+                m[i,j] = values[i]
+
+        # Test map values after calling friend methods on one pixel
+        ipix = 18
+        self.test_value(m.sqrt()[ipix,0], values[ipix]**0.5, 'Check sqrt method')
+        self.test_value(m.log()[ipix,0], math.log(values[ipix]), 'Check log method')
+        self.test_value(m.log10()[ipix,0], math.log10(values[ipix]), 'Check log10 method')
+
+        # Set the interesting value to a negative value
+        m[ipix,0] = -values[ipix]
+        self.test_value(m.abs()[ipix,0], values[ipix], 'Check abs method')
+        self.test_value(m.sign()[ipix,0], -1, 'Check sign method')
+
+        # Return
+        return
+
     # Test sky map operators
     def test_operators(self):
         """