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

Fix stepping for slices. 

The previous implementation was wrong compared to what the documentatino claimed.
This commit is contained in:
Sébastien Crozet 2017-08-13 19:52:55 +02:00 committed by Sébastien Crozet
parent 41f5231446
commit c258d13f98
8 changed files with 106 additions and 83 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CHANGELOG.md
View File

@ -7,6 +7,14 @@ This project adheres to [Semantic Versioning](http://semver.org/).
## [0.13.0]
### Breaking semantic change
* The implementation of slicing with steps now matches the documentation.
Before, step identified the number to add to pass from one column/row index
to the next one. This made 0 step invalid. Now (and on the documentation so
far), the step is the number of ignored row/columns between each
row/column. Thus, a step of 0 means that no row/column is ignored. For
example, a step of, say, 3 on previous versions should now bet set to 2.
### Modified
* The trait `Axpy` takes one additional parameter for the type of `x`.
* The alias `MatrixNM` is now deprecated. Use `MatrixMN` instead (we

8
nalgebra-lapack/tests/linalg/cholesky.rs
View File

@ -41,8 +41,8 @@ quickcheck!{
let b1 = DVector::new_random(n);
let b2 = DMatrix::new_random(n, nb);
let sol1 = chol.solve(b1.clone()).unwrap();
let sol2 = chol.solve(b2.clone()).unwrap();
let sol1 = chol.solve(&b1).unwrap();
let sol2 = chol.solve(&b2).unwrap();
return relative_eq!(&m * sol1, b1, epsilon = 1.0e-6) &&
relative_eq!(&m * sol2, b2, epsilon = 1.0e-6)
@ -59,8 +59,8 @@ quickcheck!{
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = chol.solve(b1).unwrap();
let sol2 = chol.solve(b2).unwrap();
let sol1 = chol.solve(&b1).unwrap();
let sol2 = chol.solve(&b2).unwrap();
relative_eq!(m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m * sol2, b2, epsilon = 1.0e-7)

16
nalgebra-lapack/tests/linalg/lu.rs
View File

@ -45,11 +45,11 @@ quickcheck!{
let b1 = DVector::new_random(n);
let b2 = DMatrix::new_random(n, nb);
let sol1 = lup.solve(b1.clone()).unwrap();
let sol2 = lup.solve(b2.clone()).unwrap();
let sol1 = lup.solve(&b1).unwrap();
let sol2 = lup.solve(&b2).unwrap();
let tr_sol1 = lup.solve_transpose(b1.clone()).unwrap();
let tr_sol2 = lup.solve_transpose(b2.clone()).unwrap();
let tr_sol1 = lup.solve_transpose(&b1).unwrap();
let tr_sol2 = lup.solve_transpose(&b2).unwrap();
relative_eq!(&m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(&m * sol2, b2, epsilon = 1.0e-7) &&
@ -66,10 +66,10 @@ quickcheck!{
let b1 = Vector4::new_random();
let b2 = Matrix4x3::new_random();
let sol1 = lup.solve(b1).unwrap();
let sol2 = lup.solve(b2).unwrap();
let tr_sol1 = lup.solve_transpose(b1).unwrap();
let tr_sol2 = lup.solve_transpose(b2).unwrap();
let sol1 = lup.solve(&b1).unwrap();
let sol2 = lup.solve(&b2).unwrap();
let tr_sol1 = lup.solve_transpose(&b1).unwrap();
let tr_sol2 = lup.solve_transpose(&b2).unwrap();
relative_eq!(m * sol1, b1, epsilon = 1.0e-7) &&
relative_eq!(m * sol2, b2, epsilon = 1.0e-7) &&

8
src/core/matrix.rs
View File

@ -305,8 +305,12 @@ impl<N: Scalar, R: Dim, C: Dim, S: Storage<N, R, C>> Matrix<N, R, C, S> {
/// Returns a matrix containing the result of `f` applied to each entries of `self` and
/// `rhs`.
#[inline]
pub fn zip_map<F: FnMut(N, N) -> N>(&self, rhs: &Matrix<N, R, C, S>, mut f: F) -> MatrixMN<N, R, C>
where DefaultAllocator: Allocator<N, R, C> {
pub fn zip_map<N2, N3, S2, F>(&self, rhs: &Matrix<N2, R, C, S2>, mut f: F) -> MatrixMN<N3, R, C>
where N2: Scalar,
N3: Scalar,
S2: Storage<N2, R, C>,
F: FnMut(N, N2) -> N3,
DefaultAllocator: Allocator<N3, R, C> {
let (nrows, ncols) = self.data.shape();
let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) };

101
src/core/matrix_slice.rs
View File

@ -3,7 +3,7 @@ use std::ops::{Range, RangeFrom, RangeTo, RangeFull};
use std::slice;
use core::{Scalar, Matrix};
use core::dimension::{Dim, DimName, Dynamic, DimMul, DimProd, U1};
use core::dimension::{Dim, DimName, Dynamic, U1};
use core::iter::MatrixIter;
use core::storage::{Storage, StorageMut, Owned};
use core::allocator::Allocator;
@ -258,38 +258,36 @@ macro_rules! matrix_slice_impl(
$me.$rows_generic(first_row, Dynamic::new(nrows))
}
/// Extracts from this matrix a set of consecutive rows regularly spaced by `step` rows.
/// Extracts from this matrix a set of consecutive rows regularly skipping `step` rows.
#[inline]
pub fn $rows_with_step($me: $Me, first_row: usize, nrows: usize, step: usize)
-> $MatrixSlice<N, Dynamic, C, Dynamic, S::CStride> {
$me.$rows_generic_with_step(first_row, Dynamic::new(nrows), Dynamic::new(step))
$me.$rows_generic_with_step(first_row, Dynamic::new(nrows), step)
}
/// Extracts a compile-time number of consecutive rows from this matrix.
#[inline]
pub fn $fixed_rows<RSlice>($me: $Me, first_row: usize)
-> $MatrixSlice<N, RSlice, C, S::RStride, S::CStride>
where RSlice: DimName {
pub fn $fixed_rows<RSlice: DimName>($me: $Me, first_row: usize)
-> $MatrixSlice<N, RSlice, C, S::RStride, S::CStride> {
$me.$rows_generic(first_row, RSlice::name())
}
/// Extracts from this matrix a compile-time number of rows regularly spaced by `step` rows.
/// Extracts from this matrix a compile-time number of rows regularly skipping `step`
/// rows.
#[inline]
pub fn $fixed_rows_with_step<RSlice>($me: $Me, first_row: usize, step: usize)
-> $MatrixSlice<N, RSlice, C, Dynamic, S::CStride>
where RSlice: DimName {
pub fn $fixed_rows_with_step<RSlice: DimName>($me: $Me, first_row: usize, step: usize)
-> $MatrixSlice<N, RSlice, C, Dynamic, S::CStride> {
$me.$rows_generic_with_step(first_row, RSlice::name(), Dynamic::new(step))
$me.$rows_generic_with_step(first_row, RSlice::name(), step)
}
/// Extracts from this matrix `nrows` rows regularly spaced by `step` rows. Both argument may
/// or may not be values known at compile-time.
/// Extracts from this matrix `nrows` rows regularly skipping `step` rows. Both
/// argument may or may not be values known at compile-time.
#[inline]
pub fn $rows_generic<RSlice>($me: $Me, row_start: usize, nrows: RSlice)
-> $MatrixSlice<N, RSlice, C, S::RStride, S::CStride>
where RSlice: Dim {
pub fn $rows_generic<RSlice: Dim>($me: $Me, row_start: usize, nrows: RSlice)
-> $MatrixSlice<N, RSlice, C, S::RStride, S::CStride> {
let my_shape = $me.data.shape();
$me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (1, 1));
@ -302,19 +300,18 @@ macro_rules! matrix_slice_impl(
}
}
/// Extracts from this matrix `nrows` rows regularly spaced by `step` rows. Both argument may
/// or may not be values known at compile-time.
/// Extracts from this matrix `nrows` rows regularly skipping `step` rows. Both
/// argument may or may not be values known at compile-time.
#[inline]
pub fn $rows_generic_with_step<RSlice, RStep>($me: $Me, row_start: usize, nrows: RSlice, step: RStep)
-> $MatrixSlice<N, RSlice, C, DimProd<RStep, S::RStride>, S::CStride>
where RSlice: Dim,
RStep: DimMul<S::RStride> {
pub fn $rows_generic_with_step<RSlice>($me: $Me, row_start: usize, nrows: RSlice, step: usize)
-> $MatrixSlice<N, RSlice, C, Dynamic, S::CStride>
where RSlice: Dim {
let my_shape = $me.data.shape();
let my_strides = $me.data.strides();
$me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (step.value(), 1));
$me.assert_slice_index((row_start, 0), (nrows.value(), my_shape.1.value()), (step, 1));
let strides = (step.mul(my_strides.0), my_strides.1);
let strides = (Dynamic::new((step + 1) * my_strides.0.value()), my_strides.1);
let shape = (nrows, my_shape.1);
unsafe {
@ -348,39 +345,37 @@ macro_rules! matrix_slice_impl(
$me.$columns_generic(first_col, Dynamic::new(ncols))
}
/// Extracts from this matrix a set of consecutive columns regularly spaced by `step` columns.
/// Extracts from this matrix a set of consecutive columns regularly skipping `step`
/// columns.
#[inline]
pub fn $columns_with_step($me: $Me, first_col: usize, ncols: usize, step: usize)
-> $MatrixSlice<N, R, Dynamic, S::RStride, Dynamic> {
$me.$columns_generic_with_step(first_col, Dynamic::new(ncols), Dynamic::new(step))
$me.$columns_generic_with_step(first_col, Dynamic::new(ncols), step)
}
/// Extracts a compile-time number of consecutive columns from this matrix.
#[inline]
pub fn $fixed_columns<CSlice>($me: $Me, first_col: usize)
-> $MatrixSlice<N, R, CSlice, S::RStride, S::CStride>
where CSlice: DimName {
pub fn $fixed_columns<CSlice: DimName>($me: $Me, first_col: usize)
-> $MatrixSlice<N, R, CSlice, S::RStride, S::CStride> {
$me.$columns_generic(first_col, CSlice::name())
}
/// Extracts from this matrix a compile-time number of columns regularly spaced by `step`
/// columns.
/// Extracts from this matrix a compile-time number of columns regularly skipping
/// `step` columns.
#[inline]
pub fn $fixed_columns_with_step<CSlice>($me: $Me, first_col: usize, step: usize)
-> $MatrixSlice<N, R, CSlice, S::RStride, Dynamic>
where CSlice: DimName {
pub fn $fixed_columns_with_step<CSlice: DimName>($me: $Me, first_col: usize, step: usize)
-> $MatrixSlice<N, R, CSlice, S::RStride, Dynamic> {
$me.$columns_generic_with_step(first_col, CSlice::name(), Dynamic::new(step))
$me.$columns_generic_with_step(first_col, CSlice::name(), step)
}
/// Extracts from this matrix `ncols` columns. The number of columns may or may not be
/// known at compile-time.
#[inline]
pub fn $columns_generic<CSlice>($me: $Me, first_col: usize, ncols: CSlice)
-> $MatrixSlice<N, R, CSlice, S::RStride, S::CStride>
where CSlice: Dim {
pub fn $columns_generic<CSlice: Dim>($me: $Me, first_col: usize, ncols: CSlice)
-> $MatrixSlice<N, R, CSlice, S::RStride, S::CStride> {
let my_shape = $me.data.shape();
$me.assert_slice_index((0, first_col), (my_shape.0.value(), ncols.value()), (1, 1));
@ -393,20 +388,18 @@ macro_rules! matrix_slice_impl(
}
/// Extracts from this matrix `ncols` columns regularly spaced by `step` columns. Both argument may
/// Extracts from this matrix `ncols` columns skipping `step` columns. Both argument may
/// or may not be values known at compile-time.
#[inline]
pub fn $columns_generic_with_step<CSlice, CStep>($me: $Me, first_col: usize, ncols: CSlice, step: CStep)
-> $MatrixSlice<N, R, CSlice, S::RStride, DimProd<CStep, S::CStride>>
where CSlice: Dim,
CStep: DimMul<S::CStride> {
pub fn $columns_generic_with_step<CSlice: Dim>($me: $Me, first_col: usize, ncols: CSlice, step: usize)
-> $MatrixSlice<N, R, CSlice, S::RStride, Dynamic> {
let my_shape = $me.data.shape();
let my_strides = $me.data.strides();
$me.assert_slice_index((0, first_col), (my_shape.0.value(), ncols.value()), (1, step.value()));
$me.assert_slice_index((0, first_col), (my_shape.0.value(), ncols.value()), (1, step));
let strides = (my_strides.0, step.mul(my_strides.1));
let strides = (my_strides.0, Dynamic::new((step + 1) * my_strides.1.value()));
let shape = (my_shape.0, ncols);
unsafe {
@ -444,7 +437,6 @@ macro_rules! matrix_slice_impl(
pub fn $slice_with_steps($me: $Me, start: (usize, usize), shape: (usize, usize), steps: (usize, usize))
-> $MatrixSlice<N, Dynamic, Dynamic, Dynamic, Dynamic> {
let shape = (Dynamic::new(shape.0), Dynamic::new(shape.1));
let steps = (Dynamic::new(steps.0), Dynamic::new(steps.1));
$me.$generic_slice_with_steps(start, shape, steps)
}
@ -476,7 +468,6 @@ macro_rules! matrix_slice_impl(
where RSlice: DimName,
CSlice: DimName {
let shape = (RSlice::name(), CSlice::name());
let steps = (Dynamic::new(steps.0), Dynamic::new(steps.1));
$me.$generic_slice_with_steps(start, shape, steps)
}
@ -497,21 +488,19 @@ macro_rules! matrix_slice_impl(
/// Creates a slice that may or may not have a fixed size and stride.
#[inline]
pub fn $generic_slice_with_steps<RSlice, CSlice, RStep, CStep>($me: $Me,
pub fn $generic_slice_with_steps<RSlice, CSlice>($me: $Me,
start: (usize, usize),
shape: (RSlice, CSlice),
steps: (RStep, CStep))
-> $MatrixSlice<N, RSlice, CSlice, DimProd<RStep, S::RStride>, DimProd<CStep, S::CStride>>
steps: (usize, usize))
-> $MatrixSlice<N, RSlice, CSlice, Dynamic, Dynamic>
where RSlice: Dim,
CSlice: Dim,
RStep: DimMul<S::RStride>,
CStep: DimMul<S::CStride> {
CSlice: Dim {
assert!(steps.0.value() > 0 && steps.1.value() > 0, "Matrix slicing steps must not be zero.");
$me.assert_slice_index(start, (shape.0.value(), shape.1.value()), (steps.0.value(), steps.1.value()));
$me.assert_slice_index(start, (shape.0.value(), shape.1.value()), steps);
let my_strides = $me.data.strides();
let strides = (steps.0.mul(my_strides.0), steps.1.mul(my_strides.1));
let strides = (Dynamic::new((steps.0 + 1) * my_strides.0.value()),
Dynamic::new((steps.1 + 1) * my_strides.1.value()));
unsafe {
let data = $SliceStorage::new_with_strides_unchecked($data, start, shape, strides);

8
src/linalg/svd.rs
View File

@ -344,11 +344,11 @@ impl<N: Real, R: DimMin<C>, C: Dim> SVD<N, R, C>
if b.is_upper_diagonal() {
if let Some(ref mut u) = *u {
rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::<U2>(i, k - i + 1));
rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::<U2>(i, k - i));
}
}
else if let Some(ref mut v_t) = *v_t {
rot.rotate(&mut v_t.fixed_rows_with_step_mut::<U2>(i, k - i + 1));
rot.rotate(&mut v_t.fixed_rows_with_step_mut::<U2>(i, k - i));
}
if k + 1 != end {
@ -377,11 +377,11 @@ impl<N: Real, R: DimMin<C>, C: Dim> SVD<N, R, C>
if b.is_upper_diagonal() {
if let Some(ref mut v_t) = *v_t {
rot.rotate(&mut v_t.fixed_rows_with_step_mut::<U2>(k, i + 1 - k));
rot.rotate(&mut v_t.fixed_rows_with_step_mut::<U2>(k, i - k));
}
}
else if let Some(ref mut u) = *u {
rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::<U2>(k, i + 1 - k));
rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::<U2>(k, i - k));
}
if k > 0 {

22
tests/core/matrix.rs
View File

@ -434,6 +434,28 @@ fn map() {
assert_eq!(computed, expected);
}
#[test]
fn zip_map() {
let a = Matrix3::new(
11i32, 12, 13,
21, 22, 23,
31, 32, 33);
let b = Matrix3::new(
11u32, 12, 13,
21, 22, 23,
31, 32, 33);
let expected = Matrix3::new(
22.0f32, 24.0, 26.0,
42.0, 44.0, 46.0,
62.0, 64.0, 66.0);
let computed = a.zip_map(&b, |ea, eb| ea as f32 + eb as f32);
assert_eq!(computed, expected);
}
#[test]
#[should_panic]
fn trace_panic() {

14
tests/core/matrix_slice.rs
View File

@ -13,7 +13,7 @@ fn nested_fixed_slices() {
let s1 = a.fixed_slice::<U3, U3>(0, 1); // Simple slice.
let s2 = s1.fixed_slice::<U2, U2>(1, 1); // Slice of slice.
let s3 = s1.fixed_slice_with_steps::<U2, U2>((0, 0), (2, 2)); // Slice of slice with steps.
let s3 = s1.fixed_slice_with_steps::<U2, U2>((0, 0), (1, 1)); // Slice of slice with steps.
let expected_owned_s1 = Matrix3::new(12.0, 13.0, 14.0,
22.0, 23.0, 24.0,
@ -38,7 +38,7 @@ fn nested_slices() {
let s1 = a.slice((0, 1), (3, 3));
let s2 = s1.slice((1, 1), (2, 2));
let s3 = s1.slice_with_steps((0, 0), (2, 2), (2, 2));
let s3 = s1.slice_with_steps((0, 0), (2, 2), (1, 1));
let expected_owned_s1 = DMatrix::from_row_slice(3, 3, &[ 12.0, 13.0, 14.0,
22.0, 23.0, 24.0,
@ -63,7 +63,7 @@ fn slice_mut() {
{
// We modify `a` through the mutable slice.
let mut s1 = a.slice_with_steps_mut((0, 1), (2, 2), (2, 2));
let mut s1 = a.slice_with_steps_mut((0, 1), (2, 2), (1, 1));
s1.fill(0.0);
}
@ -83,7 +83,7 @@ fn nested_row_slices() {
51.0, 52.0,
61.0, 62.0);
let s1 = a.fixed_rows::<U4>(1);
let s2 = s1.fixed_rows_with_step::<U2>(1, 2);
let s2 = s1.fixed_rows_with_step::<U2>(1, 1);
let expected_owned_s1 = Matrix4x2::new(21.0, 22.0,
31.0, 32.0,
@ -107,7 +107,7 @@ fn row_slice_mut() {
61.0, 62.0);
{
// We modify `a` through the mutable slice.
let mut s1 = a.rows_with_step_mut(1, 3, 2);
let mut s1 = a.rows_with_step_mut(1, 3, 1);
s1.fill(0.0);
}
@ -126,7 +126,7 @@ fn nested_col_slices() {
let a = Matrix2x6::new(11.0, 12.0, 13.0, 14.0, 15.0, 16.0,
21.0, 22.0, 23.0, 24.0, 25.0, 26.0);
let s1 = a.fixed_columns::<U4>(1);
let s2 = s1.fixed_columns_with_step::<U2>(1, 2);
let s2 = s1.fixed_columns_with_step::<U2>(1, 1);
let expected_owned_s1 = Matrix2x4::new(12.0, 13.0, 14.0, 15.0,
22.0, 23.0, 24.0, 25.0);
@ -145,7 +145,7 @@ fn col_slice_mut() {
{
// We modify `a` through the mutable slice.
let mut s1 = a.columns_with_step_mut(1, 3, 2);
let mut s1 = a.columns_with_step_mut(1, 3, 1);
s1.fill(0.0);
}