Add another pairwise distance algorithm

* Fix #245: return error for NaN in naive bayes
* Implement error handling for NaN values in NBayes predict:
* general behaviour has been kept unchanged according to original tests in `mod.rs`
* aka: error is returned only if all the predicted probabilities are NaN
* Add tests
* Add test with static values
* Add test for numerical stability with numpy

.github/workflows/ci.yml

@@ -36,7 +36,7 @@ jobs:
       - name: Install Rust toolchain
         uses: actions-rs/toolchain@v1
         with:
-          toolchain: stable
+          toolchain: 1.81 # 1.82 seems to break wasm32 tests https://github.com/rustwasm/wasm-bindgen/issues/4274
           target: ${{ matrix.platform.target }}
           profile: minimal
           default: true

Cargo.toml

@@ -48,7 +48,7 @@ getrandom = { version = "0.2.8", optional = true }
 wasm-bindgen-test = "0.3"
 
 [dev-dependencies]
-itertools = "0.12.0"
+itertools = "0.13.0"
 serde_json = "1.0"
 bincode = "1.3.1"
 

src/algorithm/neighbour/cover_tree.rs

@@ -124,7 +124,7 @@ impl<T: Debug + PartialEq, D: Distance<T>> CoverTree<T, D> {
         current_cover_set.push((d, &self.root));
 
         let mut heap = HeapSelection::with_capacity(k);
-        heap.add(std::f64::MAX);
+        heap.add(f64::MAX);
 
         let mut empty_heap = true;
         if !self.identical_excluded || self.get_data_value(self.root.idx) != p {
@@ -145,7 +145,7 @@ impl<T: Debug + PartialEq, D: Distance<T>> CoverTree<T, D> {
                     }
 
                     let upper_bound = if empty_heap {
-                        std::f64::INFINITY
+                        f64::INFINITY
                     } else {
                         *heap.peek()
                     };
@@ -291,7 +291,7 @@ impl<T: Debug + PartialEq, D: Distance<T>> CoverTree<T, D> {
         } else {
             let max_dist = self.max(point_set);
             let next_scale = (max_scale - 1).min(self.get_scale(max_dist));
-            if next_scale == std::i64::MIN {
+            if next_scale == i64::MIN {
                 let mut children: Vec<Node> = Vec::new();
                 let mut leaf = self.new_leaf(p);
                 children.push(leaf);
@@ -435,7 +435,7 @@ impl<T: Debug + PartialEq, D: Distance<T>> CoverTree<T, D> {
 
     fn get_scale(&self, d: f64) -> i64 {
         if d == 0f64 {
-            std::i64::MIN
+            i64::MIN
         } else {
             (self.inv_log_base * d.ln()).ceil() as i64
         }

src/algorithm/neighbour/eppstein.rs

@@ -0,0 +1,219 @@
+//! This module provides FastPair, a data-structure for efficiently tracking the dynamic
+//! closest pairs in a set of points, with an example usage in hierarchical clustering.[2][3][5]
+//!
+//! ## Purpose
+//!
+//! FastPair allows quick retrieval of the nearest neighbor for each data point by maintaining
+//! a "conga line" of closest pairs. Each point retains a link to its known nearest neighbor,
+//! and updates in the data structure propagate accordingly. This can be leveraged in
+//! agglomerative clustering steps, where merging or insertion of new points must be reflected
+//! in nearest-neighbor relationships.
+//!
+//! ## Example
+//!
+//! ```
+//! use smartcore::metrics::distance::PairwiseDistance;
+//! use smartcore::linalg::basic::matrix::DenseMatrix;
+//! use smartcore::algorithm::neighbour::fastpair::FastPair;
+//!
+//! let x = DenseMatrix::from_2d_array(&[
+//!     &[5.1, 3.5, 1.4, 0.2],
+//!     &[4.9, 3.0, 1.4, 0.2],
+//!     &[4.7, 3.2, 1.3, 0.2],
+//!     &[4.6, 3.1, 1.5, 0.2],
+//!     &[5.0, 3.6, 1.4, 0.2],
+//!     &[5.4, 3.9, 1.7, 0.4],
+//! ]).unwrap();
+//!
+//! let fastpair = FastPair::new(&x).unwrap();
+//! let closest = fastpair.closest_pair();
+//! println!("Closest pair: {:?}", closest);
+//! ```
+use std::collections::HashMap;
+
+use num::Bounded;
+
+use crate::error::{Failed, FailedError};
+use crate::linalg::basic::arrays::{Array, Array1, Array2};
+use crate::metrics::distance::euclidian::Euclidian;
+use crate::metrics::distance::PairwiseDistance;
+use crate::numbers::floatnum::FloatNumber;
+use crate::numbers::realnum::RealNumber;
+
+/// Eppstein dynamic closet-pair structure
+/// 'M' can be a matrix-like trait that provides row access
+#[derive(Debug)]
+pub struct EppsteinDCP<'a, T: RealNumber + FloatNumber, M: Array2<T>> {
+    samples: &'a M,
+    // "buckets" store, for each row, a small structure recording potential neighbors
+    neighbors: HashMap<usize, PairwiseDistance<T>>,
+}
+
+impl<'a, T: RealNumber + FloatNumber, M: Array2<T>> EppsteinDCP<'a, T, M> {
+    /// Creates a new EppsteinDCP instance with the given data
+    pub fn new(m: &'a M) -> Result<Self, Failed> {
+        if m.shape().0 < 3 {
+            return Err(Failed::because(
+                FailedError::FindFailed,
+                "min number of rows should be 3",
+            ));
+        }
+
+        let mut this = Self {
+            samples: m,
+            neighbors: HashMap::with_capacity(m.shape().0),
+        };
+        this.initialize();
+        Ok(this)
+    }
+
+    /// Build an initial "conga line" or chain of potential neighbors
+    /// akin to Eppstein’s technique[2].
+    fn initialize(&mut self) {
+        let n = self.samples.shape().0;
+        if n < 2 {
+            return;
+        }
+        // Assign each row i some large distance by default
+        for i in 0..n {
+            self.neighbors.insert(
+                i,
+                PairwiseDistance {
+                    node: i,
+                    neighbour: None,
+                    distance: Some(<T as Bounded>::max_value()),
+                },
+            );
+        }
+        // Example: link each i to the next, forming a chain
+        // (depending on the actual Eppstein approach, can refine)
+        for i in 0..(n - 1) {
+            let dist = self.compute_dist(i, i + 1);
+            self.neighbors.entry(i).and_modify(|pd| {
+                pd.neighbour = Some(i + 1);
+                pd.distance = Some(dist);
+            });
+        }
+        // Potential refinement steps omitted for brevity
+    }
+
+    /// Insert a point into the structure.
+    pub fn insert(&mut self, row_idx: usize) {
+        // Expand data, find neighbor to link with
+        // For example, link row_idx to nearest among existing
+        let mut best_neighbor = None;
+        let mut best_d = <T as Bounded>::max_value();
+        for (i, _) in &self.neighbors {
+            let d = self.compute_dist(*i, row_idx);
+            if d < best_d {
+                best_d = d;
+                best_neighbor = Some(*i);
+            }
+        }
+        self.neighbors.insert(
+            row_idx,
+            PairwiseDistance {
+                node: row_idx,
+                neighbour: best_neighbor,
+                distance: Some(best_d),
+            },
+        );
+        // For the best_neighbor, you might want to see if row_idx becomes closer
+        if let Some(kn) = best_neighbor {
+            let dist = self.compute_dist(row_idx, kn);
+            let entry = self.neighbors.get_mut(&kn).unwrap();
+            if dist < entry.distance.unwrap() {
+                entry.neighbour = Some(row_idx);
+                entry.distance = Some(dist);
+            }
+        }
+    }
+
+    /// For hierarchical clustering, discover minimal pairs, then merge
+    pub fn closest_pair(&self) -> Option<PairwiseDistance<T>> {
+        let mut min_pair: Option<PairwiseDistance<T>> = None;
+        for (_, pd) in &self.neighbors {
+            if let Some(d) = pd.distance {
+                if min_pair.is_none() || d < min_pair.as_ref().unwrap().distance.unwrap() {
+                    min_pair = Some(pd.clone());
+                }
+            }
+        }
+        min_pair
+    }
+
+    fn compute_dist(&self, i: usize, j: usize) -> T {
+        // Example: Euclidean
+        let row_i = self.samples.get_row(i);
+        let row_j = self.samples.get_row(j);
+        row_i
+            .iterator(0)
+            .zip(row_j.iterator(0))
+            .map(|(a, b)| (*a - *b) * (*a - *b))
+            .sum()
+    }
+}
+
+/// Simple usage
+#[cfg(test)]
+mod tests_eppstein {
+    use super::*;
+    use crate::linalg::basic::matrix::DenseMatrix;
+
+    #[test]
+    fn test_eppstein() {
+        let matrix =
+            DenseMatrix::from_2d_array(&[&vec![1.0, 2.0], &vec![2.0, 2.0], &vec![5.0, 3.0]])
+                .unwrap();
+        let mut dcp = EppsteinDCP::new(&matrix).unwrap();
+        dcp.insert(2);
+        let cp = dcp.closest_pair();
+        assert!(cp.is_some());
+    }
+
+    #[test]
+    fn compare_fastpair_eppstein() {
+        use crate::algorithm::neighbour::fastpair::FastPair;
+        // Assuming EppsteinDCP is implemented in a similar module
+        use crate::algorithm::neighbour::eppstein::EppsteinDCP;
+
+        // Create a static example matrix
+        let x = DenseMatrix::from_2d_array(&[
+            &[5.1, 3.5, 1.4, 0.2],
+            &[4.9, 3.0, 1.4, 0.2],
+            &[4.7, 3.2, 1.3, 0.2],
+            &[4.6, 3.1, 1.5, 0.2],
+            &[5.0, 3.6, 1.4, 0.2],
+            &[5.4, 3.9, 1.7, 0.4],
+            &[4.6, 3.4, 1.4, 0.3],
+            &[5.0, 3.4, 1.5, 0.2],
+            &[4.4, 2.9, 1.4, 0.2],
+            &[4.9, 3.1, 1.5, 0.1],
+        ])
+        .unwrap();
+
+        // Build FastPair
+        let fastpair = FastPair::new(&x).unwrap();
+        let pair_fastpair = fastpair.closest_pair();
+
+        // Build EppsteinDCP
+        let eppstein = EppsteinDCP::new(&x).unwrap();
+        let pair_eppstein = eppstein.closest_pair();
+
+        // Compare the results
+        assert_eq!(pair_fastpair.node, pair_eppstein.as_ref().unwrap().node);
+        assert_eq!(
+            pair_fastpair.neighbour.unwrap(),
+            pair_eppstein.as_ref().unwrap().neighbour.unwrap()
+        );
+
+        // Use a small epsilon for floating-point comparison
+        let epsilon = 1e-9;
+        let diff: f64 =
+            pair_fastpair.distance.unwrap() - pair_eppstein.as_ref().unwrap().distance.unwrap();
+        assert!(diff.abs() < epsilon);
+
+        println!("FastPair result: {:?}", pair_fastpair);
+        println!("EppsteinDCP result: {:?}", pair_eppstein);
+    }
+}

src/algorithm/neighbour/fastpair.rs

@@ -52,10 +52,8 @@ pub struct FastPair<'a, T: RealNumber + FloatNumber, M: Array2<T>> {
 }
 
 impl<'a, T: RealNumber + FloatNumber, M: Array2<T>> FastPair<'a, T, M> {
-    ///
     /// Constructor
-    /// Instantiate and inizialise the algorithm
-    ///
+    /// Instantiate and initialize the algorithm
     pub fn new(m: &'a M) -> Result<Self, Failed> {
         if m.shape().0 < 3 {
             return Err(Failed::because(
@@ -74,10 +72,8 @@ impl<'a, T: RealNumber + FloatNumber, M: Array2<T>> FastPair<'a, T, M> {
         Ok(init)
     }
 
-    ///
     /// Initialise `FastPair` by passing a `Array2`.
     /// Build a FastPairs data-structure from a set of (new) points.
-    ///
     fn init(&mut self) {
         // basic measures
         let len = self.samples.shape().0;
@@ -158,9 +154,7 @@ impl<'a, T: RealNumber + FloatNumber, M: Array2<T>> FastPair<'a, T, M> {
         self.neighbours = neighbours;
     }
 
-    ///
     /// Find closest pair by scanning list of nearest neighbors.
-    ///
     #[allow(dead_code)]
     pub fn closest_pair(&self) -> PairwiseDistance<T> {
         let mut a = self.neighbours[0]; // Start with first point
@@ -179,6 +173,21 @@ impl<'a, T: RealNumber + FloatNumber, M: Array2<T>> FastPair<'a, T, M> {
         }
     }
 
+    ///
+    /// Return order dissimilarities from closest to furthest
+    ///
+    #[allow(dead_code)]
+    pub fn ordered_pairs(&self) -> std::vec::IntoIter<&PairwiseDistance<T>> {
+        // improvement: implement this to return `impl Iterator<Item = &PairwiseDistance<T>>`
+        // need to implement trait `Iterator` for `Vec<&PairwiseDistance<T>>`
+        let mut distances = self
+            .distances
+            .values()
+            .collect::<Vec<&PairwiseDistance<T>>>();
+        distances.sort_by(|a, b| a.partial_cmp(b).unwrap());
+        distances.into_iter()
+    }
+
     //
     // Compute distances from input to all other points in data-structure.
     // input is the row index of the sample matrix
@@ -217,10 +226,10 @@ mod tests_fastpair {
     use super::*;
     use crate::linalg::basic::{arrays::Array, matrix::DenseMatrix};
 
-    ///
     /// Brute force algorithm, used only for comparison and testing
-    ///
-    pub fn closest_pair_brute(fastpair: &FastPair<f64, DenseMatrix<f64>>) -> PairwiseDistance<f64> {
+    pub fn closest_pair_brute(
+        fastpair: &FastPair<'_, f64, DenseMatrix<f64>>,
+    ) -> PairwiseDistance<f64> {
         use itertools::Itertools;
         let m = fastpair.samples.shape().0;
 
@@ -594,4 +603,103 @@ mod tests_fastpair {
 
         assert_eq!(closest, min_dissimilarity);
     }
+
+    #[test]
+    fn fastpair_ordered_pairs() {
+        let x = DenseMatrix::<f64>::from_2d_array(&[
+            &[5.1, 3.5, 1.4, 0.2],
+            &[4.9, 3.0, 1.4, 0.2],
+            &[4.7, 3.2, 1.3, 0.2],
+            &[4.6, 3.1, 1.5, 0.2],
+            &[5.0, 3.6, 1.4, 0.2],
+            &[5.4, 3.9, 1.7, 0.4],
+            &[4.9, 3.1, 1.5, 0.1],
+            &[7.0, 3.2, 4.7, 1.4],
+            &[6.4, 3.2, 4.5, 1.5],
+            &[6.9, 3.1, 4.9, 1.5],
+            &[5.5, 2.3, 4.0, 1.3],
+            &[6.5, 2.8, 4.6, 1.5],
+            &[4.6, 3.4, 1.4, 0.3],
+            &[5.0, 3.4, 1.5, 0.2],
+            &[4.4, 2.9, 1.4, 0.2],
+        ])
+        .unwrap();
+        let fastpair = FastPair::new(&x).unwrap();
+
+        let ordered = fastpair.ordered_pairs();
+
+        let mut previous: f64 = -1.0;
+        for p in ordered {
+            if previous == -1.0 {
+                previous = p.distance.unwrap();
+            } else {
+                let current = p.distance.unwrap();
+                assert!(current >= previous);
+                previous = current;
+            }
+        }
+    }
+
+    #[test]
+    fn test_empty_set() {
+        let empty_matrix = DenseMatrix::<f64>::zeros(0, 0);
+        let result = FastPair::new(&empty_matrix);
+        assert!(result.is_err());
+        if let Err(e) = result {
+            assert_eq!(
+                e,
+                Failed::because(FailedError::FindFailed, "min number of rows should be 3")
+            );
+        }
+    }
+
+    #[test]
+    fn test_single_point() {
+        let single_point = DenseMatrix::from_2d_array(&[&[1.0, 2.0, 3.0]]).unwrap();
+        let result = FastPair::new(&single_point);
+        assert!(result.is_err());
+        if let Err(e) = result {
+            assert_eq!(
+                e,
+                Failed::because(FailedError::FindFailed, "min number of rows should be 3")
+            );
+        }
+    }
+
+    #[test]
+    fn test_two_points() {
+        let two_points = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        let result = FastPair::new(&two_points);
+        assert!(result.is_err());
+        if let Err(e) = result {
+            assert_eq!(
+                e,
+                Failed::because(FailedError::FindFailed, "min number of rows should be 3")
+            );
+        }
+    }
+
+    #[test]
+    fn test_three_identical_points() {
+        let identical_points =
+            DenseMatrix::from_2d_array(&[&[1.0, 1.0], &[1.0, 1.0], &[1.0, 1.0]]).unwrap();
+        let result = FastPair::new(&identical_points);
+        assert!(result.is_ok());
+        let fastpair = result.unwrap();
+        let closest_pair = fastpair.closest_pair();
+        assert_eq!(closest_pair.distance, Some(0.0));
+    }
+
+    #[test]
+    fn test_result_unwrapping() {
+        let valid_matrix =
+            DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0], &[5.0, 6.0], &[7.0, 8.0]])
+                .unwrap();
+
+        let result = FastPair::new(&valid_matrix);
+        assert!(result.is_ok());
+
+        // This should not panic
+        let _fastpair = result.unwrap();
+    }
 }

src/algorithm/neighbour/linear_search.rs

@@ -61,7 +61,7 @@ impl<T, D: Distance<T>> LinearKNNSearch<T, D> {
 
         for _ in 0..k {
             heap.add(KNNPoint {
-                distance: std::f64::INFINITY,
+                distance: f64::INFINITY,
                 index: None,
             });
         }
@@ -215,7 +215,7 @@ mod tests {
         };
 
         let point_inf = KNNPoint {
-            distance: std::f64::INFINITY,
+            distance: f64::INFINITY,
             index: Some(3),
         };
 

src/algorithm/neighbour/mod.rs

@@ -41,7 +41,9 @@ use serde::{Deserialize, Serialize};
 pub(crate) mod bbd_tree;
 /// tree data structure for fast nearest neighbor search
 pub mod cover_tree;
-/// fastpair closest neighbour algorithm
+/// eppstein pairwise closest neighbour algorithm
+pub mod eppstein;
+/// fastpair pairwise closest neighbour algorithm
 pub mod fastpair;
 /// very simple algorithm that sequentially checks each element of the list until a match is found or the whole list has been searched.
 pub mod linear_search;

src/algorithm/sort/heap_select.rs

@@ -133,7 +133,7 @@ mod tests {
     #[test]
     fn test_add1() {
         let mut heap = HeapSelection::with_capacity(3);
-        heap.add(std::f64::INFINITY);
+        heap.add(f64::INFINITY);
         heap.add(-5f64);
         heap.add(4f64);
         heap.add(-1f64);
@@ -151,7 +151,7 @@ mod tests {
     #[test]
     fn test_add2() {
         let mut heap = HeapSelection::with_capacity(3);
-        heap.add(std::f64::INFINITY);
+        heap.add(f64::INFINITY);
         heap.add(0.0);
         heap.add(8.4852);
         heap.add(5.6568);

src/algorithm/sort/quick_sort.rs

@@ -3,6 +3,7 @@ use num_traits::Num;
 pub trait QuickArgSort {
     fn quick_argsort_mut(&mut self) -> Vec<usize>;
 
+    #[allow(dead_code)]
     fn quick_argsort(&self) -> Vec<usize>;
 }
 

src/cluster/kmeans.rs

@@ -96,7 +96,7 @@ impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>> PartialEq for KMeans<
                     return false;
                 }
                 for j in 0..self.centroids[i].len() {
-                    if (self.centroids[i][j] - other.centroids[i][j]).abs() > std::f64::EPSILON {
+                    if (self.centroids[i][j] - other.centroids[i][j]).abs() > f64::EPSILON {
                         return false;
                     }
                 }
@@ -270,7 +270,7 @@ impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>> KMeans<TX, TY, X, Y>
 
         let (n, d) = data.shape();
 
-        let mut distortion = std::f64::MAX;
+        let mut distortion = f64::MAX;
         let mut y = KMeans::<TX, TY, X, Y>::kmeans_plus_plus(data, parameters.k, parameters.seed);
         let mut size = vec![0; parameters.k];
         let mut centroids = vec![vec![0f64; d]; parameters.k];
@@ -331,7 +331,7 @@ impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>> KMeans<TX, TY, X, Y>
         let mut row = vec![0f64; x.shape().1];
 
         for i in 0..n {
-            let mut min_dist = std::f64::MAX;
+            let mut min_dist = f64::MAX;
             let mut best_cluster = 0;
 
             for j in 0..self.k {
@@ -361,7 +361,7 @@ impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>> KMeans<TX, TY, X, Y>
             .cloned()
             .collect();
 
-        let mut d = vec![std::f64::MAX; n];
+        let mut d = vec![f64::MAX; n];
         let mut row = vec![TX::zero(); data.shape().1];
 
         for j in 1..k {

src/ensemble/random_forest_classifier.rs

@@ -580,37 +580,6 @@ impl<TX: FloatNumber + PartialOrd, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY
         which_max(&result)
     }
 
-    /// Predict the per-class probabilties for each observation.
-    /// The probability is calculated as the fraction of trees that predicted a given class
-    pub fn predict_proba<R: Array2<f64>>(&self, x: &X) -> Result<R, Failed> {
-        let mut result: R = R::zeros(x.shape().0, self.classes.as_ref().unwrap().len());
-
-        let (n, _) = x.shape();
-
-        for i in 0..n {
-            let row_probs = self.predict_proba_for_row(x, i);
-
-            for (j, item) in row_probs.iter().enumerate() {
-                result.set((i, j), *item);
-            }
-        }
-
-        Ok(result)
-    }
-
-    fn predict_proba_for_row(&self, x: &X, row: usize) -> Vec<f64> {
-        let mut result = vec![0; self.classes.as_ref().unwrap().len()];
-
-        for tree in self.trees.as_ref().unwrap().iter() {
-            result[tree.predict_for_row(x, row)] += 1;
-        }
-
-        result
-            .iter()
-            .map(|n| *n as f64 / self.trees.as_ref().unwrap().len() as f64)
-            .collect()
-    }
-
     fn sample_with_replacement(y: &[usize], num_classes: usize, rng: &mut impl Rng) -> Vec<usize> {
         let class_weight = vec![1.; num_classes];
         let nrows = y.len();
@@ -638,7 +607,6 @@ impl<TX: FloatNumber + PartialOrd, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY
 #[cfg(test)]
 mod tests {
     use super::*;
-    use crate::linalg::basic::arrays::Array;
     use crate::linalg::basic::matrix::DenseMatrix;
     use crate::metrics::*;
 
@@ -831,69 +799,4 @@ mod tests {
 
         assert_eq!(forest, deserialized_forest);
     }
-
-    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
-    #[test]
-    fn fit_predict_probabilities() {
-        let x = DenseMatrix::<f64>::from_2d_array(&[
-            &[5.1, 3.5, 1.4, 0.2],
-            &[4.9, 3.0, 1.4, 0.2],
-            &[4.7, 3.2, 1.3, 0.2],
-            &[4.6, 3.1, 1.5, 0.2],
-            &[5.0, 3.6, 1.4, 0.2],
-            &[5.4, 3.9, 1.7, 0.4],
-            &[4.6, 3.4, 1.4, 0.3],
-            &[5.0, 3.4, 1.5, 0.2],
-            &[4.4, 2.9, 1.4, 0.2],
-            &[4.9, 3.1, 1.5, 0.1],
-            &[7.0, 3.2, 4.7, 1.4],
-            &[6.4, 3.2, 4.5, 1.5],
-            &[6.9, 3.1, 4.9, 1.5],
-            &[5.5, 2.3, 4.0, 1.3],
-            &[6.5, 2.8, 4.6, 1.5],
-            &[5.7, 2.8, 4.5, 1.3],
-            &[6.3, 3.3, 4.7, 1.6],
-            &[4.9, 2.4, 3.3, 1.0],
-            &[6.6, 2.9, 4.6, 1.3],
-            &[5.2, 2.7, 3.9, 1.4],
-        ]);
-        let y = vec![0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1];
-
-        let classifier = RandomForestClassifier::fit(
-            &x,
-            &y,
-            RandomForestClassifierParameters {
-                criterion: SplitCriterion::Gini,
-                max_depth: None,
-                min_samples_leaf: 1,
-                min_samples_split: 2,
-                n_trees: 100, // this is n_estimators in sklearn
-                m: Option::None,
-                keep_samples: false,
-                seed: 0,
-            },
-        )
-        .unwrap();
-
-        println!("{:?}", classifier.classes);
-
-        let results: DenseMatrix<f64> = classifier.predict_proba(&x).unwrap();
-        println!("{:?}", x.shape());
-        println!("{:?}", results);
-        println!("{:?}", results.shape());
-
-        assert_eq!(
-            results,
-            DenseMatrix::<f64>::new(
-                20,
-                2,
-                vec![
-                    1.0, 0.0, 0.78, 0.22, 0.95, 0.05, 0.82, 0.18, 1.0, 0.0, 0.92, 0.08, 0.99, 0.01,
-                    0.96, 0.04, 0.36, 0.64, 0.33, 0.67, 0.02, 0.98, 0.02, 0.98, 0.0, 1.0, 0.0, 1.0,
-                    0.0, 1.0, 0.0, 1.0, 0.03, 0.97, 0.05, 0.95, 0.0, 1.0, 0.02, 0.98
-                ],
-                true
-            )
-        );
-    }
 }

src/lib.rs

@@ -7,7 +7,6 @@
     clippy::approx_constant
 )]
 #![warn(missing_docs)]
-#![warn(rustdoc::missing_doc_code_examples)]
 
 //! # smartcore
 //!

src/linalg/basic/arrays.rs

@@ -265,11 +265,11 @@ pub trait ArrayView1<T: Debug + Display + Copy + Sized>: Array<T, usize> {
         if p.is_infinite() && p.is_sign_positive() {
             self.iterator(0)
                 .map(|x| x.to_f64().unwrap().abs())
-                .fold(std::f64::NEG_INFINITY, |a, b| a.max(b))
+                .fold(f64::NEG_INFINITY, |a, b| a.max(b))
         } else if p.is_infinite() && p.is_sign_negative() {
             self.iterator(0)
                 .map(|x| x.to_f64().unwrap().abs())
-                .fold(std::f64::INFINITY, |a, b| a.min(b))
+                .fold(f64::INFINITY, |a, b| a.min(b))
         } else {
             let mut norm = 0f64;
 
@@ -558,11 +558,11 @@ pub trait ArrayView2<T: Debug + Display + Copy + Sized>: Array<T, (usize, usize)
         if p.is_infinite() && p.is_sign_positive() {
             self.iterator(0)
                 .map(|x| x.to_f64().unwrap().abs())
-                .fold(std::f64::NEG_INFINITY, |a, b| a.max(b))
+                .fold(f64::NEG_INFINITY, |a, b| a.max(b))
         } else if p.is_infinite() && p.is_sign_negative() {
             self.iterator(0)
                 .map(|x| x.to_f64().unwrap().abs())
-                .fold(std::f64::INFINITY, |a, b| a.min(b))
+                .fold(f64::INFINITY, |a, b| a.min(b))
         } else {
             let mut norm = 0f64;
 
@@ -731,34 +731,34 @@ pub trait MutArrayView1<T: Debug + Display + Copy + Sized>:
 pub trait MutArrayView2<T: Debug + Display + Copy + Sized>:
     MutArray<T, (usize, usize)> + ArrayView2<T>
 {
-    ///
+    /// copy values from another array
     fn copy_from(&mut self, other: &dyn Array<T, (usize, usize)>) {
         self.iterator_mut(0)
             .zip(other.iterator(0))
             .for_each(|(s, o)| *s = *o);
     }
-    ///
+    /// update view with absolute values
     fn abs_mut(&mut self)
     where
         T: Number + Signed,
     {
         self.iterator_mut(0).for_each(|v| *v = v.abs());
     }
-    ///
+    /// update view values with opposite sign
     fn neg_mut(&mut self)
     where
         T: Number + Neg<Output = T>,
     {
         self.iterator_mut(0).for_each(|v| *v = -*v);
     }
-    ///
+    /// update view values at power `p`
     fn pow_mut(&mut self, p: T)
     where
         T: RealNumber,
     {
         self.iterator_mut(0).for_each(|v| *v = v.powf(p));
     }
-    ///
+    /// scale view values
     fn scale_mut(&mut self, mean: &[T], std: &[T], axis: u8)
     where
         T: Number,
@@ -784,27 +784,27 @@ pub trait MutArrayView2<T: Debug + Display + Copy + Sized>:
 
 /// Trait for mutable 1D-array view
 pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized + Clone {
-    ///
+    /// return a view of the array
     fn slice<'a>(&'a self, range: Range<usize>) -> Box<dyn ArrayView1<T> + 'a>;
-    ///
+    /// return a mutable view of the array
     fn slice_mut<'a>(&'a mut self, range: Range<usize>) -> Box<dyn MutArrayView1<T> + 'a>;
-    ///
+    /// fill array with a given value
     fn fill(len: usize, value: T) -> Self
     where
         Self: Sized;
-    ///
+    /// create array from iterator
     fn from_iterator<I: Iterator<Item = T>>(iter: I, len: usize) -> Self
     where
         Self: Sized;
-    ///
+    /// create array from vector
     fn from_vec_slice(slice: &[T]) -> Self
     where
         Self: Sized;
-    ///
+    /// create array from slice
     fn from_slice(slice: &'_ dyn ArrayView1<T>) -> Self
     where
         Self: Sized;
-    ///
+    /// create a zero array
     fn zeros(len: usize) -> Self
     where
         T: Number,
@@ -812,7 +812,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
     {
         Self::fill(len, T::zero())
     }
-    ///
+    /// create an array of ones
     fn ones(len: usize) -> Self
     where
         T: Number,
@@ -820,7 +820,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
     {
         Self::fill(len, T::one())
     }
-    ///
+    /// create an array of random values
     fn rand(len: usize) -> Self
     where
         T: RealNumber,
@@ -828,7 +828,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
     {
         Self::from_iterator((0..len).map(|_| T::rand()), len)
     }
-    ///
+    /// add a scalar to the array
     fn add_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -838,7 +838,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.add_scalar_mut(x);
         result
     }
-    ///
+    /// subtract a scalar from the array
     fn sub_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -848,7 +848,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.sub_scalar_mut(x);
         result
     }
-    ///
+    /// divide a scalar from the array
     fn div_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -858,7 +858,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.div_scalar_mut(x);
         result
     }
-    ///
+    /// multiply a scalar to the array
     fn mul_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -868,7 +868,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.mul_scalar_mut(x);
         result
     }
-    ///
+    /// sum of two arrays
     fn add(&self, other: &dyn Array<T, usize>) -> Self
     where
         T: Number,
@@ -878,7 +878,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.add_mut(other);
         result
     }
-    ///
+    /// subtract two arrays
     fn sub(&self, other: &impl Array1<T>) -> Self
     where
         T: Number,
@@ -888,7 +888,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.sub_mut(other);
         result
     }
-    ///
+    /// multiply two arrays
     fn mul(&self, other: &dyn Array<T, usize>) -> Self
     where
         T: Number,
@@ -898,7 +898,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.mul_mut(other);
         result
     }
-    ///
+    /// divide two arrays
     fn div(&self, other: &dyn Array<T, usize>) -> Self
     where
         T: Number,
@@ -908,7 +908,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.div_mut(other);
         result
     }
-    ///
+    /// replace values with another array
     fn take(&self, index: &[usize]) -> Self
     where
         Self: Sized,
@@ -920,7 +920,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         );
         Self::from_iterator(index.iter().map(move |&i| *self.get(i)), index.len())
     }
-    ///
+    /// create a view of the array with absolute values
     fn abs(&self) -> Self
     where
         T: Number + Signed,
@@ -930,7 +930,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.abs_mut();
         result
     }
-    ///
+    /// create a view of the array with opposite sign
     fn neg(&self) -> Self
     where
         T: Number + Neg<Output = T>,
@@ -940,7 +940,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.neg_mut();
         result
     }
-    ///
+    /// create a view of the array with values at power `p`
     fn pow(&self, p: T) -> Self
     where
         T: RealNumber,
@@ -950,7 +950,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.pow_mut(p);
         result
     }
-    ///
+    /// apply argsort to the array
     fn argsort(&self) -> Vec<usize>
     where
         T: Number + PartialOrd,
@@ -958,12 +958,12 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         let mut v = self.clone();
         v.argsort_mut()
     }
-    ///
+    /// map values of the array
     fn map<O: Debug + Display + Copy + Sized, A: Array1<O>, F: FnMut(&T) -> O>(self, f: F) -> A {
         let len = self.shape();
         A::from_iterator(self.iterator(0).map(f), len)
     }
-    ///
+    /// apply softmax to the array
     fn softmax(&self) -> Self
     where
         T: RealNumber,
@@ -973,7 +973,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result.softmax_mut();
         result
     }
-    ///
+    /// multiply array by matrix
     fn xa(&self, a_transpose: bool, a: &dyn ArrayView2<T>) -> Self
     where
         T: Number,
@@ -1003,7 +1003,7 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
         result
     }
 
-    ///
+    /// check if two arrays are approximately equal
     fn approximate_eq(&self, other: &Self, error: T) -> bool
     where
         T: Number + RealNumber,
@@ -1015,13 +1015,13 @@ pub trait Array1<T: Debug + Display + Copy + Sized>: MutArrayView1<T> + Sized +
 
 /// Trait for mutable 2D-array view
 pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized + Clone {
-    ///
+    /// fill 2d array with a given value
     fn fill(nrows: usize, ncols: usize, value: T) -> Self;
-    ///
+    /// get a view of the 2d array
     fn slice<'a>(&'a self, rows: Range<usize>, cols: Range<usize>) -> Box<dyn ArrayView2<T> + 'a>
     where
         Self: Sized;
-    ///
+    /// get a mutable view of the 2d array
     fn slice_mut<'a>(
         &'a mut self,
         rows: Range<usize>,
@@ -1029,31 +1029,31 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
     ) -> Box<dyn MutArrayView2<T> + 'a>
     where
         Self: Sized;
-    ///
+    /// create 2d array from iterator
     fn from_iterator<I: Iterator<Item = T>>(iter: I, nrows: usize, ncols: usize, axis: u8) -> Self;
-    ///
+    /// get row from 2d array
     fn get_row<'a>(&'a self, row: usize) -> Box<dyn ArrayView1<T> + 'a>
     where
         Self: Sized;
-    ///
+    /// get column from 2d array
     fn get_col<'a>(&'a self, col: usize) -> Box<dyn ArrayView1<T> + 'a>
     where
         Self: Sized;
-    ///
+    /// create a zero 2d array
     fn zeros(nrows: usize, ncols: usize) -> Self
     where
         T: Number,
     {
         Self::fill(nrows, ncols, T::zero())
     }
-    ///
+    /// create a 2d array of ones
     fn ones(nrows: usize, ncols: usize) -> Self
     where
         T: Number,
     {
         Self::fill(nrows, ncols, T::one())
     }
-    ///
+    /// create an identity matrix
     fn eye(size: usize) -> Self
     where
         T: Number,
@@ -1066,29 +1066,29 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
 
         matrix
     }
