From d8e151a10a9a510986017eb6e88a355ca739ee3a Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:28:46 -0700 Subject: [PATCH 01/14] Leave long RR artifact runs as gaps Lipponen & Tarvainen (2019), Peltola (2012). The 32-beat cap is an internal ESTIMATE rule, not a published threshold. Include synthetic input recipes and old/new measurements. --- ALGORITHMS.md | 6 + CHANGELOG.md | 12 + lib/src/onehz/foundations/rr_correction.dart | 213 ++++++---- .../foundations/rr_correction_kernel.dart | 373 ++++++++++++++++++ test/fixtures/output_changes.json | 134 +++++++ test/output_changes/rr_test.dart | 16 + tool/reproduce_outputs.dart | 62 +++ 7 files changed, 742 insertions(+), 74 deletions(-) create mode 100644 CHANGELOG.md create mode 100644 lib/src/onehz/foundations/rr_correction_kernel.dart create mode 100644 test/fixtures/output_changes.json create mode 100644 test/output_changes/rr_test.dart create mode 100644 tool/reproduce_outputs.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 593ec67a..751061d9 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -158,3 +158,9 @@ al. 2017 · Pimentel et al. · Brage et al. 2004 · Keytel et al. 2005 · Tudor- 2019 · Hopkins 2004 · Mahalanobis 1936 · Killick et al. 2012 · Wittmann & Roenneberg 2006 · Pietilä et al. 2018 · Malik et al. + + + +| Batch change | Method and limits | +|---|---| +| Leave long RR artifact runs as gaps | Lipponen & Tarvainen 2019; Peltola 2012; The 32-beat cap is an internal ESTIMATE rule with low confidence, not a published threshold. Longer runs split the series; each side is corrected independently. The synthetic 33-beat case keeps 200 NN beats and drops 33. | diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 00000000..080e707b --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,12 @@ +# Changelog + +## Split #99 output changes + +Edge must bump `kAlgoVersion` when adopting these changed batch outputs. + +### Leave long RR artifact runs as gaps + +The 32-beat cap is an internal ESTIMATE rule with low confidence, not a published threshold. Longer runs split the series; each side is corrected independently. The synthetic 33-beat case keeps 200 NN beats and drops 33. + +Method: Lipponen & Tarvainen 2019; Peltola 2012. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/foundations/rr_correction.dart b/lib/src/onehz/foundations/rr_correction.dart index 89840130..41e69169 100644 --- a/lib/src/onehz/foundations/rr_correction.dart +++ b/lib/src/onehz/foundations/rr_correction.dart @@ -1,3 +1,11 @@ +// Lipponen & Tarvainen 2019 classification with Peltola 2012 correction policy. +// The unpublished artifact-run cap is an internal rule, ESTIMATE with low +// confidence: +// over-cap runs become gaps and each side is re-corrected. Without a cap split +// the result equals uncapped correction. Shared kernel arithmetic and raw seam +// context preserve capped batch NN/classes; settled stream plus provisional tail +// equals capped batch. RrSandwich preserves source-beat absolute timestamps; +// uncapped arbitrarily long artifact runs have no exactness guarantee. // FOUNDATION — RR artifact correction. // // Lipponen & Tarvainen 2019 ("A robust algorithm for heart rate variability @@ -13,8 +21,14 @@ // This is PRV. Output = cleaned NN series + per-beat artifact mask + clean // fraction. Multi-beat gaps are never silently bridged. -import 'dart:math' as math; -import '../util.dart'; +import 'rr_correction_kernel.dart'; + +/// Longest artifact run (beats) the corrector will carry. A longer run is not a +/// correction problem but an absence: it breaks the series into independent +/// segments and its beats are left out (a gap). 32 beats leaves no run on the +/// owner's 18 days of data over the cap (longest seen: 16) and bounds the +/// streaming state; see combined-bump/brief.md. +const int kMaxArtifactRun = 32; /// Artifact label for a single beat. enum BeatClass { normal, ectopic, longShort, missed, extra } @@ -34,6 +48,13 @@ class RrCorrectionResult { /// spliced out; pass the timestamps. final List nnTimesMs; + /// For each NN, the index of the input beat it came from: the beat itself, + /// or, for a spline-corrected beat, that beat's own index. Strictly + /// increasing; the input beats not listed here were dropped. Lets a caller + /// pair each NN with that beat's own (absolute) timestamp. Empty on a + /// hand-built result that did not track it. + final List nnSourceIndex; + /// Per-input-beat classification (same length as the input rr). final List classes; @@ -53,6 +74,7 @@ class RrCorrectionResult { required this.cleanFraction, required this.droppedCount, required this.correctedCount, + this.nnSourceIndex = const [], }); } @@ -74,6 +96,10 @@ class RrCorrectionResult { /// `rec_ts*1000` and therefore whole-second — and re-anchors to the real /// timestamp whenever the wall step exceeds the interval by more than /// [reanchorGapMs], i.e. at a dropout. +/// +/// [maxArtifactRun] caps a run of artifact beats: a longer run splits the series +/// (output = correctRr(before) ++ correctRr(after), times absolute, the run's +/// beats absent). null = uncapped, the pre-cap behaviour. RrCorrectionResult correctRr( List rrMs, { List? rrTsMs, @@ -81,7 +107,55 @@ RrCorrectionResult correctRr( int windowBeats = 91, double minThresholdMs = 100, double reanchorGapMs = 1000, + int? maxArtifactRun = kMaxArtifactRun, }) { + if (maxArtifactRun != null && maxArtifactRun < 1) { + throw ArgumentError.value(maxArtifactRun, 'maxArtifactRun'); + } + if (windowBeats < 1) throw ArgumentError.value(windowBeats, 'windowBeats'); + if (maxArtifactRun != null) { + // The overwhelmingly common path has no run over the cap. Reuse the + // original vectorized correction verbatim, including its clock, and only + // evaluate split windows when a run actually requires a cut. + var hardRun = 0, hardGap = false; + for (final rr in rrMs) { + hardRun = rr < 300 || rr > 2000 ? hardRun + 1 : 0; + if (hardRun > maxArtifactRun) { hardGap = true; break; } + } + if (!hardGap) { + final ordinary = correctRr(rrMs, rrTsMs: rrTsMs, alpha: alpha, + windowBeats: windowBeats, minThresholdMs: minThresholdMs, + reanchorGapMs: reanchorGapMs, maxArtifactRun: null); + var run = 0, gap = false; + for (final cls in ordinary.classes) { + run = cls == BeatClass.normal ? 0 : run + 1; + if (run > maxArtifactRun) { gap = true; break; } + } + if (!gap) return ordinary; + } + final w = RrCapWindow.evaluate(rrMs, + ts: rrTsMs, cap: maxArtifactRun, alpha: alpha, + windowBeats: windowBeats, minThresholdMs: minThresholdMs, + reanchorGapMs: reanchorGapMs); + final nn = [], times = [], source = []; + var normal = 0, dropped = 0, corrected = 0; + for (var i = 0; i < rrMs.length; i++) { + if (w.classes[i] == BeatClass.normal) normal++; + final v = w.values[i]; + if (v == null) { + dropped++; + } else { + nn.add(v); + times.add(w.times[i]); + source.add(i); + if (w.classes[i] != BeatClass.normal) corrected++; + } + } + return RrCorrectionResult(nn: nn, nnTimesMs: times, classes: w.classes, + cleanFraction: rrMs.isEmpty ? 0 : normal / rrMs.length, + droppedCount: dropped, correctedCount: corrected, + nnSourceIndex: source); + } final n = rrMs.length; // Beat-end time for EVERY input beat, kept or not. One array, built once, so // every branch below reads a clock instead of advancing its own. @@ -105,10 +179,12 @@ RrCorrectionResult correctRr( ]; final nn = []; final times = []; + final source = []; for (var i = 0; i < n; i++) { if (classes[i] == BeatClass.normal) { nn.add(rrMs[i]); times.add(tAll[i]); + source.add(i); } } return RrCorrectionResult( @@ -118,6 +194,7 @@ RrCorrectionResult correctRr( cleanFraction: nn.length / n, droppedCount: n - nn.length, correctedCount: 0, + nnSourceIndex: source, ); } @@ -133,32 +210,13 @@ RrCorrectionResult correctRr( // medRR: rr minus local median (for missed/extra long-range tests). final med = _slidingMedian(rrMs, windowBeats); - final mRR = List.generate(n, (i) { - final d = rrMs[i] - med[i]; - return d < 0 ? d * 2 : d; // paper asymmetry weight - }); + final mRR = List.generate(n, (i) => medianDeviation(rrMs[i], med[i])); final th2 = _slidingThreshold(mRR, windowBeats, alpha, minThresholdMs); - final classes = List.filled(n, BeatClass.normal); - for (var i = 0; i < n; i++) { - final hardLong = rrMs[i] > 2000; - final hardShort = rrMs[i] < 300; - final bigJump = dRR[i].abs() > th1[i]; - final bigDev = mRR[i].abs() > th2[i]; - if (hardLong || (bigDev && mRR[i] > 0)) { - // Long interval: likely a MISSED beat (interval ~ multiple of normal). - classes[i] = (med[i] > 0 && rrMs[i] > 1.5 * med[i]) - ? BeatClass.missed - : BeatClass.longShort; - } else if (hardShort || (bigDev && mRR[i] < 0)) { - // Short interval: likely an EXTRA (spurious) beat. - classes[i] = (med[i] > 0 && rrMs[i] < 0.6 * med[i]) - ? BeatClass.extra - : BeatClass.longShort; - } else if (bigJump) { - classes[i] = BeatClass.ectopic; - } - } + final classes = [ + for (var i = 0; i < n; i++) + classifyBeat(rrMs[i], dRR[i], th1[i], med[i], mRR[i], th2[i]) + ]; // Compensatory-pair reconciliation. A single ectopic/extra beat shows up as // TWO successive dRR spikes of opposite sign (the bad beat, then the @@ -169,17 +227,13 @@ RrCorrectionResult correctRr( // single isolated artifact rather than a spurious 2-beat run. for (var k = 1; k < n; k++) { if (classes[k] != BeatClass.ectopic) continue; - final prevBad = classes[k - 1] == BeatClass.extra || - classes[k - 1] == BeatClass.missed || - classes[k - 1] == BeatClass.ectopic || - classes[k - 1] == BeatClass.longShort; - if (!prevBad) continue; - final oppositeSign = dRR[k] * dRR[k - 1] < 0; - final valueNormal = med[k] > 0 && - rrMs[k] >= 300 && - rrMs[k] <= 2000 && - (rrMs[k] - med[k]).abs() <= 0.2 * med[k]; - if (oppositeSign && valueNormal) { + final prevBad = classes[k - 1] != BeatClass.normal; + if (isRecoveryBeat( + prevIsArtifact: prevBad, + rr: rrMs[k], + dRR: dRR[k], + dRRPrev: dRR[k - 1], + med: med[k])) { classes[k] = BeatClass.normal; } } @@ -189,6 +243,7 @@ RrCorrectionResult correctRr( final isArtifact = [for (final c in classes) c != BeatClass.normal]; final nn = []; final times = []; + final source = []; var dropped = 0; var corrected = 0; var i = 0; @@ -196,6 +251,7 @@ RrCorrectionResult correctRr( if (!isArtifact[i]) { nn.add(rrMs[i]); times.add(tAll[i]); + source.add(i); i++; continue; } @@ -216,6 +272,7 @@ RrCorrectionResult correctRr( // — 1.6 % of a night at one missed beat in sixty. nn.add(corr); times.add(tAll[i]); + source.add(i); corrected++; } else { dropped++; // no anchors -> honest drop; the clock already elapsed @@ -241,6 +298,7 @@ RrCorrectionResult correctRr( cleanFraction: normalCount / n, droppedCount: dropped, correctedCount: corrected, + nnSourceIndex: source, ); } @@ -306,42 +364,63 @@ List _slidingThreshold( final n = x.length; final out = List.filled(n, 0); final half = win ~/ 2; + final w = _SortedWindow(x, half); for (var i = 0; i < n; i++) { - final lo = math.max(0, i - half); - final hi = math.min(n - 1, i + half); - final seg = []; - for (var k = lo; k <= hi; k++) { - seg.add(x[k]); - } - final q1 = percentile(seg, 25) ?? 0; - final q3 = percentile(seg, 75) ?? 0; - final qd = (q3 - q1) / 2; - // Floor keeps a gross outlier detectable on (near-)quantized clean data - // where the QD genuinely collapses to 0 (constant RR). On any series with - // real beat-to-beat variability α·QD dominates the floor, so the floor - // never governs a physiological signal. - out[i] = math.max(alpha * qd, floor); + w.advanceTo(i); + out[i] = thresholdOfSorted(w.a, alpha, floor); } return out; } +/// Median of the window around each beat EXCLUDING the beat itself (falls back +/// to the beat when the window holds nothing else). Same value as sorting the +/// window-minus-self per beat, but the window is kept sorted as it slides and +/// the self element is skipped by index arithmetic instead of copied out. List _slidingMedian(List x, int win) { final n = x.length; final out = List.filled(n, 0); final half = win ~/ 2; + final w = _SortedWindow(x, half); for (var i = 0; i < n; i++) { - final lo = math.max(0, i - half); - final hi = math.min(n - 1, i + half); - final seg = []; - for (var k = lo; k <= hi; k++) { - if (k == i) continue; - seg.add(x[k]); - } - out[i] = median(seg) ?? x[i]; + w.advanceTo(i); + out[i] = medianExcluding(w.a, x[i]) ?? x[i]; } return out; } +/// Ascending-sorted copy of `x[max(0,i-half) .. min(n-1,i+half)]`, maintained +/// incrementally as [advanceTo] steps i forward by one (O(log w) search plus a +/// w-element memmove per beat instead of an O(w log w) sort). +class _SortedWindow { + final List _x; + final int _half; + final _set = SortedMultiset(); + List get a => _set.a; + int _lo = 0; + int _hi = -1; // inclusive; -1 = empty + + _SortedWindow(this._x, this._half); + + void advanceTo(int i) { + final n = _x.length; + final newLo = i - _half < 0 ? 0 : i - _half; + final newHi = i + _half > n - 1 ? n - 1 : i + _half; + if (_hi < _lo) { + // First call: build once. + _set.resetFrom(_x.getRange(newLo, newHi + 1)); + } else { + while (_lo < newLo) { + _set.remove(_x[_lo++]); + } + while (_hi < newHi) { + _set.add(_x[++_hi]); + } + } + _lo = newLo; + _hi = newHi; + } +} + /// Catmull-Rom cubic interpolation at the artifact index using the nearest two /// NORMAL beats on each side. Returns null if anchors are unavailable. double? _splineCorrect(List rr, List isArtifact, int idx) { @@ -353,19 +432,5 @@ double? _splineCorrect(List rr, List isArtifact, int idx) { for (var k = idx + 1; k < rr.length && right.length < 2; k++) { if (!isArtifact[k]) right.add(rr[k]); } - if (left.isEmpty || right.isEmpty) return null; - final p1 = left.last; - final p2 = right.first; - final p0 = left.length >= 2 ? left.first : p1; - final p3 = right.length >= 2 ? right.last : p2; - // Catmull-Rom at t=0.5 between p1 and p2. - const t = 0.5; - final t2 = t * t; - final t3 = t2 * t; - final v = 0.5 * - ((2 * p1) + - (-p0 + p2) * t + - (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 + - (-p0 + 3 * p1 - 3 * p2 + p3) * t3); - return v; + return splineMid(left, right); } diff --git a/lib/src/onehz/foundations/rr_correction_kernel.dart b/lib/src/onehz/foundations/rr_correction_kernel.dart new file mode 100644 index 00000000..8acb72b5 --- /dev/null +++ b/lib/src/onehz/foundations/rr_correction_kernel.dart @@ -0,0 +1,373 @@ +// Lipponen & Tarvainen 2019 classification with Peltola 2012 correction policy. +// The unpublished artifact-run cap is an internal rule, ESTIMATE with low +// confidence: +// over-cap runs become gaps and each side is re-corrected. Without a cap split +// the result equals uncapped correction. Shared kernel arithmetic and raw seam +// context preserve capped batch NN/classes; settled stream plus provisional tail +// equals capped batch. RrSandwich preserves source-beat absolute timestamps; +// uncapped arbitrarily long artifact runs have no exactness guarantee. +// FOUNDATION — the per-beat pieces of Lipponen–Tarvainen correction that both +// `correctRr` (whole series) and `RrCorrector` (streaming) run. +// +// One copy of the sorted window, the threshold, the leave-one-out median, the +// class decision, the compensatory-pair rule and the spline, so the batch and +// the stream cannot drift apart: the stream's contract is bit-identity with +// `correctRr`, and sharing the arithmetic is what keeps that cheap to hold. +// +// Not exported from `onehz.dart` — internal to foundations/. + +import 'dart:math' as math; +import '../util.dart'; +import 'rr_correction.dart' show BeatClass; + +/// First index whose element is not less than [v] under `compareTo` — the same +/// total order `List.sort()` uses (-0.0 < 0.0, NaN last), so a window +/// kept in this order matches the one a per-beat sort produced. +int lowerBound(List a, double v) { + var lo = 0; + var hi = a.length; + while (lo < hi) { + final m = (lo + hi) >> 1; + if (a[m].compareTo(v) < 0) { + lo = m + 1; + } else { + hi = m; + } + } + return lo; +} + +/// An ascending multiset of doubles: O(log w) search plus a w-element memmove +/// per add / remove, instead of an O(w log w) sort per beat. +class SortedMultiset { + /// The values, ascending. Read-only to callers. + final List a = []; + + void add(double v) => a.insert(lowerBound(a, v), v); + + /// Removes one occurrence of [v], which must be present. + void remove(double v) => a.removeAt(lowerBound(a, v)); + + void resetFrom(Iterable xs) { + a + ..clear() + ..addAll(xs) + ..sort(); + } +} + +/// `alpha × QD` of an ascending window, floored: QD = (Q3 − Q1)/2. +/// +/// Floor keeps a gross outlier detectable on (near-)quantized clean data where +/// the QD genuinely collapses to 0 (constant RR). On any series with real +/// beat-to-beat variability α·QD dominates the floor. +double thresholdOfSorted(List sorted, double alpha, double floor) { + final q1 = percentileSorted(sorted, 25) ?? 0; + final q3 = percentileSorted(sorted, 75) ?? 0; + final qd = (q3 - q1) / 2; + return math.max(alpha * qd, floor); +} + +/// `median(sorted minus one occurrence of v)` — [percentile]'s linear +/// interpolation at p = 50 over the sorted list with v's slot skipped, in the +/// same operation order so the result is bit-identical. Null if nothing is left. +double? medianExcluding(List sorted, double v) { + final len = sorted.length - 1; + if (len <= 0) return null; + final skip = lowerBound(sorted, v); + double at(int j) => j < skip ? sorted[j] : sorted[j + 1]; + if (len == 1) return at(0); + final rank = (50.0 / 100) * (len - 1); + final lo = rank.floor(); + final hi = rank.ceil(); + if (lo == hi) return at(lo); + final frac = rank - lo; + return at(lo) + (at(hi) - at(lo)) * frac; +} + +/// The paper's asymmetric deviation from the local median: below-median +/// deviations count double. +double medianDeviation(double rr, double med) { + final d = rr - med; + return d < 0 ? d * 2 : d; +} + +/// First-pass class of one beat from its dRR, mRR and the two thresholds. +BeatClass classifyBeat( + double rr, double dRR, double th1, double med, double mRR, double th2) { + final hardLong = rr > 2000; + final hardShort = rr < 300; + final bigJump = dRR.abs() > th1; + final bigDev = mRR.abs() > th2; + if (hardLong || (bigDev && mRR > 0)) { + // Long interval: likely a MISSED beat (interval ~ multiple of normal). + return (med > 0 && rr > 1.5 * med) ? BeatClass.missed : BeatClass.longShort; + } else if (hardShort || (bigDev && mRR < 0)) { + // Short interval: likely an EXTRA (spurious) beat. + return (med > 0 && rr < 0.6 * med) ? BeatClass.extra : BeatClass.longShort; + } else if (bigJump) { + return BeatClass.ectopic; + } + return BeatClass.normal; +} + +/// Compensatory-pair rule: beat k was flagged ONLY by its dRR jump, the beat +/// before it was an artifact, the two dRR spikes have opposite sign and k's own +/// value is normal and close to the local median — k is just the recovery from +/// the previous beat's event, so it is demoted to normal. +bool isRecoveryBeat({ + required bool prevIsArtifact, + required double rr, + required double dRR, + required double dRRPrev, + required double med, +}) { + if (!prevIsArtifact) return false; + final oppositeSign = dRR * dRRPrev < 0; + final valueNormal = + med > 0 && rr >= 300 && rr <= 2000 && (rr - med).abs() <= 0.2 * med; + return oppositeSign && valueNormal; +} + +/// Catmull-Rom at t = 0.5 between the nearest NORMAL beats either side of an +/// isolated artifact. [left] / [right] hold up to two values each, nearest to +/// the artifact LAST in [left] and FIRST in [right]. Null without an anchor on +/// both sides. +double? splineMid(List left, List right) { + if (left.isEmpty || right.isEmpty) return null; + final p1 = left.last; + final p2 = right.first; + final p0 = left.length >= 2 ? left.first : p1; + final p3 = right.length >= 2 ? right.last : p2; + const t = 0.5; + final t2 = t * t; + final t3 = t2 * t; + return 0.5 * + ((2 * p1) + + (-p0 + p2) * t + + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 + + (-p0 + 3 * p1 - 3 * p2 + p3) * t3); +} + +/// Internal bounded-window evaluator shared by the capped batch and stream. +/// [from] has its complete left window in [rr], or is the segment's first beat. +/// Only classes at and after [from] are evaluated; the preceding class and +/// normal anchors are carried explicitly, so truncating old context cannot +/// change recovery or spline decisions. +class RrCapWindow { + final List classes; + final List values; + final List times; + final List gaps; + final List hardGaps; + RrCapWindow._(int n) + : classes = List.filled(n, BeatClass.normal), + values = List.filled(n, null), + times = List.filled(n, 0), + gaps = List.filled(n, false), + hardGaps = List.filled(n, false); + + static RrCapWindow evaluate( + List rr, { + List? ts, + required int cap, + double alpha = 5.2, + int windowBeats = 91, + double minThresholdMs = 100, + double reanchorGapMs = 1000, + int from = 0, + BeatClass previous = BeatClass.normal, + List leftNormals = const [], + bool leadingGap = false, + bool leadingHardGap = false, + int segmentStart = 0, + double clock = 0, + double? origin, + }) { + final out = RrCapWindow._(rr.length); + if (rr.isEmpty) return out; + final wall = ts != null && ts.length == rr.length; + final base = origin ?? (wall ? ts[0] - rr[0] : 0.0); + final h = windowBeats ~/ 2; + + List classify(int a, int b, int start, BeatClass prev) { + final n = b - a; + if (n < 3) { + return [ + for (var i = start; i < b; i++) + rr[i] >= 300 && rr[i] <= 2000 + ? BeatClass.normal + : BeatClass.longShort + ]; + } + final med = List.filled(n, 0); + final mrr = List.filled(n, 0); + final th1 = List.filled(n, 0); + final wd = SortedMultiset(), wr = SortedMultiset(); + var lo = a, hi = a; + double d(int i) => i == a ? 0 : rr[i] - rr[i - 1]; + for (var i = a; i < b; i++) { + final nextLo = math.max(a, i - h), nextHi = math.min(b, i + h + 1); + while (hi < nextHi) { + wd.add(d(hi)); + wr.add(rr[hi++]); + } + while (lo < nextLo) { + wd.remove(d(lo)); + wr.remove(rr[lo++]); + } + med[i - a] = medianExcluding(wr.a, rr[i]) ?? rr[i]; + mrr[i - a] = medianDeviation(rr[i], med[i - a]); + th1[i - a] = thresholdOfSorted(wd.a, alpha, minThresholdMs); + } + final wm = SortedMultiset(); + lo = a; + hi = a; + final cl = []; + for (var i = start; i < b; i++) { + final nextLo = math.max(a, i - h), nextHi = math.min(b, i + h + 1); + if (i == start) lo = nextLo; + if (i == start) hi = nextLo; + while (hi < nextHi) wm.add(mrr[hi++ - a]); + while (lo < nextLo) wm.remove(mrr[lo++ - a]); + var c = classifyBeat(rr[i], d(i), th1[i - a], med[i - a], mrr[i - a], + thresholdOfSorted(wm.a, alpha, minThresholdMs)); + if (i > a && + c == BeatClass.ectopic && + isRecoveryBeat( + prevIsArtifact: prev != BeatClass.normal, + rr: rr[i], + dRR: d(i), + dRRPrev: d(i - 1), + med: med[i - a])) { + c = BeatClass.normal; + } + cl.add(c); + prev = c; + } + return cl; + } + + // A stack avoids recursion depth proportional to the number of gaps. + final pending = <(int, int, int, BeatClass, List, bool, bool)>[ + ( + segmentStart, + rr.length, + from, + previous, + leftNormals, + leadingGap, + leadingHardGap + ) + ]; + while (pending.isNotEmpty) { + final (a, b, start, prev, anchors, inGap, hardLeading) = + pending.removeLast(); + if (start >= b) continue; + if (inGap && hardLeading && rr[start] >= 300 && rr[start] <= 2000) { + pending + .add((start, b, start, BeatClass.normal, const [], false, false)); + continue; + } + var hardCut = false; + var cutA = -1, cutB = -1; + // Physiologically impossible intervals are artifacts independently of + // every dispersion window. Cut these first: otherwise their medians can + // absorb valid beats on either side into the detected run itself. + for (var i = start; i < b;) { + if (rr[i] >= 300 && rr[i] <= 2000) { + i++; + continue; + } + var j = i + 1; + while (j < b && (rr[j] < 300 || rr[j] > 2000)) { + j++; + } + if (j - i > cap || (inGap && i == start)) { + cutA = i; + cutB = j; + hardCut = true; + break; + } + i = j; + } + final cl = cutA >= 0 ? [] : classify(a, b, start, prev); + if (cutA < 0 && inGap && cl.first == BeatClass.normal) { + pending + .add((start, b, start, BeatClass.normal, const [], false, false)); + continue; + } + if (cutA < 0) { + for (var i = start; i < b;) { + if (cl[i - start] == BeatClass.normal) { + i++; + continue; + } + var j = i + 1; + while (j < b && cl[j - start] != BeatClass.normal) { + j++; + } + if (j - i > cap || (inGap && i == start)) { + cutA = i; + cutB = j; + break; + } + i = j; + } + } + if (cutA >= 0) { + for (var i = cutA; i < cutB; i++) { + out.classes[i] = BeatClass.longShort; + out.gaps[i] = true; + out.hardGaps[i] = hardCut || (inGap && hardLeading && cutA == start); + } + // Cut BEFORE classifying either side, including their dispersion + // windows. The gap's labels are canonical, never provisional medians. + pending.add((cutB, b, cutB, BeatClass.normal, const [], false, false)); + pending.add((a, cutA, start, prev, anchors, inGap, hardLeading)); + continue; + } + for (var i = start; i < b; i++) { + out.classes[i] = cl[i - start]; + } + final normals = List.of(anchors); + for (var i = start; i < b;) { + if (out.classes[i] == BeatClass.normal) { + out.values[i] = rr[i]; + normals.add(rr[i]); + if (normals.length > 2) normals.removeAt(0); + i++; + continue; + } + var j = i + 1; + while (j < b && out.classes[j] != BeatClass.normal) { + j++; + } + if (j == i + 1) { + final right = []; + for (var k = j; k < b && right.length < 2; k++) { + if (out.classes[k] == BeatClass.normal) right.add(rr[k]); + } + out.values[i] = splineMid(normals, right); + } + i = j; + } + } + var cur = clock; + var gap = leadingGap; + for (var i = from; i < rr.length; i++) { + final reset = gap && !out.gaps[i]; + final dropout = + wall && i > 0 && ts[i] - ts[i - 1] - rr[i] > reanchorGapMs; + final anchored = wall ? ts[i] - base : 0.0; + cur = wall && reset + ? anchored + : dropout && anchored > cur + ? anchored + : cur + rr[i]; + out.times[i] = cur; + gap = out.gaps[i]; + } + return out; + } +} diff --git a/test/fixtures/output_changes.json b/test/fixtures/output_changes.json new file mode 100644 index 00000000..f3c232f5 --- /dev/null +++ b/test/fixtures/output_changes.json @@ -0,0 +1,134 @@ +{ + "provenance": "Synthetic deterministic arithmetic cases, no personal data. Formula inputs are expanded by tool/reproduce_outputs.dart. Old values measured on up/main 3a7ca34; new values measured on f8e7d3c plus #98.", + "rr": { + "normal": 900.0, + "before": 100, + "run": 33, + "artifact": 250.0, + "after": 100 + }, + "spectral": { + "count": 96, + "origin": 1800000000.0, + "step": 0.8125, + "frequency": 0.137 + }, + "cosinor": { + "count": 96, + "step": 0.7 + }, + "near_singular": { + "count": 20000, + "origin": 3.0, + "step": 1e-09 + }, + "hrv": { + "count": 1800 + }, + "load": { + "count": 90 + }, + "zones": { + "count": 120 + }, + "circadian": { + "days": 3, + "epochs": 24 + }, + "old": { + "rr": [ + 199, + 34, + 0 + ], + "spectral": 32627.260306094682, + "cosinor": [ + 60.20591440999805, + 1.9937577719142656, + 0.9902626829318173, + 2.900200339561918e-94 + ], + "near_singular": true, + "hrv": [ + 3.615265234174367, + 21.52663924532639, + 1.067383333873884, + 21.55214616158443, + 0.9849552962327551 + ], + "nonfinite": "throws", + "trimp": 32.71197112260296, + "ctl": [ + 98.49620732954641, + 103.12182624417179, + -4.625618914625377 + ], + "zones": [ + 0.0, + 120.01199999999982, + 0.0, + 0.0, + 0.0 + ], + "circadian": [ + 1.0, + 0.331267605633807, + 3.505000000000011, + 1.3600000000000079, + 14, + 0 + ], + "irregular_gap": [ + 22.360679774997894, + 83.33333333333333 + ] + }, + "new": { + "rr": [ + 200, + 33, + 0 + ], + "spectral": 32627.260527635714, + "cosinor": [ + 60.205914409998044, + 1.9937577719142632, + 0.9902626829318162, + 2.9002003395774255e-94 + ], + "near_singular": false, + "hrv": [ + 3.615265234174367, + 21.52663924532641, + 1.0673833338741983, + 21.55214616158443, + 0.984955296232754 + ], + "nonfinite": false, + "trimp": 32.711971122602975, + "ctl": [ + 98.49620732954642, + 103.12182624417177, + -4.625618914625349 + ], + "zones": [ + 0.0, + 120.01199999999977, + 0.0, + 0.0, + 0.0 + ], + "circadian": [ + 1.0, + 0.3312676056338072, + 3.505000000000011, + 1.3600000000000079, + 14, + 0 + ], + "irregular_gap": [ + 22.360679774997894, + 83.33333333333333 + ] + } +} diff --git a/test/output_changes/rr_test.dart b/test/output_changes/rr_test.dart new file mode 100644 index 00000000..1421dbd6 --- /dev/null +++ b/test/output_changes/rr_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Leave long RR artifact runs as gaps: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['rr']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/tool/reproduce_outputs.dart b/tool/reproduce_outputs.dart new file mode 100644 index 00000000..ea320f7c --- /dev/null +++ b/tool/reproduce_outputs.dart @@ -0,0 +1,62 @@ +// Synthetic reproduction cases. Run with --json to capture unrounded values. +// Uses only public batch functions, so the same file runs on up/main and #99. +import 'dart:convert'; +import 'dart:io'; +import 'dart:math' as math; +import 'package:openstrap_analytics/onehz.dart'; +import 'package:openstrap_analytics/src/onehz/util.dart' as u; + +Map reproduce(Map f) { + final out = {}; + final r = f['rr'] as Map; + final rr = [for (var i=0;i<(r['before'] as int);i++) (r['normal'] as num).toDouble(), + for (var i=0;i<(r['run'] as int);i++) (r['artifact'] as num).toDouble(), + for (var i=0;i<(r['after'] as int);i++) (r['normal'] as num).toDouble()]; + final clean = correctRr(rr); + out['rr'] = [clean.nn.length, clean.droppedCount, clean.correctedCount]; + final s=f['spectral'] as Map; + final t=[for(var i=0;i<(s['count'] as int);i++) (s['origin'] as num).toDouble()+i*(s['step'] as num)]; + final y=[for(var i=0;i; + final ct=[for(var i=0;i<(c['count'] as int);i++) i*(c['step'] as num).toDouble()]; + final cy=[for(var i=0;i; + final nt=[for(var i=0;i<(n['count'] as int);i++) (n['origin'] as num).toDouble()+i*(n['step'] as num)]; + out['near_singular']=cosinor(nt,[for(var i=0;iv.toDouble()).toList(), + nnTimesMs:[800,1700,2700,3800,14500,15350,16250,17250],minBeats:4).value!; + out['irregular_gap']=[ir.sd1,ir.pnnPct]; + return out; +} +void main(List args) { + final f=jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final values=reproduce(f); + if(args.contains('--json')) { print(jsonEncode(values)); return; } + print('| Change | up/main 3a7ca34 | Current batch |'); + print('|---|---|---|'); + final old=f['old'] as Map; + for(final entry in values.entries) { + print('| ${entry.key} | ${old[entry.key]} | ${entry.value} |'); + } +} From cdb51240b49557fe2e2d7231bceef594441d35d6 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:28:58 -0700 Subject: [PATCH 02/14] Reduce spectral error on epoch timestamps Lomb (1976), Scargle (1982), Press & Rybicki (1989). Shift the spectral origin to the first sample before forming trigonometric arguments; pin the unrounded power on synthetic epoch seconds. --- ALGORITHMS.md | 6 ++++++ CHANGELOG.md | 6 ++++++ lib/src/onehz/util.dart | 7 ++++--- test/output_changes/spectral_test.dart | 16 ++++++++++++++++ 4 files changed, 32 insertions(+), 3 deletions(-) create mode 100644 test/output_changes/spectral_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 751061d9..2ac318f0 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -164,3 +164,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Leave long RR artifact runs as gaps | Lipponen & Tarvainen 2019; Peltola 2012; The 32-beat cap is an internal ESTIMATE rule with low confidence, not a published threshold. Longer runs split the series; each side is corrected independently. The synthetic 33-beat case keeps 200 NN beats and drops 33. | + + + +| Batch change | Method and limits | +|---|---| +| Reduce spectral error on epoch timestamps | Lomb 1976; Scargle 1982; Press & Rybicki 1989; Evaluate trigonometric arguments relative to the first sample. The mathematical spectrum is time-shift invariant; floating-point band powers and threshold decisions can change for all Lomb-Scargle callers. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 080e707b..97bf1499 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -10,3 +10,9 @@ The 32-beat cap is an internal ESTIMATE rule with low confidence, not a publishe Method: Lipponen & Tarvainen 2019; Peltola 2012. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Reduce spectral error on epoch timestamps + +Evaluate trigonometric arguments relative to the first sample. The mathematical spectrum is time-shift invariant; floating-point band powers and threshold decisions can change for all Lomb-Scargle callers. + +Method: Lomb 1976; Scargle 1982; Press & Rybicki 1989. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/util.dart b/lib/src/onehz/util.dart index 634880dd..c847b5d2 100644 --- a/lib/src/onehz/util.dart +++ b/lib/src/onehz/util.dart @@ -345,6 +345,7 @@ LombScargle? lombScargle(List t, List y, List freqsHz) { } final span = tMax - tMin; if (span <= 0) return null; + final t0 = t[0]; // One-sided periodogram → density: 2·Δt, with Δt the MEAN sample interval. final psdScale = 2.0 * span / (n - 1); @@ -358,14 +359,14 @@ LombScargle? lombScargle(List t, List y, List freqsHz) { // τ: phase reference for time-shift invariance. var sin2 = 0.0, cos2 = 0.0; for (final ti in t) { - sin2 += math.sin(2 * w * ti); - cos2 += math.cos(2 * w * ti); + sin2 += math.sin(2 * w * (ti - t0)); + cos2 += math.cos(2 * w * (ti - t0)); } final tau = math.atan2(sin2, cos2) / (2 * w); var cNum = 0.0, cDen = 0.0, sNum = 0.0, sDen = 0.0; for (var i = 0; i < n; i++) { - final arg = w * (t[i] - tau); + final arg = w * ((t[i] - t0) - tau); final c = math.cos(arg); final s = math.sin(arg); cNum += yc[i] * c; diff --git a/test/output_changes/spectral_test.dart b/test/output_changes/spectral_test.dart new file mode 100644 index 00000000..01e2fce5 --- /dev/null +++ b/test/output_changes/spectral_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Reduce spectral error on epoch timestamps: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['spectral']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From b7423f8815db36808f8af118d38486c592ef78b0 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:30:14 -0700 Subject: [PATCH 03/14] Abstain when a beat has no finite time Preserve Task Force (1996) and Welch (1967) formulas for finite inputs. Timestamp validity is an implementation guard, preventing endless segment walks and invalid integer conversions. --- ALGORITHMS.md | 6 +++++ CHANGELOG.md | 6 +++++ lib/src/onehz/clinical/hrv_freq.dart | 7 +++++ lib/src/onehz/clinical/hrv_time.dart | 23 +++++++++++++++++ lib/src/onehz/respiration/cvhr_apnea.dart | 7 +++++ lib/src/onehz/respiration/resp_rate.dart | 7 +++++ test/output_changes/nonfinite_test.dart | 31 +++++++++++++++++++++++ 7 files changed, 87 insertions(+) create mode 100644 test/output_changes/nonfinite_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 2ac318f0..a7519de1 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -170,3 +170,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Reduce spectral error on epoch timestamps | Lomb 1976; Scargle 1982; Press & Rybicki 1989; Evaluate trigonometric arguments relative to the first sample. The mathematical spectrum is time-shift invariant; floating-point band powers and threshold decisions can change for all Lomb-Scargle callers. | + + + +| Batch change | Method and limits | +|---|---| +| Abstain when a beat has no finite time | Task Force 1996 HRV; Welch 1967 spectra; input-validity guard is an implementation rule; HRV frequency/time, nocturnal and session RMSSD, RSA and CVHR return absent metrics for non-finite beat times. This prevents endless segment walks and invalid integer conversions. Finite-input formulas are unchanged. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 97bf1499..51a9c454 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -16,3 +16,9 @@ Evaluate trigonometric arguments relative to the first sample. The mathematical Method: Lomb 1976; Scargle 1982; Press & Rybicki 1989. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Abstain when a beat has no finite time + +HRV frequency/time, nocturnal and session RMSSD, RSA and CVHR return absent metrics for non-finite beat times. This prevents endless segment walks and invalid integer conversions. Finite-input formulas are unchanged. + +Method: Task Force 1996 HRV; Welch 1967 spectra; input-validity guard is an implementation rule. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/clinical/hrv_freq.dart b/lib/src/onehz/clinical/hrv_freq.dart index d7592486..2190a1d3 100644 --- a/lib/src/onehz/clinical/hrv_freq.dart +++ b/lib/src/onehz/clinical/hrv_freq.dart @@ -84,6 +84,13 @@ Metric hrvFreq( // Times in seconds for Hz output; tachogram value = NN (ms). final tSec = [for (final t in nnTimesMs) t / 1000.0]; + if (tSec.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { return const Metric.absent( diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index d2b35104..cd1a5658 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -632,6 +632,15 @@ Metric hrvTime( ); } + if (nnTimesMs != null && + nnTimesMs.length == nnMs.length && + nnTimesMs.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } // RMSSD / pNN50: root mean square of SUCCESSIVE differences — successive in // TIME, not merely adjacent in the compacted list. correctRr drops multi-beat // artifact runs while advancing its clock across them, so nn[i-1] and nn[i] @@ -818,6 +827,13 @@ Metric nocturnalRmssd( note: 'too few NN intervals for windowed nocturnal RMSSD', ); } + if (nnTimesMs.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } final t0 = nnTimesMs.first; // Bucket beat INDICES by window index, so each window keeps its beat times and // the successive differences below can skip the ones that straddle a dropped @@ -1042,6 +1058,13 @@ Metric sleepSessionRmssdDetail( ); } + if (rrTsMs.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } final buckets = >{}; final bucketsTs = >{}; final inRr = []; diff --git a/lib/src/onehz/respiration/cvhr_apnea.dart b/lib/src/onehz/respiration/cvhr_apnea.dart index fe1ada73..e4046b8b 100644 --- a/lib/src/onehz/respiration/cvhr_apnea.dart +++ b/lib/src/onehz/respiration/cvhr_apnea.dart @@ -127,6 +127,13 @@ Metric cvhrApneaScreen( } final tSec = [for (final t in nnTimesMs) t / 1000.0]; + if (tSec.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } if (tSec.last - tSec.first <= 0) { return const Metric.absent( tier: Tier.high, diff --git a/lib/src/onehz/respiration/resp_rate.dart b/lib/src/onehz/respiration/resp_rate.dart index 6cd3452a..377b0ee3 100644 --- a/lib/src/onehz/respiration/resp_rate.dart +++ b/lib/src/onehz/respiration/resp_rate.dart @@ -187,6 +187,13 @@ Metric rsaRespRate( ); } final tSec = [for (final t in nnTimesMs) t / 1000.0]; + if (tSec.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { return const Metric.absent( diff --git a/test/output_changes/nonfinite_test.dart b/test/output_changes/nonfinite_test.dart new file mode 100644 index 00000000..f26f6a1d --- /dev/null +++ b/test/output_changes/nonfinite_test.dart @@ -0,0 +1,31 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import 'package:openstrap_analytics/onehz.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + for (final bad in [double.nan, double.infinity, double.negativeInfinity]) { + test('all batch timestamp guards abstain on $bad', () { + final nn = [for (var i = 0; i < 1500; i++) 900.0 + 40 * (i % 7)]; + final ts = [for (var i = 0; i < nn.length; i++) i * 1100.0]; + ts[700] = bad; + for (final metric in [hrvTime(nn, nnTimesMs: ts), hrvFreq(nn, ts, artifactFraction: 0), + rsaRespRate(nn, ts, artifactFraction: 0), cvhrApneaScreen(nn, ts, artifactFraction: 0), nocturnalRmssd(nn, ts), + sleepSessionRmssdDetail(nn, ts, startSec: 1, endSec: 7200)]) { + expect(metric.present, isFalse); + expect(metric.note, 'non-finite beat time'); + } + }); + } + + test('Abstain when a beat has no finite time: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['nonfinite']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From f54d922674a7a806e4cf461defec59b25717a599 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:30:26 -0700 Subject: [PATCH 04/14] Refuse unidentifiable cosinor fits Nelson et al. (1979), Methods for cosinor-rhythmometry; Welford (1962) centered co-moments. A relative determinant gate replaces the absolute pivot test; preserve upstream significance policy. --- ALGORITHMS.md | 6 + CHANGELOG.md | 6 + lib/src/onehz/clinical/cosinor.dart | 173 +++++++++++++------------- test/output_changes/cosinor_test.dart | 16 +++ 4 files changed, 113 insertions(+), 88 deletions(-) create mode 100644 test/output_changes/cosinor_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index a7519de1..b397fb3d 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -176,3 +176,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Abstain when a beat has no finite time | Task Force 1996 HRV; Welch 1967 spectra; input-validity guard is an implementation rule; HRV frequency/time, nocturnal and session RMSSD, RSA and CVHR return absent metrics for non-finite beat times. This prevents endless segment walks and invalid integer conversions. Finite-input formulas are unchanged. | + + + +| Batch change | Method and limits | +|---|---| +| Refuse unidentifiable cosinor fits | Nelson et al. 1979, Methods for cosinor-rhythmometry; Welford 1962 centered co-moments; Use centered least squares and a relative determinant gate of 1e-12 of the diagonal product. This is a numerical policy. Nearly collinear phase columns abstain, and ordinary fit parameters can move in their last digits. Upstream zero-amplitude significance gating is retained. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 51a9c454..46723fe2 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,3 +22,9 @@ HRV frequency/time, nocturnal and session RMSSD, RSA and CVHR return absent metr Method: Task Force 1996 HRV; Welch 1967 spectra; input-validity guard is an implementation rule. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Refuse unidentifiable cosinor fits + +Use centered least squares and a relative determinant gate of 1e-12 of the diagonal product. This is a numerical policy. Nearly collinear phase columns abstain, and ordinary fit parameters can move in their last digits. Upstream zero-amplitude significance gating is retained. + +Method: Nelson et al. 1979, Methods for cosinor-rhythmometry; Welford 1962 centered co-moments. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/clinical/cosinor.dart b/lib/src/onehz/clinical/cosinor.dart index 6d4ddfab..595ccfa4 100644 --- a/lib/src/onehz/clinical/cosinor.dart +++ b/lib/src/onehz/clinical/cosinor.dart @@ -60,14 +60,6 @@ class CosinorFit { /// reported. const double cosinorAcrophaseAlpha = 0.05; -/// Upper tail of F(2, d): P(F > f) = (1 + 2f/d)^(−d/2), exact for 2 numerator -/// degrees of freedom. -double _f2Survival(double f, int d) { - if (f.isNaN) return 1.0; - if (f.isInfinite) return 0.0; - if (f <= 0) return 1.0; - return math.pow(1 + 2 * f / d, -d / 2).toDouble(); -} /// Minimum samples for a cosinor fit. The model has 3 free parameters /// (M, β, γ); with n = 4 there is a single residual degree of freedom, so R² @@ -86,74 +78,102 @@ Metric cosinor( List y, { double periodHours = 24, int minPoints = cosinorMinPoints, +}) { + if (tHours.length != y.length) { + return const Metric.absent(tier: Tier.high, + inputs_used: ['signal_timeseries'], note: 'too few points for cosinor'); + } + final state = _CosinorState(); + final w = 2 * math.pi / periodHours; + for (var i = 0; i < y.length; i++) { + state.addCentered(math.cos(w * tHours[i]), math.sin(w * tHours[i]), y[i]); + } + return _finishCosinor(state, periodHours: periodHours, minPoints: minPoints); +} + +// Nelson et al. (1979), Methods for cosinor-rhythmometry. +// Centered least squares changes rounding. Relative 1e-12 gate is a +// numerical identifiability policy, not a physiological cutoff. +class _CosinorState { + var n = 0; + var meanC = 0.0; + var meanS = 0.0; + var meanY = 0.0; + var mCC = 0.0; + var mSS = 0.0; + var mYY = 0.0; + var mCS = 0.0; + var mCY = 0.0; + var mSY = 0.0; + + void addCentered(double c, double s, double y) { + final nextN = n + 1; + final dc = c - meanC; + final ds = s - meanS; + final dy = y - meanY; + final nextMeanC = meanC + dc / nextN; + final nextMeanS = meanS + ds / nextN; + final nextMeanY = meanY + dy / nextN; + mCC += dc * (c - nextMeanC); + mSS += ds * (s - nextMeanS); + mYY += dy * (y - nextMeanY); + mCS += dc * (s - nextMeanS); + mCY += dc * (y - nextMeanY); + mSY += ds * (y - nextMeanY); + meanC = nextMeanC; + meanS = nextMeanS; + meanY = nextMeanY; + n = nextN; + } + +} + +Metric _finishCosinor( + _CosinorState state, { + required double periodHours, + required int minPoints, }) { const inputs = ['signal_timeseries']; - final n = y.length; - if (n < math.max(4, minPoints) || tHours.length != n || periodHours <= 0) { + final n = state.n; + if (n < math.max(4, minPoints) || periodHours <= 0) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'too few points for cosinor', ); } - final w = 2 * math.pi / periodHours; - final c = [for (final t in tHours) math.cos(w * t)]; - final s = [for (final t in tHours) math.sin(w * t)]; - - // Normal equations for [M, β, γ] of y = M + β·c + γ·s. - // Build 3×3 symmetric system. - var sC = 0.0, sS = 0.0, sCC = 0.0, sSS = 0.0, sCS = 0.0; - var sY = 0.0, sYC = 0.0, sYS = 0.0; - for (var i = 0; i < n; i++) { - sC += c[i]; - sS += s[i]; - sCC += c[i] * c[i]; - sSS += s[i] * s[i]; - sCS += c[i] * s[i]; - sY += y[i]; - sYC += y[i] * c[i]; - sYS += y[i] * s[i]; - } - // Matrix A = [[n, sC, sS],[sC, sCC, sCS],[sS, sCS, sSS]], b=[sY,sYC,sYS]. - final a = [ - [n.toDouble(), sC, sS], - [sC, sCC, sCS], - [sS, sCS, sSS], - ]; - final b = [sY, sYC, sYS]; - final sol = _solve3(a, b); - if (sol == null) { + final determinant = state.mCC * state.mSS - state.mCS * state.mCS; + // The determinant is relative to the diagonal co-moments. A design needs + // at least 1e-12 of that scale after removing its means to identify both + // cosine coefficients; below that, numerical noise decides the fit. + if (determinant <= 1e-12 * state.mCC * state.mSS) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'singular cosinor design (insufficient phase coverage)', ); } - final mesor = sol[0]; - final beta = sol[1]; - final gamma = sol[2]; + final beta = + (state.mCY * state.mSS - state.mSY * state.mCS) / determinant; + final gamma = + (state.mSY * state.mCC - state.mCY * state.mCS) / determinant; + final mesor = state.meanY - beta * state.meanC - gamma * state.meanS; final amplitude = math.sqrt(beta * beta + gamma * gamma); - final acro = math.atan2(-gamma, beta); // peak phase - // Convert acrophase to clock-hours within the period (0..period). - var phaseHours = -acro / w; // time of peak relative to origin + final acro = math.atan2(-gamma, beta); + final w = 2 * math.pi / periodHours; + var phaseHours = -acro / w; phaseHours = phaseHours % periodHours; if (phaseHours < 0) phaseHours += periodHours; - // R². - final yMean = sY / n; - var ssTot = 0.0, ssRes = 0.0; - for (var i = 0; i < n; i++) { - final fit = mesor + beta * c[i] + gamma * s[i]; - ssTot += (y[i] - yMean) * (y[i] - yMean); - ssRes += (y[i] - fit) * (y[i] - fit); - } + final ssTot = state.mYY; + final ssRes = state.mYY - + 2 * beta * state.mCY - + 2 * gamma * state.mSY + + beta * beta * state.mCC + + gamma * gamma * state.mSS + + 2 * beta * gamma * state.mCS; final double r2 = ssTot == 0 ? 0.0 : (1 - ssRes / ssTot).clamp(0, 1); - // Adjusted for the 3 fitted parameters (M, β, γ). Confidence MUST come from - // the adjusted value: the raw R² of a 3-parameter fit is upward-biased - // (E[R²] = 2/(n−1) under the null), so a handful of noise points used to - // score confidence 0.95 at tier HIGH. final double r2Adj = (1 - (1 - r2) * (n - 1) / (n - 3)).clamp(0, 1); - // Zero-amplitude F-test (Nelson 1979): the acrophase is only defined once // the test rejects. A flat or rhythm-free series still yields some // atan2(−γ, β) from residual noise or float rounding. @@ -169,7 +189,6 @@ Metric cosinor( } final phaseDefined = ampP < cosinorAcrophaseAlpha; - final conf = r2Adj.clamp(0.1, 0.95); return Metric( value: CosinorFit( mesor: mesor, @@ -181,7 +200,7 @@ Metric cosinor( r2Adj: r2Adj, amplitudeP: ampP, ), - confidence: conf, + confidence: r2Adj.clamp(0.1, 0.95), tier: Tier.high, inputs_used: inputs, note: phaseDefined @@ -192,34 +211,12 @@ Metric cosinor( ); } -/// Solve a 3×3 linear system by Gaussian elimination with partial pivoting. -/// Returns null if singular. -List? _solve3(List> a, List b) { - final m = [ - [a[0][0], a[0][1], a[0][2], b[0]], - [a[1][0], a[1][1], a[1][2], b[1]], - [a[2][0], a[2][1], a[2][2], b[2]], - ]; - for (var col = 0; col < 3; col++) { - var piv = col; - for (var r = col + 1; r < 3; r++) { - if (m[r][col].abs() > m[piv][col].abs()) piv = r; - } - if (m[piv][col].abs() < 1e-12) return null; - final tmp = m[col]; - m[col] = m[piv]; - m[piv] = tmp; - final d = m[col][col]; - for (var j = col; j < 4; j++) { - m[col][j] /= d; - } - for (var r = 0; r < 3; r++) { - if (r == col) continue; - final f = m[r][col]; - for (var j = col; j < 4; j++) { - m[r][j] -= f * m[col][j]; - } - } - } - return [m[0][3], m[1][3], m[2][3]]; + +/// Upper tail of F(2, d): P(F > f) = (1 + 2f/d)^(−d/2), exact for 2 numerator +/// degrees of freedom. +double _f2Survival(double f, int d) { + if (f.isNaN) return 1.0; + if (f.isInfinite) return 0.0; + if (f <= 0) return 1.0; + return math.pow(1 + 2 * f / d, -d / 2).toDouble(); } diff --git a/test/output_changes/cosinor_test.dart b/test/output_changes/cosinor_test.dart new file mode 100644 index 00000000..6d0e6d45 --- /dev/null +++ b/test/output_changes/cosinor_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Refuse unidentifiable cosinor fits: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['cosinor', 'near_singular']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From 10c6b1c6d28f9fb9c74d346f98f589cd41afbac7 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:30:33 -0700 Subject: [PATCH 05/14] Retain small TRIMP contributions in the total Morton, Fitz-Clarke & Banister (1990); Neumaier (1974) compensated summation. Preserve weighting and anchors; pin changed last digits. --- ALGORITHMS.md | 6 ++++++ CHANGELOG.md | 6 ++++++ lib/src/onehz/clinical/load_trimp.dart | 13 ++++++++++--- test/output_changes/trimp_test.dart | 16 ++++++++++++++++ 4 files changed, 38 insertions(+), 3 deletions(-) create mode 100644 test/output_changes/trimp_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index b397fb3d..0f555bff 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -182,3 +182,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Refuse unidentifiable cosinor fits | Nelson et al. 1979, Methods for cosinor-rhythmometry; Welford 1962 centered co-moments; Use centered least squares and a relative determinant gate of 1e-12 of the diagonal product. This is a numerical policy. Nearly collinear phase columns abstain, and ordinary fit parameters can move in their last digits. Upstream zero-amplitude significance gating is retained. | + + + +| Batch change | Method and limits | +|---|---| +| Retain small TRIMP contributions in the total | Morton, Fitz-Clarke & Banister 1990; Neumaier 1974 compensated summation; Compensate the per-minute sum. The reserve and sex-specific weighting are unchanged; the synthetic total moves in its last digits. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 46723fe2..6bf6b8de 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,3 +28,9 @@ Use centered least squares and a relative determinant gate of 1e-12 of the diago Method: Nelson et al. 1979, Methods for cosinor-rhythmometry; Welford 1962 centered co-moments. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Retain small TRIMP contributions in the total + +Compensate the per-minute sum. The reserve and sex-specific weighting are unchanged; the synthetic total moves in its last digits. + +Method: Morton, Fitz-Clarke & Banister 1990; Neumaier 1974 compensated summation. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/clinical/load_trimp.dart b/lib/src/onehz/clinical/load_trimp.dart index 3c89ca4f..42ce56df 100644 --- a/lib/src/onehz/clinical/load_trimp.dart +++ b/lib/src/onehz/clinical/load_trimp.dart @@ -63,7 +63,8 @@ Metric banisterTrimp( ); } final reserve = maxHr - restingHr; - var trimp = 0.0; + // Neumaier (1974), Rundungsfehleranalyse einiger Verfahren zur Summation. + var trimp = 0.0, correction = 0.0; for (final hr in hrPerMin) { if (!hr.isFinite || hr <= 0) continue; // off-skin + non-finite guard var hrr = (hr - restingHr) / reserve; @@ -71,10 +72,16 @@ Metric banisterTrimp( if (hrr > 1) hrr = 1; // ONE Banister implementation for the whole package: the sex-specific // weighting factor y lives in [StrainScorer.banisterY]. 1 minute each. - trimp += 1.0 * hrr * StrainScorer.banisterY(hrr, female: sex == Sex.female); + final value = hrr * StrainScorer.banisterY(hrr, female: sex == Sex.female); + final next = trimp + value; + if (next.isFinite) { + correction += trimp.abs() >= value.abs() + ? (trimp - next) + value : (value - next) + trimp; + } + trimp = next; } return Metric( - value: trimp, + value: trimp + correction, confidence: 0.6, tier: Tier.estimate, inputs_used: inputs, diff --git a/test/output_changes/trimp_test.dart b/test/output_changes/trimp_test.dart new file mode 100644 index 00000000..afaf57d5 --- /dev/null +++ b/test/output_changes/trimp_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Retain small TRIMP contributions in the total: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['trimp']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From e66cbe16e5c978c1eee2e298f65f715c8956ce36 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:30:39 -0700 Subject: [PATCH 06/14] Evaluate training-load decay with affine arithmetic Banister et al. (1975); Morton et al. (1990); Neumaier (1974). Preserve EWMA formulas and bounded-prefix priming while matching the arithmetic of the later incremental path. --- ALGORITHMS.md | 6 +++++ CHANGELOG.md | 6 +++++ lib/src/onehz/clinical/load_trimp.dart | 32 +++++++++++++++++++++----- test/output_changes/ctl_test.dart | 16 +++++++++++++ 4 files changed, 54 insertions(+), 6 deletions(-) create mode 100644 test/output_changes/ctl_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 0f555bff..6f3ce74b 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -188,3 +188,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Retain small TRIMP contributions in the total | Morton, Fitz-Clarke & Banister 1990; Neumaier 1974 compensated summation; Compensate the per-minute sum. The reserve and sex-specific weighting are unchanged; the synthetic total moves in its last digits. | + + + +| Batch change | Method and limits | +|---|---| +| Evaluate training-load decay with affine arithmetic | Banister et al. 1975; Morton et al. 1990; Neumaier 1974 prefix summation; Evaluate the same EWMA recurrence as scale times seed plus offset. This matches the arithmetic used by the later incremental implementation. Last digits change; 42/7-day defaults and bounded-prefix priming remain implementation choices. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 6bf6b8de..028f5338 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -34,3 +34,9 @@ Compensate the per-minute sum. The reserve and sex-specific weighting are unchan Method: Morton, Fitz-Clarke & Banister 1990; Neumaier 1974 compensated summation. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Evaluate training-load decay with affine arithmetic + +Evaluate the same EWMA recurrence as scale times seed plus offset. This matches the arithmetic used by the later incremental implementation. Last digits change; 42/7-day defaults and bounded-prefix priming remain implementation choices. + +Method: Banister et al. 1975; Morton et al. 1990; Neumaier 1974 prefix summation. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/clinical/load_trimp.dart b/lib/src/onehz/clinical/load_trimp.dart index 42ce56df..89136ef5 100644 --- a/lib/src/onehz/clinical/load_trimp.dart +++ b/lib/src/onehz/clinical/load_trimp.dart @@ -863,13 +863,33 @@ Metric ctlAtlTsb(List dailyTrimp, final lc = 1 - math.exp(-1 / ctlDays); final la = 1 - math.exp(-1 / atlDays); final prime = math.min(math.max(primeDays, 1), dailyTrimp.length); - final seed = mean(dailyTrimp.sublist(0, prime))!; - var ctl = seed; - var atl = seed; - for (var i = prime; i < dailyTrimp.length; i++) { - ctl = ctl + lc * (dailyTrimp[i] - ctl); - atl = atl + la * (dailyTrimp[i] - atl); + // Neumaier (1974) compensation for the bounded-prefix mean. + var sum = 0.0, correction = 0.0; + for (final value in dailyTrimp.take(prime)) { + if (value.isNaN) continue; + final next = sum + value; + if (next.isFinite) { + correction += sum.abs() >= value.abs() + ? (sum - next) + value : (value - next) + sum; + } + sum = next; + } + final seed = (sum + correction) / prime; + // Banister et al. (1975) recurrence evaluated as an affine map. + // f(x) = scale*x + offset, with keep = 1-lambda. + double evaluate(double lambda) { + var scale = 1.0, offset = 0.0; + final keep = 1 - lambda; + for (var i = prime; i < dailyTrimp.length; i++) { + final value = dailyTrimp[i]; + if (value.isNaN) continue; + scale *= keep; + offset = keep * offset + lambda * value; + } + return scale * seed + offset; } + final ctl = evaluate(lc); + final atl = evaluate(la); final conf = (dailyTrimp.length / 42.0).clamp(0.3, 0.85); return Metric( value: LoadState(ctl, atl, ctl - atl), diff --git a/test/output_changes/ctl_test.dart b/test/output_changes/ctl_test.dart new file mode 100644 index 00000000..0c1d61ad --- /dev/null +++ b/test/output_changes/ctl_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Evaluate training-load decay with affine arithmetic: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['ctl']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From 58e4f1480df1133c038728b6e9e53427f3235c8e Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:30:44 -0700 Subject: [PATCH 07/14] Count repeated zone intervals before summing Karvonen (1957) supports reserve zones only. Duration credit is an internal ESTIMATE rule with low confidence; preserve cadence validation, median capping and tail credit. --- ALGORITHMS.md | 6 ++++ CHANGELOG.md | 6 ++++ lib/src/onehz/workout/hr_zones.dart | 46 +++++++++++++---------------- test/output_changes/zones_test.dart | 16 ++++++++++ 4 files changed, 49 insertions(+), 25 deletions(-) create mode 100644 test/output_changes/zones_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 6f3ce74b..ec1fc7be 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -194,3 +194,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Evaluate training-load decay with affine arithmetic | Banister et al. 1975; Morton et al. 1990; Neumaier 1974 prefix summation; Evaluate the same EWMA recurrence as scale times seed plus offset. This matches the arithmetic used by the later incremental implementation. Last digits change; 42/7-day defaults and bounded-prefix priming remain implementation choices. | + + + +| Batch change | Method and limits | +|---|---| +| Count repeated zone intervals before summing | Karvonen 1957 supports reserve zones only; interval credit is an internal ESTIMATE rule; Group equal duration credits by zone. Cadence validation, median capping and tail credit are retained. Grouped multiplication changes the last digits of zone seconds. | diff --git a/CHANGELOG.md b/CHANGELOG.md index 028f5338..b6b5230c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,3 +40,9 @@ Evaluate the same EWMA recurrence as scale times seed plus offset. This matches Method: Banister et al. 1975; Morton et al. 1990; Neumaier 1974 prefix summation. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Count repeated zone intervals before summing + +Group equal duration credits by zone. Cadence validation, median capping and tail credit are retained. Grouped multiplication changes the last digits of zone seconds. + +Method: Karvonen 1957 supports reserve zones only; interval credit is an internal ESTIMATE rule. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/workout/hr_zones.dart b/lib/src/onehz/workout/hr_zones.dart index 16c66e37..d2efec62 100644 --- a/lib/src/onehz/workout/hr_zones.dart +++ b/lib/src/onehz/workout/hr_zones.dart @@ -182,33 +182,29 @@ class HeartRateZones { HeartRateZoneSet zoneSet, ) { final sorted = [...hr]..sort((a, b) => a.tsMs.compareTo(b.tsMs)); - final zoneSeconds = List.filled(5, 0); - var below = 0.0; - - final tailSeconds = - sampleCadenceSeconds([for (final s in sorted) s.tsMs / 1000.0]); - if (tailSeconds == null) return null; - for (var i = 0; i < sorted.length; i++) { - final sample = sorted[i]; - if (!sample.valid) continue; - final durSeconds = i < sorted.length - 1 - ? _boundedGapSeconds(sorted[i + 1].tsMs - sample.tsMs, tailSeconds) - : tailSeconds; - final zone = zoneSet.zoneNumber(sample.hr); - if (zone >= 1) { - zoneSeconds[zone - 1] += durSeconds; - } else { - below += durSeconds; + // Karvonen (1957) supports reserve zones only. Duration credit is an + // internal ESTIMATE rule. Count equal gaps before summing their credit. + final cadence = sampleCadenceSeconds([for (final s in sorted) s.tsMs / 1000.0]); + if (cadence == null) return null; + final credits = List.generate(6, (_) => {}); + for (var i = 1; i < sorted.length; i++) { + final previous = sorted[i - 1]; + if (!previous.valid) continue; + final gap = (sorted[i].tsMs - previous.tsMs) / 1000.0; + final histogram = credits[zoneSet.zoneNumber(previous.hr)]; + histogram[gap] = (histogram[gap] ?? 0) + 1; + } + final seconds = List.filled(6, 0); + for (var zone = 0; zone < credits.length; zone++) { + for (final entry in credits[zone].entries) { + final gap = entry.key; + seconds[zone] += entry.value * (gap > 0 ? math.min(gap, cadence) : cadence); } } - return TimeInHeartRateZone(seconds: zoneSeconds, belowZone1: below); - } - - static double _boundedGapSeconds(double gapMs, double fallbackSeconds) { - final gapSeconds = gapMs / 1000.0; - return gapSeconds > 0 - ? math.min(gapSeconds, fallbackSeconds) - : fallbackSeconds; + if (sorted.isNotEmpty && sorted.last.valid) { + seconds[zoneSet.zoneNumber(sorted.last.hr)] += cadence; + } + return TimeInHeartRateZone(seconds: seconds.sublist(1), belowZone1: seconds[0]); } } diff --git a/test/output_changes/zones_test.dart b/test/output_changes/zones_test.dart new file mode 100644 index 00000000..c9b772dc --- /dev/null +++ b/test/output_changes/zones_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Count repeated zone intervals before summing: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['zones']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From c8e325a958c5b7455a6e48011165cf4a3af47f41 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:30:57 -0700 Subject: [PATCH 08/14] Compute circadian variance with centered moments Witting (1990); van Someren (1999); Welford (1962). Preserve batch profile sums and first-winner ties; centered variance changes last digits in IS and IV. --- ALGORITHMS.md | 6 ++++++ CHANGELOG.md | 6 ++++++ lib/src/onehz/sleep/circadian_np.dart | 17 +++++++++++++---- test/output_changes/circadian_test.dart | 16 ++++++++++++++++ 4 files changed, 41 insertions(+), 4 deletions(-) create mode 100644 test/output_changes/circadian_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index ec1fc7be..06321405 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -200,3 +200,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Count repeated zone intervals before summing | Karvonen 1957 supports reserve zones only; interval credit is an internal ESTIMATE rule; Group equal duration credits by zone. Cadence validation, median capping and tail credit are retained. Grouped multiplication changes the last digits of zone seconds. | + + + +| Batch change | Method and limits | +|---|---| +| Compute circadian variance with centered moments | Witting 1990; van Someren 1999; Welford 1962 centered variance; IS and IV may move in their last digits. Batch circular profile sums and first-winner tie selection are unchanged on the rebased head. Chunked merges can reorder profile sums and move near-tie L5/M10 start epochs; that stream-specific behavior belongs with the aggregators. | diff --git a/CHANGELOG.md b/CHANGELOG.md index b6b5230c..a81e3617 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -46,3 +46,9 @@ Group equal duration credits by zone. Cadence validation, median capping and tai Method: Karvonen 1957 supports reserve zones only; interval credit is an internal ESTIMATE rule. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Compute circadian variance with centered moments + +IS and IV may move in their last digits. Batch circular profile sums and first-winner tie selection are unchanged on the rebased head. Chunked merges can reorder profile sums and move near-tie L5/M10 start epochs; that stream-specific behavior belongs with the aggregators. + +Method: Witting 1990; van Someren 1999; Welford 1962 centered variance. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/sleep/circadian_np.dart b/lib/src/onehz/sleep/circadian_np.dart index c80c4789..d5365969 100644 --- a/lib/src/onehz/sleep/circadian_np.dart +++ b/lib/src/onehz/sleep/circadian_np.dart @@ -74,10 +74,19 @@ Metric circadianNonparametric( final lEp = l5Epochs ?? (epochsPerDay * 5 ~/ 24); final grand = mean(x)!; - var ssTot = 0.0; - for (final v in x) { - final d = v - grand; - ssTot += d * d; + // Welford (1962) centered variance; Witting (1990), van Someren (1999) + // indices retain their formulas. Circular profile sums keep input order. + var ssTot = 0.0, center = 0.0, count = 0; + for (final value in x) { + count++; + if (!value.isFinite) { + ssTot = double.nan; + center = x.take(count).fold(0, (a, b) => a + b) / count; + continue; + } + final delta = value - center; + center += delta / count; + ssTot += delta * (value - center); } if (ssTot == 0) { return const Metric.absent( diff --git a/test/output_changes/circadian_test.dart b/test/output_changes/circadian_test.dart new file mode 100644 index 00000000..480799e7 --- /dev/null +++ b/test/output_changes/circadian_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Compute circadian variance with centered moments: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['circadian']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From 580e16733ea6cd0dec81066417d2e89972bdbb54 Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:31:38 -0700 Subject: [PATCH 09/14] Compute HRV moments without subtracting a large mean Task Force (1996); Welford (1962). Port centered NN, lag and bin moments into batch code and retain shared spectral jitter policy. Rejection and breathing-line exemptions remain internal ESTIMATE rules with low confidence. --- ALGORITHMS.md | 6 + CHANGELOG.md | 6 + lib/src/onehz/clinical/hrv_time.dart | 370 ++++++++----------------- lib/src/onehz/clinical/rmssd_gate.dart | 319 +++++++++++++++++++++ test/output_changes/hrv_test.dart | 16 ++ 5 files changed, 470 insertions(+), 247 deletions(-) create mode 100644 lib/src/onehz/clinical/rmssd_gate.dart create mode 100644 test/output_changes/hrv_test.dart diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 06321405..efbd80f2 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -206,3 +206,9 @@ Malik et al. | Batch change | Method and limits | |---|---| | Compute circadian variance with centered moments | Witting 1990; van Someren 1999; Welford 1962 centered variance; IS and IV may move in their last digits. Batch circular profile sums and first-winner tie selection are unchanged on the rebased head. Chunked merges can reorder profile sums and move near-tie L5/M10 start epochs; that stream-specific behavior belongs with the aggregators. | + + + +| Batch change | Method and limits | +|---|---| +| Compute HRV moments without subtracting a large mean | Task Force 1996; Welford 1962 centered sums and products; Use origin-shifted moments for NN levels and five-minute bins, and centered lag co-moments for ACF. Last digits and decisions within rounding of jitter gates can change. Existing respiratory-line exemptions, coverage gates and absent-value rules remain. | diff --git a/CHANGELOG.md b/CHANGELOG.md index a81e3617..79caddc1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -52,3 +52,9 @@ IS and IV may move in their last digits. Batch circular profile sums and first-w Method: Witting 1990; van Someren 1999; Welford 1962 centered variance. Reproduce with `dart run tool/reproduce_outputs.dart`. +### Compute HRV moments without subtracting a large mean + +Use origin-shifted moments for NN levels and five-minute bins, and centered lag co-moments for ACF. Last digits and decisions within rounding of jitter gates can change. Existing respiratory-line exemptions, coverage gates and absent-value rules remain. + +Method: Task Force 1996; Welford 1962 centered sums and products. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index cd1a5658..96c4c09d 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -17,6 +17,7 @@ import 'dart:math' as math; import '../types.dart'; import '../util.dart'; +import 'rmssd_gate.dart'; /// Lag-1 ACF floor for the NN successive-difference series, below which RMSSD /// and pNN50 are REFUSED (SDNN / SDANN survive it and are the honest lead). @@ -121,82 +122,13 @@ const double kNnDiffNoiseShareCeiling = 0.7; /// clear of the floor, or diffs on a beat-time grid coarse against their size /// (see [_onCoarseLattice]). double? nnDiffNoiseShare(List> diffRuns, {int minSegments = 20}) { - const n = 64, kLo = 10, kHi = n ~/ 2; // kLo/n ≈ 0.15 cycles/beat - final w = [ - for (var i = 0; i < n; i++) 0.5 - 0.5 * math.cos(2 * math.pi * i / (n - 1)) - ]; - var w2 = 0.0; - for (final x in w) { - w2 += x * x; - } - final cosT = [ - for (var k = kLo; k <= kHi; k++) - [for (var i = 0; i < n; i++) math.cos(2 * math.pi * k * i / n)] - ]; - final sinT = [ - for (var k = kLo; k <= kHi; k++) - [for (var i = 0; i < n; i++) math.sin(2 * math.pi * k * i / n)] - ]; - final psd = List.filled(kHi + 1, 0.0); - var segs = 0, nd = 0; - var ssd = 0.0, wsd = 0.0, wsum = 0.0, pSum = 0.0, p2Sum = 0.0; + final sums = NoiseShareSums(); for (final r in diffRuns) { - for (final d in r) { - ssd += d * d; - nd++; - } - // Integrate the run back to RR levels (up to a constant the mean removes). - final x = List.filled(r.length + 1, 0.0); - for (var i = 0; i < r.length; i++) { - x[i + 1] = x[i] + r[i]; - } - final cov = List.filled(x.length, 0.0); // Σ w² per sample - for (var s = 0; s + n <= x.length; s += n ~/ 2) { - var m = 0.0; - for (var i = 0; i < n; i++) { - m += x[s + i]; - cov[s + i] += w[i] * w[i]; - } - m /= n; - final v = [for (var i = 0; i < n; i++) (x[s + i] - m) * w[i]]; - var p = 0.0; - for (var k = kLo; k <= kHi; k++) { - final c = cosT[k - kLo], sn = sinT[k - kLo]; - var re = 0.0, im = 0.0; - for (var i = 0; i < n; i++) { - re += v[i] * c[i]; - im -= v[i] * sn[i]; - } - psd[k] += re * re + im * im; - p += re * re + im * im; - } - pSum += p; - p2Sum += p * p; - segs++; - } - for (var i = 0; i < r.length; i++) { - final c = (cov[i] + cov[i + 1]) / 2; - wsd += c * r[i] * r[i]; - wsum += c; - } + sums.addRun(r); } - if (segs < minSegments || ssd == 0 || wsd == 0) return null; - if (_onCoarseLattice(diffRuns, ssd / nd)) return null; - // One loud stretch dominates the average: count segments by power, not number. - if (pSum * pSum / p2Sum < 0.75 * minSegments) return null; - final band = psd.sublist(kLo); - var peak = 0; - for (var k = 1; k < band.length; k++) { - if (band[k] > band[peak]) peak = k; - } - if (peak + kLo >= kHi - 1) return null; // ~0.477–0.5 band, see above - if (band.last > 0.2 * band[peak]) return null; // hump on Nyquist, see above - final med = median(band)!; - if (band[peak] < 4.5 * med) return null; // no line, no exemption - final floor = med / segs / w2; - final vetted = wsd / wsum, all = ssd / nd; - // 1 − (structured share of RMSSD²); equals 2σ²/mean(d²) when stationary. - return 1 - (vetted - 2 * floor) / math.max(vetted, all); + return noiseShareFrom(sums, + lattice: () => _onCoarseLattice(diffRuns, sums.ssd / sums.nd), + minSegments: minSegments); } /// True when every successive difference is a whole multiple of one step q and @@ -208,22 +140,8 @@ double? nnDiffNoiseShare(List> diffRuns, {int minSegments = 20}) { /// sets RMSSD. // ponytail: exact lattice only; a grid re-rounded to whole ms (7.8 ms → 7/8) // breaks the common step and is not caught. -bool _onCoarseLattice(List> diffRuns, double msd) { - var q = double.infinity; - for (final r in diffRuns) { - for (final d in r) { - if (d.abs() > 1e-6 && d.abs() < q) q = d.abs(); - } - } - if (q.isInfinite || msd >= 2 * q * q) return false; - for (final r in diffRuns) { - for (final d in r) { - final k = d.abs() / q; - if ((k - k.roundToDouble()).abs() > 1e-3) return false; - } - } - return true; -} +bool _onCoarseLattice(List> diffRuns, double msd) => + coarseLattice([for (final r in diffRuns) for (final d in r) d.abs()], msd); /// The one RMSSD jitter verdict every path shares: ACF1 below the floor AND /// the spectrum does not show a respiratory line on a low white floor. @@ -256,48 +174,6 @@ List _clearedWindows(double? nightAcf1, List rmssds, ]; } -final _hann128 = [ - for (var i = 0; i < 128; i++) 0.5 - 0.5 * math.cos(2 * math.pi * i / 127) -]; -final _hann128Sq = _hann128.fold(0.0, (a, w) => a + w * w); - -/// Welch power (128-beat Hann, 50 % overlap, averaged over segments) of one -/// window's RR, rebuilt from its difference runs, at cycles-per-beat -/// frequencies [f]. Null when no run holds a full segment. -List? _windowPsd(List> diffRuns, List f) { - const n = 128; - final p = List.filled(f.length, 0.0); - var segs = 0; - for (final r in diffRuns) { - final x = List.filled(r.length + 1, 0.0); - for (var i = 0; i < r.length; i++) { - x[i + 1] = x[i] + r[i]; - } - for (var s = 0; s + n <= x.length; s += n ~/ 2) { - var m = 0.0; - for (var i = 0; i < n; i++) { - m += x[s + i]; - } - m /= n; - for (var j = 0; j < f.length; j++) { - final c = math.cos(2 * math.pi * f[j]), sn = math.sin(2 * math.pi * f[j]); - var cr = 1.0, ci = 0.0, re = 0.0, im = 0.0; - for (var i = 0; i < n; i++) { - final v = (x[s + i] - m) * _hann128[i]; - re += v * cr; - im += v * ci; - final t = cr * c - ci * sn; - ci = cr * sn + ci * c; - cr = t; - } - p[j] += re * re + im * im; - } - segs++; - } - } - return segs == 0 ? null : [for (final v in p) v / segs]; -} - /// Last resort for a night the jitter gate refused: the windows whose /// breathing line holds still in Hz while the heart rate moves. /// @@ -324,115 +200,29 @@ List? _steadyBreathingWindows( nd++; } } - if (nd == 0 || _onCoarseLattice(all, ssd / nd)) return null; - final ref = median([for (final m in meanRrMs) if (m > 0) m]); - if (ref == null) return null; - // 0.15–0.5 cycles/beat at the night's median beat. - final fc = [for (var j = 38; j <= 128; j++) j / 256]; - final byBeat = List.filled(fc.length, 0.0); - final byHz = List.filled(fc.length, 0.0); - final cover = List.filled(fc.length, 0); - final peakHz = {}; // window -> its own peak, cycles per ms - final psd = >{}; // window -> its cycles-per-beat PSD - for (var w = 0; w < winRuns.length; w++) { - if (meanRrMs[w] <= 0) continue; - final pb = _windowPsd(winRuns[w], fc); - if (pb == null) continue; - final norm = median(pb)!; - if (norm <= 0) continue; - final k = meanRrMs[w] / ref; - final fh = [for (final f in fc) f * k]; - final ph = _windowPsd(winRuns[w], fh)!; - var best = -1; - for (var j = 0; j < fc.length; j++) { - byBeat[j] += pb[j] / norm; - if (fh[j] > 0.5) continue; - byHz[j] += ph[j] / norm; - cover[j]++; - if (best < 0 || ph[j] > ph[best]) best = j; - } - if (best >= 0) peakHz[w] = fc[best] / ref; - psd[w] = pb; - } - final nw = peakHz.length; - if (nw < 12) return null; - final bins = [for (var j = 0; j < fc.length; j++) if (cover[j] == nw) j]; - if (bins.length < 20) return null; - final hz = [for (final j in bins) byHz[j]]; - final bt = [for (final j in bins) byBeat[j]]; - var pk = 0; - for (var i = 1; i < hz.length; i++) { - if (hz[i] > hz[pk]) pk = i; - } - if (pk == 0 || pk == hz.length - 1) return null; - final f0 = fc[bins[pk]] / ref; - final brpm = f0 * 60000; - if (brpm < 8 || brpm > 30) return null; - // Same bins, same normalisation: the line must stand taller aligned in Hz - // than the beat pooling stands anywhere within ±12 bins of it. - var near = 0.0; - for (var i = math.max(0, pk - 12); i <= math.min(bt.length - 1, pk + 12); i++) { - near = math.max(near, bt[i]); - } - // The pooled noise ripple shrinks as 1/√windows, so the bar does too. - final promHz = hz[pk] / median(hz)!; - if (promHz < math.max(2.5, 1.5 + 7 / math.sqrt(nw))) return null; - if (promHz < 1.1 * near / median(bt)!) return null; - // Alternation lives at Nyquist on the beat axis, outside these bins. - if (byBeat.last / median(bt)! >= 0.5 * promHz) return null; - final keep = [ - for (final e in peakHz.entries) - if (math.log(e.value / f0).abs() <= 0.15) e.key - ]; - if (keep.length < 6) return null; - // The line proves breathing is there, not that it carries RMSSD: beat-time - // jitter loud enough to fail the gate would still be most of the number. - // Same ceiling as [nnDiffNoiseShare], over the kept windows. Jitter σ² has - // an RR spectrum σ²·(2 − 2cos ω), flat once divided by that shape, and its - // differences carry 6σ²; physiology only adds power on top. Breathing - // wanders across the night and its line spreads with it, so each window - // drops the bins of every rate the kept windows peaked at (plus a Hann - // main lobe), mapped onto its own beats. σ² is the floor of what is left, - // pooled over windows and smoothed over four resolution cells. - var lo = double.infinity, hi = 0.0; - for (final w in keep) { - lo = math.min(lo, peakHz[w]!); - hi = math.max(hi, peakHz[w]!); - } - var sq = 0.0, nAll = 0; - final acc = List.filled(fc.length, 0.0); - final wt = List.filled(fc.length, 0.0); - for (final w in keep) { - var n = 0; - for (final r in winRuns[w]) { - for (final d in r) { - sq += d * d; - n++; + return steadyBreathingFrom( + windows: winRuns.length, + meanRrMs: meanRrMs, + psdAt: (w, f) => windowPsd(winRuns[w], f), + diffs: (w) { + var n = 0; + for (final r in winRuns[w]) { + n += r.length; } - } - nAll += n; - final a = lo * meanRrMs[w] - 3 / 128, b = hi * meanRrMs[w] + 3 / 128; - for (var j = 0; j < fc.length; j++) { - if (fc[j] > a && fc[j] < b) continue; - acc[j] += n * psd[w]![j] / (2 - 2 * math.cos(2 * math.pi * fc[j])); - wt[j] += n; - } - } - var floor = double.infinity; - for (var j = 8; j < fc.length - 8; j++) { - var s = 0.0, m = 0; - for (var i = j - 8; i <= j + 8; i++) { - if (wt[i] < nAll / 2) break; - s += acc[i] / wt[i]; - m++; - } - if (m == 17) floor = math.min(floor, s / 17); - } - if (floor == double.infinity) return null; - // The lowest stretch of a noisy curve reads below its mean; 1.1 puts the - // floor back at these pooled sizes, so jitter alone is not undercounted. - final noise = 1.1 * 6 * floor / _hann128Sq * nAll; - return noise < kNnDiffNoiseShareCeiling * sq ? keep : null; + return n; + }, + addSquares: (w, acc) { + for (final r in winRuns[w]) { + for (final d in r) { + acc += d * d; + } + } + return acc; + }, + allSsd: ssd, + allNd: nd, + lattice: (msd) => _onCoarseLattice(all, msd), + ); } /// The window RMSSDs a windowed headline publishes, or null for none: the @@ -681,7 +471,29 @@ Metric hrvTime( // are far less contaminated (jitter share 1–27 % against RMSSD's 11–100 % on // the audit corpus) — they keep publishing, which is what the header has // always advised. - final acf1 = nnDiffAcf1(runs); + // Welford (1962) centered sums/products, Task Force (1996) HRV. + // Accumulate difference moments and lagged co-moments in sample order. + var dn = 0, ln = 0; + var dm = 0.0, d2 = 0.0, ma = 0.0, mb = 0.0, co = 0.0; + for (final r in runs) { + for (var i = 0; i < r.length; i++) { + final d = r[i]; + dn++; + final delta = d - dm; + dm += delta / dn; + d2 += delta * (d - dm); + if (i > 0) { + ln++; + final da = r[i - 1] - ma; + final db = d - mb; + ma += da / ln; + mb += db / ln; + co += da * (d - mb); + } + } + } + final double? acf1 = dn >= 30 && d2 > 0 + ? (co + ln * (ma - dm) * (mb - dm)) / d2 : null; final jittery = _jitterRefused(acf1, runs); // More beat-time than elapsed: the stream itself is wrong, so neither the // differences nor any breathing line in them describe one heart. @@ -695,15 +507,16 @@ Metric hrvTime( rmssd = _breathingRmssd(nnMs, nnTimesMs); } final byBreathing = jittery && rmssd != null; - final sdnn = stddev(nnMs); + final sdnn = _centeredMoments(nnMs).sd; double? sdann, sdnnIndex; if (gapAware) { final seg = _fiveMinSegments(nnMs, nnTimesMs); if (seg.length >= 2) { - final means = [for (final s in seg) mean(s)!]; - sdann = stddev(means); - final sds = [for (final s in seg) stddev(s)].whereType().toList(); + final moments = [for (final s in seg) _centeredMoments(s)]; + final means = [for (final m in moments) m.mean!]; + sdann = _centeredMoments(means).sd; + final sds = [for (final m in moments) m.sd].whereType().toList(); sdnnIndex = sds.isEmpty ? null : mean(sds); } } @@ -900,7 +713,29 @@ Metric nocturnalRmssd( meanRr.add(mean([for (final i in seg) nnMs[i]])!); rmssds.add(math.sqrt(ssd / nd)); } - final acf1 = nnDiffAcf1(runs); + // Welford (1962) centered sums/products, Task Force (1996) HRV. + // Accumulate difference moments and lagged co-moments in sample order. + var dn = 0, ln = 0; + var dm = 0.0, d2 = 0.0, ma = 0.0, mb = 0.0, co = 0.0; + for (final r in runs) { + for (var i = 0; i < r.length; i++) { + final d = r[i]; + dn++; + final delta = d - dm; + dm += delta / dn; + d2 += delta * (d - dm); + if (i > 0) { + ln++; + final da = r[i - 1] - ma; + final db = d - mb; + ma += da / ln; + mb += db / ln; + co += da * (d - mb); + } + } + } + final double? acf1 = dn >= 30 && d2 > 0 + ? (co + ln * (ma - dm) * (mb - dm)) / d2 : null; final refused = _jitterRefused(acf1, runs); if (rmssds.isEmpty && !refused) { return const Metric.absent( @@ -1166,7 +1001,29 @@ Metric sleepSessionRmssdDetail( // THE HEADLINE nightly RMSSD (→ ln_rmssd → readiness). When the differences // are noise, the honest output is no headline, not a plausible one — the // readiness composite already treats a null HRV driver as absent. - final acf1 = nnDiffAcf1(runs); + // Welford (1962) centered sums/products, Task Force (1996) HRV. + // Accumulate difference moments and lagged co-moments in sample order. + var dn = 0, ln = 0; + var dm = 0.0, d2 = 0.0, ma = 0.0, mb = 0.0, co = 0.0; + for (final r in runs) { + for (var i = 0; i < r.length; i++) { + final d = r[i]; + dn++; + final delta = d - dm; + dm += delta / dn; + d2 += delta * (d - dm); + if (i > 0) { + ln++; + final da = r[i - 1] - ma; + final db = d - mb; + ma += da / ln; + mb += db / ln; + co += da * (d - mb); + } + } + } + final double? acf1 = dn >= 30 && d2 > 0 + ? (co + ln * (ma - dm) * (mb - dm)) / d2 : null; final refused = _jitterRefused(acf1, runs); if (rmssds.isEmpty && !refused) { return Metric.absent( @@ -1305,3 +1162,22 @@ List> _cleanWindowRuns(List rr, List ts) { if (run.isNotEmpty) runs.add(run); return runs; } + +// Welford, B. P. (1962), Technometrics 4(3):419-420. +// Origin-shifted batch recurrence; skips NaN just as #99's Moments does. +({double? mean, double? sd}) _centeredMoments(List values) { + var n = 0; + var origin = 0.0, offset = 0.0, m2 = 0.0; + for (final value in values) { + if (value.isNaN) continue; + if (!value.isFinite) throw ArgumentError.value(value, 'value', 'Must be finite'); + if (n == 0) { origin = value; n = 1; continue; } + final centered = value - origin; + final delta = centered - offset; + n++; + offset += delta / n; + m2 += delta * (centered - offset); + } + return (mean: n == 0 ? null : origin + offset, + sd: n < 2 ? null : math.sqrt(m2 / (n - 1))); +} diff --git a/lib/src/onehz/clinical/rmssd_gate.dart b/lib/src/onehz/clinical/rmssd_gate.dart new file mode 100644 index 00000000..6c8972df --- /dev/null +++ b/lib/src/onehz/clinical/rmssd_gate.dart @@ -0,0 +1,319 @@ +// Rejection and breathing-line exemption thresholds are an internal rule. +// Welch 1967 spectral estimation; Hirsch & Bishop 1981 provides the +// breathing/RSA basis. This is hrvTime()'s existing internal jitter policy, +// ESTIMATE with low confidence, not a published classifier. Shared formulas +// preserve the calculation; reordered sums can change decisions within rounding +// of a threshold. +// CLINICAL — the spectral half of the RMSSD jitter gate, as pieces. +// +// `hrvTime`, `nocturnalRmssd` and `sleepSessionWindowedRmssd` (hrv_time.dart) +// read the successive-difference runs themselves when the lag-1 ACF is under +// [kNnDiffAcf1Floor]: a respiratory line on a low white floor exempts the night +// ([nnDiffNoiseShare]), and a refused night may still publish from the 5-min +// windows on a breathing line steady in Hz. Both reduce to sums over Welch +// segments, so this file holds them as SUMS plus the decisions over those +// sums, written once: the batch fills the sums from the runs in hand, and +// `hrvTimeAggregator` fills the same sums from mergeable per-chunk state. +// +// Not exported from `onehz.dart`. + +import 'dart:math' as math; + +import '../util.dart'; +import 'hrv_time.dart' show kNnDiffNoiseShareCeiling; + +/// Welch segment length (beats) of the noise-share spectrum, the first and last +/// bins of its 0.15–0.5 cycles/beat band. +const int kNoiseSegment = 64, kNoiseLoBin = 10, kNoiseHiBin = kNoiseSegment ~/ 2; + +/// The segment's Hann window (symmetric, n − 1 in the denominator). +final List noiseHann = [ + for (var i = 0; i < kNoiseSegment; i++) + 0.5 - 0.5 * math.cos(2 * math.pi * i / (kNoiseSegment - 1)) +]; + +final double _noiseHannSq = () { + var w2 = 0.0; + for (final x in noiseHann) { + w2 += x * x; + } + return w2; +}(); + +final List> _cosT = [ + for (var k = kNoiseLoBin; k <= kNoiseHiBin; k++) + [for (var i = 0; i < kNoiseSegment; i++) math.cos(2 * math.pi * k * i / kNoiseSegment)] +]; +final List> _sinT = [ + for (var k = kNoiseLoBin; k <= kNoiseHiBin; k++) + [for (var i = 0; i < kNoiseSegment; i++) math.sin(2 * math.pi * k * i / kNoiseSegment)] +]; + +/// Adds the Welch segment `x[s .. s + 63]` (RR levels integrated from the +/// differences) to [psd] (bins [kNoiseLoBin]..[kNoiseHiBin]) and returns its +/// band power. +double addNoiseSegment(List x, int s, List psd) { + const n = kNoiseSegment; + var m = 0.0; + for (var i = 0; i < n; i++) { + m += x[s + i]; + } + m /= n; + final v = [for (var i = 0; i < n; i++) (x[s + i] - m) * noiseHann[i]]; + var p = 0.0; + for (var k = kNoiseLoBin; k <= kNoiseHiBin; k++) { + final c = _cosT[k - kNoiseLoBin], sn = _sinT[k - kNoiseLoBin]; + var re = 0.0, im = 0.0; + for (var i = 0; i < n; i++) { + re += v[i] * c[i]; + im -= v[i] * sn[i]; + } + psd[k] += re * re + im * im; + p += re * re + im * im; + } + return p; +} + +/// The sums [nnDiffNoiseShare] decides on. +class NoiseShareSums { + final List psd = List.filled(kNoiseHiBin + 1, 0.0); + int segs = 0, nd = 0; + double ssd = 0.0, wsd = 0.0, wsum = 0.0, pSum = 0.0, p2Sum = 0.0; + + /// One contiguous run of differences, in the batch's own order. + void addRun(List r) { + const n = kNoiseSegment; + for (final d in r) { + ssd += d * d; + nd++; + } + // Integrate the run back to RR levels (up to a constant the mean removes). + final x = List.filled(r.length + 1, 0.0); + for (var i = 0; i < r.length; i++) { + x[i + 1] = x[i] + r[i]; + } + final cov = List.filled(x.length, 0.0); // Σ w² per sample + for (var s = 0; s + n <= x.length; s += n ~/ 2) { + for (var i = 0; i < n; i++) { + cov[s + i] += noiseHann[i] * noiseHann[i]; + } + final p = addNoiseSegment(x, s, psd); + pSum += p; + p2Sum += p * p; + segs++; + } + for (var i = 0; i < r.length; i++) { + final c = (cov[i] + cov[i + 1]) / 2; + wsd += c * r[i] * r[i]; + wsum += c; + } + } +} + +/// The white-noise share from [s] (see [nnDiffNoiseShare] in hrv_time.dart for +/// what it measures and every refusal). [lattice] says whether the differences +/// sit on a coarse beat-time grid ([coarseLattice]); it is only asked when +/// the counts allow a verdict. +double? noiseShareFrom(NoiseShareSums s, + {required bool Function() lattice, int minSegments = 20}) { + final segs = s.segs; + if (segs < minSegments || s.ssd == 0 || s.wsd == 0) return null; + if (lattice()) return null; + // One loud stretch dominates the average: count segments by power, not number. + if (s.pSum * s.pSum / s.p2Sum < 0.75 * minSegments) return null; + final band = s.psd.sublist(kNoiseLoBin); + var peak = 0; + for (var k = 1; k < band.length; k++) { + if (band[k] > band[peak]) peak = k; + } + if (peak + kNoiseLoBin >= kNoiseHiBin - 1) return null; // ~0.477–0.5 band + if (band.last > 0.2 * band[peak]) return null; // hump on Nyquist + final med = median(band)!; + if (band[peak] < 4.5 * med) return null; // no line, no exemption + final floor = med / segs / _noiseHannSq; + final vetted = s.wsd / s.wsum, all = s.ssd / s.nd; + // 1 − (structured share of RMSSD²); equals 2σ²/mean(d²) when stationary. + return 1 - (vetted - 2 * floor) / math.max(vetted, all); +} + +/// True when every successive difference is a whole multiple of one step q +/// (the smallest |d| over 1e-6) and the mean squared difference [msd] is under +/// 2q². [magnitudes] are the |d| values (repeats change nothing), in any +/// order. See hrv_time.dart `_onCoarseLattice` for why. +bool coarseLattice(Iterable magnitudes, double msd) { + var q = double.infinity; + for (final a in magnitudes) { + if (a > 1e-6 && a < q) q = a; + } + if (q.isInfinite || msd >= 2 * q * q) return false; + for (final a in magnitudes) { + final k = a / q; + if ((k - k.roundToDouble()).abs() > 1e-3) return false; + } + return true; +} + +/// The breathing-line test's Hann window (128 beats). +final List hann128 = [ + for (var i = 0; i < 128; i++) 0.5 - 0.5 * math.cos(2 * math.pi * i / 127) +]; +final double _hann128Sq = hann128.fold(0.0, (a, w) => a + w * w); + +/// The cycles-per-beat axis the breathing-line test pools on: 0.15–0.5 +/// cycles/beat at the night's median beat. +final List breathingAxis = [for (var j = 38; j <= 128; j++) j / 256]; + +/// Welch power (128-beat Hann, 50 % overlap, averaged over segments) of one +/// window's RR, rebuilt from its difference runs, at cycles-per-beat +/// frequencies [f]. Null when no run holds a full segment. +List? windowPsd(List> diffRuns, List f) { + const n = 128; + final p = List.filled(f.length, 0.0); + var segs = 0; + for (final r in diffRuns) { + final x = List.filled(r.length + 1, 0.0); + for (var i = 0; i < r.length; i++) { + x[i + 1] = x[i] + r[i]; + } + for (var s = 0; s + n <= x.length; s += n ~/ 2) { + var m = 0.0; + for (var i = 0; i < n; i++) { + m += x[s + i]; + } + m /= n; + for (var j = 0; j < f.length; j++) { + final c = math.cos(2 * math.pi * f[j]), sn = math.sin(2 * math.pi * f[j]); + var cr = 1.0, ci = 0.0, re = 0.0, im = 0.0; + for (var i = 0; i < n; i++) { + final v = (x[s + i] - m) * hann128[i]; + re += v * cr; + im += v * ci; + final t = cr * c - ci * sn; + ci = cr * sn + ci * c; + cr = t; + } + p[j] += re * re + im * im; + } + segs++; + } + } + return segs == 0 ? null : [for (final v in p) v / segs]; +} + +/// The breathing-line test over [windows] 5-min windows, each read through +/// callbacks so a caller can hold its windows as runs or as sums: +/// [psdAt] is the window's [windowPsd] at the given frequencies (null without +/// a full segment), [diffs] its difference count, [addSquares] adds its +/// squared differences to a running total (in the caller's order), [meanRrMs] +/// its mean RR. [allSsd] / [allNd] are the squared differences and their count +/// over every window, and [lattice] is [coarseLattice] over every window's +/// differences at the mean square it is given. Returns the indices of the +/// windows whose own peak sits on the line, or null. See hrv_time.dart +/// `_steadyBreathingWindows` for the method. +List? steadyBreathingFrom({ + required int windows, + required List meanRrMs, + required List? Function(int w, List f) psdAt, + required int Function(int w) diffs, + required double Function(int w, double acc) addSquares, + required double allSsd, + required int allNd, + required bool Function(double msd) lattice, +}) { + if (allNd == 0 || lattice(allSsd / allNd)) return null; + final ref = median([for (final m in meanRrMs) if (m > 0) m]); + if (ref == null) return null; + final fc = breathingAxis; + final byBeat = List.filled(fc.length, 0.0); + final byHz = List.filled(fc.length, 0.0); + final cover = List.filled(fc.length, 0); + final peakHz = {}; // window -> its own peak, cycles per ms + final psd = >{}; // window -> its cycles-per-beat PSD + for (var w = 0; w < windows; w++) { + if (meanRrMs[w] <= 0) continue; + final pb = psdAt(w, fc); + if (pb == null) continue; + final norm = median(pb)!; + if (norm <= 0) continue; + final k = meanRrMs[w] / ref; + final fh = [for (final f in fc) f * k]; + final ph = psdAt(w, fh)!; + var best = -1; + for (var j = 0; j < fc.length; j++) { + byBeat[j] += pb[j] / norm; + if (fh[j] > 0.5) continue; + byHz[j] += ph[j] / norm; + cover[j]++; + if (best < 0 || ph[j] > ph[best]) best = j; + } + if (best >= 0) peakHz[w] = fc[best] / ref; + psd[w] = pb; + } + final nw = peakHz.length; + if (nw < 12) return null; + final bins = [for (var j = 0; j < fc.length; j++) if (cover[j] == nw) j]; + if (bins.length < 20) return null; + final hz = [for (final j in bins) byHz[j]]; + final bt = [for (final j in bins) byBeat[j]]; + var pk = 0; + for (var i = 1; i < hz.length; i++) { + if (hz[i] > hz[pk]) pk = i; + } + if (pk == 0 || pk == hz.length - 1) return null; + final f0 = fc[bins[pk]] / ref; + final brpm = f0 * 60000; + if (brpm < 8 || brpm > 30) return null; + // Same bins, same normalisation: the line must stand taller aligned in Hz + // than the beat pooling stands anywhere within ±12 bins of it. + var near = 0.0; + for (var i = math.max(0, pk - 12); i <= math.min(bt.length - 1, pk + 12); i++) { + near = math.max(near, bt[i]); + } + // The pooled noise ripple shrinks as 1/√windows, so the bar does too. + final promHz = hz[pk] / median(hz)!; + if (promHz < math.max(2.5, 1.5 + 7 / math.sqrt(nw))) return null; + if (promHz < 1.1 * near / median(bt)!) return null; + // Alternation lives at Nyquist on the beat axis, outside these bins. + if (byBeat.last / median(bt)! >= 0.5 * promHz) return null; + final keep = [ + for (final e in peakHz.entries) + if (math.log(e.value / f0).abs() <= 0.15) e.key + ]; + if (keep.length < 6) return null; + // Jitter that fails the gate would still be most of the number: same + // ceiling as the noise share, over the kept windows (see hrv_time.dart). + var lo = double.infinity, hi = 0.0; + for (final w in keep) { + lo = math.min(lo, peakHz[w]!); + hi = math.max(hi, peakHz[w]!); + } + var sq = 0.0, nAll = 0; + final acc = List.filled(fc.length, 0.0); + final wt = List.filled(fc.length, 0.0); + for (final w in keep) { + sq = addSquares(w, sq); + final n = diffs(w); + nAll += n; + final a = lo * meanRrMs[w] - 3 / 128, b = hi * meanRrMs[w] + 3 / 128; + for (var j = 0; j < fc.length; j++) { + if (fc[j] > a && fc[j] < b) continue; + acc[j] += n * psd[w]![j] / (2 - 2 * math.cos(2 * math.pi * fc[j])); + wt[j] += n; + } + } + var floor = double.infinity; + for (var j = 8; j < fc.length - 8; j++) { + var s = 0.0, m = 0; + for (var i = j - 8; i <= j + 8; i++) { + if (wt[i] < nAll / 2) break; + s += acc[i] / wt[i]; + m++; + } + if (m == 17) floor = math.min(floor, s / 17); + } + if (floor == double.infinity) return null; + // The lowest stretch of a noisy curve reads below its mean; 1.1 puts the + // floor back at these pooled sizes, so jitter alone is not undercounted. + final noise = 1.1 * 6 * floor / _hann128Sq * nAll; + return noise < kNnDiffNoiseShareCeiling * sq ? keep : null; +} diff --git a/test/output_changes/hrv_test.dart b/test/output_changes/hrv_test.dart new file mode 100644 index 00000000..05a5a9d9 --- /dev/null +++ b/test/output_changes/hrv_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Compute HRV moments without subtracting a large mean: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['hrv']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} From c17280c2139a9cbff6bf859546a4ead4e430619e Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Sun, 1 Nov 2026 04:34:30 -0700 Subject: [PATCH 10/14] Make the changed batch outputs independently reproducible Record synthetic old/new observations and freeze twelve batch-only cases before aggregator adoption. Pin unchanged batch timestamp-gap exclusion and circular near-tie selection as controls. Methods and internal-policy limits are documented with the results. --- docs/OUTPUT_REPRODUCTION.md | 81 ++++++++++++++++++++ test/fixtures/output_batch_parity.json | 1 + test/output_changes/batch_controls_test.dart | 42 ++++++++++ test/output_changes/batch_parity_test.dart | 21 +++++ tool/reproduce_outputs.dart | 23 +++++- 5 files changed, 166 insertions(+), 2 deletions(-) create mode 100644 docs/OUTPUT_REPRODUCTION.md create mode 100644 test/fixtures/output_batch_parity.json create mode 100644 test/output_changes/batch_controls_test.dart create mode 100644 test/output_changes/batch_parity_test.dart diff --git a/docs/OUTPUT_REPRODUCTION.md b/docs/OUTPUT_REPRODUCTION.md new file mode 100644 index 00000000..f8eb4ff9 --- /dev/null +++ b/docs/OUTPUT_REPRODUCTION.md @@ -0,0 +1,81 @@ +# Reproducing the #99 batch output changes + +This branch replaces the batch-output portion of #99. The aggregator branch stacks +on this branch. Edge must bump `kAlgoVersion` when adopting these changed outputs. + +Run `dart pub get --offline`, then `dart run tool/reproduce_outputs.dart`. +`--json` prints unrounded machine-readable results. The script calls batch APIs +only. Its input recipes are in `test/fixtures/output_changes.json`, a small +synthetic fixture. It contains no personal measurements. Formula parameters +expand to beat sequences in memory, including 20,000 nearly collinear cosinor +samples. The fixture's `old` and `new` values are recorded observations, not +values substituted into library results. + +Old values were measured on `3a7ca34`, new values on `f8e7d3c` with #98 merged. +The output branch reproduces the new values. Each output commit pins its own +unrounded result. Tests contain no NaN list-equality expectations. + +| Calculation | up/main `3a7ca34` | Split output branch | +|---|---|---| +| RR kept / dropped / corrected | `[199, 34, 0]` | `[200, 33, 0]` | +| Lomb-Scargle PSD at 0.137 Hz | `32627.260306094682` | `32627.260527635714` | +| Cosinor mesor / amplitude / R² / p | `[60.20591440999805, 1.9937577719142656, 0.9902626829318173, 2.900200339561918e-94]` | `[60.205914409998044, 1.9937577719142632, 0.9902626829318162, 2.9002003395774255e-94]` | +| Near-singular cosinor present | `True` | `False` | +| HRV RMSSD / SDNN / SDANN / SDNN index / ACF1 | `[3.615265234174367, 21.52663924532639, 1.067383333873884, 21.55214616158443, 0.9849552962327551]` | `[3.615265234174367, 21.52663924532641, 1.0673833338741983, 21.55214616158443, 0.984955296232754]` | +| HRV with infinite final timestamp present | `throws` | `False` | +| Banister TRIMP | `32.71197112260296` | `32.711971122602975` | +| CTL / ATL / TSB | `[98.49620732954641, 103.12182624417179, -4.625618914625377]` | `[98.49620732954642, 103.12182624417177, -4.625618914625349]` | +| Zone seconds, zones 1 through 5 | `[0.0, 120.01199999999982, 0.0, 0.0, 0.0]` | `[0.0, 120.01199999999977, 0.0, 0.0, 0.0]` | +| Circadian IS / IV / M10 / L5 / start epochs | `[1.0, 0.331267605633807, 3.505000000000011, 1.3600000000000079, 14, 0]` | `[1.0, 0.3312676056338072, 3.505000000000011, 1.3600000000000079, 14, 0]` | +| Gap-aware batch irregular SD1 / pNN70 | `[22.360679774997894, 83.33333333333333]` | `[22.360679774997894, 83.33333333333333]` | + +The synthetic RR run contains 100 normal 900 ms beats, 33 artifacts at 250 ms, +and 100 normal beats. The over-cap run becomes a gap. At 32 artifacts or fewer, +the capped and uncapped results agree. The cap is an internal ESTIMATE policy, +not a published physiological threshold. + +The cosinor phase columns in the near-singular case cover only about 20 +microhours around hour 3. Main accepts this numerically unidentifiable design; +the relative determinant gate abstains. The threshold is a numerical policy. +Nelson's significance test and nullable acrophase from #101 remain intact. + +The invalid-timestamp table uses HRV because main throws promptly. The focused +guard test also calls frequency HRV, RSA, CVHR, nocturnal RMSSD and session +RMSSD with NaN and both infinities. Running old frequency HRV/RSA with an infinite +last time is unsafe because its segment walk may never terminate. The script +does not call those old endless paths. + +Main already excludes timestamp gaps in the batch irregular-rhythm screen. +The unchanged row above is a control. #99 adds the corresponding stream policy, +previous-beat timestamp checkpointing and old-checkpoint refusal in the +aggregator branch. Diagnostics and their tests also stay there. + +Main and the rebased #99 head use the same circular profile sums and strict +first-winner comparisons for batch L5/M10. Their starts agree here and on the +near-tie controls. Centered variance changes IS/IV last digits. Chunked merges +can reorder profile sums, so near-tie starts may move in the aggregator branch. +The original #99 description's broad near-tie warning should be read in that +context. These checked fixtures and tolerances are not universal error bounds. + +To independently recompute the old column, create a detached worktree at +`3a7ca34`, copy `tool/reproduce_outputs.dart` and +`test/fixtures/output_changes.json` into the same paths there, run +`dart pub get --offline`, then run the tool with `--json`. The table printer +labels its recorded baseline separately from the currently computed result. + +Methods and limits are beside the code and in `ALGORITHMS.md`: Lipponen & +Tarvainen 2019 and Peltola 2012 for RR correction; Lomb 1976, Scargle 1982 and +Press & Rybicki 1989 for spectra; Nelson et al. 1979 for cosinor; Task Force +1996 and Welford 1962 for HRV moments; Morton, Fitz-Clarke & Banister 1990 and +Banister et al. 1975 for training load; Neumaier 1974 for compensated sums; +Witting 1990 and van Someren 1999 for circadian indices. Karvonen 1957 supports +reserve zones only. Gap credit, tail credit and the RR cap remain internal rules. + +`--corpus-json` expands twelve deterministic cases, varying RR run lengths, +spectral origins and frequencies, jitter frequencies, timestamp gaps and record +lengths. `test/fixtures/output_batch_parity.json` freezes the output branch's +unrounded public results before the aggregator conversion. Both branches run +`batch_parity_test.dart` against these fixed results. All twelve cases also +matched the rebased #99 head exactly in the split audit. This supplements the +existing law tests and frozen pre-#99 references; it does not establish parity +for every possible input. Both reproduction fixtures together are about 10 KB. diff --git a/test/fixtures/output_batch_parity.json b/test/fixtures/output_batch_parity.json new file mode 100644 index 00000000..9ddfdaff --- /dev/null +++ b/test/fixtures/output_batch_parity.json @@ -0,0 +1 @@ +{"source":"580e16733ea6cd0dec81066417d2e89972bdbb54","provenance":"Public batch outputs recorded on the batch-only output commit, before the framework refactor. Inputs are reproductionCases in tool/reproduce_outputs.dart.","cases":[{"rr":[201,0,1],"spectral":2191.9035532963126,"cosinor":[60.20705778058244,1.994548567920997,0.99021822497395,3.584442855968536e-94],"near_singular":false,"hrv":[3.615265234174367,21.52663924532641,1.0673833338741983,21.55214616158443,0.984955296232754],"nonfinite":false,"trimp":10.252064196950478,"ctl":[93.59303565371741,96.28901109956863,-2.695975445851218],"zones":[0.0,50.00499999999995,0.0,0.0,0.0],"circadian":[1.0,0.02086956521739154,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,17,0],"spectral":240.47476890807894,"cosinor":[60.20546901353539,1.9944760220371391,0.9906874461342637,3.6450131402730934e-95],"near_singular":false,"hrv":[14.544681533327767,21.511404623818784,0.2902542164251966,21.539302061762694,0.7620460195717869],"nonfinite":false,"trimp":13.023941743036463,"ctl":[94.97450463635269,101.20837984038745,-6.23387520403476],"zones":[0.0,61.00609999999995,0.0,0.0,0.0],"circadian":[1.0,0.2553191489361734,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,31,0],"spectral":4.814151167014313,"cosinor":[60.20670784840657,1.9932291959984847,0.9907102525961057,3.252232588931257e-95],"near_singular":false,"hrv":[39.550764408744726,21.5122083009717,0.40991536302833925,21.55029614999938,-0.7370544026973745],"nonfinite":false,"trimp":15.28358496194033,"ctl":[94.92528070209704,96.00888145442389,-1.083600752326845],"zones":[0.0,72.00719999999993,0.0,0.0,0.0],"circadian":[1.0,0.3312676056338072,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,32,0],"spectral":358.69230942776977,"cosinor":[60.206138700452115,1.9956478767563555,0.9905970400705738,5.7121975661113537e-95],"near_singular":false,"hrv":[42.09631660353519,21.499141419377167,0.29941934343651677,21.512944620556805,-0.9682323671330902],"nonfinite":false,"trimp":18.40304984916745,"ctl":[96.58887865497881,103.10914807456258,-6.520269419583769],"zones":[0.0,83.00829999999988,0.0,0.0,0.0],"circadian":[1.0,0.3688421052631628,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,33,0],"spectral":19.577975582737867,"cosinor":[60.206354166666685,1.99256198692952,0.9903700870164812,1.7315692772094542e-94],"near_singular":false,"hrv":[3.615265234174367,21.52663924532641,1.0673833338741983,21.55214616158443,0.984955296232754],"nonfinite":false,"trimp":20.354197211570078,"ctl":[95.91769269304092,96.02635036170024,-0.10865766865931903],"zones":[0.0,94.00939999999984,0.0,0.0,0.0],"circadian":[1.0,0.02086956521739154,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[200,33,0],"spectral":11606.502913565982,"cosinor":[60.2065027398695,1.994982877334022,0.9902341170988117,3.323420791543799e-94],"near_singular":false,"hrv":[14.544681533327767,21.511404623818784,0.2902542164251966,21.539302061762694,0.7620460195717869],"nonfinite":false,"trimp":23.84820212239254,"ctl":[97.91559020097185,105.11939669789578,-7.2038064969239315],"zones":[0.0,105.0104999999998,0.0,0.0,0.0],"circadian":[1.0,0.2553191489361734,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[201,0,1],"spectral":179.12281661844955,"cosinor":[60.20613433801807,1.993285411438357,0.9902158539353855,3.625067770914454e-94],"near_singular":false,"hrv":[39.54331943857053,21.508583812826043,0.27113546544007344,21.551316007749993,-0.7371891534838562],"nonfinite":false,"trimp":25.46704673479356,"ctl":[96.68685633705937,96.44160847368228,0.24524786337708804],"zones":[0.0,116.01159999999977,0.0,0.0,0.0],"circadian":[1.0,0.3312676056338072,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,17,0],"spectral":527.7181728812988,"cosinor":[60.20664182391473,1.993388211364034,0.990292034354844,2.520360883833827e-94],"near_singular":false,"hrv":[42.11090878105368,21.502568386886722,0.29644031390575065,21.528250446057537,-0.9685070047216077],"nonfinite":false,"trimp":28.750997708819092,"ctl":[98.08372430811947,101.79002204026516,-3.7062977321456856],"zones":[0.0,127.01269999999977,0.0,0.0,0.0],"circadian":[1.0,0.3688421052631628,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,31,0],"spectral":21.91751421630071,"cosinor":[60.206102619835434,1.9948632846663268,0.9904511511664523,1.168751936831454e-94],"near_singular":false,"hrv":[3.6114616155271952,21.50300002384371,0.916280807418444,21.477133240134368,0.9853180909453777],"nonfinite":false,"trimp":30.6254910760735,"ctl":[97.31583171818741,97.22062215666931,0.09520956151810367],"zones":[0.0,138.0137999999999,0.0,0.0,0.0],"circadian":[1.0,0.02086956521739154,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,32,0],"spectral":94.441233233861,"cosinor":[60.206607219531165,1.9918709847273792,0.9905159686967243,8.514781581412862e-95],"near_singular":false,"hrv":[14.550213162597531,21.50963654798735,0.2980266237242866,21.540300395940907,0.7627458277372915],"nonfinite":false,"trimp":33.68218231936987,"ctl":[98.15600470749584,99.32476832904142,-1.1687636215455797],"zones":[0.0,149.01489999999987,0.0,0.0,0.0],"circadian":[1.0,0.2553191489361734,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,33,0],"spectral":12.677409836990822,"cosinor":[60.20622592520504,1.9957534839219526,0.9905696980951562,6.537929637620507e-95],"near_singular":false,"hrv":[39.54331943857053,21.508583812826043,0.27113546544007344,21.551316007749993,-0.7371891534838562],"nonfinite":false,"trimp":35.83311260770479,"ctl":[97.86366868742407,98.31752817052501,-0.45385948310094193],"zones":[0.0,160.01599999999982,0.0,0.0,0.0],"circadian":[1.0,0.3312676056338072,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[200,33,0],"spectral":35.03266947297398,"cosinor":[60.20646575449843,1.9913371893547758,0.9904828062401032,1.0015309381130758e-94],"near_singular":false,"hrv":[42.11090878105368,21.502568386886722,0.29644031390575065,21.528250446057537,-0.9685070047216077],"nonfinite":false,"trimp":38.64417150766156,"ctl":[98.18187926518608,97.60954398514504,0.5723352800410453],"zones":[0.0,171.01709999999983,0.0,0.0,0.0],"circadian":[1.0,0.3688421052631628,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]}]} diff --git a/test/output_changes/batch_controls_test.dart b/test/output_changes/batch_controls_test.dart new file mode 100644 index 00000000..ebde99c9 --- /dev/null +++ b/test/output_changes/batch_controls_test.dart @@ -0,0 +1,42 @@ +import 'package:openstrap_analytics/onehz.dart'; +import 'package:test/test.dart'; + +void main() { + test('RR cap boundary preserves uncapped NN for short classified runs', () { + for (final n in [30, 31, 33]) { + final rr = [...List.filled(100, 900), + ...List.filled(n, 250), ...List.filled(100, 900)]; + final capped = correctRr(rr); + if (n <= 31) { + final old = correctRr(rr, maxArtifactRun: null); + expect(capped.nn, old.nn); + expect(capped.nnTimesMs, old.nnTimesMs); + expect(capped.classes, old.classes); + } else { + expect(capped.nn, List.filled(200, 900)); + expect(capped.droppedCount, 33); + } + } + }); + test('batch circular near-ties retain strict first-winner selection', () { + // Equal 10-epoch windows on an alternating profile tie at every start. + // A 5-epoch trough starts at zero, which includes three low epochs. + final profile = [for (var i = 0; i < 24; i++) i.isEven ? 1.0 : 2.0]; + final equal = circadianNonparametric(profile, 24).value!; + expect(equal.m10StartEpoch, 0); + expect(equal.l5StartEpoch, 0); + final near = [...profile]..[0] += 1e-12; + final result = circadianNonparametric(near, 24).value!; + expect(result.m10StartEpoch, 0); + expect(result.l5StartEpoch, 2); + }); + test('batch irregular screen excludes the gap pair already excluded by main', () { + final result = irregularBeatScreen([800, 900, 1000, 1100, 700, 850, 900, 1000], + nnTimesMs: [800, 1700, 2700, 3800, 14500, 15350, 16250, 17250], + minBeats: 4).value!; + // Six successive differences remain: 100,100,100,150,50,100. + // Five exceed 70 ms. The -400 ms difference across the gap is absent. + expect(result.pnnPct, 100 * 5 / 6); + expect(result.nBeats, 8); + }); +} diff --git a/test/output_changes/batch_parity_test.dart b/test/output_changes/batch_parity_test.dart new file mode 100644 index 00000000..d81a7b1e --- /dev/null +++ b/test/output_changes/batch_parity_test.dart @@ -0,0 +1,21 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + final config = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final golden = jsonDecode(File('test/fixtures/output_batch_parity.json').readAsStringSync()) as Map; + final cases = reproductionCases(config); + test('the full synthetic corpus is available', () { + expect(cases.length, 12); + expect(golden['cases'], hasLength(12)); + }); + for (var i = 0; i < 12; i++) { + test('split-2 public batch arithmetic survives refactoring, case $i', () { + // Frozen before aggregator adoption, not computed by the current path. + // Covers high-frequency jitter, time holes, cap boundaries and epochs. + expect(reproduce(cases[i]), golden['cases'][i], reason: 'case $i'); + }); + } +} diff --git a/tool/reproduce_outputs.dart b/tool/reproduce_outputs.dart index ea320f7c..93fdc2b4 100644 --- a/tool/reproduce_outputs.dart +++ b/tool/reproduce_outputs.dart @@ -6,6 +6,24 @@ import 'dart:math' as math; import 'package:openstrap_analytics/onehz.dart'; import 'package:openstrap_analytics/src/onehz/util.dart' as u; +// Deterministic differential cases, shared by the tool and parity tests. +List> reproductionCases(Map base) => [ + for (var i = 0; i < 12; i++) (() { + final f = jsonDecode(jsonEncode(base)) as Map; + f['hrv']['count'] = i < 6 ? 1800 : 5400; + f['hrv']['frequency'] = [.17, .7, 2.4, 2.9][i % 4]; + f['hrv']['gapped'] = i.isOdd; + f['spectral']['origin'] = i.isEven ? 1800000000.0 : 0.0; + f['spectral']['frequency'] = .037 + i * .019; + f['rr']['run'] = [1, 16, 30, 31, 32, 33][i % 6]; + f['cosinor']['step'] = .31 + i * .11; + f['load']['count'] = 30 + i * 7; + f['zones']['count'] = 50 + i * 11; + f['circadian']['days'] = 1 + i % 4; + return f; + })() +]; + Map reproduce(Map f) { final out = {}; final r = f['rr'] as Map; @@ -26,9 +44,9 @@ Map reproduce(Map f) { final n=f['near_singular'] as Map; final nt=[for(var i=0;i<(n['count'] as int);i++) (n['origin'] as num).toDouble()+i*(n['step'] as num)]; out['near_singular']=cosinor(nt,[for(var i=0;i0 && i%400==0 ? 10000 : 0)]; final hv=hrvTime(nn,nnTimesMs:times).value!; out['hrv']=[hv.rmssd,hv.sdnn,hv.sdann,hv.sdnnIndex,hv.diffAcf1]; try { @@ -51,6 +69,7 @@ Map reproduce(Map f) { } void main(List args) { final f=jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + if(args.contains('--corpus-json')) { print(jsonEncode([for(final c in reproductionCases(f)) reproduce(c)])); return; } final values=reproduce(f); if(args.contains('--json')) { print(jsonEncode(values)); return; } print('| Change | up/main 3a7ca34 | Current batch |'); From 3c21387cc1e6815e74384cf8842396ab3aa290ed Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Tue, 3 Nov 2026 16:04:48 -0700 Subject: [PATCH 11/14] Refuse malformed metric inputs and label capped RR gaps explicitly --- CHANGELOG.md | 11 +++ lib/src/onehz/clinical/cosinor.dart | 7 ++ lib/src/onehz/clinical/hrv_time.dart | 27 +++++++ lib/src/onehz/clinical/load_trimp.dart | 12 ++- lib/src/onehz/foundations/rr_correction.dart | 9 ++- .../foundations/rr_correction_kernel.dart | 2 +- test/output_changes/review_inputs_test.dart | 75 +++++++++++++++++++ 7 files changed, 138 insertions(+), 5 deletions(-) create mode 100644 test/output_changes/review_inputs_test.dart diff --git a/CHANGELOG.md b/CHANGELOG.md index 79caddc1..759a3132 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,17 @@ Edge must bump `kAlgoVersion` when adopting these changed batch outputs. +### Review input and gap validity + +HRV time, nocturnal and session RMSSD abstain on non-finite NN intervals before +segmentation. Cosinor rejects non-finite or non-positive periods before computing +its design. CTL/ATL abstains on unknown/non-finite daily load; explicit observed +zero-load rest days still decay. These are input-validity policies, not new +physiological models. Capped RR runs now report `BeatClass.gap`, appended after +existing enum values, rather than claiming every omitted beat is long/short. +Finite-input formulas, corrected NN and cap accounting are unchanged. + + ### Leave long RR artifact runs as gaps The 32-beat cap is an internal ESTIMATE rule with low confidence, not a published threshold. Longer runs split the series; each side is corrected independently. The synthetic 33-beat case keeps 200 NN beats and drops 33. diff --git a/lib/src/onehz/clinical/cosinor.dart b/lib/src/onehz/clinical/cosinor.dart index 595ccfa4..b6cf3ca4 100644 --- a/lib/src/onehz/clinical/cosinor.dart +++ b/lib/src/onehz/clinical/cosinor.dart @@ -79,6 +79,13 @@ Metric cosinor( double periodHours = 24, int minPoints = cosinorMinPoints, }) { + if (!periodHours.isFinite || periodHours <= 0) { + return const Metric.absent( + tier: Tier.high, + inputs_used: ['signal_timeseries'], + note: 'invalid cosinor period', + ); + } if (tHours.length != y.length) { return const Metric.absent(tier: Tier.high, inputs_used: ['signal_timeseries'], note: 'too few points for cosinor'); diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index 96c4c09d..76e8a192 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -414,6 +414,15 @@ Metric hrvTime( RrCoverage? coverage, }) { const inputs = ['rr_cleaned']; + // An unknown interval cannot support a finite HRV estimate. Refuse the + // input before segmentation or moments rather than discarding evidence. + if (nnMs.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite NN interval', + ); + } if (nnMs.length < 2) { return const Metric.absent( tier: Tier.high, @@ -626,6 +635,15 @@ Metric nocturnalRmssd( RrCoverage? coverage, }) { const inputs = ['rr_cleaned', 'beat_times']; + // An unknown interval cannot support a finite HRV estimate. Refuse the + // input before segmentation or moments rather than discarding evidence. + if (nnMs.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite NN interval', + ); + } if (coverage != null && coverage.overCounted) { return Metric.absent( tier: Tier.high, @@ -881,6 +899,15 @@ Metric sleepSessionRmssdDetail( int minDiffsPerWindow = kMinDiffsPerRmssdWindow, }) { const inputs = ['rr_sleep_window']; + // An unknown interval cannot support a finite HRV estimate. Refuse the + // input before segmentation or moments rather than discarding evidence. + if (rrMs.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite NN interval', + ); + } if (startSec <= 0 || endSec <= startSec || rrMs.isEmpty || diff --git a/lib/src/onehz/clinical/load_trimp.dart b/lib/src/onehz/clinical/load_trimp.dart index 89136ef5..53187efa 100644 --- a/lib/src/onehz/clinical/load_trimp.dart +++ b/lib/src/onehz/clinical/load_trimp.dart @@ -836,7 +836,8 @@ const int ctlAtlMinDays = 14; /// CTL/ATL/TSB from a time-ordered daily-TRIMP series (oldest→newest). /// EWMA with time constants 42 d (CTL) and 7 d (ATL): λ = 1 − e^(−1/τ). -/// A missing day contributes a 0-load impulse (rest day) — the EWMA decays. +/// An observed rest day is 0 load and decays the EWMA. Unknown/non-finite +/// daily load abstains; it is never substituted with a rest day. /// /// SEEDING (Banister 1975 impulse-response; the load before the record started /// is UNKNOWN): both accumulators used to be seeded at `dailyTrimp.first`, @@ -853,6 +854,13 @@ Metric ctlAtlTsb(List dailyTrimp, int minDays = ctlAtlMinDays, int primeDays = 7}) { const inputs = ['daily_trimp']; + if (dailyTrimp.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.estimate, + inputs_used: inputs, + note: 'non-finite daily TRIMP', + ); + } if (dailyTrimp.length < minDays) { return Metric.absent( tier: Tier.estimate, @@ -866,7 +874,6 @@ Metric ctlAtlTsb(List dailyTrimp, // Neumaier (1974) compensation for the bounded-prefix mean. var sum = 0.0, correction = 0.0; for (final value in dailyTrimp.take(prime)) { - if (value.isNaN) continue; final next = sum + value; if (next.isFinite) { correction += sum.abs() >= value.abs() @@ -882,7 +889,6 @@ Metric ctlAtlTsb(List dailyTrimp, final keep = 1 - lambda; for (var i = prime; i < dailyTrimp.length; i++) { final value = dailyTrimp[i]; - if (value.isNaN) continue; scale *= keep; offset = keep * offset + lambda * value; } diff --git a/lib/src/onehz/foundations/rr_correction.dart b/lib/src/onehz/foundations/rr_correction.dart index 41e69169..a8e1738f 100644 --- a/lib/src/onehz/foundations/rr_correction.dart +++ b/lib/src/onehz/foundations/rr_correction.dart @@ -31,7 +31,12 @@ import 'rr_correction_kernel.dart'; const int kMaxArtifactRun = 32; /// Artifact label for a single beat. -enum BeatClass { normal, ectopic, longShort, missed, extra } +enum BeatClass { + normal, ectopic, longShort, missed, extra, + /// An over-cap artifact run, without a claim about its individual types. + /// Appended to preserve the indices of existing checkpoint labels. + gap, +} class RrCorrectionResult { /// Cleaned NN intervals (ms). Isolated artifacts spline-corrected; multi-beat @@ -56,6 +61,8 @@ class RrCorrectionResult { final List nnSourceIndex; /// Per-input-beat classification (same length as the input rr). + /// Over-cap runs use [BeatClass.gap]; other artifact labels retain their + /// inferred classification. A gap makes no physiological type claim. final List classes; /// Fraction of input beats classified normal (0..1). diff --git a/lib/src/onehz/foundations/rr_correction_kernel.dart b/lib/src/onehz/foundations/rr_correction_kernel.dart index 8acb72b5..54ca09af 100644 --- a/lib/src/onehz/foundations/rr_correction_kernel.dart +++ b/lib/src/onehz/foundations/rr_correction_kernel.dart @@ -317,7 +317,7 @@ class RrCapWindow { } if (cutA >= 0) { for (var i = cutA; i < cutB; i++) { - out.classes[i] = BeatClass.longShort; + out.classes[i] = BeatClass.gap; out.gaps[i] = true; out.hardGaps[i] = hardCut || (inGap && hardLeading && cutA == start); } diff --git a/test/output_changes/review_inputs_test.dart b/test/output_changes/review_inputs_test.dart new file mode 100644 index 00000000..e65c3af5 --- /dev/null +++ b/test/output_changes/review_inputs_test.dart @@ -0,0 +1,75 @@ +import 'package:openstrap_analytics/onehz.dart'; +import 'package:test/test.dart'; + +void main() { + for (final bad in [double.nan, double.infinity, double.negativeInfinity]) { + test('HRV entry points abstain on non-finite NN interval $bad', () { + final nn = List.filled(1200, 900); + final ts = [for (var i = 0; i < nn.length; i++) 2000.0 + i * 900]; + nn[400] = bad; + for (final metric in [ + hrvTime(nn, nnTimesMs: ts), + nocturnalRmssd(nn, ts), + sleepSessionRmssdDetail(nn, ts, startSec: 1, endSec: 1200), + sleepSessionWindowedRmssd(nn, ts, startSec: 1, endSec: 1200), + ]) { + expect(metric.present, isFalse); + expect(metric.note, 'non-finite NN interval'); + } + expect(hrvTime(nn).present, isFalse); + }); + test('unknown training load $bad abstains in prime and later history', () { + for (final index in [0, 6, 7, 13]) { + final loads = List.filled(14, 50)..[index] = bad; + final metric = ctlAtlTsb(loads); + expect(metric.present, isFalse, reason: 'day $index'); + expect(metric.note, 'non-finite daily TRIMP'); + } + }); + } + test('a five-minute segment containing only NaNs abstains without throwing', () { + final nn = List.filled(1200, 900); + for (var i = 0; i < 334; i++) { nn[i] = double.nan; } + final ts = [for (var i = 0; i < nn.length; i++) i * 900.0]; + expect(hrvTime(nn, nnTimesMs: ts).present, isFalse); + }); + for (final period in [0.0, -24.0, double.nan, double.infinity, + double.negativeInfinity]) { + test('cosinor refuses invalid period $period', () { + final metric = cosinor( + [for (var i = 0; i < 48; i++) i.toDouble()], + [for (var i = 0; i < 48; i++) 60.0 + i % 3], + periodHours: period, + ); + expect(metric.present, isFalse); + expect(metric.note, 'invalid cosinor period'); + }); + } + test('over-cap artifacts have an explicit gap label', () { + final rr = [ + ...List.filled(100, 900), + ...List.filled(33, 100), + ...List.filled(100, 900), + ]; + final result = correctRr(rr); + expect(result.classes.sublist(100, 133).map((c) => c.name), + everyElement('gap')); + expect(result.droppedCount, 33); + expect(result.nn, List.filled(200, 900)); + expect(result.nnSourceIndex, + [...List.generate(100, (i) => i), + ...List.generate(100, (i) => i + 133)]); + }); + + test('observed rest days still decay load', () { + // Seven 50-load days prime at 50. Seven observed rest days then apply + // exp(-1/tau) seven times, giving 50*exp(-7/42) and 50*exp(-7/7). + final metric = ctlAtlTsb([ + ...List.filled(7, 50), ...List.filled(7, 0), + ]); + expect(metric.value!.ctl, closeTo(42.3240862445307, 1e-12)); + expect(metric.value!.atl, closeTo(18.3939720585721, 1e-12)); + expect(ctlAtlTsb(List.filled(14, 50)).value!.ctl, + closeTo(50, 1e-12)); + }); +} From 9eea67ce4923d077cd228f308a6f42c5f1ba71ca Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Tue, 3 Nov 2026 16:16:04 -0700 Subject: [PATCH 12/14] Refuse overlapping HRV clocks without shifting recorded segment times --- CHANGELOG.md | 5 +++- docs/OUTPUT_REPRODUCTION.md | 4 ++- lib/src/onehz/clinical/hrv_time.dart | 28 +++++++++++++++++++ .../foundations/rr_correction_kernel.dart | 3 ++ lib/src/onehz/util.dart | 11 ++++++++ test/output_changes/review_inputs_test.dart | 24 ++++++++++++++++ 6 files changed, 73 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 759a3132..634e71e1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,7 +12,10 @@ its design. CTL/ATL abstains on unknown/non-finite daily load; explicit observed zero-load rest days still decay. These are input-validity policies, not new physiological models. Capped RR runs now report `BeatClass.gap`, appended after existing enum values, rather than claiming every omitted beat is long/short. -Finite-input formulas, corrected NN and cap accounting are unchanged. +HRV also abstains when retained beat times regress, including when RR +accumulation disagrees with the recorded time at a capped-gap reset. Recorded +segment anchors are preserved; timestamps are never shifted to hide the conflict. +Other finite-input formulas, corrected NN and cap accounting are unchanged. ### Leave long RR artifact runs as gaps diff --git a/docs/OUTPUT_REPRODUCTION.md b/docs/OUTPUT_REPRODUCTION.md index f8eb4ff9..606088c2 100644 --- a/docs/OUTPUT_REPRODUCTION.md +++ b/docs/OUTPUT_REPRODUCTION.md @@ -3,7 +3,9 @@ This branch replaces the batch-output portion of #99. The aggregator branch stacks on this branch. Edge must bump `kAlgoVersion` when adopting these changed outputs. -Run `dart pub get --offline`, then `dart run tool/reproduce_outputs.dart`. +Run `dart pub get` with network access on a fresh clone, then +`dart run tool/reproduce_outputs.dart`. Once dependencies are cached, +`dart pub get --offline` also works. `--json` prints unrounded machine-readable results. The script calls batch APIs only. Its input recipes are in `test/fixtures/output_changes.json`, a small synthetic fixture. It contains no personal measurements. Formula parameters diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index 76e8a192..5465a26d 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -416,6 +416,16 @@ Metric hrvTime( const inputs = ['rr_cleaned']; // An unknown interval cannot support a finite HRV estimate. Refuse the // input before segmentation or moments rather than discarding evidence. + // A reset may expose disagreement between accumulated RR and recorded + // sensor time. Do not turn overlapping segments into successive beats. + if (nnTimesMs != null && nnTimesMs.length == nnMs.length && + hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } if (nnMs.any((value) => !value.isFinite)) { return const Metric.absent( tier: Tier.high, @@ -637,6 +647,15 @@ Metric nocturnalRmssd( const inputs = ['rr_cleaned', 'beat_times']; // An unknown interval cannot support a finite HRV estimate. Refuse the // input before segmentation or moments rather than discarding evidence. + // A reset may expose disagreement between accumulated RR and recorded + // sensor time. Do not turn overlapping segments into successive beats. + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } if (nnMs.any((value) => !value.isFinite)) { return const Metric.absent( tier: Tier.high, @@ -901,6 +920,15 @@ Metric sleepSessionRmssdDetail( const inputs = ['rr_sleep_window']; // An unknown interval cannot support a finite HRV estimate. Refuse the // input before segmentation or moments rather than discarding evidence. + // A reset may expose disagreement between accumulated RR and recorded + // sensor time. Do not turn overlapping segments into successive beats. + if (hasRegressingBeatTimes(rrTsMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } if (rrMs.any((value) => !value.isFinite)) { return const Metric.absent( tier: Tier.high, diff --git a/lib/src/onehz/foundations/rr_correction_kernel.dart b/lib/src/onehz/foundations/rr_correction_kernel.dart index 54ca09af..c0eeb7cf 100644 --- a/lib/src/onehz/foundations/rr_correction_kernel.dart +++ b/lib/src/onehz/foundations/rr_correction_kernel.dart @@ -360,6 +360,9 @@ class RrCapWindow { final dropout = wall && i > 0 && ts[i] - ts[i - 1] - rr[i] > reanchorGapMs; final anchored = wall ? ts[i] - base : 0.0; + // Independent segments restart on recorded time, even if invalid RR + // accumulated ahead of it. A monotonic clamp would invent timestamps. + // HRV consumers refuse any regression in the retained beat clock. cur = wall && reset ? anchored : dropout && anchored > cur diff --git a/lib/src/onehz/util.dart b/lib/src/onehz/util.dart index c847b5d2..e9ccd4a8 100644 --- a/lib/src/onehz/util.dart +++ b/lib/src/onehz/util.dart @@ -500,3 +500,14 @@ List benjaminiHochberg(List ps) { } return out; } + +/// Whether supplied beat times move backward. Equal times remain allowed: +/// band timestamps can be quantized to whole seconds. +bool hasRegressingBeatTimes(List times) { + for (var i = 1; i < times.length; i++) { + // Non-finite timestamps have their own refusal note at the call site. + if (times[i].isFinite && times[i - 1].isFinite && + times[i] < times[i - 1]) return true; + } + return false; +} diff --git a/test/output_changes/review_inputs_test.dart b/test/output_changes/review_inputs_test.dart index e65c3af5..598f3572 100644 --- a/test/output_changes/review_inputs_test.dart +++ b/test/output_changes/review_inputs_test.dart @@ -45,6 +45,30 @@ void main() { expect(metric.note, 'invalid cosinor period'); }); } + test('a capped reset keeps source time but invalidates regressing HRV time', () { + // The RR clock exceeds elapsed sensor time before the cap. Re-anchoring + // the next independent segment must retain its actual sensor timestamp. + final rr = [...List.filled(100, 1500), + ...List.filled(33, 3000), ...List.filled(100, 1500)]; + final ts = [for (var i = 0; i < rr.length; i++) 1.7e12 + i * 1000]; + final result = correctRr(rr, rrTsMs: ts); + expect(result.nnTimesMs[99], 150000); + expect(result.nnTimesMs[100], 134500); + for (final metric in [ + hrvTime(result.nn, nnTimesMs: result.nnTimesMs), + nocturnalRmssd(result.nn, result.nnTimesMs), + sleepSessionRmssdDetail(result.nn, result.nnTimesMs, + startSec: 1, endSec: 400), + ]) { + expect(metric.present, isFalse); + expect(metric.note, 'regressing beat time'); + } + }); + + test('an entirely unknown load history abstains', () { + expect(ctlAtlTsb(List.filled(42, double.nan)).present, isFalse); + }); + test('over-cap artifacts have an explicit gap label', () { final rr = [ ...List.filled(100, 900), From d059e6e8b0083f643d9e729f017bf5ac9e63aecf Mon Sep 17 00:00:00 2001 From: MikenRafter Date: Tue, 3 Nov 2026 17:20:11 -0700 Subject: [PATCH 13/14] Refuse invalid clocks across spectral consumers and non-finite cosinor designs --- CHANGELOG.md | 4 ++ lib/src/onehz/clinical/cosinor.dart | 10 ++++ lib/src/onehz/clinical/hrv_freq.dart | 7 +++ lib/src/onehz/clinical/hrv_time.dart | 7 +++ lib/src/onehz/respiration/cvhr_apnea.dart | 7 +++ lib/src/onehz/respiration/resp_rate.dart | 7 +++ test/output_changes/review_inputs_test.dart | 59 +++++++++++++++++++++ 7 files changed, 101 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 634e71e1..be0e6d9d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,10 @@ existing enum values, rather than claiming every omitted beat is long/short. HRV also abstains when retained beat times regress, including when RR accumulation disagrees with the recorded time at a capped-gap reset. Recorded segment anchors are preserved; timestamps are never shifted to hide the conflict. +Frequency HRV, RSA and CVHR use the same regressing-clock refusal. Time HRV +also refuses supplied clock vectors whose length does not match the NN vector; +NN-only calls remain supported. Cosinor refuses non-finite design/response +moments, including overflow from a positive subnormal period or finite responses. Other finite-input formulas, corrected NN and cap accounting are unchanged. diff --git a/lib/src/onehz/clinical/cosinor.dart b/lib/src/onehz/clinical/cosinor.dart index b6cf3ca4..eaa8f5b6 100644 --- a/lib/src/onehz/clinical/cosinor.dart +++ b/lib/src/onehz/clinical/cosinor.dart @@ -149,6 +149,16 @@ Metric _finishCosinor( note: 'too few points for cosinor', ); } + // A finite period or response can still overflow the design arithmetic. + // Reject invalid moments before NaN can bypass the identifiability gate. + if (![state.meanC, state.meanS, state.meanY, state.mCC, state.mSS, + state.mYY, state.mCS, state.mCY, state.mSY].every((v) => v.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite cosinor design', + ); + } final determinant = state.mCC * state.mSS - state.mCS * state.mCS; // The determinant is relative to the diagonal co-moments. A design needs // at least 1e-12 of that scale after removing its means to identify both diff --git a/lib/src/onehz/clinical/hrv_freq.dart b/lib/src/onehz/clinical/hrv_freq.dart index 2190a1d3..1cc57995 100644 --- a/lib/src/onehz/clinical/hrv_freq.dart +++ b/lib/src/onehz/clinical/hrv_freq.dart @@ -91,6 +91,13 @@ Metric hrvFreq( note: 'non-finite beat time', ); } + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { return const Metric.absent( diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index 5465a26d..e34f0477 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -416,6 +416,13 @@ Metric hrvTime( const inputs = ['rr_cleaned']; // An unknown interval cannot support a finite HRV estimate. Refuse the // input before segmentation or moments rather than discarding evidence. + if (nnTimesMs != null && nnTimesMs.length != nnMs.length) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'beat time length mismatch', + ); + } // A reset may expose disagreement between accumulated RR and recorded // sensor time. Do not turn overlapping segments into successive beats. if (nnTimesMs != null && nnTimesMs.length == nnMs.length && diff --git a/lib/src/onehz/respiration/cvhr_apnea.dart b/lib/src/onehz/respiration/cvhr_apnea.dart index e4046b8b..348c2aff 100644 --- a/lib/src/onehz/respiration/cvhr_apnea.dart +++ b/lib/src/onehz/respiration/cvhr_apnea.dart @@ -134,6 +134,13 @@ Metric cvhrApneaScreen( note: 'non-finite beat time', ); } + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } if (tSec.last - tSec.first <= 0) { return const Metric.absent( tier: Tier.high, diff --git a/lib/src/onehz/respiration/resp_rate.dart b/lib/src/onehz/respiration/resp_rate.dart index 377b0ee3..8e57c1d9 100644 --- a/lib/src/onehz/respiration/resp_rate.dart +++ b/lib/src/onehz/respiration/resp_rate.dart @@ -194,6 +194,13 @@ Metric rsaRespRate( note: 'non-finite beat time', ); } + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { return const Metric.absent( diff --git a/test/output_changes/review_inputs_test.dart b/test/output_changes/review_inputs_test.dart index 598f3572..4e3a3dae 100644 --- a/test/output_changes/review_inputs_test.dart +++ b/test/output_changes/review_inputs_test.dart @@ -59,12 +59,71 @@ void main() { nocturnalRmssd(result.nn, result.nnTimesMs), sleepSessionRmssdDetail(result.nn, result.nnTimesMs, startSec: 1, endSec: 400), + sleepSessionRmssdDetail(result.nn, result.nnTimesMs, + startSec: 1, endSec: 100), ]) { expect(metric.present, isFalse); expect(metric.note, 'regressing beat time'); } }); + for (final entry in Function(List, List)>{ + 'frequency HRV': (nn, ts) => hrvFreq(nn, ts, artifactFraction: 0), + 'RSA': (nn, ts) => rsaRespRate(nn, ts, artifactFraction: 0), + 'CVHR': (nn, ts) => cvhrApneaScreen(nn, ts, artifactFraction: 0), + }.entries) { + test('${entry.key} refuses a regressing spectral clock', () { + final nn = List.filled(1500, 900); + final ts = [for (var i = 0; i < nn.length; i++) i * 900.0]; + ts[700] = ts[699] - 500; + final metric = entry.value(nn, ts); + expect(metric.present, isFalse); + expect(metric.note, 'regressing beat time'); + }); + } + + for (final times in >[[], [800], [800, 1700, 2500]]) { + test('time HRV refuses a supplied ${times.length}-beat clock for two NN', () { + final metric = hrvTime([800, 900], nnTimesMs: times); + expect(metric.present, isFalse); + expect(metric.note, 'beat time length mismatch'); + // NN-only remains supported, as does a matching two-beat clock. + // One 100 ms difference has RMSSD exactly 100 ms. + expect(hrvTime([800, 900]).value!.rmssd, 100); + expect(hrvTime([800, 900], nnTimesMs: [800, 1700]).value!.rmssd, 100); + }); + } + + test('cosinor refuses overflow from a subnormal positive period', () { + final metric = cosinor( + List.generate(48, (i) => i.toDouble()), + List.generate(48, (i) => 60.0 + i % 3), + periodHours: 1e-320, + ); + expect(metric.present, isFalse); + expect(metric.note, 'non-finite cosinor design'); + }); + for (final bad in [double.nan, double.infinity, double.negativeInfinity]) { + for (final axis in ['time', 'response']) { + test('cosinor refuses non-finite $axis moment $bad', () { + final ts = List.generate(48, (i) => i.toDouble()); + final values = List.generate(48, (i) => 60.0 + i % 3); + (axis == 'time' ? ts : values)[12] = bad; + final metric = cosinor(ts, values); + expect(metric.present, isFalse); + expect(metric.note, 'non-finite cosinor design'); + }); + } + } + test('cosinor refuses non-finite moments from finite overflow', () { + final metric = cosinor( + List.generate(48, (i) => i.toDouble()), + List.generate(48, (i) => i.isEven ? 1e308 : -1e308), + ); + expect(metric.present, isFalse); + expect(metric.note, 'non-finite cosinor design'); + }); + test('an entirely unknown load history abstains', () { expect(ctlAtlTsb(List.filled(42, double.nan)).present, isFalse); }); From 993b0787fc634631b007c50cb70d5b7e0e3fe961 Mon Sep 17 00:00:00 2001 From: abdulsaheel Date: Fri, 13 Nov 2026 14:30:51 +0530 Subject: [PATCH 14/14] frequency hrv: leave out beats behind a re-anchored clock instead of blanking the record a capped artifact reset puts the next segment behind the rr clock of the one before it, and refusing the whole night over that threw away the spectrum. the beats behind the latest time are left out (a few percent at most, otherwise it still refuses), nothing is interpolated. also drops a private-data comment and stops walking the series once a non-finite value makes the variance nan. --- lib/src/onehz/clinical/hrv_freq.dart | 54 ++++++++++++++++---- lib/src/onehz/foundations/rr_correction.dart | 6 +-- lib/src/onehz/sleep/circadian_np.dart | 5 +- test/output_changes/review_inputs_test.dart | 36 ++++++++++++- 4 files changed, 84 insertions(+), 17 deletions(-) diff --git a/lib/src/onehz/clinical/hrv_freq.dart b/lib/src/onehz/clinical/hrv_freq.dart index 1cc57995..c8625e06 100644 --- a/lib/src/onehz/clinical/hrv_freq.dart +++ b/lib/src/onehz/clinical/hrv_freq.dart @@ -61,6 +61,10 @@ class HrvFreq { }; } +/// Largest share of beats [hrvFreq] will leave out for sitting behind the +/// record's latest time before it refuses the clock instead. +const double kMaxRegressingBeatFraction = 0.05; + /// Frequency-domain HRV from cleaned NN. /// /// [nnMs] cleaned NN intervals (ms). [nnTimesMs] their cumulative beat times @@ -83,7 +87,7 @@ Metric hrvFreq( } // Times in seconds for Hz output; tachogram value = NN (ms). - final tSec = [for (final t in nnTimesMs) t / 1000.0]; + var tSec = [for (final t in nnTimesMs) t / 1000.0]; if (tSec.any((t) => !t.isFinite)) { return const Metric.absent( tier: Tier.high, @@ -91,12 +95,37 @@ Metric hrvFreq( note: 'non-finite beat time', ); } + // A capped artifact reset re-anchors the next segment at its source time, + // which can sit behind the cumulative RR clock of the segment before it (the + // clock ran ahead while the run was being absorbed). Refusing the whole + // record for that blanked a night over one 30 s dropout. The beats that fall + // behind the latest time already seen are the ones whose position is not + // known to be after their predecessors, so they are left out; the rest keep + // their own times and the periodogram is irregular-sampling by construction, + // so nothing is interpolated or invented. A clock that goes backward for more + // than a few percent of the record is not a seam, and still refuses. + var nn = nnMs; + var droppedBeats = 0; if (hasRegressingBeatTimes(nnTimesMs)) { - return const Metric.absent( - tier: Tier.high, - inputs_used: inputs, - note: 'regressing beat time', - ); + final keep = []; + var latest = double.negativeInfinity; + for (var i = 0; i < tSec.length; i++) { + if (tSec[i] >= latest) { + keep.add(i); + latest = tSec[i]; + } + } + droppedBeats = tSec.length - keep.length; + if (droppedBeats > tSec.length * kMaxRegressingBeatFraction || + keep.length < 16) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } + nn = [for (final i in keep) nnMs[i]]; + tSec = [for (final i in keep) tSec[i]]; } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { @@ -112,10 +141,10 @@ Metric hrvFreq( // null — absent, not a leakage-filled 0.0 (ULF used to be emitted as exactly // 0.0 for any session over 333 s, which is not a 24-h record by any reading). final ulf = - _welchBandPower(tSec, nnMs, 0.0003, 0.003, oversample: oversample); - final vlf = _welchBandPower(tSec, nnMs, 0.003, 0.04, oversample: oversample); - final lf = _welchBandPower(tSec, nnMs, 0.04, 0.15, oversample: oversample); - final hfRaw = _welchBandPower(tSec, nnMs, 0.15, 0.40, oversample: oversample); + _welchBandPower(tSec, nn, 0.0003, 0.003, oversample: oversample); + final vlf = _welchBandPower(tSec, nn, 0.003, 0.04, oversample: oversample); + final lf = _welchBandPower(tSec, nn, 0.04, 0.15, oversample: oversample); + final hfRaw = _welchBandPower(tSec, nn, 0.15, 0.40, oversample: oversample); if (lf == null && hfRaw == null) { return const Metric.absent( tier: Tier.high, @@ -184,7 +213,10 @@ Metric hrvFreq( note: hfGated ? 'HF suppressed: artifact fraction ${round6(artifactFraction)} ' '> gate — LF/VLF reported, HF/LF-HF/nu/total withheld' - : 'PRV spectrum; HF band quantization-limited at 1 Hz', + : droppedBeats > 0 + ? 'PRV spectrum; $droppedBeats beats behind a re-anchored clock ' + 'left out; HF band quantization-limited at 1 Hz' + : 'PRV spectrum; HF band quantization-limited at 1 Hz', ); } diff --git a/lib/src/onehz/foundations/rr_correction.dart b/lib/src/onehz/foundations/rr_correction.dart index a8e1738f..455a5f38 100644 --- a/lib/src/onehz/foundations/rr_correction.dart +++ b/lib/src/onehz/foundations/rr_correction.dart @@ -25,9 +25,9 @@ import 'rr_correction_kernel.dart'; /// Longest artifact run (beats) the corrector will carry. A longer run is not a /// correction problem but an absence: it breaks the series into independent -/// segments and its beats are left out (a gap). 32 beats leaves no run on the -/// owner's 18 days of data over the cap (longest seen: 16) and bounds the -/// streaming state; see combined-bump/brief.md. +/// segments and its beats are left out (a gap). 32 beats is an internal +/// ESTIMATE: comfortably above the longest artifact run seen in real +/// recordings (16 beats) while bounding the streaming state. const int kMaxArtifactRun = 32; /// Artifact label for a single beat. diff --git a/lib/src/onehz/sleep/circadian_np.dart b/lib/src/onehz/sleep/circadian_np.dart index d5365969..e37de96b 100644 --- a/lib/src/onehz/sleep/circadian_np.dart +++ b/lib/src/onehz/sleep/circadian_np.dart @@ -80,9 +80,10 @@ Metric circadianNonparametric( for (final value in x) { count++; if (!value.isFinite) { + // NaN is absorbing: every later update keeps ssTot NaN, so the running + // centre is not needed again and the rest of the series is not walked. ssTot = double.nan; - center = x.take(count).fold(0, (a, b) => a + b) / count; - continue; + break; } final delta = value - center; center += delta / count; diff --git a/test/output_changes/review_inputs_test.dart b/test/output_changes/review_inputs_test.dart index 4e3a3dae..68ec3b12 100644 --- a/test/output_changes/review_inputs_test.dart +++ b/test/output_changes/review_inputs_test.dart @@ -68,7 +68,6 @@ void main() { }); for (final entry in Function(List, List)>{ - 'frequency HRV': (nn, ts) => hrvFreq(nn, ts, artifactFraction: 0), 'RSA': (nn, ts) => rsaRespRate(nn, ts, artifactFraction: 0), 'CVHR': (nn, ts) => cvhrApneaScreen(nn, ts, artifactFraction: 0), }.entries) { @@ -82,6 +81,41 @@ void main() { }); } + group('frequency HRV and a clock that steps backward', () { + test('a capped reset leaves out the beats behind the re-anchored clock', () { + // 1,500 ms intervals for 400 beats, a 33-beat artifact run that forces + // the cap, then 400 more: the second segment is re-anchored at its source + // time, which sits behind the first segment's cumulative RR clock. + double wob(int i) => 1500.0 + (i % 7) * 12 - (i % 3) * 9; + final rr = [for (var i = 0; i < 400; i++) wob(i), + ...List.filled(33, 3000), for (var i = 0; i < 400; i++) wob(i)]; + final ts = [for (var i = 0; i < rr.length; i++) 1.7e12 + i * 1400]; + final result = correctRr(rr, rrTsMs: ts); + expect(hasRegressingBeatTimes(result.nnTimesMs), isTrue); + final metric = hrvFreq(result.nn, result.nnTimesMs, artifactFraction: 0); + expect(metric.present, isTrue); + expect(metric.note, contains('left out')); + }); + + test('one stray beat is left out, the rest keep their own times', () { + final nn = [for (var i = 0; i < 1500; i++) 900.0 + (i % 7) * 20 + (i % 3) * 15]; + final ts = [for (var i = 0; i < nn.length; i++) i * 900.0]; + ts[700] = ts[699] - 500; + final metric = hrvFreq(nn, ts, artifactFraction: 0); + expect(metric.present, isTrue); + expect(metric.note, contains('1 beats')); + }); + + test('a clock that stays behind for a tenth of the record still refuses', () { + final nn = [for (var i = 0; i < 1500; i++) 900.0 + (i % 7) * 20 + (i % 3) * 15]; + final ts = [for (var i = 0; i < nn.length; i++) i * 900.0]; + for (var i = 700; i < 900; i++) { ts[i] = ts[699] - 5000 + (i - 700); } + final metric = hrvFreq(nn, ts, artifactFraction: 0); + expect(metric.present, isFalse); + expect(metric.note, 'regressing beat time'); + }); + }); + for (final times in >[[], [800], [800, 1700, 2500]]) { test('time HRV refuses a supplied ${times.length}-beat clock for two NN', () { final metric = hrvTime([800, 900], nnTimesMs: times);