From scalars to tables#
GGDLPC.Codec(blob) derives everything the coder uses from the fitted scalars: the statistic thresholds that select a cell of the scale grid, one folded GGD probability mass function per cell, a canonical Huffman table per cell under the derived quotient split, and the zero-run mode. This chapter regenerates each table inline and checks it against the package’s build_channel output.
Parameters#
N = 128 # alphabet bound: 8-bit RGB planes centered to [-128, 127]
RANK = 2 # spatial rank (image)
SEED = 20260820
FIT_IMAGES = 16 # LSDIR_8k training images used for the fit (Kodak is test-only)
FIT_CROP = 256 # center-crop side for the fitting images
KODAK_INDEX = 22 # the test image: the same Kodak image as the FRAPPE v3 notebooks
Data and fit#
Identical to the previous chapter (same images, same passes).
import numpy as np
from datasets import load_dataset
def center_crop(img, side):
w, h = img.size
left, top = (w - side) // 2, (h - side) // 2
return img.crop((left, top, left + side, top + side))
lsdir = load_dataset('danjacobellis/LSDIR_8k', split='train')
fit_set = []
for sample in lsdir:
img = sample['image'].convert('RGB')
if min(img.size) < FIT_CROP:
continue
arr = np.asarray(center_crop(img, FIT_CROP), np.int64) - 128
fit_set.append(arr.transpose(2, 0, 1)) # (3, H, W) in [-128, 127]
if len(fit_set) == FIT_IMAGES:
break
kodak = load_dataset('danjacobellis/kodak', split='validation')
test_img = kodak[KODAK_INDEX]['image'].convert('RGB')
test_tensor = np.asarray(test_img, np.int64).transpose(2, 0, 1) - 128
import torch
import GGDLPC
def loader_factory():
batches = [np.stack(fit_set[i:i + 4]) for i in range(0, len(fit_set), 4)]
return (torch.from_numpy(b) for b in batches)
channels, provenance = GGDLPC.fit(loader_factory, N=N, rank=RANK)
blob = GGDLPC.new_blob(N, RANK, channels, provenance)
GGDLPC fit N=128 rank=2: predictor moments over 4 batches (0s)
GGDLPC fit N=128 rank=2: dense h(S, r) over 4 batches (0s)
GGDLPC fit N=128 rank=2: ch 0: mu=+0.003 a=0.8852 b=0.4220 u=0.692 v=+0.080 wq=[179, 200, -105, -7, -14] bq=-35
GGDLPC fit N=128 rank=2: ch 1: mu=-0.079 a=0.8896 b=0.4631 u=0.818 v=+0.047 wq=[173, 196, -97, -7, -14] bq=-55
GGDLPC fit N=128 rank=2: ch 2: mu=-0.017 a=0.8912 b=0.3859 u=0.628 v=+0.093 wq=[175, 199, -95, -9, -17] bq=-70
builds = [GGDLPC.build_channel(ch, N) for ch in blob['channels']]
G = GGDLPC.n_grid(N)
nsym = 2 * N + 1
print(f'G = {G} grid cells, statistic cap = {GGDLPC.stat_cap(N)}, '
f'{nsym} folded symbols per cell')
G = 49 grid cells, statistic cap = 4095, 257 folded symbols per cell
The scale grid and derived edges#
The scale grid is geometric: cell \(j\) covers residual scale \(\hat m = m_\min \tau^j\) with \(\tau = 2^{1/4}\) and \(m_\min = 1/16\). Inverting the scale law maps the cell boundaries back into integer statistic thresholds: cell \(j\) is selected when \(S\) falls between edges[j-1] and edges[j]. The thresholds are forced non-negative and non-decreasing.
TAU = 2 ** 0.25
M_MIN = 1 / 16
ch0 = blob['channels'][0]
jgrid = np.arange(G - 1)
bound = M_MIN * TAU ** (jgrid + 0.5) # scale at each cell boundary
edges = np.ceil(16.0 * (bound - ch0['b']) / max(ch0['a'], 1e-9)).astype(np.int64)
edges = np.maximum.accumulate(np.maximum(edges, 0))
assert np.array_equal(edges, builds[0].edges)
import matplotlib.pyplot as plt
plt.figure(figsize=(6, 3.5))
plt.semilogy(np.arange(G - 1), np.maximum(edges, 0.5), '.-')
plt.xlabel('cell j')
plt.ylabel('statistic threshold edges[j]')
plt.tight_layout()
print('first 16 edges:', edges[:16].tolist())
first 16 edges: [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 5, 7, 9]
Per-cell conditionals: the folded GGD#
For cell \(j\), the shape law gives \(\beta\) from the cell scale, and the mean-absolute-deviation scale \(\hat m\) converts to the GGD scale \(\alpha\) through \(c(\beta) = \Gamma(1/\beta)/\Gamma(2/\beta)\). The analytic density is integrated over unit bins across three fold periods, folded modulo \(2N + 1\) (the residual \(r \in [-2N, 2N]\) reduced to \(w = (r + N) \bmod (2N+1)\)), floored at \(10^{-12}\), and normalized.
from scipy.special import gammainc, gammaln
def ggd_cdf(x, beta):
return 0.5 + 0.5 * np.sign(x) * gammainc(1.0 / beta, np.abs(x) ** beta + 1e-300)
def folded_pmf(mu, beta, alpha, N):
nsym = 2 * N + 1
half = (3 * nsym - 1) // 2
r = np.arange(-half, half + 1)
p = np.maximum(ggd_cdf((r + 0.5 - mu) / alpha, beta)
- ggd_cdf((r - 0.5 - mu) / alpha, beta), 0.0)
q = np.zeros(nsym)
np.add.at(q, (r + N) % nsym, p)
q = np.maximum(q, 1e-12)
return q / q.sum()
def cell_pmf(ch, j, N):
m = M_MIN * TAU ** j
beta = float(np.clip(ch['u'] + ch['v'] * np.log2(m), 0.2, 4.0))
alpha = np.exp(gammaln(1 / beta) - gammaln(2 / beta)) * m
return folded_pmf(ch['mu'], beta, alpha, N)
pmfs = [cell_pmf(ch0, j, N) for j in range(G)]
plt.figure(figsize=(7, 4))
r_axis = np.arange(-N, N + 1)
for j in (6, 16, 26, 36, 46):
plt.semilogy(r_axis, pmfs[j], label=f'cell {j}')
plt.ylim(1e-9, 1)
plt.xlabel('residual r')
plt.ylabel('q(r)')
plt.legend()
plt.tight_layout()
The derived quotient split#
A single Huffman table over all 257 folded symbols would be large, and most of its codewords would cover negligible mass. Each folded symbol is therefore split into \(k\) raw low bits and a quotient \(\text{hi} = w \gg k\); a contiguous window of at most 63 hi values receives Huffman codewords, and anything outside the window goes through an escape codeword followed by hi written raw. The split \((k, \text{window}, \text{escape})\) is derived per cell by minimizing the model-expected code length over \(k\). On these laws the low-scale cells keep \(k = 0\) with the window centered on zero and an escape for the tails; the high-scale cells shift low bits out raw until the remaining hi values all fit the window, with no escape.
sp_all = builds[0].splits
fig, axs = plt.subplots(1, 2, figsize=(11, 3.5))
axs[0].step(range(G), [sp.k for sp in sp_all], where='mid')
axs[0].set_xlabel('cell j')
axs[0].set_ylabel('k (raw low bits)')
axs[1].step(range(G), [sp.esc for sp in sp_all], where='mid')
axs[1].set_xlabel('cell j')
axs[1].set_ylabel('escape present')
plt.tight_layout()
j_esc, j_demo = 12, 30 # a quiet cell with escape; a busy cell
fig, axs = plt.subplots(1, 2, figsize=(11, 4))
for ax, j in zip(axs, (j_esc, j_demo)):
sp = sp_all[j]
sp_inline = GGDLPC.derive_split(pmfs[j])
assert (sp_inline.k, sp_inline.win_lo, sp_inline.n_win, sp_inline.esc) == \
(sp.k, sp.win_lo, sp.n_win, sp.esc)
assert np.array_equal(sp_inline.lengths, sp.lengths)
hi_of_w = np.arange(nsym) >> sp.k
in_win = (hi_of_w >= sp.win_lo) & (hi_of_w < sp.win_lo + sp.n_win)
print(f'cell {j}: k = {sp.k}, hi window = [{sp.win_lo}, {sp.win_lo + sp.n_win}) '
f'of {((nsym - 1) >> sp.k) + 1} hi values, escape = {bool(sp.esc)} '
f'({sp.hi_bits} raw hi bits when taken), '
f'escape mass = {pmfs[j][~in_win].sum():.2e}')
ax.semilogy(r_axis, pmfs[j], lw=0.8, label='folded pmf')
ax.semilogy(r_axis[in_win], pmfs[j][in_win], 'r.', ms=3, label='hi window (Huffman)')
ax.set_ylim(1e-13, 1)
ax.set_xlabel('residual r')
ax.set_ylabel('q(r)')
ax.set_title(f'cell {j}')
ax.legend()
plt.tight_layout()
sp = sp_all[j_demo]
cell 12: k = 0, hi window = [97, 160) of 257 hi values, escape = True (9 raw hi bits when taken), escape mass = 4.37e-10
cell 30: k = 3, hi window = [0, 33) of 33 hi values, escape = False (0 raw hi bits when taken), escape mass = 0.00e+00
Canonical Huffman lengths#
The split’s Huffman table is built over the windowed hi alphabet (plus the escape symbol): the hi-marginal comes from cascaded pairwise sums of the pmf, the window’s masses get exact Huffman code lengths, and canonical codewords follow by ordering symbols by (length, index) and assigning consecutive codes. Only the lengths need to be stored or regenerated; the codewords are deterministic.
import heapq
def huffman_lengths(weights):
w = np.asarray(weights, np.float64)
idx = np.nonzero(w > 0)[0]
heap = [(w[i], n) for n, i in enumerate(idx)]
heapq.heapify(heap)
parent = {}
nxt = len(idx)
while len(heap) > 1:
w1, n1 = heapq.heappop(heap)
w2, n2 = heapq.heappop(heap)
parent[n1] = parent[n2] = nxt
heapq.heappush(heap, (w1 + w2, nxt))
nxt += 1
L = np.zeros_like(w)
depths = {heap[0][1]: 0}
for n in range(nxt - 1, -1, -1):
if n in parent:
depths[n] = depths[parent[n]] + 1
for n, i in enumerate(idx):
L[i] = depths[n]
L[w <= 0] = L.max()
return L
def canonical_codes(lengths):
lengths = np.asarray(lengths, np.int64)
order = sorted(range(len(lengths)), key=lambda s: (lengths[s], s))
codes = np.zeros(len(lengths), np.int64)
code, prev = 0, int(lengths[order[0]])
for sym in order:
code <<= int(lengths[sym]) - prev
prev = int(lengths[sym])
codes[sym] = code
code += 1
return codes
# hi-marginal of the cell pmf: cascaded pairwise sums, one level per low bit
qh = pmfs[j_demo].copy()
for _ in range(sp.k):
if len(qh) % 2:
qh = np.concatenate([qh, [0.0]])
qh = qh.reshape(-1, 2).sum(axis=1)
qh = qh[:((nsym - 1) >> sp.k) + 1]
if sp.esc:
csum = np.concatenate([[0.0], np.cumsum(qh)])
esc_mass = max(1.0 - float(csum[sp.win_lo + sp.n_win] - csum[sp.win_lo]), 0.0)
probs = np.concatenate([qh[sp.win_lo:sp.win_lo + sp.n_win], [esc_mass]])
else:
probs = qh
L = huffman_lengths(np.maximum(probs, 1e-12)).astype(np.int64)
assert np.array_equal(L, sp.lengths)
codes = canonical_codes(L)
order = np.argsort(L, kind='stable')
print(f'cell {j_demo}, shortest codewords (each hi covers 2^k = {1 << sp.k} residuals):')
for s in order[:8]:
if sp.esc and s == sp.n_win:
print(f' escape len {L[s]:2d} code {codes[s]:0{L[s]}b}')
continue
w_lo = (sp.win_lo + s) << sp.k
print(f' r in [{w_lo - N:+4d}, {w_lo - N + (1 << sp.k) - 1:+4d}] '
f'len {L[s]:2d} code {codes[s]:0{L[s]}b}')
plt.figure(figsize=(6, 3.5))
plt.plot(-np.log2(np.maximum(probs, 1e-12)), L, '.')
plt.plot([0, L.max()], [0, L.max()], 'r--', lw=0.8)
plt.xlabel('-log2(mass) (ideal bits)')
plt.ylabel('Huffman length (bits)')
plt.tight_layout()
cell 30, shortest codewords (each hi covers 2^k = 8 residuals):
r in [ -8, -1] len 2 code 00
r in [ +0, +7] len 2 code 01
r in [ -16, -9] len 3 code 100
r in [ +8, +15] len 3 code 101
r in [ -24, -17] len 4 code 1100
r in [ +16, +23] len 4 code 1101
r in [ -32, -25] len 5 code 11100
r in [ +24, +31] len 5 code 11101
The run mode#
When a cell’s pmf puts more than half its mass on zero, symbol-by-symbol coding costs at least one bit per zero. Starting from the lowest cell, each candidate cell’s model-expected rate under elementary Golomb zero-run coding (order \(2^g\), with an exclusion table for the symbol that interrupts a run) is compared against its symbol table; run mode extends cell by cell while the run rate is lower. The low-scale cells qualify.
fig, ax = plt.subplots(figsize=(7, 4))
for c in range(len(builds)):
p0 = [cell_pmf(blob['channels'][c], j, N)[N] for j in range(G)]
ax.plot(np.arange(G), p0, '.-', label=f'channel {c} (run_jmax = {builds[c].run_jmax})')
ax.axhline(0.5, color='k', lw=0.8, ls='--')
ax.set_xlabel('cell j')
ax.set_ylabel('q(r = 0)')
ax.legend()
plt.tight_layout()
for c in range(len(builds)):
b = builds[c]
if b.run_jmax >= 0:
print(f'channel {c}: run mode in cells 0..{b.run_jmax}, '
f'Golomb orders g = {b.run_g[:b.run_jmax + 1].tolist()}')
else:
print(f'channel {c}: run mode never triggers')
channel 0: run mode in cells 0..13, Golomb orders g = [6, 5, 5, 4, 4, 3, 3, 3, 2, 2, 1, 1, 1, 1]
channel 1: run mode in cells 0..12, Golomb orders g = [7, 6, 6, 5, 4, 4, 3, 3, 2, 2, 1, 1, 1]
channel 2: run mode in cells 0..13, Golomb orders g = [5, 5, 4, 4, 4, 3, 3, 3, 2, 2, 2, 1, 1, 1]
The next chapter encodes an image with these tables one sample at a time and verifies byte identity with the C engine.