1/**
2 * The process and its conditional Gibbs sampler: x iid N(0, σ²) → y = h*x
3 * (causal FIR) → z = round(y), and a stationary chain of exact draws of
4 * z_{M+1} given a fixed past z_1..z_M.
5 *
6 * TypeScript port of model.py from the sibling timeseries-entropy package —
7 * keep the two in step. The numpy version updates same-color coordinate
8 * blocks vectorized; here a sweep is a plain sequential scan over the
9 * coordinates, an equally valid systematic-scan Gibbs sweep. This CPU sweep
10 * is the piece a WebGPU backend would replace.
11 */
12import { ndtr, ndtri } from './normal'
13import type { Rng } from './rng'
15/** Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into the
16 * lower tail so the CDF differences keep precision. */
17export function truncatedStdNormal(lo: number, hi: number, rng: Rng): number {
18 const flip = lo > -hi // midpoint above 0 (robust to (-inf, inf) intervals)
19 const a = flip ? -hi : lo
20 const b = flip ? -lo : hi
21 const fa = ndtr(a)
22 const fb = ndtr(b)
23 const u = Math.min(Math.max(fa + (fb - fa) * rng.uniform(), 1e-300), 1 - 1e-16)
24 let x = ndtri(u)
25 if (flip) x = -x
26 return Math.min(Math.max(x, lo), hi)
27}
29/**
30 * Constructing the chain draws the past from the prior; the generating
31 * latents are themselves an exact draw from p(x | z), so the Gibbs chain
32 * starts in stationarity — no burn-in bias, only autocorrelation. Each Gibbs
33 * conditional x_i | rest is N(0, σ²) truncated to the interval read off the
34 * ≤ L constraint boxes x_i appears in. draw(k) advances the chain k steps
35 * (thin sweeps each) and returns k sampled z_{M+1} values, each marginally
36 * distributed exactly as z_{M+1} | z_1..z_M.
37 */
38export class ConditionalChain {
39 private readonly h: Float64Array
40 private readonly sigma: number
41 /** Sweeps per emitted sample. May be reassigned between draws (e.g. probe
42 * at thin = 1, then thin by the measured autocorrelation time);
43 * stationarity is unaffected. */
44 thin: number
45 private readonly rng: Rng
46 private readonly L: number
47 private readonly M: number
48 /** Latents x_0..x_{M+L-2}; boxes live in padded rows so that coordinate i
49 * sees exactly L constraint rows i..i+L-1 (rows outside the data are
50 * unconstrained), with coefficient h[j] in row i+j. */
51 private readonly x: Float64Array
52 private readonly ypad: Float64Array
53 private readonly lo: Float64Array
54 private readonly hi: Float64Array
56 constructor(kernel: Float64Array, sigma: number, past: number, rng: Rng, thin = 1) {
57 if (kernel.length === 0) throw new Error('kernel must be nonempty')
58 if (!(sigma > 0)) throw new Error('sigma must be positive')
59 if (past < 1) throw new Error('past must be >= 1')
60 this.h = kernel
61 this.sigma = sigma
62 this.thin = thin
63 this.rng = rng
64 const L = (this.L = kernel.length)
65 const M = (this.M = past)
67 const x = (this.x = new Float64Array(M + L - 1))
68 for (let i = 0; i < x.length; i++) x[i] = sigma * rng.normal()
70 const P = L - 1
71 this.ypad = new Float64Array(M + 2 * P)
72 this.lo = new Float64Array(M + 2 * P).fill(-Infinity)
73 this.hi = new Float64Array(M + 2 * P).fill(Infinity)
74 this.refreshY()
75 for (let m = 0; m < M; m++) {
76 const z = Math.floor(this.ypad[P + m] + 0.5)
77 this.lo[P + m] = z - 0.5
78 this.hi[P + m] = z + 0.5
79 }
80 }
82 /** ypad[P + m] = y_m = Σ_j h[j] x_{m+L-1-j}, recomputed to kill fp drift. */
83 private refreshY(): void {
84 const { h, x, ypad, L, M } = this
85 const P = L - 1
86 for (let m = 0; m < M; m++) {
87 let y = 0
88 for (let j = 0; j < L; j++) y += h[j] * x[m + L - 1 - j]
89 ypad[P + m] = y
90 }
91 }
93 private sweep(): void {
94 const { h, x, ypad, lo, hi, sigma, rng, L } = this
95 this.refreshY()
96 for (let i = 0; i < x.length; i++) {
97 const xi = x[i]
98 let xlo = -Infinity
99 let xhi = Infinity
100 for (let j = 0; j < L; j++) {
101 const hj = h[j]
102 if (hj === 0) continue
103 const row = i + j
104 const res = ypad[row] - hj * xi
105 const b1 = (lo[row] - res) / hj
106 const b2 = (hi[row] - res) / hj
107 if (hj > 0) {
108 if (b1 > xlo) xlo = b1
109 if (b2 < xhi) xhi = b2
110 } else {
111 if (b2 > xlo) xlo = b2
112 if (b1 < xhi) xhi = b1
113 }
114 }
115 const xn = truncatedStdNormal(xlo / sigma, xhi / sigma, rng) * sigma
116 const d = xn - xi
117 if (d !== 0) {
118 for (let j = 0; j < L; j++) ypad[i + j] += d * h[j]
119 x[i] = xn
120 }
121 }
122 }
124 /** The next k samples of z_{M+1}, continuing the chain. */
125 draw = (k: number): Int32Array => {
126 const { h, x, sigma, rng, L, M } = this
127 const out = new Int32Array(k)
128 for (let s = 0; s < k; s++) {
129 for (let t = 0; t < this.thin; t++) this.sweep()
130 // z_{M+1} = round(c + h[0] · x_free) with x_free ~ N(0, σ²) fresh.
131 let c = 0
132 for (let i = 0; i < L - 1; i++) c += h[L - 1 - i] * x[M + i]
133 out[s] = Math.floor(c + sigma * h[0] * rng.normal() + 0.5)
134 }
135 return out
136 }
137}