Estimate the reference rate R in the browser
TypeScript port of the timeseries-entropy estimator in src/entropy
(hand-synced with the package): Gibbs conditional sampler with a
sequential-scan sweep, Rhee-Glynn randomized telescoping, Cody erfc +
Halley-polished Acklam ndtri (validated against scipy to 1e-13 in the
sampler-critical range), xoshiro128** RNG (a past consumes ~5e7 draws,
too many for the 2^32-period generators used elsewhere in the app).
An estimate button runs it in a worker, one independent past at a time
until stopped - live mean +/- se in the stat row, dashed R line and
table row on the chart - resuming deterministically after a stop and
invalidating on any model change. Defaults match the CLI (M, n0, r,
reps), so browser runs and command-line runs target the same estimand;
end-to-end estimates validated against the Python package.
12 changed files+674−16
README.mdmodified+15−5View file
@@ -17,13 +17,18 @@ entropy coder cannot beat, which ANS misses by 1–2% (its symbol table plus its
1717 own arithmetic loss).
1818
1919 The reference rate R — the entropy rate of the process, the bits/sample limit
20-no lossless method can beat — comes from the companion
20+no lossless method can beat — is estimated in the browser by the method of the
21+companion
2122 [timeseries-entropy](https://github.com/concept-collection/timeseries-entropy)
2223 package: an unbiased Monte-Carlo estimator of H(z_next | a long past), by Gibbs
2324 sampling the latent Gaussian under the rounding constraints and applying
24-Rhee–Glynn randomized telescoping to the sampled chain. The app shows the exact
25-command to run it at the current settings; estimating R in the browser is
26-planned, and until then the UI shows a placeholder for it.
25+Rhee–Glynn randomized telescoping to the sampled chain. A button starts a web
26+worker that averages one independent past at a time (live mean ± se, dashed
27+line on the chart) until stopped; the app also shows the exact command to run
28+the Python original at the same settings as an independent check. The
29+in-browser code in `src/entropy/` is a hand-synced TypeScript port of that
30+package — change one, change the other. A WebGPU Gibbs sweep may replace the
31+scalar one someday; the sweep is isolated so it can be swapped.
2732
2833 ## Run it
2934
@@ -38,10 +43,15 @@ npm run dev
3843 src/model/ the latent source (fixed seeded randomness indexed by sample
3944 position, convolved zero-phase with the kernel on demand)
4045 and the FIR presets
46+src/entropy/ the unbiased reference-rate estimator: hand-synced TypeScript
47+ port of the timeseries-entropy package (Gibbs conditional
48+ sampler, Rhee–Glynn telescoping, Cody erfc / Acklam ndtri,
49+ xoshiro128** RNG)
4150 src/compress/ lossless codecs run in the browser: zlib (fflate), zstd (wasm),
4251 ans.ts (a bit-identical port of simple_ans), and FLAC-style
4352 integer LPC; borrowed from entropy-quantized-linear-transform
44-src/worker/ the codecs run off the main thread on a debounced parameter set
53+src/worker/ the codecs and the estimator run off the main thread; the
54+ estimator worker refines one past at a time until terminated
4555 src/components/ controls, filter plots, signal canvas, compression chart,
4656 and the reference-rate method note
4757 ```
src/App.tsxmodified+29−7View file
@@ -5,6 +5,7 @@ import FilterViz from './components/FilterViz'
55 import ScrollingView from './components/ScrollingView'
66 import CompressionChart from './components/CompressionChart'
77 import MethodNote from './components/MethodNote'
8+import { useReferenceRate } from './components/useReferenceRate'
89 import { DEFAULT_SPEC, clampSpec, designKernel, kernelNorm } from './model/filters'
910 import { LATENT_SEED } from './model/latent'
1011 import { DEFAULT_LPC_ORDER, LPC_ORDERS } from './compress/codecs'
@@ -129,6 +130,7 @@ export default function App() {
129130 const kernel = useMemo(() => designKernel(spec, sampleRateHz), [spec, sampleRateHz])
130131 const sigmaY = useMemo(() => sigma * kernelNorm(kernel), [kernel, sigma])
131132 const compression = useCompression(kernel, sigma, lpcOrder, blockSize)
133+ const refRate = useReferenceRate(kernel, sigma)
132134
133135 return (
134136 <div className="app">
@@ -175,12 +177,18 @@ export default function App() {
175177 <small>steps</small>
176178 </span>
177179 </div>
178- {/* Placeholder: R will be estimated in the browser by the unbiased
179- estimator; until then the command below produces it locally. */}
180180 <div className="stat">
181181 <span className="label">reference rate R</span>
182182 <span className="value">
183- — <small>bits/sample</small>
183+ {refRate.mean !== null ? refRate.mean.toFixed(2) : '—'}
184+ {refRate.se !== null && <small> ± {refRate.se.toFixed(2)}</small>}{' '}
185+ <small>bits/sample</small>
186+ </span>
187+ </div>
188+ <div className="stat">
189+ <span className="label">implied best ratio</span>
190+ <span className="value">
191+ {refRate.mean !== null && refRate.mean > 0 ? `${(16 / refRate.mean).toFixed(2)}×` : '—'}
184192 </span>
185193 </div>
186194 </div>
@@ -214,6 +222,7 @@ export default function App() {
214222 <CompressionChart
215223 results={compression.results}
216224 bounds={compression.bounds}
225+ refBits={refRate.mean}
217226 computing={compression.computing}
218227 />
219228 )}
@@ -223,12 +232,25 @@ export default function App() {
223232 coefficients). Baseline is raw int16 (16 bits/sample). The hollow bar in each group is
224233 that group's entropy limit — the order-0 entropy of the stream being coded, which no
225234 per-sample entropy coder can beat and ANS falls short of by its symbol table plus its
226- own arithmetic loss. The reference rate R — the entropy rate of the process itself,
227- the limit no lossless method can beat — is not yet computed in the browser; the
228- command below estimates it locally (see the method section at the bottom).
235+ own arithmetic loss. The dashed line, once estimated, is the reference rate R — the
236+ entropy rate of the process itself, the limit no lossless method whatsoever can beat
237+ (see the method section at the bottom).
229238 </p>
239+ <div className="estimate-row">
240+ <button onClick={refRate.running ? refRate.stop : refRate.start}>
241+ {refRate.running ? 'stop' : refRate.perPast.length > 0 ? 'refine R further' : 'estimate R in this browser'}
242+ </button>
243+ <span className="estimate-status">
244+ {refRate.running
245+ ? `${refRate.perPast.length} independent pasts averaged, M = ${refRate.past}` +
246+ (refRate.progress ? ` · ${refRate.progress}` : '')
247+ : refRate.perPast.length > 0
248+ ? `${refRate.perPast.length} independent pasts averaged, M = ${refRate.past}`
249+ : `unbiased Monte-Carlo conditioning on M = ${refRate.past} past samples; refines until stopped`}
250+ </span>
251+ </div>
230252 <CopyableCommand
231- label="reference rate R by unbiased Monte-Carlo (runs locally):"
253+ label="or cross-check R from the command line:"
232254 command={mcCommand(sigma, spec, sampleRateHz)}
233255 />
234256 </section>
src/app.cssmodified+29−0View file
@@ -443,6 +443,35 @@ body {
443443 font-weight: 500;
444444 }
445445
446+.estimate-row {
447+ display: flex;
448+ flex-wrap: wrap;
449+ align-items: center;
450+ gap: 6px 10px;
451+ margin-top: 12px;
452+ font-size: 12px;
453+}
454+
455+.estimate-row button {
456+ background: var(--surface);
457+ color: var(--ink);
458+ border: 1px solid var(--baseline);
459+ border-radius: 6px;
460+ padding: 4px 12px;
461+ font: inherit;
462+ font-weight: 600;
463+ cursor: pointer;
464+}
465+
466+.estimate-row button:hover {
467+ border-color: var(--series-1);
468+}
469+
470+.estimate-status {
471+ color: var(--muted);
472+ font-variant-numeric: tabular-nums;
473+}
474+
446475 .command-row {
447476 display: flex;
448477 flex-wrap: wrap;
src/components/CompressionChart.tsxmodified+46−2View file
@@ -89,6 +89,8 @@ function buildRows(results: CodecResult[], bounds: BoundResult[]): Row[] {
8989 export default function CompressionChart(props: {
9090 results: CodecResult[]
9191 bounds: BoundResult[]
92+ /** The browser-estimated reference rate R, once at least one past is in. */
93+ refBits: number | null
9294 computing: boolean
9395 }) {
9496 const ref = useRef<HTMLDivElement>(null)
@@ -97,7 +99,7 @@ export default function CompressionChart(props: {
9799 const [tip, setTip] = useState<Tip | null>(null)
98100 const [hovered, setHovered] = useState<string | null>(null)
99101
100- const { results, bounds } = props
102+ const { results, bounds, refBits } = props
101103 if (results.length === 0) {
102104 return <p className="card-note">Computing compression on the first block…</p>
103105 }
@@ -105,8 +107,12 @@ export default function CompressionChart(props: {
105107 const rows = buildRows(results, bounds)
106108 const value = (r: { bitsPerSample: number; ratio: number }) =>
107109 metric === 'bits' ? r.bitsPerSample : r.ratio
110+ const refValue =
111+ refBits !== null && refBits > 0 ? (metric === 'bits' ? refBits : 16 / refBits) : null
108112 const xMax =
109- metric === 'bits' ? Math.max(16, ...rows.map(value)) * 1.02 : Math.max(...rows.map(value)) * 1.1
113+ metric === 'bits'
114+ ? Math.max(16, ...rows.map(value), refValue ?? 0) * 1.02
115+ : Math.max(...rows.map(value), refValue ?? 0) * 1.1
110116
111117 const plotW = width - LABEL_W - RIGHT_PAD
112118 const height = AXIS_H + GROUPS.length * (GROUP_H + ROWS_PER_GROUP * ROW_H) + 6
@@ -124,6 +130,14 @@ export default function CompressionChart(props: {
124130 setTip({ x: e.clientX - box.left, y: e.clientY - box.top, row })
125131 }
126132
133+ const refX = refValue !== null ? xOf(refValue) : 0
134+ const refLabel =
135+ refBits !== null && refBits > 0
136+ ? metric === 'bits'
137+ ? `R = ${refBits.toFixed(2)}`
138+ : `R ⇒ ${(16 / refBits).toFixed(2)}×`
139+ : ''
140+
127141 return (
128142 <div>
129143 <div className="chart-header">
@@ -224,6 +238,28 @@ export default function CompressionChart(props: {
224238 </g>
225239 )
226240 })}
241+ {refValue !== null && (
242+ <g>
243+ <line
244+ x1={refX}
245+ x2={refX}
246+ y1={AXIS_H - 2}
247+ y2={height - 4}
248+ stroke="var(--ink-2)"
249+ strokeWidth={1.5}
250+ strokeDasharray="5 4"
251+ />
252+ <text
253+ x={refX + (refX > width - 150 ? -6 : 6)}
254+ y={AXIS_H + 10}
255+ textAnchor={refX > width - 150 ? 'end' : 'start'}
256+ className="bar-value"
257+ fill="var(--ink)"
258+ >
259+ {refLabel}
260+ </text>
261+ </g>
262+ )}
227263 </svg>
228264 {tip && (
229265 <div className="viz-tooltip" style={{ left: tip.x + 14, top: tip.y - 8 }}>
@@ -260,6 +296,14 @@ export default function CompressionChart(props: {
260296 <td>{r.ratio.toFixed(3)}</td>
261297 </tr>
262298 ))}
299+ {refBits !== null && refBits > 0 && (
300+ <tr>
301+ <td>reference rate R (Monte-Carlo)</td>
302+ <td>—</td>
303+ <td>{refBits.toFixed(3)}</td>
304+ <td>{(16 / refBits).toFixed(3)}</td>
305+ </tr>
306+ )}
263307 </tbody>
264308 </table>
265309 </details>
src/components/MethodNote.tsxmodified+4−2View file
@@ -33,8 +33,10 @@ export default function MethodNote() {
3333 honest standard error.
3434 </p>
3535 <p className="card-note">
36- Running the estimator in the browser is planned; until then, the command under the
37- compression chart runs it locally at the current settings.
36+ The estimate button under the compression chart runs exactly this method in a web
37+ worker — a TypeScript port of the package (src/entropy, hand-synced), one independent
38+ past at a time until stopped. The command line runs the Python original at the same
39+ settings for an independent check.
3840 </p>
3941 </div>
4042 )
src/components/useReferenceRate.tsadded+90−0View file
@@ -0,0 +1,90 @@
1+import { useEffect, useRef, useState } from 'react'
2+import { defaultPast } from '../entropy'
3+import type { EntropyRequest, EntropyUpdate } from '../worker/entropyWorker'
4+
5+/** Fixed base seed: resumed runs continue the same per-past seed sequence,
6+ * so a stop/start pair reproduces an uninterrupted run exactly. */
7+const BASE_SEED = 20260731
8+
9+export interface ReferenceRate {
10+ /** One unbiased estimate per independent past, in completion order. */
11+ perPast: number[]
12+ mean: number | null
13+ se: number | null
14+ running: boolean
15+ /** "past 3 · rep 5/8" while computing. */
16+ progress: string | null
17+ /** Conditioning window M. */
18+ past: number
19+ start: () => void
20+ stop: () => void
21+}
22+
23+/**
24+ * The in-browser reference-rate estimate: a worker refines it (one
25+ * independent past at a time) until stopped, and any change to the model
26+ * invalidates both the values and a run in flight.
27+ */
28+export function useReferenceRate(kernel: Float64Array, sigma: number): ReferenceRate {
29+ const [perPast, setPerPast] = useState<number[]>([])
30+ const [running, setRunning] = useState(false)
31+ const [progress, setProgress] = useState<string | null>(null)
32+ const workerRef = useRef<Worker | null>(null)
33+ const perPastRef = useRef<number[]>([])
34+ const past = defaultPast(kernel.length)
35+
36+ useEffect(() => {
37+ workerRef.current?.terminate()
38+ workerRef.current = null
39+ perPastRef.current = []
40+ setPerPast([])
41+ setRunning(false)
42+ setProgress(null)
43+ }, [kernel, sigma])
44+
45+ useEffect(() => () => workerRef.current?.terminate(), [])
46+
47+ const start = () => {
48+ if (workerRef.current) return
49+ const worker = new Worker(new URL('../worker/entropyWorker.ts', import.meta.url), {
50+ type: 'module',
51+ })
52+ worker.onmessage = (e: MessageEvent<EntropyUpdate>) => {
53+ const u = e.data
54+ if (u.type === 'past') {
55+ perPastRef.current = [...perPastRef.current, u.value]
56+ setPerPast(perPastRef.current)
57+ } else {
58+ setProgress(`past ${u.pastIndex + 1} · rep ${u.repsDone}/${u.reps}`)
59+ }
60+ }
61+ workerRef.current = worker
62+ const request: EntropyRequest = {
63+ kernel,
64+ sigma,
65+ past,
66+ seed: BASE_SEED,
67+ startPast: perPastRef.current.length,
68+ }
69+ worker.postMessage(request)
70+ setRunning(true)
71+ setProgress(null)
72+ }
73+
74+ const stop = () => {
75+ workerRef.current?.terminate()
76+ workerRef.current = null
77+ setRunning(false)
78+ setProgress(null)
79+ }
80+
81+ const n = perPast.length
82+ const mean = n > 0 ? perPast.reduce((a, b) => a + b, 0) / n : null
83+ let se: number | null = null
84+ if (mean !== null && n > 1) {
85+ const v = perPast.reduce((a, b) => a + (b - mean) * (b - mean), 0) / (n - 1)
86+ se = Math.sqrt(v / n)
87+ }
88+
89+ return { perPast, mean, se, running, progress, past, start, stop }
90+}
src/entropy/estimator.tsadded+68−0View file
@@ -0,0 +1,68 @@
1+/**
2+ * Rhee–Glynn (randomized telescoping) unbiased entropy estimation.
3+ *
4+ * TypeScript port of estimator.py from the sibling timeseries-entropy
5+ * package — keep the two in step. See that file for the full derivation;
6+ * in short: plug-in entropies of blocks whose sizes double per level form a
7+ * telescoping sum via the antithetic correction
8+ * Δ_m = h(B_m) − [h(B_m¹) + h(B_m²)]/2; truncating at a random level N with
9+ * P(N ≥ m) = 2^(−r m) and reweighting by the survival probabilities gives
10+ * an estimator whose expectation is exactly the entropy of the stationary
11+ * marginal, despite the bias of every finite-block plug-in estimate and any
12+ * autocorrelation of the draws. All entropies are in bits.
13+ */
14+import type { Rng } from './rng'
15+
16+/** draw(k) returns the next k consecutive samples of a stationary chain. */
17+export type Draw = (k: number) => Int32Array
18+
19+function entropyFromCounts(counts: Map<number, number>, n: number): number {
20+ let s = 0
21+ for (const c of counts.values()) s += c * Math.log2(c)
22+ return Math.log2(n) - s / n
23+}
24+
25+function countInto(counts: Map<number, number>, seg: Int32Array): void {
26+ for (let i = 0; i < seg.length; i++) {
27+ counts.set(seg[i], (counts.get(seg[i]) ?? 0) + 1)
28+ }
29+}
30+
31+/** h(B_0) and [Δ_1 … Δ_levels] over one growing block; counts merge upward
32+ * so the cost is linear in the n0 * 2**levels samples drawn. */
33+function telescope(draw: Draw, n0: number, levels: number): { h0: number; deltas: number[] } {
34+ const counts = new Map<number, number>()
35+ countInto(counts, draw(n0))
36+ let size = n0
37+ const h0 = entropyFromCounts(counts, size)
38+ let hPrev = h0
39+ const deltas: number[] = []
40+ for (let m = 0; m < levels; m++) {
41+ const half = new Map<number, number>()
42+ countInto(half, draw(size))
43+ const h2 = entropyFromCounts(half, size)
44+ for (const [v, c] of half) counts.set(v, (counts.get(v) ?? 0) + c)
45+ size *= 2
46+ const hFull = entropyFromCounts(counts, size)
47+ deltas.push(hFull - 0.5 * (hPrev + h2))
48+ hPrev = hFull
49+ }
50+ return { h0, deltas }
51+}
52+
53+/**
54+ * One randomized-telescoping realization of the marginal entropy (bits).
55+ * Consumes n0 * 2**N samples with P(N ≥ m) = 2^(−r m); average many
56+ * realizations (they may continue one chain back-to-back) — each has
57+ * expectation exactly H. r = 1.5 suits the typical Δ second-moment decay
58+ * of 2^(−2m); finite work needs r > 1, finite variance needs decay > r.
59+ */
60+export function unbiasedEntropy(draw: Draw, n0: number, r: number, rng: Rng): number {
61+ const rho = 2 ** -r
62+ let N = 0
63+ while (rng.uniform() < rho) N++
64+ const { h0, deltas } = telescope(draw, n0, N)
65+ let est = h0
66+ for (let m = 1; m <= N; m++) est += deltas[m - 1] * 2 ** (r * m)
67+ return est
68+}
src/entropy/index.tsadded+51−0View file
@@ -0,0 +1,51 @@
1+/**
2+ * In-browser unbiased estimation of H(z_{M+1} | z_1..z_M) — the reference
3+ * rate R — for the app's model. Hand-synced TypeScript port of the sibling
4+ * timeseries-entropy package (see the headers of estimator.ts / model.ts);
5+ * defaults match its CLI so browser runs and `timeseries-entropy` runs
6+ * target the same estimand with the same variance behavior.
7+ */
8+import { unbiasedEntropy } from './estimator'
9+import { ConditionalChain } from './model'
10+import { Rng } from './rng'
11+
12+export { unbiasedEntropy } from './estimator'
13+export { ConditionalChain, truncatedStdNormal } from './model'
14+export { Rng } from './rng'
15+export { ndtr, ndtri, erfc } from './normal'
16+
17+export const N0 = 128
18+export const R_EXPONENT = 1.5
19+export const REPS_PER_PAST = 8
20+export const THIN = 1
21+
22+/** The conditioning window M at a given kernel length, as in the CLI. */
23+export function defaultPast(kernelLength: number): number {
24+ return Math.max(512, 4 * kernelLength)
25+}
26+
27+/** The seed for one independent past, derived so past i is reproducible
28+ * whether or not the run was stopped and resumed in between. */
29+export function pastSeed(baseSeed: number, pastIndex: number): number {
30+ return (baseSeed + Math.imul(0x9e3779b9, pastIndex + 1)) >>> 0
31+}
32+
33+/** One independent past's unbiased estimate: a fresh stationary chain,
34+ * averaged over reps randomized-telescoping realizations run back-to-back
35+ * on it. Averaging these over pasts estimates H(z_{M+1} | z_1..z_M). */
36+export function estimateOnePast(
37+ kernel: Float64Array,
38+ sigma: number,
39+ past: number,
40+ seed: number,
41+ onRep?: (repsDone: number) => void,
42+): number {
43+ const rng = new Rng(seed)
44+ const chain = new ConditionalChain(kernel, sigma, past, rng, THIN)
45+ let sum = 0
46+ for (let rep = 0; rep < REPS_PER_PAST; rep++) {
47+ sum += unbiasedEntropy(chain.draw, N0, R_EXPONENT, rng)
48+ onRep?.(rep + 1)
49+ }
50+ return sum / REPS_PER_PAST
51+}
src/entropy/model.tsadded+134−0View file
@@ -0,0 +1,134 @@
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+ */
12+import { ndtr, ndtri } from './normal'
13+import type { Rng } from './rng'
14+
15+/** Standard normal truncated to [lo, hi], by inverse CDF. Mirrored into the
16+ * lower tail so the CDF differences keep precision. */
17+export 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+}
28+
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+ */
38+export class ConditionalChain {
39+ private readonly h: Float64Array
40+ private readonly sigma: number
41+ private readonly thin: number
42+ private readonly rng: Rng
43+ private readonly L: number
44+ private readonly M: number
45+ /** Latents x_0..x_{M+L-2}; boxes live in padded rows so that coordinate i
46+ * sees exactly L constraint rows i..i+L-1 (rows outside the data are
47+ * unconstrained), with coefficient h[j] in row i+j. */
48+ private readonly x: Float64Array
49+ private readonly ypad: Float64Array
50+ private readonly lo: Float64Array
51+ private readonly hi: Float64Array
52+
53+ constructor(kernel: Float64Array, sigma: number, past: number, rng: Rng, thin = 1) {
54+ if (kernel.length === 0) throw new Error('kernel must be nonempty')
55+ if (!(sigma > 0)) throw new Error('sigma must be positive')
56+ if (past < 1) throw new Error('past must be >= 1')
57+ this.h = kernel
58+ this.sigma = sigma
59+ this.thin = thin
60+ this.rng = rng
61+ const L = (this.L = kernel.length)
62+ const M = (this.M = past)
63+
64+ const x = (this.x = new Float64Array(M + L - 1))
65+ for (let i = 0; i < x.length; i++) x[i] = sigma * rng.normal()
66+
67+ const P = L - 1
68+ this.ypad = new Float64Array(M + 2 * P)
69+ this.lo = new Float64Array(M + 2 * P).fill(-Infinity)
70+ this.hi = new Float64Array(M + 2 * P).fill(Infinity)
71+ this.refreshY()
72+ for (let m = 0; m < M; m++) {
73+ const z = Math.floor(this.ypad[P + m] + 0.5)
74+ this.lo[P + m] = z - 0.5
75+ this.hi[P + m] = z + 0.5
76+ }
77+ }
78+
79+ /** ypad[P + m] = y_m = Σ_j h[j] x_{m+L-1-j}, recomputed to kill fp drift. */
80+ private refreshY(): void {
81+ const { h, x, ypad, L, M } = this
82+ const P = L - 1
83+ for (let m = 0; m < M; m++) {
84+ let y = 0
85+ for (let j = 0; j < L; j++) y += h[j] * x[m + L - 1 - j]
86+ ypad[P + m] = y
87+ }
88+ }
89+
90+ private sweep(): void {
91+ const { h, x, ypad, lo, hi, sigma, rng, L } = this
92+ this.refreshY()
93+ for (let i = 0; i < x.length; i++) {
94+ const xi = x[i]
95+ let xlo = -Infinity
96+ let xhi = Infinity
97+ for (let j = 0; j < L; j++) {
98+ const hj = h[j]
99+ if (hj === 0) continue
100+ const row = i + j
101+ const res = ypad[row] - hj * xi
102+ const b1 = (lo[row] - res) / hj
103+ const b2 = (hi[row] - res) / hj
104+ if (hj > 0) {
105+ if (b1 > xlo) xlo = b1
106+ if (b2 < xhi) xhi = b2
107+ } else {
108+ if (b2 > xlo) xlo = b2
109+ if (b1 < xhi) xhi = b1
110+ }
111+ }
112+ const xn = truncatedStdNormal(xlo / sigma, xhi / sigma, rng) * sigma
113+ const d = xn - xi
114+ if (d !== 0) {
115+ for (let j = 0; j < L; j++) ypad[i + j] += d * h[j]
116+ x[i] = xn
117+ }
118+ }
119+ }
120+
121+ /** The next k samples of z_{M+1}, continuing the chain. */
122+ draw = (k: number): Int32Array => {
123+ const { h, x, sigma, rng, L, M } = this
124+ const out = new Int32Array(k)
125+ for (let s = 0; s < k; s++) {
126+ for (let t = 0; t < this.thin; t++) this.sweep()
127+ // z_{M+1} = round(c + h[0] · x_free) with x_free ~ N(0, σ²) fresh.
128+ let c = 0
129+ for (let i = 0; i < L - 1; i++) c += h[L - 1 - i] * x[M + i]
130+ out[s] = Math.floor(c + sigma * h[0] * rng.normal() + 0.5)
131+ }
132+ return out
133+ }
134+}
src/entropy/normal.tsadded+108−0View file
@@ -0,0 +1,108 @@
1+/**
2+ * Standard normal CDF and inverse for the truncated-Gaussian sampler,
3+ * standing in for scipy.special.ndtr/ndtri in the Python package.
4+ *
5+ * Φ comes from W. J. Cody's rational erfc (SPECFUN calerf), full double
6+ * precision with correct relative accuracy deep in the lower tail — which is
7+ * what the mirrored truncated sampler relies on. Φ⁻¹ is Acklam's
8+ * approximation (~1.15e-9 relative) polished by one Halley step against this
9+ * Φ, which reaches ~1e-13 everywhere the step is taken: down to x ≈ -37.5
10+ * (below, exp(x²/2) overflows) and up to p = 1 - 1e-9 (beyond, the spacing
11+ * of doubles near 1 limits any method to about what Acklam already gives).
12+ */
13+
14+const THRESH = 0.46875
15+const SQRPI = 5.6418958354775628695e-1 // 1/√π
16+
17+const ERF_A = [3.1611237438705656, 1.13864154151050156e2, 3.77485237685302021e2, 3.20937758913846947e3, 1.85777706184603153e-1]
18+const ERF_B = [2.36012909523441209e1, 2.44024637934444173e2, 1.28261652607737228e3, 2.84423683343917062e3]
19+const ERF_C = [5.64188496988670089e-1, 8.88314979438837594, 6.61191906371416295e1, 2.98635138197400131e2, 8.8195222124176909e2, 1.71204761263407058e3, 2.05107837782607147e3, 1.23033935479799725e3, 2.15311535474403846e-8]
20+const ERF_D = [1.57449261107098347e1, 1.17693950891312499e2, 5.37181101862009858e2, 1.62138957456669019e3, 3.29079923573345963e3, 4.36261909014324716e3, 3.43936767414372164e3, 1.23033935480374942e3]
21+const ERF_P = [3.05326634961232344e-1, 3.60344899949804439e-1, 1.25781726111229246e-1, 1.60837851487422766e-2, 6.58749161529837803e-4, 1.63153871373020978e-2]
22+const ERF_Q = [2.56852019228982242, 1.87295284992346047, 5.27905102951428412e-1, 6.05183413124413191e-2, 2.33520497626869185e-3]
23+
24+/** erf(x) for |x| ≤ 0.46875. */
25+function erfSmall(x: number): number {
26+ const z = x * x
27+ let num = ERF_A[4] * z
28+ let den = z
29+ for (let i = 0; i < 3; i++) {
30+ num = (num + ERF_A[i]) * z
31+ den = (den + ERF_B[i]) * z
32+ }
33+ return (x * (num + ERF_A[3])) / (den + ERF_B[3])
34+}
35+
36+/** erfc(y) for y > 0.46875 (with the split exp(-y²) trick for accuracy). */
37+function erfcLarge(y: number): number {
38+ let result: number
39+ if (y <= 4) {
40+ let num = ERF_C[8] * y
41+ let den = y
42+ for (let i = 0; i < 7; i++) {
43+ num = (num + ERF_C[i]) * y
44+ den = (den + ERF_D[i]) * y
45+ }
46+ result = (num + ERF_C[7]) / (den + ERF_D[7])
47+ } else {
48+ const z = 1 / (y * y)
49+ let num = ERF_P[5] * z
50+ let den = z
51+ for (let i = 0; i < 4; i++) {
52+ num = (num + ERF_P[i]) * z
53+ den = (den + ERF_Q[i]) * z
54+ }
55+ result = (z * (num + ERF_P[4])) / (den + ERF_Q[4])
56+ result = (SQRPI - result) / y
57+ }
58+ const ysq = Math.trunc(y * 16) / 16
59+ const del = (y - ysq) * (y + ysq)
60+ return Math.exp(-ysq * ysq) * Math.exp(-del) * result
61+}
62+
63+/** Complementary error function, double precision over the whole line. */
64+export function erfc(x: number): number {
65+ const y = Math.abs(x)
66+ if (y <= THRESH) return 1 - erfSmall(x)
67+ const tail = y > 26.6 ? 0 : erfcLarge(y)
68+ return x < 0 ? 2 - tail : tail
69+}
70+
71+/** Standard normal CDF Φ(t), accurate relative to its size in both tails. */
72+export function ndtr(t: number): number {
73+ return 0.5 * erfc(-t / Math.SQRT2)
74+}
75+
76+const ACK_A = [-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2, 1.38357751867269e2, -3.066479806614716e1, 2.506628277459239]
77+const ACK_B = [-5.447609879822406e1, 1.615858368580409e2, -1.556989798598866e2, 6.680131188771972e1, -1.328068155288572e1]
78+const ACK_C = [-7.784894002430293e-3, -3.223964580411365e-1, -2.400758277161838, -2.549732539343734, 4.374664141464968, 2.938163982698783]
79+const ACK_D = [7.784695709041462e-3, 3.224671290700398e-1, 2.445134137142996, 3.754408661907416]
80+
81+const SQRT_2PI = Math.sqrt(2 * Math.PI)
82+
83+/** Standard normal quantile Φ⁻¹(p) for p in (0, 1). */
84+export function ndtri(p: number): number {
85+ const pLow = 0.02425
86+ let x: number
87+ if (p >= pLow && p <= 1 - pLow) {
88+ const q = p - 0.5
89+ const r = q * q
90+ x =
91+ ((((((ACK_A[0] * r + ACK_A[1]) * r + ACK_A[2]) * r + ACK_A[3]) * r + ACK_A[4]) * r + ACK_A[5]) * q) /
92+ (((((ACK_B[0] * r + ACK_B[1]) * r + ACK_B[2]) * r + ACK_B[3]) * r + ACK_B[4]) * r + 1)
93+ } else {
94+ const lower = p < pLow
95+ const q = Math.sqrt(-2 * Math.log(lower ? p : 1 - p))
96+ x =
97+ (((((ACK_C[0] * q + ACK_C[1]) * q + ACK_C[2]) * q + ACK_C[3]) * q + ACK_C[4]) * q + ACK_C[5]) /
98+ ((((ACK_D[0] * q + ACK_D[1]) * q + ACK_D[2]) * q + ACK_D[3]) * q + 1)
99+ if (!lower) x = -x
100+ }
101+ // One Halley step against the accurate Φ (see the header for the guards).
102+ if (p <= 1 - 1e-9 && x > -37.5) {
103+ const e = ndtr(x) - p
104+ const u = e * SQRT_2PI * Math.exp((x * x) / 2)
105+ x -= u / (1 + (x * u) / 2)
106+ }
107+ return x
108+}
src/entropy/rng.tsadded+64−0View file
@@ -0,0 +1,64 @@
1+/**
2+ * Seeded RNG for the entropy estimator: xoshiro128** uniforms (128-bit
3+ * state, period 2^128 − 1, seeded through splitmix32) with Box–Muller
4+ * normals. A single past consumes ~5e7 uniforms, so the 2^32-period
5+ * generators used elsewhere in the app are not enough here: their seeds are
6+ * mere offsets into one cyclic stream, and long refine-until-stopped runs
7+ * would reuse randomness across "independent" pasts. Owned by this module so
8+ * that src/entropy stays self-contained — the plan is to maintain it in step
9+ * with (and eventually extract it back into) the timeseries-entropy package.
10+ */
11+export class Rng {
12+ private s0: number
13+ private s1: number
14+ private s2: number
15+ private s3: number
16+ private spare: number | null = null
17+
18+ constructor(seed: number) {
19+ // splitmix32 stream fills the state; any seed gives a non-zero state.
20+ let z = seed >>> 0
21+ const next = () => {
22+ z = (z + 0x9e3779b9) >>> 0
23+ let t = z
24+ t = Math.imul(t ^ (t >>> 16), 0x21f0aaad)
25+ t = Math.imul(t ^ (t >>> 15), 0x735a2d97)
26+ return (t ^ (t >>> 15)) >>> 0
27+ }
28+ this.s0 = next()
29+ this.s1 = next()
30+ this.s2 = next()
31+ this.s3 = next()
32+ if ((this.s0 | this.s1 | this.s2 | this.s3) === 0) this.s0 = 1
33+ }
34+
35+ /** Uniform on [0, 1). */
36+ uniform(): number {
37+ const s1 = this.s1
38+ const x = Math.imul(s1, 5)
39+ const result = (Math.imul((x << 7) | (x >>> 25), 9) >>> 0) / 4294967296
40+ const t = s1 << 9
41+ this.s2 ^= this.s0
42+ this.s3 ^= s1
43+ this.s1 = s1 ^ this.s2
44+ this.s0 ^= this.s3
45+ this.s2 ^= t
46+ this.s3 = (this.s3 << 11) | (this.s3 >>> 21)
47+ return result
48+ }
49+
50+ /** Standard normal. */
51+ normal(): number {
52+ if (this.spare !== null) {
53+ const v = this.spare
54+ this.spare = null
55+ return v
56+ }
57+ let u = 0
58+ while (u === 0) u = this.uniform()
59+ const r = Math.sqrt(-2 * Math.log(u))
60+ const theta = 2 * Math.PI * this.uniform()
61+ this.spare = r * Math.sin(theta)
62+ return r * Math.cos(theta)
63+ }
64+}
src/worker/entropyWorker.tsadded+36−0View file
@@ -0,0 +1,36 @@
1+/**
2+ * The reference-rate estimator off the main thread: one message starts an
3+ * endless refinement loop that posts one unbiased value per completed past
4+ * (with per-rep progress in between). Stopping is worker termination — the
5+ * loop never yields, but outgoing messages still flow, and per-past values
6+ * already delivered survive on the main thread, so a later run resumes at
7+ * startPast with the same derived seeds as an uninterrupted one.
8+ */
9+import { estimateOnePast, pastSeed, REPS_PER_PAST } from '../entropy'
10+
11+export interface EntropyRequest {
12+ kernel: Float64Array
13+ sigma: number
14+ /** Conditioning window M. */
15+ past: number
16+ seed: number
17+ /** Index of the first past to compute (count already done). */
18+ startPast: number
19+}
20+
21+export type EntropyUpdate =
22+ | { type: 'progress'; pastIndex: number; repsDone: number; reps: number }
23+ | { type: 'past'; pastIndex: number; value: number }
24+
25+const post = self.postMessage as (message: EntropyUpdate) => void
26+
27+self.onmessage = (e: MessageEvent<EntropyRequest>) => {
28+ const { kernel, sigma, past, seed, startPast } = e.data
29+ for (let i = startPast; ; i++) {
30+ post({ type: 'progress', pastIndex: i, repsDone: 0, reps: REPS_PER_PAST })
31+ const value = estimateOnePast(kernel, sigma, past, pastSeed(seed, i), repsDone =>
32+ post({ type: 'progress', pastIndex: i, repsDone, reps: REPS_PER_PAST }),
33+ )
34+ post({ type: 'past', pastIndex: i, value })
35+ }
36+}