espdev
diff --git a/‎src/csaps.cpp‎
Lines changed: 41 additions & 0 deletions b/‎src/csaps.cpp‎
Lines changed: 41 additions & 0 deletions
diff --git a/‎src/csaps.h‎
Lines changed: 10 additions & 0 deletions b/‎src/csaps.h‎
Lines changed: 10 additions & 0 deletions
@@ -1,4 +1,5 @@
 #include <iostream>
+#include <cmath>
 
 #include "csaps.h"
 
@@ -86,6 +87,17 @@ void UnivariateCubicSmoothingSpline::MakeSpline()
 
 DoubleArray UnivariateCubicSmoothingSpline::Evaluate(const DoubleArray & xidata)
 {
+  const size_t pcount = m_xdata.size();
+
+  auto mesh = m_xdata.segment(1, pcount - 2);
+  DoubleArray edges(pcount);
+
+  edges(0) = -DoubleLimits::infinity();
+  edges.segment(1, pcount - 2) = mesh;
+  edges(pcount - 1) = DoubleLimits::infinity();
+
+  auto indexes = Digitize(xidata, edges);
+
   return DoubleArray();
 }
 
@@ -95,4 +107,33 @@ DoubleArray UnivariateCubicSmoothingSpline::Diff(const DoubleArray &vec)
   return vec.tail(n) - vec.head(n);
 }
 
+IndexArray UnivariateCubicSmoothingSpline::Digitize(const DoubleArray &arr, const DoubleArray &bins)
+{
+  // This code works if `arr` and `bins` are monotonically increasing
+
+  IndexArray indexes(arr.size());
+
+  // Greater or equal (a >= b)
+  auto ge = [](double a, double b)
+  {
+    return a > b || std::abs(a - b) < std::abs(std::min(a, b)) * 1.e-8;
+  };
+
+  Eigen::DenseIndex kstart = 1;
+
+  for (Eigen::DenseIndex i = 0; i < arr.size(); ++i) {
+    double val = arr(i);
+
+    for (Eigen::DenseIndex k = kstart; k < bins.size(); ++k) {
+      if (ge(val, bins(k - 1)) && val < bins(k)) {
+        indexes(i) = k;
+        kstart = k;
+        break;
+      }
+    }
+  }
+
+  return indexes;
+}
+
 } // namespace csaps
@@ -2,13 +2,18 @@
 #ifndef CSAPS_H
 #define CSAPS_H
 
+#include <limits>
+
 #include <Eigen/Dense>
 
 
 namespace csaps
 {
 
 typedef Eigen::ArrayXd DoubleArray;
+typedef Eigen::Array<Eigen::DenseIndex, Eigen::Dynamic, 1> IndexArray;
+
+typedef std::numeric_limits<double> DoubleLimits;
 
 
 class UnivariateCubicSmoothingSpline
@@ -27,8 +32,13 @@ class UnivariateCubicSmoothingSpline
 
   void MakeSpline();
   DoubleArray Evaluate(const DoubleArray &xidata);
+
+  //! Calculate the 1-th discrete difference
   static DoubleArray Diff(const DoubleArray &vec);
 
+  //! Return the indices of the bins to which each value in input array belongs
+  static IndexArray Digitize(const DoubleArray &arr, const DoubleArray &bins);
+
 protected:
   DoubleArray m_xdata;
   DoubleArray m_ydata;