mscore/algorithm/

utility.rs

use std::collections::HashMap;
use std::f64::consts::SQRT_2;
use rayon::prelude::*;
use rayon::ThreadPoolBuilder;

use std::collections::VecDeque;

fn gauss_kronrod(f: &dyn Fn(f64) -> f64, a: f64, b: f64) -> (f64, f64) {
    let nodes = [
        0.0, 0.20778495500789848, 0.40584515137739717, 0.58608723546769113,
        0.74153118559939444, 0.86486442335976907, 0.94910791234275852, 0.99145537112081264,
    ];
    let weights_gauss = [
        0.41795918367346939, 0.38183005050511894, 0.27970539148927667, 0.12948496616886969,
    ];
    let weights_kronrod = [
        0.20948214108472783, 0.20443294007529889, 0.19035057806478541, 0.16900472663926790,
        0.14065325971552592, 0.10479001032225018, 0.06309209262997855, 0.02293532201052922,
    ];

    let c1 = (b - a) / 2.0;
    let c2 = (b + a) / 2.0;

    let mut integral_gauss = 0.0;
    let mut integral_kronrod = 0.0;

    for i in 0..4 {
        let x = c1 * nodes[i] + c2;
        integral_gauss += weights_gauss[i] * (f(x) + f(2.0 * c2 - x));
    }

    for i in 0..8 {
        let x = c1 * nodes[i] + c2;
        integral_kronrod += weights_kronrod[i] * (f(x) + f(2.0 * c2 - x));
    }

    integral_gauss *= c1;
    integral_kronrod *= c1;

    (integral_kronrod, (integral_kronrod - integral_gauss).abs())
}

pub fn adaptive_integration(f: &dyn Fn(f64) -> f64, a: f64, b: f64, epsabs: f64, epsrel: f64) -> (f64, f64) {
    let mut intervals = VecDeque::new();
    intervals.push_back((a, b));

    let mut result = 0.0;
    let mut total_error = 0.0;

    while let Some((a, b)) = intervals.pop_front() {
        let (integral, error) = gauss_kronrod(f, a, b);
        if error < epsabs || error < epsrel * integral.abs() {
            result += integral;
            total_error += error;
        } else {
            let mid = (a + b) / 2.0;
            intervals.push_back((a, mid));
            intervals.push_back((mid, b));
        }
    }

    (result, total_error)
}




// Numerical integration using the trapezoidal rule
fn integrate<F>(f: F, a: f64, b: f64, n: usize) -> f64
    where
        F: Fn(f64) -> f64,
{
    let dx = (b - a) / n as f64;
    let mut sum = 0.0;
    for i in 0..n {
        let x = a + i as f64 * dx;
        sum += f(x);
    }
    sum * dx
}

// Complementary error function (erfc)
fn erfc(x: f64) -> f64 {
    1.0 - erf(x)
}

// Error function (erf)
fn erf(x: f64) -> f64 {
    let t = 1.0 / (1.0 + 0.5 * x.abs());
    let tau = t * (-x * x - 1.26551223 + t * (1.00002368 +
        t * (0.37409196 + t * (0.09678418 + t * (-0.18628806 +
            t * (0.27886807 + t * (-1.13520398 + t * (1.48851587 +
                t * (-0.82215223 + t * 0.17087277)))))))))
        .exp();
    if x >= 0.0 {
        1.0 - tau
    } else {
        tau - 1.0
    }
}

// Exponentially modified Gaussian function
fn emg(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
    let part1 = lambda / 2.0 * (-lambda * (x - mu) + lambda * lambda * sigma * sigma / 2.0).exp();
    let part2 = erfc((mu + lambda * sigma * sigma - x) / (sigma * 2.0_f64.sqrt()));
    part1 * part2
}

pub fn custom_cdf_normal(x: f64, mean: f64, std_dev: f64) -> f64 {
    let z = (x - mean) / std_dev;
    0.5 * (1.0 + erf(z / SQRT_2))
}

