Alexandria 2.31.0
SDC-CH common library for the Euclid project
Loading...
Searching...
No Matches
interpolation.icpp
Go to the documentation of this file.
1#ifndef INTERPOLATION_IMPL
2#error Please, include "MathUtils/interpolation/interpolation.h"
3#endif
4
5#include "AlexandriaKernel/index_sequence.h"
6#include "MathUtils/interpolation/GridInterpolation.h"
7
8namespace Euclid {
9namespace MathUtils {
10
11template <std::size_t, typename Seq>
12struct InterpNAdapter;
13
14/**
15 * @description
16 * GridInterpolation expects the data to follow GridContainer memory layout, but
17 * originally interpn is expected to follow a numpy memory layout.
18 * We use template dirty tricks to re-use GridInterpolation, transposing the axis at creation time,
19 * and the arguments at interpolation time
20 */
21template <std::size_t N, std::size_t... Is>
22struct InterpNAdapter<N, _index_sequence<Is...>> : NAryFunction<N> {
23 // GCC 4.8 needs this hack
24 template <std::size_t>
25 struct Doubles {
26 using type = double;
27 };
28
29 template <std::size_t>
30 struct Vectors {
31 using type = std::vector<double>;
32 };
33
34 InterpNAdapter(const Coordinates<N>& grid, const NdArray::NdArray<double>& values, InterpolationType type,
35 bool extrapolate)
36 : m_interpn(std::tuple<std::vector<typename Doubles<Is>::type>...>{grid[N - Is - 1]...}, values, extrapolate) {
37 if (type != InterpolationType::LINEAR) {
38 throw InterpolationException() << "Only linear interpolation is supported for N-dimensional grids";
39 }
40 }
41
42 double operator()(typename Doubles<Is>::type... xn) const override {
43 auto as_tuple = std::make_tuple(xn...);
44 return m_interpn(std::get<N - Is - 1>(as_tuple)...);
45 }
46
47 void operator()(const typename Vectors<Is>::type&..., std::vector<double>&) const override {
48 throw Elements::Exception() << "Not implemented";
49 }
50
51 std::unique_ptr<NAryFunction<N>> clone() const override {
52 return Euclid::make_unique<InterpNAdapter>(*this);
53 }
54
55 InterpNAdapter(const InterpNAdapter&) = default;
56
57private:
58 InterpN<typename Doubles<Is>::type...> m_interpn;
59};
60
61template <std::size_t N>
62std::unique_ptr<NAryFunction<N>> interpn(const Coordinates<N>& grid, const NdArray::NdArray<double>& values,
63 InterpolationType type, bool extrapolate) {
64 return Euclid::make_unique<InterpNAdapter<N, _make_index_sequence<N>>>(grid, values, type, extrapolate);
65}
66
67} // namespace MathUtils
68} // namespace Euclid