More tests

src/kernel.rs

@@ -10,8 +10,35 @@ pub use exponential::Exponential;
 pub use matern32::Matern32;
 pub use matern52::Matern52;
 
+#[inline]
+pub(crate) fn transition(t0: f64, t1: f64, feedback: Array2<f64>) -> Array2<f64> {
+    let a = feedback * (t1 - t0);
+
+    if a.shape() == [1, 1] {
+        array![[a[(0, 0)].exp()]]
+    } else {
+        let mut b = Array2::<f64>::zeros(a.dim());
+
+        crate::expm::expm(&a, &mut b);
+
+        b
+    }
+}
+
+pub(crate) fn distance(ts1: &[f64], ts2: &[f64]) -> Array2<f64> {
+    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 - v2).abs();
+        }
+    }
+
+    r
+}
+
 pub trait Kernel {
-    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> ArrayD<f64>;
+    fn k_mat(&self, ts1: &[f64], ts2: Option<&[f64]>) -> Array2<f64>;
     fn k_diag(&self, ts: &[f64]) -> Array1<f64>;
     fn order(&self) -> usize;
     fn state_mean(&self, t: f64) -> Array1<f64>;
@@ -28,14 +55,7 @@ pub trait Kernel {
     }
 
     fn transition(&self, t0: f64, t1: f64) -> Array2<f64> {
-        let f = self.feedback();
-
-        let a = f * (t1 - t0);
-        let mut b = Array2::<f64>::zeros(a.dim());
-
-        crate::expm::expm(&a, &mut b);
-
-        b
+        transition(t0, t1, self.feedback())
     }
 
     fn noise_cov(&self, _t0: f64, _t1: f64) -> Array2<f64> {
@@ -60,7 +80,7 @@ pub trait Kernel {
 }
 
 impl Kernel for Vec<Box<dyn Kernel>> {
-    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> ArrayD<f64> {
+    fn k_mat(&self, _ts1: &[f64], _ts2: Option<&[f64]>) -> Array2<f64> {
         unimplemented!();
     }