diff --git a/ALGORITHMS.md b/ALGORITHMS.md index 593ec67a..efbd80f2 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -158,3 +158,57 @@ al. 2017 · Pimentel et al. · Brage et al. 2004 · Keytel et al. 2005 · Tudor- 2019 · Hopkins 2004 · Mahalanobis 1936 · Killick et al. 2012 · Wittmann & Roenneberg 2006 · Pietilä et al. 2018 · Malik et al. + + + +| Batch change | Method and limits | +|---|---| +| Leave long RR artifact runs as gaps | Lipponen & Tarvainen 2019; Peltola 2012; The 32-beat cap is an internal ESTIMATE rule with low confidence, not a published threshold. Longer runs split the series; each side is corrected independently. The synthetic 33-beat case keeps 200 NN beats and drops 33. | + + + +| Batch change | Method and limits | +|---|---| +| Reduce spectral error on epoch timestamps | Lomb 1976; Scargle 1982; Press & Rybicki 1989; Evaluate trigonometric arguments relative to the first sample. The mathematical spectrum is time-shift invariant; floating-point band powers and threshold decisions can change for all Lomb-Scargle callers. | + + + +| Batch change | Method and limits | +|---|---| +| Abstain when a beat has no finite time | Task Force 1996 HRV; Welch 1967 spectra; input-validity guard is an implementation rule; HRV frequency/time, nocturnal and session RMSSD, RSA and CVHR return absent metrics for non-finite beat times. This prevents endless segment walks and invalid integer conversions. Finite-input formulas are unchanged. | + + + +| Batch change | Method and limits | +|---|---| +| Refuse unidentifiable cosinor fits | Nelson et al. 1979, Methods for cosinor-rhythmometry; Welford 1962 centered co-moments; Use centered least squares and a relative determinant gate of 1e-12 of the diagonal product. This is a numerical policy. Nearly collinear phase columns abstain, and ordinary fit parameters can move in their last digits. Upstream zero-amplitude significance gating is retained. | + + + +| Batch change | Method and limits | +|---|---| +| Retain small TRIMP contributions in the total | Morton, Fitz-Clarke & Banister 1990; Neumaier 1974 compensated summation; Compensate the per-minute sum. The reserve and sex-specific weighting are unchanged; the synthetic total moves in its last digits. | + + + +| Batch change | Method and limits | +|---|---| +| Evaluate training-load decay with affine arithmetic | Banister et al. 1975; Morton et al. 1990; Neumaier 1974 prefix summation; Evaluate the same EWMA recurrence as scale times seed plus offset. This matches the arithmetic used by the later incremental implementation. Last digits change; 42/7-day defaults and bounded-prefix priming remain implementation choices. | + + + +| Batch change | Method and limits | +|---|---| +| Count repeated zone intervals before summing | Karvonen 1957 supports reserve zones only; interval credit is an internal ESTIMATE rule; Group equal duration credits by zone. Cadence validation, median capping and tail credit are retained. Grouped multiplication changes the last digits of zone seconds. | + + + +| Batch change | Method and limits | +|---|---| +| Compute circadian variance with centered moments | Witting 1990; van Someren 1999; Welford 1962 centered variance; IS and IV may move in their last digits. Batch circular profile sums and first-winner tie selection are unchanged on the rebased head. Chunked merges can reorder profile sums and move near-tie L5/M10 start epochs; that stream-specific behavior belongs with the aggregators. | + + + +| Batch change | Method and limits | +|---|---| +| Compute HRV moments without subtracting a large mean | Task Force 1996; Welford 1962 centered sums and products; Use origin-shifted moments for NN levels and five-minute bins, and centered lag co-moments for ACF. Last digits and decisions within rounding of jitter gates can change. Existing respiratory-line exemptions, coverage gates and absent-value rules remain. | diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 00000000..be0e6d9d --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,78 @@ +# Changelog + +## Split #99 output changes + +Edge must bump `kAlgoVersion` when adopting these changed batch outputs. + +### Review input and gap validity + +HRV time, nocturnal and session RMSSD abstain on non-finite NN intervals before +segmentation. Cosinor rejects non-finite or non-positive periods before computing +its design. CTL/ATL abstains on unknown/non-finite daily load; explicit observed +zero-load rest days still decay. These are input-validity policies, not new +physiological models. Capped RR runs now report `BeatClass.gap`, appended after +existing enum values, rather than claiming every omitted beat is long/short. +HRV also abstains when retained beat times regress, including when RR +accumulation disagrees with the recorded time at a capped-gap reset. Recorded +segment anchors are preserved; timestamps are never shifted to hide the conflict. +Frequency HRV, RSA and CVHR use the same regressing-clock refusal. Time HRV +also refuses supplied clock vectors whose length does not match the NN vector; +NN-only calls remain supported. Cosinor refuses non-finite design/response +moments, including overflow from a positive subnormal period or finite responses. +Other finite-input formulas, corrected NN and cap accounting are unchanged. + + +### Leave long RR artifact runs as gaps + +The 32-beat cap is an internal ESTIMATE rule with low confidence, not a published threshold. Longer runs split the series; each side is corrected independently. The synthetic 33-beat case keeps 200 NN beats and drops 33. + +Method: Lipponen & Tarvainen 2019; Peltola 2012. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Reduce spectral error on epoch timestamps + +Evaluate trigonometric arguments relative to the first sample. The mathematical spectrum is time-shift invariant; floating-point band powers and threshold decisions can change for all Lomb-Scargle callers. + +Method: Lomb 1976; Scargle 1982; Press & Rybicki 1989. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Abstain when a beat has no finite time + +HRV frequency/time, nocturnal and session RMSSD, RSA and CVHR return absent metrics for non-finite beat times. This prevents endless segment walks and invalid integer conversions. Finite-input formulas are unchanged. + +Method: Task Force 1996 HRV; Welch 1967 spectra; input-validity guard is an implementation rule. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Refuse unidentifiable cosinor fits + +Use centered least squares and a relative determinant gate of 1e-12 of the diagonal product. This is a numerical policy. Nearly collinear phase columns abstain, and ordinary fit parameters can move in their last digits. Upstream zero-amplitude significance gating is retained. + +Method: Nelson et al. 1979, Methods for cosinor-rhythmometry; Welford 1962 centered co-moments. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Retain small TRIMP contributions in the total + +Compensate the per-minute sum. The reserve and sex-specific weighting are unchanged; the synthetic total moves in its last digits. + +Method: Morton, Fitz-Clarke & Banister 1990; Neumaier 1974 compensated summation. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Evaluate training-load decay with affine arithmetic + +Evaluate the same EWMA recurrence as scale times seed plus offset. This matches the arithmetic used by the later incremental implementation. Last digits change; 42/7-day defaults and bounded-prefix priming remain implementation choices. + +Method: Banister et al. 1975; Morton et al. 1990; Neumaier 1974 prefix summation. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Count repeated zone intervals before summing + +Group equal duration credits by zone. Cadence validation, median capping and tail credit are retained. Grouped multiplication changes the last digits of zone seconds. + +Method: Karvonen 1957 supports reserve zones only; interval credit is an internal ESTIMATE rule. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Compute circadian variance with centered moments + +IS and IV may move in their last digits. Batch circular profile sums and first-winner tie selection are unchanged on the rebased head. Chunked merges can reorder profile sums and move near-tie L5/M10 start epochs; that stream-specific behavior belongs with the aggregators. + +Method: Witting 1990; van Someren 1999; Welford 1962 centered variance. Reproduce with `dart run tool/reproduce_outputs.dart`. + +### Compute HRV moments without subtracting a large mean + +Use origin-shifted moments for NN levels and five-minute bins, and centered lag co-moments for ACF. Last digits and decisions within rounding of jitter gates can change. Existing respiratory-line exemptions, coverage gates and absent-value rules remain. + +Method: Task Force 1996; Welford 1962 centered sums and products. Reproduce with `dart run tool/reproduce_outputs.dart`. + diff --git a/docs/OUTPUT_REPRODUCTION.md b/docs/OUTPUT_REPRODUCTION.md new file mode 100644 index 00000000..606088c2 --- /dev/null +++ b/docs/OUTPUT_REPRODUCTION.md @@ -0,0 +1,83 @@ +# Reproducing the #99 batch output changes + +This branch replaces the batch-output portion of #99. The aggregator branch stacks +on this branch. Edge must bump `kAlgoVersion` when adopting these changed outputs. + +Run `dart pub get` with network access on a fresh clone, then +`dart run tool/reproduce_outputs.dart`. Once dependencies are cached, +`dart pub get --offline` also works. +`--json` prints unrounded machine-readable results. The script calls batch APIs +only. Its input recipes are in `test/fixtures/output_changes.json`, a small +synthetic fixture. It contains no personal measurements. Formula parameters +expand to beat sequences in memory, including 20,000 nearly collinear cosinor +samples. The fixture's `old` and `new` values are recorded observations, not +values substituted into library results. + +Old values were measured on `3a7ca34`, new values on `f8e7d3c` with #98 merged. +The output branch reproduces the new values. Each output commit pins its own +unrounded result. Tests contain no NaN list-equality expectations. + +| Calculation | up/main `3a7ca34` | Split output branch | +|---|---|---| +| RR kept / dropped / corrected | `[199, 34, 0]` | `[200, 33, 0]` | +| Lomb-Scargle PSD at 0.137 Hz | `32627.260306094682` | `32627.260527635714` | +| Cosinor mesor / amplitude / R² / p | `[60.20591440999805, 1.9937577719142656, 0.9902626829318173, 2.900200339561918e-94]` | `[60.205914409998044, 1.9937577719142632, 0.9902626829318162, 2.9002003395774255e-94]` | +| Near-singular cosinor present | `True` | `False` | +| HRV RMSSD / SDNN / SDANN / SDNN index / ACF1 | `[3.615265234174367, 21.52663924532639, 1.067383333873884, 21.55214616158443, 0.9849552962327551]` | `[3.615265234174367, 21.52663924532641, 1.0673833338741983, 21.55214616158443, 0.984955296232754]` | +| HRV with infinite final timestamp present | `throws` | `False` | +| Banister TRIMP | `32.71197112260296` | `32.711971122602975` | +| CTL / ATL / TSB | `[98.49620732954641, 103.12182624417179, -4.625618914625377]` | `[98.49620732954642, 103.12182624417177, -4.625618914625349]` | +| Zone seconds, zones 1 through 5 | `[0.0, 120.01199999999982, 0.0, 0.0, 0.0]` | `[0.0, 120.01199999999977, 0.0, 0.0, 0.0]` | +| Circadian IS / IV / M10 / L5 / start epochs | `[1.0, 0.331267605633807, 3.505000000000011, 1.3600000000000079, 14, 0]` | `[1.0, 0.3312676056338072, 3.505000000000011, 1.3600000000000079, 14, 0]` | +| Gap-aware batch irregular SD1 / pNN70 | `[22.360679774997894, 83.33333333333333]` | `[22.360679774997894, 83.33333333333333]` | + +The synthetic RR run contains 100 normal 900 ms beats, 33 artifacts at 250 ms, +and 100 normal beats. The over-cap run becomes a gap. At 32 artifacts or fewer, +the capped and uncapped results agree. The cap is an internal ESTIMATE policy, +not a published physiological threshold. + +The cosinor phase columns in the near-singular case cover only about 20 +microhours around hour 3. Main accepts this numerically unidentifiable design; +the relative determinant gate abstains. The threshold is a numerical policy. +Nelson's significance test and nullable acrophase from #101 remain intact. + +The invalid-timestamp table uses HRV because main throws promptly. The focused +guard test also calls frequency HRV, RSA, CVHR, nocturnal RMSSD and session +RMSSD with NaN and both infinities. Running old frequency HRV/RSA with an infinite +last time is unsafe because its segment walk may never terminate. The script +does not call those old endless paths. + +Main already excludes timestamp gaps in the batch irregular-rhythm screen. +The unchanged row above is a control. #99 adds the corresponding stream policy, +previous-beat timestamp checkpointing and old-checkpoint refusal in the +aggregator branch. Diagnostics and their tests also stay there. + +Main and the rebased #99 head use the same circular profile sums and strict +first-winner comparisons for batch L5/M10. Their starts agree here and on the +near-tie controls. Centered variance changes IS/IV last digits. Chunked merges +can reorder profile sums, so near-tie starts may move in the aggregator branch. +The original #99 description's broad near-tie warning should be read in that +context. These checked fixtures and tolerances are not universal error bounds. + +To independently recompute the old column, create a detached worktree at +`3a7ca34`, copy `tool/reproduce_outputs.dart` and +`test/fixtures/output_changes.json` into the same paths there, run +`dart pub get --offline`, then run the tool with `--json`. The table printer +labels its recorded baseline separately from the currently computed result. + +Methods and limits are beside the code and in `ALGORITHMS.md`: Lipponen & +Tarvainen 2019 and Peltola 2012 for RR correction; Lomb 1976, Scargle 1982 and +Press & Rybicki 1989 for spectra; Nelson et al. 1979 for cosinor; Task Force +1996 and Welford 1962 for HRV moments; Morton, Fitz-Clarke & Banister 1990 and +Banister et al. 1975 for training load; Neumaier 1974 for compensated sums; +Witting 1990 and van Someren 1999 for circadian indices. Karvonen 1957 supports +reserve zones only. Gap credit, tail credit and the RR cap remain internal rules. + +`--corpus-json` expands twelve deterministic cases, varying RR run lengths, +spectral origins and frequencies, jitter frequencies, timestamp gaps and record +lengths. `test/fixtures/output_batch_parity.json` freezes the output branch's +unrounded public results before the aggregator conversion. Both branches run +`batch_parity_test.dart` against these fixed results. All twelve cases also +matched the rebased #99 head exactly in the split audit. This supplements the +existing law tests and frozen pre-#99 references; it does not establish parity +for every possible input. Both reproduction fixtures together are about 10 KB. diff --git a/lib/src/onehz/clinical/cosinor.dart b/lib/src/onehz/clinical/cosinor.dart index 6d4ddfab..eaa8f5b6 100644 --- a/lib/src/onehz/clinical/cosinor.dart +++ b/lib/src/onehz/clinical/cosinor.dart @@ -60,14 +60,6 @@ class CosinorFit { /// reported. const double cosinorAcrophaseAlpha = 0.05; -/// Upper tail of F(2, d): P(F > f) = (1 + 2f/d)^(−d/2), exact for 2 numerator -/// degrees of freedom. -double _f2Survival(double f, int d) { - if (f.isNaN) return 1.0; - if (f.isInfinite) return 0.0; - if (f <= 0) return 1.0; - return math.pow(1 + 2 * f / d, -d / 2).toDouble(); -} /// 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² @@ -86,74 +78,119 @@ Metric cosinor( List y, { double periodHours = 24, int minPoints = cosinorMinPoints, +}) { + if (!periodHours.isFinite || periodHours <= 0) { + return const Metric.absent( + tier: Tier.high, + inputs_used: ['signal_timeseries'], + note: 'invalid cosinor period', + ); + } + if (tHours.length != y.length) { + return const Metric.absent(tier: Tier.high, + inputs_used: ['signal_timeseries'], note: 'too few points for cosinor'); + } + final state = _CosinorState(); + final w = 2 * math.pi / periodHours; + for (var i = 0; i < y.length; i++) { + state.addCentered(math.cos(w * tHours[i]), math.sin(w * tHours[i]), y[i]); + } + return _finishCosinor(state, periodHours: periodHours, minPoints: minPoints); +} + +// Nelson et al. (1979), Methods for cosinor-rhythmometry. +// Centered least squares changes rounding. Relative 1e-12 gate is a +// numerical identifiability policy, not a physiological cutoff. +class _CosinorState { + var n = 0; + var meanC = 0.0; + var meanS = 0.0; + var meanY = 0.0; + var mCC = 0.0; + var mSS = 0.0; + var mYY = 0.0; + var mCS = 0.0; + var mCY = 0.0; + var mSY = 0.0; + + void addCentered(double c, double s, double y) { + final nextN = n + 1; + final dc = c - meanC; + final ds = s - meanS; + final dy = y - meanY; + final nextMeanC = meanC + dc / nextN; + final nextMeanS = meanS + ds / nextN; + final nextMeanY = meanY + dy / nextN; + mCC += dc * (c - nextMeanC); + mSS += ds * (s - nextMeanS); + mYY += dy * (y - nextMeanY); + mCS += dc * (s - nextMeanS); + mCY += dc * (y - nextMeanY); + mSY += ds * (y - nextMeanY); + meanC = nextMeanC; + meanS = nextMeanS; + meanY = nextMeanY; + n = nextN; + } + +} + +Metric _finishCosinor( + _CosinorState state, { + required double periodHours, + required int minPoints, }) { const inputs = ['signal_timeseries']; - final n = y.length; - if (n < math.max(4, minPoints) || tHours.length != n || periodHours <= 0) { + final n = state.n; + if (n < math.max(4, minPoints) || periodHours <= 0) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'too few points for cosinor', ); } - final w = 2 * math.pi / periodHours; - final c = [for (final t in tHours) math.cos(w * t)]; - final s = [for (final t in tHours) math.sin(w * t)]; - - // Normal equations for [M, β, γ] of y = M + β·c + γ·s. - // Build 3×3 symmetric system. - var sC = 0.0, sS = 0.0, sCC = 0.0, sSS = 0.0, sCS = 0.0; - var sY = 0.0, sYC = 0.0, sYS = 0.0; - for (var i = 0; i < n; i++) { - sC += c[i]; - sS += s[i]; - sCC += c[i] * c[i]; - sSS += s[i] * s[i]; - sCS += c[i] * s[i]; - sY += y[i]; - sYC += y[i] * c[i]; - sYS += y[i] * s[i]; + // A finite period or response can still overflow the design arithmetic. + // Reject invalid moments before NaN can bypass the identifiability gate. + if (![state.meanC, state.meanS, state.meanY, state.mCC, state.mSS, + state.mYY, state.mCS, state.mCY, state.mSY].every((v) => v.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite cosinor design', + ); } - // Matrix A = [[n, sC, sS],[sC, sCC, sCS],[sS, sCS, sSS]], b=[sY,sYC,sYS]. - final a = [ - [n.toDouble(), sC, sS], - [sC, sCC, sCS], - [sS, sCS, sSS], - ]; - final b = [sY, sYC, sYS]; - final sol = _solve3(a, b); - if (sol == null) { + final determinant = state.mCC * state.mSS - state.mCS * state.mCS; + // The determinant is relative to the diagonal co-moments. A design needs + // at least 1e-12 of that scale after removing its means to identify both + // cosine coefficients; below that, numerical noise decides the fit. + if (determinant <= 1e-12 * state.mCC * state.mSS) { return const Metric.absent( tier: Tier.high, inputs_used: inputs, note: 'singular cosinor design (insufficient phase coverage)', ); } - final mesor = sol[0]; - final beta = sol[1]; - final gamma = sol[2]; + final beta = + (state.mCY * state.mSS - state.mSY * state.mCS) / determinant; + final gamma = + (state.mSY * state.mCC - state.mCY * state.mCS) / determinant; + final mesor = state.meanY - beta * state.meanC - gamma * state.meanS; final amplitude = math.sqrt(beta * beta + gamma * gamma); - final acro = math.atan2(-gamma, beta); // peak phase - // Convert acrophase to clock-hours within the period (0..period). - var phaseHours = -acro / w; // time of peak relative to origin + final acro = math.atan2(-gamma, beta); + final w = 2 * math.pi / periodHours; + var phaseHours = -acro / w; phaseHours = phaseHours % periodHours; if (phaseHours < 0) phaseHours += periodHours; - // R². - final yMean = sY / n; - var ssTot = 0.0, ssRes = 0.0; - for (var i = 0; i < n; i++) { - final fit = mesor + beta * c[i] + gamma * s[i]; - ssTot += (y[i] - yMean) * (y[i] - yMean); - ssRes += (y[i] - fit) * (y[i] - fit); - } + final ssTot = state.mYY; + final ssRes = state.mYY - + 2 * beta * state.mCY - + 2 * gamma * state.mSY + + beta * beta * state.mCC + + gamma * gamma * state.mSS + + 2 * beta * gamma * state.mCS; final double r2 = ssTot == 0 ? 0.0 : (1 - ssRes / ssTot).clamp(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 double r2Adj = (1 - (1 - r2) * (n - 1) / (n - 3)).clamp(0, 1); - // Zero-amplitude F-test (Nelson 1979): the acrophase is only defined once // the test rejects. A flat or rhythm-free series still yields some // atan2(−γ, β) from residual noise or float rounding. @@ -169,7 +206,6 @@ Metric cosinor( } final phaseDefined = ampP < cosinorAcrophaseAlpha; - final conf = r2Adj.clamp(0.1, 0.95); return Metric( value: CosinorFit( mesor: mesor, @@ -181,7 +217,7 @@ Metric cosinor( r2Adj: r2Adj, amplitudeP: ampP, ), - confidence: conf, + confidence: r2Adj.clamp(0.1, 0.95), tier: Tier.high, inputs_used: inputs, note: phaseDefined @@ -192,34 +228,12 @@ Metric cosinor( ); } -/// Solve a 3×3 linear system by Gaussian elimination with partial pivoting. -/// Returns null if singular. -List? _solve3(List> a, List b) { - final m = [ - [a[0][0], a[0][1], a[0][2], b[0]], - [a[1][0], a[1][1], a[1][2], b[1]], - [a[2][0], a[2][1], a[2][2], b[2]], - ]; - for (var col = 0; col < 3; col++) { - var piv = col; - for (var r = col + 1; r < 3; r++) { - if (m[r][col].abs() > m[piv][col].abs()) piv = r; - } - if (m[piv][col].abs() < 1e-12) return null; - final tmp = m[col]; - m[col] = m[piv]; - m[piv] = tmp; - final d = m[col][col]; - for (var j = col; j < 4; j++) { - m[col][j] /= d; - } - for (var r = 0; r < 3; r++) { - if (r == col) continue; - final f = m[r][col]; - for (var j = col; j < 4; j++) { - m[r][j] -= f * m[col][j]; - } - } - } - return [m[0][3], m[1][3], m[2][3]]; + +/// Upper tail of F(2, d): P(F > f) = (1 + 2f/d)^(−d/2), exact for 2 numerator +/// degrees of freedom. +double _f2Survival(double f, int d) { + if (f.isNaN) return 1.0; + if (f.isInfinite) return 0.0; + if (f <= 0) return 1.0; + return math.pow(1 + 2 * f / d, -d / 2).toDouble(); } diff --git a/lib/src/onehz/clinical/hrv_freq.dart b/lib/src/onehz/clinical/hrv_freq.dart index d7592486..1cc57995 100644 --- a/lib/src/onehz/clinical/hrv_freq.dart +++ b/lib/src/onehz/clinical/hrv_freq.dart @@ -84,6 +84,20 @@ Metric hrvFreq( // Times in seconds for Hz output; tachogram value = NN (ms). final tSec = [for (final t in nnTimesMs) t / 1000.0]; + if (tSec.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { return const Metric.absent( diff --git a/lib/src/onehz/clinical/hrv_time.dart b/lib/src/onehz/clinical/hrv_time.dart index d2b35104..e34f0477 100644 --- a/lib/src/onehz/clinical/hrv_time.dart +++ b/lib/src/onehz/clinical/hrv_time.dart @@ -17,6 +17,7 @@ import 'dart:math' as math; import '../types.dart'; import '../util.dart'; +import 'rmssd_gate.dart'; /// Lag-1 ACF floor for the NN successive-difference series, below which RMSSD /// and pNN50 are REFUSED (SDNN / SDANN survive it and are the honest lead). @@ -121,82 +122,13 @@ const double kNnDiffNoiseShareCeiling = 0.7; /// clear of the floor, or diffs on a beat-time grid coarse against their size /// (see [_onCoarseLattice]). double? nnDiffNoiseShare(List> diffRuns, {int minSegments = 20}) { - const n = 64, kLo = 10, kHi = n ~/ 2; // kLo/n ≈ 0.15 cycles/beat - final w = [ - for (var i = 0; i < n; i++) 0.5 - 0.5 * math.cos(2 * math.pi * i / (n - 1)) - ]; - var w2 = 0.0; - for (final x in w) { - w2 += x * x; - } - final cosT = [ - for (var k = kLo; k <= kHi; k++) - [for (var i = 0; i < n; i++) math.cos(2 * math.pi * k * i / n)] - ]; - final sinT = [ - for (var k = kLo; k <= kHi; k++) - [for (var i = 0; i < n; i++) math.sin(2 * math.pi * k * i / n)] - ]; - final psd = List.filled(kHi + 1, 0.0); - var segs = 0, nd = 0; - var ssd = 0.0, wsd = 0.0, wsum = 0.0, pSum = 0.0, p2Sum = 0.0; + final sums = NoiseShareSums(); for (final r in diffRuns) { - for (final d in r) { - ssd += d * d; - nd++; - } - // Integrate the run back to RR levels (up to a constant the mean removes). - final x = List.filled(r.length + 1, 0.0); - for (var i = 0; i < r.length; i++) { - x[i + 1] = x[i] + r[i]; - } - final cov = List.filled(x.length, 0.0); // Σ w² per sample - for (var s = 0; s + n <= x.length; s += n ~/ 2) { - var m = 0.0; - for (var i = 0; i < n; i++) { - m += x[s + i]; - cov[s + i] += w[i] * w[i]; - } - m /= n; - final v = [for (var i = 0; i < n; i++) (x[s + i] - m) * w[i]]; - var p = 0.0; - for (var k = kLo; k <= kHi; k++) { - final c = cosT[k - kLo], sn = sinT[k - kLo]; - var re = 0.0, im = 0.0; - for (var i = 0; i < n; i++) { - re += v[i] * c[i]; - im -= v[i] * sn[i]; - } - psd[k] += re * re + im * im; - p += re * re + im * im; - } - pSum += p; - p2Sum += p * p; - segs++; - } - for (var i = 0; i < r.length; i++) { - final c = (cov[i] + cov[i + 1]) / 2; - wsd += c * r[i] * r[i]; - wsum += c; - } + sums.addRun(r); } - if (segs < minSegments || ssd == 0 || wsd == 0) return null; - if (_onCoarseLattice(diffRuns, ssd / nd)) return null; - // One loud stretch dominates the average: count segments by power, not number. - if (pSum * pSum / p2Sum < 0.75 * minSegments) return null; - final band = psd.sublist(kLo); - var peak = 0; - for (var k = 1; k < band.length; k++) { - if (band[k] > band[peak]) peak = k; - } - if (peak + kLo >= kHi - 1) return null; // ~0.477–0.5 band, see above - if (band.last > 0.2 * band[peak]) return null; // hump on Nyquist, see above - final med = median(band)!; - if (band[peak] < 4.5 * med) return null; // no line, no exemption - final floor = med / segs / w2; - final vetted = wsd / wsum, all = ssd / nd; - // 1 − (structured share of RMSSD²); equals 2σ²/mean(d²) when stationary. - return 1 - (vetted - 2 * floor) / math.max(vetted, all); + return noiseShareFrom(sums, + lattice: () => _onCoarseLattice(diffRuns, sums.ssd / sums.nd), + minSegments: minSegments); } /// True when every successive difference is a whole multiple of one step q and @@ -208,22 +140,8 @@ double? nnDiffNoiseShare(List> diffRuns, {int minSegments = 20}) { /// sets RMSSD. // ponytail: exact lattice only; a grid re-rounded to whole ms (7.8 ms → 7/8) // breaks the common step and is not caught. -bool _onCoarseLattice(List> diffRuns, double msd) { - var q = double.infinity; - for (final r in diffRuns) { - for (final d in r) { - if (d.abs() > 1e-6 && d.abs() < q) q = d.abs(); - } - } - if (q.isInfinite || msd >= 2 * q * q) return false; - for (final r in diffRuns) { - for (final d in r) { - final k = d.abs() / q; - if ((k - k.roundToDouble()).abs() > 1e-3) return false; - } - } - return true; -} +bool _onCoarseLattice(List> diffRuns, double msd) => + coarseLattice([for (final r in diffRuns) for (final d in r) d.abs()], msd); /// The one RMSSD jitter verdict every path shares: ACF1 below the floor AND /// the spectrum does not show a respiratory line on a low white floor. @@ -256,48 +174,6 @@ List _clearedWindows(double? nightAcf1, List rmssds, ]; } -final _hann128 = [ - for (var i = 0; i < 128; i++) 0.5 - 0.5 * math.cos(2 * math.pi * i / 127) -]; -final _hann128Sq = _hann128.fold(0.0, (a, w) => a + w * w); - -/// Welch power (128-beat Hann, 50 % overlap, averaged over segments) of one -/// window's RR, rebuilt from its difference runs, at cycles-per-beat -/// frequencies [f]. Null when no run holds a full segment. -List? _windowPsd(List> diffRuns, List f) { - const n = 128; - final p = List.filled(f.length, 0.0); - var segs = 0; - for (final r in diffRuns) { - final x = List.filled(r.length + 1, 0.0); - for (var i = 0; i < r.length; i++) { - x[i + 1] = x[i] + r[i]; - } - for (var s = 0; s + n <= x.length; s += n ~/ 2) { - var m = 0.0; - for (var i = 0; i < n; i++) { - m += x[s + i]; - } - m /= n; - for (var j = 0; j < f.length; j++) { - final c = math.cos(2 * math.pi * f[j]), sn = math.sin(2 * math.pi * f[j]); - var cr = 1.0, ci = 0.0, re = 0.0, im = 0.0; - for (var i = 0; i < n; i++) { - final v = (x[s + i] - m) * _hann128[i]; - re += v * cr; - im += v * ci; - final t = cr * c - ci * sn; - ci = cr * sn + ci * c; - cr = t; - } - p[j] += re * re + im * im; - } - segs++; - } - } - return segs == 0 ? null : [for (final v in p) v / segs]; -} - /// Last resort for a night the jitter gate refused: the windows whose /// breathing line holds still in Hz while the heart rate moves. /// @@ -324,115 +200,29 @@ List? _steadyBreathingWindows( nd++; } } - if (nd == 0 || _onCoarseLattice(all, ssd / nd)) return null; - final ref = median([for (final m in meanRrMs) if (m > 0) m]); - if (ref == null) return null; - // 0.15–0.5 cycles/beat at the night's median beat. - final fc = [for (var j = 38; j <= 128; j++) j / 256]; - final byBeat = List.filled(fc.length, 0.0); - final byHz = List.filled(fc.length, 0.0); - final cover = List.filled(fc.length, 0); - final peakHz = {}; // window -> its own peak, cycles per ms - final psd = >{}; // window -> its cycles-per-beat PSD - for (var w = 0; w < winRuns.length; w++) { - if (meanRrMs[w] <= 0) continue; - final pb = _windowPsd(winRuns[w], fc); - if (pb == null) continue; - final norm = median(pb)!; - if (norm <= 0) continue; - final k = meanRrMs[w] / ref; - final fh = [for (final f in fc) f * k]; - final ph = _windowPsd(winRuns[w], fh)!; - var best = -1; - for (var j = 0; j < fc.length; j++) { - byBeat[j] += pb[j] / norm; - if (fh[j] > 0.5) continue; - byHz[j] += ph[j] / norm; - cover[j]++; - if (best < 0 || ph[j] > ph[best]) best = j; - } - if (best >= 0) peakHz[w] = fc[best] / ref; - psd[w] = pb; - } - final nw = peakHz.length; - if (nw < 12) return null; - final bins = [for (var j = 0; j < fc.length; j++) if (cover[j] == nw) j]; - if (bins.length < 20) return null; - final hz = [for (final j in bins) byHz[j]]; - final bt = [for (final j in bins) byBeat[j]]; - var pk = 0; - for (var i = 1; i < hz.length; i++) { - if (hz[i] > hz[pk]) pk = i; - } - if (pk == 0 || pk == hz.length - 1) return null; - final f0 = fc[bins[pk]] / ref; - final brpm = f0 * 60000; - if (brpm < 8 || brpm > 30) return null; - // Same bins, same normalisation: the line must stand taller aligned in Hz - // than the beat pooling stands anywhere within ±12 bins of it. - var near = 0.0; - for (var i = math.max(0, pk - 12); i <= math.min(bt.length - 1, pk + 12); i++) { - near = math.max(near, bt[i]); - } - // The pooled noise ripple shrinks as 1/√windows, so the bar does too. - final promHz = hz[pk] / median(hz)!; - if (promHz < math.max(2.5, 1.5 + 7 / math.sqrt(nw))) return null; - if (promHz < 1.1 * near / median(bt)!) return null; - // Alternation lives at Nyquist on the beat axis, outside these bins. - if (byBeat.last / median(bt)! >= 0.5 * promHz) return null; - final keep = [ - for (final e in peakHz.entries) - if (math.log(e.value / f0).abs() <= 0.15) e.key - ]; - if (keep.length < 6) return null; - // The line proves breathing is there, not that it carries RMSSD: beat-time - // jitter loud enough to fail the gate would still be most of the number. - // Same ceiling as [nnDiffNoiseShare], over the kept windows. Jitter σ² has - // an RR spectrum σ²·(2 − 2cos ω), flat once divided by that shape, and its - // differences carry 6σ²; physiology only adds power on top. Breathing - // wanders across the night and its line spreads with it, so each window - // drops the bins of every rate the kept windows peaked at (plus a Hann - // main lobe), mapped onto its own beats. σ² is the floor of what is left, - // pooled over windows and smoothed over four resolution cells. - var lo = double.infinity, hi = 0.0; - for (final w in keep) { - lo = math.min(lo, peakHz[w]!); - hi = math.max(hi, peakHz[w]!); - } - var sq = 0.0, nAll = 0; - final acc = List.filled(fc.length, 0.0); - final wt = List.filled(fc.length, 0.0); - for (final w in keep) { - var n = 0; - for (final r in winRuns[w]) { - for (final d in r) { - sq += d * d; - n++; + return steadyBreathingFrom( + windows: winRuns.length, + meanRrMs: meanRrMs, + psdAt: (w, f) => windowPsd(winRuns[w], f), + diffs: (w) { + var n = 0; + for (final r in winRuns[w]) { + n += r.length; } - } - nAll += n; - final a = lo * meanRrMs[w] - 3 / 128, b = hi * meanRrMs[w] + 3 / 128; - for (var j = 0; j < fc.length; j++) { - if (fc[j] > a && fc[j] < b) continue; - acc[j] += n * psd[w]![j] / (2 - 2 * math.cos(2 * math.pi * fc[j])); - wt[j] += n; - } - } - var floor = double.infinity; - for (var j = 8; j < fc.length - 8; j++) { - var s = 0.0, m = 0; - for (var i = j - 8; i <= j + 8; i++) { - if (wt[i] < nAll / 2) break; - s += acc[i] / wt[i]; - m++; - } - if (m == 17) floor = math.min(floor, s / 17); - } - if (floor == double.infinity) return null; - // The lowest stretch of a noisy curve reads below its mean; 1.1 puts the - // floor back at these pooled sizes, so jitter alone is not undercounted. - final noise = 1.1 * 6 * floor / _hann128Sq * nAll; - return noise < kNnDiffNoiseShareCeiling * sq ? keep : null; + return n; + }, + addSquares: (w, acc) { + for (final r in winRuns[w]) { + for (final d in r) { + acc += d * d; + } + } + return acc; + }, + allSsd: ssd, + allNd: nd, + lattice: (msd) => _onCoarseLattice(all, msd), + ); } /// The window RMSSDs a windowed headline publishes, or null for none: the @@ -624,6 +414,32 @@ Metric hrvTime( RrCoverage? coverage, }) { const inputs = ['rr_cleaned']; + // An unknown interval cannot support a finite HRV estimate. Refuse the + // input before segmentation or moments rather than discarding evidence. + if (nnTimesMs != null && nnTimesMs.length != nnMs.length) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'beat time length mismatch', + ); + } + // A reset may expose disagreement between accumulated RR and recorded + // sensor time. Do not turn overlapping segments into successive beats. + if (nnTimesMs != null && nnTimesMs.length == nnMs.length && + hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } + if (nnMs.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite NN interval', + ); + } if (nnMs.length < 2) { return const Metric.absent( tier: Tier.high, @@ -632,6 +448,15 @@ Metric hrvTime( ); } + if (nnTimesMs != null && + nnTimesMs.length == nnMs.length && + nnTimesMs.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } // RMSSD / pNN50: root mean square of SUCCESSIVE differences — successive in // TIME, not merely adjacent in the compacted list. correctRr drops multi-beat // artifact runs while advancing its clock across them, so nn[i-1] and nn[i] @@ -672,7 +497,29 @@ Metric hrvTime( // are far less contaminated (jitter share 1–27 % against RMSSD's 11–100 % on // the audit corpus) — they keep publishing, which is what the header has // always advised. - final acf1 = nnDiffAcf1(runs); + // Welford (1962) centered sums/products, Task Force (1996) HRV. + // Accumulate difference moments and lagged co-moments in sample order. + var dn = 0, ln = 0; + var dm = 0.0, d2 = 0.0, ma = 0.0, mb = 0.0, co = 0.0; + for (final r in runs) { + for (var i = 0; i < r.length; i++) { + final d = r[i]; + dn++; + final delta = d - dm; + dm += delta / dn; + d2 += delta * (d - dm); + if (i > 0) { + ln++; + final da = r[i - 1] - ma; + final db = d - mb; + ma += da / ln; + mb += db / ln; + co += da * (d - mb); + } + } + } + final double? acf1 = dn >= 30 && d2 > 0 + ? (co + ln * (ma - dm) * (mb - dm)) / d2 : null; final jittery = _jitterRefused(acf1, runs); // More beat-time than elapsed: the stream itself is wrong, so neither the // differences nor any breathing line in them describe one heart. @@ -686,15 +533,16 @@ Metric hrvTime( rmssd = _breathingRmssd(nnMs, nnTimesMs); } final byBreathing = jittery && rmssd != null; - final sdnn = stddev(nnMs); + final sdnn = _centeredMoments(nnMs).sd; double? sdann, sdnnIndex; if (gapAware) { final seg = _fiveMinSegments(nnMs, nnTimesMs); if (seg.length >= 2) { - final means = [for (final s in seg) mean(s)!]; - sdann = stddev(means); - final sds = [for (final s in seg) stddev(s)].whereType().toList(); + final moments = [for (final s in seg) _centeredMoments(s)]; + final means = [for (final m in moments) m.mean!]; + sdann = _centeredMoments(means).sd; + final sds = [for (final m in moments) m.sd].whereType().toList(); sdnnIndex = sds.isEmpty ? null : mean(sds); } } @@ -804,6 +652,24 @@ Metric nocturnalRmssd( RrCoverage? coverage, }) { const inputs = ['rr_cleaned', 'beat_times']; + // An unknown interval cannot support a finite HRV estimate. Refuse the + // input before segmentation or moments rather than discarding evidence. + // A reset may expose disagreement between accumulated RR and recorded + // sensor time. Do not turn overlapping segments into successive beats. + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } + if (nnMs.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite NN interval', + ); + } if (coverage != null && coverage.overCounted) { return Metric.absent( tier: Tier.high, @@ -818,6 +684,13 @@ Metric nocturnalRmssd( note: 'too few NN intervals for windowed nocturnal RMSSD', ); } + if (nnTimesMs.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } final t0 = nnTimesMs.first; // Bucket beat INDICES by window index, so each window keeps its beat times and // the successive differences below can skip the ones that straddle a dropped @@ -884,7 +757,29 @@ Metric nocturnalRmssd( meanRr.add(mean([for (final i in seg) nnMs[i]])!); rmssds.add(math.sqrt(ssd / nd)); } - final acf1 = nnDiffAcf1(runs); + // Welford (1962) centered sums/products, Task Force (1996) HRV. + // Accumulate difference moments and lagged co-moments in sample order. + var dn = 0, ln = 0; + var dm = 0.0, d2 = 0.0, ma = 0.0, mb = 0.0, co = 0.0; + for (final r in runs) { + for (var i = 0; i < r.length; i++) { + final d = r[i]; + dn++; + final delta = d - dm; + dm += delta / dn; + d2 += delta * (d - dm); + if (i > 0) { + ln++; + final da = r[i - 1] - ma; + final db = d - mb; + ma += da / ln; + mb += db / ln; + co += da * (d - mb); + } + } + } + final double? acf1 = dn >= 30 && d2 > 0 + ? (co + ln * (ma - dm) * (mb - dm)) / d2 : null; final refused = _jitterRefused(acf1, runs); if (rmssds.isEmpty && !refused) { return const Metric.absent( @@ -1030,6 +925,24 @@ Metric sleepSessionRmssdDetail( int minDiffsPerWindow = kMinDiffsPerRmssdWindow, }) { const inputs = ['rr_sleep_window']; + // An unknown interval cannot support a finite HRV estimate. Refuse the + // input before segmentation or moments rather than discarding evidence. + // A reset may expose disagreement between accumulated RR and recorded + // sensor time. Do not turn overlapping segments into successive beats. + if (hasRegressingBeatTimes(rrTsMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } + if (rrMs.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite NN interval', + ); + } if (startSec <= 0 || endSec <= startSec || rrMs.isEmpty || @@ -1042,6 +955,13 @@ Metric sleepSessionRmssdDetail( ); } + if (rrTsMs.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } final buckets = >{}; final bucketsTs = >{}; final inRr = []; @@ -1143,7 +1063,29 @@ Metric sleepSessionRmssdDetail( // THE HEADLINE nightly RMSSD (→ ln_rmssd → readiness). When the differences // are noise, the honest output is no headline, not a plausible one — the // readiness composite already treats a null HRV driver as absent. - final acf1 = nnDiffAcf1(runs); + // Welford (1962) centered sums/products, Task Force (1996) HRV. + // Accumulate difference moments and lagged co-moments in sample order. + var dn = 0, ln = 0; + var dm = 0.0, d2 = 0.0, ma = 0.0, mb = 0.0, co = 0.0; + for (final r in runs) { + for (var i = 0; i < r.length; i++) { + final d = r[i]; + dn++; + final delta = d - dm; + dm += delta / dn; + d2 += delta * (d - dm); + if (i > 0) { + ln++; + final da = r[i - 1] - ma; + final db = d - mb; + ma += da / ln; + mb += db / ln; + co += da * (d - mb); + } + } + } + final double? acf1 = dn >= 30 && d2 > 0 + ? (co + ln * (ma - dm) * (mb - dm)) / d2 : null; final refused = _jitterRefused(acf1, runs); if (rmssds.isEmpty && !refused) { return Metric.absent( @@ -1282,3 +1224,22 @@ List> _cleanWindowRuns(List rr, List ts) { if (run.isNotEmpty) runs.add(run); return runs; } + +// Welford, B. P. (1962), Technometrics 4(3):419-420. +// Origin-shifted batch recurrence; skips NaN just as #99's Moments does. +({double? mean, double? sd}) _centeredMoments(List values) { + var n = 0; + var origin = 0.0, offset = 0.0, m2 = 0.0; + for (final value in values) { + if (value.isNaN) continue; + if (!value.isFinite) throw ArgumentError.value(value, 'value', 'Must be finite'); + if (n == 0) { origin = value; n = 1; continue; } + final centered = value - origin; + final delta = centered - offset; + n++; + offset += delta / n; + m2 += delta * (centered - offset); + } + return (mean: n == 0 ? null : origin + offset, + sd: n < 2 ? null : math.sqrt(m2 / (n - 1))); +} diff --git a/lib/src/onehz/clinical/load_trimp.dart b/lib/src/onehz/clinical/load_trimp.dart index 3c89ca4f..53187efa 100644 --- a/lib/src/onehz/clinical/load_trimp.dart +++ b/lib/src/onehz/clinical/load_trimp.dart @@ -63,7 +63,8 @@ Metric banisterTrimp( ); } final reserve = maxHr - restingHr; - var trimp = 0.0; + // Neumaier (1974), Rundungsfehleranalyse einiger Verfahren zur Summation. + var trimp = 0.0, correction = 0.0; for (final hr in hrPerMin) { if (!hr.isFinite || hr <= 0) continue; // off-skin + non-finite guard var hrr = (hr - restingHr) / reserve; @@ -71,10 +72,16 @@ Metric banisterTrimp( if (hrr > 1) hrr = 1; // 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); + final value = hrr * StrainScorer.banisterY(hrr, female: sex == Sex.female); + final next = trimp + value; + if (next.isFinite) { + correction += trimp.abs() >= value.abs() + ? (trimp - next) + value : (value - next) + trimp; + } + trimp = next; } return Metric( - value: trimp, + value: trimp + correction, confidence: 0.6, tier: Tier.estimate, inputs_used: inputs, @@ -829,7 +836,8 @@ 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. +/// An observed rest day is 0 load and decays the EWMA. Unknown/non-finite +/// daily load abstains; it is never substituted with a rest day. /// /// SEEDING (Banister 1975 impulse-response; the load before the record started /// is UNKNOWN): both accumulators used to be seeded at `dailyTrimp.first`, @@ -846,6 +854,13 @@ Metric ctlAtlTsb(List dailyTrimp, int minDays = ctlAtlMinDays, int primeDays = 7}) { const inputs = ['daily_trimp']; + if (dailyTrimp.any((value) => !value.isFinite)) { + return const Metric.absent( + tier: Tier.estimate, + inputs_used: inputs, + note: 'non-finite daily TRIMP', + ); + } if (dailyTrimp.length < minDays) { return Metric.absent( tier: Tier.estimate, @@ -856,13 +871,31 @@ Metric ctlAtlTsb(List dailyTrimp, final lc = 1 - math.exp(-1 / ctlDays); final la = 1 - math.exp(-1 / atlDays); 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); + // Neumaier (1974) compensation for the bounded-prefix mean. + var sum = 0.0, correction = 0.0; + for (final value in dailyTrimp.take(prime)) { + final next = sum + value; + if (next.isFinite) { + correction += sum.abs() >= value.abs() + ? (sum - next) + value : (value - next) + sum; + } + sum = next; + } + final seed = (sum + correction) / prime; + // Banister et al. (1975) recurrence evaluated as an affine map. + // f(x) = scale*x + offset, with keep = 1-lambda. + double evaluate(double lambda) { + var scale = 1.0, offset = 0.0; + final keep = 1 - lambda; + for (var i = prime; i < dailyTrimp.length; i++) { + final value = dailyTrimp[i]; + scale *= keep; + offset = keep * offset + lambda * value; + } + return scale * seed + offset; } + final ctl = evaluate(lc); + final atl = evaluate(la); final conf = (dailyTrimp.length / 42.0).clamp(0.3, 0.85); return Metric( value: LoadState(ctl, atl, ctl - atl), diff --git a/lib/src/onehz/clinical/rmssd_gate.dart b/lib/src/onehz/clinical/rmssd_gate.dart new file mode 100644 index 00000000..6c8972df --- /dev/null +++ b/lib/src/onehz/clinical/rmssd_gate.dart @@ -0,0 +1,319 @@ +// Rejection and breathing-line exemption thresholds are an internal rule. +// Welch 1967 spectral estimation; Hirsch & Bishop 1981 provides the +// breathing/RSA basis. This is hrvTime()'s existing internal jitter policy, +// ESTIMATE with low confidence, not a published classifier. Shared formulas +// preserve the calculation; reordered sums can change decisions within rounding +// of a threshold. +// CLINICAL — the spectral half of the RMSSD jitter gate, as pieces. +// +// `hrvTime`, `nocturnalRmssd` and `sleepSessionWindowedRmssd` (hrv_time.dart) +// read the successive-difference runs themselves when the lag-1 ACF is under +// [kNnDiffAcf1Floor]: a respiratory line on a low white floor exempts the night +// ([nnDiffNoiseShare]), and a refused night may still publish from the 5-min +// windows on a breathing line steady in Hz. Both reduce to sums over Welch +// segments, so this file holds them as SUMS plus the decisions over those +// sums, written once: the batch fills the sums from the runs in hand, and +// `hrvTimeAggregator` fills the same sums from mergeable per-chunk state. +// +// Not exported from `onehz.dart`. + +import 'dart:math' as math; + +import '../util.dart'; +import 'hrv_time.dart' show kNnDiffNoiseShareCeiling; + +/// Welch segment length (beats) of the noise-share spectrum, the first and last +/// bins of its 0.15–0.5 cycles/beat band. +const int kNoiseSegment = 64, kNoiseLoBin = 10, kNoiseHiBin = kNoiseSegment ~/ 2; + +/// The segment's Hann window (symmetric, n − 1 in the denominator). +final List noiseHann = [ + for (var i = 0; i < kNoiseSegment; i++) + 0.5 - 0.5 * math.cos(2 * math.pi * i / (kNoiseSegment - 1)) +]; + +final double _noiseHannSq = () { + var w2 = 0.0; + for (final x in noiseHann) { + w2 += x * x; + } + return w2; +}(); + +final List> _cosT = [ + for (var k = kNoiseLoBin; k <= kNoiseHiBin; k++) + [for (var i = 0; i < kNoiseSegment; i++) math.cos(2 * math.pi * k * i / kNoiseSegment)] +]; +final List> _sinT = [ + for (var k = kNoiseLoBin; k <= kNoiseHiBin; k++) + [for (var i = 0; i < kNoiseSegment; i++) math.sin(2 * math.pi * k * i / kNoiseSegment)] +]; + +/// Adds the Welch segment `x[s .. s + 63]` (RR levels integrated from the +/// differences) to [psd] (bins [kNoiseLoBin]..[kNoiseHiBin]) and returns its +/// band power. +double addNoiseSegment(List x, int s, List psd) { + const n = kNoiseSegment; + var m = 0.0; + for (var i = 0; i < n; i++) { + m += x[s + i]; + } + m /= n; + final v = [for (var i = 0; i < n; i++) (x[s + i] - m) * noiseHann[i]]; + var p = 0.0; + for (var k = kNoiseLoBin; k <= kNoiseHiBin; k++) { + final c = _cosT[k - kNoiseLoBin], sn = _sinT[k - kNoiseLoBin]; + var re = 0.0, im = 0.0; + for (var i = 0; i < n; i++) { + re += v[i] * c[i]; + im -= v[i] * sn[i]; + } + psd[k] += re * re + im * im; + p += re * re + im * im; + } + return p; +} + +/// The sums [nnDiffNoiseShare] decides on. +class NoiseShareSums { + final List psd = List.filled(kNoiseHiBin + 1, 0.0); + int segs = 0, nd = 0; + double ssd = 0.0, wsd = 0.0, wsum = 0.0, pSum = 0.0, p2Sum = 0.0; + + /// One contiguous run of differences, in the batch's own order. + void addRun(List r) { + const n = kNoiseSegment; + for (final d in r) { + ssd += d * d; + nd++; + } + // Integrate the run back to RR levels (up to a constant the mean removes). + final x = List.filled(r.length + 1, 0.0); + for (var i = 0; i < r.length; i++) { + x[i + 1] = x[i] + r[i]; + } + final cov = List.filled(x.length, 0.0); // Σ w² per sample + for (var s = 0; s + n <= x.length; s += n ~/ 2) { + for (var i = 0; i < n; i++) { + cov[s + i] += noiseHann[i] * noiseHann[i]; + } + final p = addNoiseSegment(x, s, psd); + pSum += p; + p2Sum += p * p; + segs++; + } + for (var i = 0; i < r.length; i++) { + final c = (cov[i] + cov[i + 1]) / 2; + wsd += c * r[i] * r[i]; + wsum += c; + } + } +} + +/// The white-noise share from [s] (see [nnDiffNoiseShare] in hrv_time.dart for +/// what it measures and every refusal). [lattice] says whether the differences +/// sit on a coarse beat-time grid ([coarseLattice]); it is only asked when +/// the counts allow a verdict. +double? noiseShareFrom(NoiseShareSums s, + {required bool Function() lattice, int minSegments = 20}) { + final segs = s.segs; + if (segs < minSegments || s.ssd == 0 || s.wsd == 0) return null; + if (lattice()) return null; + // One loud stretch dominates the average: count segments by power, not number. + if (s.pSum * s.pSum / s.p2Sum < 0.75 * minSegments) return null; + final band = s.psd.sublist(kNoiseLoBin); + var peak = 0; + for (var k = 1; k < band.length; k++) { + if (band[k] > band[peak]) peak = k; + } + if (peak + kNoiseLoBin >= kNoiseHiBin - 1) return null; // ~0.477–0.5 band + if (band.last > 0.2 * band[peak]) return null; // hump on Nyquist + final med = median(band)!; + if (band[peak] < 4.5 * med) return null; // no line, no exemption + final floor = med / segs / _noiseHannSq; + final vetted = s.wsd / s.wsum, all = s.ssd / s.nd; + // 1 − (structured share of RMSSD²); equals 2σ²/mean(d²) when stationary. + return 1 - (vetted - 2 * floor) / math.max(vetted, all); +} + +/// True when every successive difference is a whole multiple of one step q +/// (the smallest |d| over 1e-6) and the mean squared difference [msd] is under +/// 2q². [magnitudes] are the |d| values (repeats change nothing), in any +/// order. See hrv_time.dart `_onCoarseLattice` for why. +bool coarseLattice(Iterable magnitudes, double msd) { + var q = double.infinity; + for (final a in magnitudes) { + if (a > 1e-6 && a < q) q = a; + } + if (q.isInfinite || msd >= 2 * q * q) return false; + for (final a in magnitudes) { + final k = a / q; + if ((k - k.roundToDouble()).abs() > 1e-3) return false; + } + return true; +} + +/// The breathing-line test's Hann window (128 beats). +final List hann128 = [ + for (var i = 0; i < 128; i++) 0.5 - 0.5 * math.cos(2 * math.pi * i / 127) +]; +final double _hann128Sq = hann128.fold(0.0, (a, w) => a + w * w); + +/// The cycles-per-beat axis the breathing-line test pools on: 0.15–0.5 +/// cycles/beat at the night's median beat. +final List breathingAxis = [for (var j = 38; j <= 128; j++) j / 256]; + +/// Welch power (128-beat Hann, 50 % overlap, averaged over segments) of one +/// window's RR, rebuilt from its difference runs, at cycles-per-beat +/// frequencies [f]. Null when no run holds a full segment. +List? windowPsd(List> diffRuns, List f) { + const n = 128; + final p = List.filled(f.length, 0.0); + var segs = 0; + for (final r in diffRuns) { + final x = List.filled(r.length + 1, 0.0); + for (var i = 0; i < r.length; i++) { + x[i + 1] = x[i] + r[i]; + } + for (var s = 0; s + n <= x.length; s += n ~/ 2) { + var m = 0.0; + for (var i = 0; i < n; i++) { + m += x[s + i]; + } + m /= n; + for (var j = 0; j < f.length; j++) { + final c = math.cos(2 * math.pi * f[j]), sn = math.sin(2 * math.pi * f[j]); + var cr = 1.0, ci = 0.0, re = 0.0, im = 0.0; + for (var i = 0; i < n; i++) { + final v = (x[s + i] - m) * hann128[i]; + re += v * cr; + im += v * ci; + final t = cr * c - ci * sn; + ci = cr * sn + ci * c; + cr = t; + } + p[j] += re * re + im * im; + } + segs++; + } + } + return segs == 0 ? null : [for (final v in p) v / segs]; +} + +/// The breathing-line test over [windows] 5-min windows, each read through +/// callbacks so a caller can hold its windows as runs or as sums: +/// [psdAt] is the window's [windowPsd] at the given frequencies (null without +/// a full segment), [diffs] its difference count, [addSquares] adds its +/// squared differences to a running total (in the caller's order), [meanRrMs] +/// its mean RR. [allSsd] / [allNd] are the squared differences and their count +/// over every window, and [lattice] is [coarseLattice] over every window's +/// differences at the mean square it is given. Returns the indices of the +/// windows whose own peak sits on the line, or null. See hrv_time.dart +/// `_steadyBreathingWindows` for the method. +List? steadyBreathingFrom({ + required int windows, + required List meanRrMs, + required List? Function(int w, List f) psdAt, + required int Function(int w) diffs, + required double Function(int w, double acc) addSquares, + required double allSsd, + required int allNd, + required bool Function(double msd) lattice, +}) { + if (allNd == 0 || lattice(allSsd / allNd)) return null; + final ref = median([for (final m in meanRrMs) if (m > 0) m]); + if (ref == null) return null; + final fc = breathingAxis; + final byBeat = List.filled(fc.length, 0.0); + final byHz = List.filled(fc.length, 0.0); + final cover = List.filled(fc.length, 0); + final peakHz = {}; // window -> its own peak, cycles per ms + final psd = >{}; // window -> its cycles-per-beat PSD + for (var w = 0; w < windows; w++) { + if (meanRrMs[w] <= 0) continue; + final pb = psdAt(w, fc); + if (pb == null) continue; + final norm = median(pb)!; + if (norm <= 0) continue; + final k = meanRrMs[w] / ref; + final fh = [for (final f in fc) f * k]; + final ph = psdAt(w, fh)!; + var best = -1; + for (var j = 0; j < fc.length; j++) { + byBeat[j] += pb[j] / norm; + if (fh[j] > 0.5) continue; + byHz[j] += ph[j] / norm; + cover[j]++; + if (best < 0 || ph[j] > ph[best]) best = j; + } + if (best >= 0) peakHz[w] = fc[best] / ref; + psd[w] = pb; + } + final nw = peakHz.length; + if (nw < 12) return null; + final bins = [for (var j = 0; j < fc.length; j++) if (cover[j] == nw) j]; + if (bins.length < 20) return null; + final hz = [for (final j in bins) byHz[j]]; + final bt = [for (final j in bins) byBeat[j]]; + var pk = 0; + for (var i = 1; i < hz.length; i++) { + if (hz[i] > hz[pk]) pk = i; + } + if (pk == 0 || pk == hz.length - 1) return null; + final f0 = fc[bins[pk]] / ref; + final brpm = f0 * 60000; + if (brpm < 8 || brpm > 30) return null; + // Same bins, same normalisation: the line must stand taller aligned in Hz + // than the beat pooling stands anywhere within ±12 bins of it. + var near = 0.0; + for (var i = math.max(0, pk - 12); i <= math.min(bt.length - 1, pk + 12); i++) { + near = math.max(near, bt[i]); + } + // The pooled noise ripple shrinks as 1/√windows, so the bar does too. + final promHz = hz[pk] / median(hz)!; + if (promHz < math.max(2.5, 1.5 + 7 / math.sqrt(nw))) return null; + if (promHz < 1.1 * near / median(bt)!) return null; + // Alternation lives at Nyquist on the beat axis, outside these bins. + if (byBeat.last / median(bt)! >= 0.5 * promHz) return null; + final keep = [ + for (final e in peakHz.entries) + if (math.log(e.value / f0).abs() <= 0.15) e.key + ]; + if (keep.length < 6) return null; + // Jitter that fails the gate would still be most of the number: same + // ceiling as the noise share, over the kept windows (see hrv_time.dart). + var lo = double.infinity, hi = 0.0; + for (final w in keep) { + lo = math.min(lo, peakHz[w]!); + hi = math.max(hi, peakHz[w]!); + } + var sq = 0.0, nAll = 0; + final acc = List.filled(fc.length, 0.0); + final wt = List.filled(fc.length, 0.0); + for (final w in keep) { + sq = addSquares(w, sq); + final n = diffs(w); + nAll += n; + final a = lo * meanRrMs[w] - 3 / 128, b = hi * meanRrMs[w] + 3 / 128; + for (var j = 0; j < fc.length; j++) { + if (fc[j] > a && fc[j] < b) continue; + acc[j] += n * psd[w]![j] / (2 - 2 * math.cos(2 * math.pi * fc[j])); + wt[j] += n; + } + } + var floor = double.infinity; + for (var j = 8; j < fc.length - 8; j++) { + var s = 0.0, m = 0; + for (var i = j - 8; i <= j + 8; i++) { + if (wt[i] < nAll / 2) break; + s += acc[i] / wt[i]; + m++; + } + if (m == 17) floor = math.min(floor, s / 17); + } + if (floor == double.infinity) return null; + // The lowest stretch of a noisy curve reads below its mean; 1.1 puts the + // floor back at these pooled sizes, so jitter alone is not undercounted. + final noise = 1.1 * 6 * floor / _hann128Sq * nAll; + return noise < kNnDiffNoiseShareCeiling * sq ? keep : null; +} diff --git a/lib/src/onehz/foundations/rr_correction.dart b/lib/src/onehz/foundations/rr_correction.dart index 89840130..a8e1738f 100644 --- a/lib/src/onehz/foundations/rr_correction.dart +++ b/lib/src/onehz/foundations/rr_correction.dart @@ -1,3 +1,11 @@ +// Lipponen & Tarvainen 2019 classification with Peltola 2012 correction policy. +// The unpublished artifact-run cap is an internal rule, ESTIMATE with low +// confidence: +// over-cap runs become gaps and each side is re-corrected. Without a cap split +// the result equals uncapped correction. Shared kernel arithmetic and raw seam +// context preserve capped batch NN/classes; settled stream plus provisional tail +// equals capped batch. RrSandwich preserves source-beat absolute timestamps; +// uncapped arbitrarily long artifact runs have no exactness guarantee. // FOUNDATION — RR artifact correction. // // Lipponen & Tarvainen 2019 ("A robust algorithm for heart rate variability @@ -13,11 +21,22 @@ // This is PRV. Output = cleaned NN series + per-beat artifact mask + clean // fraction. Multi-beat gaps are never silently bridged. -import 'dart:math' as math; -import '../util.dart'; +import 'rr_correction_kernel.dart'; + +/// Longest artifact run (beats) the corrector will carry. A longer run is not a +/// correction problem but an absence: it breaks the series into independent +/// segments and its beats are left out (a gap). 32 beats leaves no run on the +/// owner's 18 days of data over the cap (longest seen: 16) and bounds the +/// streaming state; see combined-bump/brief.md. +const int kMaxArtifactRun = 32; /// Artifact label for a single beat. -enum BeatClass { normal, ectopic, longShort, missed, extra } +enum BeatClass { + normal, ectopic, longShort, missed, extra, + /// An over-cap artifact run, without a claim about its individual types. + /// Appended to preserve the indices of existing checkpoint labels. + gap, +} class RrCorrectionResult { /// Cleaned NN intervals (ms). Isolated artifacts spline-corrected; multi-beat @@ -34,7 +53,16 @@ class RrCorrectionResult { /// spliced out; pass the timestamps. final List nnTimesMs; + /// For each NN, the index of the input beat it came from: the beat itself, + /// or, for a spline-corrected beat, that beat's own index. Strictly + /// increasing; the input beats not listed here were dropped. Lets a caller + /// pair each NN with that beat's own (absolute) timestamp. Empty on a + /// hand-built result that did not track it. + final List nnSourceIndex; + /// Per-input-beat classification (same length as the input rr). + /// Over-cap runs use [BeatClass.gap]; other artifact labels retain their + /// inferred classification. A gap makes no physiological type claim. final List classes; /// Fraction of input beats classified normal (0..1). @@ -53,6 +81,7 @@ class RrCorrectionResult { required this.cleanFraction, required this.droppedCount, required this.correctedCount, + this.nnSourceIndex = const [], }); } @@ -74,6 +103,10 @@ class RrCorrectionResult { /// `rec_ts*1000` and therefore whole-second — and re-anchors to the real /// timestamp whenever the wall step exceeds the interval by more than /// [reanchorGapMs], i.e. at a dropout. +/// +/// [maxArtifactRun] caps a run of artifact beats: a longer run splits the series +/// (output = correctRr(before) ++ correctRr(after), times absolute, the run's +/// beats absent). null = uncapped, the pre-cap behaviour. RrCorrectionResult correctRr( List rrMs, { List? rrTsMs, @@ -81,7 +114,55 @@ RrCorrectionResult correctRr( int windowBeats = 91, double minThresholdMs = 100, double reanchorGapMs = 1000, + int? maxArtifactRun = kMaxArtifactRun, }) { + if (maxArtifactRun != null && maxArtifactRun < 1) { + throw ArgumentError.value(maxArtifactRun, 'maxArtifactRun'); + } + if (windowBeats < 1) throw ArgumentError.value(windowBeats, 'windowBeats'); + if (maxArtifactRun != null) { + // The overwhelmingly common path has no run over the cap. Reuse the + // original vectorized correction verbatim, including its clock, and only + // evaluate split windows when a run actually requires a cut. + var hardRun = 0, hardGap = false; + for (final rr in rrMs) { + hardRun = rr < 300 || rr > 2000 ? hardRun + 1 : 0; + if (hardRun > maxArtifactRun) { hardGap = true; break; } + } + if (!hardGap) { + final ordinary = correctRr(rrMs, rrTsMs: rrTsMs, alpha: alpha, + windowBeats: windowBeats, minThresholdMs: minThresholdMs, + reanchorGapMs: reanchorGapMs, maxArtifactRun: null); + var run = 0, gap = false; + for (final cls in ordinary.classes) { + run = cls == BeatClass.normal ? 0 : run + 1; + if (run > maxArtifactRun) { gap = true; break; } + } + if (!gap) return ordinary; + } + final w = RrCapWindow.evaluate(rrMs, + ts: rrTsMs, cap: maxArtifactRun, alpha: alpha, + windowBeats: windowBeats, minThresholdMs: minThresholdMs, + reanchorGapMs: reanchorGapMs); + final nn = [], times = [], source = []; + var normal = 0, dropped = 0, corrected = 0; + for (var i = 0; i < rrMs.length; i++) { + if (w.classes[i] == BeatClass.normal) normal++; + final v = w.values[i]; + if (v == null) { + dropped++; + } else { + nn.add(v); + times.add(w.times[i]); + source.add(i); + if (w.classes[i] != BeatClass.normal) corrected++; + } + } + return RrCorrectionResult(nn: nn, nnTimesMs: times, classes: w.classes, + cleanFraction: rrMs.isEmpty ? 0 : normal / rrMs.length, + droppedCount: dropped, correctedCount: corrected, + nnSourceIndex: source); + } final n = rrMs.length; // Beat-end time for EVERY input beat, kept or not. One array, built once, so // every branch below reads a clock instead of advancing its own. @@ -105,10 +186,12 @@ RrCorrectionResult correctRr( ]; final nn = []; final times = []; + final source = []; for (var i = 0; i < n; i++) { if (classes[i] == BeatClass.normal) { nn.add(rrMs[i]); times.add(tAll[i]); + source.add(i); } } return RrCorrectionResult( @@ -118,6 +201,7 @@ RrCorrectionResult correctRr( cleanFraction: nn.length / n, droppedCount: n - nn.length, correctedCount: 0, + nnSourceIndex: source, ); } @@ -133,32 +217,13 @@ RrCorrectionResult correctRr( // medRR: rr minus local median (for missed/extra long-range tests). final med = _slidingMedian(rrMs, windowBeats); - final mRR = List.generate(n, (i) { - final d = rrMs[i] - med[i]; - return d < 0 ? d * 2 : d; // paper asymmetry weight - }); + final mRR = List.generate(n, (i) => medianDeviation(rrMs[i], med[i])); final th2 = _slidingThreshold(mRR, windowBeats, alpha, minThresholdMs); - final classes = List.filled(n, BeatClass.normal); - for (var i = 0; i < n; i++) { - final hardLong = rrMs[i] > 2000; - final hardShort = rrMs[i] < 300; - final bigJump = dRR[i].abs() > th1[i]; - final bigDev = mRR[i].abs() > th2[i]; - if (hardLong || (bigDev && mRR[i] > 0)) { - // Long interval: likely a MISSED beat (interval ~ multiple of normal). - classes[i] = (med[i] > 0 && rrMs[i] > 1.5 * med[i]) - ? BeatClass.missed - : BeatClass.longShort; - } else if (hardShort || (bigDev && mRR[i] < 0)) { - // Short interval: likely an EXTRA (spurious) beat. - classes[i] = (med[i] > 0 && rrMs[i] < 0.6 * med[i]) - ? BeatClass.extra - : BeatClass.longShort; - } else if (bigJump) { - classes[i] = BeatClass.ectopic; - } - } + final classes = [ + for (var i = 0; i < n; i++) + classifyBeat(rrMs[i], dRR[i], th1[i], med[i], mRR[i], th2[i]) + ]; // Compensatory-pair reconciliation. A single ectopic/extra beat shows up as // TWO successive dRR spikes of opposite sign (the bad beat, then the @@ -169,17 +234,13 @@ RrCorrectionResult correctRr( // single isolated artifact rather than a spurious 2-beat run. for (var k = 1; k < n; k++) { if (classes[k] != BeatClass.ectopic) continue; - final prevBad = classes[k - 1] == BeatClass.extra || - classes[k - 1] == BeatClass.missed || - classes[k - 1] == BeatClass.ectopic || - classes[k - 1] == BeatClass.longShort; - if (!prevBad) continue; - final oppositeSign = dRR[k] * dRR[k - 1] < 0; - final valueNormal = med[k] > 0 && - rrMs[k] >= 300 && - rrMs[k] <= 2000 && - (rrMs[k] - med[k]).abs() <= 0.2 * med[k]; - if (oppositeSign && valueNormal) { + final prevBad = classes[k - 1] != BeatClass.normal; + if (isRecoveryBeat( + prevIsArtifact: prevBad, + rr: rrMs[k], + dRR: dRR[k], + dRRPrev: dRR[k - 1], + med: med[k])) { classes[k] = BeatClass.normal; } } @@ -189,6 +250,7 @@ RrCorrectionResult correctRr( final isArtifact = [for (final c in classes) c != BeatClass.normal]; final nn = []; final times = []; + final source = []; var dropped = 0; var corrected = 0; var i = 0; @@ -196,6 +258,7 @@ RrCorrectionResult correctRr( if (!isArtifact[i]) { nn.add(rrMs[i]); times.add(tAll[i]); + source.add(i); i++; continue; } @@ -216,6 +279,7 @@ RrCorrectionResult correctRr( // — 1.6 % of a night at one missed beat in sixty. nn.add(corr); times.add(tAll[i]); + source.add(i); corrected++; } else { dropped++; // no anchors -> honest drop; the clock already elapsed @@ -241,6 +305,7 @@ RrCorrectionResult correctRr( cleanFraction: normalCount / n, droppedCount: dropped, correctedCount: corrected, + nnSourceIndex: source, ); } @@ -306,42 +371,63 @@ List _slidingThreshold( final n = x.length; final out = List.filled(n, 0); final half = win ~/ 2; + final w = _SortedWindow(x, half); for (var i = 0; i < n; i++) { - final lo = math.max(0, i - half); - final hi = math.min(n - 1, i + half); - final seg = []; - for (var k = lo; k <= hi; k++) { - 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 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); + w.advanceTo(i); + out[i] = thresholdOfSorted(w.a, alpha, floor); } return out; } +/// Median of the window around each beat EXCLUDING the beat itself (falls back +/// to the beat when the window holds nothing else). Same value as sorting the +/// window-minus-self per beat, but the window is kept sorted as it slides and +/// the self element is skipped by index arithmetic instead of copied out. List _slidingMedian(List x, int win) { final n = x.length; final out = List.filled(n, 0); final half = win ~/ 2; + final w = _SortedWindow(x, half); for (var i = 0; i < n; i++) { - final lo = math.max(0, i - half); - final hi = math.min(n - 1, i + half); - final seg = []; - for (var k = lo; k <= hi; k++) { - if (k == i) continue; - seg.add(x[k]); - } - out[i] = median(seg) ?? x[i]; + w.advanceTo(i); + out[i] = medianExcluding(w.a, x[i]) ?? x[i]; } return out; } +/// Ascending-sorted copy of `x[max(0,i-half) .. min(n-1,i+half)]`, maintained +/// incrementally as [advanceTo] steps i forward by one (O(log w) search plus a +/// w-element memmove per beat instead of an O(w log w) sort). +class _SortedWindow { + final List _x; + final int _half; + final _set = SortedMultiset(); + List get a => _set.a; + int _lo = 0; + int _hi = -1; // inclusive; -1 = empty + + _SortedWindow(this._x, this._half); + + void advanceTo(int i) { + final n = _x.length; + final newLo = i - _half < 0 ? 0 : i - _half; + final newHi = i + _half > n - 1 ? n - 1 : i + _half; + if (_hi < _lo) { + // First call: build once. + _set.resetFrom(_x.getRange(newLo, newHi + 1)); + } else { + while (_lo < newLo) { + _set.remove(_x[_lo++]); + } + while (_hi < newHi) { + _set.add(_x[++_hi]); + } + } + _lo = newLo; + _hi = newHi; + } +} + /// Catmull-Rom cubic interpolation at the artifact index using the nearest two /// NORMAL beats on each side. Returns null if anchors are unavailable. double? _splineCorrect(List rr, List isArtifact, int idx) { @@ -353,19 +439,5 @@ double? _splineCorrect(List rr, List isArtifact, int idx) { for (var k = idx + 1; k < rr.length && right.length < 2; k++) { if (!isArtifact[k]) right.add(rr[k]); } - if (left.isEmpty || right.isEmpty) return null; - final p1 = left.last; - final p2 = right.first; - final p0 = left.length >= 2 ? left.first : p1; - final p3 = right.length >= 2 ? right.last : p2; - // Catmull-Rom at t=0.5 between p1 and p2. - const t = 0.5; - final t2 = t * t; - final t3 = t2 * t; - final v = 0.5 * - ((2 * p1) + - (-p0 + p2) * t + - (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 + - (-p0 + 3 * p1 - 3 * p2 + p3) * t3); - return v; + return splineMid(left, right); } diff --git a/lib/src/onehz/foundations/rr_correction_kernel.dart b/lib/src/onehz/foundations/rr_correction_kernel.dart new file mode 100644 index 00000000..c0eeb7cf --- /dev/null +++ b/lib/src/onehz/foundations/rr_correction_kernel.dart @@ -0,0 +1,376 @@ +// Lipponen & Tarvainen 2019 classification with Peltola 2012 correction policy. +// The unpublished artifact-run cap is an internal rule, ESTIMATE with low +// confidence: +// over-cap runs become gaps and each side is re-corrected. Without a cap split +// the result equals uncapped correction. Shared kernel arithmetic and raw seam +// context preserve capped batch NN/classes; settled stream plus provisional tail +// equals capped batch. RrSandwich preserves source-beat absolute timestamps; +// uncapped arbitrarily long artifact runs have no exactness guarantee. +// FOUNDATION — the per-beat pieces of Lipponen–Tarvainen correction that both +// `correctRr` (whole series) and `RrCorrector` (streaming) run. +// +// One copy of the sorted window, the threshold, the leave-one-out median, the +// class decision, the compensatory-pair rule and the spline, so the batch and +// the stream cannot drift apart: the stream's contract is bit-identity with +// `correctRr`, and sharing the arithmetic is what keeps that cheap to hold. +// +// Not exported from `onehz.dart` — internal to foundations/. + +import 'dart:math' as math; +import '../util.dart'; +import 'rr_correction.dart' show BeatClass; + +/// First index whose element is not less than [v] under `compareTo` — the same +/// total order `List.sort()` uses (-0.0 < 0.0, NaN last), so a window +/// kept in this order matches the one a per-beat sort produced. +int lowerBound(List a, double v) { + var lo = 0; + var hi = a.length; + while (lo < hi) { + final m = (lo + hi) >> 1; + if (a[m].compareTo(v) < 0) { + lo = m + 1; + } else { + hi = m; + } + } + return lo; +} + +/// An ascending multiset of doubles: O(log w) search plus a w-element memmove +/// per add / remove, instead of an O(w log w) sort per beat. +class SortedMultiset { + /// The values, ascending. Read-only to callers. + final List a = []; + + void add(double v) => a.insert(lowerBound(a, v), v); + + /// Removes one occurrence of [v], which must be present. + void remove(double v) => a.removeAt(lowerBound(a, v)); + + void resetFrom(Iterable xs) { + a + ..clear() + ..addAll(xs) + ..sort(); + } +} + +/// `alpha × QD` of an ascending window, floored: QD = (Q3 − Q1)/2. +/// +/// 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. +double thresholdOfSorted(List sorted, double alpha, double floor) { + final q1 = percentileSorted(sorted, 25) ?? 0; + final q3 = percentileSorted(sorted, 75) ?? 0; + final qd = (q3 - q1) / 2; + return math.max(alpha * qd, floor); +} + +/// `median(sorted minus one occurrence of v)` — [percentile]'s linear +/// interpolation at p = 50 over the sorted list with v's slot skipped, in the +/// same operation order so the result is bit-identical. Null if nothing is left. +double? medianExcluding(List sorted, double v) { + final len = sorted.length - 1; + if (len <= 0) return null; + final skip = lowerBound(sorted, v); + double at(int j) => j < skip ? sorted[j] : sorted[j + 1]; + if (len == 1) return at(0); + final rank = (50.0 / 100) * (len - 1); + final lo = rank.floor(); + final hi = rank.ceil(); + if (lo == hi) return at(lo); + final frac = rank - lo; + return at(lo) + (at(hi) - at(lo)) * frac; +} + +/// The paper's asymmetric deviation from the local median: below-median +/// deviations count double. +double medianDeviation(double rr, double med) { + final d = rr - med; + return d < 0 ? d * 2 : d; +} + +/// First-pass class of one beat from its dRR, mRR and the two thresholds. +BeatClass classifyBeat( + double rr, double dRR, double th1, double med, double mRR, double th2) { + final hardLong = rr > 2000; + final hardShort = rr < 300; + final bigJump = dRR.abs() > th1; + final bigDev = mRR.abs() > th2; + if (hardLong || (bigDev && mRR > 0)) { + // Long interval: likely a MISSED beat (interval ~ multiple of normal). + return (med > 0 && rr > 1.5 * med) ? BeatClass.missed : BeatClass.longShort; + } else if (hardShort || (bigDev && mRR < 0)) { + // Short interval: likely an EXTRA (spurious) beat. + return (med > 0 && rr < 0.6 * med) ? BeatClass.extra : BeatClass.longShort; + } else if (bigJump) { + return BeatClass.ectopic; + } + return BeatClass.normal; +} + +/// Compensatory-pair rule: beat k was flagged ONLY by its dRR jump, the beat +/// before it was an artifact, the two dRR spikes have opposite sign and k's own +/// value is normal and close to the local median — k is just the recovery from +/// the previous beat's event, so it is demoted to normal. +bool isRecoveryBeat({ + required bool prevIsArtifact, + required double rr, + required double dRR, + required double dRRPrev, + required double med, +}) { + if (!prevIsArtifact) return false; + final oppositeSign = dRR * dRRPrev < 0; + final valueNormal = + med > 0 && rr >= 300 && rr <= 2000 && (rr - med).abs() <= 0.2 * med; + return oppositeSign && valueNormal; +} + +/// Catmull-Rom at t = 0.5 between the nearest NORMAL beats either side of an +/// isolated artifact. [left] / [right] hold up to two values each, nearest to +/// the artifact LAST in [left] and FIRST in [right]. Null without an anchor on +/// both sides. +double? splineMid(List left, List right) { + if (left.isEmpty || right.isEmpty) return null; + final p1 = left.last; + final p2 = right.first; + final p0 = left.length >= 2 ? left.first : p1; + final p3 = right.length >= 2 ? right.last : p2; + const t = 0.5; + final t2 = t * t; + final t3 = t2 * t; + return 0.5 * + ((2 * p1) + + (-p0 + p2) * t + + (2 * p0 - 5 * p1 + 4 * p2 - p3) * t2 + + (-p0 + 3 * p1 - 3 * p2 + p3) * t3); +} + +/// Internal bounded-window evaluator shared by the capped batch and stream. +/// [from] has its complete left window in [rr], or is the segment's first beat. +/// Only classes at and after [from] are evaluated; the preceding class and +/// normal anchors are carried explicitly, so truncating old context cannot +/// change recovery or spline decisions. +class RrCapWindow { + final List classes; + final List values; + final List times; + final List gaps; + final List hardGaps; + RrCapWindow._(int n) + : classes = List.filled(n, BeatClass.normal), + values = List.filled(n, null), + times = List.filled(n, 0), + gaps = List.filled(n, false), + hardGaps = List.filled(n, false); + + static RrCapWindow evaluate( + List rr, { + List? ts, + required int cap, + double alpha = 5.2, + int windowBeats = 91, + double minThresholdMs = 100, + double reanchorGapMs = 1000, + int from = 0, + BeatClass previous = BeatClass.normal, + List leftNormals = const [], + bool leadingGap = false, + bool leadingHardGap = false, + int segmentStart = 0, + double clock = 0, + double? origin, + }) { + final out = RrCapWindow._(rr.length); + if (rr.isEmpty) return out; + final wall = ts != null && ts.length == rr.length; + final base = origin ?? (wall ? ts[0] - rr[0] : 0.0); + final h = windowBeats ~/ 2; + + List classify(int a, int b, int start, BeatClass prev) { + final n = b - a; + if (n < 3) { + return [ + for (var i = start; i < b; i++) + rr[i] >= 300 && rr[i] <= 2000 + ? BeatClass.normal + : BeatClass.longShort + ]; + } + final med = List.filled(n, 0); + final mrr = List.filled(n, 0); + final th1 = List.filled(n, 0); + final wd = SortedMultiset(), wr = SortedMultiset(); + var lo = a, hi = a; + double d(int i) => i == a ? 0 : rr[i] - rr[i - 1]; + for (var i = a; i < b; i++) { + final nextLo = math.max(a, i - h), nextHi = math.min(b, i + h + 1); + while (hi < nextHi) { + wd.add(d(hi)); + wr.add(rr[hi++]); + } + while (lo < nextLo) { + wd.remove(d(lo)); + wr.remove(rr[lo++]); + } + med[i - a] = medianExcluding(wr.a, rr[i]) ?? rr[i]; + mrr[i - a] = medianDeviation(rr[i], med[i - a]); + th1[i - a] = thresholdOfSorted(wd.a, alpha, minThresholdMs); + } + final wm = SortedMultiset(); + lo = a; + hi = a; + final cl = []; + for (var i = start; i < b; i++) { + final nextLo = math.max(a, i - h), nextHi = math.min(b, i + h + 1); + if (i == start) lo = nextLo; + if (i == start) hi = nextLo; + while (hi < nextHi) wm.add(mrr[hi++ - a]); + while (lo < nextLo) wm.remove(mrr[lo++ - a]); + var c = classifyBeat(rr[i], d(i), th1[i - a], med[i - a], mrr[i - a], + thresholdOfSorted(wm.a, alpha, minThresholdMs)); + if (i > a && + c == BeatClass.ectopic && + isRecoveryBeat( + prevIsArtifact: prev != BeatClass.normal, + rr: rr[i], + dRR: d(i), + dRRPrev: d(i - 1), + med: med[i - a])) { + c = BeatClass.normal; + } + cl.add(c); + prev = c; + } + return cl; + } + + // A stack avoids recursion depth proportional to the number of gaps. + final pending = <(int, int, int, BeatClass, List, bool, bool)>[ + ( + segmentStart, + rr.length, + from, + previous, + leftNormals, + leadingGap, + leadingHardGap + ) + ]; + while (pending.isNotEmpty) { + final (a, b, start, prev, anchors, inGap, hardLeading) = + pending.removeLast(); + if (start >= b) continue; + if (inGap && hardLeading && rr[start] >= 300 && rr[start] <= 2000) { + pending + .add((start, b, start, BeatClass.normal, const [], false, false)); + continue; + } + var hardCut = false; + var cutA = -1, cutB = -1; + // Physiologically impossible intervals are artifacts independently of + // every dispersion window. Cut these first: otherwise their medians can + // absorb valid beats on either side into the detected run itself. + for (var i = start; i < b;) { + if (rr[i] >= 300 && rr[i] <= 2000) { + i++; + continue; + } + var j = i + 1; + while (j < b && (rr[j] < 300 || rr[j] > 2000)) { + j++; + } + if (j - i > cap || (inGap && i == start)) { + cutA = i; + cutB = j; + hardCut = true; + break; + } + i = j; + } + final cl = cutA >= 0 ? [] : classify(a, b, start, prev); + if (cutA < 0 && inGap && cl.first == BeatClass.normal) { + pending + .add((start, b, start, BeatClass.normal, const [], false, false)); + continue; + } + if (cutA < 0) { + for (var i = start; i < b;) { + if (cl[i - start] == BeatClass.normal) { + i++; + continue; + } + var j = i + 1; + while (j < b && cl[j - start] != BeatClass.normal) { + j++; + } + if (j - i > cap || (inGap && i == start)) { + cutA = i; + cutB = j; + break; + } + i = j; + } + } + if (cutA >= 0) { + for (var i = cutA; i < cutB; i++) { + out.classes[i] = BeatClass.gap; + out.gaps[i] = true; + out.hardGaps[i] = hardCut || (inGap && hardLeading && cutA == start); + } + // Cut BEFORE classifying either side, including their dispersion + // windows. The gap's labels are canonical, never provisional medians. + pending.add((cutB, b, cutB, BeatClass.normal, const [], false, false)); + pending.add((a, cutA, start, prev, anchors, inGap, hardLeading)); + continue; + } + for (var i = start; i < b; i++) { + out.classes[i] = cl[i - start]; + } + final normals = List.of(anchors); + for (var i = start; i < b;) { + if (out.classes[i] == BeatClass.normal) { + out.values[i] = rr[i]; + normals.add(rr[i]); + if (normals.length > 2) normals.removeAt(0); + i++; + continue; + } + var j = i + 1; + while (j < b && out.classes[j] != BeatClass.normal) { + j++; + } + if (j == i + 1) { + final right = []; + for (var k = j; k < b && right.length < 2; k++) { + if (out.classes[k] == BeatClass.normal) right.add(rr[k]); + } + out.values[i] = splineMid(normals, right); + } + i = j; + } + } + var cur = clock; + var gap = leadingGap; + for (var i = from; i < rr.length; i++) { + final reset = gap && !out.gaps[i]; + final dropout = + wall && i > 0 && ts[i] - ts[i - 1] - rr[i] > reanchorGapMs; + final anchored = wall ? ts[i] - base : 0.0; + // Independent segments restart on recorded time, even if invalid RR + // accumulated ahead of it. A monotonic clamp would invent timestamps. + // HRV consumers refuse any regression in the retained beat clock. + cur = wall && reset + ? anchored + : dropout && anchored > cur + ? anchored + : cur + rr[i]; + out.times[i] = cur; + gap = out.gaps[i]; + } + return out; + } +} diff --git a/lib/src/onehz/respiration/cvhr_apnea.dart b/lib/src/onehz/respiration/cvhr_apnea.dart index fe1ada73..348c2aff 100644 --- a/lib/src/onehz/respiration/cvhr_apnea.dart +++ b/lib/src/onehz/respiration/cvhr_apnea.dart @@ -127,6 +127,20 @@ Metric cvhrApneaScreen( } final tSec = [for (final t in nnTimesMs) t / 1000.0]; + if (tSec.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } if (tSec.last - tSec.first <= 0) { return const Metric.absent( tier: Tier.high, diff --git a/lib/src/onehz/respiration/resp_rate.dart b/lib/src/onehz/respiration/resp_rate.dart index 6cd3452a..8e57c1d9 100644 --- a/lib/src/onehz/respiration/resp_rate.dart +++ b/lib/src/onehz/respiration/resp_rate.dart @@ -187,6 +187,20 @@ Metric rsaRespRate( ); } final tSec = [for (final t in nnTimesMs) t / 1000.0]; + if (tSec.any((t) => !t.isFinite)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'non-finite beat time', + ); + } + if (hasRegressingBeatTimes(nnTimesMs)) { + return const Metric.absent( + tier: Tier.high, + inputs_used: inputs, + note: 'regressing beat time', + ); + } final spanSec = tSec.last - tSec.first; if (spanSec <= 0) { return const Metric.absent( diff --git a/lib/src/onehz/sleep/circadian_np.dart b/lib/src/onehz/sleep/circadian_np.dart index c80c4789..d5365969 100644 --- a/lib/src/onehz/sleep/circadian_np.dart +++ b/lib/src/onehz/sleep/circadian_np.dart @@ -74,10 +74,19 @@ Metric circadianNonparametric( final lEp = l5Epochs ?? (epochsPerDay * 5 ~/ 24); final grand = mean(x)!; - var ssTot = 0.0; - for (final v in x) { - final d = v - grand; - ssTot += d * d; + // Welford (1962) centered variance; Witting (1990), van Someren (1999) + // indices retain their formulas. Circular profile sums keep input order. + var ssTot = 0.0, center = 0.0, count = 0; + for (final value in x) { + count++; + if (!value.isFinite) { + ssTot = double.nan; + center = x.take(count).fold(0, (a, b) => a + b) / count; + continue; + } + final delta = value - center; + center += delta / count; + ssTot += delta * (value - center); } if (ssTot == 0) { return const Metric.absent( diff --git a/lib/src/onehz/util.dart b/lib/src/onehz/util.dart index 634880dd..e9ccd4a8 100644 --- a/lib/src/onehz/util.dart +++ b/lib/src/onehz/util.dart @@ -345,6 +345,7 @@ LombScargle? lombScargle(List t, List y, List freqsHz) { } final span = tMax - tMin; if (span <= 0) return null; + final t0 = t[0]; // One-sided periodogram → density: 2·Δt, with Δt the MEAN sample interval. final psdScale = 2.0 * span / (n - 1); @@ -358,14 +359,14 @@ LombScargle? lombScargle(List t, List y, List freqsHz) { // τ: phase reference for time-shift invariance. var sin2 = 0.0, cos2 = 0.0; for (final ti in t) { - sin2 += math.sin(2 * w * ti); - cos2 += math.cos(2 * w * ti); + sin2 += math.sin(2 * w * (ti - t0)); + cos2 += math.cos(2 * w * (ti - t0)); } final tau = math.atan2(sin2, cos2) / (2 * w); var cNum = 0.0, cDen = 0.0, sNum = 0.0, sDen = 0.0; for (var i = 0; i < n; i++) { - final arg = w * (t[i] - tau); + final arg = w * ((t[i] - t0) - tau); final c = math.cos(arg); final s = math.sin(arg); cNum += yc[i] * c; @@ -499,3 +500,14 @@ List benjaminiHochberg(List ps) { } return out; } + +/// Whether supplied beat times move backward. Equal times remain allowed: +/// band timestamps can be quantized to whole seconds. +bool hasRegressingBeatTimes(List times) { + for (var i = 1; i < times.length; i++) { + // Non-finite timestamps have their own refusal note at the call site. + if (times[i].isFinite && times[i - 1].isFinite && + times[i] < times[i - 1]) return true; + } + return false; +} diff --git a/lib/src/onehz/workout/hr_zones.dart b/lib/src/onehz/workout/hr_zones.dart index 16c66e37..d2efec62 100644 --- a/lib/src/onehz/workout/hr_zones.dart +++ b/lib/src/onehz/workout/hr_zones.dart @@ -182,33 +182,29 @@ class HeartRateZones { HeartRateZoneSet zoneSet, ) { final sorted = [...hr]..sort((a, b) => a.tsMs.compareTo(b.tsMs)); - final zoneSeconds = List.filled(5, 0); - var below = 0.0; - - final tailSeconds = - sampleCadenceSeconds([for (final s in sorted) s.tsMs / 1000.0]); - if (tailSeconds == null) return null; - for (var i = 0; i < sorted.length; i++) { - final sample = sorted[i]; - if (!sample.valid) continue; - final durSeconds = i < sorted.length - 1 - ? _boundedGapSeconds(sorted[i + 1].tsMs - sample.tsMs, tailSeconds) - : tailSeconds; - final zone = zoneSet.zoneNumber(sample.hr); - if (zone >= 1) { - zoneSeconds[zone - 1] += durSeconds; - } else { - below += durSeconds; + // Karvonen (1957) supports reserve zones only. Duration credit is an + // internal ESTIMATE rule. Count equal gaps before summing their credit. + final cadence = sampleCadenceSeconds([for (final s in sorted) s.tsMs / 1000.0]); + if (cadence == null) return null; + final credits = List.generate(6, (_) => {}); + for (var i = 1; i < sorted.length; i++) { + final previous = sorted[i - 1]; + if (!previous.valid) continue; + final gap = (sorted[i].tsMs - previous.tsMs) / 1000.0; + final histogram = credits[zoneSet.zoneNumber(previous.hr)]; + histogram[gap] = (histogram[gap] ?? 0) + 1; + } + final seconds = List.filled(6, 0); + for (var zone = 0; zone < credits.length; zone++) { + for (final entry in credits[zone].entries) { + final gap = entry.key; + seconds[zone] += entry.value * (gap > 0 ? math.min(gap, cadence) : cadence); } } - return TimeInHeartRateZone(seconds: zoneSeconds, belowZone1: below); - } - - static double _boundedGapSeconds(double gapMs, double fallbackSeconds) { - final gapSeconds = gapMs / 1000.0; - return gapSeconds > 0 - ? math.min(gapSeconds, fallbackSeconds) - : fallbackSeconds; + if (sorted.isNotEmpty && sorted.last.valid) { + seconds[zoneSet.zoneNumber(sorted.last.hr)] += cadence; + } + return TimeInHeartRateZone(seconds: seconds.sublist(1), belowZone1: seconds[0]); } } diff --git a/test/fixtures/output_batch_parity.json b/test/fixtures/output_batch_parity.json new file mode 100644 index 00000000..9ddfdaff --- /dev/null +++ b/test/fixtures/output_batch_parity.json @@ -0,0 +1 @@ +{"source":"580e16733ea6cd0dec81066417d2e89972bdbb54","provenance":"Public batch outputs recorded on the batch-only output commit, before the framework refactor. Inputs are reproductionCases in tool/reproduce_outputs.dart.","cases":[{"rr":[201,0,1],"spectral":2191.9035532963126,"cosinor":[60.20705778058244,1.994548567920997,0.99021822497395,3.584442855968536e-94],"near_singular":false,"hrv":[3.615265234174367,21.52663924532641,1.0673833338741983,21.55214616158443,0.984955296232754],"nonfinite":false,"trimp":10.252064196950478,"ctl":[93.59303565371741,96.28901109956863,-2.695975445851218],"zones":[0.0,50.00499999999995,0.0,0.0,0.0],"circadian":[1.0,0.02086956521739154,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,17,0],"spectral":240.47476890807894,"cosinor":[60.20546901353539,1.9944760220371391,0.9906874461342637,3.6450131402730934e-95],"near_singular":false,"hrv":[14.544681533327767,21.511404623818784,0.2902542164251966,21.539302061762694,0.7620460195717869],"nonfinite":false,"trimp":13.023941743036463,"ctl":[94.97450463635269,101.20837984038745,-6.23387520403476],"zones":[0.0,61.00609999999995,0.0,0.0,0.0],"circadian":[1.0,0.2553191489361734,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,31,0],"spectral":4.814151167014313,"cosinor":[60.20670784840657,1.9932291959984847,0.9907102525961057,3.252232588931257e-95],"near_singular":false,"hrv":[39.550764408744726,21.5122083009717,0.40991536302833925,21.55029614999938,-0.7370544026973745],"nonfinite":false,"trimp":15.28358496194033,"ctl":[94.92528070209704,96.00888145442389,-1.083600752326845],"zones":[0.0,72.00719999999993,0.0,0.0,0.0],"circadian":[1.0,0.3312676056338072,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,32,0],"spectral":358.69230942776977,"cosinor":[60.206138700452115,1.9956478767563555,0.9905970400705738,5.7121975661113537e-95],"near_singular":false,"hrv":[42.09631660353519,21.499141419377167,0.29941934343651677,21.512944620556805,-0.9682323671330902],"nonfinite":false,"trimp":18.40304984916745,"ctl":[96.58887865497881,103.10914807456258,-6.520269419583769],"zones":[0.0,83.00829999999988,0.0,0.0,0.0],"circadian":[1.0,0.3688421052631628,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,33,0],"spectral":19.577975582737867,"cosinor":[60.206354166666685,1.99256198692952,0.9903700870164812,1.7315692772094542e-94],"near_singular":false,"hrv":[3.615265234174367,21.52663924532641,1.0673833338741983,21.55214616158443,0.984955296232754],"nonfinite":false,"trimp":20.354197211570078,"ctl":[95.91769269304092,96.02635036170024,-0.10865766865931903],"zones":[0.0,94.00939999999984,0.0,0.0,0.0],"circadian":[1.0,0.02086956521739154,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[200,33,0],"spectral":11606.502913565982,"cosinor":[60.2065027398695,1.994982877334022,0.9902341170988117,3.323420791543799e-94],"near_singular":false,"hrv":[14.544681533327767,21.511404623818784,0.2902542164251966,21.539302061762694,0.7620460195717869],"nonfinite":false,"trimp":23.84820212239254,"ctl":[97.91559020097185,105.11939669789578,-7.2038064969239315],"zones":[0.0,105.0104999999998,0.0,0.0,0.0],"circadian":[1.0,0.2553191489361734,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[201,0,1],"spectral":179.12281661844955,"cosinor":[60.20613433801807,1.993285411438357,0.9902158539353855,3.625067770914454e-94],"near_singular":false,"hrv":[39.54331943857053,21.508583812826043,0.27113546544007344,21.551316007749993,-0.7371891534838562],"nonfinite":false,"trimp":25.46704673479356,"ctl":[96.68685633705937,96.44160847368228,0.24524786337708804],"zones":[0.0,116.01159999999977,0.0,0.0,0.0],"circadian":[1.0,0.3312676056338072,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,17,0],"spectral":527.7181728812988,"cosinor":[60.20664182391473,1.993388211364034,0.990292034354844,2.520360883833827e-94],"near_singular":false,"hrv":[42.11090878105368,21.502568386886722,0.29644031390575065,21.528250446057537,-0.9685070047216077],"nonfinite":false,"trimp":28.750997708819092,"ctl":[98.08372430811947,101.79002204026516,-3.7062977321456856],"zones":[0.0,127.01269999999977,0.0,0.0,0.0],"circadian":[1.0,0.3688421052631628,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,31,0],"spectral":21.91751421630071,"cosinor":[60.206102619835434,1.9948632846663268,0.9904511511664523,1.168751936831454e-94],"near_singular":false,"hrv":[3.6114616155271952,21.50300002384371,0.916280807418444,21.477133240134368,0.9853180909453777],"nonfinite":false,"trimp":30.6254910760735,"ctl":[97.31583171818741,97.22062215666931,0.09520956151810367],"zones":[0.0,138.0137999999999,0.0,0.0,0.0],"circadian":[1.0,0.02086956521739154,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,32,0],"spectral":94.441233233861,"cosinor":[60.206607219531165,1.9918709847273792,0.9905159686967243,8.514781581412862e-95],"near_singular":false,"hrv":[14.550213162597531,21.50963654798735,0.2980266237242866,21.540300395940907,0.7627458277372915],"nonfinite":false,"trimp":33.68218231936987,"ctl":[98.15600470749584,99.32476832904142,-1.1687636215455797],"zones":[0.0,149.01489999999987,0.0,0.0,0.0],"circadian":[1.0,0.2553191489361734,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[199,33,0],"spectral":12.677409836990822,"cosinor":[60.20622592520504,1.9957534839219526,0.9905696980951562,6.537929637620507e-95],"near_singular":false,"hrv":[39.54331943857053,21.508583812826043,0.27113546544007344,21.551316007749993,-0.7371891534838562],"nonfinite":false,"trimp":35.83311260770479,"ctl":[97.86366868742407,98.31752817052501,-0.45385948310094193],"zones":[0.0,160.01599999999982,0.0,0.0,0.0],"circadian":[1.0,0.3312676056338072,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]},{"rr":[200,33,0],"spectral":35.03266947297398,"cosinor":[60.20646575449843,1.9913371893547758,0.9904828062401032,1.0015309381130758e-94],"near_singular":false,"hrv":[42.11090878105368,21.502568386886722,0.29644031390575065,21.528250446057537,-0.9685070047216077],"nonfinite":false,"trimp":38.64417150766156,"ctl":[98.18187926518608,97.60954398514504,0.5723352800410453],"zones":[0.0,171.01709999999983,0.0,0.0,0.0],"circadian":[1.0,0.3688421052631628,3.505000000000011,1.3600000000000079,14,0],"irregular_gap":[22.360679774997894,83.33333333333333]}]} diff --git a/test/fixtures/output_changes.json b/test/fixtures/output_changes.json new file mode 100644 index 00000000..f3c232f5 --- /dev/null +++ b/test/fixtures/output_changes.json @@ -0,0 +1,134 @@ +{ + "provenance": "Synthetic deterministic arithmetic cases, no personal data. Formula inputs are expanded by tool/reproduce_outputs.dart. Old values measured on up/main 3a7ca34; new values measured on f8e7d3c plus #98.", + "rr": { + "normal": 900.0, + "before": 100, + "run": 33, + "artifact": 250.0, + "after": 100 + }, + "spectral": { + "count": 96, + "origin": 1800000000.0, + "step": 0.8125, + "frequency": 0.137 + }, + "cosinor": { + "count": 96, + "step": 0.7 + }, + "near_singular": { + "count": 20000, + "origin": 3.0, + "step": 1e-09 + }, + "hrv": { + "count": 1800 + }, + "load": { + "count": 90 + }, + "zones": { + "count": 120 + }, + "circadian": { + "days": 3, + "epochs": 24 + }, + "old": { + "rr": [ + 199, + 34, + 0 + ], + "spectral": 32627.260306094682, + "cosinor": [ + 60.20591440999805, + 1.9937577719142656, + 0.9902626829318173, + 2.900200339561918e-94 + ], + "near_singular": true, + "hrv": [ + 3.615265234174367, + 21.52663924532639, + 1.067383333873884, + 21.55214616158443, + 0.9849552962327551 + ], + "nonfinite": "throws", + "trimp": 32.71197112260296, + "ctl": [ + 98.49620732954641, + 103.12182624417179, + -4.625618914625377 + ], + "zones": [ + 0.0, + 120.01199999999982, + 0.0, + 0.0, + 0.0 + ], + "circadian": [ + 1.0, + 0.331267605633807, + 3.505000000000011, + 1.3600000000000079, + 14, + 0 + ], + "irregular_gap": [ + 22.360679774997894, + 83.33333333333333 + ] + }, + "new": { + "rr": [ + 200, + 33, + 0 + ], + "spectral": 32627.260527635714, + "cosinor": [ + 60.205914409998044, + 1.9937577719142632, + 0.9902626829318162, + 2.9002003395774255e-94 + ], + "near_singular": false, + "hrv": [ + 3.615265234174367, + 21.52663924532641, + 1.0673833338741983, + 21.55214616158443, + 0.984955296232754 + ], + "nonfinite": false, + "trimp": 32.711971122602975, + "ctl": [ + 98.49620732954642, + 103.12182624417177, + -4.625618914625349 + ], + "zones": [ + 0.0, + 120.01199999999977, + 0.0, + 0.0, + 0.0 + ], + "circadian": [ + 1.0, + 0.3312676056338072, + 3.505000000000011, + 1.3600000000000079, + 14, + 0 + ], + "irregular_gap": [ + 22.360679774997894, + 83.33333333333333 + ] + } +} diff --git a/test/output_changes/batch_controls_test.dart b/test/output_changes/batch_controls_test.dart new file mode 100644 index 00000000..ebde99c9 --- /dev/null +++ b/test/output_changes/batch_controls_test.dart @@ -0,0 +1,42 @@ +import 'package:openstrap_analytics/onehz.dart'; +import 'package:test/test.dart'; + +void main() { + test('RR cap boundary preserves uncapped NN for short classified runs', () { + for (final n in [30, 31, 33]) { + final rr = [...List.filled(100, 900), + ...List.filled(n, 250), ...List.filled(100, 900)]; + final capped = correctRr(rr); + if (n <= 31) { + final old = correctRr(rr, maxArtifactRun: null); + expect(capped.nn, old.nn); + expect(capped.nnTimesMs, old.nnTimesMs); + expect(capped.classes, old.classes); + } else { + expect(capped.nn, List.filled(200, 900)); + expect(capped.droppedCount, 33); + } + } + }); + test('batch circular near-ties retain strict first-winner selection', () { + // Equal 10-epoch windows on an alternating profile tie at every start. + // A 5-epoch trough starts at zero, which includes three low epochs. + final profile = [for (var i = 0; i < 24; i++) i.isEven ? 1.0 : 2.0]; + final equal = circadianNonparametric(profile, 24).value!; + expect(equal.m10StartEpoch, 0); + expect(equal.l5StartEpoch, 0); + final near = [...profile]..[0] += 1e-12; + final result = circadianNonparametric(near, 24).value!; + expect(result.m10StartEpoch, 0); + expect(result.l5StartEpoch, 2); + }); + test('batch irregular screen excludes the gap pair already excluded by main', () { + final result = irregularBeatScreen([800, 900, 1000, 1100, 700, 850, 900, 1000], + nnTimesMs: [800, 1700, 2700, 3800, 14500, 15350, 16250, 17250], + minBeats: 4).value!; + // Six successive differences remain: 100,100,100,150,50,100. + // Five exceed 70 ms. The -400 ms difference across the gap is absent. + expect(result.pnnPct, 100 * 5 / 6); + expect(result.nBeats, 8); + }); +} diff --git a/test/output_changes/batch_parity_test.dart b/test/output_changes/batch_parity_test.dart new file mode 100644 index 00000000..d81a7b1e --- /dev/null +++ b/test/output_changes/batch_parity_test.dart @@ -0,0 +1,21 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + final config = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final golden = jsonDecode(File('test/fixtures/output_batch_parity.json').readAsStringSync()) as Map; + final cases = reproductionCases(config); + test('the full synthetic corpus is available', () { + expect(cases.length, 12); + expect(golden['cases'], hasLength(12)); + }); + for (var i = 0; i < 12; i++) { + test('split-2 public batch arithmetic survives refactoring, case $i', () { + // Frozen before aggregator adoption, not computed by the current path. + // Covers high-frequency jitter, time holes, cap boundaries and epochs. + expect(reproduce(cases[i]), golden['cases'][i], reason: 'case $i'); + }); + } +} diff --git a/test/output_changes/circadian_test.dart b/test/output_changes/circadian_test.dart new file mode 100644 index 00000000..480799e7 --- /dev/null +++ b/test/output_changes/circadian_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Compute circadian variance with centered moments: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['circadian']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/cosinor_test.dart b/test/output_changes/cosinor_test.dart new file mode 100644 index 00000000..6d0e6d45 --- /dev/null +++ b/test/output_changes/cosinor_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Refuse unidentifiable cosinor fits: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['cosinor', 'near_singular']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/ctl_test.dart b/test/output_changes/ctl_test.dart new file mode 100644 index 00000000..0c1d61ad --- /dev/null +++ b/test/output_changes/ctl_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Evaluate training-load decay with affine arithmetic: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['ctl']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/hrv_test.dart b/test/output_changes/hrv_test.dart new file mode 100644 index 00000000..05a5a9d9 --- /dev/null +++ b/test/output_changes/hrv_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Compute HRV moments without subtracting a large mean: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['hrv']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/nonfinite_test.dart b/test/output_changes/nonfinite_test.dart new file mode 100644 index 00000000..f26f6a1d --- /dev/null +++ b/test/output_changes/nonfinite_test.dart @@ -0,0 +1,31 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import 'package:openstrap_analytics/onehz.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + for (final bad in [double.nan, double.infinity, double.negativeInfinity]) { + test('all batch timestamp guards abstain on $bad', () { + final nn = [for (var i = 0; i < 1500; i++) 900.0 + 40 * (i % 7)]; + final ts = [for (var i = 0; i < nn.length; i++) i * 1100.0]; + ts[700] = bad; + for (final metric in [hrvTime(nn, nnTimesMs: ts), hrvFreq(nn, ts, artifactFraction: 0), + rsaRespRate(nn, ts, artifactFraction: 0), cvhrApneaScreen(nn, ts, artifactFraction: 0), nocturnalRmssd(nn, ts), + sleepSessionRmssdDetail(nn, ts, startSec: 1, endSec: 7200)]) { + expect(metric.present, isFalse); + expect(metric.note, 'non-finite beat time'); + } + }); + } + + test('Abstain when a beat has no finite time: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['nonfinite']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/review_inputs_test.dart b/test/output_changes/review_inputs_test.dart new file mode 100644 index 00000000..4e3a3dae --- /dev/null +++ b/test/output_changes/review_inputs_test.dart @@ -0,0 +1,158 @@ +import 'package:openstrap_analytics/onehz.dart'; +import 'package:test/test.dart'; + +void main() { + for (final bad in [double.nan, double.infinity, double.negativeInfinity]) { + test('HRV entry points abstain on non-finite NN interval $bad', () { + final nn = List.filled(1200, 900); + final ts = [for (var i = 0; i < nn.length; i++) 2000.0 + i * 900]; + nn[400] = bad; + for (final metric in [ + hrvTime(nn, nnTimesMs: ts), + nocturnalRmssd(nn, ts), + sleepSessionRmssdDetail(nn, ts, startSec: 1, endSec: 1200), + sleepSessionWindowedRmssd(nn, ts, startSec: 1, endSec: 1200), + ]) { + expect(metric.present, isFalse); + expect(metric.note, 'non-finite NN interval'); + } + expect(hrvTime(nn).present, isFalse); + }); + test('unknown training load $bad abstains in prime and later history', () { + for (final index in [0, 6, 7, 13]) { + final loads = List.filled(14, 50)..[index] = bad; + final metric = ctlAtlTsb(loads); + expect(metric.present, isFalse, reason: 'day $index'); + expect(metric.note, 'non-finite daily TRIMP'); + } + }); + } + test('a five-minute segment containing only NaNs abstains without throwing', () { + final nn = List.filled(1200, 900); + for (var i = 0; i < 334; i++) { nn[i] = double.nan; } + final ts = [for (var i = 0; i < nn.length; i++) i * 900.0]; + expect(hrvTime(nn, nnTimesMs: ts).present, isFalse); + }); + for (final period in [0.0, -24.0, double.nan, double.infinity, + double.negativeInfinity]) { + test('cosinor refuses invalid period $period', () { + final metric = cosinor( + [for (var i = 0; i < 48; i++) i.toDouble()], + [for (var i = 0; i < 48; i++) 60.0 + i % 3], + periodHours: period, + ); + expect(metric.present, isFalse); + expect(metric.note, 'invalid cosinor period'); + }); + } + test('a capped reset keeps source time but invalidates regressing HRV time', () { + // The RR clock exceeds elapsed sensor time before the cap. Re-anchoring + // the next independent segment must retain its actual sensor timestamp. + final rr = [...List.filled(100, 1500), + ...List.filled(33, 3000), ...List.filled(100, 1500)]; + final ts = [for (var i = 0; i < rr.length; i++) 1.7e12 + i * 1000]; + final result = correctRr(rr, rrTsMs: ts); + expect(result.nnTimesMs[99], 150000); + expect(result.nnTimesMs[100], 134500); + for (final metric in [ + hrvTime(result.nn, nnTimesMs: result.nnTimesMs), + nocturnalRmssd(result.nn, result.nnTimesMs), + sleepSessionRmssdDetail(result.nn, result.nnTimesMs, + startSec: 1, endSec: 400), + sleepSessionRmssdDetail(result.nn, result.nnTimesMs, + startSec: 1, endSec: 100), + ]) { + expect(metric.present, isFalse); + expect(metric.note, 'regressing beat time'); + } + }); + + for (final entry in Function(List, List)>{ + 'frequency HRV': (nn, ts) => hrvFreq(nn, ts, artifactFraction: 0), + 'RSA': (nn, ts) => rsaRespRate(nn, ts, artifactFraction: 0), + 'CVHR': (nn, ts) => cvhrApneaScreen(nn, ts, artifactFraction: 0), + }.entries) { + test('${entry.key} refuses a regressing spectral clock', () { + final nn = List.filled(1500, 900); + final ts = [for (var i = 0; i < nn.length; i++) i * 900.0]; + ts[700] = ts[699] - 500; + final metric = entry.value(nn, ts); + expect(metric.present, isFalse); + expect(metric.note, 'regressing beat time'); + }); + } + + for (final times in >[[], [800], [800, 1700, 2500]]) { + test('time HRV refuses a supplied ${times.length}-beat clock for two NN', () { + final metric = hrvTime([800, 900], nnTimesMs: times); + expect(metric.present, isFalse); + expect(metric.note, 'beat time length mismatch'); + // NN-only remains supported, as does a matching two-beat clock. + // One 100 ms difference has RMSSD exactly 100 ms. + expect(hrvTime([800, 900]).value!.rmssd, 100); + expect(hrvTime([800, 900], nnTimesMs: [800, 1700]).value!.rmssd, 100); + }); + } + + test('cosinor refuses overflow from a subnormal positive period', () { + final metric = cosinor( + List.generate(48, (i) => i.toDouble()), + List.generate(48, (i) => 60.0 + i % 3), + periodHours: 1e-320, + ); + expect(metric.present, isFalse); + expect(metric.note, 'non-finite cosinor design'); + }); + for (final bad in [double.nan, double.infinity, double.negativeInfinity]) { + for (final axis in ['time', 'response']) { + test('cosinor refuses non-finite $axis moment $bad', () { + final ts = List.generate(48, (i) => i.toDouble()); + final values = List.generate(48, (i) => 60.0 + i % 3); + (axis == 'time' ? ts : values)[12] = bad; + final metric = cosinor(ts, values); + expect(metric.present, isFalse); + expect(metric.note, 'non-finite cosinor design'); + }); + } + } + test('cosinor refuses non-finite moments from finite overflow', () { + final metric = cosinor( + List.generate(48, (i) => i.toDouble()), + List.generate(48, (i) => i.isEven ? 1e308 : -1e308), + ); + expect(metric.present, isFalse); + expect(metric.note, 'non-finite cosinor design'); + }); + + test('an entirely unknown load history abstains', () { + expect(ctlAtlTsb(List.filled(42, double.nan)).present, isFalse); + }); + + test('over-cap artifacts have an explicit gap label', () { + final rr = [ + ...List.filled(100, 900), + ...List.filled(33, 100), + ...List.filled(100, 900), + ]; + final result = correctRr(rr); + expect(result.classes.sublist(100, 133).map((c) => c.name), + everyElement('gap')); + expect(result.droppedCount, 33); + expect(result.nn, List.filled(200, 900)); + expect(result.nnSourceIndex, + [...List.generate(100, (i) => i), + ...List.generate(100, (i) => i + 133)]); + }); + + test('observed rest days still decay load', () { + // Seven 50-load days prime at 50. Seven observed rest days then apply + // exp(-1/tau) seven times, giving 50*exp(-7/42) and 50*exp(-7/7). + final metric = ctlAtlTsb([ + ...List.filled(7, 50), ...List.filled(7, 0), + ]); + expect(metric.value!.ctl, closeTo(42.3240862445307, 1e-12)); + expect(metric.value!.atl, closeTo(18.3939720585721, 1e-12)); + expect(ctlAtlTsb(List.filled(14, 50)).value!.ctl, + closeTo(50, 1e-12)); + }); +} diff --git a/test/output_changes/rr_test.dart b/test/output_changes/rr_test.dart new file mode 100644 index 00000000..1421dbd6 --- /dev/null +++ b/test/output_changes/rr_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Leave long RR artifact runs as gaps: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['rr']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/spectral_test.dart b/test/output_changes/spectral_test.dart new file mode 100644 index 00000000..01e2fce5 --- /dev/null +++ b/test/output_changes/spectral_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Reduce spectral error on epoch timestamps: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['spectral']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/trimp_test.dart b/test/output_changes/trimp_test.dart new file mode 100644 index 00000000..afaf57d5 --- /dev/null +++ b/test/output_changes/trimp_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Retain small TRIMP contributions in the total: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['trimp']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/test/output_changes/zones_test.dart b/test/output_changes/zones_test.dart new file mode 100644 index 00000000..c9b772dc --- /dev/null +++ b/test/output_changes/zones_test.dart @@ -0,0 +1,16 @@ +import 'dart:convert'; +import 'dart:io'; +import 'package:test/test.dart'; +import '../../tool/reproduce_outputs.dart'; + +void main() { + test('Count repeated zone intervals before summing: unrounded synthetic batch regression', () { + // Old/new values were captured independently on 3a7ca34 and f8e7d3c. + // The fixture gives the input formulas and provenance, never personal data. + final fixture = jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + final actual = reproduce(fixture); + for (final key in ['zones']) { + expect(actual[key], fixture['new'][key], reason: key); + } + }); +} diff --git a/tool/reproduce_outputs.dart b/tool/reproduce_outputs.dart new file mode 100644 index 00000000..93fdc2b4 --- /dev/null +++ b/tool/reproduce_outputs.dart @@ -0,0 +1,81 @@ +// Synthetic reproduction cases. Run with --json to capture unrounded values. +// Uses only public batch functions, so the same file runs on up/main and #99. +import 'dart:convert'; +import 'dart:io'; +import 'dart:math' as math; +import 'package:openstrap_analytics/onehz.dart'; +import 'package:openstrap_analytics/src/onehz/util.dart' as u; + +// Deterministic differential cases, shared by the tool and parity tests. +List> reproductionCases(Map base) => [ + for (var i = 0; i < 12; i++) (() { + final f = jsonDecode(jsonEncode(base)) as Map; + f['hrv']['count'] = i < 6 ? 1800 : 5400; + f['hrv']['frequency'] = [.17, .7, 2.4, 2.9][i % 4]; + f['hrv']['gapped'] = i.isOdd; + f['spectral']['origin'] = i.isEven ? 1800000000.0 : 0.0; + f['spectral']['frequency'] = .037 + i * .019; + f['rr']['run'] = [1, 16, 30, 31, 32, 33][i % 6]; + f['cosinor']['step'] = .31 + i * .11; + f['load']['count'] = 30 + i * 7; + f['zones']['count'] = 50 + i * 11; + f['circadian']['days'] = 1 + i % 4; + return f; + })() +]; + +Map reproduce(Map f) { + final out = {}; + final r = f['rr'] as Map; + final rr = [for (var i=0;i<(r['before'] as int);i++) (r['normal'] as num).toDouble(), + for (var i=0;i<(r['run'] as int);i++) (r['artifact'] as num).toDouble(), + for (var i=0;i<(r['after'] as int);i++) (r['normal'] as num).toDouble()]; + final clean = correctRr(rr); + out['rr'] = [clean.nn.length, clean.droppedCount, clean.correctedCount]; + final s=f['spectral'] as Map; + final t=[for(var i=0;i<(s['count'] as int);i++) (s['origin'] as num).toDouble()+i*(s['step'] as num)]; + final y=[for(var i=0;i; + final ct=[for(var i=0;i<(c['count'] as int);i++) i*(c['step'] as num).toDouble()]; + final cy=[for(var i=0;i; + final nt=[for(var i=0;i<(n['count'] as int);i++) (n['origin'] as num).toDouble()+i*(n['step'] as num)]; + out['near_singular']=cosinor(nt,[for(var i=0;i0 && i%400==0 ? 10000 : 0)]; + final hv=hrvTime(nn,nnTimesMs:times).value!; + out['hrv']=[hv.rmssd,hv.sdnn,hv.sdann,hv.sdnnIndex,hv.diffAcf1]; + try { + final bad=[...times]..[times.length-1]=double.infinity; + out['nonfinite']=hrvTime(nn,nnTimesMs:bad).present; + } catch(e) { out['nonfinite']='throws'; } + final load=[for(var i=0;i<(f['load']['count'] as int);i++) 80.1+(i%13)*3.17]; + out['trimp']=banisterTrimp(load,restingHr:60,maxHr:190,sex:Sex.male).value; + final lv=ctlAtlTsb(load).value!; + out['ctl']=[lv.ctl,lv.atl,lv.tsb]; + final z=[for(var i=0;i<(f['zones']['count'] as int);i++) HrSample(i*1000.1,120.0)]; + out['zones']=HeartRateZones.timeInZone(z,HeartRateZones.zonesFromMaxHr(190))!.seconds; + final x=[for(var i=0;i<(f['circadian']['days'] as int)*(f['circadian']['epochs'] as int);i++) 1.1+(i%24)*.13+(i%3)*1e-14]; + final ci=circadianNonparametric(x,f['circadian']['epochs'] as int).value!; + out['circadian']=[ci.interdailyStability,ci.intradailyVariability,ci.m10,ci.l5,ci.m10StartEpoch,ci.l5StartEpoch]; + final ir=irregularBeatScreen([800,900,1000,1100,700,850,900,1000].map((v)=>v.toDouble()).toList(), + nnTimesMs:[800,1700,2700,3800,14500,15350,16250,17250],minBeats:4).value!; + out['irregular_gap']=[ir.sd1,ir.pnnPct]; + return out; +} +void main(List args) { + final f=jsonDecode(File('test/fixtures/output_changes.json').readAsStringSync()) as Map; + if(args.contains('--corpus-json')) { print(jsonEncode([for(final c in reproductionCases(f)) reproduce(c)])); return; } + final values=reproduce(f); + if(args.contains('--json')) { print(jsonEncode(values)); return; } + print('| Change | up/main 3a7ca34 | Current batch |'); + print('|---|---|---|'); + final old=f['old'] as Map; + for(final entry in values.entries) { + print('| ${entry.key} | ${old[entry.key]} | ${entry.value} |'); + } +}