ALFI
Advanced Library for Function Interpolation
Loading...
Searching...
No Matches
poly.h
Go to the documentation of this file.
1#pragma once
2
3#include <iostream>
4#include <cmath>
5
6#include "config.h"
7#include "util/numeric.h"
8
9namespace alfi::poly {
10 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
11 std::enable_if_t<!traits::has_size<Number>::value, Number>
12 val(const Container<Number>& coeffs, const Number& x) {
13 Number result = 0;
14 for (const Number& c : coeffs) {
15 result = result * x + c;
16 }
17 return result;
18 }
19
20 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
21 std::enable_if_t<traits::has_size<Container<Number>>::value, Container<Number>>
22 val(const Container<Number>& coeffs, const Container<Number>& xx) {
23 Container<Number> result(xx.size(), 0);
24#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
25#pragma omp parallel for
26#endif
27 for (SizeT i = 0; i < xx.size(); ++i) {
28 for (const Number& c : coeffs) {
29 result[i] = result[i] * xx[i] + c;
30 }
31 }
32 return result;
33 }
34
35 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
36 Container<Number> lagrange(const Container<Number>& X, const Container<Number>& Y) {
37 if (X.size() != Y.size()) {
38 std::cerr << "Error in function " << __FUNCTION__
39 << ": Vectors X (of size " << X.size()
40 << ") and Y (of size " << Y.size()
41 << ") are not the same size. Returning {NAN}..." << std::endl;
42 return {NAN};
43 }
44
45 if (X.empty()) {
46 std::cerr << "Error in function " << __FUNCTION__
47 << ": Vectors X and Y are empty. Cannot interpolate. Returning {NAN}..." << std::endl;
48 return {NAN};
49 }
50
51 const auto N = X.size();
52
53 Container<Number> P(N);
54 std::fill(P.begin(), P.end(), 0);
55
56 Container<Number> l(N);
57
58 for (SizeT k = 0; k < N; ++k) {
59 l.resize(1);
60 l[0] = 1;
61 for (SizeT j = 0; j < N; ++j) {
62 if (j != k) {
63 // l = conv(l, [1/(X[k]-X[j]), -X[j]/(X[k]-X[j])]);
64 l.resize(l.size() + 1);
65 l[l.size()-1] = 0;
66 for (SizeT i = l.size() - 1; i > 0; --i) {
67 l[i] = (l[i] - X[j] * l[i-1]) / (X[k] - X[j]);
68 }
69 l[0] /= X[k] - X[j];
70 }
71 }
72 for (SizeT i = 0; i < l.size(); ++i) {
73 P[i] += Y[k] * l[i];
74 }
75 }
76
77 return P;
78 }
79
80 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
81 Container<Number> lagrange_vals(const Container<Number>& X, const Container<Number>& Y, const Container<Number>& xx) {
82 const auto nn = xx.size();
83
84 if (X.size() != Y.size()) {
85 std::cerr << "Error in function " << __FUNCTION__
86 << ": Vectors X (of size " << X.size()
87 << ") and Y (of size " << Y.size()
88 << ") are not the same size. Returning an array of NaNs..." << std::endl;
89 Container<Number> yy(nn);
90 std::fill(yy.begin(), yy.end(), NAN);
91 return yy;
92 }
93
94 if (X.empty()) {
95 std::cerr << "Error in function " << __FUNCTION__
96 << ": Vectors X and Y are empty. Cannot interpolate. Returning an array of NaNs..." << std::endl;
97 Container<Number> yy(nn);
98 std::fill(yy.begin(), yy.end(), NAN);
99 return yy;
100 }
101
102 const auto N = X.size();
103
104 Container<Number> yy(nn);
105
106#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
107#pragma omp parallel for
108#endif
109 for (SizeT k = 0; k < nn; ++k) {
110 yy[k] = 0;
111 for (SizeT i = 0; i < N; ++i) {
112 Number l = 1;
113 for (SizeT j = 0; j < N; ++j) {
114 if (i != j) {
115 l *= (xx[k] - X[j]) / (X[i] - X[j]);
116 }
117 }
118 yy[k] += Y[i] * l;
119 }
120 }
121
122 return yy;
123 }
124
125 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
126 Container<Number> imp_lagrange(const Container<Number>& X, const Container<Number>& Y) {
127 if (X.size() != Y.size()) {
128 std::cerr << "Error in function " << __FUNCTION__
129 << ": Vectors X (of size " << X.size()
130 << ") and Y (of size " << Y.size()
131 << ") are not the same size. Returning {NAN}..." << std::endl;
132 return {NAN};
133 }
134
135 if (X.empty()) {
136 std::cerr << "Error in function " << __FUNCTION__
137 << ": Vectors X and Y are empty. Cannot interpolate. Returning {NAN}..." << std::endl;
138 return {NAN};
139 }
140
141 const auto N = X.size();
142
143 Container<Number> w_rev(N);
144 std::fill(w_rev.begin(), w_rev.end(), 1);
145 for (SizeT i = 0; i < N; ++i) {
146 for (SizeT j = 0; j < N; ++j) {
147 if (i != j) {
148 w_rev[i] *= (X[i] - X[j]);
149 }
150 }
151 }
152
153 Container<Number> P(N);
154 std::fill(P.begin(), P.end(), 0);
155
156 Container<Number> l(N + 1);
157
158 for (SizeT k = 0; k < N; ++k) {
159 l.resize(1);
160 l[0] = 1;
161 for (SizeT j = 0; j < N; ++j) {
162 // l = conv(l, [1, -X[j]]);
163 l.resize(l.size() + 1);
164 l[l.size()-1] = 0;
165 for (SizeT i = l.size() - 1; i > 0; --i) {
166 l[i] -= X[j] * l[i-1];
167 }
168 }
169 }
170
171 for (SizeT k = 0; k < N; ++k) {
172 Container<Number> l_cur(N);
173 l_cur[0] = l[0];
174 for (SizeT i = 1; i < N; ++i) {
175 l_cur[i] = l[i] + l_cur[i-1] * X[k];
176 }
177 for (SizeT i = 0; i < N; ++i) {
178 P[i] += Y[k] * l_cur[i] / w_rev[k];
179 }
180 }
181
182 return P;
183 }
184
185 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
186 Container<Number> imp_lagrange_vals(
187 const Container<Number>& X,
188 const Container<Number>& Y,
189 const Container<Number>& xx,
190 const Number& epsilon = std::numeric_limits<Number>::epsilon()
191 ) {
192 const auto nn = xx.size();
193
194 if (X.size() != Y.size()) {
195 std::cerr << "Error in function " << __FUNCTION__
196 << ": Vectors X (of size " << X.size()
197 << ") and Y (of size " << Y.size()
198 << ") are not the same size. Returning an array of NaNs..." << std::endl;
199 Container<Number> yy(nn);
200 std::fill(yy.begin(), yy.end(), NAN);
201 return yy;
202 }
203
204 if (X.empty()) {
205 std::cerr << "Error in function " << __FUNCTION__
206 << ": Vectors X and Y are empty. Cannot interpolate. Returning an array of NaNs..." << std::endl;
207 Container<Number> yy(nn);
208 std::fill(yy.begin(), yy.end(), NAN);
209 return yy;
210 }
211
212 const auto N = X.size();
213
214 Container<Number> w_rev(N);
215 std::fill(w_rev.begin(), w_rev.end(), 1);
216 for (SizeT i = 0; i < N; ++i) {
217 for (SizeT j = 0; j < N; ++j) {
218 if (i != j) {
219 w_rev[i] *= (X[i] - X[j]);
220 }
221 }
222 }
223
224 Container<Number> yy(nn);
225
226#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
227#pragma omp parallel for
228#endif
229 for (SizeT k = 0; k < nn; ++k) {
230 Number l = 1;
231 for (SizeT i = 0; i < N; ++i) {
232 l *= (xx[k] - X[i]);
233 }
234 yy[k] = 0;
235 for (SizeT i = 0; i < N; ++i) {
236 if (util::numeric::are_equal(xx[k], X[i], epsilon)) {
237 yy[k] = Y[i];
238 break;
239 } else {
240 yy[k] += Y[i] * l / (w_rev[i] * (xx[k] - X[i]));
241 }
242 }
243 }
244
245 return yy;
246 }
247
248 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
249 Container<Number> newton(const Container<Number>& X, const Container<Number>& Y) {
250 if (X.size() != Y.size()) {
251 std::cerr << "Error in function " << __FUNCTION__
252 << ": Vectors X (of size " << X.size()
253 << ") and Y (of size " << Y.size()
254 << ") are not the same size. Returning {NAN}..." << std::endl;
255 return {NAN};
256 }
257
258 if (X.empty()) {
259 std::cerr << "Error in function " << __FUNCTION__
260 << ": Vectors X and Y are empty. Cannot interpolate. Returning {NAN}..." << std::endl;
261 return {NAN};
262 }
263
264 const auto N = X.size();
265
266 Container<Number> F = Y;
267
268 for (SizeT i = 1; i < N; ++i) {
269 for (SizeT j = N - 1; j >= i; --j) {
270 F[j] = (F[j] - F[j-1]) / (X[j] - X[j-i]);
271 }
272 }
273
274 Container<Number> P(N);
275 std::fill(P.begin(), P.end() - 1, 0);
276 P[P.size()-1] = F[0];
277
278 Container<Number> f(N);
279
280 for (SizeT i = 0; i < N - 1; ++i) {
281 f.resize(1);
282 f[0] = F[i+1];
283 for (SizeT j = 0; j <= i; ++j) {
284 // f = conv(f, [1, -X[j]]);
285 f.resize(f.size() + 1);
286 f[f.size()-1] = 0;
287 for (SizeT k = f.size() - 1; k > 0; --k) {
288 f[k] -= X[j] * f[k-1];
289 }
290 }
291 // P = conv(P, [0, 1]); P = P + f;
292 // P is buffered, so no need in conv
293 for (SizeT j = 0; j < f.size(); ++j) {
294 P[N-1-i-1+j] += f[j];
295 }
296 }
297
298 return P;
299 }
300
301 template <typename Number = DefaultNumber, template <typename, typename...> class Container = DefaultContainer>
302 Container<Number> newton_vals(const Container<Number>& X, const Container<Number>& Y, const Container<Number>& xx) {
303 const auto nn = xx.size();
304
305 if (X.size() != Y.size()) {
306 std::cerr << "Error in function " << __FUNCTION__
307 << ": Vectors X (of size " << X.size()
308 << ") and Y (of size " << Y.size()
309 << ") are not the same size. Returning an array of NaNs..." << std::endl;
310 Container<Number> yy(nn);
311 std::fill(yy.begin(), yy.end(), NAN);
312 return yy;
313 }
314
315 if (X.empty()) {
316 std::cerr << "Error in function " << __FUNCTION__
317 << ": Vectors X and Y are empty. Cannot interpolate. Returning an array of NaNs..." << std::endl;
318 Container<Number> yy(nn);
319 std::fill(yy.begin(), yy.end(), NAN);
320 return yy;
321 }
322
323 const auto N = X.size();
324
325 Container<Number> F = Y;
326
327 for (SizeT i = 1; i < N; ++i) {
328 for (SizeT j = N - 1; j >= i; --j) {
329 F[j] = (F[j] - F[j-1]) / (X[j] - X[j-i]);
330 }
331 }
332
333 Container<Number> yy(nn);
334
335#if defined(_OPENMP) && !defined(ALFI_DISABLE_OPENMP)
336#pragma omp parallel for
337#endif
338 for (SizeT k = 0; k < nn; ++k) {
339 yy[k] = F[N-1];
340 for (SizeT iter = 0; iter < N - 1; ++iter) {
341 const auto i = N - 2 - iter;
342 yy[k] = yy[k] * (xx[k] - X[i]) + F[i];
343 }
344 }
345
346 return yy;
347 }
348}
alfi::poly
Definition poly.h:9
alfi::poly::val
std::enable_if_t<!traits::has_size< Number >::value, Number > val(const Container< Number > &coeffs, const Number &x)
Definition poly.h:12
alfi::poly::newton
Container< Number > newton(const Container< Number > &X, const Container< Number > &Y)
Definition poly.h:249
alfi::poly::lagrange_vals
Container< Number > lagrange_vals(const Container< Number > &X, const Container< Number > &Y, const Container< Number > &xx)
Definition poly.h:81
alfi::poly::newton_vals
Container< Number > newton_vals(const Container< Number > &X, const Container< Number > &Y, const Container< Number > &xx)
Definition poly.h:302
alfi::poly::lagrange
Container< Number > lagrange(const Container< Number > &X, const Container< Number > &Y)
Definition poly.h:36
alfi::poly::imp_lagrange
Container< Number > imp_lagrange(const Container< Number > &X, const Container< Number > &Y)
Definition poly.h:126
alfi::poly::imp_lagrange_vals
Container< Number > imp_lagrange_vals(const Container< Number > &X, const Container< Number > &Y, const Container< Number > &xx, const Number &epsilon=std::numeric_limits< Number >::epsilon())
Definition poly.h:186
alfi::util::numeric::are_equal
bool are_equal(const Number &a, const Number &b, const Number &epsilon=std::numeric_limits< Number >::epsilon())
Definition numeric.h:9
alfi::DefaultNumber
ALFI_DEFAULT_NUMBER DefaultNumber
Definition config.h:17
alfi::SizeT
ALFI_SIZE_TYPE SizeT
Definition config.h:22
alfi::DefaultContainer
ALFI_DEFAULT_CONTAINER< Number > DefaultContainer
Definition config.h:20