Merge pull request #21 from gjbex/development

gjbex · web-flow · commit b4fc7fd9fbf0 · 2024-12-09T15:48:45.000+01:00
Development
diff --git a/python_for_hpc.pptx b/python_for_hpc.pptx
diff --git a/source-code/convolution/cpp/CMakeLists.txt b/source-code/convolution/cpp/CMakeLists.txt
@@ -7,4 +7,6 @@ set(CMAKE_CXX_EXTENSIONS NO)
 
 add_compile_options(-Wall -Wextra -Wpedantic)
 
+enable_testing()
+
 add_subdirectory(src)
diff --git a/source-code/convolution/cpp/src/CMakeLists.txt b/source-code/convolution/cpp/src/CMakeLists.txt
@@ -1,11 +1,22 @@
 add_subdirectory(convolution)
 
-add_executable(test_convolution.exe
+find_package(Catch2)
+
+if(Catch2_FOUND)
+    add_executable(test_convolution.exe
         test_convolution.cpp)
-target_include_directories(test_convolution.exe
-    PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution)
-target_link_libraries(test_convolution.exe
-    PRIVATE convolution)
+    target_link_libraries(test_convolution.exe
+        PRIVATE Catch2::Catch2WithMain)
+    target_include_directories(test_convolution.exe
+        PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution)
+    target_link_libraries(test_convolution.exe
+        PRIVATE convolution)
+
+    add_test(NAME convolution
+        COMMAND test_convolution.exe)
+else()
+    message(STATUS "Catch2 not found, tests will not be built")
+endif()
 
 add_executable(benchmark_convolution.exe
         benchmark_convolution.cpp)
diff --git a/source-code/convolution/cpp/src/convolution/convolution.cpp b/source-code/convolution/cpp/src/convolution/convolution.cpp
@@ -1,42 +1,40 @@
 #include "convolution.h"
 #include <exception>
 
-Matrix convolve(const Matrix& image, const Matrix& kernel) {
+Matrix convolve(const Matrix& input, const Matrix& kernel) {
     if (kernel.rows() % 2 != 1 || kernel.cols() % 2 != 1) {
         throw std::invalid_argument("Only odd dimensions on kernel supported");
     }
     /*
-      s_mid and t_mid are number of pixels between the center pixel
+      kernel_row_mid and kernel_col_mid are number of pixels between the center pixel
       and the edge, ie for a 5x5 filter they will be 2.
       
-      The output size is calculated by adding s_mid, t_mid to each
-      side of the dimensions of the input image.
+      The output size is calculated by adding kernel_row_mid, kernel_col_mid to each
+      side of the dimensions of the input input.
     */
-    auto s_mid {kernel.rows()/2};
-    auto t_mid {kernel.cols()/2};
-    auto x_max {image.rows() + 2*s_mid};
-    auto y_max {image.cols() + 2*t_mid};
-    // Allocate result image.
-    Matrix new_image(x_max, y_max);
-    // Do convolution
-    for (int x = 0; x < x_max; ++x) {
-        for (int y = 0; y < y_max; ++y) {
-            // Calculate pixel value for h at (x,y). Sum one component
-            // for each pixel (s, t) of the filter kernel.
-            auto s_from {std::max(s_mid - x, -s_mid)};
-            auto s_to {std::min((x_max - x) - s_mid, s_mid + 1)};
-            auto t_from {std::max(t_mid - y, -t_mid)};
-            auto t_to {std::min((y_max - y) - t_mid, t_mid + 1)};
+    auto kernel_row_mid {kernel.rows()/2};
+    auto kernel_col_mid {kernel.cols()/2};
+    // Allocate result input.
+    Matrix output(input.rows(), input.cols());
+    // Do convolution.
+    for (int row = 0; row < input.rows(); ++row) {
+        auto min_kernel_row {std::max(0, kernel_row_mid - row)};
+        auto max_kernel_row {std::min(kernel.cols(), kernel_row_mid + input.rows() - row)};
+        for (int col = 0; col < input.cols(); ++col) {
+            auto min_kernel_col {std::max(0, kernel_col_mid - col)};
+            auto max_kernel_col {std::min(kernel.cols(), kernel_col_mid + input.cols() - col)};
             double value {0.0};
-            for (int s = s_from; s < s_to; ++s) {
-                for (int t = t_from; t < t_to; ++t) {
-                    auto v {x - s_mid + s};
-                    auto w {y - t_mid + t};
-                    value += kernel(s_mid - s, t_mid - t)*image(v, w);
+            // i and j are the kernel images.
+            for (int i = min_kernel_row; i < max_kernel_row; ++i) {
+                for (int j = min_kernel_col; j < max_kernel_col; ++j) {
+                    // k and l are the image indices.
+                    auto k {row - kernel_row_mid + i};
+                    auto l {col - kernel_col_mid + j};
+                    value += kernel(i, j)*input(k, l);
                 }
             }
-            new_image(x, y) = value;
+            output(row, col) = value;
         }
     }
-    return new_image;
+    return output;
 }
diff --git a/source-code/convolution/cpp/src/convolution/matrices.cpp b/source-code/convolution/cpp/src/convolution/matrices.cpp
@@ -19,6 +19,7 @@ Matrix& Matrix::operator=(const Matrix& other) {
     }
     return *this;
 }
+
 Matrix::Matrix(Matrix&& other) noexcept :
         rows_{other.rows_}, cols_{other.cols_},
         data_{std::move(other.data_)} {
@@ -38,6 +39,23 @@ Matrix& Matrix::operator=(Matrix&& other) noexcept {
     return *this;
 }
 
+bool Matrix::operator==(const Matrix& other) const {
+    if (this == &other) {
+        return true;
+    }
+    if (rows_ != other.rows_ || cols_ != other.cols_) {
+        return false;
+    }
+    for (int i = 0; i < rows_; ++i) {
+        for (int j = 0; j < cols_; ++j) {
+            if ((*this)(i, j) != other(i, j)) {
+                return false;
+            }   
+        }
+    }
+    return true;
+}
+
 std::ostream& operator<<(std::ostream& os, const Matrix& m) {
     for (int i = 0; i < m.rows_; ++i) {
         for (int j = 0; j < m.cols_; ++j) {
diff --git a/source-code/convolution/cpp/src/convolution/matrices.h b/source-code/convolution/cpp/src/convolution/matrices.h
@@ -12,19 +12,21 @@ struct Matrix {
         std::unique_ptr<double[]> data_;
     public:
         Matrix(int rows, int cols) :
-            rows_(rows), cols_(cols), data_(new double[rows * cols]) {}
+            rows_(rows), cols_(cols), data_(new double[rows*cols]) {}
         // copy constructor & assignment operator
         Matrix(const Matrix& other);
         Matrix& operator=(const Matrix& other);
         // move constructor & assignment operator
         Matrix(Matrix&& other) noexcept;
         Matrix& operator=(Matrix&& other) noexcept;
         // matrix indexing by row and column
-        double& operator()(int i, int j) { return data_[i * cols_ + j]; }
-        double operator()(int i, int j) const { return data_[i * cols_ + j]; }
+        double& operator()(int i, int j) { return data_[i*cols_ + j]; }
+        double operator()(int i, int j) const { return data_[i*cols_ + j]; }
         // getters for number of rows and columns
         int rows() const { return rows_; }
         int cols() const { return cols_; }
+        bool operator==(const Matrix& other) const;
+        bool operator!=(const Matrix& other) const { return !(*this == other); }
         // accessors for the data
         double* data() { return data_.get(); }
         const double* data() const { return data_.get(); }
diff --git a/source-code/convolution/cpp/src/test_convolution.cpp b/source-code/convolution/cpp/src/test_convolution.cpp
@@ -1,34 +1,134 @@
+#define CATCH_CONFIG_MAIN
+#include <catch2/catch.hpp>
 #include <convolution.h>
-#include <iostream>
-
-int main() {
-    // Create a 10x10 image
-    Matrix image(10, 10);
-    for (int i = 0; i < image.rows(); ++i) {
-        for (int j = 0; j < image.cols(); ++j) {
-            image(i, j) = i*image.cols() + j;
+
+TEST_CASE("3x3 kernel", "[convolution]") {
+    const int input_size {5};
+    const int kernel_size {3};
+
+    Matrix input(input_size, input_size);
+    for (int i = 0; i < input.rows(); ++i) {
+        for (int j = 0; j < input.cols(); ++j) {
+            input(i, j) = 1.0;
+        }
+    }
+
+    Matrix kernel(kernel_size, kernel_size);
+    for (int i = 0; i < kernel.rows(); ++i) {
+        for (int j = 0; j < kernel.cols(); ++j) {
+            kernel(i, j) = 1.0;
         }
     }
 
-    // Print the image
-    std::cout << image << std::endl;
+    Matrix target(input_size, input_size);
+    target(0, 0) = 4.0;
+    target(0, 1) = 6.0;
+    target(0, 2) = 6.0;
+    target(0, 3) = 6.0;
+    target(0, 4) = 4.0;
+    target(1, 0) = 6.0;
+    target(1, 1) = 9.0;
+    target(1, 2) = 9.0;
+    target(1, 3) = 9.0;
+    target(1, 4) = 6.0;
+    target(2, 0) = 6.0;
+    target(2, 1) = 9.0;
+    target(2, 2) = 9.0;
+    target(2, 3) = 9.0;
+    target(2, 4) = 6.0;
+    target(3, 0) = 6.0;
+    target(3, 1) = 9.0;
+    target(3, 2) = 9.0;
+    target(3, 3) = 9.0;
+    target(3, 4) = 6.0;
+    target(4, 0) = 4.0;
+    target(4, 1) = 6.0;
+    target(4, 2) = 6.0;
+    target(4, 3) = 6.0;
+    target(4, 4) = 4.0;
+
+    auto result = convolve(input, kernel);
+
+    REQUIRE(result == target);
+}
+
+TEST_CASE("5x5 kernel", "[convolution]") {
+    const int input_size {6};
+    const int kernel_size {5};
 
-    // Create a 3x3 kernel
-    Matrix kernel(3, 3);
+    Matrix input(input_size, input_size);
+    for (int i = 0; i < input.rows(); ++i) {
+        for (int j = 0; j < input.cols(); ++j) {
+            input(i, j) = 1.0;
+        }
+    }
+
+    Matrix kernel(kernel_size, kernel_size);
     for (int i = 0; i < kernel.rows(); ++i) {
         for (int j = 0; j < kernel.cols(); ++j) {
-            kernel(i, j) = 1.0/(kernel.rows()*kernel.cols());
+            kernel(i, j) = 1.0;
         }
     }
 
-    // Print the kernel
-    std::cout << kernel << std::endl;
+    Matrix target(input_size, input_size);
+    target(0, 0) = 9.0;
+    target(0, 1) = 12.0;
+    target(0, 2) = 15.0;
+    target(0, 3) = 15.0;
+    target(0, 4) = 12.0;
+    target(0, 5) = 9.0;
+    target(1, 0) = 12.0;
+    target(1, 1) = 16.0;
+    target(1, 2) = 20.0;
+    target(1, 3) = 20.0;
+    target(1, 4) = 16.0;
+    target(1, 5) = 12.0;
+    target(2, 0) = 15.0;
+    target(2, 1) = 20.0;
+    target(2, 2) = 25.0;
+    target(2, 3) = 25.0;
+    target(2, 4) = 20.0;
+    target(2, 5) = 15.0;
+    target(3, 0) = 15.0;
+    target(3, 1) = 20.0;
+    target(3, 2) = 25.0;
+    target(3, 3) = 25.0;
+    target(3, 4) = 20.0;
+    target(3, 5) = 15.0;
+    target(4, 0) = 12.0;
+    target(4, 1) = 16.0;
+    target(4, 2) = 20.0;
+    target(4, 3) = 20.0;
+    target(4, 4) = 16.0;
+    target(4, 5) = 12.0;
+    target(5, 0) = 9.0;
+    target(5, 1) = 12.0;
+    target(5, 2) = 15.0;
+    target(5, 3) = 15.0;
+    target(5, 4) = 12.0;
+    target(5, 5) = 9.0;
+    auto result = convolve(input, kernel);
+
+    REQUIRE(result == target);
+}
 
-    // Create a convolution object
-    auto new_image = convolve(image, kernel);
+TEST_CASE("even-sized kernel error", "[convolution]") {
+    const int input_size {5};
+    const int kernel_size {4};
 
-    // Print the result
-    std::cout << new_image << std::endl;
+    Matrix input(input_size, input_size);
+    for (int i = 0; i < input.rows(); ++i) {
+        for (int j = 0; j < input.cols(); ++j) {
+            input(i, j) = 1.0;
+        }
+    }
 
-    return 0;
+    Matrix kernel(kernel_size, kernel_size);
+    for (int i = 0; i < kernel.rows(); ++i) {
+        for (int j = 0; j < kernel.cols(); ++j) {
+            kernel(i, j) = 1.0;
+        }
+    }
+    REQUIRE_THROWS_AS(convolve(input, kernel), std::invalid_argument);
 }
+
diff --git a/source-code/convolution/python/convolution.py b/source-code/convolution/python/convolution.py
@@ -14,7 +14,7 @@ def convolve(input, kernel):
     Returns
     -------
     output : numpy.ndarray
-        output matrix, it is not cropped
+        output matrix, shape is the same as input
 
     Raises
     ------
@@ -23,31 +23,28 @@ def convolve(input, kernel):
     '''
     if kernel.shape[0] % 2 != 1 or kernel.shape[1] % 2 != 1:
         raise ValueError("Only odd dimensions on filter supported")
-    # s_mid and t_mid are number of pixels between the center pixel
-    # and the edge, ie for a 5x5 filter they will be 2.
-    #
-    # The output size is calculated by adding s_mid, t_mid to each
-    # side of the dimensions of the input image.
-    s_mid = kernel.shape[0] // 2
-    t_mid = kernel.shape[1] // 2
-    x_max = input.shape[0] + 2 * s_mid
-    y_max = input.shape[1] + 2 * t_mid
-    # Allocate result image.
-    output = np.zeros([x_max, y_max], dtype=input.dtype)
-    # Do convolution
-    for x in range(x_max):
-        for y in range(y_max):
-            # Calculate pixel value for h at (x,y). Sum one component
-            # for each pixel (s, t) of the filter kernel.
-            s_from = max(s_mid - x, -s_mid)
-            s_to = min((x_max - x) - s_mid, s_mid + 1)
-            t_from = max(t_mid - y, -t_mid)
-            t_to = min((y_max - y) - t_mid, t_mid + 1)
-            value = 0
-            for s in range(s_from, s_to):
-                for t in range(t_from, t_to):
-                    v = x - s_mid + s
-                    w = y - t_mid + t
-                    value += kernel[s_mid - s, t_mid - t]*input[v, w]
-            output[x, y] = value
+    # kernel_row_mid and kernel_col_mid are number of elements between the
+    # center elements # and the edge, i.e., for a 5x5 filter they will be 2.
+    kernel_row_mid = kernel.shape[0]//2
+    kernel_col_mid = kernel.shape[1]//2
+    # Allocate result input.
+    output = np.empty_like(input)
+    # Do convolution.
+    for row in range(input.shape[0]):
+        min_kernel_row = max(0, kernel_row_mid - row)
+        max_kernel_row = min(kernel.shape[0],
+                             kernel_row_mid + input.shape[0] - row)
+        for col in range(input.shape[1]):
+            min_kernel_col = max(0, kernel_col_mid - col)
+            max_kernel_col = min(kernel.shape[1],
+                                 kernel_col_mid + input.shape[1] - col)
+            value = 0.0
+            # i and j are the kernel images
+            for i in range(min_kernel_row, max_kernel_row):
+                for j in range(min_kernel_col, max_kernel_col):
+                    # k and l are the image indices
+                    k = row - kernel_row_mid + i
+                    l = col - kernel_col_mid + j
+                    value += kernel[i, j]*input[k, l]
+            output[row, col] = value
     return output
diff --git a/source-code/convolution/python/test_convolution.py b/source-code/convolution/python/test_convolution.py
