data_beans/

qc.rs

use crate::sparse_data_visitors::*;
use crate::sparse_io::*;
use crate::sparse_io_vector::*;

use indicatif::ParallelProgressIterator;
use legume_numeric::matrix::sparse_stat::{SparseColumnRunningStatistics, SparseRunningStatistics};
use legume_numeric::matrix::traits::RunningStatOps;
use legume_numeric::matrix::utils::partition_by_membership;
use log::warn;
use rayon::prelude::*;
use std::sync::{Arc, Mutex};

// The binary compiles its own copy of this module and reads only part of
// the axis.
#[allow(dead_code)]
pub mod axis;
pub use axis::Axis;

use rustc_hash::FxHashMap as HashMap;

#[derive(Clone)]
pub struct SqueezeCutoffs {
    pub row: usize,
    pub column: usize,
}

/// Whether a row or column with `nnz` non-zeros falls below `cutoff` and is
/// squeezed out. Every preview of a cutoff counts with this rule, so what it
/// reports matches what [`squeeze_by_nnz`] removes.
pub fn below_nnz_cutoff(nnz: f32, cutoff: usize) -> bool {
    (nnz as usize) < cutoff
}

/// squeeze out rows and columns with excessive zero values
///
/// The two cutoffs apply independently, each to the non-zeros counted over
/// the whole matrix before anything is removed: a row is kept by its count
/// across every column, the removed ones included, and a column likewise.
/// So a kept row can end below `cutoffs.row` among the kept columns (all
/// zero, even), and squeezing the result again can remove more.
pub fn squeeze_by_nnz(
    data: &dyn SparseIo<IndexIter = Vec<usize>>,
    cutoffs: SqueezeCutoffs,
    block_size: Option<usize>,
    preload: bool,
) -> anyhow::Result<()> {
    // Without a shape there are no counts to compare with the cutoffs.
    anyhow::ensure!(
        data.num_rows().is_some() && data.num_columns().is_some(),
        "missing shape information"
    );
    let (col_stat, row_stat) = collect_column_row_stat(data, block_size)?;

    fn nnz_index(nnz: &[f32], cutoff: usize) -> Option<Vec<usize>> {
        let ret: Vec<usize> = nnz
            .iter()
            .enumerate()
            .filter(|&(_, &x)| !below_nnz_cutoff(x, cutoff))
            .map(|(i, _)| i)
            .collect();

        (!ret.is_empty()).then_some(ret)
    }

    let row_nnz_vec = row_stat.count_positives();
    let col_nnz_vec = col_stat.count_positives();
    let row_idx = nnz_index(&row_nnz_vec, cutoffs.row);
    let col_idx = nnz_index(&col_nnz_vec, cutoffs.column);

    if row_idx.is_none() {
        warn!(
            "No rows can be kept with this cutoff {}!\n\
	     \n\
	     We will stop squeezing on the rows.\n\
	     \n",
            cutoffs.row
        );
    }

    if col_idx.is_none() {
        warn!(
            "No columns can be kept with this cutoff {}!\n\
	     \n\
	     We will stop squeezing on the columns.\n\
	     \n",
            cutoffs.column
        );
    }

    // A side on which nothing falls below its cutoff is left whole, and the
    // backend as it is when that is both.
    let row_idx = row_idx.filter(|kept| kept.len() < row_nnz_vec.len());
    let col_idx = col_idx.filter(|kept| kept.len() < col_nnz_vec.len());
    if row_idx.is_none() && col_idx.is_none() {
        return Ok(());
    }

    let file = data.get_backend_file_name();
    let backend = data.backend_type();

    let mut data = open_sparse_matrix(file, &backend)?;
    if preload {
        data.preload_columns()?;
    }
    data.subset_columns_rows(col_idx.as_ref(), row_idx.as_ref())
}