pub fn accumulated_intensity_cdf_normal(sample_start: f64, sample_end: f64, mean: f64, std_dev: f64) -> f64 {
    let cdf_start = custom_cdf_normal(sample_start, mean, std_dev);
    let cdf_end = custom_cdf_normal(sample_end, mean, std_dev);
    cdf_end - cdf_start
}

pub fn calculate_bounds_normal(mean: f64, std: f64, z_score: f64) -> (f64, f64) {
    (mean - z_score * std, mean + z_score * std)
}

pub fn emg_function(x: f64, mu: f64, sigma: f64, lambda: f64) -> f64 {
    let prefactor = lambda / 2.0 * ((lambda / 2.0) * (2.0 * mu + lambda * sigma.powi(2) - 2.0 * x)).exp();
    let erfc_part = erfc((mu + lambda * sigma.powi(2) - x) / (SQRT_2 * sigma));
    prefactor * erfc_part
}

pub fn emg_cdf_range(lower_limit: f64, upper_limit: f64, mu: f64, sigma: f64, lambda: f64, n_steps: Option<usize>) -> f64 {
    let n_steps = n_steps.unwrap_or(1000);
    integrate(|x| emg(x, mu, sigma, lambda), lower_limit, upper_limit, n_steps)
}

pub fn calculate_bounds_emg(mu: f64, sigma: f64, lambda: f64, step_size: f64, target: f64, lower_start: f64, upper_start: f64, n_steps: Option<usize>) -> (f64, f64) {
    assert!(0.0 <= target && target <= 1.0, "target must be in [0, 1]");

    // Degenerate parameters must not reach the allocation below. An infinite sigma makes
    // `(upper - lower) / step_size` infinite, `as usize` saturates it to usize::MAX, and the
    // `collect` then aborts the whole process with a capacity overflow rather than returning.
    // Hand back non-finite bounds and let the caller decide what an unusable peak means.
    if !mu.is_finite() || !sigma.is_finite() || !lambda.is_finite() || !(step_size > 0.0) {
        return (f64::NAN, f64::NAN);
    }

    let lower_initial = mu - lower_start * sigma - 2.0;
    let upper_initial = mu + upper_start * sigma;

    let span = upper_initial - lower_initial;
    if !span.is_finite() || span < 0.0 {
        return (f64::NAN, f64::NAN);
    }
    // A finite but absurd span would still allocate unboundedly: 1e8 points is already 800 MB.
    // Anything past that is a nonsensical elution window, not a slow one.
    const MAX_STEPS: f64 = 1e8;
    let steps_f = (span / step_size).round();
    if !steps_f.is_finite() || steps_f > MAX_STEPS {
        return (f64::NAN, f64::NAN);
    }

    let steps = steps_f as usize;
    let search_space: Vec<f64> = (0..=steps).map(|i| lower_initial + i as f64 * step_size).collect();

    let calc_cdf = |low: usize, high: usize| -> f64 {
        emg_cdf_range(search_space[low], search_space[high], mu, sigma, lambda, n_steps)
    };

    // Binary search for cutoff values
    let (mut low, mut high) = (0, steps);
    while low < high {
        let mid = low + (high - low) / 2;
        if calc_cdf(0, mid) < target {
            low = mid + 1;
        } else {
            high = mid;
        }
    }
    let upper_cutoff_index = low;

    low = 0;
    high = upper_cutoff_index;
    while low < high {
        let mid = high - (high - low) / 2;
        let prob_mid_to_upper = calc_cdf(mid, upper_cutoff_index);

        if prob_mid_to_upper < target {
            high = mid - 1;
        } else {
            low = mid;
        }
    }
    let lower_cutoff_index = high;

    (search_space[lower_cutoff_index], search_space[upper_cutoff_index])
}

