rsl_interpolation/types/
linear.rs

1use std::marker::PhantomData;
2
3use crate::Accelerator;
4use crate::InterpType;
5use crate::Interpolation;
6use crate::{DomainError, InterpolationError};
7
8use crate::{check_if_inbounds, check1d_data};
9
10const MIN_SIZE: usize = 2;
11
12/// Linear Interpolation type.
13///
14/// The simplest type of interpolation.
15#[doc(alias = "gsl_interp_linear")]
16pub struct Linear;
17
18impl<T> InterpType<T> for Linear
19where
20    T: crate::Num,
21{
22    type Interpolation = LinearInterp<T>;
23
24    /// Constructs a Linear Interpolator.
25    ///
26    /// # Example
27    ///
28    /// ```
29    /// # use rsl_interpolation::*;
30    /// #
31    /// # fn main() -> Result<(), InterpolationError>{
32    /// let xa = [0.0, 1.0, 2.0];
33    /// let ya = [0.0, 2.0, 4.0];
34    /// let interp = Linear.build(&xa, &ya)?;
35    /// # Ok(())
36    /// # }
37    /// ```
38    fn build(&self, xa: &[T], ya: &[T]) -> Result<LinearInterp<T>, InterpolationError> {
39        check1d_data(xa, ya, MIN_SIZE)?;
40        Ok(LinearInterp {
41            _variable_type: PhantomData,
42        })
43    }
44
45    fn name(&self) -> &str {
46        "Linear"
47    }
48
49    fn min_size(&self) -> usize {
50        MIN_SIZE
51    }
52}
53
54// ===============================================================================================
55
56/// Linear Interpolator.
57///
58/// Provides all the evaluation methods.
59///
60/// Should be constructed through the [`Linear`] type.
61pub struct LinearInterp<T> {
62    _variable_type: PhantomData<T>,
63}
64
65impl<T> Interpolation<T> for LinearInterp<T>
66where
67    T: crate::Num,
68{
69    fn eval(&self, xa: &[T], ya: &[T], x: T, acc: &mut Accelerator) -> Result<T, DomainError> {
70        check_if_inbounds(xa, x)?;
71        let index = acc.find(xa, x);
72
73        let xlo = xa[index];
74        let xhi = xa[index + 1];
75        let ylo = ya[index];
76        let yhi = ya[index + 1];
77
78        let dx = xhi - xlo;
79
80        debug_assert!(dx > T::zero());
81        Ok(ylo + (x - xlo) / dx * (yhi - ylo))
82    }
83
84    fn eval_deriv(
85        &self,
86        xa: &[T],
87        ya: &[T],
88        x: T,
89        acc: &mut Accelerator,
90    ) -> Result<T, DomainError> {
91        check_if_inbounds(xa, x)?;
92        let index = acc.find(xa, x);
93
94        let xlo = xa[index];
95        let xhi = xa[index + 1];
96        let ylo = ya[index];
97        let yhi = ya[index + 1];
98
99        let dx = xhi - xlo;
100        let dy = yhi - ylo;
101
102        debug_assert!(dx > T::zero());
103        Ok(dy / dx)
104    }
105
106    /// Always returns `0`.
107    #[allow(unused_variables)]
108    fn eval_deriv2(
109        &self,
110        xa: &[T],
111        ya: &[T],
112        x: T,
113        acc: &mut Accelerator,
114    ) -> Result<T, DomainError> {
115        check_if_inbounds(xa, x)?;
116        Ok(T::zero())
117    }
118
119    fn eval_integ(
120        &self,
121        xa: &[T],
122        ya: &[T],
123        a: T,
124        b: T,
125        acc: &mut Accelerator,
126    ) -> Result<T, DomainError> {
127        check_if_inbounds(xa, a)?;
128        check_if_inbounds(xa, b)?;
129        let index_a = acc.find(xa, a);
130        let index_b = acc.find(xa, b);
131
132        let mut result = T::zero();
133
134        for i in index_a..=index_b {
135            let xlo = xa[i];
136            let xhi = xa[i + 1];
137            let ylo = ya[i];
138            let yhi = ya[i + 1];
139
140            let dx = xhi - xlo;
141            let d = (yhi - ylo) / dx;
142            let half = T::from(0.5).unwrap();
143
144            if dx.is_zero() {
145                continue;
146            }
147
148            if (i == index_a) | (i == index_b) {
149                let x1 = if i == index_a { a } else { xlo };
150                let x2 = if i == index_b { b } else { xhi };
151                result += (x2 - x1) * (ylo + half * d * ((x2 - xlo) + (x1 - xlo)));
152            } else {
153                result += half * dx * (ylo + yhi);
154            }
155        }
156        Ok(result)
157    }
158}
rsl_interpolation/types/linear.rs

rsl_interpolation/types/
linear.rs
