1#ifndef INTERPOLATION_IMPL
2#error Please, include "MathUtils/interpolation/interpolation.h"
5#include "AlexandriaKernel/index_sequence.h"
6#include "MathUtils/interpolation/GridInterpolation.h"
11template <std::size_t, typename Seq>
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
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>
29 template <std::size_t>
31 using type = std::vector<double>;
34 InterpNAdapter(const Coordinates<N>& grid, const NdArray::NdArray<double>& values, InterpolationType type,
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";
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)...);
47 void operator()(const typename Vectors<Is>::type&..., std::vector<double>&) const override {
48 throw Elements::Exception() << "Not implemented";
51 std::unique_ptr<NAryFunction<N>> clone() const override {
52 return Euclid::make_unique<InterpNAdapter>(*this);
55 InterpNAdapter(const InterpNAdapter&) = default;
58 InterpN<typename Doubles<Is>::type...> m_interpn;
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);
67} // namespace MathUtils