append: Fix situations where we need to recompute stride

When the axis has length 0, or 1, we need to carefully compute new
strides.

Author: bluss

src/impl_owned_array.rs

@@ -261,6 +261,7 @@ impl<A, D> Array<A, D>
             return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape));
         }
 
+        let current_axis_len = self.len_of(axis);
         let remaining_shape = self.raw_dim().remove_axis(axis);
         let array_rem_shape = array.raw_dim().remove_axis(axis);
 
@@ -280,22 +281,46 @@ impl<A, D> Array<A, D>
 
         let self_is_empty = self.is_empty();
 
-        // array must be empty or have `axis` as the outermost (longest stride)
-        // axis
-        if !(self_is_empty ||
-             self.axes().max_by_key(|ax| ax.stride).map(|ax| ax.axis) == Some(axis))
-        {
-            return Err(ShapeError::from_kind(ErrorKind::IncompatibleLayout));
+        // array must be empty or have `axis` as the outermost (longest stride) axis
+        if !self_is_empty && current_axis_len > 1 {
+            // `axis` must be max stride axis or equal to its stride
+            let max_stride_axis = self.axes().max_by_key(|ax| ax.stride).unwrap();
+            if max_stride_axis.axis != axis && max_stride_axis.stride > self.stride_of(axis) {
+                return Err(ShapeError::from_kind(ErrorKind::IncompatibleLayout));
+            }
         }
 
         // array must be be "full" (have no exterior holes)
         if self.len() != self.data.len() {
             return Err(ShapeError::from_kind(ErrorKind::IncompatibleLayout));
         }
+
         let strides = if self_is_empty {
-            // recompute strides - if the array was previously empty, it could have
-            // zeros in strides.
-            res_dim.default_strides()
+            // recompute strides - if the array was previously empty, it could have zeros in
+            // strides.
+            // The new order is based on c/f-contig but must have `axis` as outermost axis.
+            if axis == Axis(self.ndim() - 1) {
+                // prefer f-contig when appending to the last axis
+                // Axis n - 1 is outermost axis
+                res_dim.fortran_strides()
+            } else {
+                // Default with modification
+                res_dim.slice_mut().swap(0, axis.index());
+                let mut strides = res_dim.default_strides();
+                res_dim.slice_mut().swap(0, axis.index());
+                strides.slice_mut().swap(0, axis.index());
+                strides
+            }
+        } else if current_axis_len == 1 {
+            // This is the outermost/longest stride axis; so we find the max across the other axes
+            let new_stride = self.axes().fold(1, |acc, ax| {
+                if ax.axis == axis { acc } else {
+                    Ord::max(acc, ax.len as isize * ax.stride)
+                }
+            });
+            let mut strides = self.strides.clone();
+            strides[axis.index()] = new_stride as usize;
+            strides
         } else {
             self.strides.clone()
         };
@@ -383,7 +408,8 @@ where
         return;
     }
     sort_axes_impl(&mut a.dim, &mut a.strides, &mut b.dim, &mut b.strides);
-    debug_assert!(a.is_standard_layout());
+    debug_assert!(a.is_standard_layout(), "not std layout dim: {:?}, strides: {:?}",
+                  a.shape(), a.strides());
 }
 
 fn sort_axes_impl<D>(adim: &mut D, astrides: &mut D, bdim: &mut D, bstrides: &mut D)

tests/append.rs

@@ -146,3 +146,52 @@ fn append_array_3d() {
                 [83, 84],
                 [87, 88]]]);
 }
+
+#[test]
+fn test_append_2d() {
+    // create an empty array and append
+    let mut a = Array::zeros((0, 4));
+    let ones = ArrayView::from(&[1.; 12]).into_shape((3, 4)).unwrap();
+    let zeros = ArrayView::from(&[0.; 8]).into_shape((2, 4)).unwrap();
+    a.try_append_array(Axis(0), ones).unwrap();
+    a.try_append_array(Axis(0), zeros).unwrap();
+    a.try_append_array(Axis(0), ones).unwrap();
+    println!("{:?}", a);
+    assert_eq!(a.shape(), &[8, 4]);
+    for (i, row) in a.rows().into_iter().enumerate() {
+        let ones = i < 3 || i >= 5;
+        assert!(row.iter().all(|&x| x == ones as i32 as f64), "failed on lane {}", i);
+    }
+
+    let mut a = Array::zeros((0, 4));
+    a = a.reversed_axes();
+    let ones = ones.reversed_axes();
+    let zeros = zeros.reversed_axes();
+    a.try_append_array(Axis(1), ones).unwrap();
+    a.try_append_array(Axis(1), zeros).unwrap();
+    a.try_append_array(Axis(1), ones).unwrap();
+    println!("{:?}", a);
+    assert_eq!(a.shape(), &[4, 8]);
+
+    for (i, row) in a.columns().into_iter().enumerate() {
+        let ones = i < 3 || i >= 5;
+        assert!(row.iter().all(|&x| x == ones as i32 as f64), "failed on lane {}", i);
+    }
+}
+
+#[test]
+fn test_append_middle_axis() {
+    // ensure we can append to Axis(1) by letting it become outermost
+    let mut a = Array::<i32, _>::zeros((3, 0, 2));
+    a.try_append_array(Axis(1), Array::from_iter(0..12).into_shape((3, 2, 2)).unwrap().view()).unwrap();
+    println!("{:?}", a);
+    a.try_append_array(Axis(1), Array::from_iter(12..24).into_shape((3, 2, 2)).unwrap().view()).unwrap();
+    println!("{:?}", a);
+
+    // ensure we can append to Axis(1) by letting it become outermost
+    let mut a = Array::<i32, _>::zeros((3, 1, 2));
+    a.try_append_array(Axis(1), Array::from_iter(0..12).into_shape((3, 2, 2)).unwrap().view()).unwrap();
+    println!("{:?}", a);
+    a.try_append_array(Axis(1), Array::from_iter(12..24).into_shape((3, 2, 2)).unwrap().view()).unwrap();
+    println!("{:?}", a);
+}