Variance along an axis

src/impl_methods.rs

@@ -124,7 +124,36 @@ impl<A, S, D> ArrayBase<S, D> where S: Data<Elem=A>, D: Dimension
         }
     }
 
-    /// Return an uniquely owned copy of the array
+    /// Return an uniquely owned copy of the array.
+    ///
+    /// If the input array is contiguous and its strides are positive, then the
+    /// output array will have the same memory layout. Otherwise, the layout of
+    /// the output array is unspecified. If you need a particular layout, you
+    /// can allocate a new array with the desired memory layout and
+    /// [`.assign()`](#method.assign) the data. Alternatively, you can collect
+    /// an iterator, like this for a result in standard layout:
+    ///
+    /// ```
+    /// # use ndarray::prelude::*;
+    /// # let arr = Array::from_shape_vec((2, 2).f(), vec![1, 2, 3, 4]).unwrap();
+    /// # let owned = {
+    /// Array::from_shape_vec(arr.raw_dim(), arr.iter().cloned().collect()).unwrap()
+    /// # };
+    /// # assert!(owned.is_standard_layout());
+    /// # assert_eq!(arr, owned);
+    /// ```
+    ///
+    /// or this for a result in column-major (Fortran) layout:
+    ///
+    /// ```
+    /// # use ndarray::prelude::*;
+    /// # let arr = Array::from_shape_vec((2, 2), vec![1, 2, 3, 4]).unwrap();
+    /// # let owned = {
+    /// Array::from_shape_vec(arr.raw_dim().f(), arr.t().iter().cloned().collect()).unwrap()
+    /// # };
+    /// # assert!(owned.t().is_standard_layout());
+    /// # assert_eq!(arr, owned);
+    /// ```
     pub fn to_owned(&self) -> Array<A, D>
         where A: Clone
     {

src/lib.rs

@@ -168,6 +168,7 @@ mod numeric_util;
 mod error;
 mod shape_builder;
 mod stacking;
+#[macro_use]
 mod zip;
 
 pub use zip::{
@@ -239,7 +240,7 @@ pub type Ixs = isize;
 ///
 /// ## `Array`
 ///
-/// [`Array`](type.Array.html) is an owned array that ows the underlying array
+/// [`Array`](type.Array.html) is an owned array that owns the underlying array
 /// elements directly (just like a `Vec`) and it is the default way to create and
 /// store n-dimensional data. `Array<A, D>` has two type parameters: `A` for
 /// the element type, and `D` for the dimensionality. A particular

src/numeric/impl_numeric.rs

@@ -115,6 +115,69 @@ impl<A, S, D> ArrayBase<S, D>
         sum / &aview0(&cnt)
     }
 
+    /// Return variance along `axis`.
+    ///
+    /// The variance is computed using the [Welford one-pass
+    /// algorithm](https://www.jstor.org/stable/1266577).
+    ///
+    /// The parameter `ddof` specifies the "delta degrees of freedom". For
+    /// example, to calculate the population variance, use `ddof = 0`, or to
+    /// calculate the sample variance, use `ddof = 1`.
+    ///
+    /// The variance is defined as:
+    ///
+    /// ```text
+    ///               1       n
+    /// variance = ――――――――   ∑ (xᵢ - x̅)²
+    ///            n - ddof  i=1
+    /// ```
+    ///
+    /// where
+    ///
+    /// ```text
+    ///     1   n
+    /// x̅ = ―   ∑ xᵢ
+    ///     n  i=1
+    /// ```
+    ///
+    /// **Panics** if `ddof` is greater equal than the length of `axis`.
+    /// **Panics** if `axis` is out of bounds or if length of `axis` is zero.
+    ///
+    /// # Example
+    ///
+    /// ```
+    /// use ndarray::{aview1, arr2, Axis};
+    ///
+    /// let a = arr2(&[[1., 2.],
+    ///                [3., 4.]]);
+    /// let var = a.var_axis(Axis(0), 0.);
+    /// assert_eq!(var, aview1(&[1., 1.]));
+    /// ```
+    pub fn var_axis(&self, axis: Axis, ddof: A) -> Array<A, D::Smaller>
+    where
+        A: Float,
+        D: RemoveAxis,
+    {
+        let mut count = A::zero();
+        let mut mean = Array::<A, _>::zeros(self.dim.remove_axis(axis));
+        let mut sum_sq = Array::<A, _>::zeros(self.dim.remove_axis(axis));
+        for subview in self.axis_iter(axis) {
+            count = count + A::one();
+            azip!(mut mean, mut sum_sq, x (subview) in {
+                let delta = x - *mean;
+                *mean = *mean + delta / count;
+                *sum_sq = *sum_sq + delta * (x - *mean);
+            });
+        }
+        if ddof >= count {
+            panic!("Ddof needs to be strictly smaller than the length \
+                    of the axis you are computing the variance for!")
+        } else {
+            let dof = count - ddof;
+            sum_sq.mapv(|s| s / dof)
+        }
+    }
+
     /// Return `true` if the arrays' elementwise differences are all within
     /// the given absolute tolerance, `false` otherwise.
     ///
@@ -137,4 +200,3 @@ impl<A, S, D> ArrayBase<S, D>
             }).is_done()
     }
 }
