dml_core/
crossfit.rs

1use ndarray::{Array1, ArrayView1, ArrayView2, Axis};
2use rand::seq::SliceRandom;
3use rand::SeedableRng;
4use rand_chacha::ChaCha8Rng;
5
6use crate::data::gather;
7use crate::error::{DmlError, Result};
8use crate::learner::{Classifier, Regressor};
9
10/// One train/test split.
11#[derive(Debug, Clone, PartialEq)]
12pub struct Fold {
13    pub train: Vec<usize>,
14    pub test: Vec<usize>,
15}
16
17/// Shuffled K-fold splitter with a fixed seed, so results are reproducible.
18#[derive(Debug, Clone)]
19pub struct KFold {
20    n_folds: usize,
21    seed: u64,
22}
23
24impl KFold {
25    pub fn new(n_folds: usize, seed: u64) -> Result<Self> {
26        if n_folds < 2 {
27            return Err(DmlError::InvalidArgument(
28                "n_folds must be at least 2".into(),
29            ));
30        }
31        Ok(Self { n_folds, seed })
32    }
33
34    pub fn n_folds(&self) -> usize {
35        self.n_folds
36    }
37
38    /// Split `0..n` into `n_folds` disjoint test sets covering every index once.
39    pub fn split(&self, n: usize) -> Result<Vec<Fold>> {
40        if n < self.n_folds {
41            return Err(DmlError::InvalidArgument(format!(
42                "{n} observations is fewer than {} folds",
43                self.n_folds
44            )));
45        }
46        let mut idx: Vec<usize> = (0..n).collect();
47        idx.shuffle(&mut ChaCha8Rng::seed_from_u64(self.seed));
48
49        let base = n / self.n_folds;
50        let extra = n % self.n_folds;
51        let mut folds = Vec::with_capacity(self.n_folds);
52        let mut start = 0;
53        for k in 0..self.n_folds {
54            let size = base + usize::from(k < extra);
55            let mut test = idx[start..start + size].to_vec();
56            let mut train: Vec<usize> = idx[..start]
57                .iter()
58                .chain(&idx[start + size..])
59                .copied()
60                .collect();
61            test.sort_unstable();
62            train.sort_unstable();
63            folds.push(Fold { train, test });
64            start += size;
65        }
66        Ok(folds)
67    }
68
69    /// Stratified split: the rows of each stratum (e.g. treated / control) are
70    /// shuffled and dealt round-robin over the folds, so every fold holds about
71    /// the same share of each stratum. Fails if a stratum has fewer rows than
72    /// there are folds.
73    pub fn split_stratified(&self, strata: &[usize]) -> Result<Vec<Fold>> {
74        let n = strata.len();
75        if n < self.n_folds {
76            return Err(DmlError::InvalidArgument(format!(
77                "{n} observations is fewer than {} folds",
78                self.n_folds
79            )));
80        }
81        let mut groups: std::collections::BTreeMap<usize, Vec<usize>> = Default::default();
82        for (i, &s) in strata.iter().enumerate() {
83            groups.entry(s).or_default().push(i);
84        }
85        let mut rng = ChaCha8Rng::seed_from_u64(self.seed);
86        let mut tests: Vec<Vec<usize>> = vec![Vec::new(); self.n_folds];
87        let mut next = 0;
88        for (stratum, mut idx) in groups {
89            if idx.len() < self.n_folds {
90                return Err(DmlError::InvalidArgument(format!(
91                    "stratum {stratum} has {} observations, fewer than {} folds",
92                    idx.len(),
93                    self.n_folds
94                )));
95            }
96            idx.shuffle(&mut rng);
97            for i in idx {
98                tests[next].push(i);
99                next = (next + 1) % self.n_folds;
100            }
101        }
102        Ok(tests
103            .into_iter()
104            .map(|mut test| {
105                test.sort_unstable();
106                let mut in_test = vec![false; n];
107                for &i in &test {
108                    in_test[i] = true;
109                }
110                let train = (0..n).filter(|&i| !in_test[i]).collect();
111                Fold { train, test }
112            })
113            .collect())
114    }
115}
116
117/// How observations are split for cross-fitting: `n_rep` independent
118/// repetitions of seeded K-fold splitting.
119#[derive(Debug, Clone, Copy, PartialEq, Eq)]
120pub struct SampleSplitting {
121    pub n_folds: usize,
122    pub n_rep: usize,
123    pub seed: u64,
124}
125
126impl Default for SampleSplitting {
127    fn default() -> Self {
128        Self {
129            n_folds: 5,
130            n_rep: 1,
131            seed: 0,
132        }
133    }
134}
135
136impl SampleSplitting {
137    /// Draw the folds for every repetition. Repetition `r` uses its own seed
138    /// derived from `seed`, so repetition 0 equals `KFold::new(n_folds, seed)`.
139    pub fn draw(&self, n: usize) -> Result<Vec<Vec<Fold>>> {
140        if self.n_rep == 0 {
141            return Err(DmlError::InvalidArgument("n_rep must be at least 1".into()));
142        }
143        (0..self.n_rep as u64)
144            .map(|r| {
145                // Golden-ratio increment keeps repetition seeds far apart.
146                let seed = self
147                    .seed
148                    .wrapping_add(r.wrapping_mul(0x9E37_79B9_7F4A_7C15));
149                KFold::new(self.n_folds, seed)?.split(n)
150            })
151            .collect()
152    }
153
154    /// Cluster-level splits for every repetition (see [`crate::cluster`]).
155    pub fn draw_clustered(&self, clusters: &[Array1<i64>]) -> Result<Vec<Vec<crate::ClusterFold>>> {
156        if self.n_rep == 0 {
157            return Err(DmlError::InvalidArgument("n_rep must be at least 1".into()));
158        }
159        (0..self.n_rep as u64)
160            .map(|r| {
161                let seed = self
162                    .seed
163                    .wrapping_add(r.wrapping_mul(0x9E37_79B9_7F4A_7C15));
164                crate::cluster::split_clusters(clusters, self.n_folds, seed)
165            })
166            .collect()
167    }
168}
169
170/// Out-of-fold predictions: for each fold, fit a fresh copy of `learner` on the
171/// training rows and predict the held-out rows. Every row is predicted by a
172/// model that never saw it — the core of DML cross-fitting.
173pub fn cross_fit_predict(
174    learner: &dyn Regressor,
175    x: ArrayView2<'_, f64>,
176    target: ArrayView1<'_, f64>,
177    folds: &[Fold],
178) -> Result<Array1<f64>> {
179    if target.len() != x.nrows() {
180        return Err(DmlError::DimensionMismatch(
181            "cross_fit_predict x/target rows".into(),
182        ));
183    }
184    cross_fit_with(x.nrows(), folds, |fold| {
185        let mut model = learner.fresh();
186        model.fit(
187            x.select(Axis(0), &fold.train).view(),
188            gather(target, &fold.train).view(),
189        )?;
190        model.predict(x.select(Axis(0), &fold.test).view())
191    })
192}
193
194/// Like [`cross_fit_predict`], but each model is trained only on the training
195/// rows where `keep[i]` is true (e.g. controls only); predictions still cover
196/// every test row.
197pub fn cross_fit_predict_subset(
198    learner: &dyn Regressor,
199    x: ArrayView2<'_, f64>,
200    target: ArrayView1<'_, f64>,
201    folds: &[Fold],
202    keep: &[bool],
203) -> Result<Array1<f64>> {
204    cross_fit_predict(learner, x, target, &subset_folds(folds, keep, x.nrows())?)
205}
206
207/// Like [`cross_fit_predict_proba`], but each model is trained only on the
208/// training rows where `keep[i]` is true; predictions cover every test row.
209pub fn cross_fit_predict_proba_subset(
210    learner: &dyn Classifier,
211    x: ArrayView2<'_, f64>,
212    labels: ArrayView1<'_, f64>,
213    folds: &[Fold],
214    keep: &[bool],
215) -> Result<Array1<f64>> {
216    cross_fit_predict_proba(learner, x, labels, &subset_folds(folds, keep, x.nrows())?)
217}
218
219/// Restrict every training set to the rows where `keep` is true.
220fn subset_folds(folds: &[Fold], keep: &[bool], n: usize) -> Result<Vec<Fold>> {
221    if keep.len() != n {
222        return Err(DmlError::DimensionMismatch("subgroup mask length".into()));
223    }
224    let subset: Vec<Fold> = folds
225        .iter()
226        .map(|f| Fold {
227            train: f.train.iter().copied().filter(|&i| keep[i]).collect(),
228            test: f.test.clone(),
229        })
230        .collect();
231    if subset.iter().any(|f| f.train.is_empty()) {
232        return Err(DmlError::InvalidArgument(
233            "a training fold has no observations in the requested subgroup".into(),
234        ));
235    }
236    Ok(subset)
237}
238
239/// Out-of-fold class-1 probabilities from a [`Classifier`].
240pub fn cross_fit_predict_proba(
241    learner: &dyn Classifier,
242    x: ArrayView2<'_, f64>,
243    labels: ArrayView1<'_, f64>,
244    folds: &[Fold],
245) -> Result<Array1<f64>> {
246    if labels.len() != x.nrows() {
247        return Err(DmlError::DimensionMismatch(
248            "cross_fit_predict_proba x/labels rows".into(),
249        ));
250    }
251    cross_fit_with(x.nrows(), folds, |fold| {
252        let mut model = learner.fresh();
253        model.fit(
254            x.select(Axis(0), &fold.train).view(),
255            gather(labels, &fold.train).view(),
256        )?;
257        model.predict_proba(x.select(Axis(0), &fold.test).view())
258    })
259}
260
261/// Run `fit_predict` on every fold and scatter each fold's test-set
262/// predictions into one length-`n` vector.
263///
264/// With the `parallel` feature (on by default) folds run on rayon's thread
265/// pool; results are identical to the sequential order because every fold
266/// writes to disjoint rows.
267pub(crate) fn cross_fit_with<F>(n: usize, folds: &[Fold], fit_predict: F) -> Result<Array1<f64>>
268where
269    F: Fn(&Fold) -> Result<Array1<f64>> + Sync,
270{
271    #[cfg(feature = "parallel")]
272    let preds: Vec<Array1<f64>> = {
273        use rayon::prelude::*;
274        folds.par_iter().map(&fit_predict).collect::<Result<_>>()?
275    };
276    #[cfg(not(feature = "parallel"))]
277    let preds: Vec<Array1<f64>> = folds.iter().map(&fit_predict).collect::<Result<_>>()?;
278
279    let mut out = Array1::<f64>::from_elem(n, f64::NAN);
280    for (fold, pred) in folds.iter().zip(preds) {
281        if pred.len() != fold.test.len() {
282            return Err(DmlError::DimensionMismatch(
283                "learner returned wrong number of predictions".into(),
284            ));
285        }
286        for (&i, &p) in fold.test.iter().zip(pred.iter()) {
287            out[i] = p;
288        }
289    }
290    if out.iter().any(|v| v.is_nan()) {
291        return Err(DmlError::InvalidArgument(
292            "folds do not cover every observation, or a learner returned NaN".into(),
293        ));
294    }
295    Ok(out)
296}
297
298#[cfg(test)]
299mod tests {
300    use super::*;
301    use crate::learner::Ols;
302    use ndarray::Array2;
303
304    #[test]
305    fn folds_partition_all_indices() {
306        let folds = KFold::new(3, 7).unwrap().split(10).unwrap();
307        let mut seen: Vec<usize> = folds.iter().flat_map(|f| f.test.clone()).collect();
308        seen.sort_unstable();
309        assert_eq!(seen, (0..10).collect::<Vec<_>>());
310        for f in &folds {
311            assert!(f.test.iter().all(|i| !f.train.contains(i)));
312            assert_eq!(f.train.len() + f.test.len(), 10);
313        }
314        let sizes: Vec<usize> = folds.iter().map(|f| f.test.len()).collect();
315        assert_eq!(sizes, vec![4, 3, 3]);
316    }
317
318    #[test]
319    fn same_seed_same_split() {
320        let a = KFold::new(5, 42).unwrap().split(50).unwrap();
321        let b = KFold::new(5, 42).unwrap().split(50).unwrap();
322        let c = KFold::new(5, 43).unwrap().split(50).unwrap();
323        assert_eq!(a, b);
324        assert_ne!(a, c);
325    }
326
327    #[test]
328    fn repetitions_differ_and_first_matches_kfold() {
329        let reps = SampleSplitting {
330            n_folds: 4,
331            n_rep: 3,
332            seed: 9,
333        }
334        .draw(40)
335        .unwrap();
336        assert_eq!(reps.len(), 3);
337        assert_eq!(reps[0], KFold::new(4, 9).unwrap().split(40).unwrap());
338        assert_ne!(reps[0], reps[1]);
339        assert_ne!(reps[1], reps[2]);
340    }
341
342    #[test]
343    fn subset_training_uses_only_kept_rows() {
344        // Target is 0 on kept rows and 100 elsewhere: a mean model trained
345        // only on kept rows predicts ~0 everywhere.
346        let x = Array2::<f64>::zeros((20, 0));
347        let keep: Vec<bool> = (0..20).map(|i| i % 2 == 0).collect();
348        let y: Array1<f64> = keep.iter().map(|&k| if k { 0.0 } else { 100.0 }).collect();
349        let folds = KFold::new(4, 1).unwrap().split(20).unwrap();
350        let pred =
351            cross_fit_predict_subset(&Ols::new(), x.view(), y.view(), &folds, &keep).unwrap();
352        assert!(pred.iter().all(|p| p.abs() < 1e-12));
353    }
354
355    #[test]
356    fn stratified_folds_balance_each_stratum() {
357        let strata: Vec<usize> = (0..53).map(|i| usize::from(i % 5 == 0)).collect(); // 11 treated
358        let folds = KFold::new(5, 3).unwrap().split_stratified(&strata).unwrap();
359        let mut seen: Vec<usize> = folds.iter().flat_map(|f| f.test.clone()).collect();
360        seen.sort_unstable();
361        assert_eq!(seen, (0..53).collect::<Vec<_>>());
362        for f in &folds {
363            let treated = f.test.iter().filter(|&&i| strata[i] == 1).count();
364            assert!((2..=3).contains(&treated), "treated per fold: {treated}");
365            assert_eq!(f.train.len() + f.test.len(), 53);
366        }
367        assert!(KFold::new(5, 0)
368            .unwrap()
369            .split_stratified(&[0, 0, 0, 0, 0, 1])
370            .is_err());
371    }
372
373    #[test]
374    fn rejects_bad_settings() {
375        assert!(KFold::new(1, 0).is_err());
376        assert!(KFold::new(5, 0).unwrap().split(3).is_err());
377    }
378
379    #[test]
380    fn cross_fit_recovers_linear_signal() {
381        let n = 100;
382        let x = Array2::from_shape_fn((n, 2), |(i, j)| ((i * 7 + j * 3) % 11) as f64);
383        let y: Array1<f64> = x
384            .rows()
385            .into_iter()
386            .map(|r| 1.0 + 0.5 * r[0] - 2.0 * r[1])
387            .collect();
388        let folds = KFold::new(5, 1).unwrap().split(n).unwrap();
389        let pred = cross_fit_predict(&Ols::new(), x.view(), y.view(), &folds).unwrap();
390        assert!((&pred - &y).iter().all(|e| e.abs() < 1e-8));
391    }
392}
dml_core/crossfit.rs

dml_core/
crossfit.rs