pub fn calculate_frame_occurrence_emg(retention_times: &[f64], rt: f64, sigma: f64, lambda_: f64, target_p: f64, step_size: f64, n_steps: Option<usize>) -> Vec<i32> {
    let (rt_min, rt_max) = calculate_bounds_emg(rt, sigma, lambda_, step_size, target_p, 20.0, 60.0, n_steps);

    // A degenerate EMG parameter upstream can make these bounds non-finite. `partial_cmp` returns
    // None for NaN, so the comparators below used to unwrap a None and panic a worker thread
    // mid-simulation. A peptide with no finite elution window occupies no frames; say that instead
    // of aborting the run.
    if !rt_min.is_finite() || !rt_max.is_finite() {
        return Vec::new();
    }

    // Closest-frame search. `total_cmp` rather than `partial_cmp().unwrap()`: the inputs are finite
    // by the guard above, but this can never panic regardless of what reaches it.
    let first_frame = retention_times.iter()
        .enumerate()
        .min_by(|(_, &a), (_, &b)| (a - rt_min).abs().total_cmp(&(b - rt_min).abs()))
        .map(|(idx, _)| idx + 1) // Rust is zero-indexed, so +1 to match Python's 1-indexing
        .unwrap_or(0); // Fallback in case of an empty slice

    // Finding the frame closest to rt_max
    let last_frame = retention_times.iter()
        .enumerate()
        .min_by(|(_, &a), (_, &b)| (a - rt_max).abs().total_cmp(&(b - rt_max).abs()))
        .map(|(idx, _)| idx + 1) // Same adjustment for 1-indexing
        .unwrap_or(0); // Fallback

    // Generating the range of frames
    (first_frame..=last_frame).map(|x| x as i32).collect()
}

pub fn calculate_frame_abundance_emg(time_map: &HashMap<i32, f64>, occurrences: &[i32], rt: f64, sigma: f64, lambda_: f64, rt_cycle_length: f64, n_steps: Option<usize>) -> Vec<f64> {
    let mut frame_abundance = Vec::new();

    for &occurrence in occurrences {
        if let Some(&time) = time_map.get(&occurrence) {
            let start = time - rt_cycle_length;
            let i = emg_cdf_range(start, time, rt, sigma, lambda_, n_steps);
            frame_abundance.push(i);
        }
    }

    frame_abundance
}

// retention_times: &[f64], rt: f64, sigma: f64, lambda_: f64
pub fn calculate_frame_occurrences_emg_par(retention_times: &[f64], rts: Vec<f64>, sigmas: Vec<f64>, lambdas: Vec<f64>, target_p: f64, step_size: f64, num_threads: usize, n_steps: Option<usize>) -> Vec<Vec<i32>> {
    let thread_pool = ThreadPoolBuilder::new().num_threads(num_threads).build().unwrap();
    let result = thread_pool.install(|| {
        rts.into_par_iter().zip(sigmas.into_par_iter()).zip(lambdas.into_par_iter())
            .map(|((rt, sigma), lambda)| {
                calculate_frame_occurrence_emg(retention_times, rt, sigma, lambda, target_p, step_size, n_steps)
            })
            .collect()
    });
    result
}

pub fn calculate_frame_abundances_emg_par(time_map: &HashMap<i32, f64>, occurrences: Vec<Vec<i32>>, rts: Vec<f64>, sigmas: Vec<f64>, lambdas: Vec<f64>, rt_cycle_length: f64, num_threads: usize, n_steps: Option<usize>) -> Vec<Vec<f64>> {
    let thread_pool = ThreadPoolBuilder::new().num_threads(num_threads).build().unwrap();
    let result = thread_pool.install(|| {
        occurrences.into_par_iter().zip(rts.into_par_iter()).zip(sigmas.into_par_iter()).zip(lambdas.into_par_iter())
            .map(|(((occurrences, rt), sigma), lambda)| {
                calculate_frame_abundance_emg(time_map, &occurrences, rt, sigma, lambda, rt_cycle_length, n_steps)
            })
            .collect()
    });
    result
}

