diff --git a/src/lib.rs b/src/lib.rs index 4ac8fb3..083b95f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -88,5 +88,7 @@ pub mod model_selection; /// Supervised neighbors-based learning methods pub mod neighbors; pub(crate) mod optimization; +/// Support Vector Machines +pub mod svm; /// Supervised tree-based learning methods pub mod tree; diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 3937b47..a637bdf 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -85,6 +85,12 @@ pub trait BaseVector: Clone + Debug { /// Create new vector of size `len` where each element is set to `value`. fn fill(len: usize, value: T) -> Self; + + /// Vector dot product + fn dot(&self, other: &Self) -> T; + + /// Returns True if matrices are element-wise equal within a tolerance `error`. + fn approximate_eq(&self, other: &Self, error: T) -> bool; } /// Generic matrix type. @@ -110,6 +116,10 @@ pub trait BaseMatrix: Clone + Debug { /// * `row` - row number fn get_row_as_vec(&self, row: usize) -> Vec; + /// Get the `row`'th row + /// * `row` - row number + fn get_row(&self, row: usize) -> Self::RowVector; + /// Copies a vector with elements of the `row`'th row into `result` /// * `row` - row number /// * `result` - receiver for the row diff --git a/src/linalg/naive/dense_matrix.rs b/src/linalg/naive/dense_matrix.rs index ae9d1d2..ab6bf43 100644 --- a/src/linalg/naive/dense_matrix.rs +++ b/src/linalg/naive/dense_matrix.rs @@ -44,6 +44,32 @@ impl BaseVector for Vec { fn fill(len: usize, value: T) -> Self { vec![value; len] } + + fn dot(&self, other: &Self) -> T { + if self.len() != other.len() { + panic!("A and B should have the same size"); + } + + let mut result = T::zero(); + for i in 0..self.len() { + result = result + self[i] * other[i]; + } + + result + } + + fn approximate_eq(&self, other: &Self, error: T) -> bool { + if self.len() != other.len() { + false + } else { + for i in 0..other.len() { + if (self[i] - other[i]).abs() > error { + return false; + } + } + true + } + } } /// Column-major, dense matrix. See [Simple Dense Matrix](../index.html). @@ -371,6 +397,16 @@ impl BaseMatrix for DenseMatrix { self.values[col * self.nrows + row] } + fn get_row(&self, row: usize) -> Self::RowVector { + let mut v = vec![T::zero(); self.ncols]; + + for c in 0..self.ncols { + v[c] = self.get(row, c); + } + + v + } + fn get_row_as_vec(&self, row: usize) -> Vec { let mut result = vec![T::zero(); self.ncols]; for c in 0..self.ncols { @@ -865,6 +901,21 @@ impl BaseMatrix for DenseMatrix { mod tests { use super::*; + #[test] + fn vec_dot() { + let v1 = vec![1., 2., 3.]; + let v2 = vec![4., 5., 6.]; + assert_eq!(32.0, BaseVector::dot(&v1, &v2)); + } + + #[test] + fn vec_approximate_eq() { + let a = vec![1., 2., 3.]; + let b = vec![1. + 1e-5, 2. + 2e-5, 3. + 3e-5]; + assert!(a.approximate_eq(&b, 1e-4)); + assert!(!a.approximate_eq(&b, 1e-5)); + } + #[test] fn from_array() { let vec = [1., 2., 3., 4., 5., 6.]; @@ -939,6 +990,12 @@ mod tests { assert_eq!(result, expected); } + #[test] + fn get_row() { + let a = DenseMatrix::from_2d_array(&[&[1., 2., 3.], &[4., 5., 6.], &[7., 8., 9.]]); + assert_eq!(vec![4., 5., 6.], a.get_row(1)); + } + #[test] fn matmul() { let a = DenseMatrix::from_2d_array(&[&[1., 2., 3.], &[4., 5., 6.]]); diff --git a/src/linalg/nalgebra_bindings.rs b/src/linalg/nalgebra_bindings.rs index 5e52d14..badd8c4 100644 --- a/src/linalg/nalgebra_bindings.rs +++ b/src/linalg/nalgebra_bindings.rs @@ -79,6 +79,20 @@ impl BaseVector for MatrixMN { m.fill(value); m } + + fn dot(&self, other: &Self) -> T { + self.dot(other) + } + + fn approximate_eq(&self, other: &Self, error: T) -> bool { + if self.shape() != other.shape() { + false + } else { + self.iter() + .zip(other.iter()) + .all(|(a, b)| (*a - *b).abs() <= error) + } + } } impl @@ -102,6 +116,10 @@ impl Self::RowVector { + self.row(row).into_owned() + } + fn copy_row_as_vec(&self, row: usize, result: &mut Vec) { let mut r = 0; for e in self.row(row).iter() { @@ -486,6 +504,21 @@ mod tests { assert_eq!(twos, RowDVector::from_vec(vec![2., 2., 2.])); } + #[test] + fn vec_dot() { + let v1 = RowDVector::from_vec(vec![1., 2., 3.]); + let v2 = RowDVector::from_vec(vec![4., 5., 6.]); + assert_eq!(32.0, BaseVector::dot(&v1, &v2)); + } + + #[test] + fn vec_approximate_eq() { + let a = RowDVector::from_vec(vec![1., 2., 3.]); + let noise = RowDVector::from_vec(vec![1e-5, 2e-5, 3e-5]); + assert!(a.approximate_eq(&(&noise + &a), 1e-4)); + assert!(!a.approximate_eq(&(&noise + &a), 1e-5)); + } + #[test] fn get_set_dynamic() { let mut m = DMatrix::from_row_slice(2, 3, &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0]); @@ -579,6 +612,12 @@ mod tests { assert_eq!(m.get_col_as_vec(1), vec!(2., 5., 8.)); } + #[test] + fn get_row() { + let a = DMatrix::from_row_slice(3, 3, &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]); + assert_eq!(RowDVector::from_vec(vec![4., 5., 6.]), a.get_row(1)); + } + #[test] fn copy_row_col_as_vec() { let m = DMatrix::from_row_slice(3, 3, &[1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0]); diff --git a/src/linalg/ndarray_bindings.rs b/src/linalg/ndarray_bindings.rs index c6324eb..9cfb6d7 100644 --- a/src/linalg/ndarray_bindings.rs +++ b/src/linalg/ndarray_bindings.rs @@ -57,7 +57,7 @@ use crate::linalg::Matrix; use crate::linalg::{BaseMatrix, BaseVector}; use crate::math::num::RealNumber; -impl BaseVector for ArrayBase, Ix1> { +impl BaseVector for ArrayBase, Ix1> { fn get(&self, i: usize) -> T { self[i] } @@ -84,6 +84,14 @@ impl BaseVector for ArrayBase, Ix1> { fn fill(len: usize, value: T) -> Self { Array::from_elem(len, value) } + + fn dot(&self, other: &Self) -> T { + self.dot(other) + } + + fn approximate_eq(&self, other: &Self, error: T) -> bool { + (self - other).iter().all(|v| v.abs() <= error) + } } impl @@ -109,6 +117,10 @@ impl Self::RowVector { + self.row(row).to_owned() + } + fn copy_row_as_vec(&self, row: usize, result: &mut Vec) { let mut r = 0; for e in self.row(row).iter() { @@ -437,6 +449,21 @@ mod tests { assert_eq!(vec![1., 2., 3.], v.to_vec()); } + #[test] + fn vec_dot() { + let v1 = arr1(&[1., 2., 3.]); + let v2 = arr1(&[4., 5., 6.]); + assert_eq!(32.0, BaseVector::dot(&v1, &v2)); + } + + #[test] + fn vec_approximate_eq() { + let a = arr1(&[1., 2., 3.]); + let noise = arr1(&[1e-5, 2e-5, 3e-5]); + assert!(a.approximate_eq(&(&noise + &a), 1e-4)); + assert!(!a.approximate_eq(&(&noise + &a), 1e-5)); + } + #[test] fn from_to_row_vec() { let vec = arr1(&[1., 2., 3.]); @@ -678,6 +705,12 @@ mod tests { assert_eq!(res, vec![4., 5., 6.]); } + #[test] + fn get_row() { + let a = arr2(&[[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]]); + assert_eq!(arr1(&[4., 5., 6.]), a.get_row(1)); + } + #[test] fn get_col_as_vector() { let a = arr2(&[[1., 2., 3.], [4., 5., 6.], [7., 8., 9.]]); diff --git a/src/math/num.rs b/src/math/num.rs index 9ffbdff..894e5a3 100644 --- a/src/math/num.rs +++ b/src/math/num.rs @@ -6,10 +6,23 @@ use num_traits::{Float, FromPrimitive}; use rand::prelude::*; use std::fmt::{Debug, Display}; use std::iter::{Product, Sum}; +use std::ops::{AddAssign, DivAssign, MulAssign, SubAssign}; /// Defines real number /// -pub trait RealNumber: Float + FromPrimitive + Debug + Display + Copy + Sum + Product { +pub trait RealNumber: + Float + + FromPrimitive + + Debug + + Display + + Copy + + Sum + + Product + + AddAssign + + SubAssign + + MulAssign + + DivAssign +{ /// Copy sign from `sign` - another real number fn copysign(self, sign: Self) -> Self; diff --git a/src/svm/mod.rs b/src/svm/mod.rs new file mode 100644 index 0000000..404b281 --- /dev/null +++ b/src/svm/mod.rs @@ -0,0 +1,25 @@ +//! # Support Vector Machines +//! + +pub mod svr; + +use serde::{Deserialize, Serialize}; + +use crate::linalg::BaseVector; +use crate::math::num::RealNumber; + +/// Kernel +pub trait Kernel> { + /// Apply kernel function to x_i and x_j + fn apply(&self, x_i: &V, x_j: &V) -> T; +} + +/// Linear Kernel +#[derive(Serialize, Deserialize, Debug)] +pub struct LinearKernel {} + +impl> Kernel for LinearKernel { + fn apply(&self, x_i: &V, x_j: &V) -> T { + x_i.dot(x_j) + } +} diff --git a/src/svm/svr.rs b/src/svm/svr.rs new file mode 100644 index 0000000..371c4ca --- /dev/null +++ b/src/svm/svr.rs @@ -0,0 +1,538 @@ +//! # Epsilon-Support Vector Regression. +//! +//! Example +//! +//! ``` +//! use smartcore::linalg::naive::dense_matrix::*; +//! use smartcore::linear::linear_regression::*; +//! use smartcore::svm::*; +//! use smartcore::svm::svr::{SVR, SVRParameters}; +//! +//! // Longley dataset (https://www.statsmodels.org/stable/datasets/generated/longley.html) +//! let x = DenseMatrix::from_2d_array(&[ +//! &[234.289, 235.6, 159.0, 107.608, 1947., 60.323], +//! &[259.426, 232.5, 145.6, 108.632, 1948., 61.122], +//! &[258.054, 368.2, 161.6, 109.773, 1949., 60.171], +//! &[284.599, 335.1, 165.0, 110.929, 1950., 61.187], +//! &[328.975, 209.9, 309.9, 112.075, 1951., 63.221], +//! &[346.999, 193.2, 359.4, 113.270, 1952., 63.639], +//! &[365.385, 187.0, 354.7, 115.094, 1953., 64.989], +//! &[363.112, 357.8, 335.0, 116.219, 1954., 63.761], +//! &[397.469, 290.4, 304.8, 117.388, 1955., 66.019], +//! &[419.180, 282.2, 285.7, 118.734, 1956., 67.857], +//! &[442.769, 293.6, 279.8, 120.445, 1957., 68.169], +//! &[444.546, 468.1, 263.7, 121.950, 1958., 66.513], +//! &[482.704, 381.3, 255.2, 123.366, 1959., 68.655], +//! &[502.601, 393.1, 251.4, 125.368, 1960., 69.564], +//! &[518.173, 480.6, 257.2, 127.852, 1961., 69.331], +//! &[554.894, 400.7, 282.7, 130.081, 1962., 70.551], +//! ]); +//! +//! let y: Vec = vec![83.0, 88.5, 88.2, 89.5, 96.2, 98.1, 99.0, +//! 100.0, 101.2, 104.6, 108.4, 110.8, 112.6, 114.2, 115.7, 116.9]; +//! +//! let svr = SVR::fit(&x, &y, +//! LinearKernel {}, +//! SVRParameters { +//! eps: 2.0, +//! c: 10.0, +//! tol: 1e-3, +//! }).unwrap(); +//! +//! let y_hat = svr.predict(&x).unwrap(); +//! ``` +use std::cell::{Ref, RefCell}; +use std::fmt::Debug; + +use serde::{Deserialize, Serialize}; + +use crate::error::Failed; +use crate::linalg::BaseVector; +use crate::linalg::Matrix; +use crate::math::num::RealNumber; +use crate::svm::Kernel; + +#[derive(Serialize, Deserialize, Debug)] + +/// SVR Parameters +pub struct SVRParameters { + /// Epsilon in the epsilon-SVR model + pub eps: T, + /// Regularization parameter. + pub c: T, + /// Tolerance for stopping criterion + pub tol: T, +} + +#[derive(Serialize, Deserialize, Debug)] +#[serde(bound( + serialize = "M::RowVector: Serialize, K: Serialize, T: Serialize", + deserialize = "M::RowVector: Deserialize<'de>, K: Deserialize<'de>, T: Deserialize<'de>", +))] + +/// Epsilon-Support Vector Regression +pub struct SVR, K: Kernel> { + kernel: K, + instances: Vec, + w: Vec, + b: T, +} + +#[derive(Serialize, Deserialize, Debug)] +struct SupportVector> { + index: usize, + x: V, + alpha: [T; 2], + grad: [T; 2], + k: T, +} + +struct Optimizer<'a, T: RealNumber, M: Matrix, K: Kernel> { + tol: T, + c: T, + svmin: usize, + svmax: usize, + gmin: T, + gmax: T, + gminindex: usize, + gmaxindex: usize, + tau: T, + sv: Vec>, + kernel: &'a K, +} + +struct Cache { + data: Vec>>>, +} + +impl Default for SVRParameters { + fn default() -> Self { + SVRParameters { + eps: T::from_f64(0.1).unwrap(), + c: T::one(), + tol: T::from_f64(1e-3).unwrap(), + } + } +} + +impl, K: Kernel> SVR { + /// Fits SVR to your data. + /// * `x` - _NxM_ matrix with _N_ observations and _M_ features in each observation. + /// * `y` - target values + /// * `kernel` - the kernel function + /// * `parameters` - optional parameters, use `Default::default()` to set parameters to default values. + pub fn fit( + x: &M, + y: &M::RowVector, + kernel: K, + parameters: SVRParameters, + ) -> Result, Failed> { + let (n, _) = x.shape(); + + if n != y.len() { + return Err(Failed::fit(&format!( + "Number of rows of X doesn't match number of rows of Y" + ))); + } + + let optimizer = Optimizer::optimize(x, y, &kernel, ¶meters); + + let (support_vectors, weight, b) = optimizer.smo(); + + Ok(SVR { + kernel: kernel, + instances: support_vectors, + w: weight, + b: b, + }) + } + + /// Predict target values from `x` + /// * `x` - _KxM_ data where _K_ is number of observations and _M_ is number of features. + pub fn predict(&self, x: &M) -> Result { + let (n, _) = x.shape(); + + let mut y_hat = M::RowVector::zeros(n); + + for i in 0..n { + y_hat.set(i, self.predict_for_row(x.get_row(i))); + } + + Ok(y_hat) + } + + pub(in crate) fn predict_for_row(&self, x: M::RowVector) -> T { + let mut f = self.b; + + for i in 0..self.instances.len() { + f += self.w[i] * self.kernel.apply(&x, &self.instances[i]); + } + + return f; + } +} + +impl, K: Kernel> PartialEq for SVR { + fn eq(&self, other: &Self) -> bool { + if self.b != other.b + || self.w.len() != other.w.len() + || self.instances.len() != other.instances.len() + { + return false; + } else { + for i in 0..self.w.len() { + if (self.w[i] - other.w[i]).abs() > T::epsilon() { + return false; + } + } + for i in 0..self.instances.len() { + if !self.instances[i].approximate_eq(&other.instances[i], T::epsilon()) { + return false; + } + } + return true; + } + } +} + +impl> SupportVector { + fn new>(i: usize, x: V, y: T, eps: T, k: &K) -> SupportVector { + let k_v = k.apply(&x, &x); + SupportVector { + index: i, + x: x, + grad: [eps + y, eps - y], + k: k_v, + alpha: [T::zero(), T::zero()], + } + } +} + +impl<'a, T: RealNumber, M: Matrix, K: Kernel> Optimizer<'a, T, M, K> { + fn optimize( + x: &M, + y: &M::RowVector, + kernel: &'a K, + parameters: &SVRParameters, + ) -> Optimizer<'a, T, M, K> { + let (n, _) = x.shape(); + + let mut support_vectors: Vec> = Vec::with_capacity(n); + + for i in 0..n { + support_vectors.push(SupportVector::new( + i, + x.get_row(i), + y.get(i), + parameters.eps, + kernel, + )); + } + + Optimizer { + tol: parameters.tol, + c: parameters.c, + svmin: 0, + svmax: 0, + gmin: T::max_value(), + gmax: T::min_value(), + gminindex: 0, + gmaxindex: 0, + tau: T::from_f64(1e-12).unwrap(), + sv: support_vectors, + kernel: kernel, + } + } + + fn minmax(&mut self) { + self.gmin = T::max_value(); + self.gmax = T::min_value(); + + for i in 0..self.sv.len() { + let v = &self.sv[i]; + let g = -v.grad[0]; + let a = v.alpha[0]; + if g < self.gmin && a > T::zero() { + self.gmin = g; + self.gminindex = 0; + self.svmin = i; + } + if g > self.gmax && a < self.c { + self.gmax = g; + self.gmaxindex = 0; + self.svmax = i; + } + + let g = v.grad[1]; + let a = v.alpha[1]; + if g < self.gmin && a < self.c { + self.gmin = g; + self.gminindex = 1; + self.svmin = i; + } + if g > self.gmax && a > T::zero() { + self.gmax = g; + self.gmaxindex = 1; + self.svmax = i; + } + } + } + + fn smo(mut self) -> (Vec, Vec, T) { + let cache: Cache = Cache::new(self.sv.len()); + + self.minmax(); + + while self.gmax - self.gmin > self.tol { + let v1 = self.svmax; + let i = self.gmaxindex; + let old_alpha_i = self.sv[v1].alpha[i]; + + let k1 = cache.get(self.sv[v1].index, || { + self.sv + .iter() + .map(|vi| self.kernel.apply(&self.sv[v1].x, &vi.x)) + .collect() + }); + + let mut v2 = self.svmin; + let mut j = self.gminindex; + let mut old_alpha_j = self.sv[v2].alpha[j]; + + let mut best = T::zero(); + let gi = if i == 0 { + -self.sv[v1].grad[0] + } else { + self.sv[v1].grad[1] + }; + for jj in 0..self.sv.len() { + let v = &self.sv[jj]; + let mut curv = self.sv[v1].k + v.k - T::two() * k1[v.index]; + if curv <= T::zero() { + curv = self.tau; + } + + let mut gj = -v.grad[0]; + if v.alpha[0] > T::zero() && gj < gi { + let gain = -((gi - gj) * (gi - gj)) / curv; + if gain < best { + best = gain; + v2 = jj; + j = 0; + old_alpha_j = self.sv[v2].alpha[0]; + } + } + + gj = v.grad[1]; + if v.alpha[1] < self.c && gj < gi { + let gain = -((gi - gj) * (gi - gj)) / curv; + if gain < best { + best = gain; + v2 = jj; + j = 1; + old_alpha_j = self.sv[v2].alpha[1]; + } + } + } + + let k2 = cache.get(self.sv[v2].index, || { + self.sv + .iter() + .map(|vi| self.kernel.apply(&self.sv[v2].x, &vi.x)) + .collect() + }); + + let mut curv = self.sv[v1].k + self.sv[v2].k - T::two() * k1[self.sv[v2].index]; + if curv <= T::zero() { + curv = self.tau; + } + + if i != j { + let delta = (-self.sv[v1].grad[i] - self.sv[v2].grad[j]) / curv; + let diff = self.sv[v1].alpha[i] - self.sv[v2].alpha[j]; + self.sv[v1].alpha[i] += delta; + self.sv[v2].alpha[j] += delta; + + if diff > T::zero() { + if self.sv[v2].alpha[j] < T::zero() { + self.sv[v2].alpha[j] = T::zero(); + self.sv[v1].alpha[i] = diff; + } + } else { + if self.sv[v1].alpha[i] < T::zero() { + self.sv[v1].alpha[i] = T::zero(); + self.sv[v2].alpha[j] = -diff; + } + } + + if diff > T::zero() { + if self.sv[v1].alpha[i] > self.c { + self.sv[v1].alpha[i] = self.c; + self.sv[v2].alpha[j] = self.c - diff; + } + } else { + if self.sv[v2].alpha[j] > self.c { + self.sv[v2].alpha[j] = self.c; + self.sv[v1].alpha[i] = self.c + diff; + } + } + } else { + let delta = (self.sv[v1].grad[i] - self.sv[v2].grad[j]) / curv; + let sum = self.sv[v1].alpha[i] + self.sv[v2].alpha[j]; + self.sv[v1].alpha[i] -= delta; + self.sv[v2].alpha[j] += delta; + + if sum > self.c { + if self.sv[v1].alpha[i] > self.c { + self.sv[v1].alpha[i] = self.c; + self.sv[v2].alpha[j] = sum - self.c; + } + } else { + if self.sv[v2].alpha[j] < T::zero() { + self.sv[v2].alpha[j] = T::zero(); + self.sv[v1].alpha[i] = sum; + } + } + + if sum > self.c { + if self.sv[v2].alpha[j] > self.c { + self.sv[v2].alpha[j] = self.c; + self.sv[v1].alpha[i] = sum - self.c; + } + } else { + if self.sv[v1].alpha[i] < T::zero() { + self.sv[v1].alpha[i] = T::zero(); + self.sv[v2].alpha[j] = sum; + } + } + } + + let delta_alpha_i = self.sv[v1].alpha[i] - old_alpha_i; + let delta_alpha_j = self.sv[v2].alpha[j] - old_alpha_j; + + let si = T::two() * T::from_usize(i).unwrap() - T::one(); + let sj = T::two() * T::from_usize(j).unwrap() - T::one(); + for v in self.sv.iter_mut() { + v.grad[0] -= si * k1[v.index] * delta_alpha_i + sj * k2[v.index] * delta_alpha_j; + v.grad[1] += si * k1[v.index] * delta_alpha_i + sj * k2[v.index] * delta_alpha_j; + } + + self.minmax(); + } + + let b = -(self.gmax + self.gmin) / T::two(); + + let mut result: Vec = Vec::new(); + let mut alpha: Vec = Vec::new(); + + for v in self.sv { + if v.alpha[0] != v.alpha[1] { + result.push(v.x); + alpha.push(v.alpha[1] - v.alpha[0]); + } + } + + (result, alpha, b) + } +} + +impl Cache { + fn new(n: usize) -> Cache { + Cache { + data: vec![RefCell::new(None); n], + } + } + + fn get Vec>(&self, i: usize, or: F) -> Ref> { + if self.data[i].borrow().is_none() { + self.data[i].replace(Some(or())); + } + Ref::map(self.data[i].borrow(), |v| v.as_ref().unwrap()) + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::linalg::naive::dense_matrix::*; + use crate::metrics::mean_squared_error; + use crate::svm::*; + + #[test] + fn svr_fit_predict() { + let x = DenseMatrix::from_2d_array(&[ + &[234.289, 235.6, 159.0, 107.608, 1947., 60.323], + &[259.426, 232.5, 145.6, 108.632, 1948., 61.122], + &[258.054, 368.2, 161.6, 109.773, 1949., 60.171], + &[284.599, 335.1, 165.0, 110.929, 1950., 61.187], + &[328.975, 209.9, 309.9, 112.075, 1951., 63.221], + &[346.999, 193.2, 359.4, 113.270, 1952., 63.639], + &[365.385, 187.0, 354.7, 115.094, 1953., 64.989], + &[363.112, 357.8, 335.0, 116.219, 1954., 63.761], + &[397.469, 290.4, 304.8, 117.388, 1955., 66.019], + &[419.180, 282.2, 285.7, 118.734, 1956., 67.857], + &[442.769, 293.6, 279.8, 120.445, 1957., 68.169], + &[444.546, 468.1, 263.7, 121.950, 1958., 66.513], + &[482.704, 381.3, 255.2, 123.366, 1959., 68.655], + &[502.601, 393.1, 251.4, 125.368, 1960., 69.564], + &[518.173, 480.6, 257.2, 127.852, 1961., 69.331], + &[554.894, 400.7, 282.7, 130.081, 1962., 70.551], + ]); + + let y: Vec = vec![ + 83.0, 88.5, 88.2, 89.5, 96.2, 98.1, 99.0, 100.0, 101.2, 104.6, 108.4, 110.8, 112.6, + 114.2, 115.7, 116.9, + ]; + + let y_hat = SVR::fit( + &x, + &y, + LinearKernel {}, + SVRParameters { + eps: 2.0, + c: 10.0, + tol: 1e-3, + }, + ) + .and_then(|lr| lr.predict(&x)) + .unwrap(); + + println!("{:?}", y_hat); + + assert!(mean_squared_error(&y_hat, &y) < 2.5); + } + + #[test] + fn svr_serde() { + let x = DenseMatrix::from_2d_array(&[ + &[234.289, 235.6, 159.0, 107.608, 1947., 60.323], + &[259.426, 232.5, 145.6, 108.632, 1948., 61.122], + &[258.054, 368.2, 161.6, 109.773, 1949., 60.171], + &[284.599, 335.1, 165.0, 110.929, 1950., 61.187], + &[328.975, 209.9, 309.9, 112.075, 1951., 63.221], + &[346.999, 193.2, 359.4, 113.270, 1952., 63.639], + &[365.385, 187.0, 354.7, 115.094, 1953., 64.989], + &[363.112, 357.8, 335.0, 116.219, 1954., 63.761], + &[397.469, 290.4, 304.8, 117.388, 1955., 66.019], + &[419.180, 282.2, 285.7, 118.734, 1956., 67.857], + &[442.769, 293.6, 279.8, 120.445, 1957., 68.169], + &[444.546, 468.1, 263.7, 121.950, 1958., 66.513], + &[482.704, 381.3, 255.2, 123.366, 1959., 68.655], + &[502.601, 393.1, 251.4, 125.368, 1960., 69.564], + &[518.173, 480.6, 257.2, 127.852, 1961., 69.331], + &[554.894, 400.7, 282.7, 130.081, 1962., 70.551], + ]); + + let y: Vec = vec![ + 83.0, 88.5, 88.2, 89.5, 96.2, 98.1, 99.0, 100.0, 101.2, 104.6, 108.4, 110.8, 112.6, + 114.2, 115.7, 116.9, + ]; + + let svr = SVR::fit(&x, &y, LinearKernel {}, Default::default()).unwrap(); + + let deserialized_svr: SVR, LinearKernel> = + serde_json::from_str(&serde_json::to_string(&svr).unwrap()).unwrap(); + + assert_eq!(svr, deserialized_svr); + } +}