/// collect row-wise sufficient statistics for Q/C
/// * `data` - `SparseIoVec` across many data matrices
/// * `block_size` - a block size for each parallelized job
pub fn collect_row_stat_across_vec(
    data: &SparseIoVec,
    block_size: Option<usize>,
) -> anyhow::Result<SparseRunningStatistics<f32>> {
    let mut row_stat = SparseRunningStatistics::new(data.num_rows());
    data.visit_columns_by_block(
        &row_stat_vec_visitor,
        &EmptyArgs {},
        &mut row_stat,
        block_size,
    )?;
    Ok(row_stat)
}

/// collect row statistics for each group of columns
/// * `data` - `SparseIo`
/// * `column_membership` - a hashmap assign columns to groups
/// * `block_size` - a block size for each parallelized job
#[allow(clippy::type_complexity)]
pub fn collect_stratified_row_stat_across_vec(
    data: &SparseIoVec,
    column_membership: &HashMap<Box<str>, Box<str>>,
    block_size: Option<usize>,
) -> anyhow::Result<(Vec<Box<str>>, Vec<SparseRunningStatistics<f32>>)> {
    let column_names = data.column_names()?;
    let default = "".to_string().into_boxed_str();
    let membership = column_names
        .into_iter()
        .map(|k| column_membership.get(&k).unwrap_or(&default).clone())
        .collect::<Vec<_>>();

    let partitions = partition_by_membership(&membership, None);
    let mut group_names = Vec::with_capacity(partitions.len());
    let mut group_stats = Vec::with_capacity(partitions.len());
    let num_features = data.num_rows();

    for (k, cols) in partitions {
        let jobs = create_jobs(cols.len(), num_features, block_size);
        let mut row_stat = SparseRunningStatistics::new(data.num_rows());
        let arc_stat = Arc::new(Mutex::new(&mut row_stat));

        jobs.par_iter()
            .progress_with(styled_progress_bar(jobs.len() as u64, "blocks"))
            .for_each(|&(lb, ub)| {
                let cols_sub = cols[lb..ub].iter().cloned();
                let csc = data
                    .read_columns_csc(cols_sub)
                    .expect("failed to read data");
                let mut stat = arc_stat.lock().expect("failed to lock row_stat");
                stat.add_csc(&csc);
            });

        group_names.push(k);
        group_stats.push(row_stat);
    }

    Ok((group_names, group_stats))
}

/// collect row-wise sufficient statistics for Q/C
/// * `data` - `SparseIo`
/// * `block_size` - a block size for each parallelized job
pub fn collect_row_stat(
    data: &dyn SparseIo<IndexIter = Vec<usize>>,
    block_size: Option<usize>,
) -> anyhow::Result<SparseRunningStatistics<f32>> {
    let nrows = data.num_rows().unwrap_or(0);
    fold_column_blocks(
        data,
        block_size,
        || SparseRunningStatistics::new(nrows),
        |stat, csc, _| stat.add_csc(csc),
        |stat, other| stat.merge(&other),
    )
}

/// Column blocks folded into one accumulator at a leaf of
/// [`fold_in_fixed_order`] by a streaming pass over `n_blocks` of them. A
/// leaf is folded on one thread and allocates one accumulator the size of
/// the feature axis, so the pass is cut into at most a few hundred leaves:
/// each block is its own leaf, and its own task, up to that many blocks,
/// and a longer pass shares an accumulator among neighbouring blocks. The
/// leaf depends on the number of blocks alone, as the fold's order must.
pub(crate) fn blocks_per_leaf(n_blocks: usize) -> usize {
    const LEAVES: usize = 256;
    n_blocks.div_ceil(LEAVES).max(1)
}

