Add convolution example

Author: gjbex

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/README.md

@@ -0,0 +1,31 @@
+# Convolution
+
+Illustration of various concepts using a naive convolution algorithm
+as an example.
+
+
+## What is it?
+
+1. `src`: source directory of C++ code.
+   a. `convolution`: directory that contains the convolution and matrix code.
+      1. `matrices.h`: definition of the Matrix class.
+      1. `matrices.cpp` definition of methods for the Matrix class.
+      1. `convolution.h`: decleartion of the convolution functions.
+      1. `convolution.cpp`: implementation of the convolution functions.
+      1. `CMakeLists.txt`: CMake file to build the shared library.
+   a. `bindings`: directory that contains the pybind11 bindings
+      for Matrix and convolve.
+      1. `matrix_bind.cpp`: pybind11 bindings for the Matrix class.
+      1. `convolve_bind.cpp`: pybind11 bindings for the convolve function.
+         This uses the buffer protocol.
+      1. `CMakeLists.txt`: CMake file to build the Python modules.
+   a. `test_convolution.cpp`: C++ application to test the C++ convolution
+      implementation.
+   a. benchmark_convolution.cpp`: C++ application to benchmark the C++ convolution
+      implementation.
+   a. `CMakeLists.txt`: CMake file to build the library, modules and C++
+      applications.
+1. `CMakeLists.txt`: CMake file to build the library, modules and C++
+   applications.
+1. `test_convolution.py`: Python script to test and benchmark the generated
+   modules.

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_subdirectory(convolution)
+add_subdirectory(bindings)
 
 add_executable(test_convolution.exe
         test_convolution.cpp)

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/bindings/CMakeLists.txt

@@ -0,0 +1,9 @@
+find_package(pybind11 REQUIRED)
+
+pybind11_add_module(matrices matrix_bind.cpp)
+target_include_directories(matrices PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../convolution)
+target_link_libraries(matrices PUBLIC convolution)
+
+pybind11_add_module(convolve convolve_bind.cpp)
+target_include_directories(convolve PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../convolution)
+target_link_libraries(convolve PUBLIC convolution)

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/bindings/convolve_bind.cpp

@@ -0,0 +1,38 @@
+#include <convolution.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+namespace py = pybind11;
+
+void convolve_func(py::buffer image_buffer, py::buffer kernel_buffer,
+        py::buffer new_image_buffer) {
+    py::buffer_info image_info = image_buffer.request();
+    if (image_info.format != py::format_descriptor<double>::format())
+        throw std::runtime_error("Incompatible format: expected a double array");
+    if (image_info.ndim != 2)
+        throw std::runtime_error("Incompatible buffer dimension");
+    py::buffer_info kernel_info = kernel_buffer.request();
+    if (kernel_info.format != py::format_descriptor<double>::format())
+        throw std::runtime_error("Incompatible format: expected a double array");
+    if (kernel_info.ndim != 2)
+        throw std::runtime_error("Incompatible buffer dimension");
+    py::buffer_info result_info = new_image_buffer.request();
+    if (result_info.format != py::format_descriptor<double>::format())
+        throw std::runtime_error("Incompatible format: expected a double array");
+    if (result_info.ndim != 2)
+        throw std::runtime_error("Incompatible buffer dimension");
+    if (result_info.shape[0] != image_info.shape[0] + kernel_info.shape[0] - 1 ||
+            result_info.shape[1] != image_info.shape[1] + kernel_info.shape[1] - 1)
+        throw std::runtime_error("Incompatible result buffer shape");
+    Matrix image(image_info.shape[0], image_info.shape[1]);
+    Matrix kernel(kernel_info.shape[0], kernel_info.shape[1]);
+    std::memcpy(image.data(), image_info.ptr, sizeof(double)*image.rows()*image.cols());
+    std::memcpy(kernel.data(), kernel_info.ptr, sizeof(double)*kernel.rows()*kernel.cols());
+    Matrix result = convolve(image, kernel);
+    std::memcpy(result_info.ptr, result.data(), sizeof(double)*result.rows()*result.cols());
+}
+
+PYBIND11_MODULE(convolve, module) {
+    module.doc() = "pybind11 wrapper module for convolution.h";
+    module.def("convolve", &convolve_func, "compute convolution of image with kernel");
+}

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/bindings/matrix_bind.cpp

@@ -0,0 +1,31 @@
+#include <matrices.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(matrices, module) {
+    module.doc() = "pybind11 wrapper module for matrices.h";
+    py::class_<Matrix>(module, "Matrix", py::buffer_protocol())
+        .def_buffer([](Matrix &m) -> py::buffer_info {
+            return py::buffer_info(
+                m.data(),                                /* Pointer to buffer */
+                sizeof(double),                          /* Size of one scalar */
+                py::format_descriptor<double>::format(), /* Python struct-style format descriptor */
+                2,                                       /* Number of dimensions */
+                {m.rows(), m.cols()},                    /* Buffer dimensions */
+                {sizeof(double)*m.cols(),                /* Strides (in bytes) for each index */
+                 sizeof(double)});
+        })
+        .def(py::init<int, int>(), "initialize matrix with given number of rows and columns")
+        .def(py::init([](const py::buffer b) {
+            py::buffer_info info = b.request();
+            if (info.format != py::format_descriptor<double>::format())
+                throw std::runtime_error("Incompatible format: expected a double array");
+            if (info.ndim != 2)
+                throw std::runtime_error("Incompatible buffer dimension!");
+            Matrix m(info.shape[0], info.shape[1]);
+            std::memcpy(m.data(), info.ptr, sizeof(double)*m.rows()*m.cols());
+            return m;
+        }), "initialize matrix from a numpy array");
+}

source-code/interfacing-c-c++-fortran/Pybind11/Convolution/src/convolution/convolution.h

@@ -3,6 +3,14 @@
 
 #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/interfacing-c-c++-fortran/Pybind11/Convolution/test_convolution.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python
+
+import argparse
+from build.src.bindings.convolve import convolve
+import numpy as np
+import timeit
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--img_size', type=int, default=10, help='number of rows in the image')
+    parser.add_argument('--kernel_size', type=int, default=3, help='number of rows in the kernel')
+    parser.add_argument('--num_iter', type=int, default=1, help='number of iterations to run')
+    args = parser.parse_args()
+
+    # Create input image and kernel
+    input_image = np.random.uniform(0.0, 1.0, size=(args.img_size, args.img_size))
+    kernel = np.random.rand(args.kernel_size, args.kernel_size)
+    result = np.empty((args.img_size + args.kernel_size - 1, args.img_size + args.kernel_size - 1))
+    time = timeit.timeit(lambda: convolve(input_image, kernel, result), number=args.num_iter)
+    print(f'Image size: {args.img_size}x{args.img_size}')
+    print(f'Kernel size: {args.kernel_size}x{args.kernel_size}')
+    print(f'Time: {time/args.num_iter:.6f} seconds')
+    print(f'Sum: {np.sum(result)}')
+
+if __name__ == '__main__':
+    main()

source-code/interfacing-c-c++-fortran/Pybind11/README.md

@@ -7,4 +7,5 @@ pybind11 is a wrapper generator for C++ code that has a lot of nice features.
 1. `Simple`: very simple illustration of wrapping C++ functions.
 1. `Spectrum`: illustration of warpping to support the buffer protocol.
 1. `Stats`: illustration of wrapping a C++ class.
+1. `Convolution`: illustration of using the buffer protocol.
 1. `environment.yml`: conda environment specification for this directory.