FEAT: Fix memory continuity judgment when stride is negative

rust-ndarray · bluss · Jan 13, 2021 · Jan 6, 2021 · Jan 8, 2021 · Jan 11, 2021
commit f0dafb39bf48c1d3a9ab40299d24993b4889699a
diff --git a/blas-tests/tests/oper.rs b/blas-tests/tests/oper.rs
@@ -173,7 +173,7 @@ where
     S2: Data<Elem = A>,
 {
     let ((m, _), k) = (lhs.dim(), rhs.dim());
-    reference_mat_mul(lhs, &rhs.to_owned().into_shape((k, 1)).unwrap())
+    reference_mat_mul(lhs, &rhs.as_standard_layout().into_shape((k, 1)).unwrap())
         .into_shape(m)
         .unwrap()
 }
@@ -186,7 +186,7 @@ where
     S2: Data<Elem = A>,
 {
     let (m, (_, n)) = (lhs.dim(), rhs.dim());
-    reference_mat_mul(&lhs.to_owned().into_shape((1, m)).unwrap(), rhs)
+    reference_mat_mul(&lhs.as_standard_layout().into_shape((1, m)).unwrap(), rhs)
         .into_shape(n)
         .unwrap()
 }

diff --git a/src/dimension/dimension_trait.rs b/src/dimension/dimension_trait.rs
@@ -286,17 +286,16 @@ pub trait Dimension:
             return true;
         }
         if dim.ndim() == 1 {
-            return false;
+            return strides[0] as isize == -1;
         }
         let order = strides._fastest_varying_stride_order();
         let strides = strides.slice();
 