/// Fold `items` into one accumulator over a fixed binary tree: a run of at
/// most `leaf` items is folded in order into a fresh accumulator; a longer
/// run is halved, the halves folded on separate rayon tasks and merged left
/// into right.
///
/// The grouping of the floating-point sums then depends on `items.len()` and
/// `leaf` alone, not on the thread count or on how rayon schedules the work,
/// so the result is the same bits run to run. (A plain rayon `fold`/`reduce`
/// splits where work stealing happens to split.)
pub(crate) fn fold_in_fixed_order<T, A, I, F, M>(
    items: &[T],
    leaf: usize,
    init: &I,
    fold: &F,
    merge: &M,
) -> anyhow::Result<A>
where
    T: Sync,
    A: Send,
    I: Fn() -> A + Sync,
    F: Fn(&mut A, &T) -> anyhow::Result<()> + Sync,
    M: Fn(&mut A, A) + Sync,
{
    if items.len() <= leaf.max(1) {
        let mut acc = init();
        for item in items {
            fold(&mut acc, item)?;
        }
        return Ok(acc);
    }
    let (left, right) = items.split_at(items.len() / 2);
    let (left, right) = rayon::join(
        || fold_in_fixed_order(left, leaf, init, fold, merge),
        || fold_in_fixed_order(right, leaf, init, fold, merge),
    );
    let mut acc = left?;
    merge(&mut acc, right?);
    Ok(acc)
}

/// Read the columns of `data` block by block in parallel and fold the blocks
/// in a fixed order ([`fold_in_fixed_order`]): `add(acc, csc, lb)` folds a
/// block, given the index of its first column, into the accumulator of its
/// leaf (a fresh `init()`), and `merge` joins two accumulators. Sums that
/// run across blocks, as a row's do, are then the same bits on any number of
/// threads. A backend without a shape has no blocks.
///
/// # Errors
/// If a block cannot be read.
fn fold_column_blocks<A, I, F, M>(
    data: &dyn SparseIo<IndexIter = Vec<usize>>,
    block_size: Option<usize>,
    init: I,
    add: F,
    merge: M,
) -> anyhow::Result<A>
where
    A: Send,
    I: Fn() -> A + Sync,
    F: Fn(&mut A, &CscMatrix<f32>, usize) + Sync,
    M: Fn(&mut A, A) + Sync,
{
    let ncols = data.num_columns().unwrap_or(0);
    let nrows = data.num_rows().unwrap_or(0);
    let jobs = create_jobs(ncols, nrows, block_size);
    let prog_bar = styled_progress_bar(jobs.len() as u64, "blocks");
    fold_in_fixed_order(
        &jobs,
        blocks_per_leaf(jobs.len()),
        &init,
        &|acc, &(lb, ub)| {
            let csc = data.read_columns_csc((lb..ub).collect())?;
            add(acc, &csc, lb);
            prog_bar.inc(1);
            Ok(())
        },
        &merge,
    )
}

/// The rows of `select_rows` (those within `nrows`) as a mask over all
/// rows, and empty statistics of `ncols` columns over them: the rows
/// selected are the denominator of the mean.
pub fn masked_column_stat(
    ncols: usize,
    nrows: usize,
    select_rows: &[usize],
) -> (Vec<bool>, SparseColumnRunningStatistics<f32>) {
    let mut mask = vec![false; nrows];
    for &r in select_rows {
        if r < nrows {
            mask[r] = true;
        }
    }
    let selected = mask.iter().filter(|x| **x).count();
    (mask, SparseColumnRunningStatistics::new(ncols, selected))
}

/// collect column-wise sufficient statistics for Q/C
/// * `data` - `SparseIoVec` across many data matrices
/// * `select_rows` - selected row indices
/// * `block_size` - a block size for each parallelized job
pub fn collect_column_stat_across_vec(
    data: &SparseIoVec,
    select_rows: Option<&[usize]>,
    block_size: Option<usize>,
) -> anyhow::Result<SparseColumnRunningStatistics<f32>> {
    let ncols = data.num_columns();
    let nrows_total = data.num_rows();

    let (row_mask, mut col_stat) = match select_rows {
        Some(sel) => {
            let (mask, stat) = masked_column_stat(ncols, nrows_total, sel);
            (Some(mask), stat)
        }
        None => (None, SparseColumnRunningStatistics::new(ncols, nrows_total)),
    };
    data.visit_columns_by_block(&col_stat_visitor, &row_mask, &mut col_stat, block_size)?;
    Ok(col_stat)
}

