M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'misc'

This commit is contained in:
Sébastien Crozet 2016-01-10 15:46:06 +01:00
parent b63e60560f 91c4b58bbb
commit bfd97beffe
16 changed files with 1200 additions and 800 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Cargo.toml
View File

@ -1,6 +1,6 @@
[package] [package]
name = "nalgebra" name = "nalgebra"
version = "0.4.0" version = "0.5.0"
authors = [ "Sébastien Crozet <developer@crozet.re>" ] # FIXME: add the contributors. authors = [ "Sébastien Crozet <developer@crozet.re>" ] # FIXME: add the contributors.
description = "Linear algebra library for computer physics, computer graphics and general low-dimensional linear algebra for Rust." description = "Linear algebra library for computer physics, computer graphics and general low-dimensional linear algebra for Rust."

4
benches/dmat.rs
View File

@ -13,7 +13,7 @@ macro_rules! bench_mul_dmat(
let a: DMat<f64> = DMat::new_random($nrows, $ncols); let a: DMat<f64> = DMat::new_random($nrows, $ncols);
let mut b: DMat<f64> = DMat::new_random($nrows, $ncols); let mut b: DMat<f64> = DMat::new_random($nrows, $ncols);
for _ in (0usize .. 1000) { for _ in 0usize .. 1000 {
// XXX: the clone here is highly undesirable! // XXX: the clone here is highly undesirable!
b = a.clone() * b; b = a.clone() * b;
} }
@ -55,7 +55,7 @@ macro_rules! bench_mul_dmat_dvec(
let m : DMat<f64> = DMat::new_random($nrows, $ncols); let m : DMat<f64> = DMat::new_random($nrows, $ncols);
let mut v : DVec<f64> = DVec::new_random($ncols); let mut v : DVec<f64> = DVec::new_random($ncols);
for _ in (0usize .. 1000) { for _ in 0usize .. 1000 {
// XXX: the clone here is highly undesirable! // XXX: the clone here is highly undesirable!
v = m.clone() * v v = m.clone() * v
} }

2
src/lib.rs
View File

@ -135,7 +135,7 @@ pub use traits::{
pub use structs::{ pub use structs::{
Identity, Identity,
DMat, DMat, DMat1, DMat2, DMat3, DMat4, DMat5, DMat6,
DVec, DVec1, DVec2, DVec3, DVec4, DVec5, DVec6, DVec, DVec1, DVec2, DVec3, DVec4, DVec5, DVec6,
Iso2, Iso3, Iso4, Iso2, Iso3, Iso4,
Mat1, Mat2, Mat3, Mat4, Mat1, Mat2, Mat3, Mat4,

837
src/structs/dmat.rs
View File

@ -3,14 +3,15 @@
#![allow(missing_docs)] // we hide doc to not have to document the $trhs double dispatch trait. #![allow(missing_docs)] // we hide doc to not have to document the $trhs double dispatch trait.
use std::cmp; use std::cmp;
use std::mem;
use std::iter::repeat; use std::iter::repeat;
use std::ops::{Add, Sub, Mul, Div, Index, IndexMut}; use std::ops::{Add, Sub, Mul, Div, Index, IndexMut};
use std::fmt::{Debug, Formatter, Result}; use std::fmt::{Debug, Formatter, Result};
use rand::{self, Rand}; use rand::{self, Rand};
use num::{Zero, One}; use num::{Zero, One};
use structs::dvec::DVec; use structs::dvec::{DVec, DVec1, DVec2, DVec3, DVec4, DVec5, DVec6};
use traits::operations::{ApproxEq, Inv, Transpose, Mean, Cov}; use traits::operations::{ApproxEq, Inv, Transpose, Mean, Cov};
use traits::structure::{Cast, Col, Row, ColSlice, RowSlice, Diag, DiagMut, Eye, Indexable, Shape, BaseNum}; use traits::structure::{Cast, Col, ColSlice, Row, RowSlice, Diag, DiagMut, Eye, Indexable, Shape, BaseNum};
#[cfg(feature="arbitrary")] #[cfg(feature="arbitrary")]
use quickcheck::{Arbitrary, Gen}; use quickcheck::{Arbitrary, Gen};
@ -38,47 +39,6 @@ impl<N> DMat<N> {
} }
} }
impl<N: Zero + Clone + Copy> DMat<N> {
/// Builds a matrix filled with zeros.
///
/// # Arguments
/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim`
/// components.
#[inline]
pub fn new_zeros(nrows: usize, ncols: usize) -> DMat<N> {
DMat::from_elem(nrows, ncols, ::zero())
}
/// Tests if all components of the matrix are zeroes.
#[inline]
pub fn is_zero(&self) -> bool {
self.mij.iter().all(|e| e.is_zero())
}
#[inline]
pub fn reset(&mut self) {
for mij in self.mij.iter_mut() {
*mij = ::zero();
}
}
}
impl<N: Rand> DMat<N> {
/// Builds a matrix filled with random values.
#[inline]
pub fn new_random(nrows: usize, ncols: usize) -> DMat<N> {
DMat::from_fn(nrows, ncols, |_, _| rand::random())
}
}
impl<N: One + Clone + Copy> DMat<N> {
/// Builds a matrix filled with a given constant.
#[inline]
pub fn new_ones(nrows: usize, ncols: usize) -> DMat<N> {
DMat::from_elem(nrows, ncols, ::one())
}
}
impl<N: Clone + Copy> DMat<N> { impl<N: Clone + Copy> DMat<N> {
/// Builds a matrix filled with a given constant. /// Builds a matrix filled with a given constant.
#[inline] #[inline]
@ -92,13 +52,35 @@ impl<N: Clone + Copy> DMat<N> {
/// Builds a matrix filled with the components provided by a vector. /// Builds a matrix filled with the components provided by a vector.
/// The vector contains the matrix data in row-major order. /// The vector contains the matrix data in row-major order.
/// Note that `from_col_vec` is a lot faster than `from_row_vec` since a `DMat` stores its data /// Note that `from_col_vec` is much faster than `from_row_vec` since a `DMat` stores its data
/// in column-major order. /// in column-major order.
/// ///
/// The vector must have at least `nrows * ncols` elements. /// The vector must have exactly `nrows * ncols` elements.
#[inline] #[inline]
pub fn from_row_vec(nrows: usize, ncols: usize, vec: &[N]) -> DMat<N> { pub fn from_row_vec(nrows: usize, ncols: usize, vec: &[N]) -> DMat<N> {
let mut res = DMat::from_col_vec(ncols, nrows, vec); DMat::from_row_iter(nrows, ncols, vec.to_vec())
}
/// Builds a matrix filled with the components provided by a vector.
/// The vector contains the matrix data in column-major order.
/// Note that `from_col_vec` is much faster than `from_row_vec` since a `DMat` stores its data
/// in column-major order.
///
/// The vector must have exactly `nrows * ncols` elements.
#[inline]
pub fn from_col_vec(nrows: usize, ncols: usize, vec: &[N]) -> DMat<N> {
DMat::from_col_iter(nrows, ncols, vec.to_vec())
}
/// Builds a matrix filled with the components provided by a source that may be moved into an iterator.
/// The source contains the matrix data in row-major order.
/// Note that `from_col_iter` is much faster than `from_row_iter` since a `DMat` stores its data
/// in column-major order.
///
/// The source must have exactly `nrows * ncols` elements.
#[inline]
pub fn from_row_iter<I: IntoIterator<Item = N>>(nrows: usize, ncols: usize, param: I) -> DMat<N> {
let mut res = DMat::from_col_iter(ncols, nrows, param);
// we transpose because the buffer is row_major // we transpose because the buffer is row_major
res.transpose_mut(); res.transpose_mut();
@ -106,26 +88,29 @@ impl<N: Clone + Copy> DMat<N> {
res res
} }
/// Builds a matrix filled with the components provided by a vector.
/// The vector contains the matrix data in column-major order. /// Builds a matrix filled with the components provided by a source that may be moved into an iterator.
/// Note that `from_col_vec` is a lot faster than `from_row_vec` since a `DMat` stores its data /// The source contains the matrix data in column-major order.
/// Note that `from_col_iter` is much faster than `from_row_iter` since a `DMat` stores its data
/// in column-major order. /// in column-major order.
/// ///
/// The vector must have at least `nrows * ncols` elements. /// The source must have exactly `nrows * ncols` elements.
#[inline] #[inline]
pub fn from_col_vec(nrows: usize, ncols: usize, vec: &[N]) -> DMat<N> { pub fn from_col_iter<I: IntoIterator<Item = N>>(nrows: usize, ncols: usize, param: I) -> DMat<N> {
assert!(nrows * ncols == vec.len()); let mij: Vec<N> = param.into_iter().collect();
assert!(nrows * ncols == mij.len(), "The ammount of data provided does not matches the matrix size.");
DMat { DMat {
nrows: nrows, nrows: nrows,
ncols: ncols, ncols: ncols,
mij: vec.to_vec() mij: mij
} }
} }
} }
impl<N> DMat<N> { impl<N> DMat<N> {
/// Builds a matrix filled with a given constant. /// Builds a matrix using an initialization function.
#[inline(always)] #[inline(always)]
pub fn from_fn<F: FnMut(usize, usize) -> N>(nrows: usize, ncols: usize, mut f: F) -> DMat<N> { pub fn from_fn<F: FnMut(usize, usize) -> N>(nrows: usize, ncols: usize, mut f: F) -> DMat<N> {
DMat { DMat {
@ -135,708 +120,102 @@ impl<N> DMat<N> {
} }
} }
/// Transforms this matrix isizeo an array. This consumes the matrix and is O(1). /// Transforms this matrix into an array. This consumes the matrix and is O(1).
/// The returned vector contains the matrix data in column-major order. /// The returned vector contains the matrix data in column-major order.
#[inline] #[inline]
pub fn to_vec(self) -> Vec<N> { pub fn into_vec(self) -> Vec<N> {
self.mij self.mij
} }
/// Gets a reference to this matrix data.
/// The returned vector contains the matrix data in column-major order.
#[inline]
pub fn as_vec(&self) -> &[N] {
&self.mij
}
/// Gets a mutable reference to this matrix data.
/// The returned vector contains the matrix data in column-major order.
#[inline]
pub fn as_mut_vec(&mut self) -> &mut [N] {
&mut self.mij[..]
}
}
// FIXME: add a function to modify the dimension (to avoid useless allocations)?
impl<N: One + Zero + Clone + Copy> Eye for DMat<N> {
/// Builds an identity matrix.
///
/// # Arguments
/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim`
/// components.
#[inline]
fn new_identity(dim: usize) -> DMat<N> {
let mut res = DMat::new_zeros(dim, dim);
for i in 0..dim {
let _1: N = ::one();
res[(i, i)] = _1;
}
res
}
}
impl<N> DMat<N> {
#[inline(always)]
fn offset(&self, i: usize, j: usize) -> usize {
i + j * self.nrows
}
}
impl<N: Copy> Indexable<(usize, usize), N> for DMat<N> {
/// Just like `set` without bounds checking.
#[inline]
unsafe fn unsafe_set(&mut self, rowcol: (usize, usize), val: N) {
let (row, col) = rowcol;
let offset = self.offset(row, col);
*self.mij[..].get_unchecked_mut(offset) = val
}
/// Just like `at` without bounds checking.
#[inline]
unsafe fn unsafe_at(&self, rowcol: (usize, usize)) -> N {
let (row, col) = rowcol;
*self.mij.get_unchecked(self.offset(row, col))
}
#[inline]
fn swap(&mut self, rowcol1: (usize, usize), rowcol2: (usize, usize)) {
let (row1, col1) = rowcol1;
let (row2, col2) = rowcol2;
let offset1 = self.offset(row1, col1);
let offset2 = self.offset(row2, col2);
let count = self.mij.len();
assert!(offset1 < count);
assert!(offset2 < count);
self.mij[..].swap(offset1, offset2);
}
}
impl<N: Copy> Col<DVec<N>> for DMat<N> {
/// The number of columns on the matrix.
#[inline]
fn ncols(&self) -> usize {
self.ncols
}
/// The `i`-th column of matrix.
#[inline]
fn col(&self, col_id: usize) -> DVec<N> {
assert!(col_id < self.ncols);
let start = self.offset(0, col_id);
let stop = self.offset(self.nrows, col_id);
DVec::from_slice(self.nrows, &self.mij[start .. stop])
}
/// Writes the `i`-th column of matrix.
#[inline]
fn set_col(&mut self, col_id: usize, col: DVec<N>) {
assert!(col_id < self.ncols);
for row_id in 0..self.nrows {
unsafe {
self.unsafe_set((row_id, col_id), col.unsafe_at(row_id));
}
}
}
}
impl<N: Copy> Row<DVec<N>> for DMat<N> {
/// The number of row on the matrix.
#[inline]
fn nrows(&self) -> usize {
self.nrows
}
/// The `i`-th row of matrix.
#[inline]
fn row(&self, row_id: usize) -> DVec<N> {
assert!(row_id < self.nrows);
let mut slice : DVec<N> = unsafe {
DVec::new_uninitialized(self.ncols)
};
for col_id in 0..self.ncols {
unsafe {
slice.unsafe_set(col_id, self.unsafe_at((row_id, col_id)));
}
}
slice
}
/// Writes the `i`-th row of matrix.
#[inline]
fn set_row(&mut self, row_id: usize, row: DVec<N>) {
assert!(row_id < self.nrows);
for col_id in 0..self.ncols {
unsafe {
self.unsafe_set((row_id, col_id), row.unsafe_at(col_id));
}
}
}
}
impl<N> Shape<(usize, usize)> for DMat<N> {
#[inline]
fn shape(&self) -> (usize, usize) {
(self.nrows, self.ncols)
}
}
impl<N> Index<(usize, usize)> for DMat<N> {
type Output = N;
fn index(&self, (i, j): (usize, usize)) -> &N {
assert!(i < self.nrows);
assert!(j < self.ncols);
unsafe {
self.mij.get_unchecked(self.offset(i, j))
}
}
}
impl<N> IndexMut<(usize, usize)> for DMat<N> {
fn index_mut(&mut self, (i, j): (usize, usize)) -> &mut N {
assert!(i < self.nrows);
assert!(j < self.ncols);
let offset = self.offset(i, j);
unsafe {
self.mij[..].get_unchecked_mut(offset)
}
}
}
impl<N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<DMat<N>> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn mul(self, right: DMat<N>) -> DMat<N> {
(&self) * (&right)
}
}
impl<'a, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<&'a DMat<N>> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn mul(self, right: &'a DMat<N>) -> DMat<N> {
(&self) * right
}
}
impl<'a, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<DMat<N>> for &'a DMat<N> {
type Output = DMat<N>;
#[inline]
fn mul(self, right: DMat<N>) -> DMat<N> {
self * (&right)
}
}
impl<'a, 'b, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<&'b DMat<N>> for &'a DMat<N> {
type Output = DMat<N>;
#[inline]
fn mul(self, right: &DMat<N>) -> DMat<N> {
assert!(self.ncols == right.nrows);
let mut res = unsafe { DMat::new_uninitialized(self.nrows, right.ncols) };
for i in 0..self.nrows {
for j in 0..right.ncols {
let mut acc: N = ::zero();
unsafe {
for k in 0..self.ncols {
acc = acc
+ self.unsafe_at((i, k)) * right.unsafe_at((k, j));
}
res.unsafe_set((i, j), acc);
}
}
}
res
}
}
impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero> Mul<DVec<N>> for DMat<N> {
type Output = DVec<N>;
fn mul(self, right: DVec<N>) -> DVec<N> {
assert!(self.ncols == right.at.len());
let mut res : DVec<N> = unsafe { DVec::new_uninitialized(self.nrows) };
for i in 0..self.nrows {
let mut acc: N = ::zero();
for j in 0..self.ncols {
unsafe {
acc = acc + self.unsafe_at((i, j)) * right.unsafe_at(j);
}
}
res.at[i] = acc;
}
res
}
}
impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero> Mul<DMat<N>> for DVec<N> {
type Output = DVec<N>;
fn mul(self, right: DMat<N>) -> DVec<N> {
assert!(right.nrows == self.at.len());
let mut res : DVec<N> = unsafe { DVec::new_uninitialized(right.ncols) };
for i in 0..right.ncols {
let mut acc: N = ::zero();
for j in 0..right.nrows {
unsafe {
acc = acc + self.unsafe_at(j) * right.unsafe_at((j, i));
}
}
res.at[i] = acc;
}
res
}
}
impl<N: BaseNum + Clone> Inv for DMat<N> {
#[inline]
fn inv(&self) -> Option<DMat<N>> {
let mut res: DMat<N> = self.clone();
if res.inv_mut() {
Some(res)
}
else {
None
}
}
fn inv_mut(&mut self) -> bool {
assert!(self.nrows == self.ncols);
let dim = self.nrows;
let mut res: DMat<N> = Eye::new_identity(dim);
// inversion using Gauss-Jordan elimination
for k in 0..dim {
// search a non-zero value on the k-th column
// FIXME: would it be worth it to spend some more time searching for the
// max instead?
let mut n0 = k; // index of a non-zero entry
while n0 != dim {
if unsafe { self.unsafe_at((n0, k)) } != ::zero() {
break;
}
n0 = n0 + 1;
}
if n0 == dim {
return false
}
// swap pivot line
if n0 != k {
for j in 0..dim {
let off_n0_j = self.offset(n0, j);
let off_k_j = self.offset(k, j);
self.mij[..].swap(off_n0_j, off_k_j);
res.mij[..].swap(off_n0_j, off_k_j);
}
} }
unsafe { dmat_impl!(DMat, DVec);
let pivot = self.unsafe_at((k, k));
for j in k..dim {
let selfval = self.unsafe_at((k, j)) / pivot;
self.unsafe_set((k, j), selfval);
}
for j in 0..dim {
let resval = res.unsafe_at((k, j)) / pivot;
res.unsafe_set((k, j), resval);
}
for l in 0..dim {
if l != k {
let normalizer = self.unsafe_at((l, k));
for j in k..dim {
let selfval = self.unsafe_at((l, j)) - self.unsafe_at((k, j)) * normalizer;
self.unsafe_set((l, j), selfval);
}
for j in 0..dim {
let resval = res.unsafe_at((l, j)) - res.unsafe_at((k, j)) * normalizer;
res.unsafe_set((l, j), resval);
}
}
}
}
}
*self = res;
true
}
}
impl<N: Clone + Copy> Transpose for DMat<N> {
#[inline]
fn transpose(&self) -> DMat<N> {
if self.nrows == self.ncols {
let mut res = self.clone();
res.transpose_mut();
res
}
else {
let mut res = unsafe { DMat::new_uninitialized(self.ncols, self.nrows) };
for i in 0..self.nrows {
for j in 0..self.ncols {
unsafe {
res.unsafe_set((j, i), self.unsafe_at((i, j)))
}
}
}
res
}
}
#[inline]
fn transpose_mut(&mut self) {
if self.nrows == self.ncols {
let n = self.nrows;
for i in 0..n - 1 {
for j in i + 1..n {
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);
}
}
}
else {
// FIXME: implement a better algorithm which does that in-place.
*self = Transpose::transpose(self);
}
}
}
impl<N: BaseNum + Cast<f64> + 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 = Cast::from(1.0f64 / self.nrows as f64);
for i in 0..self.nrows { pub struct DMat1<N> {
for j in 0..self.ncols { nrows: usize,
unsafe { ncols: usize,
let acc = res.unsafe_at(j) + self.unsafe_at((i, j)) * normalizer; mij: [N; 1 * 1],
res.unsafe_set(j, acc);
}
}
} }
res small_dmat_impl!(DMat1, DVec1, 1, 0);
} small_dmat_from_impl!(DMat1, 1, ::zero());
}
impl<N: BaseNum + Cast<f64> + Clone> 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) }; pub struct DMat2<N> {
let mean = self.mean(); nrows: usize,
ncols: usize,
// FIXME: use the rows iterator when available mij: [N; 2 * 2],
for i in 0..self.nrows {
for j in 0..self.ncols {
unsafe {
centered.unsafe_set((i, j), self.unsafe_at((i, j)) - mean.unsafe_at(j));
}
}
} }
// FIXME: return a triangular matrix? small_dmat_impl!(DMat2, DVec2, 2, 0, 1,
let fnormalizer: f64 = Cast::from(self.nrows() - 1); 2, 3);
let normalizer: N = Cast::from(fnormalizer); small_dmat_from_impl!(DMat2, 2, ::zero(), ::zero(),
// FIXME: this will do 2 allocations for temporaries! ::zero(), ::zero());
(Transpose::transpose(&centered) * centered) / normalizer
}
}
impl<N: Copy + Clone> ColSlice<DVec<N>> for DMat<N> {
fn col_slice(&self, col_id :usize, row_start: usize, row_end: usize) -> DVec<N> {
assert!(col_id < self.ncols);
assert!(row_start < row_end);
assert!(row_end <= self.nrows);
// 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_slice(row_end - row_start, &self.mij[start .. stop]);
slice
}
}
impl<N: Copy> RowSlice<DVec<N>> for DMat<N> { pub struct DMat3<N> {
fn row_slice(&self, row_id :usize, col_start: usize, col_end: usize) -> DVec<N> { nrows: usize,
assert!(row_id < self.nrows); ncols: usize,
assert!(col_start < col_end); mij: [N; 3 * 3],
assert!(col_end <= self.ncols);
let mut slice : DVec<N> = unsafe {
DVec::new_uninitialized(col_end - col_start)
};
let mut slice_idx = 0;
for col_id in col_start..col_end {
unsafe {
slice.unsafe_set(slice_idx, self.unsafe_at((row_id, col_id)));
}
slice_idx += 1;
}
slice
}
} }
impl<N: Copy + Clone + Zero> Diag<DVec<N>> for DMat<N> { small_dmat_impl!(DMat3, DVec3, 3, 0, 1, 2,
#[inline] 3, 4, 5,
fn from_diag(diag: &DVec<N>) -> DMat<N> { 6, 7, 8);
let mut res = DMat::new_zeros(diag.len(), diag.len()); small_dmat_from_impl!(DMat3, 3, ::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero());
res.set_diag(diag);
res pub struct DMat4<N> {
nrows: usize,
ncols: usize,
mij: [N; 4 * 4],
} }
#[inline] small_dmat_impl!(DMat4, DVec4, 4, 0, 1, 2, 3,
fn diag(&self) -> DVec<N> { 4, 5, 6, 7,
let smallest_dim = cmp::min(self.nrows, self.ncols); 8, 9, 10, 11,
12, 13, 14, 15);
small_dmat_from_impl!(DMat4, 4, ::zero(), ::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero(), ::zero());
let mut diag: DVec<N> = DVec::new_zeros(smallest_dim);
for i in 0..smallest_dim { pub struct DMat5<N> {
unsafe { diag.unsafe_set(i, self.unsafe_at((i, i))) } nrows: usize,
ncols: usize,
mij: [N; 5 * 5],
} }
diag small_dmat_impl!(DMat5, DVec5, 5, 0, 1, 2, 3, 4,
} 5, 6, 7, 8, 9,
} 10, 11, 12, 13, 14,
15, 16, 17, 18, 19,
impl<N: Copy + Clone + Zero> DiagMut<DVec<N>> for DMat<N> { 20, 21, 22, 23, 24);
#[inline] small_dmat_from_impl!(DMat5, 5, ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
fn set_diag(&mut self, diag: &DVec<N>) { ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
let smallest_dim = cmp::min(self.nrows, self.ncols); ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
assert!(diag.len() == smallest_dim); ::zero(), ::zero(), ::zero(), ::zero(), ::zero());
for i in 0..smallest_dim {
unsafe { self.unsafe_set((i, i), diag.unsafe_at(i)) }
}
}
}
impl<N: ApproxEq<N>> ApproxEq<N> for DMat<N> {
#[inline]
fn approx_epsilon(_: Option<DMat<N>>) -> N {
ApproxEq::approx_epsilon(None::<N>)
}
#[inline]
fn approx_ulps(_: Option<DMat<N>>) -> u32 {
ApproxEq::approx_ulps(None::<N>)
}
#[inline]
fn approx_eq_eps(&self, other: &DMat<N>, epsilon: &N) -> bool {
let mut zip = self.mij.iter().zip(other.mij.iter());
zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
}
#[inline] pub struct DMat6<N> {
fn approx_eq_ulps(&self, other: &DMat<N>, ulps: u32) -> bool { nrows: usize,
let mut zip = self.mij.iter().zip(other.mij.iter()); ncols: usize,
zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps)) mij: [N; 6 * 6],
}
} }
impl<N: Debug + Copy> Debug for DMat<N> { small_dmat_impl!(DMat6, DVec6, 6, 0, 1, 2, 3, 4, 5,
fn fmt(&self, form:&mut Formatter) -> Result { 6, 7, 8, 9, 10, 11,
for i in 0..self.nrows() { 12, 13, 14, 15, 16, 17,
for j in 0..self.ncols() { 18, 19, 20, 21, 22, 23,
let _ = write!(form, "{:?} ", self[(i, j)]); 24, 25, 26, 27, 28, 29,
} 30, 31, 32, 33, 34, 35);
let _ = write!(form, "\n"); small_dmat_from_impl!(DMat6, 6, ::zero(), ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
} ::zero(), ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
write!(form, "\n") ::zero(), ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
} ::zero(), ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
} ::zero(), ::zero(), ::zero(), ::zero(), ::zero(), ::zero(),
::zero(), ::zero(), ::zero(), ::zero(), ::zero(), ::zero());
impl<N: Copy + Mul<N, Output = N>> Mul<N> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn mul(self, right: N) -> DMat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij * right;
}
res
}
}
impl<N: Copy + Div<N, Output = N>> Div<N> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn div(self, right: N) -> DMat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij / right;
}
res
}
}
impl<N: Copy + Add<N, Output = N>> Add<N> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn add(self, right: N) -> DMat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij + right;
}
res
}
}
impl<N: Copy + Add<N, Output = N>> Add<DMat<N>> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn add(self, right: DMat<N>) -> DMat<N> {
self + (&right)
}
}
impl<'a, N: Copy + Add<N, Output = N>> Add<DMat<N>> for &'a DMat<N> {
type Output = DMat<N>;
#[inline]
fn add(self, right: DMat<N>) -> DMat<N> {
right + self
}
}
impl<'a, N: Copy + Add<N, Output = N>> Add<&'a DMat<N>> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn add(self, right: &'a DMat<N>) -> DMat<N> {
assert!(self.nrows == right.nrows && self.ncols == right.ncols,
"Unable to add matrices with different dimensions.");
let mut res = self;
for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
*mij = *mij + *right_ij;
}
res
}
}
impl<N: Copy + Sub<N, Output = N>> Sub<N> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn sub(self, right: N) -> DMat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij - right;
}
res
}
}
impl<N: Copy + Sub<N, Output = N>> Sub<DMat<N>> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn sub(self, right: DMat<N>) -> DMat<N> {
self - (&right)
}
}
impl<'a, N: Copy + Sub<N, Output = N>> Sub<DMat<N>> for &'a DMat<N> {
type Output = DMat<N>;
#[inline]
fn sub(self, right: DMat<N>) -> DMat<N> {
right - self
}
}
impl<'a, N: Copy + Sub<N, Output = N>> Sub<&'a DMat<N>> for DMat<N> {
type Output = DMat<N>;
#[inline]
fn sub(self, right: &'a DMat<N>) -> DMat<N> {
assert!(self.nrows == right.nrows && self.ncols == right.ncols,
"Unable to subtract matrices with different dimensions.");
let mut res = self;
for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
*mij = *mij - *right_ij;
}
res
}
}
#[cfg(feature="arbitrary")]
impl<N: Arbitrary> Arbitrary for DMat<N> {
fn arbitrary<G: Gen>(g: &mut G) -> DMat<N> {
DMat::from_fn(
Arbitrary::arbitrary(g), Arbitrary::arbitrary(g),
|_, _| Arbitrary::arbitrary(g)
)
}
}

