/ proof-of-concept / timeseries-compressibility
Sign in
proof-of-concept / timeseries-compressibility
Code against the prediction, not the integer residual
LPC + a memoryless coder plateaus far above the entropy rate whenever the prediction error is a fraction of a quantization step: the integer residual has already destroyed the fractional phase of the prediction, and whether it falls near a bin centre or a bin edge is worth 0.3-0.4 bits/sample. At sigma 5 with a 31-tap 300-2000 Hz bandpass, that is the whole distance between the 14.3x LPC+ANS plateau and the ~23x the rate allows. So add a tenth method that keeps it: the same Levinson predictor at full precision, each sample arithmetic-coded under a discretized Gaussian centred on the real-valued prediction. It reaches 21.7x on the same block, with the float32 coefficients charged and every reported size decoded and verified. Drop the hollow order-0 limit bars while here. They invited a comparison against R that does not hold -- they bound a per-sample coder on one prefiltered stream, not the process -- and the conditional coder sails past them, which read as a contradiction. The information they carried was really about coder overhead, so each coder now reports the model it codes against (ANS its histogram, the arithmetic coder its own predictive distribution) and a strip under the chart scores output against model: 1.4-1.9% for ANS, its symbol table plus arithmetic loss. That separates how well a coder does its job from how good its model was -- and both coders here lose about the same, so the 14x-to-22x gap is entirely the model.
Jeremy Magland <jmagland@flatironinstitute.org> committed commit 4ec013334674 parent 6b485df Browse files
8 changed files+576−154
README.mdmodified+25−9View file
@@ -7,14 +7,27 @@ steps) → FIR filter → round to integers. The app shows the filter (convoluti
77 kernel and frequency response, with cutoffs in Hz against a chosen sample
88 rate), a window of the generated integer signal (stationary by default, with a
99 play toggle to let it stream endlessly), and the measured compression of a
10-block of the generated integers under nine methods — zlib, zstd, and an rANS
11-entropy coder, each raw, delta-coded, and LPC-residual-coded — as bits per
12-sample and as ratio against raw int16 storage. The predictor order and the
10+block of the generated integers under ten methods — zlib, zstd, and an rANS
11+entropy coder, each raw, delta-coded, and LPC-residual-coded, plus the
12+conditional-Gaussian coder below — as bits per sample and as ratio against raw
13+int16 storage. The predictor order and the
1314 block size are controls, so the measurement can be pushed from 10k to a million
14-samples and LPC from order 1 to 128. Each prefilter group also carries a hollow
15-bar: the order-0 entropy of the stream being coded, the limit a per-sample
16-entropy coder cannot beat, which ANS misses by 1–2% (its symbol table plus its
17-own arithmetic loss).
15+samples and LPC from order 1 to 128. Under the chart, every coder that codes
16+against an explicit model is scored against it — ANS against the order-0
17+entropy of the stream it was handed, the arithmetic coder against its own
18+predictive distribution — which separates how well a coder does its job (1–2%
19+overhead for ANS, its symbol table plus its arithmetic loss) from how good the
20+model was in the first place.
21+
22+The tenth method is the one that can pass those hollow bars: the same LPC
23+prediction kept at full precision, each sample arithmetic-coded under a
24+discretized Gaussian centred on the real-valued prediction. When the
25+prediction error is a fraction of a quantization step (narrowband filters,
26+moderate σ), whether the prediction falls near a bin centre or a bin edge is
27+worth ~0.3–0.4 bits/sample — information the integer residual has already
28+destroyed, which is why LPC+ANS plateaus far above R there. Its only limit is
29+R itself. Like every other bar, its size is real: encoded, decoded, verified,
30+side information included.
1831
1932 The entropy rate R of the process — the bits/sample limit no lossless method
2033 can beat — is estimated in the browser by the method of the companion
@@ -47,8 +60,11 @@ src/entropy/ the unbiased entropy-rate estimator: hand-synced TypeScript
4760 sampler, Rhee–Glynn telescoping, Cody erfc / Acklam ndtri,
4861 xoshiro128** RNG)
4962 src/compress/ lossless codecs run in the browser: zlib (fflate), zstd (wasm),
50- ans.ts (a bit-identical port of simple_ans), and FLAC-style
51- integer LPC; borrowed from entropy-quantized-linear-transform
63+ ans.ts (a bit-identical port of simple_ans), FLAC-style
64+ integer LPC (borrowed from entropy-quantized-linear-transform),
65+ and conditionalGaussian.ts — real-coefficient prediction with
66+ each sample arithmetic-coded under a discretized Gaussian at
67+ the real-valued prediction
5268 src/worker/ the codecs and the estimator run off the main thread; the
5369 estimator worker refines one past at a time until terminated
5470 src/components/ controls, filter plots, signal canvas, compression chart,
src/App.tsxmodified+10−12View file
@@ -9,7 +9,7 @@ import { useEntropyRate } from './components/useEntropyRate'
99 import { DEFAULT_SPEC, clampSpec, designKernel, kernelNorm } from './model/filters'
1010 import { LATENT_SEED } from './model/latent'
1111 import { DEFAULT_LPC_ORDER, LPC_ORDERS } from './compress/codecs'
12-import type { BoundResult, CodecResult } from './compress/codecs'
12+import type { CodecResult } from './compress/codecs'
1313 import type { CompressRequest, CompressResponse } from './worker/compressWorker'
1414
1515 const BLOCK_SIZES = [10000, 20000, 50000, 100000, 200000, 500000, 1000000]
@@ -17,13 +17,12 @@ const DEFAULT_BLOCK_SIZE = 100000
1717
1818 interface CompressionState {
1919 results: CodecResult[]
20- bounds: BoundResult[]
2120 empiricalStd: number
2221 computing: boolean
2322 error: string | null
2423 }
2524
26-/** The nine codec sizes, measured in a worker on a debounced parameter set. */
25+/** The ten codec sizes, measured in a worker on a debounced parameter set. */
2726 function useCompression(
2827 kernel: Float64Array,
2928 sigma: number,
@@ -32,7 +31,6 @@ function useCompression(
3231 ): CompressionState {
3332 const [state, setState] = useState<CompressionState>({
3433 results: [],
35- bounds: [],
3634 empiricalStd: 0,
3735 computing: true,
3836 error: null,
@@ -48,7 +46,6 @@ function useCompression(
4846 if (e.data.id !== idRef.current) return
4947 setState({
5048 results: e.data.error ? [] : e.data.results,
51- bounds: e.data.error ? [] : e.data.bounds,
5249 empiricalStd: e.data.empiricalStd,
5350 computing: false,
5451 error: e.data.error ?? null,
@@ -239,7 +236,6 @@ export default function App() {
239236 ) : (
240237 <CompressionChart
241238 results={compression.results}
242- bounds={compression.bounds}
243239 rateBits={entropyRate.mean}
244240 computing={compression.computing}
245241 />
@@ -247,12 +243,14 @@ export default function App() {
247243 <p className="card-note">
248244 Measured on a {blockSize.toLocaleString()}-sample block of the same latent data the
249245 signal view shows; sizes include everything a decoder needs (ANS symbol table, LPC
250- coefficients). Baseline is raw int16 (16 bits/sample). The hollow bar in each group is
251- that group's entropy limit — the order-0 entropy of the stream being coded, which no
252- per-sample entropy coder can beat and ANS falls short of by its symbol table plus its
253- own arithmetic loss. The dashed line is the entropy rate R of the process itself,
254- estimated by the button in the stat row above — the limit no lossless method
255- whatsoever can beat (see the method section at the bottom).
246+ coefficients). Baseline is raw int16 (16 bits/sample). The dashed line is the entropy
247+ rate R of the process itself, estimated by the button in the stat row above — the one
248+ limit no lossless method whatsoever can beat (see the method section at the bottom).
249+ What separates the methods is the model each one codes against: ANS uses the histogram
250+ of whatever stream it is given, so a better prefilter is the only way it improves,
251+ while the conditional-Gaussian coder codes each sample against a prediction and can
252+ therefore approach R. The strip under the chart scores each coder against its own
253+ model, which is a question about the coder rather than about the model.
256254 </p>
257255 <CopyableCommand
258256 label="cross-check R from the command line:"
src/app.cssmodified+32−11View file
@@ -13,6 +13,7 @@
1313 --series-1: #2a78d6;
1414 --series-2: #eb6834;
1515 --series-3: #1baf7a;
16+ --series-4: #8a5cd6;
1617 }
1718 @media (prefers-color-scheme: dark) {
1819 :root {
@@ -28,6 +29,7 @@
2829 --series-1: #3987e5;
2930 --series-2: #d95926;
3031 --series-3: #199e70;
32+ --series-4: #9a6fe8;
3133 }
3234 }
3335
@@ -385,12 +387,6 @@ body {
385387 vertical-align: -1px;
386388 }
387389
388-/* Hollow: the mark for a limit, matching the outlined bars. */
389-.legend .swatch.hollow {
390- background: color-mix(in srgb, var(--muted) 12%, transparent);
391- border: 1.25px solid var(--muted);
392-}
393-
394390 .chart-body {
395391 position: relative;
396392 transition: opacity 0.15s;
@@ -426,11 +422,36 @@ body {
426422 font-variant-numeric: tabular-nums;
427423 }
428424
429-/* Limit rows read one step quieter than the measured ones. */
430-.bar-row-label.bound,
431-.bar-value.bound {
432- fill: var(--muted);
433- font-style: italic;
425+/* Coder overhead: a quiet strip under the chart, one entry per coder that
426+ has an explicit model to be scored against. */
427+.coder-efficiency {
428+ display: flex;
429+ flex-wrap: wrap;
430+ align-items: baseline;
431+ gap: 4px 18px;
432+ margin-top: 10px;
433+ font-size: 12px;
434+ color: var(--ink-2);
435+}
436+
437+.coder-efficiency-title {
438+ flex-basis: 100%;
439+ color: var(--muted);
440+}
441+
442+.coder-efficiency-item {
443+ display: inline-flex;
444+ align-items: baseline;
445+ gap: 6px;
446+ font-variant-numeric: tabular-nums;
447+}
448+
449+.coder-efficiency-bits {
450+ color: var(--muted);
451+}
452+
453+.coder-efficiency-pct {
454+ font-weight: 600;
434455 }
435456
436457 .axis-tick {
src/components/CompressionChart.tsxmodified+91−75View file
@@ -1,11 +1,11 @@
11 import { useRef, useState } from 'react'
2-import type { BoundResult, CodecResult } from '../compress/codecs'
2+import type { CodecResult } from '../compress/codecs'
33 import { useWidth } from './useWidth'
44
55 const GROUPS = ['no prefilter', 'delta', 'LPC']
66 const CODER_NAMES = ['zlib', 'zstd', 'ANS']
77 const CODER_VARS = ['var(--series-1)', 'var(--series-2)', 'var(--series-3)']
8-const BOUND_LABEL = 'entropy limit'
8+const COND_VAR = 'var(--series-4)'
99
1010 const LABEL_W = 100
1111 const RIGHT_PAD = 64
@@ -13,11 +13,10 @@ const AXIS_H = 26
1313 const GROUP_H = 20
1414 const ROW_H = 24
1515 const BAR_H = 16
16-const ROWS_PER_GROUP = 4
1716
1817 type Metric = 'bits' | 'ratio'
1918
20-/** One drawn bar: a measured codec size, or the entropy limit for its group. */
19+/** One drawn bar: a measured codec size. */
2120 interface Row {
2221 key: string
2322 label: string
@@ -26,7 +25,6 @@ interface Row {
2625 bytes: number
2726 bitsPerSample: number
2827 ratio: number
29- isBound: boolean
3028 color: string
3129 }
3230
@@ -49,14 +47,19 @@ function barPath(x0: number, y: number, len: number, h: number): string {
4947 return `M${x0},${y} h${len - r} a${r},${r} 0 0 1 ${r},${r} v${h - 2 * r} a${r},${r} 0 0 1 ${-r},${r} h${-(len - r)} z`
5048 }
5149
52-/** Codec rows then the entropy limit, group by group. */
53-function buildRows(results: CodecResult[], bounds: BoundResult[]): Row[] {
54- const rows: Row[] = []
55- GROUPS.forEach((group, g) => {
56- CODER_NAMES.forEach((coder, c) => {
57- const r = results[g * 3 + c]
58- if (!r) return
59- rows.push({
50+interface Group {
51+ label: string
52+ rows: Row[]
53+}
54+
55+/** The three prefilter groups, then the conditional coder — results[9] — as
56+ * its own group. */
57+function buildGroups(results: CodecResult[]): Group[] {
58+ const groups: Group[] = GROUPS.map((group, g) => ({
59+ label: group,
60+ rows: CODER_NAMES.map((coder, c) => ({ coder, c, r: results[g * 3 + c] }))
61+ .filter(x => x.r)
62+ .map(({ coder, c, r }) => ({
6063 key: r.codec,
6164 label: coder,
6265 name: r.codec,
@@ -64,31 +67,32 @@ function buildRows(results: CodecResult[], bounds: BoundResult[]): Row[] {
6467 bytes: r.bytes,
6568 bitsPerSample: r.bitsPerSample,
6669 ratio: r.ratio,
67- isBound: false,
6870 color: CODER_VARS[c],
69- })
71+ })),
72+ }))
73+ const extra = results[GROUPS.length * CODER_NAMES.length]
74+ if (extra) {
75+ groups.push({
76+ label: 'conditional',
77+ rows: [
78+ {
79+ key: extra.codec,
80+ label: 'Gaussian AC',
81+ name: extra.codec,
82+ note: extra.note,
83+ bytes: extra.bytes,
84+ bitsPerSample: extra.bitsPerSample,
85+ ratio: extra.ratio,
86+ color: COND_VAR,
87+ },
88+ ],
7089 })
71- const b = bounds.find(x => x.group === group)
72- if (b) {
73- rows.push({
74- key: `${group}-bound`,
75- label: BOUND_LABEL,
76- name: `${BOUND_LABEL} (${group})`,
77- note: b.note,
78- bytes: b.bytes,
79- bitsPerSample: b.bitsPerSample,
80- ratio: b.ratio,
81- isBound: true,
82- color: 'var(--muted)',
83- })
84- }
85- })
86- return rows
90+ }
91+ return groups
8792 }
8893
8994 export default function CompressionChart(props: {
9095 results: CodecResult[]
91- bounds: BoundResult[]
9296 /** The browser-estimated entropy rate R, once at least one past is in. */
9397 rateBits: number | null
9498 computing: boolean
@@ -99,12 +103,30 @@ export default function CompressionChart(props: {
99103 const [tip, setTip] = useState<Tip | null>(null)
100104 const [hovered, setHovered] = useState<string | null>(null)
101105
102- const { results, bounds, rateBits } = props
106+ const { results, rateBits } = props
103107 if (results.length === 0) {
104108 return <p className="card-note">Computing compression on the first block…</p>
105109 }
106110
107- const rows = buildRows(results, bounds)
111+ const groups = buildGroups(results)
112+ // Every coder that codes against an explicit model can be scored against it.
113+ const efficiency = results.filter(r => r.modelBitsPerSample !== undefined)
114+
115+ // Group headers and bar rows laid out top to bottom; groups may differ in
116+ // row count, so positions accumulate rather than being indexed.
117+ const groupLabels: { label: string; y: number }[] = []
118+ const placed: { row: Row; y: number }[] = []
119+ let yCursor = AXIS_H
120+ for (const g of groups) {
121+ groupLabels.push({ label: g.label, y: yCursor + 15 })
122+ yCursor += GROUP_H
123+ for (const row of g.rows) {
124+ placed.push({ row, y: yCursor })
125+ yCursor += ROW_H
126+ }
127+ }
128+
129+ const rows = placed.map(p => p.row)
108130 const value = (r: { bitsPerSample: number; ratio: number }) =>
109131 metric === 'bits' ? r.bitsPerSample : r.ratio
110132 const rateValue =
@@ -115,13 +137,8 @@ export default function CompressionChart(props: {
115137 : Math.max(...rows.map(value), rateValue ?? 0) * 1.1
116138
117139 const plotW = width - LABEL_W - RIGHT_PAD
118- const height = AXIS_H + GROUPS.length * (GROUP_H + ROWS_PER_GROUP * ROW_H) + 6
140+ const height = yCursor + 6
119141 const xOf = (v: number) => LABEL_W + (v / xMax) * plotW
120- const rowY = (i: number) =>
121- AXIS_H +
122- Math.floor(i / ROWS_PER_GROUP) * (GROUP_H + ROWS_PER_GROUP * ROW_H) +
123- GROUP_H +
124- (i % ROWS_PER_GROUP) * ROW_H
125142
126143 const fmt = (r: Row) => (metric === 'bits' ? r.bitsPerSample.toFixed(2) : `${r.ratio.toFixed(2)}×`)
127144
@@ -162,8 +179,8 @@ export default function CompressionChart(props: {
162179 </span>
163180 ))}
164181 <span>
165- <span className="swatch hollow" />
166- entropy limit (not achieved)
182+ <span className="swatch" style={{ background: COND_VAR }} />
183+ cond. Gaussian AC
167184 </span>
168185 </div>
169186 </div>
@@ -179,45 +196,20 @@ export default function CompressionChart(props: {
179196 </g>
180197 ))}
181198 <line x1={xOf(0)} x2={xOf(0)} y1={AXIS_H - 6} y2={height - 4} stroke="var(--baseline)" strokeWidth={1} />
182- {GROUPS.map((g, gi) => (
183- <text
184- key={g}
185- x={0}
186- y={AXIS_H + gi * (GROUP_H + ROWS_PER_GROUP * ROW_H) + 15}
187- className="bar-group-label"
188- >
189- {g}
199+ {groupLabels.map(g => (
200+ <text key={g.label} x={0} y={g.y} className="bar-group-label">
201+ {g.label}
190202 </text>
191203 ))}
192- {rows.map((r, i) => {
193- const y = rowY(i)
204+ {placed.map(({ row: r, y }) => {
194205 const len = Math.max(1, (value(r) / xMax) * plotW)
195206 return (
196207 <g key={r.key} opacity={hovered === null || hovered === r.key ? 1 : 0.45}>
197- <text
198- x={8}
199- y={y + BAR_H / 2 + 4}
200- className={r.isBound ? 'bar-row-label bound' : 'bar-row-label'}
201- >
208+ <text x={8} y={y + BAR_H / 2 + 4} className="bar-row-label">
202209 {r.label}
203210 </text>
204- {r.isBound ? (
205- // Hollow: a limit nobody reached, not a measured size.
206- <path
207- d={barPath(xOf(0), y + 1, len, BAR_H - 2)}
208- fill="var(--muted)"
209- fillOpacity={0.12}
210- stroke="var(--muted)"
211- strokeWidth={1.25}
212- />
213- ) : (
214- <path d={barPath(xOf(0), y, len, BAR_H)} fill={r.color} />
215- )}
216- <text
217- x={xOf(0) + len + 6}
218- y={y + BAR_H / 2 + 4}
219- className={r.isBound ? 'bar-value bound' : 'bar-value'}
220- >
211+ <path d={barPath(xOf(0), y, len, BAR_H)} fill={r.color} />
212+ <text x={xOf(0) + len + 6} y={y + BAR_H / 2 + 4} className="bar-value">
221213 {fmt(r)}
222214 </text>
223215 <rect
@@ -270,12 +262,36 @@ export default function CompressionChart(props: {
270262 <span className="tip-label">{tip.row.name}</span>
271263 </div>
272264 <div className="tip-label">
273- {tip.row.isBound ? 'equivalent to ' : ''}
274265 {tip.row.bytes.toLocaleString()} bytes · {tip.row.note}
275266 </div>
276267 </div>
277268 )}
278269 </div>
270+ {/* How much each entropy coder loses against the model it is coding
271+ against — its own overhead, separate from how good the model is. */}
272+ {efficiency.length > 0 && (
273+ <div className="coder-efficiency">
274+ <span className="coder-efficiency-title">
275+ entropy-coder overhead — output vs the model it codes against
276+ </span>
277+ {efficiency.map(r => {
278+ const model = r.modelBitsPerSample as number
279+ const over = model > 0 ? (r.bitsPerSample / model - 1) * 100 : 0
280+ return (
281+ <span key={r.codec} className="coder-efficiency-item">
282+ <span className="coder-efficiency-name">{r.codec}</span>
283+ <span className="coder-efficiency-bits">
284+ {r.bitsPerSample.toFixed(3)} / {model.toFixed(3)} bits
285+ </span>
286+ <span className="coder-efficiency-pct">
287+ {over >= 0 ? '+' : ''}
288+ {over.toFixed(1)}%
289+ </span>
290+ </span>
291+ )
292+ })}
293+ </div>
294+ )}
279295 <details className="chart-table">
280296 <summary>Table view</summary>
281297 <table>
src/compress/codecs.tsmodified+25−39View file
@@ -14,12 +14,19 @@ import { fitLpc, lpcResidual, lpcRestore, modelSize } from './lpc'
1414 // cannot be resolved.
1515 import zstdWasmUrl from '../../node_modules/@bokuweb/zstd-wasm/dist/web/zstd.wasm?url'
1616
17+/**
18+ * Compressed size in bytes, including any table the decoder needs. Codecs
19+ * that code against an explicit probability model also report that model's
20+ * ideal cost in bits, so the coder's own overhead — its tables, its
21+ * coefficients, and its arithmetic loss — can be shown against it.
22+ */
23+export type CodecSize = number | { bytes: number; modelBits: number }
24+
1725 export interface Codec {
1826 name: string
1927 /** Longer description, shown on hover. */
2028 note: string
21- /** Compressed size in bytes, including any table the decoder needs. */
22- size: (samples: Int16Array, bytes: Uint8Array) => number
29+ size: (samples: Int16Array, bytes: Uint8Array) => CodecSize
2330 }
2431
2532 /**
@@ -54,20 +61,26 @@ function asBytes(samples: Int16Array): Uint8Array {
5461 * original samples — so a reported size always belongs to an encoding that
5562 * actually round-trips. `extraBytes` covers anything else the decoder needs,
5663 * such as predictor coefficients.
64+ *
65+ * ANS codes each symbol against the histogram of `coded`, so that histogram's
66+ * order-0 entropy is exactly the model this coder is aiming at.
5767 */
5868 function ansSize(
5969 samples: Int16Array,
6070 coded: Int16Array,
6171 inverse: (coded: Int16Array) => Int16Array,
6272 extraBytes = 0,
63-): number {
73+): CodecSize {
6474 const encoded = ansEncode(coded)
6575 const decoded = inverse(ansDecode(encoded))
6676 if (decoded.length !== samples.length) throw new Error('ANS round-trip length mismatch')
6777 for (let i = 0; i < samples.length; i++) {
6878 if (decoded[i] !== samples[i]) throw new Error(`ANS round-trip mismatch at ${i}`)
6979 }
70- return encodedSize(encoded) + extraBytes
80+ return {
81+ bytes: encodedSize(encoded) + extraBytes,
82+ modelBits: order0Entropy(coded) * coded.length,
83+ }
7184 }
7285
7386 /** Predictor orders the UI offers; 32 is the FLAC default and ours. */
@@ -182,8 +195,7 @@ export const GENERAL_CODECS: Codec[] = [ZLIB, ZSTD]
182195 /**
183196 * Order-0 (memoryless) entropy of an int16 stream, in bits per sample: what a
184197 * perfect entropy coder for the sample histogram would spend, with nothing
185- * charged for describing that histogram. The ANS bars sit above this by the
186- * symbol table plus the coder's own arithmetic loss.
198+ * charged for describing that histogram — the model an ANS coder aims at.
187199 */
188200 export function order0Entropy(samples: Int16Array): number {
189201 const counts = new Int32Array(65536)
@@ -198,38 +210,6 @@ export function order0Entropy(samples: Int16Array): number {
198210 return bits
199211 }
200212
201-export interface BoundResult {
202- /** The prefilter group this bounds, matching the codec groups. */
203- group: string
204- note: string
205- bitsPerSample: number
206- ratio: number
207- bytes: number
208-}
209-
210-/** The order-0 bound for each prefilter: raw samples, delta, LPC residual. */
211-export function entropyBounds(samples: Int16Array, lpcOrder: number): BoundResult[] {
212- const streams: { group: string; what: string; data: Int16Array }[] = [
213- { group: 'no prefilter', what: 'the samples themselves', data: samples },
214- { group: 'delta', what: 'the first differences', data: delta(samples) },
215- {
216- group: 'LPC',
217- what: `the order-${lpcOrder} prediction residual`,
218- data: lpcTransform(samples, lpcOrder).residual,
219- },
220- ]
221- return streams.map(({ group, what, data }) => {
222- const bitsPerSample = order0Entropy(data)
223- return {
224- group,
225- note: `order-0 entropy of ${what} — the limit for a per-sample entropy coder, with no symbol table or coefficients charged`,
226- bitsPerSample,
227- ratio: 16 / bitsPerSample,
228- bytes: Math.ceil((bitsPerSample * samples.length) / 8),
229- }
230- })
231-}
232-
233213 let ready: Promise<void> | null = null
234214
235215 // zstd-wasm publishes the *node* build's types while the bundler resolves the
@@ -254,6 +234,9 @@ export interface CodecResult {
254234 bytes: number
255235 ratio: number
256236 bitsPerSample: number
237+ /** Ideal cost of this coder's own probability model, when it has one — the
238+ * yardstick its output is measured against. Absent for zlib and zstd. */
239+ modelBitsPerSample?: number
257240 }
258241
259242 /**
@@ -264,13 +247,16 @@ export interface CodecResult {
264247 export function compressAll(samples: Int16Array, codecs: Codec[]): CodecResult[] {
265248 const bytes = asBytes(samples)
266249 return codecs.map(codec => {
267- const size = codec.size(samples, bytes)
250+ const out = codec.size(samples, bytes)
251+ const size = typeof out === 'number' ? out : out.bytes
268252 return {
269253 codec: codec.name,
270254 note: codec.note,
271255 bytes: size,
272256 ratio: bytes.byteLength / size,
273257 bitsPerSample: (8 * size) / samples.length,
258+ modelBitsPerSample:
259+ typeof out === 'number' ? undefined : out.modelBits / samples.length,
274260 }
275261 })
276262 }
src/compress/conditionalGaussian.tsadded+380−0View file
@@ -0,0 +1,380 @@
1+/**
2+ * Conditional-Gaussian arithmetic coding: the prefilter-free way past the
3+ * order-0 entropy limit.
4+ *
5+ * LPC + a memoryless coder codes the *integer* residual with one pooled
6+ * histogram — the mixture over all fractional parts of the prediction. When
7+ * the prediction error is a fraction of a quantization step, that mixture
8+ * costs ~0.3-0.4 bits/sample more than the conditional entropy: whether the
9+ * real-valued prediction μ falls near a bin centre (~0.4 bits) or a bin edge
10+ * (~1.1 bits) is information an integer residual has already destroyed.
11+ *
12+ * This coder keeps it: the same order-p predictor as the LPC codecs, but with
13+ * float coefficients, and each sample is arithmetic-coded under the
14+ * discretized Gaussian N(μ_t, s²) — bin k gets Φ((k+½−μ)/s) − Φ((k−½−μ)/s),
15+ * tails folded into the edge bins, one escape symbol so that any int16 input
16+ * still round-trips. The decoder rebuilds μ_t from already-decoded samples
17+ * with the identical arithmetic, so both sides derive bit-identical tables.
18+ *
19+ * The header charges everything the decoder needs: the order, the float32
20+ * coefficients, the residual scale s, and the marginal std for the warm-up
21+ * samples (coded before p samples of context exist).
22+ */
23+import { ndtr } from '../entropy/normal'
24+import { fitRealCoeffs } from './lpc'
25+import type { Codec, CodecSize } from './codecs'
26+
27+const TOTAL = 1 << 16
28+const FULL = 2 ** 32 - 1
29+const HALF = 2 ** 31
30+const QUARTER = 2 ** 30
31+
32+/** Symbols cover residuals in [-kwin, kwin); beyond that the escape fires. */
33+const MAX_KWIN = 8192
34+
35+class BitWriter {
36+ bytes: number[] = []
37+ private acc = 0
38+ private n = 0
39+
40+ write(bit: number) {
41+ this.acc = this.acc * 2 + bit
42+ if (++this.n === 8) {
43+ this.bytes.push(this.acc)
44+ this.acc = 0
45+ this.n = 0
46+ }
47+ }
48+
49+ flush() {
50+ while (this.n) this.write(0)
51+ }
52+}
53+
54+class BitReader {
55+ private pos = 0
56+ constructor(private data: number[]) {}
57+
58+ read(): number {
59+ const byte = this.pos >> 3 < this.data.length ? this.data[this.pos >> 3] : 0
60+ const bit = (byte >> (7 - (this.pos & 7))) & 1
61+ this.pos++
62+ return bit
63+ }
64+}
65+
66+/** Witten–Neal–Cleary arithmetic coder, 32-bit width. All products stay
67+ * below 2^48, exact in a double, so no BigInt is needed. */
68+class ArithEncoder {
69+ private low = 0
70+ private high = FULL
71+ private pending = 0
72+ readonly out = new BitWriter()
73+
74+ private emit(bit: number) {
75+ this.out.write(bit)
76+ for (; this.pending > 0; this.pending--) this.out.write(1 - bit)
77+ }
78+
79+ encode(cum: number, freq: number, tot: number) {
80+ const span = this.high - this.low + 1
81+ this.high = this.low + Math.floor((span * (cum + freq)) / tot) - 1
82+ this.low = this.low + Math.floor((span * cum) / tot)
83+ for (;;) {
84+ if (this.high < HALF) {
85+ this.emit(0)
86+ } else if (this.low >= HALF) {
87+ this.emit(1)
88+ this.low -= HALF
89+ this.high -= HALF
90+ } else if (this.low >= QUARTER && this.high < HALF + QUARTER) {
91+ this.pending++
92+ this.low -= QUARTER
93+ this.high -= QUARTER
94+ } else {
95+ break
96+ }
97+ this.low *= 2
98+ this.high = this.high * 2 + 1
99+ }
100+ }
101+
102+ finish(): number[] {
103+ this.pending++
104+ this.emit(this.low < QUARTER ? 0 : 1)
105+ this.out.flush()
106+ return this.out.bytes
107+ }
108+}
109+
110+class ArithDecoder {
111+ private low = 0
112+ private high = FULL
113+ private value = 0
114+ private inp: BitReader
115+
116+ constructor(data: number[]) {
117+ this.inp = new BitReader(data)
118+ for (let i = 0; i < 32; i++) this.value = this.value * 2 + this.inp.read()
119+ }
120+
121+ target(tot: number): number {
122+ const span = this.high - this.low + 1
123+ return Math.floor(((this.value - this.low + 1) * tot - 1) / span)
124+ }
125+
126+ consume(cum: number, freq: number, tot: number) {
127+ const span = this.high - this.low + 1
128+ this.high = this.low + Math.floor((span * (cum + freq)) / tot) - 1
129+ this.low = this.low + Math.floor((span * cum) / tot)
130+ for (;;) {
131+ if (this.high < HALF) {
132+ // nothing to subtract
133+ } else if (this.low >= HALF) {
134+ this.low -= HALF
135+ this.high -= HALF
136+ this.value -= HALF
137+ } else if (this.low >= QUARTER && this.high < HALF + QUARTER) {
138+ this.low -= QUARTER
139+ this.high -= QUARTER
140+ this.value -= QUARTER
141+ } else {
142+ break
143+ }
144+ this.low *= 2
145+ this.high = this.high * 2 + 1
146+ this.value = this.value * 2 + this.inp.read()
147+ }
148+ }
149+}
150+
151+interface Table {
152+ kwin: number
153+ freqs: Int32Array // 2*kwin residual symbols, then the escape at index 2*kwin
154+ cum: Int32Array
155+}
156+
157+/** Frequencies (sum TOTAL) of the discretized N(d, s²) over residual bins
158+ * [-kwin, kwin), tails folded into the edge bins, every bin ≥ 1, escape = 1.
159+ * Deterministic in (d, s): encoder and decoder call it with identical
160+ * doubles, so the tables agree bit for bit. */
161+function freqTable(d: number, s: number): Table {
162+ const kwin = Math.max(8, Math.min(MAX_KWIN, Math.ceil(10 * s) + 4))
163+ const nsym = 2 * kwin + 1
164+ const freqs = new Int32Array(nsym).fill(1)
165+ const spare = TOTAL - nsym
166+ const m = Math.min(kwin - 1, Math.floor(8 * s) + 2)
167+
168+ // Bin probabilities for k in [-m, m]; the ends absorb their whole tails.
169+ const probs: number[] = []
170+ let loCdf = 0
171+ for (let k = -m; k <= m; k++) {
172+ const hiCdf = k === m ? 1 : ndtr((k + 0.5 - d) / s)
173+ probs.push(hiCdf - loCdf)
174+ loCdf = hiCdf
175+ }
176+ let assigned = 0
177+ const scaled = probs.map(p => {
178+ const f = Math.floor(p * spare)
179+ assigned += f
180+ return f
181+ })
182+ // Give the rounding deficit to the biggest bin.
183+ let imax = 0
184+ for (let i = 1; i < scaled.length; i++) if (scaled[i] > scaled[imax]) imax = i
185+ scaled[imax] += spare - assigned
186+ for (let i = 0; i < scaled.length; i++) freqs[i - m + kwin] += scaled[i]
187+
188+ const cum = new Int32Array(nsym + 1)
189+ for (let i = 0; i < nsym; i++) cum[i + 1] = cum[i] + freqs[i]
190+ return { kwin, freqs, cum }
191+}
192+
193+/**
194+ * Tables keyed by the quantized fractional phase d = μ − round(μ). The
195+ * quantization step grows with s (the phase matters less the wider the
196+ * bump: the KL cost is ~(Δd/s)²), which caps the cache at ~1k tables of
197+ * O(s)-sized windows. Encoder and decoder quantize identically, so the
198+ * cache is pure speed, not a protocol difference.
199+ */
200+function makeTableCache(s: number): (d: number) => Table {
201+ const dstep = Math.max(1 / 1024, s / 256)
202+ const cache = new Map<number, Table>()
203+ return d => {
204+ const key = Math.round(d / dstep)
205+ let table = cache.get(key)
206+ if (!table) {
207+ table = freqTable(key * dstep, s)
208+ cache.set(key, table)
209+ }
210+ return table
211+ }
212+}
213+
214+/** μ_t from the p previous samples — the one shared piece of float
215+ * arithmetic both sides must reproduce exactly, hence one function. */
216+function predict(coeffs: Float64Array, zf: Float64Array, t: number): number {
217+ let mu = 0
218+ for (let k = 0; k < coeffs.length; k++) mu += coeffs[k] * zf[t - 1 - k]
219+ return Number.isFinite(mu) ? mu : 0
220+}
221+
222+interface Model {
223+ coeffs: Float64Array // float32-rounded values, as the decoder will see them
224+ s: number
225+ stdZ: number
226+}
227+
228+/** Header: order (2) + float32 coefficients + s (4) + warm-up std (4). */
229+function headerBytes(model: Model): number {
230+ return 2 + 4 * model.coeffs.length + 4 + 4
231+}
232+
233+function fitModel(samples: Int16Array, order: number): Model {
234+ const a = fitRealCoeffs(samples, order)
235+ const coeffs = Float64Array.from(a ?? [], Math.fround)
236+
237+ let sum = 0
238+ let sumSq = 0
239+ for (let i = 0; i < samples.length; i++) {
240+ sum += samples[i]
241+ sumSq += samples[i] * samples[i]
242+ }
243+ const n = samples.length
244+ const varZ = Math.max(0, sumSq / n - (sum / n) ** 2)
245+
246+ // Residual scale: the variance of the residuals, but with outliers gated
247+ // at 8× a median-based scale first. On clean data nothing reaches 8σ and
248+ // this is exactly the residual variance; a lone spike would otherwise
249+ // inflate the variance — and with it every table's width — even though the
250+ // escape symbol already prices outliers individually.
251+ const zf = Float64Array.from(samples)
252+ const first = coeffs.length
253+ const total = Math.max(n - first, 0)
254+ const stride = Math.max(1, Math.ceil(total / (1 << 18)))
255+ const abs: number[] = []
256+ for (let t = first; t < n; t += stride) {
257+ abs.push(Math.abs(zf[t] - predict(coeffs, zf, t)))
258+ }
259+ abs.sort((x, y) => x - y)
260+ const medE = abs.length > 0 ? abs[abs.length >> 1] : Math.sqrt(varZ)
261+ const gate = Math.max(8 * 1.4826 * medE, 1)
262+ let errSq = 0
263+ let count = 0
264+ for (let t = first; t < n; t++) {
265+ const e = zf[t] - predict(coeffs, zf, t)
266+ if (Math.abs(e) <= gate) {
267+ errSq += e * e
268+ count++
269+ }
270+ }
271+ const varE = count > 0 ? errSq / count : varZ
272+ // The coding distribution round(N(μ, s²)) has variance ≈ s² + 1/12, so the
273+ // measured residual variance overshoots s² by the rounding term.
274+ const s = Math.fround(Math.max(Math.sqrt(Math.max(varE - 1 / 12, 0)), 0.02))
275+ const stdZ = Math.fround(Math.max(Math.sqrt(varZ), 0.05))
276+ return { coeffs, s, stdZ }
277+}
278+
279+const UNIFORM_TOT = 256
280+
281+/** The payload, and the ideal cost of the model it was coded against — the
282+ * sum of −log2 p over the symbols actually emitted. The gap between that and
283+ * the payload is the arithmetic coder's own rounding loss. */
284+function encodeAll(samples: Int16Array, model: Model): { payload: number[]; modelBits: number } {
285+ const { coeffs, s, stdZ } = model
286+ const order = coeffs.length
287+ const zf = Float64Array.from(samples)
288+ const enc = new ArithEncoder()
289+ const warm = freqTable(0, stdZ)
290+ const tableFor = makeTableCache(s)
291+ let modelBits = 0
292+ for (let t = 0; t < samples.length; t++) {
293+ const mu = t < order ? 0 : predict(coeffs, zf, t)
294+ const c = t < order ? 0 : Math.round(mu)
295+ const table = t < order ? warm : tableFor(mu - c)
296+ const sym = samples[t] - c + table.kwin
297+ if (sym >= 0 && sym < 2 * table.kwin) {
298+ enc.encode(table.cum[sym], table.freqs[sym], TOTAL)
299+ modelBits -= Math.log2(table.freqs[sym] / TOTAL)
300+ } else {
301+ const esc = 2 * table.kwin
302+ enc.encode(table.cum[esc], table.freqs[esc], TOTAL)
303+ const raw = samples[t] & 0xffff
304+ enc.encode(raw >> 8, 1, UNIFORM_TOT)
305+ enc.encode(raw & 0xff, 1, UNIFORM_TOT)
306+ modelBits += -Math.log2(table.freqs[esc] / TOTAL) + 16 // escape, then the raw sample
307+ }
308+ }
309+ return { payload: enc.finish(), modelBits }
310+}
311+
312+function decodeAll(payload: number[], n: number, model: Model): Int16Array {
313+ const { coeffs, s, stdZ } = model
314+ const order = coeffs.length
315+ const zf = new Float64Array(n)
316+ const out = new Int16Array(n)
317+ const dec = new ArithDecoder(payload)
318+ const warm = freqTable(0, stdZ)
319+ const tableFor = makeTableCache(s)
320+ for (let t = 0; t < n; t++) {
321+ const mu = t < order ? 0 : predict(coeffs, zf, t)
322+ const c = t < order ? 0 : Math.round(mu)
323+ const table = t < order ? warm : tableFor(mu - c)
324+ const tgt = dec.target(TOTAL)
325+ let lo = 0
326+ let hi = table.cum.length - 1
327+ while (hi - lo > 1) {
328+ const mid = (lo + hi) >> 1
329+ if (table.cum[mid] <= tgt) lo = mid
330+ else hi = mid
331+ }
332+ dec.consume(table.cum[lo], table.freqs[lo], TOTAL)
333+ let z: number
334+ if (lo < 2 * table.kwin) {
335+ z = lo - table.kwin + c
336+ } else {
337+ const hiByte = dec.target(UNIFORM_TOT)
338+ dec.consume(hiByte, 1, UNIFORM_TOT)
339+ const loByte = dec.target(UNIFORM_TOT)
340+ dec.consume(loByte, 1, UNIFORM_TOT)
341+ z = (((hiByte << 8) | loByte) << 16) >> 16
342+ }
343+ out[t] = z
344+ zf[t] = out[t]
345+ }
346+ return out
347+}
348+
349+/**
350+ * Compressed size in bytes — header plus arithmetic-coded payload — after
351+ * decoding the payload and checking it reproduces the samples exactly, so a
352+ * reported size always belongs to an encoding that round-trips.
353+ */
354+export function conditionalGaussianSize(samples: Int16Array, order: number): CodecSize {
355+ const model = fitModel(samples, order)
356+ const { payload, modelBits } = encodeAll(samples, model)
357+ const decoded = decodeAll(payload, samples.length, model)
358+ if (decoded.length !== samples.length) {
359+ throw new Error('conditional-Gaussian round-trip length mismatch')
360+ }
361+ for (let i = 0; i < samples.length; i++) {
362+ if (decoded[i] !== samples[i]) {
363+ throw new Error(`conditional-Gaussian round-trip mismatch at ${i}`)
364+ }
365+ }
366+ return { bytes: payload.length + headerBytes(model), modelBits }
367+}
368+
369+/** The codec, shaped like the others so compressAll can report it. */
370+export function conditionalGaussianCodec(order: number): Codec {
371+ return {
372+ name: `LPC(${order}) + cond. Gaussian AC`,
373+ note:
374+ `Order-${order} prediction kept at full precision; each sample arithmetic-coded ` +
375+ 'under a discretized Gaussian centred on the real-valued prediction, using the ' +
376+ 'fractional phase that integer residuals discard; the size includes the float32 ' +
377+ 'coefficients',
378+ size: samples => conditionalGaussianSize(samples, order),
379+ }
380+}
src/compress/lpc.tsmodified+5−0View file
@@ -87,6 +87,11 @@ export function fitLpc(x: Int16Array, order: number): LpcModel | null {
8787 return a ? quantiseCoefficients(a) : null
8888 }
8989
90+/** The unquantised order-p predictor, for coders that keep full precision. */
91+export function fitRealCoeffs(x: Int16Array, order: number): Float64Array | null {
92+ return levinson(autocorrelation(x, order), order)
93+}
94+
9095 /**
9196 * Prediction residual. The first `order` entries are the samples themselves
9297 * (the warm-up the decoder needs before it can predict). Differences are taken
src/worker/compressWorker.tsmodified+8−8View file
@@ -1,14 +1,13 @@
11 /**
22 * Compression measurements off the main thread, so the scrolling view never
33 * stutters while zstd -19 or the LPC fit runs. One message in (the model),
4- * one message out (the nine codec results).
4+ * one message out (the ten codec results).
55 */
66 import { LatentSource } from '../model/latent'
7+import { conditionalGaussianCodec } from '../compress/conditionalGaussian'
78 import {
89 initCodecs,
910 compressAll,
10- entropyBounds,
11- type BoundResult,
1211 ZLIB,
1312 ZSTD,
1413 ANS,
@@ -31,8 +30,6 @@ export interface CompressRequest {
3130 export interface CompressResponse {
3231 id: number
3332 results: CodecResult[]
34- /** Order-0 entropy limit for each prefilter group. */
35- bounds: BoundResult[]
3633 /** Empirical std of the quantized block, for display sanity. */
3734 empiricalStd: number
3835 error?: string
@@ -57,11 +54,14 @@ self.onmessage = async (e: MessageEvent<CompressRequest>) => {
5754 const empiricalStd = Math.sqrt(Math.max(0, sumSq / samples.length - mean * mean))
5855 post({
5956 id,
60- results: compressAll(samples, [...PLAIN_CODECS, ...lpcCodecs(lpcOrder)]),
61- bounds: entropyBounds(samples, lpcOrder),
57+ results: compressAll(samples, [
58+ ...PLAIN_CODECS,
59+ ...lpcCodecs(lpcOrder),
60+ conditionalGaussianCodec(lpcOrder),
61+ ]),
6262 empiricalStd,
6363 })
6464 } catch (err) {
65- post({ id, results: [], bounds: [], empiricalStd: 0, error: String(err) })
65+ post({ id, results: [], empiricalStd: 0, error: String(err) })
6666 }
6767 }
↑↓ move↵openescclose