/// collect column-wise and row-wise sufficient statistics for Q/C, reading
/// every block once ([`collect_column_stat`] and [`collect_row_stat`])
/// * `data` - `SparseIo`
/// * `block_size` - a block size for each parallelized job
pub fn collect_column_row_stat(
    data: &dyn SparseIo<IndexIter = Vec<usize>>,
    block_size: Option<usize>,
) -> anyhow::Result<(
    SparseColumnRunningStatistics<f32>,
    SparseRunningStatistics<f32>,
)> {
    let ncols = data.num_columns().unwrap_or(0);
    let nrows = data.num_rows().unwrap_or(0);
    // A column's statistics are written once, by the block that holds it,
    // so the order the blocks reach the lock in does not show in them.
    let col_stat = Mutex::new(SparseColumnRunningStatistics::<f32>::new(ncols, nrows));
    let row_stat = fold_column_blocks(
        data,
        block_size,
        || SparseRunningStatistics::<f32>::new(nrows),
        |row_stat, csc, lb| {
            col_stat
                .lock()
                .expect("failed to lock the column statistics")
                .add_csc(csc, lb);
            row_stat.add_csc(csc);
        },
        |row_stat, other| row_stat.merge(&other),
    )?;
    let col_stat = col_stat
        .into_inner()
        .expect("failed to unlock the column statistics");
    Ok((col_stat, row_stat))
}

/// collect column-wise sufficient statistics for Q/C
/// * `data` - `SparseIo`
/// * `block_size` - a block size for each parallelized job
pub fn collect_column_stat(
    data: &dyn SparseIo<IndexIter = Vec<usize>>,
    block_size: Option<usize>,
) -> anyhow::Result<SparseColumnRunningStatistics<f32>> {
    let ncols = data.num_columns().unwrap_or(0);
    let nrows = data.num_rows().unwrap_or(0);
    // A column's statistics are written once, by the block that holds it,
    // so the order the blocks reach the lock in does not show in them.
    let col_stat = Mutex::new(SparseColumnRunningStatistics::<f32>::new(ncols, nrows));
    fold_column_blocks(
        data,
        block_size,
        || (),
        |_, csc, lb| {
            col_stat
                .lock()
                .expect("failed to lock the column statistics")
                .add_csc(csc, lb)
        },
        |_, _| (),
    )?;
    Ok(col_stat
        .into_inner()
        .expect("failed to unlock the column statistics"))
}

struct EmptyArgs {}

fn row_stat_vec_visitor(
    job: (usize, usize),
    data: &SparseIoVec,
    _: &EmptyArgs,
    arc_stat: Arc<Mutex<&mut SparseRunningStatistics<f32>>>,
) -> anyhow::Result<()> {
    let (lb, ub) = job;
    let csc = data.read_columns_csc(lb..ub)?;

    let mut stat = arc_stat.lock().expect("failed to lock row_stat");
    stat.add_csc(&csc);
    Ok(())
}

fn col_stat_visitor(
    job: (usize, usize),
    data: &SparseIoVec,
    row_mask: &Option<Vec<bool>>,
    arc_stat: Arc<Mutex<&mut SparseColumnRunningStatistics<f32>>>,
) -> anyhow::Result<()> {
    let (lb, ub) = job;
    let csc = data.read_columns_csc(lb..ub)?;
    let mut stat = arc_stat.lock().expect("failed to lock col_stat");
    match row_mask {
        Some(mask) => stat.add_csc_masked(&csc, lb, mask),
        None => stat.add_csc(&csc, lb),
    }
    Ok(())
}

//////////////////////////////////////////////////////////////////////////////////
// Automatic nnz-cutoff selection + ASCII histogram (shared by `squeeze` and by //
// callers that want cell-calling on a per-column nnz vector).                   //
//////////////////////////////////////////////////////////////////////////////////

