Added wiener kernel

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<f64> {
     let mut r = Array2::zeros((ts1.len(), ts2.len()));
@@ -30,6 +34,7 @@ pub trait Kernel {
     fn state_cov(&self, t: f64) -> Array2<f64>;
     fn measurement_vector(&self) -> Array1<f64>;
     fn feedback(&self) -> Array2<f64>;
+    fn transition(&self, t0: f64, t1: f64) -> Array2<f64>;
 
     fn noise_effect(&self) -> Array2<f64> {
         unimplemented!();
@@ -39,10 +44,8 @@ pub trait Kernel {
         unimplemented!();
     }
 
-    fn transition(&self, t0: f64, t1: f64) -> Array2<f64>;
-
     fn noise_cov(&self, _t0: f64, _t1: f64) -> Array2<f64> {
-        todo!();
+        unimplemented!();
     }
 }
 

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<f64> {
+        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<f64> {
+        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<f64> {
+        array![0.0, 0.0]
+    }
+
+    fn state_cov(&self, t: f64) -> Array2<f64> {
+        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<f64> {
+        array![1.0, 0.0]
+    }
+
+    fn feedback(&self) -> Array2<f64> {
+        array![[0.0, 1.0], [0.0, 0.0]]
+    }
+
+    fn noise_effect(&self) -> Array2<f64> {
+        array![[0.0], [1.0]]
+    }
+
+    fn transition(&self, t0: f64, t1: f64) -> Array2<f64> {
+        array![[1.0, t1 - t0], [0.0, 1.0]]
+    }
+
+    fn noise_cov(&self, _t0: f64, _t1: f64) -> Array2<f64> {
+        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::<f64, _>(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::<f64, _>(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::<Vec<_>>();
+
+        assert_abs_diff_eq!(Array::from(vars), kernel.k_diag(&ts));
+    }
+}

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<f64> {
+        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<f64> {
+        (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<f64> {
+        array![0.0]
+    }
+
+    fn state_cov(&self, t: f64) -> Array2<f64> {
+        array![[self.var * (t - self.t0) + self.var_t0]]
+    }
+
+    fn measurement_vector(&self) -> Array1<f64> {
+        array![1.0]
+    }
+
+    fn feedback(&self) -> Array2<f64> {
+        array![[0.0]]
+    }
+
+    fn noise_effect(&self) -> Array2<f64> {
+        array![[1.0]]
+    }
+
+    fn transition(&self, _t0: f64, _t1: f64) -> Array2<f64> {
+        array![[1.0]]
+    }
+
+    fn noise_cov(&self, t0: f64, t1: f64) -> Array2<f64> {
+        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::<f64, _>(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::<f64, _>(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::<Vec<_>>();
+
+        assert_abs_diff_eq!(Array::from(vars), kernel.k_diag(&ts));
+    }
+}