-    ///
+    /// create a 2d array of random values
     fn rand(nrows: usize, ncols: usize) -> Self
     where
         T: RealNumber,
     {
         Self::from_iterator((0..nrows * ncols).map(|_| T::rand()), nrows, ncols, 0)
     }
-    ///
+    /// crate from 2d slice
     fn from_slice(slice: &dyn ArrayView2<T>) -> Self {
         let (nrows, ncols) = slice.shape();
         Self::from_iterator(slice.iterator(0).cloned(), nrows, ncols, 0)
     }
-    ///
+    /// create from row
     fn from_row(slice: &dyn ArrayView1<T>) -> Self {
         let ncols = slice.shape();
         Self::from_iterator(slice.iterator(0).cloned(), 1, ncols, 0)
     }
-    ///
+    /// create from column
     fn from_column(slice: &dyn ArrayView1<T>) -> Self {
         let nrows = slice.shape();
         Self::from_iterator(slice.iterator(0).cloned(), nrows, 1, 0)
     }
-    ///
+    /// transpose 2d array
     fn transpose(&self) -> Self {
         let (nrows, ncols) = self.shape();
         let mut m = Self::fill(ncols, nrows, *self.get((0, 0)));
@@ -1099,7 +1099,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         }
         m
     }
-    ///
+    /// change shape of 2d array
     fn reshape(&self, nrows: usize, ncols: usize, axis: u8) -> Self {
         let (onrows, oncols) = self.shape();
 
@@ -1110,7 +1110,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
 
         Self::from_iterator(self.iterator(0).cloned(), nrows, ncols, axis)
     }
