Alexandria  2.22.0
Please provide a description of the project.
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 
8 namespace Euclid {
9 namespace MathUtils {
10 
11 template <std::size_t, typename Seq>
12 struct 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  */
21 template <std::size_t N, std::size_t... Is>
22 struct 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 
57 private:
58  InterpN<typename Doubles<Is>::type...> m_interpn;
59 };
60 
61 template <std::size_t N>
62 std::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