Function _calculate_dispersion_stats

src/shared/processing/mod.rs:353–394

(
    log_dispersions: &[f64],
    bin_indices: &[usize],
    mean_bins: &[usize],
)

}

pub fn _calculate_dispersion_stats(
    log_dispersions: &[f64],
    bin_indices: &[usize],
    mean_bins: &[usize],
) -> anyhow::Result<(Vec<f64>, Vec<f64>)> {
    let n_bins = mean_bins.len();
    let mut bin_means = vec![0.0; n_bins];
    let mut bin_stds = vec![0.0; n_bins];
    let mut bin_sums = vec![0.0; n_bins];
    let mut bin_sum_squares = vec![0.0; n_bins];

    for (i, &bin_idx) in bin_indices.iter().enumerate() {
        let disp = log_dispersions[i];
        if !disp.is_nan() {
            bin_sums[bin_idx] += disp;
            bin_sum_squares[bin_idx] += disp * disp;
        }
    }

    for bin_idx in 0..n_bins {
        let count = mean_bins[bin_idx] as f64;
        if count > 0.0 {
            bin_means[bin_idx] = bin_sums[bin_idx] / count;

            if count > 1.0 {
                let variance =
                    (bin_sum_squares[bin_idx] - bin_sums[bin_idx].powi(2) / count) / (count - 1.0);

                // Add small epsilon to prevent zero standard deviation
                let min_variance = 1e-12;
                bin_stds[bin_idx] = (variance.max(min_variance)).sqrt();
            } else {
                bin_stds[bin_idx] = f64::NAN;
            }
        } else {
            bin_means[bin_idx] = f64::NAN;
            bin_stds[bin_idx] = f64::NAN;
        }
    }

    Ok((bin_means, bin_stds))
}