/// Bin width of the trough search, in natural-log units (~5% per bin).
const TROUGH_BIN: f64 = 0.05;
/// Gaussian smoothing of the binned density, in bins.
const TROUGH_SMOOTH_BINS: f64 = 3.0;
/// A cut needs the trough at most this fraction of the lower of its two peaks.
const TROUGH_MAX_DEPTH: f64 = 0.25;
/// Ambient and cell peaks sit at least a decade apart; closer modes (a doublet
/// bump, a low-complexity cell type) are not an empty↔cell boundary.
const TROUGH_MIN_PEAK_RATIO: f64 = 10.0;
/// Each side of a cut must hold at least this many columns (or 0.1% of all).
const TROUGH_MIN_SIDE: usize = 20;

/// Suggest an nnz cutoff at the **deepest trough** of the nnz distribution,
/// the gap between the ambient (empty barcode) peak and the cell peak.
/// Columns with `nnz >= cutoff` are kept.
///
/// The density lives on the log axis, where ambient and cells form two
/// peaks decades apart; in linear units the cell peak is spread so thin
/// that the trough in front of it vanishes. The bins, though, come from
/// the actual counts: each integer count `v` spreads its mass over its own
/// interval `[v - 1/2, v + 1/2)` mapped to `ln(1 + ·)`, so small counts,
/// whose log spacing exceeds a bin, leave no empty bins to pass for troughs.
///
/// A cut is proposed only when the smoothed density at the trough is at most
/// [`TROUGH_MAX_DEPTH`] of the lower of the two peaks it separates, the
/// peaks are at least [`TROUGH_MIN_PEAK_RATIO`] apart, and each side holds
/// enough columns. Unimodal data — already-called cells included — gets
/// `None`, whatever its tails look like. Deterministic, no RNG.
pub fn suggest_nnz_cutoff(nnz: &[f32]) -> Option<usize> {
    let n = nnz.len();
    let min_side = TROUGH_MIN_SIDE.max(n / 1000);
    if n < 2 * min_side {
        return None;
    }

    let mut vals: Vec<u64> = nnz.iter().map(|&x| x.max(0.0).round() as u64).collect();
    vals.sort_unstable();
    let (vmin, vmax) = (vals[0], vals[n - 1]);
    if vmin == vmax {
        return None;
    }

    // Count `v`'s interval `[v - 1/2, v + 1/2)` on the `ln(1 + ·)` axis.
    let edge = |v: u64, half: f64| (v as f64 + 1.0 + half).ln();
    let lo = edge(vmin, -0.5);
    let nbins = ((edge(vmax, 0.5) - lo) / TROUGH_BIN).ceil() as usize;
    let mut hist = vec![0.0_f64; nbins];
    let mut i = 0;
    while i < n {
        let v = vals[i];
        let mut j = i;
        while j < n && vals[j] == v {
            j += 1;
        }
        let a = (edge(v, -0.5) - lo) / TROUGH_BIN;
        let b = (edge(v, 0.5) - lo) / TROUGH_BIN;
        let mass = (j - i) as f64 / (b - a);
        for (k, h) in hist
            .iter_mut()
            .enumerate()
            .take((b.ceil() as usize).min(nbins))
            .skip(a.floor() as usize)
        {
            let overlap = b.min(k as f64 + 1.0) - a.max(k as f64);
            if overlap > 0.0 {
                *h += mass * overlap;
            }
        }
        i = j;
    }

    // Gaussian smoothing (finite ±3σ kernel, so a real gap stays exactly 0).
    let radius = (3.0 * TROUGH_SMOOTH_BINS).ceil() as isize;
    let kernel: Vec<f64> = (-radius..=radius)
        .map(|d| (-0.5 * (d as f64 / TROUGH_SMOOTH_BINS).powi(2)).exp())
        .collect();
    let smooth: Vec<f64> = (0..nbins as isize)
        .map(|k| {
            (-radius..=radius)
                .filter_map(|d| {
                    let t = k + d;
                    (0..nbins as isize)
                        .contains(&t)
                        .then(|| hist[t as usize] * kernel[(d + radius) as usize])
                })
                .sum()
        })
        .collect();

    // Mass left of each bin, and the running peaks from either end.
    let mut below = vec![0.0_f64; nbins + 1];
    for k in 0..nbins {
        below[k + 1] = below[k] + hist[k];
    }
    let mut left_peak = vec![(0.0_f64, 0usize); nbins];
    for k in 0..nbins {
        let prev = if k > 0 { left_peak[k - 1] } else { (-1.0, 0) };
        left_peak[k] = if smooth[k] > prev.0 {
            (smooth[k], k)
        } else {
            prev
        };
    }
    let mut right_peak = vec![(0.0_f64, 0usize); nbins];
    for k in (0..nbins).rev() {
        let next = if k + 1 < nbins {
            right_peak[k + 1]
        } else {
            (-1.0, k)
        };
        right_peak[k] = if smooth[k] >= next.0 {
            (smooth[k], k)
        } else {
            next
        };
    }

    let min_bins_apart = TROUGH_MIN_PEAK_RATIO.ln() / TROUGH_BIN;
    let total = below[nbins];
    let mut best: Option<(f64, usize, usize)> = None; // (depth, first, last) of the deepest run
    for t in 1..nbins.saturating_sub(1) {
        let (l_mass, r_mass) = (below[t], total - below[t + 1]);
        if l_mass < min_side as f64 || r_mass < min_side as f64 {
            continue;
        }
        let ((l_h, l_k), (r_h, r_k)) = (left_peak[t - 1], right_peak[t + 1]);
        if ((r_k - l_k) as f64) < min_bins_apart || l_h <= 0.0 || r_h <= 0.0 {
            continue;
        }
        let depth = smooth[t] / l_h.min(r_h);
        match best {
            Some((d, first, last)) if depth == d && last + 1 == t => best = Some((d, first, t)),
            Some((d, _, _)) if depth >= d => {}
            _ => best = Some((depth, t, t)),
        }
    }

    let Some((depth, first, last)) = best else {
        log::info!("nnz cell-calling: no two peaks a decade apart → unimodal, no cutoff");
        return None;
    };
    // Middle of the deepest run (an exact-zero gap is a run, not a point).
    let x = lo + ((first + last) as f64 / 2.0 + 0.5) * TROUGH_BIN;
    let cutoff = (x.exp() - 1.0).ceil().max(1.0) as usize;
    let favors_cut = depth <= TROUGH_MAX_DEPTH;
    log::info!(
        "nnz cell-calling: deepest trough at nnz {} (depth {:.3}) → {}",
        cutoff,
        depth,
        if favors_cut {
            format!("bimodal, cutoff at nnz {cutoff}")
        } else {
            "unimodal, no cutoff".to_string()
        }
    );
    favors_cut.then_some(cutoff)
}

