diff --git a/src/nalgebra.rc b/src/nalgebra.rc index 01e71d42..42111bed 100644 --- a/src/nalgebra.rc +++ b/src/nalgebra.rc @@ -39,6 +39,8 @@ pub mod dim3 /// n-dimensional linear algebra (slower). pub mod ndim { + // pub mod dmat; + pub mod dvec; pub mod nvec; pub mod nmat; } diff --git a/src/ndim/dvec.rs b/src/ndim/dvec.rs new file mode 100644 index 00000000..1d1f1a38 --- /dev/null +++ b/src/ndim/dvec.rs @@ -0,0 +1,233 @@ +use core::num::{Zero, One, Algebraic}; +use core::vec::{map_zip, map, all2, len, from_elem, all}; +use core::cmp::ApproxEq; +use traits::ring::Ring; +use traits::division_ring::DivisionRing; +use traits::dot::Dot; +use traits::sub_dot::SubDot; +use traits::norm::Norm; +use traits::translation::Translation; +use traits::workarounds::scalar_op::{ScalarMul, ScalarDiv, ScalarAdd, ScalarSub}; + +#[deriving(Eq, ToStr, Clone)] +pub struct DVec +{ + at: ~[T] +} + +pub fn zero_with_dim(dim: uint) -> DVec +{ DVec { at: from_elem(dim, Zero::zero::()) } } + +pub fn is_zero(vec: &DVec) -> bool +{ all(vec.at, |e| e.is_zero()) } + +// FIXME: is Clone needed? +impl> DVec +{ + pub fn canonical_basis_with_dim(dim: uint) -> ~[DVec] + { + let mut res : ~[DVec] = ~[]; + + for uint::range(0u, dim) |i| + { + let mut basis_element : DVec = zero_with_dim(dim); + + basis_element.at[i] = One::one(); + + res.push(basis_element); + } + + res + } + + pub fn orthogonal_subspace_basis(&self) -> ~[DVec] + { + // compute the basis of the orthogonal subspace using Gram-Schmidt + // orthogonalization algorithm + let dim = len(self.at); + let mut res : ~[DVec] = ~[]; + + for uint::range(0u, dim) |i| + { + let mut basis_element : DVec = zero_with_dim(len(self.at)); + + basis_element.at[i] = One::one(); + + if (res.len() == dim - 1) + { break; } + + let mut elt = basis_element.clone(); + + elt -= self.scalar_mul(&basis_element.dot(self)); + + for res.each |v| + { elt -= v.scalar_mul(&elt.dot(v)) }; + + if (!elt.sqnorm().approx_eq(&Zero::zero())) + { res.push(elt.normalized()); } + } + + assert!(res.len() == dim - 1); + + res + } +} + +impl> Add, DVec> for DVec +{ + fn add(&self, other: &DVec) -> DVec + { + assert!(len(self.at) == len(other.at)); + DVec { at: map_zip(self.at, other.at, | a, b | { *a + *b }) } + } +} + +impl> Sub, DVec> for DVec +{ + fn sub(&self, other: &DVec) -> DVec + { + assert!(len(self.at) == len(other.at)); + DVec { at: map_zip(self.at, other.at, | a, b | *a - *b) } + } +} + +impl> Neg> for DVec +{ + fn neg(&self) -> DVec + { DVec { at: map(self.at, |a| -a) } } +} + +impl +Dot for DVec +{ + fn dot(&self, other: &DVec) -> T + { + assert!(len(self.at) == len(other.at)); + + let mut res = Zero::zero::(); + + for uint::range(0u, len(self.at)) |i| + { res += self.at[i] * other.at[i]; } + + res + } +} + +impl SubDot for DVec +{ + fn sub_dot(&self, a: &DVec, b: &DVec) -> T + { + let mut res = Zero::zero::(); + + for uint::range(0u, len(self.at)) |i| + { res += (self.at[i] - a.at[i]) * b.at[i]; } + + res + } +} + +impl> +ScalarMul for DVec +{ + fn scalar_mul(&self, s: &T) -> DVec + { DVec { at: map(self.at, |a| a * *s) } } + + fn scalar_mul_inplace(&mut self, s: &T) + { + for uint::range(0u, len(self.at)) |i| + { self.at[i] *= *s; } + } +} + + +impl> +ScalarDiv for DVec +{ + fn scalar_div(&self, s: &T) -> DVec + { DVec { at: map(self.at, |a| a / *s) } } + + fn scalar_div_inplace(&mut self, s: &T) + { + for uint::range(0u, len(self.at)) |i| + { self.at[i] /= *s; } + } +} + +impl> +ScalarAdd for DVec +{ + fn scalar_add(&self, s: &T) -> DVec + { DVec { at: map(self.at, |a| a + *s) } } + + fn scalar_add_inplace(&mut self, s: &T) + { + for uint::range(0u, len(self.at)) |i| + { self.at[i] += *s; } + } +} + +impl> +ScalarSub for DVec +{ + fn scalar_sub(&self, s: &T) -> DVec + { DVec { at: map(self.at, |a| a - *s) } } + + fn scalar_sub_inplace(&mut self, s: &T) + { + for uint::range(0u, len(self.at)) |i| + { self.at[i] -= *s; } + } +} + +impl> Translation> for DVec +{ + fn translation(&self) -> DVec + { self.clone() } + + fn translated(&self, t: &DVec) -> DVec + { self + *t } + + fn translate(&mut self, t: &DVec) + { *self = *self + *t; } +} + +impl +Norm for DVec +{ + fn sqnorm(&self) -> T + { self.dot(self) } + + fn norm(&self) -> T + { self.sqnorm().sqrt() } + + fn normalized(&self) -> DVec + { + let mut res : DVec = self.clone(); + + res.normalize(); + + res + } + + fn normalize(&mut self) -> T + { + let l = self.norm(); + + for uint::range(0u, len(self.at)) |i| + { self.at[i] /= l; } + + l + } +} + +impl> ApproxEq for DVec +{ + fn approx_epsilon() -> T + { ApproxEq::approx_epsilon::() } + + fn approx_eq(&self, other: &DVec) -> bool + { all2(self.at, other.at, |a, b| a.approx_eq(b)) } + + fn approx_eq_eps(&self, other: &DVec, epsilon: &T) -> bool + { all2(self.at, other.at, |a, b| a.approx_eq_eps(b, epsilon)) } +} diff --git a/src/ndim/nmat.rs b/src/ndim/nmat.rs index f83d4a38..a93c25a9 100644 --- a/src/ndim/nmat.rs +++ b/src/ndim/nmat.rs @@ -109,7 +109,7 @@ RMul> for NMat for uint::range(0u, dim) |i| { for uint::range(0u, dim) |j| - { res.at[i] = res.at[i] + other.at[j] * self[(i, j)]; } + { res.at.at[i] = res.at.at[i] + other.at.at[j] * self[(i, j)]; } } res @@ -127,7 +127,7 @@ LMul> for NMat for uint::range(0u, dim) |i| { for uint::range(0u, dim) |j| - { res.at[i] = res.at[i] + other.at[j] * self[(j, i)]; } + { res.at.at[i] = res.at.at[i] + other.at.at[j] * self[(j, i)]; } } res @@ -158,8 +158,8 @@ Inv for NMat for uint::range(0u, dim) |k| { // search a non-zero value on the k-th column - // FIXME: is it worth it to spend some more time searching for the max - // instead? + // FIXME: would it be worth it to spend some more time searching for the + // max instead? // FIXME: this is kind of uggly… // … but we cannot use position_between since we are iterating on one diff --git a/src/ndim/nvec.rs b/src/ndim/nvec.rs index cb04ceb1..b887908d 100644 --- a/src/ndim/nvec.rs +++ b/src/ndim/nvec.rs @@ -1,7 +1,8 @@ -use core::num::{Zero, One, Algebraic}; +use core::num::{Zero, Algebraic}; use core::rand::{Rand, Rng, RngUtil}; -use core::vec::{map_zip, from_elem, map, all, all2}; +use core::vec::{map}; use core::cmp::ApproxEq; +use ndim::dvec::{DVec, zero_with_dim, is_zero}; use traits::basis::Basis; use traits::ring::Ring; use traits::division_ring::DivisionRing; @@ -15,14 +16,12 @@ use traits::workarounds::scalar_op::{ScalarMul, ScalarDiv, ScalarAdd, ScalarSub} // D is a phantom parameter, used only as a dimensional token. // Its allows use to encode the vector dimension at the type-level. // It can be anything implementing the Dim trait. However, to avoid confusion, -// using d0, d1, d2, d3 and d4 tokens are prefered. +// using d0, d1, d2, d3, ..., d7 (or your own dn) are prefered. // FIXME: it might be possible to implement type-level integers and use them // here? #[deriving(Eq, ToStr)] pub struct NVec -{ - at: ~[T] -} +{ at: DVec } impl Dim for NVec @@ -40,59 +39,42 @@ impl Clone for NVec impl> Add, NVec> for NVec { fn add(&self, other: &NVec) -> NVec - { NVec { at: map_zip(self.at, other.at, | a, b | { *a + *b }) } } + { NVec { at: self.at + other.at } } } impl> Sub, NVec> for NVec { fn sub(&self, other: &NVec) -> NVec - { NVec { at: map_zip(self.at, other.at, | a, b | *a - *b) } } + { NVec { at: self.at - other.at } } } impl> Neg> for NVec { fn neg(&self) -> NVec - { NVec { at: map(self.at, |a| -a) } } + { NVec { at: -self.at } } } impl Dot for NVec { fn dot(&self, other: &NVec) -> T - { - let mut res = Zero::zero::(); - - for uint::range(0u, Dim::dim::()) |i| - { res += self.at[i] * other.at[i]; } - - res - } + { self.at.dot(&other.at) } } impl SubDot for NVec { fn sub_dot(&self, a: &NVec, b: &NVec) -> T - { - let mut res = Zero::zero::(); - - for uint::range(0u, Dim::dim::()) |i| - { res += (self.at[i] - a.at[i]) * b.at[i]; } - - res - } + { self.at.sub_dot(&a.at, &b.at) } } impl> ScalarMul for NVec { fn scalar_mul(&self, s: &T) -> NVec - { NVec { at: map(self.at, |a| a * *s) } } + { NVec { at: self.at.scalar_mul(s) } } fn scalar_mul_inplace(&mut self, s: &T) - { - for uint::range(0u, Dim::dim::()) |i| - { self.at[i] *= *s; } - } + { self.at.scalar_mul_inplace(s) } } @@ -100,39 +82,30 @@ impl> ScalarDiv for NVec { fn scalar_div(&self, s: &T) -> NVec - { NVec { at: map(self.at, |a| a / *s) } } + { NVec { at: self.at.scalar_div(s) } } fn scalar_div_inplace(&mut self, s: &T) - { - for uint::range(0u, Dim::dim::()) |i| - { self.at[i] /= *s; } - } + { self.at.scalar_div_inplace(s) } } impl> ScalarAdd for NVec { fn scalar_add(&self, s: &T) -> NVec - { NVec { at: map(self.at, |a| a + *s) } } + { NVec { at: self.at.scalar_add(s) } } fn scalar_add_inplace(&mut self, s: &T) - { - for uint::range(0u, Dim::dim::()) |i| - { self.at[i] += *s; } - } + { self.at.scalar_add_inplace(s) } } impl> ScalarSub for NVec { fn scalar_sub(&self, s: &T) -> NVec - { NVec { at: map(self.at, |a| a - *s) } } + { NVec { at: self.at.scalar_sub(s) } } fn scalar_sub_inplace(&mut self, s: &T) - { - for uint::range(0u, Dim::dim::()) |i| - { self.at[i] -= *s; } - } + { self.scalar_sub_inplace(s) } } impl> Translation> for NVec @@ -166,14 +139,7 @@ Norm for NVec } fn normalize(&mut self) -> T - { - let l = self.norm(); - - for uint::range(0u, Dim::dim::()) |i| - { self.at[i] /= l; } - - l - } + { self.at.normalize() } } impl { fn canonical_basis() -> ~[NVec] - { - let dim = Dim::dim::(); - let mut res : ~[NVec] = ~[]; - - for uint::range(0u, dim) |i| - { - let mut basis_element : NVec = Zero::zero(); - - basis_element.at[i] = One::one(); - - res.push(basis_element); - } - - res - } + { map(DVec::canonical_basis_with_dim(Dim::dim::()), |&e| NVec { at: e }) } fn orthogonal_subspace_basis(&self) -> ~[NVec] - { - // compute the basis of the orthogonal subspace using Gram-Schmidt - // orthogonalization algorithm - let dim = Dim::dim::(); - let mut res : ~[NVec] = ~[]; - - for uint::range(0u, dim) |i| - { - let mut basis_element : NVec = Zero::zero(); - - basis_element.at[i] = One::one(); - - if (res.len() == dim - 1) - { break; } - - let mut elt = basis_element.clone(); - - elt -= self.scalar_mul(&basis_element.dot(self)); - - for res.each |v| - { elt -= v.scalar_mul(&elt.dot(v)) }; - - if (!elt.sqnorm().approx_eq(&Zero::zero())) - { res.push(elt.normalized()); } - } - - assert!(res.len() == dim - 1); - - res - } + { map(self.at.orthogonal_subspace_basis(), |&e| NVec { at: e }) } } // FIXME: I dont really know how te generalize the cross product int @@ -241,16 +164,10 @@ Basis for NVec impl Zero for NVec { fn zero() -> NVec - { - let _0 = Zero::zero(); - - NVec { at: from_elem(Dim::dim::(), _0) } - } + { NVec { at: zero_with_dim(Dim::dim::()) } } fn is_zero(&self) -> bool - { - all(self.at, |e| e.is_zero()) - } + { is_zero(&self.at) } } impl> ApproxEq for NVec @@ -259,10 +176,10 @@ impl> ApproxEq for NVec { ApproxEq::approx_epsilon::() } fn approx_eq(&self, other: &NVec) -> bool - { all2(self.at, other.at, |a, b| a.approx_eq(b)) } + { self.at.approx_eq(&other.at) } fn approx_eq_eps(&self, other: &NVec, epsilon: &T) -> bool - { all2(self.at, other.at, |a, b| a.approx_eq_eps(b, epsilon)) } + { self.at.approx_eq_eps(&other.at, epsilon) } } impl Rand for NVec @@ -273,7 +190,7 @@ impl Rand for NVec let mut res : NVec = Zero::zero(); for uint::range(0u, dim) |i| - { res.at[i] = rng.gen() } + { res.at.at[i] = rng.gen() } res }