Cubic Spline in Rust

main.rs

mod matrix;
mod solver;
mod spline;
use crate::spline::NaturalSpline;

use egui_plot::{Line, Plot, PlotPoints};

fn main() {
    println!("Hello, world!");

    let start = 1.0;
    let step = 0.25;
    let total_steps = 32;

    let x = [1.0, 3.0, 5.0, 8.0];
    let y = [2.0, 3.0, 9.0, 10.0];

    let spline = match NaturalSpline::init(&x, &y) {
        Ok(s) => s,
        Err(msg) => {
            println!("{msg}");
            return;
        }
    };

    let options = eframe::NativeOptions::default();
    eframe::run_simple_native("Spline", options, move |ctx, _frame| {
        egui::CentralPanel::default().show(ctx, |ui| {
            Plot::new("spline").show(ui, |plot_ui| {
                let points: PlotPoints = (0..total_steps + 1)
                    .map(|i| {
                        let x = start + i as f64 * step;
                        [x, spline.spline(x)]
                    })
                    .collect();
                plot_ui.line(Line::new(points));
            });
        });
    })
    .expect("Failed to run");
}

matrix.rs

use std::fmt;
use std::fmt::{Display, Formatter};
use std::ops::{Index, IndexMut};

pub struct Matrix {
    matrix: Vec<f64>,
    pub width: usize,
    pub height: usize,
    det: Option<f64>,
}

impl Index<(usize, usize)> for Matrix {
    type Output = f64;
    fn index(&self, (row, col): (usize, usize)) -> &f64 {
        &self.matrix[row * self.width + col]
    }
}

impl IndexMut<(usize, usize)> for Matrix {
    fn index_mut(&mut self, (row, col): (usize, usize)) -> &mut f64 {
        &mut self.matrix[row * self.width + col]
    }
}

impl Matrix {
    pub fn init(matrix: &[&[f64]], width: usize, height: usize) -> Result<Self, &'static str> {
        if width == 0 || height == 0 {
            return Err("Invalid dimensions.");
        }
        
        let mat = matrix.iter().flat_map(|row| row.iter().copied()).collect();
        Ok(Self {
            matrix: mat,
            width,
            height,
            det: None,
        })
    }

    pub fn new(width: usize, height: usize) -> Result<Self, &'static str> {
        if width == 0 || height == 0 {
            return Err("Invalid dimensions.");
        }
        
        Ok(Self {
            matrix: vec![0.0; width * height],
            width,
            height,
            det: None,
        })
    }

    pub fn minor(&self, col: usize, row: usize) -> Matrix {
        let mut m = Matrix::new(self.width - 1, self.height - 1).unwrap();
        let mut idx_i = 0;
        for i in 0..self.width {
            let mut idx_j = 0;
            if i == row {
                continue;
            }
            for j in 0..self.height {
                if j == col {
                    continue;
                }
                m[(idx_i, idx_j)] = self[(i, j)];
                idx_j += 1;
            }
            idx_i += 1;
        }
        m
    }

    pub fn determinant(&mut self) -> f64 {
        if let Some(d) = self.det {
            return d;
        }

        if self.width != self.height {
            panic!("Determinant can only be calculated for a square matrix!");
        }

        let det = if self.width == 1 {
            self[(0, 0)]
        } else if self.width == 2 {
            (self[(0, 0)] * self[(1, 1)]) - (self[(1, 0)] * self[(0, 1)])
        } else {
            let mut calculated_determinant = 0.0;
            for i in 0..self.width {
                let d = self[(0, i)] * self.minor(i, 0).determinant();
                if (i & 1) == 1 {
                    calculated_determinant -= d
                } else {
                    calculated_determinant += d;
                }
            }
            calculated_determinant
        };

        self.det = Some(det);
        det
    }

    pub fn multiply(&self, matrix: &Matrix) -> Result<Matrix, &'static str> {
        if matrix.width != self.height {
            return Err(
                "Column of this matrix is not equal to the row of the matrix being multiplied!",
            );
        }

        let mut mat = Matrix::new(self.width, matrix.height)?;

        for i in 0..self.width {
            for j in 0..matrix.height {
                let mut t = 0.0;
                for k in 0..self.width {
                    t += self[(i, k)] * matrix[(k, j)];
                }

                mat[(i, j)] = t;
            }
        }

        Ok(mat)
    }

    pub fn swap_row(&mut self, i: usize, j: usize) {
        for k in 0..self.width {
            let temp = self[(i, k)];
            self[(i, k)] = self[(j, k)];
            self[(j, k)] = temp;
        }
    }

    pub fn swap_col(&mut self, i: usize, j: usize) {
        for k in 0..self.height {
            let temp = self[(k, i)];
            self[(k, i)] = self[(k, j)];
            self[(k, j)] = temp;
        }
    }
}

