feat: vector algebra (for parallelization)

src/CMakeLists.txt

@@ -3,6 +3,7 @@ add_subdirectory(datastructures)
 add_subdirectory(algorithms)
 add_subdirectory(control)
 add_subdirectory(strings)
+add_subdirectory(math)
 
 include_directories(${CMAKE_CURRENT_SOURCE_DIR})
 
@@ -13,4 +14,5 @@ target_link_libraries(main
   algorithms_lib 
   control_lib 
   strings_lib
+  math_lib
 )

src/datastructures/vector.cc

@@ -8,6 +8,14 @@ OpsVector::PrintVector(std::vector<int> &v) {
   std::cout << std::endl;
 }
 
+void
+OpsVector::PrintVector(std::vector<double> &v) {
+  for (int i = 0; i < (int) v.size(); i += 1) {
+    std::cout << v[i] << " ";
+  }
+  std::cout << std::endl;
+}
+
 /**
  * @brief Add 1 to each element in the vector.
  *

src/datastructures/vector.h

@@ -8,6 +8,7 @@
 namespace OpsVector {
 
 void PrintVector(std::vector<int> &v);
+void PrintVector(std::vector<double> &v);
 std::vector<int> ModifyVector(std::vector<int> &v);
 std::vector<int> SearchVector(std::vector<int> &v, int n);
 std::vector<int> SortVector(std::vector<int> &v);

src/generator/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_library(generator_lib
+  genarray.cc
   genvector.cc
   genstring.cc
   genlist.cc

src/generator/genarray.cc

@@ -0,0 +1,22 @@
+#include "genarray.h"
+
+double *
+GenArray::RandomArray(int n, int m) {
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_real_distribution<> dis(0, m);
+  double *arr = new double[n];
+  for (int i = 0; i < n; i++) {
+    arr[i] = dis(gen);
+  }
+  return arr;
+}
+
+double *
+GenArray::ZeroArray(int n) {
+  double *arr = new double[n];
+  for (int i = 0; i < n; i++) {
+    arr[i] = 0;
+  }
+  return arr;
+}

src/generator/genarray.h

@@ -0,0 +1,13 @@
+#ifndef __GENERATOR_GENARRAY_H__
+#define __GENERATOR_GENARRAY_H__
+
+#include <random>
+
+namespace GenArray {
+
+double *RandomArray(int n, int m);
+double *ZeroArray(int n);
+
+};   // namespace GenArray
+
+#endif   // __GENERATOR_GENARRAY_H__

src/generator/genvector.cc

@@ -1,4 +1,6 @@
 #include "genvector.h"
+#include <random>
+#include <vector>
 
 /** @brief Generates a random vector of length n
  *
@@ -18,6 +20,27 @@ GenVector::RandomVector(int n, int m) {
   return ret;
 }
 
+/** @brief Generates a random vector of length n
+ *
+ *  @param n Length of the vector to generate
+ *  @param l Lower bound (inclusive) of the values in the vector
+ *  @param u Upper bound (non-inclusive) of the values in the vector
+ *  @return A random vector of length n
+ */
+std::vector<double>
+GenVector::RandomVectorDouble(int n, int l, int u) {
+  std::vector<double> ret = std::vector<double>(n);
+
+  std::mt19937 gen(0);
+  std::uniform_real_distribution<double> dis(l, u);
+
+  for (int i = 0; i < n; i += 1) {
+    ret[i] = dis(gen);
+  }
+
+  return ret;
+}
+
 /** @brief Generates a random square matrix of size n
  *
  *  @param n Size of the matrix to generate

src/generator/genvector.h

@@ -1,14 +1,15 @@
-#ifndef __GENERATOR_GENARRAY_H__
-#define __GENERATOR_GENARRAY_H__
+#ifndef __GENERATOR_GENVECTOR_H__
+#define __GENERATOR_GENVECTOR_H__
 
 #include <random>
 #include <vector>
 
 namespace GenVector {
 
 std::vector<int> RandomVector(int n, int m);
+std::vector<double> RandomVectorDouble(int n, int l, int u);
 std::vector<std::vector<int>> RandomSquareMatrix(int n, int m);
 
 };   // namespace GenVector
 
-#endif   // __GENERATOR_GENARRAY_H__
+#endif   // __GENERATOR_GENVECTOR_H__

src/main.cc

@@ -6,8 +6,10 @@
 #include "control/single.h"
 #include "datastructures/linkedlist.h"
 #include "datastructures/vector.h"
+#include "generator/genarray.h"
 #include "generator/genstring.h"
 #include "generator/genvector.h"
+#include "math/vectoralgebra.h"
 #include "strings/strops.h"
 
 void
@@ -159,6 +161,26 @@ _linkedlist() {
   std::cout << std::endl;
 }
 
+void
+_vectoralgebra() {
+  std::cout << "Vector Algebra" << std::endl;
+  std::cout << "--------------" << std::endl;
+
+  int n = 20;
+  std::vector<double> x = GenVector::RandomVectorDouble(n, -10, 10);
+  std::vector<double> h = GenVector::RandomVectorDouble(3, -1, 1);
+  OpsVector::PrintVector(x);
+  OpsVector::PrintVector(h);
+
+  std::vector<double> y = MathVectorAlgebra::Cubic(1, 2, 3, 4, x);
+  std::cout << "Cubic(1, 2, 3, 4, x): ";
+  OpsVector::PrintVector(y);
+
+  std::vector<double> z = MathVectorAlgebra::Convolve(1, h, x);
+  std::cout << "Convolve(1, h, x): ";
+  OpsVector::PrintVector(z);
+}
+
 int
 main(int argc, char **argv) {
   _single();
@@ -168,5 +190,6 @@ main(int argc, char **argv) {
   _sort();
   _vector();
   _linkedlist();
+  _vectoralgebra();
   return 0;
 }

src/math/CMakeLists.txt

@@ -0,0 +1,4 @@
+add_library(math_lib
+  vectoralgebra.cc
+)
+

src/math/vectoralgebra.cc

@@ -0,0 +1,56 @@
+#include "vectoralgebra.h"
+
+#include <cmath>
+
+/** @brief solve a cubic equation
+ *
+ * @param a Coefficient of x^3
+ * @param b Coefficient of x^2
+ * @param c Coefficient of x
+ * @param d Constant term
+ * @param x Vector of x values (input)
+ *
+ * @return Vector of y values (output)
+ */
+std::vector<double>
+MathVectorAlgebra::Cubic(int a, int b, int c, int d, std::vector<double> &x) {
+  int n = x.size();
+  std::vector<double> y = std::vector<double>(n);
+  for (int i = 0; i < n; i += 1) {
+    y[i] = a * x[i] * x[i] * x[i] + b * x[i] * x[i] + c * x[i] + d;
+  }
+  return y;
+}
+
+/** @brief convolve vector x with smaller vector h
+ *
+ * offset defined such that h[offset] is the center of h
+ *
+ * @param offset Offset of h
+ * @param h Pointer to array of h values (input)
+ * @param x Pointer to array of x values (input)
+ *
+ * @return Pointer to array of y values (output)
+ */
+std::vector<double>
+MathVectorAlgebra::Convolve(int offset, std::vector<double> &h, std::vector<double> &x) {
+  int m = h.size();
+  int n = x.size();
+  int sum_h = 0;
+  std::vector<double> y = std::vector<double>(n);
+  for (int i = 0; i < m; i += 1) {
+    sum_h += std::abs(h[i]);
+  }
+  if (sum_h == 0) {
+    sum_h = 1;
+  }
+  for (int i = 0; i < n; i += 1) {
+    y[i] = 0;
+    for (int j = 0; j < m; j += 1) {
+      if (i - offset + j >= 0 && i - offset + j < n) {
+        y[i] += h[j] * x[i - offset + j] / sum_h;
+      }
+    }
+  }
+  return y;
+}

