Skip to main content

fdars_core/streaming_depth/
bd.rs

1//! Streaming Band Depth (BD) estimator.
2
3use crate::iter_maybe_parallel;
4use crate::matrix::FdMatrix;
5#[cfg(feature = "parallel")]
6use rayon::iter::ParallelIterator;
7
8use super::sorted_ref::SortedReferenceState;
9use super::{c2, StreamingDepth};
10
11/// Full reference state that keeps per-curve values alongside sorted columns.
12///
13/// Required by Band Depth (BD), which checks all-or-nothing containment across
14/// ALL time points and therefore cannot decompose into per-point rank queries.
15#[derive(Debug, Clone, PartialEq)]
16pub struct FullReferenceState {
17    /// Sorted columns for rank queries (shared with MBD/FM estimators if desired).
18    pub sorted: SortedReferenceState,
19    /// `values_by_curve[j][t]` = reference curve j at time point t (row layout).
20    values_by_curve: Vec<Vec<f64>>,
21}
22
23impl FullReferenceState {
24    /// Build from a column-major reference matrix.
25    #[must_use = "expensive computation whose result should not be discarded"]
26    pub fn from_reference(data_ori: &FdMatrix) -> Self {
27        let nori = data_ori.nrows();
28        let n_points = data_ori.ncols();
29        let sorted = SortedReferenceState::from_reference(data_ori);
30        let values_by_curve: Vec<Vec<f64>> = (0..nori)
31            .map(|j| (0..n_points).map(|t| data_ori[(j, t)]).collect())
32            .collect();
33        Self {
34            sorted,
35            values_by_curve,
36        }
37    }
38}
39
40/// Streaming Band Depth estimator.
41///
42/// BD requires all-or-nothing containment across ALL time points -- it does not
43/// decompose per-point like MBD. The streaming advantage here is **reference
44/// decoupling** (no re-parsing the matrix) and **early-exit per pair** (break
45/// on first time point where x is outside the band), not an asymptotic
46/// improvement.
47#[derive(Debug, Clone, PartialEq)]
48pub struct StreamingBd {
49    state: FullReferenceState,
50}
51
52impl StreamingBd {
53    pub fn new(state: FullReferenceState) -> Self {
54        Self { state }
55    }
56
57    #[inline]
58    fn bd_one_inner(&self, curve: &[f64]) -> f64 {
59        let n = self.state.sorted.nori;
60        if n < 2 {
61            return 0.0;
62        }
63        let n_pairs = c2(n);
64        let n_points = self.state.sorted.n_points;
65
66        let mut count_in_band = 0usize;
67        for j in 0..n {
68            for k in (j + 1)..n {
69                let mut inside = true;
70                for t in 0..n_points {
71                    let x_t = curve[t];
72                    let y_j_t = self.state.values_by_curve[j][t];
73                    let y_k_t = self.state.values_by_curve[k][t];
74                    let band_min = y_j_t.min(y_k_t);
75                    let band_max = y_j_t.max(y_k_t);
76                    if x_t < band_min || x_t > band_max {
77                        inside = false;
78                        break;
79                    }
80                }
81                if inside {
82                    count_in_band += 1;
83                }
84            }
85        }
86        count_in_band as f64 / n_pairs as f64
87    }
88
89    /// Compute BD for row `row` of `data` without allocating a temporary Vec.
90    #[inline]
91    fn bd_one_from_row(&self, data: &FdMatrix, row: usize) -> f64 {
92        let n = self.state.sorted.nori;
93        if n < 2 {
94            return 0.0;
95        }
96        let n_pairs = c2(n);
97        let n_points = self.state.sorted.n_points;
98
99        let mut count_in_band = 0usize;
100        for j in 0..n {
101            for k in (j + 1)..n {
102                let mut inside = true;
103                for t in 0..n_points {
104                    let x_t = data[(row, t)];
105                    let y_j_t = self.state.values_by_curve[j][t];
106                    let y_k_t = self.state.values_by_curve[k][t];
107                    let band_min = y_j_t.min(y_k_t);
108                    let band_max = y_j_t.max(y_k_t);
109                    if x_t < band_min || x_t > band_max {
110                        inside = false;
111                        break;
112                    }
113                }
114                if inside {
115                    count_in_band += 1;
116                }
117            }
118        }
119        count_in_band as f64 / n_pairs as f64
120    }
121}
122
123impl StreamingDepth for StreamingBd {
124    fn depth_one(&self, curve: &[f64]) -> f64 {
125        self.bd_one_inner(curve)
126    }
127
128    fn depth_batch(&self, data_obj: &FdMatrix) -> Vec<f64> {
129        let nobj = data_obj.nrows();
130        let n = self.state.sorted.nori;
131        if nobj == 0 || self.state.sorted.n_points == 0 || n < 2 {
132            return vec![0.0; nobj];
133        }
134        iter_maybe_parallel!(0..nobj)
135            .map(|i| self.bd_one_from_row(data_obj, i))
136            .collect()
137    }
138
139    fn n_points(&self) -> usize {
140        self.state.sorted.n_points
141    }
142
143    fn n_reference(&self) -> usize {
144        self.state.sorted.nori
145    }
146}