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 <concepts>
10 : #include <memory>
11 :
12 : #include "DataStructures/Index.hpp"
13 : #include "Utilities/ConstantExpressions.hpp"
14 : #include "Utilities/ErrorHandling/Assert.hpp"
15 : #include "Utilities/GenerateInstantiations.hpp"
16 : #include "Utilities/Gsl.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, std::floating_point... T>
90 0 : std::array<double, NumberOfVariablesToInterpolate> interpolate(
91 : std::array<size_t, NumberOfVariablesToInterpolate>&
92 : variables_to_interpolate,
93 : const T&... target_points) const {
94 : static_assert(
95 : sizeof...(T) == Dimension,
96 : "You need to provide the correct number of interpolation points");
97 : static_assert(NumberOfVariablesToInterpolate <= NumberOfVariables,
98 : "You are trying to interpolate more variables than this "
99 : "container holds.");
100 :
101 : auto weights = get_weights(target_points...);
102 :
103 : std::array<double, NumberOfVariablesToInterpolate> interpolatedValues;
104 :
105 : for (size_t nn = 0; nn < NumberOfVariablesToInterpolate; ++nn) {
106 : interpolatedValues[nn] =
107 : interpolate(weights, variables_to_interpolate[nn]);
108 : }
109 :
110 : return interpolatedValues;
111 : }
112 :
113 : template <size_t NumberOfVariablesToInterpolate, size_t... I>
114 0 : std::array<double, NumberOfVariablesToInterpolate> interpolate(
115 : std::array<size_t, NumberOfVariablesToInterpolate>&
116 : variables_to_interpolate,
117 : std::array<double, Dimension>& target_points,
118 : std::index_sequence<I...>) const {
119 : return interpolate(variables_to_interpolate, target_points[I]...);
120 : }
121 :
122 : template <size_t NumberOfVariablesToInterpolate>
123 0 : std::array<double, NumberOfVariablesToInterpolate> interpolate(
124 : std::array<size_t, NumberOfVariablesToInterpolate>&
125 : variables_to_interpolate,
126 : std::array<double, Dimension>& target_points) const {
127 : return interpolate(variables_to_interpolate, target_points,
128 : std::make_index_sequence<Dimension>{});
129 : }
130 :
131 0 : MultiLinearSpanInterpolation() = default;
132 :
133 0 : MultiLinearSpanInterpolation(std::array<gsl::span<const double>, Dimension> x,
134 : gsl::span<const double> y,
135 : Index<Dimension> number_of_points);
136 :
137 0 : double lower_bound(const size_t which_dimension) const {
138 : return x_[which_dimension][0];
139 : }
140 :
141 0 : double upper_bound(const size_t which_dimension) const {
142 : return x_[which_dimension][number_of_points_[which_dimension] - 1];
143 : }
144 :
145 : private:
146 : /// Allow extrapolation above table bounds.
147 1 : std::array<bool, Dimension> allow_extrapolation_below_data_;
148 : /// Allow extrapolation below table bounds.
149 1 : std::array<bool, Dimension> allow_extrapolation_abov_data_;
150 : /// Inverse pacing of the table. Only used for uniform grids
151 1 : std::array<double, Dimension> inverse_spacing_;
152 : /// Spacing of the table. Only used for uniform grids
153 1 : std::array<double, Dimension> spacing_;
154 : /// Number of points per dimension
155 1 : Index<Dimension> number_of_points_;
156 :
157 0 : using DataPointer = gsl::span<const double>;
158 : /// X values of the table. Only used if allocated
159 1 : std::array<DataPointer, Dimension> x_;
160 : /// Y values of the table. Only used if allocated
161 1 : DataPointer y_;
162 :
163 : /// Determine relative index bracketing for non-uniform _spacing
164 1 : size_t find_index_general(const size_t which_dimension,
165 : const double& target_points) const;
166 :
167 : /// Determine relative index bracketing for non-uniform _spacing
168 1 : size_t find_index_uniform(const size_t which_dimension,
169 : const double& target_points) const;
170 :
171 0 : size_t find_index(const size_t which_dimension,
172 : const double& target_points) const {
173 : if constexpr (UniformSpacing) {
174 : return find_index_uniform(which_dimension, target_points);
175 : } else {
176 : return find_index_general(which_dimension, target_points);
177 : }
178 : }
179 : };
180 :
181 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
182 : size_t MultiLinearSpanInterpolation<
183 : Dimension, NumberOfVariables,
184 : UniformSpacing>::find_index_general(const size_t which_dimension,
185 : const double& target_points) const {
186 : size_t lower_index_bound = 0;
187 : size_t upper_index_bound = number_of_points_[which_dimension] - 1;
188 :
189 : ASSERT((target_points > x_[which_dimension][lower_index_bound]) or
190 : allow_extrapolation_below_data_[which_dimension],
191 : "Interpolation exceeds lower table bounds. \nwhich_dimension: "
192 : << which_dimension
193 : << "\nnumber of points: " << number_of_points_[which_dimension]
194 : << "\ntarget point: " << target_points);
195 : ASSERT((target_points < x_[which_dimension][upper_index_bound]) or
196 : allow_extrapolation_abov_data_[which_dimension],
197 : "Interpolation exceeds upper table bounds. \nwhich_dimension: "
198 : << which_dimension
199 : << "\nnumber of points: " << number_of_points_[which_dimension]
200 : << "\ntarget point: " << target_points);
201 :
202 : while (upper_index_bound > 1 + lower_index_bound) {
203 : size_t current_index =
204 : lower_index_bound + (upper_index_bound - lower_index_bound) / 2;
205 : if (target_points < x_[which_dimension][current_index]) {
206 : upper_index_bound = current_index;
207 : } else {
208 : lower_index_bound = current_index;
209 : }
210 : }
211 : return lower_index_bound;
212 : }
213 :
214 : /// Determine relative index bracketing for non-uniform _spacing
215 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
216 : size_t MultiLinearSpanInterpolation<
217 : Dimension, NumberOfVariables,
218 : UniformSpacing>::find_index_uniform(const size_t which_dimension,
219 : const double& target_points) const {
220 : // Compute coordinate relative to lowest table bound
221 : const auto relative_coordinate = (target_points - x_[which_dimension][0]);
222 :
223 : auto current_index = static_cast<size_t>(relative_coordinate *
224 : inverse_spacing_[which_dimension]);
225 :
226 : // We are exceeding the table bounds:
227 : // Use linear extrapolation based of the lowest
228 : // two points in the table
229 : ASSERT(allow_extrapolation_below_data_[which_dimension] or
230 : UNLIKELY(relative_coordinate >= 0.),
231 : "Interpolation exceeds lower table bounds.\nwhich_dimension: "
232 : << which_dimension << "\ncurrent_index: " << current_index
233 : << "\nnumber of points: " << number_of_points_[which_dimension]
234 : << "\ntarget point: " << target_points
235 : << "\nrelative coordinate: " << relative_coordinate
236 : << "\nspacing: " << inverse_spacing_[which_dimension]);
237 :
238 : // Do not trigger errors for roundoff errors in index calculation
239 : if (UNLIKELY(current_index + 1 == number_of_points_[which_dimension])) {
240 : if (relative_coordinate * inverse_spacing_[which_dimension] -
241 : static_cast<double>(current_index) <
242 : 1.e-13) {
243 : current_index -= 1;
244 : }
245 : }
246 :
247 : // We are exceeding the table bounds:
248 : // Use linear extrapolation based of the highest
249 : // two points in the table
250 :
251 : ASSERT(
252 : allow_extrapolation_abov_data_[which_dimension] or
253 : UNLIKELY(current_index + 1 < number_of_points_[which_dimension]),
254 : "Interpolation exceeds upper table bounds.\nwhich_dimension: "
255 : << which_dimension << "\ncurrent_index: " << current_index
256 : << "\nnumber of points: " << number_of_points_[which_dimension]
257 : << "\ntarget point: " << target_points
258 : << "\nrelative coordinate: " << relative_coordinate << "\nposition: "
259 : << relative_coordinate * inverse_spacing_[which_dimension] -
260 : static_cast<double>(number_of_points_[which_dimension] - 1));
261 :
262 : // Enforce index ranges
263 : current_index = std::min(number_of_points_[which_dimension] - 2,
264 : std::max(0ul, current_index));
265 :
266 : return current_index;
267 : }
268 :
269 : /// Compute interpolation weights for 1D tables
270 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
271 : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
272 : UniformSpacing>::get_weights(const double x1)
273 : const -> Weight<1> {
274 : // Relative normalized coordinate in x-direction
275 : double xx;
276 :
277 : Weight<1> weights;
278 :
279 : auto index = Index<1>(find_index(0, x1));
280 :
281 : if constexpr (UniformSpacing) {
282 : xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
283 : } else {
284 : xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
285 : }
286 :
287 : // Note: first index varies fastest
288 : weights.weights = std::array<double, 2>({
289 : (1. - xx), // 0
290 : xx, // 1
291 : });
292 :
293 : // Compute indices
294 :
295 : weights.index[0] = index[0];
296 : weights.index[1] = index[0] + 1;
297 :
298 : return weights;
299 : }
300 :
301 : /// Compute interpolation weights for 2D tables
302 :
303 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
304 : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
305 : UniformSpacing>::get_weights(const double x1,
306 : const double x2)
307 : const -> Weight<2> {
308 : Weight<2> weights;
309 :
310 : auto index = Index<2>(find_index(0, x1), find_index(1, x2));
311 :
312 : // Relative normalized coordinate in x,y-direction
313 : double xx, yy;
314 : if constexpr (UniformSpacing) {
315 : xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
316 : yy = (x2 - x_[1][index[1]]) * inverse_spacing_[1];
317 : } else {
318 : xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
319 : yy = (x2 - x_[1][index[1]]) / (x_[1][index[1] + 1] - x_[1][index[1]]);
320 : }
321 :
322 : // Note: first index varies fastest
323 : weights.weights = std::array<double, 4>({
324 : (1. - xx) * (1. - yy), // 00
325 : xx * (1. - yy), // 10
326 : (1. - xx) * yy, // 01
327 : xx * yy, // 11
328 : });
329 :
330 : // Compute indices
331 : //
332 : for (size_t j = 0; j < 2; ++j) {
333 : for (size_t i = 0; i < 2; ++i) {
334 : auto tmp_index = Index<2>(index[0] + i, index[1] + j);
335 : weights.index[i + 2 * j] = collapsed_index(tmp_index, number_of_points_);
336 : }
337 : }
338 :
339 : return weights;
340 : }
341 :
342 : /// Compute interpolation weights for 3D tables
343 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
344 : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
345 : UniformSpacing>::get_weights(const double x1,
346 : const double x2,
347 : const double x3)
348 : const -> Weight<3> {
349 : // Relative normalized coordinate in x,y,z-direction
350 : double xx, yy, zz;
351 :
352 : Weight<3> weights;
353 :
354 : auto index =
355 : Index<3>(find_index(0, x1), find_index(1, x2), find_index(2, x3));
356 :
357 : if constexpr (UniformSpacing) {
358 : xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
359 : yy = (x2 - x_[1][index[1]]) * inverse_spacing_[1];
360 : zz = (x3 - x_[2][index[2]]) * inverse_spacing_[2];
361 : } else {
362 : xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
363 : yy = (x2 - x_[1][index[1]]) / (x_[1][index[1] + 1] - x_[1][index[1]]);
364 : zz = (x3 - x_[2][index[2]]) / (x_[2][index[2] + 1] - x_[2][index[2]]);
365 : }
366 :
367 : // Note: first index varies fastest
368 : weights.weights = std::array<double, 8>({
369 : (1. - xx) * (1. - yy) * (1. - zz), // 000
370 : xx * (1. - yy) * (1. - zz), // 100
371 : (1. - xx) * yy * (1. - zz), // 010
372 : xx * yy * (1. - zz), // 110
373 : (1. - xx) * (1. - yy) * zz, // 001
374 : xx * (1. - yy) * zz, // 101
375 : (1. - xx) * yy * zz, // 011
376 : xx * yy * zz, // 111
377 : });
378 :
379 : // Compute indices
380 : //
381 : for (size_t k = 0; k < 2; ++k) {
382 : for (size_t j = 0; j < 2; ++j) {
383 : for (size_t i = 0; i < 2; ++i) {
384 : auto tmp_index = Index<3>(index[0] + i, index[1] + j, index[2] + k);
385 : weights.index[i + 2 * (j + 2 * k)] =
386 : collapsed_index(tmp_index, number_of_points_);
387 : }
388 : }
389 : }
390 :
391 : return weights;
392 : }
393 :
394 : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
395 : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, UniformSpacing>::
396 : MultiLinearSpanInterpolation(
397 : std::array<gsl::span<const double>, Dimension> x,
398 : gsl::span<const double> y, Index<Dimension> number_of_points)
399 : : number_of_points_(number_of_points), x_(x), y_(y) {
400 : for (size_t i = 0; i < Dimension; ++i) {
401 : spacing_[i] = x_[i][1] - x_[i][0];
402 : inverse_spacing_[i] = 1. / spacing_[i];
403 : allow_extrapolation_below_data_[i] = false;
404 : allow_extrapolation_abov_data_[i] = false;
405 : }
406 : }
407 :
408 : /// Multilinear span interpolation with uniform grid spacing
409 : template <size_t Dimension, size_t NumberOfVariables>
410 1 : using UniformMultiLinearSpanInterpolation =
411 : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, true>;
412 :
413 : /// Multilinear span interpolation with non-uniform grid spacing
414 : template <size_t Dimension, size_t NumberOfVariables>
415 1 : using GeneralMultiLinearSpanInterpolation =
416 : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, false>;
417 :
418 : } // namespace intrp
|
