improvise/crates/improvise-core/examples/retraction_drift.rs

//! How far a SUM maintained by insert/retract deltas drifts from the eager fold,
//! in f64. ADR 009 (improvise-felo.4), "Floating point".

use improvise_core::format::format_f64;
use improvise_core::model::cell::{CellKey, CellValue};
use improvise_core::model::types::{Aggregation, Model};

const MEASURE: &str = "_Measure";

/// SplitMix64: seeded, so a run is reproducible and a second seed is a new sample.
struct Rng(u64);

impl Rng {
    fn next(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }

    fn below(&mut self, n: u64) -> u64 {
        self.next() % n
    }
}

/// A user-entered decimal as an exact count of 1e-5 units, its magnitude spread
/// over `decades` decades. The f64 the store holds is `typed(units)`.
fn decimal(rng: &mut Rng, decades: u32) -> i64 {
    let scale = 10i64.pow(rng.below(u64::from(decades)) as u32);
    let cents = rng.below(100_000) as i64;
    let sign = if rng.below(4) == 0 { -1 } else { 1 };
    sign * cents * scale
}

fn typed(units: i64) -> f64 {
    units as f64 / 1e5
}

/// The exact decimal sum, rounded to f64 once: what ADR 005's arithmetic would show.
fn exact(units: &[i64]) -> f64 {
    units.iter().map(|&u| i128::from(u)).sum::<i128>() as f64 / 1e5
}

fn key(row: usize) -> CellKey {
    CellKey::new(vec![
        ("Row".to_string(), format!("r{row}")),
        (MEASURE.to_string(), "Amount".to_string()),
    ])
}

/// `cells` rows of one measure; the eager SUM is the one the grid draws.
fn document(rng: &mut Rng, cells: usize, decades: u32) -> (Model, Vec<i64>) {
    let mut m = Model::new("drift");
    m.add_category("Row").expect("under the limit");
    let values: Vec<i64> = (0..cells).map(|_| decimal(rng, decades)).collect();
    for (row, v) in values.iter().enumerate() {
        m.category_mut("Row")
            .expect("just added")
            .add_item(format!("r{row}"));
        m.set_cell(key(row), CellValue::Number(typed(*v)));
    }
    (m, values)
}

fn eager(m: &Model) -> f64 {
    let total = CellKey::new(vec![(MEASURE.to_string(), "Amount".to_string())]);
    m.evaluate_aggregated_f64(&total, Aggregation::OverHidden)
}

/// Neumaier's compensated running sum: the cheap mitigation, measured beside the naive one.
#[derive(Default)]
struct Compensated {
    sum: f64,
    carry: f64,
}

impl Compensated {
    fn add(&mut self, x: f64) {
        let t = self.sum + x;
        self.carry += match self.sum.abs() >= x.abs() {
            true => (self.sum - t) + x,
            false => (x - t) + self.sum,
        };
        self.sum = t;
    }

    fn value(&self) -> f64 {
        self.sum + self.carry
    }
}

fn ulps(a: f64, b: f64) -> u64 {
    (a.to_bits() as i64).abs_diff(b.to_bits() as i64)
}

/// Each running sum against the eager fold (the differential oracle) and against the
/// exact decimal sum. `flip`: the gap reaches `FLOAT_EQ_EPSILON`, so `=` can answer differently.
fn report(label: &str, cycles: usize, exact: f64, eager: f64, naive: f64, comp: f64) {
    let shown = |v: f64| format_f64(v, true, 2);
    println!(
        "{label:<22} {cycles:>8} {:>10.2e} {:>6} {:>10.2e} {:>6} {:>5} {:>10.2e} {:>10.2e} {:>10.2e}",
        (naive - eager).abs(),
        ulps(naive, eager),
        (comp - eager).abs(),
        ulps(comp, eager),
        (naive - eager).abs() >= 1e-10,
        (eager - exact).abs(),
        (naive - exact).abs(),
        (comp - exact).abs(),
    );
    if shown(naive) != shown(eager) || shown(eager) != shown(exact) {
        println!(
            "  shown at ,.2: naive {} eager {} exact {}",
            shown(naive),
            shown(eager),
            shown(exact)
        );
    }
}

/// Each cycle retracts one live value and inserts a fresh one at the same key.
fn random_updates(seed: u64, cells: usize, decades: u32, checkpoints: &[usize]) {
    let mut rng = Rng(seed);
    let (mut m, mut values) = document(&mut rng, cells, decades);
    let mut naive = eager(&m);
    let mut comp = Compensated::default();
    comp.add(naive);
    let mut done = 0;
    for &at in checkpoints {
        while done < at {
            let row = rng.below(cells as u64) as usize;
            let new = decimal(&mut rng, decades);
            naive = naive - typed(values[row]) + typed(new);
            comp.add(-typed(values[row]));
            comp.add(typed(new));
            values[row] = new;
            m.set_cell(key(row), CellValue::Number(typed(new)));
            done += 1;
        }
        let label = format!("{cells} cells, {decades} dec");
        report(&label, at, exact(&values), eager(&m), naive, comp.value());
    }
}

/// A large value arrives, small ones accumulate under it, then it is retracted.
fn absorbed_by_a_large_value() {
    let mut m = Model::new("absorb");
    m.add_category("Row").expect("under the limit");
    m.category_mut("Row").expect("just added").add_item("r0");
    let mut naive = 0.0;
    let mut comp = Compensated::default();
    m.set_cell(key(0), CellValue::Number(1e17));
    naive += 1e17;
    comp.add(1e17);
    for row in 1..=1000 {
        m.category_mut("Row")
            .expect("just added")
            .add_item(format!("r{row}"));
        m.set_cell(key(row), CellValue::Number(0.01));
        naive += 0.01;
        comp.add(0.01);
    }
    m.clear_cell(&key(0));
    naive -= 1e17;
    comp.add(-1e17);
    report("1e17, 1000 x 0.01, -1e17", 1001, 10.0, eager(&m), naive, comp.value());
}

fn main() {
    println!(
        "{:<22} {:>8} {:>10} {:>6} {:>10} {:>6} {:>5} {:>10} {:>10} {:>10}",
        "case",
        "cycles",
        "naive-eag",
        "ulp",
        "comp-eag",
        "ulp",
        "flip",
        "eag-exact",
        "nai-exact",
        "cmp-exact"
    );
    let checkpoints = [1_000, 10_000, 100_000, 1_000_000];
    for seed in [1, 2] {
        random_updates(seed, 1_000, 3, &checkpoints);
        random_updates(seed, 18_000, 3, &checkpoints);
        random_updates(seed, 18_000, 9, &checkpoints);
    }
    absorbed_by_a_large_value();
}