fdars_core/depth/
mod.rs

1//! Depth measures for functional data.
2//!
3//! This module provides various depth measures for assessing the centrality
4//! of functional observations within a reference sample.
5
6use 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
29// Re-export all public functions
30pub use band::{band, modified_band, modified_epigraph_index};
31pub use dispatch::{functional_boxplot, functional_depth, DepthMethod, FunctionalBoxplotResult};
32pub use erl::extreme_rank_length_depth;
33pub use extremal::extremal_depth;
34pub use fraiman_muniz::fraiman_muniz;
35pub use half_region::{half_region_depth, modified_half_region_depth};
36pub use hypo_epi::{epigraph_index, hypograph_index, modified_hypograph_index};
37pub use linf::linfinity_depth;
38pub use modal::modal;
39pub use random_projection::{random_projection, random_projection_1d_seeded};
40pub use random_tukey::{random_tukey, random_tukey_1d_seeded};
41pub use rpd::{rpd_depth, rpd_depth_1d_seeded};
42pub use spatial::{functional_spatial, kernel_functional_spatial};
43pub use tvd::{total_variation_depth, TvdMssResult};
44
45// Crate-internal re-exports of the underlying `_1d` primitives. These remain
46// callable within the crate (notably `dispatch.rs` and inter-depth callers)
47// after the public surface consolidated onto the `Dim`-dispatched entry points.
48pub(crate) use band::{band_1d, modified_band_1d, modified_epigraph_index_1d};
49pub(crate) use erl::extreme_rank_length_depth_1d;
50pub(crate) use extremal::extremal_depth_1d;
51pub(crate) use fraiman_muniz::fraiman_muniz_1d;
52pub(crate) use half_region::{half_region_depth_1d, modified_half_region_depth_1d};
53pub(crate) use hypo_epi::{epigraph_index_1d, hypograph_index_1d, modified_hypograph_index_1d};
54pub(crate) use linf::linfinity_depth_1d;
55pub(crate) use tvd::total_variation_depth_1d;
56// Test-only: the Cat 1 primitives whose sole in-crate caller (besides their own
57// dispatcher) is the inline `#[cfg(test)]` bit-identity module.
58#[cfg(test)]
59pub(crate) use modal::modal_1d;
60#[cfg(test)]
61pub(crate) use random_projection::random_projection_1d;
62#[cfg(test)]
63pub(crate) use random_tukey::random_tukey_1d;
64#[cfg(test)]
65pub(crate) use rpd::rpd_depth_1d;
66
67// ---------------------------------------------------------------------------
68// Shared helpers
69// ---------------------------------------------------------------------------
70
71/// Generate `nproj` unit-norm random projection vectors of dimension `m`.
72///
73/// Returns a flat buffer of length `nproj * m` where projection `p` occupies
74/// `[p*m .. (p+1)*m]`.
75///
76/// If `seed` is Some, uses a deterministic RNG seeded from the given value.
77pub(super) fn generate_random_projections(nproj: usize, m: usize, seed: Option<u64>) -> Vec<f64> {
78    let mut rng: Box<dyn RngCore> = match seed {
79        Some(s) => Box::new(StdRng::seed_from_u64(s)),
80        None => Box::new(rand::thread_rng()),
81    };
82    let mut projections = vec![0.0; nproj * m];
83    for p_idx in 0..nproj {
84        let base = p_idx * m;
85        let mut norm_sq = 0.0;
86        for t in 0..m {
87            let v: f64 = rng.sample(StandardNormal);
88            projections[base + t] = v;
89            norm_sq += v * v;
90        }
91        let inv_norm = 1.0 / norm_sq.sqrt();
92        for t in 0..m {
93            projections[base + t] *= inv_norm;
94        }
95    }
96    projections
97}
98
99/// Project each reference curve onto each projection direction and sort.
100///
101/// Returns a flat buffer of length `nproj * nori` where the sorted projections
102/// for direction `p` occupy `[p*nori .. (p+1)*nori]`.
103pub(super) fn project_and_sort_reference(
104    data_ori: &FdMatrix,
105    projections: &[f64],
106    nproj: usize,
107    nori: usize,
108    m: usize,
109) -> Vec<f64> {
110    let mut sorted = vec![0.0; nproj * nori];
111    maybe_par_chunks_mut_enumerate!(sorted, nori, |(p_idx, spo): (usize, &mut [f64])| {
112        let proj = &projections[p_idx * m..(p_idx + 1) * m];
113        for j in 0..nori {
114            let mut dot = 0.0;
115            for t in 0..m {
116                dot += data_ori[(j, t)] * proj[t];
117            }
118            spo[j] = dot;
119        }
120        spo.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
121    });
122    sorted
123}
124
125/// Shared implementation for random projection-based depth measures.
126///
127/// Generates `nproj` random projections, projects both object and reference
128/// curves to scalars, and computes univariate depth via binary-search ranking.
129/// The `aggregate` and `finalize` closures control how per-projection depths
130/// are combined (e.g. average for RP depth, minimum for Tukey depth).
131pub(super) fn random_depth_core(
132    data_obj: &FdMatrix,
133    data_ori: &FdMatrix,
134    nproj: usize,
135    seed: Option<u64>,
136    init: f64,
137    aggregate: impl Fn(f64, f64) -> f64 + Sync,
138    finalize: impl Fn(f64, usize) -> f64 + Sync,
139) -> Vec<f64> {
140    use crate::iter_maybe_parallel;
141    #[cfg(feature = "parallel")]
142    use rayon::iter::ParallelIterator;
143
144    let nobj = data_obj.nrows();
145    let nori = data_ori.nrows();
146    let m = data_obj.ncols();
147
148    if nobj == 0 || nori == 0 || m == 0 || nproj == 0 {
149        return Vec::new();
150    }
151
152    let projections = generate_random_projections(nproj, m, seed);
153    let sorted_proj_ori = project_and_sort_reference(data_ori, &projections, nproj, nori, m);
154    let denom = nori as f64 + 1.0;
155
156    iter_maybe_parallel!(0..nobj)
157        .map(|i| {
158            let mut acc = init;
159            for p_idx in 0..nproj {
160                let proj = &projections[p_idx * m..(p_idx + 1) * m];
161                let sorted_ori = &sorted_proj_ori[p_idx * nori..(p_idx + 1) * nori];
162
163                let mut proj_i = 0.0;
164                for t in 0..m {
165                    proj_i += data_obj[(i, t)] * proj[t];
166                }
167
168                let below = sorted_ori.partition_point(|&v| v < proj_i);
169                let above = nori - sorted_ori.partition_point(|&v| v <= proj_i);
170                let depth = (below.min(above) as f64 + 1.0) / denom;
171                acc = aggregate(acc, depth);
172            }
173            finalize(acc, nproj)
174        })
175        .collect()
176}
fdars_core/depth/mod.rs

fdars_core/depth/
mod.rs
