From ee0d54d3504158a3e8d29e9cf7745f43df330916 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 07:17:51 +0100
Subject: [PATCH 01/93] Added function signature.

---
 src/numeric/impl_numeric.rs | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index b1ef676cc..530c64836 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -115,6 +115,15 @@ impl<A, S, D> ArrayBase<S, D>
         sum / &aview0(&cnt)
     }
 
+
+    /// Return the qth percentile of the data along the specified axis.
+    pub fn percentile_axis<E>(&self, axis: Axis, q: E) -> Array<A, D::Smaller>
+        where E: Float,
+              D: RemoveAxis,
+    {
+        unimplemented!()
+    }
+
     /// Return `true` if the arrays' elementwise differences are all within
     /// the given absolute tolerance, `false` otherwise.
     ///

From 9b775d013e69b307e8fd6070b6aa56573b6b3b9e Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 07:36:41 +0100
Subject: [PATCH 02/93] Determine the index of the desired element if the array
 were to be sorted.

---
 src/numeric/impl_numeric.rs | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 530c64836..1a5fc36fc 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -117,10 +117,11 @@ impl<A, S, D> ArrayBase<S, D>
 
 
     /// Return the qth percentile of the data along the specified axis.
-    pub fn percentile_axis<E>(&self, axis: Axis, q: E) -> Array<A, D::Smaller>
-        where E: Float,
-              D: RemoveAxis,
+    pub fn percentile_axis(&self, axis: Axis, q: f32) -> Array<A, D::Smaller>
+        where D: RemoveAxis,
     {
+        let n = self.len_of(axis);
+        let i = ((n as f32) * q).floor() as usize;
         unimplemented!()
     }
 

From a4bec9d9058fc87fc9e04131f07cadf4242a659d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 08:08:22 +0100
Subject: [PATCH 03/93] Added macro use for slice

---
 src/lib.rs | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 9f05581b1..375d7912d 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -158,6 +158,7 @@ mod free_functions;
 pub use free_functions::*;
 pub use iterators::iter;
 
+#[macro_use]
 mod slice;
 mod layout;
 mod indexes;
@@ -339,7 +340,7 @@ pub type Ixs = isize;
 /// A *column major* (a.k.a. “f” or fortran) memory order array has
 /// columns (or, in general, the outermost axis) with contiguous elements.
 ///
-/// The logical order of any array’s elements is the row major order 
+/// The logical order of any array’s elements is the row major order
 /// (the rightmost index is varying the fastest).
 /// The iterators `.iter(), .iter_mut()` always adhere to this order, for example.
 ///
@@ -366,7 +367,7 @@ pub type Ixs = isize;
 ///
 /// These are the element iterators of arrays and they produce an element
 /// sequence in the logical order of the array, that means that the elements
-/// will be visited in the sequence that corresponds to increasing the 
+/// will be visited in the sequence that corresponds to increasing the
 /// last index first: *0, ..., 0,  0*; *0, ..., 0, 1*; *0, ...0, 2* and so on.
 ///
 /// ### `.outer_iter()` and `.axis_iter()`
@@ -683,7 +684,7 @@ pub type RcArray<A, D> = ArrayBase<OwnedRcRepr<A>, D>;
 /// `Sync` (when allowed by the element type of the array too).
 ///
 /// [**`ArrayBase`**](struct.ArrayBase.html) is used to implement both the owned
-/// arrays and the views; see its docs for an overview of all array features.  
+/// arrays and the views; see its docs for an overview of all array features.
 ///
 /// See also:
 ///
@@ -700,7 +701,7 @@ pub type ArcArray<A, D> = ArrayBase<OwnedArcRepr<A>, D>;
 /// the dimensionality.
 ///
 /// [**`ArrayBase`**](struct.ArrayBase.html) is used to implement both the owned
-/// arrays and the views; see its docs for an overview of all array features.  
+/// arrays and the views; see its docs for an overview of all array features.
 ///
 /// See also:
 ///

From 470af6cd9a65c6ee3ff5e4192a369b179b6900f0 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 08:33:37 +0100
Subject: [PATCH 04/93] Barebone implementation of randomized_select with the
 signature of randomized_partition.

---
 src/numeric/impl_numeric.rs | 23 +++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 1a5fc36fc..197a5f709 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -148,3 +148,26 @@ impl<A, S, D> ArrayBase<S, D>
     }
 }
 
+fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
+    where A: Ord + Clone
+{
+    let n = a.len(); if n == 0 {
+        (&a[0]).clone()
+    } else {
+        let (q, mut a) = randomized_partition(&a);
+        let k = q + 1;
+        if i == k {
+            (&a[q]).clone()
+        } else if i < k {
+            randomized_select(a.slice_mut(s![0..q]).to_owned(), i)
+        } else {
+            randomized_select(a.slice_mut(s![(q+1)..n]).to_owned(), i - k)
+        }
+    }
+}
+
+fn randomized_partition<A>(mut a: &Array1<A>) -> (usize, Array1<A>)
+    where A: Ord + Clone
+{
+    unimplemented!()
+}

From f2eec962c872dcc8d93c7eb9e8e28edc511cf2b1 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 08:35:53 +0100
Subject: [PATCH 05/93] Added rand as a dependency.

---
 Cargo.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Cargo.toml b/Cargo.toml
index 0557fd63a..607d4fb71 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,6 +28,7 @@ num-traits = { version = "0.1.32", default-features = false }
 num-complex = { version = "0.1.32", default-features = false }
 rustc-serialize = { version = "0.3.20", optional = true }
 itertools = { version = "0.7.0", default-features = false }
+rand = { version = "0.5.0-pre.1" }
 
 # Use via the `blas` crate feature!
 cblas-sys = { version = "0.1.4", optional = true, default-features = false }

From 92d343388cc44ddeb1224cba4ef00f4ea7dd4618 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 08:42:50 +0100
Subject: [PATCH 06/93] Added extern crate rand to lib.rs.

---
 src/lib.rs | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/lib.rs b/src/lib.rs
index 375d7912d..3e918d09b 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -100,6 +100,8 @@ extern crate matrixmultiply;
 extern crate num_traits as libnum;
 extern crate num_complex;
 
+extern crate rand;
+
 #[cfg(feature = "docs")]
 pub mod doc;
 

From 8e25ae0f4e0dd4ff7f91de77457a602ba962762b Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 14 May 2018 08:55:56 +0100
Subject: [PATCH 07/93] Basic implementation of randomized_partition first
 section. Missing partition.

---
 src/numeric/impl_numeric.rs | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 197a5f709..5defd24ce 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -13,6 +13,9 @@ use itertools::free::enumerate;
 use imp_prelude::*;
 use numeric_util;
 
+use rand::distributions::Uniform;
+use rand::thread_rng;
+
 use {
     LinalgScalar,
     FoldWhile,
@@ -154,7 +157,7 @@ fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
     let n = a.len(); if n == 0 {
         (&a[0]).clone()
     } else {
-        let (q, mut a) = randomized_partition(&a);
+        let (q, mut a) = randomized_partition(&mut a);
         let k = q + 1;
         if i == k {
             (&a[q]).clone()
@@ -166,8 +169,14 @@ fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
     }
 }
 
-fn randomized_partition<A>(mut a: &Array1<A>) -> (usize, Array1<A>)
+fn randomized_partition<A>(a: &mut Array1<A>) -> (usize, Array1<A>)
     where A: Ord + Clone
 {
+    let n = a.len();
+    let mut rng = thread_rng();
+    let i: usize = Uniform::sample_single(0, n, &mut rng);
+    let z = (&a[i]).clone();
+    a[i] = (&a[n-1]).clone();
+    a[n-1] = z;
     unimplemented!()
 }

From fc3e570884125e9c1cee6ad490424839be5b3842 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 08:47:02 +0100
Subject: [PATCH 08/93] Using swap where needed, instead of manually doing the
 swap using clone. Implemented partition and a small test, which currently
 fails due to index overflow.

---
 src/numeric/impl_numeric.rs | 41 ++++++++++++++++++++++++++++++++-----
 1 file changed, 36 insertions(+), 5 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 5defd24ce..fa89430f2 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -169,14 +169,45 @@ fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
     }
 }
 
-fn randomized_partition<A>(a: &mut Array1<A>) -> (usize, Array1<A>)
+fn randomized_partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
     where A: Ord + Clone
 {
     let n = a.len();
     let mut rng = thread_rng();
     let i: usize = Uniform::sample_single(0, n, &mut rng);
-    let z = (&a[i]).clone();
-    a[i] = (&a[n-1]).clone();
-    a[n-1] = z;
-    unimplemented!()
+    a.swap(i, n-1);
+    // let z = (&a[i]).clone();
+    // a[i] = (&a[n-1]).clone();
+    // a[n-1] = z;
+    partition(a)
+}
+
+fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
+    where A: Ord + Clone
+{
+    let n = a.len();
+    let x = (&a[n-1]).clone();
+    let mut i: usize = 0;
+    for j in 0..n {
+        println!("{:}", j);
+        if a[j] <= x {
+            i += 1;
+            a.swap(i, j)
+            // let z = (&a[i]).clone();
+            // a[i] = (&a[j]).clone();
+            // a[j] = z;
+        }
+    }
+    a.swap(i+1, n-1);
+    // let z = (&a[i+1]).clone();
+    // a[i+1] = (&a[n-1]).clone();
+    // a[n-1] = z;
+    (i+1, a)
+}
+
+#[test]
+fn test_partition() {
+    let mut a = arr1(&[1, 3, 2, 10]);
+    let (j, a) = partition(&mut a);
+    assert_eq!(j, 3);
 }

From 0f5fb47a3538834233f6fca28e31cd11795cee28 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 08:47:29 +0100
Subject: [PATCH 09/93] Cleaned comments and debugging prints.

---
 src/numeric/impl_numeric.rs | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index fa89430f2..59dd70d49 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -176,9 +176,6 @@ fn randomized_partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
     let mut rng = thread_rng();
     let i: usize = Uniform::sample_single(0, n, &mut rng);
     a.swap(i, n-1);
-    // let z = (&a[i]).clone();
-    // a[i] = (&a[n-1]).clone();
-    // a[n-1] = z;
     partition(a)
 }
 
@@ -193,15 +190,9 @@ fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
         if a[j] <= x {
             i += 1;
             a.swap(i, j)
-            // let z = (&a[i]).clone();
-            // a[i] = (&a[j]).clone();
-            // a[j] = z;
         }
     }
     a.swap(i+1, n-1);
-    // let z = (&a[i+1]).clone();
-    // a[i+1] = (&a[n-1]).clone();
-    // a[n-1] = z;
     (i+1, a)
 }
 

From 177e03b85546d8ad2d4c23446a6a81231e56e3d4 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 08:51:42 +0100
Subject: [PATCH 10/93] Fixed index overflow bug.

---
 src/numeric/impl_numeric.rs | 11 +++++------
 1 file changed, 5 insertions(+), 6 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 59dd70d49..af7d70042 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -184,16 +184,15 @@ fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
 {
     let n = a.len();
     let x = (&a[n-1]).clone();
-    let mut i: usize = 0;
-    for j in 0..n {
-        println!("{:}", j);
+    let mut i: isize = -1;
+    for j in 0..n-1 {
         if a[j] <= x {
             i += 1;
-            a.swap(i, j)
+            a.swap(i as usize, j)
         }
     }
-    a.swap(i+1, n-1);
-    (i+1, a)
+    a.swap((i as usize)+1, n-1);
+    ((i as usize)+1, a)
 }
 
 #[test]

From 811e84f01e99080cd4b1a154172b69d447f8d868 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 21:25:28 +0100
Subject: [PATCH 11/93] Fixed issues with index overflow.

---
 src/numeric/impl_numeric.rs | 21 ++++++++++++++++-----
 1 file changed, 16 insertions(+), 5 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index af7d70042..dee38892d 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -157,7 +157,7 @@ fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
     let n = a.len(); if n == 0 {
         (&a[0]).clone()
     } else {
-        let (q, mut a) = randomized_partition(&mut a);
+        let (q, a) = randomized_partition(&mut a);
         let k = q + 1;
         if i == k {
             (&a[q]).clone()
@@ -184,15 +184,15 @@ fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
 {
     let n = a.len();
     let x = (&a[n-1]).clone();
-    let mut i: isize = -1;
+    let mut i = 0;
     for j in 0..n-1 {
         if a[j] <= x {
+            a.swap(i as usize, j);
             i += 1;
-            a.swap(i as usize, j)
         }
     }
-    a.swap((i as usize)+1, n-1);
-    ((i as usize)+1, a)
+    a.swap(i, n-1);
+    (i, a)
 }
 
 #[test]
@@ -201,3 +201,14 @@ fn test_partition() {
     let (j, a) = partition(&mut a);
     assert_eq!(j, 3);
 }
+
+#[test]
+fn test_randomized_select() {
+    let a = arr1(&[1, 3, 2, 10]);
+    let j = randomized_select(a.clone(), 2);
+    assert_eq!(j, 2);
+    let j = randomized_select(a.clone(), 1);
+    assert_eq!(j, 1);
+    let j = randomized_select(a.clone(), 3);
+    assert_eq!(j, 3);
+}

From b6183773b4236659c859aa596ee8b8a9a5e16279 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 21:32:13 +0100
Subject: [PATCH 12/93] Remove useless type annotation.

---
 src/numeric/impl_numeric.rs | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index dee38892d..c3baeabde 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -154,7 +154,8 @@ impl<A, S, D> ArrayBase<S, D>
 fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
     where A: Ord + Clone
 {
-    let n = a.len(); if n == 0 {
+    let n = a.len();
+    if n == 0 {
         (&a[0]).clone()
     } else {
         let (q, a) = randomized_partition(&mut a);
@@ -187,7 +188,7 @@ fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
     let mut i = 0;
     for j in 0..n-1 {
         if a[j] <= x {
-            a.swap(i as usize, j);
+            a.swap(i, j);
             i += 1;
         }
     }
@@ -200,6 +201,9 @@ fn test_partition() {
     let mut a = arr1(&[1, 3, 2, 10]);
     let (j, a) = partition(&mut a);
     assert_eq!(j, 3);
+    let mut a = arr1(&[2, 3, 4, 1]);
+    let (j, a) = partition(&mut a);
+    assert_eq!(j, 0);
 }
 
 #[test]

From 8994eaee92a987a280cb336ca52396aa7993ae11 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 23:26:23 +0100
Subject: [PATCH 13/93] Refactored randomized_partition and partition to
 properly use mutability to modify in place the array passed as argument.
 Added some tests for partition.

---
 src/numeric/impl_numeric.rs | 86 ++++++++++++++++++++++---------------
 1 file changed, 51 insertions(+), 35 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index c3baeabde..977ba8005 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,26 +151,35 @@ impl<A, S, D> ArrayBase<S, D>
     }
 }
 
-fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
-    where A: Ord + Clone
-{
-    let n = a.len();
-    if n == 0 {
-        (&a[0]).clone()
-    } else {
-        let (q, a) = randomized_partition(&mut a);
-        let k = q + 1;
-        if i == k {
-            (&a[q]).clone()
-        } else if i < k {
-            randomized_select(a.slice_mut(s![0..q]).to_owned(), i)
-        } else {
-            randomized_select(a.slice_mut(s![(q+1)..n]).to_owned(), i - k)
-        }
-    }
-}
+// fn randomized_select<A>(mut a: Array1<A>, i: usize) -> (A, Array1<A>)
+//     where A: Ord + Clone + Copy
+// {
+//     let n = a.len();
+//     if n == 0 {
+//        ((&a[0]).clone(), a)
+//     } else {
+//         let (q, a) = randomized_partition(&mut a);
+//         let k = q + 1;
+//         if i == k {
+//             ((&a[q]).clone(), *a)
+//         } else {
+//             let (mut pre, mut after) = a.view_mut().split_at(Axis(0), q).;
+//             let mut x;
+//             if i < k {
+//                 let (z, y) = randomized_select(pre.to_owned(), i);
+//                 x = z;
+//                 pre = y.view_mut().to_owned();
+//             } else {
+//                 let (z, y) = randomized_select(after.to_owned(), i - k);
+//                 x = z;
+//                 after = y.view_mut().to_owned();
+//             }
+//             (x, stack![Axis(0), pre, after])
+//         }
+//     }
+// }
 
-fn randomized_partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
+fn randomized_partition<A>(a: &mut Array1<A>) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
@@ -180,7 +189,7 @@ fn randomized_partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
     partition(a)
 }
 
-fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
+fn partition<A>(a: &mut Array1<A>) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
@@ -193,26 +202,33 @@ fn partition<A>(a: &mut Array1<A>) -> (usize, &mut Array1<A>)
         }
     }
     a.swap(i, n-1);
-    (i, a)
+    i
 }
 
 #[test]
 fn test_partition() {
-    let mut a = arr1(&[1, 3, 2, 10]);
-    let (j, a) = partition(&mut a);
-    assert_eq!(j, 3);
+    let mut a = arr1(&[1, 3, 2, 10, 10]);
+    let j = partition(&mut a);
+    assert_eq!(j, 4);
+    for i in 0..j {
+        assert!(a[i] <= a[j]);
+    }
     let mut a = arr1(&[2, 3, 4, 1]);
-    let (j, a) = partition(&mut a);
+    let j = partition(&mut a);
     assert_eq!(j, 0);
+    let n = a.len();
+    for i in j+1..n {
+        assert!(a[i] > a[j]);
+    }
 }
 
