diff --git a/ALGORITHMS.md b/ALGORITHMS.md index f43d820..62d7b0a 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -56,7 +56,8 @@ Grouped by family (subdirectory under `lib/src/onehz/`). File paths are relative |---|---|---|---| | `hrvTime` | `clinical/hrv_time.dart` | RMSSD/SDNN/pNNx | standard time-domain HRV | | `nocturnalRmssd` | `clinical/hrv_time.dart` | median-of-5-min-window nightly RMSSD | — | -| `sleepSessionWindowedRmssd` | `clinical/hrv_time.dart` | mean-of-5-min-window RMSSD with Malik ectopic rejection | Malik et al. | +| `sleepSessionWindowedRmssd` / `sleepSessionRmssdDetail` | `clinical/hrv_time.dart` | mean-of-5-min-window RMSSD; windows need ≥ 20 clean successive differences (ultra-short RMSSD literature); Malik 20 % ectopic rejection | Malik et al.; Munoz et al. 2015; Baek et al. 2015 | +| `rrCoverage` | `clinical/hrv_time.dart` | Σ plausible RR ÷ wall-clock span; > 1.10 means duplicated or interleaved beats, and `hrvTime` / `nocturnalRmssd` / `nightHrvShape` / the sleep-session headline refuse RMSSD; the headline also drops any 5-min window holding more beat-time than its length | — (integrity check; contiguous runs measure 0.963–1.001) | | `hrvFreq` | `clinical/hrv_freq.dart` | LF/HF via Lomb-Scargle periodogram on native (unevenly-sampled) beat times | Laguna, Moody & Mark 1998; Bigger 1992 | | `decelerationCapacity` / `accelerationCapacity` | `clinical/prsa.dart` | Phase-rectified signal averaging (DC/AC) | Bauer et al. 2006 | | `nocturnalRhr` / `hrDip` | `clinical/nocturnal.dart` | nocturnal resting HR + dip classification | — | diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index b5c7438..70cd001 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -470,6 +470,113 @@ String _jitterNote(double acf1) => '$kNnDiffAcf1Floor), so RMSSD/pNN50 would measure beat-timing jitter, not ' 'vagal tone'; +/// Σ reported RR ÷ elapsed wall time above which the stream cannot be one +/// heart's beats (it banks more beat-time than time passed): a double ingest or +/// two interleaved streams. Contiguous runs measure 0.963 (gen4), 0.999 (W5), +/// 1.001 (MG) — see rr_correction.dart `_beatTimes` — so 1.10 has margin. +/// Duplicated beats add zero differences and DEFLATE RMSSD by ~1/√2; +/// interleaved streams inflate it. Every other gate here passes both. +const double kRrCoverageCeiling = 1.10; + +/// Shortest wall span [rrCoverage] will judge: whole-second stamps make a +/// shorter one meaningless. +const double kRrCoverageMinSpanSec = 600; + +/// Intervals [rrCoverage] counts as beat-time, ms. Wider than the cleaners' +/// 2000 ms on purpose: a 2000–2400 ms beat is a slow heart (25–30 bpm), still +/// real elapsed time; only a glitch (e.g. a 65 535 ms sentinel) must not be +/// summed, or it would fake an over-count on a short span. +const double _rrCoverageMinMs = 300; +const double _rrCoverageMaxMs = 2400; + +bool _rrCoverageCounts(double rrMs) => + rrMs >= _rrCoverageMinMs && rrMs <= _rrCoverageMaxMs; + +/// How much beat-time an RR stream banks against the wall clock it spans. +class RrCoverage { + /// Σ plausible RR ÷ wall span. Below 1 on any gap; above 1 is impossible + /// for one heart. + final double coverage; + final double sumRrSec; + final double spanSec; + final int beats; + + /// Intervals outside 300–2400 ms: counted, never summed. + final int implausibleBeats; + + /// Exact (ts, rr) repeats of an earlier beat, wherever they sit in the + /// input. A DIAGNOSTIC only: on whole-second stamps two equal beats in one + /// record repeat legitimately. + final int duplicateBeats; + const RrCoverage({ + required this.coverage, + required this.sumRrSec, + required this.spanSec, + required this.beats, + required this.implausibleBeats, + required this.duplicateBeats, + }); + bool get overCounted => coverage > kRrCoverageCeiling; + Map toJson() => { + 'rr_coverage': round6(coverage), + 'sum_rr_sec': round6(sumRrSec), + 'span_sec': round6(spanSec), + 'beats': beats, + 'implausible_beats': implausibleBeats, + 'duplicate_beats': duplicateBeats, + }; +} + +/// [RrCoverage] of raw RR [rrMs] against their beat-END epoch times [rrTsMs] +/// (same length; order does not matter). The span is `latest − earliest` beat +/// end plus the earliest beat's own interval (it began before its end stamp) — +/// but only when that interval is itself plausible. An implausible one is +/// never summed, so it must not stretch the denominator either (a 65 535 ms +/// glitch there hid a 12.5 % over-count). Null when fewer than 2 beats, the +/// lengths differ, or the span is under [kRrCoverageMinSpanSec]. +RrCoverage? rrCoverage(List rrMs, List rrTsMs) { + if (rrMs.length < 2 || rrMs.length != rrTsMs.length) return null; + var lo = 0, hi = 0; + for (var i = 1; i < rrTsMs.length; i++) { + if (rrTsMs[i] < rrTsMs[lo]) lo = i; + if (rrTsMs[i] > rrTsMs[hi]) hi = i; + } + final first = rrMs[lo]; + final spanSec = (rrTsMs[hi] - rrTsMs[lo] + + (_rrCoverageCounts(first) ? first : 0)) / + 1000.0; + if (!spanSec.isFinite || spanSec < kRrCoverageMinSpanSec) return null; + var sum = 0.0; + var implausible = 0; + var dup = 0; + final seen = <(double, double)>{}; + for (var i = 0; i < rrMs.length; i++) { + final v = rrMs[i]; + if (_rrCoverageCounts(v)) { + sum += v; + } else { + implausible++; + } + if (!seen.add((rrTsMs[i], v))) dup++; + } + final sumSec = sum / 1000.0; + return RrCoverage( + coverage: sumSec / spanSec, + sumRrSec: sumSec, + spanSec: spanSec, + beats: rrMs.length, + implausibleBeats: implausible, + duplicateBeats: dup, + ); +} + +/// The refusal note every RMSSD-family estimator gives an over-counted stream. +String rrOvercountNote(RrCoverage c) => _overcountNote(c); + +String _overcountNote(RrCoverage c) => + 'rr_overcount:coverage=${round6(c.coverage)} — more beat-time than ' + 'elapsed time; the RR stream holds duplicated or interleaved beats'; + class HrvTime { final double? rmssd; // ms final double? sdnn; // ms @@ -507,11 +614,14 @@ class HrvTime { /// the upstream corrector rejected (0..1), folded into confidence exactly as /// `hrvFreq` and `irregularBeatScreen` already do. Returns an absent Metric when /// there are too few beats; RMSSD/pNN50 alone go null when the successive -/// differences fail [kNnDiffAcf1Floor]. +/// differences fail [kNnDiffAcf1Floor], or when [coverage] (of the raw RR this +/// NN was cleaned from) is [RrCoverage.overCounted] — then the breathing-line +/// fallback does not run either. Without [coverage] that check is skipped. Metric hrvTime( List nnMs, { List? nnTimesMs, double artifactFraction = 0.0, + RrCoverage? coverage, }) { const inputs = ['rr_cleaned']; if (nnMs.length < 2) { @@ -564,11 +674,17 @@ Metric hrvTime( // always advised. final acf1 = nnDiffAcf1(runs); final jittery = _jitterRefused(acf1, runs); - var rmssd = (pairs > 0 && !jittery) ? math.sqrt(ssd / pairs) : null; - final pnn50 = (pairs > 0 && !jittery) ? 100.0 * nn50 / pairs : null; + // More beat-time than elapsed: the stream itself is wrong, so neither the + // differences nor any breathing line in them describe one heart. + final overCounted = coverage?.overCounted == true; + final usable = pairs > 0 && !jittery && !overCounted; + var rmssd = usable ? math.sqrt(ssd / pairs) : null; + final pnn50 = usable ? 100.0 * nn50 / pairs : null; // Refused: a long record may still show a breathing line steady in Hz // ([_steadyBreathingWindows]); then RMSSD alone comes from those 5-min windows. - if (jittery && gapAware) rmssd = _breathingRmssd(nnMs, nnTimesMs); + if (jittery && gapAware && !overCounted) { + rmssd = _breathingRmssd(nnMs, nnTimesMs); + } final byBreathing = jittery && rmssd != null; final sdnn = stddev(nnMs); @@ -590,7 +706,7 @@ Metric hrvTime( // were ~pure noise. The beat-count term is capped BEFORE the quality terms // multiply it; multiplying first let an all-night beat count (n/250 ≈ 100) // swallow any penalty and re-clamp to 0.95 regardless. - final conf = byBreathing + final conf = byBreathing || overCounted ? 0.3 : ((nnMs.length / 250.0).clamp(0.0, 1.0) // ~250 beats ≈ 5 min * @@ -610,7 +726,10 @@ Metric hrvTime( confidence: conf, tier: Tier.high, inputs_used: inputs, - note: jittery + note: overCounted + ? '${_overcountNote(coverage!)}. SDNN/SDANN survive it and are the ' + 'lead here. PRV not ECG-HRV.' + : jittery ? '${_jitterNote(acf1!)}. SDNN/SDANN survive it and are the lead here. ' '${byBreathing ? 'RMSSD only from 5-min windows on a breathing ' 'line steady in Hz, at floor confidence. ' : ''}PRV not ECG-HRV.' @@ -670,8 +789,9 @@ double? _breathingRmssd(List nn, List times) { /// true at the window's MIDPOINT second. /// /// Returns a Metric whose value is the median-of-windows RMSSD (ms). Keeps the -/// PRV-not-ECG honesty note. Absent when there are too few usable windows, or -/// when the night's successive differences fail [kNnDiffAcf1Floor]. A window +/// PRV-not-ECG honesty note. Absent when there are too few usable windows, when +/// the night's successive differences fail [kNnDiffAcf1Floor], or when +/// [coverage] (of the raw RR) is [RrCoverage.overCounted]. A window /// contributes only if it holds [minBeatsPerWindow] differences between beats /// that are ADJACENT IN TIME, not merely adjacent in the compacted NN list. Metric nocturnalRmssd( @@ -680,8 +800,16 @@ Metric nocturnalRmssd( double windowMs = 300000.0, int minBeatsPerWindow = 5, List? stageMaskPerSec, + RrCoverage? coverage, }) { const inputs = ['rr_cleaned', 'beat_times']; + if (coverage != null && coverage.overCounted) { + return Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: _overcountNote(coverage), + ); + } if (nnMs.length != nnTimesMs.length || nnMs.length < minBeatsPerWindow + 1) { return const Metric.absent( tier: Tier.high, @@ -791,19 +919,71 @@ Metric nocturnalRmssd( ); } +/// Fewest successive differences a 5-min window needs before its RMSSD joins +/// the nightly mean. RMSSD's relative sampling error is ~1/√(2n): 71 % at n = 1, +/// 32 % at n = 5, 16 % at n = 20 (cf. night_hrv_shape.dart). The +/// ultra-short-RMSSD literature (Munoz et al. 2015; Baek et al. 2015) puts the +/// shortest usable recording at roughly 10–30 s — 20 differences is 20–25 s at +/// a sleeping 48–60 bpm. The windows this drops sit where the signal is +/// disturbed (session edges, dropouts, ectopy), and those of ≤ 2 beats are the +/// ones `_cleanWindowRuns` cannot Malik-filter. `nocturnalRmssd`, a secondary +/// median estimator published under its own key, keeps its own floor of 5. +/// A floor, not a weight: weighting by n would weight by HEART RATE, letting +/// high-HR, low-RMSSD REM/arousal windows pull the mean down. +const int kMinDiffsPerRmssdWindow = 20; + +/// The nightly headline RMSSD plus the diagnostics a `Metric` cannot +/// carry. Present only when the headline is. +class SessionRmssd { + final double rmssd; // ms — the headline + final int windows; // 5-min windows that contributed + final int overCountedWindows; // windows dropped for more beat-time than time + final int thinWindows; // windows dropped for 0..floor−1 differences + final int minDiffsPerWindow; // the floor, which travels with the number + final double? diffAcf1; // pooled over the session's windows + final double? rrCoverage; // [RrCoverage.coverage] of the session's beats + const SessionRmssd({ + required this.rmssd, + required this.windows, + this.overCountedWindows = 0, + this.thinWindows = 0, + this.minDiffsPerWindow = kMinDiffsPerRmssdWindow, + this.diffAcf1, + this.rrCoverage, + }); + Map toJson() => { + 'rmssd_ms': round6(rmssd), + 'windows': windows, + 'overcounted_windows': overCountedWindows, + 'thin_windows': thinWindows, + 'min_diffs_per_window': minDiffsPerWindow, + if (diffAcf1 != null) 'diff_acf1': round6(diffAcf1!), + if (rrCoverage != null) 'rr_coverage': round6(rrCoverage!), + }; +} + /// Sleep-session nightly RMSSD (ms) as the arithmetic mean of cleaned /// consecutive 5-minute window RMSSDs. /// /// Split the detected sleep session into consecutive 5-minute windows, apply a /// simple RR cleaner (range-filter [300, 2000] ms + Malik-style ectopic /// rejection against a local median), compute RMSSD inside each valid window, -/// then return the ARITHMETIC MEAN across windows. This is intentionally +/// then return the ARITHMETIC MEAN across windows. A window joins the mean only +/// with at least [minDiffsPerWindow] clean successive differences +/// ([kMinDiffsPerRmssdWindow]): an unweighted mean let a window of a handful of +/// differences across an arousal or a dropout edge count as much as a full +/// window of ~300, and those thin windows cluster where the signal is +/// disturbed, so the bias was upward. Thin windows — 0 to +/// [minDiffsPerWindow] − 1 differences, including one the cleaner left with +/// none — are counted ([SessionRmssd.thinWindows]), not silently lost, and +/// contribute nothing to the jitter verdict either. This is intentionally /// distinct from [nocturnalRmssd], which uses cleaned NN + /// median-of-windows robustness. /// /// This is the nightly HEADLINE (→ `ln_rmssd` → readiness), so it refuses /// rather than approximates: absent when the successive differences fail -/// [kNnDiffAcf1Floor]. +/// [kNnDiffAcf1Floor], and absent when the session's beats bank more time than +/// elapsed ([kRrCoverageCeiling]). /// /// [rrMs]/[rrTsMs] are the raw RR intervals and their beat-end epoch times in /// milliseconds. [startSec]/[endSec] bound the chosen sleep session in epoch @@ -815,6 +995,34 @@ Metric sleepSessionWindowedRmssd( required int startSec, required int endSec, int windowSec = 300, + int minDiffsPerWindow = kMinDiffsPerRmssdWindow, +}) { + final m = sleepSessionRmssdDetail(rrMs, rrTsMs, + startSec: startSec, + endSec: endSec, + windowSec: windowSec, + minDiffsPerWindow: minDiffsPerWindow); + return m.present + ? Metric( + value: m.value!.rmssd, + confidence: m.confidence, + tier: m.tier, + inputs_used: m.inputs_used, + note: m.note, + ) + : Metric.absent( + tier: m.tier, inputs_used: m.inputs_used, note: m.note); +} + +/// [sleepSessionWindowedRmssd] with its diagnostics ([SessionRmssd]). The two +/// are one computation; this is the one that does it. +Metric sleepSessionRmssdDetail( + List rrMs, + List rrTsMs, { + required int startSec, + required int endSec, + int windowSec = 300, + int minDiffsPerWindow = kMinDiffsPerRmssdWindow, }) { const inputs = ['rr_sleep_window']; if (startSec <= 0 || @@ -822,7 +1030,7 @@ Metric sleepSessionWindowedRmssd( rrMs.isEmpty || rrTsMs.isEmpty || rrMs.length != rrTsMs.length) { - return const Metric.absent( + return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'invalid or empty RR session window', @@ -831,28 +1039,62 @@ Metric sleepSessionWindowedRmssd( final buckets = >{}; final bucketsTs = >{}; + final inRr = []; + final inTs = []; for (var i = 0; i < rrMs.length; i++) { final tsSec = (rrTsMs[i] / 1000.0).round(); if (tsSec < startSec || tsSec >= endSec) continue; final idx = ((tsSec - startSec) ~/ windowSec); (buckets[idx] ??= []).add(rrMs[i]); (bucketsTs[idx] ??= []).add(rrTsMs[i]); + inRr.add(rrMs[i]); + inTs.add(rrTsMs[i]); } if (buckets.isEmpty) { - return const Metric.absent( + return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'no RR beats inside the session window', ); } + final cov = rrCoverage(inRr, inTs); + if (cov != null && cov.overCounted) { + return Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: _overcountNote(cov), + ); + } final rmssds = []; final runs = >[]; // pooled jitter floor — see [nocturnalRmssd] final perWindow = >>[]; final meanRr = []; + var overCountedWindows = 0; + var thinWindows = 0; final indices = buckets.keys.toList()..sort(); for (final idx in indices) { + // The session-wide ratio above can be diluted below the ceiling by a gap + // elsewhere in the night while one stretch holds every beat twice. Each + // window is judged on its own too: its beats cannot bank more beat-time + // than the window has seconds — the session's last window may be cut + // short by [endSec] — plus its earliest beat's own interval, which began + // before its end stamp (as in [rrCoverage]). Dropped, not averaged in. + final winRr = buckets[idx]!, winTs = bucketsTs[idx]!; + final winStart = startSec + idx * windowSec; + final winSec = math.min(winStart + windowSec, endSec) - winStart; + var bankedMs = 0.0; + var lo = 0; + for (var i = 0; i < winRr.length; i++) { + if (winTs[i] < winTs[lo]) lo = i; + if (_rrCoverageCounts(winRr[i])) bankedMs += winRr[i]; + } + final overhangMs = _rrCoverageCounts(winRr[lo]) ? winRr[lo] : 0.0; + if (bankedMs > kRrCoverageCeiling * (winSec * 1000 + overhangMs)) { + overCountedWindows++; + continue; + } final rrRuns = [ for (final r in _cleanWindowRuns(buckets[idx]!, bucketsTs[idx]!)) if (r.length >= 2) r @@ -868,35 +1110,52 @@ Metric sleepSessionWindowedRmssd( nd++; } } - // Same floor as [nocturnalRmssd]'s minBeatsPerWindow: a window with one or - // two differences would otherwise weigh as much as a full one in the mean. - if (nd < 5) continue; + // A window with a handful of differences would otherwise weigh as much as + // a full one in the mean. Counted, not silently lost — the detail and the + // note both say how many. A window whose beats the cleaner reduced to no + // difference at all is the thinnest of them, and counts too. + if (nd < minDiffsPerWindow) { + thinWindows++; + continue; + } runs.addAll(diffRuns); perWindow.add(diffRuns); meanRr.add(mean([for (final r in rrRuns) ...r])!); rmssds.add(math.sqrt(ssd / nd)); } + // Every note says what was dropped before the estimate, so a blank or a + // thin headline can be told apart from a jitter refusal. + final droppedWhy = [ + if (thinWindows > 0) + '$thinWindows window(s) had fewer than $minDiffsPerWindow clean ' + 'successive differences', + if (overCountedWindows > 0) + '$overCountedWindows window(s) held more beat-time than elapsed time', + ]; + final dropped = droppedWhy.isEmpty ? '' : ' (${droppedWhy.join('; ')})'; + // 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); final refused = _jitterRefused(acf1, runs); if (rmssds.isEmpty && !refused) { - return const Metric.absent( + return Metric.absent( tier: Tier.high, inputs_used: inputs, - note: 'no valid 5-min windows for sleep-session RMSSD', + note: 'no valid 5-min windows for sleep-session RMSSD$dropped', ); } final kept = _keptWindows(acf1, refused, rmssds, perWindow, meanRr); if (kept == null) { - return Metric.absent( + return Metric.absent( tier: Tier.high, inputs_used: inputs, note: refused - ? _jitterNote(acf1!) - : '${_jitterNote(acf1!)}; no 5-min window clears it on its own', + ? '${_jitterNote(acf1!)}$dropped' + : '${_jitterNote(acf1!)}; no 5-min window clears it on its own' + '$dropped', ); } @@ -905,13 +1164,21 @@ Metric sleepSessionWindowedRmssd( ? 0.3 : ((kept.rmssds.length / 12.0).clamp(0.0, 1.0) * _acf1Quality(acf1)) .clamp(0.3, 0.95); - return Metric( - value: meanRmssd, + return Metric( + value: SessionRmssd( + rmssd: meanRmssd, + windows: kept.rmssds.length, + overCountedWindows: overCountedWindows, + thinWindows: thinWindows, + minDiffsPerWindow: minDiffsPerWindow, + diffAcf1: acf1, + rrCoverage: cov?.coverage, + ), confidence: conf, tier: Tier.high, inputs_used: inputs, - note: 'sleep-session HRV: mean RMSSD over cleaned 5-min windows.' - '${kept.note}', + note: 'sleep-session HRV: mean RMSSD over cleaned 5-min windows' + '$dropped.${kept.note}', ); } diff --git a/lib/src/onehz/sleep/night_hrv_shape.dart b/lib/src/onehz/sleep/night_hrv_shape.dart index 8d9b960..7778fba 100644 --- a/lib/src/onehz/sleep/night_hrv_shape.dart +++ b/lib/src/onehz/sleep/night_hrv_shape.dart @@ -110,13 +110,27 @@ const int kMinBeatsPerHrvBin = 300; /// [binMin] is the bin width; the IDEAS entry says to widen to 45–60 min if /// 30-min bins look wobbly on real nights, so it is a parameter, not a /// constant. +/// +/// [coverage] is [rrCoverage] of the raw RR this NN was cleaned from; when it +/// is [RrCoverage.overCounted] the shape is absent, as every RMSSD is. Metric nightHrvShape( List nnMs, List nnTimesMs, { double binMin = 30, int minBeatsPerBin = kMinBeatsPerHrvBin, + RrCoverage? coverage, }) { const inputs = ['rr_cleaned', 'beat_times']; + // Every bin is an RMSSD of the same stream: if its raw RR banks more + // beat-time than elapsed (duplicated or interleaved beats), no bin is one + // heart's variability, and a shape drawn from them is not either. + if (coverage != null && coverage.overCounted) { + return Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: rrOvercountNote(coverage), + ); + } if (nnMs.length != nnTimesMs.length || nnMs.length < minBeatsPerBin) { return Metric.absent( tier: Tier.high, diff --git a/test/onehz/clinical_test.dart b/test/onehz/clinical_test.dart index 6d2a5a7..80d1278 100644 --- a/test/onehz/clinical_test.dart +++ b/test/onehz/clinical_test.dart @@ -3,6 +3,36 @@ import 'dart:math' as math; import 'package:test/test.dart'; import 'package:openstrap_analytics/onehz.dart'; +/// Respiratory sinus arrhythmia: NN sampled at its own beat times, so the +/// breathing phase advances by the beat's real duration. Noise-free. +List rsaNn({ + required double hrBpm, + required double respBrpm, + required double ampMs, + int beats = 3000, +}) { + final base = 60000.0 / hrBpm, f = respBrpm / 60.0; + final out = []; + var tSec = 0.0; + for (var i = 0; i < beats; i++) { + final v = base + ampMs * math.sin(2 * math.pi * f * tSec); + out.add(v); + tSec += v / 1000.0; + } + return out; +} + +/// Cumulative beat-END times (ms) for [nn], starting at [t0Ms]. +List beatEnds(List nn, {double t0Ms = 0}) { + final out = []; + var t = t0Ms; + for (final v in nn) { + t += v; + out.add(t); + } + return out; +} + void main() { group('time-domain HRV (hand-computed)', () { test('RMSSD/SDNN/pNN50 on a constant-then-stepped NN series', () { @@ -469,6 +499,279 @@ void main() { }); }); + group('RR over-count gate (Σ RR ÷ wall span)', () { + test('one beat per second reads 1.0; every beat duplicated reads 2.0', () { + final rr = [for (var i = 0; i < 1200; i++) 1000.0]; // 20 min + final ts = [for (var i = 0; i < 1200; i++) 1e12 + (i + 1) * 1000.0]; + final c = rrCoverage(rr, ts)!; + expect(c.coverage, closeTo(1.0, 1e-9)); + expect(c.overCounted, isFalse); + expect(c.beats, 1200); + expect(c.duplicateBeats, 0); + final rr2 = [for (final v in rr) ...[v, v]]; + final ts2 = [for (final t in ts) ...[t, t]]; + final d = rrCoverage(rr2, ts2)!; + expect(d.coverage, closeTo(2.0, 1e-9)); + expect(d.overCounted, isTrue); + expect(d.duplicateBeats, 1200); + expect(d.toJson()['rr_coverage'], closeTo(2.0, 1e-9)); + }); + + test('under 10 min, mismatched or too few beats: null', () { + final rr = [for (var i = 0; i < 500; i++) 1000.0]; + final ts = [for (var i = 0; i < 500; i++) (i + 1) * 1000.0]; + expect(rrCoverage(rr, ts), isNull); + expect(rrCoverage([1000.0], [1000.0]), isNull); + expect(rrCoverage(rr, ts.sublist(1)), isNull); + }); + + test('an implausible interval is counted, never summed', () { + // One 65 535 ms glitch would otherwise fake an over-count. + final rr = [for (var i = 0; i < 1200; i++) 1000.0]..[600] = 65535.0; + final ts = [for (var i = 0; i < 1200; i++) (i + 1) * 1000.0]; + final c = rrCoverage(rr, ts)!; + expect(c.implausibleBeats, 1); + expect(c.coverage, closeTo(1199 / 1200, 1e-9)); + expect(c.overCounted, isFalse); + }); + + test('an implausible FIRST interval cannot widen the span', () { + // Every 8th beat duplicated reads 1.125 and is refused. A 65 535 ms + // glitch as the first interval used to stretch the span by 64.5 s while + // never being summed, and the same stream read ~1.07 and passed. + final rr = [], ts = []; + for (var i = 0; i < 1200; i++) { + final t = 1e12 + (i + 1) * 1000.0; + rr.add(1000.0); + ts.add(t); + if (i % 8 == 7) { + rr.add(1000.0); + ts.add(t); + } + } + expect(rrCoverage(rr, ts)!.coverage, closeTo(1.125, 1e-3)); + expect(rrCoverage(rr, ts)!.overCounted, isTrue); + rr[0] = 65535.0; + final c = rrCoverage(rr, ts)!; + expect(c.implausibleBeats, 1); + expect(c.overCounted, isTrue, reason: 'coverage ${c.coverage}'); + }); + + test('a gappy night reads well under 1 and is never refused for it', () { + final rr = [for (var i = 0; i < 1200; i++) 1000.0]; + final ts = [ + for (var i = 0; i < 1200; i++) (i + 1) * 1000.0 + (i >= 600 ? 1.2e6 : 0) + ]; + expect(rrCoverage(rr, ts)!.coverage, closeTo(0.5, 1e-9)); + // A real night with a 30-min strap-off gap: every estimator the + // coverage reaches still publishes. + final nn = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 3600); + final t0 = beatEnds(nn.sublist(0, 1800), t0Ms: 1e12); + final t1 = beatEnds(nn.sublist(1800), t0Ms: t0.last + 1800e3); + final t = [...t0, ...t1]; + final c = rrCoverage(nn, t)!; + expect(c.coverage, lessThan(0.7)); + final h = hrvTime(nn, nnTimesMs: t, coverage: c); + expect(h.value!.rmssd, isNotNull, reason: h.note); + final n = nocturnalRmssd(nn, t, coverage: c); + expect(n.present, isTrue, reason: n.note); + final d = sleepSessionRmssdDetail(nn, t, + startSec: (t.first / 1000).floor(), + endSec: (t.last / 1000).ceil() + 1); + expect(d.present, isTrue, reason: d.note); + expect(d.value!.rrCoverage!, lessThan(0.7)); + }); + + test('the headline refuses an over-counted session', () { + // HR 60 / 12 br/min: the jitter screen keeps it, so the over-count is + // the only thing under test. + final rr = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 25200); + final ts = beatEnds(rr, t0Ms: 1e12); + final startSec = (ts.first / 1000).floor(); + final endSec = (ts.last / 1000).ceil() + 1; + final ok = sleepSessionRmssdDetail(rr, ts, + startSec: startSec, endSec: endSec); + expect(ok.present, isTrue, reason: ok.note); + expect(ok.value!.rmssd, + sleepSessionWindowedRmssd(rr, ts, startSec: startSec, endSec: endSec) + .value); + expect(ok.value!.rrCoverage!, closeTo(1.0, 0.01)); + expect(ok.value!.windows, greaterThan(80)); + expect(ok.value!.toJson()['rr_coverage'], closeTo(1.0, 0.01)); + + final rr2 = [for (final v in rr) ...[v, v]]; + final ts2 = [for (final t in ts) ...[t, t]]; + final dup = sleepSessionRmssdDetail(rr2, ts2, + startSec: startSec, endSec: endSec); + expect(dup.present, isFalse); + expect(dup.note, startsWith('rr_overcount:coverage=')); + expect( + sleepSessionWindowedRmssd(rr2, ts2, + startSec: startSec, endSec: endSec) + .present, + isFalse); + + final c = rrCoverage(rr2, ts2)!; + final h = hrvTime(rr2, coverage: c); + expect(h.value!.rmssd, isNull); + expect(h.value!.pnn50, isNull); + expect(h.value!.sdnn, isNotNull, reason: 'SDNN survives'); + expect(h.note, startsWith('rr_overcount:coverage=')); + final n = nocturnalRmssd(rr2, ts2, coverage: c); + expect(n.present, isFalse); + expect(n.note, startsWith('rr_overcount:coverage=')); + }); + + test('an over-counted stream is not rescued by the breathing fallback', + () { + // A slow heart drifting 50 ± 8 bpm over a 3 h cycle, breathing steady at + // 18 br/min, beat-time jitter 20 ms: the pooled jitter gate refuses it + // and hrvTime keeps it through its breathing-line fallback. If the raw + // RR it was cleaned from banks more beat-time than elapsed, the stream + // is not one heart's beats: no fallback either. + final rnd = math.Random(0); + final nn = [], t = []; + var ms = 0.0, e0 = 0.0; + while (ms < 6 * 3600e3) { + final base = 60000 / (50 + 8 * math.sin(2 * math.pi * ms / 10800e3)); + final e1 = (rnd.nextDouble() - 0.5) * 20; + final v = base + 10 * math.sin(2 * math.pi * 0.3 * ms / 1000) + e1 - e0; + e0 = e1; + ms += v; + nn.add(v); + t.add(ms); + } + final rescued = hrvTime(nn, nnTimesMs: t); + expect(rescued.value!.rmssd, isNotNull); + expect(rescued.note, contains('steady in Hz'), + reason: 'the fallback path is what is under test'); + final raw2 = [for (final v in nn) ...[v, v]]; + final t2 = [for (final x in t) ...[x, x]]; + final c = rrCoverage(raw2, t2)!; + expect(c.overCounted, isTrue); + final h = hrvTime(nn, nnTimesMs: t, coverage: c); + expect(h.value!.rmssd, isNull); + expect(h.value!.sdnn, isNotNull); + expect(h.confidence, 0.3); + expect(h.note, startsWith('rr_overcount:coverage=')); + }); + + test('a duplicated stretch is caught per window even when a gap dilutes ' + 'the night', () { + // 4 h of honest beats, a 2 h strap-off gap, then 2 h with every beat + // stored twice: Σ RR equals the 8 h span, so the night-wide ratio reads + // ≈ 1.0 and passes. Each 5-min window of the duplicated stretch banks + // ~600 s in 300 s and is dropped. A window straddling the start of the + // stretch banks ≤ 110 % and may stay — the ceiling's own tolerance. + final clean = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 14400); + final ts = beatEnds(clean, t0Ms: 1e12); + final dupNn = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 7200); + final dupTs = beatEnds(dupNn, t0Ms: ts.last + 2 * 3600e3); + final rr = [...clean, for (final v in dupNn) ...[v, v]]; + final t = [...ts, for (final x in dupTs) ...[x, x]]; + final startSec = (t.first / 1000).floor(); + final endSec = (t.last / 1000).ceil() + 1; + expect(rrCoverage(rr, t)!.overCounted, isFalse, + reason: 'the gap dilutes the night-wide ratio'); + final honest = sleepSessionRmssdDetail(clean, ts, + startSec: startSec, endSec: endSec); + final d = sleepSessionRmssdDetail(rr, t, + startSec: startSec, endSec: endSec); + expect(d.present, isTrue, reason: d.note); + expect(d.value!.overCountedWindows, greaterThanOrEqualTo(23)); + expect(d.value!.toJson()['overcounted_windows'], + d.value!.overCountedWindows); + expect(d.value!.windows - honest.value!.windows, inInclusiveRange(0, 1)); + expect(d.value!.rmssd, closeTo(honest.value!.rmssd, 0.02 * honest.value!.rmssd)); + // Only the duplicated stretch: absent, and the note says why. + final only = sleepSessionRmssdDetail( + [for (final v in dupNn) ...[v, v]], [for (final x in dupTs) ...[x, x]], + startSec: (dupTs.first / 1000).floor(), + endSec: (dupTs.last / 1000).ceil() + 1); + expect(only.present, isFalse); + expect(only.note, startsWith('rr_overcount')); + }); + + test('a duplicated tail in the session\'s short last window is caught', + () { + // 4 h of honest beats, a 2 h gap, then 2 min with every beat stored + // twice, and the session ends 2 min into that last 5-min window. The + // tail banks ~240 s — under 110 % of a full window, but twice the time + // its window actually has. + final clean = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 14400); + final ts = beatEnds(clean, t0Ms: 1e12); + final startSec = (ts.first / 1000).floor(); + final tailNn = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 120); + final tailTs = beatEnds(tailNn, t0Ms: (startSec + 72 * 300 + 1) * 1000.0); + final rr = [...clean, for (final v in tailNn) ...[v, v]]; + final t = [...ts, for (final x in tailTs) ...[x, x]]; + final endSec = (t.last / 1000).ceil() + 1; + expect(endSec - (startSec + 72 * 300), lessThan(150)); + expect(rrCoverage(rr, t)!.overCounted, isFalse, + reason: 'the gap dilutes the night-wide ratio'); + final honest = sleepSessionRmssdDetail(clean, ts, + startSec: startSec, endSec: endSec); + final d = sleepSessionRmssdDetail(rr, t, + startSec: startSec, endSec: endSec); + expect(d.present, isTrue, reason: d.note); + expect(d.value!.overCountedWindows, 1); + expect(d.value!.windows, honest.value!.windows); + expect(d.value!.rmssd, honest.value!.rmssd); + // The same tail stored once is honest beat-time and stays in. + final once = sleepSessionRmssdDetail([...clean, ...tailNn], [...ts, ...tailTs], + startSec: startSec, endSec: endSec); + expect(once.value!.overCountedWindows, 0); + expect(once.value!.windows, honest.value!.windows + 1); + }); + + test('coverage does not depend on the order the beats arrive in', () { + final rr = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 1200); + final ts = beatEnds(rr, t0Ms: 1e12); + final sorted = rrCoverage(rr, ts)!; + final idx = [for (var i = 0; i < rr.length; i++) i] + ..shuffle(math.Random(1)); + final shuffled = rrCoverage( + [for (final i in idx) rr[i]], [for (final i in idx) ts[i]])!; + expect(shuffled.coverage, closeTo(sorted.coverage, 1e-12)); + expect(shuffled.spanSec, closeTo(sorted.spanSec, 1e-9)); + // Exact (ts, rr) repeats are counted wherever they sit. + final dupRr = [for (final v in rr) ...[v, v]]; + final dupTs = [for (final x in ts) ...[x, x]]; + final dIdx = [for (var i = 0; i < dupRr.length; i++) i] + ..shuffle(math.Random(2)); + final dupShuffled = rrCoverage( + [for (final i in dIdx) dupRr[i]], [for (final i in dIdx) dupTs[i]])!; + expect(rrCoverage(dupRr, dupTs)!.duplicateBeats, 1200); + expect(dupShuffled.duplicateBeats, 1200); + }); + + test('the nightly HRV shape is absent on an over-counted stream', () { + final nn = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30, beats: 25200); + final t = beatEnds(nn, t0Ms: 1e12); + expect(nightHrvShape(nn, t).present, isTrue); + final c = rrCoverage( + [for (final v in nn) ...[v, v]], [for (final x in t) ...[x, x]])!; + final s = nightHrvShape(nn, t, coverage: c); + expect(s.present, isFalse); + expect(s.note, startsWith('rr_overcount:coverage=')); + expect(nightHrvShape(nn, t, coverage: rrCoverage(nn, t)).present, isTrue); + }); + + test('no coverage handed over, or an honest one: no change', () { + final nn = rsaNn(hrBpm: 60, respBrpm: 12, ampMs: 30); + final t = beatEnds(nn, t0Ms: 1e12); + final c = rrCoverage(nn, t)!; + expect(c.overCounted, isFalse); + final a = hrvTime(nn, nnTimesMs: t); + final b = hrvTime(nn, nnTimesMs: t, coverage: c); + expect(b.value!.rmssd, a.value!.rmssd); + expect(b.confidence, a.confidence); + expect(b.note, a.note); + expect(nocturnalRmssd(nn, t, coverage: c).value, + nocturnalRmssd(nn, t).value); + }); + }); + group('frequency-domain HRV (Lomb-Scargle on beat times)', () { test('puts spectral peak at an injected RSA frequency in the HF band', () { // RR modulated at 0.25 Hz (respiration ~15 br/min) around 1000 ms. @@ -1248,6 +1551,8 @@ void main() { startSec: 1, endSec: 601, windowSec: 300, + // mechanics test: window size is not what is under test + minDiffsPerWindow: 1, ); expect(m.present, isTrue); expect(m.value, @@ -1266,7 +1571,9 @@ void main() { for (var i = 0; i < 10; i++) 1000.0 + i * 1000.0, 301000, 302000, ]; - final m = sleepSessionWindowedRmssd(rr, ts, startSec: 1, endSec: 601); + // The exclusion mechanics at a small floor; the default is tested below. + final m = sleepSessionWindowedRmssd(rr, ts, + startSec: 1, endSec: 601, minDiffsPerWindow: 5); expect(m.present, isTrue); expect(m.value, closeTo(10.0, 1e-9)); }); @@ -1274,7 +1581,9 @@ void main() { test('drops out-of-range and Malik-style ectopic beats before RMSSD', () { final rr = [1000, 1000, 1000, 1000, 200, 1000, 1000, 1000]; final ts = [for (var i = 0; i < 8; i++) 1000.0 + i * 1000.0]; - final m = sleepSessionWindowedRmssd(rr, ts, startSec: 1, endSec: 301); + // mechanics test: window size is not what is under test + final m = sleepSessionWindowedRmssd(rr, ts, + startSec: 1, endSec: 301, minDiffsPerWindow: 1); expect(m.present, isTrue); expect(m.value, closeTo(0.0, 1e-9)); }); @@ -1288,7 +1597,9 @@ void main() { // (87.7 -> 58.2, 82.9 -> 53.0, 76.9 -> 48.4 ms) and 2-13 % on gen4. final rr = [900, 900, 900, 900, 200, 1000, 1000, 1000]; final ts = [for (var i = 0; i < 8; i++) 1000.0 + i * 1000.0]; - final m = sleepSessionWindowedRmssd(rr, ts, startSec: 1, endSec: 301); + // mechanics test: window size is not what is under test + final m = sleepSessionWindowedRmssd(rr, ts, + startSec: 1, endSec: 301, minDiffsPerWindow: 1); expect(m.present, isTrue); expect(m.value, closeTo(0.0, 1e-9), reason: '2 runs of flat beats, no seam difference'); @@ -1323,18 +1634,121 @@ void main() { 1000, 2000, 3000, 4000, // flat run before the gap 123000, 124000, 125000, // flat run after a ~2 min dropout ]; - final m = sleepSessionWindowedRmssd(rr, ts, startSec: 1, endSec: 301); + // mechanics test: window size is not what is under test + final m = sleepSessionWindowedRmssd(rr, ts, + startSec: 1, endSec: 301, minDiffsPerWindow: 1); expect(m.present, isTrue); expect(m.value, closeTo(0.0, 1e-9), reason: 'two flat runs, no cross-gap difference manufactured'); }); + // Six full 5-min windows of a smooth oscillation (period 12 beats, ±16 ms; + // per-window RMSSD 2·16·sin(π/12)/√2 = 5.856 ms). 1800 beats are exactly + // 1 800 000 ms, so the last one ends at second 1800, alone in window 6. + const sixStart = 1000000000; // windows are (tsSec − sixStart) ~/ 300 + final sixRr = [ + for (var i = 0; i < 1800; i++) 1000.0 + 16.0 * math.sin(2 * math.pi * i / 12) + ]; + final sixTs = beatEnds(sixRr, t0Ms: sixStart * 1000.0); + + test('a 9-difference window with a 300 ms jump cannot move the headline', + () { + // Window 7 (seconds 2100–2399) holds 10 beats: five at 1000 ms, five at + // 1300 ms. Nine differences, one of them a 300 ms step across an + // arousal (window RMSSD 100 ms); every beat passes the range and Malik + // filters. Under a floor of 5 it counted as much as a full window. + final extra = [for (var i = 0; i < 10; i++) i < 5 ? 1000.0 : 1300.0]; + final rr2 = [...sixRr, ...extra]; + final ts2 = [...sixTs, ...beatEnds(extra, t0Ms: (sixStart + 2110) * 1000.0)]; + final base = sleepSessionWindowedRmssd(sixRr, sixTs, + startSec: sixStart, endSec: sixStart + 1800); + final thin = sleepSessionWindowedRmssd(rr2, ts2, + startSec: sixStart, endSec: sixStart + 2400); + expect(base.present && thin.present, isTrue); + expect(base.value!, closeTo(5.86, 0.05)); + expect(thin.value!, closeTo(base.value!, 1e-9), + reason: 'was (6 × 5.86 + 100) / 7 ≈ 19.3 under a floor of 5'); + final old = sleepSessionWindowedRmssd(rr2, ts2, + startSec: sixStart, endSec: sixStart + 2400, minDiffsPerWindow: 5); + expect(old.present, isTrue, reason: old.note); + expect(old.value!, closeTo((6 * base.value! + 100) / 7, 0.01)); + final d = sleepSessionRmssdDetail(rr2, ts2, + startSec: sixStart, endSec: sixStart + 2400); + expect(d.value!.windows, 6); + // Window 7, and window 6, whose one stray beat has no difference at all. + expect(d.value!.thinWindows, 2); + expect(d.value!.minDiffsPerWindow, kMinDiffsPerRmssdWindow); + expect(d.value!.toJson()['thin_windows'], 2); + expect(d.value!.toJson()['min_diffs_per_window'], 20); + expect(d.note, contains('2 window(s) had fewer than 20'), + reason: 'a published headline says what it dropped'); + }); + + test('a discarded thin window does not reach the pooled ACF1 verdict', () { + // The floor drops a thin window from the MEAN; it must drop it from the + // jitter screen's pooled differences too. A 20-beat window (19 + // differences, under the floor) alternating 910/1090 ms carries far + // more difference power than the six clean windows: pooled, it would + // drag ACF1 from ~0.866 to ~−0.78 and refuse a night it is not part of. + final extra = [for (var i = 0; i < 20; i++) i.isEven ? 910.0 : 1090.0]; + final rr2 = [...sixRr, ...extra]; + final ts2 = [...sixTs, ...beatEnds(extra, t0Ms: (sixStart + 2110) * 1000.0)]; + final base = sleepSessionRmssdDetail(sixRr, sixTs, + startSec: sixStart, endSec: sixStart + 2400); + final thin = sleepSessionRmssdDetail(rr2, ts2, + startSec: sixStart, endSec: sixStart + 2400); + expect(base.present && thin.present, isTrue, reason: thin.note); + expect(thin.value!.thinWindows, base.value!.thinWindows + 1); + expect(thin.value!.windows, base.value!.windows); + expect(thin.value!.rmssd, closeTo(base.value!.rmssd, 1e-9)); + expect(thin.value!.diffAcf1!, closeTo(0.863, 0.01)); + expect(thin.value!.diffAcf1!, closeTo(base.value!.diffAcf1!, 1e-12)); + expect(thin.confidence, base.confidence); + }); + + test('19 differences do not count, 20 do', () { + ({List rr, List ts}) window(int beats) { + final rr = [ + for (var i = 0; i < beats; i++) + 1000.0 + 16.0 * math.sin(2 * math.pi * i / 12) + ]; + return (rr: rr, ts: beatEnds(rr, t0Ms: 1000.0)); + } + + final a = window(21); // 20 Δ + final ma = sleepSessionRmssdDetail(a.rr, a.ts, startSec: 1, endSec: 301); + expect(ma.present, isTrue, reason: ma.note); + expect(ma.value!.windows, 1); + expect(ma.value!.thinWindows, 0); + final b = window(20); // 19 Δ + final mb = sleepSessionWindowedRmssd(b.rr, b.ts, startSec: 1, endSec: 301); + expect(mb.present, isFalse); + expect(mb.note, startsWith('no valid 5-min windows')); + expect(mb.note, contains('fewer than 20')); + }); + + test('the floor is a parameter', () { + final rr = [ + for (var i = 0; i < 20; i++) + 1000.0 + 16.0 * math.sin(2 * math.pi * i / 12) + ]; + final m = sleepSessionWindowedRmssd(rr, beatEnds(rr, t0Ms: 1000.0), + startSec: 1, endSec: 301, minDiffsPerWindow: 19); + expect(m.present, isTrue); + }); + + test('kMinDiffsPerRmssdWindow is 20', () { + expect(kMinDiffsPerRmssdWindow, 20); + }); + test('HRV-gap: no gap means no change (control)', () { // Same shape, but the second run starts right after the first (no real // gap) — the fix must not shrink a window that has nothing to exclude. final rr = [900, 900, 900, 1000, 1000, 1000]; final ts = [1000, 2000, 3000, 4000, 5000, 6000]; - final m = sleepSessionWindowedRmssd(rr, ts, startSec: 1, endSec: 301); + // mechanics test: window size is not what is under test + final m = sleepSessionWindowedRmssd(rr, ts, + startSec: 1, endSec: 301, minDiffsPerWindow: 1); expect(m.present, isTrue); // One real seam difference (900 -> 1000 = 100 ms) survives; RMSSD = 100 // over that single difference (the rest are 0).