Add dynamically sized matrix.

This commit is contained in:
Sébastien Crozet 2013-05-31 11:28:42 +02:00
parent 335794208d
commit f264b75ce6
5 changed files with 265 additions and 131 deletions

View File

@ -39,7 +39,7 @@ pub mod dim3
/// n-dimensional linear algebra (slower). /// n-dimensional linear algebra (slower).
pub mod ndim pub mod ndim
{ {
// pub mod dmat; pub mod dmat;
pub mod dvec; pub mod dvec;
pub mod nvec; pub mod nvec;
pub mod nmat; pub mod nmat;

244
src/ndim/dmat.rs Normal file
View File

@ -0,0 +1,244 @@
use core::num::{One, Zero};
use core::vec::{from_elem, swap, all, all2, len};
use core::cmp::ApproxEq;
use traits::inv::Inv;
use traits::transpose::Transpose;
use traits::workarounds::rlmul::{RMul, LMul};
use ndim::dvec::{DVec, zero_vec_with_dim};
#[deriving(Eq, ToStr, Clone)]
pub struct DMat<T>
{
dim: uint, // FIXME: handle more than just square matrices
mij: ~[T]
}
pub fn zero_mat_with_dim<T: Zero + Copy>(dim: uint) -> DMat<T>
{ DMat { dim: dim, mij: from_elem(dim * dim, Zero::zero()) } }
pub fn is_zero_mat<T: Zero>(mat: &DMat<T>) -> bool
{ all(mat.mij, |e| e.is_zero()) }
pub fn one_mat_with_dim<T: Copy + One + Zero>(dim: uint) -> DMat<T>
{
let mut res = zero_mat_with_dim(dim);
let _1 = One::one::<T>();
for uint::range(0u, dim) |i|
{ res.set(i, i, &_1); }
res
}
impl<T: Copy> DMat<T>
{
pub fn offset(&self, i: uint, j: uint) -> uint
{ i * self.dim + j }
pub fn set(&mut self, i: uint, j: uint, t: &T)
{
assert!(i < self.dim);
assert!(j < self.dim);
self.mij[self.offset(i, j)] = *t
}
}
impl<T: Copy> Index<(uint, uint), T> for DMat<T>
{
fn index(&self, &(i, j): &(uint, uint)) -> T
{ self.mij[self.offset(i, j)] }
}
impl<T: Copy + Mul<T, T> + Add<T, T> + Zero>
Mul<DMat<T>, DMat<T>> for DMat<T>
{
fn mul(&self, other: &DMat<T>) -> DMat<T>
{
assert!(self.dim == other.dim);
let dim = self.dim;
let mut res = zero_mat_with_dim(dim);
for uint::range(0u, dim) |i|
{
for uint::range(0u, dim) |j|
{
let mut acc: T = Zero::zero();
for uint::range(0u, dim) |k|
{ acc += self[(i, k)] * other[(k, j)]; }
res.set(i, j, &acc);
}
}
res
}
}
impl<T: Copy + Add<T, T> + Mul<T, T> + Zero>
RMul<DVec<T>> for DMat<T>
{
fn rmul(&self, other: &DVec<T>) -> DVec<T>
{
assert!(self.dim == len(other.at));
let dim = self.dim;
let mut res : DVec<T> = zero_vec_with_dim(dim);
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
}
}
impl<T: Copy + Add<T, T> + Mul<T, T> + Zero>
LMul<DVec<T>> for DMat<T>
{
fn lmul(&self, other: &DVec<T>) -> DVec<T>
{
assert!(self.dim == len(other.at));
let dim = self.dim;
let mut res : DVec<T> = zero_vec_with_dim(dim);
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
}
}
impl<T: Clone + Copy + Eq + One + Zero +
Mul<T, T> + Div<T, T> + Sub<T, T> + Neg<T>>
Inv for DMat<T>
{
fn inverse(&self) -> DMat<T>
{
let mut res : DMat<T> = self.clone();
res.invert();
res
}
fn invert(&mut self)
{
let dim = self.dim;
let mut res = one_mat_with_dim::<T>(dim);
let _0T = Zero::zero::<T>();
// inversion using Gauss-Jordan elimination
for uint::range(0u, dim) |k|
{
// 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?
// FIXME: this is kind of uggly…
// … but we cannot use position_between since we are iterating on one
// columns
let mut n0 = 0u; // index of a non-zero entry
while (n0 != dim)
{
if (self[(n0, k)] != _0T)
{ break; }
n0 += 1;
}
assert!(n0 != dim); // non inversible matrix
// swap pivot line
if (n0 != k)
{
for uint::range(0u, dim) |j|
{
let off_n0_j = self.offset(n0, j);
let off_k_j = self.offset(k, j);
swap(self.mij, off_n0_j, off_k_j);
swap(res.mij, off_n0_j, off_k_j);
}
}
let pivot = self[(k, k)];
for uint::range(k, dim) |j|
{
// FIXME: not to putting selfal exression directly on the nuction call
// is uggly but does not seem to compile any more…
let selfval = &(self[(k, j)] / pivot);
let resval = &(res[(k, j)] / pivot);
self.set(k, j, selfval);
res.set(k, j, resval);
}
for uint::range(0u, dim) |l|
{
if (l != k)
{
let normalizer = self[(l, k)] / pivot;
for uint::range(k, dim) |j|
{
let selfval = &(self[(l, j)] - self[(k, j)] * normalizer);
let resval = &(res[(l, j)] - res[(k, j)] * normalizer);
self.set(k, j, selfval);
res.set(k, j, resval);
}
}
}
}
}
}
impl<T:Copy> Transpose for DMat<T>
{
fn transposed(&self) -> DMat<T>
{
let mut res = copy *self;
res.transpose();
res
}
fn transpose(&mut self)
{
let dim = self.dim;
for uint::range(1u, dim) |i|
{
for uint::range(0u, dim - 1) |j|
{
let off_i_j = self.offset(i, j);
let off_j_i = self.offset(j, i);
swap(self.mij, off_i_j, off_j_i);
}
}
}
}
impl<T: ApproxEq<T>> ApproxEq<T> for DMat<T>
{
fn approx_epsilon() -> T
{ ApproxEq::approx_epsilon::<T, T>() }
fn approx_eq(&self, other: &DMat<T>) -> bool
{ all2(self.mij, other.mij, |a, b| a.approx_eq(b)) }
fn approx_eq_eps(&self, other: &DMat<T>, epsilon: &T) -> bool
{ all2(self.mij, other.mij, |a, b| a.approx_eq_eps(b, epsilon)) }
}

View File

@ -15,10 +15,10 @@ pub struct DVec<T>
at: ~[T] at: ~[T]
} }
pub fn zero_with_dim<T: Zero + Copy>(dim: uint) -> DVec<T> pub fn zero_vec_with_dim<T: Zero + Copy>(dim: uint) -> DVec<T>
{ DVec { at: from_elem(dim, Zero::zero::<T>()) } } { DVec { at: from_elem(dim, Zero::zero::<T>()) } }
pub fn is_zero<T: Zero>(vec: &DVec<T>) -> bool pub fn is_zero_vec<T: Zero>(vec: &DVec<T>) -> bool
{ all(vec.at, |e| e.is_zero()) } { all(vec.at, |e| e.is_zero()) }
// FIXME: is Clone needed? // FIXME: is Clone needed?
@ -30,7 +30,7 @@ impl<T: Copy + DivisionRing + Algebraic + Clone + ApproxEq<T>> DVec<T>
for uint::range(0u, dim) |i| for uint::range(0u, dim) |i|
{ {
let mut basis_element : DVec<T> = zero_with_dim(dim); let mut basis_element : DVec<T> = zero_vec_with_dim(dim);
basis_element.at[i] = One::one(); basis_element.at[i] = One::one();
@ -49,7 +49,7 @@ impl<T: Copy + DivisionRing + Algebraic + Clone + ApproxEq<T>> DVec<T>
for uint::range(0u, dim) |i| for uint::range(0u, dim) |i|
{ {
let mut basis_element : DVec<T> = zero_with_dim(len(self.at)); let mut basis_element : DVec<T> = zero_vec_with_dim(len(self.at));
basis_element.at[i] = One::one(); basis_element.at[i] = One::one();

View File

@ -1,11 +1,11 @@
use core::num::{One, Zero}; use core::num::{One, Zero};
use core::rand::{Rand, Rng, RngUtil}; use core::rand::{Rand, Rng, RngUtil};
use core::vec::{from_elem, swap, all, all2};
use core::cmp::ApproxEq; use core::cmp::ApproxEq;
use traits::dim::Dim; use traits::dim::Dim;
use traits::inv::Inv; use traits::inv::Inv;
use traits::transpose::Transpose; use traits::transpose::Transpose;
use traits::workarounds::rlmul::{RMul, LMul}; use traits::workarounds::rlmul::{RMul, LMul};
use ndim::dmat::{DMat, one_mat_with_dim, zero_mat_with_dim, is_zero_mat};
use ndim::nvec::NVec; use ndim::nvec::NVec;
// D is a phantom type parameter, used only as a dimensional token. // D is a phantom type parameter, used only as a dimensional token.
@ -14,15 +14,7 @@ use ndim::nvec::NVec;
// using d0, d1, d2, d3 and d4 tokens are prefered. // using d0, d1, d2, d3 and d4 tokens are prefered.
#[deriving(Eq, ToStr)] #[deriving(Eq, ToStr)]
pub struct NMat<D, T> pub struct NMat<D, T>
{ { mij: DMat<T> }
mij: ~[T]
}
impl<D, T: Clone> Clone for NMat<D, T>
{
fn clone(&self) -> NMat<D, T>
{ NMat{ mij: self.mij.clone() } }
}
impl<D: Dim, T: Copy> NMat<D, T> impl<D: Dim, T: Copy> NMat<D, T>
{ {
@ -30,7 +22,7 @@ impl<D: Dim, T: Copy> NMat<D, T>
{ i * Dim::dim::<D>() + j } { i * Dim::dim::<D>() + j }
fn set(&mut self, i: uint, j: uint, t: &T) fn set(&mut self, i: uint, j: uint, t: &T)
{ self.mij[NMat::offset::<D, T>(i, j)] = *t } { self.mij.set(i, j, t) }
} }
impl<D: Dim, T> Dim for NMat<D, T> impl<D: Dim, T> Dim for NMat<D, T>
@ -41,36 +33,23 @@ impl<D: Dim, T> Dim for NMat<D, T>
impl<D: Dim, T: Copy> Index<(uint, uint), T> for NMat<D, T> impl<D: Dim, T: Copy> Index<(uint, uint), T> for NMat<D, T>
{ {
fn index(&self, &(i, j): &(uint, uint)) -> T fn index(&self, &idx: &(uint, uint)) -> T
{ self.mij[NMat::offset::<D, T>(i, j)] } { self.mij[idx] }
} }
impl<D: Dim, T: Copy + One + Zero> One for NMat<D, T> impl<D: Dim, T: Copy + One + Zero> One for NMat<D, T>
{ {
fn one() -> NMat<D, T> fn one() -> NMat<D, T>
{ { NMat { mij: one_mat_with_dim(Dim::dim::<D>()) } }
let dim = Dim::dim::<D>();
let mut res = NMat{ mij: from_elem(dim * dim, Zero::zero()) };
let _1 = One::one::<T>();
for uint::range(0u, dim) |i|
{ res.set(i, i, &_1); }
res
}
} }
impl<D: Dim, T: Copy + Zero> Zero for NMat<D, T> impl<D: Dim, T: Copy + Zero> Zero for NMat<D, T>
{ {
fn zero() -> NMat<D, T> fn zero() -> NMat<D, T>
{ { NMat { mij: zero_mat_with_dim(Dim::dim::<D>()) } }
let dim = Dim::dim::<D>();
NMat{ mij: from_elem(dim * dim, Zero::zero()) }
}
fn is_zero(&self) -> bool fn is_zero(&self) -> bool
{ all(self.mij, |e| e.is_zero()) } { is_zero_mat(&self.mij) }
} }
impl<D: Dim, T: Copy + Mul<T, T> + Add<T, T> + Zero> impl<D: Dim, T: Copy + Mul<T, T> + Add<T, T> + Zero>
@ -140,87 +119,10 @@ impl<D: Dim,
Inv for NMat<D, T> Inv for NMat<D, T>
{ {
fn inverse(&self) -> NMat<D, T> fn inverse(&self) -> NMat<D, T>
{ { NMat { mij: self.mij.inverse() } }
let mut res : NMat<D, T> = self.clone();
res.invert();
res
}
fn invert(&mut self) fn invert(&mut self)
{ { self.mij.invert() }
let dim = Dim::dim::<D>();
let mut res = One::one::<NMat<D, T>>();
let _0T = Zero::zero::<T>();
// inversion using Gauss-Jordan elimination
for uint::range(0u, dim) |k|
{
// 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?
// FIXME: this is kind of uggly…
// … but we cannot use position_between since we are iterating on one
// columns
let mut n0 = 0u; // index of a non-zero entry
while (n0 != dim)
{
if (self[(n0, k)] != _0T)
{ break; }
n0 += 1;
}
assert!(n0 != dim); // non inversible matrix
// swap pivot line
if (n0 != k)
{
for uint::range(0u, dim) |j|
{
swap(self.mij,
NMat::offset::<D, T>(n0, j),
NMat::offset::<D, T>(k, j));
swap(res.mij,
NMat::offset::<D, T>(n0, j),
NMat::offset::<D, T>(k, j));
}
}
let pivot = self[(k, k)];
for uint::range(k, dim) |j|
{
// FIXME: not to putting selfal exression directly on the nuction call
// is uggly but does not seem to compile any more…
let selfval = &(self[(k, j)] / pivot);
let resval = &(res[(k, j)] / pivot);
self.set(k, j, selfval);
res.set(k, j, resval);
}
for uint::range(0u, dim) |l|
{
if (l != k)
{
let normalizer = self[(l, k)] / pivot;
for uint::range(k, dim) |j|
{
let selfval = &(self[(l, j)] - self[(k, j)] * normalizer);
let resval = &(res[(l, j)] - res[(k, j)] * normalizer);
self.set(k, j, selfval);
res.set(k, j, resval);
}
}
}
}
}
} }
impl<D: Dim, T:Copy> Transpose for NMat<D, T> impl<D: Dim, T:Copy> Transpose for NMat<D, T>
@ -235,19 +137,7 @@ impl<D: Dim, T:Copy> Transpose for NMat<D, T>
} }
fn transpose(&mut self) fn transpose(&mut self)
{ { self.mij.transpose() }
let dim = Dim::dim::<D>();
for uint::range(1u, dim) |i|
{
for uint::range(0u, dim - 1) |j|
{
swap(self.mij,
NMat::offset::<D, T>(i, j),
NMat::offset::<D, T>(j, i));
}
}
}
} }
impl<D, T: ApproxEq<T>> ApproxEq<T> for NMat<D, T> impl<D, T: ApproxEq<T>> ApproxEq<T> for NMat<D, T>
@ -256,10 +146,10 @@ impl<D, T: ApproxEq<T>> ApproxEq<T> for NMat<D, T>
{ ApproxEq::approx_epsilon::<T, T>() } { ApproxEq::approx_epsilon::<T, T>() }
fn approx_eq(&self, other: &NMat<D, T>) -> bool fn approx_eq(&self, other: &NMat<D, T>) -> bool
{ all2(self.mij, other.mij, |a, b| a.approx_eq(b)) } { self.mij.approx_eq(&other.mij) }
fn approx_eq_eps(&self, other: &NMat<D, T>, epsilon: &T) -> bool fn approx_eq_eps(&self, other: &NMat<D, T>, epsilon: &T) -> bool
{ all2(self.mij, other.mij, |a, b| a.approx_eq_eps(b, epsilon)) } { self.mij.approx_eq_eps(&other.mij, epsilon) }
} }
impl<D: Dim, T: Rand + Zero + Copy> Rand for NMat<D, T> impl<D: Dim, T: Rand + Zero + Copy> Rand for NMat<D, T>

