From 40c9915870dc7104f17b9898c1689159a485d94b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sat, 16 Aug 2014 12:16:26 +0200 Subject: [PATCH] Add the types: DVec1, .., DVec6. They are stack-allocated, dynamically sized vectors with a maximum size. This is useful for slicing small matrices, without allocation. --- Cargo.toml | 2 +- src/linalg/decompositions.rs | 12 +- src/linalg/mod.rs | 2 +- src/na.rs | 6 +- src/structs/dmat.rs | 4 +- src/structs/dvec.rs | 406 ++++------------------------ src/structs/dvec_macros.rs | 502 +++++++++++++++++++++++++++++++++++ src/structs/mat.rs | 30 +-- src/structs/mat_macros.rs | 36 +-- src/structs/mod.rs | 3 +- src/structs/vec.rs | 6 + src/structs/vec_macros.rs | 20 ++ src/tests/mat.rs | 96 +++++-- 13 files changed, 693 insertions(+), 432 deletions(-) create mode 100644 src/structs/dvec_macros.rs diff --git a/Cargo.toml b/Cargo.toml index d825a921..1874a254 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,6 +3,6 @@ name = "nalgebra" version = "0.1.0" authors = [ "Sébastien Crozet " ] # FIXME: add the contributors. -[[lib]] +[lib] name = "nalgebra" path = "src/lib.rs" diff --git a/src/linalg/decompositions.rs b/src/linalg/decompositions.rs index f666866e..44471f57 100644 --- a/src/linalg/decompositions.rs +++ b/src/linalg/decompositions.rs @@ -34,11 +34,12 @@ pub fn householder_matrix + Norm, - M: Clone + Eye + ColSlice + Transpose - + Indexable<(uint, uint), N> + Mul> - (m: &M) -> (M, M) { +pub fn qr + Norm, + M: Clone + Eye + ColSlice + Transpose + + Indexable<(uint, uint), N> + Mul> + (m: &M) + -> (M, M) { let (rows, cols) = m.shape(); assert!(rows >= cols); let mut q : M = Eye::new_identity(rows); @@ -67,4 +68,3 @@ pub fn decomp_qr Indexable<(uint, uint), N> for DMat { #[inline] unsafe fn unsafe_at(&self, rowcol: (uint, uint)) -> N { let (row, col) = rowcol; - (*self.mij.as_slice().unsafe_ref(self.offset(row, col))).clone() + (*self.mij.as_slice().unsafe_get(self.offset(row, col))).clone() } #[inline] @@ -524,7 +524,7 @@ impl ColSlice> for DMat { // we can init from slice thanks to the matrix being column major let start= self.offset(row_start, col_id); let stop = self.offset(row_end, col_id); - let slice = DVec::from_vec( + let slice = DVec::from_slice( row_end - row_start, self.mij.slice(start, stop)); slice } diff --git a/src/structs/dvec.rs b/src/structs/dvec.rs index 9b2fd38a..500147a6 100644 --- a/src/structs/dvec.rs +++ b/src/structs/dvec.rs @@ -11,97 +11,13 @@ use std::iter::FromIterator; use traits::geometry::{Dot, Norm}; use traits::structure::{Iterable, IterableMut, Indexable}; -/// Vector with a dimension unknown at compile-time. +/// Heap allocated, dynamically sized vector. #[deriving(Eq, PartialEq, Show, Clone)] pub struct DVec { /// Components of the vector. Contains as much elements as the vector dimension. pub at: Vec } -double_dispatch_binop_decl_trait!(DVec, DVecMulRhs) -double_dispatch_binop_decl_trait!(DVec, DVecDivRhs) -double_dispatch_binop_decl_trait!(DVec, DVecAddRhs) -double_dispatch_binop_decl_trait!(DVec, DVecSubRhs) - -mul_redispatch_impl!(DVec, DVecMulRhs) -div_redispatch_impl!(DVec, DVecDivRhs) -add_redispatch_impl!(DVec, DVecAddRhs) -sub_redispatch_impl!(DVec, DVecSubRhs) - -impl DVec { - /// Builds a vector filled with zeros. - /// - /// # Arguments - /// * `dim` - The dimension of the vector. - #[inline] - pub fn new_zeros(dim: uint) -> DVec { - DVec::from_elem(dim, Zero::zero()) - } - - /// Tests if all components of the vector are zeroes. - #[inline] - pub fn is_zero(&self) -> bool { - self.at.iter().all(|e| e.is_zero()) - } -} - -impl Indexable for DVec { - - fn at(&self, i: uint) -> N { - assert!(i < self.at.len()); - unsafe { - self.unsafe_at(i) - } - } - - fn set(&mut self, i: uint, val: N) { - assert!(i < self.at.len()); - unsafe { - self.unsafe_set(i, val); - } - } - - fn swap(&mut self, i: uint, j: uint) { - assert!(i < self.at.len()); - assert!(j < self.at.len()); - self.at.as_mut_slice().swap(i, j); - } - - fn shape(&self) -> uint { - self.at.len() - } - - #[inline] - unsafe fn unsafe_at(&self, i: uint) -> N { - (*self.at.as_slice().unsafe_ref(i)).clone() - } - - #[inline] - unsafe fn unsafe_set(&mut self, i: uint, val: N) { - *self.at.as_mut_slice().unsafe_mut_ref(i) = val - } - -} - -impl DVec { - /// Builds a vector filled with ones. - /// - /// # Arguments - /// * `dim` - The dimension of the vector. - #[inline] - pub fn new_ones(dim: uint) -> DVec { - DVec::from_elem(dim, One::one()) - } -} - -impl DVec { - /// Builds a vector filled with random values. - #[inline] - pub fn new_random(dim: uint) -> DVec { - DVec::from_fn(dim, |_| rand::random()) - } -} - impl DVec { /// Creates an uninitialized vec. #[inline] @@ -113,24 +29,6 @@ impl DVec { at: vec } } - - /// Gets a reference to of this vector data. - #[inline] - pub fn as_vec<'r>(&'r self) -> &'r [N] { - self.at.as_slice() - } - - /// Gets a mutable reference to of this vector data. - #[inline] - pub fn as_mut_vec<'r>(&'r mut self) -> &'r mut [N] { - self.at.as_mut_slice() - } - - /// Extracts this vector data. - #[inline] - pub fn to_vec(self) -> Vec { - self.at - } } impl DVec { @@ -144,7 +42,7 @@ impl DVec { /// /// The vector must have at least `dim` elements. #[inline] - pub fn from_vec(dim: uint, vec: &[N]) -> DVec { + pub fn from_slice(dim: uint, vec: &[N]) -> DVec { assert!(dim <= vec.len()); DVec { @@ -161,27 +59,6 @@ impl DVec { } } -impl Collection for DVec { - #[inline] - fn len(&self) -> uint { - self.at.len() - } -} - -impl Iterable for DVec { - #[inline] - fn iter<'l>(&'l self) -> Items<'l, N> { - self.at.iter() - } -} - -impl IterableMut for DVec { - #[inline] - fn mut_iter<'l>(&'l mut self) -> MutItems<'l, N> { - self.at.mut_iter() - } -} - impl FromIterator for DVec { #[inline] fn from_iter>(mut param: I) -> DVec { @@ -195,258 +72,71 @@ impl FromIterator for DVec { } } -impl + DVecMulRhs>> DVec { - /// Computes the canonical basis for the given dimension. A canonical basis is a set of - /// vectors, mutually orthogonal, with all its component equal to 0.0 except one which is equal - /// to 1.0. - pub fn canonical_basis_with_dim(dim: uint) -> Vec> { - let mut res : Vec> = Vec::new(); - for i in range(0u, dim) { - let mut basis_element : DVec = DVec::new_zeros(dim); - - *basis_element.at.get_mut(i) = One::one(); - - res.push(basis_element); - } - - res - } - - /// Computes a basis of the space orthogonal to the vector. If the input vector is of dimension - /// `n`, this will return `n - 1` vectors. - pub fn orthogonal_subspace_basis(&self) -> Vec> { - // compute the basis of the orthogonal subspace using Gram-Schmidt - // orthogonalization algorithm - let dim = self.at.len(); - let mut res : Vec> = Vec::new(); - - for i in range(0u, dim) { - let mut basis_element : DVec = DVec::new_zeros(self.at.len()); - - *basis_element.at.get_mut(i) = One::one(); - - if res.len() == dim - 1 { - break; - } - - let mut elt = basis_element.clone(); - - elt = elt - self * Dot::dot(&basis_element, self); - - for v in res.iter() { - elt = elt - v * Dot::dot(&elt, v) - }; - - if !ApproxEq::approx_eq(&Norm::sqnorm(&elt), &Zero::zero()) { - res.push(Norm::normalize_cpy(&elt)); - } - } - - assert!(res.len() == dim - 1); - - res +impl Collection for DVec { + #[inline] + fn len(&self) -> uint { + self.at.len() } } -impl> DVecAddRhs> for DVec { - #[inline] - fn binop(left: &DVec, right: &DVec) -> DVec { - assert!(left.at.len() == right.at.len()); - DVec { - at: left.at.iter().zip(right.at.iter()).map(|(a, b)| *a + *b).collect() - } - } +dvec_impl!(DVec, DVecMulRhs, DVecDivRhs, DVecAddRhs, DVecSubRhs) + +/// Stack-allocated, dynamically sized vector with a maximum size of 1. +pub struct DVec1 { + at: [N, ..1], + dim: uint } -impl> DVecSubRhs> for DVec { - #[inline] - fn binop(left: &DVec, right: &DVec) -> DVec { - assert!(left.at.len() == right.at.len()); - DVec { - at: left.at.iter().zip(right.at.iter()).map(|(a, b)| *a - *b).collect() - } - } +small_dvec_impl!(DVec1, 1, DVec1MulRhs, DVec1DivRhs, DVec1AddRhs, DVec1SubRhs, 0) +small_dvec_from_impl!(DVec1, 1, Zero::zero()) + + +/// Stack-allocated, dynamically sized vector with a maximum size of 2. +pub struct DVec2 { + at: [N, ..2], + dim: uint } -impl> Neg> for DVec { - #[inline] - fn neg(&self) -> DVec { - DVec { at: self.at.iter().map(|a| -*a).collect() } - } +small_dvec_impl!(DVec2, 2, DVec2MulRhs, DVec2DivRhs, DVec2AddRhs, DVec2SubRhs, 0, 1) +small_dvec_from_impl!(DVec2, 2, Zero::zero(), Zero::zero()) + + +/// Stack-allocated, dynamically sized vector with a maximum size of 3. +pub struct DVec3 { + at: [N, ..3], + dim: uint } -impl Dot for DVec { - #[inline] - fn dot(a: &DVec, b: &DVec) -> N { - assert!(a.at.len() == b.at.len()); +small_dvec_impl!(DVec3, 3, DVec3MulRhs, DVec3DivRhs, DVec3AddRhs, DVec3SubRhs, 0, 1, 2) +small_dvec_from_impl!(DVec3, 3, Zero::zero(), Zero::zero(), Zero::zero()) - let mut res: N = Zero::zero(); - for i in range(0u, a.at.len()) { - res = res + unsafe { a.unsafe_at(i) * b.unsafe_at(i) }; - } - - res - } - - #[inline] - fn sub_dot(a: &DVec, b: &DVec, c: &DVec) -> N { - let mut res: N = Zero::zero(); - - for i in range(0u, a.at.len()) { - res = res + unsafe { (a.unsafe_at(i) - b.unsafe_at(i)) * c.unsafe_at(i) }; - } - - res - } +/// Stack-allocated, dynamically sized vector with a maximum size of 4. +pub struct DVec4 { + at: [N, ..4], + dim: uint } -impl Norm for DVec { - #[inline] - fn sqnorm(v: &DVec) -> N { - Dot::dot(v, v) - } +small_dvec_impl!(DVec4, 4, DVec4MulRhs, DVec4DivRhs, DVec4AddRhs, DVec4SubRhs, 0, 1, 2, 3) +small_dvec_from_impl!(DVec4, 4, Zero::zero(), Zero::zero(), Zero::zero(), Zero::zero()) - #[inline] - fn norm(v: &DVec) -> N { - Norm::sqnorm(v).sqrt() - } - #[inline] - fn normalize_cpy(v: &DVec) -> DVec { - let mut res : DVec = v.clone(); - - let _ = res.normalize(); - - res - } - - #[inline] - fn normalize(&mut self) -> N { - let l = Norm::norm(self); - - for i in range(0u, self.at.len()) { - *self.at.get_mut(i) = self.at[i] / l; - } - - l - } +/// Stack-allocated, dynamically sized vector with a maximum size of 5. +pub struct DVec5 { + at: [N, ..5], + dim: uint } -impl> ApproxEq for DVec { - #[inline] - fn approx_epsilon(_: Option>) -> N { - ApproxEq::approx_epsilon(None::) - } +small_dvec_impl!(DVec5, 5, DVec5MulRhs, DVec5DivRhs, DVec5AddRhs, DVec5SubRhs, 0, 1, 2, 3, 4) +small_dvec_from_impl!(DVec5, 5, Zero::zero(), Zero::zero(), Zero::zero(), Zero::zero(), Zero::zero()) - #[inline] - fn approx_eq(a: &DVec, b: &DVec) -> bool { - let mut zip = a.at.iter().zip(b.at.iter()); - zip.all(|(a, b)| ApproxEq::approx_eq(a, b)) - } - - #[inline] - fn approx_eq_eps(a: &DVec, b: &DVec, epsilon: &N) -> bool { - let mut zip = a.at.iter().zip(b.at.iter()); - - zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon)) - } +/// Stack-allocated, dynamically sized vector with a maximum size of 6. +pub struct DVec6 { + at: [N, ..6], + dim: uint } -macro_rules! scalar_mul_impl ( - ($n: ident) => ( - impl DVecMulRhs<$n, DVec<$n>> for $n { - #[inline] - fn binop(left: &DVec<$n>, right: &$n) -> DVec<$n> { - DVec { at: left.at.iter().map(|a| a * *right).collect() } - } - } - ) -) - -macro_rules! scalar_div_impl ( - ($n: ident) => ( - impl DVecDivRhs<$n, DVec<$n>> for $n { - #[inline] - fn binop(left: &DVec<$n>, right: &$n) -> DVec<$n> { - DVec { at: left.at.iter().map(|a| a / *right).collect() } - } - } - ) -) - -macro_rules! scalar_add_impl ( - ($n: ident) => ( - impl DVecAddRhs<$n, DVec<$n>> for $n { - #[inline] - fn binop(left: &DVec<$n>, right: &$n) -> DVec<$n> { - DVec { at: left.at.iter().map(|a| a + *right).collect() } - } - } - ) -) - -macro_rules! scalar_sub_impl ( - ($n: ident) => ( - impl DVecSubRhs<$n, DVec<$n>> for $n { - #[inline] - fn binop(left: &DVec<$n>, right: &$n) -> DVec<$n> { - DVec { at: left.at.iter().map(|a| a - *right).collect() } - } - } - ) -) - -scalar_mul_impl!(f64) -scalar_mul_impl!(f32) -scalar_mul_impl!(u64) -scalar_mul_impl!(u32) -scalar_mul_impl!(u16) -scalar_mul_impl!(u8) -scalar_mul_impl!(i64) -scalar_mul_impl!(i32) -scalar_mul_impl!(i16) -scalar_mul_impl!(i8) -scalar_mul_impl!(uint) -scalar_mul_impl!(int) - -scalar_div_impl!(f64) -scalar_div_impl!(f32) -scalar_div_impl!(u64) -scalar_div_impl!(u32) -scalar_div_impl!(u16) -scalar_div_impl!(u8) -scalar_div_impl!(i64) -scalar_div_impl!(i32) -scalar_div_impl!(i16) -scalar_div_impl!(i8) -scalar_div_impl!(uint) -scalar_div_impl!(int) - -scalar_add_impl!(f64) -scalar_add_impl!(f32) -scalar_add_impl!(u64) -scalar_add_impl!(u32) -scalar_add_impl!(u16) -scalar_add_impl!(u8) -scalar_add_impl!(i64) -scalar_add_impl!(i32) -scalar_add_impl!(i16) -scalar_add_impl!(i8) -scalar_add_impl!(uint) -scalar_add_impl!(int) - -scalar_sub_impl!(f64) -scalar_sub_impl!(f32) -scalar_sub_impl!(u64) -scalar_sub_impl!(u32) -scalar_sub_impl!(u16) -scalar_sub_impl!(u8) -scalar_sub_impl!(i64) -scalar_sub_impl!(i32) -scalar_sub_impl!(i16) -scalar_sub_impl!(i8) -scalar_sub_impl!(uint) -scalar_sub_impl!(int) +small_dvec_impl!(DVec6, 6, DVec6MulRhs, DVec6DivRhs, DVec6AddRhs, DVec6SubRhs, 0, 1, 2, 3, 4, 5) +small_dvec_from_impl!(DVec6, 6, Zero::zero(), Zero::zero(), Zero::zero(), Zero::zero(), Zero::zero(), Zero::zero()) diff --git a/src/structs/dvec_macros.rs b/src/structs/dvec_macros.rs new file mode 100644 index 00000000..7c6d40f9 --- /dev/null +++ b/src/structs/dvec_macros.rs @@ -0,0 +1,502 @@ +#![macro_escape] + +macro_rules! dvec_impl( + ($dvec: ident, $mul: ident, $div: ident, $add: ident, $sub: ident) => ( + double_dispatch_binop_decl_trait!($dvec, $mul) + double_dispatch_binop_decl_trait!($dvec, $div) + double_dispatch_binop_decl_trait!($dvec, $add) + double_dispatch_binop_decl_trait!($dvec, $sub) + + mul_redispatch_impl!($dvec, $mul) + div_redispatch_impl!($dvec, $div) + add_redispatch_impl!($dvec, $add) + sub_redispatch_impl!($dvec, $sub) + + impl $dvec { + /// Builds a vector filled with zeros. + /// + /// # Arguments + /// * `dim` - The dimension of the vector. + #[inline] + pub fn new_zeros(dim: uint) -> $dvec { + $dvec::from_elem(dim, Zero::zero()) + } + + /// Tests if all components of the vector are zeroes. + #[inline] + pub fn is_zero(&self) -> bool { + self.as_slice().iter().all(|e| e.is_zero()) + } + } + + impl $dvec { + /// Slices this vector. + #[inline] + pub fn as_slice<'a>(&'a self) -> &'a [N] { + self.at.slice_to(self.len()) + } + + /// Mutably slices this vector. + #[inline] + pub fn as_mut_slice<'a>(&'a mut self) -> &'a mut [N] { + let len = self.len(); + self.at.mut_slice_to(len) + } + } + + impl Indexable for $dvec { + #[inline] + fn at(&self, i: uint) -> N { + assert!(i < self.len()); + unsafe { + self.unsafe_at(i) + } + } + + #[inline] + fn set(&mut self, i: uint, val: N) { + assert!(i < self.len()); + unsafe { + self.unsafe_set(i, val); + } + } + + #[inline] + fn swap(&mut self, i: uint, j: uint) { + assert!(i < self.len()); + assert!(j < self.len()); + self.as_mut_slice().swap(i, j); + } + + #[inline] + fn shape(&self) -> uint { + self.len() + } + + #[inline] + unsafe fn unsafe_at(&self, i: uint) -> N { + (*self.at.as_slice().unsafe_get(i)).clone() + } + + #[inline] + unsafe fn unsafe_set(&mut self, i: uint, val: N) { + *self.at.as_mut_slice().unsafe_mut_ref(i) = val + } + + } + + impl $dvec { + /// Builds a vector filled with ones. + /// + /// # Arguments + /// * `dim` - The dimension of the vector. + #[inline] + pub fn new_ones(dim: uint) -> $dvec { + $dvec::from_elem(dim, One::one()) + } + } + + impl $dvec { + /// Builds a vector filled with random values. + #[inline] + pub fn new_random(dim: uint) -> $dvec { + $dvec::from_fn(dim, |_| rand::random()) + } + } + + impl Iterable for $dvec { + #[inline] + fn iter<'l>(&'l self) -> Items<'l, N> { + self.as_slice().iter() + } + } + + impl IterableMut for $dvec { + #[inline] + fn mut_iter<'l>(&'l mut self) -> MutItems<'l, N> { + self.as_mut_slice().mut_iter() + } + } + + impl + $mul>> $dvec { + /// Computes the canonical basis for the given dimension. A canonical basis is a set of + /// vectors, mutually orthogonal, with all its component equal to 0.0 except one which is equal + /// to 1.0. + pub fn canonical_basis_with_dim(dim: uint) -> Vec<$dvec> { + let mut res : Vec<$dvec> = Vec::new(); + + for i in range(0u, dim) { + let mut basis_element : $dvec = $dvec::new_zeros(dim); + + basis_element.set(i, One::one()); + + res.push(basis_element); + } + + res + } + + /// Computes a basis of the space orthogonal to the vector. If the input vector is of dimension + /// `n`, this will return `n - 1` vectors. + pub fn orthogonal_subspace_basis(&self) -> Vec<$dvec> { + // compute the basis of the orthogonal subspace using Gram-Schmidt + // orthogonalization algorithm + let dim = self.len(); + let mut res : Vec<$dvec> = Vec::new(); + + for i in range(0u, dim) { + let mut basis_element : $dvec = $dvec::new_zeros(self.len()); + + basis_element.set(i, One::one()); + + if res.len() == dim - 1 { + break; + } + + let mut elt = basis_element.clone(); + + elt = elt - self * Dot::dot(&basis_element, self); + + for v in res.iter() { + elt = elt - v * Dot::dot(&elt, v) + }; + + if !ApproxEq::approx_eq(&Norm::sqnorm(&elt), &Zero::zero()) { + res.push(Norm::normalize_cpy(&elt)); + } + } + + assert!(res.len() == dim - 1); + + res + } + } + + impl + Zero> $add> for $dvec { + #[inline] + fn binop(left: &$dvec, right: &$dvec) -> $dvec { + assert!(left.len() == right.len()); + FromIterator::from_iter(left.as_slice().iter().zip(right.as_slice().iter()).map(|(a, b)| *a + *b)) + } + } + + impl + Zero> $sub> for $dvec { + #[inline] + fn binop(left: &$dvec, right: &$dvec) -> $dvec { + assert!(left.len() == right.len()); + FromIterator::from_iter(left.as_slice().iter().zip(right.as_slice().iter()).map(|(a, b)| *a - *b)) + } + } + + impl + Zero> Neg<$dvec> for $dvec { + #[inline] + fn neg(&self) -> $dvec { + FromIterator::from_iter(self.as_slice().iter().map(|a| -*a)) + } + } + + impl Dot for $dvec { + #[inline] + fn dot(a: &$dvec, b: &$dvec) -> N { + assert!(a.len() == b.len()); + + let mut res: N = Zero::zero(); + + for i in range(0u, a.len()) { + res = res + unsafe { a.unsafe_at(i) * b.unsafe_at(i) }; + } + + res + } + + #[inline] + fn sub_dot(a: &$dvec, b: &$dvec, c: &$dvec) -> N { + let mut res: N = Zero::zero(); + + for i in range(0u, a.len()) { + res = res + unsafe { (a.unsafe_at(i) - b.unsafe_at(i)) * c.unsafe_at(i) }; + } + + res + } + } + + impl Norm for $dvec { + #[inline] + fn sqnorm(v: &$dvec) -> N { + Dot::dot(v, v) + } + + #[inline] + fn norm(v: &$dvec) -> N { + Norm::sqnorm(v).sqrt() + } + + #[inline] + fn normalize_cpy(v: &$dvec) -> $dvec { + let mut res : $dvec = v.clone(); + + let _ = res.normalize(); + + res + } + + #[inline] + fn normalize(&mut self) -> N { + let l = Norm::norm(self); + + for n in self.as_mut_slice().mut_iter() { + *n = *n / l; + } + + l + } + } + + impl> ApproxEq for $dvec { + #[inline] + fn approx_epsilon(_: Option<$dvec>) -> N { + ApproxEq::approx_epsilon(None::) + } + + #[inline] + fn approx_eq(a: &$dvec, b: &$dvec) -> bool { + let mut zip = a.as_slice().iter().zip(b.as_slice().iter()); + + zip.all(|(a, b)| ApproxEq::approx_eq(a, b)) + } + + #[inline] + fn approx_eq_eps(a: &$dvec, b: &$dvec, epsilon: &N) -> bool { + let mut zip = a.as_slice().iter().zip(b.as_slice().iter()); + + zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon)) + } + } + + dvec_scalar_mul_impl!($dvec, f64, $mul) + dvec_scalar_mul_impl!($dvec, f32, $mul) + dvec_scalar_mul_impl!($dvec, u64, $mul) + dvec_scalar_mul_impl!($dvec, u32, $mul) + dvec_scalar_mul_impl!($dvec, u16, $mul) + dvec_scalar_mul_impl!($dvec, u8, $mul) + dvec_scalar_mul_impl!($dvec, i64, $mul) + dvec_scalar_mul_impl!($dvec, i32, $mul) + dvec_scalar_mul_impl!($dvec, i16, $mul) + dvec_scalar_mul_impl!($dvec, i8, $mul) + dvec_scalar_mul_impl!($dvec, uint, $mul) + dvec_scalar_mul_impl!($dvec, int, $mul) + + dvec_scalar_div_impl!($dvec, f64, $div) + dvec_scalar_div_impl!($dvec, f32, $div) + dvec_scalar_div_impl!($dvec, u64, $div) + dvec_scalar_div_impl!($dvec, u32, $div) + dvec_scalar_div_impl!($dvec, u16, $div) + dvec_scalar_div_impl!($dvec, u8, $div) + dvec_scalar_div_impl!($dvec, i64, $div) + dvec_scalar_div_impl!($dvec, i32, $div) + dvec_scalar_div_impl!($dvec, i16, $div) + dvec_scalar_div_impl!($dvec, i8, $div) + dvec_scalar_div_impl!($dvec, uint, $div) + dvec_scalar_div_impl!($dvec, int, $div) + + dvec_scalar_add_impl!($dvec, f64, $add) + dvec_scalar_add_impl!($dvec, f32, $add) + dvec_scalar_add_impl!($dvec, u64, $add) + dvec_scalar_add_impl!($dvec, u32, $add) + dvec_scalar_add_impl!($dvec, u16, $add) + dvec_scalar_add_impl!($dvec, u8, $add) + dvec_scalar_add_impl!($dvec, i64, $add) + dvec_scalar_add_impl!($dvec, i32, $add) + dvec_scalar_add_impl!($dvec, i16, $add) + dvec_scalar_add_impl!($dvec, i8, $add) + dvec_scalar_add_impl!($dvec, uint, $add) + dvec_scalar_add_impl!($dvec, int, $add) + + dvec_scalar_sub_impl!($dvec, f64, $sub) + dvec_scalar_sub_impl!($dvec, f32, $sub) + dvec_scalar_sub_impl!($dvec, u64, $sub) + dvec_scalar_sub_impl!($dvec, u32, $sub) + dvec_scalar_sub_impl!($dvec, u16, $sub) + dvec_scalar_sub_impl!($dvec, u8, $sub) + dvec_scalar_sub_impl!($dvec, i64, $sub) + dvec_scalar_sub_impl!($dvec, i32, $sub) + dvec_scalar_sub_impl!($dvec, i16, $sub) + dvec_scalar_sub_impl!($dvec, i8, $sub) + dvec_scalar_sub_impl!($dvec, uint, $sub) + dvec_scalar_sub_impl!($dvec, int, $sub) + ) +) + +macro_rules! dvec_scalar_mul_impl ( + ($dvec: ident, $n: ident, $mul: ident) => ( + impl $mul<$n, $dvec<$n>> for $n { + #[inline] + fn binop(left: &$dvec<$n>, right: &$n) -> $dvec<$n> { + FromIterator::from_iter(left.as_slice().iter().map(|a| a * *right)) + } + } + ) +) + +macro_rules! dvec_scalar_div_impl ( + ($dvec: ident, $n: ident, $div: ident) => ( + impl $div<$n, $dvec<$n>> for $n { + #[inline] + fn binop(left: &$dvec<$n>, right: &$n) -> $dvec<$n> { + FromIterator::from_iter(left.as_slice().iter().map(|a| a / *right)) + } + } + ) +) + +macro_rules! dvec_scalar_add_impl ( + ($dvec: ident, $n: ident, $add: ident) => ( + impl $add<$n, $dvec<$n>> for $n { + #[inline] + fn binop(left: &$dvec<$n>, right: &$n) -> $dvec<$n> { + FromIterator::from_iter(left.as_slice().iter().map(|a| a + *right)) + } + } + ) +) + +macro_rules! dvec_scalar_sub_impl ( + ($dvec: ident, $n: ident, $sub: ident) => ( + impl $sub<$n, $dvec<$n>> for $n { + #[inline] + fn binop(left: &$dvec<$n>, right: &$n) -> $dvec<$n> { + FromIterator::from_iter(left.as_slice().iter().map(|a| a - *right)) + } + } + ) +) + +macro_rules! small_dvec_impl ( + ($dvec: ident, $dim: expr, $mul: ident, $div: ident, $add: ident, $sub: ident $(,$idx: expr)*) => ( + impl Collection for $dvec { + #[inline] + fn len(&self) -> uint { + self.dim + } + } + + impl PartialEq for $dvec { + #[inline] + fn eq(&self, other: &$dvec) -> bool { + if self.len() != other.len() { + return false; // FIXME: fail instead? + } + + for (a, b) in self.as_slice().iter().zip(other.as_slice().iter()) { + if *a != *b { + return false; + } + } + + true + } + } + + impl Clone for $dvec { + fn clone(&self) -> $dvec { + let at: [N, ..$dim] = [ $( self.at[$idx].clone(), )* ]; + + $dvec { + at: at, + dim: self.dim + } + } + } + + dvec_impl!($dvec, $mul, $div, $add, $sub) + ) +) + +macro_rules! small_dvec_from_impl ( + ($dvec: ident, $dim: expr $(,$zeros: expr)*) => ( + impl $dvec { + /// Builds a vector filled with a constant. + #[inline] + pub fn from_elem(dim: uint, elem: N) -> $dvec { + assert!(dim <= $dim); + + let mut at: [N, ..$dim] = [ $( $zeros, )* ]; + + for n in at.mut_slice_to(dim).mut_iter() { + *n = elem.clone(); + } + + $dvec { + at: at, + dim: dim + } + } + } + + impl $dvec { + /// Builds a vector filled with the components provided by a vector. + /// + /// The vector must have at least `dim` elements. + #[inline] + pub fn from_slice(dim: uint, vec: &[N]) -> $dvec { + assert!(dim <= vec.len() && dim <= $dim); + + // FIXME: not safe. + let mut at: [N, ..$dim] = [ $( $zeros, )* ]; + + for (curr, other) in vec.iter().zip(at.mut_iter()) { + *other = curr.clone(); + } + + $dvec { + at: at, + dim: dim + } + } + } + + impl $dvec { + /// Builds a vector filled with the result of a function. + #[inline(always)] + pub fn from_fn(dim: uint, f: |uint| -> N) -> $dvec { + assert!(dim <= $dim); + + let mut at: [N, ..$dim] = [ $( $zeros, )* ]; + + for i in range(0, dim) { + at[i] = f(i); + } + + $dvec { + at: at, + dim: dim + } + } + } + + impl FromIterator for $dvec { + #[inline] + fn from_iter>(mut param: I) -> $dvec { + let mut at: [N, ..$dim] = [ $( $zeros, )* ]; + + let mut dim = 0; + + for n in param { + if dim == $dim { + break; + } + + at[dim] = n; + + dim = dim + 1; + } + + $dvec { + at: at, + dim: dim + } + } + } + ) +) diff --git a/src/structs/mat.rs b/src/structs/mat.rs index 248e18eb..7eeae95b 100644 --- a/src/structs/mat.rs +++ b/src/structs/mat.rs @@ -6,9 +6,9 @@ use std::mem; use std::num::{One, Zero}; use traits::operations::ApproxEq; use std::slice::{Items, MutItems}; -use structs::vec::{Vec1, Vec2, Vec3, Vec4, Vec5, Vec6, Vec1MulRhs, Vec4MulRhs, - Vec5MulRhs, Vec6MulRhs}; -use structs::dvec::DVec; +use structs::vec::{Vec1, Vec2, Vec3, Vec4, Vec5, Vec6, + Vec1MulRhs, Vec4MulRhs, Vec5MulRhs, Vec6MulRhs}; +use structs::dvec::{DVec1, DVec2, DVec3, DVec4, DVec5, DVec6}; use traits::structure::{Cast, Row, Col, Iterable, IterableMut, Dim, Indexable, Eye, ColSlice, RowSlice}; @@ -118,8 +118,8 @@ transpose_impl!(Mat1, 1) approx_eq_impl!(Mat1) row_impl!(Mat1, Vec1, 1) col_impl!(Mat1, Vec1, 1) -col_slice_impl!(Mat1, Vec1, 1) -row_slice_impl!(Mat1, Vec1, 1) +col_slice_impl!(Mat1, Vec1, DVec1, 1) +row_slice_impl!(Mat1, Vec1, DVec1, 1) to_homogeneous_impl!(Mat1, Mat2, 1, 2) from_homogeneous_impl!(Mat1, Mat2, 1, 2) outer_impl!(Vec1, Mat1) @@ -219,8 +219,8 @@ transpose_impl!(Mat2, 2) approx_eq_impl!(Mat2) row_impl!(Mat2, Vec2, 2) col_impl!(Mat2, Vec2, 2) -col_slice_impl!(Mat2, Vec2, 2) -row_slice_impl!(Mat2, Vec2, 2) +col_slice_impl!(Mat2, Vec2, DVec2, 2) +row_slice_impl!(Mat2, Vec2, DVec2, 2) to_homogeneous_impl!(Mat2, Mat3, 2, 3) from_homogeneous_impl!(Mat2, Mat3, 2, 3) outer_impl!(Vec2, Mat2) @@ -334,8 +334,8 @@ transpose_impl!(Mat3, 3) approx_eq_impl!(Mat3) // (specialized) row_impl!(Mat3, Vec3, 3) // (specialized) col_impl!(Mat3, Vec3, 3) -col_slice_impl!(Mat3, Vec3, 3) -row_slice_impl!(Mat3, Vec3, 3) +col_slice_impl!(Mat3, Vec3, DVec3, 3) +row_slice_impl!(Mat3, Vec3, DVec3, 3) to_homogeneous_impl!(Mat3, Mat4, 3, 4) from_homogeneous_impl!(Mat3, Mat4, 3, 4) outer_impl!(Vec3, Mat3) @@ -501,8 +501,8 @@ transpose_impl!(Mat4, 4) approx_eq_impl!(Mat4) row_impl!(Mat4, Vec4, 4) col_impl!(Mat4, Vec4, 4) -col_slice_impl!(Mat4, Vec4, 4) -row_slice_impl!(Mat4, Vec4, 4) +col_slice_impl!(Mat4, Vec4, DVec4, 4) +row_slice_impl!(Mat4, Vec4, DVec4, 4) to_homogeneous_impl!(Mat4, Mat5, 4, 5) from_homogeneous_impl!(Mat4, Mat5, 4, 5) outer_impl!(Vec4, Mat4) @@ -684,8 +684,8 @@ transpose_impl!(Mat5, 5) approx_eq_impl!(Mat5) row_impl!(Mat5, Vec5, 5) col_impl!(Mat5, Vec5, 5) -col_slice_impl!(Mat5, Vec5, 5) -row_slice_impl!(Mat5, Vec5, 5) +col_slice_impl!(Mat5, Vec5, DVec5, 5) +row_slice_impl!(Mat5, Vec5, DVec5, 5) to_homogeneous_impl!(Mat5, Mat6, 5, 6) from_homogeneous_impl!(Mat5, Mat6, 5, 6) outer_impl!(Vec5, Mat5) @@ -919,6 +919,6 @@ transpose_impl!(Mat6, 6) approx_eq_impl!(Mat6) row_impl!(Mat6, Vec6, 6) col_impl!(Mat6, Vec6, 6) -col_slice_impl!(Mat6, Vec6, 6) -row_slice_impl!(Mat6, Vec6, 6) +col_slice_impl!(Mat6, Vec6, DVec6, 6) +row_slice_impl!(Mat6, Vec6, DVec6, 6) outer_impl!(Vec6, Mat6) diff --git a/src/structs/mat_macros.rs b/src/structs/mat_macros.rs index 0ca91a3e..9274d292 100644 --- a/src/structs/mat_macros.rs +++ b/src/structs/mat_macros.rs @@ -226,20 +226,12 @@ macro_rules! indexable_impl( ) macro_rules! col_slice_impl( - ($t: ident, $tv: ident, $dim: expr) => ( - impl ColSlice> for $t { - fn col_slice(& self, cid: uint, rstart: uint, rend: uint) -> DVec { - assert!(cid < $dim); - assert!(rstart < rend); - assert!(rend <= $dim); + ($t: ident, $tv: ident, $slice: ident, $dim: expr) => ( + impl ColSlice<$slice> for $t { + fn col_slice(&self, cid: uint, rstart: uint, rend: uint) -> $slice { let col = self.col(cid); - let res = DVec::from_vec( - rend - rstart, - unsafe { - mem::transmute::<&$tv, & [N, ..$dim]> - (&col).slice(rstart, rend) - }); - res + + $slice::from_slice(rend - rstart, col.as_slice().slice(rstart, rend)) } } ) @@ -275,20 +267,12 @@ macro_rules! row_impl( ) macro_rules! row_slice_impl( - ($t: ident, $tv: ident, $dim: expr) => ( - impl RowSlice> for $t { - fn row_slice(& self, rid: uint, cstart: uint, cend: uint) -> DVec { - assert!(rid < $dim); - assert!(cstart < cend); - assert!(cend <= $dim); + ($t: ident, $tv: ident, $slice: ident, $dim: expr) => ( + impl RowSlice<$slice> for $t { + fn row_slice(&self, rid: uint, cstart: uint, cend: uint) -> $slice { let row = self.row(rid); - let res = DVec::from_vec( - cend - cstart, - unsafe { - mem::transmute::<&$tv, & [N, ..$dim]> - (&row).slice(cstart, cend) - }); - res + + $slice::from_slice(cend - cstart, row.as_slice().slice(cstart, cend)) } } ) diff --git a/src/structs/mod.rs b/src/structs/mod.rs index 3d33f068..0aef1e8c 100644 --- a/src/structs/mod.rs +++ b/src/structs/mod.rs @@ -1,7 +1,7 @@ //! Data structures and implementations. pub use self::dmat::DMat; -pub use self::dvec::DVec; +pub use self::dvec::{DVec, DVec1, DVec2, DVec3, DVec4, DVec5, DVec6}; pub use self::vec::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6}; pub use self::mat::{Identity, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6}; pub use self::rot::{Rot2, Rot3, Rot4}; @@ -18,6 +18,7 @@ pub use self::mat::{Mat1MulRhs, Mat2MulRhs, Mat3MulRhs, Mat4MulRhs, Mat5MulRhs, mod metal; mod dmat; +mod dvec_macros; mod dvec; mod vec_macros; mod vec; diff --git a/src/structs/vec.rs b/src/structs/vec.rs index 56bc1678..c6abd67f 100644 --- a/src/structs/vec.rs +++ b/src/structs/vec.rs @@ -38,6 +38,7 @@ new_impl!(Vec1, x) ord_impl!(Vec1, x) vec_axis_impl!(Vec1, x) vec_cast_impl!(Vec1, Vec1Cast, x) +as_slice_impl!(Vec1, 1) indexable_impl!(Vec1, 1) at_fast_impl!(Vec1, 1) new_repeat_impl!(Vec1, val, x) @@ -135,6 +136,7 @@ new_impl!(Vec2, x, y) ord_impl!(Vec2, x, y) vec_axis_impl!(Vec2, x, y) vec_cast_impl!(Vec2, Vec2Cast, x, y) +as_slice_impl!(Vec2, 2) indexable_impl!(Vec2, 2) at_fast_impl!(Vec2, 2) new_repeat_impl!(Vec2, val, x, y) @@ -234,6 +236,7 @@ new_impl!(Vec3, x, y, z) ord_impl!(Vec3, x, y, z) vec_axis_impl!(Vec3, x, y, z) vec_cast_impl!(Vec3, Vec3Cast, x, y, z) +as_slice_impl!(Vec3, 3) indexable_impl!(Vec3, 3) at_fast_impl!(Vec3, 3) new_repeat_impl!(Vec3, val, x, y, z) @@ -339,6 +342,7 @@ new_impl!(Vec4, x, y, z, w) ord_impl!(Vec4, x, y, z, w) vec_axis_impl!(Vec4, x, y, z, w) vec_cast_impl!(Vec4, Vec4Cast, x, y, z, w) +as_slice_impl!(Vec4, 4) indexable_impl!(Vec4, 4) at_fast_impl!(Vec4, 4) new_repeat_impl!(Vec4, val, x, y, z, w) @@ -442,6 +446,7 @@ new_impl!(Vec5, x, y, z, w, a) ord_impl!(Vec5, x, y, z, w, a) vec_axis_impl!(Vec5, x, y, z, w, a) vec_cast_impl!(Vec5, Vec5Cast, x, y, z, w, a) +as_slice_impl!(Vec5, 5) indexable_impl!(Vec5, 5) at_fast_impl!(Vec5, 5) new_repeat_impl!(Vec5, val, x, y, z, w, a) @@ -547,6 +552,7 @@ new_impl!(Vec6, x, y, z, w, a, b) ord_impl!(Vec6, x, y, z, w, a, b) vec_axis_impl!(Vec6, x, y, z, w, a, b) vec_cast_impl!(Vec6, Vec6Cast, x, y, z, w, a, b) +as_slice_impl!(Vec6, 6) indexable_impl!(Vec6, 6) at_fast_impl!(Vec6, 6) new_repeat_impl!(Vec6, val, x, y, z, w, a, b) diff --git a/src/structs/vec_macros.rs b/src/structs/vec_macros.rs index fd4779c8..60eea525 100644 --- a/src/structs/vec_macros.rs +++ b/src/structs/vec_macros.rs @@ -15,6 +15,26 @@ macro_rules! new_impl( ) ) +macro_rules! as_slice_impl( + ($t: ident, $dim: expr) => ( + impl $t { + /// Slices this vector. + pub fn as_slice<'a>(&'a self) -> &'a [N] { + unsafe { + mem::transmute::<&$t, &[N, ..$dim]>(self).as_slice() + } + } + + /// Mutably slices this vector. + pub fn as_mut_slice<'a>(&'a mut self) -> &'a mut [N] { + unsafe { + mem::transmute::<&mut $t, &mut [N, ..$dim]>(self).as_mut_slice() + } + } + } + ) +) + macro_rules! at_fast_impl( ($t: ident, $dim: expr) => ( impl $t { diff --git a/src/tests/mat.rs b/src/tests/mat.rs index 5a9da8e2..b5f9b677 100644 --- a/src/tests/mat.rs +++ b/src/tests/mat.rs @@ -2,7 +2,6 @@ use std::num::{Float, abs}; use std::rand::random; use na::{Vec1, Vec3, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6, Rot3, DMat, DVec, Indexable}; use na; -use na::decomp_qr; use std::cmp::{min, max}; macro_rules! test_inv_mat_impl( @@ -25,12 +24,12 @@ macro_rules! test_transpose_mat_impl( ); ) -macro_rules! test_decomp_qr_impl( +macro_rules! test_qr_impl( ($t: ty) => ( for _ in range(0u, 10000) { let randmat : $t = random(); - let (q, r) = decomp_qr(&randmat); + let (q, r) = na::qr(&randmat); let recomp = q * r; assert!(na::approx_eq(&randmat, &recomp)); @@ -38,6 +37,35 @@ macro_rules! test_decomp_qr_impl( ); ) +macro_rules! test_eigen_qr_impl( + ($t: ty) => { + for _ in range(0u, 10000) { + let randmat : $t = random(); + + let (eigenvalues, eigenvectors) = na::eigen_qr(&randmat, &Float::epsilon(), 1000); + + // FIXME: provide a method to initialize a matrix from its diagonal! + let diag: $t = na::zero(); + + for i in range(0, na::dim::<$t>()) { + diag.set((i, i), eigenvalues.at(i)); + } + + let recomp = na::transpose(&eigenvectors) * diag * eigenvectors; + + println!("mat: {}", randmat); + println!("eigenvectors: {}", eigenvectors); + println!("eigenvalues: {}", eigenvalues); + println!("recomp: {}", recomp); + + assert!(false); + fail!("what!"); + + assert!(na::approx_eq(&randmat, &recomp)); + } + } +) + #[test] fn test_transpose_mat1() { test_transpose_mat_impl!(Mat1); @@ -139,7 +167,7 @@ fn test_mean_dmat() { ] ); - assert!(na::approx_eq(&na::mean(&mat), &DVec::from_vec(3, [4.0f64, 5.0, 6.0]))); + assert!(na::approx_eq(&na::mean(&mat), &DVec::from_slice(3, [4.0f64, 5.0, 6.0]))); } #[test] @@ -223,14 +251,14 @@ fn test_dmat_from_vec() { } #[test] -fn test_decomp_qr() { +fn test_qr() { for _ in range(0u, 10) { let dim1: uint = random(); let dim2: uint = random(); let rows = min(40, max(dim1, dim2)); let cols = min(40, min(dim1, dim2)); let randmat: DMat = DMat::new_random(rows, cols); - let (q, r) = decomp_qr(&randmat); + let (q, r) = na::qr(&randmat); let recomp = q * r; assert!(na::approx_eq(&randmat, &recomp)); @@ -238,40 +266,68 @@ fn test_decomp_qr() { } #[test] -fn test_decomp_qr_mat1() { - test_decomp_qr_impl!(Mat1); +fn test_qr_mat1() { + test_qr_impl!(Mat1); } #[test] -fn test_decomp_qr_mat2() { - test_decomp_qr_impl!(Mat2); +fn test_qr_mat2() { + test_qr_impl!(Mat2); } #[test] -fn test_decomp_qr_mat3() { - test_decomp_qr_impl!(Mat3); +fn test_qr_mat3() { + test_qr_impl!(Mat3); } #[test] -fn test_decomp_qr_mat4() { - test_decomp_qr_impl!(Mat4); +fn test_qr_mat4() { + test_qr_impl!(Mat4); } #[test] -fn test_decomp_qr_mat5() { - test_decomp_qr_impl!(Mat5); +fn test_qr_mat5() { + test_qr_impl!(Mat5); } #[test] -fn test_decomp_qr_mat6() { - test_decomp_qr_impl!(Mat6); +fn test_qr_mat6() { + test_qr_impl!(Mat6); } +// #[test] +// fn test_eigen_qr_mat1() { +// test_eigen_qr_impl!(Mat1); +// } +// +// #[test] +// fn test_eigen_qr_mat2() { +// test_eigen_qr_impl!(Mat2); +// } +// +// #[test] +// fn test_eigen_qr_mat3() { +// test_eigen_qr_impl!(Mat3); +// } +// +// #[test] +// fn test_eigen_qr_mat4() { +// test_eigen_qr_impl!(Mat4); +// } +// +// #[test] +// fn test_eigen_qr_mat5() { +// test_eigen_qr_impl!(Mat5); +// } +// +// #[test] +// fn test_eigen_qr_mat6() { +// test_eigen_qr_impl!(Mat6); +// } #[test] fn test_from_fn() { - let actual: DMat = DMat::from_fn( 3, 4, - |i,j| 10*i + j); + let actual: DMat = DMat::from_fn(3, 4, |i, j| 10 * i + j); let expected: DMat = DMat::from_row_vec(3, 4, [0_0, 0_1, 0_2, 0_3, 1_0, 1_1, 1_2, 1_3,