Fix/numerical correctness

ALGORITHMS.md

@@ -61,7 +61,7 @@ Grouped by family (subdirectory under `lib/src/onehz/`). File paths are relative
 | `illnessCusum` | `clinical/illness_cusum.dart` | Online CUSUM state machine (green/yellow/red) over RHR — "NightSignal" | Alavi et al. 2022; Mishra et al. 2020 |
 | `readinessLnRmssd` | `clinical/readiness_lnrmssd.dart` | ln(RMSSD) z-scored against a rolling prior-nights baseline | Plews et al. 2013 |
 | `cosinor` | `clinical/cosinor.dart` | Cosinor rhythmometry (MESOR/amplitude/acrophase) | Halberg & Nelson 1979 |
-| `banisterTrimp` / `edwardsTrimp` | `clinical/load_trimp.dart` | Training impulse from HR-reserve | Banister 1991; Edwards 1993 |
+| `banisterTrimp` / `StrainScorer.edwardsTRIMP` | `clinical/load_trimp.dart` | Training impulse from HR-reserve | Banister 1991; Edwards 1993 |
 | `strainScoreMetric` | `clinical/load_trimp.dart` | log-squash of TRIMP onto a 0-21 scale | — |
 | `trimpStrain` | `clinical/load_trimp.dart` | TRIMP → 0-100 strain, honesty-wrapped (absent without real HRmax/RHR anchors) | — |
 | `ctlAtlTsb` | `clinical/load_trimp.dart` | Fitness-Fatigue-Form: EWMA CTL (42d) / ATL (7d) / TSB = CTL-ATL | Banister impulse-response model |

lib/onehz.dart

@@ -13,10 +13,11 @@ library onehz;
 // Layer: input types & math.
 export 'src/onehz/types.dart';
 export 'src/onehz/util.dart';
+// Per-device dispatch seam: unknown family refuses, never falls back to gen4.
+export 'src/onehz/device.dart';
 
 // Layer 0: foundations.
 export 'src/onehz/foundations/rr_correction.dart';
-export 'src/onehz/foundations/ppg_sqi.dart';
 export 'src/onehz/foundations/baseline.dart';
 export 'src/onehz/foundations/ewma_baselines.dart';
 export 'src/onehz/foundations/fusion.dart';

lib/src/onehz/clinical/hrv_freq.dart

@@ -6,9 +6,16 @@
 //   ULF < 0.003 Hz | VLF 0.003–0.04 | LF 0.04–0.15 | HF 0.15–0.40 | total
 // Normalized units: nu_lf = LF/(LF+HF)·100, nu_hf = HF/(LF+HF)·100.
 //
+// Band powers are WELCH-AVERAGED (Welch 1967): the record is split into
+// segments short enough that a modest frequency grid resolves them, each
+// segment's Lomb-Scargle PSD is integrated over the band on its own
+// resolution-matched grid, and the per-segment powers are averaged. See
+// [_welchBandPower] for why a single whole-night grid cannot work.
+//
 // HONESTY: HF is the band most corrupted by 1 Hz timing quantization and by
 // artifacts — we GATE HF (and LF/HF, nu) on the artifact fraction and report
-// reduced confidence. ULF needs a 24-h record; null on short reads.
+// reduced confidence. Every band is null unless the record is long enough to
+// RESOLVE it (ULF genuinely needs 24 h) — never a leakage-filled 0.0.
 
 import '../types.dart';
 import '../util.dart';
