diff --git a/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp b/MathLib/InterpolationAlgorithms/PiecewiseLinearInterpolation.cpp
index 79c8c0cd939bbbbbb8d5f605e171644f6b269a7c..a314121459a65664806848c7e3356575683985ce 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,31 +27,74 @@ 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);
}
}
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()) {
- interval_idx = 0;
- } else {
- if (_supp_pnts.back() <= pnt_to_interpolate) {
- 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())
+ {
+ return _values_at_supp_pnts[0];
}
+ if (_supp_pnts.back() <= pnt_to_interpolate)
+ {
+ return _values_at_supp_pnts[_supp_pnts.size() - 1];
+ }
+
+ 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;
+
// 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];
+}
+
+double PiecewiseLinearInterpolation::getDerivative(
+ double const pnt_to_interpolate) const
+{
+ if (pnt_to_interpolate < _supp_pnts.front() ||
+ _supp_pnts.back() < pnt_to_interpolate)
+ {
+ return 0;
+ }
- return m * (pnt_to_interpolate - _supp_pnts[interval_idx]) + _values_at_supp_pnts[interval_idx];
+ 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 tangent_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 tangent_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) * tangent_right + w * tangent_left;
+ }
+ 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 f0bdbb08568e9ce8dc9c00959b1e6c4d826c6b08..7a38a80a5af6e9de664cbf08ddcfcd792d480147 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,33 @@ 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;
+ /**
+ * \brief Calculates derivative using quadratic interpolation
+ * and central difference quotient.
+ * @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 a point is located between the first and second points
+ * (or last and second to last point), the derivative is calculated
+ * using linear interpolation.
+ */
+ double getDerivative(double const 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/GeoLib/TestLineSegmentIntersect2d.cpp b/Tests/GeoLib/TestLineSegmentIntersect2d.cpp
index 819ffc4fa097f8dba52bd9250d998f33a8e76bb7..62d7a5d3c0cc8505ee86f0e9687a9ac9352fd0b7 100644
--- a/Tests/GeoLib/TestLineSegmentIntersect2d.cpp
+++ b/Tests/GeoLib/TestLineSegmentIntersect2d.cpp
@@ -50,7 +50,7 @@ public:
ac::gtest_reporter gtest_reporter;
};
-#if !defined(_MSC_VER) || (_MSC_VER >= 1900)
+#if !defined(_MSC_VER) || (_MSC_VER >= 2000)
// Compilers of MVS below 2015 do not support unrestricted unions. The
// unrestricted union is used by autocheck to handle test data. The autocheck
// workaround for MVS compilers (below version 2015) contains a bug and in the
diff --git a/Tests/GeoLib/TestSortSegments.cpp b/Tests/GeoLib/TestSortSegments.cpp
index 247dab95dfc7a731e6715faa9f9e37bf19162aa8..1b2c6c76ef8ea55d1510338e865c7589f64d28f1 100644
--- a/Tests/GeoLib/TestSortSegments.cpp
+++ b/Tests/GeoLib/TestSortSegments.cpp
@@ -8,7 +8,7 @@
*/
#include <gtest/gtest.h>
-
+#include <numeric>
#include "Tests/GeoLib/AutoCheckGenerators.h"
#include "GeoLib/AnalyticalGeometry.h"
@@ -30,7 +30,7 @@ public:
ac::gtest_reporter gtest_reporter;
};
-#if !defined(_MSC_VER) || (_MSC_VER >= 1900)
+#if !defined(_MSC_VER) || (_MSC_VER >= 2000)
// Compilers of MVS below 2015 do not support unrestricted unions. The
// unrestricted union is used by autocheck to handle test data. The autocheck
// workaround for MVS compilers (below version 2015) contains a bug and in the
diff --git a/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp b/Tests/MathLib/TestPiecewiseLinearInterpolation.cpp
index 149d7338fbd7837d01db101c947acfba004e7b09..e29fc0e0fc720aa106ace9f2b6da64f69ee87e0b 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.5, 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.5, 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,22 +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());
+ // Extrapolation
+ 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)
+ {
+ ASSERT_NEAR(1 + 2 * k, interpolation.getDerivative(k + 0.5),
+ std::numeric_limits<double>::epsilon());
}
// Extrapolation
- ASSERT_NEAR(-0.5, interpolation.getValue(-0.5), std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0, interpolation.getDerivative(-1),
+ std::numeric_limits<double>::epsilon());
// Extrapolation
- ASSERT_NEAR(1.5, interpolation.getValue(size-0.5), std::numeric_limits<double>::epsilon());
+ ASSERT_NEAR(0, interpolation.getDerivative(1001),
+ std::numeric_limits<double>::epsilon());
}