diff --git a/lib/src/onehz/sleep/cardio_stager.dart b/lib/src/onehz/sleep/cardio_stager.dart index 9b25244..417eee4 100644 --- a/lib/src/onehz/sleep/cardio_stager.dart +++ b/lib/src/onehz/sleep/cardio_stager.dart @@ -18,15 +18,32 @@ // Per 30-s epoch we measure, over the in-bed window only: // motion = mean ENMO (van Hees 2013 amplitude index, g), against a LOCALLY // re-estimated 1 g reference (see below), NOT a whole-night scalar -// hr = mean valid HR (bpm) -// rmssd = RMSSD of cleaned RR beats in a ±2.5-min window (ms), or null -// then classify against baselines: +// hr = mean valid HR (bpm) hrSd = within-epoch SD of HR (bpm) +// rmssd, sdnn = over a ±2.5-min RR window (ms) +// lfhf, Rk = over a ±90-s RR window +// then classify against the night's own baselines: // WAKE : clearly elevated motion OR HR arousal above a LOCAL sleeping HR // median (a 90-min rolling window, not the whole night) -// REM : still body + RMSSD well BELOW the night's sleep RMSSD + HR ≥ a -// LOCAL p25 HR floor (see below for why p25, not median) -// DEEP : HR near the night floor + RMSSD high + very stable (NREM subtype) +// REM : weighted score over Rk / sdnn / hrSd / lfhf above a cutoff, GATED +// by atonia (no big movement) and HR ≥ a LOCAL p25 floor +// DEEP : weighted score over the same axes in the opposite direction +// (NREM subtype, low confidence) // LIGHT : remaining NREM +// +// 2026-08 — WHY SCORES, NOT CONJUNCTIONS, AND WHY NOT RMSSD. The rules above +// replace an earlier design that AND-ed boolean gates and used an RMSSD drop as +// the PRIMARY REM detector. Measured against 99 PSG-labelled wrist nights +// (DREAMT) via `tool/stager_harness.dart`, that design scored kappa 0.036 — +// deep PPV 5.7% against a 4.5% base rate, REM PPV 12.1% against 14.0%, i.e. at +// or below chance for both minority classes. Per-subject effect sizes explain +// why: RMSSD does not separate the stages at all (54% of subjects, d -0.02; +// Herzig 2017 reports the same flatness on PSG+ECG), mean HR in deep sleep runs +// HIGHER than light on the wrist (62% of subjects, d +0.31 — the old gate was +// inverted), and the two axes that DO separate (Rk d -0.53/+0.43, sdnn +// -0.41/+0.32) were respectively throttled to z>4 and not computed at all. +// AND-ing one good gate with one null and one inverted gate is why deep sleep +// emerged as isolated 30-second specks that the 3-min bout rule then deleted. +// A weighted sum degrades gracefully where a conjunction vetoes. // Post: Webster continuity rescore (bridges brief arousals into sleep — this is // what kills the over-call) + consolidateSleepStages (no single-epoch flicker). // @@ -59,20 +76,82 @@ class CardioStagerResult { const int _epochSec = 30; -/// Robust-z cutoffs for the SECONDARY REM axes (LF/HF elevation, R(k) burst). -/// Deliberately STRICT (4.0 robust-SD): the RMSSD drop is the PRIMARY REM -/// detector; LF/HF and R(k) only ADD REM where a strong, unambiguous autonomic -/// shift the RMSSD rule missed exists, so the OR-combination recovers -/// under-called REM without stealing normal light sleep. Calibrated on the real -/// 2026-07 overnight fixture (Apple-Watch GT: wake 3 / light 330 / deep 38 / -/// REM 162 min). RMSSD-only under-called REM at 134 min; a threshold sweep -/// showed z=2.5→204, 3.0→191, 3.5→182, 4.0→172 min REM. 4.0 lands closest to GT -/// (REM 172, light 304) while still adding +38 min of strong-signature REM over -/// RMSSD-only — a conservative, bounded sensitivity gain for the nights the -/// RMSSD rule under-calls (the motivating case: 64 min REM vs a ~115 signature), -/// erring toward specificity so the axes never re-open a Wake/REM over-call. -const double _remLfhfZ = 4.0; -const double _remRkZ = 4.0; +// The previous rule OR-ed three boolean axes (RMSSD drop as PRIMARY, plus +// LF/HF and R(k) throttled to z>4.0) and was calibrated against a single +// Apple-Watch-labelled night. Measurement against 99 PSG-labelled wrist nights +// showed the primary axis was at chance and the two throttled axes were the +// informative ones; see the weighted-score block in [classifyCardioEpochs]. + +// ── Stage-score axis weights ──────────────────────────────────────────────── +// These are MEASURED Cohen's d values (DREAMT, 99 PSG-labelled wrist nights, +// per-subject within-night z-scored), not tuned parameters. Under equal +// variance an LD's optimal weights are proportional to d, so the effect size +// IS the weight. Only the two decision cutoffs below are calibrated, on a +// locked development split. +const double _wRemRk = 0.43; +const double _wRemSdnn = 0.32; +const double _wRemHrSd = 0.27; +const double _wRemLfhf = 0.12; + +const double _wDeepRk = 0.53; +const double _wDeepHrSd = 0.43; +const double _wDeepSdnn = 0.41; +const double _wDeepLfhf = 0.33; + +/// Decision cutoffs on the weighted robust-z stage scores. +/// +/// The WEIGHTS above come from DREAMT (PSG ground truth) because that is what +/// an external corpus can tell us: which axes carry stage information, and in +/// which direction. The CUTOFFS cannot come from there. They sit on a score +/// whose spread depends on OUR feature extraction, and the DREAMT-optimal +/// values do not transfer: at remCut 0.8 DREAMT yields a normative ~23% REM +/// while real WHOOP captures yield 10.3%, and the same night the Apple-Watch +/// reference put at ~162 min REM came out at 62 min. Calibrating cutoffs on a +/// different sensor's feature distribution is exactly the train/serve mismatch +/// this file's header criticises Walch for. +/// +/// So the cutoffs are swept on REAL DEVICE CAPTURES +/// (`tool/whoop_proportions.dart`, 11 nights) and chosen to land the stage +/// proportions in the normative population range — REM 20-25% of TST, Deep +/// 13-23% of TST (Mitterling 2015 percentile curves; Boulos 2019 +/// meta-regression, N=5,273, AASM scoring): +/// +/// remCut median %REM of TST deepCut median %Deep of TST +/// 0.4 28.5% 0.2 17.6% +/// 0.5 23.6% <- chosen 0.3 14.9% <- chosen +/// 0.6 18.6% 0.4 11.9% +/// 0.8 ~10% (DREAMT opt) +/// +/// NOT calibrated to maximise kappa, on purpose. On the DREAMT dev split kappa +/// rises monotonically as deep is suppressed: +/// +/// deepCut %Deep called Deep sens Deep PPV kappa +/// 0.3 14.9% 50.8% 15.0% 0.125 +/// 0.6 4.8% 22.2% 20.4% 0.148 +/// 1.0 0.3% 2.3% 31.0% 0.151 <- kappa optimum +/// +/// That is arithmetic on a sleep-clinic cohort where deep is only 4.5% of +/// asleep epochs, not physiology, and shipping the kappa optimum would recreate +/// the "user sees ~0 deep sleep" defect this change exists to fix. +/// +/// Held-out result at the shipped cutoffs (49 unseen subjects): kappa 0.132, +/// Deep 56.0% sens / 10.9% PPV, REM 51.9% / 21.5%, against base rates of 4.5% +/// Deep and 14% REM — so both minority classes are called well above chance +/// where the previous rules sat at or below it. Still far short of the +/// 0.60-0.66 literature ceiling; the bulk of what remains is the WAKE rule +/// (11% sensitivity), untouched by this change. +/// +/// HONESTY: proportion-matching is distributional plausibility, NOT accuracy. +/// It says the hypnogram has a believable shape, not that any given epoch is +/// right. Only concurrent PSG on this device can settle that, and we do not +/// have it. +const double _remScoreCut = 0.5; +const double _deepScoreCut = 0.3; + +/// The shipped cutoffs, exposed so calibration tooling reports what actually +/// ships instead of re-declaring its own defaults and silently drifting. +const double kDefaultRemScoreCut = _remScoreCut; +const double kDefaultDeepScoreCut = _deepScoreCut; /// Physiologic RR gate (project rule): keep 300–2000 ms; drop successive jumps /// > 200 ms (ectopy / artifact). Used per-window before RMSSD. @@ -257,6 +336,8 @@ CardioStagerResult cardioStager( List rrTsMs = const [], int epochSec = _epochSec, SleepUserProfile? userProfile, + double remScoreCut = _remScoreCut, + double deepScoreCut = _deepScoreCut, }) { // Explicit arg wins (unit tests); else the ambient profile the edge set for // this staging pass; else null ⇒ pure per-night-local (cold-start behavior). @@ -303,6 +384,7 @@ CardioStagerResult cardioStager( final hr = List.filled(nEpoch, double.nan); final hrSd = List.filled(nEpoch, 0); final rmssd = List.filled(nEpoch, double.nan); + final sdnn = List.filled(nEpoch, double.nan); // REM autonomic features (P1): LF/HF from the RR spectrum + R(k) = rolling // mean |ΔIHR|. REM's signature at 1 Hz beat-timing is a SYMPATHETIC shift — // LF/HF rises, RR variability drops (low RMSSD), and instantaneous HR gets @@ -328,12 +410,130 @@ CardioStagerResult cardioStager( if (hv.length >= 2) hrSd[e] = stddev(hv) ?? 0; // rmssd over RR beats within a ±2.5-min window centred on the epoch rmssd[e] = _windowRmssd(rrMs, rrTsMs, accel, s, t, epochSec); + sdnn[e] = _windowSdnn(rrMs, rrTsMs, accel, s, t, epochSec); // LF/HF + R(k) over a ±90-s RR window centred on the epoch. final rem = _windowRemFeatures(rrMs, rrTsMs, accel, s, t, epochSec); if (rem.lfhf != null) lfhf[e] = rem.lfhf!; if (rem.rk != null) rk[e] = rem.rk!; } + return classifyCardioEpochs( + CardioEpochFeatures( + motion: motion, + hr: hr, + hrSd: hrSd, + rmssd: rmssd, + lfhf: lfhf, + rk: rk, + sdnn: sdnn, + ), + epochSec: epochSec, + userProfile: profile, + remScoreCut: remScoreCut, + deepScoreCut: deepScoreCut, + ); +} + +/// Per-epoch feature vectors — the ONLY thing the staging DECISION layer sees. +/// +/// Splitting these out from [cardioStager] draws the line between "turn 1 Hz +/// sensor streams into per-epoch numbers" (device-specific, hard to validate) +/// and "turn per-epoch numbers into stages" (the rules, which are what actually +/// encode our physiological assumptions). Only the second half can be scored +/// against an external PSG-labelled corpus whose raw signals we do not have, so +/// it needs to be callable on its own. All lists must be the same length; NaN +/// marks "not measurable for this epoch" and is never treated as a value. +class CardioEpochFeatures { + /// Mean positive ENMO over the epoch against a LOCAL 1 g reference (g). + final List motion; + + /// Mean valid HR over the epoch (bpm); NaN when the epoch had no valid HR. + final List hr; + + /// Within-epoch SD of HR (bpm); 0 when fewer than 2 valid samples. + final List hrSd; + + /// RMSSD over a ±2.5-min window centred on the epoch (ms); NaN when sparse. + final List rmssd; + + /// LF/HF band-power ratio over a ±90 s RR window; NaN when sparse. + final List lfhf; + + /// R(k) = mean |ΔIHR| over a ±90 s RR window (bpm); NaN when sparse. + final List rk; + + /// SDNN over the same ±2.5-min window as [rmssd] (ms); NaN when sparse. + /// Unlike RMSSD, SDNN DOES separate the stages — see the effect sizes in + /// [classifyCardioEpochs]. + final List sdnn; + + const CardioEpochFeatures({ + required this.motion, + required this.hr, + required this.hrSd, + required this.rmssd, + required this.lfhf, + required this.rk, + required this.sdnn, + }); + + /// Epoch count. All lists are required to be this long; the getter returns + /// the SHORTEST so a caller that violates the contract truncates instead of + /// indexing out of range. [classifyCardioEpochs] asserts the equality. + int get length => [ + motion.length, + hr.length, + hrSd.length, + rmssd.length, + lfhf.length, + rk.length, + sdnn.length, + ].reduce(math.min); +} + +/// The staging DECISION layer: per-epoch features → stages + deep flags. +/// +/// Extracted from [cardioStager] verbatim (no behavioural change) so the rules +/// can be scored directly against an external labelled corpus. See +/// [CardioEpochFeatures] for why the seam is here. +/// +/// [remScoreCut] / [deepScoreCut] default to the shipped constants and exist so +/// a calibration harness can sweep them on a DEVELOPMENT split and report the +/// held-out result. Production callers should not pass them. +CardioStagerResult classifyCardioEpochs( + CardioEpochFeatures f, { + int epochSec = _epochSec, + SleepUserProfile? userProfile, + double remScoreCut = _remScoreCut, + double deepScoreCut = _deepScoreCut, +}) { + final profile = userProfile ?? cardioUserProfile; + // [CardioEpochFeatures] documents that every list is the same length, and the + // in-tree caller guarantees it — but the type is public and the validation + // harness builds one from parsed external fixture data. A short secondary + // list would otherwise RangeError deep inside the classify loop, so assert in + // dev and, in release, fall back to the shortest list rather than throwing + // (the project's existing idiom: `cardioStager` already clamps to + // `min(hr1hz.length, accel.length)`). + assert( + f.hr.length == f.motion.length && + f.hrSd.length == f.motion.length && + f.rmssd.length == f.motion.length && + f.lfhf.length == f.motion.length && + f.rk.length == f.motion.length && + f.sdnn.length == f.motion.length, + 'CardioEpochFeatures lists must all be the same length', + ); + final nEpoch = f.length; + if (nEpoch < 3) return _abstain(epochSec); + final motion = f.motion; + final hr = f.hr; + final hrSd = f.hrSd; + final rmssd = f.rmssd; + final lfhf = f.lfhf; + final rk = f.rk; + final sdnn = f.sdnn; + // ── night baselines (from the LOW-MOTION epochs — the actual sleep) ──────── // motMed/motMad stay WHOLE-NIGHT scalars: `motion` is now computed against a // per-epoch LOCAL reference (above), so genuine stillness reads ~0 almost @@ -375,7 +575,6 @@ CardioStagerResult cardioStager( final hrAll = [for (final h in hr) if (!h.isNaN) h]; if (hrAll.isEmpty) return _abstain(epochSec); final hrMedGlobal = median(sleepHr) ?? mean(hrAll)!; - final hrFloor = percentile(sleepHr, 10) ?? hrMedGlobal; // ── LOCAL rolling HR baseline for the WAKE/REM autonomic gates ───────────── // A whole-night HR median/arousal threshold has the same non-stationarity @@ -444,20 +643,31 @@ CardioStagerResult cardioStager( // These three samples are fixed for the rest of the stage, so take their // median+MAD once instead of re-deriving them per epoch inside the classify // loop below. Same values, ~2,880 fewer sorts on a full night. - final rmssdScale = RobustScale.of(sleepRmssd); final lfhfScale = RobustScale.of(sleepLfhf); final rkScale = RobustScale.of(sleepRk); - final hrSdSample = [for (final s in hrSd) if (s > 0) s]; - final hrSdMed = median(hrSdSample) ?? double.infinity; - // RMSSD reference for the deep "not-elevated-HRV" gate, blended toward profile. - final rmssdRef = rmssdMed == null - ? profile?.rmssdMed - : blendP(rmssdMed, profile?.rmssdMed); - // Deep = HR in the lower half of the night's sleeping HR (the cardiac trough), - // with the floor/median blended toward the personal profile (P2). - final deepFloor = blendP(hrFloor, profile?.hrFloorP5); - final deepMed = blendP(hrMedGlobal, profile?.hrSleepMedian); - final deepHrCut = deepFloor + 0.5 * (deepMed - deepFloor); + // SDNN and within-epoch HR dispersion get the same night-relative treatment + // as the other axes, so every term in the stage scores is on one scale. + final sleepSdnn = [ + for (var e = 0; e < nEpoch; e++) + if (still(e) && !sdnn[e].isNaN) sdnn[e] + ]; + final sdnnScale = RobustScale.of(sleepSdnn); + // hrSd gets the SAME still(e) gating as the other three axes. Pooling wake + // epochs into this baseline raises its median and widens its MAD, which + // biases every hrSdZ low — and hrSd carries the second-largest deep weight + // (0.43), so that alone shifts stage proportions and invalidates the cutoffs. + // + // The `> 0` filter is not cosmetic: hrSd is 0-initialised and only assigned + // when an epoch had >= 2 valid HR samples, so 0 is the "no HR here" sentinel, + // NOT a real measurement of zero dispersion. Scoring it would make a + // data-gap epoch read as maximally quiet, i.e. as evidence FOR deep sleep. + // Measured on real captures: one night with intermittent wear had 226 of + // 1120 epochs (20%) with no HR, every one of which scored toward deep. + final sleepHrSd = [ + for (var e = 0; e < nEpoch; e++) + if (still(e) && hrSd[e] > 0) hrSd[e] + ]; + final hrSdScale = RobustScale.of(sleepHrSd); // ── classify ─────────────────────────────────────────────────────────────── final stages = List.filled(nEpoch, SleepStage.wake); @@ -476,42 +686,92 @@ CardioStagerResult cardioStager( stages[e] = SleepStage.wake; continue; } - // Asleep. REM vs NREM via autonomic signature. REM is OR-combined across - // three independent RR axes (any one suffices), THEN gated by atonia (no - // large movement) AND an HR floor (HR ≥ the local p25 — REM is not the - // quiescent cardiac trough). This recovers REM the RMSSD-only rule missed. - final rmZ = (rmssdMed != null && !rmssd[e].isNaN && sleepRmssd.length >= 4) - ? rmssdScale?.z(rmssd[e]) - : null; - final rmssdDown = rmZ != null && rmZ < -0.4; // RMSSD notably below sleep base - // LF/HF elevated (sympathetic shift) vs the night's sleeping LF/HF. - final lfhfZ = (sleepLfhf.length >= 4 && !lfhf[e].isNaN) - ? lfhfScale?.z(lfhf[e]) - : null; - final lfhfHigh = lfhfZ != null && lfhfZ > _remLfhfZ; - // R(k) burst: instantaneous-HR variability elevated vs sleeping R(k). - final rkZ = (sleepRk.length >= 4 && !rk[e].isNaN) - ? rkScale?.z(rk[e]) - : null; - final rkBurst = rkZ != null && rkZ > _remRkZ; - final remAutonomic = rmssdDown || lfhfHigh || rkBurst; - final atonia = !bigMove(e); // muscle-atonia proxy — no large movement - final hrTowardWake = - !hr[e].isNaN && hr[e] >= hrP25Local[e]; // HR up but not arousal - if (remAutonomic && atonia && hrTowardWake) { + // Asleep. REM vs NREM, and deep-within-NREM, are now scored as WEIGHTED + // SUMS of robust-z axes rather than boolean conjunctions. + // + // WHY THE CHANGE. The old rules were an AND of three booleans each, and + // measurement on 99 PSG-labelled wrist nights (DREAMT) showed two of the + // three deep gates carry no usable signal — one of them points the wrong + // way — while the REM primary is at chance. Per-subject agreement and + // Cohen's d, Deep-vs-Light and REM-vs-Light: + // + // axis Deep vs Light REM vs Light used before? + // Rk 91% subj, d -0.53 89% subj, d +0.43 throttled to z>4 + // hrSd 97% subj, d -0.43 d +0.27 deep only + // sdnn 74% subj, d -0.41 d +0.32 NOT COMPUTED + // lfhf 91% subj, d -0.33 72% subj, d +0.12 throttled to z>4 + // hr 62% subj, d +0.31 (!) d +0.18 deep gate, BACKWARDS + // rmssd d -0.13 54% subj, d -0.02 REM PRIMARY (chance) + // + // Scored as conjunctions those numbers are fatal: the shipped deep rule + // ANDs one good axis, one null axis and one inverted axis, so the three can + // only co-fire by coincidence — which is exactly why deep sleep came out as + // isolated 30-second specks that the 3-min minimum-bout rule then deleted. + // A weighted sum degrades gracefully instead: a missing or equivocal axis + // costs its weight rather than vetoing the whole decision. + // + // WEIGHTS ARE THE MEASURED EFFECT SIZES, not tuned. Under equal variance a + // linear discriminant's optimal weights are proportional to d, so using |d| + // directly is the principled choice AND leaves nothing fitted to the corpus + // beyond its sign. This also matches how the published cardiorespiratory + // stagers combine features (Long 2016, Fonseca 2015/2018 all use linear + // discriminants; Fonseca 2018 measured a generative HMM at kappa 0.26 vs + // 0.39 for a plain per-epoch LD on identical features, so the sophistication + // does not belong in the decoder). + // + // rmssd and mean HR are deliberately ABSENT from both scores. Neither + // earns a weight, and hr's deep contribution was actively harmful. + final rkZ = (sleepRk.length >= 4 && !rk[e].isNaN) ? rkScale?.z(rk[e]) : null; + final sdnnZ = + (sleepSdnn.length >= 4 && !sdnn[e].isNaN) ? sdnnScale?.z(sdnn[e]) : null; + final lfhfZ = + (sleepLfhf.length >= 4 && !lfhf[e].isNaN) ? lfhfScale?.z(lfhf[e]) : null; + // Same >=4-sample floor and unmeasurable-epoch rejection as the axes above. + // hrSd[e] == 0 means "fewer than 2 valid HR samples", never "perfectly + // steady" — see the sleepHrSd comment. An absent input must contribute + // nothing, not a favourable z. + final hrSdZ = + (sleepHrSd.length >= 4 && hrSd[e] > 0) ? hrSdScale?.z(hrSd[e]) : null; + + // Weighted mean over the axes that are actually MEASURABLE this epoch, so + // a night with no usable RR falls back to the hrSd axis alone instead of + // silently scoring every epoch as if the missing axes voted "no". + double? score(List<(double?, double)> terms) { + var num = 0.0, den = 0.0; + for (final (z, w) in terms) { + if (z == null || z.isNaN) continue; + num += z * w; + den += w.abs(); + } + return den == 0 ? null : num / den; + } + + final remScore = score([ + (rkZ, _wRemRk), + (sdnnZ, _wRemSdnn), + (hrSdZ, _wRemHrSd), + (lfhfZ, _wRemLfhf), + ]); + final deepScore = score([ + (rkZ, -_wDeepRk), + (hrSdZ, -_wDeepHrSd), + (sdnnZ, -_wDeepSdnn), + (lfhfZ, -_wDeepLfhf), + ]); + + // Atonia and the HR floor are RETAINED as gates, not folded into the score: + // they are not weak evidence for REM, they are preconditions. A large + // movement rules REM out outright (muscle atonia), and REM is not the + // quiescent cardiac trough, so HR below the local p25 rules it out too. + final atonia = !bigMove(e); + final hrTowardWake = !hr[e].isNaN && hr[e] >= hrP25Local[e]; + if (remScore != null && remScore > remScoreCut && atonia && hrTowardWake) { stages[e] = SleepStage.rem; } else { stages[e] = SleepStage.nrem; - // Deep (NREM subtype, LOW CONFIDENCE): the cardiac trough — HR in the - // lower half of the night's sleeping HR AND not HR-variable (deep sleep is - // autonomically quiet). RMSSD, when present, reinforces (high RMSSD) but - // isn't required (RR is sparse). Below NREM median, not the lowest third, - // so deep lands in a physiologic range instead of ~0. - final lowHr = !hr[e].isNaN && hr[e] <= deepHrCut; - final notHighRmssd = - rmssdRef == null || rmssd[e].isNaN || rmssd[e] >= rmssdRef * 0.9; - final stable = hrSd[e] <= hrSdMed * 1.5; - deepFlag[e] = lowHr && notHighRmssd && stable; + // Deep (NREM subtype, still LOW CONFIDENCE — this is a cardiac-quiescence + // overlay, not an EEG slow-wave measurement). + deepFlag[e] = deepScore != null && deepScore > deepScoreCut; } } @@ -637,6 +897,59 @@ double _windowRmssd(List rrMs, List rrTsMs, return math.sqrt(ss / (beats.length - 1)); } +/// SDNN (ms) of cleaned RR beats over the SAME ±2.5-min window as +/// [_windowRmssd]. Returns NaN when too few clean beats. +/// +/// Worth its own feature because SDNN and RMSSD behave completely differently +/// across sleep stages, despite both being "HRV". Measured per-subject on 99 +/// PSG-labelled wrist nights (DREAMT): SDNN separates Deep from Light in 74% of +/// subjects (Cohen's d −0.41) and REM from Light at d +0.32, while RMSSD is at +/// chance for both (54% of subjects, d −0.02). Herzig 2017 (PSG+ECG) reports +/// the same split — SDNN 53.8 / 68.5 / 105.5 ms for N3 / N2 / REM against a +/// 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); + if (beats.length < 5) return double.nan; + final m = mean(beats)!; + var ss = 0.0; + for (final v in beats) { + ss += (v - m) * (v - m); + } + 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 []; + final mid = (s + t) ~/ 2; + if (mid >= accel.length) return const []; + final centreMs = accel[mid].tsMs; + const halfWinMs = 150 * 1000; + final lo = centreMs - halfWinMs, hi = centreMs + halfWinMs; + 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; + } + return 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. diff --git a/tool/stager_harness.dart b/tool/stager_harness.dart new file mode 100644 index 0000000..4d69241 --- /dev/null +++ b/tool/stager_harness.dart @@ -0,0 +1,499 @@ +// VALIDATION HARNESS — score the SHIPPED staging decision layer against an +// external PSG-labelled corpus. +// +// Runs `classifyCardioEpochs` (the real production rules, not a +// reimplementation) over per-epoch features exported from a labelled dataset, +// and reports the metrics the literature says you must report: +// +// * Cohen's kappa, not accuracy. Canton 2026 (Sleep Adv, PMID 42333378) +// showed a trivial "easy-to-classify wake" baseline explains much of the +// reported performance of published models, and that wake-rich datasets +// inflate accuracy. A majority-class baseline is printed alongside so the +// number can be read honestly. +// * PER-CLASS sensitivity AND PPV. Fonseca 2018 (PMID 29620019) found a +// temporal model traded 8 pp of N3 recall for 13 pp of N3 precision — a +// single figure hides exactly the tradeoff we care about. +// * The PER-SUBJECT kappa distribution, not just the pooled value. Radha +// 2019 (PMID 31578345) reports kappa 0.61 +- 0.15 over 292 subjects; the +// SD is the whole story for "works for most, awful for a few", and the +// count below kappa 0.2 is the population we are trying to serve better. +// +// CAVEAT, and it is a real one: the fixture supplies the dataset's OWN feature +// extraction (window lengths, artifact handling, sensor), not ours. What this +// scores is the DECISION LAYER. A change that improves the rules here should +// improve them on our device only insofar as our features carry the same +// information — which is itself an assumption worth stating out loud. +// +// Usage: +// dart run tool/stager_harness.dart [flags] +// +// --sleep-window trim each subject to [first non-Wake, last non-Wake] +// --dev | --holdout score only one half of the locked subject split +// --sweep sweep both cutoffs on DEV, report HOLDOUT at the optimum +// --deep-curve deep cutoff vs call-rate / sens / PPV (DEV only) +// --rem-cut override the REM cutoff (default: the shipped value) +// --deep-cut override the deep cutoff (default: the shipped value) +// +// Fixture schema — every feature list must be the same length as `stages`, and +// `null` means "not measurable for this epoch" (never 0): +// {"epochSec": 30, +// "subjects": [{"subject": "S1", "stages": ["Wake"|"Light"|"Deep"|"REM", ..], +// "motion": [..], "hr": [..], "hrSd": [..], +// "rmssd": [..], "lfhf": [..], "rk": [..], "sdnn": [..]}]} + +import 'dart:convert'; +import 'dart:io'; +import 'dart:math' as math; + +import 'package:openstrap_analytics/onehz.dart'; + +/// 4-class labels, in the order used for the confusion matrix. +const kClasses = ['Wake', 'Light', 'Deep', 'REM']; + +/// Feature lists a subject record must carry, all aligned with `stages`. +const _kFeatureKeys = ['motion', 'hr', 'hrSd', 'rmssd', 'lfhf', 'rk', 'sdnn']; + +void main(List args) { + if (args.isEmpty) { + stderr.writeln('usage: dart run tool/stager_harness.dart ' + '[--sleep-window] [--dev|--holdout] [--sweep] [--deep-curve] ' + '[--rem-cut x] [--deep-cut x]'); + exitCode = 64; + return; + } + final fixture = + jsonDecode(File(args.first).readAsStringSync()) as Map; + final epochSec = (fixture['epochSec'] as num?)?.toInt() ?? 30; + var subjects = (fixture['subjects'] as List).cast>(); + + final schemaErrors = _validate(subjects); + if (schemaErrors.isNotEmpty) { + for (final e in schemaErrors.take(10)) { + stderr.writeln('fixture: $e'); + } + if (schemaErrors.length > 10) { + stderr.writeln('fixture: ... and ${schemaErrors.length - 10} more'); + } + exitCode = 65; + return; + } + + // Defaults come FROM THE LIBRARY. Re-declaring them here is how a harness + // silently stops measuring what actually ships. + double argOf(String name, double dflt) { + final i = args.indexOf(name); + if (i < 0 || i + 1 >= args.length) return dflt; + final v = double.tryParse(args[i + 1]); + if (v == null) { + stderr.writeln('warning: could not parse $name "${args[i + 1]}", ' + 'using $dflt'); + return dflt; + } + return v; + } + + final remCut = argOf('--rem-cut', kDefaultRemScoreCut); + final deepCut = argOf('--deep-cut', kDefaultDeepScoreCut); + + // In production `cardioStager` never sees a whole recording: it is handed a + // window that van Hees + AdvancedSleepStager.detectSleep have already picked + // out, so its baselines (`stillCut`, the local sleeping-HR median, the RMSSD + // reference) are estimated over mostly-sleep epochs. Scoring it over a full + // record including a long pre-sleep wake block is a domain mismatch that + // blames the rules for a windowing job they do not do. + if (args.contains('--sleep-window')) { + subjects = [ + for (final s in subjects) + if (_trimToSleepPeriod(s) case final t?) t + ]; + stdout.writeln('[--sleep-window] trimmed each subject to ' + '[first sleep, last sleep] — ${subjects.length} subjects'); + } + + // CALIBRATION PROTOCOL. Any cutoff chosen by looking at the data must be + // reported on subjects that were not looked at, so every mode that sweeps a + // cutoff runs on DEV only. + if (args.contains('--deep-curve')) { + final dev = _ofSplit(subjects, _Split.dev); + if (dev.isEmpty) { + stderr.writeln('--deep-curve needs a non-empty dev split'); + exitCode = 1; + return; + } + stdout.writeln('[--deep-curve] DEV only (${dev.length} subjects), ' + 'remCut=${remCut.toStringAsFixed(2)}'); + stdout.writeln('deepCut %Deep called Deep sens Deep PPV kappa'); + for (final d in [0.0, 0.2, 0.3, 0.4, 0.6, 0.8, 1.0]) { + final m = _deepMetrics(dev, epochSec, remCut, d); + stdout.writeln(' ${d.toStringAsFixed(1)} ' + '${(100 * m.$1).toStringAsFixed(1).padLeft(5)}% ' + '${(100 * m.$2).toStringAsFixed(1).padLeft(5)}% ' + '${(100 * m.$3).toStringAsFixed(1).padLeft(5)}% ' + '${_fmtKappa(m.$4)}'); + } + return; + } + if (args.contains('--sweep')) { + _sweep(subjects, epochSec); + return; + } + + final split = args.contains('--dev') + ? _Split.dev + : args.contains('--holdout') + ? _Split.holdout + : _Split.all; + if (split != _Split.all) { + subjects = _ofSplit(subjects, split); + stdout.writeln('[${split.name}] ${subjects.length} subjects'); + } + final shipped = + remCut == kDefaultRemScoreCut && deepCut == kDefaultDeepScoreCut; + stdout.writeln('cutoffs: rem=${remCut.toStringAsFixed(2)} ' + 'deep=${deepCut.toStringAsFixed(2)}' + '${shipped ? " (shipped)" : " (OVERRIDDEN — not what ships)"}'); + + final truth = []; + final pred = []; + final perSubjectKappa = []; + var skipped = 0; + + for (final s in subjects) { + final scored = _scoreSubject(s, epochSec, remCut, deepCut); + if (scored == null || scored.$1.length < 30) { + skipped++; + continue; + } + truth.addAll(scored.$1); + pred.addAll(scored.$2); + final k = _kappa(scored.$1, scored.$2); + if (!k.isNaN) perSubjectKappa.add(k); + } + + if (truth.isEmpty || perSubjectKappa.isEmpty) { + stderr.writeln('no subject produced a scoreable hypnogram ' + '($skipped skipped: abstained, too short, or unlabelled)'); + exitCode = 1; + return; + } + + stdout.writeln('subjects scored ${perSubjectKappa.length}' + '${skipped > 0 ? ' (skipped $skipped: abstained/too short)' : ''}' + ' epochs ${truth.length}'); + stdout.writeln(''); + + final counts = List.filled(4, 0); + for (final t in truth) { + counts[t]++; + } + final majority = counts.indexOf(counts.reduce(math.max)); + _report('BASELINE all-${kClasses[majority]}', truth, + List.filled(truth.length, majority)); + _report('SHIPPED classifyCardioEpochs', truth, pred); + + perSubjectKappa.sort(); + double q(double f) => + perSubjectKappa[(f * (perSubjectKappa.length - 1)).round()]; + final mean = + perSubjectKappa.reduce((a, b) => a + b) / perSubjectKappa.length; + var ss = 0.0; + for (final k in perSubjectKappa) { + ss += (k - mean) * (k - mean); + } + final sd = math.sqrt(ss / perSubjectKappa.length); + final bad = perSubjectKappa.where((k) => k < 0.2).length; + stdout.writeln(''); + stdout.writeln('per-subject kappa mean ${mean.toStringAsFixed(3)} ' + 'SD ${sd.toStringAsFixed(3)} ' + 'p10 ${q(0.10).toStringAsFixed(3)} ' + 'median ${q(0.50).toStringAsFixed(3)} ' + 'p90 ${q(0.90).toStringAsFixed(3)}'); + stdout.writeln('THE BAD TAIL kappa < 0.2 for $bad' + ' / ${perSubjectKappa.length} subjects ' + '(${(100 * bad / perSubjectKappa.length).toStringAsFixed(0)}%)'); + stdout.writeln(''); + stdout.writeln('reference: 4-class cardiorespiratory kappa 0.60-0.66 ' + '(Radha 2019 0.61+-0.15; Bakker 2021 0.643; Sridhar 2020 0.66)'); +} + +/// Structural check on the fixture. A malformed record must fail loudly here, +/// not produce a confident wrong kappa. +List _validate(List> subjects) { + final errs = []; + for (var i = 0; i < subjects.length; i++) { + final s = subjects[i]; + final id = s['subject'] ?? '#$i'; + if (s['subject'] is! String) errs.add('$id: missing string "subject"'); + final stages = s['stages']; + if (stages is! List) { + errs.add('$id: missing "stages" list'); + continue; + } + for (final k in _kFeatureKeys) { + final v = s[k]; + if (v is! List) { + errs.add('$id: missing "$k" list'); + } else if (v.length != stages.length) { + errs.add( + '$id: "$k" has ${v.length} entries, "stages" has ${stages.length}'); + } + } + } + return errs; +} + +enum _Split { dev, holdout, all } + +/// Deterministic, order-independent 50/50 partition by subject id (FNV-1a). +/// +/// Splits on a MIDDLE bit, not the parity of the final hash. FNV-1a's last step +/// multiplies by an odd constant, so the low bit of the result is just the low +/// bit of `prevHash ^ lastByte` — testing `isEven` would degenerate to +/// "alternate by the last character", which for ids like S001/S002/S003 gives a +/// striped split rather than a hashed one. +_Split _splitOf(String subject) { + var h = 0x811c9dc5; + for (final c in subject.codeUnits) { + h = ((h ^ c) * 0x01000193) & 0xFFFFFFFF; + } + return ((h >> 16) & 1) == 0 ? _Split.dev : _Split.holdout; +} + +List> _ofSplit( + List> subjects, _Split want) => + [ + for (final s in subjects) + if (_splitOf(s['subject'] as String) == want) s + ]; + +CardioEpochFeatures _featuresOf(Map s) => CardioEpochFeatures( + motion: _nums(s['motion'], 0), + hr: _nums(s['hr'], double.nan), + hrSd: _nums(s['hrSd'], 0), + rmssd: _nums(s['rmssd'], double.nan), + lfhf: _nums(s['lfhf'], double.nan), + rk: _nums(s['rk'], double.nan), + sdnn: _nums(s['sdnn'], double.nan), + ); + +/// (truth, prediction) label indices for one subject, or null when the stager +/// abstained. The single place reference labels are mapped onto our 3-class + +/// deep-flag output — this loop previously existed in three scoring paths, +/// which is exactly how they drift apart. +(List, List)? _scoreSubject( + Map s, int epochSec, double remCut, double deepCut) { + final stages = (s['stages'] as List).cast(); + final res = classifyCardioEpochs(_featuresOf(s), + epochSec: epochSec, remScoreCut: remCut, deepScoreCut: deepCut); + final out = res.base.stages; + if (out.isEmpty) return null; // honest abstain + final t = [], p = []; + final unknown = {}; + for (var e = 0; e < out.length && e < stages.length; e++) { + final ti = kClasses.indexOf(stages[e]); + if (ti < 0) { + unknown.add(stages[e]); + continue; + } + t.add(ti); + p.add(switch (out[e]) { + SleepStage.wake => 0, + SleepStage.rem => 3, + SleepStage.nrem => (e < res.deepFlag.length && res.deepFlag[e]) ? 2 : 1, + }); + } + if (unknown.isNotEmpty) { + stderr.writeln('warning: ${s['subject']}: dropped unrecognised stage ' + 'labels ${unknown.toList()..sort()}'); + } + return (t, p); +} + +/// Sweep both cutoffs on DEV only, then print the HOLDOUT score at the dev +/// optimum. The holdout column is the only one that means anything. +void _sweep(List> subjects, int epochSec) { + final dev = _ofSplit(subjects, _Split.dev); + final hold = _ofSplit(subjects, _Split.holdout); + stdout.writeln('dev ${dev.length} subjects | holdout ${hold.length} subjects'); + if (dev.isEmpty || hold.isEmpty) { + stderr.writeln('sweep needs a non-empty dev AND holdout split'); + exitCode = 1; + return; + } + stdout.writeln(''); + final cuts = [0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.3, 1.6, 2.0]; + var best = double.negativeInfinity; + var bestR = kDefaultRemScoreCut, bestD = kDefaultDeepScoreCut; + for (final r in cuts) { + for (final d in cuts) { + final k = _scoreSet(dev, epochSec, r, d); + if (!k.isNaN && k > best) { + best = k; + bestR = r; + bestD = d; + } + } + } + if (best == double.negativeInfinity) { + stderr.writeln('sweep produced no scoreable configuration'); + exitCode = 1; + return; + } + stdout.writeln(' remCut deepCut DEV kappa HOLDOUT kappa ' + 'holdout %REM %Deep'); + for (final r in cuts) { + final k = _scoreSet(dev, epochSec, r, bestD); + final kh = _scoreSet(hold, epochSec, r, bestD); + final rates = _callRates(hold, epochSec, r, bestD); + final mark = r == bestR ? ' <- dev optimum' : ''; + stdout.writeln(' ${r.toStringAsFixed(1)} ${bestD.toStringAsFixed(1)}' + ' ${_fmtKappa(k)} ${_fmtKappa(kh)}' + ' ${(100 * rates.$1).toStringAsFixed(1).padLeft(5)}% ' + '${(100 * rates.$2).toStringAsFixed(1).padLeft(6)}%$mark'); + } + stdout.writeln(''); + stdout.writeln('dev optimum remCut=${bestR.toStringAsFixed(1)} ' + 'deepCut=${bestD.toStringAsFixed(1)} DEV ${_fmtKappa(best)} ' + 'HOLDOUT ${_fmtKappa(_scoreSet(hold, epochSec, bestR, bestD))}'); + stdout.writeln('shipped ' + 'remCut=${kDefaultRemScoreCut.toStringAsFixed(1)} ' + 'deepCut=${kDefaultDeepScoreCut.toStringAsFixed(1)} ' + 'DEV ${_fmtKappa(_scoreSet(dev, epochSec, kDefaultRemScoreCut, kDefaultDeepScoreCut))} ' + 'HOLDOUT ${_fmtKappa(_scoreSet(hold, epochSec, kDefaultRemScoreCut, kDefaultDeepScoreCut))}'); + stdout.writeln(''); + stdout.writeln('NOTE: the dev optimum is NOT automatically the right choice. ' + 'See the cutoff doc comment in cardio_stager.dart for why kappa on a ' + 'low-deep cohort rewards suppressing deep sleep.'); +} + +String _fmtKappa(double k) => k.isNaN ? ' n/a' : k.toStringAsFixed(3); + +/// Fraction of epochs predicted REM and Deep — guards against a cutoff that +/// "wins" kappa by simply refusing to call the minority classes. +(double, double) _callRates(List> subjects, int epochSec, + double remCut, double deepCut) { + var n = 0, rem = 0, deep = 0; + for (final s in subjects) { + final scored = _scoreSubject(s, epochSec, remCut, deepCut); + if (scored == null) continue; + for (final p in scored.$2) { + n++; + if (p == 3) rem++; + if (p == 2) deep++; + } + } + return n == 0 ? (0.0, 0.0) : (rem / n, deep / n); +} + +/// (deep call rate, deep sensitivity, deep PPV, overall kappa) +(double, double, double, double) _deepMetrics( + List> subjects, + int epochSec, + double remCut, + double deepCut) { + final t = [], p = []; + for (final s in subjects) { + final scored = _scoreSubject(s, epochSec, remCut, deepCut); + if (scored == null) continue; + t.addAll(scored.$1); + p.addAll(scored.$2); + } + var called = 0, truthN = 0, hit = 0; + for (var i = 0; i < t.length; i++) { + if (p[i] == 2) called++; + if (t[i] == 2) truthN++; + if (t[i] == 2 && p[i] == 2) hit++; + } + return ( + t.isEmpty ? 0.0 : called / t.length, + truthN == 0 ? 0.0 : hit / truthN, + called == 0 ? 0.0 : hit / called, + _kappa(t, p), + ); +} + +double _scoreSet(List> subjects, int epochSec, + double remCut, double deepCut) { + final t = [], p = []; + for (final s in subjects) { + final scored = _scoreSubject(s, epochSec, remCut, deepCut); + if (scored == null) continue; + t.addAll(scored.$1); + p.addAll(scored.$2); + } + return _kappa(t, p); +} + +/// Trim a subject to the span between its first and last non-Wake epoch, i.e. +/// the sleep period. Returns null when the subject has no sleep at all or the +/// span is too short to stage. +Map? _trimToSleepPeriod(Map s) { + final stages = (s['stages'] as List).cast(); + var lo = -1, hi = -1; + for (var i = 0; i < stages.length; i++) { + if (stages[i] != 'Wake') { + if (lo < 0) lo = i; + hi = i; + } + } + if (lo < 0 || hi - lo + 1 < 60) return null; + final out = {'subject': s['subject']}; + for (final k in ['stages', ..._kFeatureKeys]) { + out[k] = (s[k] as List).sublist(lo, hi + 1); + } + return out; +} + +List _nums(dynamic raw, double whenNull) => [ + for (final v in (raw as List)) + v == null ? whenNull : (v as num).toDouble() + ]; + +void _report(String name, List t, List p) { + final cm = List.generate(4, (_) => List.filled(4, 0)); + for (var i = 0; i < t.length; i++) { + cm[t[i]][p[i]]++; + } + var correct = 0; + for (var i = 0; i < 4; i++) { + correct += cm[i][i]; + } + final acc = 100.0 * correct / t.length; + final buf = StringBuffer('${name.padRight(30)} ' + 'kappa=${_fmtKappa(_kappa(t, p))} ' + 'acc=${acc.toStringAsFixed(1).padLeft(5)}% '); + for (var i = 0; i < 4; i++) { + final rowSum = cm[i].reduce((a, b) => a + b); + var colSum = 0; + for (var r = 0; r < 4; r++) { + colSum += cm[r][i]; + } + final sens = rowSum == 0 ? 0.0 : 100.0 * cm[i][i] / rowSum; + final ppv = colSum == 0 ? 0.0 : 100.0 * cm[i][i] / colSum; + buf.write('${kClasses[i][0]} ' + '${sens.toStringAsFixed(1).padLeft(4)}/' + '${ppv.toStringAsFixed(1).padLeft(4)} '); + } + stdout.writeln(buf); +} + +double _kappa(List t, List p) { + final n = t.length; + if (n == 0) return double.nan; + final cm = List.generate(4, (_) => List.filled(4, 0)); + for (var i = 0; i < n; i++) { + cm[t[i]][p[i]]++; + } + var obs = 0.0, exp = 0.0; + for (var i = 0; i < 4; i++) { + obs += cm[i][i]; + var row = 0, col = 0; + for (var j = 0; j < 4; j++) { + row += cm[i][j]; + col += cm[j][i]; + } + exp += row * col / n; + } + if (n - exp == 0) return double.nan; + return (obs - exp) / (n - exp); +} diff --git a/tool/whoop_proportions.dart b/tool/whoop_proportions.dart new file mode 100644 index 0000000..ff29191 --- /dev/null +++ b/tool/whoop_proportions.dart @@ -0,0 +1,183 @@ +// Calibrate the stage-score CUTOFFS on our own device's feature distribution. +// +// The axis weights come from DREAMT (PSG ground truth, 99 wrist nights) and are +// what that corpus is good for: which features carry stage information and in +// which direction. The CUTOFFS are a different question — they sit on a +// robust-z score whose spread depends on how OUR feature extraction behaves, +// and the DREAMT-derived values demonstrably under-call REM on a real WHOOP +// night (62 min against a ~162 min reference). So sweep them here, on real +// device captures, and pick for NORMATIVE STAGE PROPORTIONS: +// +// REM 20-25% of TST Deep 13-23% of TST (Mitterling 2015 percentile +// curves; Boulos 2019 meta-regression, N=5,273, AASM scoring) +// +// There is no ground truth in these files — that is the point. This calibrates +// distributional plausibility only, and it is the honest limit of what we can +// do on our own device until we have concurrent PSG. +// +// Usage: dart run tool/whoop_proportions.dart +// Night JSON: {"day","start","end","onehz":[[ts,hr,ax,ay,az],..],"rr":[[tsMs,rr],..]} + +import 'dart:convert'; +import 'dart:io'; + +import 'package:openstrap_analytics/onehz.dart'; + +void main(List args) { + final dir = Directory(args.isEmpty ? '/tmp/sleepdbg' : args.first); + final files = dir + .listSync() + .whereType() + .where((f) => f.path.endsWith('.json')) + .toList() + ..sort((a, b) => a.path.compareTo(b.path)); + + final nights = <(String, List, List, List, + List, List)>[]; + var skipped = 0; + for (final f in files) { + // Recognise a night record by SHAPE, not by a hard-coded year in the + // filename. The directory also holds unrelated json (cursors, profiles), + // and a stale year filter would quietly stop matching in January. + final Object? decoded; + try { + decoded = jsonDecode(f.readAsStringSync()); + } catch (_) { + skipped++; + continue; + } + if (decoded is! Map || + decoded['onehz'] is! List || + decoded['rr'] is! List || + decoded['start'] is! int || + decoded['end'] is! int || + decoded['day'] is! String) { + skipped++; + continue; + } + final j = decoded; + final oh = j['onehz'] as List; + if (oh.length < 3600) { + skipped++; + continue; + } + final start = j['start'] as int, end = j['end'] as int; + final span = end - start; + // Mirror `_stageSessionCardio`: every slot carries ITS OWN timestamp (the + // RR windows are centred on accel[mid].tsMs, so a stale timestamp + // mis-centres them), accel is carried forward only within a bounded gap, + // and anything beyond that is UNUSABLE rather than a fabricated (0,0,1) + // "perfectly still" sample. One real night here has 20% of its epochs + // without HR, so getting this wrong silently invents stillness. + const maxCarrySec = 120; + final accel = List.generate( + span, (i) => AccelSample((start + i) * 1000.0, 0, 0, 1), + growable: false); + final hr = List.filled(span, 0); + final usable = List.filled(span, false); + var last = -1; + final byTs = >{}; + for (final r in oh) { + byTs[r[0] as int] = [ + (r[1] as num).toDouble(), + (r[2] as num).toDouble(), + (r[3] as num).toDouble(), + (r[4] as num).toDouble(), + ]; + } + for (var i = 0; i < span; i++) { + final ts = start + i; + final g = byTs[ts]; + if (g != null) { + accel[i] = AccelSample(ts * 1000.0, g[1], g[2], g[3]); + hr[i] = g[0]; + last = i; + usable[i] = true; + } else if (last >= 0 && i - last <= maxCarrySec) { + final p = accel[i - 1]; + accel[i] = AccelSample(ts * 1000.0, p.x, p.y, p.z); + usable[i] = true; + } + } + final rrMs = [], rrTs = []; + for (final r in (j['rr'] as List)) { + rrTs.add((r[0] as num).toDouble()); + rrMs.add((r[1] as num).toDouble()); + } + nights.add((j['day'] as String, hr, accel, rrMs, rrTs, usable)); + } + + stdout.writeln('${nights.length} real device nights' + '${skipped > 0 ? ' ($skipped file(s) skipped: not a night record)' : ''}'); + if (nights.isEmpty) { + stderr.writeln('no night records found in ${dir.path}'); + exitCode = 1; + return; + } + stdout.writeln('target: REM 20-25% of TST, Deep 13-23% of TST'); + stdout.writeln(''); + stdout.writeln('remCut deepCut | median %REM of TST median %Deep of TST ' + ' median TST(h)'); + for (final rc in [0.0, 0.2, 0.3, 0.4, 0.5, 0.6, 0.8]) { + for (final dc in [0.0, 0.2, 0.3, 0.4]) { + final rems = [], deeps = [], tsts = []; + for (final (_, hr, accel, rrMs, rrTs, usable) in nights) { + // Stage each USABLE run separately, exactly as the production path + // does, and aggregate — feeding a gap straight through would let + // carried-forward stillness masquerade as sleep. + var rem = 0, nrem = 0, deep = 0; + var i = 0; + while (i < usable.length) { + if (!usable[i]) { + i++; + continue; + } + var jx = i; + while (jx < usable.length && usable[jx]) { + jx++; + } + if (jx - i >= 90) { + final r = cardioStager( + hr.sublist(i, jx), accel.sublist(i, jx), + rrMs: rrMs, + rrTsMs: rrTs, + remScoreCut: rc, + deepScoreCut: dc); + for (var e = 0; e < r.base.stages.length; e++) { + switch (r.base.stages[e]) { + case SleepStage.rem: + rem++; + case SleepStage.nrem: + nrem++; + if (e < r.deepFlag.length && r.deepFlag[e]) deep++; + case SleepStage.wake: + break; + } + } + } + i = jx; + } + final tst = rem + nrem; + if (tst == 0) continue; + rems.add(100 * rem / tst); + deeps.add(100 * deep / tst); + tsts.add(tst * 30 / 3600); + } + if (rems.isEmpty) continue; + final flagR = (median(rems)! >= 20 && median(rems)! <= 25) ? ' REM✓' : ''; + final flagD = (median(deeps)! >= 13 && median(deeps)! <= 23) ? ' Deep✓' : ''; + stdout.writeln(' ${rc.toStringAsFixed(1)} ${dc.toStringAsFixed(1)} ' + '| ${median(rems)!.toStringAsFixed(1).padLeft(5)}%' + ' ${median(deeps)!.toStringAsFixed(1).padLeft(5)}%' + ' ${median(tsts)!.toStringAsFixed(1)}$flagR$flagD'); + } + } +} + +double? median(List v) { + if (v.isEmpty) return null; + final s = [...v]..sort(); + final n = s.length; + return n.isOdd ? s[n ~/ 2] : (s[n ~/ 2 - 1] + s[n ~/ 2]) / 2; +} +