perf(aura-std): O(1) Kahan-compensated running-sum SMA (was O(length) re-sum)
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
This commit is contained in:
@@ -2243,8 +2243,10 @@ mod tests {
|
|||||||
#[test]
|
#[test]
|
||||||
fn lookbacks_arity_matches_signature_inputs() {
|
fn lookbacks_arity_matches_signature_inputs() {
|
||||||
use aura_std::{Add, Sma};
|
use aura_std::{Add, Sma};
|
||||||
// every std node: one lookback per declared input
|
// every std node: one lookback per declared input. Sma keeps its window in
|
||||||
assert_eq!(Sma::new(3).lookbacks(), vec![3]);
|
// 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!(Sma::new(3).lookbacks().len(), Sma::builder().schema().inputs.len());
|
||||||
assert_eq!(Add::new().lookbacks().len(), Add::builder().schema().inputs.len());
|
assert_eq!(Add::new().lookbacks().len(), Add::builder().schema().inputs.len());
|
||||||
}
|
}
|
||||||
|
|||||||
+113
-13
@@ -1,16 +1,42 @@
|
|||||||
//! `Sma` — simple moving average over the last `length` values of one f64
|
//! `Sma` — simple moving average over the last `length` values of one f64 input.
|
||||||
//! 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
|
//! Computed as an **O(1) incremental window sum**, not a per-tick re-sum: the node
|
||||||
//! engine present (the test drives it by hand, as the sim loop later will).
|
//! 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::{
|
use aura_core::{
|
||||||
Cell, Ctx, FieldSpec, Firing, Node, NodeSchema, ParamSpec, PortSpec, PrimitiveBuilder,
|
Cell, Ctx, FieldSpec, Firing, Node, NodeSchema, ParamSpec, PortSpec, PrimitiveBuilder,
|
||||||
ScalarKind,
|
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 {
|
pub struct Sma {
|
||||||
length: usize,
|
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],
|
out: [Cell; 1],
|
||||||
}
|
}
|
||||||
|
|
||||||
@@ -18,7 +44,15 @@ impl Sma {
|
|||||||
/// Build an SMA of window `length` (must be >= 1).
|
/// Build an SMA of window `length` (must be >= 1).
|
||||||
pub fn new(length: usize) -> Self {
|
pub fn new(length: usize) -> Self {
|
||||||
assert!(length >= 1, "SMA length must be >= 1");
|
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
|
/// 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 {
|
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<usize> {
|
fn lookbacks(&self) -> Vec<usize> {
|
||||||
vec![self.length]
|
vec![1]
|
||||||
}
|
}
|
||||||
|
|
||||||
fn eval(&mut self, ctx: Ctx<'_>) -> Option<&[Cell]> {
|
fn eval(&mut self, ctx: Ctx<'_>) -> Option<&[Cell]> {
|
||||||
let w = ctx.f64_in(0);
|
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
|
return None; // not yet warmed up
|
||||||
}
|
}
|
||||||
let mut sum = 0.0;
|
|
||||||
for k in 0..self.length {
|
self.out[0] = Cell::from_f64(self.sum / self.length as f64);
|
||||||
sum += w[k]; // index 0 = newest (financial indexing)
|
|
||||||
}
|
|
||||||
self.out[0] = Cell::from_f64(sum / self.length as f64);
|
|
||||||
Some(&self.out)
|
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()));
|
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<f64> =
|
||||||
|
(0..5_000).map(|i| 1.1234 + (i as f64) * 1e-3 + ((i % 7) as f64) * 0.37).collect();
|
||||||
|
|
||||||
|
let mut hist: Vec<f64> = 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::<f64>() / 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]
|
#[test]
|
||||||
fn labels_carry_identifying_params() {
|
fn labels_carry_identifying_params() {
|
||||||
use crate::{Add, Exposure, LinComb, Recorder, SimBroker, Sub};
|
use crate::{Add, Exposure, LinComb, Recorder, SimBroker, Sub};
|
||||||
|
|||||||
Reference in New Issue
Block a user