diff --git a/Tools/SleepBench/Package.swift b/Tools/SleepBench/Package.swift index 53143e72c7..94a395c7af 100644 --- a/Tools/SleepBench/Package.swift +++ b/Tools/SleepBench/Package.swift @@ -19,5 +19,8 @@ let package = Package( ], targets: [ .executableTarget(name: "sleepbench", dependencies: ["StrandAnalytics", "WhoopProtocol"]), + // The scoring primitives are pure functions over label arrays, so they are unit-testable without a + // database. `swift test` here needs no health data and no `--db` argument. + .testTarget(name: "sleepbenchTests", dependencies: ["sleepbench"]), ] ) diff --git a/Tools/SleepBench/Sources/sleepbench/Metrics.swift b/Tools/SleepBench/Sources/sleepbench/Metrics.swift index 3a4f04e8f8..dd2af6de99 100644 --- a/Tools/SleepBench/Sources/sleepbench/Metrics.swift +++ b/Tools/SleepBench/Sources/sleepbench/Metrics.swift @@ -48,6 +48,45 @@ func onsetAndFinalWake(_ labels: [String]) -> (onset: Int?, final: Int?) { (labels.firstIndex(where: { $0 != "wake" }), labels.lastIndex(where: { $0 != "wake" })) } +// MARK: - Stage-fraction calibration + +/// Each stage's share of the night, as a percentage of the session's epochs. +/// +/// This is the quantity Cohen's kappa does NOT constrain. Kappa scores whether the predicted label at an +/// epoch matches the reference label at that epoch; a recipe can raise kappa while systematically shifting +/// how much of the night it calls wake, because the epochs it newly gets right can outnumber the epochs it +/// newly mislabels. PR #348 did exactly that (kappa up on all three benchmarks, held-out gap −0.027) and +/// PR #437 reverted it two days later for re-scoring a healthy night from 6% to 23% awake. Every stage in +/// `stageOrder` is present in the result, at 0 when the hypnogram never uses it, so a stage a recipe stops +/// emitting altogether 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%`. +/// Positive means the recipe spends MORE of the night at that stage than the human did. +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 staged sleep onset (the first non-wake epoch) to the first REM epoch. +/// +/// `nil` when the hypnogram has no sleep at all or never reaches REM — a night with no REM is a distinct +/// outcome from a night whose REM arrives at minute zero, and collapsing the two would let a recipe that +/// stops emitting REM look like one with a short latency. Callers report the `nil` count separately. +/// REM before onset cannot occur (REM is not wake, so onset is at or before it), but the guard is kept so +/// a caller passing a hand-built label array cannot produce a negative latency. +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 +} + // MARK: - Agreement /// A square confusion matrix over `classes`, rows = reference, cols = prediction. @@ -124,6 +163,17 @@ func sd(_ v: [Double]) -> Double { let m = mean(v) return (v.reduce(0) { $0 + ($1 - m) * ($1 - m) } / Double(v.count - 1)).squareRoot() } +/// Nearest-rank percentile, `p` in 0...1. Used for the p10/p90 spread of a stage-fraction or latency +/// distribution, where a mean alone hides the tail a regression shows up in first. +func percentile(_ v: [Double], _ p: Double) -> Double { + guard !v.isEmpty else { return .nan } + let s = v.sorted() + let i = Int((Double(s.count) - 1) * p + 0.5) + return s[min(s.count - 1, max(0, i))] +} +/// Mean absolute value — the unsigned companion to a signed bias, so a recipe that is unbiased on average +/// because it over- and under-calls in equal measure cannot pass as well calibrated. +func mae(_ v: [Double]) -> Double { mean(v.map { abs($0) }) } /// Pearson correlation, plus the two-sided t statistic that goes with it. func pearson(_ x: [Double], _ y: [Double]) -> (r: Double, n: Int, t: Double) { let n = min(x.count, y.count) diff --git a/Tools/SleepBench/Sources/sleepbench/main.swift b/Tools/SleepBench/Sources/sleepbench/main.swift index 8b851c992f..1207006723 100644 --- a/Tools/SleepBench/Sources/sleepbench/main.swift +++ b/Tools/SleepBench/Sources/sleepbench/main.swift @@ -6,6 +6,7 @@ import WhoopProtocol // independent reference the database carries. // // USAGE: sleepbench --db /path/to/whoop.sqlite [--device my-whoop-noop] [--pad 3600] [--csv out.csv] +// [--exclude ] // // The database is never written to and never lives in this repository. Pass a COPY; never point this at // a device's live database. @@ -20,6 +21,15 @@ struct Args { var streamDevice = "my-whoop" var pad = 3600 var csv: String? + /// Session `startTs` values held out of the stage-locked calibration set in section E. + /// + /// Section E's reference set has to stay FIXED across the two builds of a before/after comparison, or + /// the denominator moves underneath the delta. That is why it keys on `stagelock` rather than on + /// section B's recipe-dependent exclusion. But `stagelock` alone does not guarantee the row is + /// independent of V2 (see the self-comparison warning below), and the honest fix for a contaminated row + /// is to name it ONCE and pass the same `--exclude` list to every build being compared. Excluding by + /// explicit timestamp keeps the set frozen; excluding by "matches this build's V2" would not. + var exclude: Set = [] } var a = Args() var it = CommandLine.arguments.dropFirst().makeIterator() @@ -30,11 +40,18 @@ while let k = it.next() { case "--stream-device": a.streamDevice = it.next() ?? a.streamDevice case "--pad": a.pad = Int(it.next() ?? "") ?? a.pad case "--csv": a.csv = it.next() + case "--exclude": a.exclude = Set((it.next() ?? "").split(separator: ",").compactMap { Int($0) }) default: FileHandle.standardError.write(Data("unknown argument \(k)\n".utf8)); exit(2) } } guard !a.db.isEmpty else { - print("usage: sleepbench --db [--device ] [--pad ] [--csv ]") + print(""" + usage: sleepbench --db [--device ] [--stream-device ] [--pad ] [--csv ] + [--exclude ] + + --exclude session startTs values held out of section E's calibration reference set. Pass the SAME + list to every build in a before/after comparison; section E names the candidates. + """) exit(2) } @@ -43,7 +60,7 @@ let rows = try db.sessions(device: a.device) guard !rows.isEmpty else { print("no sessions for device \(a.device)"); exit(1) } // Which `userEdited` rows carry a HUMAN-AUTHORED hypnogram rather than machine stages re-derived over a // human-corrected window. Only the former can serve as a stage-level reference. -let stageLocked = try db.stageLockedStarts(device: a.device) +let stageLocked = try db.stageLockedStarts(device: a.device).subtracting(a.exclude) let editedCount = rows.filter { $0.userEdited }.count print(""" @@ -149,12 +166,16 @@ var v1Match = 0, v2Match = 0, neither = 0 /// This says nothing about HOW the row came to match — the harness reports the fact and excludes the /// night, and deliberately does not assert a mechanism it has not established. var indistinguishableFromV2: Set = [] +/// Per-night epoch agreement between the STORED hypnogram and a fresh V2 replay, by `startTs`. Section E +/// reads it to warn when its own reference set contains rows V2 cannot be scored against independently. +var v2AgreePct: [Int: Double] = [:] print("night edited epochs V1 agree V2 agree verdict") for nt in nights { let stored = epochLabels(nt.row.stages, start: nt.row.startTs, end: nt.row.endTs) let a1 = zip(stored, nt.v1).filter { $0 == $1 }.count let a2 = zip(stored, nt.v2).filter { $0 == $1 }.count let p1 = Double(a1) / Double(stored.count) * 100, p2 = Double(a2) / Double(stored.count) * 100 + v2AgreePct[nt.row.startTs] = p2 let verdict = p2 >= 99.9 ? "V2 exact" : (p1 >= 99.9 ? "V1 exact" : (p2 > p1 ? "V2 closer" : "V1 closer")) if !nt.row.userEdited { if p2 >= 99.9 { v2Match += 1 } else if p1 >= 99.9 { v1Match += 1 } else { neither += 1 } @@ -347,12 +368,187 @@ for (nm, grp) in [("mean HR >= \(Int(hiCut))", hi), ("mean HR < \(Int(hiCut))", print(" \(nm.padding(toLength: 18, withPad: " ", startingAt: 0)) n=\(e.count) mean machine wake error \(f(mean(e), 7, 1)) min median \(f(median(e), 7, 1))") } -// MARK: - E. Machine cohort split — did anything change mid-history? +// MARK: - E. Per-stage fraction calibration — the guard kappa does not provide + +/// The stage-level reference set for calibration: rows carrying a `stagelock` cursor, i.e. rows whose +/// stages a human authored directly. +/// +/// This is deliberately NOT section B's set. B excludes an edited night whose stored hypnogram is a +/// byte-exact replay of the CURRENT V2, which is the right call for B's question but is version-dependent +/// by construction: change the recipe and a night can enter or leave the exclusion, silently swapping the +/// denominator underneath a before/after comparison. The `stagelock` set comes from `cursors` and does not +/// move when the recipe moves, which is what makes a calibration delta between two builds legitimate. +/// +/// The cost of that choice, and why the warning below exists: a `stagelock` cursor proves the stages +/// arrived through the `edit_sleep_stages` path, NOT that they differ from what V2 would emit. A row whose +/// stored hypnogram replays byte-exact from V2 carries no information about V2 — it hands the incumbent a +/// free perfect night and charges every alternative for the same nights. Section B drops such rows; section +/// E cannot drop them the same way without reintroducing the version-dependence it exists to avoid. So the +/// harness NAMES them and leaves the choice to the operator, who freezes it with `--exclude`. +let calRef = nights.filter { $0.truth != nil && stageLocked.contains($0.row.startTs) } +/// The cohort PR #348 broke: a plausible full sleep opportunity. #437's field symptom was a HEALTHY night +/// re-scored 6% -> 23% awake, and an aggregate that pools short and fragmented nights hides exactly that. +let healthyMinutes = 300.0 +func isHealthy(_ nt: Night) -> Bool { nt.row.durMin >= healthyMinutes } +let calRefHealthy = calRef.filter(isHealthy) +let healthyNights = nights.filter(isHealthy) + +/// The recipes scored below, in the order the rest of the harness uses them. +let variants: [(String, (Night) -> [String])] = [ + ("V1", { $0.v1 }), ("V2", { $0.v2 }), ("V1+veto", { $0.v1Veto }), ("V2+veto", { $0.v2Veto }), +] + +print(""" + +================================================================================ +E. PER-STAGE FRACTION CALIBRATION (reference (a), stage-locked, n = \(calRef.count) nights) +================================================================================ +Cohen's kappa scores whether the label at an epoch matches; it does not constrain how much of the night a +recipe spends at each stage. PR #348 fitted the stager to DREAMT, raised kappa on all three benchmarks with +a held-out gap of -0.027, and was reverted 48 h later by PR #437 for re-scoring a healthy night from 6% to +23% awake. The revert's own words: "kappa doesn't guard stage-fraction calibration." This section is that +guard. Bias is signed, in percentage points of the night: POSITIVE = the recipe spends more of the night at +that stage than the human did. +""") + +// SELF-COMPARISON AUDIT of this section's own reference set. A stage-locked row whose stored hypnogram is a +// byte-exact V2 replay scores V2 against itself: it is a guaranteed 100% for the incumbent and a guaranteed +// loss for any alternative, so it biases a before/after delta toward "change nothing" by an amount nobody +// reading only the kappa line would see. The audit reports the contamination rather than silently absorbing +// it, and names the timestamps so the same held-out set can be pinned across both builds with `--exclude`. +let selfMatched = calRef.filter { (v2AgreePct[$0.row.startTs] ?? 0) >= 99.9 } +let nearMatched = calRef.filter { + let p = v2AgreePct[$0.row.startTs] ?? 0 + return p >= 95.0 && p < 99.9 +} +print(" E.-1 SELF-COMPARISON AUDIT of the reference set") +if selfMatched.isEmpty && nearMatched.isEmpty { + print(" clean — no stage-locked night replays from V2 at >= 95% agreement.") +} else { + print(" night stored-vs-V2 agreement verdict") + for nt in (selfMatched + nearMatched).sorted(by: { $0.row.startTs < $1.row.startTs }) { + let p = v2AgreePct[nt.row.startTs] ?? 0 + print(" \(nt.night) (startTs \(nt.row.startTs)) \(f(p, 12, 2))% \(p >= 99.9 ? "BYTE-EXACT — carries no information about V2" : "near-exact — nearly none")") + } + let pct = Double(selfMatched.count) / Double(max(1, calRef.count)) * 100 + print(""" + \(selfMatched.count) of \(calRef.count) reference nights (\(f(pct, 0, 0))%) are byte-exact V2 replays; \(nearMatched.count) more sit at >= 95%. + These inflate every V2 figure in E.0-E.4 and penalise any alternative recipe by construction. To + compare two builds honestly, pass the SAME list to BOTH: + --exclude \(selfMatched.map { String($0.row.startTs) }.joined(separator: ",")) + The harness asserts no mechanism for the match; it reports only that the row cannot serve as an + independent reference for the recipe it replays. + """) +} +if !a.exclude.isEmpty { + print(" HELD OUT by --exclude this run: \(a.exclude.sorted().map(String.init).joined(separator: ", "))") +} + +// Agreement on the SAME nights the calibration below scores, so the two can be read against each other. +// This is the contrast the section exists for: kappa here can rise while the biases below get worse. +print(" E.0 AGREEMENT ON THIS REFERENCE SET — the number that is not sufficient on its own") +for (name, get) in variants { + var c4 = Confusion(classes: stageOrder), c2 = Confusion(classes: ["wake", "sleep"]) + for nt in calRef { + guard let t = nt.truth else { continue } + c4.merge(confusion(ref: t, pred: get(nt))) + c2.merge(confusion(ref: toSleepWake(t), pred: toSleepWake(get(nt)), classes: ["wake", "sleep"])) + } + print(" \(name.padding(toLength: 9, withPad: " ", startingAt: 0)) 4-class acc \(f(c4.accuracy * 100, 5, 1))% kappa \(f(c4.kappa, 6, 3)) | sleep/wake acc \(f(c2.accuracy * 100, 5, 1))% kappa \(f(c2.kappa, 6, 3)) (epochs \(c4.total))") +} + +func calibrationTable(_ title: String, _ set: [Night]) { + print("\n \(title)") + guard !set.isEmpty else { print(" no reference nights in this stratum"); return } + for (name, get) in variants { + print(" \(name) — stage pred mean% ref mean% bias pp MAE pp worst night") + for s in stageOrder { + var pred: [Double] = [], ref: [Double] = [], err: [Double] = [] + var worst = ("", 0.0) + for nt in set { + guard let t = nt.truth else { continue } + let p = stagePercentages(get(nt))[s]!, q = stagePercentages(t)[s]! + pred.append(p); ref.append(q); err.append(p - q) + if abs(p - q) > abs(worst.1) { worst = (nt.night, p - q) } + } + print(" \(s.padding(toLength: 6, withPad: " ", startingAt: 0)) \(f(mean(pred), 12, 2)) \(f(mean(ref), 11, 2)) \(f(mean(err), 10, 2)) \(f(mae(err), 9, 2)) \(worst.0) \(f(worst.1, 7, 1))") + } + } +} +calibrationTable("E.1 AGGREGATE — every stage-locked reference night (n = \(calRef.count))", calRef) +calibrationTable(""" +E.2 HEALTHY STRATUM — in-bed >= \(Int(healthyMinutes / 60)) h (n = \(calRefHealthy.count) of \(calRef.count)) + Reported separately because the aggregate hides it: a recipe can look calibrated overall while + blowing out on exactly the nights a healthy sleeper has. +""", calRefHealthy) + +// A reference-free view of the same quantity. It needs no human labels, so it runs over every replayed +// night — a far larger n than the reference set, and enough to catch a distribution shift on its own. +print(""" + + E.3 SHIPPED POPULATION — V2 stage fractions over all replayed nights, no reference needed + stage ALL n = \(nights.count) HEALTHY >= \(Int(healthyMinutes / 60)) h n = \(healthyNights.count) + mean% median% p10 p90 mean% median% p10 p90 +""") +for s in stageOrder { + let all = nights.map { stagePercentages($0.v2)[s]! } + let hea = healthyNights.map { stagePercentages($0.v2)[s]! } + print(" \(s.padding(toLength: 8, withPad: " ", startingAt: 0)) \(f(mean(all), 8, 2)) \(f(median(all), 9, 2)) \(f(percentile(all, 0.10), 8, 2)) \(f(percentile(all, 0.90), 8, 2)) \(f(mean(hea), 8, 2)) \(f(median(hea), 9, 2)) \(f(percentile(hea, 0.10), 8, 2)) \(f(percentile(hea, 0.90), 8, 2))") +} + +print("\n E.4 PER NIGHT — V2 signed bias in pp, on the \(calRef.count) stage-locked reference nights") +print(" night inbed healthy wake light deep rem") +for nt in calRef { + // `calRef` filters on `truth != nil`; the guard keeps that guarantee local, because an empty reference + // would silently turn every predicted fraction into a bias against zero. + guard let t = nt.truth else { continue } + let b = stageBias(ref: t, pred: nt.v2) + print(" \(nt.night) \(f(nt.row.durMin, 6, 0)) \(isHealthy(nt) ? " Y" : " -") \(f(b["wake"]!, 7, 1)) \(f(b["light"]!, 7, 1)) \(f(b["deep"]!, 7, 1)) \(f(b["rem"]!, 7, 1))") +} + +// MARK: - F. First-REM latency — a stage-timing property kappa is also blind to + +print(""" + +================================================================================ +F. FIRST-REM LATENCY — minutes from staged sleep onset to the first REM epoch +================================================================================ +Fraction calibration pins how MUCH REM a recipe emits; it says nothing about WHEN. A recipe can hit the +right REM total at the wrong time — REM in the first minutes after onset is physiologically implausible in +a healthy adult, and a stager that emits it is mislabelling early light sleep whatever its kappa says. +Nights that never reach REM are counted separately rather than folded in as a zero. +""") +print(" set n median mean p10 p90 min max no-REM") +func latencyRow(_ name: String, _ set: [Night], _ get: (Night) -> [String]) { + let v = set.compactMap { firstRemLatencyMinutes(get($0)) } + guard !v.isEmpty else { print(" \(name.padding(toLength: 22, withPad: " ", startingAt: 0)) no nights reach REM"); return } + print(" \(name.padding(toLength: 22, withPad: " ", startingAt: 0)) \(String(format: "%4d", v.count)) \(f(median(v), 8, 1)) \(f(mean(v), 8, 1)) \(f(percentile(v, 0.10), 8, 1)) \(f(percentile(v, 0.90), 8, 1)) \(f(v.min()!, 8, 1)) \(f(v.max()!, 8, 1)) \(String(format: "%8d", set.count - v.count))") +} +for (name, get) in variants { latencyRow("\(name), all nights", nights, get) } +latencyRow("V2, healthy >= \(Int(healthyMinutes / 60)) h", healthyNights, { $0.v2 }) +latencyRow("V2, reference nights", calRef, { $0.v2 }) +// The human's own latency distribution on the same nights — the target the replay rows above are aiming at. +latencyRow("HUMAN reference", calRef, { $0.truth ?? [] }) + +print("\n PER NIGHT — first-REM latency in minutes, on the \(calRef.count) stage-locked reference nights") +print(" night inbed healthy human V2 error") +var latErr: [Double] = [] +for nt in calRef { + guard let t = nt.truth else { continue } + // n/a in either column means that night never reached REM, which is not an error of zero minutes. + let h = firstRemLatencyMinutes(t), p = firstRemLatencyMinutes(nt.v2) + var err = Double.nan + if let h, let p { err = p - h; latErr.append(err) } + print(" \(nt.night) \(f(nt.row.durMin, 6, 0)) \(isHealthy(nt) ? " Y" : " -") \(f(h ?? .nan, 7, 1)) \(f(p ?? .nan, 7, 1)) \(f(err, 7, 1))") +} +print(" signed error (V2 - human), nights where BOTH reach REM: n=\(latErr.count) median \(f(median(latErr), 6, 1)) mean \(f(mean(latErr), 6, 1)) MAE \(f(mae(latErr), 6, 1))") + +// MARK: - G. Machine cohort split — did anything change mid-history? print(""" ================================================================================ -E. PER-NIGHT DETAIL (all \(nights.count) sessions) +G. PER-NIGHT DETAIL (all \(nights.count) sessions) ================================================================================ """) print("night edited band inbed meanHR storedW V1 W V2 W storedEff") @@ -362,12 +558,19 @@ for nt in nights { } if let path = a.csv { - var out = "night,userEdited,hasBand,inbedMin,meanHR,storedWake,truthWake,machineWake,v1Wake,v2Wake,v1VetoWake,v2VetoWake,storedEff\n" + /// Blank rather than a sentinel wherever a night has no value: a night that never reaches REM has no + /// first-REM latency, and writing 0 there would be a different claim. + func num(_ v: Double?, _ p: Int = 1) -> String { v.map { String(format: "%.\(p)f", $0) } ?? "" } + var out = "night,userEdited,stageLocked,healthy,hasBand,inbedMin,meanHR,storedWake,truthWake,machineWake,v1Wake,v2Wake,v1VetoWake,v2VetoWake,storedEff" + for s in stageOrder { out += ",v2\(s.capitalized)Pct,truth\(s.capitalized)Pct" } + out += ",v2FirstRemLatMin,truthFirstRemLatMin\n" for nt in nights { let stored = epochLabels(nt.row.stages, start: nt.row.startTs, end: nt.row.endTs) - let tw = nt.truth.map { String(format: "%.1f", wakeMinutes($0)) } ?? "" - let mw = nt.machineWakeMin.map { String(format: "%.1f", $0) } ?? "" - out += "\(nt.night),\(nt.row.userEdited),\(!nt.band.isEmpty),\(String(format: "%.1f", nt.row.durMin)),\(String(format: "%.1f", nt.meanHR)),\(String(format: "%.1f", wakeMinutes(stored))),\(tw),\(mw),\(String(format: "%.1f", wakeMinutes(nt.v1))),\(String(format: "%.1f", wakeMinutes(nt.v2))),\(String(format: "%.1f", wakeMinutes(nt.v1Veto))),\(String(format: "%.1f", wakeMinutes(nt.v2Veto))),\(nt.row.efficiency.map { String(format: "%.4f", $0) } ?? "")\n" + let locked = stageLocked.contains(nt.row.startTs) + out += "\(nt.night),\(nt.row.userEdited),\(locked),\(isHealthy(nt)),\(!nt.band.isEmpty),\(num(nt.row.durMin)),\(num(nt.meanHR)),\(num(wakeMinutes(stored))),\(num(nt.truth.map { wakeMinutes($0) })),\(num(nt.machineWakeMin)),\(num(wakeMinutes(nt.v1))),\(num(wakeMinutes(nt.v2))),\(num(wakeMinutes(nt.v1Veto))),\(num(wakeMinutes(nt.v2Veto))),\(num(nt.row.efficiency, 4))" + let pv = stagePercentages(nt.v2), pt = nt.truth.map { stagePercentages($0) } + for s in stageOrder { out += ",\(num(pv[s], 2)),\(num(pt?[s], 2))" } + out += ",\(num(firstRemLatencyMinutes(nt.v2))),\(num(nt.truth.flatMap { firstRemLatencyMinutes($0) }))\n" } try out.write(toFile: path, atomically: true, encoding: .utf8) FileHandle.standardError.write(Data("wrote \(path)\n".utf8)) diff --git a/Tools/SleepBench/Tests/sleepbenchTests/CalibrationMetricsTests.swift b/Tools/SleepBench/Tests/sleepbenchTests/CalibrationMetricsTests.swift new file mode 100644 index 0000000000..a68907baab --- /dev/null +++ b/Tools/SleepBench/Tests/sleepbenchTests/CalibrationMetricsTests.swift @@ -0,0 +1,132 @@ +import StrandAnalytics +import XCTest + +@testable import sleepbench + +/// The stage-fraction and first-REM-latency primitives are pure functions over 30 s label arrays, so they +/// are pinned here on hand-built hypnograms — no database, no health data, no `--db` argument. +final class CalibrationMetricsTests: XCTestCase { + + /// A label array of `n` epochs all at `stage`. + private func run(_ stage: String, _ n: Int) -> [String] { [String](repeating: stage, count: n) } + + // MARK: - stagePercentages + + func testStagePercentagesSumTo100AndMatchCounts() { + let labels = run("wake", 10) + run("light", 50) + run("deep", 20) + run("rem", 20) + let p = stagePercentages(labels) + XCTAssertEqual(p["wake"]!, 10, accuracy: 1e-9) + XCTAssertEqual(p["light"]!, 50, accuracy: 1e-9) + XCTAssertEqual(p["deep"]!, 20, accuracy: 1e-9) + XCTAssertEqual(p["rem"]!, 20, accuracy: 1e-9) + XCTAssertEqual(stageOrder.reduce(0) { $0 + p[$1]! }, 100, accuracy: 1e-9) + } + + /// A stage a recipe stops emitting must read as 0%, not as a missing key — otherwise the bias against a + /// reference that DOES contain that stage would silently vanish instead of showing up as the regression + /// it is. This is the shape of the #348 failure: a recipe reallocating a whole stage. + func testStagePercentagesReportsUnusedStagesAsZero() { + let p = stagePercentages(run("wake", 4)) + XCTAssertEqual(Set(p.keys), Set(stageOrder)) + XCTAssertEqual(p["wake"]!, 100, accuracy: 1e-9) + for s in ["light", "deep", "rem"] { XCTAssertEqual(p[s]!, 0, accuracy: 1e-9) } + } + + func testStagePercentagesOnEmptyLabelsIsAllZeroNotNaN() { + let p = stagePercentages([]) + for s in stageOrder { XCTAssertEqual(p[s]!, 0, accuracy: 1e-9) } + } + + // MARK: - stageBias + + /// Sign convention is load-bearing: the whole point of the metric is telling an over-call from an + /// under-call, so POSITIVE must mean the recipe spends more of the night at that stage than the human. + func testStageBiasIsPositiveWhenTheRecipeOverCalls() { + let ref = run("wake", 10) + run("light", 90) // 10% wake + let pred = run("wake", 30) + run("light", 70) // 30% wake + let b = stageBias(ref: ref, pred: pred) + XCTAssertEqual(b["wake"]!, 20, accuracy: 1e-9) + XCTAssertEqual(b["light"]!, -20, accuracy: 1e-9) + XCTAssertEqual(b["deep"]!, 0, accuracy: 1e-9) + } + + /// The #437 field symptom, reproduced numerically: a healthy night re-scored 6% -> 23% awake is a + /// +17 pp wake bias, which is what this metric exists to surface. + func testStageBiasReproducesTheRevertScaleSymptom() { + let ref = run("wake", 6) + run("light", 94) + let pred = run("wake", 23) + run("light", 77) + XCTAssertEqual(stageBias(ref: ref, pred: pred)["wake"]!, 17, accuracy: 1e-9) + } + + /// Two nights that miss in opposite directions cancel in the signed mean but not in the MAE, which is + /// why both are reported: a mean bias near zero is not by itself evidence of good calibration. + func testMeanBiasCancelsWhereMAEDoesNot() { + let errs = [12.0, -12.0] + XCTAssertEqual(mean(errs), 0, accuracy: 1e-9) + XCTAssertEqual(mae(errs), 12, accuracy: 1e-9) + } + + // MARK: - firstRemLatencyMinutes + + /// Latency is measured from staged ONSET, not from the start of the session — a night with a long + /// awake-in-bed lead-in must not be credited with a long REM latency. + func testFirstRemLatencyIsMeasuredFromOnsetNotSessionStart() { + // 60 epochs awake in bed, then 20 light, then REM. Onset is epoch 60, REM at epoch 80. + let labels = run("wake", 60) + run("light", 20) + run("rem", 10) + XCTAssertEqual(firstRemLatencyMinutes(labels)!, 10.0, accuracy: 1e-9) // 20 epochs x 30 s + } + + func testFirstRemLatencyUsesTheThirtySecondGrid() { + let labels = run("light", 2) + run("rem", 1) + XCTAssertEqual(firstRemLatencyMinutes(labels)!, 1.0, accuracy: 1e-9) + XCTAssertEqual(epochSeconds, 30.0, accuracy: 1e-9) + } + + func testFirstRemLatencyIsZeroWhenREMIsTheFirstSleepEpoch() { + XCTAssertEqual(firstRemLatencyMinutes(run("wake", 5) + run("rem", 5))!, 0.0, accuracy: 1e-9) + } + + /// A night with no REM is a DIFFERENT outcome from a night whose REM lands at minute zero. Collapsing + /// them would let a recipe that stops emitting REM read as one with an excellent short latency. + func testFirstRemLatencyIsNilWithoutREMAndWithoutSleep() { + XCTAssertNil(firstRemLatencyMinutes(run("wake", 20) + run("light", 40))) + XCTAssertNil(firstRemLatencyMinutes(run("wake", 100))) + XCTAssertNil(firstRemLatencyMinutes([])) + } + + // MARK: - percentile + + func testPercentileSpansMinToMax() { + let v = [1.0, 2, 3, 4, 5, 6, 7, 8, 9, 10] + XCTAssertEqual(percentile(v, 0.0), 1, accuracy: 1e-9) + XCTAssertEqual(percentile(v, 1.0), 10, accuracy: 1e-9) + XCTAssertEqual(percentile(v, 0.5), 6, accuracy: 1e-9) + } + + func testPercentileIsOrderIndependentAndSafeOnEmpty() { + XCTAssertEqual(percentile([5.0, 1, 3], 0.0), 1, accuracy: 1e-9) + XCTAssertTrue(percentile([], 0.5).isNaN) + } + + // MARK: - the grid the metrics run on + + /// `stagePercentages` and `firstRemLatencyMinutes` both consume `epochLabels` output in the harness, so + /// the expansion from a `StageSegment` tiling has to land the stages where the metrics expect them. + func testMetricsAgreeWithEpochLabelExpansion() { + // 30 min wake, 30 min light, 30 min REM, on a 90 min session. + let start = 1_700_000_000 + let segs = [ + StageSegment(start: start, end: start + 1800, stage: "wake"), + StageSegment(start: start + 1800, end: start + 3600, stage: "light"), + StageSegment(start: start + 3600, end: start + 5400, stage: "rem"), + ] + let labels = epochLabels(segs, start: start, end: start + 5400) + XCTAssertEqual(labels.count, 180) + let p = stagePercentages(labels) + XCTAssertEqual(p["wake"]!, 100.0 / 3, accuracy: 1e-6) + XCTAssertEqual(p["light"]!, 100.0 / 3, accuracy: 1e-6) + XCTAssertEqual(p["rem"]!, 100.0 / 3, accuracy: 1e-6) + XCTAssertEqual(firstRemLatencyMinutes(labels)!, 30.0, accuracy: 1e-9) + XCTAssertEqual(wakeMinutes(labels), 30.0, accuracy: 1e-9) + } +}