1"""Lossless codecs measured in bits/sample on an int16 trace.
3Every entry returns the full encoded size including whatever the decoder needs.
4The generic byte compressors are run on the raw int16 buffer, on the
5int16-wrapped first difference, and on the byte-split ("shuffled") buffer that
6Blosc-style pipelines use.
7"""
8import bz2
9import lzma
10import math
11import zlib
13import numpy as np
14from scipy.linalg import solve_toeplitz
15from scipy.signal import lfilter
16from scipy.special import ndtr
18COEFF_PRECISION = 15
21# ------------------------------------------------------------- preprocessing
22def as_int16(z):
23 return np.asarray(z, dtype=np.int16)
26def delta(z):
27 d = np.empty_like(z)
28 d[0] = z[0]
29 d[1:] = (z[1:].astype(np.int32) - z[:-1].astype(np.int32)).astype(np.int16)
30 return d
33def byteshuffle(z):
34 b = z.tobytes()
35 a = np.frombuffer(b, dtype=np.uint8).reshape(-1, 2)
36 return np.concatenate([a[:, 0], a[:, 1]]).tobytes()
39# ------------------------------------------------------- generic compressors
40def _zlib(b):
41 return len(zlib.compress(b, 9))
44def _lzma(b):
45 return len(lzma.compress(b, preset=9 | lzma.PRESET_EXTREME))
48def _bz2(b):
49 return len(bz2.compress(b, 9))
52def _zstd(b):
53 import zstandard
54 return len(zstandard.ZstdCompressor(level=19).compress(b))
57def _brotli(b):
58 import brotli
59 return len(brotli.compress(b, quality=11))
62def _lz4(b):
63 import lz4.frame
64 return len(lz4.frame.compress(b, compression_level=12))
67GENERIC = {
68 "zlib": _zlib, "zstd": _zstd, "lzma": _lzma,
69 "bz2": _bz2, "brotli": _brotli, "lz4": _lz4,
70}
73# ------------------------------------------------------- lossless audio codecs
74def flac_bytes(z, rate=30000):
75 """FLAC via libsndfile: LPC + Rice-coded integer residual, the mature
76 instance of the architecture LPC+ANS also belongs to. libsndfile does not
77 expose the compression level; this is its default (level 5)."""
78 import io
79 import soundfile as sf
80 buf = io.BytesIO()
81 sf.write(buf, np.asarray(z, dtype=np.int16), int(rate),
82 format="FLAC", subtype="PCM_16")
83 return buf.getbuffer().nbytes
86# --------------------------------------------------------------------- rates
87def entropy0(v):
88 counts = np.unique(np.asarray(v), return_counts=True)[1]
89 p = counts / counts.sum()
90 return float(-(p * np.log2(p)).sum())
93# ------------------------------------------------ integer LPC (the app's own)
94def autocorr(z, order, ridge=1e-8):
95 """Autocorrelation lags 0..order, with a small ridge on lag 0 so Levinson
96 stays non-singular on degenerate (near-constant, near-empty) blocks."""
97 zf = np.asarray(z, dtype=np.float64)
98 r = np.array([zf @ zf if lag == 0 else zf[lag:] @ zf[:-lag]
99 for lag in range(order + 1)])
100 r[0] = max(r[0], 1e-12 * zf.size) * (1.0 + ridge)
101 return r
104def fit_lpc_quantized(z, order):
105 r = autocorr(z, order)
106 a = solve_toeplitz(r[:order], r[1:order + 1])
107 peak = float(np.abs(a).max())
108 if peak <= 0:
109 return np.zeros(order, dtype=np.int64), 0
110 shift = COEFF_PRECISION - 1 - int(np.floor(np.log2(peak))) - 1
111 shift = max(0, min(15, shift))
112 limit = 2 ** (COEFF_PRECISION - 1)
113 q = np.clip(np.round(a * 2.0 ** shift), -limit, limit - 1).astype(np.int64)
114 return q, shift
117def lpc_residual(z, q, shift):
118 z = np.asarray(z, dtype=np.int64)
119 order = len(q)
120 acc = np.convolve(z, q, mode="full")[:len(z)]
121 pred = np.zeros_like(z)
122 pred[1:] = acc[:-1] >> shift # floor division
123 e = z.copy()
124 e[order:] = z[order:] - pred[order:]
125 return (((e + (1 << 15)) & 0xFFFF) - (1 << 15)).astype(np.int16)
128def lpc_ans_bytes(z, order):
129 """Real encoded size of integer-LPC + rANS, side information included."""
130 import simple_ans
131 q, shift = fit_lpc_quantized(z, order)
132 e = lpc_residual(z, q, shift)
133 enc = simple_ans.ans_encode(e)
134 payload = 4 * enc.words.size
135 table = 4 * enc.symbol_counts.size + 2 * enc.symbol_values.size
136 header = 2 * order + 4 + 8 # coefficients, shift, n
137 return payload + table + header, entropy0(e)
140# ------------------------------- conditional-Gaussian coding (achievable rate)
141def conditional_gaussian_rate(z, order):
142 """Cross-entropy of the real-valued-prediction conditional-Gaussian model,
143 in bits/sample, plus header cost. This is what the arithmetic coder of
144 exploration/codec_gaussian.py achieves to within ~0.1%."""
145 zf = np.asarray(z, dtype=np.float64)
146 n = zf.size
147 r = autocorr(zf, order) / n
148 a = solve_toeplitz(r[:order], r[1:order + 1])
149 pred = lfilter(np.concatenate(([0.0], a)), [1.0], zf)
150 zt, mu = zf[order:], pred[order:]
151 s0 = math.sqrt(max((zt - mu).var() - 1.0 / 12.0, 1e-6))
152 best = np.inf
153 for s in s0 * np.linspace(0.7, 1.3, 25):
154 p = ndtr((zt + 0.5 - mu) / s) - ndtr((zt - 0.5 - mu) / s)
155 best = min(best, float(-np.log2(np.maximum(p, 1e-12)).mean()))
156 header_bits = 8 * (4 * order + 14) / n
157 return best + header_bits
160# ------------------------------------------------------------------ the suite
161def measure(z, lpc_order=32, generic=True):
162 """bits/sample for every method, as an ordered dict."""
163 z = as_int16(z)
164 n = z.size
165 out = {}
166 out["raw int16"] = 16.0
167 out["order-0 H(z)"] = entropy0(z)
168 if generic:
169 buffers = {"": z.tobytes(), "+delta": delta(z).tobytes(),
170 "+shuffle": byteshuffle(z)}
171 for name, fn in GENERIC.items():
172 for suffix, buf in buffers.items():
173 out[f"{name}{suffix}"] = 8.0 * fn(buf) / n
174 try:
175 out["FLAC"] = 8.0 * flac_bytes(z) / n
176 except Exception as exc: # pragma: no cover
177 out["FLAC"] = float("nan")
178 print(f" (FLAC unavailable: {exc})")
179 nbytes, resid_h0 = lpc_ans_bytes(z, lpc_order)
180 out[f"LPC({lpc_order})+ANS"] = 8.0 * nbytes / n
181 out[f"LPC({lpc_order}) resid H0"] = resid_h0
182 out[f"cond-Gauss({lpc_order})"] = conditional_gaussian_rate(z, lpc_order)
183 return out