draft for hadamard product

c++/nda/algorithms.hpp

@@ -24,6 +24,7 @@
 #include "./concepts.hpp"
 #include "./layout/for_each.hpp"
 #include "./traits.hpp"
+#include "./map.hpp"
 
 #include <algorithm>
 #include <cmath>
@@ -199,6 +200,72 @@ 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 nda::Array containing the elementwise product of the two input arrays.
+   */
+  template <Array A, Array B>
+    requires(nda::get_rank<A> == nda::get_rank<B>)
+  [[nodiscard]] constexpr auto hadamard(A &&a, B &&b) {
+    assert(a.shape() == b.shape());
+    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 Data type of the arrays.
+   * @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, size_t R>
+  [[nodiscard]] constexpr auto hadamard(std::array<T, R> const &a, std::array<T, R> const &b) {
+    assert(a.size() == b.size());
+    return a * b;
+  }
+
+  /**
+   * @brief Hadamard product of two std::vector objects.
+   * 
+   * @tparam T Data type of the first input vector.
+   * @tparam U Data 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>());
+    assert(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 Data type of the first input.
+   * @tparam U Data type of the second input.
+   * @param a First input.
+   * @param b Second input.
+   * @return Product of the two inputs.
+   */
+  template <typename T, typename U>
+  constexpr auto hadamard(T a, U b)
+    requires(std::is_arithmetic_v<T> && std::is_arithmetic_v<U>)
+  {
+    return a * b;
+  }
+
   /** @} */
 
 } // namespace nda

test/c++/nda_algorithms.cpp

@@ -116,3 +116,34 @@ 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<int, 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
+  EXPECT_EQ(nda::hadamard(2, 3.0), 6.0);
+}