diff --git a/src/ensemble/base_forest_regressor.rs b/src/ensemble/base_forest_regressor.rs index 550d600e..b48f8dad 100644 --- a/src/ensemble/base_forest_regressor.rs +++ b/src/ensemble/base_forest_regressor.rs @@ -5,6 +5,7 @@ use std::fmt::Debug; use serde::{Deserialize, Serialize}; use crate::error::{Failed, FailedError}; +use crate::linalg::basic::arrays::MutArrayView1; use crate::linalg::basic::arrays::{Array1, Array2}; use crate::numbers::basenum::Number; use crate::numbers::floatnum::FloatNumber; @@ -120,6 +121,14 @@ impl, Y: Array1 }) .transpose()?; + // Compute the order of each attribute once + let mut order: Vec> = Vec::with_capacity(num_attributes); + + for i in 0..num_attributes { + let mut col_i: Vec = x.get_col(i).iterator(0).copied().collect(); + order.push(col_i.argsort_mut()); + } + for tree_idx in 0..parameters.n_trees { if parameters.bootstrap { samples = BaseForestRegressor::::sample_with_replacement( @@ -139,7 +148,7 @@ impl, Y: Array1 min_samples_leaf: parameters.min_samples_leaf, min_samples_split: parameters.min_samples_split, seed: Some(parameters.seed.wrapping_add(tree_idx as u64)), // give each tree its own fixed seed - splitter: parameters.splitter.clone(), + splitter: parameters.splitter, }; // Only use sample weights on base tree if not already applied during bootstrapping let sample_weights_for_base_tree = if parameters.bootstrap { @@ -153,6 +162,7 @@ impl, Y: Array1 sample_weights_for_base_tree, samples.clone(), mtry, + &order, params, )?; trees.push(tree); diff --git a/src/ensemble/extra_trees_regressor.rs b/src/ensemble/extra_trees_regressor.rs index c0dd1a99..a5a06f77 100644 --- a/src/ensemble/extra_trees_regressor.rs +++ b/src/ensemble/extra_trees_regressor.rs @@ -277,6 +277,7 @@ impl, Y: Array1 #[cfg(test)] mod tests { use super::*; + use crate::error::FailedError; use crate::linalg::basic::matrix::DenseMatrix; use crate::metrics::mean_squared_error; @@ -470,6 +471,30 @@ mod tests { assert_eq!(yhat.err(), Some(Failed::predict(msg))); } + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn fit_with_empty_x_should_not_panic() { + let parameters = ExtraTreesRegressorParameters::default() + .with_n_trees(5) + .with_seed(42); + // (rows, columns, y): no rows, no columns, or both + let cases: Vec<(usize, usize, Vec)> = + vec![(0, 0, vec![]), (0, 3, vec![]), (3, 0, vec![1.0, 2.0, 3.0])]; + for (nrows, ncols, y) in cases { + let x: DenseMatrix = DenseMatrix::new(nrows, ncols, vec![], false) + .expect("Construction of empty x should work"); + let result = ExtraTreesRegressor::fit(&x, &y, parameters.clone()); + let expected = Failed::because( + FailedError::ParametersError, + "Training data must contain at least one sample and one feature.", + ); + assert_eq!(result.err(), Some(expected), "shape: ({nrows}, {ncols})"); + } + } + mod sklearn_parity { use super::*; // sklearn parity tests. diff --git a/src/ensemble/random_forest_regressor.rs b/src/ensemble/random_forest_regressor.rs index 9dabdf2d..364806e4 100644 --- a/src/ensemble/random_forest_regressor.rs +++ b/src/ensemble/random_forest_regressor.rs @@ -470,6 +470,7 @@ impl, Y: Array1 #[cfg(test)] mod tests { use super::*; + use crate::error::FailedError; use crate::linalg::basic::matrix::DenseMatrix; use crate::metrics::mean_absolute_error; @@ -801,6 +802,30 @@ mod tests { assert_eq!(yhat.err(), Some(Failed::predict(msg))); } + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn fit_with_empty_x_should_not_panic() { + let parameters = RandomForestRegressorParameters::default() + .with_n_trees(5) + .with_seed(42); + // (rows, columns, y): no rows, no columns, or both + let cases: Vec<(usize, usize, Vec)> = + vec![(0, 0, vec![]), (0, 3, vec![]), (3, 0, vec![1.0, 2.0, 3.0])]; + for (nrows, ncols, y) in cases { + let x: DenseMatrix = DenseMatrix::new(nrows, ncols, vec![], false) + .expect("Construction of empty x should work"); + let result = RandomForestRegressor::fit(&x, &y, parameters.clone()); + let expected = Failed::because( + FailedError::ParametersError, + "Training data must contain at least one sample and one feature.", + ); + assert_eq!(result.err(), Some(expected), "shape: ({nrows}, {ncols})"); + } + } + mod sklearn_parity { use super::*; // sklearn parity tests. diff --git a/src/tree/base_tree_regressor.rs b/src/tree/base_tree_regressor.rs index b833e345..56059824 100644 --- a/src/tree/base_tree_regressor.rs +++ b/src/tree/base_tree_regressor.rs @@ -1,5 +1,3 @@ -use std::collections::LinkedList; -use std::default::Default; use std::fmt::Debug; use std::marker::PhantomData; @@ -15,7 +13,7 @@ use crate::numbers::basenum::Number; use crate::rand_custom::get_rng_impl; #[cfg_attr(feature = "serde", derive(Serialize, Deserialize))] -#[derive(Debug, Clone, Default)] +#[derive(Debug, Clone, Copy, Default)] pub enum Splitter { Random, #[default] @@ -115,10 +113,10 @@ impl, Y: Array1> PartialE for BaseTreeRegressor { fn eq(&self, other: &Self) -> bool { - if self.depth != other.depth || self.nodes().len() != other.nodes().len() { + if self.depth != other.depth || self.nodes.len() != other.nodes.len() { false } else { - self.nodes() + self.nodes .iter() .zip(other.nodes().iter()) .all(|(a, b)| a == b) @@ -130,9 +128,8 @@ struct NodeVisitor<'a, TX: Number + PartialOrd, TY: Number, X: Array2, Y: Ar x: &'a X, y: &'a Y, node: usize, - samples: Vec, - sample_weights: Option<&'a [f64]>, - order: &'a [Vec], + start_idx: usize, // start index for the elements in `sorted_by_feature` of SplitWorkspace + end_idx: usize, // end index (exclusive) in the same vector(s) true_child_output: f64, false_child_output: f64, level: u16, @@ -145,9 +142,8 @@ impl<'a, TX: Number + PartialOrd, TY: Number, X: Array2, Y: Array1> { fn new( node_id: usize, - samples: Vec, - sample_weights: Option<&'a [f64]>, - order: &'a [Vec], + start_idx: usize, + end_idx: usize, x: &'a X, y: &'a Y, level: u16, @@ -156,9 +152,8 @@ impl<'a, TX: Number + PartialOrd, TY: Number, X: Array2, Y: Array1> x, y, node: node_id, - samples, - sample_weights, - order, + start_idx, + end_idx, true_child_output: 0f64, false_child_output: 0f64, level, @@ -166,11 +161,6 @@ impl<'a, TX: Number + PartialOrd, TY: Number, X: Array2, Y: Array1> _phantom_ty: PhantomData, } } - - /// Weighted count of sample `i`. The weight is 1.0 if no weights are given. - fn mass_of(&self, i: usize) -> f64 { - mass_of(i, &self.samples, self.sample_weights) - } } /// Weighted count of sample `i`. The weight is 1.0 if no weights are given. @@ -181,6 +171,177 @@ fn mass_of(i: usize, samples: &[usize], sample_weights: Option<&[f64]>) -> f64 { } } +/// Node elements store row indices and counts as `u32`. +fn check_row_count(n_rows: usize) -> Result<(), Failed> { + if u32::try_from(n_rows).is_err() { + return Err(Failed::fit( + "Number of rows in x must not be larger than u32::MAX", + )); + } + Ok(()) +} + +// Trait representing an element that belongs logically to a Node. Implementations of it are stored in SplitWorkspace +trait NodeElement: Copy + Default { + fn new(row: usize, count: usize, mass: f64) -> Self; + // The row index in the dataset + fn row(&self) -> usize; + // The number of times this row was sampled, should always be > 0 + fn count(&self) -> usize; + // The total mass of this element: count * sample_weight of row or count if no sample weights are used + fn mass(&self) -> f64; +} + +#[derive(Copy, Clone, Default)] +struct WeightedElement { + row_idx: u32, + count: u32, + mass: f64, +} + +impl NodeElement for WeightedElement { + fn new(row_idx: usize, count: usize, mass: f64) -> Self { + Self { + row_idx: row_idx as u32, + count: count as u32, + mass, + } + } + + #[inline(always)] + fn row(&self) -> usize { + self.row_idx as usize + } + + #[inline(always)] + fn count(&self) -> usize { + self.count as usize + } + + #[inline(always)] + fn mass(&self) -> f64 { + self.mass + } +} + +// implement CountedElement as a single u64, so it can be retrieved in a single read +#[derive(Copy, Clone, Default)] +struct CountedElement(u64); // low 32 bits: row, high 32 bits: count + +impl NodeElement for CountedElement { + fn new(row_idx: usize, count: usize, mass: f64) -> Self { + // mass() is derived from count, so sample weights must not be used + debug_assert_eq!(mass, count as f64); + Self(row_idx as u64 | ((count as u64) << 32)) + } + #[inline(always)] + fn row(&self) -> usize { + self.0 as u32 as usize + } + #[inline(always)] + fn count(&self) -> usize { + (self.0 >> 32) as usize + } + #[inline(always)] + fn mass(&self) -> f64 { + (self.0 >> 32) as f64 + } +} + +#[derive(Copy, Clone, Default)] +struct UnitElement { + row_idx: u32, +} + +impl NodeElement for UnitElement { + fn new(row_idx: usize, count: usize, mass: f64) -> Self { + // count() and mass() are always 1, so sample weights and counts > 1 must not be used. + // A count of 0 is allowed: such elements are overwritten during setup. + debug_assert!(count <= 1); + debug_assert_eq!(mass, count as f64); + Self { + row_idx: row_idx as u32, + } + } + + #[inline(always)] + fn row(&self) -> usize { + self.row_idx as usize + } + + #[inline(always)] + fn count(&self) -> usize { + 1 + } + + #[inline(always)] + fn mass(&self) -> f64 { + 1.0 + } +} + +// Checks whether the example indicated by node_element is a "true child" for this node +fn is_true_sample(node_element: &E, x: &X, node: &Node) -> bool +where + TX: Number + PartialOrd, + X: Array2, + E: NodeElement, +{ + x.get((node_element.row(), node.split_feature)) + .to_f64() + .unwrap() + <= node.split_value.unwrap_or(f64::NAN) +} + +// slice: slice that will be partitioned +// scratch: temp buffer +// is_true: is_true[idx] checks whether element with row idx equal to idx belongs to the true branch +// returns: index of first element of false branch +#[inline(never)] // inline(never) is faster as seen during profiling +fn stable_partition(slice: &mut [E], scratch: &mut [E], is_true: &[bool]) -> usize +where + E: NodeElement, +{ + // Note: this is intentionally written without an explicit if/else branch in the main loop + let n = slice.len(); + let scratch = &mut scratch[..n]; + let (mut w, mut f) = (0usize, 0usize); + for i in 0..n { + let e = slice[i]; + let t = is_true[e.row()]; + slice[w] = e; // w <= i, so this never clobbers an unread element + scratch[f] = e; + // advance only one of the pointers + w += t as usize; + f += (!t) as usize; + } + slice[w..].copy_from_slice(&scratch[..f]); + w +} + +struct SplitWorkspace { + in_true_branch: Vec, // indexed by row index in X + partition_buffer: Vec, + sorted_by_feature: Vec>, + variables: Vec, // feature indices which will be shuffled when mtry < n_features +} + +impl SplitWorkspace { + fn new( + in_true_branch: Vec, + partition_buffer: Vec, + sorted_by_feature: Vec>, + n_features: usize, + ) -> Self { + Self { + in_true_branch, + partition_buffer, + sorted_by_feature, + variables: (0..n_features).collect::>(), + } + } +} + impl, Y: Array1> BaseTreeRegressor { @@ -201,6 +362,14 @@ impl, Y: Array1> )); } + // Compute the order of each attribute once + let mut order: Vec> = Vec::new(); + + for i in 0..num_attributes { + let mut col_i: Vec = x.get_col(i).iterator(0).copied().collect(); + order.push(col_i.argsort_mut()); + } + let samples = vec![1; x_nrows]; BaseTreeRegressor::fit_weak_learner( x, @@ -208,6 +377,7 @@ impl, Y: Array1> sample_weights, samples, num_attributes, + &order, parameters, ) } @@ -218,12 +388,39 @@ impl, Y: Array1> sample_weights: Option<&[f64]>, samples: Vec, mtry: usize, + order: &[Vec], parameters: BaseTreeRegressorParameters, ) -> Result, Failed> { - let y_m = y.clone(); + check_row_count(y.shape())?; + match sample_weights { + Some(_) => Self::grow::( + x, + y, + sample_weights, + samples, + mtry, + order, + parameters, + ), + None if samples.iter().all(|&s| s <= 1) => { + // Note: if any sample has count zero, it will be filtered out, so we should + // use the smallest possible NodeElement + Self::grow::(x, y, None, samples, mtry, order, parameters) + } + None => Self::grow::(x, y, None, samples, mtry, order, parameters), + } + } - let y_ncols = y_m.shape(); - let (_, num_attributes) = x.shape(); + fn grow( + x: &X, + y: &Y, + sample_weights: Option<&[f64]>, + samples: Vec, + mtry: usize, + order: &[Vec], + parameters: BaseTreeRegressorParameters, + ) -> Result, Failed> { + let n_rows = y.shape(); let mut nodes: Vec = Vec::new(); let mut rng = get_rng_impl(parameters.seed); @@ -231,20 +428,42 @@ impl, Y: Array1> let mut sum = 0f64; let mut mass = 0f64; - for i in 0..y_ncols { + for i in 0..n_rows { let mass_i = mass_of(i, &samples, sample_weights); mass += mass_i; - sum += mass_i * y_m.get(i).to_f64().unwrap(); + sum += mass_i * y.get(i).to_f64().unwrap(); } let root = Node::new(sum / mass); nodes.push(root); - let mut order: Vec> = Vec::new(); - for i in 0..num_attributes { - let mut col_i: Vec = x.get_col(i).iterator(0).copied().collect(); - order.push(col_i.argsort_mut()); - } + let counts: Vec = samples.iter().map(|&s| s as u32).collect(); + // number of distinct rows of x in this tree + let n_kept = counts.iter().filter(|&&c| c > 0).count(); + + let sorted_by_feature: Vec> = order + .iter() + .map(|col_order| { + let mut out = vec![E::default(); n_kept + 1]; // preallocate, one additional place + let mut w = 0usize; // index to write to + for &i in col_order { + let c = counts[i]; + out[w] = E::new(i, c as usize, mass_of(i, &samples, sample_weights)); + w += (c > 0) as usize; // update w in a branchless way + } + out.truncate(n_kept); + out + }) + .collect(); + + let end_idx = sorted_by_feature[0].len(); + + let mut workspace = SplitWorkspace::new( + vec![false; x.shape().0], + vec![E::default(); end_idx], + sorted_by_feature, + x.shape().1, + ); let mut base_tree = BaseTreeRegressor { nodes, @@ -256,19 +475,18 @@ impl, Y: Array1> _phantom_y: PhantomData, }; - let mut visitor = - NodeVisitor::::new(0, samples, sample_weights, &order, x, &y_m, 1); + let mut visitor = NodeVisitor::::new(0, 0, end_idx, x, y, 1); - let mut visitor_queue: LinkedList> = LinkedList::new(); + let mut visitor_stack: Vec> = Vec::new(); - if base_tree.find_best_cutoff(&mut visitor, mtry, mass, &mut rng) { - visitor_queue.push_back(visitor); + if base_tree.find_best_cutoff(&mut visitor, mtry, mass, &mut rng, &mut workspace) { + visitor_stack.push(visitor); } let max_depth = base_tree.parameters().max_depth.unwrap_or(u16::MAX); - while let Some(node) = visitor_queue.pop_front() { + while let Some(node) = visitor_stack.pop() { if node.level < max_depth { - base_tree.split(node, mtry, &mut visitor_queue, &mut rng); + base_tree.split(node, mtry, &mut visitor_stack, &mut rng, &mut workspace); } } @@ -292,7 +510,7 @@ impl, Y: Array1> pub(crate) fn predict_for_row(&self, x: &X, row: usize) -> TY { let mut node_id = 0; loop { - let node = &self.nodes()[node_id]; + let node = &self.nodes[node_id]; let Some(true_child) = node.true_child else { return TY::from_f64(node.output).unwrap(); }; @@ -307,49 +525,60 @@ impl, Y: Array1> } } - fn find_best_cutoff( + fn find_best_cutoff( &mut self, visitor: &mut NodeVisitor<'_, TX, TY, X, Y>, mtry: usize, mass: f64, rng: &mut impl rand::Rng, + workspace: &mut SplitWorkspace, ) -> bool { let (_, n_attr) = visitor.x.shape(); - let n: usize = visitor.samples.iter().sum(); + let n: usize = workspace.sorted_by_feature[0][visitor.start_idx..visitor.end_idx] + .iter() + .map(|elem| elem.count()) + .sum(); if n < self.parameters().min_samples_split { return false; } - let sum = self.nodes()[visitor.node].output * mass; - - let mut variables = (0..n_attr).collect::>(); + let sum = self.nodes[visitor.node].output * mass; + // Note: in sklearn, the attributes are always considered in a random order if mtry < n_attr { - variables.shuffle(rng); + workspace.variables.shuffle(rng); } - let parent_gain = - mass * self.nodes()[visitor.node].output * self.nodes()[visitor.node].output; + let parent_gain = mass * self.nodes[visitor.node].output * self.nodes[visitor.node].output; - let splitter = self.parameters().splitter.clone(); + let splitter = self.parameters().splitter; - for variable in variables.iter().take(mtry) { + for variable in workspace.variables.iter().take(mtry) { match splitter { Splitter::Random => { - self.find_random_split(visitor, n, mass, sum, parent_gain, *variable, rng); + self.find_random_split( + visitor, + n, + mass, + sum, + parent_gain, + *variable, + rng, + workspace, + ); } Splitter::Best => { - self.find_best_split(visitor, n, mass, sum, parent_gain, *variable); + self.find_best_split(visitor, n, mass, sum, parent_gain, *variable, workspace); } } } - self.nodes()[visitor.node].split_score.is_some() + self.nodes[visitor.node].split_score.is_some() } - fn find_random_split( + fn find_random_split( &mut self, visitor: &mut NodeVisitor<'_, TX, TY, X, Y>, n: usize, @@ -358,46 +587,37 @@ impl, Y: Array1> parent_gain: f64, j: usize, rng: &mut impl rand::Rng, + workspace: &SplitWorkspace, ) { - let (min_val, max_val) = { - let mut min_opt = None; - let mut max_opt = None; - for &i in &visitor.order[j] { - if visitor.samples[i] > 0 { - min_opt = Some(*visitor.x.get((i, j))); - break; - } - } - for &i in visitor.order[j].iter().rev() { - if visitor.samples[i] > 0 { - max_opt = Some(*visitor.x.get((i, j))); - break; - } - } - if min_opt.is_none() { - return; - } - (min_opt.unwrap(), max_opt.unwrap()) - }; + if visitor.start_idx == visitor.end_idx { + return; + } + let first_elem = workspace.sorted_by_feature[j][visitor.start_idx]; + let min_val = visitor.x.get((first_elem.row(), j)); + let last_elem = workspace.sorted_by_feature[j][visitor.end_idx - 1]; + let max_val = visitor.x.get((last_elem.row(), j)); + + // Convert min_val and max_val to f64 first, so that values that are different in type TX + // but equal in f64 do not end up panicking the call to the rng.random_range function. + let min_val = min_val.to_f64().unwrap(); + let max_val = max_val.to_f64().unwrap(); if min_val >= max_val { return; } - let split_value = rng.random_range(min_val.to_f64().unwrap()..max_val.to_f64().unwrap()); + let split_value = rng.random_range(min_val..max_val); let mut true_sum = 0f64; let mut true_mass = 0f64; let mut true_count = 0; - for &i in &visitor.order[j] { - if visitor.samples[i] > 0 { - if visitor.x.get((i, j)).to_f64().unwrap() <= split_value { - true_sum += visitor.mass_of(i) * visitor.y.get(i).to_f64().unwrap(); - true_count += visitor.samples[i]; - true_mass += visitor.mass_of(i); - } else { - break; - } + for elem in &workspace.sorted_by_feature[j][visitor.start_idx..visitor.end_idx] { + if visitor.x.get((elem.row(), j)).to_f64().unwrap() <= split_value { + true_sum += elem.mass() * visitor.y.get(elem.row()).to_f64().unwrap(); + true_count += elem.count(); + true_mass += elem.mass(); + } else { + break; } } @@ -434,7 +654,7 @@ impl, Y: Array1> } } - fn find_best_split( + fn find_best_split( &mut self, visitor: &mut NodeVisitor<'_, TX, TY, X, Y>, n: usize, @@ -442,107 +662,126 @@ impl, Y: Array1> sum: f64, parent_gain: f64, j: usize, + workspace: &SplitWorkspace, ) { let mut true_sum = 0f64; let mut true_count = 0; let mut true_mass = 0f64; let mut prevx = Option::None; - for i in visitor.order[j].iter() { - if visitor.samples[*i] > 0 { - let x_ij = *visitor.x.get((*i, j)); + for elem in &workspace.sorted_by_feature[j][visitor.start_idx..visitor.end_idx] { + // Comparison for equality is done directly on f64 to avoid later problems + // with values being unequal in TX but equal in f64. + let x_ij = visitor.x.get((elem.row(), j)).to_f64().unwrap(); - if prevx.is_none() || x_ij == prevx.unwrap() { - prevx = Some(x_ij); - true_count += visitor.samples[*i]; - true_mass += visitor.mass_of(*i); - true_sum += visitor.mass_of(*i) * visitor.y.get(*i).to_f64().unwrap(); - continue; - } + if prevx.is_none() || x_ij == prevx.unwrap() { + prevx = Some(x_ij); + true_count += elem.count(); + true_mass += elem.mass(); + true_sum += elem.mass() * visitor.y.get(elem.row()).to_f64().unwrap(); + continue; + } - let false_count = n - true_count; + let false_count = n - true_count; - if true_count < self.parameters().min_samples_leaf - || false_count < self.parameters().min_samples_leaf - { - prevx = Some(x_ij); - true_count += visitor.samples[*i]; - true_mass += visitor.mass_of(*i); - true_sum += visitor.mass_of(*i) * visitor.y.get(*i).to_f64().unwrap(); - continue; - } + if true_count < self.parameters().min_samples_leaf + || false_count < self.parameters().min_samples_leaf + { + prevx = Some(x_ij); + true_count += elem.count(); + true_mass += elem.mass(); + true_sum += elem.mass() * visitor.y.get(elem.row()).to_f64().unwrap(); + continue; + } - let true_mean = if true_mass > 0.0 { - true_sum / true_mass - } else { - 0.0 - }; - let false_mass = mass - true_mass; - let false_mean = if false_mass > 0.0 { - (sum - true_sum) / false_mass - } else { - 0.0 - }; - - let gain = (true_mass * true_mean * true_mean - + false_mass * false_mean * false_mean) - - parent_gain; - - if self.nodes()[visitor.node].split_score.is_none() - || gain > self.nodes()[visitor.node].split_score.unwrap() - { - self.nodes[visitor.node].split_feature = j; - self.nodes[visitor.node].split_value = - Option::Some((x_ij + prevx.unwrap()).to_f64().unwrap() / 2f64); - self.nodes[visitor.node].split_score = Option::Some(gain); - - visitor.true_child_output = true_mean; - visitor.false_child_output = false_mean; - } + let true_mean = if true_mass > 0.0 { + true_sum / true_mass + } else { + 0.0 + }; + let false_mass = mass - true_mass; + let false_mean = if false_mass > 0.0 { + (sum - true_sum) / false_mass + } else { + 0.0 + }; - prevx = Some(x_ij); - true_sum += visitor.mass_of(*i) * visitor.y.get(*i).to_f64().unwrap(); - true_count += visitor.samples[*i]; - true_mass += visitor.mass_of(*i); + let gain = (true_mass * true_mean * true_mean + false_mass * false_mean * false_mean) + - parent_gain; + + if self.nodes[visitor.node].split_score.is_none() + || gain > self.nodes[visitor.node].split_score.unwrap() + { + let prev = prevx.unwrap(); + // compute the midpoint in a way that does not overflow when x_ij and prevx have the same sign + let mid = prev + (x_ij - prev) / 2.0; + + // The midpoint should put prevx in the true branch, and x_ij in the false branch. + // However, mid could round up to x_ij, in which case we would also assign x_ij to the true branch. + // In this case, we set the midpoint equal to prevx, which is always smaller + // than x_ij. + let threshold = if mid < x_ij { mid } else { prev }; + self.nodes[visitor.node].split_feature = j; + self.nodes[visitor.node].split_value = Option::Some(threshold); + self.nodes[visitor.node].split_score = Option::Some(gain); + + visitor.true_child_output = true_mean; + visitor.false_child_output = false_mean; } + + prevx = Some(x_ij); + true_sum += elem.mass() * visitor.y.get(elem.row()).to_f64().unwrap(); + true_count += elem.count(); + true_mass += elem.mass(); } } - fn split<'a>( + /// Apply the split that was found for `visitor.node`: add its two child nodes and + /// partition the node's range in each `sorted_by_feature` column into a true part + /// and a false part. Then find the best cutoff for each child. If a child can be + /// split, push it on `visitor_stack`. + /// If a branch has fewer than `min_samples_leaf` samples, clear the split and + /// keep the node as a leaf. + fn split<'a, E: NodeElement>( &mut self, - mut visitor: NodeVisitor<'a, TX, TY, X, Y>, + visitor: NodeVisitor<'a, TX, TY, X, Y>, mtry: usize, - visitor_queue: &mut LinkedList>, + visitor_stack: &mut Vec>, rng: &mut impl rand::Rng, + workspace: &mut SplitWorkspace, ) -> bool { - let (n, _) = visitor.x.shape(); - let mut tc = 0; - let mut fc = 0; + let this_node = &self.nodes[visitor.node]; + + let mut true_count = 0usize; let mut true_mass = 0f64; + let mut false_count = 0usize; let mut false_mass = 0f64; - let mut true_samples: Vec = vec![0; n]; - - for (i, true_sample) in true_samples.iter_mut().enumerate().take(n) { - if visitor.samples[i] > 0 { - if visitor - .x - .get((i, self.nodes()[visitor.node].split_feature)) - .to_f64() - .unwrap() - <= self.nodes()[visitor.node].split_value.unwrap_or(f64::NAN) - { - *true_sample = visitor.samples[i]; - tc += *true_sample; - true_mass += visitor.mass_of(i); - visitor.samples[i] = 0; - } else { - fc += visitor.samples[i]; - false_mass += visitor.mass_of(i); - } + // for each row_index, does it belong in the true branch or not? + let in_true_branch = &mut workspace.in_true_branch; + + // the column for split_feature is in order, so all the true samples should come first, all the false samples second + let split_feature = this_node.split_feature; + let col = &workspace.sorted_by_feature[split_feature][visitor.start_idx..visitor.end_idx]; + let mut n_true = 0usize; + for e in col { + if !is_true_sample(e, visitor.x, this_node) { + break; } + in_true_branch[e.row()] = true; + true_count += e.count(); + true_mass += e.mass(); + n_true += 1; + } + for e in &col[n_true..] { + in_true_branch[e.row()] = false; + false_count += e.count(); + false_mass += e.mass(); } - if tc < self.parameters().min_samples_leaf || fc < self.parameters().min_samples_leaf { + // Stop early if it is clear that there will be too few examples in the leaf + if true_count < self.parameters().min_samples_leaf + || false_count < self.parameters().min_samples_leaf + { self.nodes[visitor.node].split_feature = 0; self.nodes[visitor.node].split_value = Option::None; self.nodes[visitor.node].split_score = Option::None; @@ -550,43 +789,66 @@ impl, Y: Array1> return false; } - let true_child_idx = self.nodes().len(); + // Add the child nodes to the tree + let split_idx = visitor.start_idx + n_true; + + let true_child_idx = self.nodes.len(); self.nodes.push(Node::new(visitor.true_child_output)); - let false_child_idx = self.nodes().len(); + let false_child_idx = self.nodes.len(); self.nodes.push(Node::new(visitor.false_child_output)); self.nodes[visitor.node].true_child = Some(true_child_idx); self.nodes[visitor.node].false_child = Some(false_child_idx); - self.depth = u16::max(self.depth, visitor.level + 1); + let child_level = visitor.level + 1; + self.depth = u16::max(self.depth, child_level); + + // If the child nodes can not be split any further, there is no point in partitioning the ranges + let max_depth = self.parameters().max_depth.unwrap_or(u16::MAX); + let min_split = self.parameters().min_samples_split; + let true_can_split = child_level < max_depth && true_count >= min_split; + let false_can_split = child_level < max_depth && false_count >= min_split; + if !true_can_split && !false_can_split { + return true; // both children are leaves: no partition + } + + for j in 0..visitor.x.shape().1 { + stable_partition( + &mut workspace.sorted_by_feature[j][visitor.start_idx..visitor.end_idx], + &mut workspace.partition_buffer, + in_true_branch, + ); + } let mut true_visitor = NodeVisitor::::new( true_child_idx, - true_samples, - visitor.sample_weights, - visitor.order, + visitor.start_idx, + split_idx, visitor.x, visitor.y, - visitor.level + 1, + child_level, ); - if self.find_best_cutoff(&mut true_visitor, mtry, true_mass, rng) { - visitor_queue.push_back(true_visitor); + if true_can_split + && self.find_best_cutoff(&mut true_visitor, mtry, true_mass, rng, workspace) + { + visitor_stack.push(true_visitor); } let mut false_visitor = NodeVisitor::::new( false_child_idx, - visitor.samples, - visitor.sample_weights, - visitor.order, + split_idx, + visitor.end_idx, visitor.x, visitor.y, - visitor.level + 1, + child_level, ); - if self.find_best_cutoff(&mut false_visitor, mtry, false_mass, rng) { - visitor_queue.push_back(false_visitor); + if false_can_split + && self.find_best_cutoff(&mut false_visitor, mtry, false_mass, rng, workspace) + { + visitor_stack.push(false_visitor); } true @@ -600,6 +862,19 @@ mod tests { use crate::linalg::basic::matrix::DenseMatrix; use crate::metrics::mean_absolute_error; + #[test] + #[cfg(target_pointer_width = "64")] + fn check_row_count_rejects_more_rows_than_u32_max() { + assert!(check_row_count(0).is_ok()); + assert!(check_row_count(u32::MAX as usize).is_ok()); + assert_eq!( + check_row_count(u32::MAX as usize + 1).err(), + Some(Failed::fit( + "Number of rows in x must not be larger than u32::MAX" + )) + ); + } + #[test] fn test_fit_on_empty_data_returns_error() { // 2 rows x 2 features — values are arbitrary; only the empty-row case is under test @@ -967,4 +1242,537 @@ mod tests { let y_hat_repeated = tree_repeated.predict(&x).expect("Predict should work"); assert!(mean_absolute_error(&y_hat, &y_hat_repeated) < 1e-9); } + + // Above 2^53, f64 cannot hold every integer. Split search compares values in TX, + // but partitioning and predict compare them as f64 against an f64 threshold. + const TWO_POW_53: i64 = 1 << 53; + + fn split_value_parameters(splitter: Splitter) -> BaseTreeRegressorParameters { + BaseTreeRegressorParameters { + max_depth: None, + min_samples_leaf: 1, + min_samples_split: 2, + seed: Some(42), + splitter, + } + } + + #[test] + fn large_integer_midpoint_threshold_separates_values() { + // Both values are exact in f64, but the midpoint 2^53 + 3 is not. + // It rounds to 2^53 + 4 (ties to even), so both rows go to the true branch. + let x = + DenseMatrix::from_2d_vec(&vec![vec![TWO_POW_53 + 2], vec![TWO_POW_53 + 4]]).unwrap(); + let y = vec![0.0f64, 10.0]; + + let tree = + BaseTreeRegressor::fit_inner(&x, &y, None, split_value_parameters(Splitter::Best)) + .expect("Fit should work"); + + let threshold = tree.nodes()[0].split_value.expect("Root should split"); + assert!((TWO_POW_53 + 2) as f64 <= threshold); + assert!(((TWO_POW_53 + 4) as f64) > threshold); + + let y_hat = tree.predict(&x).expect("Predict should work"); + assert_eq!(y_hat, vec![0.0, 10.0]); + } + + #[test] + fn large_integer_values_equal_as_f64_do_not_split() { + // 2^53 and 2^53 + 1 differ as i64 but are the same f64. No f64 threshold can + // separate them, so the tree must not split and must predict the mean. + let x = DenseMatrix::from_2d_vec(&vec![vec![TWO_POW_53], vec![TWO_POW_53 + 1]]).unwrap(); + let y = vec![0.0f64, 10.0]; + + let tree = + BaseTreeRegressor::fit_inner(&x, &y, None, split_value_parameters(Splitter::Best)) + .expect("Fit should work"); + + let y_hat = tree.predict(&x).expect("Predict should work"); + assert_eq!(y_hat, vec![5.0, 5.0]); + } + + #[test] + fn large_integer_values_equal_as_f64_random_splitter() { + // min_val < max_val as i64, but the f64 range is empty. This must not panic. + let x = DenseMatrix::from_2d_vec(&vec![vec![TWO_POW_53], vec![TWO_POW_53 + 1]]).unwrap(); + let y = vec![0.0f64, 10.0]; + + let tree = + BaseTreeRegressor::fit_inner(&x, &y, None, split_value_parameters(Splitter::Random)) + .expect("Fit should work"); + + let y_hat = tree.predict(&x).expect("Predict should work"); + assert_eq!(y_hat, vec![5.0, 5.0]); + } + + #[test] + fn adjacent_f64_midpoint_threshold_separates_values() { + // a and b are adjacent f64 values. Their exact midpoint 1 + 1.5 * EPSILON is not + // representable. It is a tie, and ties round to the even mantissa, which is b. + // So the threshold becomes b, and both rows go to the true branch. + let a = 1.0 + f64::EPSILON; + let b = 1.0 + 2.0 * f64::EPSILON; + assert_eq!((a + b) / 2.0, b); // the rounding that causes the problem + + let x = DenseMatrix::from_2d_vec(&vec![vec![a], vec![b]]).unwrap(); + let y = vec![0.0f64, 10.0]; + + let tree = + BaseTreeRegressor::fit_inner(&x, &y, None, split_value_parameters(Splitter::Best)) + .expect("Fit should work"); + + let threshold = tree.nodes()[0].split_value.expect("Root should split"); + assert!(a <= threshold); + assert!(b > threshold); + + let y_hat = tree.predict(&x).expect("Predict should work"); + assert_eq!(y_hat, vec![0.0, 10.0]); + } + + #[test] + fn test_node_element_sizes() { + assert_eq!(std::mem::size_of::(), 4); + assert_eq!(std::mem::size_of::(), 8); + assert_eq!(std::mem::size_of::(), 16); + } + + type TestTree = BaseTreeRegressor, Vec>; + + fn params( + max_depth: Option, + min_samples_leaf: usize, + min_samples_split: usize, + splitter: Splitter, + ) -> BaseTreeRegressorParameters { + BaseTreeRegressorParameters { + max_depth, + min_samples_leaf, + min_samples_split, + seed: Some(42), + splitter, + } + } + + /// Row order of each column, computed in the same way as in `fit_inner`. + fn order_of(x: &DenseMatrix) -> Vec> { + (0..x.shape().1) + .map(|j| { + let mut col: Vec = x.get_col(j).iterator(0).copied().collect(); + col.argsort_mut() + }) + .collect() + } + + /// Index of the leaf that gets row `row` of `x`. Uses the same rule as `predict_for_row`. + fn leaf_of(tree: &TestTree, x: &DenseMatrix, row: usize) -> usize { + let mut node_id = 0; + while let Some(true_child) = tree.nodes()[node_id].true_child { + let node = &tree.nodes()[node_id]; + node_id = if *x.get((row, node.split_feature)) <= node.split_value.unwrap() { + true_child + } else { + node.false_child.unwrap() + }; + } + node_id + } + + /// Number of training samples in each leaf of `tree`. + /// `counts[i]` is the number of times that row `i` of `x` is in the training set. + fn leaf_counts(tree: &TestTree, x: &DenseMatrix, counts: &[usize]) -> Vec { + let mut per_node = vec![0usize; tree.nodes().len()]; + for (row, &count) in counts.iter().enumerate() { + per_node[leaf_of(tree, x, row)] += count; + } + tree.nodes() + .iter() + .zip(per_node) + .filter(|(node, _)| node.true_child.is_none()) // filter out the leaves + .map(|(_, count)| count) + .collect() + } + + /// Both trees must have the same structure, the same splits and the same outputs. + /// `split_score` is not compared, because a different summation order can change it. + fn assert_same_tree(a: &TestTree, b: &TestTree) { + assert_eq!(a.nodes().len(), b.nodes().len()); + assert_eq!(a.depth, b.depth); + for (na, nb) in a.nodes().iter().zip(b.nodes().iter()) { + assert_eq!(na.split_feature, nb.split_feature); + assert_eq!(na.split_value, nb.split_value); + assert_eq!(na.true_child, nb.true_child); + assert_eq!(na.false_child, nb.false_child); + assert!((na.output - nb.output).abs() < 1e-9); + } + } + + /// Repeat row `i` of `x`, `y` and `w` `counts[i]` times. + fn repeat_rows( + x: &[Vec], + y: &[f64], + w: &[f64], + counts: &[usize], + ) -> (DenseMatrix, Vec, Vec) { + let mut x_rep = Vec::new(); + let mut y_rep = Vec::new(); + let mut w_rep = Vec::new(); + for (i, &count) in counts.iter().enumerate() { + for _ in 0..count { + x_rep.push(x[i].clone()); + y_rep.push(y[i]); + w_rep.push(w[i]); + } + } + (DenseMatrix::from_2d_vec(&x_rep).unwrap(), y_rep, w_rep) + } + + /// Data for the element path tests: ties in both columns, integer targets. + fn element_path_data() -> (Vec>, Vec) { + let x = vec![ + vec![1.0, 4.0], + vec![2.0, 4.0], + vec![2.0, 1.0], + vec![3.0, 3.0], + vec![5.0, 2.0], + vec![6.0, 1.0], + vec![6.0, 5.0], + vec![7.0, 3.0], + ]; + let y = vec![3.0, 8.0, 1.0, 9.0, 4.0, 12.0, 7.0, 2.0]; + (x, y) + } + + #[test] + fn max_depth_boundary() { + let x = DenseMatrix::from_iterator((1..=8).map(|i| i as f64), 8, 1, 0); + let y = vec![0.0f64, 1.0, 10.0, 11.0, 100.0, 101.0, 110.0, 111.0]; + + // The root is at level 1. A node at level `max_depth` is not split. + // (max_depth, number of nodes, depth) + let cases = [ + (Some(1), 1, 0), + (Some(2), 3, 2), + (Some(3), 7, 3), + (None, 15, 4), + ]; + for (max_depth, n_nodes, depth) in cases { + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(max_depth, 1, 2, Splitter::Best)) + .expect("Fit should work"); + assert_eq!(tree.nodes().len(), n_nodes, "max_depth: {max_depth:?}"); + assert_eq!(tree.depth, depth, "max_depth: {max_depth:?}"); + } + + // With max_depth = 2 the root split separates the two halves. + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(Some(2), 1, 2, Splitter::Best)) + .expect("Fit should work"); + assert_eq!(tree.nodes()[0].split_value, Some(4.5)); + let y_hat = tree.predict(&x).expect("Predict should work"); + assert_eq!(y_hat, vec![5.5, 5.5, 5.5, 5.5, 105.5, 105.5, 105.5, 105.5]); + + // No limit is the same as the largest limit, and the leaves memorize y. + let tree_none: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 2, Splitter::Best)) + .expect("Fit should work"); + let tree_max: TestTree = BaseTreeRegressor::fit_inner( + &x, + &y, + None, + params(Some(u16::MAX), 1, 2, Splitter::Best), + ) + .expect("Fit should work"); + assert_same_tree(&tree_none, &tree_max); + assert_eq!(tree_none.predict(&x).expect("Predict should work"), y); + + for max_depth in 1..=4 { + let tree: TestTree = BaseTreeRegressor::fit_inner( + &x, + &y, + None, + params(Some(max_depth), 1, 2, Splitter::Random), + ) + .expect("Fit should work"); + assert!(tree.depth <= max_depth, "max_depth: {max_depth}"); + } + } + + #[test] + fn min_samples_split_one_child_only() { + // The root split gives a true child with 3 samples and a false child with 2 samples. + // With min_samples_split = 3, only the true child can split. + let x = DenseMatrix::from_2d_vec(&vec![ + vec![1.0_f64], + vec![2.0], + vec![3.0], + vec![4.0], + vec![5.0], + ]) + .unwrap(); + let y = vec![0.0f64, 1.0, 2.0, 100.0, 101.0]; + + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 3, Splitter::Best)) + .expect("Fit should work"); + assert_eq!(tree.nodes()[0].split_value, Some(3.5)); + let true_child = tree.nodes()[0].true_child.unwrap(); + let false_child = tree.nodes()[0].false_child.unwrap(); + assert!(tree.nodes()[true_child].true_child.is_some()); + assert!(tree.nodes()[false_child].true_child.is_none()); + assert_eq!(tree.nodes().len(), 5); + assert_eq!(tree.depth, 3); + + // If min_samples_split is larger than the number of samples, the root is a leaf. + for splitter in [Splitter::Best, Splitter::Random] { + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 6, splitter)) + .expect("Fit should work"); + assert_eq!(tree.nodes().len(), 1, "splitter: {splitter:?}"); + } + } + + #[test] + fn min_samples_leaf_boundary() { + let x = DenseMatrix::from_2d_vec(&vec![vec![1.0_f64], vec![2.0], vec![3.0], vec![4.0]]) + .unwrap(); + let y = vec![0.0f64, 0.0, 0.0, 100.0]; + + // With no limit, the best split isolates the last sample. + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 2, Splitter::Best)) + .expect("Fit should work"); + assert_eq!(tree.nodes()[0].split_value, Some(3.5)); + + // The best split is not permitted. The split with exactly 2 samples on each side is used. + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 2, 2, Splitter::Best)) + .expect("Fit should work"); + assert_eq!(tree.nodes()[0].split_value, Some(2.5)); + assert_eq!(tree.nodes().len(), 3); + assert_eq!(leaf_counts(&tree, &x, &[1; 4]), vec![2, 2]); + + // No split can give 3 samples on each side. + for splitter in [Splitter::Best, Splitter::Random] { + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 3, 2, splitter)) + .expect("Fit should work"); + assert_eq!(tree.nodes().len(), 1, "splitter: {splitter:?}"); + assert_eq!(tree.nodes()[0].output, 25.0); + } + } + + // Fit many different trees with many different options and check that in each case min_samples_leaf is correct. + #[test] + fn min_samples_leaf_holds_for_all_element_paths() { + let mut rng = get_rng_impl(Some(7)); + let n = 40; + let x_rows: Vec> = (0..n) + .map(|_| (0..3).map(|_| rng.random_range(0..6) as f64).collect()) + .collect(); + let y: Vec = (0..n).map(|_| rng.random_range(-10.0..10.0)).collect(); + let counts: Vec = (0..n).map(|_| rng.random_range(0..4usize)).collect(); + let weights: Vec = (0..n).map(|_| rng.random_range(0.1..3.0)).collect(); + let ones = vec![1usize; n]; + let x = DenseMatrix::from_2d_vec(&x_rows).unwrap(); + let order = order_of(&x); + + // (weights, counts): UnitElement, CountedElement and WeightedElement + let paths = [ + (None, &ones), + (None, &counts), + (Some(weights.as_slice()), &counts), + ]; + for (sample_weights, samples) in paths { + for splitter in [Splitter::Best, Splitter::Random] { + for min_samples_leaf in [1, 3, 5] { + for seed in 0..5 { + let parameters = BaseTreeRegressorParameters { + seed: Some(seed), + ..params(None, min_samples_leaf, 2, splitter) + }; + let tree: TestTree = BaseTreeRegressor::fit_weak_learner( + &x, + &y, + sample_weights, + samples.clone(), + 3, + &order, + parameters, + ) + .expect("Fit should work"); + let message = format!( + "weights: {}, splitter: {splitter:?}, leaf: {min_samples_leaf}, seed: {seed}", + sample_weights.is_some() + ); + assert!(tree.nodes().len() > 1, "{message}"); // make sure the tree has a split + let leaves = leaf_counts(&tree, &x, samples); + assert!( + leaves.iter().all(|&c| c >= min_samples_leaf), + "{message}, leaves: {leaves:?}" + ); + } + } + } + } + } + + #[test] + fn constant_features_do_not_split() { + // all features are constant, so no split is possible + let x = DenseMatrix::from_2d_vec(&vec![ + vec![1.0_f64, 7.0], + vec![1.0, 7.0], + vec![1.0, 7.0], + vec![1.0, 7.0], + ]) + .unwrap(); + let y = vec![1.0f64, 2.0, 3.0, 6.0]; + + for splitter in [Splitter::Best, Splitter::Random] { + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 2, splitter)) + .expect("Fit should work"); + assert_eq!(tree.nodes().len(), 1, "splitter: {splitter:?}"); + assert_eq!(tree.depth, 0, "splitter: {splitter:?}"); + let y_hat = tree.predict(&x).expect("Predict should work"); + assert_eq!(y_hat, vec![3.0; 4], "splitter: {splitter:?}"); + } + } + + #[test] + fn constant_column_is_skipped() { + // Column 0 is constant, so only column 1 can give a split. + let x = DenseMatrix::from_2d_vec(&vec![ + vec![7.0_f64, 1.0], + vec![7.0, 2.0], + vec![7.0, 3.0], + vec![7.0, 4.0], + ]) + .unwrap(); + let y = vec![0.0f64, 0.0, 10.0, 10.0]; + + for splitter in [Splitter::Best, Splitter::Random] { + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 2, splitter)) + .expect("Fit should work"); + assert!(tree.nodes().len() > 1, "splitter: {splitter:?}"); + assert!( + tree.nodes() + .iter() + .filter(|node| node.true_child.is_some()) + .all(|node| node.split_feature == 1), + "splitter: {splitter:?}" + ); + } + + let tree: TestTree = + BaseTreeRegressor::fit_inner(&x, &y, None, params(None, 1, 2, Splitter::Best)) + .expect("Fit should work"); + assert_eq!(tree.nodes()[0].split_value, Some(2.5)); + } + + #[test] + fn duplicate_columns_pick_first() { + // Both columns are identical. Their gains are equal, and the first column is kept. + // Note: this is the current behaviour. The behaviour of sklearn is different. + let values = [1.0_f64, 2.0, 2.0, 3.0, 5.0, 6.0]; + let y = vec![1.0f64, 4.0, 2.0, 9.0, 3.0, 8.0]; + let x_single = + DenseMatrix::from_2d_vec(&values.iter().map(|&v| vec![v]).collect::>()).unwrap(); + let x_double = + DenseMatrix::from_2d_vec(&values.iter().map(|&v| vec![v, v]).collect::>()) + .unwrap(); + + let parameters = params(None, 1, 2, Splitter::Best); + let tree_single: TestTree = + BaseTreeRegressor::fit_inner(&x_single, &y, None, parameters.clone()) + .expect("Fit should work"); + let tree_double: TestTree = + BaseTreeRegressor::fit_inner(&x_double, &y, None, parameters).expect("Fit should work"); + + assert!(tree_double.nodes().len() > 1); + assert_same_tree(&tree_single, &tree_double); + } + + #[test] + fn counted_and_unit_elements_match_repeated_rows() { + // Checks that repeating rows and giving them a count > 0 yields the same trees. + // Do not pose constraints on min_samples_leaf and min_samples_split. At the moment + // this should make no difference but it may do in the future. + let (x_rows, y) = element_path_data(); + let x = DenseMatrix::from_2d_vec(&x_rows).unwrap(); + let order = order_of(&x); + let n = y.len(); + + // Counts 0 or 1 use UnitElement. Counts larger than 1 use CountedElement. + // Rows with count 0 must have no effect. + let all_counts = [vec![1, 0, 1, 1, 0, 1, 1, 0], vec![2, 0, 1, 3, 1, 0, 2, 1]]; + for counts in all_counts { + let (x_rep, y_rep, _) = repeat_rows(&x_rows, &y, &vec![1.0; n], &counts); + for splitter in [Splitter::Best, Splitter::Random] { + for max_depth in [Some(2), None] { + let parameters = params(max_depth, 1, 2, splitter); + let tree: TestTree = BaseTreeRegressor::fit_weak_learner( + &x, + &y, + None, + counts.clone(), + 2, + &order, + parameters.clone(), + ) + .expect("Fit should work"); + let tree_rep: TestTree = + BaseTreeRegressor::fit_inner(&x_rep, &y_rep, None, parameters) + .expect("Fit should work"); + + assert!(tree.nodes().len() > 1); + assert_same_tree(&tree, &tree_rep); + let y_hat = tree.predict(&x).expect("Predict should work"); + let y_hat_rep = tree_rep.predict(&x).expect("Predict should work"); + assert!(mean_absolute_error(&y_hat, &y_hat_rep) < 1e-9); + } + } + } + } + + #[test] + fn weighted_element_with_counts_matches_repeated_rows() { + // The mass of a row is count * weight. Weights are chosen to be perfectly representable (power of two in the fraction denominator). + // This keeps all sums exact. + // Checks that WeightedElement handles counts different from one correctly. + let (x_rows, y) = element_path_data(); + let x = DenseMatrix::from_2d_vec(&x_rows).unwrap(); + let order = order_of(&x); + let weights = vec![0.5, 1.5, 2.0, 0.25, 1.0, 3.0, 0.75, 1.25]; + let counts = vec![2, 0, 1, 3, 1, 0, 2, 1]; + let (x_rep, y_rep, w_rep) = repeat_rows(&x_rows, &y, &weights, &counts); + + for splitter in [Splitter::Best, Splitter::Random] { + for max_depth in [Some(2), None] { + let parameters = params(max_depth, 1, 2, splitter); + let tree: TestTree = BaseTreeRegressor::fit_weak_learner( + &x, + &y, + Some(&weights), + counts.clone(), + 2, + &order, + parameters.clone(), + ) + .expect("Fit should work"); + // in this tree, there are no counts but the rows are repeated physically + let tree_rep: TestTree = + BaseTreeRegressor::fit_inner(&x_rep, &y_rep, Some(&w_rep), parameters) + .expect("Fit should work"); + + assert!(tree.nodes().len() > 1); + assert_same_tree(&tree, &tree_rep); + let y_hat = tree.predict(&x).expect("Predict should work"); + let y_hat_rep = tree_rep.predict(&x).expect("Predict should work"); + assert!(mean_absolute_error(&y_hat, &y_hat_rep) < 1e-9); + } + } + } } diff --git a/src/tree/decision_tree_regressor.rs b/src/tree/decision_tree_regressor.rs index 839cfdac..12ab72b6 100644 --- a/src/tree/decision_tree_regressor.rs +++ b/src/tree/decision_tree_regressor.rs @@ -370,6 +370,7 @@ impl, Y: Array1> #[cfg(test)] mod tests { use super::*; + use crate::error::FailedError; use crate::linalg::basic::matrix::DenseMatrix; #[test] @@ -629,6 +630,27 @@ mod tests { assert_eq!(yhat.err(), Some(Failed::predict(msg))); } + #[cfg_attr( + all(target_arch = "wasm32", not(target_os = "wasi")), + wasm_bindgen_test::wasm_bindgen_test + )] + #[test] + fn fit_with_empty_x_should_not_panic() { + // (rows, columns, y): no rows, no columns, or both + let cases: Vec<(usize, usize, Vec)> = + vec![(0, 0, vec![]), (0, 3, vec![]), (3, 0, vec![1.0, 2.0, 3.0])]; + for (nrows, ncols, y) in cases { + let x: DenseMatrix = DenseMatrix::new(nrows, ncols, vec![], false) + .expect("Construction of empty x should work"); + let result = DecisionTreeRegressor::fit(&x, &y, Default::default()); + let expected = Failed::because( + FailedError::ParametersError, + "Training data must contain at least one sample and one feature.", + ); + assert_eq!(result.err(), Some(expected), "shape: ({nrows}, {ncols})"); + } + } + #[cfg_attr( all(target_arch = "wasm32", not(target_os = "wasi")), wasm_bindgen_test::wasm_bindgen_test