Merge pull request #60 from vbarrielle/stride_constructor

Implement owning constructors taking dimensions and strides

Author: bluss

src/dimension.rs

@@ -2,6 +2,7 @@ use std::slice;
 
 use super::{Si, Ix, Ixs};
 use super::zipsl;
+use stride_error::StrideError;
 
 /// Calculate offset from `Ix` stride converting sign properly
 #[inline]
@@ -10,6 +11,88 @@ pub fn stride_offset(n: Ix, stride: Ix) -> isize
     (n as isize) * ((stride as Ixs) as isize)
 }
 
+/// Check whether `stride` is strictly positive
+#[inline]
+pub 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 (ind, &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] = ind;
+    }
+    res
+}
+
+/// Check whether the given `dim` and `stride` lead to overlapping indices
+///
+/// There is overlap if, when iterating through the dimensions in the order
+/// of maximum variation, the current stride is inferior to the sum of all
+/// preceding strides multiplied by their corresponding dimensions.
+///
+/// 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 mut prev_offset = 1;
+    for &ind in order.slice().iter() {
+        let s = strides.slice()[ind];
+        if (s as isize) < prev_offset {
+            return true;
+        }
+        prev_offset = stride_offset(dim.slice()[ind], s);
+    }
+    false
+}
+
+/// Check whether the given dimension and strides are memory safe
+/// to index the provided slice.
+///
+/// To be safe, no stride may be negative, and the offset corresponding
+/// to the last element of each dimension should be smaller than the length
+/// of the slice. Also, the strides should not allow a same element to be
+/// referenced by two different index.
+pub fn can_index_slice<A, D: Dimension>(data: &[A],
+                                        dim: &D,
+                                        strides: &D
+                                       ) -> Result<(), StrideError>
+{
+    if strides.slice().iter().cloned().all(stride_is_positive) {
+        let mut last_index = dim.clone();
+        for mut index in last_index.slice_mut().iter_mut() {
+            *index -= 1;
+        }
+        if let Some(offset) = dim.stride_offset_checked(strides, &last_index) {
+            // offset is guaranteed to be positive so no issue converting
+            // to usize here
+            if (offset as usize) >= data.len() {
+                return Err(StrideError::OutOfBounds);
+            }
+            if dim_stride_overlap(dim, strides) {
+                return Err(StrideError::Aliasing);
+            }
+        }
+        Ok(())
+    }
+    else {
+        Err(StrideError::Aliasing)
+    }
+}
+
 /// Trait for the shape and index types of arrays.
 ///
 /// `unsafe` because of the assumptions in the default methods.
@@ -69,6 +152,26 @@ pub unsafe trait Dimension : Clone + Eq {
         strides
     }
 
+    fn fortran_strides(&self) -> Self {
+        // Compute fortran array strides
+        // Shape (a, b, c) => Give strides (1, a, a * b)
+        let mut strides = self.clone();
+        {
+            let mut it = strides.slice_mut().iter_mut();
+            // Set first element to 1
+            for rs in it.by_ref() {
+                *rs = 1;
+                break;
+            }
+            let mut cum_prod = 1;
+            for (rs, dim) in it.zip(self.slice().iter()) {
+                cum_prod *= *dim;
+                *rs = cum_prod;
+            }
+        }
+        strides
+    }
+
     #[inline]
     fn first_index(&self) -> Option<Self>
     {
@@ -529,3 +632,39 @@ unsafe impl<'a> NdIndex for &'a [Ix] {
         Some(offset)
     }
 }
