harmos_signal/waveform/

reduce.rs

//! Choosing which samples to keep when there are more than anything downstream
//! can use.
//!
//! Every reducer answers the *indices* it selected, strictly increasing, and
//! never the samples themselves. A chart wants two columns, a store wants four,
//! and a caller that got indices back reads whichever of its own columns it
//! keeps — which is also what makes these compose:
//!
//! ```
//! use harmos_signal::{above_floor, lttb};
//!
//! let frequency: Vec<f64> = (0..1_000).map(|bin| f64::from(bin) * 100.0).collect();
//! let magnitude: Vec<f64> =
//!     (0..1_000).map(|bin| if bin == 100 || bin == 500 { 100.0 } else { 0.1 }).collect();
//!
//! // Drop the noise floor, then thin what is left to fifty points.
//! let kept = above_floor(&magnitude, 40.0);
//! let x: Vec<f64> = kept.iter().map(|bin| frequency[*bin]).collect();
//! let y: Vec<f64> = kept.iter().map(|bin| magnitude[*bin]).collect();
//! let picked: Vec<usize> = lttb(&x, &y, 50).iter().map(|index| kept[*index]).collect();
//!
//! assert!(picked.contains(&100) && picked.contains(&500));
//! ```

/// Thins `(x, y)` to `count` samples by largest-triangle-three-buckets.
///
/// LTTB keeps the shape a naive stride throws away: each bucket contributes the
/// sample forming the largest triangle with the previously kept sample and the
/// mean of the next bucket, so peaks and inflections survive. The first and
/// last samples are always kept.
///
/// Asking for at least as many samples as there are, or for fewer than three,
/// keeps everything.
#[must_use]
pub fn lttb(x: &[f64], y: &[f64], count: usize) -> Vec<usize> {
    let length = x.len().min(y.len());
    if length <= count || count < 3 {
        return (0..length).collect();
    }

    let mut kept = Vec::with_capacity(count);
    kept.push(0);

    let bucket = (length - 2) as f64 / (count - 2) as f64;
    let mut previous = 0;

    for index in 0..count - 2 {
        let start = 1 + (index as f64 * bucket).floor() as usize;
        let end = (1 + ((index + 1) as f64 * bucket).floor() as usize).min(length - 1);
        let ahead = if index < count - 3 {
            (1 + ((index + 2) as f64 * bucket).floor() as usize).min(length)
        } else {
            length
        };

        let (mean_x, mean_y) = mean_of(x, y, end..ahead).unwrap_or((x[length - 1], y[length - 1]));

        let picked = (start..end)
            .max_by(|left, right| {
                area(x[previous], y[previous], x[*left], y[*left], mean_x, mean_y).total_cmp(&area(
                    x[previous],
                    y[previous],
                    x[*right],
                    y[*right],
                    mean_x,
                    mean_y,
                ))
            })
            .unwrap_or(start);

        kept.push(picked);
        previous = picked;
    }

    kept.push(length - 1);
    kept
}

/// Thins `(x, y)` to at most `count` samples with logarithmically spaced buckets.
///
/// The same selection as [`lttb`], with the buckets laid out along `ln(x)` and
/// the triangle measured there too. On a spectrum drawn against a logarithmic
/// frequency axis this is the difference between the first decade getting one
/// sample and getting its share: linear buckets spend nearly every point on the
/// top decade, where the axis has the least room for them.
///
/// A bucket holding no sample contributes none, so the answer may be shorter
/// than `count`. An axis that does not increase, or holds no positive value to
/// take a logarithm of, is reduced by [`lttb`] instead of refused.
#[must_use]
pub fn lttb_log(x: &[f64], y: &[f64], count: usize) -> Vec<usize> {
    let length = x.len().min(y.len());
    if length <= count || count < 3 {
        return (0..length).collect();
    }

    let (first, last) = (positive(x[0]), x[length - 1]);
    if !last.is_finite() || last <= first {
        return lttb(x, y, count);
    }

    let (low, span) = (first.ln(), last.ln() - first.ln());
    if span <= 0.0 {
        return lttb(x, y, count);
    }

    let buckets = count - 2;
    let edge = |bucket: usize| (low + (bucket as f64 / buckets as f64) * span).exp();

    let mut kept = Vec::with_capacity(count);
    kept.push(0);
    let mut previous = 0;

    for bucket in 0..buckets {
        let start = x[..length]
            .partition_point(|at| *at < edge(bucket))
            .min(length - 1);
        let end = x[..length]
            .partition_point(|at| *at < edge(bucket + 1))
            .clamp(start + 1, length);
        let ahead = if bucket < buckets - 1 {
            x[..length]
                .partition_point(|at| *at < edge(bucket + 2).min(last))
                .clamp(end, length)
        } else {
            length
        };

        let (mean_x, mean_y) = mean_of(x, y, end..ahead).unwrap_or((x[length - 1], y[length - 1]));
        // The triangle is measured where the eye sees it: along ln(x).
        let (previous_x, mean_x) = (positive(x[previous]).ln(), positive(mean_x).ln());

        let picked = (start..end)
            .max_by(|left, right| {
                let of = |bin: usize| {
                    area(
                        previous_x,
                        y[previous],
                        positive(x[bin]).ln(),
                        y[bin],
                        mean_x,
                        mean_y,
                    )
                };
                of(*left).total_cmp(&of(*right))
            })
            .unwrap_or(start);

        kept.push(picked);
        previous = picked;
    }

    kept.push(length - 1);
    kept.dedup();
    kept
}

