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
extern crate num;
#[cfg(test)]
extern crate statrs;
pub fn spline_mov(
x_and_y:Vec<(f64, f64)>
)-> impl Fn(f64) -> f64
{
let x_and_y_diff:Vec<(f64, f64, f64)>=x_and_y.windows(2).map(|point_and_next|{
let (x_curr, y_curr)=point_and_next[0];
let (x_next, y_next)=point_and_next[1];
let x_diff=x_next-x_curr;
let y_diff=y_next-y_curr;
let ms=y_diff/x_diff;
(ms, x_diff, ms)
}).collect();
let mut c1s:Vec<(f64, f64, f64)>=vec![];
c1s.push(*x_and_y_diff.first().unwrap());
c1s.append(&mut x_and_y_diff.windows(2).map(|diff_point_and_next|{
let (_, x_curr, dy_dx_curr)=diff_point_and_next[0];
let (_, x_next, dy_dx_next)=diff_point_and_next[1];
let common=x_curr+x_next;
let c1s_val=if dy_dx_next*dy_dx_curr<=0.0 {
0.0
}
else {
3.0*common/((common+x_next)/dy_dx_curr+(common+x_curr)/dy_dx_next)
};
(c1s_val, x_curr, dy_dx_curr)
}).collect());
c1s.push(*x_and_y_diff.last().unwrap());
let c2_and_c3:Vec<(f64, f64)>=c1s.windows(2).map(|c1_diff|{
let (c1_curr, dx_curr, dy_dx_curr)=c1_diff[0];
let (c1_next, _, _)=c1_diff[1];
let inv_dx=1.0/dx_curr;
let common=c1_curr+c1_next-2.0*dy_dx_curr;
((dy_dx_curr-c1_curr-common)*inv_dx, common*inv_dx.powi(2))
}).collect();
move |x|{
let (found_index, results)=x_and_y.windows(2).enumerate().find(|(_, w)| {
let (x_curr, _)=w[0];
let (x_next, _)=w[1];
x>=x_curr && x<=x_next
}).unwrap();
let (x_next, y_next)=results[1];
let diff=x-x_next;
let diff_sq=diff.powi(2);
let (c1s_elem, _, _)=c1s[found_index];
let (c2s_elem, c3s_elem)=c2_and_c3[found_index];
y_next+c1s_elem*diff+c2s_elem*diff_sq+c3s_elem*diff*diff_sq
}
}
#[cfg(test)]
mod tests {
use monotone_spline::*;
#[test]
fn test_returns_value_at_knot(){
let x_y:Vec<(f64, f64)>=vec![(1.0, 2.0), (2.0, 2.5), (3.0, 3.0)];
let spline=spline_mov(x_y);
let x_y_2:Vec<(f64, f64)>=vec![(1.0, 2.0), (2.0, 2.5), (3.0, 3.0)];
for (x, y) in x_y_2.iter(){
assert_abs_diff_eq!(
spline(*x),
y,
epsilon=0.000000001
);
}
}
fn test_returns_in_between_value(){
let x_y:Vec<(f64, f64)>=vec![(1.0, 2.0), (2.0, 2.5), (3.0, 3.0)];
let spline=spline_mov(x_y);
let x_y_2:Vec<(f64, f64)>=vec![(1.0, 2.0), (2.0, 2.5), (3.0, 3.0)];
let test_x:Vec<f64>=vec![1.01, 1.02, 1.03, 1.4, 1.8, 1.98, 1.99, 2.01, 2.02, 2.03, 2.4, 2.8, 2.98, 2.99];
for x in test_x.iter(){
let x_d_ref=*x;
let spline_y=spline(x_d_ref);
let y_bounds=x_y_2.windows(2).find(|w|{
let (x_curr, _)=w[0];
let (x_next, _)=w[1];
x_next>x_d_ref && x_curr<x_d_ref
}).unwrap();
let (_, y_curr)=y_bounds[0];
let (_, y_next)=y_bounds[1];
assert_eq!(
spline_y>y_curr&&spline_y<y_next,
true
);
}
}
}