+
+#[cfg(test)]
+mod test {
+    use super::{Dimension};
+    use stride_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();
+        let dim = (2, 3, 2);
+        let strides = (1, 2, 6);
+        assert!(super::can_index_slice(&v, &dim, &strides).is_ok());
+
+        let strides = (2, 4, 12);
+        assert_eq!(super::can_index_slice(&v, &dim, &strides),
+                   Err(StrideError::OutOfBounds));
+    }
+
+    #[test]
+    fn overlapping_strides_dim()
+    {
+        let dim = (2, 3, 2);
+        let strides = (5, 2, 1);
+        assert!(super::dim_stride_overlap(&dim, &strides));
+        let strides = (6, 2, 1);
+        assert!(!super::dim_stride_overlap(&dim, &strides));
+    }
+}

src/lib.rs

@@ -75,10 +75,15 @@ use std::marker::PhantomData;
 
 use itertools::ZipSlices;
 
-pub use dimension::{Dimension, RemoveAxis};
+pub use dimension::{
+    Dimension,
+    RemoveAxis,
+};
+
 pub use dimension::NdIndex;
 pub use indexes::Indexes;
 pub use shape_error::ShapeError;
+pub use stride_error::StrideError;
 pub use si::{Si, S};
 
 use dimension::stride_offset;
@@ -107,6 +112,7 @@ mod linspace;
 mod numeric_util;
 mod si;
 mod shape_error;
+mod stride_error;
 
 // NOTE: In theory, the whole library should compile
 // and pass tests even if you change Ix and Ixs.
@@ -605,12 +611,43 @@ impl<S, A, D> ArrayBase<S, D>
         }
     }
 
+    /// Create an array with copies of `elem`, dimension `dim` and fortran
+    /// ordering.
+    ///
+    /// ```
+    /// use ndarray::Array;
+    /// use ndarray::arr3;
+    ///
+    /// let a = Array::from_elem_f((2, 2, 2), 1.);
+    ///
+    /// assert!(
+    ///     a == arr3(&[[[1., 1.],
+    ///                  [1., 1.]],
+    ///                 [[1., 1.],
+    ///                  [1., 1.]]])
+    /// );
+    /// assert!(a.strides() == &[1, 2, 4]);
+    /// ```
+    pub fn from_elem_f(dim: D, elem: A) -> ArrayBase<S, D> where A: Clone
+    {
+        let v = vec![elem; dim.size()];
+        unsafe {
+            Self::from_vec_dim_f(dim, v)
+        }
+    }
+
     /// Create an array with zeros, dimension `dim`.
     pub fn zeros(dim: D) -> ArrayBase<S, D> where A: Clone + libnum::Zero
     {
         Self::from_elem(dim, libnum::zero())
     }
 
+    /// Create an array with zeros, dimension `dim` and fortran ordering.
+    pub fn zeros_f(dim: D) -> ArrayBase<S, D> where A: Clone + libnum::Zero
+    {
+        Self::from_elem_f(dim, libnum::zero())
+    }
+
     /// Create an array with default values, dimension `dim`.
     pub fn default(dim: D) -> ArrayBase<S, D>
         where A: Default
@@ -634,6 +671,58 @@ impl<S, A, D> ArrayBase<S, D>
             dim: dim
         }
     }
+
+    /// Create an array from a vector (with no allocation needed),
+    /// using fortran ordering to interpret the data.
+    ///
+    /// Unsafe because dimension is unchecked, and must be correct.
+    pub unsafe fn from_vec_dim_f(dim: D, mut v: Vec<A>) -> ArrayBase<S, D>
+    {
+        debug_assert!(dim.size() == v.len());
+        ArrayBase {
+            ptr: v.as_mut_ptr(),
+            data: DataOwned::new(v),
+            strides: dim.fortran_strides(),
+            dim: dim
+        }
+    }
+
+
+    /// Create an array from a vector and interpret it according to the
+    /// provided dimensions and strides. No allocation needed.
+    ///
+    /// Unsafe because dimension and strides are unchecked.
+    pub unsafe fn from_vec_dim_stride_uchk(dim: D,
+                                           strides: D,
+                                           mut v: Vec<A>
+                                          ) -> ArrayBase<S, D>
+    {
+        ArrayBase {
+            ptr: v.as_mut_ptr(),
+            data: DataOwned::new(v),
+            strides: strides,
+            dim: dim
+        }
+    }
+
+    /// Create an array from a vector and interpret it according to the
+    /// provided dimensions and strides. No allocation needed.
+    ///
+    /// Checks whether `dim` and `strides` are compatible with the vector's
+    /// length, returning an `Err` if not compatible.
+    pub fn from_vec_dim_stride(dim: D,
+                               strides: D,
+                               v: Vec<A>
+                              ) -> Result<ArrayBase<S, D>, StrideError>
+    {
+        dimension::can_index_slice(&v,
+                                   &dim,
+                                   &strides).map(|_| {
+            unsafe {
+                Self::from_vec_dim_stride_uchk(dim, strides, v)
+            }
+        })
+    }
 }
 
 