900
src/structs/dmat_macros.rs Normal file
View File

@ -0,0 +1,900 @@
#![macro_use]
macro_rules! dmat_impl(
($dmat: ident, $dvec: ident) => (
impl<N: Zero + Clone + Copy> $dmat<N> {
/// Builds a matrix filled with zeros.
///
/// # Arguments
/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim`
/// components.
#[inline]
pub fn new_zeros(nrows: usize, ncols: usize) -> $dmat<N> {
$dmat::from_elem(nrows, ncols, ::zero())
}
/// Tests if all components of the matrix are zeroes.
#[inline]
pub fn is_zero(&self) -> bool {
self.mij.iter().all(|e| e.is_zero())
}
#[inline]
pub fn reset(&mut self) {
for mij in self.mij.iter_mut() {
*mij = ::zero();
}
}
}
impl<N: Zero + Copy + Rand> $dmat<N> {
/// Builds a matrix filled with random values.
#[inline]
pub fn new_random(nrows: usize, ncols: usize) -> $dmat<N> {
$dmat::from_fn(nrows, ncols, |_, _| rand::random())
}
}
impl<N: One + Zero + Clone + Copy> $dmat<N> {
/// Builds a matrix filled with a given constant.
#[inline]
pub fn new_ones(nrows: usize, ncols: usize) -> $dmat<N> {
$dmat::from_elem(nrows, ncols, ::one())
}
}
impl<N> $dmat<N> {
/// The number of row on the matrix.
#[inline]
pub fn nrows(&self) -> usize {
self.nrows
}
/// The number of columns on the matrix.
#[inline]
pub fn ncols(&self) -> usize {
self.ncols
}
/// Gets a reference to this matrix data.
/// The returned vector contains the matrix data in column-major order.
#[inline]
pub fn as_vec(&self) -> &[N] {
&self.mij
}
/// Gets a mutable reference to this matrix data.
/// The returned vector contains the matrix data in column-major order.
#[inline]
pub fn as_mut_vec(&mut self) -> &mut [N] {
&mut self.mij[..]
}
}
// FIXME: add a function to modify the dimension (to avoid useless allocations)?
impl<N: One + Zero + Clone + Copy> Eye for $dmat<N> {
/// Builds an identity matrix.
///
/// # Arguments
/// * `dim` - The dimension of the matrix. A `dim`-dimensional matrix contains `dim * dim`
/// components.
#[inline]
fn new_identity(dim: usize) -> $dmat<N> {
let mut res = $dmat::new_zeros(dim, dim);
for i in 0..dim {
let _1: N = ::one();
res[(i, i)] = _1;
}
res
}
}
impl<N> $dmat<N> {
#[inline(always)]
fn offset(&self, i: usize, j: usize) -> usize {
i + j * self.nrows
}
}
impl<N: Copy> Indexable<(usize, usize), N> for $dmat<N> {
/// Just like `set` without bounds checking.
#[inline]
unsafe fn unsafe_set(&mut self, rowcol: (usize, usize), val: N) {
let (row, col) = rowcol;
let offset = self.offset(row, col);
*self.mij[..].get_unchecked_mut(offset) = val
}
/// Just like `at` without bounds checking.
#[inline]
unsafe fn unsafe_at(&self, rowcol: (usize, usize)) -> N {
let (row, col) = rowcol;
*self.mij.get_unchecked(self.offset(row, col))
}
#[inline]
fn swap(&mut self, rowcol1: (usize, usize), rowcol2: (usize, usize)) {
let (row1, col1) = rowcol1;
let (row2, col2) = rowcol2;
let offset1 = self.offset(row1, col1);
let offset2 = self.offset(row2, col2);
let count = self.mij.len();
assert!(offset1 < count);
assert!(offset2 < count);
self.mij[..].swap(offset1, offset2);
}
}
impl<N> Shape<(usize, usize)> for $dmat<N> {
#[inline]
fn shape(&self) -> (usize, usize) {
(self.nrows, self.ncols)
}
}
impl<N> Index<(usize, usize)> for $dmat<N> {
type Output = N;
fn index(&self, (i, j): (usize, usize)) -> &N {
assert!(i < self.nrows);
assert!(j < self.ncols);
unsafe {
self.mij.get_unchecked(self.offset(i, j))
}
}
}
impl<N> IndexMut<(usize, usize)> for $dmat<N> {
fn index_mut(&mut self, (i, j): (usize, usize)) -> &mut N {
assert!(i < self.nrows);
assert!(j < self.ncols);
let offset = self.offset(i, j);
unsafe {
self.mij[..].get_unchecked_mut(offset)
}
}
}
impl<N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<$dmat<N>> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn mul(self, right: $dmat<N>) -> $dmat<N> {
(&self) * (&right)
}
}
impl<'a, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<&'a $dmat<N>> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn mul(self, right: &'a $dmat<N>) -> $dmat<N> {
(&self) * right
}
}
impl<'a, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<$dmat<N>> for &'a $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn mul(self, right: $dmat<N>) -> $dmat<N> {
self * (&right)
}
}
impl<'a, 'b, N: Copy + Mul<N, Output = N> + Add<N, Output = N> + Zero> Mul<&'b $dmat<N>> for &'a $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn mul(self, right: &$dmat<N>) -> $dmat<N> {
assert!(self.ncols == right.nrows);
let mut res = unsafe { $dmat::new_uninitialized(self.nrows, right.ncols) };
for i in 0..self.nrows {
for j in 0..right.ncols {
let mut acc: N = ::zero();
unsafe {
for k in 0..self.ncols {
acc = acc
+ self.unsafe_at((i, k)) * right.unsafe_at((k, j));
}
res.unsafe_set((i, j), acc);
}
}
}
res
}
}
impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero> Mul<$dvec<N>> for $dmat<N> {
type Output = $dvec<N>;
fn mul(self, right: $dvec<N>) -> $dvec<N> {
assert!(self.ncols == right.len());
let mut res : $dvec<N> = unsafe { $dvec::new_uninitialized(self.nrows) };
for i in 0..self.nrows {
let mut acc: N = ::zero();
for j in 0..self.ncols {
unsafe {
acc = acc + self.unsafe_at((i, j)) * right.unsafe_at(j);
}
}
unsafe {
res.unsafe_set(i, acc);
}
}
res
}
}
impl<N: Copy + Add<N, Output = N> + Mul<N, Output = N> + Zero> Mul<$dmat<N>> for $dvec<N> {
type Output = $dvec<N>;
fn mul(self, right: $dmat<N>) -> $dvec<N> {
assert!(right.nrows == self.len());
let mut res : $dvec<N> = unsafe { $dvec::new_uninitialized(right.ncols) };
for i in 0..right.ncols {
let mut acc: N = ::zero();
for j in 0..right.nrows {
unsafe {
acc = acc + self.unsafe_at(j) * right.unsafe_at((j, i));
}
}
unsafe {
res.unsafe_set(i, acc);
}
}
res
}
}
impl<N: BaseNum + Clone> Inv for $dmat<N> {
#[inline]
fn inv(&self) -> Option<$dmat<N>> {
let mut res: $dmat<N> = self.clone();
if res.inv_mut() {
Some(res)
}
else {
None
}
}
fn inv_mut(&mut self) -> bool {
assert!(self.nrows == self.ncols);
let dim = self.nrows;
let mut res: $dmat<N> = Eye::new_identity(dim);
// inversion using Gauss-Jordan elimination
for k in 0..dim {
// search a non-zero value on the k-th column
// FIXME: would it be worth it to spend some more time searching for the
// max instead?
let mut n0 = k; // index of a non-zero entry
while n0 != dim {
if unsafe { self.unsafe_at((n0, k)) } != ::zero() {
break;
}
n0 = n0 + 1;
}
if n0 == dim {
return false
}
// swap pivot line
if n0 != k {
for j in 0..dim {
let off_n0_j = self.offset(n0, j);
let off_k_j = self.offset(k, j);
self.mij[..].swap(off_n0_j, off_k_j);
res.mij[..].swap(off_n0_j, off_k_j);
}
}
unsafe {
let pivot = self.unsafe_at((k, k));
for j in k..dim {
let selfval = self.unsafe_at((k, j)) / pivot;
self.unsafe_set((k, j), selfval);
}
for j in 0..dim {
let resval = res.unsafe_at((k, j)) / pivot;
res.unsafe_set((k, j), resval);
}
for l in 0..dim {
if l != k {
let normalizer = self.unsafe_at((l, k));
for j in k..dim {
let selfval = self.unsafe_at((l, j)) - self.unsafe_at((k, j)) * normalizer;
self.unsafe_set((l, j), selfval);
}
for j in 0..dim {
let resval = res.unsafe_at((l, j)) - res.unsafe_at((k, j)) * normalizer;
res.unsafe_set((l, j), resval);
}
}
}
}
}
*self = res;
true
}
}
impl<N: Clone + Copy> Transpose for $dmat<N> {
#[inline]
fn transpose(&self) -> $dmat<N> {
if self.nrows == self.ncols {
let mut res = self.clone();
res.transpose_mut();
res
}
else {
let mut res = unsafe { $dmat::new_uninitialized(self.ncols, self.nrows) };
for i in 0..self.nrows {
for j in 0..self.ncols {
unsafe {
res.unsafe_set((j, i), self.unsafe_at((i, j)))
}
}
}
res
}
}
#[inline]
fn transpose_mut(&mut self) {
if self.nrows == self.ncols {
let n = self.nrows;
for i in 0..n - 1 {
for j in i + 1..n {
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);
}
}
}
else {
// FIXME: implement a better algorithm which does that in-place.
*self = Transpose::transpose(self);
}
}
}
impl<N: BaseNum + Cast<f64> + 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 = Cast::from(1.0f64 / self.nrows as f64);
for i in 0 .. self.nrows {
for j in 0 .. self.ncols {
unsafe {
let acc = res.unsafe_at(j) + self.unsafe_at((i, j)) * normalizer;
res.unsafe_set(j, acc);
}
}
}
res
}
}
impl<N: BaseNum + Cast<f64> + Clone> 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 0 .. self.nrows {
for j in 0 .. self.ncols {
unsafe {
centered.unsafe_set((i, j), self.unsafe_at((i, j)) - mean.unsafe_at(j));
}
}
}
// FIXME: return a triangular matrix?
let fnormalizer: f64 = Cast::from(self.nrows() - 1);
let normalizer: N = Cast::from(fnormalizer);
// FIXME: this will do 2 allocations for temporaries!
(Transpose::transpose(&centered) * centered) / normalizer
}
}
impl<N: Copy + Zero> Col<$dvec<N>> for $dmat<N> {
#[inline]
fn ncols(&self) -> usize {
self.ncols
}
#[inline]
fn set_col(&mut self, col_id: usize, v: $dvec<N>) {
assert!(col_id < self.ncols);
assert!(self.nrows == v.len());
for (i, e) in v[..].iter().enumerate() {
unsafe {
self.unsafe_set((i, col_id), *e);
}
}
}
#[inline]
fn col(&self, col_id: usize) -> $dvec<N> {
let mut res: $dvec<N> = unsafe {
$dvec::new_uninitialized(self.nrows)
};
for (row_id, e) in res[..].iter_mut().enumerate() {
*e = unsafe { self.unsafe_at((row_id, col_id)) };
}
res
}
}
impl<N: Copy + Clone + Zero> ColSlice<$dvec<N>> for $dmat<N> {
fn col_slice(&self, col_id :usize, row_start: usize, row_end: usize) -> $dvec<N> {
assert!(col_id < self.ncols);
assert!(row_start < row_end);
assert!(row_end <= self.nrows);
// 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_slice(row_end - row_start, &self.mij[start .. stop]);
slice
}
}
impl<N: Copy + Zero> Row<$dvec<N>> for $dmat<N> {
#[inline]
fn nrows(&self) -> usize {
self.nrows
}
#[inline]
fn set_row(&mut self, row_id: usize, v: $dvec<N>) {
assert!(row_id < self.nrows);
assert!(self.ncols == v.len());
for (i, e) in v[..].iter().enumerate() {
unsafe {
self.unsafe_set((row_id, i), *e);
}
}
}
#[inline]
fn row(&self, row_id: usize) -> $dvec<N> {
let mut res: $dvec<N> = unsafe {
$dvec::new_uninitialized(self.ncols)
};
for (col_id, e) in res[..].iter_mut().enumerate() {
*e = unsafe { self.unsafe_at((row_id, col_id)) };
}
res
}
}
impl<N: Copy> RowSlice<$dvec<N>> for $dmat<N> {
fn row_slice(&self, row_id :usize, col_start: usize, col_end: usize) -> $dvec<N> {
assert!(row_id < self.nrows);
assert!(col_start < col_end);
assert!(col_end <= self.ncols);
let mut slice : $dvec<N> = unsafe {
$dvec::new_uninitialized(col_end - col_start)
};
let mut slice_idx = 0;
for col_id in col_start..col_end {
unsafe {
slice.unsafe_set(slice_idx, self.unsafe_at((row_id, col_id)));
}
slice_idx += 1;
}
slice
}
}
impl<N: Copy + Clone + Zero> Diag<$dvec<N>> for $dmat<N> {
#[inline]
fn from_diag(diag: &$dvec<N>) -> $dmat<N> {
let mut res = $dmat::new_zeros(diag.len(), diag.len());
res.set_diag(diag);
res
}
#[inline]
fn diag(&self) -> $dvec<N> {
let smallest_dim = cmp::min(self.nrows, self.ncols);
let mut diag: $dvec<N> = $dvec::new_zeros(smallest_dim);
for i in 0..smallest_dim {
unsafe { diag.unsafe_set(i, self.unsafe_at((i, i))) }
}
diag
}
}
impl<N: Copy + Clone + Zero> DiagMut<$dvec<N>> for $dmat<N> {
#[inline]
fn set_diag(&mut self, diag: &$dvec<N>) {
let smallest_dim = cmp::min(self.nrows, self.ncols);
assert!(diag.len() == smallest_dim);
for i in 0..smallest_dim {
unsafe { self.unsafe_set((i, i), diag.unsafe_at(i)) }
}
}
}
impl<N: ApproxEq<N>> ApproxEq<N> for $dmat<N> {
#[inline]
fn approx_epsilon(_: Option<$dmat<N>>) -> N {
ApproxEq::approx_epsilon(None::<N>)
}
#[inline]
fn approx_ulps(_: Option<$dmat<N>>) -> u32 {
ApproxEq::approx_ulps(None::<N>)
}
#[inline]
fn approx_eq_eps(&self, other: &$dmat<N>, epsilon: &N) -> bool {
let mut zip = self.mij.iter().zip(other.mij.iter());
zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
}
#[inline]
fn approx_eq_ulps(&self, other: &$dmat<N>, ulps: u32) -> bool {
let mut zip = self.mij.iter().zip(other.mij.iter());
zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
}
}
impl<N: Debug + Copy> Debug for $dmat<N> {
fn fmt(&self, form:&mut Formatter) -> Result {
for i in 0..self.nrows() {
for j in 0..self.ncols() {
let _ = write!(form, "{:?} ", self[(i, j)]);
}
let _ = write!(form, "\n");
}
write!(form, "\n")
}
}
impl<N: Copy + Mul<N, Output = N>> Mul<N> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn mul(self, right: N) -> $dmat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij * right;
}
res
}
}
impl<N: Copy + Div<N, Output = N>> Div<N> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn div(self, right: N) -> $dmat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij / right;
}
res
}
}
impl<N: Copy + Add<N, Output = N>> Add<N> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn add(self, right: N) -> $dmat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij + right;
}
res
}
}
impl<N: Copy + Add<N, Output = N>> Add<$dmat<N>> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn add(self, right: $dmat<N>) -> $dmat<N> {
self + (&right)
}
}
impl<'a, N: Copy + Add<N, Output = N>> Add<$dmat<N>> for &'a $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn add(self, right: $dmat<N>) -> $dmat<N> {
right + self
}
}
impl<'a, N: Copy + Add<N, Output = N>> Add<&'a $dmat<N>> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn add(self, right: &'a $dmat<N>) -> $dmat<N> {
assert!(self.nrows == right.nrows && self.ncols == right.ncols,
"Unable to add matrices with different dimensions.");
let mut res = self;
for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
*mij = *mij + *right_ij;
}
res
}
}
impl<N: Copy + Sub<N, Output = N>> Sub<N> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn sub(self, right: N) -> $dmat<N> {
let mut res = self;
for mij in res.mij.iter_mut() {
*mij = *mij - right;
}
res
}
}
impl<N: Copy + Sub<N, Output = N>> Sub<$dmat<N>> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn sub(self, right: $dmat<N>) -> $dmat<N> {
self - (&right)
}
}
impl<'a, N: Copy + Sub<N, Output = N>> Sub<$dmat<N>> for &'a $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn sub(self, right: $dmat<N>) -> $dmat<N> {
right - self
}
}
impl<'a, N: Copy + Sub<N, Output = N>> Sub<&'a $dmat<N>> for $dmat<N> {
type Output = $dmat<N>;
#[inline]
fn sub(self, right: &'a $dmat<N>) -> $dmat<N> {
assert!(self.nrows == right.nrows && self.ncols == right.ncols,
"Unable to subtract matrices with different dimensions.");
let mut res = self;
for (mij, right_ij) in res.mij.iter_mut().zip(right.mij.iter()) {
*mij = *mij - *right_ij;
}
res
}
}
#[cfg(feature="arbitrary")]
impl<N: Copy + Zero + Arbitrary> Arbitrary for $dmat<N> {
fn arbitrary<G: Gen>(g: &mut G) -> $dmat<N> {
$dmat::from_fn(
Arbitrary::arbitrary(g), Arbitrary::arbitrary(g),
|_, _| Arbitrary::arbitrary(g)
)
}
}
)
);
macro_rules! small_dmat_impl (
($dmat: ident, $dvec: ident, $dim: expr, $($idx: expr),*) => (
impl<N: PartialEq> PartialEq for $dmat<N> {
#[inline]
fn eq(&self, other: &$dmat<N>) -> bool {
if self.nrows() != other.nrows() || self.ncols() != other.ncols() {
return false; // FIXME: fail instead?
}
for (a, b) in self.mij[0 .. self.nrows() * self.ncols()].iter().zip(
other.mij[0 .. self.nrows() * self.ncols()].iter()) {
if *a != *b {
return false;
}
}
true
}
}
impl<N: Clone> Clone for $dmat<N> {
fn clone(&self) -> $dmat<N> {
let mij: [N; $dim * $dim] = [ $( self.mij[$idx].clone(), )* ];
$dmat {
nrows: self.nrows,
ncols: self.ncols,
mij: mij,
}
}
}
dmat_impl!($dmat, $dvec);
)
);
macro_rules! small_dmat_from_impl(
($dmat: ident, $dim: expr, $($zeros: expr),*) => (
impl<N: Zero + Clone + Copy> $dmat<N> {
/// Builds a matrix filled with a given constant.
#[inline]
pub fn from_elem(nrows: usize, ncols: usize, elem: N) -> $dmat<N> {
assert!(nrows <= $dim);
assert!(ncols <= $dim);
let mut mij: [N; $dim * $dim] = [ $( $zeros, )* ];
for n in &mut mij[.. nrows * ncols] {
*n = elem;
}
$dmat {
nrows: nrows,
ncols: ncols,
mij: mij
}
}
/// Builds a matrix filled with the components provided by a vector.
/// The vector contains the matrix data in row-major order.
/// Note that `from_col_vec` is a lot faster than `from_row_vec` since a `$dmat` stores its data
/// in column-major order.
///
/// The vector must have at least `nrows * ncols` elements.
#[inline]
pub fn from_row_vec(nrows: usize, ncols: usize, vec: &[N]) -> $dmat<N> {
let mut res = $dmat::from_col_vec(ncols, nrows, vec);
// we transpose because the buffer is row_major
res.transpose_mut();
res
}
/// Builds a matrix filled with the components provided by a vector.
/// The vector contains the matrix data in column-major order.
/// Note that `from_col_vec` is a lot faster than `from_row_vec` since a `$dmat` stores its data
/// in column-major order.
///
/// The vector must have at least `nrows * ncols` elements.
#[inline]
pub fn from_col_vec(nrows: usize, ncols: usize, vec: &[N]) -> $dmat<N> {
assert!(nrows * ncols == vec.len());
let mut mij: [N; $dim * $dim] = [ $( $zeros, )* ];
for (n, val) in mij[.. nrows * ncols].iter_mut().zip(vec.iter()) {
*n = *val;
}
$dmat {
nrows: nrows,
ncols: ncols,
mij: mij
}
}
/// Builds a matrix using an initialization function.
#[inline(always)]
pub fn from_fn<F: FnMut(usize, usize) -> N>(nrows: usize, ncols: usize, mut f: F) -> $dmat<N> {
assert!(nrows <= $dim);
assert!(ncols <= $dim);
let mut mij: [N; $dim * $dim] = [ $( $zeros, )* ];
for i in 0 .. nrows {
for j in 0 .. ncols {
mij[i + j * nrows] = f(i, j)
}
}
$dmat {
nrows: nrows,
ncols: ncols,
mij: mij
}
}
}
impl<N> $dmat<N> {
#[inline]
pub unsafe fn new_uninitialized(nrows: usize, ncols: usize) -> $dmat<N> {
assert!(nrows <= $dim);
assert!(ncols <= $dim);
$dmat {
nrows: nrows,
ncols: ncols,
mij: mem::uninitialized()
}
}
}
)
);

5
src/structs/dvec.rs
View File

@ -6,11 +6,12 @@ use std::slice::{Iter, IterMut};
use std::iter::{FromIterator, IntoIterator}; use std::iter::{FromIterator, IntoIterator};
use std::iter::repeat; use std::iter::repeat;
use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut}; use std::ops::{Add, Sub, Mul, Div, Neg, Index, IndexMut};
use std::mem;
use rand::{self, Rand}; use rand::{self, Rand};
use num::{Zero, One}; use num::{Zero, One};
use traits::operations::{ApproxEq, Axpy}; use traits::operations::{ApproxEq, Axpy, Mean};
use traits::geometry::{Dot, Norm}; use traits::geometry::{Dot, Norm};
use traits::structure::{Iterable, IterableMut, Indexable, Shape, BaseFloat, BaseNum}; use traits::structure::{Iterable, IterableMut, Indexable, Shape, BaseFloat, BaseNum, Cast};
#[cfg(feature="arbitrary")] #[cfg(feature="arbitrary")]
use quickcheck::{Arbitrary, Gen}; use quickcheck::{Arbitrary, Gen};

74
src/structs/dvec_macros.rs
View File

@ -15,20 +15,20 @@ macro_rules! dvec_impl(
/// Tests if all components of the vector are zeroes. /// Tests if all components of the vector are zeroes.
#[inline] #[inline]
pub fn is_zero(&self) -> bool { pub fn is_zero(&self) -> bool {
self.as_slice().iter().all(|e| e.is_zero()) self.as_ref().iter().all(|e| e.is_zero())
} }
} }
impl<N> $dvec<N> { impl<N> $dvec<N> {
/// Slices this vector. /// Slices this vector.
#[inline] #[inline]
pub fn as_slice<'a>(&'a self) -> &'a [N] { pub fn as_ref<'a>(&'a self) -> &'a [N] {
&self.at[.. self.len()] &self.at[.. self.len()]
} }
/// Mutably slices this vector. /// Mutably slices this vector.
#[inline] #[inline]
pub fn as_mut_slice<'a>(&'a mut self) -> &'a mut [N] { pub fn as_mut<'a>(&'a mut self) -> &'a mut [N] {
let len = self.len(); let len = self.len();
&mut self.at[.. len] &mut self.at[.. len]
} }
@ -46,7 +46,7 @@ macro_rules! dvec_impl(
fn swap(&mut self, i: usize, j: usize) { fn swap(&mut self, i: usize, j: usize) {
assert!(i < self.len()); assert!(i < self.len());
assert!(j < self.len()); assert!(j < self.len());
self.as_mut_slice().swap(i, j); self.as_mut().swap(i, j);
} }
#[inline] #[inline]
@ -61,17 +61,17 @@ macro_rules! dvec_impl(
} }
impl<N> Index<usize> for $dvec<N> { impl<N, T> Index<T> for $dvec<N> where [N]: Index<T> {
type Output = N; type Output = <[N] as Index<T>>::Output;
fn index(&self, i: usize) -> &N { fn index(&self, i: T) -> &<[N] as Index<T>>::Output {
&self.as_slice()[i] &self.as_ref()[i]
} }
} }
impl<N> IndexMut<usize> for $dvec<N> { impl<N, T> IndexMut<T> for $dvec<N> where [N]: IndexMut<T> {
fn index_mut(&mut self, i: usize) -> &mut N { fn index_mut(&mut self, i: T) -> &mut <[N] as Index<T>>::Output {
&mut self.as_mut_slice()[i] &mut self.as_mut()[i]
} }
} }
@ -97,14 +97,14 @@ macro_rules! dvec_impl(
impl<N> Iterable<N> for $dvec<N> { impl<N> Iterable<N> for $dvec<N> {
#[inline] #[inline]
fn iter<'l>(&'l self) -> Iter<'l, N> { fn iter<'l>(&'l self) -> Iter<'l, N> {
self.as_slice().iter() self.as_ref().iter()
} }
} }
impl<N> IterableMut<N> for $dvec<N> { impl<N> IterableMut<N> for $dvec<N> {
#[inline] #[inline]
fn iter_mut<'l>(&'l mut self) -> IterMut<'l, N> { fn iter_mut<'l>(&'l mut self) -> IterMut<'l, N> {
self.as_mut_slice().iter_mut() self.as_mut().iter_mut()
} }
} }
@ -185,7 +185,7 @@ macro_rules! dvec_impl(
let mut res = self; let mut res = self;
for (left, right) in res.as_mut_slice().iter_mut().zip(right.as_slice().iter()) { for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
*left = *left * *right *left = *left * *right
} }
@ -202,7 +202,7 @@ macro_rules! dvec_impl(
let mut res = self; let mut res = self;
for (left, right) in res.as_mut_slice().iter_mut().zip(right.as_slice().iter()) { for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
*left = *left / *right *left = *left / *right
} }
@ -219,7 +219,7 @@ macro_rules! dvec_impl(
let mut res = self; let mut res = self;
for (left, right) in res.as_mut_slice().iter_mut().zip(right.as_slice().iter()) { for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
*left = *left + *right *left = *left + *right
} }
@ -236,7 +236,7 @@ macro_rules! dvec_impl(
let mut res = self; let mut res = self;
for (left, right) in res.as_mut_slice().iter_mut().zip(right.as_slice().iter()) { for (left, right) in res.as_mut().iter_mut().zip(right.as_ref().iter()) {
*left = *left - *right *left = *left - *right
} }
@ -249,7 +249,7 @@ macro_rules! dvec_impl(
#[inline] #[inline]
fn neg(self) -> $dvec<N> { fn neg(self) -> $dvec<N> {
FromIterator::from_iter(self.as_slice().iter().map(|a| -*a)) FromIterator::from_iter(self.as_ref().iter().map(|a| -*a))
} }
} }
@ -282,7 +282,7 @@ macro_rules! dvec_impl(
fn normalize_mut(&mut self) -> N { fn normalize_mut(&mut self) -> N {
let l = Norm::norm(self); let l = Norm::norm(self);
for n in self.as_mut_slice().iter_mut() { for n in self.as_mut().iter_mut() {
*n = *n / l; *n = *n / l;
} }
@ -290,6 +290,14 @@ macro_rules! dvec_impl(
} }
} }
impl<N: BaseFloat + Cast<f64>> Mean<N> for $dvec<N> {
#[inline]
fn mean(&self) -> N {
let normalizer = ::cast(1.0f64 / self.len() as f64);
self.iter().fold(::zero(), |acc, x| acc + *x * normalizer)
}
}
impl<N: ApproxEq<N>> ApproxEq<N> for $dvec<N> { impl<N: ApproxEq<N>> ApproxEq<N> for $dvec<N> {
#[inline] #[inline]
fn approx_epsilon(_: Option<$dvec<N>>) -> N { fn approx_epsilon(_: Option<$dvec<N>>) -> N {
@ -303,13 +311,13 @@ macro_rules! dvec_impl(
#[inline] #[inline]
fn approx_eq_eps(&self, other: &$dvec<N>, epsilon: &N) -> bool { fn approx_eq_eps(&self, other: &$dvec<N>, epsilon: &N) -> bool {
let mut zip = self.as_slice().iter().zip(other.as_slice().iter()); let mut zip = self.as_ref().iter().zip(other.as_ref().iter());
zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon)) zip.all(|(a, b)| ApproxEq::approx_eq_eps(a, b, epsilon))
} }
#[inline] #[inline]
fn approx_eq_ulps(&self, other: &$dvec<N>, ulps: u32) -> bool { fn approx_eq_ulps(&self, other: &$dvec<N>, ulps: u32) -> bool {
let mut zip = self.as_slice().iter().zip(other.as_slice().iter()); let mut zip = self.as_ref().iter().zip(other.as_ref().iter());
zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps)) zip.all(|(a, b)| ApproxEq::approx_eq_ulps(a, b, ulps))
} }
} }
@ -321,7 +329,7 @@ macro_rules! dvec_impl(
fn mul(self, right: N) -> $dvec<N> { fn mul(self, right: N) -> $dvec<N> {
let mut res = self; let mut res = self;
for e in res.as_mut_slice().iter_mut() { for e in res.as_mut().iter_mut() {
*e = *e * right *e = *e * right
} }
@ -336,7 +344,7 @@ macro_rules! dvec_impl(
fn div(self, right: N) -> $dvec<N> { fn div(self, right: N) -> $dvec<N> {
let mut res = self; let mut res = self;
for e in res.as_mut_slice().iter_mut() { for e in res.as_mut().iter_mut() {
*e = *e / right *e = *e / right
} }
@ -351,7 +359,7 @@ macro_rules! dvec_impl(
fn add(self, right: N) -> $dvec<N> { fn add(self, right: N) -> $dvec<N> {
let mut res = self; let mut res = self;
for e in res.as_mut_slice().iter_mut() { for e in res.as_mut().iter_mut() {
*e = *e + right *e = *e + right
} }
@ -366,7 +374,7 @@ macro_rules! dvec_impl(
fn sub(self, right: N) -> $dvec<N> { fn sub(self, right: N) -> $dvec<N> {
let mut res = self; let mut res = self;
for e in res.as_mut_slice().iter_mut() { for e in res.as_mut().iter_mut() {
*e = *e - right *e = *e - right
} }
@ -379,10 +387,24 @@ macro_rules! dvec_impl(
macro_rules! small_dvec_impl ( macro_rules! small_dvec_impl (
($dvec: ident, $dim: expr, $($idx: expr),*) => ( ($dvec: ident, $dim: expr, $($idx: expr),*) => (
impl<N> $dvec<N> { impl<N> $dvec<N> {
/// The number of elements of this vector.
#[inline] #[inline]
pub fn len(&self) -> usize { pub fn len(&self) -> usize {
self.dim self.dim
} }
/// Creates an uninitialized vec of dimension `dim`.
#[inline]
pub unsafe fn new_uninitialized(dim: usize) -> $dvec<N> {
assert!(dim <= $dim, "The chosen dimension is too high for that type of \
stack-allocated dynamic vector. Consider using the \
heap-allocated vector: DVec.");
$dvec {
at: mem::uninitialized(),
dim: dim
}
}
} }
impl<N: PartialEq> PartialEq for $dvec<N> { impl<N: PartialEq> PartialEq for $dvec<N> {
@ -392,7 +414,7 @@ macro_rules! small_dvec_impl (
return false; // FIXME: fail instead? return false; // FIXME: fail instead?
} }
for (a, b) in self.as_slice().iter().zip(other.as_slice().iter()) { for (a, b) in self.as_ref().iter().zip(other.as_ref().iter()) {
if *a != *b { if *a != *b {
return false; return false;
} }

8
src/structs/mat.rs
View File

@ -14,7 +14,7 @@ use structs::dvec::{DVec1, DVec2, DVec3, DVec4, DVec5, DVec6};
use traits::structure::{Cast, Row, Col, Iterable, IterableMut, Dim, Indexable, Eye, ColSlice, use traits::structure::{Cast, Row, Col, Iterable, IterableMut, Dim, Indexable, Eye, ColSlice,
RowSlice, Diag, DiagMut, Shape, BaseFloat, BaseNum, Repeat}; RowSlice, Diag, DiagMut, Shape, BaseFloat, BaseNum, Repeat};
use traits::operations::{Absolute, Transpose, Inv, Outer, EigenQR}; use traits::operations::{Absolute, Transpose, Inv, Outer, EigenQR, Mean};
use traits::geometry::{ToHomogeneous, FromHomogeneous, Orig}; use traits::geometry::{ToHomogeneous, FromHomogeneous, Orig};
use linalg; use linalg;
#[cfg(feature="arbitrary")] #[cfg(feature="arbitrary")]
@ -81,6 +81,7 @@ outer_impl!(Vec1, Mat1);
eigen_qr_impl!(Mat1, Vec1); eigen_qr_impl!(Mat1, Vec1);
arbitrary_impl!(Mat1, m11); arbitrary_impl!(Mat1, m11);
rand_impl!(Mat1, m11); rand_impl!(Mat1, m11);
mean_impl!(Mat1, Vec1, 1);
/// Square matrix of dimension 2. /// Square matrix of dimension 2.
#[repr(C)] #[repr(C)]
@ -134,6 +135,7 @@ outer_impl!(Vec2, Mat2);
eigen_qr_impl!(Mat2, Vec2); eigen_qr_impl!(Mat2, Vec2);
arbitrary_impl!(Mat2, m11, m12, m21, m22); arbitrary_impl!(Mat2, m11, m12, m21, m22);
rand_impl!(Mat2, m11, m12, m21, m22); rand_impl!(Mat2, m11, m12, m21, m22);
mean_impl!(Mat2, Vec2, 2);
/// Square matrix of dimension 3. /// Square matrix of dimension 3.
#[repr(C)] #[repr(C)]
@ -230,6 +232,7 @@ rand_impl!(Mat3,
m21, m22, m23, m21, m22, m23,
m31, m32, m33 m31, m32, m33
); );
mean_impl!(Mat3, Vec3, 3);
/// Square matrix of dimension 4. /// Square matrix of dimension 4.
#[repr(C)] #[repr(C)]
@ -349,6 +352,7 @@ rand_impl!(Mat4,
m31, m32, m33, m34, m31, m32, m33, m34,
m41, m42, m43, m44 m41, m42, m43, m44
); );
mean_impl!(Mat4, Vec4, 4);
/// Square matrix of dimension 5. /// Square matrix of dimension 5.
#[repr(C)] #[repr(C)]
@ -485,6 +489,7 @@ rand_impl!(Mat5,
m41, m42, m43, m44, m45, m41, m42, m43, m44, m45,
m51, m52, m53, m54, m55 m51, m52, m53, m54, m55
); );
mean_impl!(Mat5, Vec5, 5);
/// Square matrix of dimension 6. /// Square matrix of dimension 6.
#[repr(C)] #[repr(C)]
@ -626,3 +631,4 @@ rand_impl!(Mat6,
m51, m52, m53, m54, m55, m56, m51, m52, m53, m54, m55, m56,
m61, m62, m63, m64, m65, m66 m61, m62, m63, m64, m65, m66
); );
mean_impl!(Mat6, Vec6, 6);

23
src/structs/mat_macros.rs
View File

@ -726,3 +726,26 @@ macro_rules! eigen_qr_impl(
} }
) )
); );
macro_rules! mean_impl(
($t: ident, $v: ident, $dim: expr) => (
impl<N: BaseNum + Cast<f64> + Clone> Mean<$v<N>> for $t<N> {
fn mean(&self) -> $v<N> {
let mut res: $v<N> = ::zero();
let normalizer: N = Cast::from(1.0f64 / $dim as f64);
for i in 0 .. $dim {
for j in 0 .. $dim {
unsafe {
let acc = res.unsafe_at(j) + self.unsafe_at((i, j)) * normalizer;
res.unsafe_set(j, acc);
}
}
}
res
}
}
)
);

3
src/structs/mod.rs
View File

@ -1,6 +1,6 @@
//! Data structures and implementations. //! Data structures and implementations.
pub use self::dmat::DMat; pub use self::dmat::{DMat, DMat1, DMat2, DMat3, DMat4, DMat5, DMat6};
pub use self::dvec::{DVec, DVec1, DVec2, DVec3, DVec4, DVec5, DVec6}; 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::vec::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6};
pub use self::pnt::{Pnt0, Pnt1, Pnt2, Pnt3, Pnt4, Pnt5, Pnt6}; pub use self::pnt::{Pnt0, Pnt1, Pnt2, Pnt3, Pnt4, Pnt5, Pnt6};
@ -11,6 +11,7 @@ pub use self::persp::{Persp3, PerspMat3};
pub use self::ortho::{Ortho3, OrthoMat3}; pub use self::ortho::{Ortho3, OrthoMat3};
pub use self::quat::{Quat, UnitQuat}; pub use self::quat::{Quat, UnitQuat};
mod dmat_macros;
mod dmat; mod dmat;
mod dvec_macros; mod dvec_macros;
mod dvec; mod dvec;

8
src/structs/vec.rs
View File

@ -9,7 +9,7 @@ use std::slice::{Iter, IterMut};
use std::iter::{Iterator, FromIterator, IntoIterator}; use std::iter::{Iterator, FromIterator, IntoIterator};
use rand::{Rand, Rng}; use rand::{Rand, Rng};
use num::{Zero, One}; use num::{Zero, One};
use traits::operations::{ApproxEq, POrd, POrdering, Axpy, Absolute}; use traits::operations::{ApproxEq, POrd, POrdering, Axpy, Absolute, Mean};
use traits::geometry::{Transform, Rotate, FromHomogeneous, ToHomogeneous, Dot, Norm, use traits::geometry::{Transform, Rotate, FromHomogeneous, ToHomogeneous, Dot, Norm,
Translation, Translate}; Translation, Translate};
use traits::structure::{Basis, Cast, Dim, Indexable, Iterable, IterableMut, Shape, NumVec, use traits::structure::{Basis, Cast, Dim, Indexable, Iterable, IterableMut, Shape, NumVec,
@ -90,6 +90,7 @@ num_float_vec_impl!(Vec1);
absolute_vec_impl!(Vec1, x); absolute_vec_impl!(Vec1, x);
arbitrary_impl!(Vec1, x); arbitrary_impl!(Vec1, x);
rand_impl!(Vec1, x); rand_impl!(Vec1, x);
mean_impl!(Vec1);
/// Vector of dimension 2. /// Vector of dimension 2.
#[repr(C)] #[repr(C)]
@ -143,6 +144,7 @@ num_float_vec_impl!(Vec2);
absolute_vec_impl!(Vec2, x, y); absolute_vec_impl!(Vec2, x, y);
arbitrary_impl!(Vec2, x, y); arbitrary_impl!(Vec2, x, y);
rand_impl!(Vec2, x, y); rand_impl!(Vec2, x, y);
mean_impl!(Vec2);
/// Vector of dimension 3. /// Vector of dimension 3.
#[repr(C)] #[repr(C)]
@ -198,6 +200,7 @@ num_float_vec_impl!(Vec3);
absolute_vec_impl!(Vec3, x, y, z); absolute_vec_impl!(Vec3, x, y, z);
arbitrary_impl!(Vec3, x, y, z); arbitrary_impl!(Vec3, x, y, z);
rand_impl!(Vec3, x, y, z); rand_impl!(Vec3, x, y, z);
mean_impl!(Vec3);
/// Vector of dimension 4. /// Vector of dimension 4.
@ -256,6 +259,7 @@ num_float_vec_impl!(Vec4);
absolute_vec_impl!(Vec4, x, y, z, w); absolute_vec_impl!(Vec4, x, y, z, w);
arbitrary_impl!(Vec4, x, y, z, w); arbitrary_impl!(Vec4, x, y, z, w);
rand_impl!(Vec4, x, y, z, w); rand_impl!(Vec4, x, y, z, w);
mean_impl!(Vec4);
/// Vector of dimension 5. /// Vector of dimension 5.
#[repr(C)] #[repr(C)]
@ -315,6 +319,7 @@ num_float_vec_impl!(Vec5);
absolute_vec_impl!(Vec5, x, y, z, w, a); absolute_vec_impl!(Vec5, x, y, z, w, a);
arbitrary_impl!(Vec5, x, y, z, w, a); arbitrary_impl!(Vec5, x, y, z, w, a);
rand_impl!(Vec5, x, y, z, w, a); rand_impl!(Vec5, x, y, z, w, a);
mean_impl!(Vec5);
/// Vector of dimension 6. /// Vector of dimension 6.
#[repr(C)] #[repr(C)]
@ -374,3 +379,4 @@ num_float_vec_impl!(Vec6);
absolute_vec_impl!(Vec6, x, y, z, w, a, b); absolute_vec_impl!(Vec6, x, y, z, w, a, b);
arbitrary_impl!(Vec6, x, y, z, w, a, b); arbitrary_impl!(Vec6, x, y, z, w, a, b);
rand_impl!(Vec6, x, y, z, w, a, b); rand_impl!(Vec6, x, y, z, w, a, b);
mean_impl!(Vec6);

22
src/structs/vec_macros.rs
View File

@ -208,16 +208,16 @@ macro_rules! indexable_impl(
macro_rules! index_impl( macro_rules! index_impl(
($t: ident) => ( ($t: ident) => (
impl<N> Index<usize> for $t<N> { impl<N, T> Index<T> for $t<N> where [N]: Index<T> {
type Output = N; type Output = <[N] as Index<T>>::Output;
fn index(&self, i: usize) -> &N { fn index(&self, i: T) -> &<[N] as Index<T>>::Output {
&self.as_ref()[i] &self.as_ref()[i]
} }
} }
impl<N> IndexMut<usize> for $t<N> { impl<N, T> IndexMut<T> for $t<N> where [N]: IndexMut<T> {
fn index_mut(&mut self, i: usize) -> &mut N { fn index_mut(&mut self, i: T) -> &mut <[N] as Index<T>>::Output {
&mut self.as_mut()[i] &mut self.as_mut()[i]
} }
} }
@ -805,3 +805,15 @@ macro_rules! rand_impl(
} }
) )
); );
macro_rules! mean_impl(
($t: ident) => (
impl<N: BaseFloat + Cast<f64>> Mean<N> for $t<N> {
#[inline]
fn mean(&self) -> N {
let normalizer = ::cast(1.0f64 / self.len() as f64);
self.iter().fold(::zero(), |acc, x| acc + *x * normalizer)
}
}
)
);

2
src/traits/operations.rs
View File

@ -325,7 +325,7 @@ pub trait Cov<M> {
} }
} }
/// Trait for computing the covariance of a set of data. /// Trait for computing the mean of a set of data.
pub trait Mean<N> { pub trait Mean<N> {
/// Computes the mean of the observations stored by `v`. /// Computes the mean of the observations stored by `v`.
/// ///

72
tests/mat.rs
View File

@ -7,7 +7,7 @@ use na::{Vec1, Vec3, Mat1, Mat2, Mat3, Mat4, Mat5, Mat6, Rot2, Rot3, Persp3, Per
macro_rules! test_inv_mat_impl( macro_rules! test_inv_mat_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat : $t = random(); let randmat : $t = random();
match na::inv(&randmat) { match na::inv(&randmat) {
@ -20,7 +20,7 @@ macro_rules! test_inv_mat_impl(
macro_rules! test_transpose_mat_impl( macro_rules! test_transpose_mat_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat : $t = random(); let randmat : $t = random();
assert!(na::transpose(&na::transpose(&randmat)) == randmat); assert!(na::transpose(&na::transpose(&randmat)) == randmat);
@ -30,7 +30,7 @@ macro_rules! test_transpose_mat_impl(
macro_rules! test_qr_impl( macro_rules! test_qr_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat : $t = random(); let randmat : $t = random();
let (q, r) = na::qr(&randmat); let (q, r) = na::qr(&randmat);
@ -43,7 +43,7 @@ macro_rules! test_qr_impl(
macro_rules! test_cholesky_impl( macro_rules! test_cholesky_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
// construct symmetric positive definite matrix // construct symmetric positive definite matrix
let mut randmat : $t = random(); let mut randmat : $t = random();
@ -65,7 +65,7 @@ macro_rules! test_cholesky_impl(
macro_rules! test_hessenberg_impl( macro_rules! test_hessenberg_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat : $t = random(); let randmat : $t = random();
@ -90,7 +90,7 @@ macro_rules! test_hessenberg_impl(
macro_rules! test_eigen_qr_impl( macro_rules! test_eigen_qr_impl(
($t: ty) => { ($t: ty) => {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat : $t = random(); let randmat : $t = random();
// Make it symetric so that we can recompose the matrix to test at the end. // Make it symetric so that we can recompose the matrix to test at the end.
let randmat = na::transpose(&randmat) * randmat; let randmat = na::transpose(&randmat) * randmat;
@ -105,7 +105,7 @@ macro_rules! test_eigen_qr_impl(
assert!(na::approx_eq_eps(&randmat, &recomp, &1.0e-2)); assert!(na::approx_eq_eps(&randmat, &recomp, &1.0e-2));
} }
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat : $t = random(); let randmat : $t = random();
// Take only diagonal part // Take only diagonal part
let randmat: $t = Diag::from_diag(&randmat.diag()); let randmat: $t = Diag::from_diag(&randmat.diag());
@ -184,7 +184,7 @@ fn test_inv_mat6() {
#[test] #[test]
fn test_rotation2() { fn test_rotation2() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat: na::Rot2<f64> = na::one(); let randmat: na::Rot2<f64> = na::one();
let ang = Vec1::new(na::abs(&random::<f64>()) % <f64 as BaseFloat>::pi()); let ang = Vec1::new(na::abs(&random::<f64>()) % <f64 as BaseFloat>::pi());
@ -201,7 +201,7 @@ fn test_index_mat2() {
#[test] #[test]
fn test_inv_rotation3() { fn test_inv_rotation3() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let randmat: Rot3<f64> = na::one(); let randmat: Rot3<f64> = na::one();
let dir: Vec3<f64> = random(); let dir: Vec3<f64> = random();
let ang = na::normalize(&dir) * (na::abs(&random::<f64>()) % <f64 as BaseFloat>::pi()); let ang = na::normalize(&dir) * (na::abs(&random::<f64>()) % <f64 as BaseFloat>::pi());
@ -315,6 +315,33 @@ fn test_transpose_dmat() {
assert!(na::transpose(&na::transpose(&mat)) == mat); assert!(na::transpose(&na::transpose(&mat)) == mat);
} }
#[test]
fn test_row_dmat() {
let mat = DMat::from_row_vec(
8,
4,
&[
1u32,2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 14, 15, 16,
17, 18, 19, 20,
21, 22, 23, 24,
25, 26, 27, 28,
29, 30, 31, 32
]
);
assert_eq!(&DVec::from_slice(4, &[1u32, 2, 3, 4]), &mat.row(0));
assert_eq!(&DVec::from_slice(4, &[5u32, 6, 7, 8]), &mat.row(1));
assert_eq!(&DVec::from_slice(4, &[9u32, 10, 11, 12]), &mat.row(2));
assert_eq!(&DVec::from_slice(4, &[13u32, 14, 15, 16]), &mat.row(3));
assert_eq!(&DVec::from_slice(4, &[17u32, 18, 19, 20]), &mat.row(4));
assert_eq!(&DVec::from_slice(4, &[21u32, 22, 23, 24]), &mat.row(5));
assert_eq!(&DVec::from_slice(4, &[25u32, 26, 27, 28]), &mat.row(6));
assert_eq!(&DVec::from_slice(4, &[29u32, 30, 31, 32]), &mat.row(7));
}
#[test] #[test]
fn test_row_slice_dmat() { fn test_row_slice_dmat() {
let mat = DMat::from_row_vec( let mat = DMat::from_row_vec(
@ -335,6 +362,29 @@ fn test_row_slice_dmat() {
assert_eq!(&DVec::from_slice(2, &[19u32, 20]), &mat.row_slice(4, 2, 4)); assert_eq!(&DVec::from_slice(2, &[19u32, 20]), &mat.row_slice(4, 2, 4));
} }
#[test]
fn test_col_dmat() {
let mat = DMat::from_row_vec(
8,
4,
&[
1u32,2, 3, 4,
5, 6, 7, 8,
9, 10, 11, 12,
13, 14, 15, 16,
17, 18, 19, 20,
21, 22, 23, 24,
25, 26, 27, 28,
29, 30, 31, 32
]
);
assert_eq!(&DVec::from_slice(8, &[1u32, 5, 9, 13, 17, 21, 25, 29]), &mat.col(0));
assert_eq!(&DVec::from_slice(8, &[2u32, 6, 10, 14, 18, 22, 26, 30]), &mat.col(1));
assert_eq!(&DVec::from_slice(8, &[3u32, 7, 11, 15, 19, 23, 27, 31]), &mat.col(2));
assert_eq!(&DVec::from_slice(8, &[4u32, 8, 12, 16, 20, 24, 28, 32]), &mat.col(3));
}
#[test] #[test]
fn test_col_slice_dmat() { fn test_col_slice_dmat() {
let mat = DMat::from_row_vec( let mat = DMat::from_row_vec(
@ -608,7 +658,7 @@ fn test_dmat_set_row() {
/* FIXME: review qr decomposition to make it work with DMat. /* FIXME: review qr decomposition to make it work with DMat.
#[test] #[test]
fn test_qr() { fn test_qr() {
for _ in (0usize .. 10) { for _ in 0usize .. 10 {
let dim1: usize = random(); let dim1: usize = random();
let dim2: usize = random(); let dim2: usize = random();
let rows = min(40, max(dim1, dim2)); let rows = min(40, max(dim1, dim2));
@ -893,6 +943,6 @@ fn test_transpose_square_mat() {
let mut mat = DMat::from_col_vec(num_rows, num_cols, col_major_mat); let mut mat = DMat::from_col_vec(num_rows, num_cols, col_major_mat);
mat.transpose_mut(); mat.transpose_mut();
for i in 0..num_rows { for i in 0..num_rows {
assert_eq!(&[0, 1, 2, 3], mat.row_slice(i, 0, num_cols).as_slice()); assert_eq!(&[0, 1, 2, 3], &mat.row_slice(i, 0, num_cols)[..]);
} }
} }

10
tests/quat.rs
View File

@ -6,7 +6,7 @@ use rand::random;
#[test] #[test]
fn test_quat_as_mat() { fn test_quat_as_mat() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let axis_angle: Vec3<f64> = random(); let axis_angle: Vec3<f64> = random();
assert!(na::approx_eq(&UnitQuat::new(axis_angle).to_rot(), &Rot3::new(axis_angle))) assert!(na::approx_eq(&UnitQuat::new(axis_angle).to_rot(), &Rot3::new(axis_angle)))
@ -15,7 +15,7 @@ fn test_quat_as_mat() {
#[test] #[test]
fn test_quat_mul_vec_or_pnt_as_mat() { fn test_quat_mul_vec_or_pnt_as_mat() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let axis_angle: Vec3<f64> = random(); let axis_angle: Vec3<f64> = random();
let vec: Vec3<f64> = random(); let vec: Vec3<f64> = random();
let pnt: Pnt3<f64> = random(); let pnt: Pnt3<f64> = random();
@ -32,7 +32,7 @@ fn test_quat_mul_vec_or_pnt_as_mat() {
#[test] #[test]
fn test_quat_div_quat() { fn test_quat_div_quat() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let axis_angle1: Vec3<f64> = random(); let axis_angle1: Vec3<f64> = random();
let axis_angle2: Vec3<f64> = random(); let axis_angle2: Vec3<f64> = random();
@ -48,7 +48,7 @@ fn test_quat_div_quat() {
#[test] #[test]
fn test_quat_to_axis_angle() { fn test_quat_to_axis_angle() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let axis_angle: Vec3<f64> = random(); let axis_angle: Vec3<f64> = random();
let q = UnitQuat::new(axis_angle); let q = UnitQuat::new(axis_angle);
@ -60,7 +60,7 @@ fn test_quat_to_axis_angle() {
#[test] #[test]
fn test_quat_euler_angles() { fn test_quat_euler_angles() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let angles: Vec3<f64> = random(); let angles: Vec3<f64> = random();
let q = UnitQuat::new_with_euler_angles(angles.x, angles.y, angles.z); let q = UnitQuat::new_with_euler_angles(angles.x, angles.y, angles.z);

20
tests/vec.rs
View File

@ -6,7 +6,7 @@ use na::{Vec0, Vec1, Vec2, Vec3, Vec4, Vec5, Vec6, Mat3, Rot2, Rot3, Iterable, I
macro_rules! test_iterator_impl( macro_rules! test_iterator_impl(
($t: ty, $n: ty) => ( ($t: ty, $n: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v: $t = random(); let v: $t = random();
let mut mv: $t = v.clone(); let mut mv: $t = v.clone();
let n: $n = random(); let n: $n = random();
@ -24,7 +24,7 @@ macro_rules! test_iterator_impl(
macro_rules! test_commut_dot_impl( macro_rules! test_commut_dot_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v1 : $t = random(); let v1 : $t = random();
let v2 : $t = random(); let v2 : $t = random();
@ -35,7 +35,7 @@ macro_rules! test_commut_dot_impl(
macro_rules! test_scalar_op_impl( macro_rules! test_scalar_op_impl(
($t: ty, $n: ty) => ( ($t: ty, $n: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v1 : $t = random(); let v1 : $t = random();
let n : $n = random(); let n : $n = random();
@ -58,7 +58,7 @@ macro_rules! test_scalar_op_impl(
macro_rules! test_basis_impl( macro_rules! test_basis_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
na::canonical_basis(|e1: $t| { na::canonical_basis(|e1: $t| {
na::canonical_basis(|e2: $t| { na::canonical_basis(|e2: $t| {
assert!(e1 == e2 || na::approx_eq(&na::dot(&e1, &e2), &na::zero())); assert!(e1 == e2 || na::approx_eq(&na::dot(&e1, &e2), &na::zero()));
@ -76,7 +76,7 @@ macro_rules! test_basis_impl(
macro_rules! test_subspace_basis_impl( macro_rules! test_subspace_basis_impl(
($t: ty) => ( ($t: ty) => (
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v : $t = random(); let v : $t = random();
let v1 = na::normalize(&v); let v1 = na::normalize(&v);
@ -100,7 +100,7 @@ macro_rules! test_subspace_basis_impl(
#[test] #[test]
fn test_cross_vec3() { fn test_cross_vec3() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v1 : Vec3<f64> = random(); let v1 : Vec3<f64> = random();
let v2 : Vec3<f64> = random(); let v2 : Vec3<f64> = random();
let v3 : Vec3<f64> = na::cross(&v1, &v2); let v3 : Vec3<f64> = na::cross(&v1, &v2);
@ -321,7 +321,7 @@ fn test_outer_vec3() {
#[test] #[test]
fn test_vec3_rotation_between() { fn test_vec3_rotation_between() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v1: Vec3<f64> = random(); let v1: Vec3<f64> = random();
let mut v2: Vec3<f64> = random(); let mut v2: Vec3<f64> = random();
@ -335,7 +335,7 @@ fn test_vec3_rotation_between() {
#[test] #[test]
fn test_vec3_angle_between() { fn test_vec3_angle_between() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let vec: Vec3<f64> = random(); let vec: Vec3<f64> = random();
let other: Vec3<f64> = random(); let other: Vec3<f64> = random();
@ -353,7 +353,7 @@ fn test_vec3_angle_between() {
#[test] #[test]
fn test_vec2_rotation_between() { fn test_vec2_rotation_between() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let v1: Vec2<f64> = random(); let v1: Vec2<f64> = random();
let mut v2: Vec2<f64> = random(); let mut v2: Vec2<f64> = random();
@ -367,7 +367,7 @@ fn test_vec2_rotation_between() {
#[test] #[test]
fn test_vec2_angle_between() { fn test_vec2_angle_between() {
for _ in (0usize .. 10000) { for _ in 0usize .. 10000 {
let axis_ang: Vec1<f64> = random(); let axis_ang: Vec1<f64> = random();
let ang = na::norm(&axis_ang); let ang = na::norm(&axis_ang);