View File

@ -2,7 +2,7 @@ use core::num::{Zero, Algebraic};
use core::rand::{Rand, Rng, RngUtil}; use core::rand::{Rand, Rng, RngUtil};
use core::vec::{map}; use core::vec::{map};
use core::cmp::ApproxEq; use core::cmp::ApproxEq;
use ndim::dvec::{DVec, zero_with_dim, is_zero}; use ndim::dvec::{DVec, zero_vec_with_dim, is_zero_vec};
use traits::basis::Basis; use traits::basis::Basis;
use traits::ring::Ring; use traits::ring::Ring;
use traits::division_ring::DivisionRing; use traits::division_ring::DivisionRing;
@ -164,10 +164,10 @@ Basis for NVec<D, T>
impl<D: Dim, T: Copy + Zero> Zero for NVec<D, T> impl<D: Dim, T: Copy + Zero> Zero for NVec<D, T>
{ {
fn zero() -> NVec<D, T> fn zero() -> NVec<D, T>
{ NVec { at: zero_with_dim(Dim::dim::<D>()) } } { NVec { at: zero_vec_with_dim(Dim::dim::<D>()) } }
fn is_zero(&self) -> bool fn is_zero(&self) -> bool
{ is_zero(&self.at) } { is_zero_vec(&self.at) }
} }
impl<D, T: ApproxEq<T>> ApproxEq<T> for NVec<D, T> impl<D, T: ApproxEq<T>> ApproxEq<T> for NVec<D, T>