-#[test]
-fn test_randomized_select() {
-    let a = arr1(&[1, 3, 2, 10]);
-    let j = randomized_select(a.clone(), 2);
-    assert_eq!(j, 2);
-    let j = randomized_select(a.clone(), 1);
-    assert_eq!(j, 1);
-    let j = randomized_select(a.clone(), 3);
-    assert_eq!(j, 3);
-}
+// #[test]
+// fn test_randomized_select() {
+//     let a = arr1(&[1, 3, 2, 10]);
+//     let j = randomized_select(a.clone(), 2);
+//     assert_eq!(j, 2);
+//     let j = randomized_select(a.clone(), 1);
+//     assert_eq!(j, 1);
+//     let j = randomized_select(a.clone(), 3);
+//     assert_eq!(j, 3);
+// }

From 93591348ee0095ba6d9479635e49cd435daeb03f Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 23:30:05 +0100
Subject: [PATCH 14/93] Restored previous version of randomized_select.

---
 src/numeric/impl_numeric.rs | 45 +++++++++++++++----------------------
 1 file changed, 18 insertions(+), 27 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 977ba8005..55bebd5d6 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,33 +151,24 @@ impl<A, S, D> ArrayBase<S, D>
     }
 }
 
-// fn randomized_select<A>(mut a: Array1<A>, i: usize) -> (A, Array1<A>)
-//     where A: Ord + Clone + Copy
-// {
-//     let n = a.len();
-//     if n == 0 {
-//        ((&a[0]).clone(), a)
-//     } else {
-//         let (q, a) = randomized_partition(&mut a);
-//         let k = q + 1;
-//         if i == k {
-//             ((&a[q]).clone(), *a)
-//         } else {
-//             let (mut pre, mut after) = a.view_mut().split_at(Axis(0), q).;
-//             let mut x;
-//             if i < k {
-//                 let (z, y) = randomized_select(pre.to_owned(), i);
-//                 x = z;
-//                 pre = y.view_mut().to_owned();
-//             } else {
-//                 let (z, y) = randomized_select(after.to_owned(), i - k);
-//                 x = z;
-//                 after = y.view_mut().to_owned();
-//             }
-//             (x, stack![Axis(0), pre, after])
-//         }
-//     }
-// }
+fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
+    where A: Ord + Clone
+{
+    let n = a.len();
+    if n == 0 {
+        (&a[0]).clone()
+    } else {
+        let q = randomized_partition(&mut a);
+        let k = q + 1;
+        if i == k {
+            (&a[q]).clone()
+        } else if i < k {
+            randomized_select(a.slice_mut(s![0..q]).to_owned(), i)
+        } else {
+            randomized_select(a.slice_mut(s![(q+1)..n]).to_owned(), i - k)
+        }
+    }
+}
 
 fn randomized_partition<A>(a: &mut Array1<A>) -> usize
     where A: Ord + Clone

From 43d765edc0291f2d21ae6b24b7df876891581fcf Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 15 May 2018 23:46:02 +0100
Subject: [PATCH 15/93] Modified all function signatures to accept mutable
 views instead of owned arrays. Adapted tests as well.

---
 src/numeric/impl_numeric.rs | 36 ++++++++++++++++++------------------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 55bebd5d6..38d1953a7 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,26 +151,26 @@ impl<A, S, D> ArrayBase<S, D>
     }
 }
 
-fn randomized_select<A>(mut a: Array1<A>, i: usize) -> A
+fn randomized_select<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
     where A: Ord + Clone
 {
     let n = a.len();
     if n == 0 {
         (&a[0]).clone()
     } else {
-        let q = randomized_partition(&mut a);
+        let q = randomized_partition(&mut a.view_mut());
         let k = q + 1;
         if i == k {
             (&a[q]).clone()
         } else if i < k {
-            randomized_select(a.slice_mut(s![0..q]).to_owned(), i)
+            randomized_select(a.slice_mut(s![0..q]), i)
         } else {
-            randomized_select(a.slice_mut(s![(q+1)..n]).to_owned(), i - k)
+            randomized_select(a.slice_mut(s![(q+1)..n]), i - k)
         }
     }
 }
 
-fn randomized_partition<A>(a: &mut Array1<A>) -> usize
+fn randomized_partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
@@ -180,7 +180,7 @@ fn randomized_partition<A>(a: &mut Array1<A>) -> usize
     partition(a)
 }
 