/// `part` as a percentage of `total` (0 when there is nothing).
pub fn pct(part: usize, total: usize) -> f64 {
    if total > 0 {
        100.0 * part as f64 / total as f64
    } else {
        0.0
    }
}

/// Median of an ascending slice (0 when empty).
pub fn median_of_sorted(sorted: &[f32]) -> f32 {
    let n = sorted.len();
    match n {
        0 => 0.0,
        _ if n.is_multiple_of(2) => (sorted[n / 2 - 1] + sorted[n / 2]) / 2.0,
        _ => sorted[n / 2],
    }
}

/// One bin of a printed histogram that holds values, with the real range
/// of the values that fell into it.
struct HistBin {
    val_min: f32,
    val_max: f32,
    count: usize,
    is_cutoff: bool,
}

/// Bins of a printed histogram across its values' range.
const PRINTED_BINS: usize = 40;

/// A histogram of `values` on the log [`Axis`] every histogram shares, its
/// bins that hold values with their real ranges. Works on any non-negative
/// statistic (nnz, sum, mean, sd); the ranges stay exact `f32` so
/// count-like stats still print as integers.
fn create_log_histogram(values: &[f32], cutoff: usize) -> Vec<HistBin> {
    // The axis reaches a cutoff that is set; with none (cutoff == 0) it
    // spans the values alone.
    let axis = Axis::fit(values, true);
    let axis = if cutoff > 0 {
        axis.including(cutoff as f64)
    } else {
        axis
    };
    let mut bins: Vec<(usize, f32, f32)> =
        vec![(0, f32::INFINITY, f32::NEG_INFINITY); PRINTED_BINS];
    for &val in values.iter().filter(|v| v.is_finite()) {
        let b = &mut bins[axis.bin(val as f64, PRINTED_BINS)];
        b.0 += 1;
        b.1 = b.1.min(val);
        b.2 = b.2.max(val);
    }
    // Mark the first bin at or above the cutoff that holds values, so the
    // arrow always renders. With no cutoff (cutoff == 0, e.g. the
    // `histogram` command) nothing is marked.
    let cutoff_bin = (cutoff > 0)
        .then(|| {
            let at = axis.bin(cutoff as f64, PRINTED_BINS);
            (at..PRINTED_BINS).find(|&k| bins[k].0 > 0)
        })
        .flatten();
    bins.into_iter()
        .enumerate()
        .filter(|(_, b)| b.0 > 0)
        .map(|(k, (count, val_min, val_max))| HistBin {
            val_min,
            val_max,
            count,
            is_cutoff: Some(k) == cutoff_bin,
        })
        .collect()
}

