1use crate::matrix::FdMatrix;
7use crate::maybe_par_chunks_mut_enumerate;
8use rand::prelude::*;
9use rand_distr::StandardNormal;
10
11pub mod band;
12pub mod fraiman_muniz;
13pub mod modal;
14pub mod random_projection;
15pub mod random_tukey;
16pub mod spatial;
17
18#[cfg(test)]
19mod tests;
20
21pub use band::{band_1d, modified_band_1d, modified_epigraph_index_1d};
23pub use fraiman_muniz::{fraiman_muniz_1d, fraiman_muniz_2d};
24pub use modal::{modal_1d, modal_2d};
25pub use random_projection::{
26 random_projection_1d, random_projection_1d_seeded, random_projection_2d,
27};
28pub use random_tukey::{random_tukey_1d, random_tukey_1d_seeded, random_tukey_2d};
29pub use spatial::{
30 functional_spatial_1d, functional_spatial_2d, kernel_functional_spatial_1d,
31 kernel_functional_spatial_2d,
32};
33
34pub(super) fn generate_random_projections(nproj: usize, m: usize, seed: Option<u64>) -> Vec<f64> {
45 let mut rng: Box<dyn RngCore> = match seed {
46 Some(s) => Box::new(StdRng::seed_from_u64(s)),
47 None => Box::new(rand::thread_rng()),
48 };
49 let mut projections = vec![0.0; nproj * m];
50 for p_idx in 0..nproj {
51 let base = p_idx * m;
52 let mut norm_sq = 0.0;
53 for t in 0..m {
54 let v: f64 = rng.sample(StandardNormal);
55 projections[base + t] = v;
56 norm_sq += v * v;
57 }
58 let inv_norm = 1.0 / norm_sq.sqrt();
59 for t in 0..m {
60 projections[base + t] *= inv_norm;
61 }
62 }
63 projections
64}
65
66pub(super) fn project_and_sort_reference(
71 data_ori: &FdMatrix,
72 projections: &[f64],
73 nproj: usize,
74 nori: usize,
75 m: usize,
76) -> Vec<f64> {
77 let mut sorted = vec![0.0; nproj * nori];
78 maybe_par_chunks_mut_enumerate!(sorted, nori, |(p_idx, spo): (usize, &mut [f64])| {
79 let proj = &projections[p_idx * m..(p_idx + 1) * m];
80 for j in 0..nori {
81 let mut dot = 0.0;
82 for t in 0..m {
83 dot += data_ori[(j, t)] * proj[t];
84 }
85 spo[j] = dot;
86 }
87 spo.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
88 });
89 sorted
90}
91
92pub(super) fn random_depth_core(
99 data_obj: &FdMatrix,
100 data_ori: &FdMatrix,
101 nproj: usize,
102 seed: Option<u64>,
103 init: f64,
104 aggregate: impl Fn(f64, f64) -> f64 + Sync,
105 finalize: impl Fn(f64, usize) -> f64 + Sync,
106) -> Vec<f64> {
107 use crate::iter_maybe_parallel;
108 #[cfg(feature = "parallel")]
109 use rayon::iter::ParallelIterator;
110
111 let nobj = data_obj.nrows();
112 let nori = data_ori.nrows();
113 let m = data_obj.ncols();
114
115 if nobj == 0 || nori == 0 || m == 0 || nproj == 0 {
116 return Vec::new();
117 }
118
119 let projections = generate_random_projections(nproj, m, seed);
120 let sorted_proj_ori = project_and_sort_reference(data_ori, &projections, nproj, nori, m);
121 let denom = nori as f64 + 1.0;
122
123 iter_maybe_parallel!(0..nobj)
124 .map(|i| {
125 let mut acc = init;
126 for p_idx in 0..nproj {
127 let proj = &projections[p_idx * m..(p_idx + 1) * m];
128 let sorted_ori = &sorted_proj_ori[p_idx * nori..(p_idx + 1) * nori];
129
130 let mut proj_i = 0.0;
131 for t in 0..m {
132 proj_i += data_obj[(i, t)] * proj[t];
133 }
134
135 let below = sorted_ori.partition_point(|&v| v < proj_i);
136 let above = nori - sorted_ori.partition_point(|&v| v <= proj_i);
137 let depth = (below.min(above) as f64 + 1.0) / denom;
138 acc = aggregate(acc, depth);
139 }
140 finalize(acc, nproj)
141 })
142 .collect()
143}