From c3b86db4cb135c47678234694203e12f0b39115a Mon Sep 17 00:00:00 2001 From: vishk23 <119831996+vishk23@users.noreply.github.com> Date: Thu, 30 Jul 2026 19:01:45 -0400 Subject: [PATCH 1/4] Tools/SleepPSG: score the shipped sleep stager against PSG truth, reproducibly MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Every reference NOOP's sleep work can currently reach is either the recipe's own output handed back to it or another vendor's black box. Tools/SleepBench says so in its own header. A harness that scored SleepStagerV2 against real polysomnography existed once, produced the single most decisive result this project has had — that REM detection measures physiology and not the clock — and then vanished, because it lived in a session scratchpad. The next sleep investigation had to run on one instrument instead of two. This is that harness, in the repository, with a CI job. WHAT IT IS. `sleeppsg --dataset ` replays SleepStagerV2 over PhysioNet sleep-accel (Walch, Huang, Forger & Goldstein, SLEEP 42(12) zsz180, 2019) — 31 subjects of wrist accelerometer and heart rate recorded alongside human-scored PSG hypnograms — and scores it epoch-for-epoch. Read-only, dataset path always an argument, dataset never committed. README.md carries the download step and the attribution both the Open Data Commons licence and the MIT companion code require. THE PORT, AND WHY IT CANNOT DRIFT. Baseline numbers come from StrandAnalytics.SleepStagerV2.stageSession itself — the shipped file, no reimplementation. Variant numbers cannot: the recipe holds its constants as `static let`s and keeps Epoch/features() internal, so asking "what would this score with one transition row changed" needs the recipe restated with knobs. RecipePort.swift does that, and PortValidation measures the equivalence rather than asserting it: 48 randomised nights (off-grid starts, channel dropout, arousals, motion bursts, unsorted input, R-R present so the RSA term is exercised — the dataset cannot exercise it) plus 7 degenerate cases (all-HR-missing, zero-variance HR, 1- and 2-epoch nights, saturated motion gate, no coverage at all, motion absent). Every epoch label from both paths must agree. It runs in `swift test`, needs no dataset, and names the night and epoch when it fails. A failure there is never flakiness — it means RecipeConfig.shipped and SleepStagerV2.swift have diverged. WHAT IT REPORTS, AND WHY IN THAT SHAPE. - Stage fractions are a first-class result, not an appendix. Kappa does not constrain them: #348 improved kappa on all three of its benchmarks and #437 reverted it 48 h later for re-scoring a healthy night from 6 % to 23 % awake. A harness reporting only kappa could not have caught that, and ScoringTests pins a case where kappa rises while the wake fraction blows out, so the point survives as an assertion. - Splits are leave-one-SUBJECT-out, never epoch-level. Consecutive 30 s epochs from one night are nearly identical; a random epoch split lets any model memorise the subject. AblationTests fails if a held-out subject's rows can reach its own training fold. - Night length is a stratification variable, not a footnote. A REM-fraction-vs-duration coupling measured +0.579 on one wearer's nights; the report gives terciles and the same correlation for PSG TRUTH beside the prediction, because a coupling truth also carries is physiology and one only the prediction carries is an artefact. Wired into the `tools:` matrix that #943 added, alongside SleepBench and Backfill, so the equivalence check actually runs. Limits are stated in the README rather than buried: sleep-accel carries no beat-to-beat intervals, so the RSA respiration term is inert here (nil on every epoch, contributing exactly 0.0) and results are a lower bound on the same recipe with R-R present; the accelerometer is an Apple Watch's, not a WHOOP gravity decode, which tests rather than violates the recipe's night-relative motion thresholds; only staging is exercised, not session detection; and n = 31. --- .github/workflows/swift-packages.yml | 8 + Tools/SleepPSG/Package.resolved | 14 + Tools/SleepPSG/Package.swift | 33 ++ Tools/SleepPSG/README.md | 153 ++++++ .../SleepPSG/Sources/sleeppsg/Ablation.swift | 179 +++++++ .../Sources/sleeppsg/PortValidation.swift | 303 ++++++++++++ .../Sources/sleeppsg/RecipePort.swift | 446 +++++++++++++++++ Tools/SleepPSG/Sources/sleeppsg/Scoring.swift | 229 +++++++++ .../Sources/sleeppsg/SleepAccel.swift | 232 +++++++++ .../SleepPSG/Sources/sleeppsg/Variants.swift | 112 +++++ Tools/SleepPSG/Sources/sleeppsg/main.swift | 449 ++++++++++++++++++ .../DatasetAndAblationTests.swift | 260 ++++++++++ .../sleeppsgTests/PortValidationTests.swift | 92 ++++ .../Tests/sleeppsgTests/ScoringTests.swift | 141 ++++++ 14 files changed, 2651 insertions(+) create mode 100644 Tools/SleepPSG/Package.resolved create mode 100644 Tools/SleepPSG/Package.swift create mode 100644 Tools/SleepPSG/README.md create mode 100644 Tools/SleepPSG/Sources/sleeppsg/Ablation.swift create mode 100644 Tools/SleepPSG/Sources/sleeppsg/PortValidation.swift create mode 100644 Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift create mode 100644 Tools/SleepPSG/Sources/sleeppsg/Scoring.swift create mode 100644 Tools/SleepPSG/Sources/sleeppsg/SleepAccel.swift create mode 100644 Tools/SleepPSG/Sources/sleeppsg/Variants.swift create mode 100644 Tools/SleepPSG/Sources/sleeppsg/main.swift create mode 100644 Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift create mode 100644 Tools/SleepPSG/Tests/sleeppsgTests/PortValidationTests.swift create mode 100644 Tools/SleepPSG/Tests/sleeppsgTests/ScoringTests.swift 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..47901a09f8 --- /dev/null +++ b/Tools/SleepPSG/README.md @@ -0,0 +1,153 @@ +# 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. + +`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. + +--- + +## 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..20a1fb31b1 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift @@ -0,0 +1,179 @@ +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. + static func run(_ rows: [AblationRow]) -> [AblationResult] { + let subjects = Array(Set(rows.map { $0.subject })).sorted() + var out: [AblationResult] = [] + for model in AblationModel.allCases { + var pooledRef: [Bool] = [], pooledPred: [Bool] = [] + var perSubject: [String: Double] = [:] + for held in subjects { + let train = rows.filter { $0.subject != held } + let test = rows.filter { $0.subject == held } + guard !train.isEmpty, !test.isEmpty else { continue } + let X = train.map { model.features($0) } + let y = train.map { $0.isRem } + let w = logisticFit(X: X, y: y) + let trainP = X.map { sigmoid(dot(w, $0)) } + let thr = bestF1Threshold(scores: trainP, truth: y) + let testX = test.map { model.features($0) } + let pred = testX.map { sigmoid(dot(w, $0)) >= thr } + let ref = test.map { $0.isRem } + 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] { + let p = (X.first?.count ?? 0) + 1 + var w = [Double](repeating: 0, count: p) + guard !X.isEmpty, p > 1 else { return w } + for _ 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"]!, + ] + let pr = cyclePrior(f.clock, .infinity, cfg) + 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) + } + + 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..95cbf92257 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift @@ -0,0 +1,112 @@ +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 { + + /// PR #987 (upstream, open at the time of writing): forbid wake → deep and wake → REM outright and give + /// the freed mass to the wake self-loop. 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. + static var pr987: Variant { + var c = RecipeConfig.shipped + c.transition["awake"] = ["deep": 0.0, "rem": 0.0, "light": 0.10, "awake": 0.90] + return Variant(name: "#987 awake-transition row", + note: "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 + c.transition["awake"] = pr987.config.transition["awake"] + return Variant(name: "#348 entire (as reverted by #437)", + note: "all seven components at once", 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] + } +} diff --git a/Tools/SleepPSG/Sources/sleeppsg/main.swift b/Tools/SleepPSG/Sources/sleeppsg/main.swift new file mode 100644 index 0000000000..70982badd5 --- /dev/null +++ b/Tools/SleepPSG/Sources/sleeppsg/main.swift @@ -0,0 +1,449 @@ +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)) + """) +} + +// ───────────────────────────────────────────────────────────────────────────────────────────────────── +// 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,rem_f1,deep_f1,wake_f1,wake_pct,wake_bias_pp,deep_pct,deep_bias_pp,rem_pct,rem_bias_pp,median_rem_latency_min,rem_latency_mae_min"] + var base: (kappa: Double, remF1: Double, wakeBias: Double, deepBias: Double, latMae: Double)? + print(" variant kappa Δκ REM F1 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 }) + if base == nil { + base = (s.kappa, s.conf.prf("rem").f1, p["wake"]! - t["wake"]!, p["deep"]! - t["deep"]!, latMae) + } + 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(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))") + csvLines.append([v.name, "\(s.conf.total)", "\(s.kappa)", "\(s.accuracy)", + "\(s.conf.prf("rem").f1)", "\(s.conf.prf("deep").f1)", "\(s.conf.prf("wake").f1)", + "\(p["wake"]!)", "\(p["wake"]! - t["wake"]!)", + "\(p["deep"]!)", "\(p["deep"]! - t["deep"]!)", + "\(p["rem"]!)", "\(p["rem"]! - t["rem"]!)", + "\(median(s.latencyPred))", "\(latMae)"].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..cb07dd880d --- /dev/null +++ b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift @@ -0,0 +1,260 @@ +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 { + let root = try makeFixture(labels: "0 0\n30 2\n", + heartRate: "-600.0,70\n0.0,60\n", + motion: "-600.0 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 - 600) + } + + /// 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 + + var c = Variants.pr987.config + XCTAssertEqual(c.transition["awake"]!["deep"]!, 0.0) + XCTAssertEqual(c.transition["awake"]!["rem"]!, 0.0) + XCTAssertEqual(c.transition["awake"]!["awake"]!, 0.90) + 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) + } + + /// 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.pr987.config.transition["awake"]!) + } +} 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) + } +} From 1a29f85cb26141fd03bbabd72208ef7cffc3c365 Mon Sep 17 00:00:00 2001 From: vishk23 <119831996+vishk23@users.noreply.github.com> Date: Thu, 30 Jul 2026 19:28:05 -0400 Subject: [PATCH 2/4] SleepPSG: print a self-check against the old harness's figures, and a variant that falsifies the easy answer MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The rebuilt harness has to be able to say whether it is the SAME instrument as the one it replaces, so section 3 now prints the previous harness's reported figures beside its own. Nothing is tuned to them; the point is to expose disagreement, not to close it. Every truth-side figure reproduces exactly — 31 subjects, 26 773 PSG-scored epochs, deep 14.76 % of night, truth median first-REM latency 88.5 min. Predicted deep lands at 19.24 % against a reported 19.25 %, and four-class kappa at 0.356 against 0.349. Three independent numbers from PR #348's OWN Walch column also reproduce, which is the stronger check because they were produced by a different (Python) implementation of the same recipe: that PR reported the pre-#348 incumbent at kappa 0.356 (measured here: 0.356), mean deep 19.2 % (19.24 %) and wake recall 33 % (30.4 %), and the full #348 build at deep 9.9 % (10.63 %). The REM prediction does not reproduce: 27.0 % of night against a reported 20.8 %, REM F1 0.569 against 0.515, median first-REM 91.5 min against 142.0. There was an obvious explanation and it is WRONG. #930 replaced a `c < 0.12 ? 3.0 : 0` step in the session -fraction domain with a graded penalty in minutes from sleep onset; the old figures were reported while #930 was the candidate, so the incumbent they describe predates it, and the graded guard is the weaker one averaged over the epochs it covers — which would produce exactly this "more REM, arriving earlier" pattern. `Variants.preNine30Guard` runs that guard. It moves nothing: kappa 0.356 -> 0.356, REM F1 0.569 -> 0.569, REM 27.02 % -> 26.97 %. On 8-hour lab nights the guard barely binds in either form. The hypothesis is falsified and is kept in the tree because the variant is what falsified it. Two constraints on whatever the residual turns out to be, both readable off the same table: it is not the session window (`clock` drives the deep prior and the REM ramp alike, and predicted deep reproduces to 0.01 pp), and it is not the latency guard. One coincidence recorded as an observation, not a conclusion: the reported "REM F1 0.515" is exactly this harness's measured REM PRECISION. Also here: - wake SENSITIVITY joins the variant table, because that is the quantity #987 was landed on against the strap's band state, and a claim you cannot read off the table is a claim you cannot check; - the variant CSV carries precision/recall/F1/fraction/bias for all four stages instead of three F1s; - the ablation fits IRLS over one flat design matrix per model instead of refiltering row structs and reallocating a feature array per row per iteration — same arithmetic, ~93x fewer allocations, which is the difference between a minute and an afternoon in an unoptimised build; - the dataset reader clips to the recipe's own read window during the read. Each subject's files span DAYS around the lab night (the earliest rows sit four days before lights-out) and stageSession clips to [start-330, end+390) before doing anything, so this drops only rows the recipe would discard itself; - second bucketing is a floor, not a round, so a sample at t = 0.9 s stays in second 0 rather than shifting the motion series against the heart-rate series by up to a second. 33 tests, 0 failures. Full report runs in 28 s on the 31-subject cohort. --- Tools/SleepPSG/README.md | 38 +++++++ .../SleepPSG/Sources/sleeppsg/Ablation.swift | 101 ++++++++++++------ .../Sources/sleeppsg/RecipePort.swift | 34 +++++- .../SleepPSG/Sources/sleeppsg/Variants.swift | 23 +++- Tools/SleepPSG/Sources/sleeppsg/main.swift | 58 +++++++--- .../DatasetAndAblationTests.swift | 34 +++++- 6 files changed, 239 insertions(+), 49 deletions(-) diff --git a/Tools/SleepPSG/README.md b/Tools/SleepPSG/README.md index 47901a09f8..4495f78c62 100644 --- a/Tools/SleepPSG/README.md +++ b/Tools/SleepPSG/README.md @@ -135,6 +135,44 @@ Each row is the shipped recipe with **one** named change: PR #987's awake-transi fraction biases, because a component is an improvement only if it wins the first without wrecking the second. +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 **truth-side** figure reproduces exactly: 31 subjects, **26 773** PSG-scored epochs, deep **14.76 %** +of night, truth median first-REM latency **88.5 min**. Predicted deep lands at **19.24 %** against a +reported 19.25 %, and four-class kappa at **0.356** against a reported 0.349. + +The prediction side of REM does not match: more REM (27.0 % vs 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 %. 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, 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 diff --git a/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift b/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift index 20a1fb31b1..7e68f2ced2 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/Ablation.swift @@ -59,24 +59,38 @@ struct AblationResult { 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 = test.map { $0.isRem } + var trainIdx: [Int] = [], testIdx: [Int] = [] + trainIdx.reserveCapacity(rows.count) + for i in rows.indices { + if rows[i].subject == held { testIdx.append(i) } else { trainIdx.append(i) } + } + guard !trainIdx.isEmpty, !testIdx.isEmpty else { continue } + let w = logisticFit(X: X, p: p, rowIdx: trainIdx, y: y) + let thr = bestF1Threshold(scores: trainIdx.map { sigmoid(dotFlat(w, X, $0, p)) }, + truth: trainIdx.map { y[$0] }) + let pred = testIdx.map { sigmoid(dotFlat(w, X, $0, p)) >= thr } + let ref = testIdx.map { y[$0] } perSubject[held] = binaryF1(ref: ref, pred: pred) pooledRef.append(contentsOf: ref) pooledPred.append(contentsOf: pred) @@ -108,36 +122,63 @@ enum Ablation { return s } - /// Newton/IRLS with a small ridge penalty. The ridge is there for conditioning, not for tuning: the - /// feature counts here are 2–6 and the penalty is fixed across all three models, so it cannot advantage - /// one of them. Returns `[intercept, coefficients…]`. - static func logisticFit(X: [[Double]], y: [Bool], ridge: Double = 1e-4, iterations: Int = 30) -> [Double] { - let p = (X.first?.count ?? 0) + 1 + /// `w · X[row]` over the flat design matrix (whose column 0 is the intercept). + static func dotFlat(_ w: [Double], _ X: [Double], _ row: Int, _ p: Int) -> Double { + var s = 0.0 + let base = row * p + for j in 0..

