clean the code and add documents

MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp

@@ -15,9 +15,8 @@
 #include <limits>
 #include <utility>
 
-#include "PiecewiseLinearInterpolation.h"
-
 #include "BaseLib/quicksort.h"
+#include "PiecewiseLinearInterpolation.h"
 
 namespace MathLib
 {
@@ -28,9 +27,11 @@ PiecewiseLinearInterpolation::PiecewiseLinearInterpolation(
     : _supp_pnts(std::move(supporting_points)),
       _values_at_supp_pnts(std::move(values_at_supp_pnts))
 {
-    if (!supp_pnts_sorted) {
-        BaseLib::quicksort<double, double>(_supp_pnts, static_cast<std::size_t> (0),
-                                           _supp_pnts.size(), _values_at_supp_pnts);
+    if (!supp_pnts_sorted)
+    {
+        BaseLib::quicksort<double, double>(
+            _supp_pnts, static_cast<std::size_t>(0), _supp_pnts.size(),
+            _values_at_supp_pnts);
     }
 }
 
@@ -38,62 +39,69 @@ double PiecewiseLinearInterpolation::getValue(double pnt_to_interpolate) const
 {
     // search interval that has the point inside
     std::size_t interval_idx(std::numeric_limits<std::size_t>::max());
-    if (pnt_to_interpolate <= _supp_pnts.front()) {
+    if (pnt_to_interpolate <= _supp_pnts.front())
+    {
         interval_idx = 0;
         return _values_at_supp_pnts[0];
-    } else {
-        if (_supp_pnts.back() <= pnt_to_interpolate) {
+    }
+    else
+    {
+        if (_supp_pnts.back() <= pnt_to_interpolate)
+        {
             interval_idx = _supp_pnts.size() - 2;
             return _values_at_supp_pnts[_supp_pnts.size() - 1];
-        } else {
-            auto const& it(std::lower_bound(_supp_pnts.begin(), _supp_pnts.end(), pnt_to_interpolate));
+        }
+        else
+        {
+            auto const& it(std::lower_bound(
+                _supp_pnts.begin(), _supp_pnts.end(), pnt_to_interpolate));
             interval_idx = std::distance(_supp_pnts.begin(), it) - 1;
         }
     }
-    
+
     // compute linear interpolation polynom: y = m * (x - support[i]) + value[i]
-    const double m((_values_at_supp_pnts[interval_idx + 1] - _values_at_supp_pnts[interval_idx])
-                    / (_supp_pnts[interval_idx + 1] - _supp_pnts[interval_idx]));
+    const double m((_values_at_supp_pnts[interval_idx + 1] -
+                    _values_at_supp_pnts[interval_idx]) /
+                   (_supp_pnts[interval_idx + 1] - _supp_pnts[interval_idx]));
 
-    return m * (pnt_to_interpolate - _supp_pnts[interval_idx]) + _values_at_supp_pnts[interval_idx];
+    return m * (pnt_to_interpolate - _supp_pnts[interval_idx]) +
+           _values_at_supp_pnts[interval_idx];
 }
 
-double PiecewiseLinearInterpolation::GetCurveDerivative(double pnt_to_interpolate) const
+double PiecewiseLinearInterpolation::GetDerivative(
+    double const pnt_to_interpolate) const
 {
-    std::size_t interval_max(std::numeric_limits<std::size_t>::max());
-    std::size_t interval_idx(std::numeric_limits<std::size_t>::max());
-    interval_max = _supp_pnts.size();
-    if (pnt_to_interpolate <= _supp_pnts.front()) {
-        pnt_to_interpolate = _supp_pnts.front();
-        interval_idx = 0;
-    }
-    else {
-        if (_supp_pnts.back() <= pnt_to_interpolate) {
-            pnt_to_interpolate = _supp_pnts.back();
-            interval_idx = _supp_pnts.size() - 2;
-        }
-        else {
-            auto const& it(std::lower_bound(_supp_pnts.begin(), _supp_pnts.end(), pnt_to_interpolate));
-            interval_idx = std::distance(_supp_pnts.begin(), it) - 1;
-        }
+    if (pnt_to_interpolate < _supp_pnts.front() ||
+        _supp_pnts.back() < pnt_to_interpolate)
+    {
+        return 0;
     }
 
-    //interval_idx = interval_max - 1 - interval_idx;
-    if (interval_idx > 1 && interval_idx < interval_max - 2) {
-        double s1 = (0.5*_values_at_supp_pnts[interval_idx] - 0.5*_values_at_supp_pnts[interval_idx + 2])
-            / (0.5*_supp_pnts[interval_idx] - 0.5*_supp_pnts[interval_idx + 2]);
-        double s2 = (0.5*_values_at_supp_pnts[interval_idx - 1] - 0.5*_values_at_supp_pnts[interval_idx + 1])
-            / (0.5*_supp_pnts[interval_idx - 1] - 0.5*_supp_pnts[interval_idx + 1]);
-        double w = (pnt_to_interpolate - _supp_pnts[interval_idx + 1]) / (_supp_pnts[interval_idx] - _supp_pnts[interval_idx + 1]);
-        return (1. - w) * s1 + w * s2;
+    auto const& it(std::lower_bound(_supp_pnts.begin(), _supp_pnts.end(),
+                                    pnt_to_interpolate));
+    std::size_t const interval_idx = std::distance(_supp_pnts.begin(), it) - 1;
+
+    // interval_idx = interval_max - 1 - interval_idx;
+    if (interval_idx > 1 && interval_idx < _supp_pnts.size() - 2)
+    {
+        double const slope_right =
+            (_values_at_supp_pnts[interval_idx] -
+             _values_at_supp_pnts[interval_idx + 2]) /
+            (_supp_pnts[interval_idx] - _supp_pnts[interval_idx + 2]);
+        double const slope_left =
+            (_values_at_supp_pnts[interval_idx - 1] -
+             _values_at_supp_pnts[interval_idx + 1]) /
+            (_supp_pnts[interval_idx - 1] - _supp_pnts[interval_idx + 1]);
+        double const w =
+            (pnt_to_interpolate - _supp_pnts[interval_idx + 1]) /
+            (_supp_pnts[interval_idx] - _supp_pnts[interval_idx + 1]);
+        return (1. - w) * slope_right + w * slope_left;
     }
-    else {
-        if (std::fabs(_supp_pnts[interval_idx] - _supp_pnts[interval_idx + 1])>DBL_MIN)
-            return (_values_at_supp_pnts[interval_idx] - _values_at_supp_pnts[interval_idx + 1])
-            / (_supp_pnts[interval_idx] - _supp_pnts[interval_idx + 1]);
-        else
-            return 1 / DBL_EPSILON;
+    else
+    {
+        return (_values_at_supp_pnts[interval_idx] -
+                _values_at_supp_pnts[interval_idx + 1]) /
+               (_supp_pnts[interval_idx] - _supp_pnts[interval_idx + 1]);
     }
-
 }
-} // end MathLib
+}  // end MathLib

MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h

@@ -20,7 +20,8 @@
 namespace MathLib
 {
 /**
- * This class implements a one dimensional piecewise linear interpolation algorithm.
+ * This class implements a one dimensional piecewise linear interpolation
+ * algorithm.
  */
 class PiecewiseLinearInterpolation final
 {
@@ -41,7 +42,8 @@ public:
      * time.
      * @param supporting_points vector of supporting points
      * @param values_at_supp_pnts vector of values at the supporting points
-     * @param supp_pnts_sorted false (default), if it is sure the supporting points are sorted
+     * @param supp_pnts_sorted false (default), if it is sure the supporting
+     * points are sorted
      * one can set the switch to true
      */
     PiecewiseLinearInterpolation(std::vector<double>&& supporting_points,
@@ -51,17 +53,35 @@ public:
     /**
      * \brief Calculates the interpolation value.
      * @param pnt_to_interpolate The point should be located within the range
-     * \f$[x_{\min}, x_{\max}]\f$, where \f$x_{\min} = \min_{1 \le j \le n} x_j\f$ and
-     * \f$x_{\max} = \max_{1 \le j \le n} x_j\f$. Points outside of this interval are
-     * extrapolated.
+     * \f$[x_{\min}, x_{\max}]\f$, where \f$x_{\min} = \min_{1 \le j \le n}
+     * x_j\f$ and
+     * \f$x_{\max} = \max_{1 \le j \le n} x_j\f$. Points outside of this
+     * interval are
+     * set to x_{\min} or x_{\max}.
      * @return The interpolated value.
      */
     double getValue(double pnt_to_interpolate) const;
-	double GetCurveDerivative(double pnt_to_interpolate) const;
+
+    /**
+     * \brief Calculates the derivative value.
+     * @param pnt_to_interpolate The point should be located within the range
+     * \f$[x_{\min}, x_{\max}]\f$, where \f$x_{\min} = \min_{1 \le j \le n}
+     * x_j\f$ and
+     * \f$x_{\max} = \max_{1 \le j \le n} x_j\f$. when points are located
+     * outside of this interval
+     * the derivative is set to 0
+     * \attention if the points is located between the first and second points
+     * (or last and second to last point), the derivative is calculated by
+     * simple linear interpolation
+     * otherwise, it is calculated by second order of interpolation with central
+     * difference
+     */
+    double GetDerivative(double pnt_to_interpolate) const;
+
 private:
     std::vector<double> _supp_pnts;
     std::vector<double> _values_at_supp_pnts;
 };
-} // end namespace MathLib
+}  // end namespace MathLib
 
 #endif /* PIECEWISELINEARINTERPOLATION_H_ */

Tests/MathLib/TestPiecewiseLinearInterpolation.cpp

@@ -23,11 +23,15 @@ TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolation)
 {
     const std::size_t size(1000);
     std::vector<double> supp_pnts, values;
-    for (std::size_t k(0); k<size; ++k) {
+    for (std::size_t k(0); k < size; ++k)
+    {
         supp_pnts.push_back(static_cast<double>(k));
-        if (k % 2 == 0) {
+        if (k % 2 == 0)
+        {
             values.push_back(0.0);
-        } else {
+        }
+        else
+        {
             values.push_back(1.0);
         }
     }
@@ -35,35 +39,52 @@ TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolation)
     MathLib::PiecewiseLinearInterpolation interpolation{std::move(supp_pnts),
                                                         std::move(values)};
     // Interpolation
-    for (std::size_t k(0); k<size-1; ++k) {
-        ASSERT_NEAR(0.5, interpolation.getValue(k+0.5), std::numeric_limits<double>::epsilon());
+    for (std::size_t k(0); k < size - 1; ++k)
+    {
+        ASSERT_NEAR(0.5, interpolation.getValue(k + 0.5),
+                    std::numeric_limits<double>::epsilon());
     }
 
-    for (std::size_t k(0); k<size-1; ++k) {
-        if (k % 2 == 0) {
-            ASSERT_NEAR(0.25, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
-            ASSERT_NEAR(0.75, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
-        } else {
-            ASSERT_NEAR(0.75, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
-            ASSERT_NEAR(0.25, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
+    for (std::size_t k(0); k < size - 1; ++k)
+    {
+        if (k % 2 == 0)
+        {
+            ASSERT_NEAR(0.25, interpolation.getValue(k + 0.25),
+                        std::numeric_limits<double>::epsilon());
+            ASSERT_NEAR(0.75, interpolation.getValue(k + 0.75),
+                        std::numeric_limits<double>::epsilon());
+        }
+        else
+        {
+            ASSERT_NEAR(0.75, interpolation.getValue(k + 0.25),
+                        std::numeric_limits<double>::epsilon());
+            ASSERT_NEAR(0.25, interpolation.getValue(k + 0.75),
+                        std::numeric_limits<double>::epsilon());
         }
     }
 
     // Extrapolation
-    ASSERT_NEAR(0.0, interpolation.getValue(-0.5), std::numeric_limits<double>::epsilon());
+    ASSERT_NEAR(0.0, interpolation.getValue(-0.5),
+                std::numeric_limits<double>::epsilon());
     // Extrapolation
-    ASSERT_NEAR(1.0, interpolation.getValue(size-0.5), std::numeric_limits<double>::epsilon());
+    ASSERT_NEAR(1.0, interpolation.getValue(size - 0.5),
+                std::numeric_limits<double>::epsilon());
 }
 
-TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolationSupportPntsInReverseOrder)
+TEST(MathLibInterpolationAlgorithms,
+     PiecewiseLinearInterpolationSupportPntsInReverseOrder)
 {
     const std::size_t size(1000);
     std::vector<double> supp_pnts, values;
-    for (std::size_t k(0); k<size; ++k) {
-        supp_pnts.push_back(static_cast<double>(size-1-k));
-        if (k % 2 == 0) {
+    for (std::size_t k(0); k < size; ++k)
+    {
+        supp_pnts.push_back(static_cast<double>(size - 1 - k));
+        if (k % 2 == 0)
+        {
             values.push_back(1.0);
-        } else {
+        }
+        else
+        {
             values.push_back(0.0);
         }
     }
@@ -71,45 +92,61 @@ TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolationSupportPntsInRe
     MathLib::PiecewiseLinearInterpolation interpolation{std::move(supp_pnts),
                                                         std::move(values)};
     // Interpolation
-    for (std::size_t k(0); k<size-1; ++k) {
-        ASSERT_NEAR(0.5, interpolation.getValue(k+0.5), std::numeric_limits<double>::epsilon());
+    for (std::size_t k(0); k < size - 1; ++k)
+    {
+        ASSERT_NEAR(0.5, interpolation.getValue(k + 0.5),
+                    std::numeric_limits<double>::epsilon());
     }
 
-    for (std::size_t k(0); k<size-1; ++k) {
-        if (k % 2 == 0) {
-            ASSERT_NEAR(0.25, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
-            ASSERT_NEAR(0.75, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
-        } else {
-            ASSERT_NEAR(0.75, interpolation.getValue(k+0.25), std::numeric_limits<double>::epsilon());
-            ASSERT_NEAR(0.25, interpolation.getValue(k+0.75), std::numeric_limits<double>::epsilon());
+    for (std::size_t k(0); k < size - 1; ++k)
+    {
+        if (k % 2 == 0)
+        {
+            ASSERT_NEAR(0.25, interpolation.getValue(k + 0.25),
+                        std::numeric_limits<double>::epsilon());
+            ASSERT_NEAR(0.75, interpolation.getValue(k + 0.75),
+                        std::numeric_limits<double>::epsilon());
+        }
+        else
+        {
+            ASSERT_NEAR(0.75, interpolation.getValue(k + 0.25),
+                        std::numeric_limits<double>::epsilon());
+            ASSERT_NEAR(0.25, interpolation.getValue(k + 0.75),
+                        std::numeric_limits<double>::epsilon());
         }
     }
 
     // Extrapolation
-    ASSERT_NEAR(0.0, interpolation.getValue(-0.5), std::numeric_limits<double>::epsilon());
+    ASSERT_NEAR(0.0, interpolation.getValue(-0.5),
+                std::numeric_limits<double>::epsilon());
     // Extrapolation
-    ASSERT_NEAR(1.0, interpolation.getValue(size-0.5), std::numeric_limits<double>::epsilon());
+    ASSERT_NEAR(1.0, interpolation.getValue(size - 0.5),
+                std::numeric_limits<double>::epsilon());
 }
 
 TEST(MathLibInterpolationAlgorithms, PiecewiseLinearInterpolationDerivative)
 {
-	const std::size_t size(1000);
-	std::vector<double> supp_pnts, values;
-	for (std::size_t k(0); k<size; ++k) {
-		supp_pnts.push_back(static_cast<double>(k));
-		values.push_back(k*k);
-	}
-
-	MathLib::PiecewiseLinearInterpolation interpolation{ std::move(supp_pnts),
-		std::move(values) };
-	// Interpolation
-	for (std::size_t k(0); k<size - 1; ++k) {
+    const std::size_t size(1000);
+    std::vector<double> supp_pnts, values;
+    for (std::size_t k(0); k < size; ++k)
+    {
+        supp_pnts.push_back(static_cast<double>(k));
+        values.push_back(k * k);
+    }
 
-		ASSERT_NEAR(1 + 2 * k, interpolation.GetCurveDerivative(k + 0.5), std::numeric_limits<double>::epsilon());
-	}
+    MathLib::PiecewiseLinearInterpolation interpolation{std::move(supp_pnts),
+                                                        std::move(values)};
+    // Interpolation
+    for (std::size_t k(0); k < size - 1; ++k)
+    {
+        ASSERT_NEAR(1 + 2 * k, interpolation.GetDerivative(k + 0.5),
+                    std::numeric_limits<double>::epsilon());
+    }
 
-	// Extrapolation
-	ASSERT_NEAR(1, interpolation.GetCurveDerivative(-1), std::numeric_limits<double>::epsilon());
-	// Extrapolation
-	ASSERT_NEAR(1997, interpolation.GetCurveDerivative(1001), std::numeric_limits<double>::epsilon());
+    // Extrapolation
+    ASSERT_NEAR(0, interpolation.GetDerivative(-1),
+                std::numeric_limits<double>::epsilon());
+    // Extrapolation
+    ASSERT_NEAR(0, interpolation.GetDerivative(1001),
+                std::numeric_limits<double>::epsilon());
 }