impl Display for Matrix {
    fn fmt(&self, f: &mut Formatter<'_>) -> fmt::Result {
        for i in 0..self.height {
            for j in 0..self.width {
                write!(f, "{:^4} ", self[(i, j)])?;
            }
            writeln!(f)?;
        }
        Ok(())
    }
}

solver.rs

use crate::matrix::Matrix;

pub fn solve_systems_of_equations(
    coefficient_matrix: &Matrix,
    constant_matrix: &Matrix,
) -> Result<Matrix, &'static str> {
    
    if coefficient_matrix.width != coefficient_matrix.height {
        return Err("The coefficient matrix must be a square matrix!");
    }

    if coefficient_matrix.height != constant_matrix.height {
        return Err(
            "You need the same number of variables and equations in order to solve a systems of equation!",
        );
    }

    let mut working_matrix = Matrix::new(coefficient_matrix.width + 1, coefficient_matrix.height)?;

    for i in 0..coefficient_matrix.width {
        for j in 0..coefficient_matrix.height {
            working_matrix[(i, j)] = coefficient_matrix[(i, j)];
        }
    }

    for i in 0..coefficient_matrix.height {
        working_matrix[(i, coefficient_matrix.width)] = constant_matrix[(i, 0)];
    }

    for i in 0..coefficient_matrix.width {
        if working_matrix[(i, i)] == 0.0 {
            for k in (i + 1)..coefficient_matrix.width {
                if working_matrix[(k, i)].abs() > 0.0 {
                    working_matrix.swap_row(i, k);
                    break;
                }
            }
        }

        if working_matrix[(i, i)] == 0.0 {
            return Err("This systems of equations can not be solved!");
        }

        for j in (i + 1)..coefficient_matrix.width {
            let r = working_matrix[(j, i)] / working_matrix[(i, i)];
            for k in 0..(coefficient_matrix.width + 1) {
                working_matrix[(j, k)] -= r * working_matrix[(i, k)];
            }
        }
    }

    let mut solution_matrix = Matrix::new(1, coefficient_matrix.height)?;

    for i in (0..=working_matrix.height - 1).rev() {
        solution_matrix[(i, 0)] = working_matrix[(i, working_matrix.width - 1)];

        for j in (i + 1)..working_matrix.height {
            solution_matrix[(i, 0)] -= working_matrix[(i, j)] * solution_matrix[(j, 0)];
        }

        solution_matrix[(i, 0)] /= working_matrix[(i, i)];
    }

    Ok(solution_matrix)
}

spline.rs

use crate::matrix::Matrix;
use crate::solver::solve_systems_of_equations;

// https://www.youtube.com/watch?v=yOUst2672qo

pub struct SplineFunction {
    pub lower_range: f64,
    pub upper_range: f64,
    a: f64,
    b: f64,
    c: f64,
    d: f64,
}

impl SplineFunction {
    pub fn eval(&self, x: f64) -> f64 {
        (self.a * x * x * x) + (self.b * x * x) + (self.c * x) + self.d
    }
}

fn setup_coefficients_for_nth_equation(
    matrix: &mut Matrix,
    a: f64,
    b: f64,
    c: f64,
    d: f64,
    n: usize,
    i: usize,
) {
    matrix[(n, i * 4)] = a;
    matrix[(n, i * 4 + 1)] = b;
    matrix[(n, i * 4 + 2)] = c;
    matrix[(n, i * 4 + 3)] = d;
}