/// Project an EMG retention-time profile onto an explicit event timeline
/// (instrument-dispatch P2).
///
/// Each event carries its own `[start, end]` exposure interval (seconds), so
/// unlike [`calculate_frame_abundance_emg`] — which integrates a *fixed*
/// `[time - rt_cycle_length, time]` and is therefore wrong for unequal event
/// durations / dead time — this integrates the EMG over each event's true
/// interval. Only events overlapping the analyte's RT support (computed once via
/// [`calculate_bounds_emg`] at `target_p`) are integrated, implementing the
/// RT-support truncation policy. Returns `(event_index, abundance)` for events
/// with positive abundance, in ascending event order.
///
/// # Arguments
///
/// * `event_intervals` - per-event `[start, end]` exposure intervals (seconds),
///   indexed by event position in the run timeline.
/// * `mu`, `sigma`, `lambda` - EMG parameters of the analyte's RT profile.
/// * `target_p` - probability mass defining the RT support (e.g. 0.9999).
/// * `step_size` - search step for the bounds binary search.
/// * `n_steps` - integration steps for each interval CDF (defaults to 1000).
pub fn project_emg_over_events(
    event_intervals: &[(f64, f64)],
    mu: f64,
    sigma: f64,
    lambda: f64,
    target_p: f64,
    step_size: f64,
    n_steps: Option<usize>,
) -> Vec<(usize, f64)> {
    // RT support: the [lower, upper] window capturing `target_p` of the mass.
    // Events outside it are skipped (their CDF mass is negligible).
    let (lower, upper) = calculate_bounds_emg(mu, sigma, lambda, step_size, target_p, 20.0, 60.0, n_steps);
    event_intervals
        .iter()
        .enumerate()
        .filter(|(_, &(start, end))| end >= lower && start <= upper)
        .filter_map(|(idx, &(start, end))| {
            let abundance = emg_cdf_range(start, end, mu, sigma, lambda, n_steps);
            if abundance > 0.0 {
                Some((idx, abundance))
            } else {
                None
            }
        })
        .collect()
}

/// Parallel [`project_emg_over_events`] over many analytes sharing one event
/// timeline. `rts`/`sigmas`/`lambdas` are aligned per analyte; returns one
/// `(event_index, abundance)` list per analyte.
pub fn project_emg_over_events_par(
    event_intervals: &[(f64, f64)],
    rts: Vec<f64>,
    sigmas: Vec<f64>,
    lambdas: Vec<f64>,
    target_p: f64,
    step_size: f64,
    num_threads: usize,
    n_steps: Option<usize>,
) -> Vec<Vec<(usize, f64)>> {
    let thread_pool = ThreadPoolBuilder::new().num_threads(num_threads).build().unwrap();
    thread_pool.install(|| {
        rts.into_par_iter()
            .zip(sigmas.into_par_iter())
            .zip(lambdas.into_par_iter())
            .map(|((rt, sigma), lambda)| {
                project_emg_over_events(event_intervals, rt, sigma, lambda, target_p, step_size, n_steps)
            })
            .collect()
    })
}

/// Returns the CDF in the range [sample_start, sample_end] for a Normal(mean, std_dev).
pub fn normal_cdf_range(lower_limit: f64, upper_limit: f64, mean: f64, std_dev: f64) -> f64 {
    let cdf_start = custom_cdf_normal(lower_limit, mean, std_dev);
    let cdf_end = custom_cdf_normal(upper_limit, mean, std_dev);
    cdf_end - cdf_start
}

/// Calculate the bounding interval [lower, upper] around `mean` that captures `target` total probability
/// using a binary search across a discretized search space. This mirrors `calculate_bounds_emg`.
pub fn calculate_bounds_gaussian(
    mean: f64,
    sigma: f64,
    step_size: f64,
    target: f64,
    lower_start: f64,
    upper_start: f64
) -> (f64, f64) {
    assert!((0.0..=1.0).contains(&target), "target must be in [0, 1]");

    let lower_initial = mean - lower_start * sigma;
    let upper_initial = mean + upper_start * sigma;

    let steps = ((upper_initial - lower_initial) / step_size).ceil() as usize;
    let search_space: Vec<f64> = (0..=steps)
        .map(|i| lower_initial + i as f64 * step_size)
        .collect();

    let calc_cdf = |low: usize, high: usize| -> f64 {
        normal_cdf_range(search_space[low], search_space[high], mean, sigma)
    };

    // 1) Find upper cutoff
    let (mut low, mut high) = (0, steps);
    while low < high {
        let mid = low + (high - low) / 2;
        if calc_cdf(0, mid) < target {
            low = mid + 1;
        } else {
            high = mid;
        }
    }
    let upper_cutoff_index = low;

    // 2) Find lower cutoff
    low = 0;
    high = upper_cutoff_index;
    while low < high {
        let mid = high - (high - low) / 2;
        if calc_cdf(mid, upper_cutoff_index) < target {
            high = mid - 1;
        } else {
            low = mid;
        }
    }
    let lower_cutoff_index = high;

    (search_space[lower_cutoff_index], search_space[upper_cutoff_index])
}