@@ -666,6 +755,44 @@ impl<'a, A, D> ArrayView<'a, A, D>
         }
     }
 
+    /// Create an `ArrayView` borrowing its data from a slice.
+    ///
+    /// Checks whether `dim` and `strides` are compatible with the slice's
+    /// length, returning an `Err` if not compatible.
+    ///
+    /// ```
+    /// use ndarray::ArrayView;
+    /// use ndarray::arr3;
+    ///
+    /// let s = &[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12];
+    /// let a = ArrayView::from_slice_dim_stride((2, 3, 2),
+    ///                                          (1, 4, 2),
+    ///                                          s).unwrap();
+    ///
+    /// assert!(
+    ///     a == arr3(&[[[0, 2],
+    ///                  [4, 6],
+    ///                  [8, 10]],
+    ///                 [[1, 3],
+    ///                  [5, 7],
+    ///                  [9, 11]]])
+    /// );
+    /// assert!(a.strides() == &[1, 4, 2]);
+    /// ```
+    pub fn from_slice_dim_stride(dim: D,
+                                 strides: D,
+                                 s: &'a [A]
+                                 ) -> Result<Self, StrideError>
+    {
+        dimension::can_index_slice(s,
+                                   &dim,
+                                   &strides).map(|_| {
+            unsafe {
+                Self::new_(s.as_ptr(), dim, strides)
+            }
+        })
+    }
+
     #[inline]
     fn into_base_iter(self) -> Baseiter<'a, A, D> {
         unsafe {
@@ -723,6 +850,45 @@ impl<'a, A, D> ArrayViewMut<'a, A, D>
         }
     }
 
+    /// Create an `ArrayView` borrowing its data from a slice.
+    ///
+    /// Checks whether `dim` and `strides` are compatible with the slice's
+    /// length, returning an `Err` if not compatible.
+    ///
+    /// ```
+    /// use ndarray::ArrayViewMut;
+    /// use ndarray::arr3;
+    ///
+    /// let s = &mut [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12];
+    /// let mut a = ArrayViewMut::from_slice_dim_stride((2, 3, 2),
+    ///                                                 (1, 4, 2),
+    ///                                                 s).unwrap();
+    ///
+    /// a[[0, 0, 0]] = 1;
+    /// assert!(
+    ///     a == arr3(&[[[1, 2],
+    ///                  [4, 6],
+    ///                  [8, 10]],
+    ///                 [[1, 3],
+    ///                  [5, 7],
+    ///                  [9, 11]]])
+    /// );
+    /// assert!(a.strides() == &[1, 4, 2]);
+    /// ```
+    pub fn from_slice_dim_stride(dim: D,
+                                 strides: D,
+                                 s: &'a mut [A]
+                                 ) -> Result<Self, StrideError>
+    {
+        dimension::can_index_slice(s,
+                                   &dim,
+                                   &strides).map(|_| {
+            unsafe {
+                Self::new_(s.as_mut_ptr(), dim, strides)
+            }
+        })
+    }
+
     #[inline]
     fn into_base_iter(self) -> Baseiter<'a, A, D> {
         unsafe {