ARTS 2.5.0 (git: 9ee3ac6c)
interpolation_lagrange.h
Go to the documentation of this file.
1#ifndef interpolation_lagrange_h
2#define interpolation_lagrange_h
3
4#include "constants.h"
5#include "debug.h"
6#include "enums.h"
7#include "grids.h"
8#include "matpackVII.h"
9#include "nonstd.h"
10#include <algorithm>
11#include <array>
12#include <functional>
13#include <limits>
14#include <memory>
15#include <numeric>
16#include <type_traits>
17#include <vector>
18
19namespace Interpolation {
26constexpr Index cycler(const Index n, const Index N) noexcept { return n >= N ? n - N : n; }
27
39constexpr Numeric cyclic_clamp(const Numeric x,
40 const std::pair<Numeric, Numeric> xlim) noexcept {
41 if (x < xlim.first)
42 return cyclic_clamp(x + xlim.second - xlim.first, xlim);
43 if (x >= xlim.second)
44 return cyclic_clamp(x - xlim.second + xlim.first, xlim);
45 return x;
46}
47
54constexpr Numeric min_cyclic(const Numeric x,
55 const std::pair<Numeric, Numeric> xlim) noexcept {
56 const Numeric dx = xlim.second - xlim.first;
57 const bool lo = nonstd::abs(x) < nonstd::abs(x - dx);
58 const bool hi = nonstd::abs(x) < nonstd::abs(x + dx);
59 const bool me = nonstd::abs(x + dx) < nonstd::abs(x - dx);
60 return (lo and hi) ? x : me ? x + dx : x - dx;
61}
62
70constexpr bool full_cycle(const Numeric xo,
71 const Numeric xn,
72 const std::pair<Numeric, Numeric> xlim) noexcept {
73 return (xo == xlim.first and xn == xlim.second) or
74 (xo == xlim.second and xn == xlim.first);
75}
76
84template <class SortedVectorType>
85constexpr Index start_pos_finder(const Numeric x, const SortedVectorType& xvec) noexcept {
86 if (const Index n = xvec.size(); n > 1) {
87 const Numeric x0 = xvec[ 0];
88 const Numeric x1 = xvec[n - 1];
89 const Numeric frac = (x - x0) / (x1 - x0);
90 const auto start_pos = Index(frac * (Numeric)(n - 2));
91 return start_pos > 0 ? (start_pos < n ? start_pos : n - 1) : 0;
92 }
93 return 0;
94}
95
97
113constexpr Index IMAX(const Index a, const Index b) noexcept { return a > b ? a : b; }
114
116
132constexpr Index IMIN(const Index a, const Index b) noexcept { return a < b ? a : b; }
133
148template <class SortedVectorType>
149constexpr Index pos_finder(const Index pos0, const Numeric x, const SortedVectorType& xi,
150 const Index polyorder, const bool cyclic,
151 const bool ascending,
152 const std::pair<Numeric, Numeric> cycle = {
153 -180, 180}) noexcept {
154 if (cyclic and (x < cycle.first or x > cycle.second))
155 // We are below or above the cycle so we must redo calculations after clamping x to the cycle
156 return pos_finder(pos0, cyclic_clamp(x, cycle), xi, polyorder, cyclic, ascending, cycle);
157 const Index N = xi.size()-1;
158 Index p0=pos0;
159
160 // Loops to find the first position with a neighbor larger or smaller
161 if (ascending) {
162 while (p0 < N and xi[p0] < x) ++p0;
163 while (p0 > 0 and xi[p0] > x) --p0;
164 } else {
165 while (p0 < N and xi[p0] > x) ++p0;
166 while (p0 > 0 and xi[p0] < x) --p0;
167 }
168
169 // Adjustment for higher and lower polynominal orders so that x is in the middle
170 if (polyorder) {
171 if (cyclic) // Max N since we can overstep bounds
172 return IMIN(IMAX(p0 - polyorder / 2, 0), N);
173 // Max N-polyorder since we cannot overstep bounds
174 return IMIN(IMAX(p0 - polyorder / 2, 0), N-polyorder);
175 }
176 // In the nearest neighbor case, we mus choose the closest neighbor
177 if (p0 < N and (nonstd::abs(xi[p0] - x) >= nonstd::abs(xi[p0 + 1] - x)))
178 return p0 + 1;
179 return p0;
180}
181
193 Standard, /* 1-to-1 interpolation grid */
194 Cyclic, /* Cyclic interpolation grid */
195 Log, /* Natural logarithm interpolation grid */
196 Log10, /* 10-base logarithm interpolation grid */
197 Log2, /* 2-base logarithm interpolation grid */
198 SinDeg, /* Cosine in degrees interpolation grid, grid only defined [-90, 90] */
199 SinRad, /* Cosine in radians interpolation grid, grid only defined [-PI/2, PI/2] */
200 CosDeg, /* Cosine in degrees interpolation grid, grid only defined [0, 180] */
201 CosRad /* Cosine in radians interpolation grid, grid only defined [0, PI] */
202 );
203
213template <GridType type, class SortedVectorType>
214constexpr Numeric l(const Index p0, const Index n, const Numeric x,
215 const SortedVectorType& xi, const Index j,
216 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle = {
217 -180, 180}) noexcept {
218 Numeric val = 1.0;
219 for (Index m = 0; m < n; m++) {
220 if (m not_eq j) {
221 if constexpr (type == GridType::Log) {
222 // Natural log weights
223 val *= (std::log(x) - std::log(xi[m + p0])) /
224 (std::log(xi[j + p0]) - std::log(xi[m + p0]));
225 } else if constexpr (type == GridType::Log10) {
226 // Common log weights
227 val *= (std::log10(x) - std::log10(xi[m + p0])) /
228 (std::log10(xi[j + p0]) - std::log10(xi[m + p0]));
229 } else if constexpr (type == GridType::Log2) {
230 // Binary log weights
231 val *= (std::log2(x) - std::log2(xi[m + p0])) /
232 (std::log2(xi[j + p0]) - std::log2(xi[m + p0]));
233 } else if constexpr (type == GridType::SinDeg) {
234 // Sine in degrees weights
235 using Conversion::sind;
236 val *= (sind(x) - sind(xi[m + p0])) / (sind(xi[j + p0]) - sind(xi[m + p0]));
237 } else if constexpr (type == GridType::SinRad) {
238 // Sine in radians weights
239 using std::sin;
240 val *= (sin(x) - sin(xi[m + p0])) / (sin(xi[j + p0]) - sin(xi[m + p0]));
241 } else if constexpr (type == GridType::CosDeg) {
242 // Cosine in degrees weights (nb. order changed)
243 using Conversion::cosd;
244 val *= (cosd(xi[m + p0]) - cosd(x)) / (cosd(xi[m + p0]) - cosd(xi[j + p0]));
245 } else if constexpr (type == GridType::CosRad) {
246 // Cosine in radians weights (nb. order changed)
247 using std::cos;
248 val *= (cos(xi[m + p0]) - cos(x)) / (cos(xi[m + p0]) - cos(xi[j + p0]));
249 } else if constexpr (type == GridType::Standard) {
250 // Linear weights, simple and straightforward
251 val *= (x - xi[m + p0]) / (xi[j + p0] - xi[m + p0]);
252 } else if constexpr (type == GridType::Cyclic) {
253 // Cyclic weights
254 // We have to ensure that all weights are cyclic (e.g., 355 degrees < -6 degrees)
255 const Index N = Index(xi.size());
256 const Index m_pos = cycler(m + p0, N);
257 const Index j_pos = cycler(j + p0, N);
258 const Numeric x_val = cyclic_clamp(x, cycle);
259 if (full_cycle(xi[0], xi[N-1], cycle)) {
260 // We ignore the last point in full cycles
261 if (j_pos == N - 1)
262 return 0;
263 if (m_pos not_eq N - 1)
264 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
265 } else {
266 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
267 }
268 }
269 }
270 }
271 return val;
272}
273
288template <GridType type, class SortedVectorType,
289 class LagrangeVectorType>
290constexpr double dl_dval(
291 const Index p0, const Index n, const Numeric x, const SortedVectorType& xi,
292 [[maybe_unused]] const LagrangeVectorType& li, const Index j, const Index i,
293 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle) noexcept {
294 if constexpr (type == GridType::Standard) {
295 // Linear weights, simple and straightforward
296 if (x not_eq xi[i + p0]) {
297 // A simple case when x is not on the grid
298 return li[j] / (x - xi[i + p0]);
299 }
300 // We have to resort to the full recalculations
301 Numeric val = 1.0 / (xi[j + p0] - xi[i + p0]);
302 for (Index m = 0; m < n; m++) {
303 if (m not_eq j and m not_eq i) {
304 val *= (x - xi[m + p0]) / (xi[j + p0] - xi[m + p0]);
305 }
306 }
307 return val;
308 } else if constexpr (type == GridType::Cyclic) {
309 // Cyclic weights
310 // We have to ensure that all weights are cyclic (e.g., 355 degrees < -6 degrees)
311 const decltype(i + p0) N = xi.size();
312 const Index i_pos = cycler(i + p0, N);
313 const Index j_pos = cycler(j + p0, N);
314 const Numeric x_val = cyclic_clamp(x, cycle);
315 if (full_cycle(xi[0], xi[N-1], cycle)) {
316 // We ignore the last point in full cycles
317 if (i_pos == N - 1 or j_pos == N - 1)
318 return 0;
319 if (x_val not_eq xi[i_pos])
320 // A simple case when x is not on the grid
321 return li[j] / min_cyclic(x_val - xi[i_pos], cycle);
322
323 // We have to resort to the full recalculations
324 Numeric val = 1.0 / min_cyclic(xi[j_pos] - xi[i_pos], cycle);
325 for (Index m = 0; m < n; m++) {
326 if (m not_eq j and m not_eq i) {
327 const Index m_pos = cycler(m + p0, N);
328 if (m_pos not_eq N - 1) {
329 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
330 }
331 }
332 }
333 return val;
334 }
335 if (x_val not_eq xi[i_pos])
336 // A simple case when x is not on the grid
337 return li[j] / min_cyclic(x_val - xi[i_pos], cycle);
338 // We have to resort to the full recalculations
339 Numeric val = 1.0 / min_cyclic(xi[j_pos] - xi[i_pos], cycle);
340 for (Index m = 0; m < n; m++) {
341 if (m not_eq j and m not_eq i) {
342 const Index m_pos = cycler(m + p0, N);
343 val *= min_cyclic(x_val - xi[m_pos], cycle) / min_cyclic(xi[j_pos] - xi[m_pos], cycle);
344 }
345 }
346 return val;
347 } else /*if any other case */ {
348 // All other cases have to use full calculations because we need a linear derivative
349 Numeric val = 1.0 / (xi[j + p0] - xi[i + p0]);
350 for (Index m = 0; m < n; m++) {
351 if (m not_eq j and m not_eq i) {
352 val *= (x - xi[m + p0]) / (xi[j + p0] - xi[m + p0]);
353 }
354 }
355 return val;
356 }
357}
358
369template <GridType type, class SortedVectorType,
370 class LagrangeVectorType>
371constexpr Numeric dl(const Index p0, const Index n, const Numeric x,
372 const SortedVectorType& xi, const LagrangeVectorType& li,
373 const Index j,
374 const std::pair<Numeric, Numeric> cycle = {-180,
375 180}) noexcept {
376 Numeric dval = 0.0;
377 for (Index i = 0; i < n; i++) {
378 if (i not_eq j) {
379 dval += dl_dval<type>(p0, n, x, xi, li, j, i, cycle);
380 }
381 }
382 return dval;
383}
384
386template <class SortedVectorType>
387constexpr bool is_ascending(const SortedVectorType& xi) noexcept {
388 if (xi.nelem() > 1)
389 return xi[0] < xi[1];
390 return false;
391}
392
402template <class SortedVectorType>
403void check_lagrange_interpolation([[maybe_unused]] const SortedVectorType& xi,
404 [[maybe_unused]] const Index polyorder = 1,
405 [[maybe_unused]] const std::pair<Numeric, Numeric> x={std::numeric_limits<Numeric>::infinity(), -std::numeric_limits<Numeric>::infinity()},
406 [[maybe_unused]] const Numeric extrapol = 0.5,
407 [[maybe_unused]] const GridType type = GridType::Standard,
408 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle = {-180, 180}) {
409 const Index n = Index(xi.size());
410
411 ARTS_USER_ERROR_IF (polyorder >= n,"Interpolation setup has failed!\n"
412 "\tRequesting greater interpolation order than possible with given input grid")
413 ARTS_USER_ERROR_IF(type == GridType::Cyclic and cycle.first >= cycle.second,
414 "Interpolation setup has failed!\n" "\tBad cycle, must be [first, second)")
415 if (polyorder and extrapol > 0 and GridType::Cyclic not_eq type) {
416 const bool ascending = is_ascending(xi);
417 const Numeric xmin = ascending ? xi[0 ] - extrapol * nonstd::abs(xi[1 ] - xi[0 ]) :
418 xi[n-1] - extrapol * nonstd::abs(xi[n-2] - xi[n-1]) ;
419 const Numeric xmax = ascending ? xi[n-1] + extrapol * nonstd::abs(xi[n-2] - xi[n-1]) :
420 xi[0 ] + extrapol * nonstd::abs(xi[1 ] - xi[0 ]) ;
421 ARTS_USER_ERROR_IF(x.first < xmin or x.second > xmax,
422 "Interpolation setup has failed!\n"
423 "\tThe new grid has limits: ", x.first, ' ', x.second, '\n',
424 "\tThe old grid has limits: ", xmin, ' ', xmax)
425 }
426}
427
428
438template <class SortedVectorType>
439void check_lagrange_interpolation([[maybe_unused]] const SortedVectorType& xi,
440 [[maybe_unused]] const Index polyorder,
441 [[maybe_unused]] const Numeric x,
442 [[maybe_unused]] const Numeric extrapol = 0.5,
443 [[maybe_unused]] const GridType type = GridType::Standard,
444 [[maybe_unused]] const std::pair<Numeric, Numeric> cycle = {-180, 180}) {
445 check_lagrange_interpolation(xi, polyorder, {x, x}, extrapol, type, cycle);
446}
447
449struct Lagrange {
452
455
458
459 /* Number of weights */
460 [[nodiscard]] Index size() const noexcept { return lx.size(); }
461
463 Lagrange(Lagrange&& l) noexcept : pos(l.pos), lx(std::move(l.lx)), dlx(std::move(l.dlx)) {}
464
466 Lagrange(const Lagrange& l) noexcept : pos(l.pos), lx(l.lx), dlx(l.dlx) {}
467
470 pos = l.pos;
471 lx = std::move(l.lx);
472 dlx = std::move(l.dlx);
473 return *this;
474 }
475
477 Lagrange& operator=(const Lagrange& l) noexcept {
478 pos = l.pos;
479 lx = l.lx;
480 dlx = l.dlx;
481 return *this;
482 }
483
494 template <class SortedVectorType>
495 Lagrange(const Index pos0, const Numeric x, const SortedVectorType& xi,
496 const Index polyorder = 1,
497 const bool do_derivs = false,
498 const GridType type = GridType::Standard,
499 const std::pair<Numeric, Numeric> cycle = {-180, 180}) noexcept {
500 const Index n = xi.size();
501 if (n > 1) {
502 const Index p = polyorder + 1;
503 const bool ascending = is_ascending(xi);
504
505 pos = pos_finder(pos0, x, xi, polyorder, type == GridType::Cyclic, ascending, cycle);
506 lx = lx_finder(pos, p, x, xi, type, cycle);
507 if (do_derivs) dlx = dlx_finder(pos, p, x, xi, lx, type, cycle);
508 } else {
509 pos = 0;
510 lx = Array<Numeric>(1, 1);
511 dlx = Array<Numeric>(1, 0);
512 }
513 }
514
516 Lagrange() noexcept : lx(1, 1), dlx(1, 0) {}
517
519 friend std::ostream& operator<<(std::ostream& os, const Lagrange& l) {
520 os << "pos: " << l.pos << '\n' << "lx: ";
521 for (auto x : l.lx) os << ' ' << x;
522 os << '\n' << "dlx: ";
523 for (auto x : l.dlx) os << ' ' << x;
524 return os << '\n';
525 }
526
527 private:
537 template <class SortedVectorType>
539 const Index p0, const Index n, const Numeric x,
540 const SortedVectorType& xi, const GridType type,
541 const std::pair<Numeric, Numeric> cycle) noexcept {
542 Array<Numeric> out(n);
543 switch (type) {
544 case GridType::Standard:
545 for (Index j = 0; j < n; j++)
546 out[j] = l<GridType::Standard>(p0, n, x, xi, j);
547 break;
548 case GridType::Log:
549 for (Index j = 0; j < n; j++)
550 out[j] = l<GridType::Log>(p0, n, x, xi, j);
551 break;
552 case GridType::Log10:
553 for (Index j = 0; j < n; j++)
554 out[j] = l<GridType::Log10>(p0, n, x, xi, j);
555 break;
556 case GridType::Log2:
557 for (Index j = 0; j < n; j++)
558 out[j] = l<GridType::Log2>(p0, n, x, xi, j);
559 break;
560 case GridType::Cyclic:
561 for (Index j = 0; j < n; j++)
562 out[j] = l<GridType::Cyclic>(p0, n, x, xi, j, cycle);
563 break;
564 case GridType::SinDeg:
565 for (Index j = 0; j < n; j++)
566 out[j] = l<GridType::SinDeg>(p0, n, x, xi, j);
567 break;
568 case GridType::SinRad:
569 for (Index j = 0; j < n; j++)
570 out[j] = l<GridType::SinRad>(p0, n, x, xi, j);
571 break;
572 case GridType::CosDeg:
573 for (Index j = 0; j < n; j++)
574 out[j] = l<GridType::CosDeg>(p0, n, x, xi, j);
575 break;
576 case GridType::CosRad:
577 for (Index j = 0; j < n; j++)
578 out[j] = l<GridType::CosRad>(p0, n, x, xi, j);
579 break;
580 case GridType::FINAL: { /* pass */
581 }
582 }
583 return out;
584 }
585
596 template <class SortedVectorType>
598 const Index p0, const Index n, const Numeric x,
599 const SortedVectorType& xi, const Array<Numeric>& li,
600 const GridType type,
601 const std::pair<Numeric, Numeric> cycle) noexcept {
602 Array<Numeric> out(n);
603 switch (type) {
604 case GridType::Standard:
605 for (Index j = 0; j < n; j++)
606 out[j] = dl<GridType::Standard>(p0, n, x, xi, li, j);
607 break;
608 case GridType::Log:
609 for (Index j = 0; j < n; j++)
610 out[j] = dl<GridType::Log>(p0, n, x, xi, li, j);
611 break;
612 case GridType::Log10:
613 for (Index j = 0; j < n; j++)
614 out[j] = dl<GridType::Log10>(p0, n, x, xi, li, j);
615 break;
616 case GridType::Log2:
617 for (Index j = 0; j < n; j++)
618 out[j] = dl<GridType::Log2>(p0, n, x, xi, li, j);
619 break;
620 case GridType::Cyclic:
621 for (Index j = 0; j < n; j++)
622 out[j] = dl<GridType::Cyclic>(p0, n, x, xi, li, j, cycle);
623 break;
624 case GridType::SinDeg:
625 for (Index j = 0; j < n; j++)
626 out[j] = dl<GridType::SinDeg>(p0, n, x, xi, li, j);
627 break;
628 case GridType::SinRad:
629 for (Index j = 0; j < n; j++)
630 out[j] = dl<GridType::SinRad>(p0, n, x, xi, li, j);
631 break;
632 case GridType::CosDeg:
633 for (Index j = 0; j < n; j++)
634 out[j] = dl<GridType::CosDeg>(p0, n, x, xi, li, j);
635 break;
636 case GridType::CosRad:
637 for (Index j = 0; j < n; j++)
638 out[j] = dl<GridType::CosRad>(p0, n, x, xi, li, j);
639 break;
640 case GridType::FINAL: { /* pass */
641 }
642 }
643 return out;
644 }
645}; // Lagrange
646
648template <std::size_t PolyOrder>
652
654 std::array<Numeric, PolyOrder + 1> lx;
655
657 std::array<Numeric, PolyOrder + 1> dlx;
658
659 /* Number of weights */
660 static constexpr Index size() noexcept { return PolyOrder + 1; }
661
663 constexpr FixedLagrange(FixedLagrange&& l) noexcept : pos(l.pos), lx(std::move(l.lx)), dlx(std::move(l.dlx)) {}
664
666 constexpr FixedLagrange& operator=(FixedLagrange&& l) noexcept {
667 pos = l.pos;
668 lx = std::move(l.lx);
669 dlx = std::move(l.dlx);
670 return *this;
671 }
672
682 template <class SortedVectorType>
683 constexpr FixedLagrange(const Index p0, const Numeric x,
684 const SortedVectorType& xi,
685 const bool do_derivs = false,
686 const GridType type = GridType::Standard,
687 const std::pair<Numeric, Numeric> cycle = {-180, 180}) noexcept
688 : pos(pos_finder(p0, x, xi, PolyOrder,
689 type == GridType::Cyclic,
690 xi.size() > 1 ? xi[0] < xi[1] : false, cycle)),
691 lx(lx_finder(pos, x, xi, type, cycle)),
692 dlx(do_derivs ? dlx_finder(pos, x, xi, lx, type, cycle)
693 : std::array<Numeric, PolyOrder + 1>{}) {}
694
696 friend std::ostream& operator<<(std::ostream& os, const FixedLagrange& l) {
697 os << "pos: " << l.pos << '\n' << "lx: ";
698 for (auto x : l.lx) os << ' ' << x;
699 os << '\n' << "dlx: ";
700 for (auto x : l.dlx) os << ' ' << x;
701 return os << '\n';
702 }
703
704 private:
713 template <class SortedVectorType>
714 static constexpr std::array<Numeric, PolyOrder + 1> lx_finder(
715 const Index p0, const Numeric x, const SortedVectorType& xi,
716 const GridType type,
717 const std::pair<Numeric, Numeric> cycle) noexcept {
718 std::array<Numeric, PolyOrder + 1> out{};
719 constexpr Index n = PolyOrder + 1;
720 switch (type) {
721 case GridType::Standard:
722 for (Index j = 0; j < n; j++)
723 out[j] = l<GridType::Standard>(p0, n, x, xi, j);
724 break;
725 case GridType::Log:
726 for (Index j = 0; j < n; j++)
727 out[j] = l<GridType::Log>(p0, n, x, xi, j);
728 break;
729 case GridType::Log10:
730 for (Index j = 0; j < n; j++)
731 out[j] = l<GridType::Log10>(p0, n, x, xi, j);
732 break;
733 case GridType::Log2:
734 for (Index j = 0; j < n; j++)
735 out[j] = l<GridType::Log2>(p0, n, x, xi, j);
736 break;
737 case GridType::Cyclic:
738 for (Index j = 0; j < n; j++)
739 out[j] = l<GridType::Cyclic>(p0, n, x, xi, j, cycle);
740 break;
741 case GridType::SinDeg:
742 for (Index j = 0; j < n; j++)
743 out[j] = l<GridType::SinDeg>(p0, n, x, xi, j);
744 break;
745 case GridType::SinRad:
746 for (Index j = 0; j < n; j++)
747 out[j] = l<GridType::SinRad>(p0, n, x, xi, j);
748 break;
749 case GridType::CosDeg:
750 for (Index j = 0; j < n; j++)
751 out[j] = l<GridType::CosDeg>(p0, n, x, xi, j);
752 break;
753 case GridType::CosRad:
754 for (Index j = 0; j < n; j++)
755 out[j] = l<GridType::CosRad>(p0, n, x, xi, j);
756 break;
757 case GridType::FINAL: { /* pass */
758 }
759 }
760 return out;
761 }
762
772 template <class SortedVectorType>
773 static constexpr std::array<Numeric, PolyOrder + 1> dlx_finder(
774 const Index p0, const Numeric x, const SortedVectorType& xi,
775 const std::array<Numeric, PolyOrder + 1>& li, const GridType type,
776 const std::pair<Numeric, Numeric> cycle) noexcept {
777 std::array<Numeric, PolyOrder + 1> out{};
778 constexpr Index n = PolyOrder + 1;
779 switch (type) {
780 case GridType::Standard:
781 for (Index j = 0; j < n; j++)
782 out[j] = dl<GridType::Standard>(p0, n, x, xi, li, j);
783 break;
784 case GridType::Log:
785 for (Index j = 0; j < n; j++)
786 out[j] = dl<GridType::Log>(p0, n, x, xi, li, j);
787 break;
788 case GridType::Log10:
789 for (Index j = 0; j < n; j++)
790 out[j] = dl<GridType::Log10>(p0, n, x, xi, li, j);
791 break;
792 case GridType::Log2:
793 for (Index j = 0; j < n; j++)
794 out[j] = dl<GridType::Log2>(p0, n, x, xi, li, j);
795 break;
796 case GridType::SinDeg:
797 for (Index j = 0; j < n; j++)
798 out[j] = dl<GridType::SinDeg>(p0, n, x, xi, li, j);
799 break;
800 case GridType::SinRad:
801 for (Index j = 0; j < n; j++)
802 out[j] = dl<GridType::SinRad>(p0, n, x, xi, li, j);
803 break;
804 case GridType::CosDeg:
805 for (Index j = 0; j < n; j++)
806 out[j] = dl<GridType::CosDeg>(p0, n, x, xi, li, j);
807 break;
808 case GridType::CosRad:
809 for (Index j = 0; j < n; j++)
810 out[j] = dl<GridType::CosRad>(p0, n, x, xi, li, j);
811 break;
812 case GridType::Cyclic:
813 for (Index j = 0; j < n; j++)
814 out[j] = dl<GridType::Cyclic>(p0, n, x, xi, li, j, cycle);
815 break;
816 case GridType::FINAL: { /* pass */
817 }
818 }
819 return out;
820 }
821}; // FixedLagrange
822
826
841 const ConstVectorView& x, const ConstVectorView& xi, const Index polyorder,
842 const Numeric extrapol=0.5, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair<Numeric, Numeric> cycle={-180, 180});
843
855template <std::size_t PolyOrder, class UnsortedVectorType,
856 class SortedVectorType>
858 const UnsortedVectorType& xs, const SortedVectorType& xi,
859 const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair<Numeric, Numeric> cycle={-180, 180}) {
861 out.reserve(xs.size());
862 bool has_one = false;
863 for (auto x : xs) {
864 if (has_one) {
865 out.emplace_back(out.back().pos, x, xi, do_derivs, type, cycle);
866 } else {
867 out.emplace_back(start_pos_finder(x, xi), x, xi, do_derivs, type, cycle);
868 has_one = true;
869 }
870 }
871 return out;
872}
873
877
881
887Vector interpweights(const Lagrange& dim0);
888
895
901template <std::size_t PolyOrder>
902constexpr const std::array<Numeric, PolyOrder + 1>& interpweights(
903 const FixedLagrange<PolyOrder>& dim0) {
904 return dim0.lx;
905}
906
912template <std::size_t PolyOrder>
914 const Array<FixedLagrange<PolyOrder>>& dim0) {
915 Grid<std::array<Numeric, PolyOrder + 1>, 1> out(dim0.size());
916 for (std::size_t i = 0; i < dim0.size(); i++) out(i) = interpweights(dim0[i]);
917 return out;
918}
919
923
929Vector dinterpweights(const Lagrange& dim0);
930
937
943template <std::size_t PolyOrder>
944constexpr std::array<Numeric, PolyOrder + 1> dinterpweights(
945 const FixedLagrange<PolyOrder>& dim0) {
946 return dim0.dlx;
947}
948
954template <std::size_t PolyOrder>
956 const Array<FixedLagrange<PolyOrder>>& dim0) {
957 Grid<std::array<Numeric, PolyOrder + 1>, 1> out(dim0.size());
958 for (std::size_t i = 0; i < dim0.size(); i++)
959 out(i) = dinterpweights(dim0[i]);
960 return out;
961}
962
966
975Numeric interp(const ConstVectorView& yi, const ConstVectorView& iw,
976 const Lagrange& dim0);
977
986template <std::size_t PolyOrder, class VectorType>
987constexpr Numeric interp(const VectorType& yi,
988 const std::array<Numeric, PolyOrder + 1>& iw,
989 const FixedLagrange<PolyOrder>& dim0) {
990 Numeric out(0.0);
991 const Index I = yi.size();
992 for (Index i = 0; i < dim0.size(); i++)
993 out += iw[i] * yi[cycler(i + dim0.pos, I)];
994 return out;
995}
996
1000
1010void reinterp(VectorView out, const ConstVectorView& iy,
1011 const Grid<Vector, 1>& iw, const Array<Lagrange>& dim0);
1012
1021Vector reinterp(const ConstVectorView& iy, const Grid<Vector, 1>& iw,
1022 const Array<Lagrange>& dim0);
1023
1032template <std::size_t PolyOrder>
1034 const Grid<std::array<Numeric, PolyOrder + 1>, 1>& iw,
1035 const Array<FixedLagrange<PolyOrder>>& dim0) {
1036 Vector out(dim0.size());
1037 for (std::size_t i = 0; i < dim0.size(); i++)
1038 out[i] = interp(iy, iw(i), dim0[i]);
1039 return out;
1040}
1041
1045
1049
1056Matrix interpweights(const Lagrange& dim0, const Lagrange& dim1);
1057
1065 const Array<Lagrange>& dim1);
1072template <std::size_t PolyOrder0, std::size_t PolyOrder1>
1074 const FixedLagrange<PolyOrder0>& dim0,
1075 const FixedLagrange<PolyOrder1>& dim1) {
1077 for (Index i = 0; i < dim0.size(); i++)
1078 for (Index j = 0; j < dim1.size(); j++) out(i, j) = dim0.lx[i] * dim1.lx[j];
1079 return out;
1080}
1081
1088template <std::size_t PolyOrder0, std::size_t PolyOrder1>
1090 const Array<FixedLagrange<PolyOrder0>>& dim0,
1091 const Array<FixedLagrange<PolyOrder1>>& dim1) {
1093 dim1.size());
1094 for (std::size_t i = 0; i < dim0.size(); i++)
1095 for (std::size_t j = 0; j < dim1.size(); j++)
1096 out(i, j) = interpweights(dim0[i], dim1[j]);
1097 return out;
1098}
1099
1103
1111Matrix dinterpweights(const Lagrange& dim0, const Lagrange& dim1, Index dim);
1112
1121 const Array<Lagrange>& dim1, Index dim);
1122
1130template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1131 std::size_t DerivativeDim>
1133 const FixedLagrange<PolyOrder0>& dim0,
1134 const FixedLagrange<PolyOrder1>& dim1) {
1135 static_assert(DerivativeDim < 2);
1136
1138 for (Index i = 0; i < dim0.size(); i++)
1139 for (Index j = 0; j < dim1.size(); j++)
1140 out(i, j) = (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1141 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]);
1142 return out;
1143}
1144
1152template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1153 std::size_t DerivativeDim>
1155 const Array<FixedLagrange<PolyOrder0>>& dim0,
1156 const Array<FixedLagrange<PolyOrder1>>& dim1) {
1157 static_assert(DerivativeDim < 2);
1158
1160 dim1.size());
1161 for (std::size_t i = 0; i < dim0.size(); i++)
1162 for (std::size_t j = 0; j < dim1.size(); j++)
1163 out(i, j) = dinterpweights<PolyOrder0, PolyOrder1, DerivativeDim>(
1164 dim0[i], dim1[j]);
1165 return out;
1166}
1167
1171
1181Numeric interp(const ConstMatrixView& yi, const ConstMatrixView& iw,
1182 const Lagrange& dim0, const Lagrange& dim1);
1183
1193template <std::size_t PolyOrder0, std::size_t PolyOrder1, class MatrixType>
1194constexpr Numeric interp(
1195 const MatrixType& yi,
1197 const FixedLagrange<PolyOrder0>& dim0,
1198 const FixedLagrange<PolyOrder1>& dim1) {
1199 Numeric out(0.0);
1200 const Index I = yi.nrows();
1201 const Index J = yi.ncols();
1202 for (Index i = 0; i < dim0.size(); i++)
1203 for (Index j = 0; j < dim1.size(); j++)
1204 out += iw(i, j) * yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J));
1205 return out;
1206}
1207
1211
1222void reinterp(MatrixView out, const ConstMatrixView& iy,
1223 const Grid<Matrix, 2>& iw, const Array<Lagrange>& dim0,
1224 const Array<Lagrange>& dim1);
1225
1235Matrix reinterp(const ConstMatrixView& iy, const Grid<Matrix, 2>& iw,
1236 const Array<Lagrange>& dim0,
1237 const Array<Lagrange>& dim1);
1238
1248template <std::size_t PolyOrder0, std::size_t PolyOrder1>
1250 const ConstMatrixView& iy,
1252 const Array<FixedLagrange<PolyOrder0>>& dim0,
1253 const Array<FixedLagrange<PolyOrder1>>& dim1) {
1254 Matrix out(dim0.size(), dim1.size());
1255 for (std::size_t i = 0; i < dim0.size(); i++)
1256 for (std::size_t j = 0; j < dim1.size(); j++)
1257 out(i, j) = interp(iy, iw(i, j), dim0[i], dim1[j]);
1258 return out;
1259}
1260
1264
1268
1276Tensor3 interpweights(const Lagrange& dim0, const Lagrange& dim1,
1277 const Lagrange& dim2);
1278
1287 const Array<Lagrange>& dim1,
1288 const Array<Lagrange>& dim2);
1289
1297template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1298 std::size_t PolyOrder2>
1301 const FixedLagrange<PolyOrder1>& dim1,
1302 const FixedLagrange<PolyOrder2>& dim2) {
1304 for (Index i = 0; i < dim0.size(); i++)
1305 for (Index j = 0; j < dim1.size(); j++)
1306 for (Index k = 0; k < dim2.size(); k++)
1307 out(i, j, k) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k];
1308 return out;
1309}
1310
1318template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1319 std::size_t PolyOrder2>
1322 const Array<FixedLagrange<PolyOrder1>>& dim1,
1323 const Array<FixedLagrange<PolyOrder2>>& dim2) {
1325 out(dim0.size(), dim1.size(), dim2.size());
1326 for (std::size_t i = 0; i < dim0.size(); i++)
1327 for (std::size_t j = 0; j < dim1.size(); j++)
1328 for (std::size_t k = 0; k < dim2.size(); k++)
1329 out(i, j, k) = interpweights(dim0[i], dim1[j], dim2[k]);
1330 return out;
1331}
1332
1336
1345Tensor3 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
1346 const Lagrange& dim2, Index dim);
1347
1357 const Array<Lagrange>& dim1,
1358 const Array<Lagrange>& dim2, Index dim);
1359
1368template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1369 std::size_t PolyOrder2, std::size_t DerivativeDim>
1372 const FixedLagrange<PolyOrder1>& dim1,
1373 const FixedLagrange<PolyOrder2>& dim2) {
1374 static_assert(DerivativeDim < 3);
1375
1377 for (Index i = 0; i < dim0.size(); i++)
1378 for (Index j = 0; j < dim1.size(); j++)
1379 for (Index k = 0; k < dim2.size(); k++)
1380 out(i, j, k) = (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1381 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
1382 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]);
1383 return out;
1384}
1385
1394template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1395 std::size_t PolyOrder2, std::size_t DerivativeDim>
1398 const Array<FixedLagrange<PolyOrder1>>& dim1,
1399 const Array<FixedLagrange<PolyOrder2>>& dim2) {
1401 out(dim0.size(), dim0.size(), dim2.size());
1402 for (std::size_t i = 0; i < dim0.size(); i++)
1403 for (std::size_t j = 0; j < dim1.size(); j++)
1404 for (std::size_t k = 0; k < dim2.size(); k++)
1405 out(i, j, k) =
1406 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, DerivativeDim>(
1407 dim0[i], dim1[j], dim2[k]);
1408 return out;
1409}
1410
1414
1425Numeric interp(const ConstTensor3View& yi, const ConstTensor3View& iw,
1426 const Lagrange& dim0, const Lagrange& dim1,
1427 const Lagrange& dim2);
1428
1439template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1440 std::size_t PolyOrder2, class Tensor3Type>
1441constexpr Numeric interp(const Tensor3Type& yi,
1442 const FixedGrid<Numeric, PolyOrder0 + 1,
1443 PolyOrder1 + 1, PolyOrder2 + 1>& iw,
1444 const FixedLagrange<PolyOrder0>& dim0,
1445 const FixedLagrange<PolyOrder1>& dim1,
1446 const FixedLagrange<PolyOrder2>& dim2) {
1447 Numeric out(0.0);
1448 const Index I = yi.npages();
1449 const Index J = yi.nrows();
1450 const Index K = yi.ncols();
1451 for (Index i = 0; i < dim0.size(); i++)
1452 for (Index j = 0; j < dim1.size(); j++)
1453 for (Index k = 0; k < dim2.size(); k++)
1454 out +=
1455 iw(i, j, k) * yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
1456 cycler(k + dim2.pos, K));
1457 return out;
1458}
1459
1463
1475void reinterp(Tensor3View out, const ConstTensor3View& iy,
1476 const Grid<Tensor3, 3>& iw, const Array<Lagrange>& dim0,
1477 const Array<Lagrange>& dim1,
1478 const Array<Lagrange>& dim2);
1479
1491 const Array<Lagrange>& dim0,
1492 const Array<Lagrange>& dim1,
1493 const Array<Lagrange>& dim2);
1494
1505template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1506 std::size_t PolyOrder2>
1508 const ConstTensor3View& iy,
1509 const Grid<
1511 iw,
1512 const Array<FixedLagrange<PolyOrder0>>& dim0,
1513 const Array<FixedLagrange<PolyOrder1>>& dim1,
1514 const Array<FixedLagrange<PolyOrder2>>& dim2) {
1515 Tensor3 out(dim0.size(), dim1.size(), dim2.size());
1516 for (std::size_t i = 0; i < dim0.size(); i++)
1517 for (std::size_t j = 0; j < dim1.size(); j++)
1518 for (std::size_t k = 0; k < dim2.size(); k++)
1519 out(i, j, k) = interp(iy, iw(i, j, k), dim0[i], dim1[j], dim2[k]);
1520 return out;
1521}
1522
1526
1530
1539Tensor4 interpweights(const Lagrange& dim0, const Lagrange& dim1,
1540 const Lagrange& dim2, const Lagrange& dim3);
1541
1551 const Array<Lagrange>& dim1,
1552 const Array<Lagrange>& dim2,
1553 const Array<Lagrange>& dim3);
1554
1563template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1564 std::size_t PolyOrder2, std::size_t PolyOrder3>
1565constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1566 PolyOrder3 + 1>
1568 const FixedLagrange<PolyOrder1>& dim1,
1569 const FixedLagrange<PolyOrder2>& dim2,
1570 const FixedLagrange<PolyOrder3>& dim3) {
1571 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1572 PolyOrder3 + 1>
1573 out;
1574 for (Index i = 0; i < dim0.size(); i++)
1575 for (Index j = 0; j < dim1.size(); j++)
1576 for (Index k = 0; k < dim2.size(); k++)
1577 for (Index l = 0; l < dim3.size(); l++)
1578 out(i, j, k, l) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k] * dim3.lx[l];
1579 return out;
1580}
1581
1590template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1591 std::size_t PolyOrder2, std::size_t PolyOrder3>
1592Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1593 PolyOrder3 + 1>,
1594 4>
1596 const Array<FixedLagrange<PolyOrder1>>& dim1,
1597 const Array<FixedLagrange<PolyOrder2>>& dim2,
1598 const Array<FixedLagrange<PolyOrder3>>& dim3) {
1599 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1600 PolyOrder3 + 1>,
1601 4>
1602 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
1603 for (std::size_t i = 0; i < dim0.size(); i++)
1604 for (std::size_t j = 0; j < dim1.size(); j++)
1605 for (std::size_t k = 0; k < dim2.size(); k++)
1606 for (std::size_t l = 0; l < dim3.size(); l++)
1607 out(i, j, k, l) = interpweights(dim0[i], dim1[j], dim2[k], dim3[l]);
1608 return out;
1609}
1610
1614
1624Tensor4 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
1625 const Lagrange& dim2, const Lagrange& dim3, Index dim);
1626
1637 const Array<Lagrange>& dim1,
1638 const Array<Lagrange>& dim2,
1639 const Array<Lagrange>& dim3, Index dim);
1640
1650template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1651 std::size_t PolyOrder2, std::size_t PolyOrder3,
1652 std::size_t DerivativeDim>
1653constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1654 PolyOrder3 + 1>
1656 const FixedLagrange<PolyOrder1>& dim1,
1657 const FixedLagrange<PolyOrder2>& dim2,
1658 const FixedLagrange<PolyOrder3>& dim3) {
1659 static_assert(DerivativeDim < 4);
1660
1661 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1662 PolyOrder3 + 1>
1663 out;
1664 for (Index i = 0; i < dim0.size(); i++)
1665 for (Index j = 0; j < dim1.size(); j++)
1666 for (Index k = 0; k < dim2.size(); k++)
1667 for (Index l = 0; l < dim3.size(); l++)
1668 out(i, j, k, l) = (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
1669 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
1670 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
1671 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]);
1672 return out;
1673}
1674
1684template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1685 std::size_t PolyOrder2, std::size_t PolyOrder3,
1686 std::size_t DerivativeDim>
1687Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1688 PolyOrder3 + 1>,
1689 4>
1691 const Array<FixedLagrange<PolyOrder1>>& dim1,
1692 const Array<FixedLagrange<PolyOrder2>>& dim2,
1693 const Array<FixedLagrange<PolyOrder3>>& dim3) {
1694 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1695 PolyOrder3 + 1>,
1696 4>
1697 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
1698 for (std::size_t i = 0; i < dim0.size(); i++)
1699 for (std::size_t j = 0; j < dim1.size(); j++)
1700 for (std::size_t k = 0; k < dim2.size(); k++)
1701 for (std::size_t l = 0; l < dim3.size(); l++)
1702 out(i, j, k, l) =
1703 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, PolyOrder3,
1704 DerivativeDim>(dim0[i], dim1[j], dim2[k], dim3[l]);
1705 return out;
1706}
1707
1711
1723Numeric interp(const ConstTensor4View& yi, const ConstTensor4View& iw,
1724 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
1725 const Lagrange& dim3);
1726
1738template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1739 std::size_t PolyOrder2, std::size_t PolyOrder3, class Tensor4Type>
1740constexpr Numeric interp(
1741 const Tensor4Type& yi,
1742 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1743 PolyOrder3 + 1>& iw,
1744 const FixedLagrange<PolyOrder0>& dim0,
1745 const FixedLagrange<PolyOrder1>& dim1,
1746 const FixedLagrange<PolyOrder2>& dim2,
1747 const FixedLagrange<PolyOrder3>& dim3) {
1748 Numeric out(0.0);
1749 const Index I = yi.nbooks();
1750 const Index J = yi.npages();
1751 const Index K = yi.nrows();
1752 const Index L = yi.ncols();
1753 for (Index i = 0; i < dim0.size(); i++)
1754 for (Index j = 0; j < dim1.size(); j++)
1755 for (Index k = 0; k < dim2.size(); k++)
1756 for (Index l = 0; l < dim3.size(); l++)
1757 out += iw(i, j, k, l) *
1758 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
1759 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L));
1760 return out;
1761}
1762
1766
1779void reinterp(Tensor4View out, const ConstTensor4View& iy,
1780 const Grid<Tensor4, 4>& iw, const Array<Lagrange>& dim0,
1781 const Array<Lagrange>& dim1,
1782 const Array<Lagrange>& dim2,
1783 const Array<Lagrange>& dim3);
1784
1797 const Array<Lagrange>& dim0,
1798 const Array<Lagrange>& dim1,
1799 const Array<Lagrange>& dim2,
1800 const Array<Lagrange>& dim3);
1801
1813template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1814 std::size_t PolyOrder2, std::size_t PolyOrder3>
1816 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
1817 PolyOrder2 + 1, PolyOrder3 + 1>,
1818 4>& iw,
1819 const Array<FixedLagrange<PolyOrder0>>& dim0,
1820 const Array<FixedLagrange<PolyOrder1>>& dim1,
1821 const Array<FixedLagrange<PolyOrder2>>& dim2,
1822 const Array<FixedLagrange<PolyOrder3>>& dim3) {
1823 Tensor4 out(dim0.size(), dim1.size(), dim2.size(), dim3.size());
1824 for (std::size_t i = 0; i < dim0.size(); i++)
1825 for (std::size_t j = 0; j < dim1.size(); j++)
1826 for (std::size_t k = 0; k < dim2.size(); k++)
1827 for (std::size_t l = 0; l < dim3.size(); l++)
1828 out(i, j, k, l) =
1829 interp(iy, iw(i, j, k, l), dim0[i], dim1[j], dim2[k], dim3[l]);
1830 return out;
1831}
1832
1836
1840
1850Tensor5 interpweights(const Lagrange& dim0, const Lagrange& dim1,
1851 const Lagrange& dim2, const Lagrange& dim3,
1852 const Lagrange& dim4);
1853
1864 const Array<Lagrange>& dim1,
1865 const Array<Lagrange>& dim2,
1866 const Array<Lagrange>& dim3,
1867 const Array<Lagrange>& dim4);
1868
1878template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1879 std::size_t PolyOrder2, std::size_t PolyOrder3,
1880 std::size_t PolyOrder4>
1881constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1882 PolyOrder3 + 1, PolyOrder4 + 1>
1884 const FixedLagrange<PolyOrder1>& dim1,
1885 const FixedLagrange<PolyOrder2>& dim2,
1886 const FixedLagrange<PolyOrder3>& dim3,
1887 const FixedLagrange<PolyOrder4>& dim4) {
1888 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1889 PolyOrder3 + 1, PolyOrder4 + 1>
1890 out;
1891 for (Index i = 0; i < dim0.size(); i++)
1892 for (Index j = 0; j < dim1.size(); j++)
1893 for (Index k = 0; k < dim2.size(); k++)
1894 for (Index l = 0; l < dim3.size(); l++)
1895 for (Index m = 0; m < dim4.size(); m++)
1896 out(i, j, k, l, m) =
1897 dim0.lx[i] * dim1.lx[j] * dim2.lx[k] * dim3.lx[l] * dim4.lx[m];
1898 return out;
1899}
1900
1910template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1911 std::size_t PolyOrder2, std::size_t PolyOrder3,
1912 std::size_t PolyOrder4>
1913Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1914 PolyOrder3 + 1, PolyOrder4 + 1>,
1915 5>
1917 const Array<FixedLagrange<PolyOrder1>>& dim1,
1918 const Array<FixedLagrange<PolyOrder2>>& dim2,
1919 const Array<FixedLagrange<PolyOrder3>>& dim3,
1920 const Array<FixedLagrange<PolyOrder4>>& dim4) {
1921 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1922 PolyOrder3 + 1, PolyOrder4 + 1>,
1923 5>
1924 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
1925 for (std::size_t i = 0; i < dim0.size(); i++)
1926 for (std::size_t j = 0; j < dim1.size(); j++)
1927 for (std::size_t k = 0; k < dim2.size(); k++)
1928 for (std::size_t l = 0; l < dim3.size(); l++)
1929 for (std::size_t m = 0; m < dim4.size(); m++)
1930 out(i, j, k, l, m) =
1931 interpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m]);
1932 return out;
1933}
1934
1938
1949Tensor5 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
1950 const Lagrange& dim2, const Lagrange& dim3,
1951 const Lagrange& dim4, Index dim);
1952
1964 const Array<Lagrange>& dim1,
1965 const Array<Lagrange>& dim2,
1966 const Array<Lagrange>& dim3,
1967 const Array<Lagrange>& dim4, Index dim);
1968
1979template <std::size_t PolyOrder0, std::size_t PolyOrder1,
1980 std::size_t PolyOrder2, std::size_t PolyOrder3,
1981 std::size_t PolyOrder4, std::size_t DerivativeDim>
1982constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1983 PolyOrder3 + 1, PolyOrder4 + 1>
1985 const FixedLagrange<PolyOrder1>& dim1,
1986 const FixedLagrange<PolyOrder2>& dim2,
1987 const FixedLagrange<PolyOrder3>& dim3,
1988 const FixedLagrange<PolyOrder4>& dim4) {
1989 static_assert(DerivativeDim < 5);
1990
1991 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
1992 PolyOrder3 + 1, PolyOrder4 + 1>
1993 out;
1994 for (Index i = 0; i < dim0.size(); i++)
1995 for (Index j = 0; j < dim1.size(); j++)
1996 for (Index k = 0; k < dim2.size(); k++)
1997 for (Index l = 0; l < dim3.size(); l++)
1998 for (Index m = 0; m < dim4.size(); m++)
1999 out(i, j, k, l, m) =
2000 (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
2001 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
2002 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
2003 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
2004 (DerivativeDim == 4 ? dim4.dlx[m] : dim4.lx[m]);
2005 return out;
2006}
2007
2018template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2019 std::size_t PolyOrder2, std::size_t PolyOrder3,
2020 std::size_t PolyOrder4, std::size_t DerivativeDim>
2021Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2022 PolyOrder3 + 1, PolyOrder4 + 1>,
2023 5>
2025 const Array<FixedLagrange<PolyOrder1>>& dim1,
2026 const Array<FixedLagrange<PolyOrder2>>& dim2,
2027 const Array<FixedLagrange<PolyOrder3>>& dim3,
2028 const Array<FixedLagrange<PolyOrder4>>& dim4) {
2029 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2030 PolyOrder3 + 1, PolyOrder4 + 1>,
2031 5>
2032 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
2033 for (std::size_t i = 0; i < dim0.size(); i++)
2034 for (std::size_t j = 0; j < dim1.size(); j++)
2035 for (std::size_t k = 0; k < dim2.size(); k++)
2036 for (std::size_t l = 0; l < dim3.size(); l++)
2037 for (std::size_t m = 0; m < dim4.size(); m++)
2038 out(i, j, k, l, m) =
2039 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, PolyOrder3,
2040 PolyOrder4, DerivativeDim>(
2041 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m]);
2042 return out;
2043}
2044
2048
2061Numeric interp(const ConstTensor5View& yi, const ConstTensor5View& iw,
2062 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
2063 const Lagrange& dim3, const Lagrange& dim4);
2064
2077template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2078 std::size_t PolyOrder2, std::size_t PolyOrder3,
2079 std::size_t PolyOrder4, class Tensor5Type>
2080constexpr Numeric interp(
2081 const Tensor5Type& yi,
2082 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2083 PolyOrder3 + 1, PolyOrder4 + 1>& iw,
2084 const FixedLagrange<PolyOrder0>& dim0,
2085 const FixedLagrange<PolyOrder1>& dim1,
2086 const FixedLagrange<PolyOrder2>& dim2,
2087 const FixedLagrange<PolyOrder3>& dim3,
2088 const FixedLagrange<PolyOrder4>& dim4) {
2089 Numeric out(0.0);
2090 const Index I = yi.nshelves();
2091 const Index J = yi.nbooks();
2092 const Index K = yi.npages();
2093 const Index L = yi.nrows();
2094 const Index M = yi.ncols();
2095 for (Index i = 0; i < dim0.size(); i++)
2096 for (Index j = 0; j < dim1.size(); j++)
2097 for (Index k = 0; k < dim2.size(); k++)
2098 for (Index l = 0; l < dim3.size(); l++)
2099 for (Index m = 0; m < dim4.size(); m++)
2100 out += iw(i, j, k, l, m) *
2101 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
2102 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
2103 cycler(m + dim4.pos, M));
2104 return out;
2105}
2106
2110
2124void reinterp(Tensor5View out, const ConstTensor5View& iy,
2125 const Grid<Tensor5, 5>& iw, const Array<Lagrange>& dim0,
2126 const Array<Lagrange>& dim1,
2127 const Array<Lagrange>& dim2,
2128 const Array<Lagrange>& dim3,
2129 const Array<Lagrange>& dim4);
2130
2144 const Array<Lagrange>& dim0,
2145 const Array<Lagrange>& dim1,
2146 const Array<Lagrange>& dim2,
2147 const Array<Lagrange>& dim3,
2148 const Array<Lagrange>& dim4);
2149
2162template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2163 std::size_t PolyOrder2, std::size_t PolyOrder3,
2164 std::size_t PolyOrder4>
2166 const ConstTensor5View& iy,
2167 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
2168 PolyOrder2 + 1, PolyOrder3 + 1, PolyOrder4 + 1>,
2169 5>& iw,
2170 const FixedLagrange<PolyOrder0>& dim0,
2171 const FixedLagrange<PolyOrder1>& dim1,
2172 const FixedLagrange<PolyOrder2>& dim2,
2173 const FixedLagrange<PolyOrder3>& dim3,
2174 const FixedLagrange<PolyOrder4>& dim4) {
2175 Tensor5 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size());
2176 for (std::size_t i = 0; i < dim0.size(); i++)
2177 for (std::size_t j = 0; j < dim1.size(); j++)
2178 for (std::size_t k = 0; k < dim2.size(); k++)
2179 for (std::size_t l = 0; l < dim3.size(); l++)
2180 for (std::size_t m = 0; m < dim4.size(); m++)
2181 out(i, j, k, l, m) = interp(iy, iw(i, j, k, l, m), dim0[i], dim1[j],
2182 dim2[k], dim3[l], dim4[m]);
2183 return out;
2184}
2185
2189
2193
2204Tensor6 interpweights(const Lagrange& dim0, const Lagrange& dim1,
2205 const Lagrange& dim2, const Lagrange& dim3,
2206 const Lagrange& dim4, const Lagrange& dim5);
2207
2219 const Array<Lagrange>& dim1,
2220 const Array<Lagrange>& dim2,
2221 const Array<Lagrange>& dim3,
2222 const Array<Lagrange>& dim4,
2223 const Array<Lagrange>& dim5);
2224
2235template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2236 std::size_t PolyOrder2, std::size_t PolyOrder3,
2237 std::size_t PolyOrder4, std::size_t PolyOrder5>
2238constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2239 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2241 const FixedLagrange<PolyOrder1>& dim1,
2242 const FixedLagrange<PolyOrder2>& dim2,
2243 const FixedLagrange<PolyOrder3>& dim3,
2244 const FixedLagrange<PolyOrder4>& dim4,
2245 const FixedLagrange<PolyOrder5>& dim5) {
2246 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2247 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2248 out;
2249 for (Index i = 0; i < dim0.size(); i++)
2250 for (Index j = 0; j < dim1.size(); j++)
2251 for (Index k = 0; k < dim2.size(); k++)
2252 for (Index l = 0; l < dim3.size(); l++)
2253 for (Index m = 0; m < dim4.size(); m++)
2254 for (Index n = 0; n < dim5.size(); n++)
2255 out(i, j, k, l, m, n) = dim0.lx[i] * dim1.lx[j] * dim2.lx[k] *
2256 dim3.lx[l] * dim4.lx[m] * dim5.lx[n];
2257 return out;
2258}
2259
2270template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2271 std::size_t PolyOrder2, std::size_t PolyOrder3,
2272 std::size_t PolyOrder4, std::size_t PolyOrder5>
2273Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2274 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2275 6>
2277 const Array<FixedLagrange<PolyOrder1>>& dim1,
2278 const Array<FixedLagrange<PolyOrder2>>& dim2,
2279 const Array<FixedLagrange<PolyOrder3>>& dim3,
2280 const Array<FixedLagrange<PolyOrder4>>& dim4,
2281 const Array<FixedLagrange<PolyOrder5>>& dim5) {
2282 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2283 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2284 6>
2285 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2286 dim5.size());
2287 for (std::size_t i = 0; i < dim0.size(); i++)
2288 for (std::size_t j = 0; j < dim1.size(); j++)
2289 for (std::size_t k = 0; k < dim2.size(); k++)
2290 for (std::size_t l = 0; l < dim3.size(); l++)
2291 for (std::size_t m = 0; m < dim4.size(); m++)
2292 for (std::size_t n = 0; n < dim5.size(); n++)
2293 out(i, j, k, l, m, n) = interpweights(dim0[i], dim1[j], dim2[k],
2294 dim3[l], dim4[m], dim5[n]);
2295 return out;
2296}
2297
2301
2313Tensor6 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
2314 const Lagrange& dim2, const Lagrange& dim3,
2315 const Lagrange& dim4, const Lagrange& dim5, Index dim);
2316
2329 const Array<Lagrange>& dim1,
2330 const Array<Lagrange>& dim2,
2331 const Array<Lagrange>& dim3,
2332 const Array<Lagrange>& dim4,
2333 const Array<Lagrange>& dim5, Index dim);
2334
2346template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2347 std::size_t PolyOrder2, std::size_t PolyOrder3,
2348 std::size_t PolyOrder4, std::size_t PolyOrder5,
2349 std::size_t DerivativeDim>
2350constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2351 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2353 const FixedLagrange<PolyOrder1>& dim1,
2354 const FixedLagrange<PolyOrder2>& dim2,
2355 const FixedLagrange<PolyOrder3>& dim3,
2356 const FixedLagrange<PolyOrder4>& dim4,
2357 const FixedLagrange<PolyOrder5>& dim5) {
2358 static_assert(DerivativeDim < 6);
2359
2360 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2361 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>
2362 out;
2363 for (Index i = 0; i < dim0.size(); i++)
2364 for (Index j = 0; j < dim1.size(); j++)
2365 for (Index k = 0; k < dim2.size(); k++)
2366 for (Index l = 0; l < dim3.size(); l++)
2367 for (Index m = 0; m < dim4.size(); m++)
2368 for (Index n = 0; n < dim5.size(); n++)
2369 out(i, j, k, l, m, n) =
2370 (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
2371 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
2372 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
2373 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
2374 (DerivativeDim == 4 ? dim4.dlx[m] : dim4.lx[m]) *
2375 (DerivativeDim == 5 ? dim5.dlx[n] : dim5.lx[n]);
2376 return out;
2377}
2378
2390template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2391 std::size_t PolyOrder2, std::size_t PolyOrder3,
2392 std::size_t PolyOrder4, std::size_t PolyOrder5,
2393 std::size_t DerivativeDim>
2394Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2395 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2396 6>
2398 const Array<FixedLagrange<PolyOrder1>>& dim1,
2399 const Array<FixedLagrange<PolyOrder2>>& dim2,
2400 const Array<FixedLagrange<PolyOrder3>>& dim3,
2401 const Array<FixedLagrange<PolyOrder4>>& dim4,
2402 const Array<FixedLagrange<PolyOrder5>>& dim5) {
2403 Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2404 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>,
2405 6>
2406 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2407 dim5.size());
2408 for (std::size_t i = 0; i < dim0.size(); i++)
2409 for (std::size_t j = 0; j < dim1.size(); j++)
2410 for (std::size_t k = 0; k < dim2.size(); k++)
2411 for (std::size_t l = 0; l < dim3.size(); l++)
2412 for (std::size_t m = 0; m < dim4.size(); m++)
2413 for (std::size_t n = 0; n < dim5.size(); n++)
2414 out(i, j, k, l, m, n) =
2415 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2, PolyOrder3,
2416 PolyOrder4, PolyOrder5, DerivativeDim>(
2417 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m], dim5[n]);
2418 return out;
2419}
2420
2424
2438Numeric interp(const ConstTensor6View& yi, const ConstTensor6View& iw,
2439 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
2440 const Lagrange& dim3, const Lagrange& dim4,
2441 const Lagrange& dim5);
2442
2456template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2457 std::size_t PolyOrder2, std::size_t PolyOrder3,
2458 std::size_t PolyOrder4, std::size_t PolyOrder5, class Tensor6Type>
2459constexpr Numeric interp(
2460 const Tensor6Type& yi,
2461 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2462 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1>& iw,
2463 const FixedLagrange<PolyOrder0>& dim0,
2464 const FixedLagrange<PolyOrder1>& dim1,
2465 const FixedLagrange<PolyOrder2>& dim2,
2466 const FixedLagrange<PolyOrder3>& dim3,
2467 const FixedLagrange<PolyOrder4>& dim4,
2468 const FixedLagrange<PolyOrder5>& dim5) {
2469 Numeric out(0.0);
2470 const Index I = yi.nvitrines();
2471 const Index J = yi.nshelves();
2472 const Index K = yi.nbooks();
2473 const Index L = yi.npages();
2474 const Index M = yi.nrows();
2475 const Index N = yi.ncols();
2476 for (Index i = 0; i < dim0.size(); i++)
2477 for (Index j = 0; j < dim1.size(); j++)
2478 for (Index k = 0; k < dim2.size(); k++)
2479 for (Index l = 0; l < dim3.size(); l++)
2480 for (Index m = 0; m < dim4.size(); m++)
2481 for (Index n = 0; n < dim5.size(); n++)
2482 out += iw(i, j, k, l, m, n) *
2483 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
2484 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
2485 cycler(m + dim4.pos, M), cycler(n + dim5.pos, N));
2486 return out;
2487}
2488
2492
2507void reinterp(Tensor6View out, const ConstTensor6View& iy,
2508 const Grid<Tensor6, 6>& iw, const Array<Lagrange>& dim0,
2509 const Array<Lagrange>& dim1,
2510 const Array<Lagrange>& dim2,
2511 const Array<Lagrange>& dim3,
2512 const Array<Lagrange>& dim4,
2513 const Array<Lagrange>& dim5);
2514
2529 const Array<Lagrange>& dim0,
2530 const Array<Lagrange>& dim1,
2531 const Array<Lagrange>& dim2,
2532 const Array<Lagrange>& dim3,
2533 const Array<Lagrange>& dim4,
2534 const Array<Lagrange>& dim5);
2535
2549template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2550 std::size_t PolyOrder2, std::size_t PolyOrder3,
2551 std::size_t PolyOrder4, std::size_t PolyOrder5>
2553 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
2554 PolyOrder2 + 1, PolyOrder3 + 1,
2555 PolyOrder4 + 1, PolyOrder5 + 1>,
2556 6>& iw,
2557 const Array<FixedLagrange<PolyOrder0>>& dim0,
2558 const Array<FixedLagrange<PolyOrder1>>& dim1,
2559 const Array<FixedLagrange<PolyOrder2>>& dim2,
2560 const Array<FixedLagrange<PolyOrder3>>& dim3,
2561 const Array<FixedLagrange<PolyOrder4>>& dim4,
2562 const Array<FixedLagrange<PolyOrder5>>& dim5) {
2563 Tensor6 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2564 dim5.size());
2565 for (std::size_t i = 0; i < dim0.size(); i++)
2566 for (std::size_t j = 0; j < dim1.size(); j++)
2567 for (std::size_t k = 0; k < dim2.size(); k++)
2568 for (std::size_t l = 0; l < dim3.size(); l++)
2569 for (std::size_t m = 0; m < dim4.size(); m++)
2570 for (std::size_t n = 0; n < dim5.size(); n++)
2571 out(i, j, k, l, m, n) =
2572 interp(iy, iw(i, j, k, l, m, n), dim0[i], dim1[j], dim2[k],
2573 dim3[l], dim4[m], dim5[n]);
2574 return out;
2575}
2576
2580
2584
2596Tensor7 interpweights(const Lagrange& dim0, const Lagrange& dim1,
2597 const Lagrange& dim2, const Lagrange& dim3,
2598 const Lagrange& dim4, const Lagrange& dim5,
2599 const Lagrange& dim6);
2600
2613 const Array<Lagrange>& dim1,
2614 const Array<Lagrange>& dim2,
2615 const Array<Lagrange>& dim3,
2616 const Array<Lagrange>& dim4,
2617 const Array<Lagrange>& dim5,
2618 const Array<Lagrange>& dim6);
2619
2631template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2632 std::size_t PolyOrder2, std::size_t PolyOrder3,
2633 std::size_t PolyOrder4, std::size_t PolyOrder5,
2634 std::size_t PolyOrder6>
2635constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2636 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1,
2637 PolyOrder6 + 1>
2639 const FixedLagrange<PolyOrder1>& dim1,
2640 const FixedLagrange<PolyOrder2>& dim2,
2641 const FixedLagrange<PolyOrder3>& dim3,
2642 const FixedLagrange<PolyOrder4>& dim4,
2643 const FixedLagrange<PolyOrder5>& dim5,
2644 const FixedLagrange<PolyOrder6>& dim6) {
2645 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2646 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>
2647 out;
2648 for (Index i = 0; i < dim0.size(); i++)
2649 for (Index j = 0; j < dim1.size(); j++)
2650 for (Index k = 0; k < dim2.size(); k++)
2651 for (Index l = 0; l < dim3.size(); l++)
2652 for (Index m = 0; m < dim4.size(); m++)
2653 for (Index n = 0; n < dim5.size(); n++)
2654 for (Index o = 0; o < dim6.size(); o++)
2655 out(i, j, k, l, m, n, o) = dim0.lx[i] * dim1.lx[j] *
2656 dim2.lx[k] * dim3.lx[l] *
2657 dim4.lx[m] * dim5.lx[n] * dim6.lx[o];
2658 return out;
2659}
2660
2672template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2673 std::size_t PolyOrder2, std::size_t PolyOrder3,
2674 std::size_t PolyOrder4, std::size_t PolyOrder5,
2675 std::size_t PolyOrder6>
2676Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2677 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2678 7>
2680 const Array<FixedLagrange<PolyOrder1>>& dim1,
2681 const Array<FixedLagrange<PolyOrder2>>& dim2,
2682 const Array<FixedLagrange<PolyOrder3>>& dim3,
2683 const Array<FixedLagrange<PolyOrder4>>& dim4,
2684 const Array<FixedLagrange<PolyOrder5>>& dim5,
2685 const Array<FixedLagrange<PolyOrder6>>& dim6) {
2686 Grid<
2687 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2688 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2689 7>
2690 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2691 dim5.size(), dim6.size());
2692 for (std::size_t i = 0; i < dim0.size(); i++)
2693 for (std::size_t j = 0; j < dim1.size(); j++)
2694 for (std::size_t k = 0; k < dim2.size(); k++)
2695 for (std::size_t l = 0; l < dim3.size(); l++)
2696 for (std::size_t m = 0; m < dim4.size(); m++)
2697 for (std::size_t n = 0; n < dim5.size(); n++)
2698 for (std::size_t o = 0; o < dim6.size(); o++)
2699 out(i, j, k, l, m, n, o) =
2700 interpweights(dim0[i], dim1[j], dim2[k], dim3[l], dim4[m],
2701 dim5[n], dim6[o]);
2702 return out;
2703}
2704
2708
2721Tensor7 dinterpweights(const Lagrange& dim0, const Lagrange& dim1,
2722 const Lagrange& dim2, const Lagrange& dim3,
2723 const Lagrange& dim4, const Lagrange& dim5,
2724 const Lagrange& dim6, Index dim);
2725
2739 const Array<Lagrange>& dim1,
2740 const Array<Lagrange>& dim2,
2741 const Array<Lagrange>& dim3,
2742 const Array<Lagrange>& dim4,
2743 const Array<Lagrange>& dim5,
2744 const Array<Lagrange>& dim6, Index dim);
2745
2758template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2759 std::size_t PolyOrder2, std::size_t PolyOrder3,
2760 std::size_t PolyOrder4, std::size_t PolyOrder5,
2761 std::size_t PolyOrder6, std::size_t DerivativeDim>
2762constexpr FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2763 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1,
2764 PolyOrder6 + 1>
2766 const FixedLagrange<PolyOrder1>& dim1,
2767 const FixedLagrange<PolyOrder2>& dim2,
2768 const FixedLagrange<PolyOrder3>& dim3,
2769 const FixedLagrange<PolyOrder4>& dim4,
2770 const FixedLagrange<PolyOrder5>& dim5,
2771 const FixedLagrange<PolyOrder6>& dim6) {
2772 static_assert(DerivativeDim < 7);
2773
2774 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2775 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>
2776 out;
2777 for (Index i = 0; i < dim0.size(); i++)
2778 for (Index j = 0; j < dim1.size(); j++)
2779 for (Index k = 0; k < dim2.size(); k++)
2780 for (Index l = 0; l < dim3.size(); l++)
2781 for (Index m = 0; m < dim4.size(); m++)
2782 for (Index n = 0; n < dim5.size(); n++)
2783 for (Index o = 0; o < dim6.size(); o++)
2784 out(i, j, k, l, m, n, o) =
2785 (DerivativeDim == 0 ? dim0.dlx[i] : dim0.lx[i]) *
2786 (DerivativeDim == 1 ? dim1.dlx[j] : dim1.lx[j]) *
2787 (DerivativeDim == 2 ? dim2.dlx[k] : dim2.lx[k]) *
2788 (DerivativeDim == 3 ? dim3.dlx[l] : dim3.lx[l]) *
2789 (DerivativeDim == 4 ? dim4.dlx[m] : dim4.lx[m]) *
2790 (DerivativeDim == 5 ? dim5.dlx[n] : dim5.lx[n]) *
2791 (DerivativeDim == 6 ? dim6.dlx[o] : dim6.lx[o]);
2792 return out;
2793}
2794
2807template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2808 std::size_t PolyOrder2, std::size_t PolyOrder3,
2809 std::size_t PolyOrder4, std::size_t PolyOrder5,
2810 std::size_t PolyOrder6, std::size_t DerivativeDim>
2811Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2812 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2813 7>
2815 const Array<FixedLagrange<PolyOrder1>>& dim1,
2816 const Array<FixedLagrange<PolyOrder2>>& dim2,
2817 const Array<FixedLagrange<PolyOrder3>>& dim3,
2818 const Array<FixedLagrange<PolyOrder4>>& dim4,
2819 const Array<FixedLagrange<PolyOrder5>>& dim5,
2820 const Array<FixedLagrange<PolyOrder6>>& dim6) {
2821 Grid<
2822 FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2823 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1, PolyOrder6 + 1>,
2824 7>
2825 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
2826 dim5.size(), dim6.size());
2827 for (std::size_t i = 0; i < dim0.size(); i++)
2828 for (std::size_t j = 0; j < dim1.size(); j++)
2829 for (std::size_t k = 0; k < dim2.size(); k++)
2830 for (std::size_t l = 0; l < dim3.size(); l++)
2831 for (std::size_t m = 0; m < dim4.size(); m++)
2832 for (std::size_t n = 0; n < dim5.size(); n++)
2833 for (std::size_t o = 0; o < dim6.size(); o++)
2834 out(i, j, k, l, m, n, o) =
2835 dinterpweights<PolyOrder0, PolyOrder1, PolyOrder2,
2836 PolyOrder3, PolyOrder4, PolyOrder5,
2837 PolyOrder6, DerivativeDim>(
2838 dim0[i], dim1[j], dim2[k], dim3[l], dim4[m], dim5[n],
2839 dim6[o]);
2840 return out;
2841}
2842
2846
2861Numeric interp(const ConstTensor7View& yi, const ConstTensor7View& iw,
2862 const Lagrange& dim0, const Lagrange& dim1, const Lagrange& dim2,
2863 const Lagrange& dim3, const Lagrange& dim4, const Lagrange& dim5,
2864 const Lagrange& dim6);
2865
2880template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2881 std::size_t PolyOrder2, std::size_t PolyOrder3,
2882 std::size_t PolyOrder4, std::size_t PolyOrder5,
2883 std::size_t PolyOrder6, class Tensor7Type>
2884constexpr Numeric interp(
2885 const Tensor7Type& yi,
2886 const FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1, PolyOrder2 + 1,
2887 PolyOrder3 + 1, PolyOrder4 + 1, PolyOrder5 + 1,
2888 PolyOrder6 + 1>& iw,
2889 const FixedLagrange<PolyOrder0>& dim0,
2890 const FixedLagrange<PolyOrder1>& dim1,
2891 const FixedLagrange<PolyOrder2>& dim2,
2892 const FixedLagrange<PolyOrder3>& dim3,
2893 const FixedLagrange<PolyOrder4>& dim4,
2894 const FixedLagrange<PolyOrder5>& dim5,
2895 const FixedLagrange<PolyOrder6>& dim6) {
2896 Numeric out(0.0);
2897 const Index I = yi.nlibraries();
2898 const Index J = yi.nvitrines();
2899 const Index K = yi.nshelves();
2900 const Index L = yi.nbooks();
2901 const Index M = yi.npages();
2902 const Index N = yi.nrows();
2903 const Index O = yi.ncols();
2904 for (Index i = 0; i < dim0.size(); i++)
2905 for (Index j = 0; j < dim1.size(); j++)
2906 for (Index k = 0; k < dim2.size(); k++)
2907 for (Index l = 0; l < dim3.size(); l++)
2908 for (Index m = 0; m < dim4.size(); m++)
2909 for (Index n = 0; n < dim5.size(); n++)
2910 for (Index o = 0; o < dim6.size(); o++)
2911 out += iw(i, j, k, l, m, n, o) *
2912 yi(cycler(i + dim0.pos, I), cycler(j + dim1.pos, J),
2913 cycler(k + dim2.pos, K), cycler(l + dim3.pos, L),
2914 cycler(m + dim4.pos, M), cycler(n + dim5.pos, N),
2915 cycler(o + dim6.pos, O));
2916 return out;
2917}
2918
2922
2938void reinterp(Tensor7View out, const ConstTensor7View& iy,
2939 const Grid<Tensor7, 7>& iw, const Array<Lagrange>& dim0,
2940 const Array<Lagrange>& dim1,
2941 const Array<Lagrange>& dim2,
2942 const Array<Lagrange>& dim3,
2943 const Array<Lagrange>& dim4,
2944 const Array<Lagrange>& dim5,
2945 const Array<Lagrange>& dim6);
2946
2962 const Array<Lagrange>& dim0,
2963 const Array<Lagrange>& dim1,
2964 const Array<Lagrange>& dim2,
2965 const Array<Lagrange>& dim3,
2966 const Array<Lagrange>& dim4,
2967 const Array<Lagrange>& dim5,
2968 const Array<Lagrange>& dim6);
2969
2984template <std::size_t PolyOrder0, std::size_t PolyOrder1,
2985 std::size_t PolyOrder2, std::size_t PolyOrder3,
2986 std::size_t PolyOrder4, std::size_t PolyOrder5,
2987 std::size_t PolyOrder6>
2989 const ConstTensor7View& iy,
2990 const Grid<FixedGrid<Numeric, PolyOrder0 + 1, PolyOrder1 + 1,
2991 PolyOrder2 + 1, PolyOrder3 + 1, PolyOrder4 + 1,
2992 PolyOrder5 + 1, PolyOrder6 + 1>,
2993 7>& iw,
2994 const Array<FixedLagrange<PolyOrder0>>& dim0,
2995 const Array<FixedLagrange<PolyOrder1>>& dim1,
2996 const Array<FixedLagrange<PolyOrder2>>& dim2,
2997 const Array<FixedLagrange<PolyOrder3>>& dim3,
2998 const Array<FixedLagrange<PolyOrder4>>& dim4,
2999 const Array<FixedLagrange<PolyOrder5>>& dim5,
3000 const Array<FixedLagrange<PolyOrder6>>& dim6) {
3001 Tensor7 out(dim0.size(), dim1.size(), dim2.size(), dim3.size(), dim4.size(),
3002 dim5.size(), dim6.size());
3003 for (std::size_t i = 0; i < dim0.size(); i++)
3004 for (std::size_t j = 0; j < dim1.size(); j++)
3005 for (std::size_t k = 0; k < dim2.size(); k++)
3006 for (std::size_t l = 0; l < dim3.size(); l++)
3007 for (std::size_t m = 0; m < dim4.size(); m++)
3008 for (std::size_t n = 0; n < dim5.size(); n++)
3009 for (std::size_t o = 0; o < dim6.size(); o++)
3010 out(i, j, k, l, m, n, o) =
3011 interp(iy, iw(i, j, k, l, m, n, o), dim0[i], dim1[j],
3012 dim2[k], dim3[l], dim4[m], dim5[n], dim6[o]);
3013 return out;
3014}
3015} // namespace Interpolation
3016
3019
3022
3025
3028
3031
3033template <std::size_t N>
3035
3036#endif // interpolation_lagrange_h
This can be used to make arrays out of anything.
Definition: array.h:107
A constant view of a Matrix.
Definition: matpackI.h:1014
A constant view of a Tensor3.
Definition: matpackIII.h:132
A constant view of a Tensor4.
Definition: matpackIV.h:133
A constant view of a Tensor5.
Definition: matpackV.h:143
A constant view of a Tensor6.
Definition: matpackVI.h:149
A constant view of a Tensor7.
Definition: matpackVII.h:147
A constant view of a Vector.
Definition: matpackI.h:489
Row-major fixed grid creation.
Definition: grids.h:132
Row-major grid creation.
Definition: grids.h:52
The MatrixView class.
Definition: matpackI.h:1125
The Matrix class.
Definition: matpackI.h:1225
The Tensor3View class.
Definition: matpackIII.h:239
The Tensor3 class.
Definition: matpackIII.h:339
The Tensor4View class.
Definition: matpackIV.h:284
The Tensor4 class.
Definition: matpackIV.h:421
The Tensor5View class.
Definition: matpackV.h:333
The Tensor5 class.
Definition: matpackV.h:506
The Tensor6View class.
Definition: matpackVI.h:621
The Tensor6 class.
Definition: matpackVI.h:1088
The Tensor7View class.
Definition: matpackVII.h:1286
The Tensor7 class.
Definition: matpackVII.h:2382
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
Constants of physical expressions as constexpr.
Helper macros for debugging.
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
GridType
#define abs(x)
#define dx
Eigen::Matrix< Numeric, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > MatrixType
Definition: matpackI.h:114
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
auto sind(T x) -> decltype(std::sin(deg2rad(x)))
Returns the sine of the deg2rad of the input.
Definition: constants.h:346
auto cosd(T x) -> decltype(std::cos(deg2rad(x)))
Returns the cosine of the deg2rad of the input.
Definition: constants.h:340
char Type type
constexpr Numeric cyclic_clamp(const Numeric x, const std::pair< Numeric, Numeric > xlim) noexcept
constexpr Index cycler(const Index n, const Index N) noexcept
constexpr Index pos_finder(const Index pos0, const Numeric x, const SortedVectorType &xi, const Index polyorder, const bool cyclic, const bool ascending, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr double dl_dval(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const Index i, const std::pair< Numeric, Numeric > cycle) noexcept
Numeric interp(const ConstVectorView &yi, const ConstVectorView &iw, const Lagrange &dim0)
Array< FixedLagrange< PolyOrder > > FixedLagrangeVector(const UnsortedVectorType &xs, const SortedVectorType &xi, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180})
constexpr Index IMIN(const Index a, const Index b) noexcept
Return the minimum of two integer numbers.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr bool is_ascending(const SortedVectorType &xi) noexcept
Checks whether the Sorted Vector is ascending or not by checking its first two elements.
void check_lagrange_interpolation(const SortedVectorType &xi, const Index polyorder=1, const std::pair< Numeric, Numeric > x={std::numeric_limits< Numeric >::infinity(), -std::numeric_limits< Numeric >::infinity()}, const Numeric extrapol=0.5, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180})
void reinterp(VectorView out, const ConstVectorView &iy, const Grid< Vector, 1 > &iw, const Array< Lagrange > &dim0)
constexpr Numeric dl(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const LagrangeVectorType &li, const Index j, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
Vector dinterpweights(const Lagrange &dim0)
constexpr bool full_cycle(const Numeric xo, const Numeric xn, const std::pair< Numeric, Numeric > xlim) noexcept
constexpr Numeric min_cyclic(const Numeric x, const std::pair< Numeric, Numeric > xlim) noexcept
ENUMCLASS(GridType, char, Standard, Cyclic, Log, Log10, Log2, SinDeg, SinRad, CosDeg, CosRad)
constexpr Index IMAX(const Index a, const Index b) noexcept
Return the maximum of two integer numbers.
constexpr Index start_pos_finder(const Numeric x, const SortedVectorType &xvec) noexcept
Vector interpweights(const Lagrange &dim0)
Array< Lagrange > LagrangeVector(const ConstVectorView &xs, const ConstVectorView &xi, const Index polyorder, const Numeric extrapol, const bool do_derivs, const GridType type, const std::pair< Numeric, Numeric > cycle)
Definition: nonstd.h:5
constexpr T abs(T x) noexcept
Definition: nonstd.h:13
#define a
#define b
#define N
Definition: rng.cc:164
#define M
Definition: rng.cc:165
static constexpr Index size() noexcept
static constexpr std::array< Numeric, PolyOrder+1 > dlx_finder(const Index p0, const Numeric x, const SortedVectorType &xi, const std::array< Numeric, PolyOrder+1 > &li, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
constexpr FixedLagrange & operator=(FixedLagrange &&l) noexcept
Enusre that the move operator exists.
friend std::ostream & operator<<(std::ostream &os, const FixedLagrange &l)
std::array< Numeric, PolyOrder+1 > dlx
static constexpr std::array< Numeric, PolyOrder+1 > lx_finder(const Index p0, const Numeric x, const SortedVectorType &xi, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
constexpr FixedLagrange(FixedLagrange &&l) noexcept
Enusre that the move constructor exists.
std::array< Numeric, PolyOrder+1 > lx
constexpr FixedLagrange(const Index p0, const Numeric x, const SortedVectorType &xi, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
A Lagrange interpolation computer.
static Array< Numeric > lx_finder(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
Lagrange() noexcept
Default constructor for zero-length elements.
Lagrange(const Index pos0, const Numeric x, const SortedVectorType &xi, const Index polyorder=1, const bool do_derivs=false, const GridType type=GridType::Standard, const std::pair< Numeric, Numeric > cycle={-180, 180}) noexcept
friend std::ostream & operator<<(std::ostream &os, const Lagrange &l)
static Array< Numeric > dlx_finder(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Array< Numeric > &li, const GridType type, const std::pair< Numeric, Numeric > cycle) noexcept
Lagrange & operator=(Lagrange &&l) noexcept
Ensure that the move operator exists.
Lagrange(Lagrange &&l) noexcept
Ensure that the move constructor exists.
Lagrange(const Lagrange &l) noexcept
Ensure that the copy constructor exists.
Lagrange & operator=(const Lagrange &l) noexcept
Ensure that the copy operator exists.
Index size() const noexcept