Implement operators and add test code to illustrate smoothing of an image (#1731)

Jürgen Knödlseder
1 parent 3f24a1ab
Showing 4 changed files with 185 additions and 22 deletions
include/GFft.hpp
pyext/GFft.i
src/numerics/GFft.cpp
test/test_GNumerics.py
@@ -96,6 +96,8 @@ protected:
     void copy_members(const GFft& fft);
     void free_members(void);
     void set_data(const GNdarray& array);
+    bool has_same_shape(const GFft& fft) const;
+    void require_same_shape(const std::string& method, const GFft& fft) const;
  
     // Low-level FFT methods
     void transform(std::complex<double>* data,
@@ -114,6 +114,33 @@ public:
         }
         (*self)(index) = value;
     }
+    GFft __add__(const GFft &a) {
+        return (*self) + a;
+    }
+    GFft __sub__(const GFft &a) {
+        return (*self) - a;
+    }
+    GFft __mul__(const GFft &a) {
+        return (*self) * a;
+    }
+    // Python 2.x
+    GFft __div__(const GFft &a) {
+        return (*self) / a;
+    }
+    // Python 3.x
+    GFft __truediv__(const GFft &a) {
+        return (*self) / a;
+    }
+    // Python 2.x operator/=
+    GFft __idiv__(const GFft& a) {
+        self->operator/=(a);
+        return (*self);
+    }
+    // Python 3.x operator/=
+    GFft __itruediv__(const GFft& a) {
+        self->operator/=(a);
+        return (*self);
+    }
     GFft copy() {
         return (*self);
     }
@@ -35,6 +35,10 @@
 #include "GNdarray.hpp"
  
 /* __ Method name definitions ____________________________________________ */
+#define G_OP_ADD                                    "GFft::operator+=(GFft&)"
+#define G_OP_SUB                                    "GFft::operator-=(GFft&)"
+#define G_OP_MUL                                    "GFft::operator*=(GFft&)"
+#define G_OP_DIV                                    "GFft::operator/=(GFft&)"
 #define G_FORWARD                                  "GFft::forward(GNdarray&)"
 #define G_BACKWARD                                         "GFft::backward()"
  
@@ -148,13 +152,18 @@ GFft&amp; GFft::operator=(const GFft&amp; fft)
  * @param[in] fft Fast Fourier Transform.
  * @return Fast Fourier Transform.
  *