-

src/zip/mod.rs

@@ -6,6 +6,7 @@
 // option. This file may not be copied, modified, or distributed
 // except according to those terms.
 
+#[macro_use]
 mod zipmacro;
 
 use imp_prelude::*;

tests/array.rs

@@ -274,28 +274,37 @@ fn test_slice_inplace_with_subview_inplace() {
         *elt = i;
     }
 
-    let mut vi = arr.view();
-    vi.slice_inplace(s![1.., 2, ..;2]);
-    assert_eq!(vi.shape(), &[2, 1, 2]);
-    assert!(
-        vi.iter()
-            .zip(arr.slice(s![1.., 2..3, ..;2]).iter())
-            .all(|(a, b)| a == b)
-    );
-
-    let mut vi = arr.view();
-    vi.slice_inplace(s![1, 2, ..;2]);
-    assert_eq!(vi.shape(), &[1, 1, 2]);
-    assert!(
-        vi.iter()
-            .zip(arr.slice(s![1..2, 2..3, ..;2]).iter())
-            .all(|(a, b)| a == b)
-    );
+    {
+        let mut vi = arr.view();
+        vi.slice_inplace(s![1.., 2, ..;2]);
+        assert_eq!(vi.shape(), &[2, 1, 2]);
+        assert!(
+            vi.iter()
+                .zip(arr.slice(s![1.., 2..3, ..;2]).iter())
+                .all(|(a, b)| a == b)
+        );
+
+        let mut vi = arr.view();
+        vi.slice_inplace(s![1, 2, ..;2]);
+        assert_eq!(vi.shape(), &[1, 1, 2]);
+        assert!(
+            vi.iter()
+                .zip(arr.slice(s![1..2, 2..3, ..;2]).iter())
+                .all(|(a, b)| a == b)
+        );
+
+        let mut vi = arr.view();
+        vi.slice_inplace(s![1, 2, 3]);
+        assert_eq!(vi.shape(), &[1, 1, 1]);
+        assert_eq!(vi, Array3::from_elem((1, 1, 1), arr[(1, 2, 3)]));
+    }
 
-    let mut vi = arr.view();
+    // Do it to the RcArray itself
+    let elem = arr[(1, 2, 3)];
+    let mut vi = arr;
     vi.slice_inplace(s![1, 2, 3]);
     assert_eq!(vi.shape(), &[1, 1, 1]);
-    assert_eq!(vi, Array3::from_elem((1, 1, 1), arr[(1, 2, 3)]));
+    assert_eq!(vi, Array3::from_elem((1, 1, 1), elem));
 }
 
 #[should_panic]
@@ -807,6 +816,44 @@ fn owned_array_with_stride() {
     assert_eq!(a.strides(), &[1, 4, 2]);
 }
 
+#[test]
+fn owned_array_discontiguous() {
+    use ::std::iter::repeat;
+    let v: Vec<_> = (0..12).flat_map(|x| repeat(x).take(2)).collect();
+    let dim = (3, 2, 2);
+    let strides = (8, 4, 2);
+
+    let a = Array::from_shape_vec(dim.strides(strides), v).unwrap();
+    assert_eq!(a.strides(), &[8, 4, 2]);
+    println!("{:?}", a.iter().cloned().collect::<Vec<_>>());
+    itertools::assert_equal(a.iter().cloned(), 0..12);
+}
+
+#[test]
+fn owned_array_discontiguous_drop() {
+    use ::std::rc::Rc;
+    use ::std::cell::RefCell;
+    use ::std::collections::BTreeSet;
+
+    struct InsertOnDrop<T: Ord>(Rc<RefCell<BTreeSet<T>>>, Option<T>);
+    impl<T: Ord> Drop for InsertOnDrop<T> {
+        fn drop(&mut self) {
+            let InsertOnDrop(ref set, ref mut value) = *self;
+            set.borrow_mut().insert(value.take().expect("double drop!"));
+        }
+    }
+
+    let set = Rc::new(RefCell::new(BTreeSet::new()));
+    {
+        let v: Vec<_> = (0..12).map(|x| InsertOnDrop(set.clone(), Some(x))).collect();
+        let mut a = Array::from_shape_vec((2, 6), v).unwrap();
+        // discontiguous and non-zero offset
+        a.slice_inplace(s![.., 1..]);
+    }
+    // each item was dropped exactly once
+    itertools::assert_equal(set.borrow().iter().cloned(), 0..12);
+}
+
 macro_rules! assert_matches {
     ($value:expr, $pat:pat) => {
         match $value {
@@ -1372,6 +1419,18 @@ fn to_owned_neg_stride() {
     assert_eq!(c, co);
 }
 
+#[test]
+fn discontiguous_owned_to_owned() {
+    let mut c = arr2(&[[1, 2, 3],
+                       [4, 5, 6]]);
+    c.slice_inplace(s![.., ..;2]);
+
+    let co = c.to_owned();
+    assert_eq!(c.strides(), &[3, 2]);
+    assert_eq!(co.strides(), &[2, 1]);
+    assert_eq!(c, co);
+}
+
 #[test]
 fn map_memory_order() {
     let a = arr3(&[[[1, 2, 3],