Skip to content

Commit 62a62cc

Browse files
kylecarowkylecarow
authored
Merge pull request #3 from NREL/ndarray
NDArray & generics rewrite
2 parents 16e3bf8 + f3930b5 commit 62a62cc

18 files changed

+1147
-903
lines changed

‎.github/workflows/test.yaml

-1
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@ name: test
22

33
on:
44
push:
5-
branches: [main]
65
pull_request:
76

87
env:

‎Cargo.toml

+10-2
Original file line numberDiff line numberDiff line change
@@ -5,18 +5,26 @@ edition = "2021"
55
description = "Numerical interpolation for N-dimensional rectilinear grids"
66
repository = "https://github.com/NREL/ninterp"
77
license = "BSD-3-Clause"
8-
keywords = ["interpolation", "multidimensional", "multilinear", "numerical", "approximation"]
8+
keywords = [
9+
"interpolation",
10+
"multidimensional",
11+
"multilinear",
12+
"numerical",
13+
"approximation",
14+
]
915
categories = ["mathematics"]
1016

1117
[dependencies]
1218
itertools = "0.13.0"
1319
ndarray = "0.16.1"
20+
num = "0.4.3"
1421
serde = { version = "1.0.210", optional = true, features = ["derive"] }
1522
thiserror = "1.0.64"
1623

1724
[dev-dependencies]
1825
criterion = "0.5.1"
19-
rand = "0.9.0"
26+
ndarray-rand = "0.15.0"
27+
approx = "0.5.1"
2028

2129
[[bench]]
2230
name = "benchmark"

‎benches/benchmark.rs

+57-73
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,16 @@
11
//! Benchmarks for 0/1/2/3/N-dimensional linear interpolation
22
//! Run these with `cargo bench`
33
4-
use criterion::{criterion_group, criterion_main, Criterion};
4+
use criterion::{black_box, criterion_group, criterion_main, Criterion};
55

66
use ndarray::prelude::*;
77
use ninterp::prelude::*;
8-
use rand::{self, rngs::StdRng, Rng, SeedableRng};
8+
9+
use ndarray_rand::rand::{prelude::StdRng, Rng, SeedableRng};
10+
use ndarray_rand::rand_distr::Uniform;
11+
use ndarray_rand::RandomExt;
12+
13+
const RANDOM_SEED: u64 = 1234567890;
914

