Add feature importances to DecisionTreeRegressor, RandomForestRegressor and ExtraTreesRegressor

CHANGELOG.md

@@ -7,6 +7,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 - In the `...Regressor`s the features are always shuffled even when `mtry` equals the number of attributes. This changes the constructed trees. 
+- `tree`: `BaseTreeRegressor` (used by `DecisionTreeRegressor`, `RandomForestRegressor` and `ExtraTreesRegressor`) uses a new internal `Criterion` trait with `MSECriterion`. The MSE split improvement uses a numerically stable formula.  The selected splits do not change when a constant is added to y. Exactly tied splits can be broken in a different way, so fitted trees can differ slightly from earlier versions.
+- `tree`: in the `...Regressor`s, rows with sample weight 0 are removed before the tree grows. They no longer count toward `min_samples_leaf` and `min_samples_split`, and they no longer move a split threshold. A tree fitted with a zero-weight row is the same as a tree fitted without that row.
+- `tree`, `ensemble`: `DecisionTreeRegressor`, `RandomForestRegressor` and `ExtraTreesRegressor` get a `compute_feature_importances` method. 
+- `tree`, `ensemble`: the fitted `DecisionTreeRegressor`, `RandomForestRegressor` and `ExtraTreesRegressor` store new fields that the feature importances use. Models serialized with an earlier version cannot be deserialized. Fit these models again.
 
 ## [0.6.16]
 ### Fixed

src/ensemble/base_forest_regressor.rs

@@ -71,6 +71,7 @@ pub struct BaseForestRegressor<
     Y: Array1<TY>,
 > {
     trees: Option<Vec<BaseTreeRegressor<TX, TY, X, Y>>>,
+    num_features: Option<usize>,
     samples: Option<Vec<Vec<bool>>>,
 }
 
@@ -170,6 +171,7 @@ impl<TX: Number + FloatNumber + PartialOrd, TY: Number, X: Array2<TX>, Y: Array1
 
         Ok(BaseForestRegressor {
             trees: Some(trees),
+            num_features: Some(x.shape().1),
             samples: maybe_all_samples,
         })
     }
@@ -267,6 +269,42 @@ impl<TX: Number + FloatNumber + PartialOrd, TY: Number, X: Array2<TX>, Y: Array1
 
         samples
     }
+
+    /// Mean of the normalized tree importances, normalized again.
+    pub fn compute_feature_importances(&self) -> Result<Vec<f64>, Failed> {
+        match (&self.trees, self.num_features) {
+            (Some(trees), Some(num_features)) => {
+                let mut importances = vec![0.0f64; num_features];
+                let mut scratch = vec![0.0f64; num_features];
+
+                for tree in trees {
+                    tree.compute_feature_importances(true, &mut scratch);
+                    for (i, s) in importances.iter_mut().zip(scratch.iter()) {
+                        *i += s;
+                    }
+                }
+
+                // compute the mean of each feature importance
+                let num_trees = trees.len() as f64;
+                for i in importances.iter_mut() {
+                    *i /= num_trees;
+                }
+
+                // normalize once more
+                let sum: f64 = importances.iter().sum();
+                if sum > 0.0 {
+                    for i in importances.iter_mut() {
+                        *i /= sum;
+                    }
+                }
+
+                Ok(importances)
+            }
+            _ => Err(Failed::predict(
+                "'fit' should be called before calling 'compute_feature_importances'",
+            )),
+        }
+    }
 }
 
 #[cfg(test)]

src/ensemble/extra_trees_regressor.rs

@@ -272,6 +272,36 @@ impl<TX: Number + FloatNumber + PartialOrd, TY: Number, X: Array2<TX>, Y: Array1
             )),
         }
     }
+
+    /// Compute the feature importances for this forest.
+    /// * returns: a vector where index `j` gives the feature importances for feature `j`
+    ///
+    /// ```
+    /// use smartcore::ensemble::extra_trees_regressor::*;
+    /// use smartcore::linalg::basic::matrix::DenseMatrix;
+    ///
+    /// // Column 1 is constant, so it cannot split the data.
+    /// let rows: Vec<Vec<f64>> = (0..20).map(|i| vec![i as f64, 7.0]).collect();
+    /// let x = DenseMatrix::from_2d_vec(&rows).unwrap();
+    /// let y: Vec<f64> = (0..20).map(|i| (i * i) as f64).collect();
+    ///
+    /// let parameters = ExtraTreesRegressorParameters::default()
+    ///     .with_n_trees(10)
+    ///     .with_seed(42);
+    /// let forest = ExtraTreesRegressor::fit(&x, &y, parameters).unwrap();
+    ///
+    /// // All splits use column 0, so it gets all the importance.
+    /// // The importances of a forest are normalized: they add up to 1.
+    /// assert_eq!(forest.compute_feature_importances().unwrap(), vec![1.0, 0.0]);
+    /// ```
+    pub fn compute_feature_importances(&self) -> Result<Vec<f64>, Failed> {
+        match &self.forest_regressor {
+            Some(forest) => forest.compute_feature_importances(),
+            None => Err(Failed::predict(
+                "'fit' should be called before calling 'compute_feature_importances'",
+            )),
+        }
+    }
 }
 
 #[cfg(test)]
@@ -704,4 +734,139 @@ mod tests {
             );
         }
     }
