diff --git a/Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift b/Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift index 0f8b043496..441623be45 100644 --- a/Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift +++ b/Packages/StrandAnalytics/Sources/StrandAnalytics/SleepStagerV2.swift @@ -323,6 +323,8 @@ public enum SleepStagerV2 { } var raws: [Raw] = [] var allJerks: [Double] = [] + // The RSA transform's twiddle factors, per grid length, for this night only (see `RespDFT`). + var respDFT: [Int: RespDFT] = [:] let firstE = ((start + 29) / 30) * 30 var e = firstE while e < end { @@ -354,7 +356,7 @@ public enum SleepStagerV2 { if let vs = rrBy[s] { for v in vs { beats.append((Double(s), min(max(v, 300), 2000))) } } } beats.sort { $0.0 != $1.0 ? $0.0 < $1.0 : $0.1 < $1.1 } - let respReg = respRegularity(beats) + let respReg = respRegularity(beats, dft: &respDFT) raws.append(Raw(start: e, hr: hrMean, hrVar: hrVar, hrFlat11: hrFlat11, jerks: jerks, gapSec: max(1, gseq.count - 1), jerkMax: jerkMax, @@ -387,10 +389,43 @@ public enum SleepStagerV2 { return feats } + /// The band-limited DFT's twiddle factors for one 4 Hz grid length `n`: `cos`/`sin` of + /// `-2π·k/n · j` for every in-band bin `k` and sample `j`. + /// + /// `respRegularity` runs once per 30-second epoch and used to evaluate these ~53 × ~836 pairs every time, + /// which was 13% of a cold re-score's CPU (Time Profiler, `__sincos_stret`). They depend on `n` alone, + /// and a night's beat window is whole seconds wide, so it sees only a handful of distinct `n`. The table + /// computes each factor with the very expressions the transform used inline, so reading it back is + /// bit-identical to recomputing it; a night builds it once per `n` and drops it with the night. + /// Kotlin twin: `SleepStagerV2.RespDft`. + struct RespDFT { + let kLo: Int + let cosines: [[Double]] + let sines: [[Double]] + + init(n: Int, kLo: Int, kHi: Int) { + self.kLo = kLo + var cosines: [[Double]] = [] + var sines: [[Double]] = [] + for k in kLo...kHi { + let w = -2.0 * Double.pi * Double(k) / Double(n) + var c = [Double](repeating: 0, count: n) + var s = [Double](repeating: 0, count: n) + for j in 0.. Double? { + /// `dft` holds the twiddle factors per grid length (`RespDFT`), filled on first use; the caller keeps + /// one for a night. Kotlin twin: `SleepStagerV2.respRegularity`. + static func respRegularity(_ beats: [(Double, Double)], dft: inout [Int: RespDFT]) -> Double? { if beats.count < 12 { return nil } let t0 = beats.first!.0, tN = beats.last!.0 if tN <= t0 { return nil } @@ -414,11 +449,18 @@ public enum SleepStagerV2 { let kLo = Int(ceil(0.15 * 0.25 * Double(n))) let kHi = Int(floor(0.40 * 0.25 * Double(n))) if kHi < kLo || kLo < 0 { return nil } + let table: RespDFT + if let cached = dft[n] { + table = cached + } else { + table = RespDFT(n: n, kLo: kLo, kHi: kHi) + dft[n] = table + } var maxP = 0.0, sumP = 0.0 for k in kLo...kHi { var re = 0.0, im = 0.0 - let w = -2.0 * Double.pi * Double(k) / Double(n) - for j in 0.. maxP { maxP = p } diff --git a/Packages/StrandAnalytics/Tests/StrandAnalyticsTests/SleepStagerRespDFTTableTests.swift b/Packages/StrandAnalytics/Tests/StrandAnalyticsTests/SleepStagerRespDFTTableTests.swift new file mode 100644 index 0000000000..2a2219deb4 --- /dev/null +++ b/Packages/StrandAnalytics/Tests/StrandAnalyticsTests/SleepStagerRespDFTTableTests.swift @@ -0,0 +1,105 @@ +import XCTest +@testable import StrandAnalytics + +/// `SleepStagerV2.respRegularity` reads its twiddle factors from a per-night table (`RespDFT`) instead of +/// evaluating `cos`/`sin` for every epoch. The stage labels hang on this value, so the table must be a pure +/// speed change: the reference below is the function as it was, and the two must return the SAME Double — +/// not a close one — for every window, including the nil cases, with one table reused across windows the way +/// a night reuses it. Mirrors the Kotlin `SleepStagerRespDftTableTest`. +final class SleepStagerRespDFTTableTests: XCTestCase { + + /// The function before the table, verbatim. + private func respRegularityBeforeTable(_ 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 + } + + /// A beat window as `features()` builds one: whole-second times over `[e-90, e+120)`, several beats a second + /// or none, values clamped to 300…2000 ms, sorted by time then value. + private func windowForTest(_ rng: inout SplitMixForDFTTest, span: Int, density: Int) -> [(Double, Double)] { + var beats: [(Double, Double)] = [] + for s in 0.. 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) + } +} diff --git a/Tools/parity_twin_map.json b/Tools/parity_twin_map.json index 1c53c8c166..42ba0a5d73 100644 --- a/Tools/parity_twin_map.json +++ b/Tools/parity_twin_map.json @@ -19,15 +19,15 @@ }, "authority": { "files": {"count": 505, "sha256": "581ec0dcab0d8486ce113646cb5cbad71769353eb3aef75f45486e03591733a0"}, - "functions": {"count": 4509, "sha256": "828aa64d6abfd21f2dbb94bb552099409555b564dad5d04dbab88d8555466486"}, + "functions": {"count": 4509, "sha256": "ea93147013013ce1c99232e5a77093b38d321f231073456cb25c9cfe393de567"}, "properties": {"count": 465, "sha256": "64fcd18c85bf33a4738d154d45ec86967fdc974e78a448d16e9504d9b13609eb"}, "constants": {"count": 1969, "sha256": "0cfd548d65ff82bdaa5860c8c63a0f8ac5bacbabc9a5abeb83773c1522147652"}, - "file_pairs": {"count": 75, "sha256": "d581f60241ec523c595a7e5ff123f7987caee501d8be4ebf6df5ba105e26fe1f"}, - "function_pairs": {"count": 206, "sha256": "5b2e4542dfea5e6e75ae3132e018c66bc12fd77092d0e2c996699c32df20ce0b"}, + "file_pairs": {"count": 76, "sha256": "453da3e71534978be500547a44b9add326631fa1920c8871b83e11457f336ade"}, + "function_pairs": {"count": 207, "sha256": "8fb4bad334ea308f0a6bf7e9b94b6001afbe7e2e7cc08b7bce7cf6edaf3c5cb0"}, "property_pairs": {"count": 150, "sha256": "be7168d309c8a9cf5f3f25438c028cc60ab64202bb6e9c09a6c553e3699ad0e2"}, "constant_pairs": {"count": 685, "sha256": "a949b8f3ffc009b8054d7bc56f4131f7baf134b6f09823f323853d5cfb52c246"}, - "unpaired_files": {"count": 377, "sha256": "14e70f77d7e11b4eae3b8412cc493c1aa1bf24dfe8fbc5eca71f81b5daf5253f"}, - "unpaired_functions": {"count": 4110, "sha256": "16b46b58f2be94f5653bc25ea4099ed8d220a418b0cf193e0d408722c4910708"}, + "unpaired_files": {"count": 375, "sha256": "ba14f7aa2a3967ebc2aa21311917ae3c0c3db218186fac0a6d297f8969580f1a"}, + "unpaired_functions": {"count": 4108, "sha256": "7e83a1360a1a8dfe35ee8fefa988bf3ad5a8def29b94782dbff31d8a7d24cdc9"}, "unpaired_properties": {"count": 165, "sha256": "fe609dd25384cda6de1616d36753fed68431b1a1f9a73b6fd664c7ff35089955"}, "unpaired_constants": {"count": 599, "sha256": "b7467c65623ecc296604558a639d62eee2e00874c56d21acc7ec6d3b6264fdc7"} } diff --git a/android/app/src/main/java/com/noop/analytics/SleepStagerV2.kt b/android/app/src/main/java/com/noop/analytics/SleepStagerV2.kt index 918a6607e8..3cdf29a4d7 100644 --- a/android/app/src/main/java/com/noop/analytics/SleepStagerV2.kt +++ b/android/app/src/main/java/com/noop/analytics/SleepStagerV2.kt @@ -340,6 +340,8 @@ object SleepStagerV2 { ) val raws = ArrayList() val allJerks = ArrayList() + // The RSA transform's twiddle factors, per grid length, for this night only (see [RespDft]). + val respDft = HashMap() val firstE = ((start + 29) / 30) * 30 var e = firstE while (e < end) { @@ -373,7 +375,7 @@ object SleepStagerV2 { bs++ } beats.sortWith(compareBy({ it.first }, { it.second })) - val respReg = respRegularity(beats) + val respReg = respRegularity(beats, respDft) raws.add(Raw( start = e, hr = hrMean, hrVar = hrVar, hrFlat11 = hrFlat11, @@ -405,12 +407,35 @@ object SleepStagerV2 { return feats } + /** + * The band-limited DFT's twiddle factors for one 4 Hz grid length `n`: `cos`/`sin` of `-2π·k/n · j` for + * every in-band bin `k` and sample `j`. [respRegularity] runs once per 30-second epoch and used to evaluate + * these ~53 × ~836 pairs every time; they depend on `n` alone, and a night's beat window is whole seconds + * wide, so it sees only a handful of distinct `n`. Each factor is computed with the very expressions the + * transform used inline, so reading it back is bit-identical to recomputing it. Swift twin: + * `SleepStagerV2.RespDFT`. + */ + internal class RespDft(n: Int, val kLo: Int, kHi: Int) { + val cosines: Array = Array(kHi - kLo + 1) { DoubleArray(n) } + val sines: Array = Array(kHi - kLo + 1) { DoubleArray(n) } + + init { + for (k in kLo..kHi) { + val w = -2.0 * PI * k / n + val c = cosines[k - kLo]; val s = sines[k - kLo] + for (j in 0 until n) { val a = w * j; c[j] = cos(a); s[j] = sin(a) } + } + } + } + /** * RSA respiration regularity: tachogram → 4 Hz resample → detrend → power spectrum → peak/sum of the * 0.15–0.40 Hz (9–24 brpm) band. Returns spectral peakedness (higher = more regular breathing) or null * when there are too few beats. A direct band-limited DFT (only the ~50 in-band bins are needed). + * [dft] holds the twiddle factors per grid length ([RespDft]), filled on first use; the caller keeps one + * for a night. Swift twin: `SleepStagerV2.respRegularity`. */ - private fun respRegularity(beats: List>): Double? { + internal fun respRegularity(beats: List>, dft: MutableMap): Double? { if (beats.size < 12) return null val t0 = beats.first().first; val tN = beats.last().first if (tN <= t0) return null @@ -434,11 +459,12 @@ object SleepStagerV2 { val kLo = ceil(0.15 * 0.25 * n).toInt() val kHi = floor(0.40 * 0.25 * n).toInt() if (kHi < kLo || kLo < 0) return null + val table = dft.getOrPut(n) { RespDft(n, kLo, kHi) } var maxP = 0.0; var sumP = 0.0 for (k in kLo..kHi) { var re = 0.0; var im = 0.0 - val w = -2.0 * PI * k / n - for (j in 0 until n) { val a = w * j; re += y[j] * cos(a); im += y[j] * sin(a) } + val c = table.cosines[k - kLo]; val s = table.sines[k - kLo] + for (j in 0 until n) { re += y[j] * c[j]; im += y[j] * s[j] } val p = re * re + im * im sumP += p if (p > maxP) maxP = p diff --git a/android/app/src/test/java/com/noop/analytics/SleepStagerRespDftTableTest.kt b/android/app/src/test/java/com/noop/analytics/SleepStagerRespDftTableTest.kt new file mode 100644 index 0000000000..8e2b88fc4a --- /dev/null +++ b/android/app/src/test/java/com/noop/analytics/SleepStagerRespDftTableTest.kt @@ -0,0 +1,98 @@ +package com.noop.analytics + +import kotlin.math.PI +import kotlin.math.ceil +import kotlin.math.cos +import kotlin.math.floor +import kotlin.math.sin +import org.junit.Assert.assertEquals +import org.junit.Assert.assertTrue +import org.junit.Test + +/** + * [SleepStagerV2.respRegularity] reads its twiddle factors from a per-night table ([SleepStagerV2.RespDft]) + * instead of evaluating `cos`/`sin` for every epoch. The stage labels hang on this value, so the table must be a + * pure speed change: the reference below is the function as it was, and the two must return the SAME Double, + * bit for bit, for every window, including the null cases, with one table reused across windows the way a + * night reuses it. Mirrors the Swift `SleepStagerRespDFTTableTests`. + */ +class SleepStagerRespDftTableTest { + + /** The function before the table, verbatim. */ + private fun respRegularityBeforeTable(beats: List>): Double? { + if (beats.size < 12) return null + val t0 = beats.first().first; val tN = beats.last().first + if (tN <= t0) return null + val n = ceil((tN - t0) / 0.25 - 1e-9).toInt() + if (n < 16) return null + val y = DoubleArray(n) + var seg = 0 + for (i in 0 until n) { + val t = t0 + 0.25 * i + while (seg < beats.size - 2 && beats[seg + 1].first < t) seg++ + val ta = beats[seg].first; val tb = beats[seg + 1].first + val va = beats[seg].second; val vb = beats[seg + 1].second + y[i] = if (tb <= ta) va else va + ((t - ta) / (tb - ta)).coerceIn(0.0, 1.0) * (vb - va) + } + val mean = y.sum() / n + for (i in 0 until n) y[i] -= mean + val kLo = ceil(0.15 * 0.25 * n).toInt() + val kHi = floor(0.40 * 0.25 * n).toInt() + if (kHi < kLo || kLo < 0) return null + var maxP = 0.0; var sumP = 0.0 + for (k in kLo..kHi) { + var re = 0.0; var im = 0.0 + val w = -2.0 * PI * k / n + for (j in 0 until n) { val a = w * j; re += y[j] * cos(a); im += y[j] * sin(a) } + val p = re * re + im * im + sumP += p + if (p > maxP) maxP = p + } + if (sumP == 0.0) return null + return maxP / sumP + } + + /** A beat window as `features()` builds one: whole-second times, several beats a second or none, 300…2000 ms. */ + private fun windowForTest(rng: java.util.Random, span: Int, density: Int): List> { + val beats = ArrayList>() + for (s in 0 until span) { + if (rng.nextInt(100) >= density) continue + repeat(1 + rng.nextInt(2)) { + val breathing = 60.0 * sin(s * 2 * PI / (3 + rng.nextInt(4))) + val v = 650.0 + breathing + rng.nextInt(400) - 200 + beats.add((1_790_000_000L + s).toDouble() to v.coerceIn(300.0, 2000.0)) + } + } + beats.sortWith(compareBy>({ it.first }, { it.second })) + return beats + } + + @Test + fun tableReturnsTheSameValueAsRecomputingEveryFactor() { + val rng = java.util.Random(0x57A6E) + val dft = HashMap() + var nonNull = 0 + repeat(3_000) { + val span = if (rng.nextInt(5) == 0) 2 + rng.nextInt(208) else 205 + rng.nextInt(5) + val beats = windowForTest(rng, span, 30 + rng.nextInt(71)) + val before = respRegularityBeforeTable(beats) + val now = SleepStagerV2.respRegularity(beats, dft) + assertEquals("window of ${beats.size} beats over $span s", before?.toRawBits(), now?.toRawBits()) + if (before != null) nonNull++ + } + assertTrue("most windows must reach the transform: $nonNull", nonNull > 2_000) + assertTrue("several grid lengths must share the table: ${dft.size}", dft.size > 3) + } + + @Test + fun answerDoesNotDependOnWhatTheTableHeldBefore() { + val rng = java.util.Random(0xD1F7) + val warm = HashMap() + repeat(200) { + val beats = windowForTest(rng, 200 + rng.nextInt(10), 80) + val fromWarm = SleepStagerV2.respRegularity(beats, warm) + val fromCold = SleepStagerV2.respRegularity(beats, HashMap()) + assertEquals(fromCold?.toRawBits(), fromWarm?.toRawBits()) + } + } +}