[Double] { var w = [Double](repeating: 0, count: p) - guard !X.isEmpty, p > 1 else { return w } + 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]? { diff --git a/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift b/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift index 83a50b6762..44fe52d816 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift @@ -33,6 +33,17 @@ import WhoopProtocol /// The defaults are DUPLICATED from the shipped file rather than read from it, because `SleepStagerV2` /// keeps them internal to `StrandAnalytics`. Duplication that can drift is normally a defect; here it is /// pinned by `PortValidation`, which fails the build the moment the two disagree on any label. +/// How the REM emission is held back around sleep onset. +enum RemLatencyMode: Equatable { + /// #930, shipped: a penalty of `remLatencyPenalty` log-odds at sleep ONSET, decaying linearly to zero + /// `remLatencyMinutes` later. Onset is itself a staging output, so the recipe runs Viterbi twice — once + /// with the guard off to find onset, once with it re-based on that onset. + case gradedFromOnset + /// What #930 replaced: a hard step `c < 0.12 ? penalty : 0`, in the FRACTION domain and measured from + /// the window start rather than from sleep onset, with a single Viterbi pass. + case preNine30FractionStep +} + struct RecipeConfig: Equatable { // Population base rates, as log-priors. var priorLight: Double @@ -69,6 +80,11 @@ struct RecipeConfig: Equatable { var remLatencyMinutes: Double var onsetSustainedEpochs: Int + /// Which REM-latency guard the recipe runs. Shipped is `.gradedFromOnset` (#930). The pre-#930 step is + /// kept because the reference numbers this harness was rebuilt to reproduce were measured against it — + /// see `Variants.preNine30Guard`. It is a diagnostic, not a candidate. + var remLatencyMode: RemLatencyMode = .gradedFromOnset + /// The shipped recipe, as of `Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift`. static let shipped = RecipeConfig( priorLight: log(0.50), priorDeep: log(0.18), priorRem: log(0.22), priorAwake: log(0.10), @@ -85,7 +101,8 @@ struct RecipeConfig: Equatable { "light": ["deep": 0.06, "rem": 0.06, "light": 0.85, "awake": 0.03], "awake": ["deep": 0.01, "rem": 0.02, "light": 0.27, "awake": 0.70], ], - remLatencyPenalty: 3.0, remLatencyMinutes: 60.0, onsetSustainedEpochs: 10) + remLatencyPenalty: 3.0, remLatencyMinutes: 60.0, onsetSustainedEpochs: 10, + remLatencyMode: .gradedFromOnset) var baseLogPrior: [String: Double] { ["light": priorLight, "deep": priorDeep, "rem": priorRem, "awake": priorAwake] @@ -429,13 +446,26 @@ enum V2Recipe { "light": prior["light"]!, "awake": cfg.awakeZmv * zmvv + awakeCardiac + prior["awake"]!, ] - let pr = cyclePrior(f.clock, .infinity, cfg) + // 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 { diff --git a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift index 95cbf92257..2e26cbae8e 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift @@ -100,6 +100,27 @@ enum Variants { 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. + /// + /// 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) @@ -107,6 +128,6 @@ enum Variants { static var all: [Variant] { [incumbent, pr987, p348Priors, p348Motion, p348Emissions, p348DeepGate, - p348Deadzone, p348OtherRows, p348All] + p348Deadzone, p348OtherRows, p348All, preNine30Guard] } } diff --git a/Tools/SleepPSG/Sources/sleeppsg/main.swift b/Tools/SleepPSG/Sources/sleeppsg/main.swift index 70982badd5..133cbb872b 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/main.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/main.swift @@ -265,6 +265,33 @@ if want("baseline") { 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))") + } } // ───────────────────────────────────────────────────────────────────────────────────────────────────── @@ -412,9 +439,14 @@ if want("variants") { `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,rem_f1,deep_f1,wake_f1,wake_pct,wake_bias_pp,deep_pct,deep_bias_pp,rem_pct,rem_bias_pp,median_rem_latency_min,rem_latency_mae_min"] - var base: (kappa: Double, remF1: Double, wakeBias: Double, deepBias: Double, latMae: Double)? - print(" variant kappa Δκ REM F1 wake% bias deep% bias REM% bias latMAE") + 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 @@ -422,17 +454,17 @@ if want("variants") { 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 }) - if base == nil { - base = (s.kappa, s.conf.prf("rem").f1, p["wake"]! - t["wake"]!, p["deep"]! - t["deep"]!, latMae) - } + 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(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))") - csvLines.append([v.name, "\(s.conf.total)", "\(s.kappa)", "\(s.accuracy)", - "\(s.conf.prf("rem").f1)", "\(s.conf.prf("deep").f1)", "\(s.conf.prf("wake").f1)", - "\(p["wake"]!)", "\(p["wake"]! - t["wake"]!)", - "\(p["deep"]!)", "\(p["deep"]! - t["deep"]!)", - "\(p["rem"]!)", "\(p["rem"]! - t["rem"]!)", - "\(median(s.latencyPred))", "\(latMae)"].joined(separator: ",")) + 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(""" diff --git a/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift index cb07dd880d..f2d514c613 100644 --- a/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift +++ b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift @@ -89,13 +89,16 @@ final class SleepAccelTests: XCTestCase { /// 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: "-600.0,70\n0.0,60\n", - motion: "-600.0 0.0 0.0 1.0\n0.0 0.0 0.0 1.0\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 - 600) + 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 @@ -217,6 +220,31 @@ final class VariantsTests: XCTestCase { 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 From ac079d7ba6d1bd5378e4b616edb43c71d9d8506d Mon Sep 17 00:00:00 2001 From: vishk23 <119831996+vishk23@users.noreply.github.com> Date: Thu, 30 Jul 2026 19:31:36 -0400 Subject: [PATCH 3/4] SleepPSG: make the #987 comparison work from either side of the PR MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `RecipeConfig.shipped` has to describe `SleepStagerV2` as compiled on the branch it sits on — that is the whole basis of the port-validation check. The fork and upstream currently disagree about one row of the transition matrix: this repository's main carries PR #987 (wake→deep and wake→rem zeroed), upstream's does not until the PR merges. Hard-coding the variant as "apply #987" is therefore wrong on one of the two branches, and would have shipped a benchmark whose headline comparison is a no-op against itself. `Variants.pr987` now reads the shipped awake row and offers the OTHER one, naming itself accordingly ("#987 awake-transition row" upstream, "#987 REVERTED (pre-#987 row)" here). Same comparison, same pair of rows; only the sign of the delta and the label flip. The test asserts both directions rather than assuming one, so it passes wherever it is compiled, and it fails if the shipped row is neither side of #987 — which would mean `RecipeConfig.shipped` had gone stale for a different reason. `#348 entire` pins the awake row to the literal instead of to `pr987`, since #348's own row is the one #987 later restored and that is a fixed historical fact, not a branch-relative one. 33 tests, 0 failures. --- .../Sources/sleeppsg/RecipePort.swift | 6 ++++ .../SleepPSG/Sources/sleeppsg/Variants.swift | 34 +++++++++++++++---- .../DatasetAndAblationTests.swift | 15 +++++--- 3 files changed, 44 insertions(+), 11 deletions(-) diff --git a/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift b/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift index 44fe52d816..4c481169fb 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/RecipePort.swift @@ -99,6 +99,12 @@ struct RecipeConfig: Equatable { "deep": ["deep": 0.86, "rem": 0.007, "light": 0.126, "awake": 0.007], "rem": ["deep": 0.005, "rem": 0.88, "light": 0.10, "awake": 0.015], "light": ["deep": 0.06, "rem": 0.06, "light": 0.85, "awake": 0.03], + // THE AWAKE ROW IS THE ONE LINE THAT LEGITIMATELY DIFFERS BETWEEN BRANCHES. PR #987 zeroes + // wake→deep and wake→rem; this fork's `main` carries it and upstream's does not until the PR + // merges. `RecipeConfig.shipped` must always describe `SleepStagerV2` AS COMPILED on this + // branch — `PortValidation` is what enforces that, and it will fail loudly on the wrong value. + // `Variants.pr987` reads this row and offers the other one, so the #987 comparison works from + // either side without a second edit. "awake": ["deep": 0.01, "rem": 0.02, "light": 0.27, "awake": 0.70], ], remLatencyPenalty: 3.0, remLatencyMinutes: 60.0, onsetSustainedEpochs: 10, diff --git a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift index 2e26cbae8e..3a612c49e7 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift @@ -19,14 +19,32 @@ struct Variant { enum Variants { - /// PR #987 (upstream, open at the time of writing): forbid wake → deep and wake → REM outright and give - /// the freed mass to the wake self-loop. 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. + /// 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"] = ["deep": 0.0, "rem": 0.0, "light": 0.10, "awake": 0.90] - return Variant(name: "#987 awake-transition row", - note: "wake→deep/rem forbidden; wake self-loop 0.70→0.90", config: c) + 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 %, @@ -95,7 +113,9 @@ enum Variants { c.deepGateThresh = p348DeepGate.config.deepGateThresh c.awakeDeadzone = d.awakeDeadzone c.transition = p348OtherRows.config.transition - c.transition["awake"] = pr987.config.transition["awake"] + // #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) } diff --git a/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift index f2d514c613..25bad9b120 100644 --- a/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift +++ b/Tools/SleepPSG/Tests/sleeppsgTests/DatasetAndAblationTests.swift @@ -190,10 +190,16 @@ final class VariantsTests: XCTestCase { 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 - XCTAssertEqual(c.transition["awake"]!["deep"]!, 0.0) - XCTAssertEqual(c.transition["awake"]!["rem"]!, 0.0) - XCTAssertEqual(c.transition["awake"]!["awake"]!, 0.90) + 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") @@ -283,6 +289,7 @@ final class VariantsTests: XCTestCase { 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.pr987.config.transition["awake"]!) + XCTAssertEqual(all.transition["awake"]!, Variants.awakeRowPR987, + "#348's awake row is the one #987 later restored") } } From 8804c2c8f724ac5fb703054f18a571836bdc8e8f Mon Sep 17 00:00:00 2001 From: vishk23 <119831996+vishk23@users.noreply.github.com> Date: Thu, 30 Jul 2026 19:59:09 -0400 Subject: [PATCH 4/4] SleepPSG: say which denominator every stage fraction uses MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The tool has always printed both conventions, labelled, one table above the other. The README quoted one of them in the self-check and the other in the variant section without ever saying so, which read as two different answers to the same question: section 3, pooled over all scored epochs deep truth 13.76 % predicted 18.94 % bias +5.18 pp section 3, mean of per-subject percentages deep truth 14.76 % predicted 19.24 % bias +4.48 pp Both are right and both are measured on the same 26 773 epochs. The denominator WITHIN a night is identical either way — all scored epochs, wake included. What differs is the weighting ACROSS the 31 subjects: pooled gives a 490-minute night more weight than a 208-minute one, the per-subject mean gives every subject one vote. For deep that is a ~1 pp gap, which is large enough to look like a contradiction and small enough to go unnoticed. The self-check is stated in the per-subject mean because that is the convention the previous harness reported in, and reproducing it on any other convention would not be reproducing it. Section 6's variant table and every `bias pp` column in it are pooled, which is why the incumbent's deep bias reads +5.18 pp there. Neither changes; both are now named at the point of use. Adds a convention table to the README, labels every percentage in the self-check section, and records the same thing on `Variants.preNine30Guard`, whose docstring quotes the per-subject figures. Also un-rounds "REM 27.0 %" to 26.96 % — 27.0 was ambiguous between the two conventions (26.96 per-subject, 27.02 pooled) precisely where the distinction matters. Documentation only: no measurement, no threshold and no reported number changes. Verified against the full 31-subject cohort locally; `swift test` in Tools/SleepPSG 33 tests, 0 failures. --- Tools/SleepPSG/README.md | 39 +++++++++++++++---- .../SleepPSG/Sources/sleeppsg/Variants.swift | 6 +++ 2 files changed, 38 insertions(+), 7 deletions(-) diff --git a/Tools/SleepPSG/README.md b/Tools/SleepPSG/README.md index 4495f78c62..45b5d98f7a 100644 --- a/Tools/SleepPSG/README.md +++ b/Tools/SleepPSG/README.md @@ -100,6 +100,22 @@ improved kappa on all three of its benchmarks and was reverted 48 hours later fo 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. @@ -133,7 +149,8 @@ avoid. 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. +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 @@ -148,18 +165,25 @@ harness reported on this same dataset, so anyone can see at a glance whether the 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 night, truth median first-REM latency **88.5 min**. Predicted deep lands at **19.24 %** against a -reported 19.25 %, and four-class kappa at **0.356** against a reported 0.349. +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 (27.0 % vs 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 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 %. On 8-hour lab nights the guard barely binds either way. The +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. @@ -167,7 +191,8 @@ What the residual is remains **unresolved**, and the tool is not tuned to close 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, so `clock` is the same quantity in both harnesses. + 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 diff --git a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift index 3a612c49e7..c1ad7667ce 100644 --- a/Tools/SleepPSG/Sources/sleeppsg/Variants.swift +++ b/Tools/SleepPSG/Sources/sleeppsg/Variants.swift @@ -127,6 +127,12 @@ enum Variants { /// 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