Repository: metric-space-ai/metric (language: python)

Issue:
Add zero derivative padding for wavelet convolution matrix
Add `PaddingType` parameter to 
`template <typename T>
blaze::CompressedMatrix<T> DaubechiesMat(size_t size, int order = 4)` function



Discussion / hints:

Fix patch:
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 498fc5f5..8b2cf5d4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -51,7 +51,7 @@ Include(FetchContent)
 FetchContent_Declare(
 		Catch2
 		GIT_REPOSITORY https://github.com/catchorg/Catch2.git
-		GIT_TAG        v2.13.4)
+		GIT_TAG        v2.13.7)
 
 FetchContent_MakeAvailable(Catch2)
 
diff --git a/modules/transform/wavelet.cpp b/modules/transform/wavelet.cpp
index 015548bf..524e8121 100644
--- a/modules/transform/wavelet.cpp
+++ b/modules/transform/wavelet.cpp
@@ -666,67 +666,141 @@ typename std::enable_if<blaze::IsMatrix<Container2d>::value, Container2d>::type
 
 // ------------------------------ DWT based on matrix multiplication
 
+/**
+ * Construct convolutional matrix for wavelet transform
+ * @tparam T float type
+ * @param size
+ * @param order
+ * @param padding padding type
+ * @return
+ */
 template <typename T>
