OGS
PiecewiseLinearInterpolation.cpp
Go to the documentation of this file.
1
15
16#include <algorithm>
17#include <cmath>
18#include <limits>
19#include <range/v3/range/conversion.hpp>
20#include <utility>
21
22#include "BaseLib/Error.h"
23#include "BaseLib/quicksort.h"
24
25namespace MathLib
26{
28 std::vector<double> supporting_points,
29 std::vector<double>
30 values_at_supp_pnts,
31 bool supp_pnts_sorted)
32 : supp_pnts_(std::move(supporting_points)),
33 values_at_supp_pnts_(std::move(values_at_supp_pnts))
34{
35 if (!supp_pnts_sorted)
36 {
38 supp_pnts_, static_cast<std::size_t>(0), supp_pnts_.size(),
40 }
41
42 const auto it = std::adjacent_find(supp_pnts_.begin(), supp_pnts_.end());
43 if (it != supp_pnts_.end())
44 {
45 const std::size_t i = std::distance(supp_pnts_.begin(), it);
47 "Variable {:d} and variable {:d} are the same. Piecewise linear "
48 "interpolation is not possible\n",
49 i, i + 1);
50 }
51}
52
54 std::vector<int> const& supporting_points,
55 std::vector<double> const& values_at_supp_pnts,
56 bool supp_pnts_sorted)
57 : supp_pnts_(supporting_points | ranges::to<std::vector<double>>()),
58 values_at_supp_pnts_(values_at_supp_pnts)
59{
60 if (!supp_pnts_sorted)
61 {
63 supp_pnts_, static_cast<std::size_t>(0), supp_pnts_.size(),
65 }
66
67 const auto it = std::adjacent_find(supp_pnts_.begin(), supp_pnts_.end());
68 if (it != supp_pnts_.end())
69 {
70 const std::size_t i = std::distance(supp_pnts_.begin(), it);
72 "Variable {:d} and variable {:d} are the same. Piecewise linear "
73 "interpolation is not possible\n",
74 i, i + 1);
75 }
76}
77
78double PiecewiseLinearInterpolation::getValue(double pnt_to_interpolate) const
79{
80 // search interval that has the point inside
81 if (pnt_to_interpolate <= supp_pnts_.front())
82 {
83 return values_at_supp_pnts_[0];
84 }
85
86 if (supp_pnts_.back() <= pnt_to_interpolate)
87 {
88 return values_at_supp_pnts_[supp_pnts_.size() - 1];
89 }
90
91 auto const& it(std::lower_bound(supp_pnts_.begin(), supp_pnts_.end(),
92 pnt_to_interpolate));
93 std::size_t const interval_idx = std::distance(supp_pnts_.begin(), it) - 1;
94
95 // support points.
96 double const x = supp_pnts_[interval_idx];
97 double const x_r = supp_pnts_[interval_idx + 1];
98
99 // values.
100 double const f = values_at_supp_pnts_[interval_idx];
101 double const f_r = values_at_supp_pnts_[interval_idx + 1];
102
103 // compute linear interpolation polynom: y = m * (x - support[i]) + value[i]
104 const double t = (pnt_to_interpolate - x) / (x_r - x);
105 return std::lerp(f, f_r, t);
106}
107
109 double const pnt_to_interpolate) const
110{
111 if (pnt_to_interpolate < supp_pnts_.front() ||
112 supp_pnts_.back() < pnt_to_interpolate)
113 {
114 return 0;
115 }
116
117 auto const& it(std::lower_bound(supp_pnts_.begin(), supp_pnts_.end(),
118 pnt_to_interpolate));
119 std::size_t interval_idx = std::distance(supp_pnts_.begin(), it);
120
121 if (pnt_to_interpolate == supp_pnts_.front())
122 {
123 interval_idx = 1;
124 }
125
126 if (interval_idx > 1 && interval_idx < supp_pnts_.size() - 2)
127 {
128 // left and right support points.
129 double const x_ll = supp_pnts_[interval_idx - 2];
130 double const x_l = supp_pnts_[interval_idx - 1];
131 double const x = supp_pnts_[interval_idx];
132 double const x_r = supp_pnts_[interval_idx + 1];
133
134 // left and right values.
135 double const f_ll = values_at_supp_pnts_[interval_idx - 2];
136 double const f_l = values_at_supp_pnts_[interval_idx - 1];
137 double const f = values_at_supp_pnts_[interval_idx];
138 double const f_r = values_at_supp_pnts_[interval_idx + 1];
139
140 double const tangent_right = (f_l - f_r) / (x_l - x_r);
141 double const tangent_left = (f_ll - f) / (x_ll - x);
142 double const w = (pnt_to_interpolate - x) / (x_l - x);
143 return (1. - w) * tangent_right + w * tangent_left;
144 }
145
146 return (values_at_supp_pnts_[interval_idx] -
147 values_at_supp_pnts_[interval_idx - 1]) /
148 (supp_pnts_[interval_idx] - supp_pnts_[interval_idx - 1]);
149}
150
152{
153 assert(!supp_pnts_.empty());
154 return supp_pnts_.back();
155}
157{
158 assert(!supp_pnts_.empty());
159 return supp_pnts_.front();
160}
161} // namespace MathLib
#define OGS_FATAL(...)
Definition Error.h:26
Definition of the PiecewiseLinearInterpolation class.
PiecewiseLinearInterpolation(std::vector< double > supporting_points, std::vector< double > values_at_supp_pnts, bool supp_pnts_sorted=false)
double getDerivative(double const pnt_to_interpolate) const
Calculates derivative using quadratic interpolation and central difference quotient.
double getValue(double pnt_to_interpolate) const
Calculates the interpolation value.
void quicksort(It1 first1, It1 last1, It2 first2, Comparator compare)
Definition quicksort.h:26
static const double t
Definition of the quicksort function.