/// Returns all scan indices (0-based) that fall into the range where Normal(mean, sigma)
/// has at least `target_p` coverage.
///
/// For timsTOF data, `inverse_ion_mobility` runs backward (highest to lowest values correspond to scans).
///
/// # Arguments
///
/// - `inverse_ion_mobility`: The inverse ion mobility values for all scans (descending order).
/// - `mean`: The mean of the Gaussian distribution.
/// - `sigma`: The standard deviation of the Gaussian distribution.
/// - `target_p`: The target probability to capture.
/// - `step_size`: Step size for searching bounds.
/// - `n_lower_start`: Initial lower bound factor (relative to sigma).
/// - `n_upper_start`: Initial upper bound factor (relative to sigma).
///
/// # Returns
///
/// A `Vec<usize>` containing all scan indices (0-based) within the computed range.
///
/// # Example
///
/// ```rust
/// use mscore::algorithm::utility::calculate_scan_occurrence_gaussian;
///
/// let inverse_ion_mobility = vec![0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3];
/// let scans = calculate_scan_occurrence_gaussian(
///     &inverse_ion_mobility,
///     1.1,  // mean
///     0.001,  // sigma
///     0.9999, // target probability
///     0.01, // step size
///     3.0,  // n_lower_start
///     3.0   // n_upper_start
/// );
///
/// assert_eq!(scans, vec![2]); // Scans corresponding to 1.3 ± 1σ
/// ```
pub fn calculate_scan_occurrence_gaussian(
    inverse_ion_mobility: &[f64],
    mean: f64,
    sigma: f64,
    target_p: f64,
    step_size: f64,
    n_lower_start: f64,
    n_upper_start: f64,
) -> Vec<i32> {
    // Calculate bounds for the Gaussian
    let (ims_lower, ims_upper) = calculate_bounds_gaussian(mean, sigma, step_size, target_p, n_lower_start, n_upper_start);

    // Create a list of tuples (inverse_ion_mobility_value, index) in reverse order
    let indexed_values: Vec<(f64, usize)> = inverse_ion_mobility
        .iter()
        .rev()
        .enumerate()
        .map(|(i, &val)| (val, i))
        .collect();

    // `total_cmp`, not `partial_cmp().unwrap()`: the same NaN hazard as the retention-time search
    // above — a non-finite bound would otherwise panic a worker thread rather than degrade.
    // Find the closest index to ims_lower
    let upper_idx = indexed_values
        .iter()
        .enumerate()
        .min_by(|(_, (val_a, _)), (_, (val_b, _))| {
            (val_a - ims_lower).abs().total_cmp(&(val_b - ims_lower).abs())
        })
        .map(|(idx, _)| idx)
        .unwrap_or(0);

    // Find the closest index to ims_upper
    let lower_idx = indexed_values
        .iter()
        .enumerate()
        .min_by(|(_, (val_a, _)), (_, (val_b, _))| {
            (val_a - ims_upper).abs().total_cmp(&(val_b - ims_upper).abs())
        })
        .map(|(idx, _)| idx)
        .unwrap_or(indexed_values.len() - 1);

    // Extract the indices of the scans in the found range
    if lower_idx <= upper_idx {
        indexed_values[lower_idx..=upper_idx]
            .iter()
            .map(|&(_, idx)| idx as i32)
            .collect()
    } else {
        Vec::new()
    }
}