/// Format a statistic value compactly: whole numbers (nnz, integer counts)
/// print without a decimal point; fractional values (mean, sd) get
/// `decimals` places.
pub fn fmt_stat(v: f32, decimals: usize) -> String {
    if v.fract() == 0.0 {
        (v as i64).to_string()
    } else {
        format!("{:.*}", decimals, v)
    }
}

/// A value as a printed histogram's range shows it: whole numbers as they
/// are, others to three significant digits (0.000412, 0.0731, 2.45).
fn fmt_sig(v: f32) -> String {
    if v.fract() == 0.0 {
        return (v as i64).to_string();
    }
    let places = (2 - (v.abs() as f64).log10().floor() as i32).max(0) as usize;
    format!("{:.*}", places, v)
}

/// Print summary statistics + an ASCII log-binned histogram of a per-row or
/// per-column statistic vector (`metric` names it, e.g. "nnz", "sum", "mean").
/// A non-zero `cutoff` marks the cutoff bin and reports how much it removes;
/// an optional `suggested` value reports the trough suggestion.
///
/// Used by `data-beans histogram`, whose result it is; `qc --qc-histogram`
/// prints the same on stderr ([`eprint_nnz_summary`]).
pub fn print_nnz_summary(
    label: &str,
    metric: &str,
    values: &[f32],
    cutoff: usize,
    suggested: Option<usize>,
) {
    let _ = write_nnz_summary(
        &mut std::io::stdout().lock(),
        label,
        metric,
        values,
        cutoff,
        suggested,
    );
}

/// [`print_nnz_summary`] on stderr: for a histogram that is a diagnostic
/// beside a command's result (`qc --qc-histogram`, senna's and pinto's QC),
/// not the result itself, so stdout carries only the result.
pub fn eprint_nnz_summary(
    label: &str,
    metric: &str,
    values: &[f32],
    cutoff: usize,
    suggested: Option<usize>,
) {
    let _ = write_nnz_summary(
        &mut std::io::stderr().lock(),
        label,
        metric,
        values,
        cutoff,
        suggested,
    );
}