-blaze::CompressedMatrix<T> DaubechiesMat(size_t size, int order = 4)
-{  // Daubechies Transform matrix generator
-
+blaze::CompressedMatrix<T> DaubechiesMat(size_t size, int order = 4, Padding padding = Padding::Periodized)
+{
     assert(order % 2 == 0);
-	assert(size >= order);
+    assert(size >= order);
 
-	auto [Lo_D, Hi_D, Lo_R, Hi_R] = orthfilt(dbwavf<std::vector<T>>(order / 2, T()));
+    auto [Lo_D, Hi_D, Lo_R, Hi_R] = orthfilt(dbwavf<std::vector<T>>(order / 2, T()));
 
-/* Reverse filters for convolution */
-	std::reverse(Lo_D.begin(), Lo_D.end());
-	std::reverse(Hi_D.begin(), Hi_D.end());
+    /* Reverse filters for convolution */
+    std::reverse(Lo_D.begin(), Lo_D.end());
+    std::reverse(Hi_D.begin(), Hi_D.end());
 
+    /* Low filter part */
     auto mat = blaze::CompressedMatrix<T>(size, size);
     mat.reserve(size * Lo_D.size());
-    for (size_t i = 0; i < size / 2; ++i) {
+    if (padding == Padding::Periodized) {
+        for (size_t i = 0; i < size / 2; ++i) {
 
-        size_t ci = mat.columns() - 2 * i;
-        if (ci > Lo_D.size()) {
-          ci = 0;
+            size_t ci = mat.columns() - 2 * i;
+            if (ci > Lo_D.size()) {
+                ci = 0;
+            }
+            for (size_t a = 0; a < Lo_D.size(); ++a) {
+                if (ci >= Lo_D.size()) {
+                    ci = ci % Lo_D.size();
+                }
+
+                size_t j = i * 2 + ci;
+                if (j >= mat.columns()) {
+                    j = j % mat.columns();
+                }
+
+                mat.append(i, j, Lo_D[ci]);
+                ++ci;
+            }
+            mat.finalize(i);
         }
-        for (size_t a = 0; a < Lo_D.size(); ++a) {
-          if (ci >= Lo_D.size()) {
-            ci = ci % Lo_D.size();
-          } 
-
-          size_t j = i * 2 + ci;
-          if (j >= mat.columns()) {
-            j = j % mat.columns();
-          }
-
-          mat.append(i, j, Lo_D[ci]);
-          
-          ++ci;
+
+        /* Hi filter part */
+        for (size_t i = 0; i < size / 2; ++i) {
+            size_t ci = mat.columns() - 2 * i;
+            if (ci > Hi_D.size()) {
+                ci = 0;
+            }
+            for (size_t a = 0; a < Hi_D.size(); ++a) {
+                if (ci >= Hi_D.size()) {
+                    ci = ci % Hi_D.size();
+                }
+
+                size_t j = i * 2 + ci;
+                if (j >= mat.columns()) {
+                    j = j % mat.columns();
+                }
+
+                mat.append(size / 2 + i, j, Hi_D[ci]);
+                ++ci;
+            }
+            mat.finalize(size / 2 + i);
         }
-        mat.finalize(i);
+    } else if (padding == Padding::ZeroDerivative) {
 
-    }
+        /* Calculate padding size */
+        int paddingSize = Lo_D.size() - 1;
 
-    for (size_t i = 0; i < size / 2; ++i) {
+        /* Low filter part */
+        for (size_t i = 0; i < size / 2; ++i) {
+            int leftPadding = paddingSize / 2;
+            int j0 = -leftPadding + 2 * static_cast<int>(i);
+
+            /* Left padding */
+            T lp = 0;
+            for (int k = 0; k < -j0; ++k) {
+                lp += Lo_D[k];
+            }
 
+            /* Right padding */
+            T rp = 0;
+            int l = j0 + Lo_D.size() - size;
+            for (int k = 0; k < l; ++k) {
+                rp += Lo_D[Lo_D.size() - 1 - k];
+            }
 
-        size_t ci = mat.columns() - 2 * i;
-        if (ci > Hi_D.size()) {
-          ci = 0;
+            for (int k = 0; k < Lo_D.size(); ++k) {
+                const int j = j0 + k;
+                if (j == 0) {
+                    mat.append(i, j, lp + Lo_D[k]);
+                } else if (j == size - 1) {
+                    mat.append(i, j, rp + Lo_D[k]);
+                } else if (j > 0 and j < size - 1) {
+                    mat.append(i, j, Lo_D[k]);
+                } else {}
+            }
+            mat.finalize(i);
         }
-        for (size_t a = 0; a < Hi_D.size(); ++a) {
-          if (ci >= Hi_D.size()) {
-            ci = ci % Hi_D.size();
-          } 
-
-          size_t j = i * 2 + ci;
-          if (j >= mat.columns()) {
-            j = j % mat.columns();
-          }
-
-          mat.append(size / 2 + i, j, Hi_D[ci]);
-          
-          ++ci;
+
+        /* Hi filter part */
+        for (size_t i = 0; i < size / 2; ++i) {
+            int leftPadding = paddingSize / 2;
+            int j0 = -leftPadding + 2 * static_cast<int>(i);
+
+            /* Left padding */
+            T lp = 0;
+            for (int k = 0; k < -j0; ++k) {
+                lp += Hi_D[k];
+            }
+
+            /* Right padding */
+            T rp = 0;
+            int l = j0 + Hi_D.size() - size;
+            for (int k = 0; k < l; ++k) {
+                rp += Hi_D[Hi_D.size() - 1 - k];
+            }
+
+            for (int k = 0; k < Hi_D.size(); ++k) {
+                const int j = j0 + k;
+                if (j == 0) {
+                    mat.append(size / 2 + i, j, lp + Hi_D[k]);
+                } else if (j == size - 1) {
+                    mat.append(size / 2 + i, j, rp + Hi_D[k]);
+                } else if (j > 0 and j < size - 1) {
+                    mat.append(size / 2 + i, j, Hi_D[k]);
+                } else {}
+            }
+            mat.finalize(size / 2 + i);
         }
-        mat.finalize(size / 2 + i);
+    } else {
+        return mat;
     }
 
     return mat;
diff --git a/modules/transform/wavelet.hpp b/modules/transform/wavelet.hpp
index d2a45796..e7d3a6fc 100644
--- a/modules/transform/wavelet.hpp
+++ b/modules/transform/wavelet.hpp
@@ -90,6 +90,13 @@ using InternalBlazeMatT = typename InternalBlazeType<Container2d, SO>::matrix_ty
 
 // wavelet functions
 
+/**
+ * Padding type
+ */
+enum class Padding {
+    ZeroDerivative,
+    Periodized
+};
 
 ///**
 // * @brief valid convolution


Test patch:
diff --git a/tests/transform_test/wavelet_tests.cpp b/tests/transform_test/wavelet_tests.cpp
index 9e3e767f..93cf386b 100644
--- a/tests/transform_test/wavelet_tests.cpp
+++ b/tests/transform_test/wavelet_tests.cpp
@@ -30,6 +30,26 @@ TEMPLATE_TEST_CASE("DaubechiesMat", "[transform][wavelet]", float, double)
   }
 }
 
+TEMPLATE_TEST_CASE("DaubechiesMat()_ZeroDerivative", "[transform][wavelet]", float, double)
+{
+   const blaze::DynamicVector<TestType> data = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
+
+   const auto dmat = wavelet::DaubechiesMat<TestType>(data.size(), 4, wavelet::Padding::ZeroDerivative);
+   const auto result = dmat * data;
+
+   const blaze::DynamicVector<TestType> expected = {{1.379538}, {3.725003}, {6.553430}, {9.381957}, {12.339693},
+                                                    {-0.129410}, {-9.901802e-13}, {-9.9e-13}, {-9.9e-13}, {0.482963}};
+
+   REQUIRE(result.size() == expected.size());
+   for (size_t i = 0; i < expected.size(); ++i) {
+       if (result[i] > 1e-10) {
+           REQUIRE(expected[i] == Approx(result[i]));
+       }
+   }
+
+
+}
+
 TEMPLATE_TEST_CASE("dbwavf", "[transform][wavelet]", float, double)
 {
 	using Vector = blaze::DynamicVector<TestType>;