/// Keeps every `n`th index of a run of `length`, for at most `count` of them.
///
/// The cheap reduction, and the one that loses peaks: a stride knows nothing
/// about the samples it skips. It needs neither column, only how many there
/// are.
#[must_use]
pub fn decimate(length: usize, count: usize) -> Vec<usize> {
    if length == 0 || count == 0 {
        return Vec::new();
    }

    let stride = (length as f64 / count as f64).ceil() as usize;
    (0..length).step_by(stride.max(1)).collect()
}

/// Keeps the lowest and highest sample of each of `buckets` equal spans.
///
/// The envelope reduction: up to two indices per bucket, in sample order, so a
/// waveform drawn from the answer still touches every extreme the full run
/// touched. A bucket whose lowest and highest are the same sample contributes
/// it once.
#[must_use]
pub fn extrema(values: &[f64], buckets: usize) -> Vec<usize> {
    if values.is_empty() || buckets == 0 {
        return Vec::new();
    }

    let span = ((values.len() as f64 / buckets as f64).ceil() as usize).max(1);
    let mut kept = Vec::with_capacity(buckets * 2);

    for (bucket, samples) in values.chunks(span).enumerate() {
        // Strict comparisons on purpose: a bucket whose samples never differ
        // leaves both ends on the first one, and contributes it alone.
        let (mut lowest, mut highest) = (0, 0);
        for (index, value) in samples.iter().enumerate() {
            if *value < samples[lowest] {
                lowest = index;
            }
            if *value > samples[highest] {
                highest = index;
            }
        }

        let base = bucket * span;
        kept.push(base + lowest.min(highest));
        if lowest != highest {
            kept.push(base + lowest.max(highest));
        }
    }

    kept
}

/// Keeps the samples within `db_below_peak` decibels of the largest one.
///
/// The noise-floor cut: a spectrum's floor is what a display spends most of its
/// points drawing, and forty decibels below the peak is rarely signal. Amplitude
/// decibels, so the threshold is `peak · 10^(-db / 20)`.
///
/// A run whose largest sample is not positive keeps nothing.
#[must_use]
pub fn above_floor(values: &[f64], db_below_peak: f64) -> Vec<usize> {
    let peak = values.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    if peak.is_nan() || peak <= 0.0 {
        return Vec::new();
    }

    let floor = peak * 10.0_f64.powf(-db_below_peak / 20.0);
    values
        .iter()
        .enumerate()
        .filter(|(_, value)| **value >= floor)
        .map(|(index, _)| index)
        .collect()
}

/// The value, or the smallest positive number when it is not one.
///
/// Only the logarithmic reducer needs it, and only so that a zero or negative
/// coordinate has a logarithm rather than a `NaN` that would poison a
/// comparison.
fn positive(at: f64) -> f64 {
    at.max(f64::MIN_POSITIVE)
}

/// Twice the area of the triangle three points make.
fn area(x0: f64, y0: f64, x1: f64, y1: f64, x2: f64, y2: f64) -> f64 {
    ((x0 - x2) * (y1 - y0) - (x0 - x1) * (y2 - y0)).abs()
}

/// The mean point of a range of samples, when the range holds any.
fn mean_of(x: &[f64], y: &[f64], range: std::ops::Range<usize>) -> Option<(f64, f64)> {
    let count = range.len();
    if count == 0 {
        return None;
    }

    let sum_x: f64 = x[range.clone()].iter().sum();
    let sum_y: f64 = y[range].iter().sum();

    Some((sum_x / count as f64, sum_y / count as f64))
}

#[cfg(test)]
mod tests {
    use std::f64::consts::TAU;

    use super::*;

    /// A run of `count` samples of `f`, on a unit-spaced axis.
    fn run(count: usize, f: impl Fn(f64) -> f64) -> (Vec<f64>, Vec<f64>) {
        (0..count)
            .map(|index| (index as f64, f(index as f64)))
            .unzip()
    }

    #[test]
    fn lttb_keeps_the_endpoints_and_the_count_asked_for() {
        let (x, y) = run(1_000, |at| (at * 0.01).sin());

        let kept = lttb(&x, &y, 100);

        assert_eq!(kept.len(), 100);
        assert_eq!(kept[0], 0);
        assert_eq!(kept[99], 999);
    }

