diff --git a/lib/src/onehz/sleep/cardio_stager.dart b/lib/src/onehz/sleep/cardio_stager.dart index 7a7957f..a6e97d1 100644 --- a/lib/src/onehz/sleep/cardio_stager.dart +++ b/lib/src/onehz/sleep/cardio_stager.dart @@ -366,8 +366,9 @@ void resetCardioObservations() => _cardioObservations.clear(); /// [hr1hz] per-second HR (bpm; 0 = off-skin) over the in-bed window. /// [accel] per-second gravity vectors, SAME length/time base as [hr1hz]. /// [rrMs] / [rrTsMs] beat-to-beat RR (ms) and their ABSOLUTE times (ms), same -/// clock as `accel[i].tsMs`. Sparse/empty is fine — REM/deep just lean more on -/// HR and confidence drops. +/// clock as `accel[i].tsMs`, ASCENDING in [rrTsMs] (the per-epoch window +/// gather lower-bounds into it). Sparse/empty is fine — REM/deep just lean +/// more on HR and confidence drops. CardioStagerResult cardioStager( List hr1hz, List accel, { @@ -674,11 +675,14 @@ CardioStagerResult classifyCardioEpochs( for (var k = lo; k < hi; k++) if (still(k) && !hr[k].isNaN) hr[k] ]; - final m = median(win); + // One sort for both percentiles of this 361-epoch window (stddev doesn't + // care about order). `win` is freshly built here, so sorting it is local. + win.sort(); + final m = percentileSorted(win, 50); if (m != null) { hrMedLocal[e] = m; hrArousalLocal[e] = m + math.max(6.0, (stddev(win) ?? 6)); - hrP25Local[e] = percentile(win, 25) ?? m; + hrP25Local[e] = percentileSorted(win, 25) ?? m; } } // Blend the per-epoch LOCAL HR gates toward the sleeper's rolling profile. @@ -872,13 +876,14 @@ CardioStagerResult classifyCardioEpochs( // Only when the edge armed recording AND the staged span is a real sleep // (≥60 min) — naps must not pollute the sleeper's night profile. if (cardioRecordObservations && nEpoch >= 120 && sleepHr.isNotEmpty) { + final hrSleepMed = median(sleepHr); _cardioObservations.add(SleepNightObservation( epochs: nEpoch, hrFloorP5: percentile(sleepHr, 5), hrFloorP25: percentile(sleepHr, 25), - hrSleepMedian: median(sleepHr), - hrArousal: (median(sleepHr) ?? hrMedGlobal) + - math.max(6.0, stddev(sleepHr) ?? 6.0), + hrSleepMedian: hrSleepMed, + hrArousal: + (hrSleepMed ?? hrMedGlobal) + math.max(6.0, stddev(sleepHr) ?? 6.0), rmssdMed: rmssdMed, rmssdMad: mad(sleepRmssd), enmoStillCut: motMed + 1.5 * motMad, @@ -950,30 +955,7 @@ CardioStagerResult _abstain(int epochSec) => CardioStagerResult( /// window centred on epoch [s,t). Returns NaN when too few clean beats. double _windowRmssd(List rrMs, List rrTsMs, List accel, int s, int t, int epochSec) { - if (rrMs.isEmpty || rrTsMs.length != rrMs.length) return double.nan; - final mid = (s + t) ~/ 2; - if (mid >= accel.length) return double.nan; - final centreMs = accel[mid].tsMs; - const halfWinMs = 150 * 1000; // ±2.5 min for a stable RMSSD on sparse RR - final lo = centreMs - halfWinMs, hi = centreMs + halfWinMs; - // Gather clean beats in window (300–2000 ms, drop big successive jumps). - final beats = []; - double? prev; - for (var i = 0; i < rrMs.length; i++) { - final ts = rrTsMs[i]; - if (ts < lo || ts > hi) continue; - final v = rrMs[i]; - if (v < _rrMin || v > _rrMax) { - prev = null; - continue; - } - if (prev != null && (v - prev).abs() > _rrMaxStep) { - prev = v; - continue; - } - beats.add(v); - prev = v; - } + final beats = _cleanBeatsInWindow(rrMs, rrTsMs, accel, s, t).beats; if (beats.length < 5) return double.nan; var ss = 0.0; for (var i = 1; i < beats.length; i++) { @@ -995,7 +977,7 @@ double _windowRmssd(List rrMs, List rrTsMs, /// flat RMSSD of 67.3 / 70.7 / 79.7. double _windowSdnn(List rrMs, List rrTsMs, List accel, int s, int t, int epochSec) { - final beats = _cleanBeatsInWindow(rrMs, rrTsMs, accel, s, t); + final beats = _cleanBeatsInWindow(rrMs, rrTsMs, accel, s, t).beats; if (beats.length < 5) return double.nan; final m = mean(beats)!; var ss = 0.0; @@ -1005,22 +987,53 @@ double _windowSdnn(List rrMs, List rrTsMs, return math.sqrt(ss / (beats.length - 1)); } -/// Clean RR beats (ms) inside the ±2.5-min window centred on epoch [s,t). -/// Shared by [_windowRmssd] and [_windowSdnn] so both see exactly the same -/// beats — a divergence there would make their z-scores incomparable. -List _cleanBeatsInWindow(List rrMs, List rrTsMs, - List accel, int s, int t) { - if (rrMs.isEmpty || rrTsMs.length != rrMs.length) return const []; +/// Clean RR beats (ms) inside a ±[halfWinMs] window centred on epoch [s,t), +/// with each beat's time rebased to the window start (seconds). +/// +/// The physiologic gate (300–2000 ms, drop big successive jumps) lives here and +/// only here, so every window feature sees exactly the same beats for a given +/// window — [_windowRmssd] and [_windowSdnn] are computed over one window and +/// must agree on its contents. +/// +/// [_windowRemFeatures] passes a DIFFERENT [halfWinMs] on purpose (±90 s, not +/// ±2.5 min): its LF/HF grid and its `beats.length` abstain gate are specified +/// against that shorter window. Do not hand it the default list. +({List beats, List tsSec}) _cleanBeatsInWindow( + List rrMs, + List rrTsMs, + List accel, + int s, + int t, { + int halfWinMs = 150 * 1000, +}) { + const empty = (beats: [], tsSec: []); + if (rrMs.isEmpty || rrTsMs.length != rrMs.length) return empty; final mid = (s + t) ~/ 2; - if (mid >= accel.length) return const []; + if (mid >= accel.length) return empty; final centreMs = accel[mid].tsMs; - const halfWinMs = 150 * 1000; final lo = centreMs - halfWinMs, hi = centreMs + halfWinMs; + // [rrTsMs] is ascending, so lower-bound into it and stop at the far edge + // instead of walking every beat of the night once per epoch per feature. + // + // TRUE lower bound — the first index with ts >= lo, NOT the first index + // strictly greater than lo. `rr_ts_ms` is `rec_ts * 1000`, so beats TIE on + // the second boundary; dropping the tied beats at a window edge changes + // RMSSD / SDNN / LF-HF and with them the REM and deep epoch counts, silently. + var a = 0, b = rrTsMs.length; + while (a < b) { + final m = (a + b) >> 1; + if (rrTsMs[m] < lo) { + a = m + 1; + } else { + b = m; + } + } final beats = []; + final tsSec = []; double? prev; - for (var i = 0; i < rrMs.length; i++) { + for (var i = a; i < rrMs.length; i++) { final ts = rrTsMs[i]; - if (ts < lo || ts > hi) continue; + if (ts > hi) break; final v = rrMs[i]; if (v < _rrMin || v > _rrMax) { prev = null; @@ -1031,11 +1044,33 @@ List _cleanBeatsInWindow(List rrMs, List rrTsMs, continue; } beats.add(v); + // Rebase to the window start (lo), NOT absolute epoch ms. Lomb–Scargle is + // time-shift invariant (the τ phase reference cancels any offset), so the + // LF/HF output is unchanged — but this keeps beat times in [0, 180] s + // instead of ~1.75e9 s. Absolute epoch seconds force every sin/cos in the + // periodogram onto libm's __kernel_rem_pio2 multi-precision slow path + // (args ~9e9 rad), which — run per 30-s epoch over a full night on the + // main isolate — caused main-thread ANRs on Android (Crashlytics: libm.so + // __kernel_rem_pio2 / sin / cos, 0.9.13). + tsSec.add((ts - lo) / 1000.0); prev = v; } - return beats; + return (beats: beats, tsSec: tsSec); } +/// [_cleanBeatsInWindow]'s beats, exposed so the window-edge regression test +/// can drive the gather directly. Not a staging entry point — call +/// [cardioStager]. +List cleanBeatsInWindowForTest( + List rrMs, + List rrTsMs, + List accel, + int s, + int t, { + int halfWinMs = 150 * 1000, +}) => + _cleanBeatsInWindow(rrMs, rrTsMs, accel, s, t, halfWinMs: halfWinMs).beats; + /// Webster sleep-continuity rescore: brief wake bouts flanked by enough sleep /// are re-labelled sleep (NREM). This is the published actigraphy step that /// prevents normal in-sleep repositioning from inflating WASO. @@ -1135,41 +1170,12 @@ void _websterRescore(List sm, int epochSec) { /// nulls when too few clean beats for a stable estimate. ({double? lfhf, double? rk}) _windowRemFeatures(List rrMs, List rrTsMs, List accel, int s, int t, int epochSec) { - if (rrMs.isEmpty || rrTsMs.length != rrMs.length) { - return (lfhf: null, rk: null); - } - final mid = (s + t) ~/ 2; - if (mid >= accel.length) return (lfhf: null, rk: null); - final centreMs = accel[mid].tsMs; - const halfWinMs = 90 * 1000; // ±90 s per the REM feature spec - final lo = centreMs - halfWinMs, hi = centreMs + halfWinMs; - final beats = []; // clean RR (ms) - final beatTsSec = []; // matching beat times (s) - double? prev; - for (var i = 0; i < rrMs.length; i++) { - final ts = rrTsMs[i]; - if (ts < lo || ts > hi) continue; - final v = rrMs[i]; - if (v < _rrMin || v > _rrMax) { - prev = null; - continue; - } - if (prev != null && (v - prev).abs() > _rrMaxStep) { - prev = v; - continue; - } - beats.add(v); - // Rebase to the window start (lo), NOT absolute epoch ms. Lomb–Scargle is - // time-shift invariant (the τ phase reference cancels any offset), so the - // LF/HF output is unchanged — but this keeps beat times in [0, 180] s - // instead of ~1.75e9 s. Absolute epoch seconds force every sin/cos in the - // periodogram onto libm's __kernel_rem_pio2 multi-precision slow path - // (args ~9e9 rad), which — run per 30-s epoch over a full night on the - // main isolate — caused main-thread ANRs on Android (Crashlytics: libm.so - // __kernel_rem_pio2 / sin / cos, 0.9.13). - beatTsSec.add((ts - lo) / 1000.0); - prev = v; - } + // ±90 s per the REM feature spec — a DIFFERENT window from the RMSSD/SDNN + // one; see [_cleanBeatsInWindow]. + final win = _cleanBeatsInWindow(rrMs, rrTsMs, accel, s, t, + halfWinMs: 90 * 1000); + final beats = win.beats; // clean RR (ms) + final beatTsSec = win.tsSec; // matching beat times (s), rebased to window if (beats.length < 16) return (lfhf: null, rk: null); // spectral stability gate // R(k): mean absolute successive difference of instantaneous HR (bpm). diff --git a/lib/src/onehz/util.dart b/lib/src/onehz/util.dart index ff584fc..5b73896 100644 --- a/lib/src/onehz/util.dart +++ b/lib/src/onehz/util.dart @@ -40,10 +40,14 @@ double? stddevPop(List xs) { return math.sqrt(s / xs.length); } -/// Linear-interpolated percentile (p in [0,100]); null if empty. -double? percentile(List values, double p) { - if (values.isEmpty) return null; - final sorted = [...values]..sort(); +/// Linear-interpolated percentile (p in [0,100]) of an ALREADY-ASCENDING list; +/// null if empty. For when you hold a sorted list already, or want several +/// percentiles of one window without re-sorting it each time. +/// +/// Null on empty, never 0 — a 0.0 is a measurement, and an empty window is the +/// absence of one. +double? percentileSorted(List sorted, double p) { + if (sorted.isEmpty) return null; if (sorted.length == 1) return sorted[0]; final rank = (p / 100) * (sorted.length - 1); final lo = rank.floor(); @@ -53,6 +57,10 @@ double? percentile(List values, double p) { return sorted[lo] + (sorted[hi] - sorted[lo]) * frac; } +/// Linear-interpolated percentile (p in [0,100]); null if empty. +double? percentile(List values, double p) => + percentileSorted([...values]..sort(), p); + /// Median (linear-interpolated), or null if empty. double? median(List xs) => percentile(xs, 50); diff --git a/test/onehz/real_night_cardio_stager_test.dart b/test/onehz/real_night_cardio_stager_test.dart index c1d470c..7ac34cf 100644 --- a/test/onehz/real_night_cardio_stager_test.dart +++ b/test/onehz/real_night_cardio_stager_test.dart @@ -45,6 +45,8 @@ import 'package:test/test.dart'; import 'package:openstrap_analytics/onehz.dart'; void main() { + _windowEdgeTies(); + test('cardioStager on the real 2026-07 overnight capture matches Apple ' 'Watch ground truth within a wide band (regression: was wake=294min ' 'rem=41min pre-fix)', () { @@ -135,3 +137,37 @@ void main() { expect(tot, closeTo(stages.length * epochMin, 0.01)); }); } + +// The per-epoch RR window is found by binary search. `rr_ts_ms` is +// `rec_ts * 1000`, so beat timestamps TIE — several beats share one second, +// and a second boundary can land exactly on a window edge. The search must be +// a TRUE lower bound (first index >= lo) and the far edge must be inclusive; +// "first index strictly greater than lo" drops the tied beats at the edge, +// which moves RMSSD / SDNN / LF-HF and with them the REM and deep counts. +void _windowEdgeTies() { + test('window gather keeps tied beats at BOTH edges', () { + final accel = [ + for (var i = 0; i < 30; i++) AccelSample(i * 1000.0, 0, 0, 1) + ]; + // epoch [0,30) -> mid 15 -> centre 15000 ms; ±5 s -> window [10000, 20000]. + final ts = [ + 9000, // just before -> out + 10000, 10000, 10000, // tie ON the low edge -> in + 15000, + 20000, 20000, 20000, // tie ON the high edge -> in + 21000, // just after -> out + ]; + final rr = [for (var i = 0; i < ts.length; i++) 900.0]; + + final beats = + cleanBeatsInWindowForTest(rr, ts, accel, 0, 30, halfWinMs: 5000); + // 7, not 4 (low ties dropped) and not 5 (both edges exclusive). + expect(beats.length, 7); + + // and the beats really are the in-window ones, in order. + final marked = [900, 901, 902, 903, 904, 905, 906, 907, 908]; + final got = + cleanBeatsInWindowForTest(marked, ts, accel, 0, 30, halfWinMs: 5000); + expect(got, [901, 902, 903, 904, 905, 906, 907]); + }); +} diff --git a/test/onehz/util_test.dart b/test/onehz/util_test.dart index 177c50a..49d6e55 100644 --- a/test/onehz/util_test.dart +++ b/test/onehz/util_test.dart @@ -16,6 +16,15 @@ void main() { expect(stddev([5]), isNull); }); + test('percentileSorted matches percentile and stays null on empty', () { + final sorted = [1.0, 2.0, 3.0, 4.0]; + expect(percentileSorted(sorted, 50), percentile(sorted, 50)); + expect(percentileSorted(sorted, 25), percentile(sorted, 25)); + expect(percentileSorted([7.0], 25), 7.0); + // NOT 0 — an empty window is no measurement, not a 0 bpm floor. + expect(percentileSorted([], 25), isNull); + }); + test('MAD and robust z; MAD=0 guard on quantized data', () { // [1,1,1,1,1] -> median 1, MAD 0 => robustZ null (not div-by-zero). expect(mad([1, 1, 1, 1, 1], scaled: false), 0);