Add convolution exercise

Author: gjbex

source-code/README.md

@@ -41,6 +41,7 @@ to create it.  There is some material not covered in the presentation as well.
 1. `pypy`: code to experiment with the Pypy interpreter.
 1. `file-formats`: influcence of file formats on performance.
 1. `performance`: general considerations about performance.
+1. `convolution`: wrap up exercise to apply all techniques.
 
 **Note:** the GPU code in this repository was moved to
 [its own repository](https://github.com/gjbex/Python-on-GPUs)

source-code/convolution/README.md

@@ -0,0 +1,21 @@
+# Convolution
+
+Convolution of an image using a kernel makes a nice problem to implement
+using various HPC technologies.  It is conceptually simple enough to be used
+as an exercise, yet computationally sufficiently challenging to make it
+interesting.
+
+To get you started, you get a
+  * [Python](python) implementations, and
+  * [C++](cpp) implementation. 
+
+
+You can try to:
+
+  * use numba,
+  * use Cython,
+  * use Swig to bind the C++ implementation,
+  * use PyBind11 to bind the C++ implementation,
+  * parallelize the code using Cython,
+  * parallelize the code using multiprocessing,
+  * parallelize and run the application on multiple nodes MPI.

source-code/convolution/cpp/CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.18)
+project(algorithms LANGUAGES CXX)
+
+set(CMAKE_CXX_STANDARD 23)
+set(CMaKE_CXX_STANDARD_REQUIRED YES)
+set(CMAKE_CXX_EXTENSIONS NO)
+
+add_compile_options(-Wall -Wextra -Wpedantic)
+
+add_subdirectory(src)

source-code/convolution/cpp/README.md

@@ -0,0 +1,12 @@
+# C++ implementation
+
+This is a (pretty naive) C++ implementation of 2D convolution.
+
+
+## What is it?
+
+1. `src/convolution/convolution.h`: declaration of the convolution function.
+1. `src/convolution/convolution.cpp`: definition of the convolution function.
+1. `src/benchmark_convolution.cpp`: a C++ application to benchmark the
+   implementation.
+1. `src/test_convolution.cpp`: a C++ application to test the implementation.

source-code/convolution/cpp/src/CMakeLists.txt

@@ -0,0 +1,15 @@
+add_subdirectory(convolution)
+
+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)
+
+add_executable(benchmark_convolution.exe
+        benchmark_convolution.cpp)
+target_include_directories(benchmark_convolution.exe
+    PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/convolution)
+target_link_libraries(benchmark_convolution.exe
+    PRIVATE convolution)

source-code/convolution/cpp/src/benchmark_convolution.cpp

@@ -0,0 +1,57 @@
+#include <convolution.h>
+#include <chrono>
+#include <iostream>
+#include <numeric>
+#include <random>
+
+Matrix create_image(int rows, int cols) {
+    std::mt19937 gen(1234);
+    std::uniform_real_distribution<double> dis(0.0, 1.0);
+    Matrix image(rows, cols);
+    for (int i = 0; i < rows; ++i) {
+        for (int j = 0; j < cols; ++j) {
+            image(i, j) = dis(gen);
+        }
+    }
+    return image;
+}
+
+Matrix create_kernel(int rows, int cols) {
+    Matrix kernel(rows, cols);
+    for (int i = 0; i < rows; ++i) {
+        for (int j = 0; j < cols; ++j) {
+            kernel(i, j) = 1.0/(rows*cols);
+        }
+    }
+    return kernel;
+}
+
+double element_sum(const Matrix& matrix) {
+    return std::accumulate(matrix.data(), matrix.data() + matrix.rows()*matrix.cols(), 0.0);
+}
+
+int main(int argc, char** argv) {
+    int rows = 1000;
+    int cols = 1000;
+    int kernel_rows = 7;
+    int kernel_cols = 7;
+    if (argc > 1) {
+        rows = atoi(argv[1]);
+        cols = atoi(argv[1]);
+    }
+    if (argc > 2) {
+        kernel_rows = atoi(argv[2]);
+        kernel_cols = atoi(argv[2]);
+    }
+    std::cout << "Image size: " << rows << "x" << cols << "\n";
+    std::cout << "Kernel size: " << kernel_rows << "x" << kernel_cols << "\n";
+    Matrix image = create_image(rows, cols);
+    Matrix kernel = create_kernel(kernel_rows, kernel_cols);
+    auto start = std::chrono::high_resolution_clock::now();
+    auto result = convolve(image, kernel);
+    auto end = std::chrono::high_resolution_clock::now();
+    std::chrono::duration<double> diff = end - start;
+    std::cout << "Time: " << diff.count() << " s\n";
+    std::cout << "Sum: " << element_sum(result) << "\n";
+    return 0;
+}

source-code/convolution/cpp/src/convolution/CMakeLists.txt

@@ -0,0 +1,2 @@
+add_library(convolution SHARED
+    convolution.cpp matrices.cpp)

source-code/convolution/cpp/src/convolution/convolution.cpp

@@ -0,0 +1,42 @@
+#include "convolution.h"
+#include <exception>
+
+Matrix convolve(const Matrix& image, 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
+      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.
+    */
+    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)};
+            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);
+                }
+            }
+            new_image(x, y) = value;
+        }
+    }
+    return new_image;
+}

