From 07ad4355854c384aaa8ea4ae1e2ffaeca7f2a2d1 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Tue, 12 Aug 2025 12:39:35 +0100 Subject: [PATCH 01/19] Cubic-Splines: Implementation of the cubic splines interpolation --- .../src/interpolation/cubic_spline.rs | 417 ++++++++++++++++++ 1 file changed, 417 insertions(+) create mode 100644 crates/RustQuant_math/src/interpolation/cubic_spline.rs diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs new file mode 100644 index 00000000..1dd36e59 --- /dev/null +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -0,0 +1,417 @@ +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// RustQuant: A Rust library for quantitative finance tools. +// Copyright (C) 2023 https://github.com/avhz +// Dual licensed under Apache 2.0 and MIT. +// See: +// - LICENSE-APACHE.md +// - LICENSE-MIT.md +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +//! Module containing functionality for interpolation. + + +use crate::interpolation::{InterpolationIndex, InterpolationValue, Interpolator, IntoValue}; +use RustQuant_error::RustQuantError; + +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// STRUCTS & ENUMS +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +/// Cubic Spline Interpolator. +pub struct CubicSplineInterpolator +where + IndexType: InterpolationIndex, + ValueType: InterpolationValue, +{ + /// X-axis values for the interpolator. + pub xs: Vec, + + /// Y-axis values for the interpolator. + pub ys: Vec, + + /// Whether the interpolator has been fitted. + pub fitted: bool, + + /// The second derivative of the spline at each point. + pub second_derivatives: Vec, + + /// The time steps between each point. + pub time_steps: Vec, +} + +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// IMPLEMENTATIONS, FUNCTIONS, AND MACROS +// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +impl CubicSplineInterpolator +where + IndexType: InterpolationIndex, + ValueType: InterpolationValue, +{ + /// Create a new CubicSplineInterpolator. + /// + /// # Errors + /// - `RustQuantError::UnequalLength` if ```xs.length() != ys.length()```. + /// + /// # Panics + /// Panics if NaN is in the index. + pub fn new( + xs: Vec, + ys: Vec, + ) -> Result, RustQuantError> { + if xs.len() != ys.len() { + return Err(RustQuantError::UnequalLength); + } + + let mut tmp: Vec<_> = xs.into_iter().zip(ys).collect(); + + tmp.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap()); + + let (xs, ys): (Vec, Vec) = tmp.into_iter().unzip(); + + Ok(Self { + xs, + ys, + fitted: false, + second_derivatives: vec![], + time_steps: vec![], + }) + } + + // Compute L and D matrices in A = L * D * L^T + // Computation of L optimized by tracking sub-diagonal entries + // as diagonal entries of L are always 1 + // and all other entries outside the sub-diagonal are 0 + // + // E.g. the matrix + // [ 1 0 0 0 ] + // [ a 1 0 0 ] + // [ 0 b 1 0 ] + // [ 0 0 c 1 ] + // is represented as [a, b, c] + fn compute_lower_tri_and_diagonal_matrices( + &self, + time_steps: &[IndexType::Delta] + ) -> ( + Vec, + Vec + ) { + let mut diagonal: Vec = time_steps + .windows(2) + .map(|t| { + (t[0] + t[1]) * ValueType::from_f64(2.0).unwrap() + }).collect(); + + let mut lower_tri_matrix_sub_diag_entries: Vec = Vec::with_capacity(diagonal.len() - 1); + + for i in 1..diagonal.len() { + let prev_diag: IndexType::Delta = diagonal[i - 1]; + let t_prev: IndexType::Delta = time_steps[i - 1]; + + diagonal[i] = diagonal[i] - t_prev * (t_prev / prev_diag); + lower_tri_matrix_sub_diag_entries.push(t_prev / prev_diag); + } + + (diagonal, lower_tri_matrix_sub_diag_entries) + } + + // L * D * L^T inverts to (L^T)^-1 * D^-1 * L^-1 + // This method computes L^-1 + // + // The inverse of a lower triangular matrix + // is a lower triangular matrix + // represented as a vector of vectors + // without the known 0 entries + // i.e. + // + // [a 0 0 0] + // [b c 0 0] + // [d e f 0] + // [g h i j] + // + // is represented as + // + // [[a], [b, c], [d, e, f], [g, h, i, j]] + fn invert_lower_tri_matrix( + &self, + lower_tri_matrix_sub_diag_entries: &[ValueType] + ) -> Vec> { + + let mut temp: ValueType; + let mut inverse: Vec> = vec![vec![]; lower_tri_matrix_sub_diag_entries.len() + 1]; + + inverse[0].push(ValueType::one()); + + for i in 0..lower_tri_matrix_sub_diag_entries.len() { + temp = ValueType::one(); + for j in i..lower_tri_matrix_sub_diag_entries.len() { + temp *= - lower_tri_matrix_sub_diag_entries[j]; + inverse[j + 1].push(temp) + } + inverse[i + 1].push(ValueType::one()); + } + + inverse + } + + // Transpose a matrix + fn transpose( + &self, matrix: Vec> + ) -> Vec> { + let mut transposed_matrix: Vec> = vec![]; + + for i in 0..matrix.len() { + let mut row: Vec = vec![]; + for j in i..matrix.len() { + row.push(matrix[j][i]); + } + transposed_matrix.push(row); + } + transposed_matrix + } + + // L * D * L^T inverts to (L^T)^-1 * D^-1 * L^-1 + // This method computes (L^T)^-1 * D^-1 + // + // The inverse of an upper triangular matrix + // will be a upper triangular matrix + // represented as an vector of vectors + // without the known 0 entries + // i.e. + // + // [a b c d] + // [0 e f g] + // [0 0 h i] + // [0 0 0 j] + // + // is represented as + // + // [[a, b, c, d], [e, f, g], [h, i], [j]] + // + // Note that the method undertaken here takes + // advantage of the fact that calculating the + // inverse of an upper triangular matrix is equal + // to calculating the inverse its transpose + // (a lower triangular matrix) then transposing the result. + fn lower_tri_transpose_times_diag_inv( + &self, + lower_tri_inverse: &[Vec], + diagonal: &[IndexType::Delta], + ) -> Vec> { + let mut upper_tri_row: Vec = vec![]; + let mut product: Vec> = vec![]; + let one: ValueType = ValueType::one(); + + for i in 0..diagonal.len() { + upper_tri_row.clear(); + for j in 0..(i + 1) { + if i == j { + upper_tri_row.push(one / diagonal[j].into_value()) + } else { + upper_tri_row.push(lower_tri_inverse[i][j] / diagonal[j].into_value()); + } + } + product.push(upper_tri_row.clone()) + } + + self.transpose(product) + } + + // Multiply (L^T)^-1 * D^-1 (Upper triangular matrix) with L^-1 (Lower triangular matrix) + // to finally get the inverse of A = L * D * L^T + // The product is represented as a vector of vectors. + // + // Note that the *transpose* of the lower triangular matrix + // is used to simplify the matrix multiplication process + fn upper_tri_times_lower_tri( + &self, + upper_tri_matrix: &[Vec], + lower_tri_matrix_transpose: &[Vec] + ) -> Vec> + { + let mut product: Vec> = vec![vec![]; upper_tri_matrix.len()]; + let mut matrix_entry: ValueType; + + for i in 0..upper_tri_matrix.len() { + for j in 0..lower_tri_matrix_transpose.len() { + matrix_entry = ValueType::zero(); + let lower_diff: usize = (i as i64 - j as i64).max(0) as usize; + let upper_diff: usize = (j as i64 - i as i64).max(0) as usize; + for k in 0..(upper_tri_matrix[i].len().min(lower_tri_matrix_transpose[j].len())) { + matrix_entry = matrix_entry + upper_tri_matrix[i][k + upper_diff] * lower_tri_matrix_transpose[j][k + lower_diff]; + } + product[i].push(matrix_entry); + } + } + product + } + + // Compute the spline value at a given point. + fn spline( + &self, + point: IndexType, + time_step: IndexType::Delta, + x_coordinate_lower: IndexType, + x_coordinate_upper: IndexType, + y_coordinate_lower: ValueType, + y_coordinate_upper: ValueType, + second_derivative_lower: ValueType, + second_derivative_upper: ValueType + ) -> ValueType { + let term_1: ValueType = ( + second_derivative_lower / (ValueType::from_f64(6.0).unwrap() + * time_step.into_value())) + * (x_coordinate_upper - point).into_value() + * (x_coordinate_upper - point).into_value() + * (x_coordinate_upper - point).into_value(); + + let term_2: ValueType = ( + second_derivative_upper / (ValueType::from_f64(6.0).unwrap() + * time_step.into_value())) + * (point - x_coordinate_lower).into_value() + * (point - x_coordinate_lower).into_value() + * (point - x_coordinate_lower).into_value(); + + let term_3: ValueType = ( + (y_coordinate_upper / time_step.into_value()) - (time_step.into_value() + * second_derivative_upper / ValueType::from_f64(6.0).unwrap()) + ) * (point - x_coordinate_lower).into_value(); + + let term_4: ValueType = ( + (y_coordinate_lower / time_step.into_value()) - (time_step.into_value() + * second_derivative_lower / ValueType::from_f64(6.0).unwrap()) + ) * (x_coordinate_upper - point).into_value(); + + term_1 + term_2 + term_3 + term_4 + } +} + +impl Interpolator + for CubicSplineInterpolator +where + IndexType: InterpolationIndex, + ValueType: InterpolationValue, +{ + // The first step in fitting the cubic spline interpolator + // is to compute the 2nd derivative at each point of the spline. + // This is done by solving a system of linear equations Ax = b + // where A is a tridiagonal matrix, b is a vector of values and + // x is the vector of second derivatives at each point. + // + // The system is solved by decomposing A into L * D * L^T + // and subsequently computing the inverse of A. + // + // Once all the second derivatives are computed, we can then + // compute the spline value at a given point. + fn fit(&mut self) -> Result<(), RustQuantError> { + + self.time_steps = self.xs.windows(2) + .map(|x| (x[1] - x[0])) + .collect(); + + let (diag_matrix, lower_tri_matrix) = self.compute_lower_tri_and_diagonal_matrices(&self.time_steps); + + let mut inv_lower_tri_matrix: Vec> = self.invert_lower_tri_matrix(&lower_tri_matrix); + + let upper_tri_matrix: Vec> = self.lower_tri_transpose_times_diag_inv( + &inv_lower_tri_matrix, + &diag_matrix, + ); + + inv_lower_tri_matrix = self.transpose(inv_lower_tri_matrix); + let tridiagonal_inverse: Vec> = self.upper_tri_times_lower_tri( + &upper_tri_matrix, + &inv_lower_tri_matrix + ); + + let rhs_vector: Vec = self.xs.iter().zip(self.ys.iter()) + .collect::>() + .windows(3) + .map(|window| { + let x0: &IndexType = window[0].0; + let x1: &IndexType = window[1].0; + let x2: &IndexType = window[2].0; + + let y0: &ValueType = window[0].1; + let y1: &ValueType = window[1].1; + let y2: &ValueType = window[2].1; + + ValueType::from_f64(6.0).unwrap() * ( + (ValueType::one() / (*x2 - *x1).into_value()) * (*y2) + - ( + (ValueType::one() / (*x2 - *x1).into_value()) + + (ValueType::one() / (*x1 - *x0).into_value()) + ) * (*y1) + + (ValueType::one() / (*x1 - *x0).into_value()) * (*y0) + ) + }).collect(); + + self.second_derivatives = vec![ValueType::zero(); tridiagonal_inverse.len()]; + let mut matrix_entry: ValueType; + + for i in 0..tridiagonal_inverse.len() { + matrix_entry = self.second_derivatives[i]; + for j in 0..tridiagonal_inverse[i].len() { + matrix_entry = matrix_entry + tridiagonal_inverse[i][j] * rhs_vector[j]; + } + self.second_derivatives[i] = matrix_entry; + } + + self.second_derivatives.insert(0, ValueType::zero()); + self.second_derivatives.push(ValueType::zero()); + + self.fitted = true; + + Ok(()) + } + + fn range(&self) -> (IndexType, IndexType) { + (*self.xs.first().unwrap(), *self.xs.last().unwrap()) + } + + fn add_point(&mut self, point: (IndexType, ValueType)) { + let idx = self.xs.partition_point(|&x| x < point.0); + self.xs.insert(idx, point.0); + self.ys.insert(idx, point.1); + let _ = self.fit(); + } + + fn interpolate(&self, point: IndexType) -> Result { + let range = self.range(); + if point.partial_cmp(&range.0).unwrap() == std::cmp::Ordering::Less + || point.partial_cmp(&range.1).unwrap() == std::cmp::Ordering::Greater + { + return Err(RustQuantError::OutsideOfRange); + } + if let Ok(idx) = self + .xs + .binary_search_by(|p| p.partial_cmp(&point).expect("Cannot compare values.")) + { + return Ok(self.ys[idx]); + } + + let mut result: Option = None; + for k in 0..(self.xs.len() - 1) { + if self.xs[k] <= point && point < self.xs[k + 1] { + result = Some( + self.spline( + point, + self.time_steps[k], + self.xs[k], + self.xs[k + 1], + self.ys[k], + self.ys[k + 1], + self.second_derivatives[k], + self.second_derivatives[k + 1], + ) + ); + break; + } + } + + match result { + Some(value) => Ok(value), + None => Err(RustQuantError::OutsideOfRange), + } + } +} From 0785221709f6ce766c51116b2008694b41f12cf9 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 17 Aug 2025 21:33:40 +0100 Subject: [PATCH 02/19] Cubic-Splines: Add interpolation unit tests --- .../src/interpolation/cubic_spline.rs | 77 +++++++++++++++++++ 1 file changed, 77 insertions(+) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 1dd36e59..d3c91bbe 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -415,3 +415,80 @@ where } } } + +#[cfg(test)] +mod tests_cubic_spline_interpolation { + use super::*; + use RustQuant_utils::{assert_approx_equal, RUSTQUANT_EPSILON}; + + #[test] + fn test_natural_cubic_interpolation() { + + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + assert_approx_equal!( + 36.64909523809524, + interpolator.interpolate(2.5).unwrap(), + RUSTQUANT_EPSILON + ); + } + + #[test] + fn test_natural_cubic_interpolation_dates() { + let now: time::OffsetDateTime = time::OffsetDateTime::now_utc(); + + let xs: Vec = vec![ + now, + now + time::Duration::days(1), + now + time::Duration::days(2), + now + time::Duration::days(3), + now + time::Duration::days(4), + ]; + + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs.clone(), ys).unwrap(); + let _ = interpolator.fit(); + + assert_approx_equal!( + 36.64909523809524, + interpolator + .interpolate(xs[2] + time::Duration::hours(12)) + .unwrap(), + RUSTQUANT_EPSILON + ); + } + + #[test] + fn test_cubic_interpolation_out_of_range() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys = vec![1., 2., 3., 4., 5.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + assert!(interpolator.interpolate(6.).is_err()); + } + + #[test] + fn test_cubic_interpolation_add_range_and_refit() { + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + interpolator.add_point((5.0, 625.0)); + let _ = interpolator.fit(); + + assert_approx_equal!( + 39.966361318634604, + interpolator.interpolate(2.5).unwrap(), + RUSTQUANT_EPSILON + ); + } +} From d292c621c5102ce0d915d162bf684a5c36b0a092 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 17 Aug 2025 21:34:43 +0100 Subject: [PATCH 03/19] Cubic-Splines: Add helper functions unit tests --- .../src/interpolation/cubic_spline.rs | 129 ++++++++++++++++++ 1 file changed, 129 insertions(+) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index d3c91bbe..52cf8d01 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -492,3 +492,132 @@ mod tests_cubic_spline_interpolation { ); } } + +#[cfg(test)] +mod tests_cubic_spline_helper_functions { + use super::*; + + #[test] + fn test_compute_lower_tri_and_diagonal_matrices() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let time_steps: Vec = xs.windows(2) + .map(|x| (x[1] - x[0])) + .collect(); + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let (lower_triangular, diagonal) = interpolator.compute_lower_tri_and_diagonal_matrices(&time_steps); + + assert!(lower_triangular == vec![4.0, 3.75, 3.7333333333333334]); + assert!(diagonal == vec![0.25, 0.26666666666666666]); + } + + #[test] + fn test_invert_lower_tri_matrix() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let inv_lower_tri_matrix: Vec> = interpolator.invert_lower_tri_matrix(&[3.0, 2.0, 1.0]); + + assert!( + inv_lower_tri_matrix == vec![ + vec![1.0], + vec![-3.0, 1.0], + vec![6.0, -2.0, 1.0], + vec![-6.0, 2.0, -1.0, 1.0] + ] + ) + } + + #[test] + fn test_transpose() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let transposed_matrix = interpolator.transpose(vec![ + vec![1.0], + vec![2.0, 3.0], + vec![4.0, 5.0, 6.0], + ]); + + assert!( + transposed_matrix == vec![ + vec![1.0, 2.0, 4.0], + vec![3.0, 5.0], + vec![6.0], + ] + ); + } + + #[test] + fn lower_tri_transpose_times_diag_inv() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let lower_tri_inverse: Vec> = vec![vec![1.0], vec![2.0, 3.0], vec![4.0, 5.0, 6.0]]; + let diagonal: Vec = vec![1.0, 2.0, 3.0]; + + assert!( + interpolator.lower_tri_transpose_times_diag_inv(&lower_tri_inverse, &diagonal) == vec![ + vec![1.0, 2.0, 4.0], + vec![0.5, 2.5], + vec![0.3333333333333333] + ] + ); + } + + #[test] + fn test_upper_tri_times_lower_tri() { + let xs: Vec = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let upper_tri_matrix = vec![ + vec![1.0, 2.0, 3.0], + vec![4.0, 5.0], + vec![6.0], + ]; + let lower_tri_matrix_transpose = vec![ + vec![-1.0, -2.0, -3.0], + vec![-4.0, -5.0], + vec![-6.0], + ]; + + let product: Vec> = interpolator.upper_tri_times_lower_tri( + &upper_tri_matrix, + &lower_tri_matrix_transpose + ); + + assert!( + product == vec![ + [-14.0, -23.0, -18.0], + [-23.0, -41.0, -30.0], + [-18.0, -30.0, -36.0] + ] + ); + } + + #[test] + fn test_spline() { + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let spline_result: f64 = interpolator.spline( + 2.5, + 1.0, + 2.0, + 3.0, + 16.0, + 81.0, + 32.75733333333333, + 156.85714285714286 + ); + + assert!(spline_result == 36.64909523809524); + } +} From b4869d7cffaff19cf9c224677984b5180590257c Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 17 Aug 2025 21:37:31 +0100 Subject: [PATCH 04/19] Cubic-Splines: Modifications to InterpolationIndex and InterpolationValue traits + introduce the IntoValue trait --- .../RustQuant_math/src/interpolation/mod.rs | 90 +++++++++++++++---- 1 file changed, 72 insertions(+), 18 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/mod.rs b/crates/RustQuant_math/src/interpolation/mod.rs index d3b69d9a..a2d92274 100644 --- a/crates/RustQuant_math/src/interpolation/mod.rs +++ b/crates/RustQuant_math/src/interpolation/mod.rs @@ -7,7 +7,7 @@ // - LICENSE-MIT.md // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -use std::ops::{AddAssign, Div, Mul, Sub}; +use std::ops::{Div, Mul, MulAssign, Sub, Add, AddAssign, Neg}; use RustQuant_error::RustQuantError; pub mod linear_interpolator; @@ -19,22 +19,46 @@ pub use exponential_interpolator::*; pub mod b_splines; pub use b_splines::*; +pub mod cubic_spline; +pub use cubic_spline::*; + + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ /// Trait describing requirements to be interpolated. -pub trait InterpolationValue: - num::Num + AddAssign + std::fmt::Debug + Copy + Clone + Sized -{ +pub trait InterpolationValue: num::Num + + num::FromPrimitive + + Neg + + AddAssign + + MulAssign + + Copy + + Clone + + Sized + + std::fmt::Display + + std::fmt::Debug + + Send + + Sync {} + +/// Trait to convert a Delta type into a value of `ValueType`. +pub trait IntoValue { + /// Convert `self` into `ValueType`. + fn into_value(self) -> ValueType; } /// Trait describing requirements to be an index of interpolation. pub trait InterpolationIndex: - Sub + PartialOrd + Copy + Clone + Sized + std::fmt::Display + Sub + PartialOrd + Copy + Clone + Sized + std::fmt::Display + Send + Sync { /// Type of the difference of `Self` - `Self` type Delta: Div - + Mul; + + Mul + + Add + + Sub + + IntoValue + + Copy + + Send + + Sync; /// Type of `Delta` / `Delta` type DeltaDiv: InterpolationValue; @@ -66,10 +90,18 @@ where fn add_point(&mut self, point: (IndexType, ValueType)); } -impl InterpolationValue for T where - T: num::Num + AddAssign + std::fmt::Debug + Copy + Clone + Sized -{ -} +impl InterpolationValue for T where T: num::Num + + num::FromPrimitive + + Neg + + AddAssign + + MulAssign + + Copy + + Clone + + Sized + + std::fmt::Display + + std::fmt::Debug + + Send + + Sync {} macro_rules! impl_interpolation_index { ($a:ty, $b:ty, $c:ty) => { @@ -80,31 +112,52 @@ macro_rules! impl_interpolation_index { }; } +macro_rules! impl_num_delta_into_value { + ($b:ty, $c:ty) => { + impl IntoValue<$c> for $b { + fn into_value(self) -> $c { + self as $c + } + } + }; +} + +macro_rules! impl_time_delta_into_value { + ($b:ty, $c:ty) => { + impl IntoValue<$c> for $b { + fn into_value(self) -> $c { + self.as_seconds_f64() + } + } + }; +} + // Implement InterpolationIndex for all signed integer types. impl_interpolation_index!(i8, i8, i8); +impl_num_delta_into_value!(i8, i8); impl_interpolation_index!(i16, i16, i16); +impl_num_delta_into_value!(i16, i16); impl_interpolation_index!(i32, i32, i32); +impl_num_delta_into_value!(i32, i32); impl_interpolation_index!(i64, i64, i64); +impl_num_delta_into_value!(i64, i64); impl_interpolation_index!(i128, i128, i128); +impl_num_delta_into_value!(i128, i128); impl_interpolation_index!(isize, isize, isize); - -// Implement InterpolationIndex for all unsigned integer types. -impl_interpolation_index!(u8, u8, u8); -impl_interpolation_index!(u16, u16, u16); -impl_interpolation_index!(u32, u32, u32); -impl_interpolation_index!(u64, u64, u64); -impl_interpolation_index!(u128, u128, u128); -impl_interpolation_index!(usize, usize, usize); +impl_num_delta_into_value!(isize, isize); // Implement InterpolationIndex for all floating point types. impl_interpolation_index!(f32, f32, f32); +impl_num_delta_into_value!(f32, f32); impl_interpolation_index!(f64, f64, f64); +impl_num_delta_into_value!(f64, f64); // Implement InterpolationIndex for date/time types. impl_interpolation_index!(time::Date, time::Duration, f64); impl_interpolation_index!(time::Time, time::Duration, f64); impl_interpolation_index!(time::OffsetDateTime, time::Duration, f64); impl_interpolation_index!(time::PrimitiveDateTime, time::Duration, f64); +impl_time_delta_into_value!(time::Duration, f64); // Implement InterpolationIndex for Decimal type. impl_interpolation_index!( @@ -112,3 +165,4 @@ impl_interpolation_index!( rust_decimal::Decimal, rust_decimal::Decimal ); +impl_num_delta_into_value!(rust_decimal::Decimal, rust_decimal::Decimal); From 32f9de6fd6cd8cd9ad8d76e8f6d02c5b1e436b77 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 31 Aug 2025 00:49:18 +0100 Subject: [PATCH 05/19] Cubic-Splines: Ensure interpolator is fitted when interpolating otherwise raise error --- .../src/interpolation/cubic_spline.rs | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 52cf8d01..41b4cfea 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -377,6 +377,9 @@ where } fn interpolate(&self, point: IndexType) -> Result { + if !self.fitted { + return Err(RustQuantError::Unfitted); + } let range = self.range(); if point.partial_cmp(&range.0).unwrap() == std::cmp::Ordering::Less || point.partial_cmp(&range.1).unwrap() == std::cmp::Ordering::Greater @@ -491,6 +494,22 @@ mod tests_cubic_spline_interpolation { RUSTQUANT_EPSILON ); } + + #[test] + fn test_cubic_spline_unfitted() { + + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + + assert!( + matches!( + interpolator.interpolate(2.5), + Err(RustQuantError::Unfitted), + ) + ); + } } #[cfg(test)] From d48cd946d06ac879070d9a6a40c9a4e50dd142e0 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 31 Aug 2025 00:54:26 +0100 Subject: [PATCH 06/19] Cubic-Splines: Naming amendments for clarity + remove redundant variable --- .../src/interpolation/cubic_spline.rs | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 41b4cfea..a155c281 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -193,25 +193,24 @@ where // inverse of an upper triangular matrix is equal // to calculating the inverse its transpose // (a lower triangular matrix) then transposing the result. - fn lower_tri_transpose_times_diag_inv( + fn lower_tri_transpose_inv_times_diag_inv( &self, lower_tri_inverse: &[Vec], diagonal: &[IndexType::Delta], ) -> Vec> { - let mut upper_tri_row: Vec = vec![]; + let mut product_row: Vec = vec![]; let mut product: Vec> = vec![]; - let one: ValueType = ValueType::one(); for i in 0..diagonal.len() { - upper_tri_row.clear(); + product_row.clear(); for j in 0..(i + 1) { if i == j { - upper_tri_row.push(one / diagonal[j].into_value()) + product_row.push(ValueType::one() / diagonal[j].into_value()) } else { - upper_tri_row.push(lower_tri_inverse[i][j] / diagonal[j].into_value()); + product_row.push(lower_tri_inverse[i][j] / diagonal[j].into_value()); } } - product.push(upper_tri_row.clone()) + product.push(product_row.clone()) } self.transpose(product) @@ -313,7 +312,7 @@ where let mut inv_lower_tri_matrix: Vec> = self.invert_lower_tri_matrix(&lower_tri_matrix); - let upper_tri_matrix: Vec> = self.lower_tri_transpose_times_diag_inv( + let upper_tri_matrix: Vec> = self.lower_tri_transpose_inv_times_diag_inv( &inv_lower_tri_matrix, &diag_matrix, ); @@ -572,7 +571,7 @@ mod tests_cubic_spline_helper_functions { } #[test] - fn lower_tri_transpose_times_diag_inv() { + fn lower_tri_transpose_inv_times_diag_inv() { let xs: Vec = vec![1., 2., 3., 4., 5.]; let ys: Vec = vec![1., 2., 3., 4., 5.]; @@ -581,7 +580,7 @@ mod tests_cubic_spline_helper_functions { let diagonal: Vec = vec![1.0, 2.0, 3.0]; assert!( - interpolator.lower_tri_transpose_times_diag_inv(&lower_tri_inverse, &diagonal) == vec![ + interpolator.lower_tri_transpose_inv_times_diag_inv(&lower_tri_inverse, &diagonal) == vec![ vec![1.0, 2.0, 4.0], vec![0.5, 2.5], vec![0.3333333333333333] From 74684200b31d55d87c1e9bfb9970bb4ecf1a4edf Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Mon, 1 Sep 2025 18:57:25 +0100 Subject: [PATCH 07/19] Cubic-Spline: Resolve Clippy warnings --- .../src/interpolation/cubic_spline.rs | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index a155c281..4f64cba2 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -161,11 +161,9 @@ where let mut transposed_matrix: Vec> = vec![]; for i in 0..matrix.len() { - let mut row: Vec = vec![]; - for j in i..matrix.len() { - row.push(matrix[j][i]); - } - transposed_matrix.push(row); + transposed_matrix.push( + (i..matrix.len()).map(|j| matrix[j][i]).collect() + ); } transposed_matrix } @@ -201,6 +199,7 @@ where let mut product_row: Vec = vec![]; let mut product: Vec> = vec![]; + #[allow(clippy::needless_range_loop)] for i in 0..diagonal.len() { product_row.clear(); for j in 0..(i + 1) { @@ -232,12 +231,12 @@ where let mut matrix_entry: ValueType; for i in 0..upper_tri_matrix.len() { - for j in 0..lower_tri_matrix_transpose.len() { + for (j, lower_tri_matrix_transpose_entry) in lower_tri_matrix_transpose.iter().enumerate() { matrix_entry = ValueType::zero(); let lower_diff: usize = (i as i64 - j as i64).max(0) as usize; let upper_diff: usize = (j as i64 - i as i64).max(0) as usize; - for k in 0..(upper_tri_matrix[i].len().min(lower_tri_matrix_transpose[j].len())) { - matrix_entry = matrix_entry + upper_tri_matrix[i][k + upper_diff] * lower_tri_matrix_transpose[j][k + lower_diff]; + for k in 0..(upper_tri_matrix[i].len().min(lower_tri_matrix_transpose_entry.len())) { + matrix_entry += upper_tri_matrix[i][k + upper_diff] * lower_tri_matrix_transpose_entry[k + lower_diff]; } product[i].push(matrix_entry); } @@ -246,6 +245,7 @@ where } // Compute the spline value at a given point. + #[allow(clippy::too_many_arguments)] fn spline( &self, point: IndexType, @@ -348,10 +348,10 @@ where self.second_derivatives = vec![ValueType::zero(); tridiagonal_inverse.len()]; let mut matrix_entry: ValueType; - for i in 0..tridiagonal_inverse.len() { + for (i, tridiagonal_inverse_row) in tridiagonal_inverse.iter().enumerate() { matrix_entry = self.second_derivatives[i]; - for j in 0..tridiagonal_inverse[i].len() { - matrix_entry = matrix_entry + tridiagonal_inverse[i][j] * rhs_vector[j]; + for (j, tridiagonal_inverse_entry) in tridiagonal_inverse_row.iter().enumerate() { + matrix_entry += *tridiagonal_inverse_entry * rhs_vector[j]; } self.second_derivatives[i] = matrix_entry; } From b93660ee77660e581afd525730c5e8d48f7769f7 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Tue, 2 Sep 2025 00:45:17 +0100 Subject: [PATCH 08/19] Cubic-Spline: Amend add_point() to not fit interpolator and set fitted attribute to false --- crates/RustQuant_math/src/interpolation/cubic_spline.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 4f64cba2..2abf7f8d 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -372,7 +372,7 @@ where let idx = self.xs.partition_point(|&x| x < point.0); self.xs.insert(idx, point.0); self.ys.insert(idx, point.1); - let _ = self.fit(); + self.fitted = false; } fn interpolate(&self, point: IndexType) -> Result { From c5053a643f6a57265a1fd1d659bca357ea208f62 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Tue, 2 Sep 2025 00:45:53 +0100 Subject: [PATCH 09/19] Cubic-Spline: Unit test for adding a point and not refitting --- .../src/interpolation/cubic_spline.rs | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 2abf7f8d..dd49f274 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -494,6 +494,24 @@ mod tests_cubic_spline_interpolation { ); } + #[test] + fn test_cubic_interpolation_add_range_and_no_refit() { + let xs: Vec = vec![0., 1., 2., 3., 4.]; + let ys: Vec = vec![0., 1., 16., 81., 256.]; + + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + let _ = interpolator.fit(); + + interpolator.add_point((5.0, 625.0)); + + assert!( + matches!( + interpolator.interpolate(2.5), + Err(RustQuantError::Unfitted), + ) + ); + } + #[test] fn test_cubic_spline_unfitted() { From 855f2e8aa19b946970e10551b7f7859ac008f0f9 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Tue, 2 Sep 2025 19:23:36 +0100 Subject: [PATCH 10/19] Cubic-Spline: Amend IntoValue trait implementation for numerical data types --- crates/RustQuant_math/src/interpolation/mod.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/crates/RustQuant_math/src/interpolation/mod.rs b/crates/RustQuant_math/src/interpolation/mod.rs index a2d92274..ff90115c 100644 --- a/crates/RustQuant_math/src/interpolation/mod.rs +++ b/crates/RustQuant_math/src/interpolation/mod.rs @@ -116,7 +116,7 @@ macro_rules! impl_num_delta_into_value { ($b:ty, $c:ty) => { impl IntoValue<$c> for $b { fn into_value(self) -> $c { - self as $c + self } } }; From bd0c0acc3b29922ec9f82fd14d1cc3780d587c44 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sat, 4 Oct 2025 19:30:04 +0100 Subject: [PATCH 11/19] Cubic-Spline: Rewrite Interpolator by replacing LDLT decomposition with Cholesky Decomposition --- .../src/interpolation/cubic_spline.rs | 213 +++--------------- 1 file changed, 36 insertions(+), 177 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index dd49f274..e8e2beb8 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -78,170 +78,52 @@ where }) } - // Compute L and D matrices in A = L * D * L^T - // Computation of L optimized by tracking sub-diagonal entries - // as diagonal entries of L are always 1 - // and all other entries outside the sub-diagonal are 0 - // - // E.g. the matrix - // [ 1 0 0 0 ] - // [ a 1 0 0 ] - // [ 0 b 1 0 ] - // [ 0 0 c 1 ] - // is represented as [a, b, c] - fn compute_lower_tri_and_diagonal_matrices( - &self, - time_steps: &[IndexType::Delta] - ) -> ( - Vec, - Vec - ) { - let mut diagonal: Vec = time_steps - .windows(2) - .map(|t| { - (t[0] + t[1]) * ValueType::from_f64(2.0).unwrap() - }).collect(); + fn cholesky_decomposition(&self) -> Vec> { + let n: usize = self.time_steps.len(); + let mut lower_tri_matrix: Vec> = vec![vec![]; n - 1]; - let mut lower_tri_matrix_sub_diag_entries: Vec = Vec::with_capacity(diagonal.len() - 1); + let mut prev_diag: ValueType = (ValueType::from_f64(2.0).unwrap() * (self.time_steps[0] + self.time_steps[1]).into_value()).square_root(); + lower_tri_matrix[0].push(prev_diag); + let mut prev_sub_diag: ValueType; - for i in 1..diagonal.len() { - let prev_diag: IndexType::Delta = diagonal[i - 1]; - let t_prev: IndexType::Delta = time_steps[i - 1]; + for (i, time) in self.time_steps[1..(self.time_steps.len() - 1)].iter().enumerate() { + prev_sub_diag = time.into_value() / prev_diag; + lower_tri_matrix[i + 1].push(prev_sub_diag); - diagonal[i] = diagonal[i] - t_prev * (t_prev / prev_diag); - lower_tri_matrix_sub_diag_entries.push(t_prev / prev_diag); + prev_diag = (ValueType::from_f64(2.0).unwrap() * (*time + self.time_steps[i + 1]).into_value() - (prev_sub_diag * prev_sub_diag)).square_root(); + lower_tri_matrix[i + 1].push(prev_diag) } - (diagonal, lower_tri_matrix_sub_diag_entries) + lower_tri_matrix } - // L * D * L^T inverts to (L^T)^-1 * D^-1 * L^-1 - // This method computes L^-1 - // - // The inverse of a lower triangular matrix - // is a lower triangular matrix - // represented as a vector of vectors - // without the known 0 entries - // i.e. - // - // [a 0 0 0] - // [b c 0 0] - // [d e f 0] - // [g h i j] - // - // is represented as - // - // [[a], [b, c], [d, e, f], [g, h, i, j]] - fn invert_lower_tri_matrix( - &self, - lower_tri_matrix_sub_diag_entries: &[ValueType] - ) -> Vec> { - - let mut temp: ValueType; - let mut inverse: Vec> = vec![vec![]; lower_tri_matrix_sub_diag_entries.len() + 1]; - - inverse[0].push(ValueType::one()); - - for i in 0..lower_tri_matrix_sub_diag_entries.len() { - temp = ValueType::one(); - for j in i..lower_tri_matrix_sub_diag_entries.len() { - temp *= - lower_tri_matrix_sub_diag_entries[j]; - inverse[j + 1].push(temp) - } - inverse[i + 1].push(ValueType::one()); + fn backward_substitution(&self, matrix: &[Vec], rhs: &[ValueType]) -> Vec { + + let matrix_len: usize = matrix.len(); + let mut prev_value: ValueType = rhs[matrix_len - 1] / matrix[matrix_len - 1][1]; + let mut result: Vec = vec![prev_value]; + + for i in (1..(matrix_len - 1)).rev() { + prev_value = (rhs[i] - matrix[i + 1][0] * prev_value) / matrix[i][1]; + result.insert(0, prev_value) } - inverse - } + result.insert(0, (rhs[0] - matrix[1][0] * prev_value) / matrix[0][0]); - // Transpose a matrix - fn transpose( - &self, matrix: Vec> - ) -> Vec> { - let mut transposed_matrix: Vec> = vec![]; + result + } - for i in 0..matrix.len() { - transposed_matrix.push( - (i..matrix.len()).map(|j| matrix[j][i]).collect() - ); - } - transposed_matrix - } + fn forward_substitution(&self, matrix: &[Vec], rhs: &[ValueType]) -> Vec { - // L * D * L^T inverts to (L^T)^-1 * D^-1 * L^-1 - // This method computes (L^T)^-1 * D^-1 - // - // The inverse of an upper triangular matrix - // will be a upper triangular matrix - // represented as an vector of vectors - // without the known 0 entries - // i.e. - // - // [a b c d] - // [0 e f g] - // [0 0 h i] - // [0 0 0 j] - // - // is represented as - // - // [[a, b, c, d], [e, f, g], [h, i], [j]] - // - // Note that the method undertaken here takes - // advantage of the fact that calculating the - // inverse of an upper triangular matrix is equal - // to calculating the inverse its transpose - // (a lower triangular matrix) then transposing the result. - fn lower_tri_transpose_inv_times_diag_inv( - &self, - lower_tri_inverse: &[Vec], - diagonal: &[IndexType::Delta], - ) -> Vec> { - let mut product_row: Vec = vec![]; - let mut product: Vec> = vec![]; - - #[allow(clippy::needless_range_loop)] - for i in 0..diagonal.len() { - product_row.clear(); - for j in 0..(i + 1) { - if i == j { - product_row.push(ValueType::one() / diagonal[j].into_value()) - } else { - product_row.push(lower_tri_inverse[i][j] / diagonal[j].into_value()); - } - } - product.push(product_row.clone()) - } + let mut prev_value: ValueType = rhs[0] / matrix[0][0]; + let mut result: Vec = vec![prev_value]; - self.transpose(product) - } - - // Multiply (L^T)^-1 * D^-1 (Upper triangular matrix) with L^-1 (Lower triangular matrix) - // to finally get the inverse of A = L * D * L^T - // The product is represented as a vector of vectors. - // - // Note that the *transpose* of the lower triangular matrix - // is used to simplify the matrix multiplication process - fn upper_tri_times_lower_tri( - &self, - upper_tri_matrix: &[Vec], - lower_tri_matrix_transpose: &[Vec] - ) -> Vec> - { - let mut product: Vec> = vec![vec![]; upper_tri_matrix.len()]; - let mut matrix_entry: ValueType; - - for i in 0..upper_tri_matrix.len() { - for (j, lower_tri_matrix_transpose_entry) in lower_tri_matrix_transpose.iter().enumerate() { - matrix_entry = ValueType::zero(); - let lower_diff: usize = (i as i64 - j as i64).max(0) as usize; - let upper_diff: usize = (j as i64 - i as i64).max(0) as usize; - for k in 0..(upper_tri_matrix[i].len().min(lower_tri_matrix_transpose_entry.len())) { - matrix_entry += upper_tri_matrix[i][k + upper_diff] * lower_tri_matrix_transpose_entry[k + lower_diff]; - } - product[i].push(matrix_entry); - } + for i in 1..matrix.len() { + prev_value = (rhs[i] - matrix[i][0] * prev_value) / matrix[i][1]; + result.push(prev_value) } - product + + result } // Compute the spline value at a given point. @@ -308,22 +190,7 @@ where .map(|x| (x[1] - x[0])) .collect(); - let (diag_matrix, lower_tri_matrix) = self.compute_lower_tri_and_diagonal_matrices(&self.time_steps); - - let mut inv_lower_tri_matrix: Vec> = self.invert_lower_tri_matrix(&lower_tri_matrix); - - let upper_tri_matrix: Vec> = self.lower_tri_transpose_inv_times_diag_inv( - &inv_lower_tri_matrix, - &diag_matrix, - ); - - inv_lower_tri_matrix = self.transpose(inv_lower_tri_matrix); - let tridiagonal_inverse: Vec> = self.upper_tri_times_lower_tri( - &upper_tri_matrix, - &inv_lower_tri_matrix - ); - - let rhs_vector: Vec = self.xs.iter().zip(self.ys.iter()) + let mut rhs: Vec = self.xs.iter().zip(self.ys.iter()) .collect::>() .windows(3) .map(|window| { @@ -345,17 +212,9 @@ where ) }).collect(); - self.second_derivatives = vec![ValueType::zero(); tridiagonal_inverse.len()]; - let mut matrix_entry: ValueType; - - for (i, tridiagonal_inverse_row) in tridiagonal_inverse.iter().enumerate() { - matrix_entry = self.second_derivatives[i]; - for (j, tridiagonal_inverse_entry) in tridiagonal_inverse_row.iter().enumerate() { - matrix_entry += *tridiagonal_inverse_entry * rhs_vector[j]; - } - self.second_derivatives[i] = matrix_entry; - } - + let lower_tri_matrix: Vec> = self.cholesky_decomposition(); + rhs = self.forward_substitution(&lower_tri_matrix, &rhs); + self.second_derivatives = self.backward_substitution(&lower_tri_matrix, &rhs); self.second_derivatives.insert(0, ValueType::zero()); self.second_derivatives.push(ValueType::zero()); From 7d604809e2dfc3f27f8801252f3ef1bc8c5e0a0f Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sat, 4 Oct 2025 19:32:55 +0100 Subject: [PATCH 12/19] Cubic-Spline: Amend test reference values to reflect more accurate interpolator implementation --- crates/RustQuant_math/src/interpolation/cubic_spline.rs | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index e8e2beb8..8009f8b4 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -292,7 +292,7 @@ mod tests_cubic_spline_interpolation { let _ = interpolator.fit(); assert_approx_equal!( - 36.64909523809524, + 36.660714285714285, interpolator.interpolate(2.5).unwrap(), RUSTQUANT_EPSILON ); @@ -316,7 +316,7 @@ mod tests_cubic_spline_interpolation { let _ = interpolator.fit(); assert_approx_equal!( - 36.64909523809524, + 36.660714285714285, interpolator .interpolate(xs[2] + time::Duration::hours(12)) .unwrap(), @@ -327,7 +327,7 @@ mod tests_cubic_spline_interpolation { #[test] fn test_cubic_interpolation_out_of_range() { let xs: Vec = vec![1., 2., 3., 4., 5.]; - let ys = vec![1., 2., 3., 4., 5.]; + let ys: Vec = vec![1., 2., 3., 4., 5.]; let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); let _ = interpolator.fit(); @@ -347,7 +347,7 @@ mod tests_cubic_spline_interpolation { let _ = interpolator.fit(); assert_approx_equal!( - 39.966361318634604, + 39.97368421052632, interpolator.interpolate(2.5).unwrap(), RUSTQUANT_EPSILON ); From c742b52c4158a88b8b10de8e56cdddd4250df7c5 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sat, 4 Oct 2025 19:34:26 +0100 Subject: [PATCH 13/19] Cubic-Spline: Rewrite helper methods tests for Cholesky decomposition implementation --- .../src/interpolation/cubic_spline.rs | 113 +++++------------- 1 file changed, 33 insertions(+), 80 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 8009f8b4..8fefcd8a 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -393,107 +393,60 @@ mod tests_cubic_spline_helper_functions { use super::*; #[test] - fn test_compute_lower_tri_and_diagonal_matrices() { + fn test_cholesky_decomposition() { let xs: Vec = vec![1., 2., 3., 4., 5.]; let ys: Vec = vec![1., 2., 3., 4., 5.]; - let time_steps: Vec = xs.windows(2) + let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); + interpolator.time_steps = interpolator.xs.windows(2) .map(|x| (x[1] - x[0])) .collect(); + let result: Vec> = interpolator.cholesky_decomposition(); + let expected: &[&[f64]] = &[ + &[2.0], + &[0.5, 1.9364916731037085], + &[0.5163977794943222, 1.9321835661585918], + ]; - let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); - let (lower_triangular, diagonal) = interpolator.compute_lower_tri_and_diagonal_matrices(&time_steps); - - assert!(lower_triangular == vec![4.0, 3.75, 3.7333333333333334]); - assert!(diagonal == vec![0.25, 0.26666666666666666]); - } - - #[test] - fn test_invert_lower_tri_matrix() { - let xs: Vec = vec![1., 2., 3., 4., 5.]; - let ys: Vec = vec![1., 2., 3., 4., 5.]; - - let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); - let inv_lower_tri_matrix: Vec> = interpolator.invert_lower_tri_matrix(&[3.0, 2.0, 1.0]); - - assert!( - inv_lower_tri_matrix == vec![ - vec![1.0], - vec![-3.0, 1.0], - vec![6.0, -2.0, 1.0], - vec![-6.0, 2.0, -1.0, 1.0] - ] - ) - } - - #[test] - fn test_transpose() { - let xs: Vec = vec![1., 2., 3., 4., 5.]; - let ys: Vec = vec![1., 2., 3., 4., 5.]; - - let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); - let transposed_matrix = interpolator.transpose(vec![ - vec![1.0], - vec![2.0, 3.0], - vec![4.0, 5.0, 6.0], - ]); - - assert!( - transposed_matrix == vec![ - vec![1.0, 2.0, 4.0], - vec![3.0, 5.0], - vec![6.0], - ] - ); + assert_eq!(result, expected); } #[test] - fn lower_tri_transpose_inv_times_diag_inv() { + fn test_forward_substitution() { let xs: Vec = vec![1., 2., 3., 4., 5.]; let ys: Vec = vec![1., 2., 3., 4., 5.]; let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); - let lower_tri_inverse: Vec> = vec![vec![1.0], vec![2.0, 3.0], vec![4.0, 5.0, 6.0]]; - let diagonal: Vec = vec![1.0, 2.0, 3.0]; - - assert!( - interpolator.lower_tri_transpose_inv_times_diag_inv(&lower_tri_inverse, &diagonal) == vec![ - vec![1.0, 2.0, 4.0], - vec![0.5, 2.5], - vec![0.3333333333333333] - ] + let result: Vec = interpolator.forward_substitution( + &[ + [2.0].to_vec(), + [1.0, 2.0].to_vec(), + [1.0, 2.0].to_vec(), + ], + &[1.0, 1.0, 1.0] ); + let expected: &[f64] = &[0.5, 0.25, 0.375]; + + assert_eq!(result, expected); } #[test] - fn test_upper_tri_times_lower_tri() { + fn test_backward_substitution() { let xs: Vec = vec![1., 2., 3., 4., 5.]; let ys: Vec = vec![1., 2., 3., 4., 5.]; let interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs, ys).unwrap(); - let upper_tri_matrix = vec![ - vec![1.0, 2.0, 3.0], - vec![4.0, 5.0], - vec![6.0], - ]; - let lower_tri_matrix_transpose = vec![ - vec![-1.0, -2.0, -3.0], - vec![-4.0, -5.0], - vec![-6.0], - ]; - - let product: Vec> = interpolator.upper_tri_times_lower_tri( - &upper_tri_matrix, - &lower_tri_matrix_transpose - ); - - assert!( - product == vec![ - [-14.0, -23.0, -18.0], - [-23.0, -41.0, -30.0], - [-18.0, -30.0, -36.0] - ] + let result: Vec = interpolator.backward_substitution( + &[ + [2.0].to_vec(), + [1.0, 2.0].to_vec(), + [1.0, 2.0].to_vec(), + ], + &[1.0, 1.0, 1.0] ); + let expected: &[f64] = &[0.375, 0.25, 0.5]; + + assert_eq!(result, expected); } #[test] @@ -513,6 +466,6 @@ mod tests_cubic_spline_helper_functions { 156.85714285714286 ); - assert!(spline_result == 36.64909523809524); + assert_eq!(spline_result, 36.64909523809524); } } From 74d31c3e246a941a1fec7e09de1a73430d3bea34 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sat, 4 Oct 2025 19:35:10 +0100 Subject: [PATCH 14/19] Cubic-Spline: Create square_root trait for the Cholesky decomposition --- .../RustQuant_math/src/interpolation/mod.rs | 50 ++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/crates/RustQuant_math/src/interpolation/mod.rs b/crates/RustQuant_math/src/interpolation/mod.rs index ff90115c..b0e0048f 100644 --- a/crates/RustQuant_math/src/interpolation/mod.rs +++ b/crates/RustQuant_math/src/interpolation/mod.rs @@ -7,7 +7,8 @@ // - LICENSE-MIT.md // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -use std::ops::{Div, Mul, MulAssign, Sub, Add, AddAssign, Neg}; +use std::{i8, ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub}}; +use rust_decimal::prelude::*; use RustQuant_error::RustQuantError; pub mod linear_interpolator; @@ -32,6 +33,7 @@ pub trait InterpolationValue: num::Num + Neg + AddAssign + MulAssign + + SquareRoot + Copy + Clone + Sized @@ -46,6 +48,12 @@ pub trait IntoValue { fn into_value(self) -> ValueType; } +/// Trait to compute the square root of a value. +pub trait SquareRoot { + /// Calculate square root of a value. + fn square_root(self) -> ValueType; +} + /// Trait describing requirements to be an index of interpolation. pub trait InterpolationIndex: Sub + PartialOrd + Copy + Clone + Sized + std::fmt::Display + Send + Sync @@ -95,6 +103,7 @@ impl InterpolationValue for T where T: num::Num + Neg + AddAssign + MulAssign + + SquareRoot + Copy + Clone + Sized @@ -132,25 +141,63 @@ macro_rules! impl_time_delta_into_value { }; } +macro_rules! impl_int_square_root { + ($c:ty) => { + impl SquareRoot<$c> for $c { + fn square_root(self) -> $c { + self.isqrt() + } + } + }; +} + +macro_rules! impl_float_square_root { + ($c:ty) => { + impl SquareRoot<$c> for $c { + fn square_root(self) -> $c { + self.sqrt() + } + } + }; +} + +macro_rules! impl_decimal_square_root { + ($c:ty) => { + impl SquareRoot<$c> for $c { + fn square_root(self) -> $c { + self.sqrt().unwrap() + } + } + }; +} + // Implement InterpolationIndex for all signed integer types. impl_interpolation_index!(i8, i8, i8); impl_num_delta_into_value!(i8, i8); +impl_int_square_root!(i8); impl_interpolation_index!(i16, i16, i16); impl_num_delta_into_value!(i16, i16); +impl_int_square_root!(i16); impl_interpolation_index!(i32, i32, i32); impl_num_delta_into_value!(i32, i32); +impl_int_square_root!(i32); impl_interpolation_index!(i64, i64, i64); impl_num_delta_into_value!(i64, i64); +impl_int_square_root!(i64); impl_interpolation_index!(i128, i128, i128); impl_num_delta_into_value!(i128, i128); +impl_int_square_root!(i128); impl_interpolation_index!(isize, isize, isize); impl_num_delta_into_value!(isize, isize); +impl_int_square_root!(isize); // Implement InterpolationIndex for all floating point types. impl_interpolation_index!(f32, f32, f32); impl_num_delta_into_value!(f32, f32); impl_interpolation_index!(f64, f64, f64); impl_num_delta_into_value!(f64, f64); +impl_float_square_root!(f32); +impl_float_square_root!(f64); // Implement InterpolationIndex for date/time types. impl_interpolation_index!(time::Date, time::Duration, f64); @@ -166,3 +213,4 @@ impl_interpolation_index!( rust_decimal::Decimal ); impl_num_delta_into_value!(rust_decimal::Decimal, rust_decimal::Decimal); +impl_decimal_square_root!(rust_decimal::Decimal); From fc5edfc5c1da2bcda16931884ed1078ef2152c8d Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sat, 4 Oct 2025 19:36:04 +0100 Subject: [PATCH 15/19] Cubic-Spline: Import maths feature for rust_decimal to apply square roots needed for Cholesky decomposition --- Cargo.toml | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 6cf5ed73..c80bfa57 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -97,7 +97,7 @@ plotly = "0.10.0" # https://docs.rs/plotly/latest/plotly/ rand = "0.8.5" # https://docs.rs/rand/latest/rand/ rand_distr = "0.4.3" # https://docs.rs/rand_distr/latest/rand_distr/ rayon = "1.9.0" # https://docs.rs/rayon/latest/rayon/ -rust_decimal = "1.34.3" # https://docs.rs/rust_decimal/latest/rust_decimal/ +rust_decimal = { version = "1.34.3", features = ["maths"] } # https://docs.rs/rust_decimal/latest/rust_decimal/ statrs = "0.17.1" # https://docs.rs/statrs/latest/statrs/ thiserror = "1.0.57" # https://docs.rs/thiserror/latest/thiserror/ @@ -117,4 +117,17 @@ polars = { version = "0.44.0", features = ["docs-selection"] } serde = { version = "1.0.213", features = ["derive"] } # https://docs.rs/crate/pyo3/latest -pyo3 = {version = "0.26.0", features = ["time"] } +pyo3 = {version = "0.26.0", features = ["extension-module", "time"]} + +## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +## PYTHON BINDINGS +## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +# [lib] +# name = "RustQuant" +# crate-type = ["cdylib"] + +# [dependencies.pyo3] +# version = "0.22.0" +# features = ["extension-module"] +# features = ["abi3-py37", "extension-module"] From 46211133978cac95520f9cd0ccbfe99b1f3654a3 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 5 Oct 2025 15:11:51 +0100 Subject: [PATCH 16/19] Replace dates with floats for interpolation and curves --- Cargo.toml | 1 + crates/RustQuant_data/Cargo.toml | 4 +- crates/RustQuant_data/src/curves.rs | 117 ++++++++++-------- .../src/interpolation/b_splines.rs | 26 ---- .../src/interpolation/cubic_spline.rs | 26 ---- .../interpolation/exponential_interpolator.rs | 20 --- .../src/interpolation/linear_interpolator.rs | 48 ------- .../RustQuant_math/src/interpolation/mod.rs | 17 --- .../RustQuant_stochastics/src/curve_model.rs | 6 +- .../src/nelson_siegel_svensson.rs | 23 ++-- 10 files changed, 79 insertions(+), 209 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index c80bfa57..0e4dbcf2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -100,6 +100,7 @@ rayon = "1.9.0" # https://docs.rs/rayon/latest/rayon/ rust_decimal = { version = "1.34.3", features = ["maths"] } # https://docs.rs/rust_decimal/latest/rust_decimal/ statrs = "0.17.1" # https://docs.rs/statrs/latest/statrs/ thiserror = "1.0.57" # https://docs.rs/thiserror/latest/thiserror/ +ordered-float = "5.1.0" # https://docs.rs/ordered-float/latest/ordered_float/ # https://docs.rs/ndarray/latest/ndarray/ ndarray = { version = "0.16.1", features = ["rayon"] } diff --git a/crates/RustQuant_data/Cargo.toml b/crates/RustQuant_data/Cargo.toml index 7554bda3..0a5feef6 100644 --- a/crates/RustQuant_data/Cargo.toml +++ b/crates/RustQuant_data/Cargo.toml @@ -30,8 +30,8 @@ time = { workspace = true } plotly = { workspace = true } argmin = { workspace = true } argmin-math = { workspace = true } -pyo3 = { workspace = true } - +ordered-float = { workspace=true } +pyo3 = { workspace=true } ## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ## RUSTDOC CONFIGURATION diff --git a/crates/RustQuant_data/src/curves.rs b/crates/RustQuant_data/src/curves.rs index bea485b1..7d07a77b 100644 --- a/crates/RustQuant_data/src/curves.rs +++ b/crates/RustQuant_data/src/curves.rs @@ -35,7 +35,7 @@ use plotly::{common::Mode, Plot, Scatter}; use pyo3::{pyclass, pymethods, PyResult}; use std::collections::BTreeMap; use std::sync::Arc; -use time::Date; +use ordered_float::OrderedFloat; use RustQuant_math::interpolation::{ExponentialInterpolator, Interpolator, LinearInterpolator}; use RustQuant_stochastics::{CurveModel, NelsonSiegelSvensson}; @@ -84,7 +84,7 @@ pub enum InterpolationMethod { #[pyclass] pub struct Curve { /// The nodes of the curve. - nodes: BTreeMap, + nodes: BTreeMap, f64>, /// The type of the curve. curve_type: CurveType, @@ -93,7 +93,7 @@ pub struct Curve { interpolation_method: InterpolationMethod, /// Interpolator backend. - interpolator: Arc>, + interpolator: Arc>, /// Nelson-Siegel-Svensson curve parameters. nss: Option, @@ -104,12 +104,12 @@ impl Curve { /// Create a new curve. #[new] pub fn new( - dates: Vec, + dates: Vec, rates: Vec, curve_type: CurveType, interpolation_method: InterpolationMethod, ) -> PyResult { - let interpolator: Arc> = match interpolation_method { + let interpolator: Arc> = match interpolation_method { InterpolationMethod::Linear => { Arc::new(LinearInterpolator::new(dates.clone(), rates.clone())?) } @@ -122,10 +122,10 @@ impl Curve { InterpolationMethod::Lagrange => { todo!("Implement LagrangeInterpolator") } - }; + }; Ok(Self { - nodes: dates.into_iter().zip(rates.into_iter()).collect(), + nodes: dates.into_iter().zip(rates.into_iter()).map(|(a, b)| (OrderedFloat(a), b)).collect(), curve_type, interpolation_method, interpolator, @@ -133,38 +133,39 @@ impl Curve { }) } - /// Create a new Curve from a list of nodes. - #[staticmethod] - pub fn from_nodes( - nodes: BTreeMap, - curve_type: CurveType, - interpolation_method: InterpolationMethod, - ) -> PyResult { - let interpolator: Arc> = match interpolation_method { - InterpolationMethod::Linear => Arc::new(LinearInterpolator::new( - nodes.keys().cloned().collect(), - nodes.values().cloned().collect(), - )?), - InterpolationMethod::Exponential => Arc::new(ExponentialInterpolator::new( - nodes.keys().cloned().collect(), - nodes.values().cloned().collect(), - )?), - InterpolationMethod::CubicSpline => { - todo!("Implement CubicSplineInterpolator") - } - InterpolationMethod::Lagrange => { - todo!("Implement LagrangeInterpolator") - } - }; - Ok(Self { - nodes, - curve_type, - interpolation_method, - interpolator, - nss: None, - }) - } + // /// Create a new Curve from a list of nodes. + // #[staticmethod] + // pub fn from_nodes( + // nodes: BTreeMap, + // curve_type: CurveType, + // interpolation_method: InterpolationMethod, + // ) -> PyResult { + // let interpolator: Arc> = match interpolation_method { + // InterpolationMethod::Linear => Arc::new(LinearInterpolator::new( + // nodes.keys().cloned().collect(), + // nodes.values().cloned().collect(), + // )?), + // InterpolationMethod::Exponential => Arc::new(ExponentialInterpolator::new( + // nodes.keys().cloned().collect(), + // nodes.values().cloned().collect(), + // )?), + // InterpolationMethod::CubicSpline => { + // todo!("Implement CubicSplineInterpolator") + // } + // InterpolationMethod::Lagrange => { + // todo!("Implement LagrangeInterpolator") + // } + // }; + + // Ok(Self { + // nodes, + // curve_type, + // interpolation_method, + // interpolator, + // nss: None, + // }) + // } /// Get the interpolation method used by the curve. pub fn interpolation_method(&self) -> InterpolationMethod { @@ -172,8 +173,8 @@ impl Curve { } /// Get a rate from the curve. - pub fn get_rate(&self, date: Date) -> Option { - match self.nodes.get(&date) { + pub fn get_rate(&self, date: f64) -> Option { + match self.nodes.get(&OrderedFloat(date)) { Some(rate) => Some(*rate), None => self.interpolator.interpolate(date).ok(), } @@ -192,7 +193,7 @@ impl Curve { } /// Get multiple rates from the curve. - pub fn get_rates(&self, dates: Vec) -> Vec> { + pub fn get_rates(&self, dates: Vec) -> Vec> { dates.iter().map(|date| self.get_rate(*date)).collect() } @@ -205,30 +206,36 @@ impl Curve { } /// Set a rate in the curve. - pub fn set_rate(&mut self, date: Date, rate: f64) { - self.nodes.insert(date, rate); + pub fn set_rate(&mut self, date: f64, rate: f64) { + self.nodes.insert(OrderedFloat(date), rate); } /// Set multiple rates in the curve. - pub fn set_rates(&mut self, rates: Vec<(Date, f64)>) { + pub fn set_rates(&mut self, rates: Vec<(f64, f64)>) { for (date, rate) in rates { self.set_rate(date, rate); } } /// Get the first date in the curve. - pub fn first_date(&self) -> Option<&Date> { - self.nodes.keys().next() + pub fn first_date(&self) -> Option<&f64> { + match self.nodes.keys().next() { + Some(date) => Some(&date.0), + None => None, + } } /// Get the last date in the curve. - pub fn last_date(&self) -> Option<&Date> { - self.nodes.keys().next_back() + pub fn last_date(&self) -> Option<&f64> { + match self.nodes.keys().next_back() { + Some(date) => Some(&date.0), + None => None, + } } /// Get the dates of the curve. - pub fn dates(&self) -> Vec { - self.nodes.keys().cloned().collect() + pub fn dates(&self) -> Vec { + self.nodes.keys().map(|k| k.0).collect()//.cloned().collect() } /// Get the rates of the curve. @@ -257,7 +264,7 @@ impl Curve { } /// Get the bracketing indices for a specific index. - pub fn get_brackets(&self, index: Date) -> (Date, Date) { + pub fn get_brackets(&self, index: f64) -> (f64, f64) { let first = self.first_date().unwrap(); let last = self.last_date().unwrap(); @@ -268,10 +275,10 @@ impl Curve { return (*last, *last); } - let left = self.nodes.range(..index).next_back().unwrap().0; - let right = self.nodes.range(index..).next().unwrap().0; + let left = self.nodes.range(..OrderedFloat(index)).next_back().unwrap().0; + let right = self.nodes.range(OrderedFloat(index)..).next().unwrap().0; - return (*left, *right); + return (left.0, right.0); } /// Shift the curve by a constant value. @@ -326,7 +333,7 @@ impl CostFunction for Curve { let y_model = x .into_iter() - .map(|date| curve_function(&nss, *date)) + .map(|date| curve_function(&nss, **date)) .collect::>(); let data = std::iter::zip(y, y_model); diff --git a/crates/RustQuant_math/src/interpolation/b_splines.rs b/crates/RustQuant_math/src/interpolation/b_splines.rs index 73813319..63f8c218 100644 --- a/crates/RustQuant_math/src/interpolation/b_splines.rs +++ b/crates/RustQuant_math/src/interpolation/b_splines.rs @@ -181,32 +181,6 @@ mod tests_b_splines { ); } - #[test] - fn test_b_spline_dates() { - let now = time::OffsetDateTime::now_utc(); - let knots: Vec = vec![ - now, - now + time::Duration::days(1), - now + time::Duration::days(2), - now + time::Duration::days(3), - now + time::Duration::days(4), - now + time::Duration::days(5), - now + time::Duration::days(6), - ]; - let control_points = vec![-1.0, 2.0, 0.0, -1.0]; - - let mut interpolator = BSplineInterpolator::new(knots.clone(), control_points, 2).unwrap(); - let _ = interpolator.fit(); - - assert_approx_equal!( - 1.375, - interpolator - .interpolate(knots[2] + time::Duration::hours(12)) - .unwrap(), - RUSTQUANT_EPSILON - ); - } - #[test] fn test_b_spline_inconsistent_parameters() { let knots = vec![0.0, 1.0, 2.0, 3.0, 4.0]; diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index 8fefcd8a..a7928982 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -298,32 +298,6 @@ mod tests_cubic_spline_interpolation { ); } - #[test] - fn test_natural_cubic_interpolation_dates() { - let now: time::OffsetDateTime = time::OffsetDateTime::now_utc(); - - let xs: Vec = vec![ - now, - now + time::Duration::days(1), - now + time::Duration::days(2), - now + time::Duration::days(3), - now + time::Duration::days(4), - ]; - - let ys: Vec = vec![0., 1., 16., 81., 256.]; - - let mut interpolator: CubicSplineInterpolator = CubicSplineInterpolator::new(xs.clone(), ys).unwrap(); - let _ = interpolator.fit(); - - assert_approx_equal!( - 36.660714285714285, - interpolator - .interpolate(xs[2] + time::Duration::hours(12)) - .unwrap(), - RUSTQUANT_EPSILON - ); - } - #[test] fn test_cubic_interpolation_out_of_range() { let xs: Vec = vec![1., 2., 3., 4., 5.]; diff --git a/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs b/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs index bdc90b6d..560ab8d3 100644 --- a/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs +++ b/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs @@ -188,24 +188,4 @@ mod tests_exponential_interpolation { RUSTQUANT_EPSILON ); } - - #[test] - fn test_exponential_interpolation_dates() { - let d_1m = date!(1990 - 06 - 16); - let d_2m = date!(1990 - 07 - 17); - - let r_1m = 0.9870; - let r_2m = 0.9753; - - let dates = vec![d_1m, d_2m]; - let rates = vec![r_1m, r_2m]; - - let interpolator = ExponentialInterpolator::new(dates, rates).unwrap(); - - assert_approx_equal!( - 0.9854824711068088, - interpolator.interpolate(date!(1990 - 06 - 20)).unwrap(), - RUSTQUANT_EPSILON - ); - } } diff --git a/crates/RustQuant_math/src/interpolation/linear_interpolator.rs b/crates/RustQuant_math/src/interpolation/linear_interpolator.rs index f98dec10..d0e9a01f 100644 --- a/crates/RustQuant_math/src/interpolation/linear_interpolator.rs +++ b/crates/RustQuant_math/src/interpolation/linear_interpolator.rs @@ -129,7 +129,6 @@ where #[cfg(test)] mod tests_linear_interpolation { use super::*; - use time::macros::date; use RustQuant_utils::{assert_approx_equal, RUSTQUANT_EPSILON}; #[test] @@ -162,51 +161,4 @@ mod tests_linear_interpolation { assert!(interpolator.interpolate(6.).is_err()); } - - #[test] - fn test_linear_interpolation_dates() { - let now = time::OffsetDateTime::now_utc(); - - let xs = vec![ - now, - now + time::Duration::days(1), - now + time::Duration::days(2), - now + time::Duration::days(3), - now + time::Duration::days(4), - ]; - - let ys = vec![1., 2., 3., 4., 5.]; - - let mut interpolator = LinearInterpolator::new(xs.clone(), ys).unwrap(); - let _ = interpolator.fit(); - - assert_approx_equal!( - 2.5, - interpolator - .interpolate(xs[1] + time::Duration::hours(12)) - .unwrap(), - RUSTQUANT_EPSILON - ); - } - - #[test] - fn test_linear_interpolation_dates_textbook() { - let d_1m = date!(1990 - 06 - 16); - let d_2m = date!(1990 - 07 - 17); - - let r_1m = 0.9870; - let r_2m = 0.9753; - - let dates = vec![d_1m, d_2m]; - let rates = vec![r_1m, r_2m]; - - let interpolator = LinearInterpolator::new(dates, rates).unwrap(); - - let d = date!(1990 - 06 - 20); - assert_approx_equal!( - interpolator.interpolate(d).unwrap(), - 0.9854903225806452, - RUSTQUANT_EPSILON - ); - } } diff --git a/crates/RustQuant_math/src/interpolation/mod.rs b/crates/RustQuant_math/src/interpolation/mod.rs index b0e0048f..9c1c8ea5 100644 --- a/crates/RustQuant_math/src/interpolation/mod.rs +++ b/crates/RustQuant_math/src/interpolation/mod.rs @@ -131,16 +131,6 @@ macro_rules! impl_num_delta_into_value { }; } -macro_rules! impl_time_delta_into_value { - ($b:ty, $c:ty) => { - impl IntoValue<$c> for $b { - fn into_value(self) -> $c { - self.as_seconds_f64() - } - } - }; -} - macro_rules! impl_int_square_root { ($c:ty) => { impl SquareRoot<$c> for $c { @@ -199,13 +189,6 @@ impl_num_delta_into_value!(f64, f64); impl_float_square_root!(f32); impl_float_square_root!(f64); -// Implement InterpolationIndex for date/time types. -impl_interpolation_index!(time::Date, time::Duration, f64); -impl_interpolation_index!(time::Time, time::Duration, f64); -impl_interpolation_index!(time::OffsetDateTime, time::Duration, f64); -impl_interpolation_index!(time::PrimitiveDateTime, time::Duration, f64); -impl_time_delta_into_value!(time::Duration, f64); - // Implement InterpolationIndex for Decimal type. impl_interpolation_index!( rust_decimal::Decimal, diff --git a/crates/RustQuant_stochastics/src/curve_model.rs b/crates/RustQuant_stochastics/src/curve_model.rs index fd3f41f2..c34c20b9 100644 --- a/crates/RustQuant_stochastics/src/curve_model.rs +++ b/crates/RustQuant_stochastics/src/curve_model.rs @@ -12,11 +12,11 @@ /// A trait for curve models. pub trait CurveModel { /// Returns the forward rate for a given date. - fn forward_rate(&self, date: time::Date) -> f64; + fn forward_rate(&self, date: f64) -> f64; /// Returns the spot rate for a given date. - fn spot_rate(&self, date: time::Date) -> f64; + fn spot_rate(&self, date: f64) -> f64; /// Returns the discount factor for a given date. - fn discount_factor(&self, date: time::Date) -> f64; + fn discount_factor(&self, date: f64) -> f64; } diff --git a/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs b/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs index a3b71210..df5e0219 100644 --- a/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs +++ b/crates/RustQuant_stochastics/src/nelson_siegel_svensson.rs @@ -65,10 +65,10 @@ impl NelsonSiegelSvensson { impl CurveModel for NelsonSiegelSvensson { /// Returns the forward rate for a given date. - fn forward_rate(&self, date: Date) -> f64 { - assert!(date > today(), "Date must be in the future."); + fn forward_rate(&self, tau: f64) -> f64 { + // assert!(date > today(), "Date must be in the future."); - let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); + // let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); let term1 = f64::exp(-tau / self.lambda1); let term2 = (tau / self.lambda1) * term1; @@ -78,10 +78,10 @@ impl CurveModel for NelsonSiegelSvensson { } /// Returns the spot rate for a given date. - fn spot_rate(&self, date: Date) -> f64 { - assert!(date > today(), "Date must be in the future."); + fn spot_rate(&self, tau: f64) -> f64 { + // assert!(date > today(), "Date must be in the future."); - let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); + // let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); let term1 = self.lambda1 * (1. - f64::exp(-tau / self.lambda1)) / tau; let term2 = term1 - f64::exp(-tau / self.lambda1); @@ -91,10 +91,10 @@ impl CurveModel for NelsonSiegelSvensson { self.beta0 + self.beta1 * term1 + self.beta2 * term2 + self.beta3 * term3 } - fn discount_factor(&self, date: Date) -> f64 { - let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); + fn discount_factor(&self, tau: f64) -> f64 { + // let tau = DayCountConvention::Actual_365_25.day_count_factor(today(), date); - f64::exp(-self.spot_rate(date) * tau / 100.) + f64::exp(-self.spot_rate(tau) * tau / 100.) } } @@ -105,7 +105,6 @@ impl CurveModel for NelsonSiegelSvensson { #[cfg(test)] mod tests_nelson_siegel_svensson { use super::*; - use time::Duration; #[test] fn test_nelson_siegel_svensson() { @@ -119,8 +118,8 @@ mod tests_nelson_siegel_svensson { }; let dates = (2..365 * 30) - .map(|i| today() + Duration::days(i)) - .collect::>(); + .map(|i| 1.0 / (i as f64)) + .collect::>(); let _forward_curve = dates .iter() From a326478416cd6ec174f61c5047d1458017c9c95a Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 5 Oct 2025 15:52:52 +0100 Subject: [PATCH 17/19] Cubic-Spline: Clippy amendments --- crates/RustQuant_math/src/interpolation/cubic_spline.rs | 2 +- crates/RustQuant_math/src/interpolation/mod.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/RustQuant_math/src/interpolation/cubic_spline.rs b/crates/RustQuant_math/src/interpolation/cubic_spline.rs index a7928982..da0aea9a 100644 --- a/crates/RustQuant_math/src/interpolation/cubic_spline.rs +++ b/crates/RustQuant_math/src/interpolation/cubic_spline.rs @@ -187,7 +187,7 @@ where fn fit(&mut self) -> Result<(), RustQuantError> { self.time_steps = self.xs.windows(2) - .map(|x| (x[1] - x[0])) + .map(|x| x[1] - x[0]) .collect(); let mut rhs: Vec = self.xs.iter().zip(self.ys.iter()) diff --git a/crates/RustQuant_math/src/interpolation/mod.rs b/crates/RustQuant_math/src/interpolation/mod.rs index 9c1c8ea5..2fa2bc62 100644 --- a/crates/RustQuant_math/src/interpolation/mod.rs +++ b/crates/RustQuant_math/src/interpolation/mod.rs @@ -7,7 +7,7 @@ // - LICENSE-MIT.md // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -use std::{i8, ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub}}; +use std::ops::{Add, AddAssign, Div, Mul, MulAssign, Neg, Sub}; use rust_decimal::prelude::*; use RustQuant_error::RustQuantError; From d49ffd250c8f040cda89820487ad436e4d6b0eea Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 5 Oct 2025 18:06:45 +0100 Subject: [PATCH 18/19] Missed replacement of Date with f64 --- crates/RustQuant_data/src/curves.rs | 8 ++++---- .../src/interpolation/exponential_interpolator.rs | 1 - 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/crates/RustQuant_data/src/curves.rs b/crates/RustQuant_data/src/curves.rs index 7d07a77b..adc29b01 100644 --- a/crates/RustQuant_data/src/curves.rs +++ b/crates/RustQuant_data/src/curves.rs @@ -181,12 +181,12 @@ impl Curve { } /// Get a rate, and simultaneously add it to the nodes. - pub fn get_rate_and_insert(&mut self, date: Date) -> Option { - match self.nodes.get(&date) { + pub fn get_rate_and_insert(&mut self, date: f64) -> Option { + match self.nodes.get(&OrderedFloat(date)) { Some(rate) => Some(*rate), None => { let rate = self.interpolator.interpolate(date).ok()?; - self.nodes.insert(date, rate); + self.nodes.insert(OrderedFloat(date), rate); Some(rate) } } @@ -198,7 +198,7 @@ impl Curve { } /// Get multiple rates from the curve, and simultaneously add them to the nodes. - pub fn get_rates_and_insert(&mut self, dates: Vec) -> Vec> { + pub fn get_rates_and_insert(&mut self, dates: Vec) -> Vec> { dates .iter() .map(|date| self.get_rate_and_insert(*date)) diff --git a/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs b/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs index 560ab8d3..58d075c9 100644 --- a/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs +++ b/crates/RustQuant_math/src/interpolation/exponential_interpolator.rs @@ -173,7 +173,6 @@ where #[cfg(test)] mod tests_exponential_interpolation { use super::*; - use time::macros::date; use RustQuant_utils::{assert_approx_equal, RUSTQUANT_EPSILON}; #[test] From c61fe8e0d45a846f40be16b561095d1a47e7a393 Mon Sep 17 00:00:00 2001 From: Yasser Naji Date: Sun, 5 Oct 2025 23:54:21 +0100 Subject: [PATCH 19/19] Cubic-Spline: Further Clippy amendments --- crates/RustQuant_data/src/curves.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/crates/RustQuant_data/src/curves.rs b/crates/RustQuant_data/src/curves.rs index adc29b01..0c2b7e0d 100644 --- a/crates/RustQuant_data/src/curves.rs +++ b/crates/RustQuant_data/src/curves.rs @@ -125,7 +125,7 @@ impl Curve { }; Ok(Self { - nodes: dates.into_iter().zip(rates.into_iter()).map(|(a, b)| (OrderedFloat(a), b)).collect(), + nodes: dates.into_iter().zip(rates).map(|(a, b)| (OrderedFloat(a), b)).collect(), curve_type, interpolation_method, interpolator, @@ -278,7 +278,7 @@ impl Curve { let left = self.nodes.range(..OrderedFloat(index)).next_back().unwrap().0; let right = self.nodes.range(OrderedFloat(index)..).next().unwrap().0; - return (left.0, right.0); + (left.0, right.0) } /// Shift the curve by a constant value.