@@ -16,6 +16,8 @@ use crate::{
16
16
17
17
use super :: { Strategy , StrategyBuilder } ;
18
18
19
+ const AX0 : Axis = Axis ( 0 ) ;
20
+
19
21
#[ derive( Debug ) ]
20
22
pub struct CubicSpline ;
21
23
impl < Sd , Sx , D > StrategyBuilder < Sd , Sx , D > for CubicSpline
@@ -130,14 +132,14 @@ impl CubicSpline {
130
132
// RHS vector
131
133
let mut rhs: Array < Sd :: Elem , D > = Array :: zeros ( dim. clone ( ) ) ;
132
134
133
- let mut inner_rhs = rhs. slice_axis_mut ( Axis ( 0 ) , Slice :: from ( 1 ..-1 ) ) ;
134
- Zip :: from ( inner_rhs. axis_iter_mut ( Axis ( 0 ) ) )
135
+ let mut inner_rhs = rhs. slice_axis_mut ( AX0 , Slice :: from ( 1 ..-1 ) ) ;
136
+ Zip :: from ( inner_rhs. axis_iter_mut ( AX0 ) )
135
137
. and ( x. windows ( 3 ) )
136
- . and ( data. axis_windows ( Axis ( 0 ) , 3 ) )
138
+ . and ( data. axis_windows ( AX0 , 3 ) )
137
139
. for_each ( |rhs, x, data| {
138
- let y_left = data. index_axis ( Axis ( 0 ) , 0 ) ;
139
- let y_mid = data. index_axis ( Axis ( 0 ) , 1 ) ;
140
- let y_right = data. index_axis ( Axis ( 0 ) , 2 ) ;
140
+ let y_left = data. index_axis ( AX0 , 0 ) ;
141
+ let y_mid = data. index_axis ( AX0 , 1 ) ;
142
+ let y_right = data. index_axis ( AX0 , 2 ) ;
141
143
let x_left = x[ 0 ] ;
142
144
let x_mid = x[ 1 ] ;
143
145
let x_right = x[ 2 ] ;
@@ -152,19 +154,19 @@ impl CubicSpline {
152
154
) ;
153
155
} ) ;
154
156
155
- let rhs_0 = rhs. index_axis_mut ( Axis ( 0 ) , 0 ) ;
156
- let data_0 = data. index_axis ( Axis ( 0 ) , 0 ) ;
157
- let data_1 = data. index_axis ( Axis ( 0 ) , 1 ) ;
157
+ let rhs_0 = rhs. index_axis_mut ( AX0 , 0 ) ;
158
+ let data_0 = data. index_axis ( AX0 , 0 ) ;
159
+ let data_1 = data. index_axis ( AX0 , 1 ) ;
158
160
Zip :: from ( rhs_0)
159
161
. and ( data_0)
160
162
. and ( data_1)
161
163
. for_each ( |rhs_0, & y_0, & y_1| {
162
164
* rhs_0 = three * ( y_1 - y_0) / ( x_1 - x_0) . pow ( two) ;
163
165
} ) ;
164
166
165
- let rhs_n = rhs. index_axis_mut ( Axis ( 0 ) , len - 1 ) ;
166
- let data_n = data. index_axis ( Axis ( 0 ) , len - 1 ) ;
167
- let data_n1 = data. index_axis ( Axis ( 0 ) , len - 2 ) ;
167
+ let rhs_n = rhs. index_axis_mut ( AX0 , len - 1 ) ;
168
+ let data_n = data. index_axis ( AX0 , len - 1 ) ;
169
+ let data_n1 = data. index_axis ( AX0 , len - 2 ) ;
168
170
Zip :: from ( rhs_n)
169
171
. and ( data_n)
170
172
. and ( data_n1)
@@ -174,12 +176,12 @@ impl CubicSpline {
174
176
175
177
// now solving With Thomas algorithm
176
178
177
- let mut rhs_left = rhs. index_axis ( Axis ( 0 ) , 0 ) . into_owned ( ) ;
179
+ let mut rhs_left = rhs. index_axis ( AX0 , 0 ) . into_owned ( ) ;
178
180
for i in 1 ..len {
179
181
let w = a_low[ i] / a_mid[ i - 1 ] ;
180
182
a_mid[ i] -= w * a_up[ i - 1 ] ;
181
183
182
- let rhs = rhs. index_axis_mut ( Axis ( 0 ) , i) ;
184
+ let rhs = rhs. index_axis_mut ( AX0 , i) ;
183
185
Zip :: from ( rhs)
184
186
. and ( rhs_left. view_mut ( ) )
185
187
. for_each ( |rhs, rhs_left| {
@@ -190,17 +192,17 @@ impl CubicSpline {
190
192
}
191
193
192
194
let mut k = Array :: zeros ( dim) ;
193
- Zip :: from ( k. index_axis_mut ( Axis ( 0 ) , len - 1 ) )
194
- . and ( rhs. index_axis ( Axis ( 0 ) , len - 1 ) )
195
+ Zip :: from ( k. index_axis_mut ( AX0 , len - 1 ) )
196
+ . and ( rhs. index_axis ( AX0 , len - 1 ) )
195
197
. for_each ( |k, & rhs| {
196
198
* k = rhs / a_mid[ len - 1 ] ;
197
199
} ) ;
198
200
199
- let mut k_right = k. index_axis ( Axis ( 0 ) , len - 1 ) . into_owned ( ) ;
201
+ let mut k_right = k. index_axis ( AX0 , len - 1 ) . into_owned ( ) ;
200
202
for i in ( 0 ..len - 1 ) . rev ( ) {
201
- Zip :: from ( k. index_axis_mut ( Axis ( 0 ) , i) )
203
+ Zip :: from ( k. index_axis_mut ( AX0 , i) )
202
204
. and ( k_right. view_mut ( ) )
203
- . and ( rhs. index_axis ( Axis ( 0 ) , i) )
205
+ . and ( rhs. index_axis ( AX0 , i) )
204
206
. for_each ( |k, k_right, & rhs| {
205
207
let new_k = ( rhs - a_up[ i] * * k_right) / a_mid[ i] ;
206
208
* k = new_k;
@@ -211,12 +213,12 @@ impl CubicSpline {
211
213
let mut c_a = Array :: zeros ( a_b_dim. clone ( ) ) ;
212
214
let mut c_b = Array :: zeros ( a_b_dim) ;
213
215
for index in 0 ..len - 1 {
214
- Zip :: from ( c_a. index_axis_mut ( Axis ( 0 ) , index) )
215
- . and ( c_b. index_axis_mut ( Axis ( 0 ) , index) )
216
- . and ( k. index_axis ( Axis ( 0 ) , index) )
217
- . and ( k. index_axis ( Axis ( 0 ) , index + 1 ) )
218
- . and ( data. index_axis ( Axis ( 0 ) , index) )
219
- . and ( data. index_axis ( Axis ( 0 ) , index + 1 ) )
216
+ Zip :: from ( c_a. index_axis_mut ( AX0 , index) )
217
+ . and ( c_b. index_axis_mut ( AX0 , index) )
218
+ . and ( k. index_axis ( AX0 , index) )
219
+ . and ( k. index_axis ( AX0 , index + 1 ) )
220
+ . and ( data. index_axis ( AX0 , index) )
221
+ . and ( data. index_axis ( AX0 , index + 1 ) )
220
222
. for_each ( |c_a, c_b, & k, & k_right, & y, & y_right| {
221
223
* c_a = k * ( x[ index + 1 ] - x[ index] ) - ( y_right - y) ;
222
224
* c_b = ( y_right - y) - k_right * ( x[ index + 1 ] - x[ index] ) ;
@@ -260,8 +262,8 @@ where
260
262
let idx = interp. get_left_index ( x) ;
261
263
let ( x_left, data_left) = interp. get_point ( idx) ;
262
264
let ( x_right, data_right) = interp. get_point ( idx + 1 ) ;
263
- let a_left = self . a . index_axis ( Axis ( 0 ) , idx) ;
264
- let b_left = self . b . index_axis ( Axis ( 0 ) , idx) ;
265
+ let a_left = self . a . index_axis ( AX0 , idx) ;
266
+ let b_left = self . b . index_axis ( AX0 , idx) ;
265
267
let one: Sd :: Elem = cast ( 1.0 ) . unwrap_or_else ( || unimplemented ! ( ) ) ;
266
268
267
269
let t = ( x - x_left) / ( x_right - x_left) ;
0 commit comments