flow_peak_detection/
lib.rs

1//! KDE-based peak finding for flow cytometry histograms.
2
3use anyhow::{Result, bail};
4
5/// Configuration for peak isolation.
6#[derive(Debug, Clone, Copy)]
7pub struct PeakConfig {
8    /// Fraction of max density required to count as a peak (0–1).
9    pub threshold: f64,
10    /// Upper/lower bias within the peak (0.5 = centered 50%). `1.0` keeps the full peak.
11    pub peak_bias: f64,
12    /// Minimum events retained after isolation.
13    pub min_events: usize,
14    /// KDE grid resolution.
15    pub resolution: usize,
16}
17
18impl Default for PeakConfig {
19    fn default() -> Self {
20        Self {
21            threshold: 0.3,
22            peak_bias: 1.0,
23            min_events: 100,
24            resolution: 512,
25        }
26    }
27}
28
29/// Result of isolating a peak region on a 1-D intensity sample.
30#[derive(Debug, Clone)]
31pub struct PeakResult {
32    pub range: (f64, f64),
33    pub median: f64,
34    pub event_indices: Vec<usize>,
35    pub density: f64,
36    pub combined_score: f64,
37}
38
39/// Gaussian KDE peaks on a uniform grid (Silverman bandwidth when `bandwidth` is `None`).
40pub fn detect_peaks_kde(
41    data: &[f64],
42    bandwidth: Option<f64>,
43    resolution: usize,
44    threshold: f64,
45) -> Vec<(f64, f64)> {
46    if data.is_empty() || resolution < 8 {
47        return Vec::new();
48    }
49    let finite: Vec<f64> = data.iter().copied().filter(|v| v.is_finite()).collect();
50    if finite.len() < 2 {
51        return Vec::new();
52    }
53    let min = finite.iter().copied().fold(f64::INFINITY, f64::min);
54    let max = finite.iter().copied().fold(f64::NEG_INFINITY, f64::max);
55    if !(max > min) {
56        return Vec::new();
57    }
58    let n = finite.len() as f64;
59    let mean = finite.iter().sum::<f64>() / n;
60    let var = finite.iter().map(|x| (x - mean).powi(2)).sum::<f64>() / n;
61    let std = var.sqrt().max(f64::EPSILON);
62    let bw = bandwidth
63        .unwrap_or_else(|| 1.06 * std * n.powf(-0.2))
64        .max(f64::EPSILON);
65
66    let mut xs = Vec::with_capacity(resolution);
67    let mut dens = vec![0.0_f64; resolution];
68    for i in 0..resolution {
69        let x = min + (max - min) * (i as f64) / ((resolution - 1) as f64);
70        xs.push(x);
71        let mut d = 0.0;
72        for &v in &finite {
73            let z = (x - v) / bw;
74            d += (-0.5 * z * z).exp();
75        }
76        dens[i] = d / (n * bw * (std::f64::consts::TAU).sqrt());
77    }
78    let max_d = dens.iter().copied().fold(0.0_f64, f64::max);
79    if max_d <= 0.0 {
80        return Vec::new();
81    }
82    let cut = threshold.clamp(0.0, 1.0) * max_d;
83    let mut peaks = Vec::new();
84    for i in 1..resolution.saturating_sub(1) {
85        if dens[i] >= cut && dens[i] >= dens[i - 1] && dens[i] >= dens[i + 1] {
86            peaks.push((xs[i], dens[i]));
87        }
88    }
89    peaks
90}
91
92fn median_of(sorted: &[f64]) -> f64 {
93    if sorted.is_empty() {
94        return f64::NAN;
95    }
96    let mid = sorted.len() / 2;
97    if sorted.len() % 2 == 1 {
98        sorted[mid]
99    } else {
100        0.5 * (sorted[mid - 1] + sorted[mid])
101    }
102}
103
104#[derive(Clone, Copy)]
105enum PeakSide {
106    Positive,
107    Negative,
108}
109
110fn isolate_peak(data: &[f64], config: &PeakConfig, side: PeakSide) -> Result<PeakResult> {
111    let peaks = detect_peaks_kde(data, None, config.resolution, config.threshold);
112    if peaks.is_empty() {
113        bail!("no peaks detected above threshold {}", config.threshold);
114    }
115    let (peak_x, peak_d) = match side {
116        PeakSide::Positive => peaks
117            .iter()
118            .copied()
119            .max_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal))
120            .unwrap_or(peaks[0]),
121        PeakSide::Negative => peaks
122            .iter()
123            .copied()
124            .min_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal))
125            .unwrap_or(peaks[0]),
126    };
127
128    let finite: Vec<(usize, f64)> = data
129        .iter()
130        .enumerate()
131        .filter_map(|(i, &v)| v.is_finite().then_some((i, v)))
132        .collect();
133    if finite.is_empty() {
134        bail!("no finite events");
135    }
136    let min = finite.iter().map(|(_, v)| *v).fold(f64::INFINITY, f64::min);
137    let max = finite
138        .iter()
139        .map(|(_, v)| *v)
140        .fold(f64::NEG_INFINITY, f64::max);
141    let half = ((max - min) * 0.25).max(f64::EPSILON);
142    let mut in_window: Vec<(usize, f64)> = finite
143        .into_iter()
144        .filter(|(_, v)| (*v - peak_x).abs() <= half)
145        .collect();
146    if in_window.is_empty() {
147        bail!("no events near peak at {peak_x}");
148    }
149    in_window.sort_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));
150    let bias = config.peak_bias.clamp(0.05, 1.0);
151    let keep = ((in_window.len() as f64) * bias).ceil() as usize;
152    let keep = keep.max(1).min(in_window.len());
153    let trimmed = match side {
154        PeakSide::Positive => &in_window[in_window.len() - keep..],
155        PeakSide::Negative => &in_window[..keep],
156    };
157    if trimmed.len() < config.min_events {
158        bail!(
159            "only {} events in peak (minimum {})",
160            trimmed.len(),
161            config.min_events
162        );
163    }
164    let values: Vec<f64> = trimmed.iter().map(|(_, v)| *v).collect();
165    let lo = values[0];
166    let hi = *values.last().unwrap_or(&lo);
167    Ok(PeakResult {
168        range: (lo, hi),
169        median: median_of(&values),
170        event_indices: trimmed.iter().map(|(i, _)| *i).collect(),
171        density: peak_d,
172        combined_score: peak_d * (trimmed.len() as f64).ln_1p(),
173    })
174}
175
176/// Isolate the brightest dense peak (positive population).
177pub fn isolate_positive_peak(data: &[f64], config: &PeakConfig) -> Result<PeakResult> {
178    isolate_peak(data, config, PeakSide::Positive)
179}
180
181/// Isolate the leftmost peak (negative population).
182pub fn isolate_negative_peak(data: &[f64], config: &PeakConfig) -> Result<PeakResult> {
183    isolate_peak(data, config, PeakSide::Negative)
184}
185
186/// Boolean mask over `data` for the positive peak.
187pub fn isolate_positive_peak_mask(
188    data: &[f64],
189    threshold: f64,
190    peak_bias: f64,
191) -> Result<Vec<bool>> {
192    let config = PeakConfig {
193        threshold,
194        peak_bias,
195        min_events: 1,
196        ..PeakConfig::default()
197    };
198    let peak = isolate_positive_peak(data, &config)?;
199    let mut mask = vec![false; data.len()];
200    for i in peak.event_indices {
201        if i < mask.len() {
202            mask[i] = true;
203        }
204    }
205    Ok(mask)
206}
207
208#[cfg(test)]
209mod tests {
210    use super::*;
211
212    #[test]
213    fn detects_bimodal_and_isolates_right_peak() {
214        let mut data = Vec::new();
215        for _ in 0..400 {
216            data.push(1.0);
217        }
218        for _ in 0..400 {
219            data.push(10.0);
220        }
221        for (i, v) in data.iter_mut().enumerate() {
222            *v += (i % 7) as f64 * 0.01;
223        }
224        let peaks = detect_peaks_kde(&data, None, 256, 0.2);
225        assert!(!peaks.is_empty());
226        let pos = isolate_positive_peak(
227            &data,
228            &PeakConfig {
229                min_events: 50,
230                ..PeakConfig::default()
231            },
232        )
233        .unwrap();
234        assert!(pos.median > 5.0);
235    }
236}
flow_peak_detection/lib.rs

flow_peak_detection/
lib.rs