-    ///
+    /// multiply two 2d arrays
     fn matmul(&self, other: &dyn ArrayView2<T>) -> Self
     where
         T: Number,
@@ -1136,7 +1136,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
 
         result
     }
-    ///
+    /// matrix multiplication
     fn ab(&self, a_transpose: bool, b: &dyn ArrayView2<T>, b_transpose: bool) -> Self
     where
         T: Number,
@@ -1171,7 +1171,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
             result
         }
     }
-    ///
+    /// matrix vector multiplication
     fn ax(&self, a_transpose: bool, x: &dyn ArrayView1<T>) -> Self
     where
         T: Number,
@@ -1199,7 +1199,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         }
         result
     }
-    ///
+    /// concatenate 1d array
     fn concatenate_1d<'a>(arrays: &'a [&'a dyn ArrayView1<T>], axis: u8) -> Self {
         assert!(
             axis == 1 || axis == 0,
@@ -1237,7 +1237,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
             ),
         }
     }
-    ///
+    /// concatenate 2d array
     fn concatenate_2d<'a>(arrays: &'a [&'a dyn ArrayView2<T>], axis: u8) -> Self {
         assert!(
             axis == 1 || axis == 0,
@@ -1294,7 +1294,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
             }
         }
     }
-    ///
+    /// merge 1d arrays
     fn merge_1d<'a>(&'a self, arrays: &'a [&'a dyn ArrayView1<T>], axis: u8, append: bool) -> Self {
         assert!(
             axis == 1 || axis == 0,
@@ -1362,7 +1362,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
             }
         }
     }
-    ///
+    /// Stack arrays in sequence vertically
     fn v_stack(&self, other: &dyn ArrayView2<T>) -> Self {
         let (nrows, ncols) = self.shape();
         let (other_nrows, other_ncols) = other.shape();
@@ -1378,7 +1378,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
             0,
         )
     }
-    ///
+    /// Stack arrays in sequence horizontally
     fn h_stack(&self, other: &dyn ArrayView2<T>) -> Self {
         let (nrows, ncols) = self.shape();
         let (other_nrows, other_ncols) = other.shape();
@@ -1394,20 +1394,20 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
             1,
         )
     }
-    ///
+    /// map  array values
     fn map<O: Debug + Display + Copy + Sized, A: Array2<O>, F: FnMut(&T) -> O>(self, f: F) -> A {
         let (nrows, ncols) = self.shape();
         A::from_iterator(self.iterator(0).map(f), nrows, ncols, 0)
     }
-    ///
+    /// iter rows
     fn row_iter<'a>(&'a self) -> Box<dyn Iterator<Item = Box<dyn ArrayView1<T> + 'a>> + 'a> {
         Box::new((0..self.shape().0).map(move |r| self.get_row(r)))
     }
-    ///
+    /// iter cols
     fn col_iter<'a>(&'a self) -> Box<dyn Iterator<Item = Box<dyn ArrayView1<T> + 'a>> + 'a> {
         Box::new((0..self.shape().1).map(move |r| self.get_col(r)))
     }
-    ///
+    /// take elements from 2d array
     fn take(&self, index: &[usize], axis: u8) -> Self {
         let (nrows, ncols) = self.shape();
 
@@ -1447,7 +1447,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
     fn take_column(&self, column_index: usize) -> Self {
         self.take(&[column_index], 1)
     }
-    ///
+    /// add a scalar to the array
     fn add_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -1456,7 +1456,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.add_scalar_mut(x);
         result
     }
-    ///
+    /// subtract a scalar from the array
     fn sub_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -1465,7 +1465,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.sub_scalar_mut(x);
         result
     }
-    ///
+    /// divide a scalar from the array
     fn div_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -1474,7 +1474,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.div_scalar_mut(x);
         result
     }
-    ///
+    /// multiply a scalar to the array
     fn mul_scalar(&self, x: T) -> Self
     where
         T: Number,
@@ -1483,7 +1483,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.mul_scalar_mut(x);
         result
     }
-    ///
+    /// sum of two arrays
     fn add(&self, other: &dyn Array<T, (usize, usize)>) -> Self
     where
         T: Number,
@@ -1492,7 +1492,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.add_mut(other);
         result
     }
-    ///
+    /// subtract two arrays
     fn sub(&self, other: &dyn Array<T, (usize, usize)>) -> Self
     where
         T: Number,
@@ -1501,7 +1501,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.sub_mut(other);
         result
     }
-    ///
+    /// multiply two arrays
     fn mul(&self, other: &dyn Array<T, (usize, usize)>) -> Self
     where
         T: Number,
@@ -1510,7 +1510,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.mul_mut(other);
         result
     }