/// [`print_nnz_summary`] to `out`.
pub fn write_nnz_summary(
    out: &mut impl std::io::Write,
    label: &str,
    metric: &str,
    values: &[f32],
    cutoff: usize,
    suggested: Option<usize>,
) -> std::io::Result<()> {
    const MAX_BAR_WIDTH: usize = 50; // Maximum width for histogram bars

    let total = values.len();
    let below_cutoff = values
        .iter()
        .filter(|&&x| below_nnz_cutoff(x, cutoff))
        .count();
    let pct_removed = pct(below_cutoff, total);

    // Calculate basic statistics
    let min = values.iter().copied().fold(f32::INFINITY, f32::min);
    let max = values.iter().copied().fold(f32::NEG_INFINITY, f32::max);
    // Summed in f64: an f32 sum stops growing once it dwarfs its addends.
    let sum: f64 = values.iter().map(|&v| v as f64).sum();
    let mean = if total > 0 { sum / total as f64 } else { 0.0 };

    let mut sorted = values.to_vec();
    sorted.sort_unstable_by(f32::total_cmp);
    let median = median_of_sorted(&sorted);

    writeln!(out, "{} {} distribution:", label, metric)?;
    writeln!(out, "  Total: {}", total)?;
    writeln!(
        out,
        "  Min: {}, Max: {}, Mean: {:.2}, Median: {:.2}",
        fmt_stat(min, 2),
        fmt_stat(max, 2),
        mean,
        median
    )?;
    if cutoff > 0 {
        writeln!(
            out,
            "  Cutoff: {} (removes {} / {} = {:.2}%)",
            cutoff, below_cutoff, total, pct_removed
        )?;
    }
    if let Some(s) = suggested {
        let below_s = values.iter().filter(|&&x| below_nnz_cutoff(x, s)).count();
        let pct_s = pct(below_s, total);
        writeln!(
            out,
            "  Suggested cutoff (histogram trough of {}): {} (would remove {} / {} = {:.2}%)",
            metric, s, below_s, total, pct_s
        )?;
    }

    // Log-spaced bins, tracking the real value range per bin
    let hist = create_log_histogram(values, cutoff);

    // Scale bar width on log10(count+1) so a few outlier bins don't flatten the rest
    let max_log_count = hist
        .iter()
        .map(|b| ((b.count as f64) + 1.0).log10())
        .fold(0.0_f64, f64::max)
        .max(1e-9);

    writeln!(
        out,
        "  Histogram (x: {m} in log-spaced bins, each its range; bar: log10(count+1)):",
        m = metric
    )?;
    for b in hist {
        let marker = if b.is_cutoff { " <-- CUTOFF" } else { "" };
        let log_count = ((b.count as f64) + 1.0).log10();
        let bar_width = ((log_count / max_log_count) * MAX_BAR_WIDTH as f64).round() as usize;
        let bar_width = if b.count > 0 { bar_width.max(1) } else { 0 };
        let bar = "█".repeat(bar_width);
        let range = if b.val_min == b.val_max {
            fmt_sig(b.val_min)
        } else {
            format!("{}-{}", fmt_sig(b.val_min), fmt_sig(b.val_max))
        };
        writeln!(out, "    {:>17}: {:>6} {}{}", range, b.count, bar, marker)?;
    }
    Ok(())
}

/// Cut, in place, the rows and columns of the backend at `backend_file`
/// with fewer than `row` and `column` non-zeros: the `--squeeze` of the
/// commands that build or derive a backend. The path is one the caller
/// resolved and finalizes itself, so nothing is zipped here.
pub fn squeeze_in_place(
    backend_file: &str,
    row: usize,
    column: usize,
    block_size: Option<usize>,
) -> anyhow::Result<()> {
    log::info!("Squeeze the backend data {backend_file}");
    let (backend, file) = crate::hdf5_io::resolve_backend_file(backend_file, None)?;
    let data = crate::sparse_io::open_sparse_matrix(&file, &backend)?;
    squeeze_by_nnz(
        data.as_ref(),
        SqueezeCutoffs { row, column },
        block_size,
        true,
    )?;
    let data = crate::sparse_io::open_sparse_matrix(&file, &backend)?;
    log::info!(
        "after squeeze -- data: {} rows x {} columns",
        data.num_rows().unwrap_or_default(),
        data.num_columns().unwrap_or_default()
    );
    Ok(())
}

#[cfg(test)]
mod tests;