Skip to main content

contree/algorithms/
contree_lds.rs

1use crate::algorithms::depth2::{ConTreeDepth2, Depth2Loss};
2use crate::algorithms::interval_pruner::{Bound, IntervalsPruner};
3use crate::algorithms::shared;
4use crate::algorithms::support_feasible_splits;
5use crate::caching::{Cache, Entry, SearchedUnder};
6use crate::common::{
7    Budget, BudgetSchedule, Cost, FitOutcome, Loss, PassReport, PointSelector, ScheduleBounds,
8    ScheduleKind, SearchConfig, SearchError, SearchStatus, Skipping, Statistics,
9};
10use crate::data::view::DataView;
11use crate::data::{Dataset, Feature};
12use crate::tree::Tree;
13use rand::rngs::StdRng;
14use rand::SeedableRng;
15use std::collections::VecDeque;
16use std::time::Instant;
17
18/// Anytime ConTree search based on limited discrepancy search.
19///
20/// The search runs in passes. Each pass explores the tree space under a
21/// discrepancy budget (how far down the Gini ranking of features it may go)
22/// and, for some point selectors, a split budget (how many thresholds per
23/// feature it may try). A [`BudgetSchedule`] widens both budgets between
24/// passes, so a good tree is available early and the last complete pass
25/// proves optimality.
26pub struct ConTreeLds {
27    config: SearchConfig,
28    statistics: Statistics,
29    cache: Cache,
30    specialized: ConTreeDepth2,
31    runtime: Instant,
32    rng: StdRng,
33    max_discrepancy: usize,
34    split_budget: usize,
35    schedule_kind: ScheduleKind,
36    schedule: Option<Box<dyn BudgetSchedule>>,
37    /// Which budgets truncated the pass in progress, read by the schedule.
38    pass_report: PassReport,
39    status: SearchStatus,
40    /// Every improvement of the root incumbent, as `(seconds, error)`.
41    trajectory: Vec<(f64, Cost)>,
42    loss: Loss,
43    /// The first negative or NaN cost the loss returned during the fit.
44    invalid_cost: Option<Cost>,
45}
46
47impl ConTreeLds {
48    /// Builds a solver from a [`SearchConfig`].
49    ///
50    /// This is the preferred constructor; [`Self::new`] takes the same
51    /// settings as positional arguments.
52    pub fn with_config(config: SearchConfig) -> Self {
53        let mut solver = Self::new(
54            config.min_sup,
55            config.max_depth,
56            config.max_time,
57            config.max_error,
58            config.point_selector,
59            config.max_gap,
60            config.use_heuristic,
61            config.fast_d2,
62        );
63        solver.config = config;
64        solver
65    }
66
67    /// Builds a solver from positional settings. See [`Self::with_config`].
68    #[allow(clippy::too_many_arguments)]
69    pub fn new(
70        min_sup: usize,
71        max_depth: usize,
72        max_time: f64,
73        max_error: Cost,
74        split_selection_strategy: PointSelector,
75        max_gap: Cost,
76        use_heuristic: bool,
77        fast_d2: bool,
78    ) -> Self {
79        let config = SearchConfig::new(
80            min_sup,
81            max_depth,
82            max_time,
83            max_gap,
84            max_error,
85            use_heuristic,
86            fast_d2,
87            split_selection_strategy,
88        );
89
90        Self {
91            cache: Cache::default(),
92            config,
93            statistics: Statistics::default(),
94            specialized: ConTreeDepth2,
95            runtime: Instant::now(),
96            rng: StdRng::from_os_rng(),
97            max_discrepancy: usize::MAX,
98            split_budget: 1,
99            schedule_kind: ScheduleKind::default(),
100            schedule: None,
101            pass_report: PassReport::default(),
102            status: SearchStatus::Optimal,
103            trajectory: Vec::new(),
104            loss: Loss::default(),
105            invalid_cost: None,
106        }
107    }
108
109    /// Minimises `loss` instead of misclassification.
110    ///
111    /// A loss that reads instance ids turns the depth-2 solver off, since it
112    /// works on class counts only.
113    pub fn with_loss(mut self, loss: Loss) -> Self {
114        self.loss = loss;
115        self
116    }
117
118    /// Chooses whether thresholds are skipped; see [`Skipping`].
119    pub fn with_skipping(mut self, skipping: Skipping) -> Self {
120        self.config.skipping = skipping;
121        self
122    }
123
124    /// The first negative or NaN cost the loss returned, if any. `fit`
125    /// reports it as an error; a caller of `partial_fit` checks it here.
126    pub fn invalid_cost(&self) -> Option<Cost> {
127        self.invalid_cost
128    }
129
130    /// Seeds the random generator used by `PointSelector::Random`, making
131    /// runs reproducible. Unseeded solvers draw their seed from the OS.
132    pub fn with_random_state(mut self, seed: u64) -> Self {
133        self.rng = StdRng::seed_from_u64(seed);
134        self
135    }
136
137    /// Runs the anytime search until it stops.
138    ///
139    /// Each pass widens the discrepancy and split budgets. The search ends
140    /// when a pass completes without being truncated (the tree is optimal),
141    /// when the schedule has no larger budget to offer, or when the time
142    /// limit is reached.
143    pub fn fit(&mut self, dataset: &Dataset) -> Result<FitOutcome, SearchError> {
144        crate::algorithms::validate(&self.config, dataset)?;
145
146        let root_view = DataView::root(dataset, self.config.use_heuristic);
147        while !self.partial_fit(&root_view) {}
148        if let Some(cost) = self.invalid_cost {
149            return Err(SearchError::InvalidLeafCost { cost });
150        }
151
152        let tree = self.get_solution_tree();
153        tree.validate()?;
154
155        Ok(FitOutcome {
156            tree,
157            statistics: self.statistics,
158            status: self.status,
159        })
160    }
161
162    /// Chooses how the budget widens from pass to pass. See [`ScheduleKind`].
163    pub fn with_schedule(mut self, schedule: ScheduleKind) -> Self {
164        self.schedule_kind = schedule;
165        self
166    }
167
168    /// The budget schedule this solver uses.
169    pub fn schedule(&self) -> ScheduleKind {
170        self.schedule_kind
171    }
172
173    fn apply(&mut self, budget: Budget) {
174        self.config.budget = budget.discrepancy;
175        self.split_budget = budget.split_budget();
176    }
177
178    /// The configuration this solver was built with.
179    pub fn config(&self) -> &SearchConfig {
180        &self.config
181    }
182
183    /// The anytime profile of the last fit: `(seconds, error)` at every
184    /// improvement of the root incumbent.
185    pub fn trajectory(&self) -> &[(f64, Cost)] {
186        &self.trajectory
187    }
188
189    /// Why the last `fit` stopped.
190    pub fn status(&self) -> SearchStatus {
191        self.status
192    }
193
194    /// Runs one pass of the search and returns `true` once the search is
195    /// over. The first call initialises the cache and the budget schedule.
196    pub fn partial_fit(&mut self, root_view: &DataView<'_>) -> bool {
197        self.config.nb_runs += 1;
198        let mut root_index = 0;
199        let mut entry = Entry::default();
200        if self.config.nb_runs <= 1 {
201            self.config.cost_step = SearchConfig::cost_step_for(&self.loss);
202            self.config.fast_d2 &= !self.loss.uses_indices();
203            self.invalid_cost = None;
204            self.statistics.skipping_violations = 0;
205            self.cache = Cache::new(self.config.max_depth, root_view.total_instances);
206
207            self.max_discrepancy = self.max_discrepancy.min(Self::discrepancy_limit(
208                root_view.get_feature_number(),
209                self.config.max_depth,
210            ));
211            let mut schedule = self.schedule_kind.build(ScheduleBounds {
212                max_discrepancy: self.max_discrepancy,
213                max_splits: root_view.get_max_splits(),
214                split_applies: self.config.point_selector == PointSelector::First,
215            });
216            let budget = schedule.first();
217            self.apply(budget);
218            self.schedule = Some(schedule);
219            root_index = self.cache.init();
220
221            self.statistics.num_features = root_view.get_feature_number();
222            self.statistics.num_samples = root_view.total_instances;
223
224            let (error, label) = self.leaf(root_view);
225            entry.error = error;
226            entry.label = label;
227
228            if let Some(entry) = self.cache.get_mut(root_index) {
229                self.config.max_error = self.config.max_error.min(error);
230                entry.error = error;
231                entry.label = label;
232                entry.lower_bound = 0.0;
233                entry.is_valid = true;
234            }
235
236            self.runtime = Instant::now();
237            self.trajectory.clear();
238        }
239
240        let mut is_optimal = false;
241
242        if let Some(e) = self.cache.get(root_index) {
243            entry = *e
244        }
245        if entry.is_optimal {
246            return true;
247        }
248
249        let mut config = self.config;
250
251        config.discrepancy = 0;
252
253        self.pass_report = PassReport::default();
254        let error = entry.error;
255        let stopped = self.expand_node_with_view(root_view, &config, &mut entry, 0, true, error);
256        self.pass_report.improved = entry.error < error;
257
258        // Once the schedule has no larger budget, another pass would repeat
259        // the same truncated search, so the search stops here.
260        let next_budget = self
261            .schedule
262            .as_mut()
263            .and_then(|schedule| schedule.next(&self.pass_report));
264        let search_exhausted = next_budget.is_none();
265
266        if !stopped || !self.time_remains() || search_exhausted {
267            entry.is_optimal = true;
268            is_optimal = true;
269            self.status = if !stopped && self.statistics.skipping_violations > 0 {
270                SearchStatus::NotProven
271            } else if !stopped {
272                SearchStatus::Optimal
273            } else if !self.time_remains() {
274                SearchStatus::TimeLimit
275            } else {
276                SearchStatus::BudgetExhausted
277            };
278        }
279
280        if let Some(budget) = next_budget {
281            self.apply(budget);
282        }
283
284        self.statistics.error = entry.error;
285        self.statistics.cache_size = self.cache.len();
286        self.statistics.duration = self.elapsed_time();
287        is_optimal
288    }
289
290    /// The `(cost, label)` of the view as a single leaf.
291    fn leaf(&mut self, view: &DataView<'_>) -> (Cost, usize) {
292        shared::checked_leaf(
293            &self.loss,
294            view.get_labels_freqs(),
295            || view.instance_ids(),
296            &mut self.invalid_cost,
297        )
298    }
299
300    /// Looks a child subproblem up in the cache, seeding a new entry with the
301    /// cost of the view as a single leaf.
302    ///
303    /// `child_depth` is the child's remaining depth. Returns whether the entry
304    /// is new, its cache index, and a copy of it.
305    fn cache_child(&mut self, view: &DataView<'_>, child_depth: usize) -> (bool, usize, Entry) {
306        let (is_new, index) = self.cache.insert(&view.bitset, child_depth);
307        let depth = self.config.max_depth - child_depth;
308        let mut entry = Entry::default();
309
310        if let Some(slot) = self.cache.get_mut(index) {
311            if is_new {
312                let (error, label) = shared::checked_leaf(
313                    &self.loss,
314                    view.get_labels_freqs(),
315                    || view.instance_ids(),
316                    &mut self.invalid_cost,
317                );
318                slot.error = error;
319                slot.label = label;
320                slot.depth = depth;
321                slot.lower_bound = 0.0;
322                slot.is_valid = true;
323            }
324            entry = *slot;
325        }
326        (is_new, index, entry)
327    }
328
329    fn expand_node_with_view(
330        &mut self,
331        view: &DataView<'_>,
332        config: &SearchConfig,
333        current_best: &mut Entry,
334        parent_index: usize,
335        is_new: bool,
336        upper_bound: Cost,
337    ) -> bool {
338        if current_best.error == 0.0 || view.is_empty() {
339            if let Some(entry) = self.cache.get_mut(parent_index) {
340                entry.mark_exact();
341                *current_best = *entry;
342            }
343            return false;
344        }
345
346        // The budgets this visit runs under. The split budget only applies to
347        // the `first` selector and to the first pass of every selector.
348        let under = SearchedUnder {
349            discrepancy: config.budget.saturating_sub(config.discrepancy),
350            split_budget: if config.point_selector == PointSelector::First
351                || self.config.nb_runs <= 1
352            {
353                self.split_budget
354            } else {
355                usize::MAX
356            },
357            upper_bound,
358            truncated: false,
359        };
360
361        if !is_new {
362            if current_best.is_optimal {
363                self.statistics.cache_hits += 1;
364                return false;
365            }
366            if let Some(previous) = current_best.searched_under {
367                // A complete earlier search proved its lower bound. If that
368                // bound already reaches the parent's upper bound, no subtree
369                // here can help the parent.
370                if !previous.truncated && current_best.lower_bound >= upper_bound {
371                    self.statistics.cache_hits += 1;
372                    return false;
373                }
374                // An earlier search with at least this much budget and upper
375                // bound found everything this visit could, and was truncated
376                // exactly when this visit would be.
377                if previous.covers(&under, current_best.is_valid) {
378                    self.statistics.cache_hits += 1;
379                    return previous.truncated;
380                }
381            }
382        }
383
384        if config.max_depth == 0 {
385            current_best.is_leaf = true;
386            current_best.mark_exact();
387
388            if let Some(entry) = self.cache.get_mut(parent_index) {
389                entry.is_leaf = true;
390                entry.mark_exact();
391                *current_best = *entry;
392            }
393
394            return false;
395        }
396
397        if current_best.error <= config.max_gap || view.len() <= 1 {
398            return false;
399        }
400
401        // Only at depth exactly 2: `ConTreeDepth2` always builds two levels of
402        // tests, which would exceed a depth-1 request.
403        if config.fast_d2 && config.max_depth == 2 {
404            // Every pass revisits the same depth-2 subproblems, so each one is
405            // solved exactly once, without the parent's upper bound, and then
406            // reused. The parent still compares the result to its own bound.
407            let skipping = config.skipping.applies(&self.loss);
408            let tree = self.specialized.fit(
409                view,
410                config,
411                current_best,
412                Cost::INFINITY,
413                &mut self.statistics,
414                &mut Depth2Loss::new(&self.loss, skipping, &mut self.invalid_cost),
415            );
416            current_best.mark_exact();
417            let tree_index = self.cache.insert_tree(tree);
418            if let Some(entry) = self.cache.get_mut(parent_index) {
419                *entry = *current_best;
420                entry.tree_idx = Some(tree_index);
421            }
422            self.statistics.specialized_solver_call += 1;
423
424            return false;
425        }
426
427        let num_features = view.get_feature_number();
428        let heuristics_data = view.features_best_score();
429        debug_assert!(
430            num_features == heuristics_data.len(),
431            "Missmatch with heuristics and features number"
432        );
433        let mut stopped = false;
434        for (it, &(_, feat)) in heuristics_data.iter().take(num_features).enumerate() {
435            let feat_discrepancy = config.discrepancy + it;
436            if feat_discrepancy > config.budget {
437                // Features come in Gini order, so every later one is over the
438                // budget too. Fall through to record what this search covered.
439                self.pass_report.cut_by_discrepancy = true;
440                stopped = true;
441                break;
442            }
443
444            let mut node_config = *config;
445            node_config.discrepancy = feat_discrepancy;
446
447            stopped |= self.expand_on_feature(
448                view,
449                feat,
450                parent_index,
451                &node_config,
452                current_best,
453                upper_bound.min(current_best.error),
454            );
455
456            if current_best.error == 0.0 {
457                current_best.mark_exact();
458                if let Some(entry) = self.cache.get_mut(parent_index) {
459                    entry.mark_exact();
460                }
461
462                return false;
463            }
464
465            if !self.time_remains() {
466                return true;
467            }
468        }
469        // A result is exact only if the search was neither truncated by a
470        // budget nor cut short by the upper bound: a search completed under a
471        // tight bound says nothing about a looser one.
472        current_best.finalize_lower_bound(upper_bound);
473        current_best.searched_under = Some(SearchedUnder {
474            truncated: stopped,
475            ..under
476        });
477        if let Some(entry) = self.cache.get_mut(parent_index) {
478            entry.lower_bound = current_best.lower_bound;
479            entry.is_valid = current_best.is_valid;
480            entry.is_optimal = !stopped && current_best.is_valid;
481            entry.searched_under = current_best.searched_under;
482        }
483        current_best.is_optimal = !stopped && current_best.is_valid;
484
485        stopped
486    }
487
488    fn expand_on_feature(
489        &mut self,
490        view: &DataView<'_>,
491        feature_index: usize,
492        cache_index: usize,
493        config: &SearchConfig,
494        current_best: &mut Entry,
495        upper_bound: Cost,
496    ) -> bool {
497        let feature_column = view.get_sorted_feature(feature_index);
498        let feature_column_ids = view.get_feature_indices(feature_index);
499
500        if config.point_selector == PointSelector::First || self.config.nb_runs <= 1 {
501            let stopped = self.expand_on_feature_gini_priority(
502                view,
503                feature_index,
504                cache_index,
505                config,
506                current_best,
507                upper_bound,
508            );
509            return stopped;
510        }
511
512        let possible_index_split = view.get_possible_split_indices(feature_index);
513        if possible_index_split.is_empty() {
514            return false;
515        }
516
517        // Restrict the search to splits that meet the minimum support. The
518        // per-candidate check below skips a whole interval, so it cannot be
519        // relied on to find the feasible splits inside it.
520        let feasible =
521            support_feasible_splits(possible_index_split, view.len(), self.config.min_sup);
522        if feasible.is_empty() {
523            return false;
524        }
525
526        let mut pruner = IntervalsPruner::new(
527            possible_index_split,
528            config.max_gap,
529            config.min_sup,
530            config.skipping.applies(&self.loss),
531        );
532        let mut queue = VecDeque::new();
533        let init_bound = Bound::new(feasible.start, feasible.end - 1, None, None);
534        queue.push_back(init_bound);
535
536        let mut stopped = false;
537
538        while let Some(mut current_bound) = queue.pop_front() {
539            if !self.time_remains() {
540                return true;
541            }
542
543            // Prune against the tighter of the incumbent and the parent's upper
544            // bound: a split that cannot beat the parent's bound is of no use
545            // to the parent.
546            if pruner.subinterval_pruning(&current_bound, current_best.error.min(upper_bound)) {
547                continue;
548            }
549
550            pruner.interval_shrinking(&mut current_bound, current_best.error.min(upper_bound));
551            if !current_bound.is_valid() {
552                continue;
553            }
554
555            let selected_point = self.select_point(config, &current_bound);
556            let split_point = possible_index_split[selected_point];
557            let int_half_distance = split_point
558                .saturating_sub(possible_index_split[current_bound.left_bound])
559                .max(possible_index_split[current_bound.right_bound].saturating_sub(split_point))
560                as Cost;
561
562            let threshold_value = if selected_point > 0 {
563                let previous = feature_column_ids[possible_index_split[selected_point - 1]];
564                let point = feature_column_ids[split_point];
565                shared::threshold_between(
566                    feature_column[previous].value(),
567                    feature_column[point].value(),
568                )
569            } else {
570                let point = feature_column_ids[split_point];
571                shared::threshold_between(
572                    feature_column[feature_column_ids[0]].value(),
573                    feature_column[point].value(),
574                )
575            };
576
577            let (left_view, right_view) = view.split(feature_index, split_point);
578
579            if left_view.len() < self.config.min_sup || right_view.len() < self.config.min_sup {
580                continue;
581            }
582
583            // Search the larger child first: its error tightens the bound
584            // passed to the smaller one.
585            let process_left_first = left_view.len() >= right_view.len();
586
587            let left_config = config.derive_left();
588            let mut left_entry = Entry::default();
589            let mut right_entry = Entry::default();
590
591            // A cache index of 0 means "no child", kept by a side the search
592            // does not descend into.
593            let (mut left_index, mut right_index) = (0, 0);
594            let (mut left_is_new, mut right_is_new);
595
596            let larger_upper_bound = current_best.error.min(upper_bound + config.cost_step);
597            self.statistics.general_solver_call += 1;
598
599            if process_left_first {
600                (left_is_new, left_index, left_entry) =
601                    self.cache_child(&left_view, left_config.max_depth);
602
603                stopped |= self.expand_node_with_view(
604                    &left_view,
605                    &left_config,
606                    &mut left_entry,
607                    left_index,
608                    left_is_new,
609                    larger_upper_bound,
610                );
611                left_entry.finalize_lower_bound(larger_upper_bound);
612            } else {
613                (right_is_new, right_index, right_entry) =
614                    self.cache_child(&right_view, left_config.max_depth);
615
616                stopped |= self.expand_node_with_view(
617                    &right_view,
618                    &left_config,
619                    &mut right_entry,
620                    right_index,
621                    right_is_new,
622                    larger_upper_bound,
623                );
624                right_entry.finalize_lower_bound(larger_upper_bound);
625            }
626
627            let larger_error = if process_left_first {
628                left_entry.lower_bound
629            } else {
630                right_entry.lower_bound
631            };
632            // The smaller child's bound is what the larger child left of the
633            // incumbent, widened by the interval half-distance so that one
634            // search can prune the whole interval around this split. Only the
635            // interval pruner gains from it, so without threshold skipping the
636            // bound is not widened.
637            let budget = current_best.error.min(upper_bound + config.cost_step) - larger_error;
638            let widening = if config.skipping.applies(&self.loss) {
639                int_half_distance
640            } else {
641                0.0
642            };
643            let smaller_ub = budget + widening;
644            let smaller_upper_bound = smaller_ub.max(0.0);
645            // Stays `None` when the smaller child is not searched because the
646            // larger one alone already exceeds the upper bound.
647            let mut right_error: Option<Cost> = None;
648
649            // Search the smaller child unless the budget is negative. A budget
650            // of exactly zero is still explored, since a zero-error subtree
651            // may exist.
652            if smaller_ub > 0.0 || budget == 0.0 {
653                self.statistics.general_solver_call += 1;
654                let right_config = config.derive_right(left_config.max_gap);
655
656                if process_left_first {
657                    (right_is_new, right_index, right_entry) =
658                        self.cache_child(&right_view, right_config.max_depth);
659
660                    stopped |= self.expand_node_with_view(
661                        &right_view,
662                        &right_config,
663                        &mut right_entry,
664                        right_index,
665                        right_is_new,
666                        smaller_upper_bound,
667                    );
668                    right_entry.finalize_lower_bound(smaller_upper_bound);
669                } else {
670                    (left_is_new, left_index, left_entry) =
671                        self.cache_child(&left_view, right_config.max_depth);
672
673                    stopped |= self.expand_node_with_view(
674                        &left_view,
675                        &right_config,
676                        &mut left_entry,
677                        left_index,
678                        left_is_new,
679                        smaller_upper_bound,
680                    );
681                    left_entry.finalize_lower_bound(smaller_upper_bound);
682                }
683
684                // Use lower bounds rather than errors: when a bound cut a
685                // child search short, its `error` is only an upper bound.
686                right_error = Some(right_entry.lower_bound);
687
688                let feature_best = left_entry.lower_bound + right_entry.lower_bound;
689                if left_entry.is_valid && right_entry.is_valid && feature_best < current_best.error
690                {
691                    current_best.error = feature_best;
692                    if config.is_root {
693                        self.trajectory.push((self.elapsed_time(), feature_best));
694                    }
695                    current_best.feature = feature_index;
696                    current_best.split = threshold_value;
697                    current_best.left = left_index;
698                    current_best.right = right_index;
699
700                    if let Some(entry) = self.cache.get_mut(cache_index) {
701                        *entry = *current_best;
702                    }
703
704                    if feature_best == 0.0 {
705                        return false;
706                    }
707                }
708            }
709
710            // `None` for a child the bound invalidated: its lower bound may be
711            // zero without the subtree being error-free, and the pruner would
712            // read a zero as "nothing beyond this split can do better".
713            let left_score = left_entry.is_valid.then_some(left_entry.lower_bound);
714            pruner.add_result(selected_point, left_score, right_error);
715            if current_bound.left_bound == current_bound.right_bound {
716                continue;
717            }
718
719            let score_difference = (left_entry.lower_bound + right_error.unwrap_or(0.0)
720                - current_best.error.min(upper_bound))
721            .max(0.0);
722
723            let (left_bound, right_bound) = pruner.neighbourhood_pruning(
724                score_difference,
725                current_bound.left_bound,
726                current_bound.right_bound,
727                selected_point,
728            );
729
730            if left_bound <= current_bound.right_bound {
731                queue.push_back(Bound {
732                    left_bound,
733                    right_bound: current_bound.right_bound,
734                    last_split_left_index: Some(selected_point),
735                    last_split_right_index: current_bound.last_split_right_index,
736                });
737            }
738
739            if current_bound.left_bound <= right_bound {
740                queue.push_back(Bound {
741                    left_bound: current_bound.left_bound,
742                    right_bound,
743                    last_split_left_index: current_bound.last_split_left_index,
744                    last_split_right_index: Some(selected_point),
745                });
746            }
747        }
748        stopped
749    }
750
751    /// Tries the thresholds of one feature one at a time, in Gini order when
752    /// the heuristic is on and in position order otherwise, up to the split
753    /// budget. Returns `true` if the budget truncated the search.
754    fn expand_on_feature_gini_priority(
755        &mut self,
756        view: &DataView<'_>,
757        feature_index: usize,
758        cache_index: usize,
759        config: &SearchConfig,
760        current_best: &mut Entry,
761        upper_bound: Cost,
762    ) -> bool {
763        let feature_column = view.get_sorted_feature(feature_index);
764        let feature_column_ids = view.get_feature_indices(feature_index);
765
766        let possible_splits = view.get_possible_split_indices(feature_index);
767
768        if possible_splits.is_empty() {
769            return false;
770        }
771
772        let mut pruner = IntervalsPruner::new(
773            possible_splits,
774            config.max_gap,
775            config.min_sup,
776            config.skipping.applies(&self.loss),
777        );
778
779        // Thresholds already evaluated or ruled out by the pruner.
780        let mut pruned = vec![false; possible_splits.len()];
781        let mut stopped = false;
782
783        let local_split_budget = possible_splits.len().min(self.split_budget);
784
785        if config.use_heuristic {
786            let sorted_indices = view.ordered_possible_splits(feature_index);
787            for (idx, &split_idx) in sorted_indices.iter().enumerate() {
788                if self.config.nb_runs <= 1 && idx > 0 {
789                    self.pass_report.cut_by_split = true;
790                    return true;
791                }
792
793                if idx > local_split_budget - 1 {
794                    self.pass_report.cut_by_split = true;
795                    return true;
796                }
797
798                if !self.time_remains() {
799                    return true;
800                }
801
802                if pruned[split_idx] {
803                    continue;
804                }
805
806                stopped |= self.evaluate_split_gini_priority(
807                    view,
808                    feature_index,
809                    split_idx,
810                    possible_splits,
811                    feature_column,
812                    feature_column_ids,
813                    config,
814                    current_best,
815                    cache_index,
816                    upper_bound,
817                    &mut pruner,
818                    &mut pruned,
819                );
820
821                if current_best.error == 0.0 {
822                    return false;
823                }
824            }
825        } else {
826            for split_idx in 0..possible_splits.len() {
827                if self.config.nb_runs <= 1 && split_idx > 0 {
828                    self.pass_report.cut_by_split = true;
829                    return true;
830                }
831
832                if split_idx > local_split_budget - 1 {
833                    self.pass_report.cut_by_split = true;
834                    return true;
835                }
836
837                if !self.time_remains() {
838                    return true;
839                }
840
841                if pruned[split_idx] {
842                    continue;
843                }
844
845                stopped |= self.evaluate_split_gini_priority(
846                    view,
847                    feature_index,
848                    split_idx,
849                    possible_splits,
850                    feature_column,
851                    feature_column_ids,
852                    config,
853                    current_best,
854                    cache_index,
855                    upper_bound,
856                    &mut pruner,
857                    &mut pruned,
858                );
859
860                if current_best.error == 0.0 {
861                    return false;
862                }
863            }
864        }
865
866        if local_split_budget < possible_splits.len() {
867            self.pass_report.cut_by_split = true;
868            stopped = true;
869        }
870
871        stopped
872    }
873
874    #[inline]
875    #[allow(clippy::too_many_arguments)]
876    fn evaluate_split_gini_priority(
877        &mut self,
878        view: &DataView<'_>,
879        feature_index: usize,
880        split_idx: usize,
881        possible_splits: &[usize],
882        feature_column: &Feature,
883        feature_column_ids: &[usize],
884        config: &SearchConfig,
885        current_best: &mut Entry,
886        cache_index: usize,
887        upper_bound: Cost,
888        pruner: &mut IntervalsPruner<'_>,
889        pruned: &mut [bool],
890    ) -> bool {
891        // The interval around `split_idx` that no pruned threshold separates
892        // from it.
893        let mut current_left = split_idx;
894        while current_left > 0 && pruned[current_left - 1] {
895            current_left -= 1;
896        }
897        current_left = current_left.saturating_sub(1);
898
899        let mut current_right = split_idx;
900        while current_right < possible_splits.len() - 1 && pruned[current_right + 1] {
901            current_right += 1;
902        }
903        if current_right < possible_splits.len() - 1 {
904            current_right += 1;
905        }
906
907        let split_point = possible_splits[split_idx];
908
909        let threshold_value = if split_idx > 0 {
910            let previous = feature_column_ids[possible_splits[split_idx - 1]];
911            let point = feature_column_ids[split_point];
912            shared::threshold_between(
913                feature_column[previous].value(),
914                feature_column[point].value(),
915            )
916        } else {
917            let point = feature_column_ids[split_point];
918            shared::threshold_between(
919                feature_column[feature_column_ids[0]].value(),
920                feature_column[point].value(),
921            )
922        };
923
924        let (left_view, right_view) = view.split(feature_index, split_point);
925
926        if left_view.len() < self.config.min_sup || right_view.len() < self.config.min_sup {
927            pruned[split_idx] = true;
928            return false;
929        }
930
931        // Search the larger child first: its error tightens the bound passed
932        // to the smaller one.
933        let process_left_first = left_view.len() >= right_view.len();
934
935        let left_config = config.derive_left();
936        let mut left_entry = Entry::default();
937        let mut right_entry = Entry::default();
938
939        // A cache index of 0 means "no child", kept by a side the search does
940        // not descend into.
941        let (mut left_index, mut right_index) = (0, 0);
942        let (mut left_is_new, mut right_is_new);
943        let mut stopped = false;
944
945        let int_half_distance = split_point
946            .saturating_sub(possible_splits[0])
947            .max(possible_splits[possible_splits.len() - 1].saturating_sub(split_point))
948            as Cost;
949
950        let larger_upper_bound = current_best.error.min(upper_bound + config.cost_step);
951        self.statistics.general_solver_call += 1;
952
953        if process_left_first {
954            (left_is_new, left_index, left_entry) =
955                self.cache_child(&left_view, left_config.max_depth);
956
957            stopped |= self.expand_node_with_view(
958                &left_view,
959                &left_config,
960                &mut left_entry,
961                left_index,
962                left_is_new,
963                larger_upper_bound,
964            );
965            left_entry.finalize_lower_bound(larger_upper_bound);
966        } else {
967            (right_is_new, right_index, right_entry) =
968                self.cache_child(&right_view, left_config.max_depth);
969
970            stopped |= self.expand_node_with_view(
971                &right_view,
972                &left_config,
973                &mut right_entry,
974                right_index,
975                right_is_new,
976                larger_upper_bound,
977            );
978            right_entry.finalize_lower_bound(larger_upper_bound);
979        }
980
981        let larger_error = if process_left_first {
982            left_entry.lower_bound
983        } else {
984            right_entry.lower_bound
985        };
986        // The smaller child's bound is what the larger child left of the
987        // incumbent, widened by the interval half-distance so that one search
988        // can prune the whole interval around this split. Only the interval
989        // pruner gains from it, so without threshold skipping the bound is
990        // not widened.
991        let budget = current_best.error.min(upper_bound + config.cost_step) - larger_error;
992        let widening = if config.skipping.applies(&self.loss) {
993            int_half_distance
994        } else {
995            0.0
996        };
997        let smaller_ub = budget + widening;
998        let smaller_upper_bound = smaller_ub.max(0.0);
999        let mut right_error: Option<Cost> = None;
1000
1001        // Search the smaller child unless the budget is negative. A budget of
1002        // exactly zero is still explored, since a zero-error subtree may exist.
1003        if smaller_ub > 0.0 || budget == 0.0 {
1004            self.statistics.general_solver_call += 1;
1005            let right_config = config.derive_right(left_config.max_gap);
1006
1007            if process_left_first {
1008                (right_is_new, right_index, right_entry) =
1009                    self.cache_child(&right_view, right_config.max_depth);
1010
1011                stopped |= self.expand_node_with_view(
1012                    &right_view,
1013                    &right_config,
1014                    &mut right_entry,
1015                    right_index,
1016                    right_is_new,
1017                    smaller_upper_bound,
1018                );
1019                right_entry.finalize_lower_bound(smaller_upper_bound);
1020            } else {
1021                (left_is_new, left_index, left_entry) =
1022                    self.cache_child(&left_view, right_config.max_depth);
1023
1024                stopped |= self.expand_node_with_view(
1025                    &left_view,
1026                    &right_config,
1027                    &mut left_entry,
1028                    left_index,
1029                    left_is_new,
1030                    smaller_upper_bound,
1031                );
1032                left_entry.finalize_lower_bound(smaller_upper_bound);
1033            }
1034
1035            // Use lower bounds rather than errors: when a bound cut a child
1036            // search short, its `error` is only an upper bound.
1037            right_error = Some(right_entry.lower_bound);
1038
1039            let feature_best = left_entry.lower_bound + right_entry.lower_bound;
1040            if left_entry.is_valid && right_entry.is_valid && feature_best < current_best.error {
1041                current_best.error = feature_best;
1042                if config.is_root {
1043                    self.trajectory.push((self.elapsed_time(), feature_best));
1044                }
1045                current_best.feature = feature_index;
1046                current_best.split = threshold_value;
1047                current_best.left = left_index;
1048                current_best.right = right_index;
1049
1050                if let Some(entry) = self.cache.get_mut(cache_index) {
1051                    *entry = *current_best;
1052                }
1053
1054                if feature_best == 0.0 {
1055                    return false;
1056                }
1057            }
1058        }
1059
1060        // `None` for a child the bound invalidated: its lower bound may be
1061        // zero without the subtree being error-free, and the pruner would read
1062        // a zero as "nothing beyond this split can do better".
1063        let left_score = left_entry.is_valid.then_some(left_entry.lower_bound);
1064        pruner.add_result(split_idx, left_score, right_error);
1065        pruned[split_idx] = true;
1066
1067        let score_difference = (left_entry.lower_bound + right_error.unwrap_or(0.0)
1068            - current_best.error.min(upper_bound))
1069        .max(0.0);
1070        let (new_left_bound, new_right_bound) =
1071            pruner.neighbourhood_pruning(score_difference, current_left, current_right, split_idx);
1072
1073        // `neighbourhood_pruning` leaves two intervals to search,
1074        // `[current_left, new_right_bound]` and `[new_left_bound, current_right]`.
1075        // Only the thresholds strictly between them can be skipped.
1076        let skip_from = new_right_bound.saturating_add(1).max(current_left);
1077        if skip_from < split_idx {
1078            pruned[skip_from..split_idx].fill(true);
1079        }
1080        let skip_to = new_left_bound.min(current_right + 1);
1081        if split_idx + 1 < skip_to {
1082            pruned[split_idx + 1..skip_to].fill(true);
1083        }
1084        stopped
1085    }
1086
1087    fn discrepancy_limit(nb_candidates: usize, remaining_depth: usize) -> usize {
1088        let mut max_discrepancy = nb_candidates;
1089        for i in 1..remaining_depth {
1090            max_discrepancy += nb_candidates.saturating_sub(i);
1091        }
1092
1093        max_discrepancy
1094    }
1095
1096    pub fn statistics(&self) -> &Statistics {
1097        &self.statistics
1098    }
1099
1100    fn time_remains(&self) -> bool {
1101        shared::time_remains(&self.runtime, self.config.max_time)
1102    }
1103
1104    fn elapsed_time(&self) -> f64 {
1105        shared::elapsed_time(&self.runtime)
1106    }
1107
1108    pub fn select_point(&mut self, config: &SearchConfig, bound: &Bound) -> usize {
1109        shared::select_point(config.point_selector, &mut self.rng, bound)
1110    }
1111
1112    /// The tree found by the passes run so far.
1113    pub fn get_solution_tree(&mut self) -> Tree {
1114        shared::build_solution_tree(&self.cache)
1115    }
1116}
1117
1118#[cfg(test)]
1119mod contree_lds_test {
1120    use crate::algorithms::{ConTree, ConTreeLds};
1121    use crate::common::{Cost, PointSelector};
1122    use crate::reader::data_reader::DataReader;
1123    use crate::reader::DataReaderError;
1124    use crate::tests::fixture;
1125
1126    #[test]
1127    fn lds_is_never_better_than_the_exhaustive_optimum() -> Result<(), DataReaderError> {
1128        let reader = DataReader::default();
1129        let mut dataset = reader.read_file(&fixture("hepatitis.txt"))?;
1130        dataset.sort_features();
1131
1132        let mut exhaustive: ConTree = ConTree::new(
1133            1,
1134            2,
1135            100.0,
1136            Cost::INFINITY,
1137            PointSelector::Mid,
1138            0.0,
1139            false,
1140            true,
1141        );
1142        exhaustive.fit(&dataset).expect("fit");
1143
1144        let mut lds: ConTreeLds = ConTreeLds::new(
1145            1,
1146            2,
1147            100.0,
1148            Cost::INFINITY,
1149            PointSelector::Mid,
1150            0.0,
1151            false,
1152            true,
1153        );
1154        lds.fit(&dataset).expect("fit");
1155
1156        // LDS restricts the search, so it can only tie or lose against the
1157        // exhaustive optimum. Beating it would mean the exhaustive search is
1158        // pruning something it should not.
1159        assert!(
1160            lds.statistics.error >= exhaustive.statistics().error,
1161            "LDS reported {} against an exhaustive optimum of {}",
1162            lds.statistics.error,
1163            exhaustive.statistics().error
1164        );
1165        assert!(lds.cache.len() > 0, "the cache was never written to");
1166
1167        Ok(())
1168    }
1169
1170    #[test]
1171    fn the_anytime_loop_terminates_without_leaning_on_the_time_limit() {
1172        // With no time limit, the search must still stop once the budget
1173        // schedule runs out.
1174        let reader = DataReader::default();
1175        let mut dataset = reader.read_file(&fixture("small.txt")).unwrap();
1176        dataset.sort_features();
1177
1178        let mut lds: ConTreeLds = ConTreeLds::new(
1179            1,
1180            3,
1181            f64::INFINITY,
1182            Cost::INFINITY,
1183            PointSelector::Mid,
1184            0.0,
1185            false,
1186            false,
1187        );
1188        lds.fit(&dataset).expect("fit");
1189
1190        assert!(lds.statistics.error <= dataset.count() as Cost);
1191    }
1192}