Add .as_slice_memory_order(), improve scalar_sum, and fix bugs in from_vec_dim_stride

benches/bench1.rs

@@ -48,7 +48,7 @@ fn sum_1d_raw(bench: &mut test::Bencher)
     let a = black_box(a);
     bench.iter(|| {
         let mut sum = 0;
-        for &elt in a.raw_data() {
+        for &elt in a.as_slice_memory_order().unwrap() {
             sum += elt;
         }
         sum
@@ -93,7 +93,7 @@ fn sum_2d_raw(bench: &mut test::Bencher)
     let a = black_box(a);
     bench.iter(|| {
         let mut sum = 0;
-        for &elt in a.raw_data() {
+        for &elt in a.as_slice_memory_order().unwrap() {
             sum += elt;
         }
         sum
@@ -373,12 +373,12 @@ fn muladd_2d_f32_blas(bench: &mut test::Bencher)
 }
 
 #[bench]
-fn assign_scalar_2d_large(bench: &mut test::Bencher)
+fn assign_scalar_2d_corder(bench: &mut test::Bencher)
 {
     let a = OwnedArray::zeros((64, 64));
     let mut a = black_box(a);
     let s = 3.;
-    bench.iter(|| a.assign_scalar(&s))
+    bench.iter(move || a.assign_scalar(&s))
 }
 
 #[bench]
@@ -388,26 +388,43 @@ fn assign_scalar_2d_cutout(bench: &mut test::Bencher)
     let a = a.slice_mut(s![1..-1, 1..-1]);
     let mut a = black_box(a);
     let s = 3.;
-    bench.iter(|| a.assign_scalar(&s))
+    bench.iter(move || a.assign_scalar(&s))
 }
 
 #[bench]
-fn assign_scalar_2d_transposed_large(bench: &mut test::Bencher)
+fn assign_scalar_2d_forder(bench: &mut test::Bencher)
 {
     let mut a = OwnedArray::zeros((64, 64));
     a.swap_axes(0, 1);
     let mut a = black_box(a);
     let s = 3.;
-    bench.iter(|| a.assign_scalar(&s))
+    bench.iter(move || a.assign_scalar(&s))
 }
 
 #[bench]
-fn assign_scalar_2d_raw_large(bench: &mut test::Bencher)
+fn assign_zero_2d_corder(bench: &mut test::Bencher)
 {
     let a = OwnedArray::zeros((64, 64));
     let mut a = black_box(a);
-    let s = 3.;
-    bench.iter(|| for elt in a.raw_data_mut() { *elt = s; });
+    bench.iter(|| a.assign_scalar(&0.))
+}
+
+#[bench]
+fn assign_zero_2d_cutout(bench: &mut test::Bencher)
+{
+    let mut a = OwnedArray::zeros((66, 66));
+    let a = a.slice_mut(s![1..-1, 1..-1]);
+    let mut a = black_box(a);
+    bench.iter(|| a.assign_scalar(&0.))
+}
+
+#[bench]
+fn assign_zero_2d_forder(bench: &mut test::Bencher)
+{
+    let mut a = OwnedArray::zeros((64, 64));
+    a.swap_axes(0, 1);
+    let mut a = black_box(a);
+    bench.iter(|| a.assign_scalar(&0.))
 }
 
 #[bench]

src/dimension.rs

@@ -8,6 +8,7 @@
 use std::cmp::Ordering;
 use std::fmt::Debug;
 use std::slice;
+use itertools::free::enumerate;
 
 use super::{Si, Ix, Ixs};
 use super::zipsl;
@@ -19,30 +20,6 @@ pub fn stride_offset(n: Ix, stride: Ix) -> isize {
     (n as isize) * ((stride as Ixs) as isize)
 }
 
-/// Check whether `stride` is strictly positive
-#[inline]
-fn stride_is_positive(stride: Ix) -> bool {
-    (stride as Ixs) > 0
-}
-
-/// Return the axis ordering corresponding to the fastest variation
-///
-/// Assumes that no stride value appears twice. This cannot yield the correct
-/// result the strides are not positive.
-fn fastest_varying_order<D: Dimension>(strides: &D) -> D {
-    let mut sorted = strides.clone();
-    sorted.slice_mut().sort();
-    let mut res = strides.clone();
-    for (index, &val) in strides.slice().iter().enumerate() {
-        let sorted_ind = sorted.slice()
-                               .iter()
-                               .position(|&x| x == val)
-                               .unwrap(); // cannot panic by construction
-        res.slice_mut()[sorted_ind] = index;
-    }
-    res
-}
-
 /// Check whether the given `dim` and `stride` lead to overlapping indices
 ///
 /// There is overlap if, when iterating through the dimensions in the order
@@ -51,15 +28,19 @@ fn fastest_varying_order<D: Dimension>(strides: &D) -> D {
 ///
 /// The current implementation assumes strides to be positive
 pub fn dim_stride_overlap<D: Dimension>(dim: &D, strides: &D) -> bool {
-    let order = fastest_varying_order(strides);
+    let order = strides._fastest_varying_stride_order();
 
+    let dim = dim.slice();
+    let strides = strides.slice();
     let mut prev_offset = 1;
-    for &index in order.slice().iter() {
-        let s = strides.slice()[index];
-        if (s as isize) < prev_offset {
+    for &index in order.slice() {
+        let d = dim[index];
+        let s = strides[index];
+        // any stride is ok if dimension is 1
+        if d != 1 && (s as isize) < prev_offset {
             return true;
         }
-        prev_offset = stride_offset(dim.slice()[index], s);
+        prev_offset = stride_offset(d, s);
     }
     false
 }
@@ -74,33 +55,42 @@ pub fn dim_stride_overlap<D: Dimension>(dim: &D, strides: &D) -> bool {
 pub fn can_index_slice<A, D: Dimension>(data: &[A], dim: &D, strides: &D)
     -> Result<(), ShapeError>
 {
-    if strides.slice().iter().cloned().all(stride_is_positive) {
-        if dim.size_checked().is_none() {
-            return Err(from_kind(ErrorKind::OutOfBounds));
+    // check lengths of axes.
+    let len = match dim.size_checked() {
+        Some(l) => l,
+        None => return Err(from_kind(ErrorKind::OutOfBounds)),
+    };
+    // check if strides are strictly positive (zero ok for len 0)
+    for &s in strides.slice() {
+        let s = s as Ixs;
+        if s < 1 && (len != 0 || s < 0) {
+            return Err(from_kind(ErrorKind::Unsupported));
         }
-        let mut last_index = dim.clone();
-        for mut index in last_index.slice_mut().iter_mut() {
-            *index -= 1;
-        }
-        if let Some(offset) = stride_offset_checked_arithmetic(dim,
-                                                               strides,
-                                                               &last_index)
-        {
-            // offset is guaranteed to be positive so no issue converting
-            // to usize here
-            if (offset as usize) >= data.len() {
-                return Err(from_kind(ErrorKind::OutOfBounds));
-            }
-            if dim_stride_overlap(dim, strides) {
-                return Err(from_kind(ErrorKind::Unsupported));
-            }
-        } else {
+    }
+    if len == 0 {
+        return Ok(());
+    }
+    // check that the maximum index is in bounds
+    let mut last_index = dim.clone();
+    for mut index in last_index.slice_mut().iter_mut() {
+        *index -= 1;
+    }
+    if let Some(offset) = stride_offset_checked_arithmetic(dim,
+                                                           strides,
+                                                           &last_index)
+    {
+        // offset is guaranteed to be positive so no issue converting
+        // to usize here
+        if (offset as usize) >= data.len() {
             return Err(from_kind(ErrorKind::OutOfBounds));
         }
-        Ok(())
+        if dim_stride_overlap(dim, strides) {
+            return Err(from_kind(ErrorKind::Unsupported));
+        }
     } else {
-        return Err(from_kind(ErrorKind::Unsupported));
+        return Err(from_kind(ErrorKind::OutOfBounds));
     }
+    Ok(())
 }
 
 /// Return stride offset for this dimension and index.
@@ -335,6 +325,21 @@ pub unsafe trait Dimension : Clone + Eq + Debug + Send + Sync {
         offset
     }
 
+    /// 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.
+    #[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
+    }
 }
 
 /// Implementation-specific extensions to `Dimension`
@@ -484,6 +489,11 @@ unsafe impl Dimension for (Ix, Ix) {
         (self.1, 1)
     }
 
+    #[inline]
+    fn _fastest_varying_stride_order(&self) -> Self {
+        if self.0 as Ixs <= self.1 as Ixs { (0, 1) } else { (1, 0) }
+    }
+
     #[inline]
     fn first_index(&self) -> Option<(Ix, Ix)> {
         let (m, n) = *self;
@@ -563,6 +573,29 @@ unsafe impl Dimension for (Ix, Ix, Ix) {
         let (s, t, u) = *strides;
         stride_offset(i, s) + stride_offset(j, t) + stride_offset(k, u)
     }
+
+    #[inline]
+    fn _fastest_varying_stride_order(&self) -> Self {
+        let mut stride = *self;
+        let mut order = (0, 1, 2);
+        macro_rules! swap {
+            ($stride:expr, $order:expr, $x:expr, $y:expr) => {
+                if $stride[$x] > $stride[$y] {
+                    $stride.swap($x, $y);
+                    $order.swap($x, $y);
+                }
+            }
+        }
+        {
+            // stable sorting network for 3 elements
+            let order = order.slice_mut();
+            let strides = stride.slice_mut();
+            swap![strides, order, 1, 2];
+            swap![strides, order, 0, 1];
+            swap![strides, order, 1, 2];
+        }
+        order
+    }
 }
 
 macro_rules! large_dim {
@@ -742,13 +775,6 @@ mod test {
     use super::Dimension;
     use error::StrideError;
 
-    #[test]
-    fn fastest_varying_order() {
-        let strides = (2, 8, 4, 1);
-        let order = super::fastest_varying_order(&strides);
-        assert_eq!(order.slice(), &[3, 0, 2, 1]);
-    }
-
     #[test]
     fn slice_indexing_uncommon_strides() {
         let v: Vec<_> = (0..12).collect();

src/error.rs

@@ -87,13 +87,13 @@ impl Error for ShapeError {
 
 impl fmt::Display for ShapeError {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        self.description().fmt(f)
+        write!(f, "ShapeError/{:?}: {}", self.kind(), self.description())
     }
 }
 
 impl fmt::Debug for ShapeError {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
-        write!(f, "ShapeError {:?}: {}", self.kind(), self.description())
+        write!(f, "ShapeError/{:?}: {}", self.kind(), self.description())
     }
 }