@@ -19,6 +26,11 @@ class HrvFreq {
   final double? lf;
   final double? hf;
   final double? total;
+
+  /// Which bands [total] is the sum of. Null exactly when [total] is null. A
+  /// band the record could not resolve is ABSENT from this list, never a 0 in
+  /// the sum — so "total" always says what it totalled.
+  final List<String>? totalBands;
   final double? lfhf;
   final double? nuLf;
   final double? nuHf;
@@ -29,6 +41,7 @@ class HrvFreq {
     this.lf,
     this.hf,
     this.total,
+    this.totalBands,
     this.lfhf,
     this.nuLf,
     this.nuHf,
@@ -40,6 +53,7 @@ class HrvFreq {
         if (lf != null) 'lf': round6(lf!),
         if (hf != null) 'hf': round6(hf!),
         if (total != null) 'total': round6(total!),
+        if (totalBands != null) 'total_bands': totalBands,
         if (lfhf != null) 'lf_hf': round6(lfhf!),
         if (nuLf != null) 'nu_lf': round6(nuLf!),
         if (nuHf != null) 'nu_hf': round6(nuHf!),
@@ -57,7 +71,7 @@ Metric<HrvFreq> hrvFreq(
   List<double> nnTimesMs, {
   required double artifactFraction,
   double hfArtifactGate = 0.15,
-  int gridPoints = 600,
+  double oversample = 4.0,
 }) {
   const inputs = ['rr_cleaned', 'beat_times'];
   if (nnMs.length < 16 || nnTimesMs.length != nnMs.length) {
@@ -79,28 +93,28 @@ Metric<HrvFreq> hrvFreq(
     );
   }
 
-  // Frequency grid: from ~1/span up to the HF ceiling (0.4 Hz).
-  final loHz = (1.0 / spanSec).clamp(0.0005, 0.04);
-  final ls = lombScargle(tSec, nnMs, freqGrid(loHz, 0.4, gridPoints));
-  if (ls == null) {
+  // Each band gets its own segment length and its own resolution-matched grid.
+  // A band whose lowest frequency the record is too short to resolve returns
+  // null — absent, not a leakage-filled 0.0 (ULF used to be emitted as exactly
+  // 0.0 for any session over 333 s, which is not a 24-h record by any reading).
+  final ulf =
+      _welchBandPower(tSec, nnMs, 0.0003, 0.003, oversample: oversample);
+  final vlf = _welchBandPower(tSec, nnMs, 0.003, 0.04, oversample: oversample);
+  final lf = _welchBandPower(tSec, nnMs, 0.04, 0.15, oversample: oversample);
+  final hfRaw = _welchBandPower(tSec, nnMs, 0.15, 0.40, oversample: oversample);
+  if (lf == null && hfRaw == null) {
     return const Metric<HrvFreq>.absent(
       tier: Tier.high,
       inputs_used: inputs,
-      note: 'spectrum undefined',
+      note: 'record too short to resolve any HRV band',
     );
   }
 
-  // Only report ULF/VLF if the record is long enough to resolve them.
-  final ulf = spanSec >= 1.0 / 0.003 ? ls.bandPower(0, 0.003) : null;
-  final vlf = spanSec >= 1.0 / 0.04 ? ls.bandPower(0.003, 0.04) : null;
-  final lf = ls.bandPower(0.04, 0.15);
-  final hfRaw = ls.bandPower(0.15, 0.40);
-
   final hfGated = artifactFraction > hfArtifactGate;
   final hf = hfGated ? null : hfRaw;
 
   double? lfhf, nuLf, nuHf;
-  if (hf != null && (lf + hf) > 0) {
+  if (lf != null && hf != null && (lf + hf) > 0) {
     lfhf = hf == 0 ? null : lf / hf;
     nuLf = 100.0 * lf / (lf + hf);
     nuHf = 100.0 * hf / (lf + hf);
@@ -112,7 +126,28 @@ Metric<HrvFreq> hrvFreq(
   // out bit-identical), and dropping HF from the sum would republish a
   // different quantity under the same name. Either way it would be dishonest —
   // so total is WITHHELD alongside HF.
-  final total = hfGated ? null : (ulf ?? 0) + (vlf ?? 0) + lf + hfRaw;
+  //
+  // A band the record could not RESOLVE is a different case, and it used to be
+  // handled with `?? 0` — which is exactly the absence-as-zero this package
+  // forbids, and it is not hypothetical: ULF needs 33 333 s of span (10 cycles
+  // at 0.0003 Hz) and a night is ~28 700 s, so ULF resolved on precisely no
+  // nights and the published total was silently the ULF-less sum. Withholding
+  // total instead would delete a rendered number on every night forever. So it
+  // is published with its composition NAMED in [totalBands]: the sum is over
+  // those bands and no others.
+  final bands = <String, double?>{
+    'ulf': ulf,
+    'vlf': vlf,
+    'lf': lf,
+    'hf': hfRaw
+  };
+  final resolved = [
+    for (final e in bands.entries)
+      if (e.value != null) e.key
+  ];
+  final total = (hfGated || lf == null || hfRaw == null)
+      ? null
+      : [for (final b in resolved) bands[b]!].reduce((a, b) => a + b);
 
   // Confidence: penalize artifacts heavily; low-band-only reads still HIGH-ish.
   final conf = clamp((1 - artifactFraction) * (hfGated ? 0.6 : 0.9), 0.2, 0.9);
@@ -123,6 +158,7 @@ Metric<HrvFreq> hrvFreq(
       lf: lf,
       hf: hf,
       total: total,
+      totalBands: total == null ? null : resolved,
       lfhf: lfhf,
       nuLf: nuLf,
       nuHf: nuHf,
@@ -137,3 +173,71 @@ Metric<HrvFreq> hrvFreq(
         : 'PRV spectrum; HF band quantization-limited at 1 Hz',
   );
 }
+
+/// Welch-averaged Lomb-Scargle power in [loHz, hiHz), or null when the record
+/// is too short to RESOLVE that band.
+///
+/// Why not one grid over the whole night: an 8 h record's periodogram has peaks
+/// ~1/28800 Hz wide, so a rectangular sum over a grid coarser than that is a
+/// lucky sample of the peaks, not an integral. The shipped 600-point grid put
+/// `lf_hf` at 0.095 on a synthetic whose converged value is 2.243 — a factor of
+/// 20+, on the one spectral number that is charted, stored and fed to the coach.
+/// The resolution-matched grid for a whole night is ~50k points and measured
+/// 18 s per night on this hardware; 600 points measured 1.2 s. Neither is
+/// acceptable: one is wrong, the other unshippable.
+///
+/// So, Welch 1967: split the record into segments SHORT enough that a modest
+/// grid resolves them, integrate each segment's PSD over the band on its own
+/// resolution-matched grid, and average. Segment length is [cyclesPerSegment]
+/// periods of the band's LOWEST frequency, so the band is resolved by
+/// construction; grid spacing is the segment's resolution (1/segment) divided
+/// by [oversample]. Cost is O(beats x gridPoints) per band and independent of
+/// how many segments the record splits into — measured ~1.4 s for a full night,
+/// i.e. no worse than the wrong version it replaces, with the periodogram's
+/// large variance averaged down as a bonus.
+///
+/// Averaging per-segment powers means the estimate excludes variance slower
+/// than one segment — correct by definition for a band whose lowest frequency
+/// sets the segment length, and the reason each band is segmented separately.
+double? _welchBandPower(
+  List<double> tSec,
+  List<double> y,
+  double loHz,
+  double hiHz, {
+  double cyclesPerSegment = 10.0,
+  double oversample = 4.0,
+  int minPointsPerSegment = 16,
+}) {
+  if (loHz <= 0 || hiHz <= loHz || tSec.length < minPointsPerSegment)
+    return null;
+  final span = tSec.last - tSec.first;
+  final segSec = cyclesPerSegment / loHz;
+  if (span < segSec) return null; // band not resolvable in this record
+
+  final df = 1.0 / segSec / oversample;
+  final nGrid = ((hiHz - loHz) / df).ceil() + 1;
+  final grid = freqGrid(loHz, hiHz, nGrid);
+
+  var sum = 0.0;
+  var k = 0;
+  final step = segSec / 2; // 50 % overlap, the Welch default
+  for (var start = tSec.first; start + segSec <= tSec.last; start += step) {
+    final end = start + segSec;
+    final ts = <double>[];
+    final ys = <double>[];
+    for (var i = 0; i < tSec.length; i++) {
+      if (tSec[i] < start) continue;
+      if (tSec[i] >= end) break;
+      ts.add(tSec[i]);
+      ys.add(y[i]);
+    }
+    if (ts.length < minPointsPerSegment) continue;
+    final ls = lombScargle(ts, ys, grid);
+    if (ls == null) continue;
+    final p = ls.bandPower(loHz, hiHz);
+    if (!p.isFinite) continue;
+    sum += p;
+    k++;
+  }
+  return k == 0 ? null : sum / k;
+}