-fn partition<A>(a: &mut Array1<A>) -> usize
+fn partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
@@ -199,13 +199,13 @@ fn partition<A>(a: &mut Array1<A>) -> usize
 #[test]
 fn test_partition() {
     let mut a = arr1(&[1, 3, 2, 10, 10]);
-    let j = partition(&mut a);
+    let j = partition(&mut a.view_mut());
     assert_eq!(j, 4);
     for i in 0..j {
         assert!(a[i] <= a[j]);
     }
     let mut a = arr1(&[2, 3, 4, 1]);
-    let j = partition(&mut a);
+    let j = partition(&mut a.view_mut());
     assert_eq!(j, 0);
     let n = a.len();
     for i in j+1..n {
@@ -213,13 +213,13 @@ fn test_partition() {
     }
 }
 
-// #[test]
-// fn test_randomized_select() {
-//     let a = arr1(&[1, 3, 2, 10]);
-//     let j = randomized_select(a.clone(), 2);
-//     assert_eq!(j, 2);
-//     let j = randomized_select(a.clone(), 1);
-//     assert_eq!(j, 1);
-//     let j = randomized_select(a.clone(), 3);
-//     assert_eq!(j, 3);
-// }
+#[test]
+fn test_randomized_select() {
+    let a = arr1(&[1, 3, 2, 10]);
+    let j = randomized_select(a.clone().view_mut(), 2);
+    assert_eq!(j, 2);
+    let j = randomized_select(a.clone().view_mut(), 1);
+    assert_eq!(j, 1);
+    let j = randomized_select(a.clone().view_mut(), 3);
+    assert_eq!(j, 3);
+}

From 0dec6992314ef137b8689d8ffc086ba0cd90439e Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 20:00:31 +0200
Subject: [PATCH 16/93] Added macro_use to get azip! in numerics.

---
 src/lib.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/lib.rs b/src/lib.rs
index 3e918d09b..ac3414239 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -171,6 +171,7 @@ mod numeric_util;
 mod error;
 mod shape_builder;
 mod stacking;
+#[macro_use]
 mod zip;
 
 pub use zip::{

From 9315413a647ee600d5835b80fd89ed0e6ae891e8 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 20:38:11 +0200
Subject: [PATCH 17/93] Added macro_use to macrozip to get azip! in numeric.

---
 src/zip/mod.rs | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/zip/mod.rs b/src/zip/mod.rs
index 3991a7e79..8ffdbae0a 100644
--- a/src/zip/mod.rs
+++ b/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::*;
@@ -318,7 +319,7 @@ impl<'a, A, D: Dimension> NdProducer for ArrayView<'a, A, D>
 
     #[inline(always)]
     fn contiguous_stride(&self) -> Self::Stride { 1 }
-    
+
     #[doc(hidden)]
     fn split_at(self, axis: Axis, index: usize) -> (Self, Self) {
         self.split_at(axis, index)

From 212d1b79bc0ea8a80296f9a861605078e863a17f Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 20:46:01 +0200
Subject: [PATCH 18/93] Added an implementation of percentile_axis that
 requires a mutable reference to the array (thus renamed _mut)

---
 src/numeric/impl_numeric.rs | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 38d1953a7..68eed9097 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -120,12 +120,19 @@ impl<A, S, D> ArrayBase<S, D>
 
 
     /// Return the qth percentile of the data along the specified axis.
-    pub fn percentile_axis(&self, axis: Axis, q: f32) -> Array<A, D::Smaller>
+    pub fn percentile_axis_mut(&mut self, axis: Axis, q: f32) -> Array<A, D::Smaller>
         where D: RemoveAxis,
+              A: Ord + Clone + Zero,
+              S: DataMut,
     {
         let n = self.len_of(axis);
         let i = ((n as f32) * q).floor() as usize;
-        unimplemented!()
+        let mapping = |x| randomized_select(x, i);
+        let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
+        azip!(mut lane (self.lanes_mut(axis)), mut out in {
+            *out = mapping(lane);
+        });
+        out
     }
 
     /// Return `true` if the arrays' elementwise differences are all within

From ccc8b042a852a77226407a8e930ebee3295d3e31 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 20:58:59 +0200
Subject: [PATCH 19/93] We need to increment the index by 1 to get the correct
 result

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 68eed9097..05c5dffc8 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -127,7 +127,7 @@ impl<A, S, D> ArrayBase<S, D>
     {
         let n = self.len_of(axis);
         let i = ((n as f32) * q).floor() as usize;
-        let mapping = |x| randomized_select(x, i);
+        let mapping = |x| randomized_select(x, i+1);
         let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
         azip!(mut lane (self.lanes_mut(axis)), mut out in {
             *out = mapping(lane);

From 54f227a1cc4066e92bc3213dfe07bcbbcc721f2e Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 21:11:03 +0200
Subject: [PATCH 20/93] Mapping passes i instead of i+1 to randomized_select,
 but we use ceil to get i from q. This provides the right behaviour when
 calling the function on arrays with odd and even number of values when
 computing the median value. Added some test for percentile_axis_mut.

---
 src/numeric/impl_numeric.rs | 27 +++++++++++++++++++++++++--
 1 file changed, 25 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 05c5dffc8..99ba537c6 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -126,8 +126,8 @@ impl<A, S, D> ArrayBase<S, D>
               S: DataMut,
     {
         let n = self.len_of(axis);
-        let i = ((n as f32) * q).floor() as usize;
-        let mapping = |x| randomized_select(x, i+1);
+        let i = ((n as f32) * q).ceil() as usize;
+        let mapping = |x| randomized_select(x, i);
         let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
         azip!(mut lane (self.lanes_mut(axis)), mut out in {
             *out = mapping(lane);
@@ -230,3 +230,26 @@ fn test_randomized_select() {
     let j = randomized_select(a.clone().view_mut(), 3);
     assert_eq!(j, 3);
 }
+
+#[test]
+fn test_percentile_axis_mut() {
+    let mut a = arr2(
+        &[
+        [1, 3, 2, 10],
+        [2, 4, 3, 11],
+        [3, 5, 6, 12]
+        ]
+    );
+    let p = a.percentile_axis_mut(Axis(0), 0.5);
+    assert!(p == a.subview(Axis(0), 1));
+    let mut b = arr2(
+        &[
+        [1, 3, 2, 10],
+        [2, 4, 3, 11],
+        [3, 5, 6, 12],
+        [4, 6, 7, 13]
+        ]
+    );
+    let q = b.percentile_axis_mut(Axis(0), 0.5);
+    assert!(q == b.subview(Axis(0), 1));
+}

From 461e258a4f3047ddef5d73509a8b645e308ced42 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 21:13:55 +0200
Subject: [PATCH 21/93] Added some documentation for percentile_axis_mut.

---
 src/numeric/impl_numeric.rs | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 99ba537c6..c7dc17960 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -118,8 +118,11 @@ impl<A, S, D> ArrayBase<S, D>
         sum / &aview0(&cnt)
     }
 
-
-    /// Return the qth percentile of the data along the specified axis.
+    /// Return the qth percentile of the data along the specified axis
+    /// `q` needs to be a float between 0 and 1, included.
+    ///
+    /// **Panics** if `axis` is out of bounds.
+    /// **Panics** if `q` is strictly smaller than 0 or strictly bigger than 1.
     pub fn percentile_axis_mut(&mut self, axis: Axis, q: f32) -> Array<A, D::Smaller>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,

From 184d65d7986a862c8b5c0c070262379b16c0ac3c Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 21:15:17 +0200
Subject: [PATCH 22/93] Added an assert to prevent invalid q values.

---
 src/numeric/impl_numeric.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index c7dc17960..903565a4a 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -128,6 +128,7 @@ impl<A, S, D> ArrayBase<S, D>
               A: Ord + Clone + Zero,
               S: DataMut,
     {
+        assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = ((n as f32) * q).ceil() as usize;
         let mapping = |x| randomized_select(x, i);

From e290bff64ed611095f04e6a6de1b73261a318694 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 21:31:12 +0200
Subject: [PATCH 23/93] Changed end of recursion condition - it needs to check
 for n==1 not n==0. Implemented a test to check that percentile 0 returns the
 array minimum

---
 src/numeric/impl_numeric.rs | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 903565a4a..c6a8ef309 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -166,7 +166,7 @@ fn randomized_select<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
     where A: Ord + Clone
 {
     let n = a.len();
-    if n == 0 {
+    if n == 1 {
         (&a[0]).clone()
     } else {
         let q = randomized_partition(&mut a.view_mut());
@@ -257,3 +257,11 @@ fn test_percentile_axis_mut() {
     let q = b.percentile_axis_mut(Axis(0), 0.5);
     assert!(q == b.subview(Axis(0), 1));
 }
+
+#[test]
+fn test_percentile_axis_mut_to_get_minimum() {
+    let mut b = arr2(&[[1, 3, 22, 10]]);
+    let q = b.percentile_axis_mut(Axis(1), 0.);
+    println!("{:?}", q);
+    assert!(q == arr1(&[1]));
+}

From 911a656497d9401391b1cd569d25539e6f3f093b Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 21:35:20 +0200
Subject: [PATCH 24/93] Added a new test for the maximum, as well as testing it
 on arr1.

---
 src/numeric/impl_numeric.rs | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index c6a8ef309..a433df624 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -262,6 +262,12 @@ fn test_percentile_axis_mut() {
 fn test_percentile_axis_mut_to_get_minimum() {
     let mut b = arr2(&[[1, 3, 22, 10]]);
     let q = b.percentile_axis_mut(Axis(1), 0.);
-    println!("{:?}", q);
     assert!(q == arr1(&[1]));
 }
+
+#[test]
+fn test_percentile_axis_mut_to_get_maximum() {
+    let mut b = arr1(&[1, 3, 22, 10]);
+    let q = b.percentile_axis_mut(Axis(0), 1.);
+    assert!(q == arr0(22));
+}

From c57fb2511bfb90e26794cad79d6472215effa278 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Mon, 21 May 2018 21:36:26 +0200
Subject: [PATCH 25/93] Split test into two more meaningful tests.

---
 src/numeric/impl_numeric.rs | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index a433df624..f70eb397d 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -236,7 +236,7 @@ fn test_randomized_select() {
 }
 
 #[test]
-fn test_percentile_axis_mut() {
+fn test_percentile_axis_mut_with_odd_axis_length() {
     let mut a = arr2(
         &[
         [1, 3, 2, 10],
@@ -246,6 +246,10 @@ fn test_percentile_axis_mut() {
     );
     let p = a.percentile_axis_mut(Axis(0), 0.5);
     assert!(p == a.subview(Axis(0), 1));
+}
+
+#[test]
+fn test_percentile_axis_mut_with_even_axis_length() {
     let mut b = arr2(
         &[
         [1, 3, 2, 10],

From 2107f565b743439df0f2177ef27262f451bd888c Mon Sep 17 00:00:00 2001
From: LukeMathWalker <LukeMathWalker@users.noreply.github.com>
Date: Thu, 24 May 2018 09:50:21 +0200
Subject: [PATCH 26/93] Added more documentation

---
 src/numeric/impl_numeric.rs | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index f70eb397d..ab957f103 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -118,8 +118,23 @@ impl<A, S, D> ArrayBase<S, D>
         sum / &aview0(&cnt)
     }
 
-    /// Return the qth percentile of the data along the specified axis
-    /// `q` needs to be a float between 0 and 1, included.
+    /// Return the qth percentile of the data along the specified axis.
+    /// `q` needs to be a float between 0 and 1, bounds included.
+    /// For instance:
+    /// - `q=0.` returns the minimum along each 1-dimensional lane;
+    /// - `q=0.5` returns the median along each 1-dimensional lane;
+    /// - `q=1.` returns the maximum along each 1-dimensional lane.
+    ///
+    /// If the desired percentile lies between two data points we always
+    /// return the smaller data point.
+    ///
+    /// The array is shuffled **in place** along each 1-dimensional lane in
+    /// order to produce the required percentile without allocating a copy
+    /// of the original array. Each 1-dimensional lane is shuffled indipendently
+    /// from the others.
+    ///
+    /// The algorithm asymptotic complexity in the worst case is O(m) where
+    /// m is the number of elements in the array.
     ///
     /// **Panics** if `axis` is out of bounds.
     /// **Panics** if `q` is strictly smaller than 0 or strictly bigger than 1.

From 07bc655fb79352f229717668a582dd5ff2be1db2 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Wed, 30 May 2018 19:07:09 +0100
Subject: [PATCH 27/93] Updated to version 0.5.0 of rand.

---
 Cargo.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Cargo.toml b/Cargo.toml
index 607d4fb71..55148861d 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,7 +28,7 @@ num-traits = { version = "0.1.32", default-features = false }
 num-complex = { version = "0.1.32", default-features = false }
 rustc-serialize = { version = "0.3.20", optional = true }
 itertools = { version = "0.7.0", default-features = false }
-rand = { version = "0.5.0-pre.1" }
+rand = { version = "0.5.0" }
 
 # Use via the `blas` crate feature!
 cblas-sys = { version = "0.1.4", optional = true, default-features = false }

From 891fe95f53795548a3058487bbf026f3e8ee97b0 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Wed, 30 May 2018 19:17:37 +0100
Subject: [PATCH 28/93] Updated to use the proper way of generating random
 integers.

---
 src/numeric/impl_numeric.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index ab957f103..be6dc7ced 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -13,7 +13,7 @@ use itertools::free::enumerate;
 use imp_prelude::*;
 use numeric_util;
 
-use rand::distributions::Uniform;
+use rand::prelude::*;
 use rand::thread_rng;
 
 use {
@@ -201,7 +201,7 @@ fn randomized_partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>) -> usize
 {
     let n = a.len();
     let mut rng = thread_rng();
-    let i: usize = Uniform::sample_single(0, n, &mut rng);
+    let i: usize = rng.gen_range(0, n);
     a.swap(i, n-1);
     partition(a)
 }

From 540c495b9841e75b420005f822a4116c8f36e08d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Wed, 30 May 2018 19:26:06 +0100
Subject: [PATCH 29/93] Updated docs.

---
 src/numeric/impl_numeric.rs | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index be6dc7ced..af23b8b9a 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -120,13 +120,16 @@ impl<A, S, D> ArrayBase<S, D>
 
     /// Return the qth percentile of the data along the specified axis.
     /// `q` needs to be a float between 0 and 1, bounds included.
-    /// For instance:
+    /// The qth percentile for a 1-dimensional lane of length N is defined
+    /// as the element that would be indexed as `Nq` if the lane were to be sorted
+    /// in increasing order. If `Nq` is not an integer the desired percentile
+    /// lies between two data points: we always return the smaller data point.
+    ///
+    /// Some examples:
     /// - `q=0.` returns the minimum along each 1-dimensional lane;
     /// - `q=0.5` returns the median along each 1-dimensional lane;
     /// - `q=1.` returns the maximum along each 1-dimensional lane.
-    ///
-    /// If the desired percentile lies between two data points we always
-    /// return the smaller data point.
+    /// (`q=0` and `q=1` are considered improper percentiles)
     ///
     /// The array is shuffled **in place** along each 1-dimensional lane in
     /// order to produce the required percentile without allocating a copy
@@ -136,8 +139,8 @@ impl<A, S, D> ArrayBase<S, D>
     /// The algorithm asymptotic complexity in the worst case is O(m) where
     /// m is the number of elements in the array.
     ///
-    /// **Panics** if `axis` is out of bounds.
-    /// **Panics** if `q` is strictly smaller than 0 or strictly bigger than 1.
+    /// **Panics** if `axis` is out of bounds, if `q` is strictly smaller
+    /// than 0 or strictly bigger than 1.
     pub fn percentile_axis_mut(&mut self, axis: Axis, q: f32) -> Array<A, D::Smaller>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,

From e22f889531b245644219907ff35f6dcb65de0930 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Wed, 30 May 2018 19:27:53 +0100
Subject: [PATCH 30/93] Updated docs.

---
 src/numeric/impl_numeric.rs | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index af23b8b9a..679af447e 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -135,6 +135,8 @@ impl<A, S, D> ArrayBase<S, D>
     /// order to produce the required percentile without allocating a copy
     /// of the original array. Each 1-dimensional lane is shuffled indipendently
     /// from the others.
+    /// No assumptions should be made on the ordering of the array elements
+    /// after this computation.
     ///
     /// The algorithm asymptotic complexity in the worst case is O(m) where
     /// m is the number of elements in the array.

From 54b1878f2e0d0ac05a0a85568938a380d5c509d0 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 14:38:13 +0100
Subject: [PATCH 31/93]  Changing tests for randomized_select - we want to get
 i to be the index  of the element returned if the array we pass were to be
 sorted.

---
 src/numeric/impl_numeric.rs | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 679af447e..35d8b8349 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -248,11 +248,11 @@ fn test_partition() {
 fn test_randomized_select() {
     let a = arr1(&[1, 3, 2, 10]);
     let j = randomized_select(a.clone().view_mut(), 2);
-    assert_eq!(j, 2);
+    assert_eq!(j, 3);
     let j = randomized_select(a.clone().view_mut(), 1);
-    assert_eq!(j, 1);
+    assert_eq!(j, 2);
     let j = randomized_select(a.clone().view_mut(), 3);
-    assert_eq!(j, 3);
+    assert_eq!(j, 10);
 }
 
 #[test]

From a54b6966ccb90bf5298d8010009574f06cf24283 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 14:58:40 +0100
Subject: [PATCH 32/93] Renamed randomized_select: now nth_mut. Modified all
 functions to make it easy to switch the strategy used to choose the
 pivot_index, in order to support different selection algorithms if needed

---
 src/numeric/impl_numeric.rs | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 35d8b8349..5cbce7f7c 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,7 +151,7 @@ impl<A, S, D> ArrayBase<S, D>
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = ((n as f32) * q).ceil() as usize;
-        let mapping = |x| randomized_select(x, i);
+        let mapping = |x| nth_mut(x, i);
         let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
         azip!(mut lane (self.lanes_mut(axis)), mut out in {
             *out = mapping(lane);
@@ -182,39 +182,37 @@ impl<A, S, D> ArrayBase<S, D>
     }
 }
 
-fn randomized_select<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
+fn nth_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
     where A: Ord + Clone
 {
     let n = a.len();
     if n == 1 {
         (&a[0]).clone()
     } else {
-        let q = randomized_partition(&mut a.view_mut());
+        let pivot_index = random_pivot(n);
+        let q = partition(&mut a.view_mut(), pivot_index);
         let k = q + 1;
         if i == k {
             (&a[q]).clone()
         } else if i < k {
-            randomized_select(a.slice_mut(s![0..q]), i)
+            nth_mut(a.slice_mut(s![0..q]), i)
         } else {
-            randomized_select(a.slice_mut(s![(q+1)..n]), i - k)
+            nth_mut(a.slice_mut(s![(q+1)..n]), i - k)
         }
     }
 }
 
-fn randomized_partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>) -> usize
-    where A: Ord + Clone
+fn random_pivot(n: usize) -> usize
 {
-    let n = a.len();
     let mut rng = thread_rng();
-    let i: usize = rng.gen_range(0, n);
-    a.swap(i, n-1);
-    partition(a)
+    rng.gen_range(0, n)
 }
 
-fn partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>) -> usize
+fn partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>, pivot_index: usize) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
+    a.swap(pivot_index, n-1);
     let x = (&a[n-1]).clone();
     let mut i = 0;
     for j in 0..n-1 {

From 337d3891fe63eb81aad3df7131fd726a79cfc2ee Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:08:59 +0100
Subject: [PATCH 33/93] nth_mut renamed ith_mut, to have coherence with
 argument naming strategy. Provided documentation for ith_mut

---
 src/numeric/impl_numeric.rs | 22 ++++++++++++++++++----
 1 file changed, 18 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 5cbce7f7c..193ff1c99 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,7 +151,7 @@ impl<A, S, D> ArrayBase<S, D>
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = ((n as f32) * q).ceil() as usize;
-        let mapping = |x| nth_mut(x, i);
+        let mapping = |x| ith_mut(x, i);
         let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
         azip!(mut lane (self.lanes_mut(axis)), mut out in {
             *out = mapping(lane);
@@ -182,7 +182,21 @@ impl<A, S, D> ArrayBase<S, D>
     }
 }
 
-fn nth_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
+/// Return the `i`-th element of `a` if `a` were to be a 1-dimensional
+/// array sorted in increasing order.
+///
+/// `a` is shuffled **in place** to retrieve the desired element:
+/// no copy of the array is allocated.
+/// No assumptions should be made on the ordering of `a` elements
+/// after this computation.
+///
+/// Complexity:
+/// - average case: O(n);
+/// - worst case: O(n^2);
+/// where n is the number of elements in `a`.
+///
+/// **Panics** if `i` is greater than or equal to n.
+fn ith_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
     where A: Ord + Clone
 {
     let n = a.len();
@@ -195,9 +209,9 @@ fn nth_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
         if i == k {
             (&a[q]).clone()
         } else if i < k {
-            nth_mut(a.slice_mut(s![0..q]), i)
+            ith_mut(a.slice_mut(s![0..q]), i)
         } else {
-            nth_mut(a.slice_mut(s![(q+1)..n]), i - k)
+            ith_mut(a.slice_mut(s![(q+1)..n]), i - k)
         }
     }
 }

From 14f160547ffcf6715e138d1c46f9e776688060ec Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:17:35 +0100
Subject: [PATCH 34/93] Improved variable naming.

---
 src/numeric/impl_numeric.rs | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 193ff1c99..7f4c495ff 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -204,14 +204,14 @@ fn ith_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
         (&a[0]).clone()
     } else {
         let pivot_index = random_pivot(n);
-        let q = partition(&mut a.view_mut(), pivot_index);
-        let k = q + 1;
+        let partition_index = partition(&mut a.view_mut(), pivot_index);
+        let k = partition_index + 1;
         if i == k {
-            (&a[q]).clone()
+            (&a[partition_index]).clone()
         } else if i < k {
-            ith_mut(a.slice_mut(s![0..q]), i)
+            ith_mut(a.slice_mut(s![0..partition_index]), i)
         } else {
-            ith_mut(a.slice_mut(s![(q+1)..n]), i - k)
+            ith_mut(a.slice_mut(s![(partition_index+1)..n]), i - k)
         }
     }
 }

From 10603a29b78a057ed4604315578ab7abccb08916 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:17:56 +0100
Subject: [PATCH 35/93] partition renamed partition_mut.

---
 src/numeric/impl_numeric.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 7f4c495ff..19d75288c 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -204,7 +204,7 @@ fn ith_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
         (&a[0]).clone()
     } else {
         let pivot_index = random_pivot(n);
-        let partition_index = partition(&mut a.view_mut(), pivot_index);
+        let partition_index = partition_mut(&mut a.view_mut(), pivot_index);
         let k = partition_index + 1;
         if i == k {
             (&a[partition_index]).clone()
@@ -222,7 +222,7 @@ fn random_pivot(n: usize) -> usize
     rng.gen_range(0, n)
 }
 
-fn partition<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>, pivot_index: usize) -> usize
+fn partition_mut<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>, pivot_index: usize) -> usize
     where A: Ord + Clone
 {
     let n = a.len();

From cbcf69aacb99e1a3d317f5dae0a54ae36f4773d0 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:27:19 +0100
Subject: [PATCH 36/93] Added documentation to partition. Renamed variable for
 better readability.

---
 src/numeric/impl_numeric.rs | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 19d75288c..080c8ee17 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -222,11 +222,26 @@ fn random_pivot(n: usize) -> usize
     rng.gen_range(0, n)
 }
 
-fn partition_mut<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>, pivot_index: usize) -> usize
+/// Return the index of `a[partition_index`]` if `a` were to be sorted
+/// in increasing order.
+/// `a` elements are rearranged in such a way that `a[partition_index]`
+/// is in the position it would be in an array sorted in increasing order.
+/// All elements smaller than `a[partition_index]` are moved to its
+/// left and all elements equal or greater than `a[partition_index]`
+/// are moved to its right.
+/// The ordering of the elements in the two partitions is undefined.
+///
+/// `a` is shuffled **in place** to operate the desired partition:
+/// no copy of the array is allocated.
+///
+/// Complexity: O(n), where n is the number of elements in `a`.
+///
+/// **Panics** if `partition_index` is greater than or equal to n.
+fn partition_mut<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>, partition_index: usize) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
-    a.swap(pivot_index, n-1);
+    a.swap(partition_index, n-1);
     let x = (&a[n-1]).clone();
     let mut i = 0;
     for j in 0..n-1 {

From 4d71df751c2cb7cde0901f9f47f4ac84d2c7ab59 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:28:19 +0100
Subject: [PATCH 37/93] Shortened function signature.

---
 src/numeric/impl_numeric.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 080c8ee17..7ff65ffab 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -196,7 +196,7 @@ impl<A, S, D> ArrayBase<S, D>
 /// where n is the number of elements in `a`.
 ///
 /// **Panics** if `i` is greater than or equal to n.
-fn ith_mut<A>(mut a: ArrayViewMut<A, Dim<[Ix; 1]>>, i: usize) -> A
+fn ith_mut<A>(mut a: ArrayViewMut1<A>, i: usize) -> A
     where A: Ord + Clone
 {
     let n = a.len();
@@ -237,7 +237,7 @@ fn random_pivot(n: usize) -> usize
 /// Complexity: O(n), where n is the number of elements in `a`.
 ///
 /// **Panics** if `partition_index` is greater than or equal to n.
-fn partition_mut<A>(a: &mut ArrayViewMut<A, Dim<[Ix; 1]>>, partition_index: usize) -> usize
+fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
     where A: Ord + Clone
 {
     let n = a.len();

From 6e06f13e923fbbc6cc97c57cca709440d448e313 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:32:20 +0100
Subject: [PATCH 38/93] Using more expressive variable names in partition_mut.

---
 src/numeric/impl_numeric.rs | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 7ff65ffab..d17ae160c 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -241,17 +241,17 @@ fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
     where A: Ord + Clone
 {
     let n = a.len();
+    let partition_value = (&a[partition_index]).clone();
     a.swap(partition_index, n-1);
-    let x = (&a[n-1]).clone();
-    let mut i = 0;
+    let mut partition_boundary_index = 0;
     for j in 0..n-1 {
-        if a[j] <= x {
-            a.swap(i, j);
-            i += 1;
+        if a[j] <= partition_value {
+            a.swap(partition_boundary_index, j);
+            partition_boundary_index += 1;
         }
     }
-    a.swap(i, n-1);
-    i
+    a.swap(partition_boundary_index, n-1);
+    partition_boundary_index
 }
 
 #[test]

From 31f0f28d61d95096a0f85bc70146c7a857c02f5e Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:39:45 +0100
Subject: [PATCH 39/93] Removed unnecessary auxiliary variable, to achieve more
 meaningful naming.

---
 src/numeric/impl_numeric.rs | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index d17ae160c..48cd8816e 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -205,13 +205,12 @@ fn ith_mut<A>(mut a: ArrayViewMut1<A>, i: usize) -> A
     } else {
         let pivot_index = random_pivot(n);
         let partition_index = partition_mut(&mut a.view_mut(), pivot_index);
-        let k = partition_index + 1;
-        if i == k {
+        if i == partition_index {
             (&a[partition_index]).clone()
-        } else if i < k {
+        } else if i < partition_index {
             ith_mut(a.slice_mut(s![0..partition_index]), i)
         } else {
-            ith_mut(a.slice_mut(s![(partition_index+1)..n]), i - k)
+            ith_mut(a.slice_mut(s![(partition_index+1)..n]), i - partition_index - 1)
         }
     }
 }

From 77629cad02e1c981e7ebefea78b05dc3700ba762 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:42:24 +0100
Subject: [PATCH 40/93] Updated test_partition into test_partition_mut.

---
 src/numeric/impl_numeric.rs | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 48cd8816e..b7cd0e49d 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -254,15 +254,17 @@ fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
 }
 
 #[test]
-fn test_partition() {
+fn test_partition_mut() {
     let mut a = arr1(&[1, 3, 2, 10, 10]);
-    let j = partition(&mut a.view_mut());
-    assert_eq!(j, 4);
+    let n = a.len();
+    let j = partition(&mut a.view_mut(), n-1);
+    assert_eq!(j, 3);
     for i in 0..j {
         assert!(a[i] <= a[j]);
     }
     let mut a = arr1(&[2, 3, 4, 1]);
-    let j = partition(&mut a.view_mut());
+    let n = a.len();
+    let j = partition(&mut a.view_mut(), n-1);
     assert_eq!(j, 0);
     let n = a.len();
     for i in j+1..n {

From eee88ce6b4bb1730d4bbafbf8aa3dfbc0eec1ce9 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:43:16 +0100
Subject: [PATCH 41/93] Updated test_randomized_select into test_ith_mut.

---
 src/numeric/impl_numeric.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index b7cd0e49d..b017c290d 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -273,13 +273,13 @@ fn test_partition_mut() {
 }
 
 #[test]
-fn test_randomized_select() {
+fn test_ith_mut() {
     let a = arr1(&[1, 3, 2, 10]);
-    let j = randomized_select(a.clone().view_mut(), 2);
+    let j = ith_mut(a.clone().view_mut(), 2);
     assert_eq!(j, 3);
-    let j = randomized_select(a.clone().view_mut(), 1);
+    let j = ith_mut(a.clone().view_mut(), 1);
     assert_eq!(j, 2);
-    let j = randomized_select(a.clone().view_mut(), 3);
+    let j = ith_mut(a.clone().view_mut(), 3);
     assert_eq!(j, 10);
 }
 

From 5a29f4796fa39d7591374477f0400ebba7f49b9c Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:44:44 +0100
Subject: [PATCH 42/93] Fixed bug in test_partition_mut.

---
 src/numeric/impl_numeric.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index b017c290d..42000e36e 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -257,14 +257,14 @@ fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
 fn test_partition_mut() {
     let mut a = arr1(&[1, 3, 2, 10, 10]);
     let n = a.len();
-    let j = partition(&mut a.view_mut(), n-1);
+    let j = partition_mut(&mut a.view_mut(), n-1);
     assert_eq!(j, 3);
     for i in 0..j {
         assert!(a[i] <= a[j]);
     }
     let mut a = arr1(&[2, 3, 4, 1]);
     let n = a.len();
-    let j = partition(&mut a.view_mut(), n-1);
+    let j = partition_mut(&mut a.view_mut(), n-1);
     assert_eq!(j, 0);
     let n = a.len();
     for i in j+1..n {

From 310b53e4094b4825e3b0dd13f2049770f3861c0f Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 15:53:35 +0100
Subject: [PATCH 43/93] Fixed comparison to make documentation aligned to
 actual function behaviour.

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 42000e36e..4dfc9df67 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -244,7 +244,7 @@ fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
     a.swap(partition_index, n-1);
     let mut partition_boundary_index = 0;
     for j in 0..n-1 {
-        if a[j] <= partition_value {
+        if a[j] < partition_value {
             a.swap(partition_boundary_index, j);
             partition_boundary_index += 1;
         }

From 63516017a6f84f7bed21a2fcafd17fd6611e7450 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 16:06:45 +0100
Subject: [PATCH 44/93] Fixed edge case

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 4dfc9df67..d5471d9dd 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,7 +151,7 @@ impl<A, S, D> ArrayBase<S, D>
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = ((n as f32) * q).ceil() as usize;
-        let mapping = |x| ith_mut(x, i);
+        let mapping = |x| ith_mut(x, if i == 0 {0} else {i-1});
         let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
         azip!(mut lane (self.lanes_mut(axis)), mut out in {
             *out = mapping(lane);

From 5defbcc2c793f337056cd146160550c5e28e66fb Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 16:16:42 +0100
Subject: [PATCH 45/93] Moved tests to test module.

---
 src/numeric/impl_numeric.rs | 71 -------------------------------------
 tests/array.rs              | 71 +++++++++++++++++++++++++++++++++++++
 2 files changed, 71 insertions(+), 71 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index d5471d9dd..60de40a54 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -252,74 +252,3 @@ fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
     a.swap(partition_boundary_index, n-1);
     partition_boundary_index
 }
-
-#[test]
-fn test_partition_mut() {
-    let mut a = arr1(&[1, 3, 2, 10, 10]);
-    let n = a.len();
-    let j = partition_mut(&mut a.view_mut(), n-1);
-    assert_eq!(j, 3);
-    for i in 0..j {
-        assert!(a[i] <= a[j]);
-    }
-    let mut a = arr1(&[2, 3, 4, 1]);
-    let n = a.len();
-    let j = partition_mut(&mut a.view_mut(), n-1);
-    assert_eq!(j, 0);
-    let n = a.len();
-    for i in j+1..n {
-        assert!(a[i] > a[j]);
-    }
-}
-
-#[test]
-fn test_ith_mut() {
-    let a = arr1(&[1, 3, 2, 10]);
-    let j = ith_mut(a.clone().view_mut(), 2);
-    assert_eq!(j, 3);
-    let j = ith_mut(a.clone().view_mut(), 1);
-    assert_eq!(j, 2);
-    let j = ith_mut(a.clone().view_mut(), 3);
-    assert_eq!(j, 10);
-}
-
-#[test]
-fn test_percentile_axis_mut_with_odd_axis_length() {
-    let mut a = arr2(
-        &[
-        [1, 3, 2, 10],
-        [2, 4, 3, 11],
-        [3, 5, 6, 12]
-        ]
-    );
-    let p = a.percentile_axis_mut(Axis(0), 0.5);
-    assert!(p == a.subview(Axis(0), 1));
-}
-
-#[test]
-fn test_percentile_axis_mut_with_even_axis_length() {
-    let mut b = arr2(
-        &[
-        [1, 3, 2, 10],
-        [2, 4, 3, 11],
-        [3, 5, 6, 12],
-        [4, 6, 7, 13]
-        ]
-    );
-    let q = b.percentile_axis_mut(Axis(0), 0.5);
-    assert!(q == b.subview(Axis(0), 1));
-}
-
-#[test]
-fn test_percentile_axis_mut_to_get_minimum() {
-    let mut b = arr2(&[[1, 3, 22, 10]]);
-    let q = b.percentile_axis_mut(Axis(1), 0.);
-    assert!(q == arr1(&[1]));
-}
-
-#[test]
-fn test_percentile_axis_mut_to_get_maximum() {
-    let mut b = arr1(&[1, 3, 22, 10]);
-    let q = b.percentile_axis_mut(Axis(0), 1.);
-    assert!(q == arr0(22));
-}
diff --git a/tests/array.rs b/tests/array.rs
index 00257d159..c01408d12 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -1640,3 +1640,74 @@ fn array_macros() {
     let empty2: Array<f32, Ix2> = array![[]];
     assert_eq!(empty2, array![[]]);
 }
+
+#[test]
+fn test_partition_mut() {
+    let mut a = arr1(&[1, 3, 2, 10, 10]);
+    let n = a.len();
+    let j = partition_mut(&mut a.view_mut(), n-1);
+    assert_eq!(j, 3);
+    for i in 0..j {
+        assert!(a[i] <= a[j]);
+    }
+    let mut a = arr1(&[2, 3, 4, 1]);
+    let n = a.len();
+    let j = partition_mut(&mut a.view_mut(), n-1);
+    assert_eq!(j, 0);
+    let n = a.len();
+    for i in j+1..n {
+        assert!(a[i] > a[j]);
+    }
+}
+
+#[test]
+fn test_ith_mut() {
+    let a = arr1(&[1, 3, 2, 10]);
+    let j = ith_mut(a.clone().view_mut(), 2);
+    assert_eq!(j, 3);
+    let j = ith_mut(a.clone().view_mut(), 1);
+    assert_eq!(j, 2);
+    let j = ith_mut(a.clone().view_mut(), 3);
+    assert_eq!(j, 10);
+}
+
+#[test]
+fn test_percentile_axis_mut_with_odd_axis_length() {
+    let mut a = arr2(
+        &[
+        [1, 3, 2, 10],
+        [2, 4, 3, 11],
+        [3, 5, 6, 12]
+        ]
+    );
+    let p = a.percentile_axis_mut(Axis(0), 0.5);
+    assert!(p == a.subview(Axis(0), 1));
+}
+
+#[test]
+fn test_percentile_axis_mut_with_even_axis_length() {
+    let mut b = arr2(
+        &[
+        [1, 3, 2, 10],
+        [2, 4, 3, 11],
+        [3, 5, 6, 12],
+        [4, 6, 7, 13]
+        ]
+    );
+    let q = b.percentile_axis_mut(Axis(0), 0.5);
+    assert!(q == b.subview(Axis(0), 1));
+}
+
+#[test]
+fn test_percentile_axis_mut_to_get_minimum() {
+    let mut b = arr2(&[[1, 3, 22, 10]]);
+    let q = b.percentile_axis_mut(Axis(1), 0.);
+    assert!(q == arr1(&[1]));
+}
+
+#[test]
+fn test_percentile_axis_mut_to_get_maximum() {
+    let mut b = arr1(&[1, 3, 22, 10]);
+    let q = b.percentile_axis_mut(Axis(0), 1.);
+    assert!(q == arr0(22));
+}

From a3122e9babd5b7ff134e626c7e743ed4c4187cee Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 16:20:32 +0100
Subject: [PATCH 46/93] Fixed erroneous complexity indication. Added link to
 algorithm implementation.

---
 src/numeric/impl_numeric.rs | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 60de40a54..42bb354c8 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -138,8 +138,10 @@ impl<A, S, D> ArrayBase<S, D>
     /// No assumptions should be made on the ordering of the array elements
     /// after this computation.
     ///
-    /// The algorithm asymptotic complexity in the worst case is O(m) where
-    /// m is the number of elements in the array.
+    /// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
+    /// - average case: O(m);
+    /// - worst case: O(m^2);
+    /// where m is the number of elements in the array.
     ///
     /// **Panics** if `axis` is out of bounds, if `q` is strictly smaller
     /// than 0 or strictly bigger than 1.

From 243f12cb2b18f9b79193020f21173851d77a4427 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 16:21:13 +0100
Subject: [PATCH 47/93] Added link to algorithm description

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 42bb354c8..c98352fa0 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -192,7 +192,7 @@ impl<A, S, D> ArrayBase<S, D>
 /// No assumptions should be made on the ordering of `a` elements
 /// after this computation.
 ///
-/// Complexity:
+/// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
 /// - average case: O(n);
 /// - worst case: O(n^2);
 /// where n is the number of elements in `a`.

From 71ad41067f416ad6ca7918afcb3ad86a2e837b45 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 17:49:16 +0100
Subject: [PATCH 48/93] Fixed unclosed braces.

---
 src/numeric/impl_numeric.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 112f526fa..d7e6fcd6f 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -159,6 +159,7 @@ impl<A, S, D> ArrayBase<S, D>
             *out = mapping(lane);
         });
         out
+    }
 
     /// Return variance along `axis`.
     ///

From 517621c9dfa29b29dbbc34cf5d45fa30d06b4937 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 18:12:43 +0100
Subject: [PATCH 49/93] Properly reorganized functions to be used as method
 whenever possible

---
 src/impl_1d.rs              | 75 ++++++++++++++++++++++++++++++++++++
 src/numeric/impl_numeric.rs | 76 +------------------------------------
 2 files changed, 76 insertions(+), 75 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 4d013e7d5..5b5d1dad0 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -10,6 +10,9 @@
 //! Methods for one-dimensional arrays.
 use imp_prelude::*;
 
+use rand::prelude::*;
+use rand::thread_rng;
+
 impl<A, S> ArrayBase<S, Ix1>
     where S: Data<Elem=A>,
 {
@@ -23,5 +26,77 @@ impl<A, S> ArrayBase<S, Ix1>
             ::iterators::to_vec(self.iter().map(|x| x.clone()))
         }
     }
+
+    /// Return the `i`-th element of the array if it were to be a 1-dimensional
+    /// array sorted in increasing order.
+    ///
+    /// The array is shuffled **in place** to retrieve the desired element:
+    /// no copy of the array is allocated.
+    /// No assumptions should be made on the ordering of elements
+    /// after this computation.
+    ///
+    /// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
+    /// - average case: O(n);
+    /// - worst case: O(n^2);
+    /// where n is the number of elements in the array.
+    ///
+    /// **Panics** if `i` is greater than or equal to n.
+    pub fn ith_mut(&mut self, i: usize) -> A
+        where A: Ord + Clone,
+              S: DataMut,
+    {
+        let n = self.len();
+        if n == 1 {
+            (&self[0]).clone()
+        } else {
+            let pivot_index = random_pivot(n);
+            let partition_index = self.view_mut().partition_mut(pivot_index);
+            if i == partition_index {
+                (&self[partition_index]).clone()
+            } else if i < partition_index {
+                self.slice_mut(s![0..partition_index]).ith_mut(i)
+            } else {
+                self.slice_mut(s![(partition_index+1)..n]).ith_mut(i - partition_index - 1)
+            }
+        }
+    }
+
+    /// Return the index of `self[partition_index`]` if `self` were to be sorted
+    /// in increasing order.
+    /// `self` elements are rearranged in such a way that `self[partition_index]`
+    /// is in the position it would be in an array sorted in increasing order.
+    /// All elements smaller than `self[partition_index]` are moved to its
+    /// left and all elements equal or greater than `self[partition_index]`
+    /// are moved to its right.
+    /// The ordering of the elements in the two partitions is undefined.
+    ///
+    /// `self` is shuffled **in place** to operate the desired partition:
+    /// no copy of the array is allocated.
+    ///
+    /// Complexity: O(n), where n is the number of elements in the array.
+    ///
+    /// **Panics** if `partition_index` is greater than or equal to n.
+    pub fn partition_mut(&mut self, partition_index: usize) -> usize
+        where A: Ord + Clone,
+              S: DataMut
+    {
+        let n = self.len();
+        let partition_value = (&self[partition_index]).clone();
+        self.swap(partition_index, n-1);
+        let mut partition_boundary_index = 0;
+        for j in 0..n-1 {
+            if self[j] < partition_value {
+                self.swap(partition_boundary_index, j);
+                partition_boundary_index += 1;
+            }
+        }
+        self.swap(partition_boundary_index, n-1);
+        partition_boundary_index
+    }
 }
 
+fn random_pivot(n: usize) -> usize
+{
+    let mut rng = thread_rng();
+    rng.gen_range(0, n)
+}
diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index d7e6fcd6f..f28cd06d4 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -13,9 +13,6 @@ use itertools::free::enumerate;
 use imp_prelude::*;
 use numeric_util;
 
-use rand::prelude::*;
-use rand::thread_rng;
-
 use {
     LinalgScalar,
     FoldWhile,
@@ -153,7 +150,7 @@ impl<A, S, D> ArrayBase<S, D>
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = ((n as f32) * q).ceil() as usize;
-        let mapping = |x| ith_mut(x, if i == 0 {0} else {i-1});
+        let mapping = |mut x: ArrayViewMut1<A>| x.ith_mut(if i == 0 {0} else {i-1});
         let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
         azip!(mut lane (self.lanes_mut(axis)), mut out in {
             *out = mapping(lane);
@@ -246,74 +243,3 @@ impl<A, S, D> ArrayBase<S, D>
             }).is_done()
     }
 }
-
-/// Return the `i`-th element of `a` if `a` were to be a 1-dimensional
-/// array sorted in increasing order.
-///
-/// `a` is shuffled **in place** to retrieve the desired element:
-/// no copy of the array is allocated.
-/// No assumptions should be made on the ordering of `a` elements
-/// after this computation.
-///
-/// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
-/// - average case: O(n);
-/// - worst case: O(n^2);
-/// where n is the number of elements in `a`.
-///
-/// **Panics** if `i` is greater than or equal to n.
-fn ith_mut<A>(mut a: ArrayViewMut1<A>, i: usize) -> A
-    where A: Ord + Clone
-{
-    let n = a.len();
-    if n == 1 {
-        (&a[0]).clone()
-    } else {
-        let pivot_index = random_pivot(n);
-        let partition_index = partition_mut(&mut a.view_mut(), pivot_index);
-        if i == partition_index {
-            (&a[partition_index]).clone()
-        } else if i < partition_index {
-            ith_mut(a.slice_mut(s![0..partition_index]), i)
-        } else {
-            ith_mut(a.slice_mut(s![(partition_index+1)..n]), i - partition_index - 1)
-        }
-    }
-}
-
-fn random_pivot(n: usize) -> usize
-{
-    let mut rng = thread_rng();
-    rng.gen_range(0, n)
-}
-
-/// Return the index of `a[partition_index`]` if `a` were to be sorted
-/// in increasing order.
-/// `a` elements are rearranged in such a way that `a[partition_index]`
-/// is in the position it would be in an array sorted in increasing order.
-/// All elements smaller than `a[partition_index]` are moved to its
-/// left and all elements equal or greater than `a[partition_index]`
-/// are moved to its right.
-/// The ordering of the elements in the two partitions is undefined.
-///
-/// `a` is shuffled **in place** to operate the desired partition:
-/// no copy of the array is allocated.
-///
-/// Complexity: O(n), where n is the number of elements in `a`.
-///
-/// **Panics** if `partition_index` is greater than or equal to n.
-fn partition_mut<A>(a: &mut ArrayViewMut1<A>, partition_index: usize) -> usize
-    where A: Ord + Clone
-{
-    let n = a.len();
-    let partition_value = (&a[partition_index]).clone();
-    a.swap(partition_index, n-1);
-    let mut partition_boundary_index = 0;
-    for j in 0..n-1 {
-        if a[j] < partition_value {
-            a.swap(partition_boundary_index, j);
-            partition_boundary_index += 1;
-        }
-    }
-    a.swap(partition_boundary_index, n-1);
-    partition_boundary_index
-}

From 08ccac95a7b2bad908ac6e8e861dc2eb87de89d1 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 18:15:42 +0100
Subject: [PATCH 50/93] Updated function calls in tests.

---
 tests/array.rs | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/array.rs b/tests/array.rs
index c01408d12..ecca1c1da 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -1645,14 +1645,14 @@ fn array_macros() {
 fn test_partition_mut() {
     let mut a = arr1(&[1, 3, 2, 10, 10]);
     let n = a.len();
-    let j = partition_mut(&mut a.view_mut(), n-1);
+    let j = a.view_mut().partition_mut(n-1);
     assert_eq!(j, 3);
     for i in 0..j {
         assert!(a[i] <= a[j]);
     }
     let mut a = arr1(&[2, 3, 4, 1]);
     let n = a.len();
-    let j = partition_mut(&mut a.view_mut(), n-1);
+    let j = a.view_mut().partition_mut(n-1);
     assert_eq!(j, 0);
     let n = a.len();
     for i in j+1..n {
@@ -1663,11 +1663,11 @@ fn test_partition_mut() {
 #[test]
 fn test_ith_mut() {
     let a = arr1(&[1, 3, 2, 10]);
-    let j = ith_mut(a.clone().view_mut(), 2);
+    let j = a.clone().view_mut().ith_mut(2);
     assert_eq!(j, 3);
-    let j = ith_mut(a.clone().view_mut(), 1);
+    let j = a.clone().view_mut().ith_mut(1);
     assert_eq!(j, 2);
-    let j = ith_mut(a.clone().view_mut(), 3);
+    let j = a.clone().view_mut().ith_mut(3);
     assert_eq!(j, 10);
 }
 

From f9ed81ec97e804b939fc8b949bb0073a178e6955 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 2 Jun 2018 18:29:07 +0100
Subject: [PATCH 51/93] Improved test formatting

---
 src/impl_1d.rs              | 16 ++++++++--------
 src/numeric/impl_numeric.rs |  8 ++++----
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 5b5d1dad0..1848a9297 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -27,8 +27,8 @@ impl<A, S> ArrayBase<S, Ix1>
         }
     }
 
-    /// Return the `i`-th element of the array if it were to be a 1-dimensional
-    /// array sorted in increasing order.
+    /// Return the element that would occupy the `i`-th position if the array
+    /// were sorted in increasing order.
     ///
     /// The array is shuffled **in place** to retrieve the desired element:
     /// no copy of the array is allocated.
@@ -36,11 +36,11 @@ impl<A, S> ArrayBase<S, Ix1>
     /// after this computation.
     ///
     /// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
-    /// - average case: O(n);
-    /// - worst case: O(n^2);
+    /// - average case: O(`n`);
+    /// - worst case: O(`n`^2);
     /// where n is the number of elements in the array.
     ///
-    /// **Panics** if `i` is greater than or equal to n.
+    /// **Panics** if `i` is greater than or equal to `n`.
     pub fn ith_mut(&mut self, i: usize) -> A
         where A: Ord + Clone,
               S: DataMut,
@@ -61,7 +61,7 @@ impl<A, S> ArrayBase<S, Ix1>
         }
     }
 
-    /// Return the index of `self[partition_index`]` if `self` were to be sorted
+    /// Return the index of `self[partition_index]` if `self` were to be sorted
     /// in increasing order.
     /// `self` elements are rearranged in such a way that `self[partition_index]`
     /// is in the position it would be in an array sorted in increasing order.
@@ -73,9 +73,9 @@ impl<A, S> ArrayBase<S, Ix1>
     /// `self` is shuffled **in place** to operate the desired partition:
     /// no copy of the array is allocated.
     ///
-    /// Complexity: O(n), where n is the number of elements in the array.
+    /// Complexity: O(`n`), where `n` is the number of elements in the array.
     ///
-    /// **Panics** if `partition_index` is greater than or equal to n.
+    /// **Panics** if `partition_index` is greater than or equal to `n`.
     pub fn partition_mut(&mut self, partition_index: usize) -> usize
         where A: Ord + Clone,
               S: DataMut
diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index f28cd06d4..659d93185 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -117,7 +117,7 @@ impl<A, S, D> ArrayBase<S, D>
 
     /// Return the qth percentile of the data along the specified axis.
     /// `q` needs to be a float between 0 and 1, bounds included.
-    /// The qth percentile for a 1-dimensional lane of length N is defined
+    /// The qth percentile for a 1-dimensional lane of length `N` is defined
     /// as the element that would be indexed as `Nq` if the lane were to be sorted
     /// in increasing order. If `Nq` is not an integer the desired percentile
     /// lies between two data points: we always return the smaller data point.
@@ -136,9 +136,9 @@ impl<A, S, D> ArrayBase<S, D>
     /// after this computation.
     ///
     /// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
-    /// - average case: O(m);
-    /// - worst case: O(m^2);
-    /// where m is the number of elements in the array.
+    /// - average case: O(`m`);
+    /// - worst case: O(`m`^2);
+    /// where `m` is the number of elements in the array.
     ///
     /// **Panics** if `axis` is out of bounds, if `q` is strictly smaller
     /// than 0 or strictly bigger than 1.

From fd11be087f6dc777b1bc4438b653dcd147d53e3b Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 5 Jun 2018 22:16:39 +0100
Subject: [PATCH 52/93] Fixed doc typo.

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 659d93185..c8b3bc5a9 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -130,7 +130,7 @@ impl<A, S, D> ArrayBase<S, D>
     ///
     /// The array is shuffled **in place** along each 1-dimensional lane in
     /// order to produce the required percentile without allocating a copy
-    /// of the original array. Each 1-dimensional lane is shuffled indipendently
+    /// of the original array. Each 1-dimensional lane is shuffled independently
     /// from the others.
     /// No assumptions should be made on the ordering of the array elements
     /// after this computation.

From e4325a12a4a7a0f6370c086e51f69d41324f65c3 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 5 Jun 2018 22:18:58 +0100
Subject: [PATCH 53/93] Removed unnecessary referencing for cloning.

---
 src/impl_1d.rs | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 1848a9297..a4f3a2e4f 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -47,12 +47,12 @@ impl<A, S> ArrayBase<S, Ix1>
     {
         let n = self.len();
         if n == 1 {
-            (&self[0]).clone()
+            self[0].clone()
         } else {
             let pivot_index = random_pivot(n);
             let partition_index = self.view_mut().partition_mut(pivot_index);
             if i == partition_index {
-                (&self[partition_index]).clone()
+                self[partition_index].clone()
             } else if i < partition_index {
                 self.slice_mut(s![0..partition_index]).ith_mut(i)
             } else {
@@ -81,7 +81,7 @@ impl<A, S> ArrayBase<S, Ix1>
               S: DataMut
     {
         let n = self.len();
-        let partition_value = (&self[partition_index]).clone();
+        let partition_value = self[partition_index].clone();
         self.swap(partition_index, n-1);
         let mut partition_boundary_index = 0;
         for j in 0..n-1 {

From 35dd35bb5de7e2c022456423cd6477c8d78fa9b6 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 5 Jun 2018 22:22:09 +0100
Subject: [PATCH 54/93] Removed unnecessary call to view_mut.

---
 src/impl_1d.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index a4f3a2e4f..b3dc88918 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -50,7 +50,7 @@ impl<A, S> ArrayBase<S, Ix1>
             self[0].clone()
         } else {
             let pivot_index = random_pivot(n);
-            let partition_index = self.view_mut().partition_mut(pivot_index);
+            let partition_index = self.partition_mut(pivot_index);
             if i == partition_index {
                 self[partition_index].clone()
             } else if i < partition_index {

From 74756106b05e0d043861de2dc502047d2f593e3d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 11 Jun 2018 21:09:30 +0100
Subject: [PATCH 55/93] Added Hoare routine, with docs and proper reworking of
 ith_mut

---
 src/impl_1d.rs | 53 +++++++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 48 insertions(+), 5 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index b3dc88918..3af638b8f 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -50,11 +50,9 @@ impl<A, S> ArrayBase<S, Ix1>
             self[0].clone()
         } else {
             let pivot_index = random_pivot(n);
-            let partition_index = self.partition_mut(pivot_index);
-            if i == partition_index {
-                self[partition_index].clone()
-            } else if i < partition_index {
-                self.slice_mut(s![0..partition_index]).ith_mut(i)
+            let partition_index = self.hoare_partition_mut(pivot_index);
+            if i <= partition_index {
+                self.slice_mut(s![0..partition_index+1]).ith_mut(i)
             } else {
                 self.slice_mut(s![(partition_index+1)..n]).ith_mut(i - partition_index - 1)
             }
@@ -73,7 +71,10 @@ impl<A, S> ArrayBase<S, Ix1>
     /// `self` is shuffled **in place** to operate the desired partition:
     /// no copy of the array is allocated.
     ///
+    /// The method uses Lomuto's partition algorithm.
     /// Complexity: O(`n`), where `n` is the number of elements in the array.
+    /// Average number of element swaps: (n - 1)/2 (see
+    /// (link)[https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550])
     ///
     /// **Panics** if `partition_index` is greater than or equal to `n`.
     pub fn partition_mut(&mut self, partition_index: usize) -> usize
@@ -93,6 +94,48 @@ impl<A, S> ArrayBase<S, Ix1>
         self.swap(partition_boundary_index, n-1);
         partition_boundary_index
     }
+
+    /// Return a `partition_index`.
+    /// `self` elements are rearranged in such a way that all elements
+    /// in `self[..partition_index]` are smaller than or equal to all elements
+    /// in `self[(partition_index+1)..]`.
+    /// The ordering of the elements in the two partitions is undefined.
+    ///
+    /// `self` is shuffled **in place** to operate the desired partition:
+    /// no copy of the array is allocated.
+    ///
+    /// The method uses Hoare's partition algorithm.
+    /// Complexity: O(`n`), where `n` is the number of elements in the array.
+    /// Average number of element swaps: n/6 - 1/3 (see
+    /// (link)[https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550])
+    ///
+    /// **Panics** if `partition_index` is greater than or equal to `n`.
+    fn hoare_partition_mut(&mut self, partition_index: usize) -> usize
+        where A: Ord + Clone,
+              S: DataMut
+    {
+        let partition_value = self[partition_index].clone();
+        self.swap(partition_index, 0);
+
+        let mut i: isize = -1;
+        let mut j: isize = self.len() as isize;
+        loop {
+            loop {
+                i += 1;
+                if self[i as usize] >= partition_value { break }
+            }
+            loop {
+                j -= 1;
+                if self[j as usize] <= partition_value { break }
+            }
+            if i >= j {
+                break
+            } else {
+                self.swap(i as usize, j as usize);
+            }
+        }
+        j as usize
+    }
 }
 
 fn random_pivot(n: usize) -> usize

From 5240811d8318d33ed7adbf1b2ac1ff844e3eef5d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Mon, 11 Jun 2018 21:19:48 +0100
Subject: [PATCH 56/93] Added some tests for hoare_partition_mut

---
 src/impl_1d.rs | 21 +++++++++++++++++++++
 1 file changed, 21 insertions(+)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 3af638b8f..1f249f557 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -143,3 +143,24 @@ fn random_pivot(n: usize) -> usize
     let mut rng = thread_rng();
     rng.gen_range(0, n)
 }
+
+#[test]
+fn test_hoare_partition_mut() {
+    let mut l = vec!(
+        arr1(&[1, 3, 2, 10, 10]),
+        arr1(&[355, 453, 452, 391, 289, 343,  44, 154, 271,  44, 314, 276, 160,
+               469, 191, 138, 163, 308, 395,   3, 416, 391, 210, 354, 200]),
+        arr1(&[ 84, 192, 216, 159,  89, 296,  35, 213, 456, 278,  98,  52, 308,
+               418, 329, 173, 286, 106, 366, 129, 125, 450,  23, 463, 151]),
+    );
+
+    for a in l.iter_mut() {
+        let n = a.len();
+        let partition_index = a.partition_mut(n-1);
+        for i in 0..partition_index+1 {
+            for j in (partition_index+1)..n {
+                assert!(a[i] <= a[j]);
+            }
+        }
+    }
+}

From ca2f3f9b25a06abe0f3c278e08cd4ce2e6f033dd Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Tue, 12 Jun 2018 06:40:59 +0100
Subject: [PATCH 57/93] Improved docs and tests

---
 src/impl_1d.rs | 29 +++++++++++++++++------------
 1 file changed, 17 insertions(+), 12 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 1f249f557..ff5840712 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -96,9 +96,11 @@ impl<A, S> ArrayBase<S, Ix1>
     }
 
     /// Return a `partition_index`.
-    /// `self` elements are rearranged in such a way that all elements
-    /// in `self[..partition_index]` are smaller than or equal to all elements
-    /// in `self[(partition_index+1)..]`.
+    /// Let `pivot_value=self[pivot_index]` before the function body is executed.
+    /// Then `self` elements are rearranged in such a way that all elements
+    /// in `self[..partition_index]` are smaller than or equal to `pivot_value`
+    /// while all elements in `self[(partition_index+1)..]` are greater than
+    /// or equal to `pivot_value`.
     /// The ordering of the elements in the two partitions is undefined.
     ///
     /// `self` is shuffled **in place** to operate the desired partition:
@@ -110,23 +112,23 @@ impl<A, S> ArrayBase<S, Ix1>
     /// (link)[https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550])
     ///
     /// **Panics** if `partition_index` is greater than or equal to `n`.
-    fn hoare_partition_mut(&mut self, partition_index: usize) -> usize
+    fn hoare_partition_mut(&mut self, pivot_index: usize) -> usize
         where A: Ord + Clone,
               S: DataMut
     {
-        let partition_value = self[partition_index].clone();
-        self.swap(partition_index, 0);
+        let pivot_value = self[pivot_index].clone();
+        self.swap(pivot_index, 0);
 
         let mut i: isize = -1;
         let mut j: isize = self.len() as isize;
         loop {
             loop {
                 i += 1;
-                if self[i as usize] >= partition_value { break }
+                if self[i as usize] >= pivot_value { break }
             }
             loop {
                 j -= 1;
-                if self[j as usize] <= partition_value { break }
+                if self[j as usize] <= pivot_value { break }
             }
             if i >= j {
                 break
@@ -156,11 +158,14 @@ fn test_hoare_partition_mut() {
 
     for a in l.iter_mut() {
         let n = a.len();
-        let partition_index = a.partition_mut(n-1);
+        let pivot_index = n-1;
+        let pivot_value = a[pivot_index].clone();
+        let partition_index = a.partition_mut(pivot_index);
         for i in 0..partition_index+1 {
-            for j in (partition_index+1)..n {
-                assert!(a[i] <= a[j]);
-            }
+            assert!(a[i] <= pivot_value);
+        }
+        for j in (partition_index+1)..n {
+            assert!(pivot_value <= a[j]);
         }
     }
 }

From 39990316fdcea6cb9f6c425900950b6a7b2e82cf Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Fri, 15 Jun 2018 07:45:49 +0100
Subject: [PATCH 58/93] Using map_axis_mut to simplify percentile_axis_mut

---
 src/numeric/impl_numeric.rs | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index d2dadcc93..84de7500e 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -151,11 +151,7 @@ impl<A, S, D> ArrayBase<S, D>
         let n = self.len_of(axis);
         let i = ((n as f32) * q).ceil() as usize;
         let mapping = |mut x: ArrayViewMut1<A>| x.ith_mut(if i == 0 {0} else {i-1});
-        let mut out = Array::zeros(self.view().remove_axis(axis).raw_dim());
-        azip!(mut lane (self.lanes_mut(axis)), mut out in {
-            *out = mapping(lane);
-        });
-        out
+        self.map_axis_mut(axis, mapping)
     }
 
     /// Return variance along `axis`.

From ee4eb05c577e4db4ebfe804e3654d496342e6493 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 18:47:47 +0100
Subject: [PATCH 59/93] Remove patch number specification for rand.

---
 Cargo.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Cargo.toml b/Cargo.toml
index dacf58dd4..6ff6af170 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -28,7 +28,7 @@ num-traits = "0.2"
 num-complex = "0.2"
 rustc-serialize = { version = "0.3.20", optional = true }
 itertools = { version = "0.7.0", default-features = false }
-rand = { version = "0.5.0" }
+rand = { version = "0.5" }
 
 # Use via the `blas` crate feature!
 cblas-sys = { version = "0.1.4", optional = true, default-features = false }

From 4be21462b1d141e14fe048e21439dea9aba053c8 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 18:55:55 +0100
Subject: [PATCH 60/93] Refactored ith_mut. Removed random_pivot.

---
 src/impl_1d.rs | 17 +++++++----------
 1 file changed, 7 insertions(+), 10 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index ff5840712..eebac993a 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -49,12 +49,15 @@ impl<A, S> ArrayBase<S, Ix1>
         if n == 1 {
             self[0].clone()
         } else {
-            let pivot_index = random_pivot(n);
+            let mut rng = thread_rng();
+            let pivot_index = rng.gen_range(0, n);
             let partition_index = self.hoare_partition_mut(pivot_index);
-            if i <= partition_index {
-                self.slice_mut(s![0..partition_index+1]).ith_mut(i)
+            if i < partition_index {
+                self.slice_mut(s![..partition_index]).ith_mut(i)
+            } else if i == partition_index {
+                self[i].clone()
             } else {
-                self.slice_mut(s![(partition_index+1)..n]).ith_mut(i - partition_index - 1)
+                self.slice_mut(s![partition_index+1..]).ith_mut(i - (partition_index+1))
             }
         }
     }
@@ -140,12 +143,6 @@ impl<A, S> ArrayBase<S, Ix1>
     }
 }
 
-fn random_pivot(n: usize) -> usize
-{
-    let mut rng = thread_rng();
-    rng.gen_range(0, n)
-}
-
 #[test]
 fn test_hoare_partition_mut() {
     let mut l = vec!(

From 37e90a94d0b6303ff72e8b4ae54225e5d26ba720 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 18:58:26 +0100
Subject: [PATCH 61/93] Refactored partition_hoare_mut.

---
 src/impl_1d.rs | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index eebac993a..80c21ec36 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -122,24 +122,28 @@ impl<A, S> ArrayBase<S, Ix1>
         let pivot_value = self[pivot_index].clone();
         self.swap(pivot_index, 0);
 
-        let mut i: isize = -1;
-        let mut j: isize = self.len() as isize;
+        let n = self.len();
+        let mut i = 1;
+        let mut j = n - 1;
         loop {
             loop {
+                if i > j { break }
+                if self[i] >= pivot_value { break }
                 i += 1;
-                if self[i as usize] >= pivot_value { break }
             }
-            loop {
+            while pivot_value <= self[j] {
                 j -= 1;
-                if self[j as usize] <= pivot_value { break }
             }
             if i >= j {
                 break
             } else {
-                self.swap(i as usize, j as usize);
+                self.swap(i, j);
+                i += 1;
+                j -= 1;
             }
         }
-        j as usize
+        self.swap(0, i-1);
+        i-1
     }
 }
 

From 5145b8c1ac87ea859a6a2882e80cad8745065356 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:00:16 +0100
Subject: [PATCH 62/93] Renamed hoare_partition_mut to partition_mut.

---
 src/impl_1d.rs | 45 +++++----------------------------------------
 1 file changed, 5 insertions(+), 40 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 80c21ec36..93af0cbbc 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -51,7 +51,7 @@ impl<A, S> ArrayBase<S, Ix1>
         } else {
             let mut rng = thread_rng();
             let pivot_index = rng.gen_range(0, n);
-            let partition_index = self.hoare_partition_mut(pivot_index);
+            let partition_index = self.partition_mut(pivot_index);
             if i < partition_index {
                 self.slice_mut(s![..partition_index]).ith_mut(i)
             } else if i == partition_index {
@@ -62,46 +62,11 @@ impl<A, S> ArrayBase<S, Ix1>
         }
     }
 
-    /// Return the index of `self[partition_index]` if `self` were to be sorted
-    /// in increasing order.
-    /// `self` elements are rearranged in such a way that `self[partition_index]`
-    /// is in the position it would be in an array sorted in increasing order.
-    /// All elements smaller than `self[partition_index]` are moved to its
-    /// left and all elements equal or greater than `self[partition_index]`
-    /// are moved to its right.
-    /// The ordering of the elements in the two partitions is undefined.
-    ///
-    /// `self` is shuffled **in place** to operate the desired partition:
-    /// no copy of the array is allocated.
-    ///
-    /// The method uses Lomuto's partition algorithm.
-    /// Complexity: O(`n`), where `n` is the number of elements in the array.
-    /// Average number of element swaps: (n - 1)/2 (see
-    /// (link)[https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550])
-    ///
-    /// **Panics** if `partition_index` is greater than or equal to `n`.
-    pub fn partition_mut(&mut self, partition_index: usize) -> usize
-        where A: Ord + Clone,
-              S: DataMut
-    {
-        let n = self.len();
-        let partition_value = self[partition_index].clone();
-        self.swap(partition_index, n-1);
-        let mut partition_boundary_index = 0;
-        for j in 0..n-1 {
-            if self[j] < partition_value {
-                self.swap(partition_boundary_index, j);
-                partition_boundary_index += 1;
-            }
-        }
-        self.swap(partition_boundary_index, n-1);
-        partition_boundary_index
-    }
-
     /// Return a `partition_index`.
+    ///
     /// Let `pivot_value=self[pivot_index]` before the function body is executed.
     /// Then `self` elements are rearranged in such a way that all elements
-    /// in `self[..partition_index]` are smaller than or equal to `pivot_value`
+    /// in `self[..partition_index]` are smaller than `pivot_value`
     /// while all elements in `self[(partition_index+1)..]` are greater than
     /// or equal to `pivot_value`.
     /// The ordering of the elements in the two partitions is undefined.
@@ -115,7 +80,7 @@ impl<A, S> ArrayBase<S, Ix1>
     /// (link)[https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550])
     ///
     /// **Panics** if `partition_index` is greater than or equal to `n`.
-    fn hoare_partition_mut(&mut self, pivot_index: usize) -> usize
+    fn partition_mut(&mut self, pivot_index: usize) -> usize
         where A: Ord + Clone,
               S: DataMut
     {
@@ -148,7 +113,7 @@ impl<A, S> ArrayBase<S, Ix1>
 }
 
 #[test]
-fn test_hoare_partition_mut() {
+fn test_partition_mut() {
     let mut l = vec!(
         arr1(&[1, 3, 2, 10, 10]),
         arr1(&[355, 453, 452, 391, 289, 343,  44, 154, 271,  44, 314, 276, 160,

From 630106b2013a32fc6391571395d88c502fd23839 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:01:49 +0100
Subject: [PATCH 63/93] Reworked docs.

---
 src/impl_1d.rs | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 93af0cbbc..76d04c7e2 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -62,13 +62,14 @@ impl<A, S> ArrayBase<S, Ix1>
         }
     }
 
-    /// Return a `partition_index`.
+    /// Return the index of `self[partition_index]` if `self` were to be sorted
+    /// in increasing order.
     ///
-    /// Let `pivot_value=self[pivot_index]` before the function body is executed.
-    /// Then `self` elements are rearranged in such a way that all elements
-    /// in `self[..partition_index]` are smaller than `pivot_value`
-    /// while all elements in `self[(partition_index+1)..]` are greater than
-    /// or equal to `pivot_value`.
+    /// `self` elements are rearranged in such a way that `self[partition_index]`
+    /// is in the position it would be in an array sorted in increasing order.
+    /// All elements smaller than `self[partition_index]` are moved to its
+    /// left and all elements equal or greater than `self[partition_index]`
+    /// are moved to its right.
     /// The ordering of the elements in the two partitions is undefined.
     ///
     /// `self` is shuffled **in place** to operate the desired partition:

From d3cc1aca66d98d628b71169989f1f973f23f0022 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:05:01 +0100
Subject: [PATCH 64/93] Modified signature. Simplified index calculation

---
 src/numeric/impl_numeric.rs | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 84de7500e..2e23a6651 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -142,15 +142,15 @@ impl<A, S, D> ArrayBase<S, D>
     ///
     /// **Panics** if `axis` is out of bounds, if `q` is strictly smaller
     /// than 0 or strictly bigger than 1.
-    pub fn percentile_axis_mut(&mut self, axis: Axis, q: f32) -> Array<A, D::Smaller>
+    pub fn percentile_axis_mut(&mut self, axis: Axis, q: f64) -> Array<A, D::Smaller>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,
               S: DataMut,
     {
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
-        let i = ((n as f32) * q).ceil() as usize;
-        let mapping = |mut x: ArrayViewMut1<A>| x.ith_mut(if i == 0 {0} else {i-1});
+        let i = (((n - 1) as f64) * q).floor() as usize;
+        let mapping = |mut x: ArrayViewMut1<A>| x.ith_mut(i);
         self.map_axis_mut(axis, mapping)
     }
 

From a6fecc7ce1a4299a73c644371fcf4e1106d96282 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:05:53 +0100
Subject: [PATCH 65/93] Inlined closure.

---
 src/numeric/impl_numeric.rs | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 2e23a6651..fea7a82b9 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -150,8 +150,7 @@ impl<A, S, D> ArrayBase<S, D>
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = (((n - 1) as f64) * q).floor() as usize;
-        let mapping = |mut x: ArrayViewMut1<A>| x.ith_mut(i);
-        self.map_axis_mut(axis, mapping)
+        self.map_axis_mut(axis, |mut x| x.ith_mut(i))
     }
 
     /// Return variance along `axis`.

From 90774bc4838ab9e2d5619d5349c7d43aab139246 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:07:19 +0100
Subject: [PATCH 66/93] Modified docs to include NaN cases.

---
 src/numeric/impl_numeric.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index fea7a82b9..add5f484f 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -140,8 +140,8 @@ impl<A, S, D> ArrayBase<S, D>
     /// - worst case: O(`m`^2);
     /// where `m` is the number of elements in the array.
     ///
-    /// **Panics** if `axis` is out of bounds, if `q` is strictly smaller
-    /// than 0 or strictly bigger than 1.
+    /// **Panics** if `axis` is out of bounds or if `q` is not between
+    /// `0.` and `1.` (inclusive).
     pub fn percentile_axis_mut(&mut self, axis: Axis, q: f64) -> Array<A, D::Smaller>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,

From 3a7f31427ea7c0cdce0471e2cb9cbe777400c12e Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:10:27 +0100
Subject: [PATCH 67/93] Fixing visibility typo.

---
 src/impl_1d.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 76d04c7e2..e2811952e 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -81,7 +81,7 @@ impl<A, S> ArrayBase<S, Ix1>
     /// (link)[https://cs.stackexchange.com/questions/11458/quicksort-partitioning-hoare-vs-lomuto/11550])
     ///
     /// **Panics** if `partition_index` is greater than or equal to `n`.
-    fn partition_mut(&mut self, pivot_index: usize) -> usize
+    pub fn partition_mut(&mut self, pivot_index: usize) -> usize
         where A: Ord + Clone,
               S: DataMut
     {

From 22415cddc6cbf481d124783892c2db8e0f3ba9be Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:14:00 +0100
Subject: [PATCH 68/93] Fixed edge case

---
 src/impl_1d.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index e2811952e..9e2ca6ca6 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -98,6 +98,7 @@ impl<A, S> ArrayBase<S, Ix1>
                 i += 1;
             }
             while pivot_value <= self[j] {
+                if j == 1 { break }
                 j -= 1;
             }
             if i >= j {

From b1ce4cceda421703d0f33f7610ea3ed9d4acff23 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:17:39 +0100
Subject: [PATCH 69/93] Merged tests. Changed test assertions

---
 src/impl_1d.rs | 24 ------------------------
 tests/array.rs | 35 +++++++++++++++++++++--------------
 2 files changed, 21 insertions(+), 38 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 9e2ca6ca6..6daf9d304 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -113,27 +113,3 @@ impl<A, S> ArrayBase<S, Ix1>
         i-1
     }
 }
-
-#[test]
-fn test_partition_mut() {
-    let mut l = vec!(
-        arr1(&[1, 3, 2, 10, 10]),
-        arr1(&[355, 453, 452, 391, 289, 343,  44, 154, 271,  44, 314, 276, 160,
-               469, 191, 138, 163, 308, 395,   3, 416, 391, 210, 354, 200]),
-        arr1(&[ 84, 192, 216, 159,  89, 296,  35, 213, 456, 278,  98,  52, 308,
-               418, 329, 173, 286, 106, 366, 129, 125, 450,  23, 463, 151]),
-    );
-
-    for a in l.iter_mut() {
-        let n = a.len();
-        let pivot_index = n-1;
-        let pivot_value = a[pivot_index].clone();
-        let partition_index = a.partition_mut(pivot_index);
-        for i in 0..partition_index+1 {
-            assert!(a[i] <= pivot_value);
-        }
-        for j in (partition_index+1)..n {
-            assert!(pivot_value <= a[j]);
-        }
-    }
-}
diff --git a/tests/array.rs b/tests/array.rs
index 253acf10e..0ff1d6f89 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -1705,20 +1705,27 @@ fn array_macros() {
 
 #[test]
 fn test_partition_mut() {
-    let mut a = arr1(&[1, 3, 2, 10, 10]);
-    let n = a.len();
-    let j = a.view_mut().partition_mut(n-1);
-    assert_eq!(j, 3);
-    for i in 0..j {
-        assert!(a[i] <= a[j]);
-    }
-    let mut a = arr1(&[2, 3, 4, 1]);
-    let n = a.len();
-    let j = a.view_mut().partition_mut(n-1);
-    assert_eq!(j, 0);
-    let n = a.len();
-    for i in j+1..n {
-        assert!(a[i] > a[j]);
+    let mut l = vec!(
+        arr1(&[1, 3, 2, 10, 10]),
+        arr1(&[2, 3, 4, 1]),
+        arr1(&[355, 453, 452, 391, 289, 343,  44, 154, 271,  44, 314, 276, 160,
+               469, 191, 138, 163, 308, 395,   3, 416, 391, 210, 354, 200]),
+        arr1(&[ 84, 192, 216, 159,  89, 296,  35, 213, 456, 278,  98,  52, 308,
+               418, 329, 173, 286, 106, 366, 129, 125, 450,  23, 463, 151]),
+    );
+
+    for a in l.iter_mut() {
+        let n = a.len();
+        let pivot_index = n-1;
+        let pivot_value = a[pivot_index].clone();
+        let partition_index = a.partition_mut(pivot_index);
+        for i in 0..partition_index {
+            assert!(a[i] < pivot_value);
+        }
+        assert!(a[partition_index] == pivot_value);
+        for j in (partition_index+1)..n {
+            assert!(pivot_value <= a[j]);
+        }
     }
 }
 

From fcc9f82d164ca9ff6105b69bef4247e646d55af7 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:18:28 +0100
Subject: [PATCH 70/93] Added test for edge case

---
 tests/array.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/array.rs b/tests/array.rs
index 0ff1d6f89..4015bd5b1 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -1706,6 +1706,7 @@ fn array_macros() {
 #[test]
 fn test_partition_mut() {
     let mut l = vec!(
+        arr1(&[1, 1, 1, 1, 1]),
         arr1(&[1, 3, 2, 10, 10]),
         arr1(&[2, 3, 4, 1]),
         arr1(&[355, 453, 452, 391, 289, 343,  44, 154, 271,  44, 314, 276, 160,

From 7e1b5034676cf3c50f058614221d7925989aa445 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:23:08 +0100
Subject: [PATCH 71/93] Renamed ith_mut to sorted_get_mut

---
 src/impl_1d.rs              | 6 +++---
 src/numeric/impl_numeric.rs | 2 +-
 tests/array.rs              | 8 ++++----
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index 6daf9d304..b9bdb61f2 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -41,7 +41,7 @@ impl<A, S> ArrayBase<S, Ix1>
     /// where n is the number of elements in the array.
     ///
     /// **Panics** if `i` is greater than or equal to `n`.
-    pub fn ith_mut(&mut self, i: usize) -> A
+    pub fn sorted_get_mut(&mut self, i: usize) -> A
         where A: Ord + Clone,
               S: DataMut,
     {
@@ -53,11 +53,11 @@ impl<A, S> ArrayBase<S, Ix1>
             let pivot_index = rng.gen_range(0, n);
             let partition_index = self.partition_mut(pivot_index);
             if i < partition_index {
-                self.slice_mut(s![..partition_index]).ith_mut(i)
+                self.slice_mut(s![..partition_index]).sorted_get_mut(i)
             } else if i == partition_index {
                 self[i].clone()
             } else {
-                self.slice_mut(s![partition_index+1..]).ith_mut(i - (partition_index+1))
+                self.slice_mut(s![partition_index+1..]).sorted_get_mut(i - (partition_index+1))
             }
         }
     }
diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index add5f484f..579ac9c31 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -150,7 +150,7 @@ impl<A, S, D> ArrayBase<S, D>
         assert!((0. <= q) && (q <= 1.));
         let n = self.len_of(axis);
         let i = (((n - 1) as f64) * q).floor() as usize;
-        self.map_axis_mut(axis, |mut x| x.ith_mut(i))
+        self.map_axis_mut(axis, |mut x| x.sorted_get_mut(i))
     }
 
     /// Return variance along `axis`.
diff --git a/tests/array.rs b/tests/array.rs
index 4015bd5b1..3d45e1a99 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -1731,13 +1731,13 @@ fn test_partition_mut() {
 }
 
 #[test]
-fn test_ith_mut() {
+fn test_sorted_get_mut() {
     let a = arr1(&[1, 3, 2, 10]);
-    let j = a.clone().view_mut().ith_mut(2);
+    let j = a.clone().view_mut().sorted_get_mut(2);
     assert_eq!(j, 3);
-    let j = a.clone().view_mut().ith_mut(1);
+    let j = a.clone().view_mut().sorted_get_mut(1);
     assert_eq!(j, 2);
-    let j = a.clone().view_mut().ith_mut(3);
+    let j = a.clone().view_mut().sorted_get_mut(3);
     assert_eq!(j, 10);
 }
 

From 393b6bb076cffce476adffa8fa2c310247b40414 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:33:31 +0100
Subject: [PATCH 72/93] Added different interpolation strategies for the
 percentile index.

---
 src/numeric/impl_numeric.rs | 18 ++++++++++++++----
 1 file changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 579ac9c31..6fe06dfe4 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -19,6 +19,12 @@ use {
     Zip,
 };
 
+pub enum InterpolationStrategy {
+    Lower,
+    Nearest,
+    Highest,
+}
+
 /// Numerical methods for arrays.
 impl<A, S, D> ArrayBase<S, D>
     where S: Data<Elem=A>,
@@ -142,15 +148,19 @@ impl<A, S, D> ArrayBase<S, D>
     ///
     /// **Panics** if `axis` is out of bounds or if `q` is not between
     /// `0.` and `1.` (inclusive).
-    pub fn percentile_axis_mut(&mut self, axis: Axis, q: f64) -> Array<A, D::Smaller>
+    pub fn percentile_axis_mut(&mut self, axis: Axis, q: f64, interpolation_strategy: InterpolationStrategy) -> Array<A, D::Smaller>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,
               S: DataMut,
     {
         assert!((0. <= q) && (q <= 1.));
-        let n = self.len_of(axis);
-        let i = (((n - 1) as f64) * q).floor() as usize;
-        self.map_axis_mut(axis, |mut x| x.sorted_get_mut(i))
+        let float_percentile_index = ((self.len_of(axis) - 1) as f64) * q;
+        let percentile_index = match interpolation_strategy {
+            InterpolationStrategy::Lower => float_percentile_index.floor() as usize,
+            InterpolationStrategy::Nearest => float_percentile_index.round() as usize,
+            InterpolationStrategy::Highest => float_percentile_index.ceil() as usize,
+        };
+        self.map_axis_mut(axis, |mut x| x.sorted_get_mut(percentile_index))
     }
 
     /// Return variance along `axis`.

From a9f148d573d9f16936932cd98c7fcda7a04efdea Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 19:41:20 +0100
Subject: [PATCH 73/93] Renamed Highest to Higher, for coherence. Fixed docs,
 added docs to InterpolationStrategy Enum

---
 src/numeric/impl_numeric.rs | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 6fe06dfe4..9664c49ef 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -19,10 +19,13 @@ use {
     Zip,
 };
 
+/// Used to choose the interpolation strategy in [`percentile_axis_mut`].
+///
+/// [`percentile_axis_mut`]: struct.ArrayBase.html#method.percentile_axis_mut
 pub enum InterpolationStrategy {
     Lower,
     Nearest,
-    Highest,
+    Higher,
 }
 
 /// Numerical methods for arrays.
@@ -122,11 +125,14 @@ impl<A, S, D> ArrayBase<S, D>
     }
 
     /// Return the qth percentile of the data along the specified axis.
+    ///
     /// `q` needs to be a float between 0 and 1, bounds included.
     /// The qth percentile for a 1-dimensional lane of length `N` is defined
-    /// as the element that would be indexed as `Nq` if the lane were to be sorted
-    /// in increasing order. If `Nq` is not an integer the desired percentile
-    /// lies between two data points: we always return the smaller data point.
+    /// as the element that would be indexed as `(N-1)q` if the lane were to be sorted
+    /// in increasing order.
+    /// If `(N-1)q` is not an integer the desired percentile lies between
+    /// two data points: we return the lower, nearest or higher datapoint depending
+    /// on `interpolation_strategy`.
     ///
     /// Some examples:
     /// - `q=0.` returns the minimum along each 1-dimensional lane;
@@ -158,7 +164,7 @@ impl<A, S, D> ArrayBase<S, D>
         let percentile_index = match interpolation_strategy {
             InterpolationStrategy::Lower => float_percentile_index.floor() as usize,
             InterpolationStrategy::Nearest => float_percentile_index.round() as usize,
-            InterpolationStrategy::Highest => float_percentile_index.ceil() as usize,
+            InterpolationStrategy::Higher => float_percentile_index.ceil() as usize,
         };
         self.map_axis_mut(axis, |mut x| x.sorted_get_mut(percentile_index))
     }

From 0bd4c71bac51c7affbc269d1ffb9f74dbbd081bb Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 21:48:33 +0100
Subject: [PATCH 74/93] Fixed visibility of InterpolationStrategy

---
 src/lib.rs         | 2 +-
 src/numeric/mod.rs | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/lib.rs b/src/lib.rs
index 919ca674f..1fcd603ea 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -896,7 +896,7 @@ impl<A, S, D> ArrayBase<S, D>
 mod impl_1d;
 mod impl_2d;
 
-mod numeric;
+pub mod numeric;
 
 pub mod linalg;
 
diff --git a/src/numeric/mod.rs b/src/numeric/mod.rs
index 60c5e39f6..88ff33649 100644
--- a/src/numeric/mod.rs
+++ b/src/numeric/mod.rs
@@ -1,3 +1,3 @@
 
-mod impl_numeric;
-
+pub mod impl_numeric;
+pub use self::impl_numeric::InterpolationStrategy;

From ded4e949696cbd956fe2c9c466ff75280a52e745 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 24 Jun 2018 21:48:46 +0100
Subject: [PATCH 75/93] Updated tests to add new argument to
 percentile_axis_mut

---
 tests/array.rs | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/tests/array.rs b/tests/array.rs
index 3d45e1a99..507b21d38 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -8,6 +8,7 @@ extern crate itertools;
 
 use ndarray::{Slice, SliceInfo, SliceOrIndex};
 use ndarray::prelude::*;
+use ndarray::numeric::InterpolationStrategy;
 use ndarray::{
     rcarr2,
     arr3,
@@ -1750,7 +1751,7 @@ fn test_percentile_axis_mut_with_odd_axis_length() {
         [3, 5, 6, 12]
         ]
     );
-    let p = a.percentile_axis_mut(Axis(0), 0.5);
+    let p = a.percentile_axis_mut(Axis(0), 0.5, InterpolationStrategy::Lower);
     assert!(p == a.subview(Axis(0), 1));
 }
 
@@ -1764,20 +1765,20 @@ fn test_percentile_axis_mut_with_even_axis_length() {
         [4, 6, 7, 13]
         ]
     );
-    let q = b.percentile_axis_mut(Axis(0), 0.5);
+    let q = b.percentile_axis_mut(Axis(0), 0.5, InterpolationStrategy::Lower);
     assert!(q == b.subview(Axis(0), 1));
 }
 
 #[test]
 fn test_percentile_axis_mut_to_get_minimum() {
     let mut b = arr2(&[[1, 3, 22, 10]]);
-    let q = b.percentile_axis_mut(Axis(1), 0.);
+    let q = b.percentile_axis_mut(Axis(1), 0., InterpolationStrategy::Lower);
     assert!(q == arr1(&[1]));
 }
 
 #[test]
 fn test_percentile_axis_mut_to_get_maximum() {
     let mut b = arr1(&[1, 3, 22, 10]);
-    let q = b.percentile_axis_mut(Axis(0), 1.);
+    let q = b.percentile_axis_mut(Axis(0), 1., InterpolationStrategy::Lower);
     assert!(q == arr0(22));
 }

From 1ec4e42b87dc0438938ee86d75c8c371a83bcd4d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 1 Sep 2018 19:20:26 +0100
Subject: [PATCH 76/93] Implementing interpolation strategies as trait bounds.

---
 src/numeric/impl_numeric.rs | 53 ++++++++++++++++++++++++++++++++++---
 1 file changed, 49 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 9664c49ef..d49bb8cf6 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -22,10 +22,55 @@ use {
 /// Used to choose the interpolation strategy in [`percentile_axis_mut`].
 ///
 /// [`percentile_axis_mut`]: struct.ArrayBase.html#method.percentile_axis_mut
-pub enum InterpolationStrategy {
-    Lower,
-    Nearest,
-    Higher,
+pub trait Interpolate<T> {
+    fn needs_lower(q: f64, len: usize) -> bool;
+    fn needs_upper(q: f64, len: usize) -> bool;
+    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T;
+}
+
+struct Upper;
+struct Lower;
+struct Nearest;
+
+impl<T> Interpolate<T> for Upper {
+    fn needs_lower(q: f64, len: usize) -> bool {
+        false
+    }
+    fn needs_upper(q: f64, len: usize) -> bool {
+        true
+    }
+    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
+       upper.unwrap()
+    }
+}
+
+impl<T> Interpolate<T> for Lower {
+    fn needs_lower(q: f64, len: usize) -> bool {
+        true
+    }
+    fn needs_upper(q: f64, len: usize) -> bool {
+        false
+    }
+    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
+        lower.unwrap()
+    }
+}
+
+impl<T> Interpolate<T> for Nearest {
+    fn needs_lower(q: f64, len: usize) -> bool {
+        let float_percentile_index = ((len - 1) as f64) * q;
+        (float_percentile_index.round() - float_percentile_index) <= 0.
+    }
+    fn needs_upper(q: f64, len: usize) -> bool {
+        !Self::needs_lower(q, len)
+    }
+    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
+        if Self::needs_lower(q, len) {
+            lower.unwrap()
+        } else {
+            upper.unwrap()
+        }
+    }
 }
 
 /// Numerical methods for arrays.

From 63bbb07526ab8775b959fc28c40b7afdb2f463bf Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 1 Sep 2018 19:28:15 +0100
Subject: [PATCH 77/93] Providing some default methods.

---
 src/numeric/impl_numeric.rs | 14 ++++++++++++--
 1 file changed, 12 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index d49bb8cf6..1f05e9112 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -23,6 +23,17 @@ use {
 ///
 /// [`percentile_axis_mut`]: struct.ArrayBase.html#method.percentile_axis_mut
 pub trait Interpolate<T> {
+    fn float_percentile_index(q: f64, len: usize) -> f64 {
+        ((len - 1) as f64) * q
+    }
+
+    fn lower_index(q: f64, len: usize) -> usize {
+        Self::float_percentile_index(q, len).floor() as usize
+    }
+
+    fn upper_index(q: f64, len: usize) -> usize {
+        Self::float_percentile_index(q, len).ceil() as usize
+    }
     fn needs_lower(q: f64, len: usize) -> bool;
     fn needs_upper(q: f64, len: usize) -> bool;
     fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T;
@@ -58,8 +69,7 @@ impl<T> Interpolate<T> for Lower {
 
 impl<T> Interpolate<T> for Nearest {
     fn needs_lower(q: f64, len: usize) -> bool {
-        let float_percentile_index = ((len - 1) as f64) * q;
-        (float_percentile_index.round() - float_percentile_index) <= 0.
+        ((Self::lower_index(q, len) as f64) - Self::float_percentile_index(q, len)) <= 0.
     }
     fn needs_upper(q: f64, len: usize) -> bool {
         !Self::needs_lower(q, len)

From 3d1a196fc7cdd981dafc3c65620af12cdd2d9e53 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sat, 1 Sep 2018 19:50:44 +0100
Subject: [PATCH 78/93] Mysterious type annotation error.

---
 src/numeric/impl_numeric.rs | 31 +++++++++++++++++++++++--------
 src/numeric/mod.rs          |  1 -
 2 files changed, 23 insertions(+), 9 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 1f05e9112..0ed0bf480 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -69,7 +69,8 @@ impl<T> Interpolate<T> for Lower {
 
 impl<T> Interpolate<T> for Nearest {
     fn needs_lower(q: f64, len: usize) -> bool {
-        ((Self::lower_index(q, len) as f64) - Self::float_percentile_index(q, len)) <= 0.
+        let lower = Self::lower_index(q, len);
+        ((lower as f64) - Self::float_percentile_index(q, len)) <= 0.
     }
     fn needs_upper(q: f64, len: usize) -> bool {
         !Self::needs_lower(q, len)
@@ -209,19 +210,33 @@ impl<A, S, D> ArrayBase<S, D>
     ///
     /// **Panics** if `axis` is out of bounds or if `q` is not between
     /// `0.` and `1.` (inclusive).
-    pub fn percentile_axis_mut(&mut self, axis: Axis, q: f64, interpolation_strategy: InterpolationStrategy) -> Array<A, D::Smaller>
+    pub fn percentile_axis_mut<I>(&mut self, axis: Axis, q: f64) -> Array<A, D::Smaller>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,
               S: DataMut,
+              I: Interpolate<Array<A, D::Smaller>>,
     {
         assert!((0. <= q) && (q <= 1.));
-        let float_percentile_index = ((self.len_of(axis) - 1) as f64) * q;
-        let percentile_index = match interpolation_strategy {
-            InterpolationStrategy::Lower => float_percentile_index.floor() as usize,
-            InterpolationStrategy::Nearest => float_percentile_index.round() as usize,
-            InterpolationStrategy::Higher => float_percentile_index.ceil() as usize,
+        let mut lower = None;
+        let mut upper = None;
+        let axis_len = self.len_of(axis);
+        if I::needs_lower(q, axis_len) {
+            lower = Some(
+                self.map_axis_mut(
+                    axis,
+                    |mut x| x.sorted_get_mut(I::lower_index(q, axis_len))
+                )
+            );
         };
-        self.map_axis_mut(axis, |mut x| x.sorted_get_mut(percentile_index))
+        if I::needs_upper(q, axis_len) {
+            upper = Some(
+                self.map_axis_mut(
+                    axis,
+                    |mut x| x.sorted_get_mut(I::upper_index(q, axis_len))
+                )
+            );
+        };
+        I::interpolate(lower, upper, q, axis_len)
     }
 
     /// Return variance along `axis`.
diff --git a/src/numeric/mod.rs b/src/numeric/mod.rs
index 88ff33649..443d209ef 100644
--- a/src/numeric/mod.rs
+++ b/src/numeric/mod.rs
@@ -1,3 +1,2 @@
 
 pub mod impl_numeric;
-pub use self::impl_numeric::InterpolationStrategy;

From 10a72d0adc7bc3404953f34b1993009243a9f162 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 11:11:41 +0100
Subject: [PATCH 79/93] Fixed.

---
 src/numeric/impl_numeric.rs | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 0ed0bf480..080e06766 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -69,14 +69,14 @@ impl<T> Interpolate<T> for Lower {
 
 impl<T> Interpolate<T> for Nearest {
     fn needs_lower(q: f64, len: usize) -> bool {
-        let lower = Self::lower_index(q, len);
-        ((lower as f64) - Self::float_percentile_index(q, len)) <= 0.
+        let lower = <Self as Interpolate<T>>::lower_index(q, len);
+        ((lower as f64) - <Self as Interpolate<T>>::float_percentile_index(q, len)) <= 0.
     }
     fn needs_upper(q: f64, len: usize) -> bool {
-        !Self::needs_lower(q, len)
+        !<Self as Interpolate<T>>::needs_lower(q, len)
     }
     fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
-        if Self::needs_lower(q, len) {
+        if <Self as Interpolate<T>>::needs_lower(q, len) {
             lower.unwrap()
         } else {
             upper.unwrap()

From f8b7713c3be5f3c4451e2c5f5c7270e6fa195540 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:21:39 +0100
Subject: [PATCH 80/93] Making Lower, Upper, Nearest public.

---
 src/numeric/impl_numeric.rs | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 080e06766..3e4eea6fa 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -39,9 +39,9 @@ pub trait Interpolate<T> {
     fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T;
 }
 
-struct Upper;
-struct Lower;
-struct Nearest;
+pub struct Upper;
+pub struct Lower;
+pub struct Nearest;
 
 impl<T> Interpolate<T> for Upper {
     fn needs_lower(q: f64, len: usize) -> bool {

From e7d1f617a998e82b31d26dcf8b0172cc96f82852 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:21:59 +0100
Subject: [PATCH 81/93] Removed empty line.

---
 src/numeric/mod.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/numeric/mod.rs b/src/numeric/mod.rs
index 443d209ef..bc959a241 100644
--- a/src/numeric/mod.rs
+++ b/src/numeric/mod.rs
@@ -1,2 +1 @@
-
 pub mod impl_numeric;

From cec96cafae978b1e7d60315f97412a628086785d Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:33:27 +0100
Subject: [PATCH 82/93] Re-exporting functions, structs and traits.

---
 src/numeric/mod.rs | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/numeric/mod.rs b/src/numeric/mod.rs
index bc959a241..30024d902 100644
--- a/src/numeric/mod.rs
+++ b/src/numeric/mod.rs
@@ -1 +1,3 @@
-pub mod impl_numeric;
+pub use self::impl_numeric::*;
+
+mod impl_numeric;

From c786eacf4a7cf3bf3b4a00c6e584321f125b4ee8 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:33:37 +0100
Subject: [PATCH 83/93] Fixed import.

---
 tests/array.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/array.rs b/tests/array.rs
index 1369c1be6..51cc36209 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -8,7 +8,7 @@ extern crate itertools;
 
 use ndarray::{Slice, SliceInfo, SliceOrIndex};
 use ndarray::prelude::*;
-use ndarray::numeric::InterpolationStrategy;
+use ndarray::numeric::{Lower, Nearest, Upper};
 use ndarray::{
     rcarr2,
     arr3,

From 744a427748cede98c3556dc8feefe1ff6be7c035 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:39:32 +0100
Subject: [PATCH 84/93] Fixed test cases.

---
 tests/array.rs | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/tests/array.rs b/tests/array.rs
index 51cc36209..75a2610fa 100644
--- a/tests/array.rs
+++ b/tests/array.rs
@@ -8,7 +8,7 @@ extern crate itertools;
 
 use ndarray::{Slice, SliceInfo, SliceOrIndex};
 use ndarray::prelude::*;
-use ndarray::numeric::{Lower, Nearest, Upper};
+use ndarray::numeric::{Lower};
 use ndarray::{
     rcarr2,
     arr3,
@@ -1759,7 +1759,7 @@ fn test_percentile_axis_mut_with_odd_axis_length() {
         [3, 5, 6, 12]
         ]
     );
-    let p = a.percentile_axis_mut(Axis(0), 0.5, InterpolationStrategy::Lower);
+    let p = a.percentile_axis_mut::<Lower>(Axis(0), 0.5);
     assert!(p == a.subview(Axis(0), 1));
 }
 
@@ -1773,20 +1773,20 @@ fn test_percentile_axis_mut_with_even_axis_length() {
         [4, 6, 7, 13]
         ]
     );
-    let q = b.percentile_axis_mut(Axis(0), 0.5, InterpolationStrategy::Lower);
+    let q = b.percentile_axis_mut::<Lower>(Axis(0), 0.5);
     assert!(q == b.subview(Axis(0), 1));
 }
 
 #[test]
 fn test_percentile_axis_mut_to_get_minimum() {
     let mut b = arr2(&[[1, 3, 22, 10]]);
-    let q = b.percentile_axis_mut(Axis(1), 0., InterpolationStrategy::Lower);
+    let q = b.percentile_axis_mut::<Lower>(Axis(1), 0.);
     assert!(q == arr1(&[1]));
 }
 
 #[test]
 fn test_percentile_axis_mut_to_get_maximum() {
     let mut b = arr1(&[1, 3, 22, 10]);
-    let q = b.percentile_axis_mut(Axis(0), 1., InterpolationStrategy::Lower);
+    let q = b.percentile_axis_mut::<Lower>(Axis(0), 1.);
     assert!(q == arr0(22));
 }

From bdc8c3a8715fb19966f3ca37e6ba93ed80ee4bb3 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:47:12 +0100
Subject: [PATCH 85/93] Updated docs.

---
 src/impl_1d.rs | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/src/impl_1d.rs b/src/impl_1d.rs
index b9bdb61f2..e69162259 100644
--- a/src/impl_1d.rs
+++ b/src/impl_1d.rs
@@ -27,13 +27,18 @@ impl<A, S> ArrayBase<S, Ix1>
         }
     }
 
-    /// Return the element that would occupy the `i`-th position if the array
-    /// were sorted in increasing order.
+    /// Return the element that would occupy the `i`-th position if
+    /// the array were sorted in increasing order.
     ///
     /// The array is shuffled **in place** to retrieve the desired element:
     /// no copy of the array is allocated.
-    /// No assumptions should be made on the ordering of elements
-    /// after this computation.
+    /// After the shuffling, all elements with an index smaller than `i`
+    /// are smaller than the desired element, while all elements with
+    /// an index greater or equal than `i` are greater than or equal
+    /// to the desired element.
+    ///
+    /// No other assumptions should be made on the ordering of the
+    /// elements after this computation.
     ///
     /// Complexity ([quickselect](https://en.wikipedia.org/wiki/Quickselect)):
     /// - average case: O(`n`);

From b5264454b0027f9da62bab8e9e61da59f0320710 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:48:35 +0100
Subject: [PATCH 86/93] Fixed docs.

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 3e4eea6fa..32896035d 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -19,7 +19,7 @@ use {
     Zip,
 };
 
-/// Used to choose the interpolation strategy in [`percentile_axis_mut`].
+/// Used to provide an interpolation strategy to [`percentile_axis_mut`].
 ///
 /// [`percentile_axis_mut`]: struct.ArrayBase.html#method.percentile_axis_mut
 pub trait Interpolate<T> {

From 1b39c6133973d7aade0f9aac9c30c0951fe88450 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 12:49:50 +0100
Subject: [PATCH 87/93] Fixed warnings on unused variables.

---
 src/numeric/impl_numeric.rs | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 32896035d..6ad427fcb 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -44,25 +44,25 @@ pub struct Lower;
 pub struct Nearest;
 
 impl<T> Interpolate<T> for Upper {
-    fn needs_lower(q: f64, len: usize) -> bool {
+    fn needs_lower(_q: f64, _len: usize) -> bool {
         false
     }
-    fn needs_upper(q: f64, len: usize) -> bool {
+    fn needs_upper(_q: f64, _len: usize) -> bool {
         true
     }
-    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
+    fn interpolate(_lower: Option<T>, upper: Option<T>, _q: f64, _len: usize) -> T {
        upper.unwrap()
     }
 }
 
 impl<T> Interpolate<T> for Lower {
-    fn needs_lower(q: f64, len: usize) -> bool {
+    fn needs_lower(_q: f64, _len: usize) -> bool {
         true
     }
-    fn needs_upper(q: f64, len: usize) -> bool {
+    fn needs_upper(_q: f64, _len: usize) -> bool {
         false
     }
-    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
+    fn interpolate(lower: Option<T>, _upper: Option<T>, _q: f64, _len: usize) -> T {
         lower.unwrap()
     }
 }

From e90555e7bbc4bfda4d5a92631e2dfdf265ddae4f Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 14:41:15 +0100
Subject: [PATCH 88/93] Implemented Midpoint. Refactored code to change types
 of variables.

---
 src/numeric/impl_numeric.rs | 47 +++++++++++++++++++++++++++++++------
 1 file changed, 40 insertions(+), 7 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 6ad427fcb..fab40e406 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -6,8 +6,8 @@
 // option. This file may not be copied, modified, or distributed
 // except according to those terms.
 
-use std::ops::Add;
-use libnum::{self, Zero, Float};
+use std::ops::{Add, Div};
+use libnum::{self, Zero, Float, FromPrimitive};
 use itertools::free::enumerate;
 
 use imp_prelude::*;
@@ -36,12 +36,17 @@ pub trait Interpolate<T> {
     }
     fn needs_lower(q: f64, len: usize) -> bool;
     fn needs_upper(q: f64, len: usize) -> bool;
-    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T;
+    fn interpolate<D>(lower: Option<Array<T, D>>,
+                      upper: Option<Array<T, D>>,
+                      q: f64,
+                      len: usize) -> Array<T, D>
+        where D: Dimension;
 }
 
 pub struct Upper;
 pub struct Lower;
 pub struct Nearest;
+pub struct Midpoint;
 
 impl<T> Interpolate<T> for Upper {
     fn needs_lower(_q: f64, _len: usize) -> bool {
@@ -50,7 +55,10 @@ impl<T> Interpolate<T> for Upper {
     fn needs_upper(_q: f64, _len: usize) -> bool {
         true
     }
-    fn interpolate(_lower: Option<T>, upper: Option<T>, _q: f64, _len: usize) -> T {
+    fn interpolate<D>(_lower: Option<Array<T, D>>,
+                      upper: Option<Array<T, D>>,
+                      _q: f64,
+                      _len: usize) -> Array<T, D> {
        upper.unwrap()
     }
 }
@@ -62,7 +70,10 @@ impl<T> Interpolate<T> for Lower {
     fn needs_upper(_q: f64, _len: usize) -> bool {
         false
     }
-    fn interpolate(lower: Option<T>, _upper: Option<T>, _q: f64, _len: usize) -> T {
+    fn interpolate<D>(lower: Option<Array<T, D>>,
+                      _upper: Option<Array<T, D>>,
+                      _q: f64,
+                      _len: usize) -> Array<T, D> {
         lower.unwrap()
     }
 }
@@ -75,7 +86,10 @@ impl<T> Interpolate<T> for Nearest {
     fn needs_upper(q: f64, len: usize) -> bool {
         !<Self as Interpolate<T>>::needs_lower(q, len)
     }
-    fn interpolate(lower: Option<T>, upper: Option<T>, q: f64, len: usize) -> T {
+    fn interpolate<D>(lower: Option<Array<T, D>>,
+                      upper: Option<Array<T, D>>,
+                      q: f64,
+                      len: usize) -> Array<T, D> {
         if <Self as Interpolate<T>>::needs_lower(q, len) {
             lower.unwrap()
         } else {
@@ -84,6 +98,25 @@ impl<T> Interpolate<T> for Nearest {
     }
 }
 
+impl<T> Interpolate<T> for Midpoint
+    where T: Add<T, Output = T> + Div<T, Output = T> + Clone + FromPrimitive
+{
+    fn needs_lower(_q: f64, _len: usize) -> bool {
+        true
+    }
+    fn needs_upper(_q: f64, _len: usize) -> bool {
+        true
+    }
+    fn interpolate<D>(lower: Option<Array<T, D>>,
+                      upper: Option<Array<T, D>>,
+                      _q: f64, _len: usize) -> Array<T, D>
+        where D: Dimension
+    {
+        let denom = T::from_u8(2).unwrap();
+        (lower.unwrap() + upper.unwrap()).mapv_into(|x| x / denom.clone())
+    }
+}
+
 /// Numerical methods for arrays.
 impl<A, S, D> ArrayBase<S, D>
     where S: Data<Elem=A>,
@@ -214,7 +247,7 @@ impl<A, S, D> ArrayBase<S, D>
         where D: RemoveAxis,
               A: Ord + Clone + Zero,
               S: DataMut,
-              I: Interpolate<Array<A, D::Smaller>>,
+              I: Interpolate<A>,
     {
         assert!((0. <= q) && (q <= 1.));
         let mut lower = None;

From 1d094e8f5f966b539be9d0b5845c1586ecd251c4 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 15:01:30 +0100
Subject: [PATCH 89/93] Defined Linear.

---
 src/numeric/impl_numeric.rs | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index fab40e406..9b717bef3 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -47,6 +47,19 @@ pub struct Upper;
 pub struct Lower;
 pub struct Nearest;
 pub struct Midpoint;
+pub struct Linear<T> {
+    coeff: T,
+}
+
+impl<T> Linear<T>
+    where T: FromPrimitive + Ord
+{
+    pub fn new(coeff: T) -> Self {
+        assert!(coeff <= T::from_u8(1).unwrap());
+        assert!(coeff >= T::from_u8(0).unwrap());
+        Linear { coeff }
+    }
+}
 
 impl<T> Interpolate<T> for Upper {
     fn needs_lower(_q: f64, _len: usize) -> bool {

From f5f77eab7b1be188aa23da95113728f471b60ac4 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 15:15:58 +0100
Subject: [PATCH 90/93] Added Interpolate implementation for Linear.

---
 src/numeric/impl_numeric.rs | 44 +++++++++++++++++++++++++------------
 1 file changed, 30 insertions(+), 14 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 9b717bef3..d97098340 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -6,7 +6,7 @@
 // option. This file may not be copied, modified, or distributed
 // except according to those terms.
 
-use std::ops::{Add, Div};
+use std::ops::{Add, Sub, Div, Mul};
 use libnum::{self, Zero, Float, FromPrimitive};
 use itertools::free::enumerate;
 
@@ -34,6 +34,11 @@ pub trait Interpolate<T> {
     fn upper_index(q: f64, len: usize) -> usize {
         Self::float_percentile_index(q, len).ceil() as usize
     }
+
+    fn float_percentile_index_fraction(q: f64, len: usize) -> f64 {
+        Self::float_percentile_index(q, len) - (Self::lower_index(q, len) as f64)
+    }
+
     fn needs_lower(q: f64, len: usize) -> bool;
     fn needs_upper(q: f64, len: usize) -> bool;
     fn interpolate<D>(lower: Option<Array<T, D>>,
@@ -47,19 +52,7 @@ pub struct Upper;
 pub struct Lower;
 pub struct Nearest;
 pub struct Midpoint;
-pub struct Linear<T> {
-    coeff: T,
-}
-
-impl<T> Linear<T>
-    where T: FromPrimitive + Ord
-{
-    pub fn new(coeff: T) -> Self {
-        assert!(coeff <= T::from_u8(1).unwrap());
-        assert!(coeff >= T::from_u8(0).unwrap());
-        Linear { coeff }
-    }
-}
+pub struct Linear;
 
 impl<T> Interpolate<T> for Upper {
     fn needs_lower(_q: f64, _len: usize) -> bool {
@@ -130,6 +123,29 @@ impl<T> Interpolate<T> for Midpoint
     }
 }
 
+impl<T> Interpolate<T> for Linear
+    where T: Add<T, Output = T> + Sub<T, Output = T> + Mul<T, Output = T> + Clone + FromPrimitive
+{
+    fn needs_lower(_q: f64, _len: usize) -> bool {
+        true
+    }
+    fn needs_upper(_q: f64, _len: usize) -> bool {
+        true
+    }
+    fn interpolate<D>(lower: Option<Array<T, D>>,
+                      upper: Option<Array<T, D>>,
+                      q: f64, len: usize) -> Array<T, D>
+        where D: Dimension
+    {
+        let fraction = T::from_f64(
+            <Self as Interpolate<T>>::float_percentile_index_fraction(q, len)
+        ).unwrap();
+        let a = lower.unwrap().mapv_into(|x| x * fraction.clone());
+        let b = upper.unwrap().mapv_into(|x| x * (T::from_u8(1).unwrap() - fraction.clone()));
+        a + b
+    }
+}
+
 /// Numerical methods for arrays.
 impl<A, S, D> ArrayBase<S, D>
     where S: Data<Elem=A>,

From ff51479c42e7904abd4807fe7b85d6bc64699d9c Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 15:21:50 +0100
Subject: [PATCH 91/93] Fixed docs!

---
 src/numeric/impl_numeric.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index d97098340..e5638e249 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -249,8 +249,8 @@ impl<A, S, D> ArrayBase<S, D>
     /// as the element that would be indexed as `(N-1)q` if the lane were to be sorted
     /// in increasing order.
     /// If `(N-1)q` is not an integer the desired percentile lies between
-    /// two data points: we return the lower, nearest or higher datapoint depending
-    /// on `interpolation_strategy`.
+    /// two data points: we return the lower, nearest, higher or interpolated
+    /// value depending on the type `Interpolate` bound `I`.
     ///
     /// Some examples:
     /// - `q=0.` returns the minimum along each 1-dimensional lane;

From d4e7de73090b41e60e7e73fcd2b9247ca896f6b1 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 15:48:19 +0100
Subject: [PATCH 92/93] Removed unused constraint on A.

---
 src/numeric/impl_numeric.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index e5638e249..3bfe6347d 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -274,7 +274,7 @@ impl<A, S, D> ArrayBase<S, D>
     /// `0.` and `1.` (inclusive).
     pub fn percentile_axis_mut<I>(&mut self, axis: Axis, q: f64) -> Array<A, D::Smaller>
         where D: RemoveAxis,
-              A: Ord + Clone + Zero,
+              A: Ord + Clone,
               S: DataMut,
               I: Interpolate<A>,
     {

From caf58ee99867d0c20b0d547e71db6bdd4c721b47 Mon Sep 17 00:00:00 2001
From: LukeMathWalker <lucapalmieri1993@hotmail.com>
Date: Sun, 2 Sep 2018 16:51:13 +0100
Subject: [PATCH 93/93] Reducing the scope of search for upper_index element if
 lower_index was already found.

---
 src/numeric/impl_numeric.rs | 16 +++++++++++++---
 1 file changed, 13 insertions(+), 3 deletions(-)

diff --git a/src/numeric/impl_numeric.rs b/src/numeric/impl_numeric.rs
index 3bfe6347d..8ffdc354c 100644
--- a/src/numeric/impl_numeric.rs
+++ b/src/numeric/impl_numeric.rs
@@ -283,14 +283,24 @@ impl<A, S, D> ArrayBase<S, D>
         let mut upper = None;
         let axis_len = self.len_of(axis);
         if I::needs_lower(q, axis_len) {
+            let lower_index = I::lower_index(q, axis_len);
             lower = Some(
                 self.map_axis_mut(
                     axis,
-                    |mut x| x.sorted_get_mut(I::lower_index(q, axis_len))
+                    |mut x| x.sorted_get_mut(lower_index)
                 )
             );
-        };
-        if I::needs_upper(q, axis_len) {
+            if I::needs_upper(q, axis_len) {
+                let upper_index = I::upper_index(q, axis_len);
+                let relative_upper_index = upper_index - lower_index;
+                upper = Some(
+                    self.map_axis_mut(
+                        axis,
+                        |mut x| x.slice_mut(s![lower_index..]).sorted_get_mut(relative_upper_index)
+                    )
+                );
+            };
+        } else {
             upper = Some(
                 self.map_axis_mut(
                     axis,