Update kdtree and fix KDE radius

Cargo.toml

@@ -53,7 +53,7 @@ optional = true
 default-features = false
 
 [dependencies.kdtree]
-version = "0.7.0"
+version = "0.8.1"
 optional = true
 
 [dev-dependencies]

src/density/kde.rs

@@ -1,7 +1,7 @@
 use kdtree::distance::squared_euclidean;
 
 use crate::{
-    density::{Container, DensityError, nearest_neighbors},
+    density::{Container, DensityError, nearest_neighbors, neighborhood_radius},
     function::kernel::{Gaussian, Kernel},
 };
 
@@ -30,7 +30,7 @@ where
     if neighbors.is_empty() {
         Err(DensityError::EmptyNeighborhood)
     } else {
-        let radius = neighbors.last().unwrap().sqrt(); // safe to unwrap here since `neighbors` is not empty
+        let radius = neighborhood_radius(&neighbors).unwrap();
         let d = x.length() as i32;
         Ok((1. / (n_samples * radius.powi(d)))
             * samples

src/density/knn.rs

@@ -1,6 +1,6 @@
 use super::Container;
 use crate::{
-    density::{DensityError, nearest_neighbors},
+    density::{DensityError, nearest_neighbors, neighborhood_radius},
     function::gamma::gamma,
 };
 use core::f64::consts::PI;
@@ -30,7 +30,7 @@ where
     if neighbors.is_empty() {
         Err(DensityError::EmptyNeighborhood)
     } else {
-        let radius = neighbors.last().unwrap().sqrt();
+        let radius = neighborhood_radius(&neighbors).unwrap();
         let d = x.length() as f64;
         Ok((k / n_samples) * (gamma(d / 2. + 1.) / (PI.powf(d / 2.) * radius.powf(d))))
     }

src/density/mod.rs

@@ -79,6 +79,14 @@ impl_container!(
 );
 pub type NearestNeighbors = (Vec<f64>, f64);
 
+pub(crate) fn neighborhood_radius(neighbors: &[f64]) -> Option<f64> {
+    neighbors
+        .iter()
+        .copied()
+        .max_by(f64::total_cmp)
+        .map(f64::sqrt)
+}
+
 pub(crate) fn nearest_neighbors<S, X>(
     x: &X,
     samples: &S,
@@ -118,6 +126,11 @@ mod tests {
 
     use super::*;
 
+    #[test]
+    fn neighborhood_radius_uses_farthest_distance() {
+        assert_eq!(neighborhood_radius(&[4.0, 1.0, 9.0, 2.0]), Some(3.0));
+    }
+
     #[test]
     fn test_vec_container() {
         let v1 = vec![1.0, 2.0, 3.0];