fix(incremental): Avg rounds the exact mean once; core pins its exactness; Meter counts sum updates (improvise-hrjz)

Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01VdfEGXA7aJ6ftMsCc9PJze

crates/improvise-core/src/model/types.rs

@@ -18,9 +18,18 @@ const FLOAT_EQ_EPSILON: f64 = 1e-10;
 /// The sum correctly rounded once, whatever the order (ADR 009 clause 8), signed as IEEE
 /// signs a zero total.
 fn exact_sum(values: &[f64]) -> f64 {
+    exact(values).value()
+}
+
+/// The mean correctly rounded once: the exact total over the count, not a rounded sum divided.
+fn exact_mean(values: &[f64]) -> f64 {
+    exact(values).mean(values.len() as u64)
+}
+
+fn exact(values: &[f64]) -> numtower::ExactSum {
     let mut sum = numtower::ExactSum::new();
     values.iter().for_each(|v| sum.add(*v, 1));
-    sum.value()
+    sum
 }
 
 /// Reduce raw values. `None` on empty: an empty cell reads as **absent**, not
@@ -28,7 +37,7 @@ fn exact_sum(values: &[f64]) -> f64 {
 fn apply_agg(values: &[f64], func: &AggFunc) -> Option<f64> {
     (!values.is_empty()).then(|| match func {
         AggFunc::Sum => exact_sum(values),
-        AggFunc::Avg => exact_sum(values) / values.len() as f64,
+        AggFunc::Avg => exact_mean(values),
         AggFunc::Min => values.iter().copied().fold(f64::INFINITY, f64::min),
         AggFunc::Max => values.iter().copied().fold(f64::NEG_INFINITY, f64::max),
         AggFunc::Count => values.len() as f64,
@@ -3615,3 +3624,88 @@ mod route_then_aggregate_tests {
         );
     }
 }
+
+/// Sums and averages are exact in the model itself (ADR 009 clause 8), whatever backend reads it.
+#[cfg(test)]
+mod exact_sum_tests {
+    use super::*;
+    use crate::formula::parse_formula;
+
+    #[test]
+    fn a_sum_does_not_depend_on_the_order_of_its_addends() {
+        for values in [[1e16, 1.0, -1e16], [1e16, -1e16, 1.0], [1.0, 1e16, -1e16]] {
+            assert_eq!(apply_agg(&values, &AggFunc::Sum), Some(1.0), "{values:?}");
+        }
+    }
+
+    #[test]
+    fn an_average_rounds_the_exact_mean_once_and_does_not_overflow() {
+        assert_eq!(
+            apply_agg(&[f64::MAX, f64::MAX], &AggFunc::Avg),
+            Some(f64::MAX)
+        );
+        let values = [0.1, 0.2, 0.3];
+        let mut sum = numtower::ExactSum::new();
+        values.iter().for_each(|v| sum.add(*v, 1));
+        assert_eq!(apply_agg(&values, &AggFunc::Avg), Some(sum.mean(3)));
+    }
+
+    fn model_with(values: &[f64], formula: &str) -> Model {
+        let mut m = Model::new("m");
+        m.add_category("Region").unwrap();
+        m.category_mut("_Measure").unwrap().add_item("Revenue");
+        for (i, v) in values.iter().enumerate() {
+            let region = format!("r{i}");
+            m.category_mut("Region").unwrap().add_item(&region);
+            let key = CellKey::new(vec![
+                ("Region".into(), region),
+                ("_Measure".into(), "Revenue".into()),
+            ]);
+            m.set_cell(key, CellValue::Number(*v));
+        }
+        m.add_formula(parse_formula(formula, "_Measure").unwrap());
+        m.recompute_formulas(&["Region".to_string()]);
+        m
+    }
+
+    fn total(m: &Model) -> Option<CellValue> {
+        let at = CellKey::new(vec![("_Measure".into(), "Tot".into())]);
+        m.evaluate_aggregated(&at, Aggregation::OverHidden)
+    }
+
+    #[test]
+    fn a_formula_sum_is_exact_whatever_the_order_of_the_records() {
+        let m = model_with(&[1e16, 1.0, -1e16], "Tot = SUM(Revenue * 1)");
+        assert_eq!(total(&m), Some(CellValue::Number(1.0)));
+        let m = model_with(&[-1e16, 1e16, 1.0], "Tot = SUM(Revenue * 1)");
+        assert_eq!(total(&m), Some(CellValue::Number(1.0)));
+    }
+
+    #[test]
+    fn a_formula_sum_over_nothing_is_positive_zero() {
+        let m = model_with(&[1.0], "Tot = SUM(Revenue WHERE Region = Nowhere)");
+        let Some(CellValue::Number(n)) = total(&m) else {
+            panic!("the formula has a value");
+        };
+        assert_eq!(
+            n.to_bits(),
+            0.0f64.to_bits(),
+            "it was -0.0 under Sum for f64"
+        );
+    }
+
+    #[test]
+    fn a_stored_group_sums_exactly_and_averages_once() {
+        let mut m = model_with(&[1e16, 1.0, -1e16], "Tot = SUM(Revenue * 1)");
+        let at = CellKey::new(vec![("_Measure".into(), "Revenue".into())]);
+        m.measure_agg.insert("Revenue".into(), AggFunc::Sum);
+        let sum = m.evaluate_aggregated(&at, Aggregation::OverHidden);
+        assert_eq!(sum, Some(CellValue::Number(1.0)));
+        m.measure_agg.insert("Revenue".into(), AggFunc::Avg);
+        let Some(CellValue::Number(got)) = m.evaluate_aggregated(&at, Aggregation::OverHidden)
+        else {
+            panic!("an average")
+        };
+        assert_eq!(got.to_bits(), (1.0f64 / 3.0).to_bits());
+    }
+}

crates/improvise-incremental/src/group.rs

@@ -124,7 +124,7 @@ fn reduce(groups: &[&Group], func: Reduce, cells: i64) -> f64 {
     match func {
         Reduce::Count => cells as f64,
         Reduce::Sum => total(groups).value(),
-        Reduce::Avg => total(groups).value() / cells as f64,
+        Reduce::Avg => total(groups).mean(cells as u64),
         Reduce::Min => {
             let smallest = groups.iter().filter_map(|g| g.smallest());
             let min = smallest.fold(f64::INFINITY, f64::min);

crates/improvise-incremental/src/shapes.rs

@@ -12,6 +12,8 @@ use std::collections::BTreeMap;
 pub struct Meter {
     /// Group updates applied from a delta.
     pub updates: u64,
+    /// Exact-sum additions and retractions applied from the change log.
+    pub sum_updates: u64,
     /// Groups looked up to answer a read.
     pub reads: u64,
     /// Shapes built from the store.
@@ -33,6 +35,7 @@ impl Meter {
     pub fn plus(self, o: Meter) -> Meter {
         Meter {
             updates: self.updates + o.updates,
+            sum_updates: self.sum_updates + o.sum_updates,
             reads: self.reads + o.reads,
             builds: self.builds + o.builds,
             built_cells: self.built_cells + o.built_cells,
@@ -68,9 +71,10 @@ fn same_bits(change: &CellChange) -> bool {
 
 impl Shape {
     /// The canonical delta moves counts and multisets; the raw changes move sums and the sign of
-    /// a zero, which tell the two zeros apart. Returns the group updates the delta caused.
-    fn apply(&mut self, cats: &[String], delta: &CellDelta, changes: &[CellChange]) -> u64 {
-        let mut updates = 0;
+    /// a zero, which tell the two zeros apart. Returns the group updates the delta caused
+    /// and the sum updates the changes caused.
+    fn apply(&mut self, cats: &[String], delta: &CellDelta, changes: &[CellChange]) -> (u64, u64) {
+        let (mut updates, mut sums) = (0, 0);
         let mut touched = Vec::new();
         for ((key, value), weight) in delta.iter() {
             let Some(at) = project(key, cats) else {
@@ -92,9 +96,11 @@ impl Shape {
             };
             if let Some(CellValue::Number(old)) = &change.old {
                 group.add_number(*old, -1);
+                sums += 1;
             }
             if let Some(CellValue::Number(new)) = &change.new {
                 group.add_number(*new, 1);
+                sums += 1;
             }
             touched.push(at);
         }
@@ -103,7 +109,7 @@ impl Shape {
                 self.groups.remove(&at);
             }
         }
-        updates
+        (updates, sums)
     }
 
     /// The shape of `store`'s cells, one pass.
@@ -157,13 +163,15 @@ impl Shapes {
 
     /// Feeds a batch of changes, and the delta they make, to every shape built so far.
     pub fn apply(&mut self, delta: &CellDelta, changes: &[CellChange]) {
-        let updates: u64 = self
-            .built
-            .get_mut()
-            .iter_mut()
-            .map(|(cats, shape)| shape.apply(cats, delta, changes))
-            .sum();
-        self.bump(|m| m.updates += updates);
+        let (mut updates, mut sums) = (0, 0);
+        for (cats, shape) in self.built.get_mut().iter_mut() {
+            let (u, s) = shape.apply(cats, delta, changes);
+            (updates, sums) = (updates + u, sums + s);
+        }
+        self.bump(|m| {
+            m.updates += updates;
+            m.sum_updates += sums;
+        });
     }
 
     /// Runs `read` on the groups at `keys` in the shape of their categories, building the

crates/improvise-incremental/tests/exact.rs

@@ -68,12 +68,20 @@ fn sums_are_exact_whatever_the_order_the_cells_were_added_in() {
 }
 
 #[test]
-fn an_average_is_the_exact_sum_over_the_count() {
+fn an_average_is_the_exact_total_over_the_count_rounded_once() {
     let values = [0.1, 0.2, 0.3, 1e16, -1e16];
-    let want = Some((exact(&values) / values.len() as f64).to_bits());
+    let mut sum = ExactSum::new();
+    values.iter().for_each(|v| sum.add(*v, 1));
+    let want = Some(sum.mean(values.len() as u64).to_bits());
     assert_eq!(both(&values, AggFunc::Avg), (want, want));
 }
 
+#[test]
+fn an_average_does_not_overflow_with_its_sum_in_either_backend() {
+    let want = Some(f64::MAX.to_bits());
+    assert_eq!(both(&[f64::MAX, f64::MAX], AggFunc::Avg), (want, want));
+}
+
 #[test]
 fn a_zero_total_has_the_sign_ieee_gives_the_same_addends() {
     let cases: [(&[f64], f64); 4] = [
@@ -111,7 +119,7 @@ fn overwriting_zero_with_negative_zero_moves_the_sign_of_a_sum() {
 }
 
 #[test]
-fn a_million_overwrites_in_a_thousand_cell_group_do_not_drift() {
+fn twenty_thousand_overwrites_in_a_thousand_cell_group_do_not_drift() {
     let values: Vec<f64> = (0..1000).map(|i| 0.1 * f64::from(i % 7 + 1)).collect();
     let model = model_of(&values, AggFunc::Sum);
     let mut eager = Eager::new(model.clone(), vec![]);

crates/improvise-incremental/tests/measure.rs

@@ -61,6 +61,7 @@ fn product_partial() -> CellKey {
 fn since(before: Meter, after: Meter) -> Meter {
     Meter {
         updates: after.updates - before.updates,
+        sum_updates: after.sum_updates - before.sum_updates,
         reads: after.reads - before.reads,
         builds: after.builds - before.builds,
         built_cells: after.built_cells - before.built_cells,
@@ -94,7 +95,7 @@ fn a_group_read_comes_from_a_maintained_group() {
 
 /// The operations one overwrite costs once `shapes` category sets are maintained.
 fn overwrite_cost(n: usize, shapes: usize, func: AggFunc) -> (Meter, usize) {
-    let mut backend = Incremental::new(table(n, func), vec![]);
+    let mut backend = Incremental::new(table(n, func.clone()), vec![]);
     let partials = [region_partial(), product_partial()];
     for partial in partials.iter().take(shapes) {
         backend.raw_fold(partial);
@@ -102,7 +103,12 @@ fn overwrite_cost(n: usize, shapes: usize, func: AggFunc) -> (Meter, usize) {
     let before = backend.meter();
     backend.apply(&Edit::SetCell(cell(1, 2, 3), CellValue::Number(2000.0)));
     let cost = since(before, backend.meter());
-    assert!(backend.raw_fold(&region_partial()).is_some());
+    let read = backend.raw_fold(&region_partial());
+    let fresh = Incremental::new(backend.model().clone(), vec![]);
+    assert_eq!(read, fresh.raw_fold(&region_partial()), "{func:?}");
+    if matches!(func, AggFunc::Max) {
+        assert_eq!(read, Some(CellValue::Number(2000.0)));
+    }
     (cost, backend.shape_count())
 }
 
@@ -121,6 +127,10 @@ fn an_overwrite_costs_two_updates_per_shape_whatever_the_table_size() {
             let (large, _) = overwrite_cost(12, shapes, func.clone());
             assert!(built >= shapes);
             let updates = 2 * shapes as u64;
+            assert_eq!(
+                small.sum_updates, updates,
+                "{func:?}: and as many exact-sum updates"
+            );
             assert_eq!(
                 small.updates, updates,
                 "{func:?}: retract and insert, per shape"

crates/numtower/src/sum.rs

@@ -155,6 +155,23 @@ impl ExactSum {
         }
     }
 
+    /// The exact sum divided by `count`, rounded once: a mean that neither overflows with the
+    /// sum nor rounds twice. Non-finite sums and the sign of a zero read as `value` does.
+    /// Panics on a zero `count`.
+    pub fn mean(&self, count: u64) -> f64 {
+        assert!(count > 0, "the mean of nothing");
+        match (self.nan != 0, self.pos_inf != 0, self.neg_inf != 0) {
+            (false, false, false) => match (self.negative_zero(), self.total()) {
+                (true, total) if total.is_zero() => -0.0,
+                (_, total) => {
+                    let unit = BigInt::one() << 1074;
+                    ratio_to_f64(&total, &(unit * BigInt::from(count)))
+                }
+            },
+            _ => self.value(),
+        }
+    }
+
     /// The finite part as an integer in units of 2^-1074.
     fn total(&self) -> BigInt {
         let mut limbs = self.limbs;

crates/numtower/tests/mean.rs

@@ -0,0 +1,47 @@
+//! `ExactSum::mean`: the exact total over a count, rounded once, against exact rationals.
+//! Evidence: `context/test-evidence/numtower.md`.
+
+mod oracle;
+
+use num_bigint::BigInt;
+use num_rational::BigRational;
+use numtower::ExactSum;
+use oracle::{finite_f64, is_nearest};
+use proptest::prelude::*;
+
+fn sum_of(values: &[f64]) -> ExactSum {
+    let mut sum = ExactSum::new();
+    values.iter().for_each(|v| sum.add(*v, 1));
+    sum
+}
+
+#[test]
+fn a_mean_whose_sum_would_overflow_is_still_finite() {
+    let sum = sum_of(&[f64::MAX, f64::MAX]);
+    assert_eq!(sum.mean(2), f64::MAX, "the sum is inf, the mean is not");
+}
+
+#[test]
+fn non_finite_sums_have_non_finite_means_and_zero_keeps_its_sign() {
+    assert_eq!(sum_of(&[f64::INFINITY, 1.0]).mean(2), f64::INFINITY);
+    assert_eq!(sum_of(&[f64::NEG_INFINITY]).mean(3), f64::NEG_INFINITY);
+    assert!(sum_of(&[f64::NAN, 1.0]).mean(2).is_nan());
+    assert!(sum_of(&[f64::INFINITY, f64::NEG_INFINITY]).mean(2).is_nan());
+    assert_eq!(sum_of(&[-0.0, -0.0]).mean(2).to_bits(), (-0.0f64).to_bits());
+    assert_eq!(sum_of(&[1.0, -1.0]).mean(2).to_bits(), 0.0f64.to_bits());
+}
+
+proptest! {
+    #[test]
+    fn the_mean_is_the_nearest_f64_to_the_exact_total_over_the_count(
+        values in prop::collection::vec(finite_f64(), 1..9),
+    ) {
+        let total = values
+            .iter()
+            .map(|v| BigRational::from_float(*v).unwrap())
+            .fold(BigRational::from_integer(BigInt::from(0)), |a, b| a + b);
+        let exact = total / BigRational::from_integer(BigInt::from(values.len()));
+        let got = sum_of(&values).mean(values.len() as u64);
+        prop_assert!(is_nearest(got, &exact), "{got} for {values:?}");
+    }
+}