Dynamic Programming c++ implementation

CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 3.8)
+project(_dupin VERSION 0.0.1)
+
+set(DEFAULT_BUILD_TYPE "Release")
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE ${DEFAULT_BUILD_TYPE})
+endif()
+
+find_package(Eigen3 REQUIRED)
+find_package(TBB REQUIRED)
+find_package(pybind11 CONFIG REQUIRED)
+
+include_directories(${PROJECT_SOURCE_DIR}/src)
+add_subdirectory(src)

dupin/detect/dynp.py

@@ -0,0 +1,52 @@
+"""Implements dynamic programming class for optimal segementation algorithm."""
+import _DynP
+import numpy as np
+
+
+class DynP:
+    """Detects the change points in a time series.
+
+    Attributes
+    ----------
+    data: np.ndarray
+        Matrix storing the time series data.
+    num_bkps: int
+        Number of change points to detect.
+    jump: int
+        Interval for checking potential change points. Changing will
+        not provide optimal detection, but will reduce runtime.
+    min_size: int
+        Minimum size of a segment. Changing will not provide optimal
+        detection, but will reduce runtime.
+    """
+
+    def __init__(
+        self, data: np.ndarray, num_bkps: int, jump: int, min_size: int
+    ):
+        """Initialize the DynamicProgramming instance with given parameters."""
+        self.dynp = _DynP.DynamicProgramming(data, num_bkps, jump, min_size)
+
+    def set_num_threads(self, num_threads: int):
+        """Set the number of threads for parallelization.
+
+        Parameters
+        ----------
+        num_threads: int
+            The number of threads to use during computation. Default
+            is determined automatically.
+        """
+        self.dynp.set_threads(num_threads)
+
+    def fit(self, num_bkps: int) -> list:
+        """Calculate the cost matrix and return the breakpoints.
+
+        Parameters
+        ----------
+        num_bkps: int
+            number of change points to detect.
+
+        Returns
+        -------
+            list: A list of integers representing the breakpoints.
+        """
+        return self.dynp.fit()

src/CMakeLists.txt

@@ -0,0 +1,17 @@
+pybind11_add_module(_dupin dupininterface.cpp
+    dupin.h dupin.cpp
+)
+
+set_target_properties(_dupin PROPERTIES
+    CXX_STANDARD 17
+    CMAKE_CXX_STANDARD_REQUIRED True
+)
+
+target_include_directories(_dupin PRIVATE
+    ${EIGEN3_INCLUDE_DIR}
+    ${TBB_INCLUDE_DIRS}
+)
+
+target_link_libraries(_dupin PRIVATE TBB::tbb)
+target_compile_definitions(_dupin PRIVATE VERSION_INFO=${PROJECT_VERSION})
+target_compile_options(_dupin PRIVATE -O2 -march=native)

src/dupin.cpp