source-code/convolution/cpp/src/convolution/convolution.h

@@ -0,0 +1,16 @@
+#ifndef CONVOLUTION_HDR
+#define CONVOLUTION_HDR
+
+#include "matrices.h"
+
+/**
+ * @brief Compute the convolution of an image with a kernel.
+ * @param image The image to convolve.  This is a 2D matrix with m rows and n columns.
+ * @param kernel The kernel to convolve with. This is a 2D matrix with k rows and l columns,
+ *               where k and l ard odd integers.
+ * @return The result of the convolution. This is a 2D matrix with m + k - 1 rows
+ *         and n + l -1 columns.
+ */
+Matrix convolve(const Matrix& image, const Matrix& kernel);
+
+#endif

source-code/convolution/cpp/src/convolution/matrices.cpp

@@ -0,0 +1,49 @@
+#include "matrices.h"
+
+Matrix::Matrix(const Matrix& other) :
+        rows_(other.rows_), cols_(other.cols_),
+        data_(new double[rows_ * cols_]) {
+    for (int i = 0; i < rows_ * cols_; ++i) {
+        data_[i] = other.data_[i];
+    }
+}
+
+Matrix& Matrix::operator=(const Matrix& other) {
+    if (this != &other) {
+        rows_ = other.rows_;
+        cols_ = other.cols_;
+        data_.reset(new double[rows_ * cols_]);
+        for (int i = 0; i < rows_ * cols_; ++i) {
+            data_[i] = other.data_[i];
+        }
+    }
+    return *this;
+}
+Matrix::Matrix(Matrix&& other) noexcept :
+        rows_{other.rows_}, cols_{other.cols_},
+        data_{std::move(other.data_)} {
+    other.rows_ = 0;
+    other.cols_ = 0;
+}
+
+Matrix& Matrix::operator=(Matrix&& other) noexcept {
+    if (&other != this) {
+        rows_ = other.rows_;
+        cols_ = other.cols_;
+        data_ = std::move(other.data_);
+
+        other.rows_ = 0;
+        other.cols_ = 0;
+    }
+    return *this;
+}
+
+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) {
+            os << m(i, j) << " ";
+        }
+        os << std::endl;
+    }
+    return os;
+}

source-code/convolution/cpp/src/convolution/matrices.h

@@ -0,0 +1,37 @@
+#ifndef MATRICES_HDR
+#define MATRICES_HDR
+
+#include <iostream>
+#include <memory>
+
+
+struct Matrix {
+    private:
+        int rows_;
+        int cols_;
+        std::unique_ptr<double[]> data_;
+    public:
+        Matrix(int rows, int 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]; }
+        // getters for number of rows and columns
+        int rows() const { return rows_; }
+        int cols() const { return cols_; }
+        // accessors for the data
+        double* data() { return data_.get(); }
+        const double* data() const { return data_.get(); }
+        // destructor
+        ~Matrix() = default;
+        // textual representation of the matrix
+        friend std::ostream& operator<<(std::ostream& os, const Matrix& m);
+};
+
+#endif

source-code/convolution/cpp/src/test_convolution.cpp

@@ -0,0 +1,34 @@
+#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;
+        }
+    }
+
+    // Print the image
+    std::cout << image << std::endl;
+
+    // Create a 3x3 kernel
+    Matrix kernel(3, 3);
+    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());
+        }
+    }
+
+    // Print the kernel
+    std::cout << kernel << std::endl;
+
+    // Create a convolution object
+    auto new_image = convolve(image, kernel);
+
+    // Print the result
+    std::cout << new_image << std::endl;
+
+    return 0;
+}

source-code/convolution/python/README.md

@@ -0,0 +1,10 @@
+# Python implementation
+
+This is a (pretty naive) Python implementation of 2D convolution.
+
+
+## What is it?
+
+1. `convolution.py`: module that contains the implementation.
+1. `benchmark_convolution.py`: a Python script to benchmark the implementation.
+1. `test_convolution.py`: a Python script to test the implementation.

source-code/convolution/python/benchmark_convolution.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python
+
+import argparse
+import convolution
+import numpy as np
+import timeit
+
+def create_image(rows, cols):
+    return np.random.uniform(0.0, 1.0, size=(rows, cols))
+
+def create_kernel(rows, cols):
+    return np.ones((rows, cols))/(rows*cols)
+
+def element_sum(matrix):
+    return np.sum(matrix)
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--rows", type=int, default=1_000)
+    parser.add_argument("--cols", type=int, default=1_000)
+    parser.add_argument("--kernel_rows", type=int, default=7)
+    parser.add_argument("--kernel_cols", type=int, default=7)
+    args = parser.parse_args()
+
+    print(f'Image size: {args.rows}x{args.cols}')
+    print(f'Kernel size: {args.kernel_rows}x{args.kernel_cols}')
+
+    image = create_image(args.rows, args.cols)
+    kernel = create_kernel(args.kernel_rows, args.kernel_cols)
+
+
+    print(f'Time: {timeit.timeit(lambda: convolution.convolve(image, kernel), number=1)} s')
+    print(f'Sum: {element_sum(convolution.convolve(image, kernel))}')
+
+if __name__ == "__main__":
+    main()