diff --git a/lib/src/onehz/clinical/readiness_lnrmssd.dart b/lib/src/onehz/clinical/readiness_lnrmssd.dart index 82243c6..296dcf5 100644 --- a/lib/src/onehz/clinical/readiness_lnrmssd.dart +++ b/lib/src/onehz/clinical/readiness_lnrmssd.dart @@ -15,14 +15,17 @@ import '../util.dart'; class ReadinessLnRmssd { final double lnRmssdToday; final double rolling7Mean; + /// The window's MEDIAN — the centre z/SWC/band are taken against. + final double rolling7Median; final double cvPct; // CV of lnRMSSD over the window - final double? z; // (today - rollingMean)/rollingSD - final double? swc; // smallest worthwhile change (0.5×SD, Plews-style) + final double? z; // (today - rollingMedian)/robustSD + final double? swc; // smallest worthwhile change (0.5×robust SD, Plews-style) final String band; // 'suppressed' | 'normal' | 'elevated' final bool saturationFlag; const ReadinessLnRmssd({ required this.lnRmssdToday, required this.rolling7Mean, + required this.rolling7Median, required this.cvPct, required this.z, required this.swc, @@ -32,6 +35,7 @@ class ReadinessLnRmssd { Map toJson() => { 'ln_rmssd_today': round6(lnRmssdToday), 'rolling7_mean': round6(rolling7Mean), + 'rolling7_median': round6(rolling7Median), 'cv_pct': round6(cvPct), if (z != null) 'z': round6(z!), if (swc != null) 'swc': round6(swc!), @@ -117,14 +121,25 @@ Metric readinessLnRmssd( ); } final cv = (sd / m).abs() * 100; - final z = sd > 0 ? (today - m) / sd : null; + // ROBUST CENTRE AND SPREAD for the decision band. One illness night in a + // seven-night window both drags the mean down and inflates the SD, so a + // genuinely low night reads "normal": on a real window holding a 35.5 ms + // night, mean 4.465 / SD 0.405 put a 4.478 night (88 ms) inside the band, + // while the other six nights sat at 4.40–4.79. Median + MAD (×1.4826) are + // what the window looks like without that one night. MAD collapses only when + // more than half the window is identical; the SD stands in then. The CV + // above stays the classic SD/mean (Plews' published CV). + final centre = median(priorWindow)!; + final robustSd = mad(priorWindow)!; + final spread = robustSd > 0 ? robustSd : sd; + final z = spread > 0 ? (today - centre) / spread : null; // Plews SWC ≈ 0.5 × within-window SD (a small worthwhile change in lnRMSSD). - final swc = 0.5 * sd; + final swc = 0.5 * spread; final String band; - if (today < m - swc) { + if (today < centre - swc) { band = 'suppressed'; - } else if (today > m + swc) { + } else if (today > centre + swc) { band = 'elevated'; } else { band = 'normal'; @@ -140,6 +155,7 @@ Metric readinessLnRmssd( value: ReadinessLnRmssd( lnRmssdToday: today, rolling7Mean: m, + rolling7Median: centre, cvPct: cv, z: z, swc: swc, diff --git a/lib/src/onehz/wellness/readiness_composite.dart b/lib/src/onehz/wellness/readiness_composite.dart index ba71ce9..ba90a64 100644 --- a/lib/src/onehz/wellness/readiness_composite.dart +++ b/lib/src/onehz/wellness/readiness_composite.dart @@ -26,8 +26,12 @@ // 3. Weighted sum with disclosed weights HRV > RHR > RR > temp; weights are // RENORMALIZED over the inputs actually present (missing inputs are // dropped, never zero-imputed). +// HRV and RHR form ONE autonomic input (mean of their oriented z's, weight +// 0.70) — see [readinessAutonomicLabels]. // 4. The composite z is mapped to a 0..100 score via a logistic so typical -// days land near 50. +// days land near 50, rescaled by the user's own composite-z spread once +// 14 prior nights exist, and re-centred on their median (capped ±0.5) +// — see [calibratedReadinessScore]. // 5. ALWAYS attach the per-input contribution breakdown |w_i·z_i| ranked. // // HONESTY: glass-box index (weights disclosed); "—" when no inputs present; @@ -127,13 +131,139 @@ ReadinessInput tempInput( class Readiness { final double score; // 0..100 glass-box readiness final double compositeZ; // weighted, sign-oriented composite z - const Readiness(this.score, this.compositeZ); + + /// The personal spread σ̂ [score] was mapped with; null while calibrating + /// (fewer than [readinessCalibrationMinNights] prior composite z's). + final double? calibrationSigma; + + /// Prior composite z's the calibration had to work with. + final int calibrationNights; + + /// Median of the prior composite z's (raw) and the centre actually removed + /// (capped to ±[readinessCentreCap]); null while calibrating. + final double? calibrationCentreRaw; + final double? calibrationCentre; + const Readiness(this.score, this.compositeZ, + {this.calibrationSigma, + this.calibrationNights = 0, + this.calibrationCentreRaw, + this.calibrationCentre}); Map toJson() => { 'score': round6(score), 'composite_z': round6(compositeZ), + 'calibration': { + 'status': calibrationSigma == null ? 'calibrating' : 'calibrated', + 'nights': calibrationNights, + if (calibrationSigma != null) 'sigma': round6(calibrationSigma!), + if (calibrationCentreRaw != null) + 'centre_raw': round6(calibrationCentreRaw!), + if (calibrationCentre != null) 'centre': round6(calibrationCentre!), + }, }; } +/// The two inputs that read the SAME autonomic state. Lower HRV and higher RHR +/// move together (oriented-z correlation 0.82 over 52 real nights), so scoring +/// them as two independent inputs counted one signal twice — 70% of the weight +/// rode on it. They now form ONE autonomic input: the mean of their oriented +/// z's, carrying their combined weight. +const Set readinessAutonomicLabels = {'HRV', 'RHR'}; + +/// Combined weight of the autonomic input (HRV 0.40 + RHR 0.30). +const double readinessAutonomicWeight = 0.70; + +/// Combine sign-oriented z's into the composite. [used] is every input that +/// produced a z; [autonomicWeight] is the weight the HRV/RHR group carries +/// whenever at least one of them is present, split evenly across the members +/// present (so their contributions are the mean of their z's × the weight). +/// Every other input keeps its own weight. Returns the composite (weights +/// renormalised over what is present), the per-label contribution (summing to +/// the composite), and the weight that was present. +({double z, Map contributions, double weightSum}) + combineReadinessZ( + List<({String label, double weight, double orientedZ})> used, { + double autonomicWeight = readinessAutonomicWeight, +}) { + final nAuto = + used.where((u) => readinessAutonomicLabels.contains(u.label)).length; + var weightSum = nAuto > 0 ? autonomicWeight : 0.0; + final raw = {}; + for (final u in used) { + final auto = readinessAutonomicLabels.contains(u.label); + if (!auto) weightSum += u.weight; + raw[u.label] = (auto ? autonomicWeight / nAuto : u.weight) * u.orientedZ; + } + if (weightSum <= 0) return (z: 0.0, contributions: const {}, weightSum: 0.0); + final contributions = { + for (final e in raw.entries) e.key: e.value / weightSum, + }; + return ( + z: contributions.values.fold(0.0, (a, b) => a + b), + contributions: contributions, + weightSum: weightSum, + ); +} + +/// Nights of prior composite z before the score mapping is personalised. +const int readinessCalibrationMinNights = 14; + +/// Floor on σ̂: a run of near-identical nights must not blow a small z up. +const double readinessCalibrationSigmaFloor = 0.3; + +/// The composite-z spread the Home band cut-offs were set on (26/37/61 are the +/// 5th/20th/75th percentiles of logistic(N(0, 0.65))). +const double readinessDesignSigma = 0.65; + +/// Cap on the centre removed from z. A user whose nights sit persistently +/// below their own per-input baselines (one real history: composite-z median +/// −0.35, HRV under its trailing median on 77% of nights) otherwise reads +/// "low" every day. Re-centring on the trailing median fixes that, but an +/// uncapped centre would also absorb a real downturn within ~2 weeks. The cap +/// lets a step-down beyond 0.5 z keep showing: the part of it past the cap is +/// never re-centred away, and the rest only as the median catches up. +/// ponytail: one population constant; personalise if real step-downs say so. +const double readinessCentreCap = 0.5; + +/// Map a composite z to 0..100. +/// +/// The band cut-offs assume composite z has SD ≈ 0.65. Real users do not share +/// one spread: on one real 52-night history it was 1.26, which put 19% of +/// nights in "Rest today" against a designed 5%. So once [zHistory] (this +/// user's PRIOR composite z's) holds ≥ [readinessCalibrationMinNights] nights, +/// z is rescaled by the user's own robust spread σ̂ = MAD × 1.4826 (floored at +/// [readinessCalibrationSigmaFloor]): +/// and re-centred on the history's median ĉ, capped to ±[readinessCentreCap]: +/// score = 100 / (1 + exp(−0.65·(z − ĉ) / σ̂)) +/// Fewer nights → the uncalibrated logistic(z), flagged as calibrating. +({ + double score, + double? sigma, + int nights, + double? centreRaw, + double? centre +}) calibratedReadinessScore(double z, List zHistory) { + final n = zHistory.length; + if (n < readinessCalibrationMinNights) { + return ( + score: 100 / (1 + math.exp(-z)), + sigma: null, + nights: n, + centreRaw: null, + centre: null, + ); + } + final sigma = math.max(readinessCalibrationSigmaFloor, mad(zHistory)!); + final raw = median(zHistory)!; + final centre = raw.clamp(-readinessCentreCap, readinessCentreCap).toDouble(); + return ( + score: 100 / (1 + math.exp(-readinessDesignSigma * (z - centre) / sigma)), + sigma: sigma, + nights: n, + centreRaw: raw, + centre: centre, + ); +} + /// Compute the honest readiness composite. /// /// Each present input with a usable robust baseline contributes a sign-oriented @@ -169,15 +299,15 @@ Metric readinessComposite( int minBaseline = readinessCompositeMinBaseline, int minInputs = readinessCompositeMinInputs, double minWeightSum = readinessCompositeMinWeight, + List compositeZHistory = const [], }) { final used = []; - final drivers = []; + final parts = <({String label, double weight, double orientedZ})>[]; + final details = {}; final refusals = [ for (final inp in inputs) if (inp.refusal != null) inp.refusal! ]; - var weightSum = 0.0; - var weightedZ = 0.0; // Track the best-covered input that has a value but a too-short baseline, so // we can emit a machine-readable need_baseline note when nothing computes. var anyValuePresent = false; @@ -230,9 +360,7 @@ Metric readinessComposite( if (zr == null) continue; final oriented = inp.goodSign * zr; // + = good for readiness used.add(inp.label); - weightSum += inp.weight; - weightedZ += inp.weight * oriented; - // Driver contribution is the signed weighted z (renormalized later). + parts.add((label: inp.label, weight: inp.weight, orientedZ: oriented)); // GLASS-BOX: the disclosed method must be the method ACTUALLY used — on a // quantized baseline the MAD collapses and the mean/SD fallback above // produced this z, so saying "robust-z" there would misstate how the @@ -240,9 +368,13 @@ Metric readinessComposite( final method = rz != null ? 'robust-z (median+MAD)' : 'z (mean+SD fallback — MAD=0 on a quantized baseline)'; - drivers.add(Driver(inp.label, inp.weight * oriented, - detail: 'oriented $method=${round6(oriented)}')); + details[inp.label] = 'oriented $method=${round6(oriented)}'; } + final combined = combineReadinessZ(parts, autonomicWeight: [ + for (final inp in inputs) + if (readinessAutonomicLabels.contains(inp.label)) inp.weight + ].fold(0.0, (a, b) => a + b)); + final weightSum = combined.weightSum; final suffix = refusals.isEmpty ? '' : ' Refused: ${refusals.join('; ')}.'; if (used.length < minInputs || weightSum < minWeightSum) { // If inputs HAD values but their baselines were too short, say so in the @@ -272,33 +404,39 @@ Metric readinessComposite( '$suffix', ); } - // Renormalize weights over present inputs. - final composite = weightedZ / weightSum; - // Renormalize driver contributions by the same factor so they sum to the + // Weights renormalised over present inputs; contributions sum to the // composite z (glass-box: contributions are definitional within the formula). + final composite = combined.z; final normDrivers = [ - for (final d in drivers) - Driver(d.label, roundTo(d.contribution / weightSum, 6), detail: d.detail) + for (final e in combined.contributions.entries) + Driver(e.key, roundTo(e.value, 6), detail: details[e.key]) ]; // Rank by |contribution| (the deterministic-narrative driver ordering). normDrivers .sort((a, b) => b.contribution.abs().compareTo(a.contribution.abs())); - // Map composite z -> 0..100 via logistic; ~50 at z=0, scale so ±2 z ~ 12/88. - final score = 100 / (1 + math.exp(-composite)); + final cal = calibratedReadinessScore(composite, compositeZHistory); // Confidence scales with how many inputs were available (more = better). final conf = (0.3 + 0.15 * used.length).clamp(0.3, 0.9); return Metric( - value: Readiness(score, composite), + value: Readiness(cal.score, composite, + calibrationSigma: cal.sigma, + calibrationNights: cal.nights, + calibrationCentreRaw: cal.centreRaw, + calibrationCentre: cal.centre), confidence: conf, tier: Tier.estimate, inputs_used: used, drivers: normDrivers, note: - 'GLASS-BOX readiness: disclosed weights HRV>RHR>RR>temp, renormalized ' - 'over present inputs. Drivers are definitional within the formula ' - '(correction, not inferred cause).$suffix', + 'GLASS-BOX readiness: disclosed weights autonomic (HRV+RHR, ' + 'mean of their z) 0.70 > RR 0.20 > temp 0.10, renormalized over present ' + 'inputs. Drivers are definitional within the formula (correction, not ' + 'inferred cause). ' + '${cal.sigma == null ? 'calibrating:have=${cal.nights},need=$readinessCalibrationMinNights' : 'calibrated:sigma=${round6(cal.sigma!)},n=${cal.nights},' + 'centre=${round6(cal.centre!)},centre_raw=${round6(cal.centreRaw!)}'}.' + '$suffix', ); } diff --git a/test/onehz/readiness_calibration_test.dart b/test/onehz/readiness_calibration_test.dart new file mode 100644 index 0000000..074c988 --- /dev/null +++ b/test/onehz/readiness_calibration_test.dart @@ -0,0 +1,181 @@ +import 'dart:math' as math; + +import 'package:openstrap_analytics/onehz.dart'; +import 'package:test/test.dart'; + +List around(double c, double step) => + [for (var i = 0; i < 14; i++) c + (i.isEven ? step : -step)]; + +void main() { + group('autonomic input (HRV + RHR count once)', () { + test('HRV and RHR share 0.70 as the mean of their oriented z', () { + final m = readinessComposite([ + hrvInput(4.0, around(4.5, 0.1)), // well below baseline + rhrInput(60.0, around(52.0, 2.0)), // well above baseline + respInput(14.0, around(14.0, 0.5)), // at baseline + tempInput(2000.0, around(2000.0, 2.0), settledFraction: 0.97), + ]); + final c = {for (final d in m.drivers!) d.label: d.contribution}; + double oz(String l) => double.parse( + m.drivers!.firstWhere((d) => d.label == l).detail!.split('=').last); + final zH = oz('HRV'), zR = oz('RHR'); + // RR and temp sit on their medians (z 0), so the composite is the + // autonomic input alone: 0.70 × mean(zH, zR) over a total weight of 1.0. + expect(m.value!.compositeZ, closeTo(0.70 * (zH + zR) / 2, 1e-5)); + expect(c['HRV'], closeTo(0.35 * zH, 1e-5)); + expect(c['RHR'], closeTo(0.35 * zR, 1e-5)); + }); + + test('one autonomic member alone carries the whole 0.70', () { + final r = combineReadinessZ([ + (label: 'HRV', weight: 0.40, orientedZ: -1.0), + (label: 'RR', weight: 0.20, orientedZ: 1.0), + ]); + expect(r.weightSum, closeTo(0.90, 1e-12)); + expect(r.z, closeTo((0.70 * -1.0 + 0.20 * 1.0) / 0.90, 1e-12)); + }); + + test('a correlated HRV/RHR dip no longer outweighs everything else', () { + // Old weights: (0.4·-2 + 0.3·-2 + 0.2·2 + 0.1·2)/1.0 = -0.8. + // Autonomic group: (0.7·-2 + 0.2·2 + 0.1·2)/1.0 = -0.8 too — the merge + // changes nothing when the two agree; it matters when they disagree. + final agree = combineReadinessZ([ + (label: 'HRV', weight: 0.4, orientedZ: -2.0), + (label: 'RHR', weight: 0.3, orientedZ: -2.0), + (label: 'RR', weight: 0.2, orientedZ: 2.0), + (label: 'temp', weight: 0.1, orientedZ: 2.0), + ]); + expect(agree.z, closeTo(-0.8, 1e-12)); + // Disagreeing: HRV -2, RHR +2 → old -0.2 (HRV wins by weight), now 0. + final disagree = combineReadinessZ([ + (label: 'HRV', weight: 0.4, orientedZ: -2.0), + (label: 'RHR', weight: 0.3, orientedZ: 2.0), + ]); + expect(disagree.z, closeTo(0.0, 1e-12)); + }); + }); + + group('personal spread calibration', () { + test('< 14 prior nights: uncalibrated logistic, marked calibrating', () { + final r = calibratedReadinessScore(-1.0, List.filled(13, 0.5)); + expect(r.sigma, isNull); + expect(r.score, closeTo(100 / (1 + math.exp(1.0)), 1e-9)); + final m = readinessComposite([ + hrvInput(4.5, around(4.5, 0.1)), + rhrInput(52.0, around(52.0, 2.0)), + ], compositeZHistory: List.filled(5, 0.0)); + expect(m.value!.toJson()['calibration'], + {'status': 'calibrating', 'nights': 5}); + expect(m.note, contains('calibrating:have=5,need=14')); + }); + + test('≥ 14 nights: z rescaled by 0.65/σ̂, σ̂ and n disclosed', () { + final hist = [for (var i = 0; i < 20; i++) (i.isEven ? 1.0 : -1.0)]; + final r = calibratedReadinessScore(-1.0, hist); + final sigma = 1.0 * 1.4826; // MAD 1 × 1.4826 + expect(r.sigma, closeTo(sigma, 1e-9)); + expect(r.nights, 20); + expect(r.score, closeTo(100 / (1 + math.exp(0.65 / sigma)), 1e-9)); + }); + + test('σ̂ is floored so a flat history cannot blow a small z up', () { + final r = calibratedReadinessScore(-0.3, List.filled(20, 0.0)); + expect(r.sigma, readinessCalibrationSigmaFloor); + expect(r.score, closeTo(100 / (1 + math.exp(0.65)), 1e-9)); + }); + + test('high-variance user: "Rest today" share falls to ~the designed 5%', + () { + // Composite z ~ N(0, 1.26) — the spread measured on a real 52-night + // history (the cut-offs assume 0.65). Deterministic seed. + final rnd = math.Random(7); + double gauss() => + math.sqrt(-2 * math.log(1 - rnd.nextDouble())) * + math.cos(2 * math.pi * rnd.nextDouble()); + final zs = [for (var i = 0; i < 600; i++) 1.26 * gauss()]; + var restOld = 0, restNew = 0, scored = 0; + for (var i = 28; i < zs.length; i++) { + final hist = zs.sublist(i - 28, i); + scored++; + if (100 / (1 + math.exp(-zs[i])) < 26) restOld++; + if (calibratedReadinessScore(zs[i], hist).score < 26) restNew++; + } + expect(restOld / scored, greaterThan(0.12)); // the bug: ~15-20% + expect(restNew / scored, inInclusiveRange(0.02, 0.09)); + }); + }); + + group('capped re-centring', () { + final rnd = math.Random(11); + double gauss() => + math.sqrt(-2 * math.log(1 - rnd.nextDouble())) * + math.cos(2 * math.pi * rnd.nextDouble()); + + test('a history centred at -0.35 scores a median of ~50', () { + final zs = [for (var i = 0; i < 400; i++) -0.35 + 0.9 * gauss()]; + final scores = [ + for (var i = 28; i < zs.length; i++) + calibratedReadinessScore(zs[i], zs.sublist(i - 28, i)).score + ]..sort(); + expect(scores[scores.length ~/ 2], closeTo(50, 4)); + }); + + test('a sustained -1.0 step-down reads low at first, converges slowly', () { + // 28 nights around 0, then every night 1.0 lower. + final zs = [ + for (var i = 0; i < 28; i++) 0.6 * math.sin(i * 2.4), // spread, median ~0 + for (var i = 0; i < 40; i++) -1.0 + 0.6 * math.sin((i + 28) * 2.4), + ]; + final s = [ + for (var i = 28; i < zs.length; i++) + calibratedReadinessScore(zs[i], zs.sublist(i - 28, i)) + ]; + double meanScore(Iterable<({double score, double? sigma, int nights, + double? centreRaw, double? centre})> xs) => + xs.map((r) => r.score).reduce((a, b) => a + b) / xs.length; + // First week: the centre has barely moved, so the step shows in full — + // below the "Take it easy" cut-off (37) on average. + final first = meanScore(s.take(7)); + expect(first, lessThan(37)); + // Converges slowly: the last week reads higher, but the cap keeps the + // half of the step past -0.5 from ever being re-centred away. + final last = meanScore(s.skip(s.length - 7)); + expect(last, greaterThan(first)); + expect(last, lessThan(50)); + expect(s.last.centreRaw!, lessThan(-0.5)); + expect(s.last.centre, -readinessCentreCap); + }); + + test('cap boundary: centre is the median inside ±0.5, clamped outside', () { + final at = calibratedReadinessScore(0, List.filled(20, 0.5)); + expect(at.centreRaw, 0.5); + expect(at.centre, 0.5); + final past = calibratedReadinessScore(0, List.filled(20, -0.8)); + expect(past.centreRaw, -0.8); + expect(past.centre, -0.5); + // Flat history → σ̂ floor 0.3: (0 − (−0.5)) · 0.65 / 0.3. + expect(past.score, closeTo(100 / (1 + math.exp(-0.65 * 0.5 / 0.3)), 1e-9)); + final m = readinessComposite([ + hrvInput(4.5, around(4.5, 0.1)), + rhrInput(52.0, around(52.0, 2.0)), + ], compositeZHistory: List.filled(20, -0.8)); + final cal = m.value!.toJson()['calibration'] as Map; + expect(cal['centre_raw'], -0.8); + expect(cal['centre'], -0.5); + expect(m.note, contains('centre=-0.5,centre_raw=-0.8')); + }); + }); + + group('readinessLnRmssd robust centre', () { + test('one illness night cannot make a low night "normal"', () { + // The real window: six nights 4.40–4.79 plus one 35.5 ms (3.568) night, + // then an 88 ms (4.478) night. + final hist = [4.56, 3.568, 4.623, 4.397, 4.662, 4.794, 4.653, 4.478]; + final m = readinessLnRmssd(hist); + expect(m.value!.rolling7Mean, closeTo(4.4653, 1e-3)); + expect(m.value!.rolling7Median, closeTo(4.623, 1e-9)); + expect(m.value!.band, 'suppressed'); // mean/SD said 'normal' + expect(m.value!.z!, lessThan(-1)); + }); + }); +}