/// Compute the abundance in each occurrence frame by looking at
/// the probability of Normal(mean, sigma) within `[time - rt_cycle_length, time]`.
pub fn calculate_abundance_gaussian(
    time_map: &HashMap<i32, f64>,
    occurrences: &[i32],
    mean: f64,
    sigma: f64,
    cycle_length: f64,
) -> Vec<f64> {
    let mut frame_abundance = Vec::new();

    for &occurrence in occurrences {
        if let Some(&time) = time_map.get(&occurrence) {
            let start = time - cycle_length;
            let val = normal_cdf_range(start, time, mean, sigma);
            frame_abundance.push(val);
        }
    }

    frame_abundance
}

pub fn calculate_scan_occurrences_gaussian_par(
    times: &[f64],
    means: Vec<f64>,
    sigmas: Vec<f64>,
    target_p: f64,
    step_size: f64,
    n_lower_start: f64,
    n_upper_start: f64,
    num_threads: usize
) -> Vec<Vec<i32>> {
    let thread_pool = ThreadPoolBuilder::new()
        .num_threads(num_threads)
        .build()
        .unwrap();

    thread_pool.install(|| {
        means.into_par_iter()
            .zip(sigmas.into_par_iter())
            .map(|(m, s)| {
                calculate_scan_occurrence_gaussian(
                    times,
                    m,
                    s,
                    target_p,
                    step_size,
                    n_lower_start,
                    n_upper_start
                )
            })
            .collect()
    })
}

/// Parallel version for multiple (mean, sigma) pairs to get abundance
pub fn calculate_scan_abundances_gaussian_par(
    time_map: &HashMap<i32, f64>,
    occurrences: Vec<Vec<i32>>,
    means: Vec<f64>,
    sigmas: Vec<f64>,
    cycle_length: f64,
    num_threads: usize
) -> Vec<Vec<f64>> {
    let thread_pool = ThreadPoolBuilder::new()
        .num_threads(num_threads)
        .build()
        .unwrap();

    thread_pool.install(|| {
        occurrences.into_par_iter()
            .zip(means.into_par_iter())
            .zip(sigmas.into_par_iter())
            .map(|((occ, m), s)| {
                calculate_abundance_gaussian(
                    time_map,
                    &occ,
                    m,
                    s,
                    cycle_length
                )
            })
            .collect()
    })
}

#[cfg(test)]
mod tests {
    use super::*;

    fn approx_eq(a: f64, b: f64, epsilon: f64) -> bool {
        (a - b).abs() < epsilon
    }

    #[test]
    fn test_project_emg_reduces_to_legacy_under_uniform_intervals() {
        // A uniform frame timeline: frames at t = 1..=N seconds, each frame's
        // exposure interval being [t - cycle, t] (the legacy assumption).
        let (mu, sigma, lambda) = (25.0, 1.5, 0.3);
        let cycle = 1.0_f64;
        let times: Vec<f64> = (1..=60).map(|i| i as f64).collect();
        let intervals: Vec<(f64, f64)> = times.iter().map(|&t| (t - cycle, t)).collect();

        // Legacy path: occurrence frames + per-frame abundance over [t-cycle, t].
        let occ = calculate_frame_occurrence_emg(&times, mu, sigma, lambda, 0.9999, 0.01, None);
        let mut time_map = std::collections::HashMap::new();
        for (i, &t) in times.iter().enumerate() {
            time_map.insert((i + 1) as i32, t); // 1-indexed frame ids
        }
        let legacy_abund =
            calculate_frame_abundance_emg(&time_map, &occ, mu, sigma, lambda, cycle, None);

        // New event-interval path.
        let projected = project_emg_over_events(&intervals, mu, sigma, lambda, 0.9999, 0.01, None);

        // The non-negligible projected events must reproduce the legacy
        // per-frame abundances (frame id == event_index + 1).
        for (&frame_id, &abund) in occ.iter().zip(legacy_abund.iter()) {
            if abund <= 0.0 {
                continue;
            }
            let event_idx = (frame_id - 1) as usize;
            let found = projected.iter().find(|(idx, _)| *idx == event_idx);
            assert!(found.is_some(), "event {event_idx} (frame {frame_id}) missing from projection");
            let (_, p_abund) = found.unwrap();
            assert!(
                (p_abund - abund).abs() < 1e-9,
                "abundance mismatch at frame {frame_id}: legacy {abund} vs projected {p_abund}"
            );
        }
    }

