From 02ada43fbfcdff1e0d27a238c140f97f1506eb5e Mon Sep 17 00:00:00 2001
From: Dominik Kiese <40888446+dominikkiese@users.noreply.github.com>
Date: Thu, 26 Sep 2024 11:21:15 -0400
Subject: [PATCH] Hadamard product (#79)

Add hadamard function to calculate element-wise product between nda::Arrays, std::array, std::vector and Scalars.

Co-authored-by: Nils Wentzell <Wentzell@users.noreply.github.com>
Co-authored-by: Thomas Hahn <thomas.hahn01@gmail.com>
---
 c++/nda/algorithms.hpp      | 65 +++++++++++++++++++++++++++++++++++++
 c++/nda/stdutil/array.hpp   | 11 ++++---
 test/c++/nda_algorithms.cpp | 35 ++++++++++++++++++++
 3 files changed, 107 insertions(+), 4 deletions(-)

diff --git a/c++/nda/algorithms.hpp b/c++/nda/algorithms.hpp
index cee3b0885..6f2b1fa8b 100644
--- a/c++/nda/algorithms.hpp
+++ b/c++/nda/algorithms.hpp
@@ -23,14 +23,19 @@
 
 #include "./concepts.hpp"
 #include "./layout/for_each.hpp"
+#include "./layout/range.hpp"
+#include "./macros.hpp"
+#include "./map.hpp"
 #include "./traits.hpp"
 
 #include <algorithm>
+#include <array>
 #include <cmath>
 #include <cstdlib>
 #include <functional>
 #include <type_traits>
 #include <utility>
+#include <vector>
 
 namespace nda {
 
@@ -199,6 +204,66 @@ namespace nda {
     return fold(std::multiplies<>{}, a, get_value_t<A>{1});
   }
 
+  /**
+   * @brief Hadamard product of two nda::Array objects.
+   *
+   * @tparam A nda::Array type.
+   * @tparam B nda::Array type.
+   * @param a nda::Array object.
+   * @param b nda::Array object.
+   * @return A lazy nda::expr_call object representing the elementwise product of the two input objects.
+   */
+  template <Array A, Array B>
+    requires(nda::get_rank<A> == nda::get_rank<B>)
+  [[nodiscard]] constexpr auto hadamard(A &&a, B &&b) {
+    return nda::map([](auto const &x, auto const &y) { return x * y; })(std::forward<A>(a), std::forward<B>(b));
+  }
+
+  /**
+   * @brief Hadamard product of two std::array objects.
+   *
+   * @tparam T Value type of the first array.
+   * @tparam U Value type of the second array.
+   * @tparam R Size of the arrays.
+   * @param a std::array object.
+   * @param b std::array object.
+   * @return std::array containing the elementwise product of the two input arrays.
+   */
+  template <typename T, typename U, size_t R>
+  [[nodiscard]] constexpr auto hadamard(std::array<T, R> const &a, std::array<U, R> const &b) {
+    return a * b;
+  }
+
+  /**
+   * @brief Hadamard product of two std::vector objects.
+   *
+   * @tparam T Value type of the first input vector.
+   * @tparam U Value type of the second input vector.
+   * @param a std::vector object.
+   * @param b std::vector object.
+   * @return std::vector containing the elementwise product of the two input vectors.
+   */
+  template <typename T, typename U>
+  [[nodiscard]] constexpr auto hadamard(std::vector<T> const &a, std::vector<U> const &b) {
+    using TU = decltype(std::declval<T>() * std::declval<U>());
+    EXPECTS(a.size() == b.size());
+
+    std::vector<TU> c(a.size());
+    for (auto i : range(c.size())) c[i] = a[i] * b[i];
+    return c;
+  }
+
+  /**
+   * @brief Hadamard product of two arithmetic types.
+   *
+   * @tparam T nda::Scalar type of the first input.
+   * @tparam U nda::Scalar type of the second input.
+   * @param a First input.
+   * @param b Second input.
+   * @return Product of the two inputs.
+   */
+  constexpr auto hadamard(nda::Scalar auto a, nda::Scalar auto b) { return a * b; }
+
   /** @} */
 
 } // namespace nda
diff --git a/c++/nda/stdutil/array.hpp b/c++/nda/stdutil/array.hpp
index 61f26a0c6..9c355bca0 100644
--- a/c++/nda/stdutil/array.hpp
+++ b/c++/nda/stdutil/array.hpp
@@ -105,15 +105,18 @@ namespace std { // has to be in the right namespace for ADL
   /**
    * @brief Multiply two std::array objects element-wise.
    *
-   * @tparam T Value type of the arrays.
+   * @tparam T Value type of the first array.
+   * @tparam U Value type of the second array.
    * @tparam R Size of the arrays.
    * @param lhs Left hand side std::array operand.
    * @param rhs Right hand side std::array operand.
    * @return std::array containing the element-wise product.
    */
-  template <typename T, size_t R>
-  constexpr std::array<T, R> operator*(std::array<T, R> const &lhs, std::array<T, R> const &rhs) {
-    std::array<T, R> res;
+  template <typename T, typename U, size_t R>
+  constexpr auto operator*(std::array<T, R> const &lhs, std::array<U, R> const &rhs) {
+    using TU = decltype(std::declval<T>() * std::declval<U>());
+
+    std::array<TU, R> res;
     for (int i = 0; i < R; ++i) res[i] = lhs[i] * rhs[i];
     return res;
   }
diff --git a/test/c++/nda_algorithms.cpp b/test/c++/nda_algorithms.cpp
index bdb2ba6a8..69bfe2db9 100644
--- a/test/c++/nda_algorithms.cpp
+++ b/test/c++/nda_algorithms.cpp
@@ -116,3 +116,38 @@ TEST_F(NDAAlgorithm, CombineAlgorithmsWithArithmeticOps) {
   EXPECT_EQ(nda::sum(B), -18);
   EXPECT_EQ(nda::sum(A + B), 18);
 }
+
+TEST_F(NDAAlgorithm, HadamardProduct) {
+  nda::array<int, 2> A(3, 3), B(3, 3);
+  for (int i = 0; i < 3; ++i) {
+    for (int j = 0; j < 3; ++j) {
+      A(i, j) = i + 2 * j + 1;
+      B(i, j) = i - 3 * j;
+    }
+  }
+
+  // test with nda::array
+  auto Amat = nda::matrix<int>(A);
+  auto Bmat = nda::matrix<int>(B);
+  EXPECT_EQ(nda::hadamard(A, B), A * B);
+  EXPECT_EQ(nda::hadamard(Amat, Bmat), A * B);
+  EXPECT_EQ(nda::hadamard(Amat, B), A * B);
+
+  // test with std::array
+  auto arr1 = std::array<int, 3>{1, 2, 3};
+  auto arr2 = std::array<long, 3>{4, 5, 6};
+  EXPECT_EQ(nda::hadamard(arr1, arr2), arr1 * arr2);
+
+  // test with std::vector
+  auto vec1 = std::vector<long>{1, 2, 3};
+  auto vec2 = std::vector<double>{4, 5, 6};
+  auto vec3 = std::vector<double>{4, 10, 18};
+  EXPECT_EQ(nda::hadamard(vec1, vec2), vec3);
+
+  // test with arithmetic types or complex
+  EXPECT_EQ(nda::hadamard(2, 3.0), 6.0);
+
+  auto x = std::complex<long>{1, 2};
+  auto y = x * 3;
+  EXPECT_EQ(nda::hadamard(x, 3), y);
+}
\ No newline at end of file