diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
index ee83e888ce9e4d37ef1f22064ba1d723904db5cb..d5a3db04232e465a852816b066bb9b8195b86072 100644
--- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
+++ b/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
diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h
index e16254241311ebdc32cfed90c015cb09c7179c6d..5338582b1d17e4b7178c307649ecb1eded48aa76 100644
--- a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.h
+++ b/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_ */
diff --git a/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp b/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp
index ce7e525a1c267f40442fd841621c2be4666035fd..26438d514c2d67c855c47c6ed1d274a4ad5d7c47 100644
--- a/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp
+++ b/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());
}