    #[test]
    fn test_project_emg_respects_unequal_event_durations() {
        // Two events covering the SAME [start, end] as one wide event must split
        // its mass (area-preserving) — the legacy fixed-cycle kernel can't do this.
        let (mu, sigma, lambda) = (10.0, 1.0, 0.2);
        let wide = vec![(8.0, 12.0)];
        let split = vec![(8.0, 10.0), (10.0, 12.0)];
        let wide_p = project_emg_over_events(&wide, mu, sigma, lambda, 0.9999, 0.01, None);
        let split_p = project_emg_over_events(&split, mu, sigma, lambda, 0.9999, 0.01, None);
        let wide_sum: f64 = wide_p.iter().map(|(_, a)| a).sum();
        let split_sum: f64 = split_p.iter().map(|(_, a)| a).sum();
        // Mass is conserved up to trapezoidal-quadrature discretization (each
        // interval is integrated with n_steps; the split is in fact finer).
        assert!((wide_sum - split_sum).abs() < 1e-3, "split must conserve mass: {wide_sum} vs {split_sum}");
        assert_eq!(split_p.len(), 2, "both sub-intervals should carry mass");
    }

    #[test]
    fn test_normal_cdf_range() {
        let mean = 0.0;
        let std_dev = 1.0;

        // For a standard normal, nearly all probability is within ~[-10, 10].
        // So normal_cdf_range(-10, 10, 0, 1) should be ~1.0
        let prob_all = normal_cdf_range(-10.0, 10.0, mean, std_dev);
        assert!(approx_eq(prob_all, 1.0, 1e-6),
                "CDF range from -10 to 10 should capture nearly all probability, got {prob_all}");

        // Check an interval around mean ± 1σ -> about 68% of mass
        let prob_1sigma = normal_cdf_range(-1.0, 1.0, mean, std_dev);
        assert!(
            (prob_1sigma - 0.68).abs() < 0.02,
            "Expected ~0.68 within ±1σ, got {prob_1sigma}"
        );
    }

    #[test]
    fn test_calculate_bounds_gaussian() {
        let mean = 0.0;
        let sigma = 1.0;
        let target = 0.68;
        // We'll discretize in steps of 0.1
        let (low, high) = calculate_bounds_gaussian(mean, sigma, 0.01, target, 5.0, 5.0);

        // Check that the coverage is close to 0.68
        let coverage = normal_cdf_range(low, high, mean, sigma);
        assert!(
            (coverage - target).abs() < 0.1,
            "Expected coverage ~0.68, got {coverage} for interval [{low}, {high}]"
        );
    }

    #[test]
    fn test_calculate_frame_occurrence_gaussian() {
        // Suppose we have 10 frames of retention times from 0.0 to 9.0
        let retention_times: Vec<f64> = (0..10).map(|x| x as f64).collect();
        let mean = 5.0;   // Centered around the 5th second
        let sigma = 1.0;
        let target_p = 0.68;
        let step_size = 0.1;

        // This should capture frames near t=5.0, about ±1.0 in the "most probable" sense.
        // "lower_start" and "upper_start" here are up to you; let's do 5.0 each side
        let frames = calculate_scan_occurrence_gaussian(
            &retention_times,
            mean,
            sigma,
            target_p,
            step_size,
            5.0,
            5.0
        );

        // Expect frames near 4, 5, 6
        // Because those times (4.0, 5.0, 6.0) are the main chunk of ±1σ around 5.0
        assert!(
            !frames.is_empty(),
            "We expect at least a few frames around 5.0"
        );
        assert!(
            frames.contains(&5),
            "We definitely expect the central frame (index=5 in 1-based indexing) to be included"
        );
    }