-    ///
+    /// divide two arrays
     fn div(&self, other: &dyn Array<T, (usize, usize)>) -> Self
     where
         T: Number,
@@ -1519,7 +1519,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.div_mut(other);
         result
     }
-    ///
+    /// absolute values of the array
     fn abs(&self) -> Self
     where
         T: Number + Signed,
@@ -1528,7 +1528,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.abs_mut();
         result
     }
-    ///
+    /// negation of the array
     fn neg(&self) -> Self
     where
         T: Number + Neg<Output = T>,
@@ -1537,7 +1537,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         result.neg_mut();
         result
     }
-    ///
+    /// values at power `p`
     fn pow(&self, p: T) -> Self
     where
         T: RealNumber,
@@ -1575,7 +1575,7 @@ pub trait Array2<T: Debug + Display + Copy + Sized>: MutArrayView2<T> + Sized +
         }
     }
 
-    /// appriximate equality of the elements of a matrix according to a given error
+    /// approximate equality of the elements of a matrix according to a given error
     fn approximate_eq(&self, other: &Self, error: T) -> bool
     where
         T: Number + RealNumber,
@@ -1631,8 +1631,8 @@ mod tests {
         let v = vec![3., -2., 6.];
         assert_eq!(v.norm(1.), 11.);
         assert_eq!(v.norm(2.), 7.);
-        assert_eq!(v.norm(std::f64::INFINITY), 6.);
-        assert_eq!(v.norm(std::f64::NEG_INFINITY), 2.);
+        assert_eq!(v.norm(f64::INFINITY), 6.);
+        assert_eq!(v.norm(f64::NEG_INFINITY), 2.);
     }
 
     #[test]

src/linalg/basic/matrix.rs

@@ -91,7 +91,7 @@ impl<'a, T: Debug + Display + Copy + Sized> DenseMatrixView<'a, T> {
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> fmt::Display for DenseMatrixView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> fmt::Display for DenseMatrixView<'_, T> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         writeln!(
             f,
@@ -142,7 +142,7 @@ impl<'a, T: Debug + Display + Copy + Sized> DenseMatrixMutView<'a, T> {
         }
     }
 
