Line data Source code
1 0 : // Distributed under the MIT License.
2 : // See LICENSE.txt for details.
3 :
4 : #pragma once
5 :
6 : #include <algorithm>
7 : #include <array>
8 : #include <cmath>
9 : #include <memory>
10 :
11 : #include "DataStructures/Index.hpp"
12 : #include "Utilities/ConstantExpressions.hpp"
13 : #include "Utilities/ErrorHandling/Assert.hpp"
14 : #include "Utilities/GenerateInstantiations.hpp"
15 : #include "Utilities/Gsl.hpp"
16 : #include "Utilities/Requires.hpp"
17 : #include "Utilities/TMPL.hpp"
18 :
19 : namespace intrp {
20 :
21 : /*!
22 : * \brief Performs linear interpolation in arbitrary dimensions.
23 : * The class is non-owning and expects a C-ordered array, (n, x, y, z).
24 : * The variable index, n, varies fastest in memory.
25 : * Note that this class is intentionally non-pupable.
26 : *
27 : * \tparam Dimension dimensionality of the table
28 : * \tparam NumberOfVariables number of variables stored in the table
29 : * \tparam UniformSpacing indicated whether the table has uniform _spacing or
30 : * not. This is useful for performance reasons in finding the appropriate table
31 : * index.
32 : */
33 :
34 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
35 1 : class MultiLinearSpanInterpolation {
36 :
37 : public:
38 : template <size_t ThisDimension>
39 0 : struct Weight {
40 0 : std::array<double, two_to_the(ThisDimension)> weights;
41 0 : Index<two_to_the(ThisDimension)> index;
42 : };
43 :
44 0 : size_t extents(const size_t which_dimension) const {
45 : return number_of_points_[which_dimension];
46 : }
47 :
48 0 : void extrapolate_above_data(const size_t which_dimension, const bool value) {
49 : allow_extrapolation_abov_data_[which_dimension] = value;
50 : }
51 :
52 0 : void extrapolate_below_data(const size_t which_dimension, const bool value) {
53 : allow_extrapolation_below_data_[which_dimension] = value;
54 : }
55 :
56 : /// Compute interpolation weights for 1D tables
57 1 : Weight<1> get_weights(const double x1) const;
58 :
59 : /// Compute interpolation weights for 2D tables
60 1 : Weight<2> get_weights(const double x1, const double x2) const;
61 :
62 : /// Compute interpolation weights for 3D tables
63 1 : Weight<3> get_weights(const double x1, const double x2,
64 : const double x3) const;
65 :
66 0 : double interpolate(const Weight<Dimension>& weights,
67 : const size_t which_variable = 0) const {
68 : double result = 0;
69 :
70 : for (size_t nn = 0; nn < weights.weights.size(); ++nn) {
71 : result += weights.weights[nn] *
72 : y_[which_variable + NumberOfVariables * weights.index[nn]];
73 : }
74 :
75 : return result;
76 : }
77 :
78 : template <size_t... variables_to_interpolate>
79 0 : std::array<double, sizeof...(variables_to_interpolate)> interpolate(
80 : const Weight<Dimension>& weights) const {
81 : static_assert(sizeof...(variables_to_interpolate) <= NumberOfVariables,
82 : "You are trying to interpolate more variables than this "
83 : "container holds.");
84 :
85 : return std::array<double, sizeof...(variables_to_interpolate)>{
86 : interpolate(weights, variables_to_interpolate)...};
87 : }
88 :
89 : template <size_t NumberOfVariablesToInterpolate, typename... T,
90 : Requires<(std::is_floating_point_v<typename std::remove_cv_t<T>> and
91 : ...)> = nullptr>
92 0 : std::array<double, NumberOfVariablesToInterpolate> interpolate(
93 : std::array<size_t, NumberOfVariablesToInterpolate>&
94 : variables_to_interpolate,
95 : const T&... target_points) const {
96 : static_assert(
97 : sizeof...(T) == Dimension,
98 : "You need to provide the correct number of interpolation points");
99 : static_assert(NumberOfVariablesToInterpolate <= NumberOfVariables,
100 : "You are trying to interpolate more variables than this "
101 : "container holds.");
102 :
103 : auto weights = get_weights(target_points...);
104 :
105 : std::array<double, NumberOfVariablesToInterpolate> interpolatedValues;
106 :
107 : for (size_t nn = 0; nn < NumberOfVariablesToInterpolate; ++nn) {
108 : interpolatedValues[nn] =
109 : interpolate(weights, variables_to_interpolate[nn]);
110 : }
111 :
112 : return interpolatedValues;
113 : }
114 :
115 : template <size_t NumberOfVariablesToInterpolate, size_t... I>
116 0 : std::array<double, NumberOfVariablesToInterpolate> interpolate(
117 : std::array<size_t, NumberOfVariablesToInterpolate>&
118 : variables_to_interpolate,
119 : std::array<double, Dimension>& target_points,
120 : std::index_sequence<I...>) const {
121 : return interpolate(variables_to_interpolate, target_points[I]...);
122 : }
123 :
124 : template <size_t NumberOfVariablesToInterpolate>
125 0 : std::array<double, NumberOfVariablesToInterpolate> interpolate(
126 : std::array<size_t, NumberOfVariablesToInterpolate>&
127 : variables_to_interpolate,
128 : std::array<double, Dimension>& target_points) const {
129 : return interpolate(variables_to_interpolate, target_points,
130 : std::make_index_sequence<Dimension>{});
131 : }
132 :
133 0 : MultiLinearSpanInterpolation() = default;
134 :
135 0 : MultiLinearSpanInterpolation(
136 : std::array<gsl::span<const double>, Dimension> x_,
137 : gsl::span<const double> y_, Index<Dimension> number_of_points__);
138 :
139 0 : double lower_bound(const size_t which_dimension) const {
140 : return x_[which_dimension][0];
141 : }
142 :
143 0 : double upper_bound(const size_t which_dimension) const {
144 : return x_[which_dimension][number_of_points_[which_dimension] - 1];
145 : }
146 :
147 : private:
148 : /// Allow extrapolation above table bounds.
149 1 : std::array<bool, Dimension> allow_extrapolation_below_data_;
150 : /// Allow extrapolation below table bounds.
151 1 : std::array<bool, Dimension> allow_extrapolation_abov_data_;
152 : /// Inverse pacing of the table. Only used for uniform grids
153 1 : std::array<double, Dimension> inverse_spacing_;
154 : /// Spacing of the table. Only used for uniform grids
155 1 : std::array<double, Dimension> spacing_;
156 : /// Number of points per dimension
157 1 : Index<Dimension> number_of_points_;
158 :
159 0 : using DataPointer = gsl::span<const double>;
160 : /// X values of the table. Only used if allocated
161 1 : std::array<DataPointer, Dimension> x_;
162 : /// Y values of the table. Only used if allocated
163 1 : DataPointer y_;
164 :
165 : /// Determine relative index bracketing for non-uniform _spacing
166 1 : size_t find_index_general(const size_t which_dimension,
167 : const double& target_points) const;
168 :
169 : /// Determine relative index bracketing for non-uniform _spacing
170 1 : size_t find_index_uniform(const size_t which_dimension,
171 : const double& target_points) const;
172 :
173 0 : size_t find_index(const size_t which_dimension,
174 : const double& target_points) const {
175 : if constexpr (UniformSpacing) {
176 : return find_index_uniform(which_dimension, target_points);
177 : } else {
178 : return find_index_general(which_dimension, target_points);
179 : }
180 : }
181 : };
182 :
183 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
184 : size_t MultiLinearSpanInterpolation<
185 : Dimension, NumberOfVariables,
186 : UniformSpacing>::find_index_general(const size_t which_dimension,
187 : const double& target_points) const {
188 : size_t lower_index_bound = 0;
189 : size_t upper_index_bound = number_of_points_[which_dimension] - 1;
190 :
191 : ASSERT((target_points > x_[which_dimension][lower_index_bound]) or
192 : allow_extrapolation_below_data_[which_dimension],
193 : "Interpolation exceeds lower table bounds. \nwhich_dimension: "
194 : << which_dimension
195 : << "\nnumber of points: " << number_of_points_[which_dimension]
196 : << "\ntarget point: " << target_points);
197 : ASSERT((target_points < x_[which_dimension][upper_index_bound]) or
198 : allow_extrapolation_abov_data_[which_dimension],
199 : "Interpolation exceeds upper table bounds. \nwhich_dimension: "
200 : << which_dimension
201 : << "\nnumber of points: " << number_of_points_[which_dimension]
202 : << "\ntarget point: " << target_points);
203 :
204 : while (upper_index_bound > 1 + lower_index_bound) {
205 : size_t current_index =
206 : lower_index_bound + (upper_index_bound - lower_index_bound) / 2;
207 : if (target_points < x_[which_dimension][current_index]) {
208 : upper_index_bound = current_index;
209 : } else {
210 : lower_index_bound = current_index;
211 : }
212 : }
213 : return lower_index_bound;
214 : }
215 :
216 : /// Determine relative index bracketing for non-uniform _spacing
217 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
218 : size_t MultiLinearSpanInterpolation<
219 : Dimension, NumberOfVariables,
220 : UniformSpacing>::find_index_uniform(const size_t which_dimension,
221 : const double& target_points) const {
222 : // Compute coordinate relative to lowest table bound
223 : const auto relative_coordinate = (target_points - x_[which_dimension][0]);
224 :
225 : auto current_index = static_cast<size_t>(relative_coordinate *
226 : inverse_spacing_[which_dimension]);
227 :
228 : // We are exceeding the table bounds:
229 : // Use linear extrapolation based of the lowest
230 : // two points in the table
231 : ASSERT(allow_extrapolation_below_data_[which_dimension] or
232 : UNLIKELY(relative_coordinate >= 0.),
233 : "Interpolation exceeds lower table bounds.\nwhich_dimension: "
234 : << which_dimension << "\ncurrent_index: " << current_index
235 : << "\nnumber of points: " << number_of_points_[which_dimension]
236 : << "\ntarget point: " << target_points
237 : << "\nrelative coordinate: " << relative_coordinate);
238 :
239 : // We are exceeding the table bounds:
240 : // Use linear extrapolation based of the highest
241 : // two points in the table
242 :
243 : ASSERT(allow_extrapolation_abov_data_[which_dimension] or
244 : UNLIKELY(current_index + 1 < number_of_points_[which_dimension]),
245 : "Interpolation exceeds upper table bounds.\nwhich_dimension: "
246 : << which_dimension << "\ncurrent_index: " << current_index
247 : << "\nnumber of points: " << number_of_points_[which_dimension]
248 : << "\ntarget point: " << target_points
249 : << "\nrelative coordinate: " << relative_coordinate);
250 :
251 : // Enforce index ranges
252 : current_index = std::min(number_of_points_[which_dimension] - 2,
253 : std::max(0ul, current_index));
254 :
255 : return current_index;
256 : }
257 :
258 : /// Compute interpolation weights for 1D tables
259 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
260 : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
261 : UniformSpacing>::get_weights(const double x1)
262 : const -> Weight<1> {
263 : // Relative normalized coordinate in x-direction
264 : double xx;
265 :
266 : Weight<1> weights;
267 :
268 : auto index = Index<1>(find_index(0, x1));
269 :
270 : if constexpr (UniformSpacing) {
271 : xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
272 : } else {
273 : xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
274 : }
275 :
276 : // Note: first index varies fastest
277 : weights.weights = std::array<double, 2>({
278 : (1. - xx), // 0
279 : xx, // 1
280 : });
281 :
282 : // Compute indices
283 :
284 : weights.index[0] = index[0];
285 : weights.index[1] = index[0] + 1;
286 :
287 : return weights;
288 : }
289 :
290 : /// Compute interpolation weights for 2D tables
291 :
292 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
293 : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
294 : UniformSpacing>::get_weights(const double x1,
295 : const double x2)
296 : const -> Weight<2> {
297 : Weight<2> weights;
298 :
299 : auto index = Index<2>(find_index(0, x1), find_index(1, x2));
300 :
301 : // Relative normalized coordinate in x,y-direction
302 : double xx, yy;
303 : if constexpr (UniformSpacing) {
304 : xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
305 : yy = (x2 - x_[1][index[1]]) * inverse_spacing_[1];
306 : } else {
307 : xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
308 : yy = (x2 - x_[1][index[1]]) / (x_[1][index[1] + 1] - x_[1][index[1]]);
309 : }
310 :
311 : // Note: first index varies fastest
312 : weights.weights = std::array<double, 4>({
313 : (1. - xx) * (1. - yy), // 00
314 : xx * (1. - yy), // 10
315 : (1. - xx) * yy, // 01
316 : xx * yy, // 11
317 : });
318 :
319 : // Compute indices
320 : //
321 : for (size_t j = 0; j < 2; ++j) {
322 : for (size_t i = 0; i < 2; ++i) {
323 : auto tmp_index = Index<2>(index[0] + i, index[1] + j);
324 : weights.index[i + 2 * j] = collapsed_index(tmp_index, number_of_points_);
325 : }
326 : }
327 :
328 : return weights;
329 : }
330 :
331 : /// Compute interpolation weights for 3D tables
332 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
333 : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
334 : UniformSpacing>::get_weights(const double x1,
335 : const double x2,
336 : const double x3)
337 : const -> Weight<3> {
338 : // Relative normalized coordinate in x,y,z-direction
339 : double xx, yy, zz;
340 :
341 : Weight<3> weights;
342 :
343 : auto index =
344 : Index<3>(find_index(0, x1), find_index(1, x2), find_index(2, x3));
345 :
346 : if constexpr (UniformSpacing) {
347 : xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
348 : yy = (x2 - x_[1][index[1]]) * inverse_spacing_[1];
349 : zz = (x3 - x_[2][index[2]]) * inverse_spacing_[2];
350 : } else {
351 : xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
352 : yy = (x2 - x_[1][index[1]]) / (x_[1][index[1] + 1] - x_[1][index[1]]);
353 : zz = (x3 - x_[2][index[2]]) / (x_[2][index[2] + 1] - x_[2][index[2]]);
354 : }
355 :
356 : // Note: first index varies fastest
357 : weights.weights = std::array<double, 8>({
358 : (1. - xx) * (1. - yy) * (1. - zz), // 000
359 : xx * (1. - yy) * (1. - zz), // 100
360 : (1. - xx) * yy * (1. - zz), // 010
361 : xx * yy * (1. - zz), // 110
362 : (1. - xx) * (1. - yy) * zz, // 001
363 : xx * (1. - yy) * zz, // 101
364 : (1. - xx) * yy * zz, // 011
365 : xx * yy * zz, // 111
366 : });
367 :
368 : // Compute indices
369 : //
370 : for (size_t k = 0; k < 2; ++k) {
371 : for (size_t j = 0; j < 2; ++j) {
372 : for (size_t i = 0; i < 2; ++i) {
373 : auto tmp_index = Index<3>(index[0] + i, index[1] + j, index[2] + k);
374 : weights.index[i + 2 * (j + 2 * k)] =
375 : collapsed_index(tmp_index, number_of_points_);
376 : }
377 : }
378 : }
379 :
380 : return weights;
381 : }
382 :
383 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
384 : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, UniformSpacing>::
385 : MultiLinearSpanInterpolation(
386 : std::array<gsl::span<double const>, Dimension> x,
387 : gsl::span<double const> y, Index<Dimension> number_of_points)
388 : : number_of_points_(number_of_points), x_(x), y_(y) {
389 : for (size_t i = 0; i < Dimension; ++i) {
390 : spacing_[i] = x_[i][1] - x_[i][0];
391 : inverse_spacing_[i] = 1. / spacing_[i];
392 : allow_extrapolation_below_data_[i] = false;
393 : allow_extrapolation_abov_data_[i] = false;
394 : }
395 : }
396 :
397 : /// Multilinear span interpolation with uniform grid spacing
398 : template <size_t Dimension, size_t NumberOfVariables>
399 1 : using UniformMultiLinearSpanInterpolation =
400 : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, true>;
401 :
402 : /// Multilinear span interpolation with non-uniform grid spacing
403 : template <size_t Dimension, size_t NumberOfVariables>
404 1 : using GeneralMultiLinearSpanInterpolation =
405 : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, false>;
406 :
407 : } // namespace intrp
|