@@ -0,0 +1,165 @@
+#include <iostream>
+#include <iomanip>
+#include <limits>
+#include <unordered_map>
+#include <vector>
+#include <Eigen/Dense>
+#include <tbb/blocked_range2d.h>
+#include <tbb/global_control.h>
+#include <tbb/parallel_for.h>
+#include "dupin.h"
+
+using namespace std;
+using namespace Eigen;
+
+DynamicProgramming::DynamicProgramming()
+    : num_bkps(1), num_parameters(0), num_timesteps(0), jump(1), min_size(3) {}
+
+DynamicProgramming::DynamicProgramming(const Eigen::MatrixXd &data, int num_bkps_,
+                                        int jump_, int min_size_)
+                                        : data(data), num_bkps(num_bkps_),
+                                        jump(jump_), min_size(min_size_) {
+  num_timesteps = data.rows();
+  num_parameters = data.cols();
+}
+
+void DynamicProgramming::scale_data() {
+  Eigen::VectorXd min_val = data.colwise().minCoeff();
+  Eigen::VectorXd max_val = data.colwise().maxCoeff();
+  Eigen::VectorXd range = max_val - min_val;
+
+  for (int j = 0; j < num_parameters; ++j) {
+    if (range(j) == 0.0) {
+      data.col(j).setZero();
+    } else {
+      data.col(j) = (data.col(j).array() - min_val(j)) / range(j);
+    }
+  }
+}
+void DynamicProgramming::regression_setup(linear_fit_struct &lfit) {
+  lfit.x = Eigen::VectorXd::LinSpaced(num_timesteps, 0, num_timesteps - 1) /
+           (num_timesteps - 1);
+  lfit.y = data;
+}
+
+Eigen::VectorXd DynamicProgramming::regression_line(int start, int end, int dim,
+                                             linear_fit_struct &lfit) {
+  int n = end - start;
+  Eigen::VectorXd x = lfit.x.segment(start, n);
+  Eigen::VectorXd y = lfit.y.col(dim).segment(start, n);
+
+  double x_mean = x.mean();
+  double y_mean = y.mean();
+
+  Eigen::VectorXd x_centered = x.array() - x_mean;
+  Eigen::VectorXd y_centered = y.array() - y_mean;
+
+  double slope = x_centered.dot(y_centered) / x_centered.squaredNorm();
+  double intercept = y_mean - slope * x_mean;
+
+  return x.unaryExpr(
+      [slope, intercept](double xi) { return slope * xi + intercept; });
+}
+
+double DynamicProgramming::l2_cost(Eigen::MatrixXd &predicted_y, int start, int end) {
+  Eigen::MatrixXd diff = predicted_y.block(start, 0, end - start, num_parameters) -
+                  data.block(start, 0, end - start, num_parameters);
+  return std::sqrt(diff.array().square().sum());
+}
+
+Eigen::MatrixXd DynamicProgramming::predicted(int start, int end,
+                                       linear_fit_struct &lfit) {
+  Eigen::MatrixXd predicted_y(num_timesteps, num_parameters);
+  for (int i = 0; i < num_parameters; ++i) {
+    predicted_y.block(start, i, end - start, 1) =
+        regression_line(start, end, i, lfit);
+  }
+  return predicted_y;
+}
+
+double DynamicProgramming::cost_function(int start, int end) {
+  linear_fit_struct lfit;
+  regression_setup(lfit);
+  Eigen::MatrixXd predicted_y = predicted(start, end, lfit);
+  return l2_cost(predicted_y, start, end);
+}
+
+void DynamicProgramming::initialize_cost_matrix() {
+  scale_data();
+  cost_matrix.initialize(num_timesteps);
+  tbb::parallel_for(tbb::blocked_range<int>(0, num_timesteps),
+                    [&](const tbb::blocked_range<int> &r) {
+                      for (int i = r.begin(); i < r.end(); ++i) {
+                        for (int j = i + min_size; j < num_timesteps; ++j) {
+                          cost_matrix(i, j) = cost_function(i, j);
+                        }
+                      }
+                    });
+  cost_computed = true;
+}
+
+std::pair<double, std::vector<int>> DynamicProgramming::seg(int start, int end,
+                                                  int num_bkps) {
+  MemoKey key = {start, end, num_bkps};
+  auto it = memo.find(key);
+  if (it != memo.end()) {
+    return it->second;
+  }
+  if (num_bkps == 0) {
+    return {cost_matrix(start, end), {end}};
+  }
+
+  std::pair<double, std::vector<int>> best = {std::numeric_limits<double>::infinity(), {}};
+
+  for (int bkp = start + min_size; bkp < end; bkp++) {
+    if ((bkp - start) >= min_size && (end - bkp) >= min_size) {
+      auto left = seg(start, bkp, num_bkps - 1);
+      auto right = seg(bkp, end, 0);
+      double cost = left.first + right.first;
+      if (cost < best.first) {
+        best.first = cost;
+        best.second = left.second;
+        best.second.push_back(bkp);
+        best.second.insert(best.second.end(), right.second.begin(),
+                           right.second.end());
+      }
+    }
+  }
+
+  memo[key] = best;
+  return best;
+}
+
+std::vector<int> DynamicProgramming::compute_breakpoints() {
+  auto result = seg(0, num_timesteps - 1, num_bkps);
+  std::vector<int> breakpoints = result.second;
+  std::sort(breakpoints.begin(), breakpoints.end());
+  breakpoints.erase(std::unique(breakpoints.begin(), breakpoints.end()),
+                    breakpoints.end());
+  return breakpoints;
+}
+
+std::vector<int> DynamicProgramming::fit(int num_bkps_in){
+  num_bkps = num_bkps_in;
+  if (!cost_computed){
+  initialize_cost_matrix();
+  }
+  return compute_breakpoints();
+}
+
+void set_parallelization(int num_threads) {
+  static tbb::global_control gc(tbb::global_control::max_allowed_parallelism,
+                                num_threads);
+}
+
+DynamicProgramming::UpperTriangularMatrix &
+DynamicProgramming::getCostMatrix() {
+  return cost_matrix;
+}
+
+void DynamicProgramming::setCostMatrix(
+    const DynamicProgramming::UpperTriangularMatrix &value) {
+  cost_matrix = value;
+}
+
+int main() { return 0; }

src/dupin.h

