criterion_stats/bivariate/
mod.rs

1//! Bivariate analysis
2
3mod bootstrap;
4mod resamples;
5
6pub mod regression;
7
8use std::cmp;
9
10use float::Float;
11use num_cpus;
12use thread_scoped as thread;
13
14use bivariate::resamples::Resamples;
15use tuple::{Tuple, TupledDistributionsBuilder};
16use univariate::Sample;
17
18/// Bivariate `(X, Y)` data
19///
20/// Invariants:
21///
22/// - No `NaN`s in the data
23/// - At least two data points in the set
24pub struct Data<'a, X, Y>(&'a [X], &'a [Y])
25where
26    X: 'a,
27    Y: 'a;
28
29impl<'a, X, Y> Copy for Data<'a, X, Y> {}
30
31#[cfg_attr(feature = "cargo-clippy", allow(clippy::expl_impl_clone_on_copy))]
32impl<'a, X, Y> Clone for Data<'a, X, Y> {
33    fn clone(&self) -> Data<'a, X, Y> {
34        *self
35    }
36}
37
38impl<'a, X, Y> Data<'a, X, Y> {
39    /// Returns the length of the data set
40    pub fn len(&self) -> usize {
41        self.0.len()
42    }
43
44    /// Checks whether the data set is empty
45    pub fn is_empty(&self) -> bool {
46        self.len() == 0
47    }
48
49    /// Iterate over the data set
50    pub fn iter(&self) -> Pairs<'a, X, Y> {
51        Pairs {
52            data: *self,
53            state: 0,
54        }
55    }
56}
57
58impl<'a, X, Y> Data<'a, X, Y>
59where
60    X: Float,
61    Y: Float,
62{
63    /// Creates a new data set from two existing slices
64    pub fn new(xs: &'a [X], ys: &'a [Y]) -> Data<'a, X, Y> {
65        assert!(
66            xs.len() == ys.len()
67                && xs.len() > 1
68                && xs.iter().all(|x| !x.is_nan())
69                && ys.iter().all(|y| !y.is_nan())
70        );
71
72        Data(xs, ys)
73    }
74
75    // TODO Remove the `T` parameter in favor of `S::Output`
76    /// Returns the bootstrap distributions of the parameters estimated by the `statistic`
77    ///
78    /// - Multi-threaded
79    /// - Time: `O(nresamples)`
80    /// - Memory: `O(nresamples)`
81    pub fn bootstrap<T, S>(&self, nresamples: usize, statistic: S) -> T::Distributions
82    where
83        S: Fn(Data<X, Y>) -> T,
84        S: Sync,
85        T: Tuple,
86        T::Distributions: Send,
87        T::Builder: Send,
88    {
89        let ncpus = num_cpus::get();
90
91        unsafe {
92            // TODO need some sensible threshold to trigger the multi-threaded path
93            if ncpus > 1 && nresamples > self.0.len() {
94                let granularity = nresamples / ncpus + 1;
95                let statistic = &statistic;
96
97                let chunks = (0..ncpus)
98                    .map(|i| {
99                        let mut sub_distributions: T::Builder =
100                            TupledDistributionsBuilder::new(granularity);
101                        let mut resamples = Resamples::new(*self);
102                        let offset = i * granularity;
103
104                        thread::scoped(move || {
105                            for _ in offset..cmp::min(offset + granularity, nresamples) {
106                                sub_distributions.push(statistic(resamples.next()))
107                            }
108                            sub_distributions
109                        })
110                    })
111                    .collect::<Vec<_>>();
112
113                let mut builder: T::Builder = TupledDistributionsBuilder::new(nresamples);
114                for chunk in chunks {
115                    builder.extend(&mut (chunk.join()));
116                }
117                builder.complete()
118            } else {
119                let mut distributions: T::Builder = TupledDistributionsBuilder::new(nresamples);
120                let mut resamples = Resamples::new(*self);
121
122                for _ in 0..nresamples {
123                    distributions.push(statistic(resamples.next()));
124                }
125
126                distributions.complete()
127            }
128        }
129    }
130
131    /// Returns a view into the `X` data
132    pub fn x(&self) -> &'a Sample<X> {
133        Sample::new(&self.0)
134    }
135
136    /// Returns a view into the `Y` data
137    pub fn y(&self) -> &'a Sample<Y> {
138        Sample::new(&self.1)
139    }
140}
141
142/// Iterator over `Data`
143pub struct Pairs<'a, X: 'a, Y: 'a> {
144    data: Data<'a, X, Y>,
145    state: usize,
146}
147
148impl<'a, X, Y> Iterator for Pairs<'a, X, Y> {
149    type Item = (&'a X, &'a Y);
150
151    fn next(&mut self) -> Option<(&'a X, &'a Y)> {
152        if self.state < self.data.len() {
153            let i = self.state;
154            self.state += 1;
155
156            // This is safe because i will always be < self.data.{0,1}.len()
157            debug_assert!(i < self.data.0.len());
158            debug_assert!(i < self.data.1.len());
159            unsafe { Some((self.data.0.get_unchecked(i), self.data.1.get_unchecked(i))) }
160        } else {
161            None
162        }
163    }
164}