From c258d13f98ad2eea96ea7982e4d53d1338c624d3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?S=C3=A9bastien=20Crozet?= Date: Sun, 13 Aug 2017 19:52:55 +0200 Subject: [PATCH] Fix stepping for slices. The previous implementation was wrong compared to what the documentatino claimed. --- CHANGELOG.md | 8 ++ nalgebra-lapack/tests/linalg/cholesky.rs | 8 +- nalgebra-lapack/tests/linalg/lu.rs | 16 ++-- src/core/matrix.rs | 8 +- src/core/matrix_slice.rs | 105 ++++++++++------------- src/linalg/svd.rs | 8 +- tests/core/matrix.rs | 22 +++++ tests/core/matrix_slice.rs | 14 +-- 8 files changed, 106 insertions(+), 83 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d04a3bdd..eebda429 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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 diff --git a/nalgebra-lapack/tests/linalg/cholesky.rs b/nalgebra-lapack/tests/linalg/cholesky.rs index 9f945c59..c0102af5 100644 --- a/nalgebra-lapack/tests/linalg/cholesky.rs +++ b/nalgebra-lapack/tests/linalg/cholesky.rs @@ -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) diff --git a/nalgebra-lapack/tests/linalg/lu.rs b/nalgebra-lapack/tests/linalg/lu.rs index 667bc860..ff36735f 100644 --- a/nalgebra-lapack/tests/linalg/lu.rs +++ b/nalgebra-lapack/tests/linalg/lu.rs @@ -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) && diff --git a/src/core/matrix.rs b/src/core/matrix.rs index eef100e8..fe14cc4a 100644 --- a/src/core/matrix.rs +++ b/src/core/matrix.rs @@ -305,8 +305,12 @@ impl> Matrix { /// Returns a matrix containing the result of `f` applied to each entries of `self` and /// `rhs`. #[inline] - pub fn zip_map N>(&self, rhs: &Matrix, mut f: F) -> MatrixMN - where DefaultAllocator: Allocator { + pub fn zip_map(&self, rhs: &Matrix, mut f: F) -> MatrixMN + where N2: Scalar, + N3: Scalar, + S2: Storage, + F: FnMut(N, N2) -> N3, + DefaultAllocator: Allocator { let (nrows, ncols) = self.data.shape(); let mut res = unsafe { MatrixMN::new_uninitialized_generic(nrows, ncols) }; diff --git a/src/core/matrix_slice.rs b/src/core/matrix_slice.rs index ec47615c..4fa518d9 100644 --- a/src/core/matrix_slice.rs +++ b/src/core/matrix_slice.rs @@ -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 { - $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($me: $Me, first_row: usize) - -> $MatrixSlice - where RSlice: DimName { + pub fn $fixed_rows($me: $Me, first_row: usize) + -> $MatrixSlice { $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($me: $Me, first_row: usize, step: usize) - -> $MatrixSlice - where RSlice: DimName { + pub fn $fixed_rows_with_step($me: $Me, first_row: usize, step: usize) + -> $MatrixSlice { - $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($me: $Me, row_start: usize, nrows: RSlice) - -> $MatrixSlice - where RSlice: Dim { + pub fn $rows_generic($me: $Me, row_start: usize, nrows: RSlice) + -> $MatrixSlice { 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($me: $Me, row_start: usize, nrows: RSlice, step: RStep) - -> $MatrixSlice, S::CStride> - where RSlice: Dim, - RStep: DimMul { + pub fn $rows_generic_with_step($me: $Me, row_start: usize, nrows: RSlice, step: usize) + -> $MatrixSlice + 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 { - $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($me: $Me, first_col: usize) - -> $MatrixSlice - where CSlice: DimName { + pub fn $fixed_columns($me: $Me, first_col: usize) + -> $MatrixSlice { $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($me: $Me, first_col: usize, step: usize) - -> $MatrixSlice - where CSlice: DimName { + pub fn $fixed_columns_with_step($me: $Me, first_col: usize, step: usize) + -> $MatrixSlice { - $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($me: $Me, first_col: usize, ncols: CSlice) - -> $MatrixSlice - where CSlice: Dim { + pub fn $columns_generic($me: $Me, first_col: usize, ncols: CSlice) + -> $MatrixSlice { 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($me: $Me, first_col: usize, ncols: CSlice, step: CStep) - -> $MatrixSlice> - where CSlice: Dim, - CStep: DimMul { + pub fn $columns_generic_with_step($me: $Me, first_col: usize, ncols: CSlice, step: usize) + -> $MatrixSlice { 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 { 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($me: $Me, - start: (usize, usize), - shape: (RSlice, CSlice), - steps: (RStep, CStep)) - -> $MatrixSlice, DimProd> + pub fn $generic_slice_with_steps($me: $Me, + start: (usize, usize), + shape: (RSlice, CSlice), + steps: (usize, usize)) + -> $MatrixSlice where RSlice: Dim, - CSlice: Dim, - RStep: DimMul, - CStep: DimMul { + 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); diff --git a/src/linalg/svd.rs b/src/linalg/svd.rs index 5489840d..fd4bd0dc 100644 --- a/src/linalg/svd.rs +++ b/src/linalg/svd.rs @@ -344,11 +344,11 @@ impl, C: Dim> SVD if b.is_upper_diagonal() { if let Some(ref mut u) = *u { - rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::(i, k - i + 1)); + rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::(i, k - i)); } } else if let Some(ref mut v_t) = *v_t { - rot.rotate(&mut v_t.fixed_rows_with_step_mut::(i, k - i + 1)); + rot.rotate(&mut v_t.fixed_rows_with_step_mut::(i, k - i)); } if k + 1 != end { @@ -377,11 +377,11 @@ impl, C: Dim> SVD if b.is_upper_diagonal() { if let Some(ref mut v_t) = *v_t { - rot.rotate(&mut v_t.fixed_rows_with_step_mut::(k, i + 1 - k)); + rot.rotate(&mut v_t.fixed_rows_with_step_mut::(k, i - k)); } } else if let Some(ref mut u) = *u { - rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::(k, i + 1 - k)); + rot.inverse().rotate_rows(&mut u.fixed_columns_with_step_mut::(k, i - k)); } if k > 0 { diff --git a/tests/core/matrix.rs b/tests/core/matrix.rs index 6ef8722a..e4557a30 100644 --- a/tests/core/matrix.rs +++ b/tests/core/matrix.rs @@ -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() { diff --git a/tests/core/matrix_slice.rs b/tests/core/matrix_slice.rs index 516b5668..9313a3db 100644 --- a/tests/core/matrix_slice.rs +++ b/tests/core/matrix_slice.rs @@ -13,7 +13,7 @@ fn nested_fixed_slices() { let s1 = a.fixed_slice::(0, 1); // Simple slice. let s2 = s1.fixed_slice::(1, 1); // Slice of slice. - let s3 = s1.fixed_slice_with_steps::((0, 0), (2, 2)); // Slice of slice with steps. + let s3 = s1.fixed_slice_with_steps::((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::(1); - let s2 = s1.fixed_rows_with_step::(1, 2); + let s2 = s1.fixed_rows_with_step::(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::(1); - let s2 = s1.fixed_columns_with_step::(1, 2); + let s2 = s1.fixed_columns_with_step::(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); }