src/math/vectoralgebra.h

@@ -0,0 +1,13 @@
+#ifndef __MATH_VECTORALGEBRA_H__
+#define __MATH_VECTORALGEBRA_H__
+
+#include <vector>
+
+namespace MathVectorAlgebra {
+
+std::vector<double> Cubic(int a, int b, int c, int d, std::vector<double> &x);
+std::vector<double> Convolve(int offset, std::vector<double> &h, std::vector<double> &x);
+
+};   // namespace MathVectorAlgebra
+
+#endif   // !__MATH_VECTORALGEBRA_H__

tests/CMakeLists.txt

@@ -17,6 +17,7 @@ add_executable(tests
   datastructures/test_vector.cc
   datastructures/test_linkedlist.cc
   datastructures/test_list.cc
+  math/test_vectoralgebra.cc
 )
 target_link_libraries(tests 
   PRIVATE 
@@ -27,4 +28,5 @@ target_link_libraries(tests
     algorithms_lib 
     control_lib 
     strings_lib
+    math_lib
 )

tests/math/test_vectoralgebra.cc

@@ -0,0 +1,67 @@
+#include <catch2/benchmark/catch_benchmark.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include "generator/genvector.h"
+#include "math/vectoralgebra.h"
+
+TEST_CASE("VectorAlgebra::Cubic::Test", "[test][math][vectoralgebra]") {
+  std::vector<double> x0 = {1, 2, 3, 4, 5};
+  std::vector<double> x1 = {0, 1, 2, 3, 4, 10};
+  std::vector<double> y0 = MathVectorAlgebra::Cubic(1, -1, -2, 4, x0);
+  std::vector<double> y1 = MathVectorAlgebra::Cubic(1, 0, 0, 10, x1);
+  std::vector<double> y2 = MathVectorAlgebra::Cubic(0, 5, -3, -4, x1);
+  std::vector<double> y3 = MathVectorAlgebra::Cubic(1, 3, -4, 4, x1);
+  std::vector<double> y0ref = {2, 4, 16, 44, 94};
+  std::vector<double> y1ref = {10, 11, 18, 37, 74, 1010};
+  std::vector<double> y2ref = {-4, -2, 10, 32, 64, 466};
+  for (int i = 0; i < 5; i++) {
+    REQUIRE(y0[i] == y0ref[i]);
+  }
+  for (int i = 0; i < 5; i++) {
+    REQUIRE(y1[i] == y1ref[i]);
+  }
+  for (int i = 0; i < 5; i++) {
+    REQUIRE(y2[i] == y2ref[i]);
+  }
+  for (int i = 0; i < 5; i++) {
+    REQUIRE(y3[i] != y2ref[i]);
+  }
+}
+
+TEST_CASE("VectorAlgebra::Cubic::Benchmark", "[benchmark][math][vectoralgebra]") {
+  std::vector<double> x0 = GenVector::RandomVectorDouble(1000, -10, 10);
+  BENCHMARK("Cubic(1000)") { return MathVectorAlgebra::Cubic(1, -1, -2, 4, x0); };
+  std::vector<double> x1 = GenVector::RandomVectorDouble(10000, -10, 10);
+  BENCHMARK("Cubic(10000)") { return MathVectorAlgebra::Cubic(-2, 1, 3, -8, x1); };
+  std::vector<double> x2 = GenVector::RandomVectorDouble(100000, -10, 10);
+  BENCHMARK("Cubic(10000)") { return MathVectorAlgebra::Cubic(3, -1, 1, 3, x2); };
+}
+
+TEST_CASE("VectorAlgebra::Convolve::Test", "[test][math][vectoralgebra]") {
+  std::vector<double> x0 = {1, 2, 3, 4, 5};
+  std::vector<double> h0 = {-1, 5, -1};
+  std::vector<double> h1 = {-1, 4, -1};
+
+  std::vector<double> y0 = MathVectorAlgebra::Convolve(1, h0, x0);
+  std::vector<double> y1 = MathVectorAlgebra::Convolve(1, h1, x0);
+  std::vector<double> y0ref = {0.428, 0.857, 1.286, 1.714, 3.000};
+  std::vector<double> y1ref = {0.333, 0.667, 1.000, 1.333, 2.667};
+  for (int i = 0; i < 5; i++) {
+    REQUIRE(std::abs(y0[i] - y0ref[i]) < 0.002);
+  }
+  for (int i = 0; i < 5; i++) {
+    REQUIRE(std::abs(y1[i] - y1ref[i]) < 0.002);
+  }
+}
+
+TEST_CASE("VectorAlgebra::Convolve::Benchmark", "[benchmark][math][vectoralgebra]") {
+  std::vector<double> x0 = GenVector::RandomVectorDouble(1000, -10, 10);
+  std::vector<double> h0 = GenVector::RandomVectorDouble(3, -1, 1);
+  BENCHMARK("Convolve(1000, 3)") { return MathVectorAlgebra::Convolve(1, h0, x0); };
+  std::vector<double> x1 = GenVector::RandomVectorDouble(10000, -10, 10);
+  std::vector<double> h1 = GenVector::RandomVectorDouble(3, -1, 1);
+  BENCHMARK("Convolve(10000, 3)") { return MathVectorAlgebra::Convolve(1, h1, x1); };
+  std::vector<double> x2 = GenVector::RandomVectorDouble(10000, -10, 10);
+  std::vector<double> h2 = GenVector::RandomVectorDouble(15, -1, 1);
+  BENCHMARK("Convolve(10000, 15)") { return MathVectorAlgebra::Convolve(7, h2, x2); };
+}