From 8ad42194eeb699f77f2e6950318f884b54a14d04 Mon Sep 17 00:00:00 2001 From: abdulsaheel Date: Sun, 26 Jul 2026 18:01:47 +0530 Subject: [PATCH 1/2] fix: abstain instead of fabricating across the metric families MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit An audit of the whole package found a recurring shape: when an input was absent or a dispersion estimate was degenerate, the code substituted a constant instead of returning an absent Metric. Every fix below restores the package's core contract — absent input yields null, never a number — and each is covered by a regression test that fails against the old code. The two worst were in sleep, where a night with no data reported as a perfect night: - advanced_stager: a window with ZERO accelerometer samples returned one full-length 'light' segment, so a strap on the nightstand produced tst=28800 at 100% efficiency. The unbounded accel carry-forward did the same for a 6 h dropout. Windows are now split at long gaps, the carry-forward is bounded to 60 s, and unusable seconds stay wake. The carry-forward also copied the stale timestamp, mis-centring every RMSSD/LF-HF window inside a gap. - cardio_stager: with HR entirely absent the baseline fell back to a literal 60 bpm, after which every HR-relative gate was dead and every epoch defaulted to NREM. It now abstains. - Webster rescoring read the list it was mutating, so bridged wake runs inflated the context for the next one and cascaded; a fragmented night came back with waso=0 and 100% efficiency. Context is now snapshotted (same bug in stager.dart). Degenerate-dispersion fallbacks that fired instead of abstaining: - anomaly: a constant quantized column was clamped to a 1e-6 scale, so a 0.4 change became a ~1e6 z and flagged an illness anomaly every time. - illness_cusum: a max(1.0, SD) floor latched a sustained red from a 5 bpm one-night bump on a quantized baseline. - stress_si: a 1 ms RR range divided through to SI 48780 ('high'). - cpc: a zero or NaN-poisoned spectrum published a 999.0 sentinel, and a non-finite band power sailed past the variance guard entirely. Method and unit errors: - rr_correction: the Lipponen-Tarvainen threshold was taken over |dRR| rather than the signed series, collapsing the quartile deviation ~15x into a fixed 100 ms cutoff. An artifact-free 400-beat record had 32 beats flagged and spline-replaced. Verified against the real capture. - circadian_lifestyle: social jetlag medianed raw clock hours with no wrap, so a 23:50 weekday vs 01:10 weekend midsleep reported 22.55 h. Now circular throughout, with chronotype's bands on the same axis. - load_trimp: two disagreeing Banister implementations (neither with the published female scale), a sample duration taken from only the first two timestamps (8.08 vs 46.85 strain on the same data), an unclamped strain reaching 107.8, and CTL/ATL seeded from a single day. - nocturnal: the "lowest 30-min mean" never slid, silently becoming the mean of everything, and compacted across off-skin gaps. - hrv_freq: HF was withheld but still summed into total. - van_hees: the final win-1 seconds shared one verdict; now per-second, with an explicit undecidable tail rather than a guess. - resp_rate: brpm came from the median grid while peak_hz came from the highest-power grid, so peak_hz*60 disagreed with brpm. Honesty envelope: - physiologicalAge returned a present Metric with zero physiological inputs while claiming six in inputs_used. - workout detection let a missing HRmax DISABLE the zone-2 gate, then computed strain against a hidden 220-age anchor next to hrmax: null, and derived resting HR without filtering off-skin samples. - vo2maxEstimate emitted Infinity for restingHr == 0. - readiness_composite disclosed "robust-z" even when the mean/SD fallback produced the value. The fallback itself is unchanged. Also: gapAwareEwma extrapolated on non-positive dt, prsa threw RangeError for l=1, cosinor scored confidence on unadjusted R-squared, and cycles truncated pre-onset beats into minute 0. Outputs change. Consumers must bump their algorithm version to recompute. --- lib/src/onehz/clinical/cosinor.dart | 30 +- lib/src/onehz/clinical/hrv_freq.dart | 11 +- lib/src/onehz/clinical/illness_cusum.dart | 26 +- lib/src/onehz/clinical/load_trimp.dart | 175 ++++-- lib/src/onehz/clinical/nocturnal.dart | 80 ++- lib/src/onehz/clinical/prsa.dart | 11 + lib/src/onehz/clinical/readiness_lnrmssd.dart | 39 +- lib/src/onehz/clinical/stress_si.dart | 24 +- lib/src/onehz/foundations/baseline.dart | 9 + lib/src/onehz/foundations/rr_correction.dart | 18 +- lib/src/onehz/human/circadian_lifestyle.dart | 163 +++++- lib/src/onehz/human/coaching.dart | 139 ++++- lib/src/onehz/respiration/resp_rate.dart | 37 +- lib/src/onehz/sleep/advanced_stager.dart | 136 +++-- lib/src/onehz/sleep/cardio_stager.dart | 67 ++- lib/src/onehz/sleep/cpc.dart | 40 +- lib/src/onehz/sleep/cycles.dart | 7 +- lib/src/onehz/sleep/segment.dart | 58 +- lib/src/onehz/sleep/stager.dart | 25 +- lib/src/onehz/sleep/van_hees.dart | 82 ++- lib/src/onehz/wellness/anomaly.dart | 64 +- .../onehz/wellness/readiness_composite.dart | 12 +- lib/src/onehz/workout/workout_detect.dart | 111 +++- test/onehz/clinical_test.dart | 309 +++++++++- test/onehz/coaching_test.dart | 191 ++++++ test/onehz/foundations_test.dart | 75 +++ test/onehz/human_test.dart | 59 ++ test/onehz/respiration_test.dart | 36 ++ test/onehz/sleep_honesty_test.dart | 546 ++++++++++++++++++ test/onehz/wellness_test.dart | 80 +++ test/onehz/workout_test.dart | 218 +++++++ 31 files changed, 2596 insertions(+), 282 deletions(-) create mode 100644 test/onehz/sleep_honesty_test.dart diff --git a/lib/src/onehz/clinical/cosinor.dart b/lib/src/onehz/clinical/cosinor.dart index 64948a5..e3a25ab 100644 --- a/lib/src/onehz/clinical/cosinor.dart +++ b/lib/src/onehz/clinical/cosinor.dart @@ -20,7 +20,12 @@ class CosinorFit { final double acrophaseRad; // phase of the peak (rad), in [-π, π] final double acrophaseHours; // acrophase expressed as clock-hours of [period] final double periodHours; - final double r2; // goodness of fit (0..1) + final double r2; // raw goodness of fit (0..1) + + /// R² adjusted for the 3 fitted parameters (M, β, γ): + /// R²adj = 1 − (1−R²)·(n−1)/(n−3). This is the honest fit quality — the raw + /// R² of a 3-parameter fit is upward-biased and approaches 1 as n → 3. + final double r2Adj; const CosinorFit({ required this.mesor, required this.amplitude, @@ -28,6 +33,7 @@ class CosinorFit { required this.acrophaseHours, required this.periodHours, required this.r2, + required this.r2Adj, }); Map toJson() => { 'mesor': round6(mesor), @@ -36,21 +42,31 @@ class CosinorFit { 'acrophase_hours': round6(acrophaseHours), 'period_hours': round6(periodHours), 'r2': round6(r2), + 'r2_adj': round6(r2Adj), }; } +/// Minimum samples for a cosinor fit. The model has 3 free parameters +/// (M, β, γ); with n = 4 there is a single residual degree of freedom, so R² +/// is essentially an interpolation score (4 RANDOM points routinely fit at +/// r² 0.76–0.99). Halberg's zero-amplitude test needs real residual dof, so we +/// require ≥8 samples over the period and report the ADJUSTED R². +const int cosinorMinPoints = 8; + /// Single-component cosinor fit. /// /// [tHours] sample times in hours (any origin). [y] sample values. [periodHours] -/// the rhythm period (default 24). Needs ≥4 points and non-degenerate design. +/// the rhythm period (default 24). Needs ≥[minPoints] points and a +/// non-degenerate design. Metric cosinor( List tHours, List y, { double periodHours = 24, + int minPoints = cosinorMinPoints, }) { const inputs = ['signal_timeseries']; final n = y.length; - if (n < 4 || tHours.length != n || periodHours <= 0) { + if (n < math.max(4, minPoints) || tHours.length != n || periodHours <= 0) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, @@ -109,8 +125,13 @@ Metric cosinor( ssRes += (y[i] - fit) * (y[i] - fit); } final r2 = ssTot == 0 ? 0.0 : clamp(1 - ssRes / ssTot, 0, 1); + // Adjusted for the 3 fitted parameters (M, β, γ). Confidence MUST come from + // the adjusted value: the raw R² of a 3-parameter fit is upward-biased + // (E[R²] = 2/(n−1) under the null), so a handful of noise points used to + // score confidence 0.95 at tier HIGH. + final r2Adj = clamp(1 - (1 - r2) * (n - 1) / (n - 3), 0, 1); - final conf = clamp(r2, 0.1, 0.95); + final conf = clamp(r2Adj, 0.1, 0.95); return Metric( value: CosinorFit( mesor: mesor, @@ -119,6 +140,7 @@ Metric cosinor( acrophaseHours: phaseHours, periodHours: periodHours, r2: r2, + r2Adj: r2Adj, ), confidence: conf, tier: Tier.high, diff --git a/lib/src/onehz/clinical/hrv_freq.dart b/lib/src/onehz/clinical/hrv_freq.dart index d25f48c..b9ebf09 100644 --- a/lib/src/onehz/clinical/hrv_freq.dart +++ b/lib/src/onehz/clinical/hrv_freq.dart @@ -105,7 +105,14 @@ Metric hrvFreq( nuLf = 100.0 * lf / (lf + hf); nuHf = 100.0 * hf / (lf + hf); } - final total = (ulf ?? 0) + (vlf ?? 0) + lf + hfRaw; + // TOTAL POWER is by definition the sum over ALL bands (Task Force 1996). + // When HF is suppressed we do not have a trustworthy HF term, so a "total" + // is not computable: silently summing the withheld hfRaw back in republished + // exactly the quantity the gate withheld (the gated and ungated totals came + // 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; // 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); @@ -126,7 +133,7 @@ Metric hrvFreq( inputs_used: inputs, note: hfGated ? 'HF suppressed: artifact fraction ${round6(artifactFraction)} ' - '> gate — LF/VLF reported, HF/LF-HF/nu withheld' + '> gate — LF/VLF reported, HF/LF-HF/nu/total withheld' : 'PRV spectrum; HF band quantization-limited at 1 Hz', ); } diff --git a/lib/src/onehz/clinical/illness_cusum.dart b/lib/src/onehz/clinical/illness_cusum.dart index 01182a7..3b115b0 100644 --- a/lib/src/onehz/clinical/illness_cusum.dart +++ b/lib/src/onehz/clinical/illness_cusum.dart @@ -22,6 +22,13 @@ enum IllnessState { green, yellow, red } /// Required minimum valid baseline nights before the CUSUM can flag. const int illnessCusumMinBaseline = 7; +/// Machine-readable note attached to a night whose trailing baseline is long +/// enough but has ZERO dispersion (MAD and SD both 0 — a fully constant, +/// quantized baseline). The standardized deviation is undefined, so the night +/// is held green and NOT accumulated rather than standardized against a +/// fabricated scale. +const String degenerateBaselineNote = 'degenerate_baseline:scale=0'; + class IllnessDay { final String date; final IllnessState state; @@ -91,10 +98,21 @@ List illnessCusum( final med = median(window)!; var scale = mad(window) ?? 0; if (scale <= 0) { - // Quantized/constant baseline: fall back to a small physiological floor - // (1 bpm) so we can still standardize, but flag low confidence by never - // letting tiny noise trip the alarm. - scale = math.max(1.0, (stddev(window) ?? 1.0)); + // Quantized baseline (whole-bpm RHR) can collapse the MAD to 0. Fall + // back to the ordinary SD so a usable baseline still standardizes — + // the same convention as wellness/readiness_composite.dart and + // wellness/changepoint.dart. + scale = stddev(window) ?? 0; + } + if (scale <= 0 || !scale.isFinite) { + // Truly constant baseline: there is NO dispersion to standardize + // against, so z is undefined. Substituting a magic 1 bpm floor here + // turned a 5 bpm one-night bump into z = 5 and latched the alarm red — + // exactly the fabrication this package forbids. ABSTAIN instead: hold + // green, do not accumulate, and say why. + out.add(IllnessDay(dates[i], IllnessState.green, null, null, + need: degenerateBaselineNote)); + continue; } final z = (r - med) / scale; // One-sided upper CUSUM on elevation (RHR up = potential illness). diff --git a/lib/src/onehz/clinical/load_trimp.dart b/lib/src/onehz/clinical/load_trimp.dart index 2a9f02e..96413e6 100644 --- a/lib/src/onehz/clinical/load_trimp.dart +++ b/lib/src/onehz/clinical/load_trimp.dart @@ -17,8 +17,16 @@ import '../util.dart'; /// Banister TRIMP over a series of per-minute mean HRs. /// +/// TRIMP = Σ Δt(min) · ΔHRr · y(ΔHRr) with the PUBLISHED sex-specific weighting +/// factor (Banister 1991; Morton 1990): +/// men y = 0.64 · e^(1.92·x) +/// women y = 0.86 · e^(1.67·x) +/// Delegates to [StrainScorer.banisterY] so this file has exactly ONE Banister +/// implementation (it previously dropped the 0.64/0.86 coefficient entirely, +/// disagreeing with [StrainScorer.banisterTRIMP] by a factor of 1.5625). +/// /// [hrPerMin] mean HR for each worn minute (bpm; pass only valid minutes). -/// [restingHr], [maxHr] the personal anchors. [sex] selects the b constant. +/// [restingHr], [maxHr] the personal anchors. [sex] selects the coefficients. /// Returns absent if anchors are missing/degenerate (no fabrication). Metric banisterTrimp( List hrPerMin, { @@ -37,7 +45,6 @@ Metric banisterTrimp( note: 'Banister TRIMP needs measured RHR and HRmax (HRmax>RHR)', ); } - final b = sex == Sex.male ? 1.92 : 1.67; final reserve = maxHr - restingHr; var trimp = 0.0; for (final hr in hrPerMin) { @@ -45,7 +52,9 @@ Metric banisterTrimp( var hrr = (hr - restingHr) / reserve; if (hrr < 0) hrr = 0; if (hrr > 1) hrr = 1; - trimp += 1.0 * hrr * math.exp(b * hrr); // 1 minute each + // ONE Banister implementation for the whole package: the sex-specific + // weighting factor y lives in [StrainScorer.banisterY]. 1 minute each. + trimp += 1.0 * hrr * StrainScorer.banisterY(hrr, female: sex == Sex.female); } return Metric( value: trimp, @@ -125,14 +134,17 @@ Metric edwardsTrimp(List zoneMinutes) { // // 1. Heart-Rate Reserve (Karvonen): HRR = HRmax − RHR. // 2. Per-sample intensity %HRR = (HR − RHR) / HRR × 100, clamped 0..100. -// 3. TRIMP over the window: +// 3. TRIMP over the window (each sample carries its OWN duration, measured +// from the real timestamps and capped at the stream's median cadence so a +// gap is never counted as effort): // a. Edwards 5-zone (default): sample contributes its zone weight (1..5 at -// 50/60/70/80/90 %HRR cut-offs) × per-sample duration (min). -// b. Banister exponential: sample contributes dur × x × 0.64 × e^(b·x). +// 50/60/70/80/90 %HRR cut-offs) × that sample's duration (min). +// b. Banister exponential: sample contributes dur × x × y(x), with +// y = 0.64·e^(1.92x) (men) / 0.86·e^(1.67x) (women). // 4. strain = 100 × ln(TRIMP + 1) / ln(D), D = 7201 (TRIMP 7200 ≈ max). // -// References: Karvonen 1957; Edwards 1993; Banister 1991 (b = 1.92 M / 1.67 F); -// Tanaka 2001 (HRmax = 208 − 0.7·age). +// References: Karvonen 1957; Edwards 1993; Banister 1991 (y = 0.64·e^(1.92x) +// men / 0.86·e^(1.67x) women); Tanaka 2001 (HRmax = 208 − 0.7·age). // // NOTE (steps/active-energy floor): strain is PURELY HR-derived // (Edwards/Banister TRIMP → log map). Steps and active calories are computed as @@ -171,11 +183,27 @@ class StrainScorer { /// Upper percentile for the observed-HRmax estimate. static const double hrmaxPercentile = 99.5; - /// Banister coefficients. - static const double banisterScale = 0.64; + /// Banister 1991 weighting factor y = c · e^(b·x), x = fractional %HRR. + /// PUBLISHED coefficients are sex-specific in BOTH terms: + /// men c = 0.64, b = 1.92 + /// women c = 0.86, b = 1.67 + /// (Applying the male c = 0.64 to women — as this class used to — understates + /// female TRIMP by ~26 %.) + static const double banisterScaleMen = 0.64; + static const double banisterScaleWomen = 0.86; static const double banisterBMen = 1.92; static const double banisterBWomen = 1.67; + /// Deprecated alias for [banisterScaleMen]; kept so existing call sites keep + /// resolving. Prefer [banisterY], which pairs c and b correctly by sex. + static const double banisterScale = banisterScaleMen; + + /// The sex-specific Banister weighting factor y(x) — the single source of + /// truth for Banister weighting in this package. + static double banisterY(double x, {required bool female}) => + (female ? banisterScaleWomen : banisterScaleMen) * + math.exp((female ? banisterBWomen : banisterBMen) * x); + /// Edwards zone cut-offs as (%HRR threshold, weight), highest-first. static const List> edwardsZones = [ [90.0, 5], @@ -242,29 +270,70 @@ class StrainScorer { // ── TRIMP accumulation ────────────────────────────────────────────────────── - /// Per-sample duration (minutes) from the first two timestamps (seconds). - /// Falls back to 1 s when <2 samples or coincident timestamps. - static double sampleDurationMinutes(List tsSec) { - if (tsSec.length < 2) return fallbackSampleMin; - final deltaS = (tsSec[1] - tsSec[0]).abs(); - return deltaS > 0 ? deltaS / 60.0 : fallbackSampleMin; + /// Median inter-sample interval (seconds) of a time-ordered stream, ignoring + /// non-positive steps and pathological (>[maxPlausibleGapSec]) ones. Floored + /// at [fallbackSampleMin] minutes' worth. Mirrors the convention already used + /// by `HeartRateZones.timeInZone`. + static double medianIntervalSeconds(List tsSec, + {double maxPlausibleGapSec = 300.0}) { + final gaps = []; + for (var i = 1; i < tsSec.length; i++) { + final g = tsSec[i] - tsSec[i - 1]; + if (g > 0 && g <= maxPlausibleGapSec) gaps.add(g); + } + if (gaps.isEmpty) return fallbackSampleMin * 60.0; + gaps.sort(); + return math.max(gaps[gaps.length ~/ 2], fallbackSampleMin * 60.0); } - static double edwardsTRIMP(List bpm, double restingHR, double hrReserve, - double sampleDurationMin) { - var weighted = 0; - for (final s in bpm) { - weighted += zoneWeight(s, restingHR, hrReserve); + /// PER-SAMPLE effort durations (minutes) from the ACTUAL timestamps. + /// + /// Sample i is credited with the interval to sample i+1; the tail sample gets + /// the stream's median cadence. Every interval is CAPPED at that median, so a + /// hole in the stream can never be counted as sustained effort (the same gap + /// policy as `HeartRateZones.timeInZone`). + /// + /// This replaces the old `sampleDurationMinutes`, which read ONE interval + /// (the first two timestamps) and applied it to every sample — catastrophic + /// on exactly the sparse/irregular streams [minSparseReadings] admits: 21 + /// samples over 20 min with the first two 1 s apart scored strain 8.08 + /// instead of ~47, and a 1 Hz stream with a 5-min leading gap scored 104. + static List sampleDurationsMinutes(List tsSec) { + final n = tsSec.length; + if (n == 0) return const []; + if (n == 1) return [fallbackSampleMin]; + final capSec = medianIntervalSeconds(tsSec); + final out = List.filled(n, capSec / 60.0); + for (var i = 0; i < n - 1; i++) { + final g = tsSec[i + 1] - tsSec[i]; + out[i] = (g > 0 ? math.min(g, capSec) : capSec) / 60.0; } - return weighted * sampleDurationMin; + return out; } - static double banisterTRIMP(List bpm, double restingHR, double hrReserve, - double sampleDurationMin, double b) { + /// Edwards 5-zone TRIMP: Σ zoneWeight(sample) × that sample's duration (min). + static double edwardsTRIMP(List bpm, double restingHR, + double hrReserve, List durationsMin) { var acc = 0.0; - for (final s in bpm) { - final x = pctHRR(s, restingHR, hrReserve) / 100.0; - if (x > 0) acc += sampleDurationMin * x * banisterScale * math.exp(b * x); + for (var i = 0; i < bpm.length; i++) { + final dur = i < durationsMin.length + ? durationsMin[i] + : (durationsMin.isEmpty ? fallbackSampleMin : durationsMin.last); + acc += zoneWeight(bpm[i], restingHR, hrReserve) * dur; + } + return acc; + } + + /// Banister exponential TRIMP: Σ duration(min) × x × y(x), y per [banisterY]. + static double banisterTRIMP(List bpm, double restingHR, + double hrReserve, List durationsMin, {bool female = false}) { + var acc = 0.0; + for (var i = 0; i < bpm.length; i++) { + final dur = i < durationsMin.length + ? durationsMin[i] + : (durationsMin.isEmpty ? fallbackSampleMin : durationsMin.last); + final x = pctHRR(bpm[i], restingHR, hrReserve) / 100.0; + if (x > 0) acc += dur * x * banisterY(x, female: female); } return acc; } @@ -272,11 +341,14 @@ class StrainScorer { // ── Logarithmic map ───────────────────────────────────────────────────────── /// Map accumulated TRIMP onto [0, 100] via 100 × ln(TRIMP+1) / ln(D), 2 dp. - /// TRIMP ≤ 0 → 0. + /// TRIMP ≤ 0 → 0; above the D−1 ceiling the score is CLAMPED at [maxStrain] + /// (it used to run off the top of its own documented range: TRIMP 14400 → + /// 107.8, while the sibling [strainScore] clamped correctly). static double trimpToStrain(double trimp, {double denominator = strainDenominator}) { if (trimp <= 0) return 0; final value = maxStrain * math.log(trimp + 1.0) / math.log(denominator); - return (value * 100).roundToDouble() / 100; + final clamped = math.min(maxStrain, math.max(0.0, value)); + return (clamped * 100).roundToDouble() / 100; } // ── TRIMP method ────────────────────────────────────────────────────────────── @@ -316,15 +388,15 @@ class StrainScorer { } if (!enoughData || effMax <= restingHR) return null; - final sampleDur = sampleDurationMinutes(tsSec); + final durations = sampleDurationsMinutes(tsSec); final hrReserve = effMax - restingHR; final double trimp; if (edwards) { - trimp = edwardsTRIMP(bpm, restingHR, hrReserve, sampleDur); + trimp = edwardsTRIMP(bpm, restingHR, hrReserve, durations); } else { - final b = female ? banisterBWomen : banisterBMen; - trimp = banisterTRIMP(bpm, restingHR, hrReserve, sampleDur, b); + trimp = banisterTRIMP(bpm, restingHR, hrReserve, durations, + female: female); } return trimpToStrain(trimp, denominator: denominator); } @@ -396,24 +468,44 @@ class LoadState { }; } +/// Minimum days of daily-TRIMP history before CTL/ATL/TSB are reported. +/// Two weeks: enough for the 7-day ATL to be converged and for the CTL prime +/// to rest on a real week of load rather than a single day. +const int ctlAtlMinDays = 14; + /// CTL/ATL/TSB from a time-ordered daily-TRIMP series (oldest→newest). /// EWMA with time constants 42 d (CTL) and 7 d (ATL): λ = 1 − e^(−1/τ). /// A missing day contributes a 0-load impulse (rest day) — the EWMA decays. +/// +/// SEEDING (Banister 1975 impulse-response; the load before the record started +/// is UNKNOWN): both accumulators used to be seeded at `dailyTrimp.first`, +/// which asserted that a single observed day had already been sustained for the +/// full 42-day chronic window — `ctlAtlTsb([500])` returned ctl 500 / atl 500 / +/// tsb 0, a fully-adapted, perfectly-fresh athlete conjured from one workout. +/// Now: ABSTAIN below [minDays] with the standard need_baseline note, and prime +/// both accumulators with the MEAN of the first [primeDays] observed days +/// (never a single day, never future days) before running the EWMA over the +/// remainder. Metric ctlAtlTsb(List dailyTrimp, - {double ctlDays = 42, double atlDays = 7}) { + {double ctlDays = 42, + double atlDays = 7, + int minDays = ctlAtlMinDays, + int primeDays = 7}) { const inputs = ['daily_trimp']; - if (dailyTrimp.isEmpty) { - return const Metric.absent( + if (dailyTrimp.length < minDays) { + return Metric.absent( tier: Tier.estimate, inputs_used: inputs, - note: 'no daily TRIMP history', + note: needBaselineNote(have: dailyTrimp.length, need: minDays), ); } final lc = 1 - math.exp(-1 / ctlDays); final la = 1 - math.exp(-1 / atlDays); - var ctl = dailyTrimp.first; - var atl = dailyTrimp.first; - for (var i = 1; i < dailyTrimp.length; i++) { + final prime = math.min(math.max(primeDays, 1), dailyTrimp.length); + final seed = mean(dailyTrimp.sublist(0, prime))!; + var ctl = seed; + var atl = seed; + for (var i = prime; i < dailyTrimp.length; i++) { ctl = ctl + lc * (dailyTrimp[i] - ctl); atl = atl + la * (dailyTrimp[i] - atl); } @@ -423,6 +515,7 @@ Metric ctlAtlTsb(List dailyTrimp, confidence: conf, tier: Tier.estimate, inputs_used: inputs, - note: 'Banister CTL(42d)/ATL(7d)/TSB; descriptive load, not injury risk', + note: 'Banister CTL(42d)/ATL(7d)/TSB, primed from the first $prime observed ' + 'days; descriptive load, not injury risk', ); } diff --git a/lib/src/onehz/clinical/nocturnal.dart b/lib/src/onehz/clinical/nocturnal.dart index 162ec85..9cf46ec 100644 --- a/lib/src/onehz/clinical/nocturnal.dart +++ b/lib/src/onehz/clinical/nocturnal.dart @@ -24,29 +24,63 @@ class NocturnalRhr { /// Nocturnal resting HR from a night of 1 Hz HR samples. /// -/// [hr] 1 Hz HR samples (bpm; 0 = off-skin, excluded). [windowSec] rolling mean -/// window (default 30 min). Assumes ~1 Hz spacing; uses a sample-count window. -Metric nocturnalRhr(List hr, {int windowSamples = 1800}) { +/// [hr] 1 Hz HR samples (bpm; 0 = off-skin, excluded). [windowSamples] rolling +/// window length in SAMPLE POSITIONS (default 1800 = 30 min at 1 Hz). +/// [minCoverage] fraction of a window's positions that must carry a valid +/// on-skin sample for the window to count. +/// +/// The window slides over WALL-CLOCK POSITIONS, not over the compacted valid +/// stream: an off-skin gap must not be closed up, or the "30-min" window can +/// silently span the whole night. A night with no window meeting [minCoverage] +/// yields an ABSENT metric — we never relabel the whole-night mean as a +/// lowest-30-min trough. +Metric nocturnalRhr(List hr, + {int windowSamples = 1800, double minCoverage = 0.9}) { const inputs = ['hr_1hz']; final valid = hr.where((h) => h > 0).toList(); - if (valid.length < windowSamples ~/ 2) { + if (windowSamples < 1 || hr.length < windowSamples) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'insufficient valid (on-skin) HR for nocturnal RHR', ); } - // Lowest rolling mean over the valid stream. - final w = windowSamples > valid.length ? valid.length : windowSamples; + // Lowest rolling mean over CONTIGUOUS wall-clock windows. A window is only + // eligible when at least [minCoverage] of its positions are on-skin; its + // mean is taken over the valid samples inside it. + final needValid = (minCoverage * windowSamples).ceil(); var sum = 0.0; - for (var i = 0; i < w; i++) { - sum += valid[i]; + var count = 0; + for (var i = 0; i < windowSamples; i++) { + if (hr[i] > 0) { + sum += hr[i]; + count++; + } + } + double? best; + if (count >= needValid) best = sum / count; + for (var i = windowSamples; i < hr.length; i++) { + if (hr[i] > 0) { + sum += hr[i]; + count++; + } + final out = hr[i - windowSamples]; + if (out > 0) { + sum -= out; + count--; + } + if (count >= needValid) { + final m = sum / count; + if (best == null || m < best) best = m; + } } - var best = sum / w; - for (var i = w; i < valid.length; i++) { - sum += valid[i] - valid[i - w]; - final m = sum / w; - if (m < best) best = m; + if (best == null) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'no contiguous on-skin window long enough for a nocturnal RHR ' + 'trough (off-skin gaps are never compacted away)', + ); } final p1 = percentile(valid, 1)!; final conf = clamp(valid.length / 7200.0, 0.4, 0.95); // ~2 h coverage => high @@ -73,13 +107,31 @@ class HrDip { }; } +/// Minimum valid on-skin samples required on EACH side (day, night) before a +/// dip % is reported: 300 samples ≈ 5 min at 1 Hz. Below that a "daytime mean" +/// and a "nocturnal mean" are single moments, not periods. +const int hrDipMinSamples = 300; + /// Nocturnal HR dip %. [dayHr] and [nightHr] are 1 Hz HR samples for the waking /// and sleeping periods respectively (0 excluded). Bands follow the BP-dip /// convention applied to HR: ≥10% dipper, 0–10% non-dipper, <0 riser. -Metric hrDip(List dayHr, List nightHr) { +/// +/// Both sides need at least [minSamples] valid samples — one day sample and one +/// night sample can produce any dip % at all, so a dip band computed from a +/// handful of samples is fabrication, not measurement. +Metric hrDip(List dayHr, List nightHr, + {int minSamples = hrDipMinSamples}) { const inputs = ['hr_1hz_day', 'hr_1hz_night']; final dv = dayHr.where((h) => h > 0).toList(); final nv = nightHr.where((h) => h > 0).toList(); + if (dv.length < minSamples || nv.length < minSamples) { + return Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'HR dip needs ≥$minSamples valid samples on each of day and night ' + '(have day=${dv.length}, night=${nv.length})', + ); + } final dm = mean(dv); final nm = mean(nv); if (dm == null || nm == null || dm <= 0) { diff --git a/lib/src/onehz/clinical/prsa.dart b/lib/src/onehz/clinical/prsa.dart index 730bb1b..65cf415 100644 --- a/lib/src/onehz/clinical/prsa.dart +++ b/lib/src/onehz/clinical/prsa.dart @@ -68,6 +68,17 @@ Metric _prsa( final kind = deceleration ? 'DC' : 'AC'; final inputs = const ['rr_cleaned']; final n = nnMs.length; + if (l < 2) { + // Bauer 2006 quantifies PRSA with the Haar contrast at wavelet scale s = 2: + // DC = [X(0) + X(1) − X(−1) − X(−2)]/4, which needs TWO profile points on + // each side of the anchor. With l < 2 the profile is too short (l = 1 used + // to index profile[l−2] = profile[−1] and throw a RangeError). + return Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'PRSA ($kind) needs l≥2 for the s=2 Haar contrast', + ); + } if (n < 2 * l + 4) { return Metric.absent( tier: Tier.high, diff --git a/lib/src/onehz/clinical/readiness_lnrmssd.dart b/lib/src/onehz/clinical/readiness_lnrmssd.dart index c14e26b..d2696cb 100644 --- a/lib/src/onehz/clinical/readiness_lnrmssd.dart +++ b/lib/src/onehz/clinical/readiness_lnrmssd.dart @@ -82,29 +82,38 @@ Metric readinessLnRmssd( } final m = mean(priorWindow)!; final sd = stddev(priorWindow); - final cv = (m != 0 && sd != null) ? (sd / m).abs() * 100 : 0.0; - final z = (sd != null && sd > 0) ? (today - m) / sd : null; + if (sd == null || m == 0) { + // The published outputs of this stack (CV, SWC, band) are all defined + // RELATIVE to the baseline's dispersion. With a single prior night the SD + // is UNDEFINED — and the metric used to fill in cvPct 0.0 and band + // 'normal', asserting "tonight is typical" on the strength of nothing at + // all. Abstain instead (a zero-but-DEFINED SD is a different case: CV is + // genuinely 0 and SWC is genuinely 0, so that still computes below). + return Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'lnRMSSD baseline dispersion undefined (needs ≥2 prior nights) — ' + 'CV/z/SWC/band are undefined', + ); + } + final cv = (sd / m).abs() * 100; + final z = sd > 0 ? (today - m) / sd : null; // Plews SWC ≈ 0.5 × within-window SD (a small worthwhile change in lnRMSSD). - final swc = sd != null ? 0.5 * sd : null; + final swc = 0.5 * sd; - String band; - if (swc != null) { - if (today < m - swc) { - band = 'suppressed'; - } else if (today > m + swc) { - band = 'elevated'; - } else { - band = 'normal'; - } + final String band; + if (today < m - swc) { + band = 'suppressed'; + } else if (today > m + swc) { + band = 'elevated'; } else { band = 'normal'; } // LnRMSSD:RR saturation guard: when mean NN is long (low HR, high vagal tone) // and lnRMSSD is at the top of the personal range, the metric saturates. - final saturation = meanNnTodayMs != null && - meanNnTodayMs > 1100 && - (sd != null && sd > 0 && (today - m) / sd > 1.0); + final saturation = + meanNnTodayMs != null && meanNnTodayMs > 1100 && z != null && z > 1.0; final conf = clamp(priorWindow.length / windowDays.toDouble(), 0.3, 0.9); return Metric( diff --git a/lib/src/onehz/clinical/stress_si.dart b/lib/src/onehz/clinical/stress_si.dart index fdd4df3..3073ee2 100644 --- a/lib/src/onehz/clinical/stress_si.dart +++ b/lib/src/onehz/clinical/stress_si.dart @@ -45,7 +45,15 @@ class StressIndex { /// variation range (MxDMn) is maximal, which drives SI to ~0 (a bug we hit on /// real data: whole-night SI read ~7). So we slide ~5-min windows (~256 beats), /// compute SI per window, and report the MEDIAN — the robust resting SI. -Metric baevskyStressIndex(List nnMs) { +/// +/// [minRangeMs] guards the MxDMn denominator. Baevsky's SI is defined on an RR +/// histogram whose variation range is a physiological quantity (0.15–0.5 s at +/// rest); a 5-min window whose whole RR range is a few ms is beat-timing +/// QUANTIZATION, not physiology, and 1/MxDMn then explodes (300 beats +/// alternating 1000/1001 ms → SI ≈ 48 780, reported as 'high'). Windows below +/// the guard are dropped; if no window survives, the metric is ABSENT. +Metric baevskyStressIndex(List nnMs, + {double minRangeMs = 20.0}) { const inputs = ['rr_cleaned']; final nn = nnMs.where((v) => v >= 300 && v <= 2000).toList(); if (nn.length < 30) { @@ -63,7 +71,7 @@ Metric baevskyStressIndex(List nnMs) { for (var start = 0; start + 30 <= nn.length; start += step) { final seg = nn.sublist(start, math.min(start + win, nn.length)); if (seg.length < 30) break; - final r = _siOfSegment(seg); + final r = _siOfSegment(seg, minRangeMs); if (r != null) { sis.add(r[0]); modes.add(r[1]); @@ -104,8 +112,9 @@ Metric baevskyStressIndex(List nnMs) { ); } -/// SI of one short NN segment → [si, modeS, amoPct, mxdmnS], or null if degenerate. -List? _siOfSegment(List seg) { +/// SI of one short NN segment → [si, modeS, amoPct, mxdmnS], or null if +/// degenerate (mode undefined, or a variation range below [minRangeMs]). +List? _siOfSegment(List seg, double minRangeMs) { const binMs = 50.0; final counts = {}; for (final v in seg) { @@ -121,7 +130,10 @@ List? _siOfSegment(List seg) { }); final modeS = ((modeBin + 0.5) * binMs) / 1000.0; final amoPct = 100.0 * modeCount / seg.length; - final mxdmnS = (seg.reduce(math.max) - seg.reduce(math.min)) / 1000.0; - if (modeS <= 0 || mxdmnS <= 0) return null; + final mxdmnMs = seg.reduce(math.max) - seg.reduce(math.min); + final mxdmnS = mxdmnMs / 1000.0; + // MxDMn sits in the denominator: a range at or below the beat-timing + // resolution is not a measurement of autonomic tension, it is quantization. + if (modeS <= 0 || mxdmnMs < minRangeMs) return null; return [amoPct / (2.0 * modeS * mxdmnS), modeS, amoPct, mxdmnS]; } diff --git a/lib/src/onehz/foundations/baseline.dart b/lib/src/onehz/foundations/baseline.dart index 022329b..adc0343 100644 --- a/lib/src/onehz/foundations/baseline.dart +++ b/lib/src/onehz/foundations/baseline.dart @@ -97,6 +97,15 @@ List gapAwareEwma( continue; } final dt = timesMs[i] - lastT!; + if (dt <= 0 || halfLifeMs <= 0 || !dt.isFinite) { + // Duplicate or non-monotonic timestamp: NO time has elapsed, so the + // estimate earns no new weight (λ = 0). Letting λ = 1 − 2^(−dt/H) go + // negative here EXTRAPOLATES the estimate outside the data (Roberts 1959 + // requires λ ∈ [0,1]). We also keep lastT at the latest time already + // seen so an out-of-order sample cannot corrupt the next dt. + out.add(EwmaPoint(est, false)); + continue; + } // Time-aware λ: half-life expressed in ms => decay over the elapsed dt. var lambda = 1 - math.pow(2, -dt / halfLifeMs).toDouble(); final isGap = maxGapMs > 0 && dt > maxGapMs; diff --git a/lib/src/onehz/foundations/rr_correction.dart b/lib/src/onehz/foundations/rr_correction.dart index c20f653..c76f609 100644 --- a/lib/src/onehz/foundations/rr_correction.dart +++ b/lib/src/onehz/foundations/rr_correction.dart @@ -212,7 +212,14 @@ RrCorrectionResult correctRr( ); } -/// Time-varying threshold: alpha × (QD = (Q3−Q1)/2) of |x| in a sliding window. +/// Time-varying threshold: alpha × QD of the SIGNED series x in a sliding +/// window, where QD = (Q3 − Q1)/2 (Lipponen & Tarvainen 2019, eq. for th1/th2: +/// "th = α · quartile deviation of dRR over the 91-beat window", α = 5.2). +/// +/// The quartile deviation MUST be taken on the SIGNED series. Taking it on |x| +/// folds the symmetric ±dRR distribution onto one side, collapsing QD by ~an +/// order of magnitude; the threshold then sinks to [floor] and the detector +/// degenerates into a fixed 100 ms cut-off that flags ordinary RSA as ectopy. List _slidingThreshold( List x, int win, double alpha, double floor) { final n = x.length; @@ -223,14 +230,15 @@ List _slidingThreshold( final hi = math.min(n - 1, i + half); final seg = []; for (var k = lo; k <= hi; k++) { - seg.add(x[k].abs()); + seg.add(x[k]); } final q1 = percentile(seg, 25) ?? 0; final q3 = percentile(seg, 75) ?? 0; final qd = (q3 - q1) / 2; - // Floor keeps a gross outlier detectable even on (near-)quantized clean - // data where the QD collapses to 0 — but the floor sits well above normal - // beat-to-beat HRV wobble so it never flags the healthy signal. + // Floor keeps a gross outlier detectable on (near-)quantized clean data + // where the QD genuinely collapses to 0 (constant RR). On any series with + // real beat-to-beat variability α·QD dominates the floor, so the floor + // never governs a physiological signal. out[i] = math.max(alpha * qd, floor); } return out; diff --git a/lib/src/onehz/human/circadian_lifestyle.dart b/lib/src/onehz/human/circadian_lifestyle.dart index c997bf5..2a8425a 100644 --- a/lib/src/onehz/human/circadian_lifestyle.dart +++ b/lib/src/onehz/human/circadian_lifestyle.dart @@ -2,28 +2,100 @@ // Catalog §A: Social jetlag [PUB Wittmann/Roenneberg 2006] (ship first) and // chronotype label via MSFsc [PUB] (HR-acrophase variant is HEUR). // -// Mid-sleep is the clock-time midpoint of a sleep episode (decimal hours, may -// exceed 24 for after-midnight midpoints — we keep it on a continuous axis so -// circular wrap doesn't corrupt the mean). Social jetlag = signed difference -// between free-day (weekend) and work-day (weekday) mid-sleep. +// Mid-sleep is the clock-time midpoint of a sleep episode, supplied by the +// caller as LOCAL TIME-OF-DAY in [0, 24). That is a CIRCULAR quantity: 23.9 h +// and 0.1 h are 12 minutes apart, not 23.8 h apart. Every statistic here is +// therefore computed with circular methods (Fisher, *Statistical Analysis of +// Circular Data*, 1993, §2.2–2.3): the summary point is the circular median +// (ordinary median on the sample unwrapped about its circular mean) and the +// difference between two midpoints is the SHORTEST signed arc on the 24 h +// circle, in (−12, +12]. Taking a plain median/subtraction on raw clock-hours +// — what this file used to do — reported a ~23:50 weekday vs ~01:10 weekend +// midsleep as 22.5 h of social jetlag instead of +1.4 h. +// +// Social jetlag = signed shortest-arc difference between free-day (weekend) and +// work-day (weekday) mid-sleep [PUB Wittmann, Dinich, Merrow & Roenneberg, +// *Chronobiol Int* 2006;23(1-2):497–509]. // // MSFsc (sleep-corrected mid-sleep on free days) corrects the free-day midpoint // for oversleep relative to the weekly average, the standard MCTQ chronotype -// proxy. We DO NOT print absolute MSFsc minutes (catalog honesty rule) — only a -// coarse type label + a percentile-of-you for stability. +// proxy [PUB Roenneberg et al., *Curr Biol* 2004; *Sleep Med Rev* 2007]. We DO +// NOT print absolute MSFsc minutes (catalog honesty rule) — only a coarse type +// label + a percentile-of-you for stability. The label bands and the stability +// percentile are evaluated on the MCTQ band axis: clock-hours unwrapped about +// 06:00, i.e. (−6, +18], so that a 23:30 mid-sleep reads as −0.5 h (early type) +// rather than 23.5 h (which the old [0,24) bands mislabelled "evening type"). // // HONESTY: report STATE (your weekend runs later) not a clinical chronotype // diagnosis; gate ≥minDays with ≥2 free days; "—" when insufficient. +import 'dart:math' as math; + import '../types.dart'; import '../util.dart'; import 'percentile_of_you.dart'; +// ── Circular helpers on the 24 h clock ────────────────────────────────────── + +/// Wrap a clock-hour into [0, 24). +double _wrap24(double h) { + final r = h % 24.0; + return r < 0 ? r + 24.0 : r; +} + +/// Shortest signed arc `a − b` on the 24 h circle, in (−12, +12]. +double _circDiffH(double a, double b) { + var d = (a - b) % 24.0; + if (d < 0) d += 24.0; + if (d > 12.0) d -= 24.0; + return d; +} + +/// Unwrap [h] onto the continuous axis centred on [centre] — i.e. the +/// representative of [h] nearest [centre]. Result lies in (centre−12, centre+12]. +double _unwrapAround(double centre, double h) => centre + _circDiffH(h, centre); + +/// Circular mean direction of clock-hours (Fisher 1993 §2.2): the angle of the +/// resultant of the unit vectors. Null when the resultant vanishes (an +/// antipodal/uniform sample has no defined mean direction) or the sample is +/// empty. +double? _circMeanH(List hs) { + if (hs.isEmpty) return null; + var sx = 0.0, sy = 0.0; + for (final h in hs) { + final a = h * math.pi / 12.0; // 24 h ↔ 2π + sx += math.cos(a); + sy += math.sin(a); + } + if (math.sqrt(sx * sx + sy * sy) < 1e-9) return null; // no resultant + return _wrap24(math.atan2(sy, sx) * 12.0 / math.pi); +} + +/// Circular median of clock-hours (Fisher 1993 §2.3): the ordinary median taken +/// on the sample unwrapped about its circular mean, wrapped back to [0, 24). +/// Robust to the odd late night AND correct across midnight. Null when the +/// sample is empty or has no defined mean direction. +double? _circMedianH(List hs) { + if (hs.isEmpty) return null; + final anchor = _circMeanH(hs); + if (anchor == null) return null; + final unwrapped = [for (final h in hs) _unwrapAround(anchor, h)]; + return _wrap24(median(unwrapped)!); +} + +/// Anchor for the MCTQ chronotype band axis: 06:00. Mid-sleep on free days sits +/// in the small hours for essentially everyone, so unwrapping about 06:00 puts +/// early types at negative hours and evening types at 5–8 h — a monotone axis +/// the label bands can be read off directly. +const double _mctqBandAnchorH = 6.0; + class SocialJetlag { - final double sjlHours; // signed: free-day midsleep − work-day midsleep - final double absHours; // |sjlHours|, the headline "jet zones" magnitude - final double midSleepFree; // mean free-day mid-sleep (decimal h) - final double midSleepWork; // mean work-day mid-sleep (decimal h) + /// Signed SHORTEST ARC free-day − work-day mid-sleep on the 24 h circle, + /// in (−12, +12]. Positive => weekends run later. + final double sjlHours; + final double absHours; // |sjlHours| ≤ 12, the headline "jet zones" magnitude + final double midSleepFree; // circular median free-day mid-sleep, [0,24) + final double midSleepWork; // circular median work-day mid-sleep, [0,24) final int nFree; final int nWork; const SocialJetlag(this.sjlHours, this.absHours, this.midSleepFree, @@ -41,9 +113,12 @@ class SocialJetlag { /// Social jetlag from per-night mid-sleep clock-hours, split into free-day /// (e.g. weekend / unconstrained) and work-day midpoints. /// -/// [freeMidSleepH] / [workMidSleepH] are decimal clock-hours of the sleep -/// midpoint for each night in the respective category. We use the MEDIAN -/// (robust to the odd late night). Positive SJL => weekends run LATER. +/// [freeMidSleepH] / [workMidSleepH] are decimal LOCAL clock-hours in [0, 24) +/// of the sleep midpoint for each night in the respective category. We use the +/// CIRCULAR MEDIAN (robust to the odd late night, and correct across midnight) +/// and the SHORTEST SIGNED ARC between the two midpoints, so |SJL| can never +/// exceed 12 h [PUB Wittmann/Roenneberg 2006]. Positive SJL => weekends run +/// LATER. Metric socialJetlag( List freeMidSleepH, List workMidSleepH, { @@ -57,9 +132,17 @@ Metric socialJetlag( note: 'need ≥2 free-day and ≥2 work-day nights to compare', ); } - final msf = median(freeMidSleepH)!; - final msw = median(workMidSleepH)!; - final sjl = msf - msw; // signed + final msf = _circMedianH(freeMidSleepH); + final msw = _circMedianH(workMidSleepH); + if (msf == null || msw == null) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'mid-sleep times have no resultant direction on the 24 h circle ' + '— no meaningful midpoint (never imputed)', + ); + } + final sjl = _circDiffH(msf, msw); // signed shortest arc, (−12, +12] final conf = clamp( (freeMidSleepH.length + workMidSleepH.length) / 14.0, 0.3, 0.9); return Metric( @@ -68,13 +151,16 @@ Metric socialJetlag( confidence: conf, tier: Tier.high, inputs_used: inputs, - note: 'signed weekend−weekday mid-sleep drift (≈ jet zones), behavioral', + note: 'signed weekend−weekday mid-sleep drift on the 24 h circle ' + '(≈ jet zones, |SJL| ≤ 12 h), behavioral', ); } class Chronotype { - /// Sleep-corrected free-day mid-sleep, decimal clock-hours. INTERNAL — never - /// surfaced as an absolute number; drives the label + percentile only. + /// Sleep-corrected free-day mid-sleep, decimal clock-hours wrapped to [0,24). + /// INTERNAL — never surfaced as an absolute number; the label and the + /// stability percentile are derived from it on the MCTQ band axis (unwrapped + /// about 06:00), so both agree with this value. final double msfScHours; final String typeLabel; // coarse early/intermediate/late label final Metric? stability; // percentile-of-you for steadiness @@ -114,23 +200,36 @@ Metric chronotype( note: 'chronotype needs ≥14 days with ≥2 free days', ); } - final msf = median(freeMidSleepH)!; - final sdFree = median(freeSleepDurH)!; + // CIRCULAR median (see the helpers at the top): a free-day mid-sleep set + // straddling midnight must not be averaged on a linear axis. + final msfCirc = _circMedianH(freeMidSleepH); + if (msfCirc == null) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'free-day mid-sleep has no resultant direction on the 24 h circle', + ); + } + final sdFree = median(freeSleepDurH)!; // a DURATION — linear, not circular // Oversleep correction only when free-day sleep exceeds the weekly average. + // MSFsc = MSF − (SD_free − SD_week)/2 [PUB Roenneberg, MCTQ]. final correction = sdFree > avgWeekSleepDurH ? (sdFree - avgWeekSleepDurH) / 2.0 : 0.0; - final msfSc = msf - correction; + final msfSc = _wrap24(msfCirc - correction); - // Coarse label off MSFsc clock-hour (population anchors are only for the - // label band — never printed as a number). + // Coarse label off MSFsc, read on the MCTQ BAND AXIS — clock-hours unwrapped + // about 06:00, so a 23:30 mid-sleep is −0.5 h (early) rather than 23.5 h, + // which the old [0,24) bands mislabelled as the latest possible type. The + // population anchors are only for the label band — never printed as a number. + final band = _unwrapAround(_mctqBandAnchorH, msfSc); String label; - if (msfSc < 3.0) { + if (band < 3.0) { label = 'early type'; - } else if (msfSc < 4.0) { + } else if (band < 4.0) { label = 'moderate early type'; - } else if (msfSc < 5.0) { + } else if (band < 5.0) { label = 'intermediate type'; - } else if (msfSc < 6.0) { + } else if (band < 6.0) { label = 'moderate evening type'; } else { label = 'evening type'; @@ -138,7 +237,13 @@ Metric chronotype( Metric? stability; if (history.length >= 14) { - stability = percentileOfYou(msfSc, history, minN: 14); + // Compare on the same continuous band axis, so a history straddling + // midnight ranks by real chronotype distance instead of by clock wrap. + stability = percentileOfYou( + band, + [for (final h in history) _unwrapAround(_mctqBandAnchorH, _wrap24(h))], + minN: 14, + ); } return Metric( diff --git a/lib/src/onehz/human/coaching.dart b/lib/src/onehz/human/coaching.dart index 5c13f76..fee25d7 100644 --- a/lib/src/onehz/human/coaching.dart +++ b/lib/src/onehz/human/coaching.dart @@ -288,10 +288,22 @@ Metric vo2maxEstimate({ required Sex sex, required double? age, }) { - if (restingHr == null || maxHr == null || maxHr <= restingHr) { + // Uth-Sørensen-Overgaard-Pedersen 2004: VO2max ≈ 15.3 · (HRmax / HRrest). + // The ratio is only defined for a STRICTLY POSITIVE resting HR — a 0 (the + // package's off-skin sentinel, see types.dart HrSample) divides to Infinity, + // which Metric.toJson emits raw and jsonEncode then throws on. `maxHr <= + // restingHr` does not catch it, so guard the denominator explicitly and + // abstain. Non-finite inputs abstain for the same reason. + if (restingHr == null || + maxHr == null || + !restingHr.isFinite || + !maxHr.isFinite || + restingHr <= 0 || + maxHr <= restingHr) { return const Metric.absent( tier: Tier.estimate, inputs_used: ['resting_hr', 'max_hr'], + note: 'VO2max needs a positive resting HR below HRmax — "—" (never imputed)', ); } final vo2 = 15.3 * (maxHr / restingHr); @@ -324,6 +336,29 @@ Metric physiologicalAge({ required double? sleepEfficiency, required double? dailySteps, }) { + // ABSTAIN when NOTHING physiological was supplied. `score` starts at the + // chronological age and only the blocks below move it, so with every + // physiological input null this used to return a PRESENT metric reading + // "physioAge == your age, delta 0" — a fabricated result — while claiming six + // inputs it never saw. A physiological age with no physiology in it is not an + // estimate, it is the birth date restated. + final used = [ + if (vo2max != null) 'vo2max', + if (restingHr != null) 'resting_hr', + if (rmssd != null) 'rmssd', + if (sleepDurationH != null) 'sleep_duration', + if (sleepEfficiency != null) 'sleep_efficiency', + if (dailySteps != null) 'steps', + ]; + if (used.isEmpty) { + return const Metric.absent( + tier: Tier.estimate, + inputs_used: ['profile'], + note: 'no physiological input present — "—" (never imputed; ' + 'chronological age alone is not a physiological age)', + ); + } + var score = chronologicalAge; if (vo2max != null) { score -= ((vo2max - 35.0) / 5.0).clamp(-8.0, 8.0); @@ -349,20 +384,28 @@ Metric physiologicalAge({ score = score.clamp(18.0, 95.0); return Metric( value: PhysioAge(physioAge: score, deltaYears: score - chronologicalAge), - confidence: 0.35, + // Confidence tracks how much physiology actually went in: one input is a + // hint, all six is the intended estimate. + confidence: (0.15 + 0.035 * used.length).clamp(0.15, 0.35), tier: Tier.estimate, - inputs_used: const [ - 'profile', - 'vo2max', - 'resting_hr', - 'rmssd', - 'sleep', - 'steps', - ], - note: 'directional physiological-age estimate', + // inputs_used reports what was ACTUALLY used, never the full menu. + inputs_used: ['profile', ...used], + note: 'directional physiological-age estimate from ${used.length}/6 ' + 'physiological inputs', ); } +/// Unbiased (n−1) sample variance about a known mean. 0 for n < 2. +double _sampleVar(List xs, double m) { + if (xs.length < 2) return 0.0; + var s = 0.0; + for (final x in xs) { + final dx = x - m; + s += dx * dx; + } + return s / (xs.length - 1); +} + class JournalDay { final String date; final Set tags; @@ -378,6 +421,15 @@ class JournalEffect { final int nUntagged; final bool insufficient; final bool meaningful; + + /// Standardized effect size — Cohen's d = delta / pooled SD (Cohen 1988). + /// Null when the pooled within-group SD is 0 (both sides constant) or the + /// comparison was not run at all. Disclosed so the "meaningful" verdict is + /// auditable rather than a bare percentage. + final double? cohensD; + + /// Pooled within-group SD used for [cohensD]; null when not computed. + final double? pooledSd; const JournalEffect({ required this.outcome, required this.delta, @@ -387,6 +439,8 @@ class JournalEffect { required this.nUntagged, required this.insufficient, required this.meaningful, + this.cohensD, + this.pooledSd, }); } @@ -396,10 +450,26 @@ class JournalTagCorrelation { const JournalTagCorrelation(this.tag, this.effects); } +/// Per-tag effect of a journal entry on each outcome series. +/// +/// [outcomes] values must be POSITIONALLY ALIGNED to [dates] (same length); a +/// series of a different length cannot be attributed to dates at all, so it is +/// reported as insufficient rather than silently truncated or index-crashed. +/// +/// [minEffectPct] and [minCohensD] set the "meaningful" bar. A percentage +/// difference of means alone is NOT evidence: with 2 days per side, two +/// noisy series routinely differ by several percent. The verdict therefore also +/// requires a standardized effect size (Cohen's d = delta / pooled SD ≥ 0.5, +/// Cohen's conventional "medium" effect) so within-group spread is accounted +/// for. When both sides are exactly constant (pooled SD = 0) d is undefined and +/// we require [minNForZeroSpread] observations per side before calling it. List journalCorrelations({ required List journal, required List dates, required Map> outcomes, + double minEffectPct = 3.0, + double minCohensD = 0.5, + int minNForZeroSpread = 3, }) { final allTags = {for (final j in journal) ...j.tags}; final tagByDate = {for (final j in journal) j.date: j.tags}; @@ -407,6 +477,25 @@ List journalCorrelations({ for (final tag in allTags) { final effects = []; for (final entry in outcomes.entries) { + // LENGTH GUARD: `entry.value[i]` used to be indexed by dates.length with + // no check, so any outcome list shorter than `dates` threw RangeError. + // A misaligned series is not partially usable — we cannot know which + // dates the values belong to — so abstain for this outcome. + if (entry.value.length != dates.length) { + effects.add( + JournalEffect( + outcome: entry.key, + delta: 0, + pctChange: null, + higherSide: 'neither', + nTagged: 0, + nUntagged: 0, + insufficient: true, + meaningful: false, + ), + ); + continue; + } final tagged = []; final untagged = []; for (var i = 0; i < dates.length; i++) { @@ -437,6 +526,30 @@ List journalCorrelations({ final pct = untaggedMean.abs() < 1e-9 ? null : (delta / untaggedMean.abs()) * 100.0; + + // DISPERSION TEST. Pooled within-group SD (Cohen 1988): + // sp = sqrt( ((n1-1)·s1² + (n2-1)·s2²) / (n1+n2-2) ), d = delta / sp. + // Without it, "3% difference of two means" was reported as a meaningful + // journal effect off 2 days per side — a difference smaller than the + // day-to-day noise of either side. + final st = _sampleVar(tagged, taggedMean); + final su = _sampleVar(untagged, untaggedMean); + final dof = tagged.length + untagged.length - 2; + final pooledVar = dof > 0 + ? ((tagged.length - 1) * st + (untagged.length - 1) * su) / dof + : 0.0; + final pooledSd = pooledVar > 0 ? math.sqrt(pooledVar) : 0.0; + final d = pooledSd > 0 ? delta / pooledSd : null; + + final bigEnough = pct != null && pct.abs() >= minEffectPct; + final separated = d != null + ? d.abs() >= minCohensD + // Both sides exactly constant: d is undefined. Only trust it with a + // real number of observations behind each constant. + : (delta.abs() > 0 && + tagged.length >= minNForZeroSpread && + untagged.length >= minNForZeroSpread); + effects.add( JournalEffect( outcome: entry.key, @@ -446,7 +559,9 @@ List journalCorrelations({ nTagged: tagged.length, nUntagged: untagged.length, insufficient: false, - meaningful: pct != null && pct.abs() >= 3.0, + meaningful: bigEnough && separated, + cohensD: d, + pooledSd: pooledSd, ), ); } diff --git a/lib/src/onehz/respiration/resp_rate.dart b/lib/src/onehz/respiration/resp_rate.dart index 7a5fc10..e1f15cc 100644 --- a/lib/src/onehz/respiration/resp_rate.dart +++ b/lib/src/onehz/respiration/resp_rate.dart @@ -94,8 +94,8 @@ Metric rsaRespRate( // peak to be stable across them. A respiratory peak is sharp & resolution- // invariant; spurious HRV structure or artifact is not. final peaks = []; // br/min - double? bestPower; - double? bestPeakHz; + final peakHz = []; // the same peaks in Hz, index-aligned + final peakPwr = []; // their spectral power, index-aligned for (final grid in const [300, 450, 700]) { final ls = lombScargle(tSec, nnMs, freqGrid(rsaLoHz, rsaHiHz, grid)); if (ls == null) continue; @@ -103,14 +103,11 @@ Metric rsaRespRate( if (pk == null) continue; // Reject aliasing: a peak at/above Nyquist is not a real breathing rate. if (pk >= respHiHz) continue; - final pwr = _powerAt(ls, pk); peaks.add(pk * 60.0); - if (bestPower == null || pwr > bestPower) { - bestPower = pwr; - bestPeakHz = pk; - } + peakHz.add(pk); + peakPwr.add(_powerAt(ls, pk)); } - if (peaks.length < 2 || bestPeakHz == null) { + if (peaks.length < 2) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, @@ -127,7 +124,23 @@ Metric rsaRespRate( '(spread ${round6(spread)} br/min) — withheld', ); } - final brpm = median(peaks)!; + // ONE SOURCE for the reported triple. `brpm` used to be median(peaks) while + // `peakHz`/`power` came from the highest-POWER grid, so `peak_hz * 60` and + // `brpm` could disagree by up to [agreeBrpm] inside a single RespEstimate. + // We now pick the MEDOID grid — the grid whose peak is closest to the median + // across grids — and report its rate, frequency and power together, so + // `brpm == peakHz * 60` holds exactly and `power` is the power measured AT + // the reported frequency. (With three grids the medoid IS the median; with + // two, the nearer of the pair. Robustness still comes from the agreement gate + // above, which already rejected any disagreeing set.) + final medBrpm = median(peaks)!; + var best = 0; + for (var i = 1; i < peaks.length; i++) { + if ((peaks[i] - medBrpm).abs() < (peaks[best] - medBrpm).abs()) best = i; + } + final brpm = peaks[best]; + final bestPeakHz = peakHz[best]; + final bestPower = peakPwr[best]; // Confidence: high when clean & resolution-stable; penalize artifacts and // wide spread. Cap below 1 (PRV ceiling). final conf = clamp( @@ -140,8 +153,10 @@ Metric rsaRespRate( confidence: conf, tier: Tier.high, inputs_used: inputs, - note: 'RSA HF-peak respiratory rate (Lomb-Scargle on native beat times); ' - 'PRV-derived; 1 Hz Nyquist caps rate at 30 br/min', + note: 'RSA HF-peak respiratory rate (Lomb-Scargle on native beat times, ' + 'medoid of ${peaks.length} spectral resolutions — brpm, peak_hz and ' + 'power all come from that one grid); PRV-derived; 1 Hz Nyquist caps ' + 'rate at 30 br/min', ); } diff --git a/lib/src/onehz/sleep/advanced_stager.dart b/lib/src/onehz/sleep/advanced_stager.dart index ef78443..8f12676 100644 --- a/lib/src/onehz/sleep/advanced_stager.dart +++ b/lib/src/onehz/sleep/advanced_stager.dart @@ -1302,66 +1302,134 @@ class AdvancedSleepStager { return out; } + /// Longest accelerometer dropout (s) the dense-array build will carry the + /// last-known gravity vector across. + /// + /// A carry-forward is honest for a BRIEF dropout (a missed 1 Hz sample is far + /// more likely than a real posture change, and ENMO off a stale-but-still + /// vector reads "no movement" rather than fabricating motion). UNBOUNDED it + /// is fabrication in the opposite direction: hours with no accelerometer at + /// all become a perfectly still — i.e. perfectly asleep — stretch, so a + /// forced 8 h window holding 2 h of data scored TST 8 h / efficiency 100%. + /// + /// 60 s = two staging epochs, and the shortest wake bout the Webster/ + /// Cole-Kripke continuity rules treat as bridgeable (Webster et al. 1982; + /// Cole et al. 1992) — i.e. the finest resolution at which this pipeline + /// claims to separate sleep from wake. A carry-forward bounded by it can + /// therefore never manufacture a scorable sleep bout on its own. Seconds past + /// the bound are UNSTAGED: they are left out of staging entirely and reported + /// as wake. + static const int maxAccelCarryForwardSec = 60; + /// DEFAULT staging path — delegates to `cardioStager` (cardio_stager.dart). /// See the file header for why this is the default. `cardioStager` expects /// per-SECOND-indexed [AccelSample]/HR arrays (index i == second i from /// [start]), not the sparse timestamped [GravTs]/[HrTs] lists this file /// otherwise uses — so this builds that dense array explicitly: accel gaps - /// carry the last-known vector forward (a brief gap is far more likely a - /// missed sample than genuine movement — and ENMO from a stale-but-still - /// vector reads as "no movement", never fabricating motion that didn't - /// happen); HR gaps fill with 0, which `cardioStager` already documents as - /// its own "off-skin" contract — no fabrication either way. + /// carry the last-known vector forward for at most + /// [maxAccelCarryForwardSec]; HR gaps fill with 0, which `cardioStager` + /// already documents as its own "off-skin" contract — no fabrication either + /// way. + /// + /// HONESTY: seconds with no usable accelerometer are NOT staged. The window + /// is split at those gaps and each contiguous usable RUN is staged on its own + /// (so a dropout cannot pollute the neighbouring run's night baselines + /// either); the gap seconds, runs too short to stage, and a window where + /// `cardioStager` itself abstains all come back as WAKE — the "stay unstaged" + /// contract [stageWindow] documents. They must NEVER come back as 'light', + /// which is what a zero-data window used to report for its entire length. static List _stageSessionCardio(int start, int end, List grav, List hr, List rr) { final span = end - start; - if (span < 3 * epochS.round()) return [StageSegment(start, end, 'light')]; + if (span <= 0) return const []; + final epSec = epochS.round(); + final minStageableSec = 3 * epSec; + if (span < minStageableSec) return [StageSegment(start, end, 'wake')]; + final gByTs = {for (final g in grav) if (g.ts >= start && g.ts < end) g.ts: g}; final hByTs = {for (final h in hr) if (h.ts >= start && h.ts < end) h.ts: h}; final accel = List.filled( span, AccelSample(start * 1000.0, 0, 0, 1.0)); final hr1hz = List.filled(span, 0.0); - var haveGrav = false; + // usable[i] — second i has a real accel sample or a BOUNDED carry-forward. + final usable = List.filled(span, false); + var lastRealIdx = -1; for (var i = 0; i < span; i++) { final ts = start + i; final g = gByTs[ts]; if (g != null) { accel[i] = AccelSample(ts * 1000.0, g.x, g.y, g.z); - haveGrav = true; - } else if (haveGrav) { - accel[i] = accel[i - 1]; // carry-forward — see doc comment above. + lastRealIdx = i; + usable[i] = true; + } else if (lastRealIdx >= 0 && (i - lastRealIdx) <= maxAccelCarryForwardSec) { + // Bounded carry-forward — see [maxAccelCarryForwardSec]. The TIMESTAMP + // is this second's, not the stale sample's: cardioStager centres its RR + // windows on `accel[mid].tsMs`, so copying the old sample wholesale + // mis-centred every RMSSD/LF-HF window inside a gap. + final p = accel[i - 1]; + accel[i] = AccelSample(ts * 1000.0, p.x, p.y, p.z); + usable[i] = true; } hr1hz[i] = hByTs[ts]?.bpm ?? 0.0; // 0 = off-skin, cardioStager's own contract. } - if (!haveGrav) return [StageSegment(start, end, 'light')]; + final rSeg = [for (final r in rr) if (r.ts >= start && r.ts < end) r]; final rrMs = [for (final r in rSeg) r.rrMs]; final rrTsMs = [for (final r in rSeg) r.ts * 1000.0]; - final result = cardioStager(hr1hz, accel, rrMs: rrMs, rrTsMs: rrTsMs); - final nEpoch = result.base.stages.length; - if (nEpoch == 0) return [StageSegment(start, end, 'light')]; - final labels = List.generate(nEpoch, (i) { - switch (result.base.stages[i]) { - case SleepStage.wake: - return 'wake'; - case SleepStage.rem: - return 'rem'; - case SleepStage.nrem: - return (i < result.deepFlag.length && result.deepFlag[i]) - ? 'deep' - : 'light'; + // Everything not staged below stays 'wake' — the honest default. + final perSec = List.filled(span, 'wake'); + var i = 0; + while (i < span) { + if (!usable[i]) { + i++; + continue; } - }); - // Reuse the same edges/segment-building [_stageSession] uses, so callers - // get byte-identical StageSegment semantics regardless of which staging - // method produced them. - final edges = [ - for (var i = 0; i <= nEpoch; i++) start + i * epochS, - ]; - edges[nEpoch] = math.max(edges[nEpoch], end.toDouble()); - final grid = _EpochGrid(edges, nEpoch, [], [], [], [], [], []); - return _buildSegments(labels, grid, end); + var j = i; + while (j < span && usable[j]) { + j++; + } + if ((j - i) >= minStageableSec) { + final result = cardioStager( + hr1hz.sublist(i, j), + accel.sublist(i, j), + rrMs: rrMs, + rrTsMs: rrTsMs, + ); + final nEpoch = result.base.stages.length; + for (var e = 0; e < nEpoch; e++) { + final label = switch (result.base.stages[e]) { + SleepStage.wake => 'wake', + SleepStage.rem => 'rem', + SleepStage.nrem => (e < result.deepFlag.length && result.deepFlag[e]) + ? 'deep' + : 'light', + }; + final lo = i + e * epSec; + // The last epoch absorbs the run's sub-epoch remainder, exactly as + // [_buildSegments] extends the final segment to the window end. + final hi = e == nEpoch - 1 ? j : math.min(j, lo + epSec); + for (var k = lo; k < hi; k++) { + perSec[k] = label; + } + } + } + i = j; + } + + // Coalesce the per-second labels into contiguous [StageSegment]s tiling + // [start, end) — same segment semantics [_buildSegments] produces. + final segments = []; + var k = 0; + while (k < span) { + var m = k; + while (m < span && perSec[m] == perSec[k]) { + m++; + } + segments.add(StageSegment(start + k, start + m, perSec[k])); + k = m; + } + return segments; } static List _stageSession(int start, int end, List grav, diff --git a/lib/src/onehz/sleep/cardio_stager.dart b/lib/src/onehz/sleep/cardio_stager.dart index cbdd1dd..27c8cf1 100644 --- a/lib/src/onehz/sleep/cardio_stager.dart +++ b/lib/src/onehz/sleep/cardio_stager.dart @@ -263,18 +263,7 @@ CardioStagerResult cardioStager( final profile = userProfile ?? cardioUserProfile; final n = math.min(hr1hz.length, accel.length); final nEpoch = n ~/ epochSec; - if (nEpoch < 3) { - return CardioStagerResult( - const StagerResult( - stages: [], - epochSec: _epochSec, - wakePct: 0, - nremPct: 0, - remPct: 0), - const [], - 0, - ); - } + if (nEpoch < 3) return _abstain(epochSec); // ── per-second ENMO (motion) against a LOCALLY-ADAPTIVE 1 g reference ────── // A single whole-night gravity-magnitude reference (the old approach) is @@ -374,8 +363,18 @@ CardioStagerResult cardioStager( for (var e = 0; e < nEpoch; e++) if (still(e) && !hr[e].isNaN) hr[e] ]; - final hrMedGlobal = - median(sleepHr) ?? (mean([for (final h in hr) if (!h.isNaN) h]) ?? 60); + // HONESTY GATE — no HR anywhere ⇒ ABSTAIN, never substitute a nominal resting + // HR. Every decision below is HR-RELATIVE (the wake/arousal gate, the REM + // p25 HR floor, the deep cardiac-trough cut), so with `hr[e]` NaN for every + // epoch none of them can ever fire and the classifier falls through to NREM + // for the whole window — i.e. a band sitting on the nightstand would be + // reported as a perfect night of sleep. The old `?? 60` fallback made that + // fabrication look like a real baseline. A window with SOME HR still gets the + // whole-window mean as the baseline when no epoch qualifies as `still` (that + // is a real measurement, just not a quiet one). + 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 ───────────── @@ -580,6 +579,21 @@ CardioStagerResult cardioStager( ); } +/// Honest "cannot stage this window" result: NO epochs, no deep flags, zero +/// confidence. Callers must treat an empty [StagerResult.stages] as UNSTAGED +/// (see `advanced_stager._stageSessionCardio`), never as sleep. +CardioStagerResult _abstain(int epochSec) => CardioStagerResult( + StagerResult( + stages: const [], + epochSec: epochSec, + wakePct: 0, + nremPct: 0, + remPct: 0, + ), + const [], + 0, + ); + /// RMSSD (ms) of cleaned RR beats whose absolute time falls within a ±2.5-min /// window centred on epoch [s,t). Returns NaN when too few clean beats. double _windowRmssd(List rrMs, List rrTsMs, @@ -620,13 +634,28 @@ double _windowRmssd(List rrMs, List rrTsMs, /// 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. +/// +/// The flanking-sleep CONTEXT is measured against an immutable SNAPSHOT of the +/// hypnogram taken before the pass — Webster's rule scores each wake bout +/// against the ORIGINAL surrounding sleep (Webster et al. 1982; Cole et al. +/// 1992), not against sleep this same pass just manufactured. Reading the list +/// being mutated let every bridged bout count as context for the next one, so +/// bridging CASCADED: a fragmented night of short sleep bouts separated by long +/// wake bouts collapsed into one continuous sleep block (WASO 0, efficiency +/// 100%). Exposed (non-private) so the regression test can drive the rule +/// directly; not part of the package's public barrel. +void websterRescoreCardio(List sm, int epochSec) => + _websterRescore(sm, epochSec); + void _websterRescore(List sm, int epochSec) { bool isSleep(SleepStage s) => s != SleepStage.wake; final n = sm.length; double minToEp(double m) => m * 60.0 / epochSec; + // Immutable context snapshot — see the doc comment above. + final snap = List.of(sm); var onset = -1, lastSleep = -1; for (var i = 0; i < n; i++) { - if (isSleep(sm[i])) { + if (isSleep(snap[i])) { if (onset < 0) onset = i; lastSleep = i; } @@ -634,7 +663,7 @@ void _websterRescore(List sm, int epochSec) { if (onset < 0) return; int runBefore(int i) { var c = 0, k = i - 1; - while (k >= onset && isSleep(sm[k])) { + while (k >= onset && isSleep(snap[k])) { c++; k--; } @@ -642,7 +671,7 @@ void _websterRescore(List sm, int epochSec) { } int runAfter(int i) { var c = 0, k = i + 1; - while (k <= lastSleep && isSleep(sm[k])) { + while (k <= lastSleep && isSleep(snap[k])) { c++; k++; } @@ -659,12 +688,12 @@ void _websterRescore(List sm, int epochSec) { ]; var i = onset; while (i <= lastSleep) { - if (isSleep(sm[i])) { + if (isSleep(snap[i])) { i++; continue; } var j = i; - while (j <= lastSleep && !isSleep(sm[j])) { + while (j <= lastSleep && !isSleep(snap[j])) { j++; } final wakeLen = (j - i).toDouble(); diff --git a/lib/src/onehz/sleep/cpc.dart b/lib/src/onehz/sleep/cpc.dart index 2890298..c009060 100644 --- a/lib/src/onehz/sleep/cpc.dart +++ b/lib/src/onehz/sleep/cpc.dart @@ -111,7 +111,43 @@ Metric cardiopulmonaryCoupling( final vlfc = couplingLs.bandPower(0.001, 0.01); final lfc = couplingLs.bandPower(0.01, 0.1); final hfc = couplingLs.bandPower(0.1, 0.4); - final ratio = lfc > 0 ? hfc / lfc : (hfc > 0 ? double.infinity : 0.0); + + // A non-finite band power means the input itself was poisoned (a NaN/±inf RR + // or beat time survives every arithmetic step, and `variance <= 0` does not + // catch NaN). Nothing measured here is real, so nothing is reported. + if (!hfc.isFinite || !lfc.isFinite || !vlfc.isFinite) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'coupling spectrum non-finite (NaN/inf in the NN series or beat ' + 'times) — no coupling was measured', + ); + } + + // Thomas et al. 2005 defines the stability index as the RATIO of high- to + // low-frequency coupling power. With no LFC power at all the ratio is a + // DIVISION BY ZERO — mathematically undefined, not "infinitely stable" (and + // 0/0 is not "maximally unstable" either). The published method says nothing + // about a spectrum with an empty 0.01–0.1 Hz band, so we abstain rather than + // report a number. (Pre-2026-07 this emitted the sentinel 999.0 inside a live + // Metric at confidence up to 0.85, which downstream read as a real, + // extraordinarily stable night.) + if (lfc <= 0) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'no low-frequency (0.01-0.1 Hz) coupling power — the Thomas 2005 ' + 'HFC/LFC stability ratio is undefined here, not "perfectly stable"', + ); + } + final ratio = hfc / lfc; + if (!ratio.isFinite) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'HFC/LFC stability ratio overflowed — not a measurement', + ); + } // Confidence grows with record length; capped (a screen, not a diagnosis). final conf = clamp(n / 3600.0, 0.3, 0.85); @@ -120,7 +156,7 @@ Metric cardiopulmonaryCoupling( hfc: hfc, lfc: lfc, vlfc: vlfc, - cpcRatio: ratio.isFinite ? ratio : 999.0, + cpcRatio: ratio, dominantHz: domF, ), confidence: conf, diff --git a/lib/src/onehz/sleep/cycles.dart b/lib/src/onehz/sleep/cycles.dart index 81d61ed..6ef0a72 100644 --- a/lib/src/onehz/sleep/cycles.dart +++ b/lib/src/onehz/sleep/cycles.dart @@ -83,7 +83,12 @@ SleepCyclesResult detectSleepCycles( final perMin = List.filled(nMin, null); final bins = List>.generate(nMin, (_) => []); for (var k = 0; k < rrMs.length; k++) { - final m = (rrTsMs[k] ~/ 1000 - onsetSec) ~/ 60; + // FLOOR-divide, don't truncate. Dart's `~/` rounds toward ZERO, so a beat + // 1–59 s BEFORE onset produced bin 0 instead of a negative bin and slipped + // past the `m < 0` guard below — up to 59 s of pre-onset (often awake, + // high-variability) beats polluted minute 0 of the cycle series. + final relSec = (rrTsMs[k] / 1000.0).floor() - onsetSec; + final m = (relSec / 60.0).floor(); if (m < 0 || m >= nMin) continue; final v = rrMs[k]; if (v >= _rrMin && v <= _rrMax) bins[m].add(v); diff --git a/lib/src/onehz/sleep/segment.dart b/lib/src/onehz/sleep/segment.dart index b64f5ff..6c7ff3b 100644 --- a/lib/src/onehz/sleep/segment.dart +++ b/lib/src/onehz/sleep/segment.dart @@ -318,6 +318,11 @@ SleepSegmentation segmentSleep( offsetMs: chosen.end * 1000.0, immobile: fallbackWindow?.immobile ?? List.filled(trimmedAccel.length, false), + // Forward the undecidable-second mask too. Dropping it here silently + // downgraded "we could not tell" into "not immobile", which is the + // conservative direction but costs the caller `unresolvedTailSec` — the + // one signal that says a night ran past the end of the record. + immobileUnknown: fallbackWindow?.immobileUnknown ?? const [], zAngleDeg: fallbackWindow?.zAngleDeg ?? List.filled(trimmedAccel.length, 0.0), sptSec: inBed, @@ -342,8 +347,19 @@ class _SleepGroup { final int start; final int end; final double asleepMin; + + /// SUM of the bridged sessions' durations — the bridge GAPS are excluded, so + /// this is NOT `end - start` for a multi-session group. Do not use it for + /// time-of-day math: `start + inBedSec ~/ 2` lands half the total gap EARLY + /// (a single 50-min bridge ⇒ 25 min early). Use [midsleepSec] for that. final int inBedSec; + /// The group's circadian centre: the midpoint of its ACTUAL SPAN. This is + /// what a midsleep anchor is defined against (the middle of the sleep period, + /// gaps included — Roenneberg's MSF/mid-sleep convention), and it is what + /// [_pickMainSleepGroup] must compare with a habitual-midsleep anchor. + int get midsleepSec => start + (end - start) ~/ 2; + const _SleepGroup({ required this.sessions, required this.start, @@ -432,8 +448,10 @@ _SleepGroup? _pickMainSleepGroup( } double alignmentBonusFor(_SleepGroup g) { - final mid = g.start + (g.inBedSec ~/ 2); - final dist = circularDistanceSec(localSecOfDay(mid), targetMidsleepSec); + // Midsleep = the middle of the SPAN, never `start + inBedSec/2` — see + // [_SleepGroup.inBedSec]/[_SleepGroup.midsleepSec]. + final dist = + circularDistanceSec(localSecOfDay(g.midsleepSec), targetMidsleepSec); if (dist <= fullWindowSec) return alignmentBonusMin; if (dist >= zeroWindowSec) return 0.0; final frac = (zeroWindowSec - dist) / (zeroWindowSec - fullWindowSec); @@ -454,12 +472,36 @@ _SleepGroup? _pickMainSleepGroup( return winner; } +/// Habitual midsleep anchor (local second-of-day) from ≥[minDays] of history. +/// +/// Timezone conversion is PER TIMESTAMP, with exactly the same precedence +/// [segmentSleep] uses for its own local-time-of-day math — because the anchor +/// this returns is compared against `_pickMainSleepGroup`'s per-timestamp +/// conversion, and the two must agree: +/// * [tzOffsetResolver] (if given) wins — inject a deterministic ts→offset map. +/// * else a fixed [tzOffsetSeconds] (if given) — LEGACY/deterministic only. +/// * else the machine's offset in effect AT EACH timestamp (DST-correct). +/// +/// A single frozen offset applied to every history block is a DST bypass: the +/// ≥14 days this function requires will regularly straddle a transition, so +/// roughly half the days convert with the wrong offset and the circular-mean +/// anchor is biased by up to ~30 min against the DST-correct comparison it +/// feeds. Prefer passing NOTHING (machine-local, DST-correct) or a resolver; +/// pass [tzOffsetSeconds] only when a caller genuinely wants one frozen offset. int? habitualMidsleepSecFromHistory( List<({int startSec, int endSec, String dayKey})> history, { - required int tzOffsetSeconds, + int? tzOffsetSeconds, + int Function(int tsSec)? tzOffsetResolver, int minDays = 14, }) { if (history.isEmpty) return null; + final int Function(int tsSec) tzAt = tzOffsetResolver ?? + (tzOffsetSeconds != null + ? (int _) => tzOffsetSeconds + : (int tsSec) => DateTime.fromMillisecondsSinceEpoch( + tsSec * 1000, + isUtc: false, + ).timeZoneOffset.inSeconds); final longestByDay = {}; for (final block in history) { final cur = longestByDay[block.dayKey]; @@ -474,10 +516,12 @@ int? habitualMidsleepSecFromHistory( if (longestByDay.length < minDays) return null; final mids = [ for (final block in longestByDay.values) - _localSecOfDay( - block.startSec + ((block.endSec - block.startSec) ~/ 2), - tzOffsetSeconds, - ), + // Convert with the offset in effect AT THAT BLOCK'S midsleep instant, not + // one offset frozen for the whole history — see the doc comment above. + () { + final mid = block.startSec + ((block.endSec - block.startSec) ~/ 2); + return _localSecOfDay(mid, tzAt(mid)); + }(), ]; return _circularMeanSec(mids); } diff --git a/lib/src/onehz/sleep/stager.dart b/lib/src/onehz/sleep/stager.dart index 46a0cbb..e75d49b 100644 --- a/lib/src/onehz/sleep/stager.dart +++ b/lib/src/onehz/sleep/stager.dart @@ -215,6 +215,11 @@ Metric autonomicStager( ); } +/// Test seam for [_websterRescore] — exposed (non-private) so the regression +/// test can drive the continuity rule directly. Not part of the public barrel. +void websterRescoreAutonomic(List sm, int epochSec) => + _websterRescore(sm, epochSec); + /// Webster/Cole-Kripke sleep-continuity rescoring (in place). /// /// After sleep onset (first sleep epoch), a contiguous run of WAKE epochs that @@ -227,16 +232,26 @@ Metric autonomicStager( /// sleep-before ≥15 min AND wake-run ≤ 4 min → sleep /// "sleep-before" is satisfied by sustained sleep on EITHER bracketing side, /// so an arousal sandwiched between two long sleep bouts is bridged. +/// +/// The flanking-sleep CONTEXT is measured against an immutable SNAPSHOT of the +/// hypnogram taken before the pass — Webster's rule scores each wake bout +/// against the ORIGINAL surrounding sleep (Webster et al. 1982; Cole et al. +/// 1992), not against sleep this same pass just manufactured. Reading the list +/// being mutated let every bridged bout count as context for the next one, so +/// bridging CASCADED and a genuinely fragmented night collapsed into one +/// continuous sleep block (WASO 0, efficiency 100%). void _websterRescore(List sm, int epochSec) { bool isSleep(SleepStage s) => s != SleepStage.wake; final n = sm.length; // epochs-per-minute (rounded; epochSec=30 → 2/min). double minToEp(double m) => m * 60.0 / epochSec; + // Immutable context snapshot — see the doc comment above. + final snap = List.of(sm); // Locate sleep onset & final sleep epoch — only rescore within the sleep body. var onset = -1, lastSleep = -1; for (var i = 0; i < n; i++) { - if (isSleep(sm[i])) { + if (isSleep(snap[i])) { if (onset < 0) onset = i; lastSleep = i; } @@ -247,7 +262,7 @@ void _websterRescore(List sm, int epochSec) { int sleepRunBefore(int i) { var c = 0; var k = i - 1; - while (k >= onset && isSleep(sm[k])) { + while (k >= onset && isSleep(snap[k])) { c++; k--; } @@ -258,7 +273,7 @@ void _websterRescore(List sm, int epochSec) { int sleepRunAfter(int i) { var c = 0; var k = i + 1; - while (k <= lastSleep && isSleep(sm[k])) { + while (k <= lastSleep && isSleep(snap[k])) { c++; k++; } @@ -279,13 +294,13 @@ void _websterRescore(List sm, int epochSec) { var i = onset; while (i <= lastSleep) { - if (isSleep(sm[i])) { + if (isSleep(snap[i])) { i++; continue; } // WAKE run [i, j). var j = i; - while (j <= lastSleep && !isSleep(sm[j])) { + while (j <= lastSleep && !isSleep(snap[j])) { j++; } final wakeLen = j - i; diff --git a/lib/src/onehz/sleep/van_hees.dart b/lib/src/onehz/sleep/van_hees.dart index 3e33909..36ba7a6 100644 --- a/lib/src/onehz/sleep/van_hees.dart +++ b/lib/src/onehz/sleep/van_hees.dart @@ -37,9 +37,25 @@ class SleepWindow { final double? onsetMs; final double? offsetMs; - /// Per-second immobility mask (true = "no movement" second), full length. + /// Per-second immobility mask (true = ASSERTED "no movement" second), full + /// length. `false` means "not asserted immobile" — which is either a real + /// movement OR an UNDECIDABLE second (see [immobileUnknown]). It is never a + /// guess: the van Hees rule needs `sustainedMin` of FUTURE angle data, and the + /// last `sustainedMin` of any record does not have it. final List immobile; + /// Per-second "undecidable" mask, full length and index-aligned with + /// [immobile]. `true` marks a second whose forward sustained-inactivity + /// window is truncated by the end of the record AND which shows no movement + /// in the data we do have — i.e. it could still resolve either way once the + /// next samples arrive. Such seconds are `false` in [immobile] (not claimed + /// as rest) and are excluded from the detected sleep period. + /// + /// Consumers that only read [immobile] therefore degrade conservatively (an + /// unresolved tail second reads as "not asserted immobile"), never + /// optimistically. Empty when unknown. + final List immobileUnknown; + /// Per-second smoothed z-angle (deg), full length — reused downstream. final List zAngleDeg; @@ -54,14 +70,26 @@ class SleepWindow { required this.immobile, required this.zAngleDeg, required this.sptSec, + this.immobileUnknown = const [], }); + /// Seconds at the end of the record whose immobility is genuinely + /// undecidable (see [immobileUnknown]). + int get unresolvedTailSec { + var c = 0; + for (final u in immobileUnknown) { + if (u) c++; + } + return c; + } + Map toJson() => { 'onset_idx': onsetIdx, 'offset_idx': offsetIdx, if (onsetMs != null) 'onset_ms': onsetMs, if (offsetMs != null) 'offset_ms': offsetMs, 'spt_sec': sptSec, + if (unresolvedTailSec > 0) 'unresolved_tail_sec': unresolvedTailSec, }; } @@ -112,32 +140,42 @@ Metric vanHeesSleepWindow( } final win = sustainedMin * 60; final immobile = List.filled(n, false); - // A second is "no movement" if the MAX absolute angle change over the - // surrounding `win` seconds stays below the threshold. + final immobileUnknown = List.filled(n, false); + // A second is "no movement" if the MAX absolute angle change over the `win` + // seconds STARTING AT IT stays below the threshold (van Hees 2015/2018: the + // sustained-inactivity rule is a property of the block that follows the + // second, so it is evaluated PER SECOND on that second's own window). + // + // The final `win-1` seconds have a TRUNCATED forward window, which splits + // them into two honest cases — never into one shared verdict (pre-2026-07 the + // tail took a single global `[n-win, n)` answer and stamped it on every tail + // second, so one twitch anywhere in the last 5 min flipped the whole tail to + // mobile, and stillness observed BEFORE a tail second was used to certify a + // full sustained window that second never actually had): + // * movement seen in [i, n) → the ≥sustainedMin rule already FAILS on the + // data in hand, whatever comes next ⇒ decided, not immobile. + // * no movement seen, window short ⇒ UNDECIDABLE. We do not assert rest + // (that would extrapolate stillness past the end of the record) and we + // record it in `immobileUnknown` rather than guessing either way. for (var i = 0; i < n; i++) { - final lo = i; final hi = math.min(n, i + win); - if (hi - lo < win) { - // Tail shorter than a full window: fall back to the trailing window. - final lo2 = math.max(0, n - win); - var maxd = 0.0; - for (var k = lo2 + 1; k < n; k++) { - if (dAng[k] > maxd) maxd = dAng[k]; - } - immobile[i] = maxd < angleThresholdDeg && (n - lo2) >= win; - continue; - } var maxd = 0.0; - for (var k = lo + 1; k < hi; k++) { + for (var k = i + 1; k < hi; k++) { if (dAng[k] > maxd) { maxd = dAng[k]; if (maxd >= angleThresholdDeg) break; } } - immobile[i] = maxd < angleThresholdDeg; + final still = maxd < angleThresholdDeg; + final fullWindow = hi - i >= win; + immobile[i] = still && fullWindow; + immobileUnknown[i] = still && !fullWindow; } - // 4. longest immobile block, bridging brief gaps. + // 4. longest immobile block, bridging brief gaps. Only ASSERTED immobile + // seconds extend a block, so a night still running when the record ends is + // reported up to the last second we can actually certify — the undecidable + // tail is left out rather than annexed on the assumption it stayed still. final bridge = bridgeGapMin * 60; var bestStart = -1, bestEnd = -1, bestLen = 0; var i = 0; @@ -185,6 +223,11 @@ Metric vanHeesSleepWindow( final offsetMs = hasTs ? accel[math.min(bestEnd, n - 1)].tsMs : null; + var unresolved = 0; + for (final u in immobileUnknown) { + if (u) unresolved++; + } + // Confidence grows with the detected SPT length up to a typical night. final conf = clamp(bestLen / (7 * 3600), 0.3, 0.95); return Metric( @@ -194,6 +237,7 @@ Metric vanHeesSleepWindow( onsetMs: onsetMs, offsetMs: offsetMs, immobile: immobile, + immobileUnknown: immobileUnknown, zAngleDeg: ang, sptSec: bestLen, ), @@ -201,7 +245,9 @@ Metric vanHeesSleepWindow( tier: Tier.high, inputs_used: inputs, note: 'van Hees angle-based REST window (5°/${sustainedMin}min); ' - 'a rest period, not PSG sleep', + 'a rest period, not PSG sleep' + '${unresolved > 0 ? '; last ${unresolved}s of the record are ' + 'undecidable (forward window truncated) and are excluded' : ''}', ); } diff --git a/lib/src/onehz/wellness/anomaly.dart b/lib/src/onehz/wellness/anomaly.dart index f162c04..ee9c31b 100644 --- a/lib/src/onehz/wellness/anomaly.dart +++ b/lib/src/onehz/wellness/anomaly.dart @@ -14,7 +14,10 @@ // (HRV is negated: a DROP in HRV is the illness direction). // // HONESTY: this is a complement, not a diagnosis. Missing features reduce the -// vector dimension (we never impute). Persistence + a conservative chi-square +// vector dimension (we never impute), and so does a feature whose baseline has +// no dispersion at all (MAD = 0 AND SD = 0) — there is no scale to standardize +// against, so it is DROPPED rather than floored to an epsilon. Persistence + a +// conservative chi-square // gate keep the false-positive rate honest, and we report the per-feature // contributions so a flag is explainable. @@ -42,9 +45,13 @@ class AnomalyDay { final bool candidate; // distance crossed gate THIS night (pre-persistence) final List drivers; // per-feature signed contribution - /// Machine-readable "need_baseline:have=H,need=N" note set on nights that - /// could not be evaluated for lack of baseline coverage (H = best per-feature - /// baseline count available, N = required minimum). Null when evaluated. + /// Machine-readable note set on nights that could not be evaluated: + /// * "need_baseline:have=H,need=N" — insufficient baseline coverage (H = + /// best per-feature baseline count available, N = required minimum). + /// * "degenerate_baseline:no_dispersion" — the surviving baseline columns + /// are exactly constant (MAD = 0 AND SD = 0), so there is no scale to + /// standardize against and we abstain rather than invent one. + /// Null when the night WAS evaluated. final String? need; const AnomalyDay(this.date, this.mahalanobis, this.flagged, this.candidate, this.drivers, {this.need}); @@ -125,20 +132,41 @@ List multivariateAnomaly( run = 0; continue; } - // Robust center (median) + scale (MAD) per available feature. - final center = [for (final f in idx) median(cols[f])!]; - final scale = [ - for (final f in idx) - () { - final s = mad(cols[f]) ?? 0; - return s <= 0 ? (stddev(cols[f]) ?? 1.0).clamp(1e-6, 1e9) : s; - }() - ]; + // Robust center (median) + scale (MAD, ordinary SD as the coarser fallback) + // per available feature. + // + // ABSTAIN, NEVER FLOOR: a feature whose trailing baseline has NO dispersion + // at all (MAD == 0 AND SD == 0 — an exactly-constant, fully-quantized column + // such as a skin-temp z that reads 0.0 every night) has no scale to + // standardize against. Clamping the scale to an epsilon (the old `.clamp( + // 1e-6, 1e9)`) turned any deviation into a ~1e6 z, so d² blew past the χ² + // gate unconditionally and a 0.4-unit change surfaced as an illness anomaly. + // The sibling modules already refuse this case — readiness_composite's + // `robustZ(v, base) ?? z(v, base)` yields null when SD is also 0, and + // changepoint guards zero variance — so we match them: DROP the degenerate + // feature from the vector, and if fewer than 2 features survive, abstain. + final keep = []; + final center = []; + final scale = []; + for (final f in idx) { + final m = mad(cols[f]) ?? 0; + final sc = m > 0 ? m : (stddev(cols[f]) ?? 0); + if (!sc.isFinite || sc <= 0) continue; // no dispersion → not standardizable + keep.add(f); + center.add(median(cols[f])!); + scale.add(sc); + } + if (keep.length < 2) { + out.add(AnomalyDay(dates[i], null, false, false, const [], + need: 'degenerate_baseline:no_dispersion')); + run = 0; + continue; + } // Standardized current vector. - final zc = [for (var a = 0; a < idx.length; a++) (cur[idx[a]]! - center[a]) / scale[a]]; + final zc = [for (var a = 0; a < keep.length; a++) (cur[keep[a]]! - center[a]) / scale[a]]; // Robust correlation matrix from aligned rows (standardized), regularized. - final cov = _robustCorr(rows, idx, center, scale, ridge); + final cov = _robustCorr(rows, keep, center, scale, ridge); final inv = _invert(cov); double d2; if (inv == null) { @@ -157,13 +185,13 @@ List multivariateAnomaly( // Per-feature contribution to d² (diagonal share), for the "why". final drivers = []; - for (var a = 0; a < idx.length; a++) { - drivers.add(Driver(_featLabels[idx[a]], round6(zc[a]), + for (var a = 0; a < keep.length; a++) { + drivers.add(Driver(_featLabels[keep[a]], round6(zc[a]), detail: 'standardized deviation')); } drivers.sort((x, y) => y.contribution.abs().compareTo(x.contribution.abs())); - final gate = chiSqGate ?? _chiSq999(idx.length); + final gate = chiSqGate ?? _chiSq999(keep.length); final candidate = d2 > gate; if (candidate) { run++; diff --git a/lib/src/onehz/wellness/readiness_composite.dart b/lib/src/onehz/wellness/readiness_composite.dart index 0e9ad1f..6337ec9 100644 --- a/lib/src/onehz/wellness/readiness_composite.dart +++ b/lib/src/onehz/wellness/readiness_composite.dart @@ -115,15 +115,23 @@ Metric readinessComposite( // nights, absent others, even with sleep detected). Fall back to an ordinary // mean/SD z so a usable input still contributes; only skip when SD is ALSO // zero (a truly constant baseline with no dispersion to normalize against). - final zr = robustZ(v, base) ?? z(v, base); + final rz = robustZ(v, base); + final zr = rz ?? z(v, base); if (zr == null) continue; final oriented = inp.goodSign * zr; // + = good for readiness used.add(inp.label); weightSum += inp.weight; weightedZ += inp.weight * oriented; // Driver contribution is the signed weighted z (renormalized later). + // GLASS-BOX: the disclosed method must be the method ACTUALLY used — on a + // quantized baseline the MAD collapses and the mean/SD fallback above + // produced this z, so saying "robust-z" there would misstate how the + // contribution was computed. + final method = rz != null + ? 'robust-z (median+MAD)' + : 'z (mean+SD fallback — MAD=0 on a quantized baseline)'; drivers.add(Driver(inp.label, inp.weight * oriented, - detail: 'oriented robust-z=${round6(oriented)}')); + detail: 'oriented $method=${round6(oriented)}')); } if (used.isEmpty || weightSum == 0) { // If inputs HAD values but their baselines were too short, say so in the diff --git a/lib/src/onehz/workout/workout_detect.dart b/lib/src/onehz/workout/workout_detect.dart index 238f26b..79f4556 100644 --- a/lib/src/onehz/workout/workout_detect.dart +++ b/lib/src/onehz/workout/workout_detect.dart @@ -45,6 +45,15 @@ class ExerciseSession { final double? caloriesKcal; final double? caloriesKJ; + /// True when [caloriesKcal]/[caloriesKJ] were computed against a FABRICATED + /// anchor — [Calories.estimateBoutCalories] falls back to a flat + /// `hrmax = 220` / `restingHr = 60` when either is null, and returns + /// `usedDefaultAnchors` precisely so that number can be caveated instead of + /// shown as if it were personal. The flag used to be computed and thrown + /// away; it is now carried through to [toJson]. False when no calories were + /// computed at all. + final bool caloriesUsedDefaultAnchors; + /// Sport label from the classifier seam ("detected" by default). final String sport; @@ -61,6 +70,7 @@ class ExerciseSession { required this.hrmaxSource, required this.caloriesKcal, required this.caloriesKJ, + this.caloriesUsedDefaultAnchors = false, this.sport = defaultSportLabel, }); @@ -77,6 +87,7 @@ class ExerciseSession { 'hrmax_source': hrmaxSource, 'calories_kcal': caloriesKcal == null ? null : round6(caloriesKcal!), 'calories_kj': caloriesKJ == null ? null : round6(caloriesKJ!), + 'calories_used_default_anchors': caloriesUsedDefaultAnchors, 'sport': sport, }; } @@ -134,12 +145,22 @@ class WorkoutDetector { return out; } - /// Day resting-HR baseline = nearest-rank RESTING_PERCENTILE of bpm values. - /// Derive resting HR from a sorted series. [bpmSorted] must already be sorted ascending. - static double _deriveRestingHR(List bpmSorted) { - final rank = - math.max(1, (restingPercentile / 100.0 * bpmSorted.length).ceil()); - return bpmSorted[rank - 1]; + /// Day resting-HR baseline = nearest-rank RESTING_PERCENTILE of the ON-SKIN + /// bpm values. + /// + /// OFF-SKIN FILTER: the package convention (types.dart, `HrSample`) is that + /// `hr == 0` means the sensor was off the wrist, NEVER bradycardia. Taking the + /// 10th percentile of the RAW stream on a day with ≥10 % dropout returned + /// restHR = 0, which dragged hrFloor down to 15 bpm and inflated every + /// downstream %HRR — an ordinary 120 bpm walk read as 63 % HRR (zone 2) + /// instead of 48 % (zone 0), so it cleared the zone-2 workout gate. Zeros are + /// dropped before the percentile; null when nothing on-skin remains (we + /// abstain rather than invent a resting HR). + static double? _deriveRestingHR(List bpm) { + final onSkin = [for (final b in bpm) if (b > 0 && b.isFinite) b]..sort(); + if (onSkin.isEmpty) return null; + final rank = math.max(1, (restingPercentile / 100.0 * onSkin.length).ceil()); + return onSkin[rank - 1]; } /// Value whose ts is nearest [target] within [tol] s, else null. Ties → later @@ -286,7 +307,10 @@ class WorkoutDetector { final motion = activitySeries(gravTs, gx, gy, gz); if (motion.isEmpty) return const []; - final restHR = restingHR ?? _deriveRestingHR([...sBpm]..sort()); + final restHR = restingHR ?? _deriveRestingHR(sBpm); + // No caller RHR and no on-skin sample to derive one from → no baseline, so + // no honest HR gate. Abstain rather than gate against a fabricated floor. + if (restHR == null) return const []; final hrFloor = restHR + hrMarginBPM; final double? effMaxHR; @@ -355,14 +379,21 @@ class WorkoutDetector { } // Intensity qualification: require ≥ MIN_INTENSITY_Z2PLUS in zone 2+. - if (zonePct.isNotEmpty) { - var z2plus = 0.0; - for (var z = 2; z <= 5; z++) { - z2plus += zonePct[z] ?? 0.0; - } - z2plus /= 100.0; - if (z2plus < minIntensityZ2Plus) continue; + // + // AN UNEVALUABLE GATE BLOCKS, IT DOES NOT PASS. `zonePct` is empty exactly + // when there is no usable HRmax anchor (no caller HRmax, no age for + // Tanaka, <600 samples for an observed estimate → estimateHRmax returns + // ("unknown", 0)), or when the anchor is at/below resting HR. The old code + // skipped the whole gate in that case, so a 6-minute walk at RHR+16 bpm + // was emitted as a durable workout. With no zone breakdown we cannot know + // whether the bout qualified, so we drop it. + if (zonePct.isEmpty) continue; + var z2plus = 0.0; + for (var z = 2; z <= 5; z++) { + z2plus += zonePct[z] ?? 0.0; } + z2plus /= 100.0; + if (z2plus < minIntensityZ2Plus) continue; // OVERLAP-DEDUP: drop a detected bout overlapping a saved/manual span. if (savedSpans.any((s) => _overlaps(start, end, s.startSec, s.endSec))) { @@ -370,6 +401,11 @@ class WorkoutDetector { } double? kcal, kj; + // Carry [Calories.estimateBoutCalories]'s usedDefaultAnchors through to + // the session — it exists so a calorie number built on the flat + // `hrmax ?? 220` / `restingHr ?? 60` fallback can be caveated, and it used + // to be computed and dropped on the floor here. + var calUsedDefaultAnchors = false; if (profile != null) { final winBpmInt = [for (final b in winBpm) b]; final cal = Calories.estimateBoutCalories( @@ -382,17 +418,29 @@ class WorkoutDetector { ); kcal = cal.kcal; kj = cal.kj; + calUsedDefaultAnchors = cal.usedDefaultAnchors; } final avg = winBpm.reduce((a, b) => a + b) / winBpm.length; final peak = winBpm.reduce(math.max).round(); // Strain via the existing StrainScorer (reused, NOT re-derived). - final strain = StrainScorer.strain( - winBpm, - [for (final t in winTs) t.toDouble()], - maxHR: effMaxHR, - restingHR: restHR, - ); + // + // NO HIDDEN ANCHOR: StrainScorer.strain silently substitutes + // `defaultMaxHR() = 220 − 30 = 190` when maxHR is null, so a session used + // to report `hrmax: null, hrmax_source: 'unknown'` next to a concrete + // strain scored against an invented 190. We ABSTAIN instead — no anchor, + // no strain. (The zone gate above already drops such bouts; this keeps the + // guarantee local so it survives any future change to that gate. The + // fallback itself lives in clinical/load_trimp.dart and is not ours to + // change.) + final strain = effMaxHR == null + ? null + : StrainScorer.strain( + winBpm, + [for (final t in winTs) t.toDouble()], + maxHR: effMaxHR, + restingHR: restHR, + ); // HYBRID SEAM: type the bout. final bout = WorkoutBout( @@ -428,6 +476,7 @@ class WorkoutDetector { hrmaxSource: hrmaxSource, caloriesKcal: kcal, caloriesKJ: kj, + caloriesUsedDefaultAnchors: calUsedDefaultAnchors, sport: sport, )); } @@ -465,12 +514,28 @@ Metric> detectWorkouts({ savedSpans: savedSpans, classify: classify, ); + // Distinguish "no workouts today" from "the zone-2 qualification gate could + // not be evaluated because there is no HRmax anchor" — with no anchor the + // detector correctly emits nothing, and the caller deserves to know why. + final noAnchor = + maxHR == null && StrainScorer.estimateHRmax(hrBpm, age).$1 == 0.0; return Metric>( value: list, confidence: list.isEmpty ? 0.0 : 0.6, tier: Tier.estimate, - inputs_used: const ['hr_1hz', 'gravity_1hz', 'profile'], - note: 'detected workouts (HR + motion gated, ≥5 min, ≥50% time in zone 2+); ' - 'wrist-HR ESTIMATE, not medical advice', + inputs_used: [ + 'hr_1hz', + 'gravity_1hz', + if (profile != null) 'profile', + if (maxHR != null) 'max_hr', + if (age != null) 'age', + if (restingHR != null) 'resting_hr', + ], + note: noAnchor + ? 'no HRmax anchor (no caller HRmax, no age, too few HR samples for an ' + 'observed estimate) — the ≥50% time-in-zone-2+ qualification gate ' + 'cannot be evaluated, so NO workout is emitted (never guessed)' + : 'detected workouts (HR + motion gated, ≥5 min, ≥50% time in zone 2+); ' + 'wrist-HR ESTIMATE, not medical advice', ); } diff --git a/test/onehz/clinical_test.dart b/test/onehz/clinical_test.dart index aea022e..606feb1 100644 --- a/test/onehz/clinical_test.dart +++ b/test/onehz/clinical_test.dart @@ -57,6 +57,29 @@ void main() { expect(m.value!.hf, isNull); // HF withheld honestly expect(m.value!.lf, isNotNull); // LF still reported }); + + test('REGRESSION: a gated HF is not republished through `total`', () { + // total used to sum hfRaw back in, so the gated and ungated totals were + // bit-identical (0.49944 in both) — the suppression was cosmetic. + final rr = []; + final times = []; + var t = 0.0; + for (var i = 0; i < 400; i++) { + final v = 1000 + 40 * math.sin(2 * math.pi * 0.25 * (t / 1000)); + rr.add(v); + t += v; + times.add(t); + } + final clean = hrvFreq(rr, times, artifactFraction: 0.0); + final gated = hrvFreq(rr, times, artifactFraction: 0.5); + expect(clean.value!.total, isNotNull); + expect(gated.value!.hfGated, isTrue); + expect(gated.value!.total, isNull, + reason: 'total power is a sum over ALL bands; with HF withheld it is ' + 'not computable'); + expect(gated.value!.toJson().containsKey('total'), isFalse); + expect(gated.note, contains('total')); + }); }); group('PRSA DC/AC (Bauer 2006)', () { @@ -78,6 +101,20 @@ void main() { test('absent without enough beats', () { expect(decelerationCapacity([1000, 1010, 990]).present, isFalse); }); + test('REGRESSION: l=1 is refused, not a RangeError', () { + // The Haar contrast at wavelet scale s=2 reads X(-2) = profile[l-2], + // i.e. profile[-1] for l=1 — it used to throw RangeError. + final rr = [ + for (var i = 0; i < 200; i++) 1000 + 20 * math.sin(2 * math.pi * i / 20) + ]; + final dc = decelerationCapacity(rr, l: 1); + expect(dc.present, isFalse); + expect(dc.confidence, 0); + expect(dc.note, contains('l≥2')); + expect(accelerationCapacity(rr, l: 1).present, isFalse); + // l=2 (the default) still works on the same series. + expect(decelerationCapacity(rr, l: 2).present, isTrue); + }); }); group('nocturnal RHR + dip', () { @@ -96,15 +133,56 @@ void main() { expect(m.value!.p1, closeTo(50, 1.0)); }); test('dip band classification', () { - final day = [for (var i = 0; i < 200; i++) 80]; - final night = [for (var i = 0; i < 200; i++) 60]; + final day = [for (var i = 0; i < 400; i++) 80]; + final night = [for (var i = 0; i < 400; i++) 60]; final m = hrDip(day, night); expect(m.value!.dipPct, closeTo(25, 1e-9)); // (80-60)/80 expect(m.value!.band, 'dipper'); - // riser case - final r = hrDip([60, 60, 60], [70, 70, 70]); + // riser case. NOTE: this used to use 3 samples a side; a 3-sample "day" + // and "night" is exactly the fabrication hrDipMinSamples now refuses, so + // the case is expressed with a real (5-min) period of coverage instead. + final r = hrDip([for (var i = 0; i < 400; i++) 60], + [for (var i = 0; i < 400; i++) 70]); expect(r.value!.band, 'riser'); }); + + test('REGRESSION: lowest-30-min mean needs a REAL contiguous 30-min window', + () { + // 900 valid samples ramping 100 -> 50 bpm. The gate admitted + // length >= windowSamples~/2 and then w = min(1800, 900) made the sliding + // loop never execute, so low30Mean was the WHOLE-STREAM mean (75.0) — + // published as a "lowest-30-min" trough with confidence 0.4. + final ramp = [for (var i = 0; i < 900; i++) 100 - i * 50 / 900]; + final m = nocturnalRhr(ramp); + expect(m.present, isFalse); + expect(m.value, isNull); + expect(m.confidence, 0); + }); + + test('REGRESSION: off-skin gaps are not compacted into a fake window', () { + // 1800 valid samples scattered one-per-16-s across 8 h. Compacting the + // valid stream made this a "30-min" window that actually spans 8 hours. + final scattered = [ + for (var i = 0; i < 28800; i++) (i % 16 == 0) ? 60.0 : 0.0 + ]; + expect(nocturnalRhr(scattered).present, isFalse); + // A genuinely contiguous 30-min block of the same samples DOES resolve. + final contiguous = [ + ...List.filled(1000, 0), + ...List.filled(1800, 60.0), + ...List.filled(1000, 0), + ]; + final ok = nocturnalRhr(contiguous); + expect(ok.present, isTrue); + expect(ok.value!.low30Mean, closeTo(60, 1e-9)); + }); + + test('REGRESSION: hrDip refuses a 1-sample day and a 1-sample night', () { + final m = hrDip([70], [60]); + expect(m.present, isFalse); + expect(m.confidence, 0); + expect(m.note, contains('${hrDipMinSamples}')); + }); }); group('illness CUSUM FSM (NightSignal)', () { @@ -139,6 +217,39 @@ void main() { expect(out.every((d) => d.state == IllnessState.green), isTrue); expect(out.every((d) => d.cusum == null), isTrue); }); + test('REGRESSION: a DEGENERATE (zero-dispersion) baseline abstains instead ' + 'of standardizing against a fabricated 1 bpm scale', () { + // 9 identical quantized nights, then a 5 bpm bump. MAD = 0 AND SD = 0, so + // there is no dispersion at all. The old `max(1.0, SD)` fallback made + // scale = 1 bpm => z = 5 => cusum 4.5 > h=4 => yellow, red by night 12: + // a one-night bump latching a sustained "illness" red. + final rhr = [...List.filled(9, 55.0), 60.0, 58.0, 58.0]; + final dates = [for (var i = 0; i < rhr.length; i++) 'd$i']; + final out = illnessCusum(dates, rhr); + // Before: green×9, then yellow, then RED, RED. + expect(out.every((d) => d.state == IllnessState.green), isTrue, + reason: 'no alarm can be raised without a dispersion estimate'); + // Night 9 has a long-enough baseline that is perfectly constant + // (MAD = 0 AND SD = 0) => abstained, not standardized against 1 bpm. + expect(out[9].cusum, isNull); + expect(out[9].z, isNull); + expect(out[9].need, degenerateBaselineNote); + // Nights 10-11 gain a real SD once the 60 enters the window, so they are + // evaluated — but from an honest scale, and they never trip the alarm. + expect(out[10].z, isNotNull); + expect(out[10].cusum!, lessThan(4.0)); + }); + test('a merely QUANTIZED baseline (MAD=0 but SD>0) still evaluates', () { + // MAD collapses on this baseline but SD does not — same convention as + // wellness/readiness_composite.dart: fall back to SD, only abstain when + // BOTH are zero. + final rhr = [...List.filled(8, 55.0), 56.0, 60.0]; + final dates = [for (var i = 0; i < rhr.length; i++) 'd$i']; + final out = illnessCusum(dates, rhr); + expect(out.last.z, isNotNull); + expect(out.last.cusum, isNotNull); + expect(out.last.need, isNull); + }); }); group('lnRMSSD readiness stack', () { @@ -162,6 +273,23 @@ void main() { expect(m.present, isTrue); expect(m.value!.rolling7Mean, 4.0); }); + test('REGRESSION: an UNDEFINED baseline SD abstains instead of emitting ' + 'cvPct 0.0 / band "normal"', () { + // One prior night => stddev() is null => CV, SWC and the band are + // undefined. The metric used to publish cvPct 0.0, swc null and + // band 'normal' anyway: "tonight is typical", asserted from nothing. + final m = readinessLnRmssd([4.0, 2.0], minNights: 2); + expect(m.present, isFalse); + expect(m.value, isNull); + expect(m.confidence, 0); + expect(m.note, contains('dispersion undefined')); + // A DEFINED (even zero) dispersion still computes — CV really is 0 there. + final flat = readinessLnRmssd([4.0, 4.0, 4.0, 2.0]); + expect(flat.present, isTrue); + expect(flat.value!.cvPct, 0.0); + expect(flat.value!.z, isNull); // z is undefined at SD = 0 + expect(flat.value!.band, 'suppressed'); + }); }); group('cosinor', () { @@ -187,6 +315,39 @@ void main() { final m = cosinor(t, y); expect(m.value!.r2, lessThan(0.2)); }); + test('REGRESSION: 4 random points never score a confident circadian fit', + () { + // A 3-parameter fit (M, β, γ) on 4 points has ONE residual degree of + // freedom: 4 random points scored raw r² 0.76–0.99 and were published at + // confidence 0.95, tier HIGH. Now: refused outright (< cosinorMinPoints). + final rnd = math.Random(20260726); + for (var trial = 0; trial < 200; trial++) { + final t = [0, 6, 12, 18]; + final y = [for (var i = 0; i < 4; i++) rnd.nextDouble()]; + final m = cosinor(t, y); + expect(m.present, isFalse); + expect(m.confidence, 0); + } + }); + test('REGRESSION: confidence comes from the ADJUSTED R² (3 fitted params)', + () { + // 8 points of noise: raw R² is upward-biased, adjusted R² is not. + final rnd = math.Random(4242); + final t = [for (var i = 0; i < 8; i++) i * 3.0]; + final y = [for (var i = 0; i < 8; i++) rnd.nextDouble()]; + final m = cosinor(t, y); + expect(m.present, isTrue); + expect(m.value!.r2Adj, lessThan(m.value!.r2)); + expect(m.confidence, closeTo(m.value!.r2Adj.clamp(0.1, 0.95), 1e-12)); + // A genuine 24-h rhythm still earns full confidence. + final tt = [for (var h = 0; h < 48; h++) h.toDouble()]; + final yy = [ + for (var h = 0; h < 48; h++) 60 + 10 * math.cos(2 * math.pi * h / 24) + ]; + final good = cosinor(tt, yy); + expect(good.value!.r2Adj, closeTo(1.0, 1e-6)); + expect(good.confidence, 0.95); + }); }); group('TRIMP + CTL/ATL/TSB', () { @@ -197,6 +358,38 @@ void main() { restingHr: 50, maxHr: 190, sex: Sex.male); expect(m.value!, greaterThan(0)); }); + + test('REGRESSION: ONE Banister implementation, matching the published ' + 'sex-specific y = c·e^(b·x)', () { + // The two implementations in this file disagreed by 1.5625× (the + // top-level one dropped the 0.64/0.86 coefficient entirely) and the + // StrainScorer one applied the MALE 0.64 to women (−26% on female load). + const x = (150.0 - 50.0) / (190.0 - 50.0); // 0.714286 %HRR + final expectedMale = x * 0.64 * math.exp(1.92 * x); + final expectedFemale = x * 0.86 * math.exp(1.67 * x); + + final male = + banisterTrimp([150], restingHr: 50, maxHr: 190, sex: Sex.male); + final female = + banisterTrimp([150], restingHr: 50, maxHr: 190, sex: Sex.female); + expect(male.value!, closeTo(expectedMale, 1e-9)); + expect(female.value!, closeTo(expectedFemale, 1e-9)); + // 1 min at 150 bpm (RHR 50, HRmax 190): 1.8016 male / 2.0250 female. + // Before: 2.8150 from the top-level fn (no c at all) and 1.5070 from + // StrainScorer for a woman (male c on the female b). + expect(male.value!, closeTo(1.801589, 1e-5)); + expect(female.value!, closeTo(2.024984, 1e-5)); + + // The StrainScorer path agrees exactly with the top-level one. + expect(StrainScorer.banisterTRIMP([150], 50, 140, [1.0], female: false), + closeTo(male.value!, 1e-12)); + expect(StrainScorer.banisterTRIMP([150], 50, 140, [1.0], female: true), + closeTo(female.value!, 1e-12)); + expect(StrainScorer.banisterY(x, female: false), + closeTo(0.64 * math.exp(1.92 * x), 1e-12)); + expect(StrainScorer.banisterY(x, female: true), + closeTo(0.86 * math.exp(1.67 * x), 1e-12)); + }); test('Edwards zone-sum is the weighted dot product', () { // zones [10,5,0,0,0] -> 10*1 + 5*2 = 20 final m = edwardsTrimp([10, 5, 0, 0, 0]); @@ -214,6 +407,35 @@ void main() { expect(s.value!.atl, greaterThan(s.value!.ctl)); expect(s.value!.tsb, lessThan(0)); }); + + test('REGRESSION: one training day does NOT fabricate 42 days of chronic ' + 'load', () { + // ctlAtlTsb([500]) used to seed BOTH accumulators at dailyTrimp.first → + // ctl 500, atl 500, tsb 0.0: a fully-adapted, perfectly-fresh athlete + // conjured from a single workout. + final one = ctlAtlTsb([500.0]); + expect(one.present, isFalse); + expect(one.value, isNull); + expect(one.confidence, 0); + expect(one.note, 'need_baseline:have=1,need=$ctlAtlMinDays'); + // Still absent one day short of the minimum... + expect(ctlAtlTsb(List.filled(ctlAtlMinDays - 1, 50.0)).present, + isFalse); + // ...and present at the minimum. + expect( + ctlAtlTsb(List.filled(ctlAtlMinDays, 50.0)).present, isTrue); + }); + + test('REGRESSION: the seed is a week of observed load, not day one', () { + // A single huge opening day must not become the chronic baseline. + final hist = [600.0, for (var i = 0; i < 20; i++) 0.0]; + final m = ctlAtlTsb(hist); + expect(m.present, isTrue); + // Prime = mean of the first 7 days = 600/7 ≈ 85.7, then 14 rest days + // decay it — nowhere near the old ctl≈600 anchor. + expect(m.value!.ctl, lessThan(90)); + expect(m.value!.atl, lessThan(m.value!.ctl)); + }); }); group('display heart-rate zones', () { @@ -388,6 +610,54 @@ void main() { expect((v * 100).round() / 100, v); }); + test('REGRESSION: trimpToStrain is CLAMPED to maxStrain', () { + // Docstring says "Map accumulated TRIMP onto [0, 100]" but nothing + // clamped: 14400 → 107.8. (The sibling strainScore() always clamped.) + expect(StrainScorer.trimpToStrain(14400), 100.0); + expect(StrainScorer.trimpToStrain(1e9), StrainScorer.maxStrain); + expect(StrainScorer.trimpToStrain(7200), closeTo(100.0, 1e-9)); + // Below the ceiling nothing changed. + expect(StrainScorer.trimpToStrain(335), lessThan(100.0)); + }); + + test('REGRESSION: strain integrates PER-SAMPLE durations, not the first ' + 'inter-sample gap applied to everything', () { + const bpmv = 150.0; + // (a) 21 samples over 20 min whose FIRST two are 1 s apart — exactly the + // sparse stream minSparseReadings admits. sampleDuration was 1 s for all + // 21 samples → strain 8.08 instead of ~47. + final tsIrregular = [0, 1, for (var i = 1; i < 20; i++) 1 + i * 63.1]; + final bpm21 = List.filled(21, bpmv); + final irregular = + StrainScorer.strain(bpm21, tsIrregular, maxHR: 190, restingHR: 50)!; + final uniform = StrainScorer.strain( + bpm21, [for (var i = 0; i < 21; i++) i * 60.0], + maxHR: 190, restingHR: 50)!; + expect(irregular, greaterThan(40.0)); + expect(irregular, closeTo(uniform, 3.0), + reason: 'same HR over the same wall-clock span → similar strain'); + + // (b) The inverse: 700 samples at 1 Hz behind a 300 s leading gap. The + // 5-min first gap became every sample's duration → strain 104.25. + final tsGap = [0, for (var i = 0; i < 699; i++) 300.0 + i]; + final gapped = StrainScorer.strain( + List.filled(700, bpmv), tsGap, + maxHR: 190, restingHR: 50)!; + final dense = StrainScorer.strain(List.filled(700, bpmv), + [for (var i = 0; i < 700; i++) i.toDouble()], + maxHR: 190, restingHR: 50)!; + expect(gapped, lessThanOrEqualTo(StrainScorer.maxStrain)); + expect(gapped, closeTo(dense, 1.0), + reason: 'a hole in the stream is not elapsed effort'); + expect(gapped, lessThan(60.0)); + + // Per-sample durations: gaps capped at the median cadence, tail gets it. + final durs = StrainScorer.sampleDurationsMinutes(tsGap); + expect(durs.length, 700); + expect(durs.first, closeTo(1 / 60.0, 1e-12)); + expect(durs.reduce(math.max), closeTo(1 / 60.0, 1e-12)); + }); + test('Edwards zone weight at %HRR boundaries (RHR=0,reserve=100 → bpm=%HRR)', () { int w(double pct) => StrainScorer.zoneWeight(pct, 0, 100); expect(w(49), 0); @@ -404,12 +674,11 @@ void main() { final lo = List.filled(30, 100.0); final hi = List.filled(30, 150.0); final ts = [for (var i = 0; i < 30; i++) i.toDouble()]; - final tLo = StrainScorer.banisterTRIMP( - lo, 50, 150, StrainScorer.sampleDurationMinutes(ts), - StrainScorer.banisterBMen); - final tHi = StrainScorer.banisterTRIMP( - hi, 50, 150, StrainScorer.sampleDurationMinutes(ts), - StrainScorer.banisterBMen); + // API change: TRIMP now integrates PER-SAMPLE durations, and the + // Banister sex is selected by name (b and its scale must stay paired). + final durs = StrainScorer.sampleDurationsMinutes(ts); + final tLo = StrainScorer.banisterTRIMP(lo, 50, 150, durs); + final tHi = StrainScorer.banisterTRIMP(hi, 50, 150, durs); expect(tHi, greaterThan(tLo)); }); @@ -536,6 +805,26 @@ void main() { expect(m.present, isFalse); expect(m.value, isNull); }); + + test('REGRESSION: a NEAR-degenerate RR range abstains instead of reporting ' + 'SI 48780 / "high"', () { + // 300 beats alternating 1000/1001 ms — plausible 1 Hz beat-timing + // quantization at a steady sleeping HR. The guard was only mxdmnS <= 0, + // so MxDMn = 0.001 s blew the 1/MxDMn denominator up to si 48780, + // level 'high'. + final nn = [ + for (var i = 0; i < 300; i++) i.isEven ? 1000.0 : 1001.0 + ]; + final m = baevskyStressIndex(nn); + expect(m.present, isFalse); + expect(m.value, isNull); + expect(m.confidence, 0); + // A genuinely varying series of the same length still computes. + final ok = baevskyStressIndex([ + for (var i = 0; i < 300; i++) 900.0 + 15.0 * math.sin(i.toDouble()) + ]); + expect(ok.present, isTrue); + }); }); group('cardiac coherence (McCraty & Zayas 2014)', () { diff --git a/test/onehz/coaching_test.dart b/test/onehz/coaching_test.dart index dfa2abc..c0e3cd9 100644 --- a/test/onehz/coaching_test.dart +++ b/test/onehz/coaching_test.dart @@ -1,5 +1,7 @@ // Coaching surface — synthetic known-answer tests, incl. a regression for the // physiological-age oversleep bug. Covers PR #11's untested coaching API. +import 'dart:convert'; + import 'package:test/test.dart'; import 'package:openstrap_analytics/src/onehz/types.dart'; import 'package:openstrap_analytics/src/onehz/human/coaching.dart'; @@ -467,4 +469,193 @@ void main() { expect(m.value, isEmpty); }); }); + + // ------------------------------------------------------------------------- + // REGRESSION: physiologicalAge must ABSTAIN with no physiology, and must + // report the inputs it ACTUALLY used. + // ------------------------------------------------------------------------- + group('physiologicalAge — honesty envelope (regression)', () { + test('every physiological input null => ABSENT, not "your age"', () { + // PRE-FIX: score started at chronologicalAge, nothing moved it, and the + // function returned a PRESENT metric (physioAge 30, delta 0, conf 0.35) + // claiming six inputs it had never seen. + final m = physiologicalAge( + chronologicalAge: 30, + sex: Sex.male, + vo2max: null, + restingHr: null, + rmssd: null, + sleepDurationH: null, + sleepEfficiency: null, + dailySteps: null, + ); + expect(m.present, isFalse); + expect(m.value, isNull); + expect(m.confidence, 0); + expect(m.toJson()['value'], '—'); + expect(m.inputs_used, ['profile']); + }); + + test('inputs_used lists only the inputs actually supplied', () { + // PRE-FIX this was a hardcoded six-entry list in EVERY partial case. + final m = physiologicalAge( + chronologicalAge: 30, + sex: Sex.male, + vo2max: null, + restingHr: 55, + rmssd: null, + sleepDurationH: 7.5, + sleepEfficiency: null, + dailySteps: null, + ); + expect(m.present, isTrue); + expect(m.inputs_used, ['profile', 'resting_hr', 'sleep_duration']); + expect(m.inputs_used, isNot(contains('vo2max'))); + expect(m.inputs_used, isNot(contains('rmssd'))); + expect(m.inputs_used, isNot(contains('steps'))); + }); + + test('confidence scales with how much physiology went in', () { + Metric build(int n) => physiologicalAge( + chronologicalAge: 40, + sex: Sex.male, + vo2max: n >= 1 ? 50 : null, + restingHr: n >= 2 ? 48 : null, + rmssd: n >= 3 ? 60 : null, + sleepDurationH: n >= 4 ? 7.5 : null, + sleepEfficiency: n >= 5 ? 94 : null, + dailySteps: n >= 6 ? 12000 : null, + ); + expect(build(6).confidence, greaterThan(build(1).confidence)); + expect(build(6).inputs_used, hasLength(7)); // profile + 6 + }); + }); + + // ------------------------------------------------------------------------- + // REGRESSION: vo2maxEstimate must not divide by a zero resting HR. + // ------------------------------------------------------------------------- + group('vo2maxEstimate — zero resting HR (regression)', () { + test('restingHr == 0 (the off-skin sentinel) ABSTAINS, never Infinity', () { + // PRE-FIX `maxHr <= restingHr` did not catch it: 15.3 * (190/0) produced + // value: Infinity, which Metric.toJson emits raw and jsonEncode throws on. + final m = vo2maxEstimate(restingHr: 0, maxHr: 190, sex: Sex.male, age: 30); + expect(m.present, isFalse); + expect(m.value, isNull); + expect(() => jsonEncode(m.toJson()), returnsNormally); + }); + + test('a negative or non-finite resting HR also abstains', () { + expect( + vo2maxEstimate(restingHr: -5, maxHr: 190, sex: Sex.male, age: 30) + .present, + isFalse); + expect( + vo2maxEstimate( + restingHr: double.nan, maxHr: 190, sex: Sex.male, age: 30) + .present, + isFalse); + }); + + test('a valid pair still computes', () { + final m = vo2maxEstimate(restingHr: 50, maxHr: 190, sex: Sex.male, age: 30); + expect(m.present, isTrue); + expect(m.value!.isFinite, isTrue); + expect(() => jsonEncode(m.toJson()), returnsNormally); + }); + }); + + // ------------------------------------------------------------------------- + // REGRESSION: journalCorrelations needs a dispersion test, and must not + // index an outcome list by dates.length without checking. + // ------------------------------------------------------------------------- + group('journalCorrelations — dispersion + length guard (regression)', () { + test('a 3% mean gap swamped by within-group spread is NOT meaningful', () { + // tagged [50,80] mean 65 vs untagged [40,86] mean 63 => +3.17%, which + // PRE-FIX cleared the bare `pct.abs() >= 3.0` bar. Each side spans 30–46 + // points, so Cohen's d is ~0.07: this is noise, not a journal effect. + final journal = [ + const JournalDay('d0', {'coffee'}), + const JournalDay('d1', {'coffee'}), + const JournalDay('d2', {}), + const JournalDay('d3', {}), + ]; + final out = journalCorrelations( + journal: journal, + dates: const ['d0', 'd1', 'd2', 'd3'], + outcomes: const { + 'recovery': [50, 80, 40, 86] + }, + ); + final eff = out.firstWhere((c) => c.tag == 'coffee').effects.single; + expect(eff.insufficient, isFalse); + expect(eff.pctChange!.abs(), greaterThanOrEqualTo(3.0), + reason: 'the old percentage bar IS cleared'); + expect(eff.cohensD, isNotNull); + expect(eff.cohensD!.abs(), lessThan(0.5)); + expect(eff.meaningful, isFalse, + reason: 'dispersion test must veto it (d=${eff.cohensD})'); + }); + + test('a large, well-separated effect is still meaningful', () { + final out = journalCorrelations( + journal: const [ + JournalDay('d0', {'alcohol'}), + JournalDay('d1', {'alcohol'}), + JournalDay('d2', {}), + JournalDay('d3', {}), + ], + dates: const ['d0', 'd1', 'd2', 'd3'], + outcomes: const { + 'recovery': [40, 42, 80, 82] + }, + ); + final eff = out.firstWhere((c) => c.tag == 'alcohol').effects.single; + expect(eff.meaningful, isTrue); + expect(eff.cohensD!.abs(), greaterThan(0.5)); + }); + + test('two constant sides with only 2 days each are NOT meaningful', () { + // Pooled SD is 0 so Cohen's d is undefined; refuse to call it. + final out = journalCorrelations( + journal: const [ + JournalDay('d0', {'x'}), + JournalDay('d1', {'x'}), + JournalDay('d2', {}), + JournalDay('d3', {}), + ], + dates: const ['d0', 'd1', 'd2', 'd3'], + outcomes: const { + 'recovery': [60, 60, 70, 70] + }, + ); + final eff = out.firstWhere((c) => c.tag == 'x').effects.single; + expect(eff.cohensD, isNull); + expect(eff.meaningful, isFalse); + }); + + test('an outcome list shorter than dates is guarded, not a RangeError', () { + // PRE-FIX `entry.value[i]` was indexed by dates.length => RangeError. + late final List out; + expect( + () => out = journalCorrelations( + journal: const [ + JournalDay('d0', {'x'}), + JournalDay('d1', {'x'}), + JournalDay('d2', {}), + JournalDay('d3', {}), + ], + dates: const ['d0', 'd1', 'd2', 'd3'], + outcomes: const { + 'recovery': [60, 62] // misaligned: 2 values for 4 dates + }, + ), + returnsNormally, + ); + final eff = out.firstWhere((c) => c.tag == 'x').effects.single; + expect(eff.insufficient, isTrue); + expect(eff.meaningful, isFalse); + expect(eff.nTagged, 0); + expect(eff.nUntagged, 0); + }); + }); } diff --git a/test/onehz/foundations_test.dart b/test/onehz/foundations_test.dart index 44a839b..1199437 100644 --- a/test/onehz/foundations_test.dart +++ b/test/onehz/foundations_test.dart @@ -1,4 +1,5 @@ // Item 2 — FOUNDATIONS. Synthetic, known-answer tests. +import 'dart:math' as math; import 'package:test/test.dart'; import 'package:openstrap_analytics/onehz.dart'; @@ -95,7 +96,64 @@ void main() { final r = correctRr(rr); expect(r.cleanFraction, closeTo(1.0, 1e-9)); expect(r.droppedCount, 0); + expect(r.correctedCount, 0); expect(r.nn.length, 60); + // Nothing was substituted: the cleaned series IS the input. + for (var i = 0; i < rr.length; i++) { + expect(r.nn[i], closeTo(rr[i], 1e-12)); + } + }); + + // 400 beats of ORDINARY resting variability: RSA at ~13 beats/breath + + // a slow LF wave + a little jitter. RR 928-1172 ms, max |dRR| 56 ms, + // ZERO injected artifacts. + List cleanRsa() { + final rnd = math.Random(11); + return [ + for (var i = 0; i < 400; i++) + 1050 + + 95 * math.sin(2 * math.pi * i / 13.0) + + 25 * math.sin(2 * math.pi * i / 61.0) + + (rnd.nextDouble() - 0.5) * 10 + ]; + } + + test('REGRESSION: a clean physiological RSA record is NOT flagged — the ' + 'quartile deviation is taken on the SIGNED dRR series ' + '(Lipponen-Tarvainen 2019)', () { + // Taking the QD of |dRR| folds the symmetric ±dRR distribution onto one + // side, collapsing the dispersion so far that the threshold sinks to the + // minThresholdMs floor and the detector degenerates into a fixed 100 ms + // cut-off. On THIS artifact-free record that flagged 32 of 400 healthy + // beats (cleanFraction 0.92) and shrank SDNN 69.82 -> 64.78 (-7%). + final rr = cleanRsa(); + var maxAbsDrr = 0.0; + for (var i = 1; i < rr.length; i++) { + final d = (rr[i] - rr[i - 1]).abs(); + if (d > maxAbsDrr) maxAbsDrr = d; + } + expect(maxAbsDrr, lessThan(100), + reason: 'sanity: every beat-to-beat step is below the 100 ms floor'); + + final r = correctRr(rr); + expect(r.cleanFraction, 1.0); + expect(r.correctedCount, 0); + expect(r.droppedCount, 0); + expect(r.nn.length, 400); + // The HRV of a clean record must survive correction untouched. + final before = hrvTime(rr).value!; + final after = hrvTime(r.nn).value!; + expect(after.rmssd!, closeTo(before.rmssd!, 1e-9)); + expect(after.sdnn!, closeTo(before.sdnn!, 1e-9)); + }); + + test('a genuine gross outlier is still caught on that same record', () { + // The signed-QD threshold must not have blinded the detector. + final rr = cleanRsa(); + rr[200] = 350; // impossible beat + final r = correctRr(rr); + expect(r.classes[200], isNot(BeatClass.normal)); + expect(r.cleanFraction, lessThan(1.0)); }); test('flags EXACTLY one injected isolated ectopic and spline-corrects it', () { @@ -182,6 +240,23 @@ void main() { expect(e.last.value, lessThan(20)); expect(e.last.value, greaterThan(10)); }); + test('REGRESSION: gap-aware EWMA never extrapolates outside the data on a ' + 'duplicate or non-monotonic timestamp', () { + // dt <= 0 made lambda = 1 - 2^(-dt/H) NEGATIVE, so the update ran + // BACKWARDS: [0,1000,500] with values [10,10,20] produced 5.857 — below + // every input (all >= 10). + final e = gapAwareEwma([0, 1000, 500], [10, 10, 20], halfLifeMs: 1000); + expect(e.length, 3); + for (final p in e) { + expect(p.value, greaterThanOrEqualTo(10.0)); + expect(p.value, lessThanOrEqualTo(20.0)); + } + // No elapsed time => no new weight => the estimate does not move. + expect(e.last.value, closeTo(10.0, 1e-12)); + // Duplicate timestamps behave the same way. + final dup = gapAwareEwma([0, 0, 0], [10, 50, 50], halfLifeMs: 1000); + expect(dup.map((p) => p.value), everyElement(closeTo(10.0, 1e-12))); + }); test('MDC gate: small change suppressed, large surfaced', () { final b = robustBaseline([10, 11, 9, 10, 12, 8, 10, 11, 9, 10]); expect(changeExceedsMdc(0.1, b), isFalse); diff --git a/test/onehz/human_test.dart b/test/onehz/human_test.dart index ebf91ea..a8399df 100644 --- a/test/onehz/human_test.dart +++ b/test/onehz/human_test.dart @@ -373,4 +373,63 @@ void main() { expect(p.value!.percentile, inInclusiveRange(0, 100)); }); }); + + // ------------------------------------------------------------------------- + // REGRESSION: mid-sleep is CIRCULAR. A midnight-straddling weekday midpoint + // must not report ~22.5 h of social jetlag for a ~1.4 h drift. + // [PUB Wittmann/Roenneberg, Chronobiol Int 2006; MCTQ MSFsc] + // ------------------------------------------------------------------------- + group('social jetlag / chronotype — circular clock (regression)', () { + test('midnight-straddling weekday mid-sleep gives the SHORT arc', () { + // ~23:50 weekday mid-sleep vs ~01:10 weekend => a real +1.4 h drift. + // PRE-FIX: median() on raw clock-hours gave msw 23.7, msf 1.15 and a + // headline sjl of −22.55 h. + final work = [23.7, 0.3, 23.9, 0.1, 23.8]; + final free = [1.0, 1.3]; + final m = socialJetlag(free, work); + expect(m.present, isTrue, reason: m.note); + final v = m.value!; + expect(v.sjlHours, greaterThan(0), reason: 'weekend runs later'); + expect(v.sjlHours, closeTo(1.25, 0.4)); + expect(v.absHours, lessThan(2.0)); + }); + + test('|SJL| can never exceed 12 h (shortest arc on the 24 h circle)', () { + // ~23:06 free vs ~01:06 work: 2 h EARLIER, not 22 h later. + final m = socialJetlag([23.0, 23.2], [1.0, 1.2, 1.1]); + expect(m.present, isTrue, reason: m.note); + expect(m.value!.absHours, closeTo(2.0, 0.3)); + expect(m.value!.absHours, lessThanOrEqualTo(12.0)); + expect(m.value!.sjlHours, lessThan(0), reason: 'weekend runs EARLIER'); + }); + + test('the ordinary (non-wrapping) case is unchanged', () { + final m = socialJetlag([5.4, 5.6, 5.5], [3.4, 3.6, 3.5, 3.5, 3.4]); + expect(m.value!.sjlHours, closeTo(2.0, 0.15)); + }); + + test('chronotype: a 23:30 free-day mid-sleep is EARLY, not evening', () { + // PRE-FIX msfSc = 23.5 fell through every band to "evening type" — the + // exact opposite of the truth. The label bands now read the MCTQ band + // axis (clock-hours unwrapped about 06:00). + final m = chronotype( + [23.4, 23.6, 23.5], + [8.0, 8.0, 8.0], + avgWeekSleepDurH: 8.0, + totalDaysObserved: 21, + ); + expect(m.present, isTrue, reason: m.note); + expect(m.value!.typeLabel, 'early type'); + }); + + test('chronotype: a late free-day mid-sleep is still an evening type', () { + final m = chronotype( + [5.3, 5.6, 5.5], + [8.5, 9.0, 8.8], + avgWeekSleepDurH: 8.7, + totalDaysObserved: 21, + ); + expect(m.value!.typeLabel, contains('evening')); + }); + }); } diff --git a/test/onehz/respiration_test.dart b/test/onehz/respiration_test.dart index 246c106..2e71824 100644 --- a/test/onehz/respiration_test.dart +++ b/test/onehz/respiration_test.dart @@ -378,4 +378,40 @@ void main() { } }); }); + + // ------------------------------------------------------------------------- + // REGRESSION: brpm, peak_hz and power inside one RespEstimate must come from + // ONE source. brpm used to be median(peaks) across the three spectral grids + // while peakHz/power came from the highest-POWER grid, so `peak_hz * 60` and + // `brpm` disagreed inside a single result. + // ------------------------------------------------------------------------- + group('rsaRespRate — internally consistent output (regression)', () { + for (final hz in const [0.25, 0.20, 0.30]) { + test('modHz=$hz: peak_hz * 60 == brpm exactly', () { + final s = syntheticRsaRr(modHz: hz, beats: 500); + final corr = correctRr(s.rr); + final m = rsaRespRate(corr.nn, corr.nnTimesMs, + artifactFraction: 1 - corr.cleanFraction); + expect(m.present, isTrue, reason: m.note); + final v = m.value!; + expect(v.brpm, isNotNull); + expect(v.peakHz, isNotNull); + expect(v.power, isNotNull); + expect(v.peakHz! * 60.0, closeTo(v.brpm!, 1e-9), + reason: 'brpm ${v.brpm} vs peak_hz*60 ${v.peakHz! * 60}'); + // The reported rate must still be the right one. + expect(v.brpm!, closeTo(hz * 60.0, 1.5)); + }); + } + + test('the JSON pair round-trips consistently', () { + final s = syntheticRsaRr(modHz: 0.25, beats: 500); + final corr = correctRr(s.rr); + final m = rsaRespRate(corr.nn, corr.nnTimesMs, + artifactFraction: 1 - corr.cleanFraction); + final j = m.value!.toJson(); + expect((j['peak_hz'] as double) * 60.0, + closeTo(j['brpm'] as double, 1e-4)); + }); + }); } diff --git a/test/onehz/sleep_honesty_test.dart b/test/onehz/sleep_honesty_test.dart new file mode 100644 index 0000000..fdb1abb --- /dev/null +++ b/test/onehz/sleep_honesty_test.dart @@ -0,0 +1,546 @@ +// SLEEP — HONESTY REGRESSIONS (2026-07). +// +// The core contract: absent input yields null / unstaged, NEVER a fabricated +// value. `AdvancedSleepStager.stageWindow`'s own docstring promises "Seconds +// with no data ... simply stay unstaged (wake) — honest about gaps, never +// fabricated". Every test here pins a place where the sleep code broke that +// promise and reported a perfect night out of an empty or fragmented signal. +// +// Each test FAILS against the pre-fix behavior; the pre-fix number is stated in +// the test so a future reader can tell a regression from a re-tune. + +import 'dart:math' as math; +import 'package:test/test.dart'; +import 'package:openstrap_analytics/onehz.dart'; + +/// Fixed absolute epoch second so every fixture is deterministic. +const int _t0 = 1700000000; + +/// Midnight (UTC) of _t0's day — fixtures anchor local clock times off this and +/// always pass `tzOffsetSec: 0`, so "local" == UTC and nothing depends on the +/// machine timezone. +final int _midnight = _t0 - (_t0 % 86400); + +void main() { + // ═══════════════════════════════════════════════════════════════════════════ + // (1) A window with NO accelerometer at all must not be scored as sleep. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — no accelerometer in the window', () { + test('a forced window holding zero samples yields NO sleep ' + '(was: 8 h of light sleep, efficiency 100%)', () { + // 4 h of perfectly good data, then a user-asserted 8 h window ~13 h later + // that contains not one sample. + final accel = []; + final hr = []; + for (var i = 0; i < 4 * 3600; i++) { + accel.add(AccelSample((_t0 + i) * 1000.0, 0.0, 0.0, 1.0)); + hr.add(55); + } + final onset = _t0 + 50000; + final s = segmentSleep(accel, hr, + forcedWindow: (onsetSec: onset, offsetSec: onset + 8 * 3600)); + + // The window itself is the user's word, so it is honored... + expect(s.present, isTrue); + expect(s.inBedSec, 8 * 3600); + // ...but NOTHING inside it may be claimed as sleep. + expect(s.tstSec, 0, reason: 'no data ⇒ no sleep (pre-fix: 28800)'); + expect(s.lightSec, 0, reason: 'pre-fix: the whole window read "light"'); + expect(s.deepSec, 0); + expect(s.remSec, 0); + expect(s.wakeSec, 8 * 3600); + expect(s.efficiencyPct, 0.0, reason: 'pre-fix: 100.0'); + expect(s.stages.every((x) => x == SleepStage.wake), isTrue); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (2) HR entirely absent ⇒ abstain, never a fabricated 60 bpm baseline. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — heart rate entirely absent', () { + test('cardioStager abstains rather than defaulting the HR baseline to 60', + () { + // 2 h of perfectly still accel, HR off-skin (0) for every second. + final accel = [ + for (var i = 0; i < 2 * 3600; i++) + AccelSample((_t0 + i) * 1000.0, 0.0, 0.0, 1.0) + ]; + final hr = List.filled(2 * 3600, 0.0); + final r = cardioStager(hr, accel); + // Pre-fix: 240 epochs, every one NREM (the `?? 60` baseline meant the + // wake / REM / deep gates could never fire, so everything fell through). + expect(r.base.stages, isEmpty, reason: 'pre-fix: 240 NREM epochs'); + expect(r.confidence, 0); + // A window that DOES have HR still stages normally. + final withHr = List.filled(2 * 3600, 52.0); + expect(cardioStager(withHr, accel).base.stages, isNotEmpty); + }); + + test('the strap-on-the-nightstand night is not reported as a perfect sleep ' + '(was: TST 7h54, efficiency 100%)', () { + final accel = []; + final hr = []; + var i = 0; + final start = _midnight + 21 * 3600; // 21:00 + void seg(int secs, {required bool active}) { + for (var k = 0; k < secs; k++, i++) { + final x = active ? (k.isEven ? 0.0 : 0.3) : 0.005; + accel.add(AccelSample( + (start + i) * 1000.0, x, 0.0, active ? 0.95 : 1.0)); + hr.add(0.0); // NO heart rate at all — the band is off the wrist. + } + } + + seg(2 * 3600, active: true); + seg(8 * 3600, active: false); // 8 h of "perfect stillness" + seg(2 * 3600, active: true); + + final s = segmentSleep(accel, hr, tzOffsetSec: 0); + expect(s.tstSec, 0, reason: 'pre-fix: 28438 s of sleep from zero HR'); + expect(s.efficiencyPct, 0.0, reason: 'pre-fix: 100.0'); + expect(s.remSec, 0); + expect(s.deepSec, 0); + expect(s.lightSec, 0); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (3) An accelerometer dropout must not be carried forward into "stillness". + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — bounded accel carry-forward', () { + test('a 6 h dropout inside an 8 h window stays unstaged ' + '(was: TST 8 h, WASO 0, efficiency 100%)', () { + final accel = []; + final hr = []; + void block(int fromSec, int secs) { + for (var k = 0; k < secs; k++) { + accel.add(AccelSample((_t0 + fromSec + k) * 1000.0, 0.005, 0.0, 1.0)); + hr.add(52); + } + } + + block(0, 3600); // hour 1: real data + // hours 2-7: NOTHING + block(7 * 3600, 3600); // hour 8: real data + + final s = segmentSleep(accel, hr, + forcedWindow: (onsetSec: _t0, offsetSec: _t0 + 8 * 3600)); + expect(s.present, isTrue); + expect(s.inBedSec, 8 * 3600); + // Only the ~2 h that actually has data can be staged; the 6 h hole is + // wake. Allow the bounded 60 s carry-forward tail on each real block. + expect(s.tstSec!, lessThan(2 * 3600 + 2 * 61), + reason: 'pre-fix: 28800 — the whole window read as sleep'); + expect(s.wakeSec!, greaterThan(5 * 3600 + 45 * 60), + reason: 'pre-fix: 0'); + expect(s.efficiencyPct!, lessThan(30.0), reason: 'pre-fix: 100.0'); + }); + + test('a SHORT dropout is still carried forward (the bound is 60 s, not 0)', + () { + // Same 8 h window, but the hole is only 45 s — a plausible missed-sample + // burst, not a data outage. It must stay staged, or every real capture's + // packet loss would be punched out of the night. + final accel = []; + final hr = []; + for (var k = 0; k < 4 * 3600; k++) { + if (k >= 3600 && k < 3600 + 45) continue; // 45 s hole + accel.add(AccelSample((_t0 + k) * 1000.0, 0.005, 0.0, 1.0)); + hr.add(52); + } + final s = segmentSleep(accel, hr, + forcedWindow: (onsetSec: _t0, offsetSec: _t0 + 4 * 3600)); + expect(s.tstSec!, greaterThan((0.9 * 4 * 3600).round())); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (4) Webster rescore must score against an IMMUTABLE snapshot. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — Webster continuity rescore does not cascade', () { + /// `leadMin` of sleep, then [reps] × (`wakeMin` wake + 1 min sleep), then a + /// trailing wake block. At 30 s epochs. + List fragmented({ + required int leadMin, + required int wakeMin, + required int reps, + }) { + final out = []; + void push(SleepStage s, int mins) { + for (var k = 0; k < mins * 2; k++) { + out.add(s); + } + } + + push(SleepStage.nrem, leadMin); + for (var r = 0; r < reps; r++) { + push(SleepStage.wake, wakeMin); + push(SleepStage.nrem, 1); + } + push(SleepStage.wake, 20); + return out; + } + + int wakeInBody(List sm) { + var first = -1, last = -1; + for (var i = 0; i < sm.length; i++) { + if (sm[i] != SleepStage.wake) { + if (first < 0) first = i; + last = i; + } + } + if (first < 0) return 0; + var c = 0; + for (var i = first; i <= last; i++) { + if (sm[i] == SleepStage.wake) c++; + } + return c; + } + + test('cardio_stager: only the bout with REAL flanking sleep is bridged', () { + // 15 min sleep, then 5 × (10 min wake + 1 min sleep). The cardio rule + // table bridges ≤10 min of wake given ≥15 min of flanking sleep, so bout + // #1 legitimately bridges. Bouts #2-#5 are flanked by ONE minute of sleep + // and must survive — pre-fix each bridged bout was counted as context for + // the next, so all five collapsed and WASO went to 0. + final sm = fragmented(leadMin: 15, wakeMin: 10, reps: 5); + websterRescoreCardio(sm, 30); + expect(wakeInBody(sm), 4 * 20, + reason: 'four 10-min bouts survive; pre-fix: 0 (full cascade)'); + }); + + test('stager: only the bout with REAL flanking sleep is bridged', () { + // Same shape, sized to the classic Webster table this file uses + // (≥15 min context bridges ≤5 min of wake). + final sm = fragmented(leadMin: 15, wakeMin: 4, reps: 5); + websterRescoreAutonomic(sm, 30); + expect(wakeInBody(sm), 4 * 8, + reason: 'four 4-min bouts survive; pre-fix: 0 (full cascade)'); + }); + + test('a genuinely fragmented night keeps its WASO end-to-end', () { + // 5.1 h forced window: 60 min sleep, 6 × (8 min wake + 3 min sleep), + // 180 min sleep. Bouts #1 and #6 are legitimately bridgeable (≥15 min of + // REAL flanking sleep); bouts #2-#5 are flanked by 3 min and must stay + // wake. Measured: 1920 s WASO / 89.5% efficiency with the snapshot, + // 0 s WASO / 100.0% efficiency with the mutating context. + final accel = []; + final hr = []; + var i = 0; + void push(int secs, {required bool awake}) { + for (var k = 0; k < secs; k++, i++) { + final ph = math.sin(k * 0.5); + accel.add(awake + ? AccelSample( + (_t0 + i) * 1000.0, 0.35 * ph, 0.3, 0.9 * (1 - 0.2 * ph)) + : AccelSample((_t0 + i) * 1000.0, 0.005, 0.0, 1.0)); + hr.add(awake ? 88.0 : 50.0); + } + } + + push(60 * 60, awake: false); + for (var r = 0; r < 6; r++) { + push(8 * 60, awake: true); + push(3 * 60, awake: false); + } + push(180 * 60, awake: false); + + final s = segmentSleep(accel, hr, + forcedWindow: (onsetSec: _t0, offsetSec: _t0 + 306 * 60)); + expect(s.wasoSec!, greaterThan(25 * 60), + reason: 'pre-fix: 0 — the cascade swallowed every wake bout'); + expect(s.efficiencyPct!, lessThan(95.0), reason: 'pre-fix: 100.0'); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (5) Main-sleep selection uses the group's SPAN midpoint, not the + // gap-excluding summed duration. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — bridged-group midsleep', () { + test('a fragmented night is scored at its true circadian centre', () { + // Two candidate nights in one capture: + // A — a single 6 h block, 18:00-00:00 (span == duration, unaffected). + // B — three 2 h blocks 01:30-09:00 bridged across two 45-min gaps. + // B's detected span is 01:32→08:57 (true midsleep 05:14), but its summed + // session duration excludes the bridges, so the old + // `start + inBedSec ~/ 2` midpoint read 04:24 — 3043 s (~51 min) early. + // With the anchor below that error is exactly enough to flip which night + // is picked as the main sleep. + final accel = []; + final hr = []; + var i = 0; + final start = _midnight + 16 * 3600; // 16:00 + + void active(int secs, {double bpm = 76}) { + for (var k = 0; k < secs; k++, i++) { + final ph = math.sin(k * 0.5); + accel.add(AccelSample( + (start + i) * 1000.0, 0.3 * ph, 0.3, 0.9 * (1 - 0.2 * ph))); + hr.add(bpm + 2 * math.sin(k / 600.0)); + } + } + + void sleep(int secs, {double bpm = 50}) { + for (var k = 0; k < secs; k++, i++) { + accel.add(AccelSample((start + i) * 1000.0, 0.02, 0.02, 1.0)); + hr.add(bpm + 1.5 * math.sin(k / 1800.0)); + } + } + + active(2 * 3600); // 16:00-18:00 + sleep(6 * 3600); // 18:00-00:00 → night A + active(90 * 60, bpm: 90); // 00:00-01:30 (>60 min ⇒ A and B never bridge) + sleep(2 * 3600); // 01:30-03:30 ┐ + active(45 * 60, bpm: 90); // │ night B — three fragments bridged + sleep(2 * 3600); // 04:15-06:15 │ across two <60 min gaps + active(45 * 60, bpm: 90); // │ + sleep(2 * 3600); // 07:00-09:00 ┘ + active(2 * 3600); // 09:00-11:00 + + const anchor = 33840; // 09:24 habitual midsleep + final s = segmentSleep(accel, hr, + hrBaseline: List.filled(200, 76), + tzOffsetSec: 0, + habitualMidsleepSec: anchor); + + expect(s.present, isTrue); + final onsetSod = (s.window!.onsetMs! ~/ 1000) % 86400; + // Night B wins: onset 01:32, in-bed span 7 h 24 min (gaps INCLUDED). + expect(onsetSod, closeTo(1 * 3600 + 32 * 60, 120), + reason: 'pre-fix the early midpoint made night A (18:02) win'); + expect(s.inBedSec!, greaterThan(7 * 3600)); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (6) The habitual-midsleep anchor must convert PER TIMESTAMP (DST). + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — habitual midsleep across a DST transition', () { + test('per-timestamp resolver keeps the anchor at the true local clock time', + () { + const stdOffset = -5 * 3600; // e.g. EST + const dstOffset = -4 * 3600; // e.g. EDT + const localMidsleep = 3 * 3600; // the sleeper is dead-on 03:00 every night + const switchDay = 8; + // The transition instant: between day 7's and day 8's midsleep. + final switchTs = _midnight + switchDay * 86400; + + int offsetAt(int ts) => ts < switchTs ? stdOffset : dstOffset; + + final history = <({int startSec, int endSec, String dayKey})>[]; + for (var d = 0; d < 16; d++) { + // The UTC instant whose LOCAL clock reads 03:00 on day d. + final off = d < switchDay ? stdOffset : dstOffset; + final mid = _midnight + d * 86400 + localMidsleep - off; + history.add(( + startSec: mid - 4 * 3600, + endSec: mid + 4 * 3600, + dayKey: 'd$d', + )); + } + + final dstCorrect = + habitualMidsleepSecFromHistory(history, tzOffsetResolver: offsetAt); + // Every night really was at 03:00 local, so the anchor is 03:00 local. + expect(dstCorrect, isNotNull); + expect(dstCorrect!, closeTo(localMidsleep, 2)); + + // The old behavior — ONE frozen offset for the whole history — splits the + // days across two local clock times and lands the anchor half an hour out. + final frozen = habitualMidsleepSecFromHistory(history, + tzOffsetSeconds: dstOffset); + expect(frozen!, closeTo(localMidsleep + 1800, 60)); + expect((frozen - dstCorrect).abs(), greaterThan(1500)); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (7) Sleep-cycle minute bins must FLOOR, so pre-onset beats are dropped. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — sleep-cycle minute binning', () { + test('beats 1-59 s BEFORE onset do not land in minute 0', () { + const onset = _t0; + const offset = _t0 + 3 * 3600; // 180 min + final rrMs = []; + final rrTs = []; + + // 40 beats in the 50 s immediately BEFORE onset. `~/` truncates toward + // zero, so pre-fix these all binned to minute 0 and slipped past the + // `m < 0` guard. + for (var k = 0; k < 40; k++) { + rrTs.add((onset - 50) * 1000.0 + k * 1200.0); + rrMs.add(k.isEven ? 520.0 : 660.0); + } + // The night's ACTUAL beats start at minute 30 (nothing at all in minutes + // 0-29), so with correct binning minutes 0-19 have no data at all — not + // even after the ±10 min smoothing — and no series point exists there. + for (var m = 30; m < 170; m++) { + for (var b = 0; b < 40; b++) { + rrTs.add((onset + m * 60) * 1000.0 + b * 1400.0); + rrMs.add(1000.0 + 30.0 * math.sin(m / 14.0) + (b.isEven ? 6.0 : -6.0)); + } + } + + final r = detectSleepCycles(rrMs, rrTs, onset, offset); + expect(r.series, isNotEmpty); + final firstT = r.series.first['t'] as int; + expect(firstT, greaterThanOrEqualTo(onset + 19 * 60), + reason: 'pre-fix the pre-onset beats created a minute-0 point at ' + 't == onsetSec'); + expect(r.series.any((p) => p['t'] == onset), isFalse); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (8) van Hees: the last `sustainedMin` of a record is UNDECIDABLE, and every + // second in it is judged on ITS OWN forward window — not on one global + // trailing verdict stamped across the whole tail. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — van Hees undecidable tail', () { + const n = 3600; // 1 h at 1 Hz + const win = 300; // sustainedMin (5) × 60, the GGIR default + + /// A dead-still hour (z-angle 90°) with one 10 s reorientation burst + /// starting at [moveAt] (pass a negative value for no movement at all). + List record({required int moveAt, int moveLen = 10}) { + return [ + for (var i = 0; i < n; i++) + if (i >= moveAt && i < moveAt + moveLen) + // Flip between 90° and 30° every second — well over the 5° floor, + // and long enough to survive the 5 s rolling median. + AccelSample((_t0 + i) * 1000.0, i.isEven ? 0.87 : 0.0, 0, + i.isEven ? 0.5 : 1.0) + else + AccelSample((_t0 + i) * 1000.0, 0, 0, 1), + ]; + } + + test('a record still to its last sample does not certify the final 5 min ' + '(pre-fix: spt_sec 3600 and immobile.last == true)', () { + // Movement at 3280 sits just BEFORE the tail, so the pre-fix global + // trailing window [n-win, n) never saw it and declared all 299 tail + // seconds "no movement" — certifying sustained stillness for seconds + // whose 5 min of following data do not exist. + final m = vanHeesSleepWindow(record(moveAt: 3280)); + expect(m.present, isTrue); + final w = m.value!; + + // van Hees 2015/2018: a second qualifies as "no movement" only when the + // |Δ z-angle| stays under threshold for the FULL sustained window that + // FOLLOWS it. The last win-1 seconds of any record have no such window. + expect(w.immobile.last, isFalse, + reason: 'pre-fix: true — certified from data BEFORE it'); + expect(w.sptSec, 3301, + reason: 'pre-fix: 3600 — the rest period annexed the whole ' + 'uncertifiable tail'); + expect(w.offsetIdx, lessThanOrEqualTo(n - win + 1)); + + // ...and the tail is reported as UNDECIDABLE, not as movement. + expect(w.unresolvedTailSec, win - 1, reason: 'pre-fix: 0 (no such state)'); + expect(w.immobileUnknown.length, n); + expect(w.immobileUnknown.sublist(n - win + 1).every((u) => u), isTrue); + for (var i = 0; i < n - win + 1; i++) { + expect(w.immobileUnknown[i], isFalse, + reason: 'second $i has a full forward window — it is decided'); + } + expect(w.toJson()['unresolved_tail_sec'], win - 1); + }); + + test('a move inside the tail is resolved PER SECOND ' + '(pre-fix: one shared verdict for all 299 tail seconds)', () { + final m = vanHeesSleepWindow(record(moveAt: 3400)); + expect(m.present, isTrue); + final w = m.value!; + + // A tail second whose own forward window CONTAINS the burst: the + // sustained-inactivity rule already fails on the data in hand, whatever + // comes after the record ends ⇒ decided, and decided "moving". + expect(w.immobile[3350], isFalse); + expect(w.immobileUnknown[3350], isFalse, + reason: 'the move at 3400 decides second 3350'); + + // A tail second AFTER the burst sees nothing but stillness, but its + // window is truncated ⇒ undecidable. Pre-fix it inherited the tail-wide + // "moving" verdict from a burst that had already ended. + expect(w.immobile[3500], isFalse); + expect(w.immobileUnknown[3500], isTrue, + reason: 'pre-fix: false — a brief move ANYWHERE in the last 5 min ' + 'flipped every remaining tail second to "moving"'); + + // The two seconds must not share one answer — that is the whole bug. + expect(w.immobileUnknown[3350] == w.immobileUnknown[3500], isFalse); + expect(w.unresolvedTailSec, 193, reason: 'pre-fix: 0 (no such state)'); + }); + }); + + // ═══════════════════════════════════════════════════════════════════════════ + // (9) CPC: the Thomas 2005 HFC/LFC ratio is a MEASUREMENT or it is nothing — + // never a sentinel, never a NaN dressed up as a number. + // ═══════════════════════════════════════════════════════════════════════════ + group('honesty — cardiopulmonary coupling ratio', () { + test('a spectrum with zero coupling power abstains ' + '(pre-fix: present, hfc=0 lfc=0 cpc_ratio=0.0)', () { + // Every beat carries the SAME timestamp (a stuck clock / corrupt offload) + // and the NN deviations cancel exactly, so the Lomb-Scargle power is 0.0 + // at every frequency and both bands integrate to exactly zero. HFC/LFC is + // then 0/0 — undefined. Pre-fix it shipped as cpc_ratio 0.0, i.e. "the + // least stable sleep measurable". + final nn = List.filled(64, 1000.0); + nn[62] = 900.0; + nn[63] = 1100.0; + final ts = List.filled(64, 0.0); + + final m = cardiopulmonaryCoupling(nn, ts); + expect(m.present, isFalse, reason: 'pre-fix: present with cpc_ratio 0.0'); + expect(m.value, isNull); + expect(m.confidence, 0); + expect(m.note, contains('undefined')); + }); + + test('a NaN anywhere in the NN series abstains ' + '(pre-fix: present at confidence 0.85 with NaN band powers)', () { + // A full hour of beats → the length-driven confidence pins at its 0.85 + // cap, which is exactly what made the pre-fix output so convincing. + final nn = []; + final ts = []; + var t = 0.0; + for (var i = 0; i < 3600; i++) { + final rr = 1000 + 40 * math.sin(2 * math.pi * 0.25 * (t / 1000.0)); + t += rr; + nn.add(rr); + ts.add(t); + } + nn[1000] = double.nan; // one poisoned beat + + final m = cardiopulmonaryCoupling(nn, ts); + // `variance <= 0` does not catch NaN, so the whole spectrum came out NaN + // and was published as hfc/lfc/vlfc = NaN with cpc_ratio 0.0. + expect(m.present, isFalse, + reason: 'pre-fix: present, conf 0.85, hfc/lfc/vlfc = NaN'); + expect(m.value, isNull); + expect(m.confidence, 0); + expect(m.note, contains('non-finite')); + }); + + test('a healthy RSA record still reports a real ratio', () { + // Guard against over-abstention: the control case must survive. + final nn = []; + final ts = []; + var t = 0.0; + for (var i = 0; i < 300; i++) { + final rr = 1000 + 40 * math.sin(2 * math.pi * 0.25 * (t / 1000.0)); + t += rr; + nn.add(rr); + ts.add(t); + } + final m = cardiopulmonaryCoupling(nn, ts); + expect(m.present, isTrue); + final c = m.value!; + expect(c.cpcRatio.isFinite, isTrue); + expect(c.cpcRatio, greaterThan(0)); + expect(c.cpcRatio, isNot(999.0)); + expect(c.lfc, greaterThan(0)); + }); + }); +} diff --git a/test/onehz/wellness_test.dart b/test/onehz/wellness_test.dart index 92ae610..fcf93ea 100644 --- a/test/onehz/wellness_test.dart +++ b/test/onehz/wellness_test.dart @@ -418,4 +418,84 @@ void main() { expect(days[tempIllnessMinBaseline].need, isNull); }); }); + + // ------------------------------------------------------------------------- + // REGRESSION: degenerate (zero-dispersion) baseline columns must be dropped, + // not floored to an epsilon scale. + // ------------------------------------------------------------------------- + group('multivariateAnomaly — degenerate baseline (regression)', () { + test('an exactly-constant baseline column is DROPPED, never floored to 1e-6', + () { + // Ten baseline nights whose skin-temp z is an exactly-constant quantized + // 0.0 (MAD == 0 AND SD == 0), alongside a real HRV column, then a night + // where temp moves by a physiologically trivial 0.4. + // + // PRE-FIX the scale was `(stddev ?? 1.0).clamp(1e-6, 1e9)` => 1e-6, so + // zc = 4e5, d2 ~ 1.6e11 >> the chi-square(2) gate of 13.82 and the night + // surfaced as an illness anomaly candidate off a 0.4 change. + final feats = [ + for (var i = 0; i < 10; i++) + AnomalyFeatures(hrv: 40.0 + (i.isEven ? 1.0 : -1.0), temp: 0.0), + const AnomalyFeatures(hrv: 40.0, temp: 0.4), + const AnomalyFeatures(hrv: 40.0, temp: 0.4), + ]; + final dates = [for (var i = 0; i < feats.length; i++) 'd$i']; + final out = + multivariateAnomaly(dates, feats, minBaseline: 10, persistDays: 2); + + expect(out[10].candidate, isFalse, + reason: 'a 0.4 move on a scale-less feature is not an anomaly'); + expect(out[10].mahalanobis, isNull, + reason: 'only one feature survives => no distance is computable'); + expect(out[10].need, 'degenerate_baseline:no_dispersion'); + expect(out.where((d) => d.flagged), isEmpty); + }); + + test('a feature WITH dispersion still computes normally', () { + // Same shape, but temp now varies => both features are standardizable. + final feats = [ + for (var i = 0; i < 12; i++) + AnomalyFeatures( + hrv: 40.0 + (i.isEven ? 1.0 : -1.0), + temp: i.isEven ? 0.1 : -0.1), + ]; + final dates = [for (var i = 0; i < feats.length; i++) 'd$i']; + final out = + multivariateAnomaly(dates, feats, minBaseline: 10, persistDays: 2); + expect(out[11].mahalanobis, isNotNull); + expect(out[11].need, isNull); + expect(out[11].drivers, hasLength(2)); + }); + }); + + // ------------------------------------------------------------------------- + // REGRESSION: the glass-box driver detail must name the method ACTUALLY used. + // ------------------------------------------------------------------------- + group('readinessComposite — disclosed method matches the method used', () { + test('quantized baseline (MAD=0) discloses the mean/SD fallback', () { + // Whole-bpm RHR pinned at 55 for most of the window: MAD collapses to 0, + // robustZ abstains and the deliberate `?? z(v, base)` fallback (#26) + // produced the contribution. PRE-FIX the detail still said "robust-z". + final base = [55, 55, 55, 55, 55, 55, 55, 58]; + final m = readinessComposite([rhrInput(60, base)]); + expect(m.present, isTrue, reason: m.note); + final d = m.drivers!.single; + expect(d.detail, contains('mean+SD fallback')); + expect(d.detail, isNot(contains('robust-z'))); + }); + + test('a dispersed baseline still discloses robust-z (median+MAD)', () { + final base = [50, 52, 54, 56, 58, 60, 62]; + final m = readinessComposite([rhrInput(70, base)]); + expect(m.present, isTrue, reason: m.note); + expect(m.drivers!.single.detail, contains('robust-z (median+MAD)')); + }); + + test('a fully constant baseline (MAD=0 AND SD=0) still abstains', () { + // The fallback is NOT a licence to score against zero dispersion. + final m = readinessComposite([rhrInput(60, List.filled(8, 55.0))]); + expect(m.present, isFalse); + expect(m.toJson()['value'], '—'); + }); + }); } diff --git a/test/onehz/workout_test.dart b/test/onehz/workout_test.dart index a1e0b24..699c26e 100644 --- a/test/onehz/workout_test.dart +++ b/test/onehz/workout_test.dart @@ -332,4 +332,222 @@ void main() { expect(kcal, lessThan(10)); // resting-only over 5 min is tiny (~5–6 kcal) }); }); + + // --------------------------------------------------------------------------- + // REGRESSION: unevaluable gates must BLOCK, hidden anchors must not exist, + // off-skin samples must not become the resting-HR baseline, and the + // fabricated-anchor calorie flag must reach the output. + // --------------------------------------------------------------------------- + group('WorkoutDetector — abstain-over-fabricate (regression)', () { + /// Build a day of [restS] still/low-HR seconds, then [workS] seconds of + /// sustained motion at [workBpm], then [restS] still seconds again. + /// [offSkinS] leading seconds report hr == 0 (the off-skin sentinel) with + /// static gravity. + ({ + List hrTs, + List hrBpm, + List gTs, + List gx, + List gy, + List gz + }) day({ + required int workS, + required double workBpm, + required double restBpm, + int restS = 60, + int offSkinS = 0, + }) { + final hrTs = []; + final hrBpm = []; + final gTs = []; + final gx = []; + final gy = []; + final gz = []; + var t = 0; + void still(int n, double bpm) { + for (var i = 0; i < n; i++, t++) { + hrTs.add(t); + hrBpm.add(bpm); + gTs.add(t); + gx.add(0); + gy.add(0); + gz.add(1.0); // static gravity -> ~0 motion intensity + } + } + + still(offSkinS, 0); // OFF-SKIN: hr == 0 (types.dart HrSample convention) + still(restS, restBpm); + for (var i = 0; i < workS; i++, t++) { + hrTs.add(t); + hrBpm.add(workBpm); + gTs.add(t); + gx.add(i.isEven ? 0.5 : -0.5); // |delta| = 1.0 >> motionThreshold + gy.add(0); + gz.add(0.8); + } + still(restS, restBpm); + return (hrTs: hrTs, hrBpm: hrBpm, gTs: gTs, gx: gx, gy: gy, gz: gz); + } + + test('no HRmax anchor => the zone-2 gate is unevaluable => NO workout', () { + // age null + maxHR null + <600 HR samples => estimateHRmax returns + // (0.0, "unknown") => effMaxHR null => zonePct empty. PRE-FIX the whole + // ">=50% time in zone 2+" gate was SKIPPED, so this 6.7-minute walk at + // RHR+16 bpm was emitted as a durable workout. + final d = day(workS: 400, workBpm: 76, restBpm: 60, restS: 60); + expect(d.hrTs.length, lessThan(600), + reason: 'must stay below estimateHRmax observed-sample minimum'); + + final out = WorkoutDetector.detect( + hrTs: d.hrTs, + hrBpm: d.hrBpm, + gravTs: d.gTs, + gx: d.gx, + gy: d.gy, + gz: d.gz, + restingHR: 60, + maxHR: null, + age: null, + ); + expect(out, isEmpty); + }); + + test('detectWorkouts SAYS the gate could not be evaluated', () { + final d = day(workS: 400, workBpm: 76, restBpm: 60, restS: 60); + final m = detectWorkouts( + hrTs: d.hrTs, + hrBpm: d.hrBpm, + gravTs: d.gTs, + gx: d.gx, + gy: d.gy, + gz: d.gz, + restingHR: 60, + ); + expect(m.value, isEmpty); + expect(m.note, contains('no HRmax anchor')); + // inputs_used must reflect what was actually supplied. + expect(m.inputs_used, contains('resting_hr')); + expect(m.inputs_used, isNot(contains('max_hr'))); + expect(m.inputs_used, isNot(contains('age'))); + expect(m.inputs_used, isNot(contains('profile'))); + }); + + test('an emitted session NEVER reports strain without an HRmax anchor', () { + // PRE-FIX StrainScorer.strain(maxHR: null) silently used + // defaultMaxHR() = 220 - 30 = 190, so a session shipped a concrete strain + // alongside `hrmax: null, hrmax_source: "unknown"`. + final scenarios = >[ + for (final anchor in [null, 190]) + WorkoutDetector.detect( + hrTs: day(workS: 400, workBpm: 160, restBpm: 60).hrTs, + hrBpm: day(workS: 400, workBpm: 160, restBpm: 60).hrBpm, + gravTs: day(workS: 400, workBpm: 160, restBpm: 60).gTs, + gx: day(workS: 400, workBpm: 160, restBpm: 60).gx, + gy: day(workS: 400, workBpm: 160, restBpm: 60).gy, + gz: day(workS: 400, workBpm: 160, restBpm: 60).gz, + restingHR: 60, + maxHR: anchor, + ), + ]; + // With an anchor a real bout is emitted; without one, nothing is. + expect(scenarios[1], isNotEmpty); + expect(scenarios[0], isEmpty); + for (final list in scenarios) { + for (final s in list) { + if (s.hrmax == null) { + expect(s.strain, isNull, + reason: 'strain must abstain without a real HRmax'); + } + } + } + }); + + test('off-skin (hr==0) samples are excluded from the resting-HR percentile', + () { + // 200 s off-skin (hr == 0) then a 400 s walk at 120 bpm on a 55 bpm day. + // PRE-FIX the 10th percentile of the RAW stream was 0, so restHR = 0, + // hrFloor = 15 and %HRR for the walk was (120-0)/190 = 63% => zone 2, so + // ordinary walking cleared the >=50%-in-zone-2+ gate. On-skin only, the + // 10th percentile is 55 and %HRR is (120-55)/135 = 48% => zone 0. + final d = day( + workS: 400, workBpm: 120, restBpm: 55, restS: 400, offSkinS: 200); + final out = WorkoutDetector.detect( + hrTs: d.hrTs, + hrBpm: d.hrBpm, + gravTs: d.gTs, + gx: d.gx, + gy: d.gy, + gz: d.gz, + maxHR: 190, // isolate the resting-HR derivation + ); + expect(out, isEmpty, + reason: 'a 120 bpm walk is zone 0 against a real 55 bpm resting HR'); + }); + + test('an all-off-skin day derives no resting HR and abstains', () { + final d = day(workS: 400, workBpm: 0, restBpm: 0, restS: 60); + final out = WorkoutDetector.detect( + hrTs: d.hrTs, + hrBpm: d.hrBpm, + gravTs: d.gTs, + gx: d.gx, + gy: d.gy, + gz: d.gz, + maxHR: 190, + ); + expect(out, isEmpty); + }); + + test('the fabricated-anchor calorie flag reaches the session JSON', () { + // Calories.estimateBoutCalories returns usedDefaultAnchors precisely so a + // number built on the flat hrmax 220 / restingHr 60 fallback can be + // caveated. PRE-FIX it was computed and dropped, and ExerciseSession had + // no field or JSON key for it at all. + final flagged = Calories.estimateBoutCalories( + [0, 60, 120], + [140, 145, 150], + profile: const WorkoutUserProfile(), + hrmax: null, + restingHr: null, + ); + expect(flagged.usedDefaultAnchors, isTrue); + + final d = day(workS: 400, workBpm: 160, restBpm: 60); + final out = WorkoutDetector.detect( + hrTs: d.hrTs, + hrBpm: d.hrBpm, + gravTs: d.gTs, + gx: d.gx, + gy: d.gy, + gz: d.gz, + restingHR: 60, + maxHR: 190, + profile: const WorkoutUserProfile( + weightKg: 75, heightCm: 178, age: 30, sex: 'male'), + ); + expect(out, isNotEmpty); + final j = out.first.toJson(); + expect(j.containsKey('calories_used_default_anchors'), isTrue); + expect(j['calories_used_default_anchors'], isFalse, + reason: 'both anchors were real here'); + + // And the flag is genuinely carried, not hardcoded false. + const caveated = ExerciseSession( + start: 0, + end: 400, + avgHR: 150, + peakHR: 160, + strain: null, + durationS: 400, + zoneTimePct: {}, + avgHRRPct: null, + hrmax: null, + hrmaxSource: 'unknown', + caloriesKcal: 123.0, + caloriesKJ: 514.6, + caloriesUsedDefaultAnchors: true, + ); + expect(caveated.toJson()['calories_used_default_anchors'], isTrue); + }); + }); } From 8b1aa4e99490b8cc15397154bc548ec3fc32c44a Mon Sep 17 00:00:00 2001 From: abdulsaheel Date: Sun, 26 Jul 2026 19:02:49 +0530 Subject: [PATCH 2/2] fix: derive activity from a calibration-invariant feature, not ENMO MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Diagnosed against a real 152 MB user database: dailyStepEstimate reported 39,384 steps for a day whose true value was ~2,000. Root cause, proven by counterfactual on that raw data. ENMO is `mean(max(0, |a| - gRef))`, and gRef is auto-calibrated per day as the median of the stillest samples. On the bad day gRef came out at 0.9797; every other day it was ~1.032. The wrist rests in a different orientation during the long sleep block, and because sleep is the stillest and longest stretch it dominates that median — hourly median |a| was 0.94-0.99 asleep vs 1.03-1.06 awake, the sensor having the usual few-percent per-axis gain error. ENMO subtracts gRef from EVERY sample, so a reference 0.05 g low adds 0.05 g to every minute of the day, which is exactly the 0.05 g walking floor. Ordinary sitting cleared the gate for hours. Sweeping gRef over the identical samples: 0.97 -> 42,155 steps; 1.00 -> 13,035; 1.02 -> 0; 1.03 -> 0. There is no stable regime. The signal and the calibration error are both ~0.05 g, so SNR is ~1 by construction and no choice of threshold fixes it — the feature is wrong. Three changes, none of them a threshold tweak: FEATURE. Gravity is a constant VECTOR in the sensor frame over short windows; motion is AC. So high-pass each axis and take the magnitude of the dynamic vector, instead of estimating a scalar |g| and subtracting it. Any per-axis offset or gain error lands in the DC term and cancels exactly. Verified end to end: injecting +0.05 g of axis offset and +5% of axis gain leaves the active-minute count bit-identical. ANCHOR. Two anchors were tested and both fail. An absolute g constant (the old 0.05) is calibration-fragile, per above. A same-day relative baseline — which the old docstring falsely claimed was implemented — fails in mirror image: on a quiet day the baseline collapses, the floor collapses with it, and everything passes (it produced 16,610 on that same day). The stable anchor is a multi-day PERSONAL floor. personalDynFloor takes the pooled trailing minutes; personalDynFloorFromDailySummaries takes one persisted value per day for callers that prune their raw substrate, and uses the median across days so a single anomalous day cannot move it. With no history the estimator ABSTAINS — it does not fall back to a constant, because falling back to a constant is the bug. METRIC. Nyquist is not negotiable: gait is 1.4-2.5 Hz and 2.0 Hz (120 spm, the most common cadence) aliases exactly to DC at 1 Hz. Steps are not resolvable from this substrate; ambulatory MINUTES are, and minutes of moderate activity is the unit public activity guidance is written in. So activeMinutes becomes the primary quantity and steps are reported as a RANGE over the Tudor-Locke free-living cadence band, narrowed only when Tier A has measured this user's real cadence. The `cadence = 85 + 220*ENMO + 0.4*hrExcess` regression is deleted: it claimed to resolve cadence from a signal that provably cannot resolve cadence. Tier A (the 100 Hz AN-2554 pedometer) is untouched — it is a real, ground-truth-calibrated count and it works. Docstrings now match the implementation. The previous ones described the per-day baseline defence as present and as the reason the estimator "CANNOT inflate"; it was neither implemented nor sufficient. Outputs change for every day. Consumers must bump their algorithm version. --- lib/src/onehz/motion/enmo.dart | 124 ++++++++- lib/src/onehz/motion/motion.dart | 15 +- lib/src/onehz/motion/steps.dart | 459 ++++++++++++++++++++++-------- test/onehz/motion_test.dart | 110 ++++++++ test/onehz/steps_test.dart | 463 ++++++++++++++++++++++++++----- 5 files changed, 962 insertions(+), 209 deletions(-) diff --git a/lib/src/onehz/motion/enmo.dart b/lib/src/onehz/motion/enmo.dart index 2804e11..23c0836 100644 --- a/lib/src/onehz/motion/enmo.dart +++ b/lib/src/onehz/motion/enmo.dart @@ -1,11 +1,40 @@ -// MOTION / ACTIVITY — ENMO + MAD per-minute amplitude index. +// MOTION / ACTIVITY — per-minute amplitude indices on a 1 Hz accel stream. // -// van Hees 2013 (ENMO = Euclidean Norm Minus One) / Vähä-Ypyä 2015 (MAD, -// Mean Amplitude Deviation). The foundational 24/7 motion index on a 1 Hz -// gravity-vector accel stream. +// THREE per-minute features live here, computed in one pass: // -// ENMO_i = max(0, ‖a_i‖ − g_ref) per sample, then aggregated/min -// MAD_min = mean_i( |‖a_i‖ − mean_min(‖a‖)| ) per minute +// ENMO_i = max(0, ‖a_i‖ − g_ref) per sample, mean/minute +// (van Hees 2013, Euclidean Norm Minus One) +// MAD_min = mean_i( |‖a_i‖ − mean_min(‖a‖)| ) per minute +// (Vähä-Ypyä 2015, Mean Amplitude Deviation) +// dynAmp = mean_i( ‖a_i − rollingMean(a)‖ ) per minute ← see below +// +// WHY dynAmp EXISTS (calibration invariance). +// ENMO subtracts a SCALAR gravity reference `g_ref` from every sample. That +// reference has to be estimated from the data, and on a wrist the estimate is +// orientation-dependent: a consumer MEMS accel carries a few percent of +// per-axis gain and offset error, so the measured ‖a‖ of a motionless wrist +// differs by several 0.01 g between, say, a sleeping posture and a sitting +// posture. Auto-calibration keys off the stillest epochs, which on a 24 h day +// is the sleep block — so g_ref is biased toward the sleep orientation and is +// then subtracted from the whole waking day. The resulting bias is the SAME +// ORDER as the walking signal itself (both ~0.05 g), i.e. SNR ≈ 1 by +// construction. No threshold on ENMO can be stable against that. +// +// Vähä-Ypyä 2015 makes exactly this argument for preferring an amplitude +// measure that does not depend on knowing ‖g‖. dynAmp takes it one step +// further and removes gravity as a VECTOR rather than as a scalar norm: +// gravity is constant in the SENSOR frame over a short window, motion is AC, +// so a per-axis high-pass leaves only the dynamic component. +// +// dx_i = a_i − rollingMean(a over the last `highPassWindowS` seconds) +// dyn_i = ‖dx_i‖ , dynAmp(minute) = mean(dyn over the minute) +// +// Under a per-axis affine sensor error a' = G·a + b (G diagonal gain, b +// offset), the rolling mean maps to G·mean + b, so the OFFSET b CANCELS +// EXACTLY — every constant per-axis bias, including whatever gravity happens +// to project onto each axis in the current posture, lands in the DC term and +// is removed. A residual gain G only SCALES dyn, so any threshold expressed +// as a quantile of the user's own dynAmp distribution is invariant to it too. // // HONESTY (catalog §"what 1 Hz accel CANNOT do", Nyquist): // * 1 Hz accel gives an AMPLITUDE index only. NO steps, NO cadence, NO gait, @@ -14,7 +43,10 @@ // * Intensity bands here are RELATIVE (within-user, percentile-of-you), // NOT absolute METs — wrist 1 Hz cannot calibrate energy in MET units. // * The 1 g reference is AUTO-CALIBRATED from the data's own still epochs, -// since the sensor's zero-g offset/gain drift. +// since the sensor's zero-g offset/gain drift. [calibrateGRef] and [enmo] +// are RETAINED for the callers that legitimately want a norm-based index +// (van Hees sleep detection, Brage energy fusion) — but any decision that +// needs a STABLE absolute cut-point should use [MotionMinute.dynAmp]. import 'dart:math' as math; import '../types.dart'; @@ -27,12 +59,22 @@ class MotionMinute { final double enmo; // mean ENMO over the minute (g), ≥0 final double mad; // mean amplitude deviation over the minute (g), ≥0 final double meanMag; // mean ‖a‖ over the minute (g) — for diagnostics + + /// Mean magnitude of the GRAVITY-REMOVED (per-axis high-passed) accel vector + /// over the minute (g), ≥0. CALIBRATION-INVARIANT: any constant per-axis + /// offset — sensor bias, or the projection of gravity in the current posture + /// — cancels exactly, and a per-axis gain only rescales it. This is the + /// feature to threshold when the cut-point must be stable across days; see + /// the file header for the derivation. + final double dynAmp; + const MotionMinute( this.tsMinStartMs, this.nSamples, this.enmo, this.mad, this.meanMag, + this.dynAmp, ); Map toJson() => { 'ts_min_start_ms': tsMinStartMs, @@ -40,6 +82,7 @@ class MotionMinute { 'enmo_g': round6(enmo), 'mad_g': round6(mad), 'mean_mag_g': round6(meanMag), + 'dyn_amp_g': round6(dynAmp), }; } @@ -76,18 +119,41 @@ double calibrateGRef(List mags) { return g > 0 ? g : 1.0; } -/// Compute the ENMO + MAD per-minute motion index over a 1 Hz accel series. +/// Default high-pass window (s) that separates GRAVITY from MOTION. +/// +/// Gravity is a constant vector in the sensor frame over short windows, so +/// everything slower than ~1/[defaultGravityWindowS] Hz is treated as gravity +/// (plus per-axis sensor bias) and removed; everything faster is motion. +/// 15 s ≈ a 0.067 Hz corner: well below any voluntary movement, well above the +/// timescale on which a wrist changes posture. +const double defaultGravityWindowS = 15.0; + +/// Compute the per-minute motion indices (ENMO + MAD + dynAmp) over a 1 Hz +/// accel series, in a single pass. /// /// [samples] need not be exactly 1 Hz nor perfectly contiguous — minutes are -/// bucketed by wall-clock `tsMs`. Invalid (off-wrist) samples are dropped. -/// [gRef] overrides auto-calibration when a personal/static reference is known. -/// [minSamplesPerMinute] gates a minute as covered (default 30 = ≥50% @1 Hz). +/// bucketed by wall-clock `tsMs`, and samples are sorted by `tsMs` first so the +/// high-pass window is causal in real time. Invalid (off-wrist) samples are +/// dropped. [gRef] overrides auto-calibration when a personal/static reference +/// is known. [minSamplesPerMinute] gates a minute as covered (default 30 = +/// ≥50% @1 Hz). [gravityWindowS] is THE GRAVITY BAND: the trailing window whose +/// per-axis mean is treated as gravity + constant sensor bias and subtracted +/// from each axis before taking the magnitude (see [MotionMinute.dynAmp]). +/// +/// Edge/gap handling is honest, not padded: the trailing mean is taken over +/// whatever samples actually fall inside the window, so the first samples of a +/// series average over fewer points (the very first sample has dyn = 0 by +/// construction, since it is its own mean), and a recording gap longer than +/// [gravityWindowS] simply empties the window rather than carrying a stale +/// gravity estimate across the gap. EnmoResult enmoSeries( List samples, { double? gRef, int minSamplesPerMinute = 30, + double gravityWindowS = defaultGravityWindowS, }) { - final valid = samples.where((s) => s.valid).toList(); + final valid = samples.where((s) => s.valid).toList() + ..sort((a, b) => a.tsMs.compareTo(b.tsMs)); if (valid.isEmpty) return EnmoResult(gRef ?? 1.0, const [], 0.0); final mags = [ @@ -95,6 +161,33 @@ EnmoResult enmoSeries( ]; final ref = gRef ?? calibrateGRef(mags); + // ── per-axis high-pass → dynamic-vector magnitude (calibration-invariant) ── + // dx = a − trailingMean(a, gravityWindowS). A constant per-axis offset (and + // the constant gravity projection of the current posture) appears identically + // in a and in its mean, so it cancels EXACTLY. Window is bounded by TIME, not + // sample count, so it is sample-rate agnostic and gap-safe. + final dyn = List.filled(valid.length, 0.0); + final windowMs = gravityWindowS * 1000.0; + var lo = 0; + var sx = 0.0, sy = 0.0, sz = 0.0; + for (var i = 0; i < valid.length; i++) { + final s = valid[i]; + sx += s.x; + sy += s.y; + sz += s.z; + while (lo < i && (s.tsMs - valid[lo].tsMs) >= windowMs) { + sx -= valid[lo].x; + sy -= valid[lo].y; + sz -= valid[lo].z; + lo++; + } + final n = i - lo + 1; + final dx = s.x - sx / n; + final dy = s.y - sy / n; + final dz = s.z - sz / n; + dyn[i] = math.sqrt(dx * dx + dy * dy + dz * dz); + } + // bucket sample indices by minute final buckets = >{}; for (var i = 0; i < valid.length; i++) { @@ -122,6 +215,12 @@ EnmoResult enmoSeries( madSum += (m - meanMag).abs(); } final mad = madSum / magsMin.length; + // dynAmp: mean gravity-removed vector magnitude over the minute. + var dynSum = 0.0; + for (final i in idxs) { + dynSum += dyn[i]; + } + final dynAmp = dynSum / idxs.length; if (idxs.length >= minSamplesPerMinute) covered++; minutes.add(MotionMinute( k * 60000.0, @@ -129,6 +228,7 @@ EnmoResult enmoSeries( enmo, mad, meanMag, + dynAmp, )); } final coverage = minutes.isEmpty ? 0.0 : covered / minutes.length; diff --git a/lib/src/onehz/motion/motion.dart b/lib/src/onehz/motion/motion.dart index 1294188..d0f55bc 100644 --- a/lib/src/onehz/motion/motion.dart +++ b/lib/src/onehz/motion/motion.dart @@ -6,7 +6,10 @@ /// /// * ENMO + MAD per-minute amplitude index (van Hees 2013 / Vähä-Ypyä 2015), /// with auto-calibrated 1 g reference and RELATIVE (not MET) intensity -/// bands. → enmo.dart +/// bands — PLUS `dynAmp`, a per-axis high-passed dynamic amplitude that is +/// invariant to per-axis sensor offset (and hence to which posture the +/// gravity reference happened to be calibrated in). Threshold dynAmp, not +/// ENMO, whenever a cut-point has to hold across days. → enmo.dart /// * Static gravity-tilt orientation → sleep-position proxy (low-motion /// epochs). → orientation.dart /// * Branched HR + accel energy fusion (Brage 2004) — RELATIVE EE index, or @@ -19,9 +22,13 @@ /// 1.4–2.5 Hz), so: /// * [livePedometer] counts REAL steps on the ~100 Hz foreground accel /// (R10 / 0x2B) — adaptive-threshold peak detection (AN-2554 family). -/// * [dailyStepEstimate] gives a 24/7 ESTIMATE from the 1 Hz substrate -/// (ambulatory-minutes × cadence) — an estimate, never a count. -/// * [calibrateCadence] lets the live path personalize the 1 Hz estimate. +/// * [dailyStepEstimate] reports ACTIVE MINUTES from the 1 Hz substrate — +/// the quantity that IS resolvable there — and derives a step RANGE from +/// the free-living cadence band. Its threshold is a multi-day personal +/// reference ([personalDynFloor]); with too little history it ABSTAINS +/// rather than substituting a constant. +/// * [calibrateCadence] lets the live path narrow that range to the user's +/// own measured cadence. /// Still genuinely impossible / not faked: dynamic-orientation limb tracking, /// frequency-domain activity TYPE classification (walk vs run vs cycle). /// At 1 Hz only an AMPLITUDE index + STATIC orientation are recoverable, and diff --git a/lib/src/onehz/motion/steps.dart b/lib/src/onehz/motion/steps.dart index c1befa1..3792df5 100644 --- a/lib/src/onehz/motion/steps.dart +++ b/lib/src/onehz/motion/steps.dart @@ -17,18 +17,40 @@ // waving/typing/handling and reads 0 at rest. Directly testable: walk N // steps with the app open and compare. // -// TIER B — [dailyStepEstimate]: a 24/7 ESTIMATE from the 1 Hz substrate. We -// cannot count steps, but we CAN detect ambulatory MINUTES (ENMO in the -// walking band, optionally confirmed by HR elevation) and multiply by a -// cadence (steps/min). This is the only method the Nyquist ceiling permits: -// bout-duration × cadence (Tudor-Locke 2011: free-living walking ≈ 100–120 -// steps/min). It is an ESTIMATE, never a count — tier is ESTIMATE. +// TIER B — [dailyStepEstimate]: a 24/7 estimate from the 1 Hz substrate. We +// cannot count steps, so the PRIMARY quantity we report is the one that IS +// resolvable at 1 Hz: ACTIVE (ambulatory) MINUTES. Steps are then reported +// as a RANGE, minutes × the free-living cadence band (Tudor-Locke 2011, +// ~100–130 steps/min), never as a single fabricated-precision number. +// +// Three things make the minute detector stable, and all three matter: +// 1. The feature is [MotionMinute.dynAmp] — per-axis high-passed dynamic +// amplitude — NOT ENMO. ENMO depends on a scalar gravity reference +// whose per-day estimate moves by about the same amount as the signal +// being measured, so an absolute cut-point on ENMO is not stable +// across days. dynAmp removes gravity as a vector and is invariant to +// per-axis offset. (Vähä-Ypyä 2015 argues for calibration-robust +// amplitude measures over ENMO for exactly this reason.) +// 2. The cut-point is a MULTI-DAY PERSONAL REFERENCE ([personalDynFloor], +// a quantile of the user's POOLED trailing dynAmp minutes) — neither +// an absolute g constant (calibration-fragile) nor a same-day relative +// baseline (which collapses on a quiet day and then passes +// everything). One floor, computed over enough history to be stable, +// applied to every day. +// 3. Corroboration + duration: HR must be lifted off rest when HR is +// available, and a minute only counts inside a run of consecutive +// ambulatory minutes. +// +// With no personal reference the estimator ABSTAINS (absent Metric with a +// `need_baseline:` note). It never substitutes a constant — substituting a +// constant IS the failure mode this design exists to prevent. // // CALIBRATION — [StepCalibration] / [calibrateCadence]: Tier A is also Tier // B's teacher. When live walking data exists we measure THIS user's real -// cadence and the ENMO level it occurred at, and feed that back so the 24/7 -// estimate is personally tuned. The live path both stands alone AND makes -// the estimate "kinda accurate" per-user. +// cadence, and that measured cadence NARROWS the reported step band. It is +// used only where it is real (cadence, from a 100 Hz count); it is never +// extrapolated into a per-minute cadence regression, because 1 Hz cannot +// resolve cadence at all (gait 1.4–2.5 Hz; 2.0 Hz aliases exactly to DC). // // Pure: dart:math only. No I/O, no clock, no randomness. @@ -232,10 +254,14 @@ PedometerResult livePedometer( // ──────────────────────── CALIBRATION: live teaches the estimate ──────────── -/// A personal cadence model learned from live (100 Hz) walking, used to scale -/// the 1 Hz daily estimate. [cadenceSpm] is the user's measured walking cadence; -/// [refEnmo] is the 1 Hz ENMO level (g) observed during that same walking; [n] -/// counts the live windows folded in (more = more trusted). +/// A personal cadence model learned from live (100 Hz) walking. +/// +/// [cadenceSpm] is the user's measured walking cadence — the one quantity Tier A +/// genuinely measures and the only one Tier B consumes (to narrow its reported +/// step band; see [dailyStepEstimate]). [refEnmo] is the concurrent 1 Hz ENMO +/// level (g), retained as a diagnostic of what the norm-based index read during +/// known walking; it is NOT part of any threshold. [n] counts the live windows +/// folded in (more = more trusted). class StepCalibration { final double cadenceSpm; final double refEnmo; @@ -295,82 +321,251 @@ StepCalibration? calibrateCadence( ); } + // ───────────────────────── TIER B: 1 Hz daily estimate ────────────────────── +// +// WHY THIS LOOKS THE WAY IT DOES — the two anchors that DON'T work: +// +// (a) An ABSOLUTE g cut-point on ENMO. ENMO = max(0, ‖a‖ − gRef), and gRef is +// estimated per-day from the stillest samples. On a wrist those are the +// sleep block, whose orientation differs from the waking day; with a few +// percent of per-axis gain/offset error the still ‖a‖ can differ by +// ~0.05 g between postures. That bias is added to EVERY waking minute, +// and it is the same size as the walking signal itself. Sweeping gRef +// over one real day's raw data moved the daily total from ~42 000 steps +// to 0 with no stable plateau in between — SNR ≈ 1 by construction. +// +// (b) A SAME-DAY RELATIVE baseline (e.g. today's p20 + k·MAD). This fails in +// mirror image: on a genuinely quiet day the baseline collapses toward +// zero, the floor collapses with it, and ordinary sedentary minutes clear +// it. Same real day, same failure magnitude, opposite direction. +// +// What DOES work is a MULTI-DAY PERSONAL REFERENCE on a calibration-invariant +// feature: one floor derived from the POOLED distribution of the user's +// [MotionMinute.dynAmp] minutes across trailing history, applied to every day. +// It is stable because it is estimated from thousands of minutes, and it is +// invariant to sensor drift because dynAmp is (see enmo.dart). +// +// And when there is not enough history to estimate that floor, we ABSTAIN. + +/// Free-living walking cadence band (steps/min), Tudor-Locke 2011 (and the +/// cadence-band literature that follows it): purposeful adult ambulation in +/// free living sits around 100 steps/min, with normal walking spanning roughly +/// 100–130. We report the BAND, not a point, because 1 Hz accel cannot resolve +/// cadence at all — see [dailyStepEstimate]. +const double freeLivingCadenceLowSpm = 100.0; +const double freeLivingCadenceHighSpm = 130.0; + +/// Physiological clamp for any personally-measured cadence used to narrow the +/// band. Outside this, the "measurement" is not walking. +const double cadenceClampLowSpm = 60.0; +const double cadenceClampHighSpm = 180.0; + +/// Half-width (fraction) of the band placed around a personally MEASURED +/// cadence. Tier A measures cadence over a bout; ±10% covers the ordinary +/// within-person spread between strolling and purposeful walking. +const double personalCadenceBandFrac = 0.10; + +/// Quantile of the POOLED trailing dynAmp minutes used as the ambulatory floor. +/// +/// p90 means "the top decile of your minutes is where ambulation lives", which +/// is both a sane prior for a wrist-worn 24/7 stream (most minutes of most days +/// are sedentary or asleep) and a hard structural cap: at most ~10% of pooled +/// minutes can clear it, so no calibration excursion can produce a 1000× +/// swing in the daily total. +const double personalDynFloorQuantile = 0.90; + +/// Minimum pooled trailing minutes before a personal floor is trustworthy. +/// 2000 minutes ≈ 1.5 days of continuous wear, in practice several partial +/// days — enough that the quantile is not dominated by one posture or one day. +const int personalDynFloorMinMinutes = 2000; + +/// Minimum trailing DAYS for [personalDynFloorFromDailySummaries]. +const int personalDynFloorMinDays = 5; + +/// Multiple of the personal floor above which a minute is VIGOROUS/non-walking +/// arm motion (shaking, lifting, sport) rather than ambulation. Expressed as a +/// RATIO so it inherits the floor's calibration-invariance — an absolute g +/// ceiling would reintroduce exactly the fragility this design removes. +const double defaultVigorousCeilingRatio = 3.0; + +/// Minimum covered minutes in a day before an estimate is meaningful at all. +const int dailyStepMinCoveredMinutes = 4; + +/// Derive the PERSONAL ambulatory floor (g, in dynAmp units) from pooled +/// trailing history. +/// +/// [pooledMinuteDynAmps] is every [MotionMinute.dynAmp] the caller has for this +/// user over its trailing window — pooled ACROSS days, deliberately: a single +/// day's distribution is not a stable anchor (see the section header). Returns +/// the [quantile] of that pooled distribution, or `null` when there is not +/// enough history ([minMinutes]) or the distribution is degenerate (a +/// non-positive quantile would pass every minute). +/// +/// This package is pure — no I/O, no clock — so it cannot read history itself. +/// The caller supplies the pool; this function only decides what a floor IS. +double? personalDynFloor( + List pooledMinuteDynAmps, { + double quantile = personalDynFloorQuantile, + int minMinutes = personalDynFloorMinMinutes, +}) { + final xs = [ + for (final v in pooledMinuteDynAmps) + if (v.isFinite && v >= 0) v + ]; + if (xs.length < minMinutes) return null; + final q = percentile(xs, clamp(quantile, 0.0, 1.0) * 100.0); + if (q == null || !q.isFinite || q <= 0) return null; + return q; +} + +/// The same personal floor, derived from PER-DAY summaries instead of the raw +/// pooled minutes. +/// +/// [personalDynFloor] is the definition, but it needs every trailing minute — +/// and a caller that prunes its raw substrate (as the on-device pipeline does, +/// within days) cannot re-read them later. Persisting one high-quantile value +/// per day is ~1400× cheaper and is what a storage-bound caller can actually +/// keep, so this variant takes that: [dailyHighQuantiles] is each trailing +/// day's own [personalDynFloorQuantile] of `dynAmp`. +/// +/// It returns the MEDIAN across days rather than a quantile-of-quantiles. That +/// is the deliberate choice: the median is robust to a single anomalous day — +/// a day spent travelling, or one where the wrist sat in an unusual posture — +/// which is exactly the single-day sensitivity that makes a same-day threshold +/// unusable in the first place. Pooling the raw minutes would let one very long +/// day dominate; the median weights every day equally. +/// +/// Returns `null` below [minDays] of history, or when the result is degenerate. +double? personalDynFloorFromDailySummaries( + List dailyHighQuantiles, { + int minDays = personalDynFloorMinDays, +}) { + final xs = [ + for (final v in dailyHighQuantiles) + if (v.isFinite && v > 0) v + ]; + if (xs.length < minDays) return null; + final m = median(xs); + if (m == null || !m.isFinite || m <= 0) return null; + return m; +} + +/// The per-day value a caller should persist to feed +/// [personalDynFloorFromDailySummaries] — this day's own high quantile of +/// `dynAmp` over its covered minutes. Returns `null` when the day is too thin +/// to summarise, so the caller stores nothing rather than a fabricated level. +double? dailyDynSummary( + List motion, { + double minSamplesPerMinute = 30, + double quantile = personalDynFloorQuantile, + int minCoveredMinutes = 60, +}) { + final xs = [ + for (final m in motion) + if (m.nSamples >= minSamplesPerMinute && m.dynAmp.isFinite && m.dynAmp >= 0) + m.dynAmp + ]; + if (xs.length < minCoveredMinutes) return null; + final q = percentile(xs, clamp(quantile, 0.0, 1.0) * 100.0); + if (q == null || !q.isFinite || q <= 0) return null; + return q; +} -/// Daily step ESTIMATE from the 1 Hz substrate (never a count — see file head). +/// Daily ACTIVITY estimate from the 1 Hz substrate. +/// +/// [activeMinutes] is the PRIMARY, honest quantity: minutes spent ambulatory. +/// It is what a 1 Hz accel stream can actually support, and it is the unit +/// public activity guidance is written in (minutes of moderate activity). +/// +/// Steps are reported as the RANGE [stepsLow]–[stepsHigh] = activeMinutes × +/// the cadence band. [steps] is the midpoint, provided only so callers that +/// must render one scalar can; it carries no more information than the range +/// and should be shown with the range wherever there is room. class DailyStepEstimate { - final int steps; - final int ambulatoryMinutes; - final double cadenceUsed; // representative steps/min applied + final int activeMinutes; // primary quantity + final int stepsLow; // activeMinutes × cadenceLowSpm + final int stepsHigh; // activeMinutes × cadenceHighSpm + final int steps; // midpoint of the range (back-compat scalar) + final double cadenceLowSpm; + final double cadenceHighSpm; + final double dynFloorG; // personal floor actually applied (g) final double coverage; // fraction of the day with valid motion data - final bool calibrated; // personal cadence model was used - const DailyStepEstimate( - this.steps, - this.ambulatoryMinutes, - this.cadenceUsed, - this.coverage, - this.calibrated, - ); + final bool calibrated; // a personally MEASURED cadence narrowed the band + + const DailyStepEstimate({ + required this.activeMinutes, + required this.stepsLow, + required this.stepsHigh, + required this.steps, + required this.cadenceLowSpm, + required this.cadenceHighSpm, + required this.dynFloorG, + required this.coverage, + required this.calibrated, + }); Map toJson() => { + 'active_min': activeMinutes, + 'steps_low': stepsLow, + 'steps_high': stepsHigh, 'steps': steps, - 'ambulatory_min': ambulatoryMinutes, - 'cadence_used_spm': round6(cadenceUsed), + 'cadence_low_spm': round6(cadenceLowSpm), + 'cadence_high_spm': round6(cadenceHighSpm), + 'cadence_source': calibrated ? 'personal_measured' : 'population_band', + 'dyn_floor_g': round6(dynFloorG), 'coverage': round6(coverage), 'calibrated': calibrated, }; } -/// Default uncalibrated free-living walking cadence (Tudor-Locke 2011). -const double defaultCadenceSpm = 110.0; - -/// Default ENMO (g) we associate with that default cadence (wrist walking band). -const double defaultRefEnmoG = 0.06; - -/// Minute ENMO ceiling (g): above this is vigorous/non-walking arm motion -/// (shaking, lifting, sport) — counted toward activity elsewhere, not steps. -const double ambulatoryEnmoCeilingG = 0.40; - -/// Default FIXED movement gate (g) when uncalibrated — above resting 1 Hz noise -/// (~0.05) but below typical walking. Calibration replaces it with refEnmo·0.5. -const double defaultWalkFloorG = 0.05; - -/// Per-minute cadence regression coefficients (steps/min), literature ballpark -/// (Tudor-Locke baseline + movement & HR terms). `cadence = C0 + Cm·ENMO_g + -/// Chr·(HR−RHR)`, clamped to a physiological band. Calibration re-centres C0. -const double kStepC0 = 85.0; -const double kStepCm = 220.0; -const double kStepChr = 0.40; - -/// 1 Hz STEP ESTIMATE (the only step method that survives the Nyquist ceiling). +/// 1 Hz ACTIVE-MINUTES estimate, with steps as a derived RANGE. /// -/// We can't peak-count gait at 1 Hz (1.4–2.5 Hz aliases past 0.5 Hz), so we -/// detect WALKING minutes from the accel amplitude and multiply by a cadence — -/// the standard sub-Nyquist pedometry method. Walking detection is self-calibrated -/// + HR-corroborated + bout-gated so it CANNOT inflate (the old fixed 0.02 g floor -/// counted resting noise → ~100k/day): -/// • a minute is "ambulatory" only if its ENMO clears the day's OWN sedentary -/// baseline (p30 + 2·MAD, floored at +0.015 g) and ≤ the vigorous ceiling, -/// • AND (when HR is present) sits above the day's resting HR + [hrMarginBpm] -/// (resting = supplied RHR, else the day's 10th-percentile HR), -/// • AND belongs to a run of ≥[minBoutMin] consecutive ambulatory minutes. -/// Then steps = Σ ambulatory-minutes × cadence, where cadence is the personal -/// model ([calib]) or a default, scaled gently by intensity and clamped to a -/// physiological band. Tier is always ESTIMATE. +/// NYQUIST, stated plainly: gait is 1.4–2.5 Hz and 2.0 Hz — 120 steps/min, the +/// most common adult cadence — aliases exactly to DC on a 1 Hz stream. Steps +/// are therefore NOT resolvable here and neither is cadence. What IS resolvable +/// is whether a minute contained sustained whole-body movement. So this +/// function detects AMBULATORY MINUTES and converts them to a step RANGE using +/// a cadence band (Tudor-Locke 2011), never a per-minute cadence estimate. +/// +/// A covered minute is ambulatory when ALL of these hold: +/// • its [MotionMinute.dynAmp] is above [personalDynFloorG] and at or below +/// `personalDynFloorG × [vigorousCeilingRatio]` (above the ceiling is +/// vigorous/non-ambulatory arm motion, counted as activity elsewhere); +/// • when HR is supplied, its HR is at least `restingHr + [hrMarginBpm]` +/// ([restingHr] if given, else the day's 10th-percentile HR); +/// • it belongs to a run of at least [minBoutMin] CONSECUTIVE ambulatory +/// minutes, where consecutive means adjacent in ORIGINAL minute index — a +/// coverage gap breaks the run and cannot stitch two short stretches into +/// one qualifying bout. +/// +/// [personalDynFloorG] is REQUIRED and MAY BE NULL. Null means "not enough +/// history to know this user's movement scale", and the honest answer to that +/// is an ABSENT metric carrying a `need_baseline:have=…,need=…` note — build +/// the floor with [personalDynFloor] and pass [pooledMinutesAvailable] so the +/// note can report progress. There is deliberately NO constant fallback: a +/// constant absolute floor is precisely the failure this design removes. +/// +/// [calib] (a personally MEASURED Tier A cadence) narrows the reported band to +/// ±[personalCadenceBandFrac] around that cadence; otherwise the population +/// band is used. Tier is always ESTIMATE. /// /// IMPORTANT (no double-count): the caller must pass ONLY minutes NOT covered by /// the live 100 Hz pedometer — 100 Hz steps are real and always preferred for the /// time they cover. This function never sees those minutes. Metric dailyStepEstimate( List motion, { + required double? personalDynFloorG, List? hrPerMin, double? restingHr, StepCalibration? calib, double hrMarginBpm = 8.0, double minSamplesPerMinute = 30, int minBoutMin = 3, + double vigorousCeilingRatio = defaultVigorousCeilingRatio, + int pooledMinutesAvailable = 0, }) { - const inputs = ['enmo_per_min', 'hr_per_min', 'cadence_calibration']; + const inputs = ['dyn_amp_per_min', 'hr_per_min', 'personal_dyn_floor']; if (motion.isEmpty) { return const Metric.absent( tier: Tier.estimate, @@ -379,79 +574,89 @@ Metric dailyStepEstimate( ); } - final baseCadence = calib?.cadenceSpm ?? defaultCadenceSpm; - final refEnmo = - (calib != null && calib.refEnmo > 0) ? calib.refEnmo : defaultRefEnmoG; + // COLD START: no personal movement scale → abstain. Never substitute a + // constant floor; a constant floor is the bug this rewrite exists to fix. + final floor = personalDynFloorG; + if (floor == null || !floor.isFinite || floor <= 0) { + return Metric.absent( + tier: Tier.estimate, + inputs_used: inputs, + note: needBaselineNote( + have: pooledMinutesAvailable, + need: personalDynFloorMinMinutes, + ), + ); + } + final ceiling = floor * math.max(vigorousCeilingRatio, 1.0); // Covered minutes only — sparse minutes can't be judged. final idx = []; - final enmos = []; + final dyns = []; for (var i = 0; i < motion.length; i++) { if (motion[i].nSamples >= minSamplesPerMinute) { idx.add(i); - enmos.add(motion[i].enmo); + dyns.add(motion[i].dynAmp); } } final covered = idx.length; final coverage = covered / motion.length; - final calibrated = calib != null && calib.n >= 3; - if (covered < 4) { - return Metric( - value: DailyStepEstimate(0, 0, baseCadence, coverage, calibrated), - confidence: 0.15, + if (covered < dailyStepMinCoveredMinutes) { + return Metric.absent( tier: Tier.estimate, inputs_used: inputs, - note: 'too few covered minutes to estimate steps', + note: 'too few covered minutes to estimate activity ' + '(have=$covered, need=$dailyStepMinCoveredMinutes)', ); } - // FIXED gate + CONTINUOUS cadence (the ChatGPT-style multi-signal regression). - // We do NOT peak-count (Nyquist) and we do NOT use a per-day relative threshold - // (that self-suppressed on active days → the over/under whipsaw). A minute is - // "walking" by a STABLE gate — movement above a fixed floor AND HR lifted off - // rest — and within walking minutes the cadence scales continuously with - // movement + HR excess. Calibration tightens the floor + re-centres cadence to - // the user; uncalibrated runs a sensible ballpark (refined after a real walk). - final moveFloor = - (calibrated && refEnmo > 0) ? refEnmo * 0.5 : defaultWalkFloorG; - + // Cadence band. A personally MEASURED cadence (Tier A, 100 Hz, real counts) + // narrows the band; otherwise we use the free-living population band. We do + // NOT model per-minute cadence — 1 Hz cannot resolve it. + final measured = calib != null && + calib.n >= 3 && + calib.cadenceSpm >= cadenceClampLowSpm && + calib.cadenceSpm <= cadenceClampHighSpm + ? calib.cadenceSpm + : null; + final calibrated = measured != null; + final cadLow = calibrated + ? clamp(measured * (1 - personalCadenceBandFrac), cadenceClampLowSpm, + cadenceClampHighSpm) + : freeLivingCadenceLowSpm; + final cadHigh = calibrated + ? clamp(measured * (1 + personalCadenceBandFrac), cadenceClampLowSpm, + cadenceClampHighSpm) + : freeLivingCadenceHighSpm; + + // HR corroboration: HR must be lifted off rest for a minute to count. final useHr = hrPerMin != null && hrPerMin.length == motion.length; double restHr = restingHr ?? 0; if (useHr && restingHr == null) { final hrs = [for (final h in hrPerMin) if (h > 0) h]; if (hrs.length >= 10) restHr = percentile(hrs, 10)!; } - final hrGate = restHr + hrMarginBpm; // HR must be lifted off rest to count - - // Re-centre the cadence intercept on the personal cadence when calibrated - // (default 110 → C0 85, the literature baseline). - final c0 = calibrated ? (baseCadence - 25.0) : kStepC0; + final hrGate = restHr + hrMarginBpm; - // pass 1: does each covered minute pass the per-minute gate on its own. - // this part already existed - the bit that was missing (despite the doc - // comment above claiming it existed) is pass 2 below. + // pass 1 — per-minute gate: movement inside the ambulatory band, and (when + // HR is available) HR lifted off rest. final gateOk = List.filled(idx.length, false); - final minuteCadence = List.filled(idx.length, 0.0); for (var k = 0; k < idx.length; k++) { - final i = idx[k]; - final m = enmos[k]; - if (m <= moveFloor || m > ambulatoryEnmoCeilingG) continue; // not walking - final hr = useHr ? hrPerMin[i] : 0.0; - if (useHr && restHr > 0 && hr > 0 && hr < hrGate) continue; // HR at rest → skip - final hrExcess = (useHr && hr > 0) ? math.max(hr - restHr, 0.0) : 0.0; + final d = dyns[k]; + if (d <= floor || d > ceiling) continue; // sedentary, or vigorous non-gait + if (useHr && restHr > 0) { + final hr = hrPerMin[idx[k]]; + if (hr > 0 && hr < hrGate) continue; // HR says still at rest + } gateOk[k] = true; - minuteCadence[k] = clamp(c0 + kStepCm * m + kStepChr * hrExcess, 70.0, 170.0); } - // pass 2: only credit minutes inside a run of >= minBoutMin CONSECUTIVE - // gate-passing minutes - a scattered single minute here or there (a brief - // HR blip mid-turnover in bed, say) doesn't count on its own anymore. runs - // are broken by original minute index (idx[k]), not just position in the + // pass 2 — bout gate: only credit minutes inside a run of >= minBoutMin + // CONSECUTIVE gate-passing minutes, so a scattered single minute (a brief + // movement/HR blip mid-turnover in bed, say) never becomes phantom activity. + // Runs are broken by ORIGINAL minute index (idx[k]), not position in the // covered-minutes array, so a coverage gap can't stitch two separate // stretches into one fake long bout. - var steps = 0.0; - var ambMin = 0; - final cadences = []; + var activeMin = 0; var k = 0; while (k < idx.length) { if (!gateOk[k]) { @@ -464,33 +669,45 @@ Metric dailyStepEstimate( idx[end + 1] == idx[end] + 1) { end++; } - if (end - k + 1 >= minBoutMin) { - for (var j = k; j <= end; j++) { - steps += minuteCadence[j]; - cadences.add(minuteCadence[j]); - ambMin++; - } - } + if (end - k + 1 >= minBoutMin) activeMin += end - k + 1; k = end + 1; } - final cadenceUsed = mean(cadences) ?? baseCadence; + final stepsLow = (activeMin * cadLow).round(); + final stepsHigh = (activeMin * cadHigh).round(); + final stepsMid = ((stepsLow + stepsHigh) / 2).round(); + + // Confidence reflects (a) how much of the day we could actually judge and + // (b) whether the cadence band is this user's or the population's. It never + // reflects the step number itself — that number is a band by construction. final conf = clamp( - (calibrated ? 0.55 : 0.35) * clamp(coverage / 0.6, 0.3, 1.0), + (calibrated ? 0.45 : 0.30) * clamp(coverage / 0.6, 0.3, 1.0), 0.1, 0.7, ); return Metric( value: DailyStepEstimate( - steps.round(), ambMin, cadenceUsed, coverage, calibrated), + activeMinutes: activeMin, + stepsLow: stepsLow, + stepsHigh: stepsHigh, + steps: stepsMid, + cadenceLowSpm: cadLow, + cadenceHighSpm: cadHigh, + dynFloorG: floor, + coverage: coverage, + calibrated: calibrated, + ), confidence: conf, tier: Tier.estimate, inputs_used: inputs, note: calibrated - ? 'ESTIMATE: per-minute cadence (movement + HR) over walking minutes, ' - 'personalized — 1 Hz cannot count steps directly' - : 'ESTIMATE: per-minute cadence (movement + HR); walk with the app open ' - 'on open ground to calibrate to your stride', + ? 'ESTIMATE: active minutes from gravity-removed 1 Hz amplitude vs your ' + 'personal movement floor; steps = minutes × your measured cadence ' + '(${cadLow.round()}–${cadHigh.round()} spm) — 1 Hz cannot count steps' + : 'ESTIMATE: active minutes from gravity-removed 1 Hz amplitude vs your ' + 'personal movement floor; steps = minutes × the free-living cadence ' + 'band (${cadLow.round()}–${cadHigh.round()} spm). Walk with the app ' + 'open to measure your own cadence and narrow the range', ); } diff --git a/test/onehz/motion_test.dart b/test/onehz/motion_test.dart index 7ed6a73..564936b 100644 --- a/test/onehz/motion_test.dart +++ b/test/onehz/motion_test.dart @@ -89,6 +89,116 @@ void main() { }); }); + group('dynAmp — calibration-invariant dynamic amplitude', () { + /// 1 Hz series: gravity on z, plus an alternating ±[amp] perturbation on + /// axis [axis] during the minutes selected by [active]. + List _series( + int minutes, { + required bool Function(int) active, + double amp = 0.2, + int axis = 0, + double gz = 1.0, + }) { + final out = []; + for (var m = 0; m < minutes; m++) { + for (var s = 0; s < 60; s++) { + final i = m * 60 + s; + final p = active(m) ? (i.isEven ? amp : -amp) : 0.0; + out.add(AccelSample( + i * 1000.0, + axis == 0 ? p : 0.0, + axis == 1 ? p : 0.0, + gz + (axis == 2 ? p : 0.0), + )); + } + } + return out; + } + + test('a still wrist reads dynAmp ≈ 0 at ANY gravity reading, while ENMO ' + 'moves with the reference', () { + // Same physical stillness, two orientations the sensor reads differently + // (this ~0.05 g spread between postures is the real, measured fault). + final high = enmoSeries(_series(2, active: (_) => false, gz: 1.03), + gRef: 1.0); + final low = enmoSeries(_series(2, active: (_) => false, gz: 0.94), + gRef: 1.0); + for (final m in high.minutes) { + expect(m.dynAmp, closeTo(0.0, 1e-12)); + } + for (final m in low.minutes) { + expect(m.dynAmp, closeTo(0.0, 1e-12)); + } + // ENMO, by contrast, is entirely determined by the reference offset. + expect(high.minutes.first.enmo, closeTo(0.03, 1e-9)); + expect(low.minutes.first.enmo, closeTo(0.0, 1e-9)); + }); + + test('known alternating motion → dynAmp ≈ the perturbation amplitude', () { + final r = enmoSeries(_series(2, active: (_) => true, amp: 0.2)); + // A 15 s trailing-mean high-pass on a ±A square wave leaves A·(14/15) = + // 0.1867 (the odd-length window keeps one sample of DC leakage). + for (final m in r.minutes) { + expect(m.dynAmp, inInclusiveRange(0.17, 0.20)); + } + }); + + test('PROPERTY: a constant per-axis offset cancels EXACTLY', () { + // The exact fault: per-axis bias (and the constant gravity projection of + // whatever posture the wrist is in) is DC, so a per-axis high-pass removes + // it identically. Offsets chosen larger than the walking signal itself. + final base = _series(6, active: (m) => m.isOdd, amp: 0.3, axis: 1); + const bx = 0.05, by = -0.03, bz = 0.11; + final shifted = [ + for (final s in base) AccelSample(s.tsMs, s.x + bx, s.y + by, s.z + bz), + ]; + final a = enmoSeries(base, gRef: 1.0); + final b = enmoSeries(shifted, gRef: 1.0); + expect(b.minutes.length, a.minutes.length); + for (var i = 0; i < a.minutes.length; i++) { + expect(b.minutes[i].dynAmp, closeTo(a.minutes[i].dynAmp, 1e-9), + reason: 'dynAmp must not see a constant per-axis offset'); + } + // ENMO does see it — which is why it cannot carry an absolute cut-point. + final dEnmo = (b.minutes.first.enmo - a.minutes.first.enmo).abs(); + expect(dEnmo, greaterThan(0.05), + reason: 'ENMO shifts by the offset; that is the bug being fixed'); + }); + + test('a recording gap does not manufacture a dynAmp spike', () { + // Wrist still in posture A, an hour off-stream, then still in a totally + // different posture B. Nothing moved while we were recording. + final out = []; + for (var i = 0; i < 300; i++) { + out.add(AccelSample(i * 1000.0, 0, 0, 1.0)); + } + for (var i = 0; i < 300; i++) { + out.add(AccelSample(3600000.0 + i * 1000.0, 1.0, 0, 0)); + } + final r = enmoSeries(out); + for (final m in r.minutes) { + expect(m.dynAmp, closeTo(0.0, 1e-9), + reason: 'the gravity window must empty across a gap, not carry ' + 'a stale orientation across it'); + } + }); + + test('unsorted input is ordered before the gravity window runs', () { + final ordered = _series(2, active: (_) => true, amp: 0.25); + final shuffled = [...ordered.reversed]; + final a = enmoSeries(ordered); + final b = enmoSeries(shuffled); + for (var i = 0; i < a.minutes.length; i++) { + expect(b.minutes[i].dynAmp, closeTo(a.minutes[i].dynAmp, 1e-9)); + } + }); + + test('dynAmp is exported in toJson', () { + final r = enmoSeries(_series(1, active: (_) => true, amp: 0.2)); + expect(r.minutes.single.toJson()['dyn_amp_g'], isA()); + }); + }); + group('Static gravity-tilt → sleep position', () { test('z-up gravity (lying flat, watch face up) → supine', () { final m = staticTilt(_still(30, 0, 0, 1.0)); diff --git a/test/onehz/steps_test.dart b/test/onehz/steps_test.dart index 3f497a4..c48eea3 100644 --- a/test/onehz/steps_test.dart +++ b/test/onehz/steps_test.dart @@ -160,125 +160,444 @@ void main() { }); }); - group('Tier B — 1 Hz step estimate', () { - // Build per-minute motion rows directly (bypass enmoSeries). - List rows(List enmos) => [ - for (var i = 0; i < enmos.length; i++) - MotionMinute(i * 60000.0, 60, enmos[i], enmos[i], 1.0 + enmos[i]), + group('Tier B — personalDynFloor', () { + test('insufficient pooled history → null (never a constant)', () { + expect(personalDynFloor(List.filled(1999, 0.5)), isNull); + expect(personalDynFloor(List.filled(2000, 0.5)), isNotNull); + expect(personalDynFloor(const []), isNull); + }); + + test('a degenerate (all-zero) pool → null, not a floor of 0', () { + // A floor of 0 would pass every minute — abstaining is the honest answer. + expect(personalDynFloor(List.filled(3000, 0.0)), isNull); + }); + + test('returns the requested quantile of the pooled distribution', () { + final pool = [for (var i = 0; i < 3000; i++) i / 3000.0]; + expect(personalDynFloor(pool), closeTo(0.9, 0.01)); + expect(personalDynFloor(pool, quantile: 0.5), closeTo(0.5, 0.01)); + }); + + test('the minimum-history requirement is a named, overridable constant', () { + expect(personalDynFloorMinMinutes, 2000); + expect( + personalDynFloor(List.filled(50, 0.4), minMinutes: 10), + closeTo(0.4, 1e-9)); + }); + }); + + // The storage-bound variant. A caller that prunes its raw substrate within + // days cannot re-read trailing minutes, so it persists ONE value per day. + group('Tier B — personalDynFloorFromDailySummaries', () { + test('too few trailing days → null (never a constant)', () { + expect( + personalDynFloorFromDailySummaries( + List.filled(personalDynFloorMinDays - 1, 0.44)), + isNull); + expect( + personalDynFloorFromDailySummaries( + List.filled(personalDynFloorMinDays, 0.44)), + isNotNull); + expect(personalDynFloorFromDailySummaries(const []), isNull); + }); + + test('a single anomalous day cannot move the floor (median, not mean)', () { + // The whole point of a multi-day anchor: one day spent travelling, or + // with the wrist in an odd posture, must not drag the threshold. + final normal = [0.44, 0.43, 0.45, 0.44, 0.46, 0.43, 0.45]; + final withOutlier = [...normal, 9.0]; + final a = personalDynFloorFromDailySummaries(normal)!; + final b = personalDynFloorFromDailySummaries(withOutlier)!; + expect((a - b).abs(), lessThan(0.02), + reason: 'a 20x outlier day must barely move a median-based floor'); + }); + + test('degenerate/non-positive day summaries are dropped, not averaged in', + () { + expect(personalDynFloorFromDailySummaries(List.filled(8, 0.0)), + isNull); + final mixed = [0.44, 0.0, 0.45, -1.0, 0.43, 0.44, 0.46, 0.45]; + // Only the 6 positive days survive, which still clears the minimum. + expect(personalDynFloorFromDailySummaries(mixed), closeTo(0.445, 0.01)); + }); + }); + + group('Tier B — dailyDynSummary (what the caller persists)', () { + List mins(List dyn, {int n = 60}) => [ + for (var i = 0; i < dyn.length; i++) + MotionMinute(i * 60000.0, n, 0.055, 0.02, 1.055, dyn[i]), + ]; + + test('a day too thin to summarise yields null, not a fabricated level', () { + expect(dailyDynSummary(mins(List.filled(59, 0.4))), isNull); + expect(dailyDynSummary(mins(List.filled(60, 0.4))), isNotNull); + expect(dailyDynSummary(const []), isNull); + }); + + test('uncovered minutes do not count toward the summary', () { + // 200 rows but all sparse → below the covered-minute floor → null. + expect(dailyDynSummary(mins(List.filled(200, 0.4), n: 5)), + isNull); + }); + + test('summarises this day at the same quantile the floor is defined on', + () { + final day = [for (var i = 0; i < 1000; i++) i / 1000.0]; + expect(dailyDynSummary(mins(day)), closeTo(0.9, 0.01)); + }); + + test('round-trips: per-day summaries feed the multi-day floor', () { + // End-to-end of the persistence path: summarise each day, pool the + // summaries, get a floor — the exact sequence the caller performs. + final summaries = [ + for (var d = 0; d < 7; d++) + dailyDynSummary(mins([for (var i = 0; i < 500; i++) i / 1000.0]))! + ]; + final floor = personalDynFloorFromDailySummaries(summaries); + expect(floor, isNotNull); + expect(floor!, greaterThan(0)); + }); + }); + + group('Tier B — 1 Hz active-minutes estimate', () { + // Per-minute motion rows built directly (bypassing enmoSeries). ENMO/MAD/ + // meanMag are filled with DELIBERATELY MISLEADING values: every sedentary + // minute carries an ENMO of 0.055 g, just above the absolute 0.05 g walking + // floor the old estimator used. If anything ever re-introduces an ENMO-based + // decision path, these tests break loudly instead of shipping 39k steps. + List rows(List dyn) => [ + for (var i = 0; i < dyn.length; i++) + MotionMinute(i * 60000.0, 60, 0.055, 0.02, 1.055, dyn[i]), ]; - // A day = `sed` sedentary minutes (low ENMO) + `walk` walking minutes. - List day(int sed, int walk, - {double sedE = 0.006, double walkE = 0.06}) => - rows([...List.filled(sed, sedE), ...List.filled(walk, walkE)]); - // A learned personal walking signature — the estimate is calibration-gated. + const sedDyn = 0.02; // a sedentary minute's dynamic amplitude (g) + const walkDyn = 0.60; // an ambulatory minute's (g) + const floorG = 0.375; // the kind of value personalDynFloor yields in practice + // A day = `sed` sedentary minutes then `walk` ambulatory minutes. + List day(int sed, int walk) => rows([ + ...List.filled(sed, sedDyn), + ...List.filled(walk, walkDyn), + ]); + // A measured personal cadence from Tier A (100 Hz, real counts). const cal = StepCalibration(cadenceSpm: 110, refEnmo: 0.06, n: 10); - test('uncalibrated → a bounded ballpark (no whipsaw)', () { - // A walking block: even uncalibrated it gives a believable number (fixed - // gate, continuous cadence), not 0 and not absurd. - final m = dailyStepEstimate(day(120, 30)); // no calib - expect(m.present, isTrue); - expect(m.value!.calibrated, isFalse); - expect(m.value!.steps, inInclusiveRange(1500, 5000)); // ~30 walking min + test('COLD START: no personal floor → ABSTAIN with a need_baseline note', () { + final m = dailyStepEstimate(day(120, 30), + personalDynFloorG: null, pooledMinutesAvailable: 640); + expect(m.present, isFalse, reason: 'no constant fallback is permitted'); + expect(m.confidence, 0); + expect(m.tier, Tier.estimate); + expect(m.note, 'need_baseline:have=640,need=$personalDynFloorMinMinutes'); }); - test('calibrated sedentary day → 0 steps', () { - final m = dailyStepEstimate(rows(List.filled(600, 0.006)), calib: cal); - expect(m.value!.steps, 0); + test('a non-positive floor is treated as absent, not as "pass everything"', + () { + final m = dailyStepEstimate(day(120, 30), personalDynFloorG: 0.0); + expect(m.present, isFalse); + expect(m.note, startsWith('need_baseline:')); }); - test('a quiet day with HR at rest → 0 steps', () { - // Below-floor movement OR resting HR → nothing counts (no whipsaw to huge). - final m = dailyStepEstimate(rows(List.filled(300, 0.02)), - calib: cal, - hrPerMin: List.filled(300, 58.0), - restingHr: 58); + test('REGRESSION: a day whose sedentary minutes sit just above an ABSOLUTE ' + '0.05 g floor produces no active minutes', () { + // The measured failure shape: a calibration excursion lifted every + // sedentary minute of one day above the old absolute 0.05 g gate, and the + // day reported 39,384 steps against a true ~2,000. The personal floor is + // a multi-day reference, so a whole quiet day sitting at 0.055 g simply + // sits far below it — there is nothing for a drift to push it over. + final drifted = rows(List.filled(1400, 0.055)); + final m = dailyStepEstimate(drifted, personalDynFloorG: floorG); + expect(m.value!.activeMinutes, 0); expect(m.value!.steps, 0); + expect(m.value!.stepsHigh, 0); }); - test('a walking block over a sedentary baseline → minutes × cadence', () { - final m = dailyStepEstimate(day(120, 30), calib: cal); - expect(m.value!.ambulatoryMinutes, inInclusiveRange(28, 30)); - expect(m.value!.steps, inInclusiveRange(2800, 3600)); // ~30 × ~110 + test('REGRESSION: a quiet day cannot collapse its own threshold', () { + // The mirror-image failure of a SAME-DAY relative baseline (day p20 + + // 4·MAD): on a quiet day the baseline collapses and everything passes. + // The floor here comes from history, so a quiet day stays quiet. + final quiet = rows(List.filled(1400, sedDyn)); + final m = dailyStepEstimate(quiet, personalDynFloorG: floorG); + expect(m.value!.activeMinutes, 0); }); - test('more walking → more steps', () { - final few = dailyStepEstimate(day(120, 10), calib: cal); - final many = dailyStepEstimate(day(120, 40), calib: cal); + test('an ambulatory block over a sedentary day → active minutes × the ' + 'cadence band', () { + final m = dailyStepEstimate(day(120, 30), personalDynFloorG: floorG); + expect(m.present, isTrue); + expect(m.value!.activeMinutes, 30); + expect(m.value!.calibrated, isFalse); + // population band 100–130 spm (Tudor-Locke 2011) + expect(m.value!.stepsLow, 3000); + expect(m.value!.stepsHigh, 3900); + expect(m.value!.steps, 3450); // midpoint, the back-compat scalar + expect(m.value!.steps, + inInclusiveRange(m.value!.stepsLow, m.value!.stepsHigh)); + expect(m.value!.dynFloorG, closeTo(floorG, 1e-12)); + expect(m.tier, Tier.estimate); + }); + + test('more ambulatory minutes → more steps, and the range scales with them', + () { + final few = dailyStepEstimate(day(120, 10), personalDynFloorG: floorG); + final many = dailyStepEstimate(day(120, 40), personalDynFloorG: floorG); + expect(many.value!.activeMinutes, + greaterThan(few.value!.activeMinutes)); expect(many.value!.steps, greaterThan(few.value!.steps)); + expect(many.value!.stepsHigh - many.value!.stepsLow, + greaterThan(few.value!.stepsHigh - few.value!.stepsLow)); + }); + + test('VIGOROUS CEILING: motion far above the floor is not walking', () { + // 30 minutes of violent arm motion (shaking / a lifting set): well over + // floor × vigorousCeilingRatio, so it is activity but not ambulation. + final m = dailyStepEstimate( + rows([ + ...List.filled(120, sedDyn), + ...List.filled(30, floorG * 5), + ]), + personalDynFloorG: floorG); + expect(m.value!.activeMinutes, 0); + // …and raising the ceiling lets the same minutes through, proving the + // ceiling (not some other gate) was what rejected them. + final loose = dailyStepEstimate( + rows([ + ...List.filled(120, sedDyn), + ...List.filled(30, floorG * 5), + ]), + personalDynFloorG: floorG, + vigorousCeilingRatio: 10.0); + expect(loose.value!.activeMinutes, 30); }); - test('HR (soft veto) suppresses walking-amplitude minutes at rest HR', () { + test('HR GATE: ambulatory-amplitude minutes at resting HR do not count', () { final m = dailyStepEstimate(day(120, 30), - calib: cal, + personalDynFloorG: floorG, hrPerMin: [ ...List.filled(120, 58.0), ...List.filled(30, 58.0) ], restingHr: 58); - expect(m.value!.ambulatoryMinutes, 0, reason: 'HR at rest → not walking'); + expect(m.value!.activeMinutes, 0, reason: 'HR at rest → not walking'); }); - test('HR elevated over the walking block → it counts', () { + test('HR GATE: HR elevated over the block → it counts', () { final m = dailyStepEstimate(day(120, 30), - calib: cal, + personalDynFloorG: floorG, hrPerMin: [ ...List.filled(120, 58.0), ...List.filled(30, 95.0) ], restingHr: 58); - expect(m.value!.ambulatoryMinutes, greaterThan(20)); + expect(m.value!.activeMinutes, 30); }); - test('an isolated elevated minute does not count on its own (bout gate)', () { - final e = List.filled(60, 0.006); - e[30] = 0.20; // one elevated minute, surrounded by sedentary ones - final m = dailyStepEstimate(rows(e), calib: cal); - // a single minute alone isn't a real walk (a brief HR/movement blip - // shouldn't turn into phantom steps) - needs minBoutMin=3 in a row. - expect(m.value!.ambulatoryMinutes, 0); + test('HR GATE: resting HR falls back to the day p10 when none is supplied', + () { + final m = dailyStepEstimate(day(120, 30), + personalDynFloorG: floorG, + hrPerMin: [ + ...List.filled(120, 55.0), + ...List.filled(30, 57.0), // < p10 + 8 bpm + ]); + expect(m.value!.activeMinutes, 0); + }); + + test('BOUT GATE: an isolated elevated minute does not count on its own', () { + final d = List.filled(60, sedDyn); + d[30] = walkDyn; + final m = dailyStepEstimate(rows(d), personalDynFloorG: floorG); + expect(m.value!.activeMinutes, 0); expect(m.value!.steps, 0); }); - test('exactly at the bout boundary: 3 in a row counts, 2 does not', () { - final e2 = List.filled(60, 0.006); - e2[30] = 0.20; - e2[31] = 0.20; - final justTwo = dailyStepEstimate(rows(e2), calib: cal); - expect(justTwo.value!.ambulatoryMinutes, 0); + test('BOUT GATE: exactly at the boundary — 3 in a row counts, 2 does not', + () { + final d2 = List.filled(60, sedDyn); + d2[30] = walkDyn; + d2[31] = walkDyn; + expect(dailyStepEstimate(rows(d2), personalDynFloorG: floorG) + .value! + .activeMinutes, + 0); + + final d3 = List.filled(60, sedDyn); + d3[30] = walkDyn; + d3[31] = walkDyn; + d3[32] = walkDyn; + expect(dailyStepEstimate(rows(d3), personalDynFloorG: floorG) + .value! + .activeMinutes, + 3); + }); - final e3 = List.filled(60, 0.006); - e3[30] = 0.20; - e3[31] = 0.20; - e3[32] = 0.20; - final justThree = dailyStepEstimate(rows(e3), calib: cal); - expect(justThree.value!.ambulatoryMinutes, 3); + test('BOUT GATE: a coverage gap breaks the run rather than stitching two ' + 'short bouts together', () { + // 4 elevated minutes total but never 3 adjacent, so none of it counts. + final d = List.filled(60, sedDyn); + d[10] = walkDyn; + d[11] = walkDyn; + d[20] = walkDyn; + d[21] = walkDyn; + final m = dailyStepEstimate(rows(d), personalDynFloorG: floorG); + expect(m.value!.activeMinutes, 0); }); - test('a coverage gap breaks the run instead of stitching two short bouts together', () { - // 2 elevated minutes, a sedentary gap, then 2 more elevated minutes - - // 4 elevated minutes total but never 3 in a row, so none of it counts. - final e = List.filled(60, 0.006); - e[10] = 0.20; - e[11] = 0.20; - e[20] = 0.20; - e[21] = 0.20; - final m = dailyStepEstimate(rows(e), calib: cal); - expect(m.value!.ambulatoryMinutes, 0); + test('a MEASURED personal cadence narrows the reported range', () { + final pop = dailyStepEstimate(day(120, 30), personalDynFloorG: floorG); + final personal = dailyStepEstimate(day(120, 30), + personalDynFloorG: floorG, calib: cal); + expect(personal.value!.calibrated, isTrue); + expect(personal.value!.activeMinutes, pop.value!.activeMinutes); + expect(personal.value!.stepsHigh - personal.value!.stepsLow, + lessThan(pop.value!.stepsHigh - pop.value!.stepsLow), + reason: 'a measured cadence is a narrower band than the population'); + // ±10% around a measured 110 spm + expect(personal.value!.cadenceLowSpm, closeTo(99.0, 1e-9)); + expect(personal.value!.cadenceHighSpm, closeTo(121.0, 1e-9)); + expect(personal.confidence, greaterThan(pop.confidence)); }); - test('a higher personal cadence lifts the count', () { - final slow = dailyStepEstimate(day(120, 30), calib: cal); + test('a faster measured cadence lifts the count', () { + final slow = + dailyStepEstimate(day(120, 30), personalDynFloorG: floorG, calib: cal); final fast = dailyStepEstimate(day(120, 30), + personalDynFloorG: floorG, calib: const StepCalibration(cadenceSpm: 135, refEnmo: 0.06, n: 10)); - expect(fast.value!.cadenceUsed, greaterThan(slow.value!.cadenceUsed)); expect(fast.value!.steps, greaterThan(slow.value!.steps)); }); + test('a thin calibration (n < 3) does not claim a personal cadence', () { + final m = dailyStepEstimate(day(120, 30), + personalDynFloorG: floorG, + calib: const StepCalibration(cadenceSpm: 110, refEnmo: 0.06, n: 1)); + expect(m.value!.calibrated, isFalse); + expect(m.value!.cadenceLowSpm, freeLivingCadenceLowSpm); + }); + test('empty motion → absent ESTIMATE', () { - final m = dailyStepEstimate(const []); + final m = dailyStepEstimate(const [], personalDynFloorG: floorG); expect(m.present, isFalse); expect(m.tier, Tier.estimate); + expect(m.note, 'no motion minutes'); + }); + + test('too few covered minutes → absent, not a fabricated zero', () { + final sparse = [ + for (var i = 0; i < 3; i++) + MotionMinute(i * 60000.0, 60, 0.055, 0.02, 1.055, walkDyn), + ]; + final m = dailyStepEstimate(sparse, personalDynFloorG: floorG); + expect(m.present, isFalse); + // uncovered minutes are excluded before the count, too + final uncovered = [ + for (var i = 0; i < 100; i++) + MotionMinute(i * 60000.0, 5, 0.055, 0.02, 1.055, walkDyn), + ]; + expect(dailyStepEstimate(uncovered, personalDynFloorG: floorG).present, + isFalse); + }); + + test('toJson leads with active minutes and carries the range', () { + final j = dailyStepEstimate(day(120, 30), personalDynFloorG: floorG) + .value! + .toJson(); + expect(j['active_min'], 30); + expect(j['steps_low'], 3000); + expect(j['steps_high'], 3900); + expect(j['steps'], 3450); + expect(j['cadence_source'], 'population_band'); + }); + }); + + group('Tier B — PROPERTY: calibration invariance end to end', () { + // Build a synthetic 1 Hz day, run the whole pipeline (enmoSeries → + // personalDynFloor → dailyStepEstimate), then run it again through a + // corrupted sensor — a constant per-axis OFFSET plus a per-axis GAIN — and + // require the same number of active minutes. This is the exact fault that + // produced 39,384 steps on a real day, so it is the regression that matters + // most: the answer must not depend on the sensor's calibration state. + List synth( + int minutes, { + required double Function(int) amp, + double gx = 1.0, + double gy = 1.0, + double gz = 1.0, + double bx = 0.0, + double by = 0.0, + double bz = 0.0, + }) { + final out = []; + for (var m = 0; m < minutes; m++) { + final a = amp(m); + // Motion DIRECTION rotates between axes so per-axis gains cannot cancel + // by a trivial common factor — the invariance being tested is real. + final axis = m % 3; + for (var s = 0; s < 60; s++) { + final i = m * 60 + s; + final p = i.isEven ? a : -a; + final x = (axis == 0 ? p : 0.0); + final y = (axis == 1 ? p : 0.0) - 0.05; + final z = (axis == 2 ? p : 0.0) + 1.03; + out.add(AccelSample( + i * 1000.0, gx * x + bx, gy * y + by, gz * z + bz)); + } + } + return out; + } + + // History pool: a CONTINUOUS, low-skewed spread of minute intensities, the + // shape a real 24/7 wrist stream has — most minutes near-still, a long + // right tail. Its p90 is what personalDynFloor will pick up. + double poolAmp(int m) { + final u = (m % 200) / 200.0; + return 0.01 + 0.5 * u * u; + } + + // The day under test: near-still, except one 30-minute walk. + double dayAmp(int m) => + (m >= 600 && m < 630) ? 0.60 : 0.01 + 0.0004 * (m % 40); + + ({int active, double floor}) run({ + double gx = 1.0, + double gy = 1.0, + double gz = 1.0, + double bx = 0.0, + double by = 0.0, + double bz = 0.0, + }) { + final pool = enmoSeries(synth(2400, + amp: poolAmp, gx: gx, gy: gy, gz: gz, bx: bx, by: by, bz: bz)); + final today = enmoSeries(synth(720, + amp: dayAmp, gx: gx, gy: gy, gz: gz, bx: bx, by: by, bz: bz)); + final floor = + personalDynFloor([for (final m in pool.minutes) m.dynAmp]); + expect(floor, isNotNull, reason: '2400 pooled minutes is enough history'); + final est = + dailyStepEstimate(today.minutes, personalDynFloorG: floor); + expect(est.present, isTrue); + return (active: est.value!.activeMinutes, floor: floor!); + } + + test('the clean sensor finds the walk', () { + final r = run(); + expect(r.active, 30, reason: 'the one 30-minute walk, and only it'); + }); + + test('a constant per-axis OFFSET changes nothing', () { + expect(run(bx: 0.05, by: -0.04, bz: 0.06).active, run().active); + }); + + test('a per-axis GAIN error changes nothing', () { + expect(run(gx: 1.05, gy: 0.96, gz: 1.03).active, run().active); + }); + + test('OFFSET + GAIN together change nothing', () { + final clean = run(); + final corrupt = + run(gx: 1.05, gy: 0.96, gz: 1.03, bx: 0.05, by: -0.04, bz: 0.06); + expect(corrupt.active, clean.active); + // The floor itself DOES move — it rides the same gain the feature does, + // which is precisely why the decision does not move. + expect(corrupt.floor, isNot(closeTo(clean.floor, 1e-12))); }); });