pub struct NaturalSpline {
    functions: Vec<SplineFunction>,
}

impl NaturalSpline {
    pub fn init(x: &[f64], y: &[f64]) -> Result<Self, &'static str> {
        if x.len() != y.len() {
            return Err("Data points are not of equal length!");
        }

        if x.is_empty() || y.is_empty() {
            return Err("Give me some data points!");
        }

        let number_of_equations = x.len() - 1;
        let mut coefficient_matrix = Matrix::new(number_of_equations * 4, number_of_equations * 4)?;
        let mut constant_matrix = Matrix::new(1, number_of_equations * 4)?;

        let mut idx = 0;
        for (i, v) in y.iter().enumerate().take((number_of_equations + 1)) {
            constant_matrix[(0, idx)] = y[i];
            if i > 0 && i != number_of_equations {
                idx += 1;
                constant_matrix[(0, idx)] = *v;
            }
            idx += 1;
        }

        let mut offset = 0;
        let mut n = 0;
        idx = 0;
        for (i, v) in x.iter().enumerate().take(number_of_equations + 1) {
            let coeff = *v;
            setup_coefficients_for_nth_equation(
                &mut coefficient_matrix,
                coeff * coeff * coeff,
                coeff * coeff,
                coeff,
                1.0,
                idx,
                n,
            );
            if i != 0 && i != number_of_equations {
                idx += 1;
                n += 1;
                setup_coefficients_for_nth_equation(
                    &mut coefficient_matrix,
                    coeff * coeff * coeff,
                    coeff * coeff,
                    coeff,
                    1.0,
                    idx,
                    n,
                );
            }
            idx += 1;
            offset = idx - 1;
        }

        n = 0;
        for v in x.iter().take(number_of_equations).skip(1) {
            let coeff = *v;
            offset += 1;
            setup_coefficients_for_nth_equation(
                &mut coefficient_matrix,
                -3.0 * coeff * coeff,
                -2.0 * coeff,
                -1.0,
                0.0,
                offset,
                n,
            );
            n += 1;
            setup_coefficients_for_nth_equation(
                &mut coefficient_matrix,
                3.0 * coeff * coeff,
                2.0 * coeff,
                1.0,
                0.0,
                offset,
                n,
            );
        }

        n = 0;
        for v in x.iter().take(number_of_equations).skip(1) {
            let coeff = *v;
            offset += 1;
            setup_coefficients_for_nth_equation(
                &mut coefficient_matrix,
                -6.0 * coeff,
                -2.0,
                0.0,
                0.0,
                offset,
                n,
            );
            n += 1;
            setup_coefficients_for_nth_equation(
                &mut coefficient_matrix,
                6.0 * coeff,
                2.0,
                0.0,
                0.0,
                offset,
                n,
            );
        }

        offset += 1;
        setup_coefficients_for_nth_equation(
            &mut coefficient_matrix,
            6.0 * x[0],
            2.0,
            0.0,
            0.0,
            offset,
            0,
        );
        
        offset += 1;
        setup_coefficients_for_nth_equation(
            &mut coefficient_matrix,
            6.0 * x.last().unwrap(),
            2.0,
            0.0,
            0.0,
            offset,
            number_of_equations - 1,
        );

        let solution = solve_systems_of_equations(&coefficient_matrix, &constant_matrix)?;

        let mut functions: Vec<SplineFunction> = Vec::new();
        for i in 0..number_of_equations {
            functions.push(SplineFunction {
                lower_range: x[i],
                upper_range: x[i + 1],
                a: solution[(0, i * 4)],
                b: solution[(0, i * 4 + 1)],
                c: solution[(0, i * 4 + 2)],
                d: solution[(0, i * 4 + 3)],
            })
        }

        Ok(Self { functions })
    }

    pub fn spline(&self, x: f64) -> f64 {
        for f in &self.functions {
            if x >= f.lower_range && x <= f.upper_range {
                return f.eval(x);
            }
        }

        if x < self.functions[0].lower_range {
            return self.functions[0].eval(x);
        }

        self.functions.last().unwrap().eval(x)
    }
}