-        // FIXME: Negative strides
         let dim_slice = dim.slice();
         let mut cstride = 1;
         for &i in order.slice() {
             // a dimension of length 1 can have unequal strides
-            if dim_slice[i] != 1 && strides[i] != cstride {
+            if dim_slice[i] != 1 && (strides[i] as isize).abs() as usize != cstride {
                 return false;
             }
             cstride *= dim_slice[i];
@@ -307,16 +306,16 @@ pub trait Dimension:
     /// Return the axis ordering corresponding to the fastest variation
     /// (in ascending order).
     ///
-    /// Assumes that no stride value appears twice. This cannot yield the correct
-    /// result the strides are not positive.
+    /// Assumes that no stride value appears twice.
     #[doc(hidden)]
     fn _fastest_varying_stride_order(&self) -> Self {
         let mut indices = self.clone();
         for (i, elt) in enumerate(indices.slice_mut()) {
             *elt = i;
         }
         let strides = self.slice();
-        indices.slice_mut().sort_by_key(|&i| strides[i]);
+        indices.slice_mut()
+            .sort_by_key(|&i| (strides[i] as isize).abs());
         indices
     }
 
@@ -645,7 +644,7 @@ impl Dimension for Dim<[Ix; 2]> {
 
     #[inline]
     fn _fastest_varying_stride_order(&self) -> Self {
-        if get!(self, 0) as Ixs <= get!(self, 1) as Ixs {
+        if (get!(self, 0) as Ixs).abs() <= (get!(self, 1) as Ixs).abs() {
             Ix2(0, 1)
         } else {
             Ix2(1, 0)
@@ -805,7 +804,7 @@ impl Dimension for Dim<[Ix; 3]> {
         let mut order = Ix3(0, 1, 2);
         macro_rules! swap {
             ($stride:expr, $order:expr, $x:expr, $y:expr) => {
-                if $stride[$x] > $stride[$y] {
+                if ($stride[$x] as isize).abs() > ($stride[$y] as isize).abs() {
                     $stride.swap($x, $y);
                     $order.ixm().swap($x, $y);
                 }

diff --git a/src/dimension/mod.rs b/src/dimension/mod.rs
@@ -46,15 +46,12 @@ pub fn stride_offset(n: Ix, stride: Ix) -> isize {
 /// There is overlap if, when iterating through the dimensions in order of
 /// increasing stride, the current stride is less than or equal to the maximum
 /// possible offset along the preceding axes. (Axes of length ≤1 are ignored.)
-///
-/// The current implementation assumes that strides of axes with length > 1 are
-/// nonnegative. Additionally, it does not check for overflow.
 pub fn dim_stride_overlap<D: Dimension>(dim: &D, strides: &D) -> bool {
     let order = strides._fastest_varying_stride_order();
     let mut sum_prev_offsets = 0;
     for &index in order.slice() {
         let d = dim[index];
-        let s = strides[index] as isize;
+        let s = (strides[index] as isize).abs();
         match d {
             0 => return false,
             1 => {}
@@ -210,8 +207,7 @@ where
 ///
 /// 2. The product of non-zero axis lengths must not exceed `isize::MAX`.
 ///
-/// 3. For axes with length > 1, the stride must be nonnegative. This is
-///    necessary to make sure the pointer cannot move backwards outside the
+/// 3. For axes with length > 1, the pointer cannot move outside the
 ///    slice. For axes with length ≤ 1, the stride can be anything.
 ///
 /// 4. If the array will be empty (any axes are zero-length), the difference
@@ -257,14 +253,6 @@ fn can_index_slice_impl<D: Dimension>(
         return Err(from_kind(ErrorKind::OutOfBounds));
     }
 
-    // Check condition 3.
-    for (&d, &s) in izip!(dim.slice(), strides.slice()) {
-        let s = s as isize;
-        if d > 1 && s < 0 {
-            return Err(from_kind(ErrorKind::Unsupported));
-        }
-    }
-
     // Check condition 5.
     if !is_empty && dim_stride_overlap(dim, strides) {
         return Err(from_kind(ErrorKind::Unsupported));
@@ -394,6 +382,19 @@ fn to_abs_slice(axis_len: usize, slice: Slice) -> (usize, usize, isize) {
     (start, end, step)
 }
 
+/// This function computes the offset from the logically first element to the first element in
+/// memory of the array. The result is always <= 0.
+pub fn offset_from_ptr_to_memory(dim: &[Ix], strides: &[Ix]) -> isize {
+    let offset = izip!(dim, strides).fold(0, |_offset, (d, s)| {
+        if (*s as isize) < 0 {
+            _offset + *s as isize * (*d as isize - 1)
+        } else {
+            _offset
+        }
+    });
+    offset
+}
+
 /// Modify dimension, stride and return data pointer offset
 ///
 /// **Panics** if stride is 0 or if any index is out of bounds.
@@ -693,13 +694,21 @@ mod test {
         let dim = (2, 3, 2).into_dimension();
         let strides = (5, 2, 1).into_dimension();
         assert!(super::dim_stride_overlap(&dim, &strides));
+        let strides = (-5isize as usize, 2, -1isize as usize).into_dimension();
+        assert!(super::dim_stride_overlap(&dim, &strides));
         let strides = (6, 2, 1).into_dimension();
         assert!(!super::dim_stride_overlap(&dim, &strides));
+        let strides = (6, -2isize as usize, 1).into_dimension();
+        assert!(!super::dim_stride_overlap(&dim, &strides));
         let strides = (6, 0, 1).into_dimension();
         assert!(super::dim_stride_overlap(&dim, &strides));
+        let strides = (-6isize as usize, 0, 1).into_dimension();
+        assert!(super::dim_stride_overlap(&dim, &strides));
         let dim = (2, 2).into_dimension();
         let strides = (3, 2).into_dimension();
         assert!(!super::dim_stride_overlap(&dim, &strides));
+        let strides = (3, -2isize as usize).into_dimension();
+        assert!(!super::dim_stride_overlap(&dim, &strides));
     }
 
     #[test]
@@ -736,7 +745,7 @@ mod test {
         can_index_slice::<i32, _>(&[1], &Ix1(2), &Ix1(1)).unwrap_err();
         can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(0)).unwrap_err();
         can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(1)).unwrap();
-        can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(-1isize as usize)).unwrap_err();
+        can_index_slice::<i32, _>(&[1, 2], &Ix1(2), &Ix1(-1isize as usize)).unwrap();
     }
 
     #[test]

diff --git a/src/impl_constructors.rs b/src/impl_constructors.rs
@@ -16,6 +16,7 @@ use num_traits::{Float, One, Zero};
 use std::mem::MaybeUninit;
 
 use crate::dimension;
+use crate::dimension::offset_from_ptr_to_memory;
 use crate::error::{self, ShapeError};
 use crate::extension::nonnull::nonnull_from_vec_data;
 use crate::imp_prelude::*;
@@ -24,6 +25,7 @@ use crate::indices;
 use crate::iterators::{to_vec, to_vec_mapped};
 use crate::StrideShape;
 use crate::{geomspace, linspace, logspace};
+use rawpointer::PointerExt;
 
 /// # Constructor Methods for Owned Arrays
 ///
@@ -431,7 +433,8 @@ where
     ///
     /// 2. The product of non-zero axis lengths must not exceed `isize::MAX`.
     ///
-    /// 3. For axes with length > 1, the stride must be nonnegative.
+    /// 3. For axes with length > 1, the pointer cannot move outside the
+    ///    slice.
     ///
     /// 4. If the array will be empty (any axes are zero-length), the
     ///    difference between the least address and greatest address accessible
@@ -457,7 +460,8 @@ where
         // debug check for issues that indicates wrong use of this constructor
         debug_assert!(dimension::can_index_slice(&v, &dim, &strides).is_ok());
         ArrayBase {
-            ptr: nonnull_from_vec_data(&mut v),
+            ptr: nonnull_from_vec_data(&mut v)
+                .offset(offset_from_ptr_to_memory(dim.slice(), strides.slice()).abs()),
             data: DataOwned::new(v),
             strides,
             dim,
@@ -483,7 +487,7 @@ where
     ///
     /// This constructor is limited to elements where `A: Copy` (no destructors)
     /// to avoid users shooting themselves too hard in the foot.
-    /// 
+    ///
     /// (Also note that the constructors `from_shape_vec` and
     /// `from_shape_vec_unchecked` allow the user yet more control, in the sense
     /// that Arrays can be created from arbitrary vectors.)

diff --git a/src/impl_methods.rs b/src/impl_methods.rs
@@ -17,7 +17,8 @@ use crate::arraytraits;
 use crate::dimension;
 use crate::dimension::IntoDimension;
 use crate::dimension::{
-    abs_index, axes_of, do_slice, merge_axes, size_of_shape_checked, stride_offset, Axes,
+    abs_index, axes_of, do_slice, merge_axes, offset_from_ptr_to_memory, size_of_shape_checked,
+    stride_offset, Axes,
 };
 use crate::error::{self, ErrorKind, ShapeError};
 use crate::itertools::zip;
@@ -1280,9 +1281,6 @@ where
     }
 
     /// Return true if the array is known to be contiguous.
-    ///
-    /// Will detect c- and f-contig arrays correctly, but otherwise
-    /// There are some false negatives.
     pub(crate) fn is_contiguous(&self) -> bool {
         D::is_contiguous(&self.dim, &self.strides)
     }
@@ -1404,14 +1402,18 @@ where
     ///
     /// If this function returns `Some(_)`, then the elements in the slice
     /// have whatever order the elements have in memory.
-    ///
-    /// Implementation notes: Does not yet support negatively strided arrays.
     pub fn as_slice_memory_order(&self) -> Option<&[A]>
     where
         S: Data,
     {
         if self.is_contiguous() {
-            unsafe { Some(slice::from_raw_parts(self.ptr.as_ptr(), self.len())) }
+            let offset = offset_from_ptr_to_memory(self.dim.slice(), self.strides.slice());
+            unsafe {
+                Some(slice::from_raw_parts(
+                    self.ptr.offset(offset).as_ptr(),
+                    self.len(),
+                ))
+            }
         } else {
             None
         }
@@ -1425,7 +1427,13 @@ where
     {
         if self.is_contiguous() {
             self.ensure_unique();
-            unsafe { Some(slice::from_raw_parts_mut(self.ptr.as_ptr(), self.len())) }
+            let offset = offset_from_ptr_to_memory(self.dim.slice(), self.strides.slice());
+            unsafe {
+                Some(slice::from_raw_parts_mut(
+                    self.ptr.offset(offset).as_ptr(),
+                    self.len(),
+                ))
+            }
         } else {
             None
         }

diff --git a/tests/dimension.rs b/tests/dimension.rs
@@ -118,19 +118,36 @@ fn fastest_varying_order() {
     let order = strides._fastest_varying_stride_order();
     assert_eq!(order.slice(), &[3, 0, 2, 1]);
 
+    let strides = Dim([-2isize as usize, 8, -4isize as usize, -1isize as usize]);
+    let order = strides._fastest_varying_stride_order();
+    assert_eq!(order.slice(), &[3, 0, 2, 1]);
+
     assert_eq!(Dim([1, 3])._fastest_varying_stride_order(), Dim([0, 1]));
+    assert_eq!(
+        Dim([1, -3isize as usize])._fastest_varying_stride_order(),
+        Dim([0, 1])
+    );
     assert_eq!(Dim([7, 2])._fastest_varying_stride_order(), Dim([1, 0]));
+    assert_eq!(
+        Dim([-7isize as usize, 2])._fastest_varying_stride_order(),
+        Dim([1, 0])
+    );
     assert_eq!(
         Dim([6, 1, 3])._fastest_varying_stride_order(),
         Dim([1, 2, 0])
     );
+    assert_eq!(
+        Dim([-6isize as usize, 1, -3isize as usize])._fastest_varying_stride_order(),
+        Dim([1, 2, 0])
+    );
 
     // it's important that it produces distinct indices. Prefer the stable order
     // where 0 is before 1 when they are equal.
     assert_eq!(Dim([2, 2])._fastest_varying_stride_order(), [0, 1]);
     assert_eq!(Dim([2, 2, 1])._fastest_varying_stride_order(), [2, 0, 1]);
     assert_eq!(
-        Dim([2, 2, 3, 1, 2])._fastest_varying_stride_order(),
+        Dim([-2isize as usize, -2isize as usize, 3, 1, -2isize as usize])
+            ._fastest_varying_stride_order(),
         [3, 0, 1, 4, 2]
     );
 }

diff --git a/tests/oper.rs b/tests/oper.rs
@@ -707,7 +707,7 @@ fn gen_mat_vec_mul() {
         S2: Data<Elem = A>,
     {
         let ((m, _), k) = (lhs.dim(), rhs.dim());
-        reference_mat_mul(lhs, &rhs.to_owned().into_shape((k, 1)).unwrap())
+        reference_mat_mul(lhs, &rhs.as_standard_layout().into_shape((k, 1)).unwrap())
             .into_shape(m)
             .unwrap()
     }
@@ -772,7 +772,7 @@ fn vec_mat_mul() {
         S2: Data<Elem = A>,
     {
         let (m, (_, n)) = (lhs.dim(), rhs.dim());
-        reference_mat_mul(&lhs.to_owned().into_shape((1, m)).unwrap(), rhs)
+        reference_mat_mul(&lhs.as_standard_layout().into_shape((1, m)).unwrap(), rhs)
             .into_shape(n)
             .unwrap()
     }