1015
#[allow(non_snake_case)]
1116
/// 0-D interpolation (hardcoded)
@@ -18,130 +23,109 @@ fn benchmark_0D() {
1823
/// 0-D interpolation (multilinear interpolator)
1924
fn benchmark_0D_multi() {
2025
let interp_0d_multi = InterpND::new(
21-
vec![vec![]],
26+
vec![array![]],
2227
array![0.5].into_dyn(),
2328
Linear,
2429
Extrapolate::Error,
2530
)
2631
.unwrap();
27-
interp_0d_multi.interpolate(&[]).unwrap();
32+
interp_0d_multi.interpolate(black_box(&[])).unwrap();
2833
}
2934

3035
#[allow(non_snake_case)]
3136
/// 1-D interpolation (hardcoded)
3237
fn benchmark_1D() {
33-
let seed = 1234567890;
34-
let mut rng = StdRng::seed_from_u64(seed);
35-
let grid_data: Vec<f64> = (0..100).map(|x| x as f64).collect();
38+
let mut rng = StdRng::seed_from_u64(RANDOM_SEED);
39+
let grid_data: Array1<f64> = (0..100).map(|x| x as f64).collect();
3640
// Generate interpolator data (same as N-D benchmark)
37-
let values_data: Vec<f64> = (0..100).map(|_| rng.random::<f64>()).collect();
41+
let values_data = Array1::random_using(100, Uniform::new(0., 1.), &mut rng);
3842
// Create a 1-D interpolator with 100 data points
3943
let interp_1d = Interp1D::new(grid_data, values_data, Linear, Extrapolate::Error).unwrap();
4044
// Sample 1,000 points
41-
let points: Vec<f64> = (0..1_000).map(|_| rng.random::<f64>() * 99.).collect();
45+
let points: Vec<f64> = (0..1_000).map(|_| rng.gen::<f64>() * 99.).collect();
4246
for point in points {
43-
interp_1d.interpolate(&[point]).unwrap();
47+
interp_1d.interpolate(black_box(&[point])).unwrap();
4448
}
4549
}
4650

4751
#[allow(non_snake_case)]
4852
/// 1-D interpolation (multilinear interpolator)
4953
fn benchmark_1D_multi() {
50-
let seed = 1234567890;
51-
let mut rng = StdRng::seed_from_u64(seed);
54+
let mut rng = StdRng::seed_from_u64(RANDOM_SEED);
5255
// Generate interpolator data (same as hardcoded benchmark)
53-
let grid_data: Vec<f64> = (0..100).map(|x| x as f64).collect();
54-
let values_data: Vec<f64> = (0..100).map(|_| rng.random::<f64>()).collect();
56+
let grid_data: Array1<f64> = (0..100).map(|x| x as f64).collect();
57+
let values_data = Array1::random_using(100, Uniform::new(0., 1.), &mut rng).into_dyn();
5558
// Create an N-D interpolator with 100x100 data (10,000 points)
56-
let interp_1d_multi = InterpND::new(
57-
vec![grid_data],
58-
ArrayD::from_shape_vec(IxDyn(&[100]), values_data).unwrap(),
59-
Linear,
60-
Extrapolate::Error,
61-
)
62-
.unwrap();
59+
let interp_1d_multi =
60+
InterpND::new(vec![grid_data], values_data, Linear, Extrapolate::Error).unwrap();
6361
// Sample 1,000 points
64-
let points: Vec<f64> = (0..1_000).map(|_| rng.random::<f64>() * 99.).collect();
62+
let points: Vec<f64> = (0..1_000).map(|_| rng.gen::<f64>() * 99.).collect();
6563
for point in points {
66-
interp_1d_multi.interpolate(&[point]).unwrap();
64+
interp_1d_multi.interpolate(black_box(&[point])).unwrap();
6765
}
6866
}
6967

7068
#[allow(non_snake_case)]
7169
/// 2-D interpolation (hardcoded)
7270
fn benchmark_2D() {
73-
let seed = 1234567890;
74-
let mut rng = StdRng::seed_from_u64(seed);
75-
let grid_data: Vec<f64> = (0..100).map(|x| x as f64).collect();
76-
// Generate interpolator data (same as N-D benchmark) and arrange into `Vec<Vec<f64>>`
77-
let values_data: Vec<f64> = (0..10_000).map(|_| rng.random::<f64>()).collect();
78-
let values_data: Vec<Vec<f64>> = (0..100)
79-
.map(|x| values_data[(100 * x)..(100 + 100 * x)].into())
80-
.collect();
71+
let mut rng = StdRng::seed_from_u64(RANDOM_SEED);
72+
let grid_data: Array1<f64> = (0..100).map(|x| x as f64).collect();
73+
let values_data = Array2::random_using((100, 100), Uniform::new(0., 1.), &mut rng);
8174
// Create a 2-D interpolator with 100x100 data (10,000 points)
8275
let interp_2d = Interp2D::new(
83-
grid_data.clone(),
84-
grid_data.clone(),
85-
values_data,
76+
grid_data.view(),
77+
grid_data.view(),
78+
values_data.view(),
8679
Linear,
8780
Extrapolate::Error,
8881
)
8982
.unwrap();
9083
// Sample 1,000 points
9184
let points: Vec<Vec<f64>> = (0..1_000)
92-
.map(|_| vec![rng.random::<f64>() * 99., rng.random::<f64>() * 99.])
85+
.map(|_| vec![rng.gen::<f64>() * 99., rng.gen::<f64>() * 99.])
9386
.collect();
9487
for point in points {
95-
interp_2d.interpolate(&point).unwrap();
88+
interp_2d.interpolate(black_box(&point)).unwrap();
9689
}
9790
}
9891

9992
#[allow(non_snake_case)]
10093
/// 2-D interpolation (multilinear interpolator)
10194
fn benchmark_2D_multi() {
102-
let seed = 1234567890;
103-
let mut rng = StdRng::seed_from_u64(seed);
95+
let mut rng = StdRng::seed_from_u64(RANDOM_SEED);
10496
// Generate interpolator data (same as hardcoded benchmark)
105-
let grid_data: Vec<f64> = (0..100).map(|x| x as f64).collect();
106-
let values_data: Vec<f64> = (0..10_000).map(|_| rng.random::<f64>()).collect();
97+
let grid_data: Array1<f64> = (0..100).map(|x| x as f64).collect();
98+
let values_data = Array2::random_using((100, 100), Uniform::new(0., 1.), &mut rng).into_dyn();
10799
// Create an N-D interpolator with 100x100 data (10,000 points)
108100
let interp_2d_multi = InterpND::new(
109-
vec![grid_data.clone(), grid_data.clone()],
110-
ArrayD::from_shape_vec(IxDyn(&[100, 100]), values_data).unwrap(),
101+
vec![grid_data.view(), grid_data.view()],
102+
values_data.view(),
111103
Linear,
112104
Extrapolate::Error,
113105
)
114106
.unwrap();
115107
// Sample 1,000 points
116108
let points: Vec<Vec<f64>> = (0..1_000)
117-
.map(|_| vec![rng.random::<f64>() * 99., rng.random::<f64>() * 99.])
109+
.map(|_| vec![rng.gen::<f64>() * 99., rng.gen::<f64>() * 99.])
118110
.collect();
119111
for point in points {
120-
interp_2d_multi.interpolate(&point).unwrap();
112+
interp_2d_multi.interpolate(black_box(&point)).unwrap();
121113
}
122114
}
123115

124116
#[allow(non_snake_case)]
125117
/// 3-D interpolation (hardcoded)
126118
fn benchmark_3D() {
127-
let seed = 1234567890;
128-
let mut rng = StdRng::seed_from_u64(seed);
129-
let grid_data: Vec<f64> = (0..100).map(|x| x as f64).collect();
119+
let mut rng = StdRng::seed_from_u64(RANDOM_SEED);
120+
let grid_data: Array1<f64> = (0..100).map(|x| x as f64).collect();
130121
// Generate interpolator data (same as N-D benchmark) and arrange into `Vec<Vec<Vec<f64>>>`
131-
let values_data: Vec<f64> = (0..1_000_000).map(|_| rng.random::<f64>()).collect();
132-
let values_data: Vec<Vec<Vec<f64>>> = (0..100)
133-
.map(|x| {
134-
(0..100)
135-
.map(|y| values_data[(100 * (y + 100 * x))..(100 + 100 * (y + 100 * x))].into())
136-
.collect()
137-
})
138-
.collect();
122+
let values_data = Array3::random_using((100, 100, 100), Uniform::new(0., 1.), &mut rng);
139123
// Create a 3-D interpolator with 100x100x100 data (1,000,000 points)
140124
let interp_3d = Interp3D::new(
141-
grid_data.clone(),
142-
grid_data.clone(),
143-
grid_data.clone(),
144-
values_data,
125+
grid_data.view(),
126+
grid_data.view(),
127+
grid_data.view(),
128+
values_data.view(),
145129
Linear,
146130
Extrapolate::Error,
147131
)
@@ -150,29 +134,29 @@ fn benchmark_3D() {
150134
let points: Vec<Vec<f64>> = (0..1_000)
151135
.map(|_| {
152136
vec![
153-
rng.random::<f64>() * 99.,
154-
rng.random::<f64>() * 99.,
155-
rng.random::<f64>() * 99.,
137+
rng.gen::<f64>() * 99.,
138+
rng.gen::<f64>() * 99.,
139+
rng.gen::<f64>() * 99.,
156140
]
157141
})
158142
.collect();
159143
for point in points {
160-
interp_3d.interpolate(&point).unwrap();
144+
interp_3d.interpolate(black_box(&point)).unwrap();
161145
}
162146
}
163147

164148
#[allow(non_snake_case)]
165149
/// 3-D interpolation (multilinear interpolator)
166150
fn benchmark_3D_multi() {
167-
let seed = 1234567890;
168-
let mut rng = StdRng::seed_from_u64(seed);
151+
let mut rng = StdRng::seed_from_u64(RANDOM_SEED);
169152
// Generate interpolator data (same as hardcoded benchmark)
170-
let grid_data: Vec<f64> = (0..100).map(|x| x as f64).collect();
171-
let values_data: Vec<f64> = (0..1_000_000).map(|_| rng.random::<f64>()).collect();
153+
let grid_data: Array1<f64> = (0..100).map(|x| x as f64).collect();
154+
let values_data =
155+
Array3::random_using((100, 100, 100), Uniform::new(0., 1.), &mut rng).into_dyn();
172156
// Create an N-D interpolator with 100x100x100 data (1,000,000 points)
173157
let interp_3d_multi = InterpND::new(
174-
vec![grid_data.clone(), grid_data.clone(), grid_data.clone()],
175-
ArrayD::from_shape_vec(IxDyn(&[100, 100, 100]), values_data).unwrap(),
158+
vec![grid_data.view(), grid_data.view(), grid_data.view()],
159+
values_data.view(),
176160
Linear,
177161
Extrapolate::Error,
178162
)
@@ -181,14 +165,14 @@ fn benchmark_3D_multi() {
181165
let points: Vec<Vec<f64>> = (0..1_000)
182166
.map(|_| {
183167
vec![
184-
rng.random::<f64>() * 99.,
185-
rng.random::<f64>() * 99.,
186-
rng.random::<f64>() * 99.,
168+
rng.gen::<f64>() * 99.,
169+
rng.gen::<f64>() * 99.,
170+
rng.gen::<f64>() * 99.,
187171
]
188172
})
189173
.collect();
190174
for point in points {
191-
interp_3d_multi.interpolate(&point).unwrap();
175+
interp_3d_multi.interpolate(black_box(&point)).unwrap();
192176
}
193177
}
194178

‎examples/custom_strategy.rs

+19-8
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,29 @@
1+
use ndarray::prelude::*;
2+
13
use ninterp::prelude::*;
24
use ninterp::strategy::*;
35

46
// Debug must be derived for custom strategies
57
#[derive(Debug)]
68
struct CustomStrategy;
79

8-
// Implement strategy for 2-D interpolation
9-
impl Strategy2D for CustomStrategy {
10+
// Implement strategy for 2-D f32 interpolation
11+
impl<D> Strategy2D<D> for CustomStrategy
12+
where
13+
// Implement for any 2-D interpolator where the contained type is `f32`
14+
// e.g. `Array2<f32>`, `ArrayView2<f32>`, `CowArray<<'a, f32>, Ix2>`, etc.
15+
// For a more generic bound, consider introducing a bound for D::Elem
16+
// e.g. D::Elem: num_traits::Num + PartialOrd
17+
D: ndarray::Data<Elem = f32>,
18+
{
1019
fn interpolate(
1120
&self,
12-
_data: &Data2D,
13-
point: &[f64; 2],
14-
) -> Result<f64, ninterp::error::InterpolateError> {
21+
_data: &InterpData2D<D>,
22+
point: &[f32; 2],
23+
) -> Result<f32, ninterp::error::InterpolateError> {
1524
// Dummy interpolation strategy, product of all point components
25+
// Here we could access the `InterpData2D` instead,
26+
// but this is just an example.
1627
Ok(point.iter().fold(1., |acc, x| acc * x))
1728
}
1829

@@ -29,9 +40,9 @@ impl Strategy2D for CustomStrategy {
2940

3041
fn main() {
3142
let interp = Interp2D::new(
32-
vec![0., 2., 4.],
33-
vec![0., 4., 8.],
34-
vec![vec![0., 0., 0.], vec![0., 0., 0.], vec![0., 0., 0.]],
43+
array![0., 2., 4.],
44+
array![0., 4., 8.],
45+
array![[0., 0., 0.], [0., 0., 0.], [0., 0., 0.]],
3546
CustomStrategy,
3647
Extrapolate::Error,
3748
)

‎examples/dynamic_interpolator.rs

+8-6
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,14 @@
1+
use ndarray::prelude::*;
2+
13
use ninterp::prelude::*;
24

35
fn main() {
46
// Create `Interpolator` trait object
5-
let mut boxed: Box<dyn Interpolator> = Box::new(
7+
let mut boxed: Box<dyn Interpolator<_>> = Box::new(
68
Interp2D::new(
7-
vec![0., 1.],
8-
vec![0., 1.],
9-
vec![vec![2., 4.], vec![4., 16.]],
9+
array![0., 1.],
10+
array![0., 1.],
11+
array![[2., 4.], [4., 16.]],
1012
Linear,
1113
Extrapolate::Enable,
1214
)
@@ -16,8 +18,8 @@ fn main() {
1618
// Change underlying interpolator
1719
boxed = Box::new(
1820
Interp1D::new(
19-
vec![0., 1., 2.],
20-
vec![0., 4., 8.],
21+
array![0., 1., 2.],
22+
array![0., 4., 8.],
2123
Nearest,
2224
Extrapolate::Error,
2325
)

‎examples/dynamic_strategy.rs

+5-3
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,15 @@
1+
use ndarray::prelude::*;
2+
13
use ninterp::prelude::*;
24
use ninterp::strategy::Strategy1D;
35

46
fn main() {
57
// Create mutable interpolator
68
let mut interp = Interp1D::new(
7-
vec![0., 1., 2.],
8-
vec![0., 3., 6.],
9+
array![0., 1., 2.],
10+
array![0., 3., 6.],
911
// Provide the strategy as a trait object
10-
Box::new(Linear) as Box<dyn Strategy1D>,
12+
Box::new(Linear) as Box<dyn Strategy1D<_>>,
1113
Extrapolate::Error,
1214
)
1315
.unwrap();

0 commit comments

Comments
 (0)