sleep: compute the RSA transform's twiddle factors once per night

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..<n { let a = w * Double(j); c[j] = cos(a); s[j] = sin(a) }
+                cosines.append(c)
+                sines.append(s)
+            }
+            self.cosines = cosines
+            self.sines = sines
+        }
+    }
+
     /// 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 nil
     /// when there are too few beats. A direct band-limited DFT (only the ~50 in-band bins are needed).
-    static func respRegularity(_ beats: [(Double, Double)]) -> 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..<n { let a = w * Double(j); re += y[j] * cos(a); im += y[j] * sin(a) }
+            let c = table.cosines[k - kLo], s = table.sines[k - kLo]
+            for j in 0..<n { re += y[j] * c[j]; im += y[j] * s[j] }
             let p = re * re + im * im
             sumP += p
             if p > maxP { maxP = p }

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..<n {
+            let t = t0 + 0.25 * Double(i)
+            while seg < beats.count - 2 && beats[seg + 1].0 < t { seg += 1 }
+            let ta = beats[seg].0, tb = beats[seg + 1].0
+            let va = beats[seg].1, vb = beats[seg + 1].1
+            y[i] = tb <= ta ? va : va + min(max((t - ta) / (tb - ta), 0), 1) * (vb - va)
+        }
+        let mean = y.reduce(0, +) / Double(n)
+        for i in 0..<n { y[i] -= mean }
+        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 }
+        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..<n { let a = w * Double(j); re += y[j] * cos(a); im += y[j] * sin(a) }
+            let p = re * re + im * im
+            sumP += p
+            if p > 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..<span where Int(rng.next() % 100) < density {
+            for _ in 0..<(1 + Int(rng.next() % 2)) {
+                let breathing = 60.0 * sin(Double(s) * 2 * .pi / Double(3 + rng.next() % 4))
+                let v = 650.0 + breathing + Double(rng.next() % 400) - 200
+                beats.append((Double(1_790_000_000 + s), min(max(v, 300), 2000)))
+            }
+        }
+        beats.sort { $0.0 != $1.0 ? $0.0 < $1.0 : $0.1 < $1.1 }
+        return beats
+    }
+
+    func testTheTableReturnsTheSameValueAsRecomputingEveryFactor() {
+        var rng = SplitMixForDFTTest(seed: 0x57A6E)
+        var dft: [Int: SleepStagerV2.RespDFT] = [:]
+        var compared = 0, nonNil = 0
+        for _ in 0..<3_000 {
+            // Mostly full-width windows (as mid-night), some short or sparse ones (session edges, dropout).
+            let span = rng.next() % 5 == 0 ? 2 + Int(rng.next() % 208) : 205 + Int(rng.next() % 5)
+            let beats = windowForTest(&rng, span: span, density: 30 + Int(rng.next() % 71))
+            let before = respRegularityBeforeTable(beats)
+            let now = SleepStagerV2.respRegularity(beats, dft: &dft)
+            XCTAssertEqual(now?.bitPattern, before?.bitPattern, "window of \(beats.count) beats over \(span) s")
+            compared += 1
+            if before != nil { nonNil += 1 }
+        }
+        XCTAssertEqual(compared, 3_000)
+        XCTAssertGreaterThan(nonNil, 2_000, "most windows must reach the transform, or equality proves little")
+        XCTAssertGreaterThan(dft.count, 3, "several grid lengths must share the table")
+    }
+
+    /// A fresh table and a warm one give the same answer: what a window returns cannot depend on which
+    /// windows came before it in the night.
+    func testTheAnswerDoesNotDependOnWhatTheTableHeldBefore() {
+        var rng = SplitMixForDFTTest(seed: 0xD1F7)
+        var warm: [Int: SleepStagerV2.RespDFT] = [:]
+        for _ in 0..<200 {
+            let beats = windowForTest(&rng, span: 200 + Int(rng.next() % 10), density: 80)
+            var cold: [Int: SleepStagerV2.RespDFT] = [:]
+            let fromWarm = SleepStagerV2.respRegularity(beats, dft: &warm)
+            let fromCold = SleepStagerV2.respRegularity(beats, dft: &cold)
+            XCTAssertEqual(fromWarm?.bitPattern, fromCold?.bitPattern)
+        }
+    }
+}
+
+/// Deterministic, so a failure names a reproducible window.
+struct SplitMixForDFTTest {
+    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)
+    }
+}

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"}
   }

android/app/src/main/java/com/noop/analytics/SleepStagerV2.kt

@@ -340,6 +340,8 @@ object SleepStagerV2 {
         )
         val raws = ArrayList<Raw>()
         val allJerks = ArrayList<Double>()
+        // The RSA transform's twiddle factors, per grid length, for this night only (see [RespDft]).
+        val respDft = HashMap<Int, RespDft>()
         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<DoubleArray> = Array(kHi - kLo + 1) { DoubleArray(n) }
+        val sines: Array<DoubleArray> = 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<Pair<Double, Double>>): Double? {
+    internal fun respRegularity(beats: List<Pair<Double, Double>>, dft: MutableMap<Int, RespDft>): 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

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<Pair<Double, Double>>): 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<Pair<Double, Double>> {
+        val beats = ArrayList<Pair<Double, Double>>()
+        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<Pair<Double, Double>>({ it.first }, { it.second }))
+        return beats
+    }
+
+    @Test
+    fun tableReturnsTheSameValueAsRecomputingEveryFactor() {
+        val rng = java.util.Random(0x57A6E)
+        val dft = HashMap<Int, SleepStagerV2.RespDft>()
+        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<Int, SleepStagerV2.RespDft>()
+        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())
+        }
+    }
+}