diff --git a/.github/workflows/swift-packages.yml b/.github/workflows/swift-packages.yml index 43b47c0f06..c46bb25628 100644 --- a/.github/workflows/swift-packages.yml +++ b/.github/workflows/swift-packages.yml @@ -13,6 +13,7 @@ on: paths: - 'Packages/**' - 'Tools/SleepBench/**' + - 'Tools/SleepPSG/**' - 'Tools/Backfill/**' - '.github/workflows/swift-packages.yml' push: @@ -20,6 +21,7 @@ on: paths: - 'Packages/**' - 'Tools/SleepBench/**' + - 'Tools/SleepPSG/**' - 'Tools/Backfill/**' - '.github/workflows/swift-packages.yml' workflow_dispatch: @@ -61,6 +63,11 @@ jobs: # it — and #935 added 13 tests that no job ran. A tool whose numbers decide whether a staging change # lands has to keep compiling against the packages it measures, or it bitrots silently and the next # person to reach for it finds it broken. Backfill has no tests yet, so it is build-only. + # + # SleepPSG has a second, sharper reason to be here. Its variant scores run through a knob-for-knob PORT of + # `SleepStagerV2`, and its test suite is what proves that port still reproduces the shipped stager + # label-for-label. Without this job the port drifts the first time someone edits a coefficient, and the + # tool keeps printing confident numbers about a recipe that no longer exists. The check needs no dataset. tools: runs-on: macos-15 strategy: @@ -68,6 +75,7 @@ jobs: matrix: tool: - SleepBench + - SleepPSG - Backfill steps: - name: Checkout diff --git a/Tools/SleepPSG/Package.resolved b/Tools/SleepPSG/Package.resolved new file mode 100644 index 0000000000..fd559f8d7e --- /dev/null +++ b/Tools/SleepPSG/Package.resolved @@ -0,0 +1,14 @@ +{ + "pins" : [ + { + "identity" : "grdb.swift", + "kind" : "remoteSourceControl", + "location" : "https://github.com/groue/GRDB.swift.git", + "state" : { + "revision" : "2cf6c756e1e5ef6901ebae16576a7e4e4b834622", + "version" : "6.29.3" + } + } + ], + "version" : 2 +} diff --git a/Tools/SleepPSG/Package.swift b/Tools/SleepPSG/Package.swift new file mode 100644 index 0000000000..a2718485f9 --- /dev/null +++ b/Tools/SleepPSG/Package.swift @@ -0,0 +1,33 @@ +// swift-tools-version:5.9 +import PackageDescription + +// SleepPSG — score `SleepStagerV2` against POLYSOMNOGRAPHY truth. +// +// `Tools/SleepBench` replays the stagers against the references a wearer's own database carries: their +// manual restages, the strap's band `sleep_state`, Apple Health. Every one of those is either the recipe's +// own output fed back to it, or another vendor's black box. None of them is truth. +// +// This tool supplies the missing reference. It replays the SHIPPED `SleepStagerV2` over PhysioNet +// `sleep-accel` (Walch et al., SLEEP 2019) — 31 subjects of wrist accelerometer + heart rate with +// concurrent, human-scored PSG hypnograms — and scores it epoch-for-epoch against those hypnograms. +// +// The dataset is NEVER committed and is ALWAYS a command-line argument (`--dataset`). See README.md for +// the download step and the licence attribution both the dataset and its companion code require. +// +// Nothing here reads or writes a NOOP database, and no health data of any wearer's lives in this repo. +let package = Package( + name: "sleeppsg", + platforms: [.macOS(.v13)], + dependencies: [ + .package(path: "../../Packages/StrandAnalytics"), + .package(path: "../../Packages/WhoopProtocol"), + ], + targets: [ + .executableTarget(name: "sleeppsg", dependencies: ["StrandAnalytics", "WhoopProtocol"]), + // The port-equivalence check and every scoring primitive are pure functions — they need neither the + // dataset nor a database, so `swift test` here runs anywhere the packages build, with no arguments + // and no downloads. That is deliberate: the check that keeps this harness honest (the recipe port + // reproducing the shipped stager label-for-label) has to run in CI, where the dataset never exists. + .testTarget(name: "sleeppsgTests", dependencies: ["sleeppsg"]), + ] +) diff --git a/Tools/SleepPSG/README.md b/Tools/SleepPSG/README.md new file mode 100644 index 0000000000..45b5d98f7a --- /dev/null +++ b/Tools/SleepPSG/README.md @@ -0,0 +1,216 @@ +# SleepPSG — scoring `SleepStagerV2` against polysomnography + +`Tools/SleepBench` replays NOOP's sleep stagers against the references a wearer's own database happens to +carry: their manual restages, the strap's band `sleep_state`, Apple Health. Every one of those is either the +recipe's own output handed back to it, or another vendor's black box. None of them is truth, and the sleep +work has repeatedly had to say so in its own commit messages. + +This tool supplies the missing reference. It replays the **shipped** `SleepStagerV2` over PhysioNet +`sleep-accel` — wrist accelerometer and heart rate recorded alongside **human-scored PSG hypnograms** — and +scores it epoch-for-epoch against those hypnograms. + +Everything here is read-only. The dataset is never committed, never modified, and is always a command-line +argument. No NOOP database is opened and no wearer's health data is involved at any point. + +--- + +## Get the dataset + +~577 MB compressed. Put it anywhere **outside** this repository. + +```bash +mkdir -p ~/datasets && cd ~/datasets +curl -L -o sleep-accel.zip \ + https://physionet.org/static/published-projects/sleep-accel/motion-and-heart-rate-from-a-wrist-worn-wearable-and-labeled-sleep-from-polysomnography-1.0.0.zip +unzip -q sleep-accel.zip +``` + +The archive expands to a long-named directory containing `labels/`, `heart_rate/`, `motion/` and `steps/`. +Pass either that directory or its parent — `--dataset` finds the level with `labels/` in it. + +## Run it + +```bash +cd Tools/SleepPSG +swift run -c release sleeppsg --dataset ~/datasets/motion-and-heart-rate-*-1.0.0 +``` + +| flag | meaning | +|---|---| +| `--dataset PATH` | the extracted `sleep-accel` root. Required for every section except `port`. | +| `--section S` | `all` (default), or one of `port`, `baseline`, `strata`, `rem`, `variants`. | +| `--subjects N` | score only the first N subject ids — a fast smoke run. | +| `--csv PATH` | write the per-variant table as CSV. | +| `--seed N` | the port-validation corpus seed. | + +`--section port` needs no dataset at all: it is the equivalence check described below, and it is also what +`swift test` runs. + +--- + +## Attribution — required by both licences + +**Dataset.** Walch O, Huang Y, Forger D, Goldstein C. *Sleep stage prediction with raw acceleration and +photoplethysmography heart rate data derived from a consumer wearable device.* **SLEEP** 42(12), zsz180 +(2019). Distributed on PhysioNet as `sleep-accel` v1.0.0 under the **Open Data Commons Attribution License +v1.0**, which permits use and redistribution and *requires* attribution. + +**Companion code.** `ojwalch/sleep_classifiers`, **MIT**. Not vendored here; cited because it is the +reference implementation of the dataset's own baselines, and because anyone reproducing this work will want +it. The comparison models in section 5 of the report are written from scratch in this tool and are not +derived from that code. + +**PhysioNet.** Goldberger AL et al. *PhysioBank, PhysioToolkit, and PhysioNet.* Circulation 101(23):e215 +(2000). + +--- + +## What the report contains + +### 1. Port validation — why you can trust the variant rows + +Every **baseline** number comes from `StrandAnalytics.SleepStagerV2.stageSession` directly: the file that +ships in the app, compiled from `Packages/`, with no reimplementation in between. + +The **variant** rows cannot work that way. `SleepStagerV2` holds its constants as `static let`s and keeps +`Epoch`/`features()` internal to `StrandAnalytics`, so asking "what would this recipe score with one +transition row changed?" means either rebuilding the app for each question or restating the recipe with its +knobs exposed. This tool does the second, in `RecipePort.swift`. + +That is the classic way a benchmark quietly starts measuring something else — someone edits the shipped +file, nobody edits the port. So equivalence is not reviewed, it is **measured, on every build**: + +- 48 randomised nights (varying length, off-grid window starts, channel dropout, arousals, motion bursts, + unsorted input, R-R present on most so the RSA term is exercised — the PhysioNet data cannot exercise it) +- plus 7 degenerate cases: all-HR-missing, zero-variance HR, 1-epoch, 2-epoch, saturated motion gate, + no-coverage-at-all, and motion-absent +- both paths stage every night; **every epoch label must agree** + +`swift test` runs this, so CI runs it, with no dataset present. When `RecipeConfig.shipped` and +`SleepStagerV2.swift` diverge, the test fails and names the night and epoch. **A failure here is never +flakiness.** It means the two have drifted and the fix is to update `RecipeConfig.shipped`. + +### 2–3. Dataset summary, and the shipped recipe against PSG + +Four-class agreement (wake / light / deep / REM, with N1+N2 → light and N3 → deep), per-stage +precision/recall/F1, and **stage fractions as a first-class result**. + +Stage fractions are not an appendix here, for a specific reason: **kappa does not constrain them.** PR #348 +improved kappa on all three of its benchmarks and was reverted 48 hours later for re-scoring a healthy night +from 6 % to 23 % awake — *"kappa doesn't guard stage-fraction calibration."* A harness that reported only +kappa could not have seen that regression, and this one is built so it cannot repeat the mistake. + +**Two conventions, both printed, and a percentage here is meaningless without naming which.** Every stage +fraction in this tool divides by **all scored epochs of the night, wake included** — that part never varies. +What varies is how the 31 subjects are combined: + +| convention | what it is | deep, shipped recipe | +|---|---|---| +| **pooled over all scored epochs** | one denominator of 26 773 epochs for the whole cohort; a long night carries more weight than a short one | predicted **18.94 %**, truth **13.76 %**, bias **+5.18 pp** | +| **mean of per-subject percentages** | each subject's own % of their own night, then averaged over 31 subjects; every subject counts once | predicted **19.24 %**, truth **14.76 %**, bias **+4.48 pp** | + +Section 3 prints both tables, one above the other, each labelled. They are the same epochs and the same +denominator *within* a night; the ~1 pp gap between them is the weighting across subjects, nothing else. +The pooled convention is what a cohort-level "% of night" means and is the one section 6's variant table +and every bias column use. The per-subject mean is what a clinician comparing individuals would use, and it +is the one the self-check below is stated in, because it is the convention the previous harness reported. +**Quote either, never both interchangeably, and always say which.** + +`SleepStagerV2` fits nothing to data — every coefficient is fixed a priori — so there is no train/test split +to make for the recipe itself. The leave-one-subject-out machinery exists for the fitted comparison models +in section 5, which do need it. + +### 4. Stratified by night length + +Pooling hides couplings. `SleepStagerV2` has been observed to emit more REM as a fraction of sleep the +longer the night runs (corr +0.579 on one wearer's nights), traced to the REM-latency guard being a fixed +60-minute penalty — a third of a short night, a tenth of a long one. Subjects are split into terciles by +scored night length, and the correlation is reported for the prediction **and for PSG truth**: a coupling +truth also carries is physiology, one only the prediction carries is an artefact. + +### 5. Is REM detection physiology, or is it the clock? + +REM rises through the night in every population, so a model given nothing but "how far into the night is it" +already scores respectably. If NOOP's REM output were mostly that, the honest description of the feature +would change. + +Three models — clock features alone, physiology features alone, both — fit **leave-one-SUBJECT-out**. + +**Never an epoch-level split.** Consecutive 30 s epochs from one night are about as independent as +consecutive frames of a film; a random epoch split puts a subject's 02:14 epoch in train and their 02:14:30 +epoch in test, and every model then scores near-perfectly by memorising the subject. Each fold holds out one +whole subject. The decision threshold is chosen on the *training* folds and applied unchanged to the +held-out subject — REM is a ~20 % minority class, so a fixed 0.5 would under-predict it unequally across +three models with different calibration, and tuning on the test fold would be the leak the design exists to +avoid. + +### 6. Recipe variants + +Each row is the shipped recipe with **one** named change: PR #987's awake-transition row, and each of PR +#348's seven components measured alone, plus all seven together. Reported with kappa **and** the stage- +fraction biases, because a component is an improvement only if it wins the first without wrecking the +second. **Every percentage and every `bias pp` in this table is pooled over all 26 773 scored epochs** — +so the incumbent's deep bias reads +5.18 pp here, not the +4.48 pp the per-subject mean gives. + +The last row is not a candidate. `pre-#930 REM guard (provenance)` runs the REM-latency guard as it stood +before PR #930 — a hard `c < 0.12 ? 3.0 : 0.0` step in the session-fraction domain — and exists to explain +the self-check below. + +--- + +## The self-check, and what it found + +This tool replaces one that was deleted. Section 3 prints its measured values beside the figures that +harness reported on this same dataset, so anyone can see at a glance whether the rebuild is the same +instrument. **Nothing here is tuned to those numbers**, and where they disagree the disagreement is +reported rather than closed. + +**Every percentage in this section is the mean of per-subject percentages**, because that is the convention +the previous harness reported in; the pooled figures for the same quantities are ~1 pp away and are printed +directly above it in section 3. The tool labels both tables, and the self-check block repeats the +convention in its own header so a figure copied out of it carries its denominator with it. + +Every **truth-side** figure reproduces exactly: 31 subjects, **26 773** PSG-scored epochs, deep **14.76 %** +of each subject's own night averaged over the 31 (pooled: 13.76 %), truth median first-REM latency +**88.5 min**. Predicted deep lands at **19.24 %** on the same convention against a reported 19.25 % +(pooled: 18.94 %), and four-class kappa at **0.356** against a reported 0.349. + +The prediction side of REM does not match: more REM (**26.96 %** per-subject mean, against a reported +20.8 %), REM F1 0.569 vs 0.515, and a first REM period arriving much earlier (median 91.5 min vs 142.0). + +**The obvious explanation was tested and is wrong.** PR #930 replaced a fraction-domain `c < 0.12 ? 3.0 : 0` +step with a graded penalty measured in minutes from sleep onset, and the old figures were reported while +#930 was the candidate — so the incumbent they describe is the recipe *before* it. That predicts exactly +this pattern. The `pre-#930` variant runs that guard, and it moves nothing: kappa 0.356 → 0.356, REM F1 +0.569 → 0.569, REM 27.02 % → 26.97 % (pooled — section 6's table, like every variant row, is pooled). +On 8-hour lab nights the guard barely binds either way. The +hypothesis is falsified, and it is recorded here rather than deleted because the variant is what falsified +it. + +What the residual is remains **unresolved**, and the tool is not tuned to close it. Two constraints on any +future explanation, both from the table above: + +- It cannot be the session window. `clock` drives the deep prior and the REM ramp alike; predicted deep + reproduces to 0.01 pp (19.24 % measured against 19.25 % reported, per-subject mean), so `clock` is the + same quantity in both harnesses. +- It is REM-specific and it is not the latency guard, which the variant just ruled out. + +One suggestive coincidence, offered as an observation and not a conclusion: the previously reported "REM +F1 0.515" is exactly this harness's measured REM **precision** (0.515), against a measured F1 of 0.569. + +--- + +## What this benchmark is not entitled to claim + +- **The R-R stream is empty.** `sleep-accel` carries no beat-to-beat intervals, so `respRegularity` returns + nil on every epoch, `zrg` collapses to a constant 0, and the recipe's RSA respiration term contributes + exactly 0.0 to every emission. Five of the recipe's six inputs are live here; the sixth is silent. Results + are a lower bound on the same recipe with an R-R stream present. +- **Different hardware.** This is an Apple Watch's ~50 Hz raw accelerometer collapsed to a per-second mean, + not a WHOOP gravity decode. The absolute scale differs — which is precisely the thing `SleepStagerV2` was + built not to depend on, since every motion threshold is a multiple of the night's own median per-second + jerk. A benchmark on foreign hardware is a test of that claim, not a violation of it. +- **Staging only.** The session window handed to the stager is the labelled span. NOOP's own session + detection is not exercised, and nothing here says anything about it. +- **n = 31 subjects**, in a sleep-lab setting, wearing a device on the wrist alongside PSG leads. That is + more independent subjects than any other reference NOOP has, and it is still 31 people. diff --git a/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift b/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift new file mode 100644 index 0000000000..7e68f2ced2 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift @@ -0,0 +1,220 @@ +import Foundation + +// Ablation — the leave-one-SUBJECT-out model comparison that answers "is REM detection reading physiology, +// or is it reading the clock?" +// +// The question is not rhetorical. REM rises through the night in every population, so a model given nothing +// but "how far into the night is it" already scores respectably. If NOOP's REM output were mostly that, the +// per-epoch cardiac and movement terms would be decoration and the honest description of the feature would +// change. The test is to fit three models on the same epochs and compare: clock features alone, physiology +// features alone, and both. +// +// LEAVE-ONE-SUBJECT-OUT, never epoch-level. Consecutive 30 s epochs from one night are about as +// independent as consecutive frames of a film: a random epoch split puts a subject's 02:14 epoch in train +// and their 02:14:30 epoch in test, and every model then scores near-perfectly by memorising the subject. +// Each fold here holds out one whole subject, fits on the other 30, and is scored only on the held-out one. +// +// The decision threshold is chosen ON THE TRAINING FOLDS (the value maximising training F1) and then +// applied unchanged to the held-out subject. A fixed 0.5 would not be neutral: REM is ~20 % of a night, so +// 0.5 systematically under-predicts a minority class, and it would do so unequally across three models with +// different calibration. Tuning it on the test fold would be the leak this whole design exists to avoid. + +/// One epoch's evidence, already reduced to the columns the recipe itself sees. +struct AblationRow { + let subject: String + let isRem: Bool + /// Time-of-night fraction and its square — enough for a clock model to express a rise, a fall, or a + /// hump, so "the clock" is not strawmanned as a straight line. + let clock: [Double] + /// Per-night z-scored heart rate, heart-rate variability, movement, and the HR-flatness percentile — + /// the same normalised quantities `stageEpochs` builds its emissions from. + let physiology: [Double] +} + +enum AblationModel: String, CaseIterable { + case clock = "clock-only" + case physiology = "physiology-only" + case both = "both" + + func features(_ r: AblationRow) -> [Double] { + switch self { + case .clock: return r.clock + case .physiology: return r.physiology + case .both: return r.clock + r.physiology + } + } +} + +struct AblationResult { + let model: AblationModel + /// F1 over every held-out epoch pooled — one number for the cohort. + let pooledF1: Double + /// Mean of the per-subject F1s — one number per subject, so a long night cannot outvote a short one. + let meanSubjectF1: Double + let subjectF1: [String: Double] + let pooledPrecision: Double + let pooledRecall: Double +} + +enum Ablation { + + /// Fit and score all three models leave-one-subject-out. + /// + /// The design matrix is built ONCE per model as a flat `[Double]` with the intercept column baked in, + /// and each fold indexes into it. Filtering the row structs per fold instead — and rebuilding a feature + /// array per row per IRLS iteration — is tens of millions of allocations across 3 models × 31 folds × + /// 30 iterations, which in an unoptimised build is the difference between a minute and an afternoon. + /// The arithmetic is unchanged; only the bookkeeping is. + static func run(_ rows: [AblationRow]) -> [AblationResult] { + let subjects = Array(Set(rows.map { $0.subject })).sorted() + guard let first = rows.first else { return [] } + let y = rows.map { $0.isRem } + var out: [AblationResult] = [] + for model in AblationModel.allCases { + let p = model.features(first).count + 1 // + intercept column + var X = [Double](repeating: 1, count: rows.count * p) + for (i, r) in rows.enumerated() { + let fs = model.features(r) + for j in 0..= thr } + let ref = testIdx.map { y[$0] } + perSubject[held] = binaryF1(ref: ref, pred: pred) + pooledRef.append(contentsOf: ref) + pooledPred.append(contentsOf: pred) + } + var tp = 0, fp = 0, fn = 0 + for i in 0.. Double { 1.0 / (1.0 + exp(-max(-40, min(40, z)))) } + + static func dot(_ w: [Double], _ x: [Double]) -> Double { + var s = w[0] // intercept + for i in 0.. Double { + var s = 0.0 + let base = row * p + for j in 0..

