From 0345052a66fd271602c5b3b9458dd8bfd8332577 Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Mon, 4 Jan 2021 13:38:55 +0100 Subject: [PATCH] Added wiener kernel --- src/kernel.rs | 9 ++- src/kernel/affine.rs | 152 +++++++++++++++++++++++++++++++++++++++++++ src/kernel/wiener.rs | 143 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 301 insertions(+), 3 deletions(-) create mode 100644 src/kernel/affine.rs create mode 100644 src/kernel/wiener.rs diff --git a/src/kernel.rs b/src/kernel.rs index 34fe9b9..20571b0 100644 --- a/src/kernel.rs +++ b/src/kernel.rs @@ -1,14 +1,18 @@ use ndarray::prelude::*; +mod affine; mod constant; mod exponential; mod matern32; mod matern52; +mod wiener; +pub use affine::Affine; pub use constant::Constant; pub use exponential::Exponential; pub use matern32::Matern32; pub use matern52::Matern52; +pub use wiener::Wiener; pub(crate) fn distance(ts1: &[f64], ts2: &[f64]) -> Array2 { let mut r = Array2::zeros((ts1.len(), ts2.len())); @@ -30,6 +34,7 @@ pub trait Kernel { fn state_cov(&self, t: f64) -> Array2; fn measurement_vector(&self) -> Array1; fn feedback(&self) -> Array2; + fn transition(&self, t0: f64, t1: f64) -> Array2; fn noise_effect(&self) -> Array2 { unimplemented!(); @@ -39,10 +44,8 @@ pub trait Kernel { unimplemented!(); } - fn transition(&self, t0: f64, t1: f64) -> Array2; - fn noise_cov(&self, _t0: f64, _t1: f64) -> Array2 { - todo!(); + unimplemented!(); } } diff --git a/src/kernel/affine.rs b/src/kernel/affine.rs new file mode 100644 index 0000000..c904932 --- /dev/null +++ b/src/kernel/affine.rs @@ -0,0 +1,152 @@ +use std::iter::FromIterator; + +use ndarray::prelude::*; + +use super::Kernel; + +pub struct Affine { + var_offset: f64, + var_slope: f64, + t0: f64, +} + +impl Affine { + pub fn new(var_offset: f64, var_slope: f64, t0: f64) -> Self { + Affine { + var_offset, + var_slope, + t0, + } + } +} + +impl Kernel for Affine { + fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> Array2 { + let ts2 = ts2.unwrap_or(ts1); + + let mut r = Array2::zeros((ts1.len(), ts2.len())); + + for (i, v1) in ts1.iter().enumerate() { + for (j, v2) in ts2.iter().enumerate() { + r[(i, j)] = (v1 - self.t0) * (v2 - self.t0); + } + } + + (r * self.var_slope) + self.var_offset + } + + fn k_diag(&self, ts: &[f64]) -> Array1 { + let r = Array1::from_iter(ts.iter().map(|v| (v - self.t0).powi(2))); + + (r * self.var_slope) + self.var_offset + } + + fn order(&self) -> usize { + 2 + } + + fn state_mean(&self, _t: f64) -> Array1 { + array![0.0, 0.0] + } + + fn state_cov(&self, t: f64) -> Array2 { + let t = t - self.t0; + + array![[t.powi(2), t], [t, 1.0]] * self.var_slope + + array![[1.0, 0.0], [0.0, 0.0]] * self.var_offset + } + + fn measurement_vector(&self) -> Array1 { + array![1.0, 0.0] + } + + fn feedback(&self) -> Array2 { + array![[0.0, 1.0], [0.0, 0.0]] + } + + fn noise_effect(&self) -> Array2 { + array![[0.0], [1.0]] + } + + fn transition(&self, t0: f64, t1: f64) -> Array2 { + array![[1.0, t1 - t0], [0.0, 1.0]] + } + + fn noise_cov(&self, _t0: f64, _t1: f64) -> Array2 { + array![[0.0, 0.0], [0.0, 0.0]] + } +} + +#[cfg(test)] +mod tests { + extern crate intel_mkl_src; + + use approx::assert_abs_diff_eq; + use rand::{distributions::Standard, thread_rng, Rng}; + + use super::*; + + #[test] + fn test_kernel_matrix() { + let kernel = Affine::new(0.0, 2.0, 0.0); + + let ts = [1.26, 1.46, 2.67]; + + assert_abs_diff_eq!( + kernel.k_mat(&ts, None), + array![ + [3.1752000000000002, 3.6792, 6.7284], + [3.6792, 4.263199999999999, 7.796399999999999], + [6.7284, 7.796399999999999, 14.2578] + ] + ); + } + + #[test] + fn test_kernel_diag() { + let kernel = Affine::new(0.0, 2.0, 0.0); + + let ts: Vec<_> = thread_rng() + .sample_iter::(Standard) + .take(10) + .map(|x| x * 10.0) + .collect(); + + assert_eq!(kernel.k_mat(&ts, None).diag(), kernel.k_diag(&ts)); + } + + #[test] + fn test_kernel_order() { + let kernel = Affine::new(0.0, 2.0, 0.0); + + let m = kernel.order(); + + assert_eq!(kernel.state_mean(0.0).shape(), &[m]); + assert_eq!(kernel.state_cov(0.0).shape(), &[m, m]); + assert_eq!(kernel.measurement_vector().shape(), &[m]); + assert_eq!(kernel.feedback().shape(), &[m, m]); + assert_eq!(kernel.noise_effect().shape()[0], m); + assert_eq!(kernel.transition(0.0, 1.0).shape(), &[m, m]); + assert_eq!(kernel.noise_cov(0.0, 1.0).shape(), &[m, m]); + } + + #[test] + fn test_ssm_variance() { + let kernel = Affine::new(0.0, 2.0, 0.0); + + let ts: Vec<_> = thread_rng() + .sample_iter::(Standard) + .take(10) + .map(|x| x * 10.0) + .collect(); + + let h = kernel.measurement_vector(); + + let vars = ts + .iter() + .map(|t| h.dot(&kernel.state_cov(*t)).dot(&h)) + .collect::>(); + + assert_abs_diff_eq!(Array::from(vars), kernel.k_diag(&ts)); + } +} diff --git a/src/kernel/wiener.rs b/src/kernel/wiener.rs new file mode 100644 index 0000000..ada3721 --- /dev/null +++ b/src/kernel/wiener.rs @@ -0,0 +1,143 @@ +use std::iter::FromIterator; + +use ndarray::prelude::*; + +use super::Kernel; + +pub struct Wiener { + var: f64, + t0: f64, + var_t0: f64, +} + +impl Wiener { + pub fn new(var: f64, t0: f64, var_t0: f64) -> Self { + Wiener { var, t0, var_t0 } + } +} + +impl Kernel for Wiener { + fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> Array2 { + let ts2 = ts2.unwrap_or(ts1); + + let mut r = Array2::zeros((ts1.len(), ts2.len())); + + for (i, v1) in ts1.iter().enumerate() { + for (j, v2) in ts2.iter().enumerate() { + r[(i, j)] = if v1 < v2 { v1 - self.t0 } else { v2 - self.t0 }; + } + } + + (r * self.var) + self.var_t0 + } + + fn k_diag(&self, ts: &[f64]) -> Array1 { + (Array1::from_iter(ts.iter().map(|v| v - self.t0)) * self.var) + self.var_t0 + } + + fn order(&self) -> usize { + 1 + } + + fn state_mean(&self, _t: f64) -> Array1 { + array![0.0] + } + + fn state_cov(&self, t: f64) -> Array2 { + array![[self.var * (t - self.t0) + self.var_t0]] + } + + fn measurement_vector(&self) -> Array1 { + array![1.0] + } + + fn feedback(&self) -> Array2 { + array![[0.0]] + } + + fn noise_effect(&self) -> Array2 { + array![[1.0]] + } + + fn transition(&self, _t0: f64, _t1: f64) -> Array2 { + array![[1.0]] + } + + fn noise_cov(&self, t0: f64, t1: f64) -> Array2 { + array![[self.var * (t1 - t0)]] + } +} + +#[cfg(test)] +mod tests { + extern crate intel_mkl_src; + + use approx::assert_abs_diff_eq; + use rand::{distributions::Standard, thread_rng, Rng}; + + use super::*; + + #[test] + fn test_kernel_matrix() { + let kernel = Wiener::new(1.2, 0.0, 0.0); + + let ts = [1.26, 1.46, 2.67]; + + assert_abs_diff_eq!( + kernel.k_mat(&ts, None), + array![ + [1.512, 1.512, 1.512], + [1.512, 1.752, 1.752], + [1.512, 1.752, 3.2039999999999997] + ] + ); + } + + #[test] + fn test_kernel_diag() { + let kernel = Wiener::new(1.2, 0.0, 0.0); + + let ts: Vec<_> = thread_rng() + .sample_iter::(Standard) + .take(10) + .map(|x| x * 10.0) + .collect(); + + assert_eq!(kernel.k_mat(&ts, None).diag(), kernel.k_diag(&ts)); + } + + #[test] + fn test_kernel_order() { + let kernel = Wiener::new(1.2, 0.0, 0.0); + + let m = kernel.order(); + + assert_eq!(kernel.state_mean(0.0).shape(), &[m]); + assert_eq!(kernel.state_cov(0.0).shape(), &[m, m]); + assert_eq!(kernel.measurement_vector().shape(), &[m]); + assert_eq!(kernel.feedback().shape(), &[m, m]); + assert_eq!(kernel.noise_effect().shape()[0], m); + assert_eq!(kernel.transition(0.0, 1.0).shape(), &[m, m]); + assert_eq!(kernel.noise_cov(0.0, 1.0).shape(), &[m, m]); + } + + #[test] + fn test_ssm_variance() { + let kernel = Wiener::new(1.2, 0.0, 0.0); + + let ts: Vec<_> = thread_rng() + .sample_iter::(Standard) + .take(10) + .map(|x| x * 10.0) + .collect(); + + let h = kernel.measurement_vector(); + + let vars = ts + .iter() + .map(|t| h.dot(&kernel.state_cov(*t)).dot(&h)) + .collect::>(); + + assert_abs_diff_eq!(Array::from(vars), kernel.k_diag(&ts)); + } +}