-    fn iter_mut<'b>(&'b mut self, axis: u8) -> Box<dyn Iterator<Item = &mut T> + 'b> {
+    fn iter_mut<'b>(&'b mut self, axis: u8) -> Box<dyn Iterator<Item = &'b mut T> + 'b> {
         let column_major = self.column_major;
         let stride = self.stride;
         let ptr = self.values.as_mut_ptr();
@@ -169,7 +169,7 @@ impl<'a, T: Debug + Display + Copy + Sized> DenseMatrixMutView<'a, T> {
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> fmt::Display for DenseMatrixMutView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> fmt::Display for DenseMatrixMutView<'_, T> {
     fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
         writeln!(
             f,
@@ -493,7 +493,7 @@ impl<T: Number + RealNumber> EVDDecomposable<T> for DenseMatrix<T> {}
 impl<T: Number + RealNumber> LUDecomposable<T> for DenseMatrix<T> {}
 impl<T: Number + RealNumber> SVDDecomposable<T> for DenseMatrix<T> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> Array<T, (usize, usize)> for DenseMatrixView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> Array<T, (usize, usize)> for DenseMatrixView<'_, T> {
     fn get(&self, pos: (usize, usize)) -> &T {
         if self.column_major {
             &self.values[pos.0 + pos.1 * self.stride]
@@ -515,7 +515,7 @@ impl<'a, T: Debug + Display + Copy + Sized> Array<T, (usize, usize)> for DenseMa
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> Array<T, usize> for DenseMatrixView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> Array<T, usize> for DenseMatrixView<'_, T> {
     fn get(&self, i: usize) -> &T {
         if self.nrows == 1 {
             if self.column_major {
@@ -553,11 +553,11 @@ impl<'a, T: Debug + Display + Copy + Sized> Array<T, usize> for DenseMatrixView<
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView2<T> for DenseMatrixView<'a, T> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView2<T> for DenseMatrixView<'_, T> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView1<T> for DenseMatrixView<'a, T> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView1<T> for DenseMatrixView<'_, T> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> Array<T, (usize, usize)> for DenseMatrixMutView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> Array<T, (usize, usize)> for DenseMatrixMutView<'_, T> {
     fn get(&self, pos: (usize, usize)) -> &T {
         if self.column_major {
             &self.values[pos.0 + pos.1 * self.stride]
@@ -579,9 +579,7 @@ impl<'a, T: Debug + Display + Copy + Sized> Array<T, (usize, usize)> for DenseMa
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, (usize, usize)>
-    for DenseMatrixMutView<'a, T>
-{
+impl<T: Debug + Display + Copy + Sized> MutArray<T, (usize, usize)> for DenseMatrixMutView<'_, T> {
     fn set(&mut self, pos: (usize, usize), x: T) {
         if self.column_major {
             self.values[pos.0 + pos.1 * self.stride] = x;
@@ -595,15 +593,16 @@ impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, (usize, usize)>
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> MutArrayView2<T> for DenseMatrixMutView<'a, T> {}
+impl<T: Debug + Display + Copy + Sized> MutArrayView2<T> for DenseMatrixMutView<'_, T> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView2<T> for DenseMatrixMutView<'a, T> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView2<T> for DenseMatrixMutView<'_, T> {}
 
 impl<T: RealNumber> MatrixStats<T> for DenseMatrix<T> {}
 
 impl<T: RealNumber> MatrixPreprocessing<T> for DenseMatrix<T> {}
 
 #[cfg(test)]
+#[warn(clippy::reversed_empty_ranges)]
 mod tests {
     use super::*;
     use approx::relative_eq;

src/linalg/basic/vector.rs

@@ -119,7 +119,7 @@ impl<T: Debug + Display + Copy + Sized> Array1<T> for Vec<T> {
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> Array<T, usize> for VecMutView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> Array<T, usize> for VecMutView<'_, T> {
     fn get(&self, i: usize) -> &T {
         &self.ptr[i]
     }
@@ -138,7 +138,7 @@ impl<'a, T: Debug + Display + Copy + Sized> Array<T, usize> for VecMutView<'a, T
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, usize> for VecMutView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> MutArray<T, usize> for VecMutView<'_, T> {
     fn set(&mut self, i: usize, x: T) {
         self.ptr[i] = x;
     }
@@ -149,10 +149,10 @@ impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, usize> for VecMutView<'a
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView1<T> for VecMutView<'a, T> {}
-impl<'a, T: Debug + Display + Copy + Sized> MutArrayView1<T> for VecMutView<'a, T> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView1<T> for VecMutView<'_, T> {}
+impl<T: Debug + Display + Copy + Sized> MutArrayView1<T> for VecMutView<'_, T> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> Array<T, usize> for VecView<'a, T> {
+impl<T: Debug + Display + Copy + Sized> Array<T, usize> for VecView<'_, T> {
     fn get(&self, i: usize) -> &T {
         &self.ptr[i]
     }
@@ -171,7 +171,7 @@ impl<'a, T: Debug + Display + Copy + Sized> Array<T, usize> for VecView<'a, T> {
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView1<T> for VecView<'a, T> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView1<T> for VecView<'_, T> {}
 
 #[cfg(test)]
 mod tests {

src/linalg/ndarray/matrix.rs

@@ -68,7 +68,7 @@ impl<T: Debug + Display + Copy + Sized> ArrayView2<T> for ArrayBase<OwnedRepr<T>
 
 impl<T: Debug + Display + Copy + Sized> MutArrayView2<T> for ArrayBase<OwnedRepr<T>, Ix2> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, (usize, usize)> for ArrayView<'a, T, Ix2> {
+impl<T: Debug + Display + Copy + Sized> BaseArray<T, (usize, usize)> for ArrayView<'_, T, Ix2> {
     fn get(&self, pos: (usize, usize)) -> &T {
         &self[[pos.0, pos.1]]
     }
@@ -144,11 +144,9 @@ impl<T: Number + RealNumber> EVDDecomposable<T> for ArrayBase<OwnedRepr<T>, Ix2>
 impl<T: Number + RealNumber> LUDecomposable<T> for ArrayBase<OwnedRepr<T>, Ix2> {}
 impl<T: Number + RealNumber> SVDDecomposable<T> for ArrayBase<OwnedRepr<T>, Ix2> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView2<T> for ArrayView<'a, T, Ix2> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView2<T> for ArrayView<'_, T, Ix2> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, (usize, usize)>
-    for ArrayViewMut<'a, T, Ix2>
-{
+impl<T: Debug + Display + Copy + Sized> BaseArray<T, (usize, usize)> for ArrayViewMut<'_, T, Ix2> {
     fn get(&self, pos: (usize, usize)) -> &T {
         &self[[pos.0, pos.1]]
     }
@@ -175,9 +173,7 @@ impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, (usize, usize)>
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, (usize, usize)>
-    for ArrayViewMut<'a, T, Ix2>
-{
+impl<T: Debug + Display + Copy + Sized> MutArray<T, (usize, usize)> for ArrayViewMut<'_, T, Ix2> {
     fn set(&mut self, pos: (usize, usize), x: T) {
         self[[pos.0, pos.1]] = x
     }
@@ -195,9 +191,9 @@ impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, (usize, usize)>
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> MutArrayView2<T> for ArrayViewMut<'a, T, Ix2> {}
+impl<T: Debug + Display + Copy + Sized> MutArrayView2<T> for ArrayViewMut<'_, T, Ix2> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView2<T> for ArrayViewMut<'a, T, Ix2> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView2<T> for ArrayViewMut<'_, T, Ix2> {}
 
 #[cfg(test)]
 mod tests {

src/linalg/ndarray/vector.rs

@@ -41,7 +41,7 @@ impl<T: Debug + Display + Copy + Sized> ArrayView1<T> for ArrayBase<OwnedRepr<T>
 
 impl<T: Debug + Display + Copy + Sized> MutArrayView1<T> for ArrayBase<OwnedRepr<T>, Ix1> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, usize> for ArrayView<'a, T, Ix1> {
+impl<T: Debug + Display + Copy + Sized> BaseArray<T, usize> for ArrayView<'_, T, Ix1> {
     fn get(&self, i: usize) -> &T {
         &self[i]
     }
@@ -60,9 +60,9 @@ impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, usize> for ArrayView<'a
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView1<T> for ArrayView<'a, T, Ix1> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView1<T> for ArrayView<'_, T, Ix1> {}
 
-impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, usize> for ArrayViewMut<'a, T, Ix1> {
+impl<T: Debug + Display + Copy + Sized> BaseArray<T, usize> for ArrayViewMut<'_, T, Ix1> {
     fn get(&self, i: usize) -> &T {
         &self[i]
     }
@@ -81,7 +81,7 @@ impl<'a, T: Debug + Display + Copy + Sized> BaseArray<T, usize> for ArrayViewMut
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, usize> for ArrayViewMut<'a, T, Ix1> {
+impl<T: Debug + Display + Copy + Sized> MutArray<T, usize> for ArrayViewMut<'_, T, Ix1> {
     fn set(&mut self, i: usize, x: T) {
         self[i] = x;
     }
@@ -92,8 +92,8 @@ impl<'a, T: Debug + Display + Copy + Sized> MutArray<T, usize> for ArrayViewMut<
     }
 }
 
-impl<'a, T: Debug + Display + Copy + Sized> ArrayView1<T> for ArrayViewMut<'a, T, Ix1> {}
-impl<'a, T: Debug + Display + Copy + Sized> MutArrayView1<T> for ArrayViewMut<'a, T, Ix1> {}
+impl<T: Debug + Display + Copy + Sized> ArrayView1<T> for ArrayViewMut<'_, T, Ix1> {}
+impl<T: Debug + Display + Copy + Sized> MutArrayView1<T> for ArrayViewMut<'_, T, Ix1> {}
 
 impl<T: Debug + Display + Copy + Sized> Array1<T> for ArrayBase<OwnedRepr<T>, Ix1> {
     fn slice<'a>(&'a self, range: Range<usize>) -> Box<dyn ArrayView1<T> + 'a> {

src/linalg/traits/evd.rs

@@ -841,7 +841,7 @@ mod tests {
         ));
         for (i, eigen_values_i) in eigen_values.iter().enumerate() {
             assert!((eigen_values_i - evd.d[i]).abs() < 1e-4);
-            assert!((0f64 - evd.e[i]).abs() < std::f64::EPSILON);
+            assert!((0f64 - evd.e[i]).abs() < f64::EPSILON);
         }
     }
     #[cfg_attr(
@@ -875,7 +875,7 @@ mod tests {
         ));
         for (i, eigen_values_i) in eigen_values.iter().enumerate() {
             assert!((eigen_values_i - evd.d[i]).abs() < 1e-4);
-            assert!((0f64 - evd.e[i]).abs() < std::f64::EPSILON);
+            assert!((0f64 - evd.e[i]).abs() < f64::EPSILON);
         }
     }
     #[cfg_attr(

src/linalg/traits/stats.rs

@@ -142,7 +142,6 @@ pub trait MatrixPreprocessing<T: RealNumber>: MutArrayView2<T> + Clone {
     ///
     /// assert_eq!(a, expected);
     /// ```
-
     fn binarize_mut(&mut self, threshold: T) {
         let (nrows, ncols) = self.shape();
         for row in 0..nrows {
@@ -217,8 +216,8 @@ mod tests {
         let expected_0 = vec![0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
         let expected_1 = vec![1.25, 1.25];
 
-        assert!(m.var(0).approximate_eq(&expected_0, std::f64::EPSILON));
-        assert!(m.var(1).approximate_eq(&expected_1, std::f64::EPSILON));
+        assert!(m.var(0).approximate_eq(&expected_0, f64::EPSILON));
+        assert!(m.var(1).approximate_eq(&expected_1, f64::EPSILON));
         assert_eq!(
             m.mean(0),
             vec![0.0, 0.25, 0.25, 1.25, 1.5, 1.75, 2.75, 3.25]

src/linalg/traits/svd.rs

@@ -48,11 +48,9 @@ pub struct SVD<T: Number + RealNumber, M: SVDDecomposable<T>> {
     pub V: M,
     /// Singular values of the original matrix
     pub s: Vec<T>,
-    ///
     m: usize,
-    ///
     n: usize,
-    ///
+    /// Tolerance
     tol: T,
 }
 

src/linear/bg_solver.rs

@@ -27,9 +27,9 @@ use crate::error::Failed;
 use crate::linalg::basic::arrays::{Array, Array1, Array2, ArrayView1, MutArrayView1};
 use crate::numbers::floatnum::FloatNumber;
 
-///
+/// Trait for Biconjugate Gradient Solver
 pub trait BiconjugateGradientSolver<'a, T: FloatNumber, X: Array2<T>> {
-    ///
+    /// Solve Ax = b
     fn solve_mut(
         &self,
         a: &'a X,
@@ -109,7 +109,7 @@ pub trait BiconjugateGradientSolver<'a, T: FloatNumber, X: Array2<T>> {
         Ok(err)
     }
 
-    ///
+    /// solve preconditioner
     fn solve_preconditioner(&self, a: &'a X, b: &[T], x: &mut [T]) {
         let diag = Self::diag(a);
         let n = diag.len();
@@ -133,7 +133,7 @@ pub trait BiconjugateGradientSolver<'a, T: FloatNumber, X: Array2<T>> {
         y.copy_from(&x.xa(true, a));
     }
 
-    ///
+    /// Extract the diagonal from a matrix
     fn diag(a: &X) -> Vec<T> {
         let (nrows, ncols) = a.shape();
         let n = nrows.min(ncols);

src/linear/lasso_optimizer.rs

@@ -16,7 +16,7 @@ use crate::linalg::basic::arrays::{Array1, Array2, ArrayView1, MutArray, MutArra
 use crate::linear::bg_solver::BiconjugateGradientSolver;
 use crate::numbers::floatnum::FloatNumber;
 
-///
+/// Interior Point Optimizer
 pub struct InteriorPointOptimizer<T: FloatNumber, X: Array2<T>> {
     ata: X,
     d1: Vec<T>,
@@ -25,9 +25,8 @@ pub struct InteriorPointOptimizer<T: FloatNumber, X: Array2<T>> {
     prs: Vec<T>,
 }
 
-///
 impl<T: FloatNumber, X: Array2<T>> InteriorPointOptimizer<T, X> {
-    ///
+    /// Initialize a new Interior Point Optimizer
     pub fn new(a: &X, n: usize) -> InteriorPointOptimizer<T, X> {
         InteriorPointOptimizer {
             ata: a.ab(true, a, false),
@@ -38,7 +37,7 @@ impl<T: FloatNumber, X: Array2<T>> InteriorPointOptimizer<T, X> {
         }
     }
 
-    ///
+    /// Run the optimization
     pub fn optimize(
         &mut self,
         x: &X,
@@ -101,7 +100,7 @@ impl<T: FloatNumber, X: Array2<T>> InteriorPointOptimizer<T, X> {
 
             // CALCULATE DUALITY GAP
             let xnu = nu.xa(false, x);
-            let max_xnu = xnu.norm(std::f64::INFINITY);
+            let max_xnu = xnu.norm(f64::INFINITY);
             if max_xnu > lambda_f64 {
                 let lnu = T::from_f64(lambda_f64 / max_xnu).unwrap();
                 nu.mul_scalar_mut(lnu);
@@ -208,7 +207,6 @@ impl<T: FloatNumber, X: Array2<T>> InteriorPointOptimizer<T, X> {
         Ok(w)
     }
 
-    ///
     fn sumlogneg(f: &X) -> T {
         let (n, _) = f.shape();
         let mut sum = T::zero();
@@ -220,11 +218,9 @@ impl<T: FloatNumber, X: Array2<T>> InteriorPointOptimizer<T, X> {
     }
 }
 
-///
 impl<'a, T: FloatNumber, X: Array2<T>> BiconjugateGradientSolver<'a, T, X>
     for InteriorPointOptimizer<T, X>
 {
-    ///
     fn solve_preconditioner(&self, a: &'a X, b: &[T], x: &mut [T]) {
         let (_, p) = a.shape();
 
@@ -234,7 +230,6 @@ impl<'a, T: FloatNumber, X: Array2<T>> BiconjugateGradientSolver<'a, T, X>
         }
     }
 
-    ///
     fn mat_vec_mul(&self, _: &X, x: &Vec<T>, y: &mut Vec<T>) {
         let (_, p) = self.ata.shape();
         let x_slice = Vec::from_slice(x.slice(0..p).as_ref());
@@ -246,7 +241,6 @@ impl<'a, T: FloatNumber, X: Array2<T>> BiconjugateGradientSolver<'a, T, X>
         }
     }
 
-    ///
     fn mat_t_vec_mul(&self, a: &X, x: &Vec<T>, y: &mut Vec<T>) {
         self.mat_vec_mul(a, x, y);
     }

src/linear/logistic_regression.rs

@@ -183,14 +183,11 @@ pub struct LogisticRegression<
 }
 
 trait ObjectiveFunction<T: Number + FloatNumber, X: Array2<T>> {
-    ///
     fn f(&self, w_bias: &[T]) -> T;
 
-    ///
     #[allow(clippy::ptr_arg)]
     fn df(&self, g: &mut Vec<T>, w_bias: &Vec<T>);
 
-    ///
     #[allow(clippy::ptr_arg)]
     fn partial_dot(w: &[T], x: &X, v_col: usize, m_row: usize) -> T {
         let mut sum = T::zero();
@@ -261,8 +258,8 @@ impl<TX: Number + FloatNumber + RealNumber, TY: Number + Ord, X: Array2<TX>, Y:
     }
 }
 
-impl<'a, T: Number + FloatNumber, X: Array2<T>> ObjectiveFunction<T, X>
-    for BinaryObjectiveFunction<'a, T, X>
+impl<T: Number + FloatNumber, X: Array2<T>> ObjectiveFunction<T, X>
+    for BinaryObjectiveFunction<'_, T, X>
 {
     fn f(&self, w_bias: &[T]) -> T {
         let mut f = T::zero();
@@ -316,8 +313,8 @@ struct MultiClassObjectiveFunction<'a, T: Number + FloatNumber, X: Array2<T>> {
     _phantom_t: PhantomData<T>,
 }
 
-impl<'a, T: Number + FloatNumber + RealNumber, X: Array2<T>> ObjectiveFunction<T, X>
-    for MultiClassObjectiveFunction<'a, T, X>
+impl<T: Number + FloatNumber + RealNumber, X: Array2<T>> ObjectiveFunction<T, X>
+    for MultiClassObjectiveFunction<'_, T, X>
 {
     fn f(&self, w_bias: &[T]) -> T {
         let mut f = T::zero();
@@ -629,11 +626,11 @@ mod tests {
         objective.df(&mut g, &vec![1., 2., 3., 4., 5., 6., 7., 8., 9.]);
         objective.df(&mut g, &vec![1., 2., 3., 4., 5., 6., 7., 8., 9.]);
 
-        assert!((g[0] + 33.000068218163484).abs() < std::f64::EPSILON);
+        assert!((g[0] + 33.000068218163484).abs() < f64::EPSILON);
 
         let f = objective.f(&[1., 2., 3., 4., 5., 6., 7., 8., 9.]);
 
-        assert!((f - 408.0052230582765).abs() < std::f64::EPSILON);
+        assert!((f - 408.0052230582765).abs() < f64::EPSILON);
 
         let objective_reg = MultiClassObjectiveFunction {
             x: &x,
@@ -689,13 +686,13 @@ mod tests {
         objective.df(&mut g, &vec![1., 2., 3.]);
         objective.df(&mut g, &vec![1., 2., 3.]);
 
-        assert!((g[0] - 26.051064349381285).abs() < std::f64::EPSILON);
-        assert!((g[1] - 10.239000702928523).abs() < std::f64::EPSILON);
-        assert!((g[2] - 3.869294270156324).abs() < std::f64::EPSILON);
+        assert!((g[0] - 26.051064349381285).abs() < f64::EPSILON);
+        assert!((g[1] - 10.239000702928523).abs() < f64::EPSILON);
+        assert!((g[2] - 3.869294270156324).abs() < f64::EPSILON);
 
         let f = objective.f(&[1., 2., 3.]);
 
-        assert!((f - 59.76994756647412).abs() < std::f64::EPSILON);
+        assert!((f - 59.76994756647412).abs() < f64::EPSILON);
 
         let objective_reg = BinaryObjectiveFunction {
             x: &x,
@@ -916,7 +913,7 @@ mod tests {
         let x: DenseMatrix<f32> = DenseMatrix::rand(52181, 94);
         let y1: Vec<i32> = vec![1; 2181];
         let y2: Vec<i32> = vec![0; 50000];
-        let y: Vec<i32> = y1.into_iter().chain(y2.into_iter()).collect();
+        let y: Vec<i32> = y1.into_iter().chain(y2).collect();
 
         let lr = LogisticRegression::fit(&x, &y, Default::default()).unwrap();
         let lr_reg = LogisticRegression::fit(
@@ -938,12 +935,12 @@ mod tests {
         let x: &DenseMatrix<f64> = &DenseMatrix::rand(52181, 94);
         let y1: Vec<u32> = vec![1; 2181];
         let y2: Vec<u32> = vec![0; 50000];
-        let y: &Vec<u32> = &(y1.into_iter().chain(y2.into_iter()).collect());
+        let y: &Vec<u32> = &(y1.into_iter().chain(y2).collect());
         println!("y vec height: {:?}", y.len());
         println!("x matrix shape: {:?}", x.shape());
 
         let lr = LogisticRegression::fit(x, y, Default::default()).unwrap();
-        let y_hat = lr.predict(&x).unwrap();
+        let y_hat = lr.predict(x).unwrap();
 
         println!("y_hat shape: {:?}", y_hat.shape());
 

src/naive_bayes/bernoulli.rs

@@ -427,6 +427,7 @@ impl<TX: Number + PartialOrd, TY: Number + Ord + Unsigned, X: Array2<TX>, Y: Arr
 
     /// Estimates the class labels for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with class estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         if let Some(threshold) = self.binarize {

src/naive_bayes/categorical.rs

@@ -95,7 +95,7 @@ impl<T: Number + Unsigned> PartialEq for CategoricalNBDistribution<T> {
                         return false;
                     }
                     for (a_i_j, b_i_j) in a_i.iter().zip(b_i.iter()) {
-                        if (*a_i_j - *b_i_j).abs() > std::f64::EPSILON {
+                        if (*a_i_j - *b_i_j).abs() > f64::EPSILON {
                             return false;
                         }
                     }
@@ -375,6 +375,7 @@ impl<T: Number + Unsigned, X: Array2<T>, Y: Array1<T>> CategoricalNB<T, X, Y> {
 
     /// Estimates the class labels for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with class estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         self.inner.as_ref().unwrap().predict(x)

src/naive_bayes/gaussian.rs

@@ -328,6 +328,7 @@ impl<TX: Number + RealNumber, TY: Number + Ord + Unsigned, X: Array2<TX>, Y: Arr
 
     /// Estimates the class labels for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with class estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         self.inner.as_ref().unwrap().predict(x)

src/naive_bayes/mod.rs

@@ -40,7 +40,7 @@ use crate::linalg::basic::arrays::{Array1, Array2, ArrayView1};
 use crate::numbers::basenum::Number;
 #[cfg(feature = "serde")]
 use serde::{Deserialize, Serialize};
-use std::{cmp::Ordering, marker::PhantomData};
+use std::marker::PhantomData;
 
 /// Distribution used in the Naive Bayes classifier.
 pub(crate) trait NBDistribution<X: Number, Y: Number>: Clone {
@@ -89,45 +89,46 @@ impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>, D: NBDistribution<TX,
 
     /// Estimates the class labels for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with class estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         let y_classes = self.distribution.classes();
-        let predictions = x
-            .row_iter()
-            .map(|row| {
-                y_classes
-                    .iter()
-                    .enumerate()
-                    .map(|(class_index, class)| {
-                        (
-                            class,
-                            self.distribution.log_likelihood(class_index, &row)
-                                + self.distribution.prior(class_index).ln(),
-                        )
-                    })
-                    // For some reason, the max_by method cannot use NaNs for finding the maximum value, it panics.
-                    // NaN must be considered as minimum values,
-                    // therefore it's like NaNs would not be considered for choosing the maximum value.
-                    // So we need to handle this case for avoiding panicking by using `Option::unwrap`.
-                    .max_by(|(_, p1), (_, p2)| match p1.partial_cmp(p2) {
-                        Some(ordering) => ordering,
-                        None => {
-                            if p1.is_nan() {
-                                Ordering::Less
-                            } else if p2.is_nan() {
-                                Ordering::Greater
+
+        if y_classes.is_empty() {
+            return Err(Failed::predict("Failed to predict, no classes available"));
+        }
+
+        let (rows, _) = x.shape();
+        let mut predictions = Vec::with_capacity(rows);
+        let mut all_probs_nan = true;
+
+        for row_index in 0..rows {
+            let row = x.get_row(row_index);
+            let mut max_log_prob = f64::NEG_INFINITY;
+            let mut max_class = None;
+
+            for (class_index, class) in y_classes.iter().enumerate() {
+                let log_likelihood = self.distribution.log_likelihood(class_index, &row);
+                let log_prob = log_likelihood + self.distribution.prior(class_index).ln();
+
+                if !log_prob.is_nan() && log_prob > max_log_prob {
+                    max_log_prob = log_prob;
+                    max_class = Some(*class);
+                    all_probs_nan = false;
+                }
+            }
+
+            predictions.push(max_class.unwrap_or(y_classes[0]));
+        }
+
+        if all_probs_nan {
+            Err(Failed::predict(
+                "Failed to predict, all probabilities were NaN",
+            ))
         } else {
-                                Ordering::Equal
+            Ok(Y::from_vec_slice(&predictions))
         }
     }
-                    })
-                    .map(|(prediction, _probability)| *prediction)
-                    .ok_or_else(|| Failed::predict("Failed to predict, there is no result"))
-            })
-            .collect::<Result<Vec<TY>, Failed>>()?;
-        let y_hat = Y::from_vec_slice(&predictions);
-        Ok(y_hat)
-    }
 }
 pub mod bernoulli;
 pub mod categorical;
@@ -146,7 +147,7 @@ mod tests {
     #[derive(Debug, PartialEq, Clone)]
     struct TestDistribution<'d>(&'d Vec<i32>);
 
-    impl<'d> NBDistribution<i32, i32> for TestDistribution<'d> {
+    impl NBDistribution<i32, i32> for TestDistribution<'_> {
         fn prior(&self, _class_index: usize) -> f64 {
             1.
         }
@@ -163,7 +164,7 @@ mod tests {
         }
 
         fn classes(&self) -> &Vec<i32> {
-            &self.0
+            self.0
         }
     }
 
@@ -176,7 +177,7 @@ mod tests {
             Ok(_) => panic!("Should return error in case of empty classes"),
             Err(err) => assert_eq!(
                 err.to_string(),
-                "Predict failed: Failed to predict, there is no result"
+                "Predict failed: Failed to predict, no classes available"
             ),
         }
 
@@ -192,4 +193,441 @@ mod tests {
             Err(_) => panic!("Should success in normal case without NaNs"),
         }
     }
+
+    // A simple test distribution using float
+    #[derive(Debug, PartialEq, Clone)]
+    struct TestDistributionAgain {
+        classes: Vec<u32>,
+        probs: Vec<f64>,
+    }
+
+    impl NBDistribution<f64, u32> for TestDistributionAgain {
+        fn classes(&self) -> &Vec<u32> {
+            &self.classes
+        }
+        fn prior(&self, class_index: usize) -> f64 {
+            self.probs[class_index]
+        }
+        fn log_likelihood<'a>(
+            &'a self,
+            class_index: usize,
+            _j: &'a Box<dyn ArrayView1<f64> + 'a>,
+        ) -> f64 {
+            self.probs[class_index].ln()
+        }
+    }
+
+    type TestNB = BaseNaiveBayes<f64, u32, DenseMatrix<f64>, Vec<u32>, TestDistributionAgain>;
+
+    #[test]
+    fn test_predict_empty_classes() {
+        let dist = TestDistributionAgain {
+            classes: vec![],
+            probs: vec![],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        assert!(nb.predict(&x).is_err());
+    }
+
+    #[test]
+    fn test_predict_single_class() {
+        let dist = TestDistributionAgain {
+            classes: vec![1],
+            probs: vec![1.0],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        let result = nb.predict(&x).unwrap();
+        assert_eq!(result, vec![1, 1]);
+    }
+
+    #[test]
+    fn test_predict_multiple_classes() {
+        let dist = TestDistributionAgain {
+            classes: vec![1, 2, 3],
+            probs: vec![0.2, 0.5, 0.3],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0], &[5.0, 6.0]]).unwrap();
+        let result = nb.predict(&x).unwrap();
+        assert_eq!(result, vec![2, 2, 2]);
+    }
+
+    #[test]
+    fn test_predict_with_nans() {
+        let dist = TestDistributionAgain {
+            classes: vec![1, 2],
+            probs: vec![f64::NAN, 0.5],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        let result = nb.predict(&x).unwrap();
+        assert_eq!(result, vec![2, 2]);
+    }
+
+    #[test]
+    fn test_predict_all_nans() {
+        let dist = TestDistributionAgain {
+            classes: vec![1, 2],
+            probs: vec![f64::NAN, f64::NAN],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        assert!(nb.predict(&x).is_err());
+    }
+
+    #[test]
+    fn test_predict_extreme_probabilities() {
+        let dist = TestDistributionAgain {
+            classes: vec![1, 2],
+            probs: vec![1e-300, 1e-301],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        let result = nb.predict(&x).unwrap();
+        assert_eq!(result, vec![1, 1]);
+    }
+
+    #[test]
+    fn test_predict_with_infinity() {
+        let dist = TestDistributionAgain {
+            classes: vec![1, 2, 3],
+            probs: vec![f64::INFINITY, 1.0, 2.0],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        let result = nb.predict(&x).unwrap();
+        assert_eq!(result, vec![1, 1]);
+    }
+
+    #[test]
+    fn test_predict_with_negative_infinity() {
+        let dist = TestDistributionAgain {
+            classes: vec![1, 2, 3],
+            probs: vec![f64::NEG_INFINITY, 1.0, 2.0],
+        };
+        let nb = TestNB::fit(dist).unwrap();
+        let x = DenseMatrix::from_2d_array(&[&[1.0, 2.0], &[3.0, 4.0]]).unwrap();
+        let result = nb.predict(&x).unwrap();
+        assert_eq!(result, vec![3, 3]);
+    }
+
+    #[test]
+    fn test_gaussian_naive_bayes_numerical_stability() {
+        #[derive(Debug, PartialEq, Clone)]
+        struct GaussianTestDistribution {
+            classes: Vec<u32>,
+            means: Vec<Vec<f64>>,
+            variances: Vec<Vec<f64>>,
+            priors: Vec<f64>,
+        }
+
+        impl NBDistribution<f64, u32> for GaussianTestDistribution {
+            fn classes(&self) -> &Vec<u32> {
+                &self.classes
+            }
+
+            fn prior(&self, class_index: usize) -> f64 {
+                self.priors[class_index]
+            }
+
+            fn log_likelihood<'a>(
+                &'a self,
+                class_index: usize,
+                j: &'a Box<dyn ArrayView1<f64> + 'a>,
+            ) -> f64 {
+                let means = &self.means[class_index];
+                let variances = &self.variances[class_index];
+                j.iterator(0)
+                    .enumerate()
+                    .map(|(i, &xi)| {
+                        let mean = means[i];
+                        let var = variances[i] + 1e-9; // Small smoothing for numerical stability
+                        let coeff = -0.5 * (2.0 * std::f64::consts::PI * var).ln();
+                        let exponent = -(xi - mean).powi(2) / (2.0 * var);
+                        coeff + exponent
+                    })
+                    .sum()
+            }
+        }
+
+        fn train_distribution(x: &DenseMatrix<f64>, y: &[u32]) -> GaussianTestDistribution {
+            let mut classes: Vec<u32> = y
+                .iter()
+                .cloned()
+                .collect::<std::collections::HashSet<u32>>()
+                .into_iter()
+                .collect();
+            classes.sort();
+            let n_classes = classes.len();
+            let n_features = x.shape().1;
+
+            let mut means = vec![vec![0.0; n_features]; n_classes];
+            let mut variances = vec![vec![0.0; n_features]; n_classes];
+            let mut class_counts = vec![0; n_classes];
+
+            // Calculate means and count samples per class
+            for (sample, &class) in x.row_iter().zip(y.iter()) {
+                let class_idx = classes.iter().position(|&c| c == class).unwrap();
+                class_counts[class_idx] += 1;
+                for (i, &value) in sample.iterator(0).enumerate() {
+                    means[class_idx][i] += value;
+                }
+            }
+
+            // Normalize means
+            for (class_idx, mean) in means.iter_mut().enumerate() {
+                for value in mean.iter_mut() {
+                    *value /= class_counts[class_idx] as f64;
+                }
+            }
+
+            // Calculate variances
+            for (sample, &class) in x.row_iter().zip(y.iter()) {
+                let class_idx = classes.iter().position(|&c| c == class).unwrap();
+                for (i, &value) in sample.iterator(0).enumerate() {
+                    let diff = value - means[class_idx][i];
+                    variances[class_idx][i] += diff * diff;
+                }
+            }
+
+            // Normalize variances and add small epsilon to avoid zero variance
+            let epsilon = 1e-9;
+            for (class_idx, variance) in variances.iter_mut().enumerate() {
+                for value in variance.iter_mut() {
+                    *value = *value / class_counts[class_idx] as f64 + epsilon;
+                }
+            }
+
+            // Calculate priors
+            let total_samples = y.len() as f64;
+            let priors: Vec<f64> = class_counts
+                .iter()
+                .map(|&count| count as f64 / total_samples)
+                .collect();
+
+            GaussianTestDistribution {
+                classes,
+                means,
+                variances,
+                priors,
+            }
+        }
+
+        type TestNBGaussian =
+            BaseNaiveBayes<f64, u32, DenseMatrix<f64>, Vec<u32>, GaussianTestDistribution>;
+
+        // Create a constant training dataset
+        let n_samples = 1000;
+        let n_features = 5;
+        let n_classes = 4;
+
+        let mut x_data = Vec::with_capacity(n_samples * n_features);
+        let mut y_data = Vec::with_capacity(n_samples);
+
+        for i in 0..n_samples {
+            for j in 0..n_features {
+                x_data.push((i * j) as f64 % 10.0);
+            }
+            y_data.push((i % n_classes) as u32);
+        }
+
+        let x = DenseMatrix::new(n_samples, n_features, x_data, true).unwrap();
+        let y = y_data;
+
+        // Train the model
+        let dist = train_distribution(&x, &y);
+        let nb = TestNBGaussian::fit(dist).unwrap();
+
+        // Create constant test data
+        let n_test_samples = 100;
+        let mut test_x_data = Vec::with_capacity(n_test_samples * n_features);
+        for i in 0..n_test_samples {
+            for j in 0..n_features {
+                test_x_data.push((i * j * 2) as f64 % 15.0);
+            }
+        }
+        let test_x = DenseMatrix::new(n_test_samples, n_features, test_x_data, true).unwrap();
+
+        // Make predictions
+        let predictions = nb
+            .predict(&test_x)
+            .map_err(|e| format!("Prediction failed: {}", e))
+            .unwrap();
+
+        // Check numerical stability
+        assert_eq!(
+            predictions.len(),
+            n_test_samples,
+            "Number of predictions should match number of test samples"
+        );
+
+        // Check that all predictions are valid class labels
+        for &pred in predictions.iter() {
+            assert!(pred < n_classes as u32, "Predicted class should be valid");
+        }
+
+        // Check consistency of predictions
+        let repeated_predictions = nb
+            .predict(&test_x)
+            .map_err(|e| format!("Repeated prediction failed: {}", e))
+            .unwrap();
+        assert_eq!(
+            predictions, repeated_predictions,
+            "Predictions should be consistent when repeated"
+        );
+
+        // Check extreme values
+        let extreme_x =
+            DenseMatrix::new(2, n_features, vec![f64::MAX; n_features * 2], true).unwrap();
+        let extreme_predictions = nb.predict(&extreme_x);
+        assert!(
+            extreme_predictions.is_err(),
+            "Extreme value input should result in an error"
+        );
+        assert_eq!(
+            extreme_predictions.unwrap_err().to_string(),
+            "Predict failed: Failed to predict, all probabilities were NaN",
+            "Incorrect error message for extreme values"
+        );
+
+        // Check for NaN handling
+        let nan_x = DenseMatrix::new(2, n_features, vec![f64::NAN; n_features * 2], true).unwrap();
+        let nan_predictions = nb.predict(&nan_x);
+        assert!(
+            nan_predictions.is_err(),
+            "NaN input should result in an error"
+        );
+
+        // Check for very small values
+        let small_x =
+            DenseMatrix::new(2, n_features, vec![f64::MIN_POSITIVE; n_features * 2], true).unwrap();
+        let small_predictions = nb
+            .predict(&small_x)
+            .map_err(|e| format!("Small value prediction failed: {}", e))
+            .unwrap();
+        for &pred in small_predictions.iter() {
+            assert!(
+                pred < n_classes as u32,
+                "Predictions for very small values should be valid"
+            );
+        }
+
+        // Check for values close to zero
+        let near_zero_x =
+            DenseMatrix::new(2, n_features, vec![1e-300; n_features * 2], true).unwrap();
+        let near_zero_predictions = nb
+            .predict(&near_zero_x)
+            .map_err(|e| format!("Near-zero value prediction failed: {}", e))
+            .unwrap();
+        for &pred in near_zero_predictions.iter() {
+            assert!(
+                pred < n_classes as u32,
+                "Predictions for near-zero values should be valid"
+            );
+        }
+
+        println!("All numerical stability checks passed!");
+    }
+
+    #[test]
+    fn test_gaussian_naive_bayes_numerical_stability_random_data() {
+        #[derive(Debug)]
+        struct MySimpleRng {
+            state: u64,
+        }
+
+        impl MySimpleRng {
+            fn new(seed: u64) -> Self {
+                MySimpleRng { state: seed }
+            }
+
+            /// Get the next u64 in the sequence.
+            fn next_u64(&mut self) -> u64 {
+                // LCG parameters; these are somewhat arbitrary but commonly used.
+                // Feel free to tweak the multiplier/adder etc.
+                self.state = self.state.wrapping_mul(6364136223846793005).wrapping_add(1);
+                self.state
+            }
+
+            /// Get an f64 in the range [min, max).
+            fn next_f64(&mut self, min: f64, max: f64) -> f64 {
+                let fraction = (self.next_u64() as f64) / (u64::MAX as f64);
+                min + fraction * (max - min)
+            }
+
+            /// Get a usize in the range [min, max). This floors the floating result.
+            fn gen_range_usize(&mut self, min: usize, max: usize) -> usize {
+                let v = self.next_f64(min as f64, max as f64);
+                // Truncate into the integer range. Because of floating inexactness,
+                // ensure we also clamp.
+                let int_v = v.floor() as isize;
+                // simple clamp to avoid any float rounding out of range
+                let clamped = int_v.max(min as isize).min((max - 1) as isize);
+                clamped as usize
+            }
+        }
+        use crate::naive_bayes::gaussian::GaussianNB;
+        // We will generate random data in a reproducible way (using a fixed seed).
+        // We will generate random data in a reproducible way:
+        let mut rng = MySimpleRng::new(42);
+
+        let n_samples = 1000;
+        let n_features = 5;
+        let n_classes = 4;
+
+        // Our feature matrix and label vector
+        let mut x_data = Vec::with_capacity(n_samples * n_features);
+        let mut y_data = Vec::with_capacity(n_samples);
+
+        // Fill x_data with random values and y_data with random class labels.
+        for _i in 0..n_samples {
+            for _j in 0..n_features {
+                // We’ll pick random values in [-10, 10).
+                x_data.push(rng.next_f64(-10.0, 10.0));
+            }
+            let class = rng.gen_range_usize(0, n_classes) as u32;
+            y_data.push(class);
+        }
+
+        // Create DenseMatrix from x_data
+        let x = DenseMatrix::new(n_samples, n_features, x_data, true).unwrap();
+
+        // Train GaussianNB
+        let gnb = GaussianNB::fit(&x, &y_data, Default::default())
+            .expect("Fitting GaussianNB with random data failed.");
+
+        // Predict on the same training data to verify no numerical instability
+        let predictions = gnb.predict(&x).expect("Prediction on random data failed.");
+
+        // Basic sanity checks
+        assert_eq!(
+            predictions.len(),
+            n_samples,
+            "Prediction size must match n_samples"
+        );
+        for &pred_class in &predictions {
+            assert!(
+                (pred_class as usize) < n_classes,
+                "Predicted class {} is out of range [0..n_classes).",
+                pred_class
+            );
+        }
+
+        // If you want to compare with scikit-learn, you can do something like:
+        // println!("X = {:?}", &x);
+        // println!("Y = {:?}", &y_data);
+        // println!("predictions = {:?}", &predictions);
+        // and then in Python:
+        //    import numpy as np
+        //    from sklearn.naive_bayes import GaussianNB
+        //    X = np.reshape(np.array(x), (1000, 5), order='F')
+        //    Y = np.array(y)
+        //    gnb = GaussianNB().fit(X, Y)
+        //    preds = gnb.predict(X)
+        //    expected = np.array(predictions)
+        //    assert expected == preds
+        // They should match closely (or exactly) depending on floating rounding.
+    }
 }

src/naive_bayes/multinomial.rs

@@ -358,6 +358,7 @@ impl<TX: Number + Unsigned, TY: Number + Ord + Unsigned, X: Array2<TX>, Y: Array
 
     /// Estimates the class labels for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with class estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         self.inner.as_ref().unwrap().predict(x)

src/neighbors/knn_classifier.rs

@@ -261,6 +261,7 @@ impl<TX: Number, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>, D: Distance<Vec
 
     /// Estimates the class labels for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with class estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         let mut result = Y::zeros(x.shape().0);

src/neighbors/knn_regressor.rs

@@ -88,25 +88,21 @@ pub struct KNNRegressor<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>, D:
 impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>, D: Distance<Vec<TX>>>
     KNNRegressor<TX, TY, X, Y, D>
 {
-    ///
     fn y(&self) -> &Y {
         self.y.as_ref().unwrap()
     }
 
-    ///
     fn knn_algorithm(&self) -> &KNNAlgorithm<TX, D> {
         self.knn_algorithm
             .as_ref()
             .expect("Missing parameter: KNNAlgorithm")
     }
 
-    ///
     fn weight(&self) -> &KNNWeightFunction {
         self.weight.as_ref().expect("Missing parameter: weight")
     }
 
     #[allow(dead_code)]
-    ///
     fn k(&self) -> usize {
         self.k.unwrap()
     }
@@ -250,6 +246,7 @@ impl<TX: Number, TY: Number, X: Array2<TX>, Y: Array1<TY>, D: Distance<Vec<TX>>>
 
     /// Predict the target for the provided data.
     /// * `x` - data of shape NxM where N is number of data points to estimate and M is number of features.
+    ///
     /// Returns a vector of size N with estimates.
     pub fn predict(&self, x: &X) -> Result<Y, Failed> {
         let mut result = Y::zeros(x.shape().0);
@@ -312,7 +309,7 @@ mod tests {
         let y_hat = knn.predict(&x).unwrap();
         assert_eq!(5, Vec::len(&y_hat));
         for i in 0..y_hat.len() {
-            assert!((y_hat[i] - y_exp[i]).abs() < std::f64::EPSILON);
+            assert!((y_hat[i] - y_exp[i]).abs() < f64::EPSILON);
         }
     }
 

src/optimization/first_order/gradient_descent.rs

@@ -1,5 +1,3 @@
-// TODO: missing documentation
-
 use std::default::Default;
 
 use crate::linalg::basic::arrays::Array1;
@@ -8,30 +6,27 @@ use crate::optimization::first_order::{FirstOrderOptimizer, OptimizerResult};
 use crate::optimization::line_search::LineSearchMethod;
 use crate::optimization::{DF, F};
 
-///
+/// Gradient Descent optimization algorithm
 pub struct GradientDescent {
-    ///
+    /// Maximum number of iterations
     pub max_iter: usize,
-    ///
+    /// Relative tolerance for the gradient norm
     pub g_rtol: f64,
-    ///
+    /// Absolute tolerance for the gradient norm
     pub g_atol: f64,
 }
 
-///
 impl Default for GradientDescent {
     fn default() -> Self {
         GradientDescent {
             max_iter: 10000,
-            g_rtol: std::f64::EPSILON.sqrt(),
-            g_atol: std::f64::EPSILON,
+            g_rtol: f64::EPSILON.sqrt(),
+            g_atol: f64::EPSILON,
         }
     }
 }
 
-///
 impl<T: FloatNumber> FirstOrderOptimizer<T> for GradientDescent {
-    ///
     fn optimize<'a, X: Array1<T>, LS: LineSearchMethod<T>>(
         &self,
         f: &'a F<'_, T, X>,

src/optimization/first_order/lbfgs.rs

@@ -11,31 +11,29 @@ use crate::optimization::first_order::{FirstOrderOptimizer, OptimizerResult};
 use crate::optimization::line_search::LineSearchMethod;
 use crate::optimization::{DF, F};
 
-///
+/// Limited-memory BFGS optimization algorithm
 pub struct LBFGS {
-    ///
+    /// Maximum number of iterations
     pub max_iter: usize,
-    ///
+    /// TODO: Add documentation
     pub g_rtol: f64,
-    ///
+    /// TODO: Add documentation
     pub g_atol: f64,
-    ///
+    /// TODO: Add documentation
     pub x_atol: f64,
-    ///
+    /// TODO: Add documentation
     pub x_rtol: f64,
-    ///
+    /// TODO: Add documentation
     pub f_abstol: f64,
-    ///
+    /// TODO: Add documentation
     pub f_reltol: f64,
-    ///
+    /// TODO: Add documentation
     pub successive_f_tol: usize,
-    ///
+    /// TODO: Add documentation
     pub m: usize,
 }
 
-///
 impl Default for LBFGS {
-    ///
     fn default() -> Self {
         LBFGS {
             max_iter: 1000,
@@ -51,9 +49,7 @@ impl Default for LBFGS {
     }
 }
 
-///
 impl LBFGS {
-    ///
     fn two_loops<T: FloatNumber + RealNumber, X: Array1<T>>(&self, state: &mut LBFGSState<T, X>) {
         let lower = state.iteration.max(self.m) - self.m;
         let upper = state.iteration;
@@ -95,7 +91,6 @@ impl LBFGS {
         state.s.mul_scalar_mut(-T::one());
     }
 
-    ///
     fn init_state<T: FloatNumber + RealNumber, X: Array1<T>>(&self, x: &X) -> LBFGSState<T, X> {
         LBFGSState {
             x: x.clone(),
@@ -119,7 +114,6 @@ impl LBFGS {
         }
     }
 
-    ///
     fn update_state<'a, T: FloatNumber + RealNumber, X: Array1<T>, LS: LineSearchMethod<T>>(
         &self,
         f: &'a F<'_, T, X>,
@@ -161,7 +155,6 @@ impl LBFGS {
         df(&mut state.x_df, &state.x);
     }
 
-    ///
     fn assess_convergence<T: FloatNumber, X: Array1<T>>(
         &self,
         state: &mut LBFGSState<T, X>,
@@ -173,7 +166,7 @@ impl LBFGS {
         }
 
         if state.x.max_diff(&state.x_prev)
-            <= T::from_f64(self.x_rtol * state.x.norm(std::f64::INFINITY)).unwrap()
+            <= T::from_f64(self.x_rtol * state.x.norm(f64::INFINITY)).unwrap()
         {
             x_converged = true;
         }
@@ -188,14 +181,13 @@ impl LBFGS {
             state.counter_f_tol += 1;
         }
 
-        if state.x_df.norm(std::f64::INFINITY) <= self.g_atol {
+        if state.x_df.norm(f64::INFINITY) <= self.g_atol {
             g_converged = true;
         }
 
         g_converged || x_converged || state.counter_f_tol > self.successive_f_tol
     }
 
-    ///
     fn update_hessian<T: FloatNumber, X: Array1<T>>(
         &self,
         _: &DF<'_, X>,
@@ -212,7 +204,6 @@ impl LBFGS {
     }
 }
 
-///
 #[derive(Debug)]
 struct LBFGSState<T: FloatNumber, X: Array1<T>> {
     x: X,
@@ -234,9 +225,7 @@ struct LBFGSState<T: FloatNumber, X: Array1<T>> {
     alpha: T,
 }
 
-///
 impl<T: FloatNumber + RealNumber> FirstOrderOptimizer<T> for LBFGS {
-    ///
     fn optimize<'a, X: Array1<T>, LS: LineSearchMethod<T>>(
         &self,
         f: &F<'_, T, X>,
@@ -248,7 +237,7 @@ impl<T: FloatNumber + RealNumber> FirstOrderOptimizer<T> for LBFGS {
 
         df(&mut state.x_df, x0);
 
-        let g_converged = state.x_df.norm(std::f64::INFINITY) < self.g_atol;
+        let g_converged = state.x_df.norm(f64::INFINITY) < self.g_atol;
         let mut converged = g_converged;
         let stopped = false;
 
@@ -299,7 +288,7 @@ mod tests {
 
         let result = optimizer.optimize(&f, &df, &x0, &ls);
 
-        assert!((result.f_x - 0.0).abs() < std::f64::EPSILON);
+        assert!((result.f_x - 0.0).abs() < f64::EPSILON);
         assert!((result.x[0] - 1.0).abs() < 1e-8);
         assert!((result.x[1] - 1.0).abs() < 1e-8);
         assert!(result.iterations <= 24);

src/optimization/first_order/mod.rs

@@ -1,6 +1,6 @@
-///
+/// Gradient descent optimization algorithm
 pub mod gradient_descent;
-///
+/// Limited-memory BFGS optimization algorithm
 pub mod lbfgs;
 
 use std::clone::Clone;
@@ -11,9 +11,9 @@ use crate::numbers::floatnum::FloatNumber;
 use crate::optimization::line_search::LineSearchMethod;
 use crate::optimization::{DF, F};
 
-///
+/// First-order optimization is a class of algorithms that use the first derivative of a function to find optimal solutions.
 pub trait FirstOrderOptimizer<T: FloatNumber> {
-    ///
+    /// run first order optimization
     fn optimize<'a, X: Array1<T>, LS: LineSearchMethod<T>>(
         &self,
         f: &F<'_, T, X>,
@@ -23,13 +23,13 @@ pub trait FirstOrderOptimizer<T: FloatNumber> {
     ) -> OptimizerResult<T, X>;
 }
 
-///
+/// Result of optimization
 #[derive(Debug, Clone)]
 pub struct OptimizerResult<T: FloatNumber, X: Array1<T>> {
-    ///
+    /// Solution
     pub x: X,
-    ///
+    /// f(x) value
     pub f_x: T,
-    ///
+    /// number of iterations
     pub iterations: usize,
 }

src/optimization/line_search.rs

@@ -1,11 +1,9 @@
-// TODO: missing documentation
-
 use crate::optimization::FunctionOrder;
 use num_traits::Float;
 
-///
+/// Line search optimization.
 pub trait LineSearchMethod<T: Float> {
-    ///
+    /// Find alpha that satisfies strong Wolfe conditions.
     fn search(
         &self,
         f: &(dyn Fn(T) -> T),
@@ -16,32 +14,31 @@ pub trait LineSearchMethod<T: Float> {
     ) -> LineSearchResult<T>;
 }
 
-///
+/// Line search result
 #[derive(Debug, Clone)]
 pub struct LineSearchResult<T: Float> {
-    ///
+    /// Alpha value
     pub alpha: T,
-    ///
+    /// f(alpha) value
     pub f_x: T,
 }
 
-///
+/// Backtracking line search method.
 pub struct Backtracking<T: Float> {
-    ///
+    /// TODO: Add documentation
     pub c1: T,
-    ///
+    /// Maximum number of iterations for Backtracking single run
     pub max_iterations: usize,
-    ///
+    /// TODO: Add documentation
     pub max_infinity_iterations: usize,
-    ///
+    /// TODO: Add documentation
     pub phi: T,
-    ///
+    /// TODO: Add documentation
     pub plo: T,
-    ///
+    /// function order
     pub order: FunctionOrder,
 }
 
-///
 impl<T: Float> Default for Backtracking<T> {
     fn default() -> Self {
         Backtracking {
@@ -55,9 +52,7 @@ impl<T: Float> Default for Backtracking<T> {
     }
 }
 
-///
 impl<T: Float> LineSearchMethod<T> for Backtracking<T> {
-    ///
     fn search(
         &self,
         f: &(dyn Fn(T) -> T),

src/optimization/mod.rs

@@ -1,21 +1,19 @@
-// TODO: missing documentation
-
-///
+/// first order optimization algorithms
 pub mod first_order;
-///
+/// line search algorithms
 pub mod line_search;
 
-///
+/// Function f(x) = y
 pub type F<'a, T, X> = dyn for<'b> Fn(&'b X) -> T + 'a;
-///
+/// Function df(x)
 pub type DF<'a, X> = dyn for<'b> Fn(&'b mut X, &'b X) + 'a;
 
-///
+/// Function order
 #[allow(clippy::upper_case_acronyms)]
 #[derive(Debug, PartialEq, Eq)]
 pub enum FunctionOrder {
-    ///
+    /// Second order
     SECOND,
-    ///
+    /// Third order
     THIRD,
 }

src/preprocessing/numerical.rs

@@ -172,18 +172,14 @@ where
     T: Number + RealNumber,
     M: Array2<T>,
 {
-    if let Some(output_matrix) = columns.first().cloned() {
-        return Some(
+    columns.first().cloned().map(|output_matrix| {
         columns
             .iter()
             .skip(1)
             .fold(output_matrix, |current_matrix, new_colum| {
                 current_matrix.h_stack(new_colum)
-                }),
-        );
-    } else {
-        None
-    }
+            })
+    })
 }
 
 #[cfg(test)]

src/readers/csv.rs

@@ -30,7 +30,7 @@ pub struct CSVDefinition<'a> {
     /// What seperates the fields in your csv-file?
     field_seperator: &'a str,
 }
-impl<'a> Default for CSVDefinition<'a> {
+impl Default for CSVDefinition<'_> {
     fn default() -> Self {
         Self {
             n_rows_header: 1,

src/svm/mod.rs

@@ -292,7 +292,7 @@ mod tests {
             .unwrap()
             .abs();
 
-        assert!((4913f64 - result) < std::f64::EPSILON);
+        assert!((4913f64 - result).abs() < f64::EPSILON);
     }
 
     #[cfg_attr(

src/svm/svc.rs

@@ -360,8 +360,8 @@ impl<'a, TX: Number + RealNumber, TY: Number + Ord, X: Array2<TX> + 'a, Y: Array
     }
 }
 
-impl<'a, TX: Number + RealNumber, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>> PartialEq
-    for SVC<'a, TX, TY, X, Y>
+impl<TX: Number + RealNumber, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>> PartialEq
+    for SVC<'_, TX, TY, X, Y>
 {
     fn eq(&self, other: &Self) -> bool {
         if (self.b.unwrap().sub(other.b.unwrap())).abs() > TX::epsilon() * TX::two()
@@ -1110,7 +1110,7 @@ mod tests {
         let svc = SVC::fit(&x, &y, &params).unwrap();
 
         // serialization
-        let deserialized_svc: SVC<f64, i32, _, _> =
+        let deserialized_svc: SVC<'_, f64, i32, _, _> =
             serde_json::from_str(&serde_json::to_string(&svc).unwrap()).unwrap();
 
         assert_eq!(svc, deserialized_svc);

src/svm/svr.rs

@@ -281,8 +281,8 @@ impl<'a, T: Number + FloatNumber + PartialOrd, X: Array2<T>, Y: Array1<T>> SVR<'
     }
 }
 
-impl<'a, T: Number + FloatNumber + PartialOrd, X: Array2<T>, Y: Array1<T>> PartialEq
-    for SVR<'a, T, X, Y>
+impl<T: Number + FloatNumber + PartialOrd, X: Array2<T>, Y: Array1<T>> PartialEq
+    for SVR<'_, T, X, Y>
 {
     fn eq(&self, other: &Self) -> bool {
         if (self.b - other.b).abs() > T::epsilon() * T::two()
@@ -702,7 +702,7 @@ mod tests {
 
         let svr = SVR::fit(&x, &y, &params).unwrap();
 
-        let deserialized_svr: SVR<f64, DenseMatrix<f64>, _> =
+        let deserialized_svr: SVR<'_, f64, DenseMatrix<f64>, _> =
             serde_json::from_str(&serde_json::to_string(&svr).unwrap()).unwrap();
 
         assert_eq!(svr, deserialized_svr);

src/tree/decision_tree_classifier.rs

@@ -77,7 +77,9 @@ use serde::{Deserialize, Serialize};
 
 use crate::api::{Predictor, SupervisedEstimator};
 use crate::error::Failed;
+use crate::linalg::basic::arrays::MutArray;
 use crate::linalg::basic::arrays::{Array1, Array2, MutArrayView1};
+use crate::linalg::basic::matrix::DenseMatrix;
 use crate::numbers::basenum::Number;
 use crate::rand_custom::get_rng_impl;
 
@@ -197,12 +199,12 @@ impl PartialEq for Node {
         self.output == other.output
             && self.split_feature == other.split_feature
             && match (self.split_value, other.split_value) {
-                (Some(a), Some(b)) => (a - b).abs() < std::f64::EPSILON,
+                (Some(a), Some(b)) => (a - b).abs() < f64::EPSILON,
                 (None, None) => true,
                 _ => false,
             }
             && match (self.split_score, other.split_score) {
-                (Some(a), Some(b)) => (a - b).abs() < std::f64::EPSILON,
+                (Some(a), Some(b)) => (a - b).abs() < f64::EPSILON,
                 (None, None) => true,
                 _ => false,
             }
@@ -613,7 +615,7 @@ impl<TX: Number + PartialOrd, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>>
             visitor_queue.push_back(visitor);
         }
 
-        while tree.depth() < tree.parameters().max_depth.unwrap_or(std::u16::MAX) {
+        while tree.depth() < tree.parameters().max_depth.unwrap_or(u16::MAX) {
             match visitor_queue.pop_front() {
                 Some(node) => tree.split(node, mtry, &mut visitor_queue, &mut rng),
                 None => break,
@@ -650,7 +652,7 @@ impl<TX: Number + PartialOrd, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>>
                     if node.true_child.is_none() && node.false_child.is_none() {
                         result = node.output;
                     } else if x.get((row, node.split_feature)).to_f64().unwrap()
-                        <= node.split_value.unwrap_or(std::f64::NAN)
+                        <= node.split_value.unwrap_or(f64::NAN)
                     {
                         queue.push_back(node.true_child.unwrap());
                     } else {
@@ -803,9 +805,7 @@ impl<TX: Number + PartialOrd, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>>
                     .get((i, self.nodes()[visitor.node].split_feature))
                     .to_f64()
                     .unwrap()
-                    <= self.nodes()[visitor.node]
-                        .split_value
-                        .unwrap_or(std::f64::NAN)
+                    <= self.nodes()[visitor.node].split_value.unwrap_or(f64::NAN)
                 {
                     *true_sample = visitor.samples[i];
                     tc += *true_sample;
@@ -889,11 +889,77 @@ impl<TX: Number + PartialOrd, TY: Number + Ord, X: Array2<TX>, Y: Array1<TY>>
         }
         importances
     }
+
+    /// Predict class probabilities for the input samples.
+    ///
+    /// # Arguments
+    ///
+    /// * `x` - The input samples as a matrix where each row is a sample and each column is a feature.
+    ///
+    /// # Returns
+    ///
+    /// A `Result` containing a `DenseMatrix<f64>` where each row corresponds to a sample and each column
+    /// corresponds to a class. The values represent the probability of the sample belonging to each class.
+    ///
+    /// # Errors
+    ///
+    /// Returns an error if at least one row prediction process fails.
+    pub fn predict_proba(&self, x: &X) -> Result<DenseMatrix<f64>, Failed> {
+        let (n_samples, _) = x.shape();
+        let n_classes = self.classes().len();
+        let mut result = DenseMatrix::<f64>::zeros(n_samples, n_classes);
+
+        for i in 0..n_samples {
+            let probs = self.predict_proba_for_row(x, i)?;
+            for (j, &prob) in probs.iter().enumerate() {
+                result.set((i, j), prob);
+            }
+        }
+
+        Ok(result)
+    }
+
+    /// Predict class probabilities for a single input sample.
+    ///
+    /// # Arguments
+    ///
+    /// * `x` - The input matrix containing all samples.
+    /// * `row` - The index of the row in `x` for which to predict probabilities.
+    ///
+    /// # Returns
+    ///
+    /// A vector of probabilities, one for each class, representing the probability
+    /// of the input sample belonging to each class.
+    fn predict_proba_for_row(&self, x: &X, row: usize) -> Result<Vec<f64>, Failed> {
+        let mut node = 0;
+
+        while let Some(current_node) = self.nodes().get(node) {
+            if current_node.true_child.is_none() && current_node.false_child.is_none() {
+                // Leaf node reached
+                let mut probs = vec![0.0; self.classes().len()];
+                probs[current_node.output] = 1.0;
+                return Ok(probs);
+            }
+
+            let split_feature = current_node.split_feature;
+            let split_value = current_node.split_value.unwrap_or(f64::NAN);
+
+            if x.get((row, split_feature)).to_f64().unwrap() <= split_value {
+                node = current_node.true_child.unwrap();
+            } else {
+                node = current_node.false_child.unwrap();
+            }
+        }
+
+        // This should never happen if the tree is properly constructed
+        Err(Failed::predict("Nodes iteration did not reach leaf"))
+    }
 }
 
 #[cfg(test)]
 mod tests {
     use super::*;
+    use crate::linalg::basic::arrays::Array;
     use crate::linalg::basic::matrix::DenseMatrix;
 
     #[test]
@@ -925,17 +991,62 @@ mod tests {
     )]
     #[test]
     fn gini_impurity() {
-        assert!((impurity(&SplitCriterion::Gini, &[7, 3], 10) - 0.42).abs() < std::f64::EPSILON);
+        assert!((impurity(&SplitCriterion::Gini, &[7, 3], 10) - 0.42).abs() < f64::EPSILON);
         assert!(
             (impurity(&SplitCriterion::Entropy, &[7, 3], 10) - 0.8812908992306927).abs()
-                < std::f64::EPSILON
+                < f64::EPSILON
         );
         assert!(
             (impurity(&SplitCriterion::ClassificationError, &[7, 3], 10) - 0.3).abs()
-                < std::f64::EPSILON
+                < f64::EPSILON
         );
     }
 
+    #[cfg_attr(
+        all(target_arch = "wasm32", not(target_os = "wasi")),
+        wasm_bindgen_test::wasm_bindgen_test
+    )]
+    #[test]
+    fn test_predict_proba() {
+        let x: DenseMatrix<f64> = DenseMatrix::from_2d_array(&[
+            &[5.1, 3.5, 1.4, 0.2],
+            &[4.9, 3.0, 1.4, 0.2],
+            &[4.7, 3.2, 1.3, 0.2],
+            &[4.6, 3.1, 1.5, 0.2],
+            &[5.0, 3.6, 1.4, 0.2],
+            &[7.0, 3.2, 4.7, 1.4],
+            &[6.4, 3.2, 4.5, 1.5],
+            &[6.9, 3.1, 4.9, 1.5],
+            &[5.5, 2.3, 4.0, 1.3],
+            &[6.5, 2.8, 4.6, 1.5],
+        ])
+        .unwrap();
+        let y: Vec<usize> = vec![0, 0, 0, 0, 0, 1, 1, 1, 1, 1];
+
+        let tree = DecisionTreeClassifier::fit(&x, &y, Default::default()).unwrap();
+        let probabilities = tree.predict_proba(&x).unwrap();
+
+        assert_eq!(probabilities.shape(), (10, 2));
+
+        for row in 0..10 {
+            let row_sum: f64 = probabilities.get_row(row).sum();
+            assert!(
+                (row_sum - 1.0).abs() < 1e-6,
+                "Row probabilities should sum to 1"
+            );
+        }
+
+        // Check if the first 5 samples have higher probability for class 0
+        for i in 0..5 {
+            assert!(probabilities.get((i, 0)) > probabilities.get((i, 1)));
+        }
+
+        // Check if the last 5 samples have higher probability for class 1
+        for i in 5..10 {
+            assert!(probabilities.get((i, 1)) > probabilities.get((i, 0)));
+        }
+    }
+
     #[cfg_attr(
         all(target_arch = "wasm32", not(target_os = "wasi")),
         wasm_bindgen_test::wasm_bindgen_test

src/tree/decision_tree_regressor.rs

@@ -311,15 +311,15 @@ impl Node {
 
 impl PartialEq for Node {
     fn eq(&self, other: &Self) -> bool {
-        (self.output - other.output).abs() < std::f64::EPSILON
+        (self.output - other.output).abs() < f64::EPSILON
             && self.split_feature == other.split_feature
             && match (self.split_value, other.split_value) {
-                (Some(a), Some(b)) => (a - b).abs() < std::f64::EPSILON,
+                (Some(a), Some(b)) => (a - b).abs() < f64::EPSILON,
                 (None, None) => true,
                 _ => false,
             }
             && match (self.split_score, other.split_score) {
-                (Some(a), Some(b)) => (a - b).abs() < std::f64::EPSILON,
+                (Some(a), Some(b)) => (a - b).abs() < f64::EPSILON,
                 (None, None) => true,
                 _ => false,
             }
@@ -478,7 +478,7 @@ impl<TX: Number + PartialOrd, TY: Number, X: Array2<TX>, Y: Array1<TY>>
             visitor_queue.push_back(visitor);
         }
 
-        while tree.depth() < tree.parameters().max_depth.unwrap_or(std::u16::MAX) {
+        while tree.depth() < tree.parameters().max_depth.unwrap_or(u16::MAX) {
             match visitor_queue.pop_front() {
                 Some(node) => tree.split(node, mtry, &mut visitor_queue, &mut rng),
                 None => break,
@@ -515,7 +515,7 @@ impl<TX: Number + PartialOrd, TY: Number, X: Array2<TX>, Y: Array1<TY>>
                     if node.true_child.is_none() && node.false_child.is_none() {
                         result = node.output;
                     } else if x.get((row, node.split_feature)).to_f64().unwrap()
-                        <= node.split_value.unwrap_or(std::f64::NAN)
+                        <= node.split_value.unwrap_or(f64::NAN)
                     {
                         queue.push_back(node.true_child.unwrap());
                     } else {
@@ -640,9 +640,7 @@ impl<TX: Number + PartialOrd, TY: Number, X: Array2<TX>, Y: Array1<TY>>
                     .get((i, self.nodes()[visitor.node].split_feature))
                     .to_f64()
                     .unwrap()
-                    <= self.nodes()[visitor.node]
-                        .split_value
-                        .unwrap_or(std::f64::NAN)
+                    <= self.nodes()[visitor.node].split_value.unwrap_or(f64::NAN)
                 {
                     *true_sample = visitor.samples[i];
                     tc += *true_sample;