@@ -0,0 +1,131 @@
+#pragma once
+
+#include <algorithm>
+#include <iostream>
+#include <limits>
+#include <unordered_map>
+#include <vector>
+#include <Eigen/Dense>
+
+
+// Calculates optimal breakpoints in time-series data using memoization
+class DynamicProgramming {
+private:
+
+  //stores upper triangular cost matrix efficiently
+  class UpperTriangularMatrix {
+  private:
+    std::vector<double> matrix;
+    std::vector<int> row_indices;
+    int length;
+
+    int index(int row, int col) const {
+        return row_indices[row] + col - row;
+    }
+
+  public:
+    UpperTriangularMatrix() : length(0) {}
+
+    void initialize(int n) {
+        length = n;
+        matrix.resize(n * (n + 1) / 2, 0.0);
+        row_indices.resize(n);
+        for (int row = 0; row < n; ++row) {
+            row_indices[row] = row * (2 * length - row + 1) / 2;
+        }
+    }
+
+    double &operator()(int row, int col) {
+        return matrix[index(row, col)];
+    }
+    int getSize() const { return length; }
+};
+  // Struct for memoization key, combining start, end, and number of
+  // breakpoints.
+  struct MemoKey {
+    int start;
+    int end;
+    int num_bkps;
+
+    // Comparison operator for MemoKey.
+    bool operator==(const MemoKey &other) const {
+      return start == other.start && end == other.end &&
+             num_bkps == other.num_bkps;
+    }
+  };
+
+  // Custom XOR-bit hash function for MemoKey, avoids clustering of data in
+  // unordered map to improve efficiency.
+  struct MemoKeyHash {
+    std::size_t operator()(const MemoKey &key) const {
+      return ((std::hash<int>()(key.start) ^
+               (std::hash<int>()(key.end) << 1)) >>
+              1) ^
+             std::hash<int>()(key.num_bkps);
+    }
+  };
+
+  // Memoization map to store the cost and partition for given parameters.
+  std::unordered_map<MemoKey, std::pair<double, std::vector<int>>, MemoKeyHash>
+      memo;
+
+  int num_bkps;          // Number of breakpoints to detect.
+  int num_parameters;    // Number of features in the dataset.
+  int num_timesteps;     // Number of data points (time steps).
+  int jump;              // Interval for checking potential breakpoints.
+  int min_size;          // Minimum size of a segment.
+
+  Eigen::MatrixXd data; // Matrix storing the dataset.
+  UpperTriangularMatrix cost_matrix; //Matrix storing costs
+  bool cost_computed = false;
+  // Structure for storing linear regression parameters.
+  struct linear_fit_struct {
+    Eigen::MatrixXd y; // Dependent variable (labels).
+    Eigen::VectorXd x; // z Independent variable (time steps).
+  };
+ // Scales the dataset using min-max normalization.
+  void scale_data();
+
+  // Prepares data for linear regression.
+  void regression_setup(linear_fit_struct &lfit);
+
+  // Calculates the regression line for a given data segment.
+  Eigen::VectorXd regression_line(int start, int end, int dim,
+                                  linear_fit_struct &lfit);
+
+  // Generates predicted values based on the linear regression model.
+  Eigen::MatrixXd predicted(int start, int end, linear_fit_struct &lfit);
+
+  // Calculates L2 cost (Euclidean distance) between predicted and actual data.
+  double l2_cost(Eigen::MatrixXd &predicted_y, int start, int end);
+
+  // Computes the cost of a specific data segment using linear regression.
+  double cost_function(int start, int end);
+
+    // Recursive function for dynamic programming segmentation.
+  std::pair<double, std::vector<int>> seg(int start, int end, int num_bkps);
+
+// Initializes and fills the cost matrix for all data segments.
+  void initialize_cost_matrix();
+
+  // Returns the optimal set of breakpoints after segmentation.
+  std::vector<int> compute_breakpoints();
+
+public:
+  // Default constructor.
+  DynamicProgramming();
+
+  // Parameterized constructor.
+  DynamicProgramming(const Eigen::MatrixXd &data, int num_bkps_, int jump_,
+                     int min_size_);
+
+  //Sets number of threads for parallelization
+  void set_parallelization(int num_threads);
+
+  // Calculates optimal breakpoints with given number of points.
+  std::vector<int> fit(int num_bkps_in);
+
+  // Getter functions for cost matrix.
+  DynamicProgramming::UpperTriangularMatrix &getCostMatrix();
+  void setCostMatrix(const DynamicProgramming::UpperTriangularMatrix &value);
+};

src/dupininterface.cpp

@@ -0,0 +1,15 @@
+#include "dupin.h"
+#include <pybind11/eigen.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(_DynP, m) {
+  py::class_<DynamicProgramming>(m, "DynamicProgramming")
+      .def(py::init<>())
+      .def_property("cost_matrix", &DynamicProgramming::getCostMatrix,
+                    &DynamicProgramming::setCostMatrix)
+      .def("fit", &DynamicProgramming::fit)
+      .def("set_threads", &DynamicProgramming::set_parallelization);
+}