rsl_interpolation/types/
akima.rs

1use std::collections::VecDeque;
2
3use crate::types::utils::integ_eval;
4use crate::types::utils::{check1d_data, check_if_inbounds};
5use crate::Accelerator;
6use crate::DomainError;
7use crate::InterpType;
8use crate::Interpolation;
9use crate::InterpolationError;
10
11const MIN_SIZE: usize = 5;
12
13/// Akima Interpolation type.
14///
15/// Non-rounded Akima spline with natural boundary conditions. This method uses the non-rounded
16/// corner algorithm of Wodicka.
17#[doc(alias = "gsl_interp_akima")]
18pub struct Akima;
19
20impl<T> InterpType<T> for Akima
21where
22    T: crate::Num,
23{
24    type Interpolation = AkimaInterp<T>;
25
26    /// Constructs an Akima Interpolator.
27    ///
28    /// ## Example
29    ///
30    /// ```
31    /// # use rsl_interpolation::*;
32    /// #
33    /// # fn main() -> Result<(), InterpolationError>{
34    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
35    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
36    /// let interp = Akima.build(&xa, &ya)?;
37    /// # Ok(())
38    /// # }
39    /// ```
40    fn build(&self, xa: &[T], ya: &[T]) -> Result<AkimaInterp<T>, InterpolationError> {
41        check1d_data(xa, ya, MIN_SIZE)?;
42
43        let size = xa.len();
44        let two = T::from(2.0).unwrap();
45        let three = T::from(3.0).unwrap();
46
47        // All m indices are shifted by +2
48        let mut m = VecDeque::<T>::with_capacity(size);
49        for i in 0..=size - 2 {
50            m.push_back((ya[i + 1] - ya[i]) / (xa[i + 1] - xa[i]));
51        }
52
53        // Non-periodic boundary conditions
54        m.push_front(two * m[0] - m[1]);
55        m.push_front(three * m[1] - two * m[2]);
56        m.push_back(two * m[size] - m[size - 1]);
57        m.push_back(three * m[size] - two * m[size - 1]);
58        let m = m.make_contiguous().to_vec();
59
60        let (b, c, d) = akima_calc(xa, &m);
61
62        let state = AkimaInterp { b, c, d, m };
63        Ok(state)
64    }
65
66    fn name(&self) -> &str {
67        "Akima"
68    }
69
70    fn min_size(&self) -> usize {
71        MIN_SIZE
72    }
73}
74
75// ===============================================================================================
76
77/// Akima Interpolator.
78///
79/// Provides all the evaluation methods.
80///
81/// Should be constructed through the [`Akima`] type.
82#[allow(dead_code)]
83#[doc(alias = "gsl_akima_interp")]
84pub struct AkimaInterp<T>
85where
86    T: crate::Num,
87{
88    b: Vec<T>,
89    c: Vec<T>,
90    d: Vec<T>,
91    m: Vec<T>,
92}
93
94impl<T> Interpolation<T> for AkimaInterp<T>
95where
96    T: crate::Num,
97{
98    fn eval(&self, xa: &[T], ya: &[T], x: T, acc: &mut Accelerator) -> Result<T, DomainError> {
99        akima_eval(xa, ya, (&self.b, &self.c, &self.d), x, acc)
100    }
101
102    #[allow(unused_variables)]
103    fn eval_deriv(
104        &self,
105        xa: &[T],
106        ya: &[T],
107        x: T,
108        acc: &mut Accelerator,
109    ) -> Result<T, DomainError> {
110        akima_eval_deriv(xa, (&self.b, &self.c, &self.d), x, acc)
111    }
112
113    #[allow(unused_variables)]
114    fn eval_deriv2(
115        &self,
116        xa: &[T],
117        ya: &[T],
118        x: T,
119        acc: &mut Accelerator,
120    ) -> Result<T, DomainError> {
121        akima_eval_deriv2(xa, (&self.c, &self.d), x, acc)
122    }
123
124    fn eval_integ(
125        &self,
126        xa: &[T],
127        ya: &[T],
128        a: T,
129        b: T,
130        acc: &mut Accelerator,
131    ) -> Result<T, DomainError> {
132        akima_eval_integ(xa, ya, (&self.b, &self.c, &self.d), a, b, acc)
133    }
134}
135
136//=================================================================================================
137
138/// Akima Periodic Interpolation type.
139///
140/// Non-rounded Akima spline with natural boundary conditions. This method uses the non-rounded
141/// corner algorithm of Wodicka.
142pub struct AkimaPeriodic;
143
144impl<T> InterpType<T> for AkimaPeriodic
145where
146    T: crate::Num,
147{
148    type Interpolation = AkimaPeriodicInterp<T>;
149
150    /// Constructs an Akima Periodic Interpolator.
151    ///
152    /// ## Example
153    ///
154    /// ```
155    /// # use rsl_interpolation::*;
156    /// #
157    /// # fn main() -> Result<(), InterpolationError>{
158    /// let xa = [0.0, 1.0, 2.0, 3.0, 4.0];
159    /// let ya = [0.0, 2.0, 4.0, 6.0, 8.0];
160    /// let interp = AkimaPeriodic.build(&xa, &ya)?;
161    /// # Ok(())
162    /// # }
163    /// ```
164    fn build(&self, xa: &[T], ya: &[T]) -> Result<AkimaPeriodicInterp<T>, InterpolationError> {
165        check1d_data(xa, ya, MIN_SIZE)?;
166
167        let size = xa.len();
168
169        // All m indices are shifted by +2
170        let mut m = VecDeque::<T>::with_capacity(size + 3);
171        for i in 0..=size - 2 {
172            m.push_back((ya[i + 1] - ya[i]) / (xa[i + 1] - xa[i]));
173        }
174
175        // Periodic boundary conditions
176        m.push_front(m[size - 1 - 1]);
177        m.push_front(m[size - 1 - 1]);
178        m.push_back(m[2]);
179        m.push_back(m[3]);
180        let m = m.make_contiguous().to_vec();
181
182        let (b, c, d) = akima_calc(xa, &m);
183
184        let state = AkimaPeriodicInterp { b, c, d, m };
185        Ok(state)
186    }
187
188    fn name(&self) -> &str {
189        "Akima Periodic"
190    }
191
192    fn min_size(&self) -> usize {
193        MIN_SIZE
194    }
195}
196
197// ===============================================================================================
198
199/// Akima Interpolator.
200///
201/// Provides all the evaluation methods.
202///
203/// Should be constructed through the [`Akima`] type.
204#[allow(dead_code)]
205#[doc(alias = "gsl_interp_akima_periodic")]
206pub struct AkimaPeriodicInterp<T>
207where
208    T: crate::Num,
209{
210    c: Vec<T>,
211    b: Vec<T>,
212    d: Vec<T>,
213    m: Vec<T>,
214}
215
216impl<T> Interpolation<T> for AkimaPeriodicInterp<T>
217where
218    T: crate::Num,
219{
220    fn eval(&self, xa: &[T], ya: &[T], x: T, acc: &mut Accelerator) -> Result<T, DomainError> {
221        akima_eval(xa, ya, (&self.b, &self.c, &self.d), x, acc)
222    }
223
224    #[allow(unused_variables)]
225    fn eval_deriv(
226        &self,
227        xa: &[T],
228        ya: &[T],
229        x: T,
230        acc: &mut Accelerator,
231    ) -> Result<T, DomainError> {
232        akima_eval_deriv(xa, (&self.b, &self.c, &self.d), x, acc)
233    }
234
235    #[allow(unused_variables)]
236    fn eval_deriv2(
237        &self,
238        xa: &[T],
239        ya: &[T],
240        x: T,
241        acc: &mut Accelerator,
242    ) -> Result<T, DomainError> {
243        akima_eval_deriv2(xa, (&self.c, &self.d), x, acc)
244    }
245
246    fn eval_integ(
247        &self,
248        xa: &[T],
249        ya: &[T],
250        a: T,
251        b: T,
252        acc: &mut Accelerator,
253    ) -> Result<T, DomainError> {
254        akima_eval_integ(xa, ya, (&self.b, &self.c, &self.d), a, b, acc)
255    }
256}
257
258//=================================================================================================
259
260fn akima_eval<T>(
261    xa: &[T],
262    ya: &[T],
263    state: (&[T], &[T], &[T]),
264    x: T,
265    acc: &mut Accelerator,
266) -> Result<T, DomainError>
267where
268    T: crate::Num,
269{
270    check_if_inbounds(xa, x)?;
271    let index = acc.find(xa, x);
272
273    let xlo = xa[index];
274    let delx = x - xlo;
275    let b = state.0[index];
276    let c = state.1[index];
277    let d = state.2[index];
278
279    Ok(ya[index] + delx * (b + delx * (c + d * delx)))
280}
281
282fn akima_eval_deriv<T>(
283    xa: &[T],
284    state: (&[T], &[T], &[T]),
285    x: T,
286    acc: &mut Accelerator,
287) -> Result<T, DomainError>
288where
289    T: crate::Num,
290{
291    check_if_inbounds(xa, x)?;
292    let two = T::from(2).unwrap();
293    let three = T::from(3).unwrap();
294
295    let index = acc.find(xa, x);
296
297    let xlo = xa[index];
298    let delx = x - xlo;
299    let b = state.0[index];
300    let c = state.1[index];
301    let d = state.2[index];
302
303    Ok(b + delx * (two * c + three * d * delx))
304}
305
306#[inline(always)]
307fn akima_eval_deriv2<T>(
308    xa: &[T],
309    state: (&[T], &[T]),
310    x: T,
311    acc: &mut Accelerator,
312) -> Result<T, DomainError>
313where
314    T: crate::Num,
315{
316    check_if_inbounds(xa, x)?;
317    let two = T::from(2).unwrap();
318    let six = T::from(6).unwrap();
319
320    let index = acc.find(xa, x);
321
322    let xlo = xa[index];
323    let delx = x - xlo;
324    let c = state.0[index];
325    let d = state.1[index];
326
327    Ok(two * c + six * d * delx)
328}
329
330fn akima_eval_integ<T>(
331    xa: &[T],
332    ya: &[T],
333    state: (&[T], &[T], &[T]),
334    a: T,
335    b: T,
336    acc: &mut Accelerator,
337) -> Result<T, DomainError>
338where
339    T: crate::Num,
340{
341    check_if_inbounds(xa, a)?;
342    check_if_inbounds(xa, b)?;
343    let index_a = acc.find(xa, a);
344    let index_b = acc.find(xa, b);
345    let bs = state.0;
346    let cs = state.1;
347    let ds = state.2;
348
349    let quarter = T::from(0.25).unwrap();
350    let half = T::from(0.5).unwrap();
351    let third = T::from(1.0 / 3.0).unwrap();
352    let mut result = T::zero();
353
354    for i in index_a..=index_b {
355        let xlo = xa[i];
356        let xhi = xa[i + 1];
357        let ylo = ya[i];
358
359        let dx = xhi - xlo;
360
361        // If two x points are the same
362        if dx.is_zero() {
363            continue;
364        }
365
366        if (i == index_a) | (i == index_b) {
367            let x1 = if i == index_a { a } else { xlo };
368            let x2 = if i == index_b { b } else { xhi };
369            result += integ_eval(ylo, bs[i], cs[i], ds[i], xlo, x1, x2);
370        } else {
371            result += dx * (ylo + dx * (half * bs[i] + dx * (third * cs[i] + quarter * ds[i] * dx)))
372        }
373    }
374    Ok(result)
375}
376
377/// Common Calculation
378fn akima_calc<T>(xa: &[T], m: &[T]) -> (Vec<T>, Vec<T>, Vec<T>)
379where
380    T: crate::Num,
381{
382    let size = xa.len();
383    let two = T::from(2.0).unwrap();
384    let three = T::from(3.0).unwrap();
385    let mut b = Vec::<T>::with_capacity(size - 1);
386    let mut c = Vec::<T>::with_capacity(size - 1);
387    let mut d = Vec::<T>::with_capacity(size - 1);
388
389    for i in 0..size - 1 {
390        let ne = (m[i + 3] - m[i + 2]).abs() + (m[i + 1] - m[i]).abs();
391        if ne.is_zero() {
392            b.push(m[i + 2]);
393            c.push(T::zero());
394            d.push(T::zero());
395        } else {
396            let hi = xa[i + 1] - xa[i];
397            let nenext = (m[i + 4] - m[i + 3]).abs() + (m[i + 2] - m[i + 1]).abs();
398            let ai = (m[i + 1] - m[i]).abs() / ne;
399            let ai_plus1: T;
400            let tli_plus1: T;
401            if nenext.is_zero() {
402                tli_plus1 = m[i + 2];
403            } else {
404                ai_plus1 = (m[i + 2] - m[i + 1]).abs() / nenext;
405                tli_plus1 = (T::one() - ai_plus1) * m[i + 2] + ai_plus1 * m[i + 3];
406            }
407            b.push((T::one() - ai) * m[i + 1] + ai * m[i + 2]);
408            c.push((three * m[i + 2] - two * b[i] - tli_plus1) / hi);
409            d.push((b[i] + tli_plus1 - two * m[i + 2]) / hi.powi(2));
410        }
411    }
412    (b, c, d)
413}
rsl_interpolation/types/akima.rs

rsl_interpolation/types/
akima.rs
