1use crate::matrix::FdMatrix;
7use crate::maybe_par_chunks_mut_enumerate;
8use rand::prelude::*;
9use rand_distr::StandardNormal;
10
11pub mod band;
12pub mod dispatch;
13pub mod erl;
14pub mod extremal;
15pub mod fraiman_muniz;
16pub mod half_region;
17pub mod hypo_epi;
18pub mod linf;
19pub mod modal;
20pub mod random_projection;
21pub mod random_tukey;
22pub mod rpd;
23pub mod spatial;
24pub mod tvd;
25
26#[cfg(test)]
27mod tests;
28
29pub use band::{band_1d, modified_band_1d, modified_epigraph_index_1d};
31pub use dispatch::{functional_boxplot, functional_depth, DepthMethod, FunctionalBoxplotResult};
32pub use erl::extreme_rank_length_depth_1d;
33pub use extremal::extremal_depth_1d;
34#[allow(deprecated)]
37pub use fraiman_muniz::{fraiman_muniz, fraiman_muniz_1d, fraiman_muniz_2d};
38pub use half_region::{half_region_depth_1d, modified_half_region_depth_1d};
39pub use hypo_epi::{epigraph_index_1d, hypograph_index_1d, modified_hypograph_index_1d};
40pub use linf::linfinity_depth_1d;
41#[allow(deprecated)]
42pub use modal::{modal, modal_1d, modal_2d};
43#[allow(deprecated)]
44pub use random_projection::{
45 random_projection, random_projection_1d, random_projection_1d_seeded, random_projection_2d,
46};
47#[allow(deprecated)]
48pub use random_tukey::{random_tukey, random_tukey_1d, random_tukey_1d_seeded, random_tukey_2d};
49pub use rpd::{rpd_depth_1d, rpd_depth_1d_seeded};
50pub use spatial::{
51 functional_spatial_1d, functional_spatial_2d, kernel_functional_spatial_1d,
52 kernel_functional_spatial_2d,
53};
54pub use tvd::{total_variation_depth_1d, TvdMssResult};
55
56pub(super) fn generate_random_projections(nproj: usize, m: usize, seed: Option<u64>) -> Vec<f64> {
67 let mut rng: Box<dyn RngCore> = match seed {
68 Some(s) => Box::new(StdRng::seed_from_u64(s)),
69 None => Box::new(rand::thread_rng()),
70 };
71 let mut projections = vec![0.0; nproj * m];
72 for p_idx in 0..nproj {
73 let base = p_idx * m;
74 let mut norm_sq = 0.0;
75 for t in 0..m {
76 let v: f64 = rng.sample(StandardNormal);
77 projections[base + t] = v;
78 norm_sq += v * v;
79 }
80 let inv_norm = 1.0 / norm_sq.sqrt();
81 for t in 0..m {
82 projections[base + t] *= inv_norm;
83 }
84 }
85 projections
86}
87
88pub(super) fn project_and_sort_reference(
93 data_ori: &FdMatrix,
94 projections: &[f64],
95 nproj: usize,
96 nori: usize,
97 m: usize,
98) -> Vec<f64> {
99 let mut sorted = vec![0.0; nproj * nori];
100 maybe_par_chunks_mut_enumerate!(sorted, nori, |(p_idx, spo): (usize, &mut [f64])| {
101 let proj = &projections[p_idx * m..(p_idx + 1) * m];
102 for j in 0..nori {
103 let mut dot = 0.0;
104 for t in 0..m {
105 dot += data_ori[(j, t)] * proj[t];
106 }
107 spo[j] = dot;
108 }
109 spo.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
110 });
111 sorted
112}
113
114pub(super) fn random_depth_core(
121 data_obj: &FdMatrix,
122 data_ori: &FdMatrix,
123 nproj: usize,
124 seed: Option<u64>,
125 init: f64,
126 aggregate: impl Fn(f64, f64) -> f64 + Sync,
127 finalize: impl Fn(f64, usize) -> f64 + Sync,
128) -> Vec<f64> {
129 use crate::iter_maybe_parallel;
130 #[cfg(feature = "parallel")]
131 use rayon::iter::ParallelIterator;
132
133 let nobj = data_obj.nrows();
134 let nori = data_ori.nrows();
135 let m = data_obj.ncols();
136
137 if nobj == 0 || nori == 0 || m == 0 || nproj == 0 {
138 return Vec::new();
139 }
140
141 let projections = generate_random_projections(nproj, m, seed);
142 let sorted_proj_ori = project_and_sort_reference(data_ori, &projections, nproj, nori, m);
143 let denom = nori as f64 + 1.0;
144
145 iter_maybe_parallel!(0..nobj)
146 .map(|i| {
147 let mut acc = init;
148 for p_idx in 0..nproj {
149 let proj = &projections[p_idx * m..(p_idx + 1) * m];
150 let sorted_ori = &sorted_proj_ori[p_idx * nori..(p_idx + 1) * nori];
151
152 let mut proj_i = 0.0;
153 for t in 0..m {
154 proj_i += data_obj[(i, t)] * proj[t];
155 }
156
157 let below = sorted_ori.partition_point(|&v| v < proj_i);
158 let above = nori - sorted_ori.partition_point(|&v| v <= proj_i);
159 let depth = (below.min(above) as f64 + 1.0) / denom;
160 acc = aggregate(acc, depth);
161 }
162 finalize(acc, nproj)
163 })
164 .collect()
165}