1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
use super::*;
use strategy::*;
use itertools::Itertools;
pub fn get_index_permutations(shape: &[usize]) -> Vec<Vec<usize>> {
if shape.is_empty() {
return vec![vec![]];
}
shape
.iter()
.map(|&len| 0..len)
.multi_cartesian_product()
.collect()
}
impl<D> StrategyND<D> for Linear
where
D: Data + RawDataClone + Clone,
D::Elem: Num + PartialOrd + Copy + Debug,
{
fn interpolate(
&self,
data: &InterpDataND<D>,
point: &[D::Elem],
) -> Result<D::Elem, InterpolateError> {
// Dimensionality
let mut n = data.values.ndim();
// Point can share up to N values of a grid point, which reduces the problem dimensionality
// i.e. the point shares one of three values of a 3-D grid point, then the interpolation becomes 2-D at that slice
// or if the point shares two of three values of a 3-D grid point, then the interpolation becomes 1-D
let mut point = point.to_owned();
let mut grid: Vec<_> = data.grid.iter().map(|arr| arr.view()).collect();
let mut values_view = data.values.view();
for dim in (0..n).rev() {
// Range is reversed so that removal doesn't affect indexing
if let Some(pos) = grid[dim]
.iter()
.position(|&grid_point| grid_point == point[dim])
{
point.remove(dim);
grid.remove(dim);
values_view.index_axis_inplace(Axis(dim), pos);
}
}
if values_view.len() == 1 {
// Supplied point is coincident with a grid point, so just return the value
return Ok(values_view.first().copied().unwrap());
}
// Simplified dimensionality
n = values_view.ndim();
// Extract the lower and upper indices for each dimension,
// as well as the fraction of how far the supplied point is between the surrounding grid points
let mut lower_idxs = Vec::with_capacity(n);
let mut interp_diffs = Vec::with_capacity(n);
for dim in 0..n {
// Extrapolation is checked previously in Interpolator::interpolate,
// meaning by now, point is within grid bounds or extrapolation is enabled
let lower_idx = if &point[dim] < grid[dim].first().unwrap() {
0
} else if &point[dim] > grid[dim].last().unwrap() {
grid[dim].len() - 2
} else {
find_nearest_index(grid[dim].view(), &point[dim])
};
let interp_diff = (point[dim] - grid[dim][lower_idx])
/ (grid[dim][lower_idx + 1] - grid[dim][lower_idx]);
lower_idxs.push(lower_idx);
interp_diffs.push(interp_diff);
}
// `interp_vals` contains all values surrounding the point of interest, starting with shape (2, 2, ...) in N dimensions
// this gets mutated and reduces in dimension each iteration, filling with the next values to interpolate with
// this ends up as a 0-dimensional array containing only the final interpolated value
let mut interp_vals = values_view
.slice_each_axis(|ax| {
let lower = lower_idxs[ax.axis.0];
ndarray::Slice::from(lower..=lower + 1)
})
.to_owned();
let mut index_permutations = get_index_permutations(interp_vals.shape());
// This loop interpolates in each dimension sequentially
// each outer loop iteration the dimensionality reduces by 1
// `interp_vals` ends up as a 0-dimensional array containing only the final interpolated value
for (dim, diff) in interp_diffs.into_iter().enumerate() {
let next_dim = n - 1 - dim;
let next_shape = vec![2; next_dim];
// Indeces used for saving results of this dimensions interpolation results
// assigned to `index_permutations` at end of loop to be used for indexing in next iteration
let next_idxs = get_index_permutations(&next_shape);
let mut intermediate_arr = Array::zeros(next_shape);
for i in 0..next_idxs.len() {
// `next_idxs` is always half the length of `index_permutations`
let l = index_permutations[i].as_slice();
let u = index_permutations[next_idxs.len() + i].as_slice();
// This calculation happens 2^(n-1) times in the first iteration of the outer loop,
// 2^(n-2) times in the second iteration, etc.
intermediate_arr[next_idxs[i].as_slice()] =
interp_vals[l] * (D::Elem::one() - diff) + interp_vals[u] * diff;
}
index_permutations = next_idxs;
interp_vals = intermediate_arr;
}
// return the only value contained within the 0-dimensional array
Ok(interp_vals.first().copied().unwrap())
}
/// Returns `true`.
fn allow_extrapolate(&self) -> bool {
true
}
}
impl<D> StrategyND<D> for Nearest
where
D: Data + RawDataClone + Clone,
D::Elem: Num + PartialOrd + Copy + Debug,
{
fn interpolate(
&self,
data: &InterpDataND<D>,
point: &[D::Elem],
) -> Result<D::Elem, InterpolateError> {
// Dimensionality
let mut n = data.values.ndim();
// Point can share up to N values of a grid point, which reduces the problem dimensionality
// i.e. the point shares one of three values of a 3-D grid point, then the interpolation becomes 2-D at that slice
// or if the point shares two of three values of a 3-D grid point, then the interpolation becomes 1-D
let mut point = point.to_owned();
let mut grid: Vec<_> = data.grid.iter().map(|arr| arr.view()).collect();
let mut values_view = data.values.view();
for dim in (0..n).rev() {
// Range is reversed so that removal doesn't affect indexing
if let Some(pos) = grid[dim]
.iter()
.position(|&grid_point| grid_point == point[dim])
{
point.remove(dim);
grid.remove(dim);
values_view.index_axis_inplace(Axis(dim), pos);
}
}
if values_view.len() == 1 {
// Supplied point is coincident with a grid point, so just return the value
return Ok(values_view.first().copied().unwrap());
}
// Simplified dimensionality
n = values_view.ndim();
// Extract the lower and upper indices for each dimension,
// as well as the fraction of how far the supplied point is between the surrounding grid points
let mut lower_idxs = Vec::with_capacity(n);
let mut lower_closers = Vec::with_capacity(n);
for dim in 0..n {
let lower_idx = find_nearest_index(grid[dim].view(), &point[dim]);
let lower_closer =
point[dim] - grid[dim][lower_idx] < grid[dim][lower_idx + 1] - point[dim];
lower_idxs.push(lower_idx);
lower_closers.push(lower_closer);
}
// `interp_vals` contains all values surrounding the point of interest, starting with shape (2, 2, ...) in N dimensions
// this gets mutated and reduces in dimension each iteration, filling with the next values to interpolate with
// this ends up as a 0-dimensional array containing only the final interpolated value
let mut interp_vals = values_view
.slice_each_axis(|ax| {
let lower = lower_idxs[ax.axis.0];
ndarray::Slice::from(lower..=lower + 1)
})
.to_owned();
let mut index_permutations = get_index_permutations(interp_vals.shape());
// This loop interpolates in each dimension sequentially
// each outer loop iteration the dimensionality reduces by 1
// `interp_vals` ends up as a 0-dimensional array containing only the final interpolated value
for (dim, lower_closer) in lower_closers.into_iter().enumerate() {
let next_dim = n - 1 - dim;
let next_shape = vec![2; next_dim];
// Indeces used for saving results of this dimensions interpolation results
// assigned to `index_permutations` at end of loop to be used for indexing in next iteration
let next_idxs = get_index_permutations(&next_shape);
let mut intermediate_arr = Array::zeros(next_shape);
for i in 0..next_idxs.len() {
// `next_idxs` is always half the length of `index_permutations`
let l = index_permutations[i].as_slice();
let u = index_permutations[next_idxs.len() + i].as_slice();
// This calculation happens 2^(n-1) times in the first iteration of the outer loop,
// 2^(n-2) times in the second iteration, etc.
intermediate_arr[next_idxs[i].as_slice()] = if lower_closer {
interp_vals[l]
} else {
interp_vals[u]
};
}
index_permutations = next_idxs;
interp_vals = intermediate_arr;
}
// return the only value contained within the 0-dimensional array
Ok(interp_vals.first().copied().unwrap())
}
/// Returns `false`.
fn allow_extrapolate(&self) -> bool {
false
}
}