    #[test]
    fn test_calculate_frame_abundance_gaussian() {
        // Set up a mock time map: frame_index -> time
        // We'll pretend each frame index i runs from i-1 to i in real-time
        let mut time_map = HashMap::new();
        for i in 1..=5 {
            time_map.insert(i as i32, i as f64);
        }

        // Suppose we only have two frames to check
        let occurrences = vec![1, 3];
        let mean = 3.0;
        let sigma = 1.0;
        let im_cycle_length = 1.0;

        let abundances = calculate_abundance_gaussian(
            &time_map,
            &occurrences,
            mean,
            sigma,
            im_cycle_length
        );

        // We'll do a basic sanity check:
        // - For frame 1, it integrates from time=0 to time=1.
        // - For frame 3, from 2 to 3.
        assert_eq!(abundances.len(), 2, "We should have 2 abundance values");
        let (a1, a2) = (abundances[0], abundances[1]);

        // The second abundance (covering [2,3]) should be bigger,
        // because it's closer to mean=3.0
        assert!(
            a2 > a1,
            "Expected frame near t=3 to have higher abundance than t=1"
        );
    }

    #[test]
    fn test_parallel_functions() {
        // Just a quick sanity check
        let retention_times: Vec<f64> = (0..10).map(|x| x as f64).collect();
        let means = vec![3.0, 5.0];
        let sigmas = vec![1.0, 1.5];

        let target_p = 0.68;
        let step_size = 0.1;
        let num_threads = 2;

        let res_occurrences = calculate_scan_occurrences_gaussian_par(
            &retention_times,
            means.clone(),
            sigmas.clone(),
            target_p,
            step_size,
            5.0,
            5.0,
            num_threads
        );
        assert_eq!(res_occurrences.len(), 2, "Should produce 2 sets of occurrences");

        // Mock time_map for abundances
        let mut time_map = HashMap::new();
        for i in 1..=10 {
            time_map.insert(i, i as f64);
        }

        let res_abundances = calculate_scan_abundances_gaussian_par(
            &time_map,
            res_occurrences,
            means,
            sigmas,
            1.0,          // rt_cycle_length
            num_threads
        );
        assert_eq!(res_abundances.len(), 2, "Should produce 2 sets of abundances");
    }
}

#[cfg(test)]
mod emg_degenerate_tests {
    use super::*;

    /// A degenerate EMG parameter upstream used to make the retention-time bounds non-finite, and
    /// the closest-frame search then unwrapped a `None` from `partial_cmp` and panicked a worker
    /// thread mid-simulation. It must degrade to "this peptide occupies no frames" instead.
    #[test]
    fn non_finite_bounds_return_no_frames_instead_of_panicking() {
        let rts: Vec<f64> = (0..100).map(|i| i as f64 * 0.5).collect();
        for (sigma, lambda_) in [(f64::NAN, 1.0), (1.0, f64::NAN), (0.0, 0.0), (f64::INFINITY, 1.0)] {
            let frames = calculate_frame_occurrence_emg(&rts, 10.0, sigma, lambda_, 0.999, 0.0001, None);
            assert!(
                frames.is_empty() || frames.iter().all(|&f| f >= 0),
                "sigma={sigma}, lambda={lambda_} produced {frames:?}"
            );
        }
    }

    #[test]
    fn well_behaved_parameters_still_return_a_contiguous_frame_window() {
        let rts: Vec<f64> = (0..100).map(|i| i as f64 * 0.5).collect();
        let frames = calculate_frame_occurrence_emg(&rts, 10.0, 0.5, 0.5, 0.999, 0.0001, None);
        assert!(!frames.is_empty(), "a valid peak must occupy at least one frame");
        for w in frames.windows(2) {
            assert_eq!(w[1], w[0] + 1, "frame window must be contiguous");
        }
    }
}