Add mean and covariance computation for DMat.

This also fixes the transposition for rectangular DMat, and add scalar
addition/subtraction/multiplication/division for DMat.
This commit is contained in:
Sébastien Crozet 2013-09-22 14:22:17 +02:00
parent d12d23d2c7
commit 4ffe274b03
7 changed files with 341 additions and 36 deletions

View File

@ -1,5 +1,7 @@
//! Matrix with dimensions unknown at compile-time.
#[doc(hidden)]; // we hide doc to not have to document the $trhs double dispatch trait.
use std::rand::Rand;
use std::rand;
use std::num::{One, Zero};
@ -7,28 +9,28 @@ use std::vec;
use std::cmp::ApproxEq;
use std::util;
use dvec::{DVec, DVecMulRhs};
use traits::operations::{Inv, Transpose};
use traits::operations::{Inv, Transpose, Mean, Cov};
#[doc(hidden)]
mod metal;
/// Matrix with dimensions unknown at compile-time.
#[deriving(Eq, ToStr, Clone)]
#[deriving(Eq, Clone)]
pub struct DMat<N> {
priv nrows: uint,
priv ncols: uint,
priv mij: ~[N]
}
/// Trait of object `o` which can be multiplied by a `DMat` `d`: `d * o`.
pub trait DMatMulRhs<N, Res> {
/// Multiplies a `DMat` by `Self`.
fn binop(left: &DMat<N>, right: &Self) -> Res;
}
double_dispatch_binop_decl_trait!(DMat, DMatMulRhs)
double_dispatch_binop_decl_trait!(DMat, DMatDivRhs)
double_dispatch_binop_decl_trait!(DMat, DMatAddRhs)
double_dispatch_binop_decl_trait!(DMat, DMatSubRhs)
impl<N, Rhs: DMatMulRhs<N, Res>, Res> Mul<Rhs, Res> for DMat<N> {
#[inline(always)]
fn mul(&self, other: &Rhs) -> Res {
DMatMulRhs::binop(self, other)
}
}
mul_redispatch_impl!(DMat, DMatMulRhs)
div_redispatch_impl!(DMat, DMatDivRhs)
add_redispatch_impl!(DMat, DMatAddRhs)
sub_redispatch_impl!(DMat, DMatSubRhs)
impl<N> DMat<N> {
/// Creates an uninitialized matrix.
@ -89,6 +91,20 @@ impl<N: Clone> DMat<N> {
mij: vec::from_elem(nrows * ncols, val)
}
}
/// Builds a matrix filled with the components provided by a vector.
///
/// The vector must have at least `nrows * ncols` elements.
#[inline]
pub fn from_vec(nrows: uint, ncols: uint, vec: &[N]) -> DMat<N> {
assert!(nrows * ncols <= vec.len());
DMat {
nrows: nrows,
ncols: ncols,
mij: vec.slice_to(nrows * ncols).to_owned()
}
}
}
impl<N> DMat<N> {
@ -116,7 +132,7 @@ impl<N> DMat<N> {
/// Transforms this matrix into an array. This consumes the matrix and is O(1).
#[inline]
pub fn to_array(self) -> ~[N] {
pub fn to_vec(self) -> ~[N] {
self.mij
}
}
@ -350,24 +366,88 @@ Inv for DMat<N> {
impl<N: Clone> Transpose for DMat<N> {
#[inline]
fn transposed(&self) -> DMat<N> {
let mut res = self.clone();
if self.nrows == self.ncols {
let mut res = self.clone();
res.transpose();
res.transpose();
res
res
}
else {
let mut res = unsafe { DMat::new_uninitialized(self.ncols, self.nrows) };
for i in range(0u, self.nrows) {
for j in range(0u, self.ncols) {
unsafe {
res.set_fast(j, i, self.at_fast(i, j))
}
}
}
res
}
}
#[inline]
fn transpose(&mut self) {
for i in range(1u, self.nrows) {
for j in range(0u, self.ncols - 1) {
let off_i_j = self.offset(i, j);
let off_j_i = self.offset(j, i);
if self.nrows == self.ncols {
for i in range(1u, self.nrows) {
for j in range(0u, self.ncols - 1) {
let off_i_j = self.offset(i, j);
let off_j_i = self.offset(j, i);
self.mij.swap(off_i_j, off_j_i);
self.mij.swap(off_i_j, off_j_i);
}
}
util::swap(&mut self.nrows, &mut self.ncols);
}
else {
// FIXME: implement a better algorithm which does that in-place.
*self = self.transposed();
}
}
}
impl<N: Num + NumCast + Clone> Mean<DVec<N>> for DMat<N> {
fn mean(&self) -> DVec<N> {
let mut res: DVec<N> = DVec::new_zeros(self.ncols);
let normalizer: N = NumCast::from(1.0f64 / NumCast::from(self.nrows));
for i in range(0u, self.nrows) {
for j in range(0u, self.ncols) {
unsafe {
let acc = res.at_fast(j) + self.at_fast(i, j) * normalizer;
res.set_fast(j, acc);
}
}
}
util::swap(&mut self.nrows, &mut self.ncols);
res
}
}
impl<N: Clone + Num + NumCast + DMatDivRhs<N, DMat<N>> + ToStr > Cov<DMat<N>> for DMat<N> {
// FIXME: this could be heavily optimized, removing all temporaries by merging loops.
fn cov(&self) -> DMat<N> {
assert!(self.nrows > 1);
let mut centered = unsafe { DMat::new_uninitialized(self.nrows, self.ncols) };
let mean = self.mean();
// FIXME: use the rows iterator when available
for i in range(0u, self.nrows) {
for j in range(0u, self.ncols) {
unsafe {
centered.set_fast(i, j, self.at_fast(i, j) - mean.at_fast(j));
}
}
}
// FIXME: return a triangular matrix?
let normalizer: N = NumCast::from(self.nrows() - 1);
// FIXME: this will do 2 allocations for temporaries!
(centered.transposed() * centered) / normalizer
}
}
@ -398,3 +478,137 @@ impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
}
}
}
macro_rules! scalar_mul_impl (
($n: ident) => (
impl DMatMulRhs<$n, DMat<$n>> for $n {
#[inline]
fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
DMat {
nrows: left.nrows,
ncols: left.ncols,
mij: left.mij.iter().map(|a| a * *right).collect()
}
}
}
)
)
macro_rules! scalar_div_impl (
($n: ident) => (
impl DMatDivRhs<$n, DMat<$n>> for $n {
#[inline]
fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
DMat {
nrows: left.nrows,
ncols: left.ncols,
mij: left.mij.iter().map(|a| a / *right).collect()
}
}
}
)
)
macro_rules! scalar_add_impl (
($n: ident) => (
impl DMatAddRhs<$n, DMat<$n>> for $n {
#[inline]
fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
DMat {
nrows: left.nrows,
ncols: left.ncols,
mij: left.mij.iter().map(|a| a + *right).collect()
}
}
}
)
)
macro_rules! scalar_sub_impl (
($n: ident) => (
impl DMatSubRhs<$n, DMat<$n>> for $n {
#[inline]
fn binop(left: &DMat<$n>, right: &$n) -> DMat<$n> {
DMat {
nrows: left.nrows,
ncols: left.ncols,
mij: left.mij.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!(float)
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!(float)
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!(float)
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!(float)
scalar_sub_impl!(uint)
scalar_sub_impl!(int)
impl<N: ToStr + Clone> ToStr for DMat<N> {
fn to_str(&self) -> ~str {
let mut res = ~"DMat ";
res = res + self.nrows.to_str() + " " + self.ncols.to_str() + " {\n";
for i in range(0u, self.nrows) {
for j in range(0u, self.ncols) {
res = res + " " + unsafe { self.at_fast(i, j).to_str() };
}
res = res + "\n";
}
res = res + "}";
res
}
}

View File

@ -12,6 +12,7 @@ use std::iter::FromIterator;
use traits::geometry::{Dot, Norm, Translation};
use traits::structure::{Iterable, IterableMut};
#[doc(hidden)]
mod metal;
/// Vector with a dimension unknown at compile-time.
@ -85,6 +86,16 @@ impl<N> DVec<N> {
at: vec
}
}
#[inline]
pub unsafe fn set_fast(&mut self, i: uint, val: N) {
*self.at.unsafe_mut_ref(i) = val
}
#[inline]
pub fn to_vec(self) -> ~[N] {
self.at
}
}
impl<N: Clone> DVec<N> {
@ -93,6 +104,18 @@ impl<N: Clone> DVec<N> {
pub fn from_elem(dim: uint, elem: N) -> DVec<N> {
DVec { at: vec::from_elem(dim, elem) }
}
/// Builds a vector filled with the components provided by a vector.
///
/// The vector must have at least `dim` elements.
#[inline]
pub fn from_vec(dim: uint, vec: &[N]) -> DVec<N> {
assert!(dim <= vec.len());
DVec {
at: vec.slice_to(dim).to_owned()
}
}
}
impl<N> DVec<N> {

View File

@ -10,7 +10,7 @@ use vec::*;
// traits
pub use traits::structure::{Mat, Dim, Indexable, Iterable, IterableMut, MatCast, Row, Col};
pub use traits::operations::{Absolute, ScalarSub, ScalarAdd, Inv, RMul, Transpose};
pub use traits::operations::{Absolute, ScalarSub, ScalarAdd, Inv, RMul, Transpose, Mean, Cov};
pub use traits::geometry::{Rotation, RotationMatrix, Rotate, Transformation, Transform,
Translation, Translate, ToHomogeneous, FromHomogeneous,
RotationWithTranslation, AbsoluteRotate};

View File

@ -1,5 +1,7 @@
#[macro_escape];
#[doc(hidden)]; // we hide doc to not have to document the $trhs double dispatch trait.
// Create the traits needed to do fancy operator oveloading.
// This is a meta version of
// http://smallcultfollowing.com/babysteps/blog/2012/10/04/refining-traits-slash-impls/

View File

@ -87,12 +87,12 @@ fn test_inv_mat6() {
#[test]
fn test_rotation2() {
do 10000.times {
let randmat: Rotmat<Mat2<f64>> = One::one();
let ang = &Vec1::new(abs::<f64>(random()) % Real::pi());
do 10000.times {
let randmat: Rotmat<Mat2<f64>> = One::one();
let ang = &Vec1::new(abs::<f64>(random()) % Real::pi());
assert!(randmat.rotated(ang).rotation().approx_eq(ang));
}
assert!(randmat.rotated(ang).rotation().approx_eq(ang));
}
}
#[test]
@ -104,12 +104,54 @@ fn test_index_mat2() {
#[test]
fn test_inv_rotation3() {
do 10000.times {
let randmat: Rotmat<Mat3<f64>> = One::one();
let dir: Vec3<f64> = random();
let ang = &(dir.normalized() * (abs::<f64>(random()) % Real::pi()));
let rot = randmat.rotated(ang);
do 10000.times {
let randmat: Rotmat<Mat3<f64>> = One::one();
let dir: Vec3<f64> = random();
let ang = &(dir.normalized() * (abs::<f64>(random()) % Real::pi()));
let rot = randmat.rotated(ang);
assert!((rot.transposed() * rot).approx_eq(&One::one()));
}
assert!((rot.transposed() * rot).approx_eq(&One::one()));
}
}
#[test]
fn test_mean_dmat() {
let mat = DMat::from_vec(
3,
3,
[
1.0f64, 2.0, 3.0,
4.0f64, 5.0, 6.0,
7.0f64, 8.0, 9.0,
]
);
assert!(mat.mean().approx_eq(&DVec::from_vec(3, [4.0f64, 5.0, 6.0])));
}
#[test]
fn test_cov_dmat() {
let mat = DMat::from_vec(
5,
3,
[
4.0, 2.0, 0.60,
4.2, 2.1, 0.59,
3.9, 2.0, 0.58,
4.3, 2.1, 0.62,
4.1, 2.2, 0.63
]
);
let expected = DMat::from_vec(
3,
3,
[
0.025, 0.0075, 0.00175,
0.0075, 0.007, 0.00135,
0.00175, 0.00135, 0.00043
]
);
assert!(mat.cov().approx_eq(&expected));
}

View File

@ -31,6 +31,22 @@ pub trait Outer<V, M> {
fn outer(&self, other: &V) -> M;
}
/// Trait for computing the covariance of a set of data.
pub trait Cov<M> {
/// Computes the covariance of the obsevations stored by `self`:
/// * for matrices, observations are stored in its rows.
/// * for vectors, observations are stored in its components (thus are 1-dimensional).
fn cov(&self) -> M;
}
/// Trait for computing the covariance of a set of data.
pub trait Mean<N> {
/// Computes the mean of the observations stored by `self`.
/// * for matrices, observations are stored in its rows.
/// * for vectors, observations are stored in its components (thus are 1-dimensional).
fn mean(&self) -> N;
}
// XXX: those two traits should not exist since there is generalized operator overloading of Add
// and Sub.
// However, using the same trait multiple time as a trait bound (ex: impl<T: Add<N, V> + Add<V, V>)

View File

@ -105,16 +105,24 @@ pub trait Row<R> {
fn row(&self, i: uint) -> R;
/// Writes the `i`-th row of `self`.
fn set_row(&mut self, i: uint, R);
// FIXME: add iterators on rows: this could be a very good way to generalize _and_ optimize
// a lot of operations.
}
/// Trait to access columns of a matrix or vector.
pub trait Col<C> {
/// The number of column of this matrix or vector.
fn num_cols(&self) -> uint;
/// Reads the `i`-th column of `self`.
fn col(&self, i: uint) -> C;
/// Writes the `i`-th column of `self`.
fn set_col(&mut self, i: uint, C);
// FIXME: add iterators on columns: this could be a very good way to generalize _and_ optimize
// a lot of operations.
}
/// Trait of objects having a spacial dimension known at compile time.