[Double] { + var w = [Double](repeating: 0, count: p) + guard !rowIdx.isEmpty, p > 0 else { return w } + var H = [Double](repeating: 0, count: p * p) + var g = [Double](repeating: 0, count: p) + for _ in 0.. [Double] { + guard let first = X.first else { return [0] } + let p = first.count + 1 + var flat = [Double](repeating: 1, count: X.count * p) + for (i, row) in X.enumerated() { for j in 0.. [Double]? { + var a = a0, b = b0 + let n = b.count + for c in 0.. abs(a[piv][c]) { piv = r } + if abs(a[piv][c]) < 1e-12 { return nil } + if piv != c { a.swapAt(piv, c); b.swapAt(piv, c) } + for r in (c + 1).. Double { + var best = 0.5, bestF1 = -1.0 + var t = 0.02 + while t < 1.0 { + let pred = scores.map { $0 >= t } + let f1 = binaryF1(ref: truth, pred: pred) + if f1 > bestF1 { bestF1 = f1; best = t } + t += 0.02 + } + return best + } +} diff --git a/Tools/SleepPSG/Sources/sleeppsg/PortValidation.swift b/Tools/SleepPSG/Sources/sleeppsg/PortValidation.swift new file mode 100644 index 0000000000..8cd2877ba7 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/PortValidation.swift @@ -0,0 +1,303 @@ +import Foundation +import StrandAnalytics +import WhoopProtocol + +// PortValidation — the check that entitles this harness to say its variant numbers describe NOOP's shipped +// recipe rather than a lookalike. +// +// `RecipePort.swift` re-states `SleepStagerV2` with its constants exposed. That is the only way to ask +// "what would this recipe score with one constant moved", and it is also the classic way a benchmark +// quietly starts measuring something else: someone edits the shipped file, nobody edits the port, and every +// number the tool prints afterwards is about code that never ran on anyone's wrist. +// +// So the port is not reviewed for equivalence — it is MEASURED for it, on every build. Both paths stage the +// same nights and every epoch label must agree. The nights are generated here rather than loaded, for two +// reasons: the check has to run in CI where no dataset and no health data exist, and generated nights can +// be aimed at the paths a real night rarely reaches (a stream with no heart rate at all, a night whose +// motion gate fires on every epoch, a two-epoch nap). +// +// A failure here is never "the test is flaky". It means `RecipeConfig.shipped` and +// `Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift` have diverged; the report names +// the night and the first epoch that differs. + +/// SplitMix64 — a tiny deterministic PRNG so the generated corpus is identical on every machine and every +/// run. A benchmark whose validation set moves between runs cannot tell a real divergence from a reroll. +struct SplitMix64: RandomNumberGenerator { + private var state: UInt64 + init(seed: UInt64) { state = seed } + mutating func next() -> UInt64 { + state &+= 0x9E37_79B9_7F4A_7C15 + var z = state + z = (z ^ (z >> 30)) &* 0xBF58_476D_1CE4_E5B9 + z = (z ^ (z >> 27)) &* 0x94D0_49BB_1331_11EB + return z ^ (z >> 31) + } + mutating func double(_ lo: Double, _ hi: Double) -> Double { Double.random(in: lo.. Int { Int.random(in: lo...hi, using: &self) } + mutating func chance(_ p: Double) -> Bool { Double.random(in: 0..<1, using: &self) < p } +} + +/// One synthetic night: the exact arguments a `stageSession` call takes, plus a name for the report. +struct SyntheticNight { + let name: String + let start: Int + let end: Int + let grav: [GravitySample] + let hr: [HRSample] + let rr: [RRInterval] + let resp: [RespSample] +} + +enum PortValidation { + + /// The pinned corpus seed. Fixed rather than time- or machine-derived so that a divergence is always a + /// code change and never a reroll — and so a failure reported by CI reproduces exactly on a laptop. + static let defaultSeed: UInt64 = 0x5EED_5_1EE9_0001 + + /// The corpus: 48 randomised nights plus the degenerate cases. Deterministic given `seed`. + static func corpus(seed: UInt64 = PortValidation.defaultSeed, randomNights: Int = 48) -> [SyntheticNight] { + if seed == defaultSeed && randomNights == 48 { return pinnedCorpus } + return build(seed: seed, randomNights: randomNights) + } + + /// The default corpus, built once. Generating it is not free — a few million RNG draws per night in a + /// debug build — and several tests want the same nights, so building it per call would put minutes of + /// pure setup into `swift test`. A `static let` is initialised lazily and exactly once, thread-safely. + static let pinnedCorpus = build(seed: defaultSeed, randomNights: 48) + + private static func build(seed: UInt64, randomNights: Int) -> [SyntheticNight] { + var g = SplitMix64(seed: seed) + var out: [SyntheticNight] = [] + for i in 0.. SyntheticNight { + // ~20 min to ~3.7 h. Sized so the whole corpus lands in the same order of magnitude as a night's + // worth of epochs per subject without making `swift test` slow: the RSA term runs a band-limited + // DFT per epoch, and this check stages every night TWICE. + let nEpochs = g.int(40, 440) + // Deliberately misalign some windows from the 30 s wall-clock grid: `features` starts at the first + // multiple of 30 at or after `start`, and an off-grid start is where a tiling bug would show. + let start = 1_700_000_000 + g.int(0, 5000) * 30 + (g.chance(0.4) ? g.int(1, 29) : 0) + let end = start + nEpochs * 30 + (g.chance(0.5) ? g.int(0, 29) : 0) + + let haveHR = g.chance(0.9) + let haveRR = g.chance(0.75) + let haveGrav = g.chance(0.95) + let hrDropout = g.double(0.0, 0.55) + let gravDropout = g.double(0.0, 0.45) + + // A slow HR baseline with a couple of arousals, so the 5- and 11-min flatness windows see real + // structure rather than white noise. + let hrBase = g.double(44, 68) + let hrDrift = g.double(-6, 6) + let arousals = (0..= ibi { + rrCarry -= ibi + rr.append(RRInterval(ts: t, rrMs: Int((ibi * 1000).rounded()))) + } + } + } + } + // Hand some nights in unsorted: both paths sort defensively, and "defensively" is worth pinning. + if g.chance(0.25) { hr.shuffle(using: &g) } + if g.chance(0.25) { grav.shuffle(using: &g) } + return SyntheticNight(name: "random-\(index)", start: start, end: end, + grav: grav, hr: hr, rr: rr, resp: []) + } + + // MARK: - Degenerate nights + + /// The seven shapes a real night reaches rarely or never, each of which lands on a branch with its own + /// early return, fallback constant, or divide-by-zero guard. + static func degenerateNights() -> [SyntheticNight] { + let base = 1_700_000_000 + var out: [SyntheticNight] = [] + + // 1. All HR missing — every HR-derived feature is nil, the prefix grid is empty, and `stdOfSeconds` + // returns nil for every window. Only motion may speak. + do { + let n = 120 + var grav: [GravitySample] = [] + for s in 0..<(n * 30) { + let a = Double(s) * 0.0007 + grav.append(GravitySample(ts: base + s, x: cos(a), y: sin(a), z: 0.05)) + } + out.append(SyntheticNight(name: "degenerate-all-hr-missing", start: base, end: base + n * 30, + grav: grav, hr: [], rr: [], resp: [])) + } + + // 2. Zero-variance HR — `zfun`'s std is 0, which it must treat as 1.0 rather than dividing by it, + // and every flatness window ties, so the percentile rank is degenerate too. + do { + let n = 120 + var grav: [GravitySample] = [], hr: [HRSample] = [] + for s in 0..<(n * 30) { + hr.append(HRSample(ts: base + s, bpm: 55)) + let a = Double(s) * 0.0011 + grav.append(GravitySample(ts: base + s, x: cos(a), y: sin(a), z: 0.02)) + } + out.append(SyntheticNight(name: "degenerate-zero-variance-hr", start: base, end: base + n * 30, + grav: grav, hr: hr, rr: [], resp: [])) + } + + // 3. One epoch — the Viterbi lattice never takes a transition, and the segment tiler's first and + // last epoch are the same epoch. + do { out.append(shortNight(name: "degenerate-1-epoch", base: base, epochs: 1)) } + + // 4. Two epochs — exactly one transition, the smallest lattice that uses the transition matrix. + do { out.append(shortNight(name: "degenerate-2-epoch", base: base, epochs: 2)) } + + // 5. Saturated motion gate — every epoch's peak jerk clears the night-relative gate, so the wake + // boost is applied everywhere and the night-relative floor is at its most fragile. + do { + let n = 90 + var grav: [GravitySample] = [], hr: [HRSample] = [] + for s in 0..<(n * 30) { + hr.append(HRSample(ts: base + s, bpm: 58 + (s % 7))) + // Alternate poles: consecutive per-second jerks are ~2 g, orders above any median floor. + let f = s % 2 == 0 ? 1.0 : -1.0 + grav.append(GravitySample(ts: base + s, x: f, y: 0, z: 0)) + } + out.append(SyntheticNight(name: "degenerate-saturated-motion", start: base, end: base + n * 30, + grav: grav, hr: hr, rr: [], resp: [])) + } + + // 6. No coverage at all — `features` returns empty and the recipe must fall back to a single + // light-labelled segment spanning the window rather than returning nothing. + do { + out.append(SyntheticNight(name: "degenerate-no-coverage", start: base, end: base + 3600, + grav: [], hr: [], rr: [], resp: [])) + } + + // 7. Motion absent, heart rate only — no jerks anywhere, so the quiescent floor falls back to its + // 1e-6 epsilon and `motionQuiescent` is true on every epoch. + do { + let n = 100 + var hr: [HRSample] = [] + for s in 0..<(n * 30) { hr.append(HRSample(ts: base + s, bpm: Int(52 + 9 * sin(Double(s) * 0.002)))) } + out.append(SyntheticNight(name: "degenerate-no-motion", start: base, end: base + n * 30, + grav: [], hr: hr, rr: [], resp: [])) + } + return out + } + + private static func shortNight(name: String, base: Int, epochs: Int) -> SyntheticNight { + var grav: [GravitySample] = [], hr: [HRSample] = [] + for s in 0..<(epochs * 30) { + hr.append(HRSample(ts: base + s, bpm: 56 + (s % 5))) + let a = Double(s) * 0.02 + grav.append(GravitySample(ts: base + s, x: cos(a), y: sin(a), z: 0.1)) + } + return SyntheticNight(name: name, start: base, end: base + epochs * 30, + grav: grav, hr: hr, rr: [], resp: []) + } + + // MARK: - The check + + struct Divergence { + let night: String + let epochIndex: Int + let shipped: String + let port: String + } + + struct Result { + var nights = 0 + var randomNights = 0 + var degenerateNights = 0 + var epochs = 0 + var matchingEpochs = 0 + var divergences: [Divergence] = [] + var ok: Bool { divergences.isEmpty && epochs > 0 && matchingEpochs == epochs } + } + + /// Stage every night through BOTH paths and compare epoch labels on the 30 s grid. + /// + /// Epoch labels rather than the raw `[StageSegment]` array, because that is the unit every score in this + /// harness is computed on — an equivalence that held only after run-length collapsing would not + /// guarantee the numbers agree. (The segment arrays do also match; comparing labels is the stricter + /// statement for our purposes and the weaker one to state, so it is what the report claims.) + static func run(seed: UInt64 = PortValidation.defaultSeed, randomNights: Int = 48) -> Result { + let nights = corpus(seed: seed, randomNights: randomNights) + // Nights are staged concurrently: staging is pure, and the shipped stager's memo cache is itself + // lock-guarded (`AnalyticsMemoCache`), so the only shared state is safe. Serially this check runs + // for minutes in a debug build — the recipe's RSA term is a per-epoch DFT and every night is staged + // twice — and a validation nobody wants to wait for is a validation that gets skipped. + var perNight = [(epochs: Int, matching: Int, divergences: [Divergence])]( + repeating: (0, 0, []), count: nights.count) + perNight.withUnsafeMutableBufferPointer { buf in + let out = buf + DispatchQueue.concurrentPerform(iterations: nights.count) { idx in + let n = nights[idx] + let shipped = SleepStagerV2.stageSession(start: n.start, end: n.end, grav: n.grav, + hr: n.hr, rr: n.rr, resp: n.resp) + let port = V2Recipe.stageSession(start: n.start, end: n.end, grav: n.grav, + hr: n.hr, rr: n.rr, resp: n.resp, cfg: .shipped) + let a = epochLabels(shipped, start: n.start, end: n.end) + let b = epochLabels(port, start: n.start, end: n.end) + var epochs = 0, matching = 0 + var divs: [Divergence] = [] + for i in 0.. [StageSegment] { + let gravW = clipToWindow(grav, lo: start - padLo, hi: end + padHi, ts: { $0.ts }) + let hrW = clipToWindow(hr, lo: start - padLo, hi: end + padHi, ts: { $0.ts }) + let rrW = clipToWindow(rr, lo: start - padLo, hi: end + padHi, ts: { $0.ts }) + + let gravS = gravW.sorted { $0.ts < $1.ts } + let hrS = hrW.sorted { $0.ts < $1.ts } + let rrS = rrW.sortedByTsStable() + + let feats = features(start: start, end: end, grav: gravS, hr: hrS, rr: rrS, cfg: cfg) + if feats.isEmpty { return [StageSegment(start: start, end: end, stage: "light")] } + let labels = stageEpochs(feats, cfg: cfg) + + var segments: [StageSegment] = [] + for (i, f) in feats.enumerated() { + let stage = labels[i] == "awake" ? "wake" : labels[i] + let segStart = i == 0 ? start : f.start + let segEnd = i == feats.count - 1 ? end : feats[i + 1].start + if let last = segments.last, last.stage == stage { + segments[segments.count - 1].end = segEnd + } else { + segments.append(StageSegment(start: segStart, end: segEnd, stage: stage)) + } + } + return segments + } + + /// Stage `[start, end]` and return one label per FEATURISED epoch, alongside those epochs. The scoring + /// paths want the epoch grid the recipe actually built (an epoch with no coverage at all is skipped by + /// `features`, and a scorer that assumed a dense grid would silently shift every later label). + static func stageEpochsDetailed(start: Int, end: Int, grav: [GravitySample], hr: [HRSample], + rr: [RRInterval], resp: [RespSample], + cfg: RecipeConfig = .shipped) -> (epochs: [Epoch], labels: [String]) { + let gravW = clipToWindow(grav, lo: start - padLo, hi: end + padHi, ts: { $0.ts }) + let hrW = clipToWindow(hr, lo: start - padLo, hi: end + padHi, ts: { $0.ts }) + let rrW = clipToWindow(rr, lo: start - padLo, hi: end + padHi, ts: { $0.ts }) + let feats = features(start: start, end: end, + grav: gravW.sorted { $0.ts < $1.ts }, + hr: hrW.sorted { $0.ts < $1.ts }, + rr: rrW.sortedByTsStable(), cfg: cfg) + return (feats, stageEpochs(feats, cfg: cfg)) + } + + private static func clipToWindow(_ samples: [T], lo: Int, hi: Int, ts: (T) -> Int) -> [T] { + if samples.isEmpty { return samples } + if ts(samples[0]) >= lo && ts(samples[samples.count - 1]) < hi { return samples } + var l = 0, h = samples.count + while l < h { let m = (l + h) / 2; if ts(samples[m]) < lo { l = m + 1 } else { h = m } } + let start = l + l = start; h = samples.count + while l < h { let m = (l + h) / 2; if ts(samples[m]) < hi { l = m + 1 } else { h = m } } + if start == 0 && l == samples.count { return samples } + return Array(samples[start.. [Epoch] { + if end <= start { return [] } + let span = Double(max(1, end - start)) + + var hrSum = [Int: Double](), hrCnt = [Int: Int]() + for s in hr { hrSum[s.ts, default: 0] += Double(s.bpm); hrCnt[s.ts, default: 0] += 1 } + var secHR = [Int: Double](); secHR.reserveCapacity(hrSum.count) + for (k, v) in hrSum { secHR[k] = v / Double(hrCnt[k]!) } + + var gxSum = [Int: Double](), gySum = [Int: Double](), gzSum = [Int: Double](), gCnt = [Int: Int]() + for g in grav { + gxSum[g.ts, default: 0] += g.x; gySum[g.ts, default: 0] += g.y + gzSum[g.ts, default: 0] += g.z; gCnt[g.ts, default: 0] += 1 + } + var secG = [Int: (Double, Double, Double)](); secG.reserveCapacity(gCnt.count) + for (k, c) in gCnt { let d = Double(c); secG[k] = (gxSum[k]! / d, gySum[k]! / d, gzSum[k]! / d) } + + var rrBy = [Int: [Double]]() + for r in rr { rrBy[r.ts, default: []].append(Double(r.rrMs)) } + + let gridLo = secHR.keys.min() + let gridHi = secHR.keys.max() + var pCnt = [Int](), pSum = [Double](), pSq = [Double]() + if let g0 = gridLo, let g1 = gridHi { + let n = g1 - g0 + 1 + pCnt = [Int](repeating: 0, count: n + 1) + pSum = [Double](repeating: 0, count: n + 1) + pSq = [Double](repeating: 0, count: n + 1) + for i in 0.. Double? { + guard let g0 = gridLo, let g1 = gridHi else { return nil } + let a = max(lo, g0) - g0 + let b = min(hi, g1 + 1) - g0 + if b <= a { return nil } + let cnt = pCnt[b] - pCnt[a] + if cnt < 2 { return nil } + let n = Double(cnt) + let sv = pSum[b] - pSum[a] + let sq = pSq[b] - pSq[a] + let m = sv / n + let v = (sq - 2 * m * sv + n * m * m) / n + return (v < 0 ? 0 : v).squareRoot() + } + + struct Raw { + let start: Int; let hr: Double?; let hrVar: Double?; let hrFlat11: Double? + let jerks: [Double]; let gapSec: Int; let jerkMax: Double; let respReg: Double?; let clock: Double + let minutes: Double + } + var raws: [Raw] = [] + var allJerks: [Double] = [] + let firstE = ((start + 29) / 30) * 30 + var e = firstE + while e < end { + var hrs: [Double] = [] + var gseq: [(Double, Double, Double)] = [] + for s in e..<(e + 30) { + if let h = secHR[s] { hrs.append(h) } + if let g = secG[s] { gseq.append(g) } + } + if hrs.isEmpty && gseq.isEmpty { e += 30; continue } + + var jerks: [Double] = [] + for i in 1.. moveThr ? 1 : 0) } + feats.append(Epoch( + start: r.start, hr: r.hr, hrVar: r.hrVar, hrFlat11: r.hrFlat11, + moveFrac: Double(moves) / Double(r.gapSec), jerkMax: r.jerkMax, respReg: r.respReg, + clock: r.clock, jerkScale: jerkScale, minutesSinceOnset: r.minutes)) + } + return feats + } + + static func respRegularity(_ beats: [(Double, Double)]) -> Double? { + if beats.count < 12 { return nil } + let t0 = beats.first!.0, tN = beats.last!.0 + if tN <= t0 { return nil } + let n = Int(ceil((tN - t0) / 0.25 - 1e-9)) + if n < 16 { return nil } + + var y = [Double](repeating: 0, count: n) + var seg = 0 + for i in 0.. maxP { maxP = p } + } + if sumP == 0 { return nil } + return maxP / sumP + } + + // MARK: - Staging + + /// Awake-emission dead-zone. Identity when `awakeDeadzone <= 0`, which is the shipped state. + static func dz(_ z: Double, _ deadzone: Double) -> Double { + if deadzone <= 0.0 { return z } + if z > deadzone { return z - deadzone } + if z < -deadzone { return z + deadzone } + return 0.0 + } + + static func motionQuiescent(_ f: Epoch, _ cfg: RecipeConfig) -> Bool { + f.moveFrac <= 0.0 && f.jerkMax <= f.jerkScale * cfg.jerkFloorGateMult + } + + static func remLatencyGuard(_ minutesSinceOnset: Double, _ cfg: RecipeConfig) -> Double { + cfg.remLatencyPenalty * min(1.0, max(0.0, 1.0 - minutesSinceOnset / cfg.remLatencyMinutes)) + } + + static func cyclePrior(_ c: Double, _ minutesSinceOnset: Double, _ cfg: RecipeConfig) -> [String: Double] { + ["deep": 1.2 * max(0.0, 1.0 - c / 0.55), + "rem": 1.0 * c - remLatencyGuard(minutesSinceOnset, cfg), + "light": 0.0, "awake": 0.0] + } + + static func sustainedSleepOnset(_ labels: [String], _ cfg: RecipeConfig) -> Int? { + var run = 0 + for i in labels.indices { + if labels[i] == "awake" { run = 0; continue } + run += 1 + if run >= cfg.onsetSustainedEpochs { return i - cfg.onsetSustainedEpochs + 1 } + } + return nil + } + + static func viterbi(_ emSeq: [[String: Double]], _ cfg: RecipeConfig) -> [String] { + if emSeq.isEmpty { return [] } + let logT = cfg.transition.mapValues { row in row.mapValues { log(max($0, 1e-9)) } } + var V = emSeq[0] + var back: [[String: String]] = [] + for t in 1.. bestVal { bestVal = val; bestPrev = p } + } + newV[s] = bestVal + emSeq[t][s]! + bp[s] = bestPrev + } + V = newV; back.append(bp) + } + var last = stageNames[0], lastV = V[last]! + for s in stageNames.dropFirst() where V[s]! > lastV { lastV = V[s]!; last = s } + var path = [last] + for bp in back.reversed() { last = bp[last]!; path.append(last) } + return path.reversed() + } + + static func stageEpochs(_ feats: [Epoch], cfg: RecipeConfig = .shipped) -> [String] { + if feats.isEmpty { return [] } + let prior = cfg.baseLogPrior + + func zfun(_ vals: [Double?]) -> (Double?) -> Double { + let present = vals.compactMap { $0 } + if present.isEmpty { return { _ in 0.0 } } + let m = present.reduce(0, +) / Double(present.count) + let sd0 = (present.reduce(0.0) { $0 + ($1 - m) * ($1 - m) } / Double(present.count)).squareRoot() + let sd = sd0 == 0 ? 1.0 : sd0 + return { v in v == nil ? 0.0 : (v! - m) / sd } + } + let zhr = zfun(feats.map { $0.hr }) + let zhv = zfun(feats.map { $0.hrVar }) + let zmv = zfun(feats.map { Optional($0.moveFrac) }) + let zrg = zfun(feats.map { $0.respReg }) + + let fsorted = feats.compactMap { $0.hrFlat11 }.sorted() + func fpct(_ v: Double?) -> Double { + guard let v = v, !fsorted.isEmpty else { return 0.5 } + var lo = 0, hi = fsorted.count + while lo < hi { let mid = (lo + hi) / 2; if fsorted[mid] <= v { lo = mid + 1 } else { hi = mid } } + return Double(lo) / Double(fsorted.count) + } + + var seq: [[String: Double]] = [] + seq.reserveCapacity(feats.count) + for f in feats { + let zhrv = zhr(f.hr), zhvv = zhv(f.hrVar), zmvv = zmv(f.moveFrac) + let gate = cfg.deepGateSlope * max(0.0, fpct(f.hrFlat11) - cfg.deepGateThresh) + let awakeCardiac0 = cfg.awakeZhv * dz(zhvv, cfg.awakeDeadzone) + + cfg.awakeZhr * dz(zhrv, cfg.awakeDeadzone) + let awakeCardiac = motionQuiescent(f, cfg) ? min(0.0, awakeCardiac0) : awakeCardiac0 + var em: [String: Double] = [ + "deep": cfg.deepZhv * zhvv + cfg.deepZhr * zhrv + cfg.deepZmv * zmvv - gate + prior["deep"]!, + "rem": cfg.remZhv * zhvv + cfg.remZmv * zmvv + cfg.remZhr * zhrv + prior["rem"]!, + "light": prior["light"]!, + "awake": cfg.awakeZmv * zmvv + awakeCardiac + prior["awake"]!, + ] + // In `.gradedFromOnset` the guard is DISABLED here (`.infinity` ⇒ guard = 0) and added by pass 2 + // once onset is known. In `.preNine30FractionStep` there is no onset to wait for — the step is + // a function of the session fraction alone — so it goes in now and there is no second pass. + let pr: [String: Double] + switch cfg.remLatencyMode { + case .gradedFromOnset: + pr = cyclePrior(f.clock, .infinity, cfg) + case .preNine30FractionStep: + pr = ["deep": 1.2 * max(0.0, 1.0 - f.clock / 0.55), + "rem": 1.0 * f.clock - (f.clock < 0.12 ? cfg.remLatencyPenalty : 0.0), + "light": 0.0, "awake": 0.0] + } + for s in stageNames { em[s]! += pr[s]! } + if f.jerkMax > f.jerkScale * cfg.jerkFloorGateMult { em["awake"]! += cfg.motionGateBoost } + if let rg = f.respReg { let z = zrg(rg); em["deep"]! += cfg.respWeight * z; em["rem"]! -= cfg.respWeight * z } + seq.append(em) + } + + if cfg.remLatencyMode == .preNine30FractionStep { return viterbi(seq, cfg) } + + let provisional = viterbi(seq, cfg) + let originMin = sustainedSleepOnset(provisional, cfg).map { feats[$0].minutesSinceOnset } ?? 0.0 + for i in feats.indices { + seq[i]["rem"]! -= remLatencyGuard(feats[i].minutesSinceOnset - originMin, cfg) + } + return viterbi(seq, cfg) + } +} diff --git a/Tools/SleepPSG/Sources/sleeppsg/Scoring.swift b/Tools/SleepPSG/Sources/sleeppsg/Scoring.swift new file mode 100644 index 0000000000..31ab5b22c8 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/Scoring.swift @@ -0,0 +1,229 @@ +import Foundation +import StrandAnalytics + +// Scoring — the agreement and calibration primitives, over plain label arrays. +// +// The definitions here match `Tools/SleepBench/Sources/sleepbench/Metrics.swift` deliberately: the two +// harnesses score the same recipe against different references, and a kappa that meant something slightly +// different in each would make their numbers incomparable exactly when someone needs to compare them. +// `ScoringTests` pins the shared ones against hand-computed values so a divergence is caught rather than +// discovered. +// +// The one primitive SleepBench does not carry is per-stage F1. It belongs here because the question this +// dataset exists to answer is about REM specifically, and REM is a minority class: on a night that is 20 % +// REM, a recipe that emits no REM at all still scores 80 % accuracy and a respectable four-class kappa. +// Precision and recall for that one class, and their harmonic mean, are what move when REM detection does. + +let stageOrder = ["wake", "light", "deep", "rem"] +let epochSeconds = 30.0 + +func epochCount(start: Int, end: Int) -> Int { + max(1, Int(ceil(Double(end - start) / epochSeconds))) +} + +/// Expand a `[StageSegment]` tiling into a per-epoch label array on the 30 s grid, taking each epoch's +/// label from the segment covering its START instant — the same convention `SleepBench` uses. +func epochLabels(_ stages: [StageSegment], start: Int, end: Int) -> [String] { + let n = epochCount(start: start, end: end) + var out = [String](repeating: "wake", count: n) + guard !stages.isEmpty else { return out } + let sorted = stages.sorted { $0.start < $1.start } + var si = 0 + for i in 0..= sorted[sorted.count - 1].end { out[i] = sorted[sorted.count - 1].stage } + } + return out +} + +// MARK: - Stage-fraction calibration + +/// Each stage's share of the epochs, as a percentage. +/// +/// This is the quantity Cohen's kappa does NOT constrain, and the reason it is reported beside kappa +/// everywhere in this tool rather than in an appendix. PR #348 raised kappa on all three of its benchmarks +/// and was reverted 48 h later for re-scoring a healthy night from 6 % to 23 % awake: a recipe can win more +/// epochs than it loses while systematically shifting how much of the night it spends in a stage. Every +/// stage in `stageOrder` is present at 0 when unused, so a stage a recipe stops emitting shows up as a bias +/// rather than as a missing row. +func stagePercentages(_ labels: [String]) -> [String: Double] { + var pct = [String: Double](uniqueKeysWithValues: stageOrder.map { ($0, 0.0) }) + guard !labels.isEmpty else { return pct } + for l in labels { pct[l, default: 0] += 1 } + for k in pct.keys { pct[k]! = pct[k]! / Double(labels.count) * 100 } + return pct +} + +/// Signed per-stage calibration error in percentage points, `predicted% − reference%`. +func stageBias(ref: [String], pred: [String]) -> [String: Double] { + let r = stagePercentages(ref), p = stagePercentages(pred) + return [String: Double](uniqueKeysWithValues: stageOrder.map { ($0, p[$0]! - r[$0]!) }) +} + +/// Minutes from sleep onset (first non-wake epoch) to the first REM epoch. `nil` when the hypnogram has no +/// sleep or never reaches REM — a night with no REM is a distinct outcome from one whose REM arrives at +/// minute zero, and collapsing them would let a recipe that stops emitting REM look fast. +func firstRemLatencyMinutes(_ labels: [String]) -> Double? { + guard let onset = labels.firstIndex(where: { $0 != "wake" }), + let rem = labels.firstIndex(of: "rem"), rem >= onset else { return nil } + return Double(rem - onset) * epochSeconds / 60.0 +} + +/// The same latency over a grid with unscored epochs in it. The nils are skipped when looking for onset +/// and for REM, but they still occupy their slot in the index arithmetic — so the answer stays wall-clock +/// minutes. Collapsing the array first and counting the survivors would shorten the latency by however +/// much of the night the PSG technician left unscored, silently and only for the truth column. +func firstRemLatencyMinutes(_ labels: [String?]) -> Double? { + guard let onset = labels.firstIndex(where: { $0 != nil && $0 != "wake" }), + let rem = labels.firstIndex(where: { $0 == "rem" }), rem >= onset else { return nil } + return Double(rem - onset) * epochSeconds / 60.0 +} + +// MARK: - Agreement + +/// A square confusion matrix over `classes`, rows = reference, cols = prediction. +struct Confusion { + let classes: [String] + var m: [[Int]] + init(classes: [String] = stageOrder) { + self.classes = classes + m = Array(repeating: Array(repeating: 0, count: classes.count), count: classes.count) + } + mutating func add(ref: String, pred: String) { + guard let r = classes.firstIndex(of: ref), let c = classes.firstIndex(of: pred) else { return } + m[r][c] += 1 + } + mutating func merge(_ o: Confusion) { + for i in 0.. 0 else { return .nan } + let po = accuracy + var pe = 0.0 + for i in 0..= 1.0 ? .nan : (po - pe) / (1 - pe) + } + /// One-vs-rest precision, recall and F1 for `cls`. `support` is the reference count, so a class the + /// recipe never emits is distinguishable from a class the reference never contains. + func prf(_ cls: String) -> (precision: Double, recall: Double, f1: Double, support: Int) { + guard let k = classes.firstIndex(of: cls) else { return (.nan, .nan, .nan, 0) } + var tp = 0, fn = 0, fp = 0 + for i in 0.. (sens: Double, spec: Double, n: Int) { + guard let k = classes.firstIndex(of: cls) else { return (.nan, .nan, 0) } + var tp = 0, fn = 0, fp = 0, tn = 0 + for i in 0.. Confusion { + var c = Confusion(classes: classes) + for i in 0.. [String] { labels.map { $0 == "wake" ? "wake" : "sleep" } } + +/// Binary F1 for a single positive class, over aligned label arrays. The REM ablation's headline number. +func binaryF1(ref: [Bool], pred: [Bool]) -> Double { + var tp = 0, fp = 0, fn = 0 + for i in 0.. Double { v.isEmpty ? .nan : v.reduce(0, +) / Double(v.count) } +func median(_ v: [Double]) -> Double { + guard !v.isEmpty else { return .nan } + let s = v.sorted() + return s.count % 2 == 1 ? s[s.count / 2] : (s[s.count / 2 - 1] + s[s.count / 2]) / 2 +} +func sd(_ v: [Double]) -> Double { + guard v.count > 1 else { return .nan } + let m = mean(v) + return (v.reduce(0) { $0 + ($1 - m) * ($1 - m) } / Double(v.count - 1)).squareRoot() +} +func mae(_ v: [Double]) -> Double { mean(v.map { abs($0) }) } + +/// Pearson correlation and its two-sided t statistic. +func pearson(_ x: [Double], _ y: [Double]) -> (r: Double, n: Int, t: Double) { + let n = min(x.count, y.count) + guard n > 2 else { return (.nan, n, .nan) } + let mx = mean(Array(x.prefix(n))), my = mean(Array(y.prefix(n))) + var sxy = 0.0, sxx = 0.0, syy = 0.0 + for i in 0.. 0, syy > 0 else { return (.nan, n, .nan) } + let r = sxy / (sxx * syy).squareRoot() + let t = r * (Double(n - 2) / max(1e-12, 1 - r * r)).squareRoot() + return (r, n, t) +} + +/// Spearman rank correlation — reported beside Pearson because a monotone-but-curved coupling (which is +/// what a saturating penalty produces) shows up in the ranks before it shows up in the linear fit. +func spearman(_ x: [Double], _ y: [Double]) -> Double { + let n = min(x.count, y.count) + guard n > 2 else { return .nan } + func ranks(_ v: [Double]) -> [Double] { + let idx = v.indices.sorted { v[$0] < v[$1] } + var r = [Double](repeating: 0, count: v.count) + var i = 0 + while i < idx.count { + var j = i + while j + 1 < idx.count && v[idx[j + 1]] == v[idx[i]] { j += 1 } + let avg = Double(i + j) / 2 + 1 + for k in i...j { r[idx[k]] = avg } + i = j + 1 + } + return r + } + return pearson(ranks(Array(x.prefix(n))), ranks(Array(y.prefix(n)))).r +} + +func f(_ v: Double, _ w: Int = 7, _ p: Int = 1) -> String { + v.isNaN ? String(repeating: " ", count: max(0, w - 3)) + "n/a" : String(format: "%\(w).\(p)f", v) +} diff --git a/Tools/SleepPSG/Sources/sleeppsg/SleepAccel.swift b/Tools/SleepPSG/Sources/sleeppsg/SleepAccel.swift new file mode 100644 index 0000000000..12e2f8e244 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/SleepAccel.swift @@ -0,0 +1,232 @@ +import Foundation +import WhoopProtocol + +// SleepAccel — read PhysioNet `sleep-accel` v1.0.0 and present each subject as the streams a NOOP stager +// takes, plus the PSG hypnogram to score against. +// +// Walch O, Huang Y, Forger D, Goldstein C. Sleep stage prediction with raw acceleration and photo- +// plethysmography heart rate data derived from a consumer wearable device. SLEEP 42(12), zsz180 (2019). +// +// Dataset licence: Open Data Commons Attribution License v1.0. Attribution is a condition of use; see +// README.md. The dataset is NEVER committed here and is always supplied by `--dataset`. +// +// ── The mapping, and why each choice is the honest one ──────────────────────────────────────────────── +// +// This dataset is an Apple Watch: a ~50 Hz raw accelerometer and a sparse PPG heart rate. NOOP's stager +// takes a 1 Hz gravity vector, a heart-rate stream, and R-R intervals. Three decisions bridge that, and +// each one can bias a result, so each is stated rather than buried: +// +// 1. GRAVITY ← per-second mean of the raw accelerometer. `SleepStagerV2.features` already collapses its +// gravity stream to one value per second by averaging, so pre-averaging here is not an extra filter — +// it is the same operation done earlier, and it keeps 31 subjects of 50 Hz data in memory. The absolute +// SCALE still differs from a WHOOP gravity decode, and that is exactly the property the recipe was +// built not to depend on: every motion threshold is a multiple of the night's OWN median per-second +// jerk (`jerkScale`), never an absolute g. A benchmark on foreign hardware is a test of that claim. +// +// 2. R-R ← EMPTY. The dataset carries no beat-to-beat intervals, so `respRegularity` returns nil on every +// epoch, `zrg` collapses to the constant 0, and the RSA term contributes exactly 0.0 to every emission. +// The term is therefore INERT here, not corrupted: the recipe runs with one of its six inputs silent. +// Any result below is a lower bound on what the same recipe does with an R-R stream present. +// +// 3. TIME ← the dataset's seconds-from-PSG-start, offset by a fixed multiple of 30 into positive unix +// territory. Label epochs land on multiples of 30, so the recipe's wall-clock 30 s grid coincides with +// the PSG scoring grid exactly and no epoch is compared against a label it only half overlaps. +// +// The session window handed to the stager is the LABELLED span — first scored epoch to last, plus 30 s. +// NOOP's own session detection is not exercised: this measures staging, which is the only thing +// `SleepStagerV2` does (detection still comes from V1). + +/// One subject: the streams, the truth hypnogram on the same 30 s grid, and the window they share. +struct PSGSubject { + let id: String + /// Window start — a multiple of 30, aligned with truth epoch 0. + let start: Int + let end: Int + let grav: [GravitySample] + let hr: [HRSample] + /// One entry per 30 s epoch of `[start, end)`. `nil` = the epoch carries no PSG score and is excluded + /// from every metric rather than silently counted as wake. + let truth: [String?] + + var scoredEpochs: Int { truth.reduce(0) { $0 + ($1 == nil ? 0 : 1) } } + /// Duration of the scored night in minutes — the stratification variable. + var scoredMinutes: Double { Double(scoredEpochs) * 30.0 / 60.0 } +} + +enum SleepAccel { + + /// A multiple of 30 near the present, so `t = 0` maps onto a 30 s wall-clock boundary and no timestamp + /// is negative (the recipe's `((start + 29) / 30) * 30` truncates toward zero, which is only the + /// intended floor for non-negative input). + static let timeBase = 1_699_999_980 + + /// Which whole second a fractional dataset timestamp belongs to. FLOOR, not round: the recipe's + /// per-second aggregation buckets by integer second, so a sample at t = 0.9 s is inside second 0, and + /// rounding would push it into second 1 and shift the motion series against the heart-rate series by + /// up to a second. Floor is also the correct answer for the dataset's large negative times. + static func wholeSecond(_ t: Double) -> Int { Int(t.rounded(.down)) } + + /// AASM stage codes as scored in `labels/*_labeled_sleep.txt`, mapped onto NOOP's four classes. + /// N1+N2 → light and N3(+N4) → deep is the standard four-class collapse; anything else — including the + /// dataset's unscored marker — becomes `nil` and is dropped from the denominator. + static func stageForCode(_ code: Int) -> String? { + switch code { + case 0: return "wake" + case 1, 2: return "light" + case 3, 4: return "deep" + case 5: return "rem" + default: return nil // -1 unscored, and any code this collapse does not claim to handle + } + } + + /// Locate the dataset root. Accepts either the extracted directory itself or its parent, so a user who + /// unzipped without `-j` does not have to guess which level to pass. + static func resolveRoot(_ path: String) -> String? { + let fm = FileManager.default + var candidates = [path] + if let subs = try? fm.contentsOfDirectory(atPath: path) { + candidates.append(contentsOf: subs.map { (path as NSString).appendingPathComponent($0) }) + } + for c in candidates { + var isDir: ObjCBool = false + let labels = (c as NSString).appendingPathComponent("labels") + if fm.fileExists(atPath: labels, isDirectory: &isDir), isDir.boolValue { return c } + } + return nil + } + + /// Subject ids, from the label directory — a subject with no PSG hypnogram cannot be scored at all. + static func subjectIDs(root: String) -> [String] { + let dir = (root as NSString).appendingPathComponent("labels") + let files = (try? FileManager.default.contentsOfDirectory(atPath: dir)) ?? [] + return files.compactMap { name -> String? in + guard name.hasSuffix("_labeled_sleep.txt") else { return nil } + return String(name.dropLast("_labeled_sleep.txt".count)) + }.sorted() + } + + /// Load one subject. Returns nil when a required file is missing or the hypnogram has no scored epoch. + static func load(root: String, id: String) -> PSGSubject? { + let labelRows = readRows((root as NSString) + .appendingPathComponent("labels/\(id)_labeled_sleep.txt"), expecting: 2) + guard !labelRows.isEmpty else { return nil } + + // The truth grid: first scored epoch to last. Epoch times in the file are multiples of 30. + var byEpoch: [Int: String] = [:] + var minT = Int.max, maxT = Int.min + for r in labelRows { + let t = Int((r[0] / 30.0).rounded()) * 30 + guard let stage = stageForCode(Int(r[1].rounded())) else { continue } + byEpoch[t] = stage + minT = min(minT, t); maxT = max(maxT, t) + } + guard minT != Int.max else { return nil } + + let start = timeBase + minT + let end = timeBase + maxT + 30 + let nEpochs = (end - start) / 30 + var truth = [String?](repeating: nil, count: nEpochs) + for (t, s) in byEpoch { + let i = (t - minT) / 30 + if i >= 0 && i < nEpochs { truth[i] = s } + } + + // Each subject's motion and heart-rate files cover DAYS around the lab night — the earliest rows sit + // four days before lights-out. Rows outside the recipe's own read window are dropped during the + // read. This is provably loss-free rather than a convenience: `SleepStagerV2.stageSession` clips its + // inputs to exactly `[start − padLo, end + padHi)` before doing anything else, so these are rows the + // recipe would discard itself. Doing it here keeps 31 subjects of daytime 50 Hz data out of memory. + let loSec = minT - V2Recipe.padLo + let hiSec = maxT + 30 + V2Recipe.padHi + + let hr = readRows((root as NSString) + .appendingPathComponent("heart_rate/\(id)_heartrate.txt"), expecting: 2) + .compactMap { r -> HRSample? in + let bpm = Int(r[1].rounded()) + let sec = wholeSecond(r[0]) + guard bpm > 0, bpm < 300, sec >= loSec, sec < hiSec else { return nil } + return HRSample(ts: timeBase + sec, bpm: bpm) + } + .sorted { $0.ts < $1.ts } + + let grav = readAccelerationPerSecond((root as NSString) + .appendingPathComponent("motion/\(id)_acceleration.txt"), loSec: loSec, hiSec: hiSec) + + return PSGSubject(id: id, start: start, end: end, grav: grav, hr: hr, truth: truth) + } + + // MARK: - Parsing + + /// Read a whitespace/comma-separated numeric table, keeping rows that yield at least `expecting` + /// values. Byte-level because the motion files run to millions of lines each and a per-line `String` + /// allocation dominates the whole run. + static func readRows(_ path: String, expecting: Int) -> [[Double]] { + guard let data = FileManager.default.contents(atPath: path) else { return [] } + var out: [[Double]] = [] + forEachRow(data, maxValues: expecting) { vals, n in + if n >= expecting { out.append(Array(vals[0.. [GravitySample] { + guard let data = FileManager.default.contents(atPath: path) else { return [] } + var sx = [Int: Double](), sy = [Int: Double](), sz = [Int: Double](), cnt = [Int: Int]() + forEachRow(data, maxValues: 4) { vals, n in + guard n >= 4 else { return } + let t = wholeSecond(vals[0]) + if t < loSec || t >= hiSec { return } + sx[t, default: 0] += vals[1]; sy[t, default: 0] += vals[2]; sz[t, default: 0] += vals[3] + cnt[t, default: 0] += 1 + } + var out: [GravitySample] = [] + out.reserveCapacity(cnt.count) + for (t, c) in cnt { + let d = Double(c) + out.append(GravitySample(ts: timeBase + t, x: sx[t]! / d, y: sy[t]! / d, z: sz[t]! / d)) + } + return out.sorted { $0.ts < $1.ts } + } + + /// Walk `data` line by line, parsing up to `maxValues` numbers per line with `strtod`, and hand each + /// row to `body`. Separators are commas and spaces/tabs; a line whose numbers run out early is passed + /// with its count so the caller can drop it. + private static func forEachRow(_ data: Data, maxValues: Int, _ body: ([Double], Int) -> Void) { + var bytes = [UInt8](data) + bytes.append(0) + var vals = [Double](repeating: 0, count: maxValues) + bytes.withUnsafeBufferPointer { buf in + guard let raw = buf.baseAddress else { return } + let base = UnsafeRawPointer(raw).assumingMemoryBound(to: CChar.self) + let count = buf.count - 1 + var i = 0 + while i < count { + var lineEnd = i + while lineEnd < count && buf[lineEnd] != 0x0A { lineEnd += 1 } + var p = i + var n = 0 + while p < lineEnd && n < maxValues { + let c = buf[p] + // Only start a parse on something that can begin a number; otherwise `strtod` would + // skip trailing whitespace and swallow the FIRST number of the following line. + let numeric = (c >= 0x30 && c <= 0x39) || c == 0x2D || c == 0x2B || c == 0x2E + if !numeric { p += 1; continue } + var endPtr: UnsafeMutablePointer? + let v = strtod(base + p, &endPtr) + guard let e = endPtr else { break } + let consumed = UnsafePointer(e) - (base + p) + if consumed <= 0 { p += 1; continue } + vals[n] = v + n += 1 + p += consumed + } + body(vals, n) + i = lineEnd + 1 + } + } + } +} diff --git a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift new file mode 100644 index 0000000000..c1ad7667ce --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift @@ -0,0 +1,159 @@ +import Foundation + +// Variants — the recipe builds this harness compares, each one a single named edit to the shipped config. +// +// The set is not arbitrary. PR #348 changed seven things about `SleepStagerV2` at once, PR #437 reverted all +// seven 48 h later, and PR #987 re-litigated the revert component-by-component and restored exactly one of +// them. That whole argument was settled on one wearer's 36 nights, against the strap's own band state — +// an independent reference, but a proprietary black box, not truth. Each component is reproduced here as +// its own variant so PSG can be asked the same question in the same shape. +// +// Every variant differs from `RecipeConfig.shipped` in the fewest fields that express it. A variant that +// bundled two changes could not be attributed. + +struct Variant { + let name: String + let note: String + let config: RecipeConfig +} + +enum Variants { + + /// The awake transition row before PR #987, and after it. #987 forbids wake → deep and wake → REM + /// outright and gives the freed mass to the wake self-loop. It was landed on band-state evidence — band + /// kappa 0.105 → 0.118, wake sensitivity 16.0 % → 17.6 %, first-REM latency MAE 53.9 → 41.6 min — and + /// never checked against PSG, which is one of the two questions this harness exists to answer. + static let awakeRowPre987: [String: Double] = ["deep": 0.01, "rem": 0.02, "light": 0.27, "awake": 0.70] + static let awakeRowPR987: [String: Double] = ["deep": 0.0, "rem": 0.0, "light": 0.10, "awake": 0.90] + + /// The #987 comparison, whichever side of it this branch is on. + /// + /// This matters because the fork and upstream disagree about the incumbent: this repository's `main` + /// carries #987 already, upstream's does not until the PR merges, and `RecipeConfig.shipped` must track + /// whatever `SleepStagerV2` says on the branch it is compiled against (that is exactly what + /// `PortValidation` enforces). Rather than hard-code one direction and be wrong half the time, the + /// variant reads the shipped awake row and offers the OTHER one. The comparison is the same either way; + /// only the sign of the delta and the row's name flip. + static var pr987: Variant { + let shippedRow = RecipeConfig.shipped.transition["awake"] ?? awakeRowPre987 + let alreadyHas987 = shippedRow == awakeRowPR987 + var c = RecipeConfig.shipped + c.transition["awake"] = alreadyHas987 ? awakeRowPre987 : awakeRowPR987 + return Variant( + name: alreadyHas987 ? "#987 REVERTED (pre-#987 row)" : "#987 awake-transition row", + note: alreadyHas987 + ? "wake→deep/rem restored to 0.01/0.02; self-loop 0.90→0.70" + : "wake→deep/rem forbidden; wake self-loop 0.70→0.90", + config: c) + } + + /// #348 component 1 — base priors. Measured alone on band state: healthy wake fraction 9.43 % → 17.76 %, + /// i.e. the #437 blow-out reproduced. + static var p348Priors: Variant { + var c = RecipeConfig.shipped + c.priorDeep = log(0.15); c.priorAwake = log(0.34) + return Variant(name: "#348-A base priors", + note: "deep 0.18→0.15, awake 0.10→0.34", config: c) + } + + /// #348 component 2 — motion gate. A second wake channel: alone, healthy wake 9.43 % → 15.92 %. + static var p348Motion: Variant { + var c = RecipeConfig.shipped + c.jerkFloorMoveMult = 75.0; c.jerkFloorGateMult = 35.0; c.motionGateBoost = 4.0 + return Variant(name: "#348-B motion gate", + note: "move 38→75, gate 55→35, boost 2→4", config: c) + } + + /// #348 component 3 — emission coefficients. Alone: band kappa 0.105 → 0.094, but the closest call of + /// the six (it cut healthy deep 22.09 % → 15.30 % and latency MAE 53.9 → 32.1 min). + static var p348Emissions: Variant { + var c = RecipeConfig.shipped + c.deepZhv = -0.8; c.deepZhr = 0.5; c.deepZmv = -0.1 + c.remZhv = 0.8; c.remZmv = -0.4; c.remZhr = 0.4 + c.awakeZmv = 1.0; c.awakeZhv = 0.5; c.awakeZhr = 0.6 + return Variant(name: "#348-C emission coefficients", + note: "deep/rem/awake log-emission weights", config: c) + } + + /// #348 component 4 — deep gate. Alone: healthy deep 22.09 % → 25.47 %, worsening an over-call. + static var p348DeepGate: Variant { + var c = RecipeConfig.shipped + c.deepGateThresh = 0.40 + return Variant(name: "#348-D deep gate", note: "deepGateThresh 0.25→0.40", config: c) + } + + /// #348 component 5 — awake cardiac dead-zone. Alone: band kappa 0.105 → 0.099. + static var p348Deadzone: Variant { + var c = RecipeConfig.shipped + c.awakeDeadzone = 0.30 + return Variant(name: "#348-E awake dead-zone", note: "awakeDeadzone off→0.30", config: c) + } + + /// #348 component 6 — the other three transition rows. Alone: band kappa 0.105 → 0.101, healthy REM + /// 29.58 % → 33.11 %. + static var p348OtherRows: Variant { + var c = RecipeConfig.shipped + c.transition["deep"] = ["deep": 0.76, "rem": 0.012, "light": 0.216, "awake": 0.012] + c.transition["rem"] = ["deep": 0.00333, "rem": 0.92, "light": 0.06667, "awake": 0.01] + c.transition["light"] = ["deep": 0.08, "rem": 0.08, "light": 0.80, "awake": 0.04] + return Variant(name: "#348-F deep/rem/light rows", + note: "self-loops retuned, off-diagonals renormalised", config: c) + } + + /// All seven together — the build #437 reverted. Reported so the component rows can be checked for + /// additivity rather than assumed to be. + static var p348All: Variant { + var c = p348Priors.config + let m = p348Motion.config, e = p348Emissions.config, d = p348Deadzone.config + c.jerkFloorMoveMult = m.jerkFloorMoveMult; c.jerkFloorGateMult = m.jerkFloorGateMult + c.motionGateBoost = m.motionGateBoost + c.deepZhv = e.deepZhv; c.deepZhr = e.deepZhr; c.deepZmv = e.deepZmv + c.remZhv = e.remZhv; c.remZmv = e.remZmv; c.remZhr = e.remZhr + c.awakeZmv = e.awakeZmv; c.awakeZhv = e.awakeZhv; c.awakeZhr = e.awakeZhr + c.deepGateThresh = p348DeepGate.config.deepGateThresh + c.awakeDeadzone = d.awakeDeadzone + c.transition = p348OtherRows.config.transition + // #348's own awake row IS the one #987 later restored — pinned to the literal rather than to + // `pr987`, which flips direction depending on which side of #987 the branch is on. + c.transition["awake"] = awakeRowPR987 + return Variant(name: "#348 entire (as reverted by #437)", + note: "all seven components at once", config: c) + } + + /// NOT a candidate — a provenance check, and the reason it is here is worth stating. + /// + /// This harness was rebuilt to reproduce a set of numbers a previous, deleted harness measured on this + /// same dataset: kappa 0.349, REM F1 0.515, REM 20.8 % of night, median first-REM latency 142.0 min. + /// The rebuild reproduces every TRUTH-side figure exactly (26 773 scored epochs, deep 14.76 % of night, + /// truth first-REM 88.5 min) but lands the prediction side elsewhere — more REM, arriving earlier. + /// + /// Those percentages are the MEAN OF PER-SUBJECT PERCENTAGES, which is the convention the old figures + /// were reported in; pooled over all 26 773 epochs the same truth deep fraction is 13.76 % and the same + /// predicted deep fraction is 18.94 %. `main.swift` prints both tables, labelled, for exactly this + /// reason — a stage fraction quoted without its convention is a number two people will disagree about + /// while both being right. Section 6's variant table, and every `bias pp` in it, is pooled. + /// + /// One mechanism explains that whole pattern, and it is dated rather than mysterious: PR #930 replaced a + /// hard `c < 0.12 ? 3.0 : 0.0` step in the SESSION-FRACTION domain with a graded penalty measured in + /// MINUTES from sleep onset. Averaged over the epochs it covers, the graded guard is the weaker one, so + /// after #930 the recipe admits REM earlier and admits more of it. The old numbers were measured while + /// #930 was the CANDIDATE — so the incumbent they describe is the recipe as it stood BEFORE it. + /// + /// This variant runs that pre-#930 guard so the claim is checked rather than asserted. + static var preNine30Guard: Variant { + var c = RecipeConfig.shipped + c.remLatencyMode = .preNine30FractionStep + return Variant(name: "pre-#930 REM guard (provenance)", + note: "c<0.12 step in the fraction domain, single Viterbi pass", config: c) + } + + /// The shipped recipe itself, so a table row exists to difference against. + static var incumbent: Variant { + Variant(name: "incumbent (shipped)", note: "SleepStagerV2 as it ships today", config: .shipped) + } + + static var all: [Variant] { + [incumbent, pr987, p348Priors, p348Motion, p348Emissions, p348DeepGate, + p348Deadzone, p348OtherRows, p348All, preNine30Guard] + } +} diff --git a/Tools/SleepPSG/Sources/sleeppsg/main.swift b/Tools/SleepPSG/Sources/sleeppsg/main.swift new file mode 100644 index 0000000000..133cbb872b --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/main.swift @@ -0,0 +1,481 @@ +import Foundation +import StrandAnalytics +import WhoopProtocol + +// sleeppsg — score NOOP's shipped sleep stager against polysomnography. +// +// USAGE: sleeppsg --dataset /path/to/sleep-accel-1.0.0 [--section all] [--csv out.csv] [--subjects N] +// sleeppsg --section port (no dataset needed — the port check alone) +// +// The dataset lives outside this repository and is always an argument. Nothing is written to it. + +struct Args { + var dataset: String? + var section = "all" + var csv: String? + var subjects: Int? + var seed: UInt64 = PortValidation.defaultSeed +} +var a = Args() +var it = CommandLine.arguments.dropFirst().makeIterator() +while let k = it.next() { + switch k { + case "--dataset": a.dataset = it.next() + case "--section": a.section = it.next() ?? a.section + case "--csv": a.csv = it.next() + case "--subjects": a.subjects = Int(it.next() ?? "") + case "--seed": a.seed = UInt64(it.next() ?? "") ?? a.seed + case "-h", "--help": + print(""" + usage: sleeppsg --dataset [--section all|port|baseline|strata|rem|variants] + [--subjects N] [--csv ] [--seed ] + + --dataset PhysioNet sleep-accel v1.0.0, extracted. See README.md for the download step and the + attribution its licence requires. Never committed; read-only here. + --section which report to print. `port` needs no dataset. + --subjects score only the first N subject ids (a fast smoke run; full cohort is the default). + """) + exit(0) + default: FileHandle.standardError.write(Data("unknown argument \(k)\n".utf8)); exit(2) + } +} + +let wantAll = a.section == "all" +func want(_ s: String) -> Bool { wantAll || a.section == s } + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// 1. PORT VALIDATION +// ───────────────────────────────────────────────────────────────────────────────────────────────────── + +if want("port") { + print(""" + + ================================================================================ + 1. PORT VALIDATION — does this harness run the SHIPPED recipe? + ================================================================================ + Every baseline number below comes from `StrandAnalytics.SleepStagerV2.stageSession` directly: the file + that ships in the app, compiled from `Packages/`, called with no reimplementation in between. The + VARIANT rows cannot work that way — `SleepStagerV2` holds its constants as `static let`s — so they run + through `V2Recipe`, a knob-for-knob port in this tool. This section is the check that the port and the + shipped stager are the same recipe. It runs in `swift test`, and therefore in CI, with no dataset. + """) + let r = PortValidation.run(seed: a.seed) + print(""" + nights \(r.nights) (\(r.randomNights) randomised + \(r.degenerateNights) degenerate) + epoch labels \(r.matchingEpochs)/\(r.epochs) identical + verdict \(r.ok ? "PASS — the port reproduces the shipped stager exactly" : "FAIL") + """) + if !r.ok { + for d in r.divergences.prefix(20) { + print(" divergence \(d.night) epoch \(d.epochIndex): shipped=\(d.shipped) port=\(d.port)") + } + FileHandle.standardError.write(Data("port validation FAILED — variant numbers are not trustworthy\n".utf8)) + exit(1) + } +} + +if a.section == "port" { exit(0) } + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// Dataset +// ───────────────────────────────────────────────────────────────────────────────────────────────────── + +guard let datasetArg = a.dataset else { + FileHandle.standardError.write(Data(""" + --dataset is required for every section except `port`. + See Tools/SleepPSG/README.md for the download step and the attribution the licence requires. + + """.utf8)) + exit(2) +} +guard let root = SleepAccel.resolveRoot(datasetArg) else { + FileHandle.standardError.write(Data("no `labels/` directory under \(datasetArg) — is that the sleep-accel root?\n".utf8)) + exit(2) +} + +var ids = SleepAccel.subjectIDs(root: root) +if let n = a.subjects { ids = Array(ids.prefix(n)) } +// Loaded concurrently — the motion files run to tens of megabytes of text each and the read dominates the +// whole run. Order is restored from `ids` afterwards so the report is identical however the threads finish. +var loaded = [PSGSubject?](repeating: nil, count: ids.count) +loaded.withUnsafeMutableBufferPointer { buf in + let out = buf + DispatchQueue.concurrentPerform(iterations: ids.count) { i in + out[i] = SleepAccel.load(root: root, id: ids[i]) + } +} +let subjects: [PSGSubject] = loaded.compactMap { $0 } +guard !subjects.isEmpty else { + FileHandle.standardError.write(Data("no scorable subjects under \(root)\n".utf8)); exit(1) +} + +let totalScored = subjects.reduce(0) { $0 + $1.scoredEpochs } +print(""" + +================================================================================ +2. DATASET +================================================================================ +PhysioNet sleep-accel v1.0.0 — Walch, Huang, Forger & Goldstein, SLEEP 42(12) zsz180 (2019). +Licence: Open Data Commons Attribution v1.0. Root: \(root) + +subjects \(subjects.count) +PSG-scored epochs \(totalScored) +scored night length median \(f(median(subjects.map { $0.scoredMinutes }), 6, 1)) min \ +range \(f(subjects.map { $0.scoredMinutes }.min() ?? .nan, 5, 1))–\(f(subjects.map { $0.scoredMinutes }.max() ?? .nan, 5, 1)) min + +The R-R stream is EMPTY for every subject — this dataset carries no beat-to-beat intervals — so the +recipe's RSA respiration term returns nil on every epoch and contributes exactly 0.0 to every emission. +Five of the recipe's six inputs are live here; the sixth is silent, not wrong. +""") + +// Truth / prediction on the shared 30 s grid, restricted to PSG-scored epochs. +struct Scored { + let subject: PSGSubject + let truth: [String] + let pred: [String] + /// Full-grid arrays (including unscored epochs) — latency is a property of the hypnogram, not of the + /// scored subset, so it is measured before masking. + let fullPred: [String] + let fullTruth: [String?] +} + +func score(_ s: PSGSubject, using stage: (PSGSubject) -> [StageSegment]) -> Scored { + let segs = stage(s) + let pred = epochLabels(segs, start: s.start, end: s.end) + var t: [String] = [], p: [String] = [] + for i in 0.. [StageSegment] { + SleepStagerV2.stageSession(start: s.start, end: s.end, grav: s.grav, hr: s.hr, rr: [], resp: []) +} +func stageVariant(_ cfg: RecipeConfig) -> (PSGSubject) -> [StageSegment] { + { s in V2Recipe.stageSession(start: s.start, end: s.end, grav: s.grav, hr: s.hr, + rr: [], resp: [], cfg: cfg) } +} + +let shipped = subjects.map { score($0, using: stageShipped) } + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// Reporting helpers +// ───────────────────────────────────────────────────────────────────────────────────────────────────── + +struct Summary { + var conf = Confusion() + var kappa: Double { conf.kappa } + var accuracy: Double { conf.accuracy } + var predPct: [String: Double] = [:] + var truthPct: [String: Double] = [:] + var latencyPred: [Double] = [] + var latencyTruth: [Double] = [] + var latencyPairsPred: [Double] = [] + var latencyPairsTruth: [Double] = [] + /// Per-subject kappa, excluding any subject whose scored night carries only one class (kappa is + /// undefined there, and letting a NaN into the mean would quietly void the whole column). + var subjectKappa: [Double] = [] + var noRemNights = 0 +} + +func summarise(_ rows: [Scored]) -> Summary { + var s = Summary() + var predAcc = [String: Double](), truthAcc = [String: Double]() + for r in rows { + var c = Confusion() + for i in 0.. [String]) -> [String: Double] { + var all: [String] = [] + for r in rows { all.append(contentsOf: pick(r)) } + return stagePercentages(all) +} + +let sh = summarise(shipped) +let pooledPred = pooledPct(shipped) { $0.pred } +let pooledTruth = pooledPct(shipped) { $0.truth } + +if want("baseline") { + print(""" + + ================================================================================ + 3. SHIPPED SleepStagerV2 vs PSG TRUTH + ================================================================================ + Four-class agreement over every PSG-scored epoch. `SleepStagerV2` fits no parameters to data — every + coefficient is fixed a priori — so there is no train/test split to make here: the recipe sees each + subject exactly once and has never seen any of them. The leave-one-subject-out machinery in section 6 + exists for the fitted comparison models, which do need it. + + epochs scored \(sh.conf.total) + accuracy \(f(sh.accuracy * 100, 6, 2)) % + Cohen's kappa \(f(sh.kappa, 6, 3)) (per-subject mean \(f(mean(sh.subjectKappa), 6, 3)) ± \(f(sd(sh.subjectKappa), 5, 3))) + + per stage precision recall F1 support + """) + for st in stageOrder { + let p = sh.conf.prf(st) + print(" \(st.padding(toLength: 18, withPad: " ", startingAt: 0))\(f(p.precision, 8, 3)) \(f(p.recall, 8, 3)) \(f(p.f1, 8, 3)) \(p.support)") + } + print(""" + + STAGE FRACTIONS — reported as a first-class result, not an appendix. + Kappa does not constrain them: PR #348 raised kappa on all three of its benchmarks and was reverted + 48 h later for re-scoring a healthy night from 6 % to 23 % awake. + + stage predicted % truth % bias pp (pooled over all scored epochs) + """) + for st in stageOrder { + print(" \(st.padding(toLength: 10, withPad: " ", startingAt: 0))\(f(pooledPred[st]!, 10, 2))\(f(pooledTruth[st]!, 10, 2))\(f(pooledPred[st]! - pooledTruth[st]!, 12, 2))") + } + print("\n stage predicted % truth % bias pp (mean of per-subject percentages)") + for st in stageOrder { + print(" \(st.padding(toLength: 10, withPad: " ", startingAt: 0))\(f(sh.predPct[st]!, 10, 2))\(f(sh.truthPct[st]!, 10, 2))\(f(sh.predPct[st]! - sh.truthPct[st]!, 12, 2))") + } + print(""" + + FIRST-REM LATENCY — minutes from staged sleep onset to the first REM epoch. + nights with no REM in the prediction: \(sh.noRemNights) of \(shipped.count) + + predicted truth + median \(f(median(sh.latencyPred), 10, 1))\(f(median(sh.latencyTruth), 10, 1)) min + mean \(f(mean(sh.latencyPred), 10, 1))\(f(mean(sh.latencyTruth), 10, 1)) min + MAE (paired) \(f(mae(zip(sh.latencyPairsPred, sh.latencyPairsTruth).map { $0 - $1 }), 10, 1)) min (n = \(sh.latencyPairsPred.count)) + """) + + // A harness that replaces a deleted one has to say whether it is the same instrument. These are the + // figures the previous harness reported on this dataset. They are printed as a self-check, NOT as a + // target: nothing in this tool is tuned to hit them, and where they disagree the disagreement is the + // finding. See `Variants.preNine30Guard` for what explains the prediction-side gap. + let ref: [(String, Double, Double, Int)] = [ + ("subjects", Double(subjects.count), 31, 0), + ("PSG-scored epochs", Double(sh.conf.total), 26773, 0), + ("kappa (4-class)", sh.kappa, 0.349, 3), + ("REM F1", sh.conf.prf("rem").f1, 0.515, 3), + ("REM % of night", sh.predPct["rem"]!, 20.8, 2), + ("deep % predicted", sh.predPct["deep"]!, 19.25, 2), + ("deep % truth", sh.truthPct["deep"]!, 14.76, 2), + ("first-REM predicted, min", median(sh.latencyPred), 142.0, 1), + ("first-REM truth, min", median(sh.latencyTruth), 88.5, 1), + ] + print(""" + + SELF-CHECK against the previous harness's reported figures + Printed to expose disagreement, not to be matched. Nothing here is tuned to these numbers. Stage + fractions use the per-subject-mean convention, which is the one those figures were reported on. + + quantity measured previously delta + """) + for (name, got, want, dp) in ref { + print(" \(name.padding(toLength: 28, withPad: " ", startingAt: 0))\(f(got, 9, dp))\(f(want, 13, dp))\(f(got - want, 11, dp))") + } +} + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// Night-length stratification +// ───────────────────────────────────────────────────────────────────────────────────────────────────── + +if want("strata") { + print(""" + + ================================================================================ + 4. STRATIFIED BY NIGHT LENGTH + ================================================================================ + Pooling hides a coupling. On one wearer's own nights, `SleepStagerV2` emits more REM as a fraction of + sleep the longer the night runs — corr(total sleep, REM %) = +0.579 — traced to the REM-latency guard + being a fixed 60-minute penalty, which therefore covers a third of a short night and a tenth of a long + one. Whether that reproduces against PSG is a question a pooled table cannot be asked. + + Subjects split into terciles by SCORED night length. + """) + let sortedRows = shipped.sorted { $0.subject.scoredMinutes < $1.subject.scoredMinutes } + let third = max(1, sortedRows.count / 3) + let strata: [(String, [Scored])] = [ + ("short", Array(sortedRows.prefix(third))), + ("mid", Array(sortedRows.dropFirst(third).dropLast(sortedRows.count - 2 * third))), + ("long", Array(sortedRows.suffix(sortedRows.count - 2 * third))), + ] + print(" stratum n night min kappa REM% pred REM% truth deep% pred deep% truth wake% pred wake% truth") + for (name, rows) in strata where !rows.isEmpty { + let p = pooledPct(rows) { $0.pred }, t = pooledPct(rows) { $0.truth } + var c = Confusion() + for r in rows { for i in 0.. Double? { + let sleep = labels.filter { $0 != "wake" }.count + guard sleep > 0 else { return nil } + return Double(labels.filter { $0 == "rem" }.count) / Double(sleep) * 100 + } + var lens: [Double] = [], remPred: [Double] = [], remTruth: [Double] = [] + for r in shipped { + guard let rp = remPctOfSleep(r.pred), let rt = remPctOfSleep(r.truth) else { continue } + lens.append(r.subject.scoredMinutes); remPred.append(rp); remTruth.append(rt) + } + let cp = pearson(lens, remPred), ct = pearson(lens, remTruth) + print(""" + + corr(scored night length, REM % of SLEEP), n = \(cp.n) + predicted Pearson \(f(cp.r, 7, 3)) Spearman \(f(spearman(lens, remPred), 7, 3)) + PSG truth Pearson \(f(ct.r, 7, 3)) Spearman \(f(spearman(lens, remTruth), 7, 3)) + A coupling the truth also carries is physiology; one only the prediction carries is an artefact. + """) +} + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// The REM question +// ───────────────────────────────────────────────────────────────────────────────────────────────────── + +if want("rem") { + print(""" + + ================================================================================ + 5. IS REM DETECTION PHYSIOLOGY, OR IS IT THE CLOCK? + ================================================================================ + Three models fit on the same epochs, leave-one-SUBJECT-out, scored only on the held-out subject. + Clock features are the time-of-night fraction and its square; physiology features are the per-night + z-scored heart rate, HR-variability and movement plus the HR-flatness percentile — the same normalised + quantities the recipe's own emissions are built from. Decision threshold tuned on the training folds. + """) + var rows: [AblationRow] = [] + for s in subjects { + let (feats, _) = V2Recipe.stageEpochsDetailed(start: s.start, end: s.end, grav: s.grav, + hr: s.hr, rr: [], resp: []) + // Per-night normalisation, matching `stageEpochs`. + func zfun(_ vals: [Double?]) -> (Double?) -> Double { + let present = vals.compactMap { $0 } + if present.isEmpty { return { _ in 0.0 } } + let m = present.reduce(0, +) / Double(present.count) + let sd0 = (present.reduce(0.0) { $0 + ($1 - m) * ($1 - m) } / Double(present.count)).squareRoot() + let sdv = sd0 == 0 ? 1.0 : sd0 + return { v in v == nil ? 0.0 : (v! - m) / sdv } + } + let zhr = zfun(feats.map { $0.hr }), zhv = zfun(feats.map { $0.hrVar }) + let zmv = zfun(feats.map { Optional($0.moveFrac) }) + let fsorted = feats.compactMap { $0.hrFlat11 }.sorted() + func fpct(_ v: Double?) -> Double { + guard let v = v, !fsorted.isEmpty else { return 0.5 } + var lo = 0, hi = fsorted.count + while lo < hi { let mid = (lo + hi) / 2; if fsorted[mid] <= v { lo = mid + 1 } else { hi = mid } } + return Double(lo) / Double(fsorted.count) + } + for e in feats { + let idx = (e.start - s.start) / 30 + guard idx >= 0, idx < s.truth.count, let tv = s.truth[idx] else { continue } + rows.append(AblationRow( + subject: s.id, isRem: tv == "rem", + clock: [e.clock, e.clock * e.clock], + physiology: [zhr(e.hr), zhv(e.hrVar), zmv(e.moveFrac), fpct(e.hrFlat11)])) + } + } + print(" epochs in the comparison: \(rows.count) REM prevalence \(f(Double(rows.filter { $0.isRem }.count) / Double(max(1, rows.count)) * 100, 5, 1)) %\n") + let res = Ablation.run(rows) + print(" model pooled F1 pooled P pooled R mean per-subject F1") + for r in res { + print(" \(r.model.rawValue.padding(toLength: 18, withPad: " ", startingAt: 0))\(f(r.pooledF1, 8, 3)) \(f(r.pooledPrecision, 9, 3)) \(f(r.pooledRecall, 9, 3)) \(f(r.meanSubjectF1, 12, 3))") + } + if let clock = res.first(where: { $0.model == .clock }), + let phys = res.first(where: { $0.model == .physiology }), + let both = res.first(where: { $0.model == .both }) { + let subs = clock.subjectF1.keys.sorted() + let dPhys = subs.compactMap { s -> Double? in + guard let c = clock.subjectF1[s], let p = phys.subjectF1[s] else { return nil } + return p - c + } + let dBoth = subs.compactMap { s -> Double? in + guard let c = clock.subjectF1[s], let b = both.subjectF1[s] else { return nil } + return b - c + } + print(""" + + what physiology adds over the clock, per subject + physiology-only − clock-only mean \(f(mean(dPhys), 7, 3)) positive in \(dPhys.filter { $0 > 0 }.count)/\(dPhys.count) subjects + both − clock-only mean \(f(mean(dBoth), 7, 3)) positive in \(dBoth.filter { $0 > 0 }.count)/\(dBoth.count) subjects + pooled: physiology-only − clock-only \(f(phys.pooledF1 - clock.pooledF1, 7, 3)), both − clock-only \(f(both.pooledF1 - clock.pooledF1, 7, 3)) + + For reference, the SHIPPED recipe's own REM F1 on these epochs is \(f(sh.conf.prf("rem").f1, 6, 3)) — it fits + nothing, so it is not comparable to a fitted model's held-out score and is printed only for scale. + """) + } +} + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// Variants +// ───────────────────────────────────────────────────────────────────────────────────────────────────── + +if want("variants") { + print(""" + + ================================================================================ + 6. RECIPE VARIANTS AGAINST PSG + ================================================================================ + Every row is the shipped recipe with ONE named change. Deltas are against the incumbent row. + `wake bias` and `deep bias` are predicted minus truth in percentage points — the #437 guard, which + kappa does not carry. + """) + var csvLines = [(["variant", "epochs", "kappa", "accuracy"] + + stageOrder.flatMap { ["\($0)_precision", "\($0)_recall", "\($0)_f1", "\($0)_pct", "\($0)_bias_pp"] } + + ["median_rem_latency_min", "rem_latency_mae_min"]).joined(separator: ",")] + var base: (kappa: Double, remF1: Double, wakeSens: Double)? + // Wake SENSITIVITY is in the table beside kappa because it is the quantity #987 was landed on against + // the strap's band state (16.0 % → 17.6 %). Printing it here is what makes that claim checkable against + // truth rather than merely restated. + print(" variant kappa Δκ REM F1 wakeSens wake% bias deep% bias REM% bias latMAE") + for v in Variants.all { + let rows = v.name.hasPrefix("incumbent") + ? shipped + : subjects.map { score($0, using: stageVariant(v.config)) } + let s = summarise(rows) + let p = pooledPct(rows) { $0.pred }, t = pooledPct(rows) { $0.truth } + let latMae = mae(zip(s.latencyPairsPred, s.latencyPairsTruth).map { $0 - $1 }) + let wakeSens = s.conf.prf("wake").recall * 100 + if base == nil { base = (s.kappa, s.conf.prf("rem").f1, wakeSens) } + let b = base! + print(" \(v.name.padding(toLength: 32, withPad: " ", startingAt: 0))\(f(s.kappa, 6, 3)) \(f(s.kappa - b.kappa, 6, 3)) \(f(s.conf.prf("rem").f1, 6, 3)) \(f(wakeSens, 6, 2))\(f(wakeSens - b.wakeSens, 6, 2)) \(f(p["wake"]!, 6, 2))\(f(p["wake"]! - t["wake"]!, 6, 2)) \(f(p["deep"]!, 6, 2))\(f(p["deep"]! - t["deep"]!, 6, 2)) \(f(p["rem"]!, 6, 2))\(f(p["rem"]! - t["rem"]!, 6, 2)) \(f(latMae, 6, 1))") + var row = [v.name, "\(s.conf.total)", "\(s.kappa)", "\(s.accuracy)"] + for st in stageOrder { + let m = s.conf.prf(st) + row += ["\(m.precision)", "\(m.recall)", "\(m.f1)", "\(p[st]!)", "\(p[st]! - t[st]!)"] + } + row += ["\(median(s.latencyPred))", "\(latMae)"] + csvLines.append(row.joined(separator: ",")) + } + print(""" + + A component is an improvement here only if it raises kappa AND does not blow out a stage fraction. + That conjunction is the lesson of #348 → #437: kappa alone certified a build that mis-called a + healthy night's wake fraction by 17 percentage points. + """) + if let path = a.csv { + try? csvLines.joined(separator: "\n").write(toFile: path, atomically: true, encoding: .utf8) + print(" per-variant CSV written to \(path)") + } +} + +print("") diff --git a/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift new file mode 100644 index 0000000000..25bad9b120 --- /dev/null +++ b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift @@ -0,0 +1,295 @@ +import XCTest +import WhoopProtocol +@testable import sleeppsg + +/// The dataset reader and the leave-one-subject-out machinery, on synthetic fixtures written to a temp +/// directory. No PhysioNet download, no health data, no network — these run in CI. +final class SleepAccelTests: XCTestCase { + + private func makeFixture(labels: String, heartRate: String, motion: String) throws -> String { + let root = NSTemporaryDirectory() + "sleeppsg-fixture-" + UUID().uuidString + let fm = FileManager.default + for d in ["labels", "heart_rate", "motion"] { + try fm.createDirectory(atPath: root + "/" + d, withIntermediateDirectories: true) + } + try labels.write(toFile: root + "/labels/9001_labeled_sleep.txt", atomically: true, encoding: .utf8) + try heartRate.write(toFile: root + "/heart_rate/9001_heartrate.txt", atomically: true, encoding: .utf8) + try motion.write(toFile: root + "/motion/9001_acceleration.txt", atomically: true, encoding: .utf8) + addTeardownBlock { try? fm.removeItem(atPath: root) } + return root + } + + /// The AASM collapse. N1 and N2 both become light and N3 becomes deep, which is the standard + /// four-class reduction; anything the collapse does not claim to handle becomes nil and leaves the + /// denominator rather than being counted as wake. + func testStageCodeCollapse() { + XCTAssertEqual(SleepAccel.stageForCode(0), "wake") + XCTAssertEqual(SleepAccel.stageForCode(1), "light") + XCTAssertEqual(SleepAccel.stageForCode(2), "light") + XCTAssertEqual(SleepAccel.stageForCode(3), "deep") + XCTAssertEqual(SleepAccel.stageForCode(4), "deep") + XCTAssertEqual(SleepAccel.stageForCode(5), "rem") + XCTAssertNil(SleepAccel.stageForCode(-1)) + XCTAssertNil(SleepAccel.stageForCode(7)) + } + + /// The load path: label grid, time base, and the per-second collapse of the raw accelerometer. + func testLoadAlignsTruthEpochsToTheStagerGrid() throws { + // Four scored epochs, one of them unscored (-1), and a REM epoch last. + let labels = "0 0\n30 1\n60 -1\n90 5\n" + // Heart rate is comma-separated in this dataset; motion is space-separated. + let hr = "0.0,60\n1.5,61\n119.0,58\n" + // Two raw samples inside second 0 and one inside second 1, so the per-second mean is checkable. + // 0.90 belongs to second 0 — the bucket is a FLOOR, so a sample late in a second is not pushed + // into the next one, which would shift motion against heart rate by up to a second. + let motion = "0.10 0.0 0.0 1.0\n0.90 1.0 0.0 0.0\n1.20 0.0 1.0 0.0\n" + let root = try makeFixture(labels: labels, heartRate: hr, motion: motion) + + XCTAssertEqual(SleepAccel.resolveRoot(root), root) + XCTAssertEqual(SleepAccel.subjectIDs(root: root), ["9001"]) + let s = try XCTUnwrap(SleepAccel.load(root: root, id: "9001")) + + XCTAssertEqual(s.start, SleepAccel.timeBase) + XCTAssertEqual(s.end, SleepAccel.timeBase + 120) + XCTAssertEqual(s.truth, ["wake", "light", nil, "rem"]) + XCTAssertEqual(s.scoredEpochs, 3) + XCTAssertEqual(s.scoredMinutes, 1.5, accuracy: 1e-12) + + // The window start is a multiple of 30, so the recipe's wall-clock epoch grid coincides with the + // PSG scoring grid and no epoch is compared against a label it only half overlaps. + XCTAssertEqual(s.start % 30, 0) + + XCTAssertEqual(s.hr.map { $0.ts }, [SleepAccel.timeBase, SleepAccel.timeBase + 1, SleepAccel.timeBase + 119]) + XCTAssertEqual(s.hr.map { $0.bpm }, [60, 61, 58]) + + XCTAssertEqual(s.grav.count, 2) + XCTAssertEqual(s.grav[0].ts, SleepAccel.timeBase) + XCTAssertEqual(s.grav[0].x, 0.5, accuracy: 1e-12, "second 0 averages its two raw samples") + XCTAssertEqual(s.grav[0].z, 0.5, accuracy: 1e-12) + XCTAssertEqual(s.grav[1].ts, SleepAccel.timeBase + 1) + XCTAssertEqual(s.grav[1].y, 1.0, accuracy: 1e-12) + } + + /// Rows outside the recipe's own read window are dropped during the read. `SleepStagerV2.stageSession` + /// clips its inputs to `[start − 330, end + 390)` before doing anything, so this discards only rows the + /// recipe provably never touches — which matters because each subject's files span DAYS around the lab + /// night and the daytime 50 Hz motion would otherwise dominate memory. + func testSamplesOutsideTheRecipesReadWindowAreDropped() throws { + // One labelled epoch at t = 0, so the read window is [-330, 420). + let root = try makeFixture( + labels: "0 2\n", + heartRate: "-400.0,70\n-100.0,60\n0.0,61\n500.0,80\n", + motion: "-400.0 0 0 1\n-100.0 0 0 1\n0.0 0 0 1\n500.0 0 0 1\n") + let s = try XCTUnwrap(SleepAccel.load(root: root, id: "9001")) + XCTAssertEqual(s.hr.map { $0.ts - SleepAccel.timeBase }, [-100, 0]) + XCTAssertEqual(s.grav.map { $0.ts - SleepAccel.timeBase }, [-100, 0]) + } + + /// Negative times (the dataset records before PSG lights-out) must not produce negative timestamps: + /// the recipe's `((start + 29) / 30) * 30` truncates toward zero, which is only the intended floor on + /// non-negative input. + func testPreRecordingSamplesStayPositive() throws { + // t = −300.5 is before lights-out but still inside the recipe's read window (−330 onward), so it + // survives clipping and exercises the negative-time path. Floor puts it in second −301. + let root = try makeFixture(labels: "0 0\n30 2\n", + heartRate: "-300.5,70\n0.0,60\n", + motion: "-300.5 0.0 0.0 1.0\n0.0 0.0 0.0 1.0\n") + let s = try XCTUnwrap(SleepAccel.load(root: root, id: "9001")) + XCTAssertTrue(s.hr.allSatisfy { $0.ts > 0 }) + XCTAssertTrue(s.grav.allSatisfy { $0.ts > 0 }) + XCTAssertEqual(s.hr.first?.ts, SleepAccel.timeBase - 301) + XCTAssertEqual(s.grav.first?.ts, SleepAccel.timeBase - 301) + } + + /// The byte-level row parser: mixed separators, blank lines, a trailing line with no newline, and a + /// short row that must be dropped rather than half-read. + func testRowParserToleratesRealFileShapes() throws { + let root = try makeFixture(labels: "0 0\n30 5\n", + heartRate: "0.0,60\n\n10.0 61\n20.0,62", + motion: "0.0 1.0 0.0 0.0\n1.0 0.0\n2.0 0.0 0.0 1.0") + let s = try XCTUnwrap(SleepAccel.load(root: root, id: "9001")) + XCTAssertEqual(s.hr.map { $0.bpm }, [60, 61, 62], "blank line skipped, both separators accepted") + XCTAssertEqual(s.grav.count, 2, "the 2-value motion row is dropped, not padded with garbage") + } + + /// Passing the directory that CONTAINS the extracted dataset must work too — the published zip expands + /// into a long-named subdirectory nobody types by hand. + func testResolveRootAcceptsTheParentDirectory() throws { + let root = try makeFixture(labels: "0 0\n", heartRate: "0.0,60\n", motion: "0.0 0 0 1\n") + let parent = (root as NSString).deletingLastPathComponent + XCTAssertEqual(SleepAccel.resolveRoot(parent + "/" + (root as NSString).lastPathComponent), root) + } +} + +final class AblationTests: XCTestCase { + + /// A subject's own epochs must never appear in the fold that scores it. This is the leak the whole + /// design exists to prevent — consecutive 30 s epochs from one night are nearly identical, so an + /// epoch-level split lets any model memorise the subject and score near-perfectly. + func testHeldOutSubjectIsNeverInItsOwnTrainingFold() { + // Two subjects with OPPOSITE physiology→REM mappings and no clock signal at all. A model that + // leaked the held-out subject's rows would score well; one that generalises cannot, because the + // training fold teaches it exactly the wrong sign. + var rows: [AblationRow] = [] + for i in 0..<200 { + let x = Double(i % 2) + rows.append(AblationRow(subject: "A", isRem: x > 0.5, + clock: [0.5, 0.25], physiology: [x, 0, 0, 0.5])) + rows.append(AblationRow(subject: "B", isRem: x < 0.5, + clock: [0.5, 0.25], physiology: [x, 0, 0, 0.5])) + } + let res = Ablation.run(rows) + let phys = res.first { $0.model == .physiology }! + XCTAssertLessThan(phys.pooledF1, 0.55, + "physiology F1 near 1.0 here would mean the held-out subject leaked into training") + } + + /// The fitter has to actually fit: a cleanly separable signal must be recovered, or a null result from + /// the clock model would be indistinguishable from a broken optimiser. + func testLogisticFitRecoversASeparableSignal() { + var X: [[Double]] = [], y: [Bool] = [] + for i in 0..<400 { + let v = Double(i) / 400.0 + X.append([v]); y.append(v > 0.5) + } + let w = Ablation.logisticFit(X: X, y: y) + XCTAssertGreaterThan(w[1], 1.0, "the coefficient must be positive and substantial") + XCTAssertLessThan(Ablation.sigmoid(Ablation.dot(w, [0.1])), 0.2) + XCTAssertGreaterThan(Ablation.sigmoid(Ablation.dot(w, [0.9])), 0.8) + } + + /// The threshold is tuned on the training folds, which matters because REM is a minority class: a + /// fixed 0.5 systematically under-predicts it, and unequally across models with different calibration. + func testThresholdTuningBeatsAFixedHalfOnAMinorityClass() { + // 10 % positives whose scores never reach 0.5 — a fixed cut scores 0, a tuned cut recovers them. + var scores: [Double] = [], truth: [Bool] = [] + for i in 0..<100 { + let pos = i < 10 + scores.append(pos ? 0.40 : 0.05); truth.append(pos) + } + let thr = Ablation.bestF1Threshold(scores: scores, truth: truth) + XCTAssertLessThanOrEqual(thr, 0.40) + XCTAssertEqual(binaryF1(ref: truth, pred: scores.map { $0 >= 0.5 }), 0.0, accuracy: 1e-12) + XCTAssertEqual(binaryF1(ref: truth, pred: scores.map { $0 >= thr }), 1.0, accuracy: 1e-12) + } + + func testLinearSolveMatchesAKnownSystem() { + // [[2,1],[1,3]] x = [5,10] → x = [1,3]. + let x = Ablation.solve([[2, 1], [1, 3]], [5, 10])! + XCTAssertEqual(x[0], 1.0, accuracy: 1e-9) + XCTAssertEqual(x[1], 3.0, accuracy: 1e-9) + XCTAssertNil(Ablation.solve([[1, 2], [2, 4]], [1, 2]), "a singular system returns nil, not NaNs") + } +} + +final class VariantsTests: XCTestCase { + + /// Each variant must differ from the shipped recipe in exactly the fields its name claims. A variant + /// that quietly moved a second constant could not be attributed, which is the entire point of measuring + /// #348's components one at a time. + func testEachVariantIsASingleNamedChange() { + let s = RecipeConfig.shipped + + // #987 is a two-sided comparison: this repository's main carries it, upstream's does not, and the + // variant offers whichever row the branch is NOT on. Both directions are asserted so the test does + // not have to know which branch it is compiled on. + var c = Variants.pr987.config + let shippedRow = s.transition["awake"]! + XCTAssertTrue(shippedRow == Variants.awakeRowPre987 || shippedRow == Variants.awakeRowPR987, + "the shipped awake row is neither side of #987 — RecipeConfig.shipped is stale") + XCTAssertEqual(c.transition["awake"]!, + shippedRow == Variants.awakeRowPR987 ? Variants.awakeRowPre987 : Variants.awakeRowPR987, + "the #987 variant must offer the row the branch is NOT on") + c.transition = s.transition + XCTAssertEqual(c, s, "#987 must change nothing but the awake transition row") + + var p = Variants.p348Priors.config + XCTAssertEqual(p.priorDeep, log(0.15), accuracy: 1e-12) + XCTAssertEqual(p.priorAwake, log(0.34), accuracy: 1e-12) + p.priorDeep = s.priorDeep; p.priorAwake = s.priorAwake + XCTAssertEqual(p, s) + + var m = Variants.p348Motion.config + m.jerkFloorMoveMult = s.jerkFloorMoveMult; m.jerkFloorGateMult = s.jerkFloorGateMult + m.motionGateBoost = s.motionGateBoost + XCTAssertEqual(m, s) + + var d = Variants.p348DeepGate.config + d.deepGateThresh = s.deepGateThresh + XCTAssertEqual(d, s) + + var z = Variants.p348Deadzone.config + z.awakeDeadzone = s.awakeDeadzone + XCTAssertEqual(z, s) + + var o = Variants.p348OtherRows.config + XCTAssertEqual(o.transition["awake"]!, s.transition["awake"]!, "the awake row is #987's, not this one's") + o.transition = s.transition + XCTAssertEqual(o, s) + + var g = Variants.preNine30Guard.config + XCTAssertEqual(g.remLatencyMode, .preNine30FractionStep) + g.remLatencyMode = s.remLatencyMode + XCTAssertEqual(g, s, "the provenance variant must change nothing but the guard's shape") + } + + /// The shipped mode is the graded one, and the pre-#930 step must actually reach the lattice — a + /// provenance variant that staged identically to the incumbent would explain nothing about the gap it + /// is there to explain. + func testPreNine30GuardActuallyChangesTheHypnogram() { + XCTAssertEqual(RecipeConfig.shipped.remLatencyMode, .gradedFromOnset) + var differed = false + for n in PortValidation.corpus().prefix(6) { + let a = V2Recipe.stageSession(start: n.start, end: n.end, grav: n.grav, hr: n.hr, + rr: n.rr, resp: n.resp, cfg: .shipped) + let b = V2Recipe.stageSession(start: n.start, end: n.end, grav: n.grav, hr: n.hr, + rr: n.rr, resp: n.resp, cfg: Variants.preNine30Guard.config) + XCTAssertFalse(a.isEmpty) + if epochLabels(a, start: n.start, end: n.end) != epochLabels(b, start: n.start, end: n.end) { + differed = true + break + } + } + XCTAssertTrue(differed, "the pre-#930 guard staged identically on every night — mode is not wired") + } + + /// The shipped dead-zone is off, and `dz` must then be the identity — that is what makes + /// `RecipeConfig.shipped` able to reproduce an emission the shipped file writes without any dead-zone + /// at all. + func testDeadzoneIsTheIdentityWhenDisabled() { + XCTAssertEqual(RecipeConfig.shipped.awakeDeadzone, 0.0) + for z in [-3.0, -0.1, 0.0, 0.29, 5.0] { + XCTAssertEqual(V2Recipe.dz(z, 0.0), z, accuracy: 1e-12) + } + XCTAssertEqual(V2Recipe.dz(0.2, 0.3), 0.0, accuracy: 1e-12) + XCTAssertEqual(V2Recipe.dz(0.5, 0.3), 0.2, accuracy: 1e-12) + XCTAssertEqual(V2Recipe.dz(-0.5, 0.3), -0.2, accuracy: 1e-12) + } + + /// Transition rows must stay probability distributions, in every variant. A row that no longer sums to + /// 1 would silently reweight the whole lattice rather than express the change its name claims. + func testEveryVariantsTransitionRowsSumToOne() { + for v in Variants.all { + for (from, row) in v.config.transition { + XCTAssertEqual(row.values.reduce(0, +), 1.0, accuracy: 1e-9, + "\(v.name): row \(from) does not sum to 1") + } + } + } + + /// `#348 entire` must be the union of the components, not an independently typed eighth build. + func testEntire348IsTheUnionOfItsComponents() { + let all = Variants.p348All.config + XCTAssertEqual(all.priorDeep, Variants.p348Priors.config.priorDeep) + XCTAssertEqual(all.priorAwake, Variants.p348Priors.config.priorAwake) + XCTAssertEqual(all.jerkFloorMoveMult, Variants.p348Motion.config.jerkFloorMoveMult) + XCTAssertEqual(all.motionGateBoost, Variants.p348Motion.config.motionGateBoost) + XCTAssertEqual(all.deepZhv, Variants.p348Emissions.config.deepZhv) + XCTAssertEqual(all.awakeZhr, Variants.p348Emissions.config.awakeZhr) + XCTAssertEqual(all.deepGateThresh, Variants.p348DeepGate.config.deepGateThresh) + XCTAssertEqual(all.awakeDeadzone, Variants.p348Deadzone.config.awakeDeadzone) + XCTAssertEqual(all.transition["light"]!, Variants.p348OtherRows.config.transition["light"]!) + XCTAssertEqual(all.transition["awake"]!, Variants.awakeRowPR987, + "#348's awake row is the one #987 later restored") + } +} diff --git a/Tools/SleepPSG/Tests/sleeppsgTests/PortValidationTests.swift b/Tools/SleepPSG/Tests/sleeppsgTests/PortValidationTests.swift new file mode 100644 index 0000000000..2065380b68 --- /dev/null +++ b/Tools/SleepPSG/Tests/sleeppsgTests/PortValidationTests.swift @@ -0,0 +1,92 @@ +import XCTest +import StrandAnalytics +import WhoopProtocol +@testable import sleeppsg + +/// The test that keeps this harness honest. +/// +/// `V2Recipe` re-states `SleepStagerV2` with its constants exposed, because that is the only way to score a +/// variant of the recipe without rebuilding the app. The failure mode of every such port is silent drift: +/// someone changes a coefficient in `Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift`, +/// nobody changes `RecipeConfig.shipped`, and the benchmark keeps printing confident numbers about a recipe +/// that no longer exists. +/// +/// So equivalence is asserted here, in CI, on every PR — not verified once by reading the two files side by +/// side. Nothing in this file needs the PhysioNet dataset, a database, or any health data. +final class PortValidationTests: XCTestCase { + + /// Every epoch label, on every night in the corpus, from both paths. + func testPortReproducesShippedStagerExactly() { + let r = PortValidation.run() + XCTAssertEqual(r.nights, 55, "corpus size changed — 48 randomised + 7 degenerate expected") + XCTAssertEqual(r.randomNights, 48) + XCTAssertEqual(r.degenerateNights, 7) + XCTAssertGreaterThan(r.epochs, 10_000, "corpus too small to be evidence of anything") + if !r.divergences.isEmpty { + let detail = r.divergences.prefix(5) + .map { "\($0.night) epoch \($0.epochIndex): shipped=\($0.shipped) port=\($0.port)" } + .joined(separator: "; ") + XCTFail(""" + RecipeConfig.shipped no longer reproduces SleepStagerV2 — \ + \(r.epochs - r.matchingEpochs) of \(r.epochs) epoch labels differ. \ + Update RecipeConfig.shipped to match the shipped file. First divergences: \(detail) + """) + } + XCTAssertEqual(r.matchingEpochs, r.epochs) + } + + /// The corpus must actually reach the branches it claims to. A degenerate case that silently produced + /// no epochs would make the headline "identical on N epochs" true and meaningless. + func testDegenerateNightsAreNonTrivial() { + for n in PortValidation.degenerateNights() where n.name != "degenerate-no-coverage" { + let segs = SleepStagerV2.stageSession(start: n.start, end: n.end, grav: n.grav, + hr: n.hr, rr: n.rr, resp: n.resp) + XCTAssertFalse(segs.isEmpty, "\(n.name) staged to nothing") + let labels = epochLabels(segs, start: n.start, end: n.end) + XCTAssertGreaterThan(labels.count, 0, "\(n.name) produced no epochs") + } + // The no-coverage night is the one that MUST collapse to the single-segment fallback. + let empty = PortValidation.degenerateNights().first { $0.name == "degenerate-no-coverage" }! + let segs = SleepStagerV2.stageSession(start: empty.start, end: empty.end, grav: [], + hr: [], rr: [], resp: []) + XCTAssertEqual(segs.count, 1) + XCTAssertEqual(segs.first?.stage, "light") + } + + /// The corpus is regenerated from a fixed seed, so a divergence is always a code change and never a + /// reroll. Compared on the GENERATED nights rather than on staged output: staging the corpus is the + /// expensive part, and determinism is a property of the generator. + func testCorpusIsDeterministic() { + func signature(_ ns: [SyntheticNight]) -> [String] { + ns.map { "\($0.name)|\($0.start)|\($0.end)|\($0.grav.count)|\($0.hr.count)|\($0.rr.count)|" + + "\($0.hr.first?.bpm ?? -1)|\($0.hr.last?.bpm ?? -1)|\($0.rr.last?.rrMs ?? -1)" } + } + XCTAssertEqual(signature(PortValidation.corpus()), signature(PortValidation.corpus())) + XCTAssertNotEqual(signature(PortValidation.corpus()), signature(PortValidation.corpus(seed: 12345)), + "a different seed must produce a different corpus") + } + + /// Equivalence must hold on ANY corpus, not only the pinned one — otherwise the seed is load-bearing + /// and the check is a coincidence. A small off-seed corpus makes that point without paying for a + /// second full run. + func testPortIsEquivalentOnAnUnpinnedCorpusToo() { + let r = PortValidation.run(seed: 0xA11CE, randomNights: 6) + XCTAssertGreaterThan(r.epochs, 500) + XCTAssertEqual(r.matchingEpochs, r.epochs, "\(r.divergences.first.map(String.init(describing:)) ?? "")") + } + + /// The randomised nights have to exercise the R-R / RSA path, which the PhysioNet dataset cannot: + /// sleep-accel carries no beat-to-beat intervals, so without this the respiration term would be + /// unvalidated in the port everywhere it matters. + func testCorpusExercisesTheRespirationTerm() { + let withRR = PortValidation.corpus().filter { !$0.rr.isEmpty } + XCTAssertGreaterThan(withRR.count, 10, "too few nights carry R-R for the RSA term to be pinned") + let n = withRR[0] + let feats = V2Recipe.features(start: n.start, end: n.end, + grav: n.grav.sorted { $0.ts < $1.ts }, + hr: n.hr.sorted { $0.ts < $1.ts }, + rr: n.rr.sortedByTsStable()) + XCTAssertTrue(feats.contains { $0.respReg != nil }, + "no epoch produced a respiration-regularity value — the RSA term is never reached") + } +} diff --git a/Tools/SleepPSG/Tests/sleeppsgTests/ScoringTests.swift b/Tools/SleepPSG/Tests/sleeppsgTests/ScoringTests.swift new file mode 100644 index 0000000000..ebe1b9abc7 --- /dev/null +++ b/Tools/SleepPSG/Tests/sleeppsgTests/ScoringTests.swift @@ -0,0 +1,141 @@ +import XCTest +import StrandAnalytics +@testable import sleeppsg + +/// Scoring primitives, pinned against values computed by hand. +/// +/// These definitions are shared, in spirit, with `Tools/SleepBench/Sources/sleepbench/Metrics.swift`: the +/// two harnesses score the same recipe against different references, and a kappa that meant something +/// slightly different in each would make their numbers incomparable precisely when someone needs to compare +/// them. Pinning the arithmetic is how that stays true without the two files being one file. +final class ScoringTests: XCTestCase { + + // MARK: - Agreement + + func testKappaOnAHandComputedMatrix() { + // 10 epochs: 4 wake (3 right), 6 light (5 right, 1 called wake). + let ref = ["wake", "wake", "wake", "wake", "light", "light", "light", "light", "light", "light"] + let pred = ["wake", "wake", "wake", "light", "light", "light", "light", "light", "light", "wake"] + let c = confusion(ref: ref, pred: pred) + // po = 8/10. Marginals: ref wake .4 light .6; pred wake .4 light .6 → pe = .16 + .36 = .52. + // kappa = (.8 − .52) / (1 − .52) = .28 / .48 = 0.58333… + XCTAssertEqual(c.accuracy, 0.8, accuracy: 1e-12) + XCTAssertEqual(c.kappa, 0.28 / 0.48, accuracy: 1e-12) + } + + func testPerfectAndChanceAgreement() { + let labels = ["wake", "light", "deep", "rem", "light", "light"] + XCTAssertEqual(confusion(ref: labels, pred: labels).kappa, 1.0, accuracy: 1e-12) + // Every epoch called light against a 50/50 reference: accuracy = pe, so kappa is exactly 0. + let ref = ["light", "wake", "light", "wake"] + let all = ["light", "light", "light", "light"] + XCTAssertEqual(confusion(ref: ref, pred: all).kappa, 0.0, accuracy: 1e-12) + } + + func testPerStageF1IsHarmonicMeanOfPrecisionAndRecall() { + // REM: 2 true REM epochs, 1 recovered; 1 false REM. precision 1/2, recall 1/2, F1 1/2. + let ref = ["rem", "rem", "light", "light"] + let pred = ["rem", "light", "rem", "light"] + let p = confusion(ref: ref, pred: pred).prf("rem") + XCTAssertEqual(p.precision, 0.5, accuracy: 1e-12) + XCTAssertEqual(p.recall, 0.5, accuracy: 1e-12) + XCTAssertEqual(p.f1, 0.5, accuracy: 1e-12) + XCTAssertEqual(p.support, 2) + } + + /// A recipe that stops emitting a stage must score 0 on it, not n/a — a silent nil would let a + /// collapsed stage disappear from a summary table instead of announcing itself. + func testStageNeverEmittedScoresZeroNotNaN() { + let p = confusion(ref: ["rem", "rem", "light"], pred: ["light", "light", "light"]).prf("rem") + XCTAssertTrue(p.precision.isNaN, "precision is undefined when nothing was predicted") + XCTAssertEqual(p.recall, 0.0, accuracy: 1e-12) + XCTAssertEqual(p.f1, 0.0, accuracy: 1e-12) + XCTAssertEqual(p.support, 2) + } + + func testBinaryF1MatchesTheConfusionMatrixVersion() { + let ref = ["rem", "rem", "light", "light", "rem"] + let pred = ["rem", "light", "rem", "light", "rem"] + XCTAssertEqual(binaryF1(ref: ref.map { $0 == "rem" }, pred: pred.map { $0 == "rem" }), + confusion(ref: ref, pred: pred).prf("rem").f1, accuracy: 1e-12) + } + + // MARK: - Stage fractions + + func testStagePercentagesAlwaysCarryEveryStage() { + let pct = stagePercentages(["light", "light", "wake", "wake"]) + XCTAssertEqual(Set(pct.keys), Set(stageOrder)) + XCTAssertEqual(pct["light"]!, 50, accuracy: 1e-12) + XCTAssertEqual(pct["wake"]!, 50, accuracy: 1e-12) + XCTAssertEqual(pct["deep"]!, 0, accuracy: 1e-12) + XCTAssertEqual(pct["rem"]!, 0, accuracy: 1e-12) + } + + /// The #348 → #437 lesson, as an assertion: a change can raise kappa AND wreck a stage fraction, so a + /// harness that reports only kappa cannot see the regression that reverted that PR. + func testKappaCanRiseWhileWakeFractionBlowsOut() { + let ref = Array(repeating: "light", count: 90) + Array(repeating: "wake", count: 10) + // Incumbent: calls only 6 epochs wake, gets 4 of them right. + var before = Array(repeating: "light", count: 100) + for i in 90..<94 { before[i] = "wake" } + for i in 0..<2 { before[i] = "wake" } + // Candidate: calls 23 epochs wake, gets 9 of the 10 real ones — kappa up, wake fraction ruined. + var after = Array(repeating: "light", count: 100) + for i in 90..<99 { after[i] = "wake" } + for i in 0..<14 { after[i] = "wake" } + let kBefore = confusion(ref: ref, pred: before).kappa + let kAfter = confusion(ref: ref, pred: after).kappa + XCTAssertGreaterThan(kAfter, kBefore, "the candidate must win on kappa for this test to mean anything") + XCTAssertEqual(stageBias(ref: ref, pred: before)["wake"]!, -4, accuracy: 1e-9) + XCTAssertEqual(stageBias(ref: ref, pred: after)["wake"]!, 13, accuracy: 1e-9) + } + + // MARK: - Latency + + func testFirstRemLatencyIsMeasuredFromOnsetNotWindowStart() { + // 4 wake epochs, then light, then REM at index 6 → 2 epochs after onset → 1.0 min. + let labels = ["wake", "wake", "wake", "wake", "light", "light", "rem", "light"] + XCTAssertEqual(firstRemLatencyMinutes(labels)!, 1.0, accuracy: 1e-12) + } + + func testNoRemIsNilRatherThanZero() { + XCTAssertNil(firstRemLatencyMinutes(["wake", "light", "light", "deep"])) + XCTAssertNil(firstRemLatencyMinutes(["wake", "wake"])) + } + + /// Unscored epochs keep their slot in the arithmetic. Collapsing them first would shorten the truth + /// latency by however much of the night the technician left unscored — and only the truth column, + /// which is exactly the kind of asymmetry that makes a prediction look better than it is. + func testLatencyOverAGridWithUnscoredEpochs() { + let sparse: [String?] = ["wake", nil, nil, "light", nil, "light", "rem"] + // onset at index 3, REM at index 6 → 3 epochs → 1.5 min, unscored slots included. + XCTAssertEqual(firstRemLatencyMinutes(sparse)!, 1.5, accuracy: 1e-12) + // Collapsing first would have given 2 epochs → 1.0 min. Pin that this is NOT what happens. + XCTAssertNotEqual(firstRemLatencyMinutes(sparse.compactMap { $0 })!, 1.5) + XCTAssertNil(firstRemLatencyMinutes([nil, nil, "wake"] as [String?])) + } + + // MARK: - Epoch expansion + + func testEpochLabelsTileTheWindowFromSegments() { + let segs = [StageSegment(start: 0, end: 60, stage: "wake"), + StageSegment(start: 60, end: 150, stage: "light"), + StageSegment(start: 150, end: 180, stage: "rem")] + XCTAssertEqual(epochLabels(segs, start: 0, end: 180), + ["wake", "wake", "light", "light", "light", "rem"]) + } + + // MARK: - Descriptive statistics + + func testPearsonAndSpearmanOnAMonotoneButCurvedRelation() { + let x = [1.0, 2, 3, 4, 5] + let y = [1.0, 4, 9, 16, 25] + XCTAssertEqual(spearman(x, y), 1.0, accuracy: 1e-12, "ranks are perfectly monotone") + XCTAssertLessThan(pearson(x, y).r, 1.0, "the linear fit is not perfect") + XCTAssertGreaterThan(pearson(x, y).r, 0.95) + } + + func testSpearmanHandlesTies() { + XCTAssertEqual(spearman([1.0, 2, 2, 3], [1.0, 2, 2, 3]), 1.0, accuracy: 1e-12) + } +}