From e20e9ca6e01b1b076e9386e2ac9d0e771db0e8bc Mon Sep 17 00:00:00 2001 From: Volodymyr Orlov Date: Fri, 5 Jun 2020 10:40:17 -0700 Subject: [PATCH] feat: adds new distance measures + LU decomposition --- benches/distance.rs | 4 +- src/algorithm/neighbour/cover_tree.rs | 14 +- src/algorithm/neighbour/linear_search.rs | 4 +- src/linalg/evd.rs | 4 +- src/linalg/lu.rs | 254 ++++++++++++++++++ src/linalg/mod.rs | 8 +- src/linalg/naive/dense_matrix.rs | 38 +++ src/linalg/nalgebra_bindings.rs | 7 + src/linalg/ndarray_bindings.rs | 7 + src/math/distance/euclidian.rs | 12 +- src/math/distance/hamming.rs | 45 ++++ src/math/distance/mahalanobis.rs | 97 +++++++ src/math/distance/manhattan.rs | 43 +++ src/math/distance/minkowski.rs | 63 +++++ src/math/distance/mod.rs | 18 +- .../first_order/gradient_descent.rs | 4 +- 16 files changed, 594 insertions(+), 28 deletions(-) create mode 100644 src/linalg/lu.rs create mode 100644 src/math/distance/hamming.rs create mode 100644 src/math/distance/mahalanobis.rs create mode 100644 src/math/distance/manhattan.rs create mode 100644 src/math/distance/minkowski.rs diff --git a/benches/distance.rs b/benches/distance.rs index d9020f3..b9512e8 100644 --- a/benches/distance.rs +++ b/benches/distance.rs @@ -4,12 +4,12 @@ extern crate smartcore; use criterion::Criterion; use criterion::black_box; -use smartcore::math::distance::euclidian::*; +use smartcore::math::distance::*; fn criterion_benchmark(c: &mut Criterion) { let a = vec![1., 2., 3.]; - c.bench_function("Euclidean Distance", move |b| b.iter(|| Euclidian::distance(black_box(&a), black_box(&a)))); + c.bench_function("Euclidean Distance", move |b| b.iter(|| Distances::euclidian().distance(black_box(&a), black_box(&a)))); } criterion_group!(benches, criterion_benchmark); diff --git a/src/algorithm/neighbour/cover_tree.rs b/src/algorithm/neighbour/cover_tree.rs index dd73d18..2df810e 100644 --- a/src/algorithm/neighbour/cover_tree.rs +++ b/src/algorithm/neighbour/cover_tree.rs @@ -44,7 +44,7 @@ impl> CoverTree } else { let mut parent: Option = Option::None; let mut p_i = 0; - let mut qi_p_ds = vec!((self.root(), D::distance(&p, &self.root().data))); + let mut qi_p_ds = vec!((self.root(), self.distance.distance(&p, &self.root().data))); let mut i = self.max_level; loop { let i_d = self.base.powf(F::from(i).unwrap()); @@ -84,7 +84,7 @@ impl> CoverTree } pub fn find(&self, p: &T, k: usize) -> Vec{ - let mut qi_p_ds = vec!((self.root(), D::distance(&p, &self.root().data))); + let mut qi_p_ds = vec!((self.root(), self.distance.distance(&p, &self.root().data))); for i in (self.min_level..self.max_level+1).rev() { let i_d = self.base.powf(F::from(i).unwrap()); let mut q_p_ds = self.get_children_dist(&p, &qi_p_ds, i); @@ -115,7 +115,7 @@ impl> CoverTree let p = &self.nodes.get(p_id.index).unwrap().data; let mut i = 0; while i != s.len() { - let d = D::distance(p, &s[i]); + let d = self.distance.distance(p, &s[i]); if d <= r { my_near.0.push(s.remove(i)); } else if d > r && d <= F::two() * r{ @@ -169,7 +169,7 @@ impl> CoverTree let q: Vec<&Node> = qi_p_ds.iter().flat_map(|(n, _)| self.get_child(n, i)).collect(); - children.extend(q.into_iter().map(|n| (n, D::distance(&n.data, &p)))); + children.extend(q.into_iter().map(|n| (n, self.distance.distance(&n.data, &p)))); children @@ -219,7 +219,7 @@ impl> CoverTree let mut p_selected: Vec<&Node> = Vec::new(); for p in next_nodes { for q in nodes { - if D::distance(&p.data, &q.data) <= tree.base.powf(F::from(i).unwrap()) { + if tree.distance.distance(&p.data, &q.data) <= tree.base.powf(F::from(i).unwrap()) { p_selected.push(*p); } } @@ -233,7 +233,7 @@ impl> CoverTree for p in nodes { for q in nodes { if p != q { - assert!(D::distance(&p.data, &q.data) > tree.base.powf(F::from(i).unwrap())); + assert!(tree.distance.distance(&p.data, &q.data) > tree.base.powf(F::from(i).unwrap())); } } } @@ -280,7 +280,7 @@ mod tests { struct SimpleDistance{} impl Distance for SimpleDistance { - fn distance(a: &i32, b: &i32) -> f64 { + fn distance(&self, a: &i32, b: &i32) -> f64 { (a - b).abs() as f64 } } diff --git a/src/algorithm/neighbour/linear_search.rs b/src/algorithm/neighbour/linear_search.rs index aea64aa..bca3d60 100644 --- a/src/algorithm/neighbour/linear_search.rs +++ b/src/algorithm/neighbour/linear_search.rs @@ -38,7 +38,7 @@ impl> LinearKNNSearch { for i in 0..self.data.len() { - let d = D::distance(&from, &self.data[i]); + let d = self.distance.distance(&from, &self.data[i]); let datum = heap.peek_mut(); if d < datum.distance { datum.distance = d; @@ -81,7 +81,7 @@ mod tests { struct SimpleDistance{} impl Distance for SimpleDistance { - fn distance(a: &i32, b: &i32) -> f64 { + fn distance(&self, a: &i32, b: &i32) -> f64 { (a - b).abs() as f64 } } diff --git a/src/linalg/evd.rs b/src/linalg/evd.rs index a74e186..1f6c772 100644 --- a/src/linalg/evd.rs +++ b/src/linalg/evd.rs @@ -829,9 +829,7 @@ mod tests { &[0.6952105, 0.43984484, -0.7036135] ]); - let evd = A.evd(false); - - println!("{}", &evd.V.abs()); + let evd = A.evd(false); assert!(eigen_vectors.abs().approximate_eq(&evd.V.abs(), 1e-4)); for i in 0..eigen_values.len() { diff --git a/src/linalg/lu.rs b/src/linalg/lu.rs new file mode 100644 index 0000000..b2f1919 --- /dev/null +++ b/src/linalg/lu.rs @@ -0,0 +1,254 @@ +#![allow(non_snake_case)] + +use std::fmt::Debug; +use std::marker::PhantomData; + +use crate::math::num::FloatExt; +use crate::linalg::BaseMatrix; + +#[derive(Debug, Clone)] +pub struct LU> { + LU: M, + pivot: Vec, + pivot_sign: i8, + singular: bool, + phantom: PhantomData +} + +impl> LU { + pub fn new(LU: M, pivot: Vec, pivot_sign: i8) -> LU { + + let (_, n) = LU.shape(); + + let mut singular = false; + for j in 0..n { + if LU.get(j, j) == T::zero() { + singular = true; + break; + } + } + + LU { + LU: LU, + pivot: pivot, + pivot_sign: pivot_sign, + singular: singular, + phantom: PhantomData + } + } + + pub fn L(&self) -> M { + let (n_rows, n_cols) = self.LU.shape(); + let mut L = M::zeros(n_rows, n_cols); + + for i in 0..n_rows { + for j in 0..n_cols { + if i > j { + L.set(i, j, self.LU.get(i, j)); + } else if i == j { + L.set(i, j, T::one()); + } else { + L.set(i, j, T::zero()); + } + } + } + + L + } + + pub fn U(&self) -> M { + let (n_rows, n_cols) = self.LU.shape(); + let mut U = M::zeros(n_rows, n_cols); + + for i in 0..n_rows { + for j in 0..n_cols { + if i <= j { + U.set(i, j, self.LU.get(i, j)); + } else { + U.set(i, j, T::zero()); + } + } + } + + U + } + + pub fn pivot(&self) -> M { + let (_, n) = self.LU.shape(); + let mut piv = M::zeros(n, n); + + for i in 0..n { + piv.set(i, self.pivot[i], T::one()); + } + + piv + } + + pub fn inverse(&self) -> M { + let (m, n) = self.LU.shape(); + + if m != n { + panic!("Matrix is not square: {}x{}", m, n); + } + + let mut inv = M::zeros(n, n); + + for i in 0..n { + inv.set(i, i, T::one()); + } + + inv = self.solve(inv); + return inv; + } + + fn solve(&self, mut b: M) -> M { + let (m, n) = self.LU.shape(); + let (b_m, b_n) = b.shape(); + + if b_m != m { + panic!("Row dimensions do not agree: A is {} x {}, but B is {} x {}", m, n, b_m, b_n); + } + + if self.singular { + panic!("Matrix is singular."); + } + + let mut X = M::zeros(b_m, b_n); + + for j in 0..b_n { + for i in 0..m { + X.set(i, j, b.get(self.pivot[i], j)); + } + } + + for k in 0..n { + for i in k+1..n { + for j in 0..b_n { + X.sub_element_mut(i, j, X.get(k, j) * self.LU.get(i, k)); + } + } + } + + for k in (0..n).rev() { + for j in 0..b_n { + X.div_element_mut(k, j, self.LU.get(k, k)); + } + + for i in 0..k { + for j in 0..b_n { + X.sub_element_mut(i, j, X.get(k, j) * self.LU.get(i, k)); + } + } + } + + for j in 0..b_n { + for i in 0..m { + b.set(i, j, X.get(i, j)); + } + } + + b + + } + +} + +pub trait LUDecomposableMatrix: BaseMatrix { + + fn lu(&self) -> LU { + self.clone().lu_mut() + } + + fn lu_mut(mut self) -> LU { + + let (m, n) = self.shape(); + + let mut piv = vec![0; m]; + for i in 0..m { + piv[i] = i; + } + + let mut pivsign = 1; + let mut LUcolj = vec![T::zero(); m]; + + for j in 0..n { + + for i in 0..m { + LUcolj[i] = self.get(i, j); + } + + for i in 0..m { + let kmax = usize::min(i, j); + let mut s = T::zero(); + for k in 0..kmax { + s = s + self.get(i, k) * LUcolj[k]; + } + + LUcolj[i] = LUcolj[i] - s; + self.set(i, j, LUcolj[i]); + } + + let mut p = j; + for i in j+1..m { + if LUcolj[i].abs() > LUcolj[p].abs() { + p = i; + } + } + if p != j { + for k in 0..n { + let t = self.get(p, k); + self.set(p, k, self.get(j, k)); + self.set(j, k, t); + } + let k = piv[p]; + piv[p] = piv[j]; + piv[j] = k; + pivsign = -pivsign; + } + + if j < m && self.get(j, j) != T::zero() { + for i in j+1..m { + self.div_element_mut(i, j, self.get(j, j)); + } + } + } + + LU::new(self, piv, pivsign) + + } + + fn lu_solve_mut(self, b: Self) -> Self { + + self.lu_mut().solve(b) + + } +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::linalg::naive::dense_matrix::*; + + #[test] + fn decompose() { + + let a = DenseMatrix::from_array(&[&[1., 2., 3.], &[0., 1., 5.], &[5., 6., 0.]]); + let expected_L = DenseMatrix::from_array(&[&[1. , 0. , 0. ], &[0. , 1. , 0. ], &[0.2, 0.8, 1. ]]); + let expected_U = DenseMatrix::from_array(&[&[ 5., 6., 0.], &[ 0., 1., 5.], &[ 0., 0., -1.]]); + let expected_pivot = DenseMatrix::from_array(&[&[0., 0., 1.], &[0., 1., 0.], &[1., 0., 0.]]); + let lu = a.lu(); + assert!(lu.L().approximate_eq(&expected_L, 1e-4)); + assert!(lu.U().approximate_eq(&expected_U, 1e-4)); + assert!(lu.pivot().approximate_eq(&expected_pivot, 1e-4)); + } + + #[test] + fn inverse() { + + let a = DenseMatrix::from_array(&[&[1., 2., 3.], &[0., 1., 5.], &[5., 6., 0.]]); + let expected = DenseMatrix::from_array(&[&[-6.0, 3.6, 1.4], &[5.0, -3.0, -1.0], &[-1.0, 0.8, 0.2]]); + let a_inv = a.lu().inverse(); + println!("{}", a_inv); + assert!(a_inv.approximate_eq(&expected, 1e-4)); + } +} \ No newline at end of file diff --git a/src/linalg/mod.rs b/src/linalg/mod.rs index 9a64a0f..dff103a 100644 --- a/src/linalg/mod.rs +++ b/src/linalg/mod.rs @@ -2,6 +2,7 @@ pub mod naive; pub mod qr; pub mod svd; pub mod evd; +pub mod lu; pub mod ndarray_bindings; pub mod nalgebra_bindings; @@ -13,6 +14,7 @@ use crate::math::num::FloatExt; use svd::SVDDecomposableMatrix; use evd::EVDDecomposableMatrix; use qr::QRDecomposableMatrix; +use lu::LUDecomposableMatrix; pub trait BaseMatrix: Clone + Debug { @@ -172,11 +174,13 @@ pub trait BaseMatrix: Clone + Debug { fn argmax(&self) -> Vec; - fn unique(&self) -> Vec; + fn unique(&self) -> Vec; + + fn cov(&self) -> Self; } -pub trait Matrix: BaseMatrix + SVDDecomposableMatrix + EVDDecomposableMatrix + QRDecomposableMatrix + PartialEq + Display {} +pub trait Matrix: BaseMatrix + SVDDecomposableMatrix + EVDDecomposableMatrix + QRDecomposableMatrix + LUDecomposableMatrix + PartialEq + Display {} pub fn row_iter>(m: &M) -> RowIter { RowIter{ diff --git a/src/linalg/naive/dense_matrix.rs b/src/linalg/naive/dense_matrix.rs index dd38361..062db5f 100644 --- a/src/linalg/naive/dense_matrix.rs +++ b/src/linalg/naive/dense_matrix.rs @@ -13,6 +13,7 @@ pub use crate::linalg::BaseMatrix; use crate::linalg::svd::SVDDecomposableMatrix; use crate::linalg::evd::EVDDecomposableMatrix; use crate::linalg::qr::QRDecomposableMatrix; +use crate::linalg::lu::LUDecomposableMatrix; use crate::math::num::FloatExt; #[derive(Debug, Clone)] @@ -188,6 +189,8 @@ impl EVDDecomposableMatrix for DenseMatrix {} impl QRDecomposableMatrix for DenseMatrix {} +impl LUDecomposableMatrix for DenseMatrix {} + impl Matrix for DenseMatrix {} impl PartialEq for DenseMatrix { @@ -679,6 +682,34 @@ impl BaseMatrix for DenseMatrix { result } + fn cov(&self) -> Self { + + let (m, n) = self.shape(); + + let mu = self.column_mean(); + + let mut cov = Self::zeros(n, n); + + for k in 0..m { + for i in 0..n { + for j in 0..=i { + cov.add_element_mut(i, j, (self.get(k, i) - mu[i]) * (self.get(k, j) - mu[j])); + } + } + } + + let m_t = T::from(m - 1).unwrap(); + + for i in 0..n { + for j in 0..=i { + cov.div_element_mut(i, j, m_t); + cov.set(j, i, cov.get(i, j)); + } + } + + cov + } + } #[cfg(test)] @@ -887,4 +918,11 @@ mod tests { assert_eq!(format!("{}", a), "[[0.9, 0.4, 0.7], [0.4, 0.5, 0.3], [0.7, 0.3, 0.8]]"); } + #[test] + fn cov() { + let a = DenseMatrix::from_array(&[&[64.0, 580.0, 29.0], &[66.0, 570.0, 33.0], &[68.0, 590.0, 37.0], &[69.0, 660.0, 46.0], &[73.0, 600.0, 55.0]]); + let expected = DenseMatrix::from_array(&[&[11.5, 50.0, 34.75], &[50.0, 1250.0, 205.0], &[34.75, 205.0, 110.0]]); + assert_eq!(a.cov(), expected); + } + } diff --git a/src/linalg/nalgebra_bindings.rs b/src/linalg/nalgebra_bindings.rs index 474bd09..ba8665b 100644 --- a/src/linalg/nalgebra_bindings.rs +++ b/src/linalg/nalgebra_bindings.rs @@ -9,6 +9,7 @@ use crate::linalg::Matrix as SmartCoreMatrix; use crate::linalg::svd::SVDDecomposableMatrix; use crate::linalg::evd::EVDDecomposableMatrix; use crate::linalg::qr::QRDecomposableMatrix; +use crate::linalg::lu::LUDecomposableMatrix; impl BaseMatrix for Matrix> { @@ -318,6 +319,10 @@ impl Self { + panic!("Not implemented"); + } + } impl SVDDecomposableMatrix for Matrix> {} @@ -326,6 +331,8 @@ impl QRDecomposableMatrix for Matrix> {} +impl LUDecomposableMatrix for Matrix> {} + impl SmartCoreMatrix for Matrix> {} #[cfg(test)] diff --git a/src/linalg/ndarray_bindings.rs b/src/linalg/ndarray_bindings.rs index a22c335..2fc88d9 100644 --- a/src/linalg/ndarray_bindings.rs +++ b/src/linalg/ndarray_bindings.rs @@ -14,6 +14,7 @@ use crate::linalg::Matrix; use crate::linalg::svd::SVDDecomposableMatrix; use crate::linalg::evd::EVDDecomposableMatrix; use crate::linalg::qr::QRDecomposableMatrix; +use crate::linalg::lu::LUDecomposableMatrix; impl BaseMatrix for ArrayBase, Ix2> @@ -286,6 +287,10 @@ impl Self { + panic!("Not implemented"); + } + } impl SVDDecomposableMatrix for ArrayBase, Ix2> {} @@ -294,6 +299,8 @@ impl QRDecomposableMatrix for ArrayBase, Ix2> {} +impl LUDecomposableMatrix for ArrayBase, Ix2> {} + impl Matrix for ArrayBase, Ix2> {} #[cfg(test)] diff --git a/src/math/distance/euclidian.rs b/src/math/distance/euclidian.rs index 0a6815b..49d35d9 100644 --- a/src/math/distance/euclidian.rs +++ b/src/math/distance/euclidian.rs @@ -21,17 +21,13 @@ impl Euclidian { sum } - - pub fn distance(x: &Vec, y: &Vec) -> T { - Euclidian::squared_distance(x, y).sqrt() - } } impl Distance, T> for Euclidian { - fn distance(x: &Vec, y: &Vec) -> T { - Self::distance(x, y) + fn distance(&self, x: &Vec, y: &Vec) -> T { + Euclidian::squared_distance(x, y).sqrt() } } @@ -46,9 +42,9 @@ mod tests { let a = vec![1., 2., 3.]; let b = vec![4., 5., 6.]; - let d_arr: f64 = Euclidian::distance(&a, &b); + let l2: f64 = Euclidian{}.distance(&a, &b); - assert!((d_arr - 5.19615242).abs() < 1e-8); + assert!((l2 - 5.19615242).abs() < 1e-8); } } \ No newline at end of file diff --git a/src/math/distance/hamming.rs b/src/math/distance/hamming.rs new file mode 100644 index 0000000..882d7b8 --- /dev/null +++ b/src/math/distance/hamming.rs @@ -0,0 +1,45 @@ +use serde::{Serialize, Deserialize}; + +use crate::math::num::FloatExt; + +use super::Distance; + +#[derive(Serialize, Deserialize, Debug)] +pub struct Hamming { +} + +impl Distance, F> for Hamming { + + fn distance(&self, x: &Vec, y: &Vec) -> F { + if x.len() != y.len() { + panic!("Input vector sizes are different"); + } + + let mut dist = 0; + for i in 0..x.len() { + if x[i] != y[i]{ + dist += 1; + } + } + + F::from_i64(dist).unwrap() / F::from_usize(x.len()).unwrap() + } + +} + + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn minkowski_distance() { + let a = vec![1, 0, 0, 1, 0, 0, 1]; + let b = vec![1, 1, 0, 0, 1, 0, 1]; + + let h: f64 = Hamming{}.distance(&a, &b); + + assert!((h - 0.42857142).abs() < 1e-8); + } + +} \ No newline at end of file diff --git a/src/math/distance/mahalanobis.rs b/src/math/distance/mahalanobis.rs new file mode 100644 index 0000000..e0aeaf7 --- /dev/null +++ b/src/math/distance/mahalanobis.rs @@ -0,0 +1,97 @@ +#![allow(non_snake_case)] + +use std::marker::PhantomData; + +use serde::{Serialize, Deserialize}; + +use crate::math::num::FloatExt; + +use super::Distance; +use crate::linalg::Matrix; + +#[derive(Serialize, Deserialize, Debug)] +pub struct Mahalanobis> { + pub sigma: M, + pub sigmaInv: M, + t: PhantomData +} + +impl> Mahalanobis { + pub fn new(data: &M) -> Mahalanobis { + let sigma = data.cov(); + let sigmaInv = sigma.lu().inverse(); + Mahalanobis { + sigma: sigma, + sigmaInv: sigmaInv, + t: PhantomData + } + } + + pub fn new_from_covariance(cov: &M) -> Mahalanobis { + let sigma = cov.clone(); + let sigmaInv = sigma.lu().inverse(); + Mahalanobis { + sigma: sigma, + sigmaInv: sigmaInv, + t: PhantomData + } + } +} + +impl> Distance, T> for Mahalanobis { + + fn distance(&self, x: &Vec, y: &Vec) -> T { + let (nrows, ncols) = self.sigma.shape(); + if x.len() != nrows { + panic!("Array x[{}] has different dimension with Sigma[{}][{}].", x.len(), nrows, ncols); + } + + if y.len() != nrows { + panic!("Array y[{}] has different dimension with Sigma[{}][{}].", y.len(), nrows, ncols); + } + + println!("{}", self.sigmaInv); + + let n = x.len(); + let mut z = vec![T::zero(); n]; + for i in 0..n { + z[i] = x[i] - y[i]; + } + + // np.dot(np.dot((a-b),VI),(a-b).T) + let mut s = T::zero(); + for j in 0..n { + for i in 0..n { + s = s + self.sigmaInv.get(i, j) * z[i] * z[j]; + } + } + + s.sqrt() + } + +} + + +#[cfg(test)] +mod tests { + use super::*; + use crate::linalg::naive::dense_matrix::*; + + #[test] + fn mahalanobis_distance() { + let data = DenseMatrix::from_array(&[ + &[ 64., 580., 29.], + &[ 66., 570., 33.], + &[ 68., 590., 37.], + &[ 69., 660., 46.], + &[ 73., 600., 55.]]); + + let a = data.column_mean(); + let b = vec![66., 640., 44.]; + + let mahalanobis = Mahalanobis::new(&data); + + println!("{}", mahalanobis.distance(&a, &b)); + } + +} \ No newline at end of file diff --git a/src/math/distance/manhattan.rs b/src/math/distance/manhattan.rs new file mode 100644 index 0000000..673958d --- /dev/null +++ b/src/math/distance/manhattan.rs @@ -0,0 +1,43 @@ +use serde::{Serialize, Deserialize}; + +use crate::math::num::FloatExt; + +use super::Distance; + +#[derive(Serialize, Deserialize, Debug)] +pub struct Manhattan { +} + +impl Distance, T> for Manhattan { + + fn distance(&self, x: &Vec, y: &Vec) -> T { + if x.len() != y.len() { + panic!("Input vector sizes are different"); + } + + let mut dist = T::zero(); + for i in 0..x.len() { + dist = dist + (x[i] - y[i]).abs(); + } + + dist + } + +} + + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn manhattan_distance() { + let a = vec![1., 2., 3.]; + let b = vec![4., 5., 6.]; + + let l1: f64 = Manhattan{}.distance(&a, &b); + + assert!((l1 - 9.0).abs() < 1e-8); + } + +} \ No newline at end of file diff --git a/src/math/distance/minkowski.rs b/src/math/distance/minkowski.rs new file mode 100644 index 0000000..96d9e13 --- /dev/null +++ b/src/math/distance/minkowski.rs @@ -0,0 +1,63 @@ +use serde::{Serialize, Deserialize}; + +use crate::math::num::FloatExt; + +use super::Distance; + +#[derive(Serialize, Deserialize, Debug)] +pub struct Minkowski { + pub p: T +} + +impl Distance, T> for Minkowski { + + fn distance(&self, x: &Vec, y: &Vec) -> T { + if x.len() != y.len() { + panic!("Input vector sizes are different"); + } + if self.p < T::one() { + panic!("p must be at least 1"); + } + + let mut dist = T::zero(); + for i in 0..x.len() { + let d = (x[i] - y[i]).abs(); + dist = dist + d.powf(self.p); + } + + dist.powf(T::one()/self.p) + } + +} + + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn minkowski_distance() { + let a = vec![1., 2., 3.]; + let b = vec![4., 5., 6.]; + + let l1: f64 = Minkowski{p: 1.0}.distance(&a, &b); + let l2: f64 = Minkowski{p: 2.0}.distance(&a, &b); + let l3: f64 = Minkowski{p: 3.0}.distance(&a, &b); + + assert!((l1 - 9.0).abs() < 1e-8); + assert!((l2 - 5.19615242).abs() < 1e-8); + assert!((l3 - 4.32674871).abs() < 1e-8); + } + + #[test] + #[should_panic(expected = "p must be at least 1")] + fn minkowski_distance_negative_p() { + let a = vec![1., 2., 3.]; + let b = vec![4., 5., 6.]; + + let _: f64 = Minkowski{p: 0.0}.distance(&a, &b); + } + + + +} \ No newline at end of file diff --git a/src/math/distance/mod.rs b/src/math/distance/mod.rs index 9e29063..e359574 100644 --- a/src/math/distance/mod.rs +++ b/src/math/distance/mod.rs @@ -1,9 +1,13 @@ pub mod euclidian; +pub mod minkowski; +pub mod manhattan; +pub mod hamming; +pub mod mahalanobis; use crate::math::num::FloatExt; pub trait Distance{ - fn distance(a: &T, b: &T) -> F; + fn distance(&self, a: &T, b: &T) -> F; } pub struct Distances{ @@ -13,4 +17,16 @@ impl Distances { pub fn euclidian() -> euclidian::Euclidian{ euclidian::Euclidian {} } + + pub fn minkowski(p: T) -> minkowski::Minkowski{ + minkowski::Minkowski {p: p} + } + + pub fn manhattan() -> manhattan::Manhattan{ + manhattan::Manhattan {} + } + + pub fn hamming() -> hamming::Hamming{ + hamming::Hamming {} + } } \ No newline at end of file diff --git a/src/optimization/first_order/gradient_descent.rs b/src/optimization/first_order/gradient_descent.rs index b642534..9598e2b 100644 --- a/src/optimization/first_order/gradient_descent.rs +++ b/src/optimization/first_order/gradient_descent.rs @@ -102,9 +102,7 @@ mod tests { ls.order = FunctionOrder::THIRD; let optimizer: GradientDescent = Default::default(); - let result = optimizer.optimize(&f, &df, &x0, &ls); - - println!("{:?}", result); + let result = optimizer.optimize(&f, &df, &x0, &ls); assert!((result.f_x - 0.0).abs() < 1e-5); assert!((result.x.get(0, 0) - 1.0).abs() < 1e-2);