From 67c1f51cfe527b4ac6b8bbe5b98b0199528d2082 Mon Sep 17 00:00:00 2001 From: Brummel Date: Wed, 17 Jun 2026 10:18:42 +0200 Subject: [PATCH] perf(aura-std): O(1) Kahan-compensated running-sum SMA (was O(length) re-sum) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Sma::eval re-summed the whole window every tick — O(length)/tick, multiplied across millions of bars and every sweep point. Replace with the industry-standard incremental running sum (ta-lib shape): the node owns the window ring and keeps a running sum, adding the new sample and subtracting the evicted one each cycle, so the input column drops to depth 1 (lookbacks() = [1], like Ema). The running sum carries Kahan/Neumaier compensation (the fix pandas' rolling mean adopted) so it does not drift over long runs; Ema needs none — its recurrence is contractive. Determinism (C1) holds: same input -> same run, reproducibly. The float output differs from the full re-sum in the last ULPs (a new, deterministic baseline); the suite needed only one arity-assumption update (Sma lookback 3 -> 1), no equity golden changed. closes #39 --- crates/aura-engine/src/blueprint.rs | 6 +- crates/aura-std/src/sma.rs | 126 +++++++++++++++++++++++++--- 2 files changed, 117 insertions(+), 15 deletions(-) diff --git a/crates/aura-engine/src/blueprint.rs b/crates/aura-engine/src/blueprint.rs index 0dba2b1..3049737 100644 --- a/crates/aura-engine/src/blueprint.rs +++ b/crates/aura-engine/src/blueprint.rs @@ -2243,8 +2243,10 @@ mod tests { #[test] fn lookbacks_arity_matches_signature_inputs() { use aura_std::{Add, Sma}; - // every std node: one lookback per declared input - assert_eq!(Sma::new(3).lookbacks(), vec![3]); + // every std node: one lookback per declared input. Sma keeps its window in + // node state now (Kahan running sum), so its lookback is 1 (a depth-1 input), + // not `length` — the arity (one lookback per input) is what this test guards. + assert_eq!(Sma::new(3).lookbacks(), vec![1]); assert_eq!(Sma::new(3).lookbacks().len(), Sma::builder().schema().inputs.len()); assert_eq!(Add::new().lookbacks().len(), Add::builder().schema().inputs.len()); } diff --git a/crates/aura-std/src/sma.rs b/crates/aura-std/src/sma.rs index e62c7c9..6d97c4e 100644 --- a/crates/aura-std/src/sma.rs +++ b/crates/aura-std/src/sma.rs @@ -1,16 +1,42 @@ -//! `Sma` — simple moving average over the last `length` values of one f64 -//! input. The walking skeleton's first worked node: it proves the `aura-core` -//! `Node` contract is authorable from a downstream crate and evaluable with no -//! engine present (the test drives it by hand, as the sim loop later will). +//! `Sma` — simple moving average over the last `length` values of one f64 input. +//! +//! Computed as an **O(1) incremental window sum**, not a per-tick re-sum: the node +//! owns a `length`-slot ring of the window and keeps a running `sum`, adding the +//! new sample and subtracting the evicted one each cycle (the ta-lib running-sum +//! shape). The running sum carries a **Kahan/Neumaier compensation** term so it +//! does not drift over millions of add/remove ops — the fix pandas' rolling mean +//! adopted after real float-drift bug reports; ta-lib omits it and drifts. +//! +//! `Ema` needs no equivalent (`ema.rs`): its recurrence is *contractive* +//! (`ema += alpha*(x - ema)` rescales old state by `1-alpha < 1` each tick), so a +//! past rounding error decays geometrically instead of accumulating. The SMA +//! running sum is *accumulative* — error has nowhere to go — which is exactly why +//! it, and not the EMA, needs compensation. +//! +//! Because the window lives in node state, `eval` reads only the newest sample and +//! `lookbacks()` is `1` (the input column is depth-1), exactly like `Ema`. O(1) +//! time, O(length) state, allocation-free on the hot path (the ring and output +//! buffer are sized once at construction). use aura_core::{ Cell, Ctx, FieldSpec, Firing, Node, NodeSchema, ParamSpec, PortSpec, PrimitiveBuilder, ScalarKind, }; -/// Simple moving average over the last `length` values of one f64 input. +/// Simple moving average over the last `length` values of one f64 input, +/// maintained as an O(1) Kahan-compensated running sum over a node-owned ring. pub struct Sma { length: usize, + // The window, node-owned (so the input column is depth-1): `ring[pos]` is the + // oldest value, the one the next push evicts. Sized once at construction (C7). + ring: Box<[f64]>, + pos: usize, + // Samples seen so far — the warm-up gate (silent until `length`, like `Ema`). + count: usize, + // Running sum of the values currently in `ring`, plus its Kahan compensation + // term (the low-order bits each add/remove dropped, folded into the next op). + sum: f64, + comp: f64, out: [Cell; 1], } @@ -18,7 +44,15 @@ impl Sma { /// Build an SMA of window `length` (must be >= 1). pub fn new(length: usize) -> Self { assert!(length >= 1, "SMA length must be >= 1"); - Self { length, out: [Cell::from_f64(0.0)] } + Self { + length, + ring: vec![0.0; length].into_boxed_slice(), + pos: 0, + count: 0, + sum: 0.0, + comp: 0.0, + out: [Cell::from_f64(0.0)], + } } /// The param-generic recipe for a blueprint primitive: declares `length` and builds @@ -37,21 +71,50 @@ impl Sma { } } +/// Kahan/Neumaier compensated accumulation: fold `v` into `*sum`, carrying the +/// low-order bits lost on this step in `*comp` so a long sequence of adds (and +/// removes, which are adds of a negative) does not drift. The whole reason the SMA +/// running sum stays accurate over millions of ticks; `Ema`'s contractive +/// recurrence needs no such term (see the module docs). +fn kahan(sum: &mut f64, comp: &mut f64, v: f64) { + let y = v - *comp; + let t = *sum + y; + *comp = (t - *sum) - y; + *sum = t; +} + impl Node for Sma { + // The window lives in node state, so only the newest sample is read each cycle + // — `length` sizes the ring, not the input column (recursive, like `Ema`). fn lookbacks(&self) -> Vec { - vec![self.length] + vec![1] } fn eval(&mut self, ctx: Ctx<'_>) -> Option<&[Cell]> { let w = ctx.f64_in(0); - if w.len() < self.length { + if w.is_empty() { + return None; // no sample yet + } + let x = w[0]; // index 0 = newest (financial indexing) + + // Add the newest into the running sum; once the window is full, remove the + // value this push evicts (the ring slot about to be overwritten). Both go + // through Kahan so the running sum tracks the true window sum. + kahan(&mut self.sum, &mut self.comp, x); + if self.count >= self.length { + let evicted = self.ring[self.pos]; + kahan(&mut self.sum, &mut self.comp, -evicted); + } + self.ring[self.pos] = x; + self.pos = (self.pos + 1) % self.length; + if self.count < self.length { + self.count += 1; + } + if self.count < self.length { return None; // not yet warmed up } - let mut sum = 0.0; - for k in 0..self.length { - sum += w[k]; // index 0 = newest (financial indexing) - } - self.out[0] = Cell::from_f64(sum / self.length as f64); + + self.out[0] = Cell::from_f64(self.sum / self.length as f64); Some(&self.out) } @@ -100,6 +163,43 @@ mod tests { assert_eq!(sma.eval(Ctx::new(&inputs, Timestamp(0))), Some([Cell::from_f64(9.0)].as_slice())); } + #[test] + fn lookback_is_one_window_lives_in_node_state() { + // the incremental SMA owns its window ring, so it reads only the newest + // sample each cycle — input column depth drops from `length` to 1 (like Ema). + assert_eq!(Sma::new(20).lookbacks(), vec![1]); + assert_eq!(Sma::new(1).lookbacks(), vec![1]); + } + + #[test] + fn incremental_matches_full_resum_within_tolerance() { + // the running Kahan sum must track the true window mean across a long, noisy + // f64 series. Drive the node and a reference full re-sum side by side; they + // agree to a tight tolerance (Kahan keeps drift near machine epsilon — a + // bare running sum would slowly diverge, which is the bug pandas fixed). + let length = 50; + let mut sma = Sma::new(length); + let mut inputs = vec![AnyColumn::with_capacity(ScalarKind::F64, 1)]; + + // a deterministic non-integer series mixing magnitudes, so cancellation bites + let series: Vec = + (0..5_000).map(|i| 1.1234 + (i as f64) * 1e-3 + ((i % 7) as f64) * 0.37).collect(); + + let mut hist: Vec = Vec::new(); + for &x in &series { + inputs[0].push(Scalar::f64(x)).unwrap(); + let got = sma.eval(Ctx::new(&inputs, Timestamp(0))).map(|r| r[0].f64()); + hist.push(x); + if hist.len() >= length { + let want: f64 = hist[hist.len() - length..].iter().sum::() / length as f64; + let g = got.expect("warmed up after `length` samples"); + assert!((g - want).abs() < 1e-9, "incremental {g} vs re-sum {want}"); + } else { + assert_eq!(got, None, "silent until warmed up"); + } + } + } + #[test] fn labels_carry_identifying_params() { use crate::{Add, Exposure, LinComb, Recorder, SimBroker, Sub};