diff --git a/ALGORITHMS.md b/ALGORITHMS.md index fd8f3e0..ec584b5 100644 --- a/ALGORITHMS.md +++ b/ALGORITHMS.md @@ -51,6 +51,7 @@ Every metric in this package is a **published, peer-reviewed algorithm** compute | 24 | HRV stability | `calcHrvStability` | **CV** = SD/mean of nightly RMSSD over a window | HRV reliability/CV (Plews/Flatt) | RMSSD series → CV % (HIGH) | | 25 | Irregular-beat screen | `calcIrregular` | **Poincaré** SD1/SD2 + ectopic-rejection fraction + pNN50 | Brennan et al., *IEEE TBME* 2001; AF-screening literature | RR stream → flag + SD1/SD2 (ESTIMATE) | | 26 | Readiness (composite) | `calcReadinessIndex` | Transparent weighted blend: recovery·0.5 + sleep·0.2 + dip·0.15 + calm·0.15 | Composite (documented weights; abstains w/o HRV) | recovery/sleep/dip/arousal → 0–100 (ESTIMATE) | +| 27 | Steps | `calcSteps` / `pedometer` (`steps.ts`); runner `runStepsImu` (backend) | **AN-2554** wrist pedometer: dynamic-threshold peak pairs + 8-step confirm, ×gain | Analog Devices AN-2554 (2023); Zhao, *Analog Dialogue* 2010 | wrist IMU accel (R10+0x33, ~100 Hz) → steps (ESTIMATE) | | – | Max HR helper | `resolveMaxHr` | measured session max → 220−age → observed → 190 | Fox 1971 (age fallback) | minutes/baseline/age → HRmax + source | > `buildCoach` (deterministic plan) and `buildNotifications` (nudges) are rule engines over @@ -234,6 +235,31 @@ the old heuristic readiness was retired.) > from the session's minute window, not in this package — they're descriptive aggregates of > the same inputs, not new physiology. +### 27. Steps — `calcSteps` / `pedometer` (`steps.ts`); runner `runStepsImu` (backend) +**Where:** the pure AN-2554 math (`pedometer`, `calcSteps`, calibration gain) lives here in the +analytics package like every other metric; the backend `steps_imu.ts` is a thin **runner** that +re-decodes the IMU frames from R2 (`frameAccel` in `decode.ts`), dedups + groups them per minute, +and feeds the signals to `calcSteps` — mirroring the HRV/resp runners. +**Why:** the WHOOP 4.0 exposes **no step counter** over Bluetooth (the official app computes +steps in the cloud from raw accelerometer + ML; even the most complete community client falls +back to phone steps on 4.0). So we derive them ourselves from the wrist accelerometer with +Analog Devices' **AN-2554** time-domain pedometer — ~97% accurate on steady wrist gait. +**How:** the high-rate IMU (~100 Hz) arrives on two live channels — **R10** (pkt 0x2B, 100 +accel samples/axis) and the **0x33 IMU stream** (10 accel + 10 gyro samples/frame: X/Y/Z +contiguous from byte 24, frame-index @14, scale 1/4096 g) — re-decoded from the raw frames in +R2, **deduped by (ts, frame-index)** (upload windows overlap), assembled into a contiguous +per-minute signal. Per minute: `sum(|x|+|y|+|z|)` → 4-tap low-pass → centered **33-sample +window** max/min peak detection → a **dynamic threshold** (running mean of recent max/min +midpoints) with a **0.1 g sensitivity** dead-zone — a `max > thr+s/2` paired with a +`min < thr−s/2` is a *possible step*, and only after **8 consecutive** possible steps does it +start counting (the regularity gate that rejects waving/typing/handling — verified to read 0 +at rest; lowering it re-introduces false positives). A per-device **calibration gain (×1.11)** +corrects the typical ~10 % wrist undercount, locked against a 100-step ground-truth walk (raw +90 → 100). ESTIMATE. **Steps only accrue while the strap is connected** (the IMU is live-only; +the historical 1 Hz record carries no usable IMU). Owned by the cron (hourly: today+yesterday; +nightly: 2 days), written to `daily.steps` *after* analytics so the IMU value is authoritative. +*(Implemented in the backend, not this package, but documented here as the algorithm of record.)* + --- ## What we deliberately do **not** do @@ -258,4 +284,6 @@ Banister 1991 · Morton/Fitz-Clarke/Banister, *J Appl Physiol* 1990 · Cole & Kr et al., *Sports Med* 2013 · Baevsky & Berseneva 2008 · Mahalanobis 1936 · Mishra et al., *Nat Biomed Eng* 2020 · Smarr et al., *Sci Rep* 2020 · Radin et al., *Lancet Digit Health* 2020 · Fox 1971 · Uth, Sørensen, Overgaard & Pedersen, *Eur J Appl Physiol* 2004 · Foster, -*Med Sci Sports Exerc* 1998 · Brennan, Palaniswami & Kamen (Poincaré HRV), *IEEE TBME* 2001. +*Med Sci Sports Exerc* 1998 · Brennan, Palaniswami & Kamen (Poincaré HRV), *IEEE TBME* 2001 · +Analog Devices, *AN-2554: Step Counting Using the ADXL367* 2023 · Zhao, "Full-Featured +Pedometer Design Realized with 3-Axis Digital Accelerometer," *Analog Dialogue* 2010. diff --git a/README.md b/README.md index 0cb01e7..81f1de2 100644 --- a/README.md +++ b/README.md @@ -55,14 +55,15 @@ so it just doesn't. |----------|------|--------------| | `calcRestingHR` | `resting.ts` | 5th percentile of heart rate across your sleep window. Falls back to your quietest 30 minutes if there's no sleep yet. | | `calcStrain` | `strain.ts` | Banister TRIMP over heart-rate reserve, `ratio·0.64·e^(1.92·ratio)` summed per minute, squashed onto a 0–21 scale. | -| `calcHrZones` | `zones.ts` | Minutes spent in five zones by percent of max HR. | +| `calcHrZones` | `zones.ts` | Minutes spent in five zones by % of max HR. (Karvonen %HRR is more individualized in theory, but with an age-predicted max it adds no real accuracy and empties light-day zones, so %HRmax is kept deliberately.) | | `calcCalories` | `calories.ts` | Keytel (2005), the active-kcal-per-minute equation, summed. Different formula for men and women; averages the two if it doesn't know. | | `calcSleep` | `sleep.ts` | Cole-Kripke scores each epoch awake or asleep from motion, then I nudge it with the overnight HR dip. Gives onset, wake, efficiency, and a beta stage estimate. | -| `calcSleepRegularity` | `regularity.ts` | The Sleep Regularity Index, 0–100, from how much your bed and wake times wander night to night. | +| `calcSleepRegularity` | `regularity.ts` | Sleep-timing regularity, 0–100, from how much your bed and wake times wander night to night (circular variance of onset/wake clock-times). Honest scope: this is *not* the Phillips epoch-agreement Sleep Regularity Index — we don't plumb minute-level sleep/wake state across days, so we don't claim that name. | | `detectSessions` | `sessions.ts` | Finds workouts: sustained stretches above 40% heart-rate reserve, then classifies them roughly as cardio, strength, or a walk. | -| `calcHrRecovery` | `recovery.ts` | HRR60, how many beats your heart drops in the minute after a peak. A real fitness signal. | -| `calcLoad`, `calcFitnessTrend` | `trends.ts` | ACWR (last 7 days over last 28) for load, and regression slopes on resting HR and HRR for whether you're getting fitter. | -| `calcReadiness`, `calcAnomaly` | `readiness.ts` | Readiness from resting-HR deviation, sleep debt, and sleep quality. Plus a flag for "your resting HR has been up two days, are you getting sick?" | +| `timeDomainHrv`, `freqDomainHrv` | `hrv.ts` | HRV from the beat-to-beat R-R stream: RMSSD/SDNN/pNN50 and LF/HF (Lomb–Scargle, gated to the Task Force 1996 window minimums — HF ≥~60 s, LF ≥~250 s — so short windows don't report spectral noise), plus the Baevsky stress index. | +| `calcRecovery`, `calcHrRecovery` | `recovery.ts` | Recovery from nightly HRV — ln-RMSSD z-scored against your own baseline (Plews). Plus HRR60, the beats your heart drops in the minute after a peak. | +| `calcLoad`, `calcFitnessTrend` | `trends.ts` | EWMA acute:chronic workload ratio (Williams 2017 — 7/28-day exponentially-weighted, fixes the rolling-average coupling), and regression slopes on resting HR and HRR for whether you're getting fitter. | +| `calcReadinessIndex`, `calcAnomaly` | `readiness_index.ts`, `readiness.ts` | An HRV-led readiness composite: recovery blended with sleep, the nocturnal dip, and arousal (abstains until there's HRV). Plus a flag for "your resting HR has been up two days, are you getting sick?" | | `calcBaselines` | `baselines.ts` | Rolling 30-day medians, the anchors everything else compares against. | | `calcStress`, `classifyArousal` | `stress.ts` | Arousal from heart rate sitting above resting while you're not moving. If you're moving it's exercise, not stress, so it's gated out. | | `calcNocturnalHeart` | `nocturnal.ts` | Your sleeping HR, its low point, how far it dipped from daytime, and a flag if it's running high. | @@ -71,15 +72,20 @@ so it just doesn't. A couple of things worth calling out so you don't go looking for them: -**There's no HRV in here, and there won't be.** WHOOP builds recovery on heart-rate -variability, the beat-to-beat timing. The band doesn't hand that to us over the wire, at -least not in any form I've been able to recover. So readiness here is built from resting -HR, sleep debt, and sleep quality instead, and it says so right on the label. The -`hrv.ts` file is empty on purpose, as a reminder of what we don't have. +**HRV is in now** — this section used to say it never would be. WHOOP builds recovery on +heart-rate variability, the beat-to-beat timing, and for a long time it looked like the +band never handed that over. It turns out the R-R intervals are sitting right there in the +1 Hz historical (V24) records; they just don't ride the live stream, so the +[backend](https://github.com/OpenStrap/backend) re-decodes them from the raw bytes off the +request path and feeds them in. `hrv.ts` does the time- and frequency-domain measures, +`recovery.ts` turns nightly ln-RMSSD into a recovery z-score, and readiness is now an +HRV-led composite. It's labelled beta because recovering the field from the bytes is +empirical — but it's the real beat-to-beat signal, the same substrate WHOOP uses. **Max heart rate** falls back gracefully: a real measured peak from your workouts if I've -seen one, otherwise the highest I've observed, otherwise `220 − age`, otherwise 190. The -worse the source, the lower the confidence on anything that depends on it. +seen one, otherwise the highest I've observed, otherwise Tanaka `208 − 0.7·age` (more +accurate than the old `220 − age`), otherwise 190. The worse the source, the lower the +confidence on anything that depends on it. ## Tests @@ -98,5 +104,5 @@ Write a function that takes minutes (or history) plus the baseline and profile, returns a `Metric`. Keep it pure, no side effects. Derive the confidence from coverage and completeness like the others do, and return `null` with `0` confidence when you don't have the inputs. Put the name of the method you used in a comment so the next -person can check your work. And if your idea needs HRV or some signal the band doesn't +person can check your work. And if your idea needs a signal the band genuinely doesn't expose, it doesn't belong here, that's the line. diff --git a/src/__tests__/_harness.ts b/src/__tests__/_harness.ts index 1312457..98d1acc 100644 --- a/src/__tests__/_harness.ts +++ b/src/__tests__/_harness.ts @@ -30,7 +30,7 @@ export function min( ts: number, hr: number, activity = 0, - opts: Partial<{ steps: number; wrist_on: boolean; hr_max: number }> = {} + opts: Partial<{ steps: number; wrist_on: boolean; hr_max: number; act_class: import('../types').ActivityClass }> = {} ) { return { ts, @@ -41,5 +41,6 @@ export function min( activity, steps: opts.steps ?? 0, wrist_on: opts.wrist_on ?? hr > 0, + ...(opts.act_class ? { act_class: opts.act_class } : {}), }; } diff --git a/src/__tests__/analytics.test.ts b/src/__tests__/analytics.test.ts index 5807cfc..f9c5a7d 100644 --- a/src/__tests__/analytics.test.ts +++ b/src/__tests__/analytics.test.ts @@ -7,19 +7,30 @@ import { calcRestingHR } from '../resting'; import { calcStrain } from '../strain'; import { calcHrZones } from '../zones'; import { calcCalories } from '../calories'; -import { calcSleep } from '../sleep'; +import { calcSleep, calcSleepPeriods } from '../sleep'; import { calcSleepRegularity } from '../regularity'; import { detectSessions } from '../sessions'; import { calcHrRecovery, calcRecovery } from '../recovery'; import { calcLoad, calcFitnessTrend } from '../trends'; import { calcAnomaly } from '../readiness'; -import { calcBaselines } from '../baselines'; +import { calcBaselines, seedBaselines } from '../baselines'; import { calcStress } from '../stress'; +import { calcSpo2Index } from '../spo2'; import { calcSleepStress } from '../arousal'; import { calcNocturnalHeart } from '../nocturnal'; import { calcIllness } from '../illness'; import { timeDomainHrv, freqDomainHrv, baevskyStressIndex, cleanRr } from '../hrv'; +import { pedometer, calcSteps, STEP_PARAMS } from '../steps'; import { resolveMaxHr } from '../util'; +import { calcCircadian, stageSleep } from '../circadian'; +import { detectSleepCycles } from '../cycles'; +import { detectWakeState, peekRecentState } from '../wake'; +import { calcCycle } from '../cycle'; +import { extractHarFeatures, classifyActivityWindow, segmentWorkout, DB10_LO, dwtDetailEnergies } from '../har'; +import type { ClassVote } from '../har'; +import { calcRestlessness } from '../restlessness'; +import { calcDaytimeHrv } from '../hrv'; +import { calcDesaturation } from '../spo2'; const baseline: Baseline = { resting_hr: 50, @@ -84,7 +95,7 @@ console.log('--- §2 calcStrain ---'); // ── §3 calcHrZones ─────────────────────────────────────────────────────────── console.log('--- §3 calcHrZones ---'); { - // maxHR 190: 100bpm=52.6%→z1, 120=63%→z2, 150=78.9%→z3, 160=84%→z4, 175=92%→z5 + // maxHR 190 (%HRmax): 100bpm=52.6%→z1, 120=63%→z2, 150=78.9%→z3, 160=84%→z4, 175=92%→z5 const mins: Minute[] = [100, 120, 150, 160, 175].map((h, i) => min(i * 60, h)); const z = calcHrZones(mins, baseline); assert(z.zone1_min === 1, '100bpm → z1'); @@ -106,7 +117,7 @@ console.log('--- §3 calcHrZones ---'); { age: 40 } ); assert(ageOnly.max_hr_source === 'age' && ageOnly.max_hr_used === 180, - 'age fallback maxHR = 220−age when no measured HR'); + 'age fallback maxHR = Tanaka 208−0.7·age (=180 at age 40) when no measured HR'); } // ── §4 calcCalories ────────────────────────────────────────────────────────── @@ -206,6 +217,29 @@ console.log('--- §5 calcSleep ---'); assert(ow.in_bed_min <= 210, `off-wrist: off-wrist stretch not counted as in-bed (got ${ow.in_bed_min})`); + // Regression (v3): the band's flash record (R24 — the entire overnight) carries NO + // actigraphy, so `activity` is ~0 for every sleep minute and Cole-Kripke is inert; the + // call is HR-driven. A night's HR legitimately runs ABOVE the 5th-pctile RHR floor, so + // the old fixed `1.15*rhr` awake-override flagged the whole night awake → 1-min "sleep" + // (observed in prod). The window-relative reference must detect it. rhr=50 floor, night + // HR ~62 (well above 1.15*50=57.5), activity 0, bracketed by waking HR ~95. + const aboveFloor: Minute[] = []; + for (let i = 0; i < 30; i++) aboveFloor.push(min(i * 60, 95, 0)); // waking evening + for (let i = 30; i < 430; i++) aboveFloor.push(min(i * 60, 62, 0)); // ~6.7h night, HR > RHR floor + for (let i = 430; i < 470; i++) aboveFloor.push(min(i * 60, 95, 0)); // waking morning + const af = calcSleep(aboveFloor, baseline); + assert(af.duration_min >= 380, + `above-floor night (HR>RHR, activity inert) is detected, not clipped to ~1min (got ${af.duration_min})`); + + // Regression (v3): the OTHER failure mode — with activity inert, a flat + ELEVATED + // sedentary-awake window must NOT be manufactured into a full night. The absolute + // backstop (HR > 1.5*RHR → awake) clips it. Flat 92 bpm for 5h, rhr 50 → 92 > 75. + const flatAwake: Minute[] = []; + for (let i = 0; i < 300; i++) flatAwake.push(min(i * 60, 92, 0)); + const fa = calcSleep(flatAwake, baseline); + assert(fa.duration_min <= 30, + `flat elevated sedentary-awake window is not mis-read as sleep (got ${fa.duration_min})`); + // Plausibility guard: even if (pathologically) almost the whole 18h window // reads low-activity, a single main-sleep period must never exceed ~14h. const giant: Minute[] = []; @@ -215,6 +249,49 @@ console.log('--- §5 calcSleep ---'); `plausibility: main-sleep period capped at ~14h (got ${g.in_bed_min})`); } +// ── §5b calcSleepPeriods (multi-period; naps = shorter sleeps) ─────────────── +console.log('--- §5b calcSleepPeriods ---'); +{ + // A main night (200 min) + a separate afternoon nap (40 min), split by a long + // awake daytime stretch. v2 must return TWO periods, longest flagged is_main. + const day: Minute[] = []; + for (let i = 0; i < 5; i++) day.push(min(i * 60, 70, 2000)); // awake before + for (let i = 5; i < 205; i++) day.push(min(i * 60, 45, 50)); // 200-min night (asleep) + for (let i = 205; i < 305; i++) day.push(min(i * 60, 75, 2500)); // 100-min awake daytime (long gap) + for (let i = 305; i < 345; i++) day.push(min(i * 60, 48, 50)); // 40-min afternoon nap (asleep) + for (let i = 345; i < 350; i++) day.push(min(i * 60, 72, 2000)); // awake after + const v2 = calcSleepPeriods(day, baseline); + assert(v2.periods.length === 2, `two sleep periods detected (got ${v2.periods.length})`); + const mainP = v2.periods[v2.main_idx ?? -1]; + assert(mainP != null && mainP.is_main === true, 'main period flagged is_main'); + assert(mainP.duration_min >= 195, `main period ≈ 200 min (got ${mainP?.duration_min})`); + const napP = v2.periods.find((p) => !p.is_main); + assert(napP != null && napP.duration_min >= 30 && napP.duration_min <= 45, + `nap treated as a shorter sleep ≈ 40 min, edge-trimmed (got ${napP?.duration_min})`); + assert(v2.total_asleep_min >= 228, `total asleep sums both periods (got ${v2.total_asleep_min})`); + assert(v2.periods.every((p) => p.confidence >= 0 && p.confidence <= 1), 'per-period confidence in [0,1]'); + + // Single night → exactly one period (backward-consistent with calcSleep). + const oneNight: Minute[] = []; + for (let i = 0; i < 5; i++) oneNight.push(min(i * 60, 70, 2000)); + for (let i = 5; i < 205; i++) oneNight.push(min(i * 60, 45, 50)); + for (let i = 205; i < 210; i++) oneNight.push(min(i * 60, 72, 2000)); + const one = calcSleepPeriods(oneNight, baseline); + assert(one.periods.length === 1 && one.periods[0].is_main, 'single night → one main period'); + + // Micro-doze (< 15 min) is discarded, not surfaced as a period. + const micro: Minute[] = []; + for (let i = 0; i < 5; i++) micro.push(min(i * 60, 70, 2000)); + for (let i = 5; i < 13; i++) micro.push(min(i * 60, 45, 50)); // 8-min doze + for (let i = 13; i < 20; i++) micro.push(min(i * 60, 72, 2000)); + const m2 = calcSleepPeriods(micro, baseline); + assert(m2.periods.length === 0 && m2.confidence === 0, 'micro-doze (<15 min) discarded'); + + // Empty input → no periods, conf 0. + const ep = calcSleepPeriods([], baseline); + assert(ep.periods.length === 0 && ep.confidence === 0, 'no data → no periods + conf 0'); +} + // ── §6 calcSleepRegularity ─────────────────────────────────────────────────── console.log('--- §6 calcSleepRegularity ---'); { @@ -257,12 +334,28 @@ console.log('--- §7 detectSessions ---'); assert(ses.type === 'run/cardio', 'high act + high HR → run/cardio'); assert(ses.strain > 0 && ses.kcal > 0, 'session carries strain + calories'); - // a <5min bout is discarded + // a 2-min sustained bout now qualifies (threshold lowered 3/5 → 2 min). + const short: Minute[] = []; + for (let i = 0; i < 10; i++) short.push(min(i * 60, 55, 1)); + for (let i = 10; i < 12; i++) short.push(min(i * 60, 150, 100)); // 2 min + for (let i = 12; i < 20; i++) short.push(min(i * 60, 55, 1)); + assert(detectSessions(short, baseline).length === 1, '2-min bout now detected'); + + // a 1-min blip is still discarded. const tiny: Minute[] = []; for (let i = 0; i < 10; i++) tiny.push(min(i * 60, 55, 1)); - for (let i = 10; i < 13; i++) tiny.push(min(i * 60, 150, 100)); // only 3 min - for (let i = 13; i < 20; i++) tiny.push(min(i * 60, 55, 1)); - assert(detectSessions(tiny, baseline).length === 0, 'bout <5 min discarded'); + tiny.push(min(10 * 60, 150, 100)); // 1 min only + for (let i = 11; i < 20; i++) tiny.push(min(i * 60, 55, 1)); + assert(detectSessions(tiny, baseline).length === 0, '1-min blip discarded'); + + // Per-minute HAR class → motion-based workout type + confidence (not the HR heuristic). + const cyc: Minute[] = []; + for (let i = 0; i < 10; i++) cyc.push(min(i * 60, 55, 1)); + for (let i = 10; i < 20; i++) cyc.push(min(i * 60, 140, 100, { hr_max: 150, act_class: 'cycle' })); + for (let i = 20; i < 30; i++) cyc.push(min(i * 60, 55, 1)); + const cs = detectSessions(cyc, baseline)[0]; + assert(cs.type === 'cycle' && cs.type_confidence > 0.4, `motion class → cycle (got ${cs.type}/${cs.type_confidence})`); + assert(cs.detected_type === 'cycle', 'detected_type recorded for calibration ledger'); } // ── §8 calcHrRecovery ──────────────────────────────────────────────────────── @@ -293,13 +386,16 @@ console.log('--- §9 calcLoad / calcFitnessTrend ---'); approx(load.acwr ?? -1, 1.0, 1e-9, 'steady strain → ACWR 1.0'); assert(load.band === 'optimal', 'ACWR 1.0 → optimal band'); - // acute spike: last 7 days at 20, prior 21 at 10 → acute 20, chronic ~12.5 → acwr 1.6 + // acute spike (EWMA, Williams 2017): prior 21 days at 10, last 7 at 20. + // EWMA acute (λ=0.25) rises to ~18.7, chronic (λ≈0.069) to ~13.9 → acwr ~1.34 + // (caution). EWMA is deliberately smoother than the rolling-average ratio. const spike = Array.from({ length: 28 }, (_, i) => ({ ts: i * 86400, strain: i >= 21 ? 20 : 10, })); const sl = calcLoad(spike); - assert(sl.band === 'high-risk', `acute spike → high-risk (acwr ${sl.acwr})`); + assert(sl.band === 'caution' && (sl.acwr ?? 0) > 1.3, + `acute spike → caution (acwr ${sl.acwr})`); // <7 days → unknown assert(calcLoad(steady.slice(0, 5)).band === 'unknown', '<7 days → unknown band'); @@ -337,6 +433,10 @@ console.log('--- §10 calcRecovery / calcAnomaly / calcIllness ---'); const an = calcAnomaly({ recent_rhr: [50, 51, 55, 56] }, baseline); assert(an.signal === true && an.triggers.includes('rhr_elevated_2d'), 'two elevated RHR days → signal'); assert(an.note === 'signal, not a diagnosis', 'anomaly non-diagnostic note'); + // §4 cycle gate: luteal phase suppresses the pure-RHR-elevation rule (expected rise). + const anLuteal = calcAnomaly({ recent_rhr: [50, 51, 55, 56] }, baseline, { cyclePhase: 'luteal' }); + assert(anLuteal.signal === false && /cycle/i.test(anLuteal.note), + 'luteal phase suppresses pure-RHR anomaly with a cycle note'); // illness (Mahalanobis): RHR↑ + RMSSD↓ + temp↑ vs baseline → signal. const hist = { @@ -349,6 +449,26 @@ console.log('--- §10 calcRecovery / calcAnomaly / calcIllness ---'); assert(sick.note === 'a signal, not a diagnosis', 'illness non-diagnostic note'); const well = calcIllness({ resting_hr: 56, rmssd: 76, skin_temp: 34.05 }, hist); assert(well.signal === false, 'normal vector → no illness signal'); + + // §5 respiratory rate as a 4th Mahalanobis feature: RMSSD↓ + resp↑ fires + lists 'resp'. + const histR = { ...hist, resp_rate: Array.from({ length: 20 }, (_, i) => 14 + (i % 2)) }; + const sickResp = calcIllness({ resting_hr: 56, rmssd: 45, skin_temp: 34.05, resp_rate: 19 }, histR); + assert(sickResp.signal === true && sickResp.triggers.includes('resp'), + 'elevated respiratory rate drives the illness signal'); + assert(sickResp.inputs_used.includes('resp_rate'), 'resp_rate listed in inputs_used'); + + // §4 cycle gating: temp+RHR rise ALONE is phase-expected → suppressed in luteal, + // but still fires when no cycle context is supplied. + const cycIn = { resting_hr: 64, rmssd: 76, skin_temp: 35.0 }; + const noCyc = calcIllness(cycIn, hist); + assert(noCyc.signal === true && noCyc.triggers.includes('rhr') && noCyc.triggers.includes('temp'), + 'temp+RHR rise → illness signal with no cycle context'); + const luteal = calcIllness(cycIn, hist, { cyclePhase: 'luteal' }); + assert(luteal.signal === false, 'luteal phase suppresses temp/RHR-only illness signal'); + assert(/cycle/i.test(luteal.note), 'suppressed signal explains the cycle phase'); + // …but HRV/resp deviations are NOT explained by the cycle → still fires. + const lutealReal = calcIllness({ resting_hr: 64, rmssd: 45, skin_temp: 35.0, resp_rate: 19 }, histR, { cyclePhase: 'luteal' }); + assert(lutealReal.signal === true, 'HRV/resp shift still fires even in luteal phase'); } // ── §11 calcBaselines ──────────────────────────────────────────────────────── @@ -362,17 +482,25 @@ console.log('--- §11 calcBaselines ---'); session_hr_max: 180 + (i % 10), zone_min: [10, 20, 15, 5, 2] as [number, number, number, number, number], })); - const bl = calcBaselines(hist); + // Observed daily peak (189) CLEARS the age-predicted floor (age 40 → Tanaka 180), + // so it's a genuine effort → trusted as the measured max. + const bl = calcBaselines(hist, { age: 40 }); assert(bl.resting_hr !== null && bl.resting_hr! >= 50 && bl.resting_hr! <= 52, 'RHR median in range'); - assert(bl.max_hr === 189 && bl.max_hr_source === 'measured', 'maxHR = max observed session'); + assert(bl.max_hr === 189 && bl.max_hr_source === 'measured', 'observed peak above age floor → measured'); assert(bl.chronic_strain !== null, 'chronic strain computed'); assert(bl.zone_min !== null && bl.zone_min![0] === 10, 'per-zone medians present'); approx(bl.confidence, 1, 1e-9, '30 days → confidence 1'); + // Guard: a quiet daily peak (≤ age floor) must NOT be promoted to a measured max — + // it would under-state HRmax and inflate zones/strain. Age floor wins instead. + const quiet = hist.map((d) => ({ ...d, session_hr_max: 150 })); + const blQuiet = calcBaselines(quiet, { age: 30 }); // Tanaka 187 > 150 + assert(blQuiet.max_hr === 187 && blQuiet.max_hr_source === 'age', 'quiet daily peak → age floor, not measured'); + // age fallback for maxHR when no sessions const noSess = hist.map((d) => ({ ...d, session_hr_max: undefined })); const bl2 = calcBaselines(noSess, { age: 30 }); - assert(bl2.max_hr === 190 && bl2.max_hr_source === 'age', 'no sessions + age → 220−age maxHR'); + assert(bl2.max_hr === 187 && bl2.max_hr_source === 'age', 'no sessions + age → Tanaka 208−0.7·age maxHR'); // seed period (3 days) → low confidence const seed = calcBaselines(hist.slice(0, 3)); @@ -441,7 +569,7 @@ console.log('--- §HRV time/freq/SI ---'); // Respiratory peak: RR modulated at 0.25 Hz (15 brpm) → resp_rate ≈ 15. const t: number[] = []; let acc = 0; const resp: number[] = []; - for (let i = 0; i < 200; i++) { + for (let i = 0; i < 320; i++) { // ≥250 s span so LF (and LF/HF) are valid per Task Force 1996 const rr = 900 + 60 * Math.sin(2 * Math.PI * 0.25 * (acc / 1000)); resp.push(Math.round(rr)); acc += rr; } @@ -477,6 +605,28 @@ console.log('--- §12 calcStress (HRV) ---'); assert(JSON.stringify(calcStress(rr, baseSI)) === JSON.stringify(calcStress(rr, baseSI)), 'stress deterministic'); } +// ── §SpO₂ relative index (red/IR ratio) ─────────────────────────────────────── +console.log('--- §calcSpo2Index ---'); +{ + // Too few minutes → null, conf 0. + assert(calcSpo2Index([0.85, 0.86, 0.85], 0.85).index === null, 'spo2: <30 min → null'); + assert(calcSpo2Index([0.85, 0.86], 0.85).confidence === 0, 'spo2: too few → conf 0'); + // No baseline yet → seed night_ratio, null index. + const stable = Array.from({ length: 200 }, () => 0.850); + const seed = calcSpo2Index(stable, null); + assert(seed.index === null, 'spo2: no baseline → null index'); + approx(seed.night_ratio!, 0.85, 0.001, 'spo2: no baseline → seed night_ratio'); + // Stable clean night vs baseline → high confidence; lower ratio than baseline → positive index. + const better = calcSpo2Index(Array.from({ length: 200 }, () => 0.840), 0.850); + assert(better.index !== null && better.index > 0, 'spo2: lower ratio than baseline → positive index'); + assert(better.confidence > 0.8, 'spo2: stable + plenty of samples → high confidence'); + // Noisy night (high intra-night CV) → low confidence even with a baseline. + const noisy = calcSpo2Index(Array.from({ length: 200 }, (_, i) => 0.85 + (i % 2 ? 0.08 : -0.08)), 0.850); + assert(noisy.confidence < 0.3, 'spo2: high intra-night CV → low confidence'); + // Plausibility gate drops garbage ratios. + assert(calcSpo2Index(Array.from({ length: 200 }, () => 3.0), 0.85).index === null, 'spo2: implausible ratios → null'); +} + // ── §sleep-stress / nocturnal arousal ───────────────────────────────────────── console.log('--- §calcSleepStress ---'); { @@ -595,15 +745,91 @@ console.log('--- regression: sleep stage proportions ---'); assert(st.light_min >= st.rem_min, `light ≥ REM (light ${st.light_min} vs REM ${st.rem_min})`); } +// ── regression: stageSleep detects REM from RR variability on a flat-HR night ── +console.log('--- regression: RR-driven REM staging ---'); +{ + // A calm night where HR is nearly flat (≈60 bpm) so HR LEVEL alone CANNOT separate + // REM from light (the real bug: REM read 0%). REM is encoded the physiological way — + // parasympathetic withdrawal → REDUCED beat-to-beat variability (low RMSSD), with HR + // only mildly above light. Deep = high RMSSD + lowest HR. stageSleep must use the RR + // autonomic axis to recover a physiological REM share (15–25%), not 0. + const ONSET = 0, WAKE = 280 * 60; + // Build a minute's RR stream with a target RMSSD: alternate ±d around the mean RR so + // successive |Δ| ≈ 2d ⇒ RMSSD ≈ 2d (kept < 200 ms so cleanRr doesn't drop beats). + const rrFor = (hr: number, rmssdTarget: number): number[] => { + const meanRr = Math.round(60000 / hr); + const d = Math.min(95, Math.round(rmssdTarget / 2)); + const out: number[] = []; + for (let j = 0; j < 48; j++) out.push(meanRr + (j % 2 === 0 ? d : -d)); + return out; + }; + type SM = { ts: number; hr_avg: number; rr?: number[] }; + const night: SM[] = []; + const push = (a: number, b: number, hr: number, rmssd: number) => { + for (let i = a; i < b; i++) night.push({ ts: i * 60, hr_avg: hr, rr: rrFor(hr, rmssd) }); + }; + push(0, 30, 60, 50); // light (medium variability) + push(30, 95, 56, 90); // deep (high RMSSD, lowest HR) + push(95, 150, 60, 50); // light + push(150, 215, 64, 16); // REM (low RMSSD, mildly elevated HR) + push(215, 280, 60, 50); // light + const ss = stageSleep(night, ONSET, WAKE, /*mesor*/ 90); + const tot = ss.light_min + ss.deep_min + ss.rem_min; + assert(tot > 0, 'RR-staged night produced stages'); + const remPct = (100 * ss.rem_min) / tot, deepPct = (100 * ss.deep_min) / tot; + assert(remPct >= 12 && remPct <= 35, `REM recovered from RR, physiological share (got ${remPct.toFixed(0)}%)`); + assert(deepPct >= 8, `deep detected from high-RMSSD block (got ${deepPct.toFixed(0)}%)`); + assert(ss.light_min >= ss.rem_min, 'light remains dominant'); + // 0 short(<20 min) awake flaps in the hypnogram. + let flaps = 0; + for (let i = 0; i < ss.hypnogram.length;) { + let j = i; while (j < ss.hypnogram.length && ss.hypnogram[j].stage === ss.hypnogram[i].stage) j++; + if (ss.hypnogram[i].stage === 'awake' && (j - i) < 20) flaps++; + i = j; + } + assert(flaps === 0, `no short awake flaps (got ${flaps})`); + // Without RR, the SAME flat-HR night cannot resolve REM → graceful HR-only fallback + // (must not throw, must not fabricate a REM-dominated night). + const noRr = night.map((m) => ({ ts: m.ts, hr_avg: m.hr_avg })); + const fb = stageSleep(noRr, ONSET, WAKE, 90); + assert((fb.light_min + fb.deep_min + fb.rem_min) > 0, 'HR-only fallback still stages'); +} + +// ── §Sleep cycles (fractal-cycle method on HRV) ─────────────────────────────── +console.log('--- §detectSleepCycles ---'); +{ + // RMSSD oscillating with a ~80-min ultradian period over a 320-min night → the + // findpeaks(20min, 0.9z) detector should recover ~4 peaks ⇒ ~3 cycles near 80 min. + // RR is built so each minute's RMSSD ≈ target: alternate ±d ⇒ RMSSD ≈ 2d. + const rrFor = (rmssd: number): number[] => { + const d = Math.max(2, Math.round(rmssd / 2)); + return Array.from({ length: 40 }, (_, j) => 900 + (j % 2 ? d : -d)); + }; + const mins: { ts: number; rr: number[] }[] = []; + for (let i = 0; i < 320; i++) { + const rmssd = 50 + 30 * Math.sin((2 * Math.PI * i) / 80); // ~80-min cycle + mins.push({ ts: i * 60, rr: rrFor(rmssd) }); + } + const c = detectSleepCycles(mins, 0, 319 * 60); + assert(c.n >= 2 && c.n <= 6, `cycles: ~3-4 ultradian cycles found (got ${c.n})`); + assert(c.mean_duration_min != null && c.mean_duration_min >= 55 && c.mean_duration_min <= 110, + `cycles: mean duration near the ~80-min period (got ${c.mean_duration_min})`); + assert(c.series.length > 0, 'cycles: emits a z-series for plotting'); + // No RR → abstain cleanly (no fabricated cycles). + const noRr = detectSleepCycles(Array.from({ length: 200 }, (_, i) => ({ ts: i * 60 })), 0, 199 * 60); + assert(noRr.n === 0 && noRr.cycles.length === 0, 'cycles: no RR → abstains'); +} + // ── regression: resolveMaxHr doesn't promote a quiet within-day peak ────────── console.log('--- regression: resolveMaxHr source ---'); { // No baseline max, age present, day peaks only at 110 bpm (a quiet day). Must - // use the age-predicted max (191) as the denominator, NOT call 110 "measured". + // use the age-predicted max (Tanaka: 208−0.7·29 ≈ 188) as the denominator, + // NOT call 110 "measured". const quiet: Minute[] = []; for (let i = 0; i < 60; i++) quiet.push(min(i * 60, 95 + (i % 5), 100, { hr_max: 110 })); const r1 = resolveMaxHr(quiet, { max_hr: 0 }, { age: 29 }); - assert(r1.source === 'age' && r1.maxHr === 191, + assert(r1.source === 'age' && r1.maxHr === 188, `quiet-day peak not promoted to measured (got ${r1.maxHr}/${r1.source})`); // A genuine hard effort above age-max IS taken as measured. const effort: Minute[] = []; @@ -616,4 +842,271 @@ console.log('--- regression: resolveMaxHr source ---'); assert(r3.source === 'measured' && r3.maxHr === 185, 'baseline max_hr wins'); } +// ── §Steps pedometer (AN-2554) ─────────────────────────────────────────────── +console.log('--- §Steps pedometer ---'); +{ + // Rest: flat 1g magnitude → the CONFIRM gate must read exactly 0. + const rest: number[] = new Array(3000).fill(1.0); + assert(pedometer(rest) === 0, 'rest signal → 0 steps (CONFIRM gate rejects non-gait)'); + + // Walk: clean ~1.8 Hz cadence over 30 s @100 Hz ≈ 54 gait cycles. Magnitude + // swings well outside the ±SENS dead-zone, so each cycle is a step. + const walk: number[] = []; + for (let i = 0; i < 3000; i++) walk.push(1.0 + 0.3 * Math.sin(2 * Math.PI * 1.8 * (i / 100))); + const raw = pedometer(walk); + assert(raw >= 45 && raw <= 60, `walk ~54 cycles → plausible step count (got ${raw})`); + + // calcSteps groups per-minute signals, sums, applies the calibration gain. + assert(calcSteps([rest]) === 0, 'calcSteps all-rest → 0'); + approx(calcSteps([walk]), Math.round(raw * STEP_PARAMS.GAIN), 0.001, + 'calcSteps applies the ×GAIN calibration to the summed raw count'); + assert(STEP_PARAMS.GAIN === 1.11, 'locked calibration gain = 1.11'); +} + +// ── §Circadian — CircaCP cosinor + bounded change-point ────────────────────── +console.log('--- §Circadian calcCircadian ---'); +{ + // 2 days of 1-min HR: asleep (hr≈55) hours [1,8), awake (hr≈80) otherwise. + // Onset ≈ 01:00, wake ≈ 08:00 each day; sharp transitions. + const mins: Minute[] = []; + for (let i = 0; i < 2 * 1440; i++) { + const ts = i * 60; + const hod = Math.floor(ts / 3600) % 24; + const asleep = hod >= 1 && hod < 8; + const hr = (asleep ? 55 : 80) + (i % 5) - 2; // tiny deterministic jitter + mins.push(min(ts, hr)); + } + const c = calcCircadian(mins); + assert(c.amplitude !== null && c.amplitude > 5, `circadian amplitude detected (got ${c.amplitude})`); + // day-2 onset ≈ 90000s (25:00 → 01:00 day 2), wake ≈ 115200s (32:00 → 08:00 day 2) + assert(c.onset_ts !== null && Math.abs(c.onset_ts - 90000) <= 1800, `onset ≈ 01:00 day2 (got ${c.onset_ts})`); + assert(c.wake_ts !== null && Math.abs(c.wake_ts - 115200) <= 1800, `wake ≈ 08:00 day2 (got ${c.wake_ts})`); + assert(c.settled === true, 'completed night marked settled'); + assert(c.confidence > 0.5, `confidence high on clean rhythm (got ${c.confidence})`); + + // flat HR (no rhythm) → abstain + const flat: Minute[] = []; + for (let i = 0; i < 2 * 1440; i++) flat.push(min(i * 60, 70)); + const cf = calcCircadian(flat); + assert(cf.onset_ts === null && cf.confidence < 0.3, 'flat HR → abstains (no fabricated boundary)'); +} + +// ── detectWakeState (sleep/wake ensemble) ───────────────────────────────────── +{ + const bl: Baseline = { resting_hr: 50, max_hr: 190, sleep_need_min: 480 }; + // 8h sleep (low HR, still) then N min awake (elevated HR, moving + steps). + const build = (awakeMin: number): Minute[] => { + const out: Minute[] = []; + let t = 0; + for (let i = 0; i < 480; i++, t += 60) out.push(min(t, 50, 0.01, { wrist_on: true })); + for (let i = 0; i < awakeMin; i++, t += 60) out.push(min(t, 72, 0.4, { steps: 20, wrist_on: true })); + return out; + }; + + const woke = detectWakeState({ minutes: build(15), baseline: bl }); + assert(woke.state === 'awake', `ensemble: state awake after waking (got ${woke.state})`); + assert(woke.wake_ts != null && Math.abs(woke.wake_ts - 480 * 60) <= 180, `ensemble: wake_ts ≈ sleep→wake boundary ±3min (got ${woke.wake_ts})`); + assert(woke.awake_min >= 12 && woke.awake_min <= 19, `ensemble: sustained awake ~15 min ±detector fuzz (got ${woke.awake_min})`); + assert(woke.asleep_min >= 90, `ensemble: main sleep ≥90 min (got ${woke.asleep_min})`); + + const tooSoon = detectWakeState({ minutes: build(5), baseline: bl }); + assert(tooSoon.wake_ts === null, 'ensemble: <10 min awake → no premature wake fire'); + + const stillAsleep = detectWakeState({ minutes: build(0), baseline: bl }); + assert(stillAsleep.state === 'asleep' && stillAsleep.wake_ts === null, 'ensemble: mid-sleep → asleep, no wake_ts'); + + const movingTail = [min(0, 72, 0.4, { steps: 20, wrist_on: true }), min(60, 73, 0.5, { steps: 25, wrist_on: true }), min(120, 71, 0.3, { steps: 10, wrist_on: true })]; + assert(peekRecentState(movingTail, bl) === 'awake', 'peek: moving + HR up → awake'); +} + +// ── regression: QUIET sedentary wake (HR up, NO motion, RR present) must fire ── +// The real-world bug: a user awake but still (on the phone in bed) has elevated HR +// and RR but ~zero motion. The OLD flat 2-of-3 majority let the two motion voters +// (blind to quiet wake) outvote cardiac → "asleep" → close never fired → no recovery. +// The ≥2 consensus must let the autonomic pair (cardiac + hrvArousal) carry the wake. +{ + const bl: Baseline = { resting_hr: 55, max_hr: 190, sleep_need_min: 480 }; + const minutes: Minute[] = []; + const rrByMin = new Map(); + let t = 0; + const rr = (meanMs: number, sd: number, n = 40) => Array.from({ length: n }, (_, j) => meanMs + (j % 2 ? sd : -sd)); + for (let i = 0; i < 480; i++, t += 60) { minutes.push(min(t, 52, 0.01, { wrist_on: true })); rrByMin.set(t, rr(1150, 12)); } // sleep: low HR, still, low RR-SD + for (let i = 0; i < 30; i++, t += 60) { minutes.push(min(t, 74, 0.01, { wrist_on: true })); rrByMin.set(t, rr(810, 60)); } // QUIET wake: HR up, NO motion, high RR-SD + + const ws = detectWakeState({ minutes, baseline: bl, rrByMin }); + assert(ws.state === 'awake', `quiet sedentary wake detected without motion (got ${ws.state})`); + assert(ws.wake_ts != null && Math.abs(ws.wake_ts - 480 * 60) <= 180, `quiet wake_ts at the boundary (got ${ws.wake_ts})`); + assert(ws.asleep_min >= 90, `main sleep preserved (got ${ws.asleep_min})`); + assert(ws.votes.cardiac === 'awake' && ws.votes.hrvArousal === 'awake', 'autonomic pair both vote awake at wake'); + + // Honest degradation: same still-but-awake tail with NO RR → only 1 signal (cardiac) + // → below the ≥2 bar → stays asleep rather than guess. + const noRr = detectWakeState({ minutes, baseline: bl }); + assert(noRr.wake_ts === null, `no-RR + no-motion quiet wake cannot be confirmed (honest) (got ${noRr.wake_ts})`); +} + +// ── menstrual cycle (log-anchored calendar method) ─────────────────────────── +{ + // No logs → empty/abstain. + const none = calcCycle([], '2026-06-20'); + assert(none.confidence === 0 && none.phase === 'unknown' && none.predicted_next === null, + 'cycle: no logs → abstain'); + + // Three regular 28-day starts → median 28, prediction = last + 28. + const starts = ['2026-04-04', '2026-05-02', '2026-05-30']; + const c = calcCycle(starts, '2026-06-06'); // day 8 of the cycle that began 05-30 + assert(c.mean_length === 28, `cycle: median length 28 (got ${c.mean_length})`); + assert(c.length_history.length === 2, `cycle: 2 observed lengths (got ${c.length_history.length})`); + assert(c.cycle_day === 8, `cycle: cycle day 8 (got ${c.cycle_day})`); + assert(c.predicted_next === '2026-06-27', `cycle: next period 06-27 (got ${c.predicted_next})`); + assert(c.ovulation_est === '2026-06-13', `cycle: ovulation = next−14 (got ${c.ovulation_est})`); + assert(c.fertile_start === '2026-06-08' && c.fertile_end === '2026-06-14', + `cycle: fertile window ov−5..ov+1 (got ${c.fertile_start}..${c.fertile_end})`); + assert(c.phase === 'follicular', `cycle: day 8 pre-ovulation → follicular (got ${c.phase})`); + assert(c.confidence > 0.5, `cycle: confidence grows with cycles (got ${c.confidence})`); + + // Menstruation window (day ≤ 5). + const m = calcCycle(starts, '2026-05-31'); // day 2 + assert(m.phase === 'menstruation', `cycle: day 2 → menstruation (got ${m.phase})`); + + // Luteal: after the fertile window, before next period. + const l = calcCycle(starts, '2026-06-20'); // day 22 + assert(l.phase === 'luteal', `cycle: day 22 → luteal (got ${l.phase})`); + + // Very overdue → prediction unreliable, abstain on phase + low confidence. + const od = calcCycle(starts, '2026-07-25'); // ~56 days since last start + assert(od.phase === 'unknown' && od.confidence <= 0.2, `cycle: very overdue → unknown/low conf (got ${od.phase}/${od.confidence})`); + + // Single log → 28-day default, low-but-nonzero confidence. + const one = calcCycle(['2026-06-10'], '2026-06-15'); + assert(one.mean_length === null && one.predicted_next === '2026-07-08' && one.confidence > 0, + `cycle: single log uses 28d default (got ${one.predicted_next}/${one.confidence})`); +} + +// ── §HAR — activity recognition (Mannini features + classifier + segmentation) ── +console.log('--- §HAR activity recognition ---'); +{ + // db10 orthonormality invariants — catch any coefficient transcription error. + const sumLo = DB10_LO.reduce((s, v) => s + v, 0); + const sumSq = DB10_LO.reduce((s, v) => s + v * v, 0); + assert(DB10_LO.length === 20, 'db10 has 20 taps'); + approx(sumLo, Math.SQRT2, 1e-6, 'db10 Σh = √2'); + approx(sumSq, 1, 1e-6, 'db10 Σh² = 1'); + + // Synthetic tri-axial window: gravity on Z + a sinusoidal swing at f0 on X. + const fs = 100, secs = 4, n = fs * secs; + // Oscillate the magnitude (gravity axis) at f0 so SMV ≈ 1 + amp·sin(2π f0 t) — this + // matches how the accel-vector magnitude actually varies with gait (avoids the sin² + // frequency-doubling artifact you get from a single off-axis sinusoid). + const mk = (f0: number, amp: number, noise = 0.004) => { + const x: number[] = [], y: number[] = [], z: number[] = []; + for (let i = 0; i < n; i++) { + const t = i / fs; + z.push(1 + amp * Math.sin(2 * Math.PI * f0 * t) + (((i * 7919) % 991) / 991 - 0.5) * noise); + x.push((((i * 1103515245 + 12345) % 1000) / 1000 - 0.5) * noise); + y.push((((i * 1103) % 997) / 997 - 0.5) * noise); + } + return { x, y, z }; + }; + + // Frequency detection: a 2.0 Hz swing → dom1_freq ≈ 2.0 (within bin resolution). + const w2 = mk(2.0, 0.5); + const f2 = extractHarFeatures(w2.x, w2.y, w2.z, fs); + approx(f2.dom1_freq, 2.0, 0.3, `HAR dom1_freq ≈ 2.0 (got ${f2.dom1_freq.toFixed(2)})`); + assert(f2.dom1_ratio > 0.25, 'HAR strong sine → periodic (high dom1_ratio)'); + + // Classification: flat (gravity only, tiny noise) → sedentary. + const flat = mk(1.0, 0.0); + assert(classifyActivityWindow(extractHarFeatures(flat.x, flat.y, flat.z, fs)).cls === 'sedentary', + 'HAR flat signal → sedentary'); + + // 2 Hz strong swing → a locomotion class (walk), not sedentary/other. + const cw = classifyActivityWindow(f2); + assert(cw.cls === 'walk', `HAR 2 Hz → walk (got ${cw.cls})`); + + // 2.8 Hz strong swing → run. + const w3 = mk(2.8, 0.6); + assert(classifyActivityWindow(extractHarFeatures(w3.x, w3.y, w3.z, fs)).cls === 'run', + 'HAR 2.8 Hz → run'); + + // wavelet detail energies present (6 levels), non-negative. + const we = dwtDetailEnergies(w2.x, 6); + assert(we.length === 6 && we.every((e) => e >= 0), 'db10 detail energies: 6 levels, ≥0'); + + // Segmentation: 5 min walk → 5 min run (one continuous bout) → two phases, primary = either. + const votes: ClassVote[] = []; + for (let t = 0; t < 300; t += 4) votes.push({ ts: 1000 + t, cls: 'walk', conf: 0.7 }); + for (let t = 300; t < 600; t += 4) votes.push({ ts: 1000 + t, cls: 'run', conf: 0.7 }); + const seg = segmentWorkout(votes); + assert(seg.segments.length === 2, `HAR segment: walk→run → 2 phases (got ${seg.segments.length})`); + assert(seg.segments[0].type === 'walk' && seg.segments[1].type === 'run', 'HAR phases ordered walk then run'); + + // A single-window blip inside a long run is smoothed away (no spurious phase). + const blip: ClassVote[] = []; + for (let t = 0; t < 600; t += 4) blip.push({ ts: 2000 + t, cls: t === 300 ? 'cycle' : 'run', conf: 0.7 }); + assert(segmentWorkout(blip).segments.length === 1, 'HAR single-window blip smoothed → one phase'); +} + +// ── §Restlessness / §Daytime HRV / §Desaturation ──────────────────────────── +console.log('--- §restlessness / daytime HRV / desaturation ---'); +{ + // Restlessness: a still night with a movement spike every 30 min → bouts detected. + const sleepMin: Minute[] = []; + for (let i = 0; i < 240; i++) { + const moving = i % 30 === 0; + sleepMin.push({ ts: 1000 + i * 60, hr_avg: 55, hr_min: 54, hr_max: 56, hr_n: 60, activity: moving ? 0.5 : 0.01, steps: 0, wrist_on: true }); + } + const rest = calcRestlessness(sleepMin); + assert(rest.score !== null && rest.movement_bouts >= 5, `restlessness: detects bouts (got ${rest.movement_bouts})`); + assert(rest.longest_still_min > 0 && rest.mobility_pct !== null, 'restlessness: still stretch + mobility'); + assert(calcRestlessness(sleepMin.slice(0, 5)).score === null, 'restlessness: <20 min → null'); + + // Daytime HRV: 60 min of RR bucketed into 5-min windows → per-window RMSSD series. + const byMin: { ts: number; rr: number[] }[] = []; + for (let i = 0; i < 60; i++) { + const rr: number[] = []; + for (let k = 0; k < 12; k++) rr.push(850 + ((i + k) % 5) * 15); + byMin.push({ ts: 1000 + i * 60, rr }); + } + const dh = calcDaytimeHrv(byMin, 300); + assert(dh.rmssd_median !== null && dh.n_windows >= 10, `daytime HRV: windows (got ${dh.n_windows})`); + assert(dh.series.length === dh.n_windows && dh.lowest_ts !== null, 'daytime HRV: series + lowest window'); + assert(calcDaytimeHrv([], 300).rmssd_median === null, 'daytime HRV: no RR → null'); + + // Desaturation: a 2-min dip (R↑ above baseline) every 20 min → events counted. + const ratios: number[] = []; + for (let i = 0; i < 120; i++) ratios.push(i % 20 < 2 ? 0.86 : 0.79); + const des = calcDesaturation(ratios, 0.80); + assert(des.events >= 4 && des.odi !== null, `desaturation: counts dips (got ${des.events})`); + assert(des.deepest_pct !== null && des.deepest_pct > 0, 'desaturation: reports deepest dip'); + const desNoBase = calcDesaturation(ratios, null); + assert(desNoBase.events === 0 && desNoBase.confidence === 0, 'desaturation: no baseline → abstain'); +} + +// ── §11b seedBaselines (profile priors) ────────────────────────────────────── +console.log('--- §11b seedBaselines ---'); +{ + const m = seedBaselines({ age: 30, sex: 'm' }); + assert(m.resting_hr === 62, `seed: male RHR prior (got ${m.resting_hr})`); + assert(m.max_hr === 187, `seed: Tanaka 208-0.7*30 (got ${m.max_hr})`); + assert(m.sleep_need_min === 480, 'seed: adult 8h sleep need'); + assert(m.sleeping_hr === 55, `seed: sleeping HR = resting-7 (got ${m.sleeping_hr})`); + approx(m.hrv_rmssd, 42, 0.5, `seed: RMSSD age-prior 60-0.6*30 (got ${m.hrv_rmssd})`); + + const f = seedBaselines({ age: 66, sex: 'f' }); + assert(f.resting_hr === 66, 'seed: female RHR prior'); + assert(f.max_hr === 162, `seed: Tanaka @66 (got ${f.max_hr})`); + assert(f.sleep_need_min === 450, 'seed: ≥65 → 450'); + assert(f.hrv_rmssd === 20, `seed: RMSSD @66 (got ${f.hrv_rmssd})`); + + const young = seedBaselines({ age: 20 }); + assert(young.sleep_need_min === 510, 'seed: <25 → 510'); + assert(young.hrv_rmssd === 48, `seed: RMSSD @20 (got ${young.hrv_rmssd})`); + + const none = seedBaselines(); + assert(none.resting_hr === 64, 'seed: no sex → neutral RHR 64'); + assert(none.max_hr === Math.round(208 - 0.7 * 35), 'seed: no age → default age 35 maxHR'); + assert(none.resp_rate === 15 && none.hrv_si === 90, 'seed: resp + SI constants'); +} + summary('analytics'); diff --git a/src/baselines.ts b/src/baselines.ts index 036d485..466f35c 100644 --- a/src/baselines.ts +++ b/src/baselines.ts @@ -1,6 +1,56 @@ // §11 Baselines — rolling 30-day medians feeding everything. import type { DayHistory, Profile, Metric, BaselinesValue } from './types'; -import { median, round } from './util'; +import { median, round, clamp } from './util'; + +/** Population/profile-prior baselines for a brand-new user (see seedBaselines). */ +export interface SeededBaselines { + resting_hr: number; + max_hr: number | null; + sleep_need_min: number; + hrv_rmssd: number; + sleeping_hr: number; + resp_rate: number; + hrv_si: number; +} + +/** + * seedBaselines(profile?) + * + * Population/profile PRIORS for a user with no measured history yet. These are + * low-confidence anchors, NOT measured truth — they let strain/HR-zones use a + * personal resting/max HR from day 1, and warm-start the EWMA-rolled autonomic + * baselines (RMSSD, sleeping-HR, SI) from physiological norms instead of cold, + * self-referential seeding on night 1. Every value is overwritten by real data as + * it accrues: calcBaselines' rolling medians, the nightly EWMA rolls, and the + * ≥5-night gated scores (recovery/stress) — which still require real nights and + * are never fabricated by these priors. + * + * Population medians (healthy adults): + * resting HR — ~62 (m) / ~66 (f) (HUNT3 / NHANES adult medians). + * max HR — Tanaka 208 − 0.7·age (Tanaka et al., JACC 2001). + * sleep need — 8 h adult default (NSF 2015), nudged at age extremes. + * RMSSD — age-declining (Umetani et al., JACC 1998): 60 − 0.6·age, 18–60 ms. + * sleeping HR — resting − ~7 (typical nocturnal dip). + * resp rate — 15 br/min (normal resting adult). + * Baevsky SI — ~90 (normal resting sympathetic tone). + */ +export function seedBaselines(profile?: Profile): SeededBaselines { + const age = profile?.age && profile.age > 0 ? profile.age : 35; + const sex = profile?.sex; + const restingHr = sex === 'f' ? 66 : sex === 'm' ? 62 : 64; + const maxHr = Math.round(208 - 0.7 * age); + const sleepNeed = age < 25 ? 510 : age >= 65 ? 450 : 480; + const rmssd = clamp(round(60 - 0.6 * age, 0), 18, 60); + return { + resting_hr: restingHr, + max_hr: maxHr > 0 ? maxHr : null, + sleep_need_min: sleepNeed, + hrv_rmssd: rmssd, + sleeping_hr: restingHr - 7, + resp_rate: 15, + hrv_si: 90, + }; +} /** * calcBaselines(history, profile?) @@ -55,17 +105,34 @@ export function calcBaselines( ? (zoneCols.map((c) => median(c) ?? 0) as [number, number, number, number, number]) : null; - // maxHR = max observed session hr_max; else 220−age. + // maxHR: the highest per-day peak observed across the window. BUT a daily peak + // on a sedentary day is just a quiet high (e.g. 120 bpm), NOT a true HRmax — + // trusting it as the zone/strain denominator under-states the scale and inflates + // those metrics. So we only treat the observed peak as a real 'measured' max when + // it EXCEEDS the age-predicted Tanaka floor (208−0.7·age, JACC 2001 — a genuine + // hard effort); otherwise the age floor wins. Mirrors resolveMaxHr (util.ts) so + // the baseline and the per-day denominator agree. const observedMax = window .map((d) => d.session_hr_max) .filter((x): x is number => x != null); - let maxHr: number | null = null; - let maxHrSource: 'measured' | 'age' = 'measured'; - if (observedMax.length > 0) { - maxHr = Math.max(...observedMax); - maxHrSource = 'measured'; - } else if (profile?.age && profile.age > 0) { - maxHr = 220 - profile.age; + const observedPeak = observedMax.length > 0 ? Math.max(...observedMax) : 0; + const ageMax = + profile?.age && profile.age > 0 ? Math.round(208 - 0.7 * profile.age) : null; + let maxHr: number | null; + let maxHrSource: 'measured' | 'age'; + if (ageMax != null) { + if (observedPeak > ageMax) { + maxHr = observedPeak; // genuine above-age effort → trust it + maxHrSource = 'measured'; + } else { + maxHr = ageMax; // a quiet daily peak can't shrink the scale + maxHrSource = 'age'; + } + } else if (observedPeak > 0) { + // No age to floor against: the observed peak is the best available, but flag it + // 'age' (an unverified within-window peak, not a true measured HRmax) so callers + // down-weight confidence. + maxHr = observedPeak; maxHrSource = 'age'; } else { maxHr = null; // honest: no measured max, no age → cannot derive @@ -79,7 +146,7 @@ export function calcBaselines( if (sleeps.length) inputs_used.push('sleep_duration_min'); if (temps.length) inputs_used.push('skin_temp'); if (strains.length) inputs_used.push('daily_strain'); - if (observedMax.length) inputs_used.push('session_hr_max'); + if (maxHrSource === 'measured') inputs_used.push('session_hr_max'); else if (profile?.age) inputs_used.push('profile.age'); return { diff --git a/src/circadian.ts b/src/circadian.ts new file mode 100644 index 0000000..80d9255 --- /dev/null +++ b/src/circadian.ts @@ -0,0 +1,389 @@ +// §Circadian — CircaCP (Chen & Sun 2021, arXiv:2111.14960): robust cosinor model +// of the circadian rhythm + a single bounded change-point per cycle for sleep +// onset / wake. Validated on the dev account: wake within ~2 min, onset within +// ~2 min of the Cole-Kripke detector on real WHOOP 1 Hz HR. +// +// WHY HR (not actigraphy): the paper uses minute actigraphy, but the method is +// signal-agnostic and HR carries a strong, clean circadian rhythm (MESOR≈79, +// amplitude≈13 bpm on the dev account) with a sharp wake transition (HR rises +// ~25 bpm on waking) — a better wake discriminator than our near-zero activity +// rollup. Published method, no training, no fabrication; pure & deterministic. +// +// Division of labour: this returns the circadian phase + the MAIN-sleep +// onset/wake boundary of the most-recent completed cycle. In-sleep staging, +// efficiency and naps stay with calcSleep / calcSleepPeriods over that window. +import type { Minute, Metric, CircadianValue } from './types'; +import { isHrUsable, clamp, round, percentile } from './util'; +import { cleanRr } from './hrv'; + +const DAY = 86400; +const W = (2 * Math.PI) / DAY; // circadian angular frequency (rad/sec) + +interface Pt { t: number; y: number } + +/** + * Robust cosinor: y ≈ M + b1·cos(ωt) + b2·sin(ωt), fit by IRLS with a Tukey + * biweight so the active-phase HR spikes (exercise) don't drag the rhythm. + * Returns MESOR M, and (b1,b2) → amplitude = hypot, phase φ = atan2(b2,b1) + * with b1·cos+b2·sin = amp·cos(ωt − φ). + */ +function fitCosinor(pts: Pt[]): { mesor: number; b1: number; b2: number; amp: number; phi: number } | null { + const n = pts.length; + if (n < 120) return null; // need at least ~2h of points for a stable fit + const rows = pts.map((p) => ({ c: Math.cos(W * p.t), s: Math.sin(W * p.t), y: p.y })); + let w = new Array(n).fill(1); + let M = 0, b1 = 0, b2 = 0; + for (let iter = 0; iter < 8; iter++) { + // weighted normal equations for design [1, cos, sin] (3×3) + const A = [[0, 0, 0], [0, 0, 0], [0, 0, 0]]; + const bv = [0, 0, 0]; + for (let i = 0; i < n; i++) { + const x = [1, rows[i].c, rows[i].s]; + const wi = w[i]; + for (let r = 0; r < 3; r++) { + bv[r] += wi * x[r] * rows[i].y; + for (let cc = 0; cc < 3; cc++) A[r][cc] += wi * x[r] * x[cc]; + } + } + const sol = solve3(A, bv); + if (!sol) return null; + [M, b1, b2] = sol; + // update biweight weights from residuals (c = 4.685·1.4826·MAD) + const res = rows.map((r) => r.y - (M + b1 * r.c + b2 * r.s)); + const absr = res.map(Math.abs).sort((a, b) => a - b); + const mad = absr[absr.length >> 1] || 1; + const cc = 4.685 * 1.4826 * mad; + w = res.map((e) => (Math.abs(e) < cc ? (1 - (e / cc) ** 2) ** 2 : 0)); + } + return { mesor: M, b1, b2, amp: Math.hypot(b1, b2), phi: Math.atan2(b2, b1) }; +} + +/** Solve 3×3 A·x = b by Gaussian elimination with partial pivot. null if singular. */ +function solve3(A: number[][], b: number[]): [number, number, number] | null { + const m = A.map((row, i) => [...row, b[i]]); + for (let col = 0; col < 3; col++) { + let piv = col; + for (let r = col + 1; r < 3; r++) if (Math.abs(m[r][col]) > Math.abs(m[piv][col])) piv = r; + if (Math.abs(m[piv][col]) < 1e-12) return null; + [m[col], m[piv]] = [m[piv], m[col]]; + const pv = m[col][col]; + for (let k = col; k < 4; k++) m[col][k] /= pv; + for (let r = 0; r < 3; r++) { + if (r === col) continue; + const f = m[r][col]; + for (let k = col; k < 4; k++) m[r][k] -= f * m[col][k]; + } + } + return [m[0][3], m[1][3], m[2][3]]; +} + +/** Per-minute RMSSD (ms) from a (raw) RR stream; null if too few beats to be + * meaningful. Cleans physiological + ectopic artifacts first (shared `cleanRr`). */ +function minuteRmssd(rr: number[] | undefined): number | null { + if (!rr || rr.length < 12) return null; + const c = cleanRr(rr); + if (c.length < 10) return null; + let s = 0; + for (let i = 1; i < c.length; i++) { const d = c[i] - c[i - 1]; s += d * d; } + return Math.sqrt(s / (c.length - 1)); +} + +/** Median of finite values; null if none. */ +function medOf(xs: (number | null)[]): number | null { + const a = xs.filter((x): x is number => x != null && Number.isFinite(x)).sort((p, q) => p - q); + return a.length ? a[a.length >> 1] : null; +} + +/** Median filter (±k) to denoise HR before change-point search. */ +function smooth(ys: number[], k: number): number[] { + const n = ys.length; + const out = new Array(n); + for (let i = 0; i < n; i++) { + const lo = Math.max(0, i - k), hi = Math.min(n, i + k + 1); + const seg = ys.slice(lo, hi).sort((a, b) => a - b); + out[i] = seg[seg.length >> 1]; + } + return out; +} + +/** + * Single mean-shift change-point in a window, constrained to a direction ('drop' = + * onset, HR falls). Maximises the SSE reduction of a 2-segment split — i.e. the + * STRONGEST drop, which for the pre-bathyphase window is the evening→sleep onset + * (bigger than any quiet-evening dip). null if the window is too thin. + */ +function changePoint(pts: Pt[], want: 'drop' | 'rise'): number | null { + const MIN = 15; + if (pts.length < 2 * MIN) return null; + const ys = smooth(pts.map((p) => p.y), 5); + const n = ys.length; + const pre = new Array(n + 1).fill(0); + const pre2 = new Array(n + 1).fill(0); + for (let i = 0; i < n; i++) { pre[i + 1] = pre[i] + ys[i]; pre2[i + 1] = pre2[i] + ys[i] * ys[i]; } + const sse = (a: number, b: number): number => { + const cnt = b - a; if (cnt <= 0) return 0; + const sum = pre[b] - pre[a]; + return (pre2[b] - pre2[a]) - (sum * sum) / cnt; + }; + const total = sse(0, n); + let best: { gain: number; tau: number } | null = null; + for (let tau = MIN; tau < n - MIN; tau++) { + const d = (pre[n] - pre[tau]) / (n - tau) - pre[tau] / tau; // +ve = rise + if (want === 'rise' && d <= 0) continue; + if (want === 'drop' && d >= 0) continue; + const gain = total - (sse(0, tau) + sse(tau, n)); + if (!best || gain > best.gain) best = { gain, tau }; + } + return best ? pts[best.tau].t : null; +} + +/** + * Main sleep period. ONSET = the strongest HR drop in [bath−8h, bath] (the + * evening→sleep fall; robust to quiet sedentary evenings that a below-mesor run + * would wrongly absorb). WAKE = the end of the consolidated below-MESOR run + * extending forward from the bathyphase, bridging interior >mesor bumps ≤ BRIDGE + * min (REM/arousal) and ending at the morning rise that persists (daytime). Mesor + * is the cosinor midline — sleep below, wake above — so it doesn't cut through + * high-sleeping-HR users' light sleep. Returns onset/wake, or null. + */ +function mainSleepPeriod(pts: Pt[], bath: number, mesor: number): { onset: number; wake: number } | null { + const n = pts.length; + if (n < 30) return null; + const ys = smooth(pts.map((p) => p.y), 5); + const asleep = ys.map((v) => v < mesor); + // anchor = index nearest the bathyphase + let a = 0; + for (let i = 1; i < n; i++) if (Math.abs(pts[i].t - bath) < Math.abs(pts[a].t - bath)) a = i; + const BRIDGE = 60 * 60; + // wake: walk forward from anchor, bridging >mesor gaps ≤ BRIDGE; stop at a long one. + let end = a; + for (let i = a + 1; i < n;) { + if (asleep[i]) { end = i; i++; continue; } + let k = i; while (k < n && !asleep[k]) k++; + if (pts[(k < n ? k : n) - 1].t - pts[i].t > BRIDGE) break; + i = k; + } + // onset: strongest evening→sleep drop; fall back to the below-mesor run start. + let start = a; + for (let i = a - 1; i >= 0;) { + if (asleep[i]) { start = i; i--; continue; } + let k = i; while (k >= 0 && !asleep[k]) k--; + if (pts[i].t - pts[k + 1].t > BRIDGE) break; + i = k; + } + const onsetCp = changePoint(pts.filter((p) => p.t >= bath - 8 * 3600 && p.t <= bath), 'drop'); + const onset = onsetCp != null && onsetCp >= pts[start].t ? onsetCp : pts[start].t; + return { onset, wake: pts[end].t }; +} + +export interface SleepStaging { + in_bed_min: number; asleep_min: number; efficiency: number; + awake_min: number; light_min: number; deep_min: number; rem_min: number; + hypnogram: { t: number; stage: 'awake' | 'light' | 'deep' | 'rem' }[]; +} + +/** + * The ONE sleep stager (ESTIMATE/beta). Classifies every minute in [onset, wake] + * into awake / light / deep / rem, then returns BOTH the per-minute hypnogram AND the + * reconciled totals — so the graph and the stage breakdown can never disagree (the + * previous bug: two different classifiers). + * + * TWO AXES: + * • HR LEVEL — sets WAKE (sustained ≥20 min elevation, or off-wrist) and the + * overall depth ordering. + * • RR/HRV AUTONOMIC (when per-minute `rr` is present) — separates DEEP from REM, + * which HR level CANNOT do on calm nights (REM-HR ≈ light-HR). Deep/SWS = high + * parasympathetic tone → HIGH beat-to-beat variability (RMSSD) + lowest HR; REM = + * parasympathetic withdrawal → REDUCED RMSSD + mildly elevated HR (verified on real + * nights: corr(HR,RMSSD) ≈ −0.4, and a pure-HR stager read REM ≈ 0%). Thresholds are + * RATIOS to the night's asleep-RMSSD median, so they self-calibrate per user: + * deep = RMSSD ≥ 1.15·median AND HR ≤ asleep-HR-median + * rem = RMSSD ≤ 0.88·median (asleep, below wake) + * light = everything else (the broad majority). + * • Without usable RR we fall back to the legacy HR-only bands (graceful, honest — + * such nights just can't resolve REM, rather than fabricate it). + * asleep = light+deep+rem; efficiency = asleep / in-bed. + */ +export function stageSleep( + minutes: { ts: number; hr_avg: number; rr?: number[] }[], onset: number, wake: number, mesor: number, +): SleepStaging { + const inBed = Math.max(1, Math.round((wake - onset) / 60)); + const win = minutes.filter((m) => m.ts >= onset && m.ts <= wake).sort((a, b) => a.ts - b.ts); + const empty: SleepStaging = { + in_bed_min: inBed, asleep_min: 0, efficiency: 0, awake_min: inBed, + light_min: 0, deep_min: 0, rem_min: 0, hypnogram: [], + }; + const worn = win.filter((m) => m.hr_avg > 0); + if (worn.length < 5) return empty; + const hrs = worn.map((m) => m.hr_avg); + const floor = percentile(hrs, 10) ?? Math.min(...hrs); + const span = Math.max(1, mesor - floor); + // Light-dominant bands (HR-only stages_beta): deep = only the lowest HR, rem = a + // narrow elevated band just below wake, light = the broad middle (most of sleep). + const tAwake = Math.max(floor + 10, floor + 0.70 * span); + // REM = the broad "elevated-but-asleep" band below wake. Anchored at 0.40·span (was + // 0.60, which sat just under tAwake → a 2-3 bpm sliver that almost never triggered, + // so REM read ~0 on calm nights). REM's HR is only modestly above light, not 60% of + // the way to the daytime mesor — 0.40 puts the band where REM actually lives. + const tRem = floor + 0.40 * span; + const tDeep = floor + 0.12 * span; + + // smoothed HR aligned to `win` (off-wrist minutes carry hr 0 → forced awake) + const ys = smooth(win.map((m) => (m.hr_avg > 0 ? m.hr_avg : tAwake + 50)), 5); + + // ── Autonomic axis: per-minute RMSSD (smoothed ±2), and the night's asleep + // RMSSD/HR medians as self-calibrating references. Only used if a real fraction of + // the asleep minutes carry enough RR — else deep/rem stay on the HR-only fallback. ── + const rmRaw = win.map((m) => minuteRmssd(m.rr)); + const rmS = rmRaw.map((_, i) => { + const seg: (number | null)[] = []; + for (let j = Math.max(0, i - 2); j < Math.min(win.length, i + 3); j++) seg.push(rmRaw[j]); + return medOf(seg); + }); + const asleepI = win.map((_, i) => i).filter((i) => win[i].hr_avg > 0 && ys[i] < tAwake); + const rmRef = medOf(asleepI.map((i) => rmS[i])); + const hrRef = medOf(asleepI.map((i) => ys[i])); + const rrUsable = rmRef != null && hrRef != null + && asleepI.filter((i) => rmS[i] != null).length >= Math.max(20, Math.floor(asleepI.length * 0.4)); + const DEEP_R = 1.15, REM_R = 0.88; // RMSSD ratios vs asleep median (tuned on real RR: REM≈21%, deep≈13%) + + const stage: ('awake' | 'light' | 'deep' | 'rem')[] = new Array(win.length).fill('light'); + // pass 1: provisional per-minute classes + for (let k = 0; k < win.length; k++) { + if (win[k].hr_avg <= 0) { stage[k] = 'awake'; continue; } + const v = ys[k]; + if (v >= tAwake) { stage[k] = 'awake'; continue; } + if (rrUsable && rmS[k] != null) { + // RR present → autonomic deep/rem split (the REM-detecting path). + const rm = rmS[k]!; + stage[k] = (rm >= DEEP_R * rmRef! && v <= hrRef!) ? 'deep' + : (rm <= REM_R * rmRef!) ? 'rem' + : 'light'; + } else { + // No usable RR this minute/night → legacy HR-level bands. + stage[k] = v < tDeep ? 'deep' : v >= tRem ? 'rem' : 'light'; + } + } + // pass 2: an "awake" minute only counts as awake if part of a SUSTAINED (≥20 min) + // run; otherwise it's a brief REM/arousal bump → reclassify as rem. + let k = 0; + while (k < win.length) { + if (stage[k] === 'awake' && win[k].hr_avg > 0) { + let j = k; while (j < win.length && stage[j] === 'awake' && win[j].hr_avg > 0) j++; + if ((win[j - 1].ts - win[k].ts) / 60 < 20) for (let x = k; x < j; x++) stage[x] = 'rem'; + k = j; + } else k++; + } + // pass 3: BOUT-SMOOTHING — the per-minute classes flicker when HR hovers near a + // threshold (sawtooth hypnogram). Merge any run shorter than MIN_BOUT minutes into + // its longer neighbour, repeatedly, so the hypnogram reads as stable sleep bouts + // (real stages last many minutes). Awake keeps a higher floor so a genuine short + // awakening survives; sub-floor awakes are noise → absorbed into surrounding sleep. + const MIN_BOUT = 6, MIN_AWAKE_BOUT = 10; + for (let iter = 0; iter < 4; iter++) { + const runs: { s: number; e: number }[] = []; + for (let i = 0; i < win.length;) { let j = i; while (j < win.length && stage[j] === stage[i]) j++; runs.push({ s: i, e: j - 1 }); i = j; } + if (runs.length <= 1) break; + let changed = false; + for (let r = 0; r < runs.length; r++) { + const { s, e } = runs[r]; + const lenMin = e - s + 1; + const floorMin = stage[s] === 'awake' ? MIN_AWAKE_BOUT : MIN_BOUT; + if (lenMin >= floorMin) continue; + const prev = r > 0 ? runs[r - 1] : null; + const next = r < runs.length - 1 ? runs[r + 1] : null; + let target: typeof stage[number] | null = null; + if (prev && next) target = (prev.e - prev.s) >= (next.e - next.s) ? stage[prev.s] : stage[next.s]; + else if (prev) target = stage[prev.s]; + else if (next) target = stage[next.s]; + if (target) { for (let x = s; x <= e; x++) stage[x] = target; changed = true; } + } + if (!changed) break; + } + let light = 0, deep = 0, rem = 0, awake = 0; + for (const s of stage) { if (s === 'awake') awake++; else if (s === 'deep') deep++; else if (s === 'rem') rem++; else light++; } + const asleep = light + deep + rem; + return { + in_bed_min: inBed, asleep_min: asleep, efficiency: clamp(asleep / inBed, 0, 1), + awake_min: awake, light_min: light, deep_min: deep, rem_min: rem, + hypnogram: win.map((m, idx) => ({ t: m.ts, stage: stage[idx] })), + }; +} + +export interface CircadianOpts { + now?: number; // unix s "now" (default = latest minute ts); for determinism + settleSec?: number; // a wake older than this = night complete (default 600s) + anchorTs?: number; // target a specific cycle's bathyphase nearest this ts +} + +/** + * calcCircadian(minutes, opts) + * + * Fits the circadian cosinor over the HR series, then for the most-recent + * completed cycle (bathyphase = HR trough) finds main-sleep onset (HR drop in + * the 8h before the trough) and wake (HR rise in the 8h after). Returns the + * circadian phase + that boundary. + * + * Confidence = rhythm strength (amplitude) × whether a clean onset+wake pair was + * found. Abstains (nulls, confidence 0) when the rhythm is too weak or data too + * thin — never fabricates a boundary. + */ +export function calcCircadian(minutes: Minute[], opts: CircadianOpts = {}): Metric { + const usable = minutes.filter(isHrUsable).map((m) => ({ t: m.ts, y: m.hr_avg })).sort((a, b) => a.t - b.t); + const settle = opts.settleSec ?? 600; + + const empty = (): Metric => ({ + mesor: null, amplitude: null, acrophase_ts: null, bathyphase_ts: null, + onset_ts: null, wake_ts: null, in_bed_min: 0, settled: false, + confidence: 0, tier: 'HIGH', inputs_used: [], + }); + if (usable.length < 120) return empty(); + + const now = opts.now ?? usable[usable.length - 1].t; + const fit = fitCosinor(usable); + if (!fit) return empty(); + // No detectable circadian rhythm (flat HR) → abstain, never fabricate a boundary. + if (fit.amp < 3) return empty(); + + // bathyphase (HR trough): amp·cos(ωt − φ) minimal → ωt − φ = π. + const bathBase = (fit.phi + Math.PI) / W; + // acrophase (HR peak): ωt − φ = 0. + const acroBase = fit.phi / W; + const nearest = (base: number, ref: number) => base + Math.round((ref - base) / DAY) * DAY; + + // Pick the cycle: the most-recent bathyphase IN THE PAST. The wake is found within + // [bath, bath+8h] clipped to available data, and `settled` (wake older than the + // settle window) decides whether the night is complete. We deliberately do NOT + // require bath+8h to be in the past — that wrongly stepped back a whole night for + // an early riser checking soon after waking (the "previous night / short sleep" bug). + let bath = nearest(bathBase, opts.anchorTs ?? now); + if (bath > now - 3600) bath -= DAY; // bathyphase must be ~in the past + const acro = nearest(acroBase, now); + + const inWin = (lo: number, hi: number) => usable.filter((p) => p.t >= lo && p.t <= hi); + // Consolidated main-sleep period around the bathyphase (handles REM bumps + daytime). + const period = mainSleepPeriod(inWin(bath - 8 * 3600, bath + 10 * 3600), bath, fit.mesor); + const onset_ts = period ? period.onset : null; + const wake_ts = period ? period.wake : null; + const in_bed_min = onset_ts != null && wake_ts != null ? Math.round((wake_ts - onset_ts) / 60) : 0; + const settled = wake_ts != null && wake_ts <= now - settle; + + // confidence: rhythm strength (amp 2→0 … 10→1) gated on a clean onset+wake pair. + const rhythm = clamp((fit.amp - 2) / 8, 0, 1); + const paired = onset_ts != null && wake_ts != null ? 1 : 0.3; + const confidence = round(rhythm * paired, 2); + + return { + mesor: round(fit.mesor, 1), + amplitude: round(fit.amp, 1), + acrophase_ts: acro, + bathyphase_ts: bath, + onset_ts, + wake_ts, + in_bed_min, + settled, + confidence, + tier: 'HIGH', + inputs_used: ['hr'], + }; +} diff --git a/src/cycle.ts b/src/cycle.ts new file mode 100644 index 0000000..83bd122 --- /dev/null +++ b/src/cycle.ts @@ -0,0 +1,100 @@ +// cycle.ts — menstrual cycle estimation. +// +// LOG-ANCHORED + calendar method: the user logs period-start dates; prediction is +// driven by those, never inferred from biometrics alone (we don't claim to detect +// ovulation from temperature — that would be fabrication). The luteal phase is +// physiologically stable at ~14 days (Wilcox 2000), so ovulation ≈ next-period − 14 +// and the fertile window is the 5 days before ovulation + ovulation day. +// +// Biometric overlay (skin-temp / RHR / HRV shifts across the cycle) is rendered by +// the caller from stored `daily` values — it ENRICHES the view but is descriptive, +// never the basis of the prediction. Honest: an ESTIMATE, not medical or +// contraceptive guidance. + +import type { Metric } from './types' +import { median } from './util' + +const DAY_MS = 86400000 +const toMs = (d: string) => Date.parse(d + 'T00:00:00Z') +const toDate = (ms: number) => new Date(ms).toISOString().slice(0, 10) +const daysBetween = (a: string, b: string) => Math.round((toMs(b) - toMs(a)) / DAY_MS) + +export type CyclePhase = 'menstruation' | 'follicular' | 'ovulation' | 'luteal' | 'unknown' + +export interface CycleValue { + cycle_day: number | null // 1-based day within the current cycle (day 1 = start) + phase: CyclePhase + mean_length: number | null // median observed cycle length (days); null if <2 starts + length_history: number[] // observed consecutive-start gaps (days) + last_start: string | null // most recent logged period start (YYYY-MM-DD) + predicted_next: string | null // predicted next period start + days_until_next: number | null + ovulation_est: string | null // predicted_next − 14d + fertile_start: string | null // ovulation − 5d + fertile_end: string | null // ovulation + 1d + note: string +} + +const DEFAULT_LEN = 28 // population default until the user has ≥2 logged periods +const LUTEAL = 14 // stable luteal length → ovulation = next_period − LUTEAL +const MENSES = 5 // assumed menses length when no explicit end is logged + +/** Estimate the current cycle position + next-period / fertile-window prediction + * from a list of logged period-START dates. `today` is supplied (pure, no clock). */ +export function calcCycle(startsRaw: string[], today: string): Metric { + const empty = (note: string): Metric => ({ + cycle_day: null, phase: 'unknown', mean_length: null, length_history: [], + last_start: null, predicted_next: null, days_until_next: null, + ovulation_est: null, fertile_start: null, fertile_end: null, note, + confidence: 0, tier: 'ESTIMATE', inputs_used: ['period_log'], + }) + + const starts = Array.from(new Set(startsRaw)) + .filter((d) => /^\d{4}-\d{2}-\d{2}$/.test(d) && toMs(d) <= toMs(today)) + .sort() + if (starts.length === 0) return empty('Log a period to start tracking your cycle.') + + // Observed cycle lengths between consecutive starts (keep physiological 15–60d). + const lengths: number[] = [] + for (let i = 1; i < starts.length; i++) { + const len = daysBetween(starts[i - 1], starts[i]) + if (len >= 15 && len <= 60) lengths.push(len) + } + const med = lengths.length ? median(lengths) : null + const meanLen = med == null ? null : Math.round(med) + const useLen = meanLen ?? DEFAULT_LEN + + const last = starts[starts.length - 1] + const cycleDay = daysBetween(last, today) + 1 // day 1 = the start date itself + + const nextMs = toMs(last) + useLen * DAY_MS + const predictedNext = toDate(nextMs) + const daysUntil = daysBetween(today, predictedNext) + const ovMs = nextMs - LUTEAL * DAY_MS + const ovulation = toDate(ovMs) + const fertileStart = toDate(ovMs - 5 * DAY_MS) + const fertileEnd = toDate(ovMs + 1 * DAY_MS) + + // Phase by calendar method. + const todayMs = toMs(today) + let phase: CyclePhase + if (cycleDay <= MENSES) phase = 'menstruation' + else if (todayMs >= toMs(fertileStart) && todayMs <= toMs(fertileEnd)) phase = 'ovulation' + else if (todayMs < ovMs) phase = 'follicular' + else phase = 'luteal' + + // Confidence grows with the number of observed cycles; collapses if very overdue + // (a missed/late period makes the calendar prediction unreliable — say so). + let conf = lengths.length === 0 ? 0.3 : Math.min(0.9, 0.4 + 0.15 * lengths.length) + if (cycleDay > useLen * 1.6) { phase = 'unknown'; conf = Math.min(conf, 0.2) } + + return { + cycle_day: cycleDay, phase, mean_length: meanLen, length_history: lengths, + last_start: last, predicted_next: predictedNext, days_until_next: daysUntil, + ovulation_est: ovulation, fertile_start: fertileStart, fertile_end: fertileEnd, + note: lengths.length === 0 + ? 'Based on one logged period and a 28-day default — accuracy improves as you log more.' + : `Based on ${lengths.length + 1} logged periods (median ${useLen}-day cycle).`, + confidence: conf, tier: 'ESTIMATE', inputs_used: ['period_log'], + } +} diff --git a/src/cycles.ts b/src/cycles.ts new file mode 100644 index 0000000..042d5ad --- /dev/null +++ b/src/cycles.ts @@ -0,0 +1,108 @@ +// §Sleep cycles — ultradian NREM↔REM cycle detection, adapted from Rosenblum et al. +// 2024 (eLife 13:RP96784, "Fractal cycles of sleep"). +// +// The paper detects sleep cycles as the interval between successive PEAKS of the +// smoothed, z-normalized EEG *aperiodic (fractal) slope* time series — MATLAB +// findpeaks(MinPeakDistance = 20 min, MinPeakProminence = 0.9 z). Peaks coincide with +// REM, troughs with non-REM; ~4–6 cycles/night, ~90 min each. It's a continuous, +// parameter-light definition with NO arbitrary per-night knob — exactly the property +// we want. +// +// We have no EEG. But heart-rate VARIABILITY carries the SAME ultradian rhythm — RMSSD +// is high in deep non-REM and low in REM — so we run the IDENTICAL peak algorithm on +// the smoothed z-normalized per-minute RMSSD series. Cycle boundaries anchor on the +// deep-sleep RMSSD peaks (the paper anchors on REM peaks: a cosmetic half-cycle phase +// shift — the cycle COUNT and DURATION, the actual outputs, are unchanged). +// +// Pure & deterministic. Tier ESTIMATE (it's a wrist-HRV proxy for an EEG method). +import { cleanRr } from './hrv'; + +export interface SleepCycle { start_ts: number; end_ts: number; duration_min: number } +export interface SleepCyclesValue { + cycles: SleepCycle[]; + mean_duration_min: number | null; + n: number; + series: { t: number; z: number }[]; // smoothed z-RMSSD, for plotting under the hypnogram +} + +const SMOOTH_MIN = 10; // ±10-min moving average → exposes the ~90-min envelope +const MIN_PEAK_DIST = 20; // min minutes between cycle boundaries (paper: 20 min / 40×30s) +const MIN_PROMINENCE = 0.9; // z — paper's MinPeakProminence (the |0.9| z descent gate) + +/** Per-minute RMSSD (ms) from a raw RR stream; null if too few clean beats. */ +function minuteRmssd(rr: number[] | undefined): number | null { + if (!rr || rr.length < 12) return null; + const c = cleanRr(rr); + if (c.length < 10) return null; + let s = 0; + for (let i = 1; i < c.length; i++) { const d = c[i] - c[i - 1]; s += d * d; } + return Math.sqrt(s / (c.length - 1)); +} + +/** Local maxima with topographic prominence ≥ minProm, then enforce a minimum spacing + * (keep the highest in each cluster) — a faithful port of MATLAB findpeaks' two gates. */ +function findPeaks(y: (number | null)[], minDist: number, minProm: number): number[] { + const n = y.length; + const cand: { i: number; v: number }[] = []; + for (let i = 1; i < n - 1; i++) { + const yi = y[i]; if (yi == null) continue; + const a = y[i - 1] ?? -Infinity, b = y[i + 1] ?? -Infinity; + if (!(yi >= a && yi > b)) continue; + // prominence: walk out each side until a higher sample; the higher of the two + // in-between minima is the reference. prom = peak − that reference. + let l = i; while (l > 0 && (y[l - 1] ?? -Infinity) < yi) l--; + let r = i; while (r < n - 1 && (y[r + 1] ?? -Infinity) < yi) r++; + let lmin = yi, rmin = yi; + for (let k = l; k <= i; k++) { const v = y[k]; if (v != null && v < lmin) lmin = v; } + for (let k = i; k <= r; k++) { const v = y[k]; if (v != null && v < rmin) rmin = v; } + if (yi - Math.max(lmin, rmin) >= minProm) cand.push({ i, v: yi }); + } + cand.sort((p, q) => q.v - p.v); // tallest first + const kept: number[] = []; + for (const c of cand) if (kept.every((k) => Math.abs(c.i - k) >= minDist)) kept.push(c.i); + return kept.sort((a, b) => a - b); +} + +/** + * detectSleepCycles(minutes, onset, wake) + * minutes: per-minute records over (at least) the sleep window, each with optional `rr`. + * Returns the ultradian cycles (peak-to-peak), their mean duration, and the smoothed + * z-RMSSD series. Abstains (empty) when there isn't enough clean RR to resolve cycles. + */ +export function detectSleepCycles( + minutes: { ts: number; rr?: number[] }[], onset: number, wake: number, +): SleepCyclesValue { + const none: SleepCyclesValue = { cycles: [], mean_duration_min: null, n: 0, series: [] }; + const win = minutes.filter((m) => m.ts >= onset && m.ts <= wake).sort((a, b) => a.ts - b.ts); + if (win.length < 60) return none; // need ~1 h to resolve even one cycle + + const raw = win.map((m) => minuteRmssd(m.rr)); + // ±SMOOTH_MIN moving average (ignoring gaps). + const sm = raw.map((_, i) => { + let s = 0, c = 0; + for (let j = Math.max(0, i - SMOOTH_MIN); j <= Math.min(raw.length - 1, i + SMOOTH_MIN); j++) { + const v = raw[j]; if (v != null) { s += v; c++; } + } + return c ? s / c : null; + }); + const vals = sm.filter((x): x is number => x != null); + if (vals.length < 60) return none; + const mean = vals.reduce((a, b) => a + b, 0) / vals.length; + const sd = Math.sqrt(vals.reduce((a, b) => a + (b - mean) ** 2, 0) / vals.length) || 1; + const z = sm.map((x) => (x == null ? null : (x - mean) / sd)); + + const peaks = findPeaks(z, MIN_PEAK_DIST, MIN_PROMINENCE); + const cycles: SleepCycle[] = []; + for (let i = 0; i + 1 < peaks.length; i++) { + const start_ts = win[peaks[i]].ts, end_ts = win[peaks[i + 1]].ts; + cycles.push({ start_ts, end_ts, duration_min: Math.round((end_ts - start_ts) / 60) }); + } + const mean_duration_min = cycles.length + ? Math.round(cycles.reduce((s, c) => s + c.duration_min, 0) / cycles.length) : null; + const series = win + .map((m, i) => ({ t: m.ts, z: z[i] })) + .filter((p): p is { t: number; z: number } => p.z != null) + .map((p) => ({ t: p.t, z: Math.round(p.z * 1000) / 1000 })); + + return { cycles, mean_duration_min, n: cycles.length, series }; +} diff --git a/src/har.ts b/src/har.ts new file mode 100644 index 0000000..642527a --- /dev/null +++ b/src/har.ts @@ -0,0 +1,325 @@ +// §HAR — Human Activity Recognition from wrist accelerometer windows. +// Method: Mannini et al., Med Sci Sports Exerc 2013 ("Activity recognition using +// a single accelerometer placed at the wrist or ankle"), wrist configuration. +// Per ~4 s window we extract Mannini's feature set (time-domain + frequency-domain +// + db10 wavelet energies) and classify with TRANSPARENT cadence/power thresholds +// (ESTIMATE tier) — no trained model yet (we have no labelled data). The feature +// vector is the same one a future SVM would consume, so the upgrade path is clean. +// +// HONEST SCOPE: this only ever runs on LIVE high-rate stream data (R10/0x33, +// ~100 Hz) — flash-drained history is 1 Hz and CANNOT be motion-classified. +// Pure + synchronous (no I/O, no wasm): a db10 DWT on a 512-sample window is +// microseconds in plain TS; a wasm dependency would only puncture the pure model. + +import { mean, stddev } from './util'; +import type { ActivityClass } from './types'; + +export type { ActivityClass }; + +export interface HarFeatures { + smv_mean: number; + smv_std: number; + smv_min: number; + smv_max: number; + total_power: number; // spectral power in 0.3–15 Hz + dom1_freq: number; // strongest spectral peak in 0.3–15 Hz + dom1_pow: number; + dom2_freq: number; // second peak + dom2_pow: number; + cad_freq: number; // dominant peak inside the locomotion band 0.6–2.5 Hz + cad_pow: number; + dom1_ratio: number; // dom1_pow / total_power (spectral peakiness → periodicity) + freq_ratio_prev: number; // dom1_freq / previous window's dom1_freq (transition cue) + wav_e5: number; // db10 detail energy at level 5 + wav_e6: number; // db10 detail energy at level 6 +} + +// ── db10 decomposition low-pass coefficients (PyWavelets `pywt.Wavelet('db10').dec_lo`). +// Validated in tests via Σh=√2 and Σh²=1. Source: wavelets.pybytes.com / PyWavelets. +export const DB10_LO: number[] = [ + 2.667005790055555358661744877130858277192498290851289932779975e-02, + 1.881768000776914890208929736790939942702546758640393484348595e-01, + 5.272011889317255864817448279595081924981402680840223445318549e-01, + 6.884590394536035657418717825492358539771364042407339537279681e-01, + 2.811723436605774607487269984455892876243888859026150413831543e-01, + -2.498464243273153794161018979207791000564669737132073715013121e-01, + -1.959462743773770435042992543190981318766776476382778474396781e-01, + 1.273693403357932600826772332014009770786177480422245995563097e-01, + 9.305736460357235116035228983545273226942917998946925868063974e-02, + -7.139414716639708714533609307605064767292611983702150917523756e-02, + -2.945753682187581285828323760141839199388200516064948779769654e-02, + 3.321267405934100173976365318215912897978337413267096043323351e-02, + 3.606553566956169655423291417133403299517350518618994762730612e-03, + -1.073317548333057504431811410651364448111548781143923213370333e-02, + 1.395351747052901165789318447957707567660542855688552426721117e-03, + 1.992405295185056117158742242640643211762555365514105280067936e-03, + -6.858566949597116265613709819265714196625043336786920516211903e-04, + -1.164668551292854509514809710258991891527461854347597362819235e-04, + 9.358867032006959133405013034222854399688456215297276443521873e-05, + -1.326420289452124481243667531226683305749240960605829756400674e-05, +]; + +/** Quadrature-mirror high-pass: g[k] = (−1)^k · h[N−1−k]. */ +function db10Hi(): number[] { + const N = DB10_LO.length; + return DB10_LO.map((_, k) => (k % 2 === 0 ? 1 : -1) * DB10_LO[N - 1 - k]); +} +const DB10_HI = db10Hi(); + +/** One DWT level with periodic ('wrap') boundary → {approx, detail} (downsampled by 2). */ +function dwtStep(sig: number[], lo: number[], hi: number[]): { a: number[]; d: number[] } { + const n = sig.length, L = lo.length; + const half = Math.floor(n / 2); + const a = new Array(half).fill(0); + const d = new Array(half).fill(0); + for (let i = 0; i < half; i++) { + let sa = 0, sd = 0; + for (let k = 0; k < L; k++) { + const idx = (2 * i + k) % n; // periodic extension + sa += lo[k] * sig[idx]; + sd += hi[k] * sig[idx]; + } + a[i] = sa; d[i] = sd; + } + return { a, d }; +} + +/** Detail-coefficient energy (Σ d²) at each level 1..levels via db10 DWT. */ +export function dwtDetailEnergies(signal: number[], levels: number): number[] { + let a = signal.slice(); + const out: number[] = []; + for (let lvl = 1; lvl <= levels; lvl++) { + if (a.length < 2) { out.push(0); continue; } + const { a: na, d } = dwtStep(a, DB10_LO, DB10_HI); + out.push(d.reduce((s, v) => s + v * v, 0)); + a = na; + } + return out; +} + +// ── 4th-order Butterworth low-pass = two RBJ biquad sections (Butterworth Qs). ── +function biquadLP(sig: number[], fs: number, fc: number, q: number): number[] { + const w0 = (2 * Math.PI * fc) / fs; + const cosw = Math.cos(w0), sinw = Math.sin(w0); + const alpha = sinw / (2 * q); + const a0 = 1 + alpha; + const b0 = ((1 - cosw) / 2) / a0, b1 = (1 - cosw) / a0, b2 = ((1 - cosw) / 2) / a0; + const a1 = (-2 * cosw) / a0, a2 = (1 - alpha) / a0; + const out = new Array(sig.length); + // Initialise to steady state at the first sample so the ~1 g DC gravity offset + // doesn't ring up a 0→1 startup transient (which would inflate power/variance). + const s0 = sig.length ? sig[0] : 0; + let x1 = s0, x2 = s0, y1 = s0, y2 = s0; + for (let i = 0; i < sig.length; i++) { + const x0 = sig[i]; + const y0 = b0 * x0 + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2; + x2 = x1; x1 = x0; y2 = y1; y1 = y0; out[i] = y0; + } + return out; +} +/** 15 Hz 4th-order Butterworth low-pass (Mannini preprocessing). */ +function butterLP4(sig: number[], fs: number, fc = 15): number[] { + return biquadLP(biquadLP(sig, fs, fc, 0.54119610), fs, fc, 1.30656296); +} + +function nextPow2(n: number): number { let p = 1; while (p < n) p <<= 1; return p; } + +/** Iterative radix-2 FFT power spectrum (length N/2+1), N = next pow2 of input. */ +function powerSpectrum(sig: number[]): number[] { + const N = nextPow2(sig.length); + const re = new Float64Array(N), im = new Float64Array(N); + for (let i = 0; i < sig.length; i++) re[i] = sig[i]; + // bit-reversal permutation + for (let i = 1, j = 0; i < N; i++) { + let bit = N >> 1; + for (; j & bit; bit >>= 1) j ^= bit; + j ^= bit; + if (i < j) { const tr = re[i]; re[i] = re[j]; re[j] = tr; const ti = im[i]; im[i] = im[j]; im[j] = ti; } + } + for (let len = 2; len <= N; len <<= 1) { + const ang = (-2 * Math.PI) / len; + const wr = Math.cos(ang), wi = Math.sin(ang); + for (let i = 0; i < N; i += len) { + let cr = 1, ci = 0; + for (let k = 0; k < len / 2; k++) { + const ur = re[i + k], ui = im[i + k]; + const vr = re[i + k + len / 2] * cr - im[i + k + len / 2] * ci; + const vi = re[i + k + len / 2] * ci + im[i + k + len / 2] * cr; + re[i + k] = ur + vr; im[i + k] = ui + vi; + re[i + k + len / 2] = ur - vr; im[i + k + len / 2] = ui - vi; + const ncr = cr * wr - ci * wi; ci = cr * wi + ci * wr; cr = ncr; + } + } + } + const half = N / 2; + const pow = new Array(half + 1); + for (let i = 0; i <= half; i++) pow[i] = (re[i] * re[i] + im[i] * im[i]) / N; + return pow; +} + +/** + * extractHarFeaturesFromSmv(smvRaw, fs, prevDomFreq?) — Mannini feature set from the + * accel-vector MAGNITUDE signal (= SMV). This is the primary entry point: the band's + * decoders give us per-sample |accel| (frameAccel `mags`), which IS the SMV, so the + * ingest path feeds it directly. The 15 Hz Butterworth LP is applied to the magnitude. + */ +export function extractHarFeaturesFromSmv(smvRaw: number[], fs: number, prevDomFreq = 0): HarFeatures { + const smv = butterLP4(smvRaw, fs); // 15 Hz LP on the magnitude + const smv_mean = mean(smv), smv_std = stddev(smv); + const smv_min = Math.min(...smv), smv_max = Math.max(...smv); + + // Spectrum of the DC-removed SMV. + const ac = smv.map((v) => v - smv_mean); + const pow = powerSpectrum(ac); + const N = nextPow2(ac.length); + const binHz = fs / N; + const idxOf = (f: number) => Math.round(f / binHz); + const loBin = Math.max(1, idxOf(0.3)), hiBin = Math.min(pow.length - 1, idxOf(15)); + let total = 0; + for (let i = loBin; i <= hiBin; i++) total += pow[i]; + + // top-2 peaks in 0.3–15 Hz + let d1i = loBin, d2i = loBin; + for (let i = loBin; i <= hiBin; i++) if (pow[i] > pow[d1i]) d1i = i; + for (let i = loBin; i <= hiBin; i++) if (i !== d1i && pow[i] > pow[d2i]) d2i = i; + // dominant peak in the locomotion band 0.6–2.5 Hz + const cLo = Math.max(1, idxOf(0.6)), cHi = Math.min(pow.length - 1, idxOf(2.5)); + let ci = cLo; + for (let i = cLo; i <= cHi; i++) if (pow[i] > pow[ci]) ci = i; + + const dom1_freq = d1i * binHz, dom1_pow = pow[d1i]; + const wav = dwtDetailEnergies(smv, 6); + + return { + smv_mean, smv_std, smv_min, smv_max, + total_power: total, + dom1_freq, dom1_pow, + dom2_freq: d2i * binHz, dom2_pow: pow[d2i], + cad_freq: ci * binHz, cad_pow: pow[ci], + dom1_ratio: total > 0 ? dom1_pow / total : 0, + freq_ratio_prev: prevDomFreq > 0 ? dom1_freq / prevDomFreq : 1, + wav_e5: wav[4] ?? 0, + wav_e6: wav[5] ?? 0, + }; +} + +/** + * extractHarFeatures(x, y, z, fs, prevDomFreq?) — convenience wrapper for tri-axial + * input (mainly tests/synthetic data): computes the SMV magnitude and delegates. + */ +export function extractHarFeatures( + x: number[], y: number[], z: number[], fs: number, prevDomFreq = 0, +): HarFeatures { + const n = Math.min(x.length, y.length, z.length); + const smvRaw = new Array(n); + for (let i = 0; i < n; i++) smvRaw[i] = Math.sqrt(x[i] * x[i] + y[i] * y[i] + z[i] * z[i]); + return extractHarFeaturesFromSmv(smvRaw, fs, prevDomFreq); +} + +// Tunable thresholds (documented, ESTIMATE tier — replace with trained model later). +// We classify off dom1_freq (strongest 0.3–15 Hz peak = locomotion fundamental); the +// Mannini cad_freq feature (0.6–2.5 Hz band, capped) is kept for the future SVM. +const SED_POWER = 0.02; // below this spectral power → not moving rhythmically +const SED_STD = 0.04; // g +const PERIODIC = 0.25; // dom1_ratio above this → clearly periodic (walk/run/cycle) +const RUN_HZ = 2.4; // dominant ≥2.4 Hz → running cadence +const WALK_HZ = 1.3; // 1.3–2.4 Hz → walking cadence +const CYCLE_HZ_LO = 0.6; // 0.6–1.3 Hz with smooth motion → cycling pedal cadence + +/** + * classifyActivityWindow(f) — transparent threshold classifier over Mannini features. + * Returns the class + a confidence from spectral peakiness & band power. ESTIMATE tier. + */ +export function classifyActivityWindow(f: HarFeatures): { cls: ActivityClass; confidence: number } { + // Sedentary: low spectral power AND low motion variance. + if (f.total_power < SED_POWER && f.smv_std < SED_STD) { + return { cls: 'sedentary', confidence: 0.6 }; + } + const periodic = f.dom1_ratio >= PERIODIC; + const peakConf = Math.min(0.9, 0.4 + f.dom1_ratio); // peakier spectrum → higher confidence + const cad = f.dom1_freq; + + if (periodic) { + if (cad >= RUN_HZ) return { cls: 'run', confidence: peakConf }; + if (cad >= WALK_HZ) return { cls: 'walk', confidence: peakConf }; + if (cad >= CYCLE_HZ_LO && f.smv_std < SED_STD * 4) + return { cls: 'cycle', confidence: peakConf * 0.9 }; // smooth low-cadence motion + } + // Elevated but non-periodic motion (irregular, no clean cadence) → resistance/lifting. + if (f.smv_std >= SED_STD && !periodic) return { cls: 'lift', confidence: 0.45 }; + return { cls: 'other', confidence: 0.4 }; +} + +// ── Segmentation: per-window votes → smoothed phases (graceful activity switches) ── +export interface ClassVote { ts: number; cls: ActivityClass; conf: number } +export interface WorkoutSegment { start_ts: number; end_ts: number; type: ActivityClass; confidence: number } +export interface SegmentResult { primary: ActivityClass; segments: WorkoutSegment[]; type_confidence: number } + +/** Mode of a class array (ties → highest mean confidence handled by caller). */ +function modeClass(window: ClassVote[]): ActivityClass { + const c: Record = {}; + for (const v of window) c[v.cls] = (c[v.cls] ?? 0) + 1; + let best: ActivityClass = window[0].cls, bestN = -1; + for (const k of Object.keys(c)) if (c[k] > bestN) { bestN = c[k]; best = k as ActivityClass; } + return best; +} + +/** + * segmentWorkout(votes, opts) — smooth the per-window class timeline (median filter to + * kill blips) and run-length-encode into phases ≥ minPhaseSec, so a run→cycle or + * lift→cardio switch becomes ONE workout with labelled phases (hysteresis, not flip-flop). + * primary = longest phase; 'mixed' when no phase dominates (circuit/brick). + */ +export function segmentWorkout( + votes: ClassVote[], + opts: { smoothWin?: number; minPhaseSec?: number } = {}, +): SegmentResult { + const smoothWin = opts.smoothWin ?? 7; // ~window count for median smoothing + const minPhaseSec = opts.minPhaseSec ?? 180; // a phase must persist ≥3 min + if (votes.length === 0) return { primary: 'other', segments: [], type_confidence: 0 }; + const sorted = [...votes].sort((a, b) => a.ts - b.ts); + + // Median-class smoothing over a sliding window (hysteresis against momentary blips). + const smoothed: ClassVote[] = sorted.map((v, i) => { + const half = Math.floor(smoothWin / 2); + const win = sorted.slice(Math.max(0, i - half), Math.min(sorted.length, i + half + 1)); + return { ts: v.ts, cls: modeClass(win), conf: v.conf }; + }); + + // Run-length encode into raw phases. + const raw: WorkoutSegment[] = []; + for (const v of smoothed) { + const last = raw[raw.length - 1]; + if (last && last.type === v.cls) { + last.end_ts = v.ts; + last.confidence = (last.confidence + v.conf) / 2; + } else { + raw.push({ start_ts: v.ts, end_ts: v.ts, type: v.cls, confidence: v.conf }); + } + } + + // Merge phases shorter than minPhaseSec into the neighbour (drop noise blips). + const phases: WorkoutSegment[] = []; + for (const seg of raw) { + const dur = seg.end_ts - seg.start_ts; + if (dur < minPhaseSec && phases.length > 0) { + phases[phases.length - 1].end_ts = seg.end_ts; // absorb blip into prior phase + } else if (dur < minPhaseSec && phases.length === 0) { + phases.push({ ...seg }); // keep until a real phase appears + } else { + if (phases.length && phases[phases.length - 1].type === seg.type) + phases[phases.length - 1].end_ts = seg.end_ts; + else phases.push({ ...seg }); + } + } + + // Primary = longest-duration phase; 'mixed' if the top phase is <50% of total. + const totalDur = phases.reduce((s, p) => s + (p.end_ts - p.start_ts), 0) || 1; + let top = phases[0]; + for (const p of phases) if ((p.end_ts - p.start_ts) > (top.end_ts - top.start_ts)) top = p; + const topShare = (top.end_ts - top.start_ts) / totalDur; + const primary: ActivityClass = topShare >= 0.5 ? top.type : 'other'; + const type_confidence = Math.round(top.confidence * topShare * 100) / 100; + + return { primary, segments: phases, type_confidence }; +} diff --git a/src/hrv.ts b/src/hrv.ts index 84cb222..3d3e04d 100644 --- a/src/hrv.ts +++ b/src/hrv.ts @@ -14,7 +14,7 @@ // RR is decoded from type-24 records (parse_r24 rr_intervals_ms) and validated on // real hardware (99.7% physiological; p50≈860 ms ≈ 70 bpm). This is NOT the parked // R17/R11 live-optical path — it's the historical RR that ships in every record. -import { round, mean, stddev } from './util'; +import { round, mean, stddev, median } from './util'; import type { Metric, HrvStabilityValue, IrregularValue } from './types'; /** Standard frequency bands (Hz) per the HRV Task Force (1996). */ @@ -126,26 +126,34 @@ export function freqDomainHrv(rrRaw: number[]): FreqDomainHrv { const mean = rr.reduce((a, b) => a + b, 0) / rr.length; const x = rr.map((r) => r - mean); const span = t[t.length - 1] - t[0]; - if (span < 60) return none; // need ≥1 min of beats for stable spectral estimates - + // Task Force 1996 rule: a band needs ≥10× the wavelength of its lower bound. + // HF lower bound 0.15 Hz → ≥~67 s; LF lower bound 0.04 Hz → ≥250 s (~4 min). + // So HF/resp are valid from ~1 min, but LF (and therefore LF/HF) are NOT + // trustworthy below ~250 s — report them null rather than ship spectral noise. + if (span < 60) return none; + const HF_MIN_SPAN = 60; + const LF_MIN_SPAN = 250; const df = 0.005; // 5 mHz grid - const lf = lombScargleBand(t, x, LF_BAND[0], LF_BAND[1], df).power; + const hfBand = lombScargleBand(t, x, HF_BAND[0], HF_BAND[1], df); - const vlf = lombScargleBand(t, x, VLF_BAND[0], VLF_BAND[1], df).power; - const total = vlf + lf + hfBand.power; + const lfValid = span >= LF_MIN_SPAN; + const lf = lfValid ? lombScargleBand(t, x, LF_BAND[0], LF_BAND[1], df).power : null; + const vlf = lfValid ? lombScargleBand(t, x, VLF_BAND[0], VLF_BAND[1], df).power : null; + const total = lf != null && vlf != null ? vlf + lf + hfBand.power : null; // Respiratory rate = HF peak frequency × 60 (breaths/min). Confidence = how - // dominant that peak is vs the mean HF power (prominence). + // dominant that peak is vs the mean HF power (prominence). Valid from ~1 min. + const hfValid = span >= HF_MIN_SPAN; const meanHf = hfBand.power / ((HF_BAND[1] - HF_BAND[0]) / df); const prominence = meanHf > 0 ? hfBand.peakPower / meanHf : 0; - const respConf = Math.max(0, Math.min(1, (prominence - 1) / 4)); // ~1×→0, ~5×→1 + const respConf = hfValid ? Math.max(0, Math.min(1, (prominence - 1) / 4)) : 0; // ~1×→0, ~5×→1 const respRate = hfBand.peakFreq * 60; return { - lf: round(lf, 1), + lf: lf == null ? null : round(lf, 1), hf: round(hfBand.power, 1), - lf_hf: hfBand.power > 0 ? round(lf / hfBand.power, 3) : null, - total_power: round(total, 1), + lf_hf: lf != null && hfBand.power > 0 ? round(lf / hfBand.power, 3) : null, + total_power: total == null ? null : round(total, 1), resp_rate: respConf >= 0.3 ? round(respRate, 1) : null, resp_conf: round(respConf, 3), }; @@ -237,3 +245,52 @@ export function calcIrregular(rrRaw: number[]): Metric { inputs_used: ['rr_intervals'], }; } + +/** §Daytime HRV — waking-hours RMSSD timeline (ultradian autonomic rhythm). */ +export interface DaytimeHrvValue { + rmssd_median: number | null; // median of per-window RMSSD across the day + series: { ts: number; rmssd: number }[]; // per-window RMSSD (ultradian rhythm / stress timeline) + lowest_ts: number | null; // window with the lowest RMSSD (most-stressed point) + n_windows: number; +} + +/** + * calcDaytimeHrv(byMinute, bucketSec=300) — HRV across the WAKING day (not just sleep). + * `byMinute` is per-minute RR arrays (ms) over the daytime window. We bucket RR into + * fixed windows (default 5 min) and run the standard time-domain RMSSD per window with + * enough beats → an ultradian HRV / daytime-stress timeline. Tier HIGH (published + * RMSSD); abstains (null) when too few windows have usable RR — no fabrication. + */ +export function calcDaytimeHrv( + byMinute: { ts: number; rr: number[] }[], + bucketSec = 300, +): Metric { + const buckets = new Map(); + for (const m of byMinute) { + if (!m.rr || m.rr.length === 0) continue; + const key = Math.floor(m.ts / bucketSec); + const b = buckets.get(key) ?? { ts: key * bucketSec, rr: [] }; + for (const v of m.rr) b.rr.push(v); + buckets.set(key, b); + } + const series: { ts: number; rmssd: number }[] = []; + for (const b of [...buckets.values()].sort((a, c) => a.ts - c.ts)) { + const td = timeDomainHrv(b.rr); + if (td.rmssd != null) series.push({ ts: b.ts, rmssd: td.rmssd }); + } + if (series.length < 3) { + return { rmssd_median: null, series, lowest_ts: null, n_windows: series.length, confidence: 0, tier: 'HIGH', inputs_used: ['rr_intervals'] }; + } + const vals = series.map((s) => s.rmssd); + let lowest = series[0]; + for (const s of series) if (s.rmssd < lowest.rmssd) lowest = s; + return { + rmssd_median: round(median(vals) ?? 0, 1), + series, + lowest_ts: lowest.ts, + n_windows: series.length, + confidence: round(Math.min(1, series.length / 24), 3), // ~2 h of windows → full + tier: 'HIGH', + inputs_used: ['rr_intervals'], + }; +} diff --git a/src/illness.ts b/src/illness.ts index de03f18..e429fe3 100644 --- a/src/illness.ts +++ b/src/illness.ts @@ -1,46 +1,77 @@ -// §Illness — multivariate under-recovery / illness signal. Combines the three -// signals that co-move when the body is fighting something: resting HR ↑, -// nocturnal HRV (RMSSD) ↓, skin temperature ↑. We compute the MAHALANOBIS distance -// (Mahalanobis 1936) of today's 3-vector from the user's personal baseline mean + -// covariance — one scalar that accounts for how these normally co-vary, rather -// than three independent flags. Validated approach for wearable illness detection -// (Mishra et al., Nat Biomed Eng 2020; Smarr et al., Sci Rep 2020). +// §Illness — multivariate under-recovery / illness signal. Combines the signals +// that co-move when the body is fighting something: resting HR ↑, nocturnal HRV +// (RMSSD) ↓, skin temperature ↑, and respiratory rate ↑. We compute the +// MAHALANOBIS distance (Mahalanobis 1936) of today's vector from the user's +// personal baseline mean + covariance — one scalar that accounts for how these +// normally co-vary, rather than independent flags. Validated approach for +// wearable illness detection (Mishra et al., Nat Biomed Eng 2020; Smarr et al., +// Sci Rep 2020); elevated nocturnal respiratory rate is among the earliest +// infection signals (Mishra 2020; Natarajan, Lancet Digit Health 2020). // A SIGNAL, NOT A DIAGNOSIS — and only fires when deviations are in the illness // direction (not just "unusual"). -import type { Metric, IllnessValue, Driver, MetricRef } from './types'; +// +// Cycle-aware (honest, no fabrication): for users tracking their menstrual cycle, +// a rise in skin temperature and resting HR through the luteal phase is normal +// physiology (Wilcox 2000), not illness. When the only deviations are those two +// phase-expected signals, we do NOT flag illness — we still fire if HRV or +// respiratory rate also shift (signals not explained by the cycle). +import type { Metric, IllnessValue, Driver } from './types'; +import type { CyclePhase } from './cycle'; import { mean, stddev, round } from './util'; export interface IllnessToday { resting_hr: number | null; rmssd: number | null; // nocturnal RMSSD (ms) skin_temp: number | null; // RELATIVE temp index + resp_rate?: number | null; // nocturnal respiratory rate (breaths/min) } export interface IllnessHistory { resting_hr: number[]; rmssd: number[]; skin_temp: number[]; + resp_rate?: number[]; +} +export interface IllnessOpts { + /** Current menstrual-cycle phase, if the user tracks it — gates phase-expected + * RHR/temp rises so they don't masquerade as illness. */ + cyclePhase?: CyclePhase | null; } -/** Invert a symmetric 3×3 matrix; null if singular. */ -function inv3(m: number[][]): number[][] | null { - const [a, b, c] = m[0], [d, e, f] = m[1], [g, h, i] = m[2]; - const A = e * i - f * h, B = -(d * i - f * g), C = d * h - e * g; - const det = a * A + b * B + c * C; - if (Math.abs(det) < 1e-12) return null; - const id = 1 / det; - return [ - [A * id, (c * h - b * i) * id, (b * f - c * e) * id], - [B * id, (a * i - c * g) * id, (c * d - a * f) * id], - [C * id, (b * g - a * h) * id, (a * e - b * d) * id], - ]; +/** Invert a symmetric N×N matrix via Gauss-Jordan with partial pivoting; null if + * singular. Handles the 2-, 3- and 4-feature cases uniformly. */ +function invMatrix(m: number[][]): number[][] | null { + const n = m.length; + const a = m.map((row, i) => [...row, ...Array.from({ length: n }, (_, j) => (i === j ? 1 : 0))]); + for (let col = 0; col < n; col++) { + let piv = col; + for (let r = col + 1; r < n; r++) if (Math.abs(a[r][col]) > Math.abs(a[piv][col])) piv = r; + if (Math.abs(a[piv][col]) < 1e-12) return null; + [a[col], a[piv]] = [a[piv], a[col]]; + const d = a[col][col]; + for (let j = 0; j < 2 * n; j++) a[col][j] /= d; + for (let r = 0; r < n; r++) { + if (r === col) continue; + const f = a[r][col]; + for (let j = 0; j < 2 * n; j++) a[r][j] -= f * a[col][j]; + } + } + return a.map((row) => row.slice(n)); } +// Signals whose elevation is expected through the luteal/menstrual phase — not illness. +const CYCLE_EXPECTED = new Set(['rhr', 'temp']); + /** - * calcIllness(today, history). Needs ≥7 days of baseline per feature. Distance - * threshold 2.5 (≈ χ²₃ 95th pct √7.8 ≈ 2.79; we use 2.5 to be a touch sensitive). - * Tier ESTIMATE; "a signal, not a diagnosis." + * calcIllness(today, history, opts?). Needs ≥7 days of baseline per feature and + * ≥2 present features. Distance threshold 2.5 (≈ χ² 95th pct). Tier ESTIMATE; + * "a signal, not a diagnosis." When opts.cyclePhase is luteal/menstruation and the + * only deviating signals are RHR/temp, the signal is suppressed (phase-expected). */ -export function calcIllness(today: IllnessToday, history: IllnessHistory): Metric { +export function calcIllness( + today: IllnessToday, + history: IllnessHistory, + opts?: IllnessOpts, +): Metric { const NOTE = 'a signal, not a diagnosis'; // Build aligned feature set from whatever is present today + has ≥7 baseline. type Feat = { key: string; label: string; today: number; hist: number[]; dir: 1 | -1 }; @@ -51,6 +82,8 @@ export function calcIllness(today: IllnessToday, history: IllnessHistory): Metri cand.push({ key: 'rmssd', label: 'HRV (RMSSD)', today: today.rmssd, hist: history.rmssd, dir: -1 }); if (today.skin_temp != null && history.skin_temp.length >= 7) cand.push({ key: 'temp', label: 'Skin temperature', today: today.skin_temp, hist: history.skin_temp, dir: 1 }); + if (today.resp_rate != null && (history.resp_rate?.length ?? 0) >= 7) + cand.push({ key: 'resp', label: 'Respiratory rate', today: today.resp_rate, hist: history.resp_rate!, dir: 1 }); const none = (): Metric => ({ signal: false, distance: null, triggers: [], note: NOTE, @@ -64,15 +97,12 @@ export function calcIllness(today: IllnessToday, history: IllnessHistory): Metri return sd > 0 ? f.dir * (f.today - mu) / sd : 0; }); - // Covariance of the (illness-oriented) baseline features for Mahalanobis. - const n = cand[0].hist.length; const lens = cand.map((f) => f.hist.length); const minLen = Math.min(...lens); let distance: number; const drivers: Driver[] = []; const dim = cand.length; - // Standardized deviation vector. - const dvec = z; + const dvec = z; // standardized deviation vector if (dim >= 2 && minLen >= 7) { // Correlation matrix over the overlapping tail (standardized → correlation). const tail = cand.map((f) => f.hist.slice(-minLen)); @@ -87,11 +117,7 @@ export function calcIllness(today: IllnessToday, history: IllnessHistory): Metri } } // Mahalanobis² = dvecᵀ · C⁻¹ · dvec (using correlation since dvec is already z). - let inv = dim === 3 ? inv3(corr) : null; - if (dim === 2) { - const det = corr[0][0] * corr[1][1] - corr[0][1] * corr[1][0]; - inv = Math.abs(det) < 1e-12 ? null : [[corr[1][1] / det, -corr[0][1] / det], [-corr[1][0] / det, corr[0][0] / det]]; - } + const inv = invMatrix(corr); if (inv) { let d2 = 0; for (let a = 0; a < dim; a++) for (let b = 0; b < dim; b++) d2 += dvec[a] * inv[a][b] * dvec[b]; @@ -103,7 +129,10 @@ export function calcIllness(today: IllnessToday, history: IllnessHistory): Metri distance = Math.sqrt(dvec.reduce((s, v) => s + v * v, 0)); } - const metricFor = (key: string): string => key === 'rmssd' ? 'hrv' : key === 'rhr' ? 'rhr' : 'temp'; + const metricFor = (key: string): string => + key === 'rmssd' ? 'hrv' : key === 'rhr' ? 'rhr' : key === 'resp' ? 'resp' : 'temp'; + const inputName = (key: string): string => + key === 'rmssd' ? 'hrv_rmssd' : key === 'rhr' ? 'resting_hr' : key === 'resp' ? 'resp_rate' : 'skin_temp'; const triggers: string[] = []; cand.forEach((f, k) => { if (z[k] > 0.75) { @@ -116,13 +145,25 @@ export function calcIllness(today: IllnessToday, history: IllnessHistory): Metri }); // Signal: distance past threshold AND ≥2 features deviating toward illness. - const signal = distance > 2.5 && triggers.length >= 2; - const confidence = Math.min(0.6, (minLen / 30) * (cand.length / 3)); + let signal = distance > 2.5 && triggers.length >= 2; + let note = NOTE; + // Cycle-aware gate: through the luteal/menstrual phase, a temp & RHR rise is + // expected. If those are the ONLY deviations, don't flag illness — but still + // fire when HRV or respiration also shift (not explained by the cycle). + const inCyclePhase = opts?.cyclePhase === 'luteal' || opts?.cyclePhase === 'menstruation'; + if (signal && inCyclePhase) { + const corroborating = triggers.filter((t) => !CYCLE_EXPECTED.has(t)); + if (corroborating.length === 0) { + signal = false; + note = `${NOTE} (a rise in temperature & resting HR can be expected in this phase of your cycle)`; + } + } + const confidence = Math.min(0.6, (minLen / 30) * (cand.length / 4)); return { - signal, distance: round(distance, 2), triggers, note: NOTE, + signal, distance: round(distance, 2), triggers, note, confidence: round(confidence, 4), tier: 'ESTIMATE', - inputs_used: cand.map((f) => f.key === 'rmssd' ? 'hrv_rmssd' : f.key === 'rhr' ? 'resting_hr' : 'skin_temp'), + inputs_used: cand.map((f) => inputName(f.key)), drivers, }; } diff --git a/src/index.ts b/src/index.ts index c1ed946..01a29f6 100644 --- a/src/index.ts +++ b/src/index.ts @@ -22,6 +22,20 @@ export { calcCalories } from './calories'; // §5 Sleep export { calcSleep } from './sleep'; +// §5b Sleep v2 — multi-period (naps = shorter sleeps). Additive; calcSleep unchanged. +export { calcSleepPeriods } from './sleep'; +// Per-minute asleep/awake mask (calcSleep's boundary) — reconciles the hypnogram's awake. +export { sleepAwakeMask } from './sleep'; +// v1-method per-minute hypnogram (Cole-Kripke + HR-percentile bands) — single source. +export { stageHypnogram } from './sleep'; +export type { NightHypnogram } from './sleep'; +// §Sleep cycles — ultradian NREM↔REM cycles (Rosenblum 2024 fractal-cycle method, HRV-adapted). +export { detectSleepCycles } from './cycles'; +export type { SleepCycle, SleepCyclesValue } from './cycles'; + +// §Menstrual cycle — log-anchored calendar method + fertile window (Wilcox 2000). +export { calcCycle } from './cycle'; +export type { CycleValue, CyclePhase } from './cycle'; // §6 Sleep regularity (SRI) export { calcSleepRegularity } from './regularity'; @@ -35,10 +49,10 @@ export { calcHrRecovery, calcRecovery } from './recovery'; // §HRV — RMSSD/SDNN/pNN50, Lomb–Scargle LF/HF, Baevsky SI, RSA respiratory rate export { timeDomainHrv, freqDomainHrv, baevskyStressIndex, cleanRr, - calcHrvStability, calcIrregular, + calcHrvStability, calcIrregular, calcDaytimeHrv, VLF_BAND, LF_BAND, HF_BAND, } from './hrv'; -export type { TimeDomainHrv, FreqDomainHrv } from './hrv'; +export type { TimeDomainHrv, FreqDomainHrv, DaytimeHrvValue } from './hrv'; // §9 Training load / fitness trend export { calcLoad, calcFitnessTrend } from './trends'; @@ -46,6 +60,28 @@ export { calcLoad, calcFitnessTrend } from './trends'; // §Fitness — VO₂max (Uth–Sørensen), Banister fitness/fatigue/form, Foster monotony export { calcVo2Max, calcFitnessModel, calcMonotony } from './fitness'; +// §Steps — AN-2554 wrist pedometer (pure math; backend re-decodes the IMU + runs it) +export { calcSteps, pedometer, STEP_PARAMS } from './steps'; + +// §HAR — wrist activity recognition (Mannini 2013): per-window features + classifier +// + workout segmentation. LIVE high-rate stream only (flash is 1 Hz). +export { + extractHarFeatures, extractHarFeaturesFromSmv, classifyActivityWindow, segmentWorkout, + dwtDetailEnergies, DB10_LO, +} from './har'; +export type { + HarFeatures, ClassVote, WorkoutSegment, SegmentResult, +} from './har'; + +// §Circadian — CircaCP cosinor + bounded change-point (physiological-day anchor) +export { calcCircadian, stageSleep } from './circadian'; +export type { CircadianOpts, SleepStaging } from './circadian'; + +// §Sleep/wake ENSEMBLE — pluggable voters (Cole-Kripke + cardiac/CPD + van Hees); +// drives the demand-driven day-close trigger. detectWakeState + cheap peekRecentState. +export { detectWakeState, peekRecentState, coleKripke, cardiac, inactivity, DEFAULT_VOTERS } from './wake'; +export type { WakeContext, WakeState, WakeLabel, Voter } from './wake'; + // §Composite Readiness — weighted HRV + sleep blend (abstains without HRV) export { calcReadinessIndex } from './readiness_index'; export type { ReadinessInputs } from './readiness_index'; @@ -53,12 +89,12 @@ export type { ReadinessInputs } from './readiness_index'; // §10 Anomaly + illness (Mahalanobis). calcReadiness REMOVED (heuristic) — use // calcRecovery (HRV) instead. export { calcAnomaly } from './readiness'; -export type { AnomalyInputs } from './readiness'; +export type { AnomalyInputs, AnomalyOpts } from './readiness'; export { calcIllness } from './illness'; -export type { IllnessToday, IllnessHistory } from './illness'; +export type { IllnessToday, IllnessHistory, IllnessOpts } from './illness'; // §11 Baselines -export { calcBaselines } from './baselines'; +export { calcBaselines, seedBaselines, type SeededBaselines } from './baselines'; // Activity metric (steps / active-sedentary) REMOVED in v0 — see activity.ts. @@ -71,9 +107,17 @@ export type { // §12 Stress — HRV-based (Baevsky Stress Index + LF/HF, personal-relative). export { calcStress } from './stress'; +// §SpO₂ — RELATIVE blood-oxygen index + overnight desaturation screen. +export { calcSpo2Index, calcDesaturation } from './spo2'; +export type { Spo2Value, DesaturationValue } from './spo2'; + // §Sleep stress / nocturnal arousal (HR surge + motion during sleep). export { calcSleepStress } from './arousal'; +// §Restlessness — nocturnal movement fragmentation from per-minute actigraphy. +export { calcRestlessness } from './restlessness'; +export type { RestlessnessValue } from './restlessness'; + // §13 Nocturnal Heart (sleeping-HR dynamics + dip + elevated-overnight flag). export { calcNocturnalHeart } from './nocturnal'; export type { NocturnalValue } from './nocturnal'; diff --git a/src/readiness.ts b/src/readiness.ts index 03346f8..e41d40a 100644 --- a/src/readiness.ts +++ b/src/readiness.ts @@ -3,6 +3,7 @@ // keeps only the simple, rule-based RHR anomaly (Radin et al., Lancet Digit Health // 2020); the richer multivariate illness signal lives in illness.ts (Mahalanobis). import type { Baseline, Metric, AnomalyValue } from './types'; +import type { CyclePhase } from './cycle'; import { round } from './util'; export interface AnomalyInputs { @@ -12,6 +13,10 @@ export interface AnomalyInputs { sleep_efficiency?: number | null; // today baseline_sleep_efficiency?: number | null; // prior typical efficiency } +export interface AnomalyOpts { + /** menstrual-cycle phase if tracked — gates phase-expected RHR rises */ + cyclePhase?: CyclePhase | null; +} /** * calcAnomaly(inputs, baseline) @@ -25,9 +30,10 @@ export interface AnomalyInputs { */ export function calcAnomaly( inputs: AnomalyInputs, - baseline: Baseline + baseline: Baseline, + opts?: AnomalyOpts ): Metric { - const NOTE = 'signal, not a diagnosis'; + let NOTE = 'signal, not a diagnosis'; const triggers: string[] = []; const used: string[] = []; @@ -66,7 +72,14 @@ export function calcAnomaly( const ruleB = rhrUp && tempUp && effDown; if (ruleB) triggers.push('rhr_temp_efficiency'); - const signal = ruleA || ruleB; + // Cycle-aware gate: through the luteal/menstrual phase a RHR (and temp) rise is + // normal physiology. Suppress the pure-RHR rule A then; rule B keeps firing because + // its sleep-efficiency drop is NOT explained by the cycle. (Wilcox 2000.) + const inCyclePhase = opts?.cyclePhase === 'luteal' || opts?.cyclePhase === 'menstruation'; + const signal = (ruleA && !inCyclePhase) || ruleB; + if (ruleA && inCyclePhase && !ruleB) { + NOTE = 'signal, not a diagnosis (an elevated resting HR can be expected in this phase of your cycle)'; + } // Confidence scales with how many corroborating inputs were evaluable. const evaluable = [ diff --git a/src/regularity.ts b/src/regularity.ts index 41f5b40..7dedcfc 100644 --- a/src/regularity.ts +++ b/src/regularity.ts @@ -1,4 +1,12 @@ -// §6 Sleep regularity (SRI). Tier HIGH (≥3 nights). +// §6 Sleep timing regularity. Tier HIGH (≥3 nights). +// +// HONEST SCOPE: this is the circular variability of sleep-ONSET and WAKE clock-times +// (how much your bed/wake times wander night to night), scaled 0–100. It is a +// legitimate, useful regularity measure — but it is NOT the Phillips et al. 2017 +// Sleep Regularity Index (Sci Rep 7:3216), which is an epoch-by-epoch probability +// that two time points 24 h apart are in the same sleep/wake state. We don't plumb +// minute-level sleep/wake state across consecutive days, so we don't claim the SRI +// name. If that state is ever wired in, implement the real SRI and relabel. import type { NightSummary, Metric, SleepRegularityValue } from './types'; import { round } from './util'; @@ -41,10 +49,13 @@ function circularStdMin(minutesOfDay: number[]): number { * calcSleepRegularity(nights[]) * * For each night compute onset/wake minute-of-day. - * SRI = max(0, 100 − (avg(circ_std_onset, circ_std_wake)/120)*100). + * score = max(0, 100 − (avg(circ_std_onset, circ_std_wake)/120)*100). * Uses CIRCULAR std (see circularStdMin) so midnight-straddling bedtimes don't * spuriously read as irregular. confidence = 0 if <3 nights, else 0.7. * + * NOTE: the returned `sri` field is this onset/wake timing-regularity score, NOT + * the Phillips epoch-agreement Sleep Regularity Index — see the file header. + * * Confidence formula: pinned 0.7 once ≥3 nights with both onset & wake present; * 0 below that threshold (input incompleteness). */ diff --git a/src/restlessness.ts b/src/restlessness.ts new file mode 100644 index 0000000..eb9a0b9 --- /dev/null +++ b/src/restlessness.ts @@ -0,0 +1,76 @@ +// §Restlessness — nocturnal MOVEMENT fragmentation from per-minute actigraphy. +// +// Distinct from calcSleepStress (which detects sympathetic HR-surge arousals): this +// is purely about physical movement during the sleep window — how often you shifted, +// how fragmented the night was, your longest still stretch. Built ONLY from the +// per-minute `activity` aggregate we already store (works on 1 Hz flash R24 — no +// high-rate stream, no new storage). Tier ESTIMATE. +// +// "Movement" is defined relative to the night's own dynamic range (a robust floor + +// a fraction of the spread), so it adapts per-user and per-night rather than using a +// fabricated absolute threshold. +import type { Minute, Metric, Driver } from './types'; +import { percentile, round } from './util'; + +export interface RestlessnessValue { + score: number | null; // 0..100 (higher = more restless) + restless_min: number; // minutes with movement above the night floor + movement_bouts: number; // transitions still→moving (tosses/turns) + mobility_pct: number | null; // fraction of the window spent moving (0..1) + longest_still_min: number; // longest contiguous still stretch (sleep continuity) +} + +/** + * calcRestlessness(sleepMinutes) — sleepMinutes are the worn minutes within the main + * sleep period (onset..wake). A minute is "movement" when its activity exceeds a + * per-night threshold = p10 + 0.4·(p90 − p10). Bouts = still→move transitions. + */ +export function calcRestlessness(sleepMinutes: Minute[]): Metric { + const m = sleepMinutes + .filter((x) => x.wrist_on !== false && x.activity != null) + .sort((a, b) => a.ts - b.ts); + const empty = (): Metric => ({ + score: null, restless_min: 0, movement_bouts: 0, mobility_pct: null, longest_still_min: 0, + confidence: 0, tier: 'ESTIMATE', inputs_used: [], + }); + if (m.length < 20) return empty(); + + const acts = m.map((x) => x.activity); + const p10 = percentile(acts, 10) ?? 0, p90 = percentile(acts, 90) ?? 0; + const thresh = p10 + 0.4 * (p90 - p10); + + let restless = 0, bouts = 0, longestStill = 0, curStill = 0; + let moving = false; + for (const x of m) { + const isMove = x.activity > thresh && x.activity > 0; + if (isMove) { + restless++; + if (!moving) bouts++; // entering a movement bout (a toss/turn) + moving = true; + if (curStill > longestStill) longestStill = curStill; + curStill = 0; + } else { + moving = false; + curStill++; + } + } + if (curStill > longestStill) longestStill = curStill; + + const total = m.length; + const mobility = restless / total; + const hours = Math.max(0.5, total / 60); + const boutsPerHour = bouts / hours; + // Score: movement-bout density + mobility fraction, mapped to 0..100. + const score = Math.max(0, Math.min(100, Math.round(boutsPerHour * 6 + mobility * 100 * 0.5))); + + const drivers: Driver[] = [ + { label: 'Movement bouts', contribution: bouts, detail: `${bouts} shifts (${round(boutsPerHour, 1)}/h)`, ref: { metric: 'activity', scale: 'day' } }, + { label: 'Mobility', contribution: round(mobility * 100, 1), detail: `${restless}/${total} min moving`, ref: { metric: 'activity', scale: 'day' } }, + ]; + return { + score, restless_min: restless, movement_bouts: bouts, + mobility_pct: round(mobility, 4), longest_still_min: longestStill, + confidence: round(Math.min(1, total / 240), 4), tier: 'ESTIMATE', + inputs_used: ['activity'], drivers, + }; +} diff --git a/src/sessions.ts b/src/sessions.ts index 2060229..c2e47e4 100644 --- a/src/sessions.ts +++ b/src/sessions.ts @@ -1,10 +1,12 @@ // §7 Auto-workout detection. Tier HIGH (event) / ESTIMATE (type). -import type { Minute, Baseline, Profile, Metric, SessionValue } from './types'; +import type { Minute, Baseline, Profile, Metric, SessionValue, ActivityClass } from './types'; import { isHrUsable, resolveMaxHr, median, mean, round } from './util'; import { calcStrain } from './strain'; import { calcCalories } from './calories'; import { calcHrZones } from './zones'; import { calcHrRecovery } from './recovery'; +import { segmentWorkout } from './har'; +import type { ClassVote } from './har'; /** * detectSessions(minutes, baseline, profile?) @@ -64,10 +66,10 @@ export function detectSessions( i = lastAboveIdx + 1; } - // 2. Require ≥3 min sustained start AND mean activity > daily median. + // 2. Require ≥2 min sustained start AND mean activity > daily median. const qualified = segs.filter((s) => { const slice = worn.slice(s.startIdx, s.endIdx + 1); - if (slice.length < 3) return false; + if (slice.length < 2) return false; const meanAct = mean(slice.map((m) => m.activity)); return meanAct > dailyMedianAct; }); @@ -88,12 +90,12 @@ export function detectSessions( } } - // 4. Discard sessions <5 min total; build outputs. + // 4. Discard sessions <2 min total; build outputs. const out: Metric[] = []; for (const s of merged) { const slice = worn.slice(s.startIdx, s.endIdx + 1); const durationMin = (slice[slice.length - 1].ts - slice[0].ts) / 60 + 1; - if (durationMin < 5) continue; + if (durationMin < 2) continue; const hrs = slice.map((m) => m.hr_avg); const avgHr = mean(hrs); @@ -107,7 +109,28 @@ export function detectSessions( const zones = calcHrZones(slice, baseline, profile); const hrr = calcHrRecovery(slice, baseline, profile); - const type = classifyType(meanAct, dailyMedianAct, avgHr, rhr, maxHr); + // Type: prefer the motion-based HAR classes carried per-minute from ingest (live + // high-rate stream → Mannini classifier). Run segmentWorkout over the bout's + // per-minute act_class to get the primary type + graceful phase breakdown. If the + // bout has no classified minutes (flash-drained, 1 Hz, no motion texture), fall + // back to the crude HR/activity heuristic — honestly low-confidence. + const votes: ClassVote[] = slice + .filter((m) => m.act_class) + .map((m) => ({ ts: m.ts, cls: m.act_class as ActivityClass, conf: 1 })); + let type: string; + let typeConf: number; + let segments: SessionValue['segments']; + if (votes.length >= 2) { + const seg = segmentWorkout(votes, { minPhaseSec: 120 }); + type = seg.primary; + // Cap at 0.75: motion-classified is far better than the HR heuristic but still + // ESTIMATE (threshold classifier, no trained model yet). + typeConf = Math.min(0.75, Math.max(0.4, seg.type_confidence)); + segments = seg.segments.length > 1 ? seg.segments : undefined; + } else { + type = classifyType(meanAct, dailyMedianAct, avgHr, rhr, maxHr); + typeConf = 0.4; + } out.push({ start_ts: slice[0].ts, @@ -131,7 +154,9 @@ export function detectSessions( mean_activity: round(meanAct, 4), peak_activity: round(peakAct, 4), type, - type_confidence: 0.4, + type_confidence: round(typeConf, 2), + segments, + detected_type: type, confidence: 0.8, // event detection, HIGH tier: 'HIGH', inputs_used: ['hr_avg', 'hr_max', 'activity', 'baseline.resting_hr'], diff --git a/src/sleep.ts b/src/sleep.ts index a8d1272..5b45673 100644 --- a/src/sleep.ts +++ b/src/sleep.ts @@ -1,11 +1,30 @@ // §5 Sleep — Cole-Kripke actigraphy + HR-dip fusion. Tier HIGH (duration/eff), // ESTIMATE (stages). One-minute epochs. -import type { Minute, Baseline, Metric, SleepValue, SleepStages } from './types'; +import type { + Minute, Baseline, Metric, SleepValue, SleepStages, + SleepPeriod, SleepPeriodsValue, +} from './types'; import { isHrUsable, mean, round } from './util'; +import { cleanRr } from './hrv'; // Cole-Kripke weights over window [-4..+2]. const CK_W = [1.06, 0.54, 0.58, 0.76, 2.3, 0.74, 0.67]; +/** Per-minute RMSSD (ms) from a raw RR stream; null if too few clean beats. */ +function minuteRmssd(rr: number[] | undefined): number | null { + if (!rr || rr.length < 12) return null; + const c = cleanRr(rr); + if (c.length < 10) return null; + let s = 0; + for (let i = 1; i < c.length; i++) { const d = c[i] - c[i - 1]; s += d * d; } + return Math.sqrt(s / (c.length - 1)); +} +function medOfNullable(xs: (number | null)[]): number | null { + const a = xs.filter((x): x is number => x != null && Number.isFinite(x)).sort((p, q) => p - q); + return a.length ? a[a.length >> 1] : null; +} +const REM_RMSSD_FACTOR = 0.90; // high-HR minute with smoothed RMSSD < 0.90×asleep-median = REM, not wake + /** * calcSleep(minutes, baseline) * @@ -47,6 +66,42 @@ export function calcSleep(minutes: Minute[], baseline: Baseline): Metric m.wrist_on && m.hr_avg > 0) + .map((m) => m.hr_avg) + .sort((a, b) => a - b); + const pctl = (p: number) => + wornHr.length ? wornHr[Math.min(wornHr.length - 1, Math.floor(p * wornHr.length))] : rhr; + const sleepHr = Math.max(rhr, pctl(0.10)); // the quiet sleep-HR level for this window + const ASLEEP_HI = 1.05; // ≤ this × sleepHr → strong dip, nudge asleep + const AWAKE_HI = 1.20; // > this × sleepHr → clearly elevated, nudge awake (clears REM) + // Absolute backstop: HR at/above this × the RHR floor reads awake REGARDLESS of the + // window. Without it, a flat + motion-inert window anchors sleepHr to its own level and + // every minute falls under ASLEEP_HI*sleepHr → a sedentary-awake stretch (TV, desk) + // would be mis-read as a full night. Real sleep, including REM, rarely SUSTAINS >1.5× + // the 5th-percentile RHR floor, so this clips the false-positive without clipping sleep. + // (The proper tie-breaker is actigraphy; R24 carries none today — see decode note.) + const ABS_WAKE = 1.5; // 1. Cole-Kripke score + HR-dip fusion → boolean asleep per epoch. const asleep: boolean[] = new Array(n).fill(false); @@ -70,8 +125,9 @@ export function calcSleep(minutes: Minute[], baseline: Baseline): Metric 0) { - if (m.hr_avg < 0.95 * rhr) isAsleep = true; // strong dip → asleep - else if (m.hr_avg > 1.15 * rhr) isAsleep = false; // clearly elevated → awake + if (m.hr_avg > ABS_WAKE * rhr) isAsleep = false; // absolute backstop → awake regardless + else if (m.hr_avg <= ASLEEP_HI * sleepHr) isAsleep = true; // at/below the sleep level → asleep + else if (m.hr_avg > AWAKE_HI * sleepHr) isAsleep = false; // clearly elevated → awake } asleep[i] = isAsleep; } @@ -202,6 +258,283 @@ export function calcSleep(minutes: Minute[], baseline: Baseline): Metric { + const sorted = [...minutes].sort((a, b) => a.ts - b.ts); + const n = sorted.length; + const rhr = baseline.resting_hr; + + const empty = (): Metric => ({ + periods: [], + total_asleep_min: 0, + main_idx: null, + stages_beta: true, + confidence: 0, + tier: 'HIGH', + inputs_used: [], + }); + if (n === 0) return empty(); + if (rhr == null || rhr <= 0) return empty(); // need a resting-HR baseline (see calcSleep) + + // 1. Per-epoch asleep — identical scorer to calcSleep (duplicated here on + // purpose so v1 stays byte-for-byte untouched). + const asleep: boolean[] = new Array(n).fill(false); + for (let i = 0; i < n; i++) { + let s = 0; + for (let k = 0; k < CK_W.length; k++) { + const off = k - 4; + const idx = i + off; + if (idx >= 0 && idx < n) s += CK_W[k] * sorted[idx].activity; + } + s *= 0.001; + const m = sorted[i]; + if (!m.wrist_on) { asleep[i] = false; continue; } + let isAsleep = s < 1; + if (m.hr_avg > 0) { + if (m.hr_avg < 0.95 * rhr) isAsleep = true; + else if (m.hr_avg > 1.15 * rhr) isAsleep = false; + } + asleep[i] = isAsleep; + } + + // 2. Collect ALL consolidated periods (same ≤20-min interior-gap rule as the + // main-sleep detector). Each is trimmed to its first/last asleep epoch. + const MAX_GAP_MIN = 20; + const MAX_SLEEP_MIN = 14 * 60; // same plausibility ceiling, applied per period + const MIN_PERIOD_MIN = 15; // shorter than this isn't a sleep period + + const raw: { start: number; end: number; asleepN: number }[] = []; + let pf = -1, pl = -1, pa = 0, gap = 0; + const close = () => { + if (pf >= 0 && pa > 0) raw.push({ start: pf, end: pl, asleepN: pa }); + pf = -1; pl = -1; pa = 0; gap = 0; + }; + for (let i = 0; i < n; i++) { + if (asleep[i]) { if (pf < 0) pf = i; pl = i; pa++; gap = 0; } + else if (pf >= 0) { if (++gap > MAX_GAP_MIN) close(); } + } + close(); + + const periods: SleepPeriod[] = []; + for (const p of raw) { + let startIdx = p.start; + let endIdx = p.end; + if (endIdx - startIdx + 1 > MAX_SLEEP_MIN) endIdx = startIdx + MAX_SLEEP_MIN - 1; + + const span = sorted.slice(startIdx, endIdx + 1); + const in_bed_min = span.length; + let duration_min = 0; + for (let i = startIdx; i <= endIdx; i++) if (asleep[i]) duration_min++; + if (duration_min < MIN_PERIOD_MIN) continue; + + const efficiency = in_bed_min > 0 ? duration_min / in_bed_min : 0; + const sleepEpochs = span.filter((_, i) => asleep[startIdx + i]); + const stages = estimateStages(sleepEpochs, rhr); + + const hasHr = span.some((m) => m.wrist_on && m.hr_avg > 0); + const hasActivity = span.some((m) => m.activity > 0); + const hasTemp = baseline.skin_temp != null; + const inputCompleteness = [hasHr, hasActivity, hasTemp].filter(Boolean).length / 3; + const coverage = Math.min(1, in_bed_min / 90); + + periods.push({ + onset_ts: sorted[startIdx].ts, + wake_ts: sorted[endIdx].ts, + duration_min, + in_bed_min, + efficiency: round(efficiency, 4), + stages, + is_main: false, + confidence: round(inputCompleteness * coverage, 4), + }); + } + + if (periods.length === 0) return empty(); + + // 3. Flag the longest period as the main one (UI hint only). + let main_idx = 0; + for (let i = 1; i < periods.length; i++) { + if (periods[i].duration_min > periods[main_idx].duration_min) main_idx = i; + } + periods[main_idx].is_main = true; + + const total_asleep_min = periods.reduce((a, p) => a + p.duration_min, 0); + + const inputs_used: string[] = ['activity']; + if (sorted.some((m) => m.wrist_on && m.hr_avg > 0)) inputs_used.push('hr_avg'); + if (baseline.skin_temp != null) inputs_used.push('baseline.skin_temp'); + + return { + periods, + total_asleep_min, + main_idx, + stages_beta: true, + confidence: periods[main_idx].confidence, + tier: 'HIGH', + inputs_used, + }; +} + +/** + * sleepAwakeMask(minutes, baseline, rrByMin?) → Map per minute. + * Cole-Kripke actigraphy + HR-dip — the authoritative asleep/awake boundary — PLUS an + * RR tiebreaker for the override's blind spot: + * + * The HR-dip rule "HR > 1.15·rhr ⇒ awake" is right for genuine wake but it ALSO catches + * REM (REM HR legitimately runs ~15% above the sleeping floor) — and on a calm wrist with + * a near-dead activity signal there's no movement to tell them apart, so REM gets called + * "awake". Fix: for those high-HR minutes, look at beat-to-beat RR. REM = parasympathetic + * withdrawal ⇒ LOW RMSSD; so a high-HR minute whose smoothed RMSSD is below 0.90× the + * night's asleep-RMSSD median is REM → stays ASLEEP, not awake. (rrByMin optional; without + * it the legacy HR-only override stands.) + * + * The day-detail uses this to drive the hypnogram + breakdown from one source. Empty map + * when there's no resting-HR baseline. calcSleep/calcSleepPeriods scorers stay untouched. + */ +export function sleepAwakeMask( + minutes: Minute[], baseline: Baseline, rrByMin?: Map, +): Map { + const out = new Map(); + const rhr = baseline.resting_hr; + if (rhr == null || rhr <= 0) return out; + const sorted = [...minutes].sort((a, b) => a.ts - b.ts); + const n = sorted.length; + + // Per-minute smoothed RMSSD + the asleep-RMSSD median → the REM cut for the tiebreaker. + let rms: (number | null)[] = []; + let remCut: number | null = null; + if (rrByMin && rrByMin.size) { + const raw = sorted.map((m) => minuteRmssd(rrByMin.get(m.ts))); + rms = raw.map((_, i) => medOfNullable(raw.slice(Math.max(0, i - 2), Math.min(n, i + 3)))); + const asleepRms = sorted.map((m, i) => (m.hr_avg > 0 ? rms[i] : null)); + const med = medOfNullable(asleepRms); + if (med != null) remCut = REM_RMSSD_FACTOR * med; + } + + for (let i = 0; i < n; i++) { + let s = 0; + for (let k = 0; k < CK_W.length; k++) { + const idx = i + (k - 4); + if (idx >= 0 && idx < n) s += CK_W[k] * sorted[idx].activity; + } + s *= 0.001; + const m = sorted[i]; + if (!m.wrist_on) { out.set(m.ts, false); continue; } // off-wrist = not asleep (boundary) + let isAsleep = s < 1; + if (m.hr_avg > 0) { + if (m.hr_avg < 0.95 * rhr) isAsleep = true; + else if (m.hr_avg > 1.15 * rhr) { + // RR tiebreaker: high HR + low RMSSD = REM (asleep); else genuine wake. + const remLike = remCut != null && rms[i] != null && rms[i]! < remCut; + isAsleep = remLike; + } + } + out.set(m.ts, isAsleep); + } + return out; +} + +/** Merge per-minute stage runs shorter than the floor into the larger neighbour + * (awake keeps a higher floor) — consolidates the per-minute classifier's flicker + * into stable bouts so the hypnogram doesn't sawtooth. */ +function boutSmoothStage(labels: string[], minRun = 5, minAwakeRun = 7, passes = 6): string[] { + const s = [...labels]; + for (let p = 0; p < passes; p++) { + const runs: { a: number; b: number }[] = []; + for (let i = 0; i < s.length;) { let j = i; while (j < s.length && s[j] === s[i]) j++; runs.push({ a: i, b: j - 1 }); i = j; } + if (runs.length <= 1) break; + let changed = false; + for (let r = 0; r < runs.length; r++) { + const { a, b } = runs[r]; + const floor = s[a] === 'awake' ? minAwakeRun : minRun; + if (b - a + 1 >= floor) continue; + const prev = r > 0 ? runs[r - 1] : null; + const next = r < runs.length - 1 ? runs[r + 1] : null; + let tgt: string | null = null; + if (prev && next) tgt = (prev.b - prev.a) >= (next.b - next.a) ? s[prev.a] : s[next.a]; + else if (prev) tgt = s[prev.a]; + else if (next) tgt = s[next.a]; + if (tgt) { for (let x = a; x <= b; x++) s[x] = tgt; changed = true; } + } + if (!changed) break; + } + return s; +} + +export interface NightHypnogram { + hypnogram: { t: number; stage: 'awake' | 'light' | 'deep' | 'rem' }[]; + light_min: number; deep_min: number; rem_min: number; awake_min: number; asleep_min: number; +} + +/** + * stageHypnogram(minutes, onset, wake, baseline) — the v1 staging method, made + * per-minute and consistent. ONE source for the whole hypnogram + breakdown: + * • asleep/awake from calcSleep's Cole-Kripke + HR-dip mask (the proven detector), + * • deep/light/rem within the asleep minutes from the SAME HR-percentile bands as + * estimateStages (deep = bottom ~22% of sleeping HR, REM = top ~21%, else light), + * • bout-smoothed so it reads as stable bouts, not minute flicker. + * No RR, no circadian, no second stager — exactly what worked in v1, now driving both + * the graph and the totals (so they can never disagree). null if no resting-HR baseline. + */ +export function stageHypnogram( + minutes: Minute[], onset: number, wake: number, baseline: Baseline, rrByMin?: Map, +): NightHypnogram | null { + const rhr = baseline.resting_hr; + if (rhr == null || rhr <= 0) return null; + const mask = sleepAwakeMask(minutes, baseline, rrByMin); // ts → asleep (REM tiebreaker via RR) + const win = minutes.filter((m) => m.ts >= onset && m.ts <= wake).sort((a, b) => a.ts - b.ts); + if (win.length < 5) return null; + + // HR-percentile bands over the night's OWN sleeping HR (same as estimateStages). + const sleepHr = win.filter((m) => mask.get(m.ts) !== false && m.hr_avg > 0).map((m) => m.hr_avg); + const hrs = sleepHr.length ? sleepHr : win.filter((m) => m.hr_avg > 0).map((m) => m.hr_avg); + const sortedHr = [...hrs].sort((a, b) => a - b); + const meanHr = hrs.length ? hrs.reduce((a, b) => a + b, 0) / hrs.length : rhr; + const q = (p: number): number => sortedHr.length ? sortedHr[Math.min(sortedHr.length - 1, Math.floor(p * sortedHr.length))] : meanHr; + const deepEdge = q(0.22), remEdge = q(0.79); + const bigJump = Math.max(6, (hrs.length ? Math.max(1, q(0.9) - q(0.1)) : 1) * 0.6); + const acts = win.map((m) => m.activity); + const meanAct = acts.reduce((a, b) => a + b, 0) / (acts.length || 1); + + const raw: string[] = win.map((m, i) => { + if (mask.get(m.ts) === false || m.hr_avg <= 0) return 'awake'; + const hr = m.hr_avg; + const prev = i > 0 && win[i - 1].hr_avg > 0 ? win[i - 1].hr_avg : hr; + const next = i + 1 < win.length && win[i + 1].hr_avg > 0 ? win[i + 1].hr_avg : hr; + const hrJump = Math.max(Math.abs(hr - prev), Math.abs(hr - next)); + const lowAct = m.activity <= meanAct; + if (lowAct && hr <= deepEdge) return 'deep'; + if (lowAct && hr >= remEdge) return 'rem'; + if (lowAct && hrJump > bigJump) return 'rem'; + return 'light'; + }); + const sm = boutSmoothStage(raw); + let light = 0, deep = 0, rem = 0, awake = 0; + for (const s of sm) { if (s === 'awake') awake++; else if (s === 'deep') deep++; else if (s === 'rem') rem++; else light++; } + return { + hypnogram: win.map((m, i) => ({ t: m.ts, stage: sm[i] as 'awake' | 'light' | 'deep' | 'rem' })), + light_min: light, deep_min: deep, rem_min: rem, awake_min: awake, asleep_min: light + deep + rem, + }; +} + /** * BETA stage estimator. Splits asleep epochs into deep/REM/light using activity * + HR relative to that night's own distribution. Honest heuristic, not clinical. diff --git a/src/spo2.ts b/src/spo2.ts new file mode 100644 index 0000000..5876070 --- /dev/null +++ b/src/spo2.ts @@ -0,0 +1,120 @@ +// §SpO₂ — RELATIVE blood-oxygen index from the red/IR reflectance ratio. +// +// Pulse oximetry is grounded in the ratio R = red/IR: dividing the two LED +// channels cancels common-mode (perfusion, skin tone, contact, motion) that a +// single channel conflates. We do NOT have factory LED/photodiode calibration, +// and the 1 Hz historical record can't resolve the cardiac AC for a true +// ratio-of-ratios — so this is a RELATIVE index (tonight's median R vs the +// user's own baseline R), never an absolute %. Validated on 4 user datasets: +// the ratio is more stable than red-alone on clean signal (CV 0.79–4.81% vs +// 0.99–5.53%); on a noisy band it isn't, which is why the confidence below +// is gated on intra-night ratio stability + sample count — a noisy night +// yields low confidence (or null), never a misleading number. +import type { Metric, Driver, MetricRef } from './types'; +import { median, mean, stddev, round, clamp } from './util'; + +export interface Spo2Value { + /** Relative deviation vs personal baseline, in % of baseline ratio. + * Signed so that POSITIVE = better-oxygenated than your baseline (lower R). */ + index: number | null; + /** Tonight's median red/IR ratio (raw) — kept so the caller can roll the baseline. */ + night_ratio: number | null; +} + +const MIN_MINUTES = 30; // need ≥30 one-minute ratios for a defensible night +const PLAUSIBLE = (r: number) => r > 0.4 && r < 1.5; // reflectance-ratio sanity band +const CV_FLOOR = 0.08; // intra-night CV at/above which confidence → 0 + +/** + * calcSpo2Index(ratios, baselineRatio) + * + * ratios: per-minute red/IR ratios over the sleep window (one per usable minute). + * baselineRatio: the user's rolling baseline ratio, or null if not yet established. + * + * Returns a RELATIVE index + a computed confidence. With no baseline yet it reports + * night_ratio only (to seed the baseline) and a null index — honest, not guessed. + */ +export function calcSpo2Index(ratios: number[], baselineRatio: number | null): Metric { + const r = ratios.filter(PLAUSIBLE); + const none = (conf = 0): Metric => ({ + index: null, night_ratio: null, confidence: conf, tier: 'RELATIVE', inputs_used: [], + }); + if (r.length < MIN_MINUTES) return none(); + + const med = median(r); + if (med == null) return none(); + const nightR = round(med, 4); + const m = mean(r); + const cv = m > 0 ? stddev(r) / m : 1; + // confidence: more minutes (cap at 3 h) AND a stable within-night ratio. + const conf = round(clamp(Math.min(r.length / 180, 1) * Math.max(0, 1 - cv / CV_FLOOR), 0, 1), 3); + const inputs_used = ['spo2_red_raw', 'spo2_ir_raw']; + const ref: MetricRef = { metric: 'spo2', scale: 'day' }; + + // No personal baseline yet → seed value only, no fabricated index. + if (baselineRatio == null || !(baselineRatio > 0)) { + return { index: null, night_ratio: nightR, confidence: round(conf * 0.5, 3), tier: 'RELATIVE', inputs_used }; + } + + // + = lower ratio than baseline = more oxygenated than your norm. + const index = round(((baselineRatio - nightR) / baselineRatio) * 100, 2); + const drivers: Driver[] = [ + { label: 'Blood-oxygen vs baseline', contribution: index, detail: `R ${nightR} vs baseline ${round(baselineRatio, 4)}`, ref }, + ]; + return { index, night_ratio: nightR, confidence: conf, tier: 'RELATIVE', inputs_used, drivers }; +} + +export interface DesaturationValue { + /** count of relative-desaturation dips overnight (R rising ≥ DESAT_REL above baseline) */ + events: number; + /** desaturation index: events per hour of usable signal (apnea-screening ODI analogue) */ + odi: number | null; + /** deepest relative dip seen, in % of baseline ratio */ + deepest_pct: number | null; + note: string; +} + +const DESAT_REL = 0.04; // R rising ≥4% above baseline ≈ a relative desaturation dip +const DESAT_MINUTES = 1; // a dip must persist ≥1 min to count (collapse consecutive) + +/** + * calcDesaturation(ratios, baselineRatio) — nocturnal RELATIVE desaturation screen. + * + * A higher red/IR ratio than baseline = relatively LESS oxygenated. We count clustered + * dips where R rises ≥ DESAT_REL above the personal baseline and report an ODI-style + * events/hour. RELATIVE + SCREENING ONLY — no absolute %SpO₂, no diagnosis (we lack + * factory calibration and 1 Hz can't resolve true ratio-of-ratios). Confidence-gated: + * a noisy night or no baseline → low confidence / null, never a misleading number. + */ +export function calcDesaturation(ratios: number[], baselineRatio: number | null): Metric { + const NOTE = 'a screen, not a diagnosis'; + const r = ratios.filter(PLAUSIBLE); + const none = (conf = 0): Metric => ({ + events: 0, odi: null, deepest_pct: null, note: NOTE, confidence: conf, tier: 'RELATIVE', inputs_used: [], + }); + if (r.length < MIN_MINUTES || baselineRatio == null || !(baselineRatio > 0)) return none(); + + const thresh = baselineRatio * (1 + DESAT_REL); + let events = 0, run = 0, deepest = 0; + for (const v of r) { + if (v >= thresh) { + run++; + const dipPct = ((v - baselineRatio) / baselineRatio) * 100; + if (dipPct > deepest) deepest = dipPct; + if (run === DESAT_MINUTES) events++; // count once per sustained dip + } else { + run = 0; + } + } + const hours = Math.max(0.5, r.length / 60); + const m = mean(r); + const cv = m > 0 ? stddev(r) / m : 1; + const conf = round(clamp(Math.min(r.length / 180, 1) * Math.max(0, 1 - cv / CV_FLOOR), 0, 1), 3); + const drivers: Driver[] = [ + { label: 'Desaturation dips', contribution: events, detail: `${events} dips (${round(events / hours, 1)}/h)`, ref: { metric: 'spo2', scale: 'day' } }, + ]; + return { + events, odi: round(events / hours, 1), deepest_pct: round(deepest, 1), note: NOTE, + confidence: conf, tier: 'RELATIVE', inputs_used: ['spo2_red_raw', 'spo2_ir_raw'], drivers, + }; +} diff --git a/src/steps.ts b/src/steps.ts new file mode 100644 index 0000000..e62011b --- /dev/null +++ b/src/steps.ts @@ -0,0 +1,92 @@ +// §Steps — wrist pedometer (AN-2554, Analog Devices ADXL367 reference). PURE math +// over an already-decoded accelerometer-magnitude signal. The I/O around it — re- +// decoding the IMU frames from R2, ordering + per-minute grouping, persistence — +// lives in the backend runner (steps_imu.ts), exactly like the HRV/resp runners. +// +// Pipeline (per contiguous signal): +// sum-of-abs accel → low-pass moving average → centered-window max/min peak +// detection → dynamic threshold ± sensitivity/2 → CONFIRM consecutive "possible +// steps" to confirm (the regularity gate that rejects waving/typing/handling — +// validated to read 0 at rest). Params scaled to our ~100 Hz IMU; a calibration +// gain corrects the normal ~10% wrist undercount (locked vs a 100-step ground- +// truth walk: raw 90 → ×1.11 ≈ 100). + +// ── locked AN-2554 parameters (calibrated on a 100-step ground-truth walk) ── +const FS = 100 // assembled IMU sample rate (Hz) +const FILTER = 8 // low-pass moving-average taps +const WINDOW = 33 // centered peak window (~0.33 s @100 Hz) +const SENS = 0.10 // g — dead-zone around the dynamic threshold +const THR_ORDER = 4 // dynamic-threshold smoothing buffer +const CONFIRM = 8 // consecutive possible steps before counting (rejects non-gait) +const MAXMIN_TIMEOUT = 120 // samples to find a min after a max (~1.2 s) +const GAIN = 1.11 // calibration: raw 90 → ~100 on the ground-truth walk + +/** The locked parameters, exposed for documentation/tests. */ +export const STEP_PARAMS = { + FS, FILTER, WINDOW, SENS, THR_ORDER, CONFIRM, MAXMIN_TIMEOUT, GAIN, +} as const + +/** + * pedometer(sig) — AN-2554 time-domain step count over ONE contiguous + * accelerometer-magnitude signal (g, ~100 Hz). Raw count, no calibration gain. + */ +export function pedometer(sig: number[]): number { + const n = sig.length + if (n < WINDOW) return 0 + // low-pass: trailing moving average + const lp = new Array(n) + let acc = 0 + for (let i = 0; i < n; i++) { + acc += sig[i] + if (i >= FILTER) acc -= sig[i - FILTER] + lp[i] = acc / Math.min(i + 1, FILTER) + } + const half = WINDOW >> 1 + // centered-window extrema candidates + const cand: { i: number; max: boolean; v: number }[] = [] + for (let i = half; i < n - half; i++) { + let isMax = true, isMin = true + const v = lp[i] + for (let j = i - half; j <= i + half; j++) { + if (lp[j] > v) isMax = false + if (lp[j] < v) isMin = false + if (!isMax && !isMin) break + } + if (isMax) cand.push({ i, max: true, v }) + else if (isMin) cand.push({ i, max: false, v }) + } + // dynamic threshold + CONFIRM-step regularity + const dyn: number[] = [] + let dynVal = sig.reduce((s, v) => s + v, 0) / n + let steps = 0, poss = 0, regulation = false + let state: 'max' | 'min' = 'max' + let curMax = 0, curMaxIdx = -1 + for (const c of cand) { + if (state === 'max') { + if (c.max) { curMax = c.v; curMaxIdx = c.i; state = 'min' } + } else { + if (c.max) { if (c.v > curMax) { curMax = c.v; curMaxIdx = c.i } continue } + if (c.i - curMaxIdx > MAXMIN_TIMEOUT) { state = 'max'; poss = 0; regulation = false; continue } + const mx = curMax, mn = c.v + if (mx > dynVal + SENS / 2 && mn < dynVal - SENS / 2) { + if (mx - mn > SENS) { dyn.push((mx + mn) / 2); if (dyn.length > THR_ORDER) dyn.shift(); dynVal = dyn.reduce((s, v) => s + v, 0) / dyn.length } + poss++ + if (regulation) steps++ + else if (poss >= CONFIRM) { steps += poss; regulation = true } + } else { poss = 0; regulation = false } + state = 'max' + } + } + return steps +} + +/** + * calcSteps(minuteSignals) — run the pedometer over each per-minute contiguous + * signal, sum, and apply the calibration gain. `minuteSignals[m]` is the ordered + * (ts, sub-frame) magnitude samples for minute m. Returns the calibrated daily total. + */ +export function calcSteps(minuteSignals: number[][]): number { + let total = 0 + for (const sig of minuteSignals) total += pedometer(sig) + return Math.round(total * GAIN) +} diff --git a/src/trends.ts b/src/trends.ts index c00f3f1..a4cd33f 100644 --- a/src/trends.ts +++ b/src/trends.ts @@ -34,10 +34,18 @@ export function calcLoad(dailyStrain: DailyStrain[]): Metric { }; } - const last7 = sorted.slice(-7).map((d) => d.strain); - const last28 = sorted.slice(-28).map((d) => d.strain); - const acute = mean(last7); - const chronic = mean(last28); + // EWMA acute/chronic (Williams et al. 2017, BJSM 51:209) — exponentially- + // weighted, λ = 2/(N+1) with N=7 (acute) / N=28 (chronic). Decays older load + // smoothly instead of the rolling-average "cliff", and is more sensitive to + // when load actually occurred. Seed both at the first day's strain. + const LAMBDA_ACUTE = 2 / (7 + 1); // 0.25 + const LAMBDA_CHRONIC = 2 / (28 + 1); // ≈0.069 + let acute = sorted[0].strain; + let chronic = sorted[0].strain; + for (let i = 1; i < sorted.length; i++) { + acute = sorted[i].strain * LAMBDA_ACUTE + acute * (1 - LAMBDA_ACUTE); + chronic = sorted[i].strain * LAMBDA_CHRONIC + chronic * (1 - LAMBDA_CHRONIC); + } const acwr = chronic > 0 ? acute / chronic : null; let band: LoadValue['band'] = 'unknown'; diff --git a/src/types.ts b/src/types.ts index 5386b0c..a58e9ca 100644 --- a/src/types.ts +++ b/src/types.ts @@ -14,8 +14,18 @@ export interface Minute { activity: number; // motion magnitude (actigraphy count) steps: number; // detected steps this minute (accel peak-count; ESTIMATE) wrist_on: boolean; + // Dominant HAR activity class for this minute, classified at ingest from the live + // high-rate accel (see har.ts). Present ONLY for live-streamed minutes (flash is + // 1 Hz → no motion texture). Drives workout typing + segmentation in detectSessions. + act_class?: ActivityClass; } +/** Wrist activity-recognition class (Mannini 2013). Canonical home for the shared type. */ +export type ActivityClass = 'sedentary' | 'walk' | 'run' | 'cycle' | 'lift' | 'other'; + +/** One labelled phase of a workout (graceful activity switches). */ +export interface SessionSegment { start_ts: number; end_ts: number; type: string; confidence: number } + /** User profile. Fields are frequently absent — algorithms must degrade honestly. */ export interface Profile { age?: number; @@ -145,6 +155,27 @@ export interface SleepValue { stages_beta: boolean; } +// ── Sleep v2 (multi-period) — naps are just shorter sleeps ─────────────────── +/** One consolidated sleep period. Same breakdown as SleepValue, per-period. */ +export interface SleepPeriod { + onset_ts: number; + wake_ts: number; + duration_min: number; // asleep minutes + in_bed_min: number; + efficiency: number; // 0..1 + stages: SleepStages | null; // BETA/ESTIMATE + is_main: boolean; // longest period of the day (UI hint only; data is uniform) + confidence: number; // per-period detection confidence (0..1) +} + +/** All sleep periods detected in a window (one card each in the UI). */ +export interface SleepPeriodsValue { + periods: SleepPeriod[]; // chronological + total_asleep_min: number; // sum across periods + main_idx: number | null; // index of the main (longest) period, or null + stages_beta: boolean; +} + export interface SleepRegularityValue { sri: number; // 0..100 onset_std_min: number; @@ -165,8 +196,10 @@ export interface SessionValue { hrr60: number | null; mean_activity: number; peak_activity: number; - type: 'walk' | 'run/cardio' | 'strength/other'; - type_confidence: number; // ESTIMATE 0.4 + type: string; // ActivityClass when motion-classified, else legacy 'walk'|'run/cardio'|'strength/other' + type_confidence: number; // ESTIMATE + segments?: SessionSegment[]; // labelled phases (multi-activity workouts) + detected_type?: string; // the model's call at detection (for the calibration ledger) } export interface HrRecoveryValue { @@ -271,6 +304,20 @@ export interface BaselinesValue { days_used: number; } +// CircaCP circadian rhythm + main-sleep boundary (see circadian.ts). Cosinor +// phase is timezone-free (the physiological day anchor); onset/wake are the +// main-sleep boundary of the most-recent completed cycle. +export interface CircadianValue { + mesor: number | null; // bpm, rhythm-adjusted mean HR + amplitude: number | null; // bpm, cosinor amplitude (half peak-to-trough) + acrophase_ts: number | null; // unix s, HR peak (active-phase center) near window end + bathyphase_ts: number | null; // unix s, HR trough (rest-phase center) of the detected cycle + onset_ts: number | null; // main-sleep onset (HR drop), or null + wake_ts: number | null; // main-sleep wake (HR rise), or null + in_bed_min: number; // wake − onset, minutes + settled: boolean; // wake older than the settle window → night complete +} + // ── history shapes for trend/baseline fns ─────────────────────────────────── /** One night's sleep summary (for SRI). */ diff --git a/src/util.ts b/src/util.ts index d264ab8..e4260d4 100644 --- a/src/util.ts +++ b/src/util.ts @@ -92,7 +92,10 @@ export function resolveMaxHr( // 2. Age-predicted max (floor), taking a genuine above-age effort if present. if (profile?.age && profile.age > 0) { - const ageMax = 220 - profile.age; + // Tanaka, Monahan & Seals 2001 (JACC 37:153) — validated on 18,712 subjects, + // more accurate than 220−age (which over-estimates in the young, under- + // estimates in the old; they converge near age 40). + const ageMax = Math.round(208 - 0.7 * profile.age); if (observed > ageMax) return { maxHr: observed, source: 'measured' }; return { maxHr: ageMax, source: 'age' }; } diff --git a/src/wake.ts b/src/wake.ts new file mode 100644 index 0000000..af05d13 --- /dev/null +++ b/src/wake.ts @@ -0,0 +1,326 @@ +// wake.ts — sleep/wake state ENSEMBLE for the demand-driven day-close trigger. +// +// Purpose: from an assembled window of minute rollups, answer "is the user asleep +// or awake right now, and if they just woke, when?" — cheaply, deterministically, +// and HONESTLY (abstain to 'unknown' rather than fabricate). This is the gate that +// fires the once-per-day Tier-3/4 derivation at the user's REAL wake, not a clock. +// +// Method: a PLUGGABLE registry of per-minute voters, each emitting asleep|awake| +// unknown for every minute. We majority-vote per minute, consolidate into bouts, +// then report the current state + the most-recent sleep→wake boundary. Voters are +// swappable; the initial set is grounded in published, NON-TRAINED methods that fit +// our exact per-minute signals (hr_avg, activity, and optional RR): +// +// • coleKripke — Cole-Kripke 1992 actigraphy weighted-sum (benchmark-competitive +// among non-DL methods; ActiGraph/Newcastle PSG eval 2023). +// • cardiac — CPD's core finding (Cakmak 2020, Sleep): cardiac change-points +// carry the WAKE signal actigraphy misses. Uses HR level + change +// vs the night's own trough, and RR-SD (HRV) when RR is supplied. +// • inactivity — van Hees-style sustained-low-movement = sleep (non-trained). +// +// Thresholds are SELF-CALIBRATING off the window's own distribution (our `activity` +// is stddev-of-|accel|, not actigraph counts, so absolute published thresholds don't +// transfer — we use the published SHAPE with adaptive scaling). Tier = ESTIMATE; we +// do NOT claim the papers' trained-model accuracy. Validate against ground truth. + +import type { Minute, Baseline } from './types'; +import { isHrUsable } from './util'; + +export type WakeLabel = 'asleep' | 'awake' | 'unknown'; + +export interface WakeContext { + minutes: Minute[]; // time-ordered window (should span the night) + baseline: Baseline; // resting_hr anchor + rrByMin?: Map; // optional per-minute RR (ms) — enables the HRV arm of the cardiac voter + now?: number; // evaluation instant (default = last minute ts) +} + +/** A voter labels EVERY minute in the window. Pure; no I/O. */ +export type Voter = (ctx: WakeContext) => WakeLabel[]; + +export interface WakeState { + state: WakeLabel; // current state at `now` + wake_ts: number | null; // start of the most-recent sustained wake bout (the day boundary) + onset_ts: number | null; // start of the main sleep bout preceding that wake + awake_min: number; // sustained awake minutes ending at `now` + asleep_min: number; // minutes asleep in the main bout + votes: Record; // each voter's verdict at `now` (transparency) + confidence: number; // 0..1 (voter agreement × data coverage) +} + +// ── helpers ────────────────────────────────────────────────────────────────── +const MIN = 60; +const MIN_MAIN_SLEEP_MIN = 90; // a real night, not a catnap, must precede a "wake" +const SUSTAINED_WAKE_MIN = 10; // your rule: awake must hold ≥10 min before we fire + +function median(xs: number[]): number { + if (!xs.length) return 0; + const s = [...xs].sort((a, b) => a - b); + const m = s.length >> 1; + return s.length % 2 ? s[m] : (s[m - 1] + s[m]) / 2; +} + +// ── Voter 1: Cole-Kripke (actigraphy weighted-sum, published shape) ─────────── +// D = P·Σ wᵢ·Aᵢ over [−4..+2] minutes; sleep if D < 1. Our `activity` is rescaled +// to the window so the unit-dependent constants behave (self-calibrating). +const CK_W = [106, 54, 58, 76, 230, 74, 67]; // weights for A[-4..+2] +const CK_P = 0.001; +export const coleKripke: Voter = ({ minutes }) => { + const act = minutes.map((m) => (m.wrist_on ? m.activity : 0)); + // normalize activity to a counts-like scale: median wake-ish activity ≈ unit. + const nz = act.filter((a) => a > 0); + const scale = median(nz) || 1; // counts ≈ activity/scale + const n = minutes.length; + return minutes.map((_, i) => { + if (!isHrUsable(minutes[i]) && minutes[i].activity === 0) return 'unknown'; + let d = 0; + for (let k = -4; k <= 2; k++) { + const j = i + k; + if (j < 0 || j >= n) continue; + d += CK_W[k + 4] * (act[j] / scale); + } + d *= CK_P; + return d < 1 ? 'asleep' : 'awake'; + }); +}; + +// ── Voter 2: Cardiac (CPD core — the wake signal actigraphy misses) ─────────── +// Awake when SMOOTHED HR sits above the night's cardiac trough by a margin. This is +// the band's single most reliable wake signal (quiet wake has elevated HR but no +// motion — exactly where the actigraphy voters go blind). Smoothing (±2 min) means a +// single REM/arousal HR spike can't flip a minute to 'awake'; only a sustained rise +// does. The RR/HRV arm now lives in its OWN voter (hrvArousal) so it isn't +// double-counted. Trough is self-calibrating off the window, floored by RHR. +export const cardiac: Voter = ({ minutes, baseline }) => { + const usable = minutes.filter(isHrUsable).map((m) => m.hr_avg); + const sorted = [...usable].sort((a, b) => a - b); + const p10 = sorted.length ? sorted[Math.floor(sorted.length * 0.1)] : baseline.resting_hr; + const trough = Math.max(baseline.resting_hr || 0, p10 || 0) || (sorted[0] ?? 0); + const wakeMargin = 8; // bpm above the sleeping trough → cardiac wake (separates wake from REM) + // ±2-min median-smoothed HR (NaN where unusable) so spikes don't fragment the vote. + const hr = minutes.map((m) => (isHrUsable(m) ? m.hr_avg : NaN)); + const hs = hr.map((_, i) => { + const seg = hr.slice(Math.max(0, i - 2), Math.min(hr.length, i + 3)) + .filter((v) => Number.isFinite(v)).sort((a, b) => a - b); + return seg.length ? seg[seg.length >> 1] : NaN; + }); + return minutes.map((_, i) => { + if (!Number.isFinite(hs[i])) return 'unknown'; + return hs[i] > trough + wakeMargin ? 'awake' : 'asleep'; + }); +}; + +// ── Voter 3: van Hees-style sustained inactivity ────────────────────────────── +// Sleep = movement held below an adaptive low threshold across a rolling window. +export const inactivity: Voter = ({ minutes }) => { + const act = minutes.map((m) => (m.wrist_on ? m.activity : NaN)); + const worn = act.filter((a) => Number.isFinite(a)) as number[]; + const s = [...worn].sort((a, b) => a - b); + const pct = (p: number) => (s.length ? s[Math.min(s.length - 1, Math.floor(s.length * p))] : 0); + const p10 = pct(0.1), p90 = pct(0.9); + // "Still" = at/below the window's low (sleep) floor + a margin. Anchoring on the + // floor (not the median) is robust when the window is sleep-dominated, where the + // median itself sits inside the sleep mode. ABS_MOVE (~0.05 g RMS) is the still→move + // boundary in our `activity` units; the range term adapts when the window has daytime. + const ABS_MOVE = 0.05; + const thr = p10 + Math.max(ABS_MOVE, 0.3 * (p90 - p10)); + const W = 5; // minutes of sustained stillness + const n = minutes.length; + return minutes.map((_, i) => { + if (!Number.isFinite(act[i])) return 'unknown'; + let still = 0, seen = 0; + for (let j = Math.max(0, i - W); j <= Math.min(n - 1, i + W); j++) { + if (!Number.isFinite(act[j])) continue; + seen++; + if (act[j] <= thr) still++; + } + if (seen === 0) return 'unknown'; + return still / seen >= 0.7 ? 'asleep' : 'awake'; + }); +}; + +// ── Voter 4: HRV / RR autonomic arousal ─────────────────────────────────────── +// Beat-to-beat RR spread (SD) climbs with autonomic arousal toward and at wake. A +// minute is 'awake' when its short-window RR-SD is elevated AND HR sits above the +// sleeping trough. Returns 'unknown' when the minute has no RR — so a night without +// RR simply drops this voter rather than biasing it. This is the SECOND autonomic +// signal: paired with `cardiac` it lets the ≥2 consensus fire on quiet wake (HR up, +// no motion) without the two blind motion voters being able to veto it. +export const hrvArousal: Voter = ({ minutes, baseline, rrByMin }) => { + const usable = minutes.filter(isHrUsable).map((m) => m.hr_avg); + const sorted = [...usable].sort((a, b) => a - b); + const p10 = sorted.length ? sorted[Math.floor(sorted.length * 0.1)] : baseline.resting_hr; + const trough = Math.max(baseline.resting_hr || 0, p10 || 0) || (sorted[0] ?? 0); + const RR_SD_WAKE = 45; // ms — elevated beat-to-beat spread (validated on real RR) + // Per-minute RR-SD (NaN when <4 beats)… + const sdRaw = minutes.map((m) => { + const rr = rrByMin?.get(Math.floor(m.ts / MIN) * MIN); + if (!rr || rr.length < 4) return NaN; + const mean = rr.reduce((s, v) => s + v, 0) / rr.length; + return Math.sqrt(rr.reduce((s, v) => s + (v - mean) ** 2, 0) / rr.length); + }); + // …then ±2-min median-smoothed. Raw minute RR-SD bounces across the threshold, which + // fragments the morning wake so cardiac+hrv only intermittently reach the ≥2 bar; + // smoothing makes the autonomic signal continuous (mirrors the cardiac HR smoothing). + const sd = sdRaw.map((_, i) => { + const seg = sdRaw.slice(Math.max(0, i - 2), Math.min(sdRaw.length, i + 3)) + .filter((v) => Number.isFinite(v)).sort((a, b) => a - b); + return seg.length ? seg[seg.length >> 1] : NaN; + }); + return minutes.map((m, i) => { + if (!Number.isFinite(sd[i]) || !isHrUsable(m)) return 'unknown'; + return sd[i] > RR_SD_WAKE && m.hr_avg > trough ? 'awake' : 'asleep'; + }); +}; + +// Pluggable registry — swap/extend these as better algorithms are validated. +// FOUR voters, two families: motion (coleKripke, inactivity) + autonomic (cardiac, +// hrvArousal). The consensus rule (≥2) + bout-smoothing is HR-led by construction — +// motion is blind to quiet wake, so the autonomic pair must be able to carry a wake. +export const DEFAULT_VOTERS: { name: string; fn: Voter }[] = [ + { name: 'coleKripke', fn: coleKripke }, + { name: 'cardiac', fn: cardiac }, + { name: 'inactivity', fn: inactivity }, + { name: 'hrvArousal', fn: hrvArousal }, +]; + +/** + * Consensus label per minute: AWAKE if ≥`minAwake` voters say awake. NOT a majority — + * a 2–2 split counts as AWAKE on purpose. Two of the four voters are motion-based and + * physically blind to quiet wakefulness (lying still, HR up); a flat majority lets them + * veto a correct cardiac+HRV wake (and ties → 'unknown' → the close never fires). The + * ≥2 rule means the autonomic pair (cardiac, hrvArousal) can carry a wake on their own, + * while still requiring corroboration (one voter alone never flips it). Else asleep if + * any voter had a known read; else unknown. + */ +function consensusPerMinute(labels: WakeLabel[][], n: number, minAwake = 2): WakeLabel[] { + const out: WakeLabel[] = []; + for (let i = 0; i < n; i++) { + let awake = 0, known = 0; + for (const arr of labels) { + const l = arr[i]; + if (l === 'awake') { awake++; known++; } + else if (l === 'asleep') known++; + } + out.push(awake >= minAwake ? 'awake' : known ? 'asleep' : 'unknown'); + } + return out; +} + +/** + * Merge per-minute label runs shorter than `minRun` minutes into their larger + * neighbour (repeated to a fixed point), so intermittent voter agreement reads as one + * stable bout instead of a sawtooth. The HRV voter flickers in/out, which fragments a + * real continuous wake into sub-threshold pieces; smoothing bridges them. Same + * technique the sleep stager uses on its hypnogram. + */ +function boutSmooth(labels: WakeLabel[], minRun = 10, passes = 4): WakeLabel[] { + const s = [...labels]; + for (let p = 0; p < passes; p++) { + const runs: { a: number; b: number }[] = []; + for (let i = 0; i < s.length;) { let j = i; while (j < s.length && s[j] === s[i]) j++; runs.push({ a: i, b: j - 1 }); i = j; } + if (runs.length <= 1) break; + let changed = false; + for (let r = 0; r < runs.length; r++) { + const { a, b } = runs[r]; + if (b - a + 1 >= minRun) continue; + const prev = r > 0 ? runs[r - 1] : null; + const next = r < runs.length - 1 ? runs[r + 1] : null; + let tgt: WakeLabel | null = null; + if (prev && next) tgt = (prev.b - prev.a) >= (next.b - next.a) ? s[prev.a] : s[next.a]; + else if (prev) tgt = s[prev.a]; + else if (next) tgt = s[next.a]; + if (tgt) { for (let x = a; x <= b; x++) s[x] = tgt; changed = true; } + } + if (!changed) break; + } + return s; +} + +/** + * Run the ensemble over the window. Returns the current state + the most-recent + * sustained sleep→wake boundary (the physiological-day anchor). Requires ≥`minVote` + * agreement (default majority of the registry) AND a main sleep bout ≥90 min AND a + * sustained-awake stretch ≥10 min before reporting state='awake' with a wake_ts. + */ +export function detectWakeState( + ctx: WakeContext, + voters = DEFAULT_VOTERS, +): WakeState { + const minutes = ctx.minutes; + const n = minutes.length; + const now = ctx.now ?? (n ? minutes[n - 1].ts : 0); + const empty: WakeState = { state: 'unknown', wake_ts: null, onset_ts: null, awake_min: 0, asleep_min: 0, votes: {}, confidence: 0 }; + if (n < SUSTAINED_WAKE_MIN) return empty; + + const perVoter = voters.map((v) => v.fn(ctx)); + // HR-led consensus (≥2 awake, ties→awake) then bout-smoothed into stable bouts. + const labels = boutSmooth(consensusPerMinute(perVoter, n)); + + // votes at `now` (last minute) for transparency. + const votes: Record = {}; + voters.forEach((v, k) => { votes[v.name] = perVoter[k][n - 1] ?? 'unknown'; }); + + // Find the most-recent sleep bout and the wake that ends it. + // Walk backward: the current trailing run is "awake" if we just woke. + let i = n - 1; + // trailing awake run + let awakeRunStart = n; + while (i >= 0 && labels[i] === 'awake') { awakeRunStart = i; i--; } + // skip any unknown gap between wake and the sleep bout + while (i >= 0 && labels[i] === 'unknown') i--; + // the sleep bout + let sleepEnd = i; + while (i >= 0 && labels[i] !== 'awake') i--; // include asleep + interior unknowns + let sleepStart = i + 1; + + const sleepBoutMin = sleepEnd >= sleepStart && sleepStart >= 0 + ? Math.round((minutes[sleepEnd].ts - minutes[sleepStart].ts) / MIN) + 1 : 0; + const awakeMin = awakeRunStart < n + ? Math.round((minutes[n - 1].ts - minutes[awakeRunStart].ts) / MIN) + 1 : 0; + + // coverage = fraction of the window with a known label. + const known = labels.filter((l) => l !== 'unknown').length; + const coverage = n ? known / n : 0; + const agree = (() => { + // mean per-minute voter agreement over known minutes (for confidence). + let acc = 0, c = 0; + for (let k = 0; k < n; k++) { + let a = 0, w = 0; + for (const arr of perVoter) { if (arr[k] === 'asleep') a++; else if (arr[k] === 'awake') w++; } + const tot = a + w; if (!tot) continue; + acc += Math.max(a, w) / tot; c++; + } + return c ? acc / c : 0; + })(); + const confidence = Math.round(coverage * agree * 100) / 100; + + // Current state. + const current: WakeLabel = labels[n - 1]; + const justWoke = current === 'awake' + && awakeMin >= SUSTAINED_WAKE_MIN + && sleepBoutMin >= MIN_MAIN_SLEEP_MIN; + + return { + state: current, + wake_ts: justWoke ? minutes[awakeRunStart].ts : null, + onset_ts: justWoke && sleepStart >= 0 && sleepStart <= sleepEnd ? minutes[sleepStart].ts : null, + awake_min: awakeMin, + asleep_min: sleepBoutMin, + votes, + confidence, + }; +} + +/** + * Cheap recent-window check for the cron's per-tick ladder (ladder-step 2): is the + * user plausibly awake in the last few minutes? Liberal (high-recall) on purpose — + * a 'maybe' only triggers the full `detectWakeState`; it never suppresses a wake. + */ +export function peekRecentState(recent: Minute[], baseline: Baseline): WakeLabel { + const worn = recent.filter((m) => m.wrist_on); + if (worn.length < 3) return 'unknown'; + const hrUp = worn.filter(isHrUsable).some((m) => m.hr_avg > (baseline.resting_hr || 0) + 6); + const moving = worn.some((m) => m.activity > 0 && m.steps > 0); + return hrUp || moving ? 'awake' : 'asleep'; +} diff --git a/src/zones.ts b/src/zones.ts index 73fb27a..dabcce7 100644 --- a/src/zones.ts +++ b/src/zones.ts @@ -21,6 +21,10 @@ export function calcHrZones( const { maxHr, source } = resolveMaxHr(minutes, baseline, profile); const worn = minutes.filter(isHrUsable); + // %HRmax zones (standard, familiar). NOTE: Karvonen %HRR is more individualized + // in principle, but it requires a TRUSTWORTHY measured max — here maxHR is + // usually age-predicted, so %HRR adds no real accuracy and makes light days + // read as "no zones." Validated on real data; %HRmax kept deliberately. const z = [0, 0, 0, 0, 0]; for (const m of worn) { const pct = (m.hr_avg / maxHr) * 100;