    #[test]
    fn lttb_answers_strictly_increasing_indices() {
        let (x, y) = run(1_000, |at| (at * 0.05).sin());

        let kept = lttb(&x, &y, 137);

        assert!(kept.windows(2).all(|pair| pair[0] < pair[1]));
    }

    #[test]
    fn lttb_keeps_a_spike_a_stride_would_miss() {
        let (x, y) = run(100, |at| if (at - 53.0).abs() < 0.5 { 100.0 } else { 0.0 });

        let kept = lttb(&x, &y, 20);

        assert!(
            kept.contains(&53),
            "the spike at 53 was thinned away: {kept:?}"
        );
        assert!(
            !decimate(100, 20).contains(&53),
            "a stride of five steps straight over it"
        );
    }

    #[test]
    fn asking_for_more_than_there_is_keeps_everything() {
        let (x, y) = run(10, |at| at);

        assert_eq!(lttb(&x, &y, 50), (0..10).collect::<Vec<_>>());
        assert_eq!(lttb(&x, &y, 2), (0..10).collect::<Vec<_>>());
        assert_eq!(lttb_log(&x, &y, 50), (0..10).collect::<Vec<_>>());
        assert!(lttb(&[], &[], 10).is_empty());
    }

    #[test]
    fn log_buckets_spend_points_on_the_low_decades() {
        // 1 kHz to 10 MHz: four decades, of which the top one holds nine tenths
        // of the samples.
        let x: Vec<f64> = (1..=10_000).map(|bin| f64::from(bin) * 1_000.0).collect();
        let y: Vec<f64> = x.iter().map(|at| (at.ln() * 0.5).sin()).collect();

        let linear = lttb(&x, &y, 100);
        let logarithmic = lttb_log(&x, &y, 100);

        let below = |kept: &[usize]| kept.iter().filter(|bin| x[**bin] < 1e6).count();

        assert!(below(&logarithmic) > below(&linear));
        assert!(logarithmic.len() <= 100);
        assert_eq!(logarithmic[0], 0);
        assert_eq!(logarithmic[logarithmic.len() - 1], 9_999);
    }

    #[test]
    fn log_buckets_answer_increasing_indices() {
        let x: Vec<f64> = (1..=5_000).map(|bin| f64::from(bin) * 1_000.0).collect();
        let y: Vec<f64> = x.iter().map(|at| 1_000.0 / at.sqrt()).collect();

        let kept = lttb_log(&x, &y, 50);

        assert!(kept.windows(2).all(|pair| pair[0] < pair[1]), "{kept:?}");
    }

    #[test]
    fn an_axis_with_no_logarithmic_span_falls_back_to_linear_buckets() {
        let x = vec![5.0; 100];
        let y: Vec<f64> = (0..100).map(f64::from).collect();

        assert_eq!(lttb_log(&x, &y, 10), lttb(&x, &y, 10));
    }

    #[test]
    fn decimation_strides_and_never_overruns() {
        assert_eq!(decimate(10, 5), vec![0, 2, 4, 6, 8]);
        assert_eq!(decimate(100, 10).len(), 10);
        assert!(decimate(0, 10).is_empty());
        assert!(decimate(10, 0).is_empty());
    }

    #[test]
    fn extrema_keep_both_ends_of_every_bucket() {
        let (_, y) = run(100, |at| (at * 0.1 * TAU).sin());

        let kept = extrema(&y, 10);

        assert!(kept.len() <= 20);
        assert!(kept.windows(2).all(|pair| pair[0] < pair[1]), "{kept:?}");

        let highest = y.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        let lowest = y.iter().copied().fold(f64::INFINITY, f64::min);
        assert!(kept.iter().any(|index| (y[*index] - highest).abs() < 1e-12));
        assert!(kept.iter().any(|index| (y[*index] - lowest).abs() < 1e-12));
    }

    #[test]
    fn a_flat_bucket_contributes_one_sample() {
        assert_eq!(extrema(&[1.0, 1.0, 1.0, 1.0], 2), vec![0, 2]);
        assert!(extrema(&[], 4).is_empty());
        assert!(extrema(&[1.0], 0).is_empty());
    }

    #[test]
    fn the_floor_keeps_what_is_within_the_stated_decibels() {
        let mut values = vec![1.0; 100];
        values[50] = 100.0;

        // Forty decibels below a peak of 100 is 1.0, which the floor keeps.
        assert_eq!(above_floor(&values, 40.0).len(), 100);
        // Twenty decibels below it is 10.0, which only the peak clears.
        assert_eq!(above_floor(&values, 20.0), vec![50]);
    }

    #[test]
    fn a_run_with_no_positive_sample_has_no_floor() {
        assert!(above_floor(&[0.0, 0.0], 40.0).is_empty());
        assert!(above_floor(&[-1.0, -2.0], 40.0).is_empty());
        assert!(above_floor(&[], 40.0).is_empty());
    }
}