hologram 0.1.4

Interpolation library with multipurpose Radial Basis Function (RBF).
Documentation 
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
109
110
111

use std::f64::consts::PI;

use crate::Interpolator;

/// Represents a Lagrangian interpolator, used to perform polynomial interpolation based on a set of
/// input and output points. NB: Only supports x: 1d and y: 3d.
///
/// # Fields
/// * `x`: A vector of input points where interpolation is performed.
/// * `y`: A vector of output points corresponding to the input points,
///        where each output point is a 3-element array [x, y, z].
/// * `denominators`: Precomputed denominators used in the Lagrangian interpolation formula to
///                   optimise the calculation.
pub struct Lagrangian {
    x: Vec<f64>,
    y: Vec<[f64; 3]>,
    denominators: Vec<f64>,
}

impl Lagrangian {
    pub fn new(x: Vec<f64>, y: Vec<[f64; 3]>, chebyshev: bool) -> Self {
        assert_eq!(x.len(), y.len(), "x and y must have the same length");

        // Chebyshev nodes
        let n = x.len();
        let x = if chebyshev {
            let a = *x.iter().min_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
            let b = *x.iter().max_by(|a, b| a.partial_cmp(b).unwrap()).unwrap();
            (0..n)
                .map(|i| {
                    let cos_value = ((2 * i + 1) as f64 * PI / (2.0 * n as f64)).cos();
                    0.5 * ((b - a) * cos_value + (b + a))
                })
                .collect::<Vec<f64>>()
        } else {
            x
        };

        // Compute denominators in advance
        let mut denominators = vec![1.0; n];

        for i in 0..n {
            for j in 0..n {
                if i != j {
                    denominators[i] *= x[i] - x[j]; //.max(std::f64::EPSILON)
                }
            }
            if denominators[i] == 0.0 {
                panic!("Division by zero in basis function calculation.");
            }
        }

        Lagrangian { x, y, denominators }
    }

    fn compute_basis(&self, x_n: f64, i: usize) -> f64 {
        let mut basis = 1.0;
        for j in 0..self.x.len() {
            if i != j {
                basis *= x_n - self.x[j];
            }
        }
        basis / self.denominators[i]
    }
}

impl Interpolator<f64, [f64; 3]> for Lagrangian {
    fn predict(&self, x_new: &[f64]) -> Result<Vec<[f64; 3]>, String> {
        let mut results = Vec::with_capacity(x_new.len());

        for &x_n in x_new {
            let mut result = [0.0; 3];

            for i in 0..self.x.len() {
                let basis = self.compute_basis(x_n, i);
                for k in 0..3 {
                    result[k] += self.y[i][k] * basis;
                }
            }

            if result.iter().any(|&val| val.is_nan()) {
                return Err(String::from("NaN values in output"));
            }

            results.push(result);
        }

        Ok(results.into_iter().rev().collect())
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_lagrangian_interpolator() {
        let x = vec![0.0, 1.0, 2.0];
        let y = vec![[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [2.0, 3.0, 1.0]];

        let interpolator = Lagrangian::new(x, y, false);

        let x_new = vec![0.5, 1.5];

        let result = interpolator.predict(&x_new).expect("Failed to predict");

        let expected = vec![[2.875, 2.375, 0.75], [2.375, 1.875, 1.75]];

        assert_eq!(result, expected);
    }
}