+
+    /// Data where only column 0 carries information: `y` is a step function of column 0.
+    /// Columns 1 and 2 are permutations of the row index, so they do not separate the two steps.
+    fn step_data() -> (DenseMatrix<f64>, Vec<f64>) {
+        let n = 50;
+        let rows: Vec<Vec<f64>> = (0..n)
+            .map(|i| vec![i as f64, ((i * 37) % n) as f64, ((i * 17) % n) as f64])
+            .collect();
+        let y = (0..n).map(|i| if i < n / 2 { 0.0 } else { 10.0 }).collect();
+        (DenseMatrix::from_2d_vec(&rows).unwrap(), y)
+    }
+
+    #[test]
+    fn feature_importances_unfitted_returns_error() {
+        let et =
+            <ExtraTreesRegressor<f64, f64, DenseMatrix<f64>, Vec<f64>> as SupervisedEstimator<
+                DenseMatrix<f64>,
+                Vec<f64>,
+                ExtraTreesRegressorParameters,
+            >>::new();
+        let result = et.compute_feature_importances();
+        assert!(result.is_err());
+        assert_eq!(result.unwrap_err().error(), FailedError::PredictFailed);
+    }
+
+    #[test]
+    fn feature_importances_sum_to_one() {
+        let (x, y) = step_data();
+        for m in [1, 2, 3] {
+            let et = ExtraTreesRegressor::fit(
+                &x,
+                &y,
+                ExtraTreesRegressorParameters::default()
+                    .with_n_trees(50)
+                    .with_m(m)
+                    .with_seed(42),
+            )
+            .unwrap();
+            let importances = et.compute_feature_importances().unwrap();
+            assert_eq!(importances.len(), 3);
+            assert!(importances.iter().all(|&v| v >= 0.0 && v.is_finite()));
+            assert!((importances.iter().sum::<f64>() - 1.0).abs() < 1e-12);
+        }
+    }
+
+    #[test]
+    fn feature_importances_constant_column_gets_zero() {
+        // The random splitter skips a constant column, so every split is on column 0,
+        // whatever the random thresholds are.
+        let x = DenseMatrix::from_2d_vec(&(0..20).map(|i| vec![i as f64, 7.0]).collect::<Vec<_>>())
+            .unwrap();
+        let y: Vec<f64> = (0..20).map(|i| (i * i) as f64).collect();
+        let et = ExtraTreesRegressor::fit(
+            &x,
+            &y,
+            ExtraTreesRegressorParameters::default()
+                .with_n_trees(10)
+                .with_seed(42),
+        )
+        .unwrap();
+        assert_eq!(et.compute_feature_importances().unwrap(), vec![1.0, 0.0]);
+    }
+
+    #[test]
+    fn feature_importances_find_informative_feature() {
+        let (x, y) = step_data();
+        for m in [1, 3] {
+            let et = ExtraTreesRegressor::fit(
+                &x,
+                &y,
+                ExtraTreesRegressorParameters::default()
+                    .with_n_trees(100)
+                    .with_m(m)
+                    .with_seed(42),
+            )
+            .unwrap();
+            let importances = et.compute_feature_importances().unwrap();
+            assert!(importances[0] > importances[1], "m = {m}: {importances:?}");
+            assert!(importances[0] > importances[2], "m = {m}: {importances:?}");
+        }
+    }
+
+    #[test]
+    fn feature_importances_do_not_depend_on_weight_scale() {
+        // There is no bootstrap, so the weights go directly to the trees.
+        // Multiplying all weights by 4 (exact in f64) must give the same importances.
+        let (x, y) = step_data();
+        let weights: Vec<f64> = (0..y.len()).map(|i| 0.5 + (i % 4) as f64 * 0.25).collect();
+        let scaled: Vec<f64> = weights.iter().map(|w| 4.0 * w).collect();
+        let parameters = ExtraTreesRegressorParameters::default()
+            .with_n_trees(20)
+            .with_m(2)
+            .with_seed(42);
+
+        let a = ExtraTreesRegressor::fit_with_weights(&x, &y, &weights, parameters.clone())
+            .unwrap()
+            .compute_feature_importances()
+            .unwrap();
+        let b = ExtraTreesRegressor::fit_with_weights(&x, &y, &scaled, parameters)
+            .unwrap()
+            .compute_feature_importances()
+            .unwrap();
+        for (ai, bi) in a.iter().zip(b.iter()) {
+            assert!((ai - bi).abs() < 1e-12);
+        }
+    }
+
+    #[test]
+    fn feature_importances_are_zero_without_improvement() {
+        let (x, y) = step_data();
+
+        // max_depth 1: every tree is only a root.
+        let et = ExtraTreesRegressor::fit(
+            &x,
+            &y,
+            ExtraTreesRegressorParameters::default()
+                .with_n_trees(5)
+                .with_max_depth(1)
+                .with_seed(42),
+        )
+        .unwrap();
+        assert_eq!(et.compute_feature_importances().unwrap(), vec![0.0; 3]);
+
+        // Constant y: trees can split, but no split decreases the impurity.
+        let y_const = vec![3.0; y.len()];
+        let et = ExtraTreesRegressor::fit(
+            &x,
+            &y_const,
+            ExtraTreesRegressorParameters::default()
+                .with_n_trees(5)
+                .with_seed(42),
+        )
+        .unwrap();
+        assert_eq!(et.compute_feature_importances().unwrap(), vec![0.0; 3]);
+    }
 }