- * Adds a Fast Fourier Transform to another Fast Fourier Transform.
- *
- * @todo Implement operator
+ * Adds a Fast Fourier Transform to another Fast Fourier Transform. The
+ * method requires that both FFTs have the same dimension.
  ***************************************************************************/
 GFft& GFft::operator+=(const GFft& fft)
 {
-    // TODO: Implement operator
+    // Throw an exception if the FFTs have not the same shape
+    require_same_shape(G_OP_ADD, fft);
+    
+    // Add elements
+    for (int i = 0; i < m_size; ++i) {
+        *(m_data+i) += *(fft.m_data+i);
+    }
  
     // Return Fast Fourier Transform
     return *this;
@@ -168,12 +177,17 @@ GFft&amp; GFft::operator+=(const GFft&amp; fft)
  * @return Fast Fourier transform.
  *
  * Subtracts a Fast Fourier Transform from another Fast Fourier Transform.
- *
- * @todo Implement operator
+ * The method requires that both FFTs have the same dimension.
  ***************************************************************************/
 GFft& GFft::operator-=(const GFft& fft)
 {
-    // TODO: Implement operator
+    // Throw an exception if the FFTs have not the same shape
+    require_same_shape(G_OP_SUB, fft);
+
+    // Subtract elements
+    for (int i = 0; i < m_size; ++i) {
+        *(m_data+i) -= *(fft.m_data+i);
+    }
  
     // Return Fast Fourier Transform
     return *this;
@@ -187,12 +201,17 @@ GFft&amp; GFft::operator-=(const GFft&amp; fft)
  * @return Fast Fourier transform.
  *
  * Multiplies a Fast Fourier Transform to another Fast Fourier Transform.
- *
- * @todo Implement operator
+ * The method requires that both FFTs have the same dimension.
  ***************************************************************************/
 GFft& GFft::operator*=(const GFft& fft)
 {
-    // TODO: Implement operator
+    // Throw an exception if the FFTs have not the same shape
+    require_same_shape(G_OP_MUL, fft);
+
+    // Multiply elements
+    for (int i = 0; i < m_size; ++i) {
+        *(m_data+i) *= *(fft.m_data+i);
+    }
  
     // Return Fast Fourier Transform
     return *this;
@@ -206,12 +225,17 @@ GFft&amp; GFft::operator*=(const GFft&amp; fft)
  * @return Fast Fourier transform.
  *
  * Divides a Fast Fourier Transform by another Fast Fourier Transform.
- *
- * @todo Implement operator
+ * The method requires that both FFTs have the same dimension.
  ***************************************************************************/
 GFft& GFft::operator/=(const GFft& fft)
 {
-    // TODO: Implement operator
+    // Throw an exception if the FFTs have not the same shape
+    require_same_shape(G_OP_DIV, fft);
+
+    // Divide elements
+    for (int i = 0; i < m_size; ++i) {
+        *(m_data+i) /= *(fft.m_data+i);
+    }
  
     // Return Fast Fourier Transform
     return *this;
@@ -614,6 +638,75 @@ void GFft::set_data(const GNdarray&amp; array)
  
  
 /***********************************************************************//**
+ * @brief Check if FFT has the same shape
+ *
+ * @param[in] fft Fast Fourier Transform.
+ * @return True if the FFT has the same shape.
+ ***************************************************************************/
+bool GFft::has_same_shape(const GFft& fft) const
+{
+    // Check if the array has the same dimension
+    bool identity = m_shape.size() == fft.m_shape.size();
+
+    // Check if the shape is identical. Break as soon as one shape value
+    // differs.
+    if (identity) {
+        for (int i = 0; i < m_shape.size(); ++i) {
+            if (m_shape[i] != fft.m_shape[i]) {
+                identity = false;
+                break;
+            }
+        }
+    }
+
+    // Return result
+    return identity;
+}
+
+
+/***********************************************************************//**
+ * @brief Throw exception if FFT shapes differ
+ *
+ * @param[in] method Method that throws exception.
+ * @param[in] fft Fast Fourier Transform.
+ *
+ * @exception GException::invalid_argument
+ *            FFT shapes differ.
+ ***************************************************************************/
+void GFft::require_same_shape(const std::string& method,
+                              const GFft&        fft) const
+{
+    // If the shape differs then throw an exception
+    if (!has_same_shape(fft)) {
+
+        // Compose exception message
+        std::string msg = "Incompatible FFT dimensions (";
+        for (int i = 0; i < m_shape.size(); ++i) {
+            if (i > 0) {
+                msg += ", ";
+            }
+            msg += gammalib::str(m_shape[i]);
+        }
+        msg += ") <=> (";
+        for (int i = 0; i < fft.m_shape.size(); ++i) {
+            if (i > 0) {
+                msg += ", ";
+            }
+            msg += gammalib::str(fft.m_shape[i]);
+        }
+        msg += ").";
+
+        // Throw exception
+        throw GException::invalid_argument(method, msg);
+
+    } // endif: array shapes differed
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Perform Fast Fourier Transform
  *
  * @param[in] data Pointer to complex array to be transformed.
 # ==========================================================================
 # This module performs unit tests for the GammaLib numerics module.
 #
-# Copyright (C) 2012-2015 Juergen Knoedlseder
+# Copyright (C) 2012-2016 Juergen Knoedlseder
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -25,12 +25,12 @@ import gammalib
 # ======================================= #
 class Test(gammalib.GPythonTestSuite):
     """
-    Test class for GammaLib numerics module.
+    Test class for GammaLib numerics module
     """
     # Constructor
     def __init__(self):
         """
-        Constructor.
+        Constructor
         """
         # Call base class constructor
         gammalib.GPythonTestSuite.__init__(self)
@@ -41,21 +41,62 @@ class Test(gammalib.GPythonTestSuite):
     # Set test functions
     def set(self):
         """
-        Set all test functions.
+        Set all test functions
         """
         # Set test name
-        self.name("numerics")
+        self.name('numerics')
  
         # Append tests
-        self.append(self.test, "Numerics module dummy test")
+        self.append(self.test_fft, 'Test GFft')
  
         # Return
         return
  
-    # Test function
-    def test(self):
+    # Test Fast Fourier Transform
+    def test_fft(self):
         """
-        Test function.
+        Test function
         """
+        # Allocate and set 2-d array
+        array = gammalib.GNdarray(10, 10)
+        for i in range(3,7):
+            array[i,i] = 1.0
+        ref = array.sum()
+        
+        # Allocate and set 2-d kernel. The kernel needs to be normalized
+        # to unity and the centre of the kernel needs to be at pixel [0,0],
+        # and it needs to be wraped around to negative indices
+        kernel = gammalib.GNdarray(10, 10)
+        kernel[0,0] = 0.4
+        # 0.4
+        kernel[0,1] = 0.1
+        kernel[1,0] = 0.1
+        kernel[0,9] = 0.1
+        kernel[9,0] = 0.1
+        # 0.2
+        kernel[1,1] = 0.05
+        kernel[9,1] = 0.05
+        kernel[9,9] = 0.05
+        kernel[1,9] = 0.05
+
+        # Smooth 2-d array using FFT
+        fft_array  = gammalib.GFft(array)
+        fft_kernel = gammalib.GFft(kernel)
+        fft_smooth = fft_array * fft_kernel
+
+        # Backtransform
+        smooth = fft_smooth.backward()
+
+        # Test sum
+        sum = smooth.sum()
+        self.test_value(sum, ref)
+
+        # Store in sky map
+        map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,10,10)
+        for iy in range(10):
+            for ix in range(10):
+                map[ix+iy*10] = smooth[ix,iy]
+        map.save('test_fft.fits', True)
+        
         # Return
         return