03 - Codifica di sorgente: Huffman e primo teorema di Shannon

Vuoi scrivere i simboli di una sorgente in binario con meno bit possibile, e rileggere il messaggio senza ambiguità. Tre domande, in ordine: quali lunghezze di parola sono possibili (la disuguaglianza di Kraft), quanto può scendere la lunghezza media (mai sotto l'entropia H; un codice per simbolo può restare sotto H + 1) e quale codice raggiunge il minimo (Huffman). In mezzo c'è una distinzione da tenere ferma: il limite per i codici per simbolo e il teorema di codifica di sorgente di Shannon sono due enunciati diversi, il secondo su blocchi lunghi e con N che tende all'infinito. Ogni listato è stato eseguito due volte con uscita identica (Python 3.14.4, numpy 2.3.5, semi fissi, nessuna rete).

Si presuppongono l'articolo 01 (entropia) e l'articolo 02 (divergenza di Kullback-Leibler): il limite inferiore discende dal fatto che una divergenza non è negativa.

Codici prefissi e decodifica istantanea

Un codice assegna a ogni simbolo una parola binaria; la sua lunghezza media è L = Σ p(x) · |c(x)|.

Un codice è prefisso se nessuna parola è l'inizio di un'altra, e questo permette di decodificare leggendo da sinistra a destra senza guardare avanti.

Con A = 0, B = 10, C = 11, la sequenza 10110 si legge 10 | 11 | 0, cioè B, C, A: appena la stringa letta coincide con una parola, quella parola è finita, perché nessun'altra la può estendere.

Un codice è univocamente decodificabile (UD) se due sequenze di simboli diverse non danno mai la stessa stringa. I codici prefissi sono UD, ma non è vero il contrario, ed è la ragione per cui conta chi enuncia la disuguaglianza di Kraft. MacKay (§5.2) riferisce che Kraft dimostrò che, date lunghezze che soddisfano la disuguaglianza, esiste un codice prefisso con quelle lunghezze, e che il viceversa, cioè che l'univoca decodificabilità implica la disuguaglianza, è di McMillan (1956): per questo la si chiama anche disuguaglianza di Kraft-McMillan.

Il listato enumera tutti i codici con da 1 a 4 parole lunghe da 1 a 3 bit (1470 insiemi) e per ognuno calcola la somma di Kraft Σ 2^−l con frazioni esatte, se è prefisso, e se è UD.

Per l'univoca decodificabilità usa il test di Sardinas-Patterson, che qui è un'implementazione mia, non verificata su una fonte primaria; come secondo controllo indipendente cerca per forza bruta due sequenze di al più 6 parole che producano la stessa stringa.

from fractions import Fraction
from itertools import combinations, product


def kraft(parole):
    return sum(Fraction(1, 2 ** len(w)) for w in parole)


def senza_prefissi(parole):
    return not any(a != b and b.startswith(a) for a in parole for b in parole)


def univocamente_decodificabile(parole):
    """Test di Sardinas-Patterson: il codice è UD se nessun 'resto' diventa una parola."""
    C = set(parole)
    resti = {b[len(a):] for a in C for b in C if a != b and b.startswith(a)}
    visti = set()
    while resti:
        if resti & C:
            return False
        visti |= resti
        nuovi = set()
        for r in resti:
            for c in C:
                if c.startswith(r) and c != r:
                    nuovi.add(c[len(r):])
                if r.startswith(c) and r != c:
                    nuovi.add(r[len(c):])
        resti = nuovi - visti
    return True


def ambiguo_per_forza_bruta(parole, max_parole=6):
    """Secondo controllo, indipendente: due sequenze di parole diverse che danno la stessa stringa."""
    visti = {}
    for k in range(1, max_parole + 1):
        for seq in product(parole, repeat=k):
            stringa = "".join(seq)
            if stringa in visti and visti[stringa] != seq:
                return True
            visti.setdefault(stringa, seq)
    return False


tutte = ["".join(b) for n in (1, 2, 3) for b in product("01", repeat=n)]   # 14 parole
tot = ud = prefix = ud_non_prefix = 0
violazioni_kraft = discordanze = 0
for n in range(1, 5):
    for C in combinations(tutte, n):
        tot += 1
        p, u = senza_prefissi(C), univocamente_decodificabile(C)
        prefix += p
        ud += u
        ud_non_prefix += u and not p
        violazioni_kraft += (p or u) and kraft(C) > 1
        discordanze += (not u) != ambiguo_per_forza_bruta(C)
print(f"codici con 1-4 parole lunghe 1-3 bit: {tot}")
print(f"  senza prefissi: {prefix}   univocamente decodificabili: {ud}   UD ma con prefissi: {ud_non_prefix}")
print(f"  codici UD o prefissi con somma di Kraft > 1: {violazioni_kraft}")
print(f"  discordanze fra Sardinas-Patterson e la forza bruta (fino a 6 parole): {discordanze}")
for C in (["0", "10", "11"], ["0", "01", "11"], ["0", "01", "10"]):
    print(C, "Kraft =", kraft(C), "| senza prefissi:", senza_prefissi(C), "| UD:", univocamente_decodificabile(C))

# Output:
# codici con 1-4 parole lunghe 1-3 bit: 1470
#   senza prefissi: 458   univocamente decodificabili: 686   UD ma con prefissi: 228
#   codici UD o prefissi con somma di Kraft > 1: 0
#   discordanze fra Sardinas-Patterson e la forza bruta (fino a 6 parole): 0
# ['0', '10', '11'] Kraft = 1 | senza prefissi: True | UD: True
# ['0', '01', '11'] Kraft = 1 | senza prefissi: False | UD: True
# ['0', '01', '10'] Kraft = 1 | senza prefissi: False | UD: False

Su 1470 codici, 458 sono prefissi e 686 UD: 228 sono UD pur avendo parole che iniziano con altre. Nessuno dei codici UD o prefissi ha somma di Kraft maggiore di 1, e i due metodi di decisione dell'UD non discordano mai.

Le ultime tre righe sono istruttive. {0, 10, 11} è un codice prefisso. {0, 01, 11} non lo è, ma è UD (lette al contrario, le parole 0, 10, 11 formano un codice prefisso); per decodificarlo, però, bisogna talvolta guardare avanti, quindi non è istantaneo. {0, 01, 10} ha anch'esso somma di Kraft 1, e non è UD: la stringa 010 si legge 0 | 10 oppure 01 | 0.

La disuguaglianza è dunque una condizione sulle lunghezze, non sul codice: non basta calcolare la somma per sapere che un dato codice funziona.

La disuguaglianza di Kraft, e perché serve

La disuguaglianza dice che le lunghezze l_i di un codice prefisso binario soddisfano Σ 2^−l_i ≤ 1.

L'idea si vede su un albero binario completo di profondità L = max l_i, con 2^L foglie: una parola di lunghezza l «occupa» le 2^L−l foglie che stanno sotto di lei, e perché il codice sia prefisso queste foglie non si sovrappongono. Quindi Σ 2^L−l_i ≤ 2^L, e dividendo per 2^L si ottiene la disuguaglianza.

Ogni parola corta costa molto spazio (una parola di 1 bit occupa metà albero), e questo è il vincolo che impedisce di dare lunghezze piccole a tutti: tre parole da 1 bit hanno somma 3/2 > 1 e non esistono.

Il senso pratico è che la disuguaglianza trasforma un problema di codici in un problema di numeri: ogni scelta di lunghezze con somma al più 1 è realizzabile da un codice prefisso (MacKay, §5.2), e cercare il codice migliore significa cercare le lunghezze intere che minimizzano L sotto quel vincolo.

Il limite: l'entropia, e un bit di margine

Sia z = Σ 2^−l_i ≤ 1 e sia q_i = 2^−l_i/z, che è una distribuzione.

Dalla definizione, l_i = log₂(1/q_i) − log₂ z, e sostituendo nella lunghezza media: L = H(p) + D(p||q) + log₂(1/z).

La divergenza è non negativa (articolo 02) e log₂(1/z) lo è perché z ≤ 1, quindi L ≥ H, con uguaglianza solo se p_i = 2^−l_i per ogni i: le probabilità devono essere potenze di 1/2. È la dimostrazione di MacKay (§5.3), che la enuncia per ogni codice univocamente decodificabile e usa proprio la disuguaglianza di Kraft-McMillan.

Per il limite superiore si scelgono lunghezze un po' più grandi di quelle ideali, l_i = ⌈log₂(1/p_i)⌉.

Soddisfano Kraft, perché 2^−⌈x⌉ ≤ 2^−x e Σ p_i = 1, e superano la lunghezza ideale di meno di un bit.

Ne segue il teorema 5.1 di MacKay, per i codici per simbolo: esiste un codice prefisso con H(X) ≤ L < H(X) + 1. Il «+1» è enorme in proporzione quando l'entropia è piccola, come si vede fra poco.

L'algoritmo di Huffman

MacKay (§5.5) descrive l'algoritmo in due frasi: prendi i due simboli meno probabili, che avranno le parole più lunghe, di uguale lunghezza e diverse solo nell'ultimo bit; fondili in un unico simbolo e ripeti. Ogni passo riduce di uno l'alfabeto, quindi servono |A| − 1 passi.

Il listato lo realizza con una coda di priorità (heapq): ogni elemento è un peso, un contatore di spareggio e il dizionario delle parole dei simboli fusi finora, e a ogni fusione si antepone 0 al ramo del primo nodo estratto e 1 al secondo.

Il contatore ordine rende il risultato deterministico quando due pesi sono uguali, ed è la sola scelta arbitraria: un altro spareggio dà parole diverse, ma la stessa lunghezza media.

Le altre due funzioni calcolano L con frazioni esatte e l'entropia in bit. Il blocco non stampa nulla: è la definizione usata dai listati seguenti.

import heapq
from fractions import Fraction
from itertools import count
from math import log2


def huffman(pesi):
    """pesi: dict simbolo -> peso (int, Fraction o float). Ritorna dict simbolo -> parola."""
    if len(pesi) == 1:
        return {s: "0" for s in pesi}
    ordine = count()                         # spareggio deterministico fra pesi uguali
    heap = [(w, next(ordine), {s: ""}) for s, w in pesi.items()]
    heapq.heapify(heap)
    while len(heap) > 1:
        w1, _, c1 = heapq.heappop(heap)      # i due nodi meno probabili...
        w2, _, c2 = heapq.heappop(heap)
        fuso = {s: "0" + p for s, p in c1.items()} | {s: "1" + p for s, p in c2.items()}
        heapq.heappush(heap, (w1 + w2, next(ordine), fuso))   # ...diventano un nodo solo
    return heap[0][2]


def lunghezza_media(pesi, codice):
    tot = sum(pesi.values())
    return sum(Fraction(pesi[s]) * len(codice[s]) for s in pesi) / tot


def entropia(pesi):
    tot = sum(pesi.values())
    return sum(float(w / tot) * log2(tot / w) for w in pesi.values())

Il listato successivo applica la funzione a due sorgenti: una con probabilità potenze di 1/2 (0,5, 0,25, 0,125, 0,125) e una a cinque simboli (0,35, 0,25, 0,20, 0,12, 0,08), scelte per l'esempio.

Per ognuna stampa le parole, L come frazione esatta, H e la verifica di H ≤ L < H + 1.

R = Fraction
sorgenti = {
    "potenze di 2": {"A": R(1, 2), "B": R(1, 4), "C": R(1, 8), "D": R(1, 8)},
    "cinque simboli": {"A": R(35, 100), "B": R(25, 100), "C": R(20, 100), "D": R(12, 100), "E": R(8, 100)},
}
for nome, pesi in sorgenti.items():
    codice = huffman(pesi)
    L, H = lunghezza_media(pesi, codice), entropia(pesi)
    print(nome)
    for s in sorted(codice, key=lambda s: (len(codice[s]), s)):
        print(f"   {s} p={float(pesi[s]):.3f}  {codice[s]}")
    print(f"   L = {L} = {float(L):.4f} bit   H = {H:.4f} bit   L - H = {float(L) - H:.4f}   (H <= L < H+1: {H <= L < H + 1})")

# Output:
# potenze di 2
#    A p=0.500  0
#    B p=0.250  10
#    C p=0.125  110
#    D p=0.125  111
#    L = 7/4 = 1.7500 bit   H = 1.7500 bit   L - H = 0.0000   (H <= L < H+1: True)
# cinque simboli
#    A p=0.350  11
#    B p=0.250  10
#    C p=0.200  00
#    D p=0.120  011
#    E p=0.080  010
#    L = 11/5 = 2.2000 bit   H = 2.1531 bit   L - H = 0.0469   (H <= L < H+1: True)

Con probabilità potenze di 1/2 il codice raggiunge l'entropia: L = 7/4 = H = 1,75 bit, e il codice a lunghezze 1, 2, 3, 3 è proprio quello di cui parla il caso di uguaglianza. Con cinque simboli L = 11/5 = 2,2 bit contro H = 2,1531: lo scarto è 0,0469 bit.

Un altro spareggio fra i due pesi uguali (0,20 e 0,20) darebbe parole diverse, ma la stessa lunghezza media.

Che il codice sia ottimale non si deve credere sulla parola. Il listato seguente lo verifica in due modi.

Prima confronta, su 3000 distribuzioni casuali con da 2 a 6 simboli (pesi interi, quindi aritmetica esatta), la L di Huffman con il minimo ottenuto enumerando tutte le liste di lunghezze (da 1 a n−1, perché un albero binario completo con n foglie non è più profondo) che soddisfano Kraft.

Poi codifica e decodifica tutte le stringhe fino a 6 simboli sui cinque simboli dell'esempio (19.530) e altre 2000 stringhe casuali da 300 simboli, e misura la lunghezza media osservata.

import random
from itertools import product

rng = random.Random(2026)


# Per ogni n, tutte le liste di lunghezze (da 1 a n-1) che soddisfano Kraft, in aritmetica intera
ammesse = {
    n: [l for l in product(range(1, n), repeat=n)
        if sum(2 ** (n - 1 - x) for x in l) <= 2 ** (n - 1)]
    for n in range(2, 7)
}

# 1) Ottimalità: confronto con la ricerca esaustiva fra TUTTE le lunghezze che soddisfano Kraft
casi, diversi = 3000, 0
for _ in range(casi):
    n = rng.randint(2, 6)
    pesi = {i: rng.randint(1, 50) for i in range(n)}
    tot = sum(pesi.values())
    L_huffman = lunghezza_media(pesi, huffman(pesi))
    minimo = min(sum(pesi[i] * l[i] for i in range(n)) for l in ammesse[n])
    diversi += L_huffman != Fraction(minimo, tot)
print(f"distribuzioni con L(Huffman) diversa dal minimo esaustivo: {diversi} su {casi}")
print("liste di lunghezze ammesse da Kraft, per n=2..6 simboli:", [len(ammesse[n]) for n in range(2, 7)])

# 2) Andata e ritorno: tutte le stringhe fino a 6 simboli, più 2000 stringhe casuali da 300
pesi = {"A": 35, "B": 25, "C": 20, "D": 12, "E": 8}
codice = huffman(pesi)
inverso = {p: s for s, p in codice.items()}


def decodifica(bits):
    out, corrente = [], ""
    for b in bits:
        corrente += b
        if corrente in inverso:              # prefisso: la parola è finita non appena la riconosci
            out.append(inverso[corrente])
            corrente = ""
    assert corrente == ""
    return "".join(out)


controllate = 0
for n in range(1, 7):
    for s in product("ABCDE", repeat=n):
        testo = "".join(s)
        assert decodifica("".join(codice[c] for c in testo)) == testo
        controllate += 1
medie = []
for _ in range(2000):
    testo = "".join(rng.choices("ABCDE", weights=[35, 25, 20, 12, 8], k=300))
    codificato = "".join(codice[c] for c in testo)
    assert decodifica(codificato) == testo
    medie.append(len(codificato) / len(testo))
    controllate += 1
media = sum(medie) / len(medie)
errore = (sum((m - media) ** 2 for m in medie) / (len(medie) - 1)) ** 0.5 / len(medie) ** 0.5
print(f"stringhe codificate e decodificate senza errori: {controllate}")
print(f"lunghezza media osservata su 600000 simboli: {media:.4f} +/- {errore:.4f} bit (teorica 2.2000)")

# Output:
# distribuzioni con L(Huffman) diversa dal minimo esaustivo: 0 su 3000
# liste di lunghezze ammesse da Kraft, per n=2..6 simboli: [1, 4, 32, 392, 6154]
# stringhe codificate e decodificate senza errori: 21530
# lunghezza media osservata su 600000 simboli: 2.1987 +/- 0.0005 bit (teorica 2.2000)

Su 3000 distribuzioni, 0 differiscono dal minimo esaustivo; le liste di lunghezze ammesse da Kraft crescono in fretta (1, 4, 32, 392, 6154), il che spiega perché si cerca un algoritmo e non una enumerazione.

Le 21.530 stringhe tornano identiche dopo la codifica e la decodifica, e su 600.000 simboli la lunghezza media osservata è 2,1987 ± 0,0005 bit contro i 2,2000 teorici: è a 2,6 errori standard, più di quanto mi aspettassi. Con altri otto semi, in una prova a parte non riportata nel listato, la media su 600.000 simboli va da 2,1993 a 2,2007, quindi lo scarto è una fluttuazione di questo campione e non un difetto del codice, che l'andata e ritorno su 21.530 stringhe esclude.

L'ottimalità riguarda i codici binari per simbolo; MacKay (soluzione dell'esercizio 5.32) descrive la modifica necessaria per codificare con un alfabeto di q simboli.

Dove Huffman perde, e il teorema che lo supera

Un codice per simbolo non può usare meno di un bit per simbolo, anche quando il simbolo è quasi certo. Per una sorgente con probabilità (0,99, 0,01) l'entropia è 0,0808 bit, ma Huffman assegna 1 bit a ciascuno dei due simboli e la lunghezza media è 1: dodici volte l'entropia.

MacKay (§5.6) lo esprime dicendo che i codici di Huffman, benché presentati come «ottimali», hanno molti difetti pratici: sono ottimali fra i codici per simbolo, ma in pratica non vogliamo un codice per simbolo.

Il rimedio è codificare blocchi di N simboli come se fossero un simbolo solo. Il listato lo prova per N da 1 a 12: la probabilità di un blocco dipende solo dal numero di simboli rari, e usando il peso intero 99^N−k (la probabilità vera è il peso diviso 100^N) la lunghezza media per simbolo è esatta.

Accanto stampa il limite H + 1/N, che segue dal teorema 5.1 applicato alla sorgente a blocchi, la cui entropia è N volte quella del singolo simbolo.

from itertools import product
from math import log2

H = 0.99 * log2(1 / 0.99) + 0.01 * log2(1 / 0.01)
print(f"entropia di un simbolo con p = (0.99, 0.01): H = {H:.4f} bit")
print("  N   parole   L/N (bit per simbolo)   H + 1/N")
for N in (1, 2, 4, 6, 8, 10, 12):
    # blocchi di N simboli: la probabilità dipende solo dal numero k di simboli rari
    # peso intero 99^(N-k): la probabilità vera è peso / 100^N, quindi il calcolo è esatto
    pesi = {b: 99 ** (N - sum(b)) for b in product((0, 1), repeat=N)}
    L = lunghezza_media(pesi, huffman(pesi))
    print(f"{N:3d}  {len(pesi):6d}   {float(L) / N:20.4f}   {H + 1 / N:8.4f}")

# Output:
# entropia di un simbolo con p = (0.99, 0.01): H = 0.0808 bit
#   N   parole   L/N (bit per simbolo)   H + 1/N
#   1       2                 1.0000     1.0808
#   2       4                 0.5150     0.5808
#   4      16                 0.2727     0.3308
#   6      64                 0.1953     0.2475
#   8     256                 0.1572     0.2058
#  10    1024                 0.1362     0.1808
#  12    4096                 0.1222     0.1641

La lunghezza per simbolo scende da 1 a 0,1222 con blocchi di 12 simboli, sempre sotto il limite H + 1/N, e tende all'entropia, ma lentamente: a N = 12 è ancora il 51% sopra 0,0808, e il codice ha 4096 parole.

Questa è la forma operativa del teorema di codifica di sorgente. MacKay (teorema 4.1, p. 78) lo enuncia così: N variabili casuali indipendenti e identicamente distribuite, ciascuna con entropia H(X), si possono comprimere in più di N H(X) bit con rischio trascurabile di perdere informazione quando N → ∞; al contrario, se si comprime in meno di N H(X) bit è praticamente certo che si perda informazione.

Che cosa garantisce, e che cosa no. Garantisce che H è il confine: sopra si può arrivare, sotto no.

Non garantisce niente per un N finito, non costruisce il codice, e presuppone che si conosca la distribuzione e che la sorgente sia i.i.d.; Shannon, nel suo articolo del 1948 (teorema 3, sulle sequenze tipiche), lavora con sorgenti ergodiche, che sono più generali. Non dice nemmeno che il costo in memoria e in tempo del codice sia accettabile: il codice a blocchi di 12 ha 4096 parole, a blocchi di 30 avrebbe 2³⁰ parole possibili, oltre un miliardo.

Codifica aritmetica: l'intera sequenza come un numero

La codifica aritmetica abolisce l'obbligo di un numero intero di bit per simbolo. MacKay (cap. 6) ne attribuisce l'invenzione a Elias, Rissanen e Pasco, e la resa pratica a Witten e altri (1987). L'idea è assegnare a ogni sequenza un sotto-intervallo di [0, 1) la cui larghezza è la probabilità della sequenza: ogni simbolo restringe l'intervallo in proporzione alla propria probabilità.

Poi si sceglie una frazione binaria dentro l'intervallo: se la larghezza è P, bastano ⌈log₂(1/P)⌉ bit.

Il listato implementa il codificatore e il decodificatore con frazioni esatte, quindi senza errori di arrotondamento. intervallo costruisce la tabella delle probabilità cumulate e restringe low e width simbolo per simbolo; codifica cerca il più piccolo k tale che 2^−k ≤ width e prende il più piccolo multiplo di 2^−k che non sia inferiore a low, che per costruzione sta nell'intervallo. Il decodificatore deve conoscere il numero N di simboli (o serve un simbolo di fine).

Lo provo su tutte le sequenze fino a 10 simboli e su 300 sequenze casuali da 1000, con P(A) = 0,99 e P(B) = 0,01.

import random
from math import ceil, log2


def intervallo(seq, prob):
    """Intervallo [low, low+width) di [0,1) associato alla sequenza, con frazioni esatte."""
    cum, c = {}, Fraction(0)
    for s in sorted(prob):
        cum[s], c = c, c + prob[s]
    low, width = Fraction(0), Fraction(1)
    for s in seq:
        low, width = low + width * cum[s], width * prob[s]
    return low, width, cum


def codifica(seq, prob):
    low, width, _ = intervallo(seq, prob)
    k = 0
    while Fraction(1, 2 ** k) > width:        # più piccolo k con 2^-k <= larghezza
        k += 1
    m = -((-low.numerator * 2 ** k) // low.denominator)       # ceil(low * 2^k)
    return format(m, f"0{k}b") if k else ""


def decodifica(bits, n, prob):
    x = Fraction(int(bits, 2), 2 ** len(bits)) if bits else Fraction(0)
    _, _, cum = intervallo("", prob)
    low, width, out = Fraction(0), Fraction(1), []
    for _ in range(n):
        for s in sorted(prob):
            a = low + width * cum[s]
            if a <= x < a + width * prob[s]:
                out.append(s)
                low, width = a, width * prob[s]
                break
    return "".join(out)


prob = {"A": Fraction(99, 100), "B": Fraction(1, 100)}
rng = random.Random(5)
# Tutte le sequenze fino a 10 simboli e 300 sequenze casuali da 1000: andata e ritorno esatti
for n in range(1, 11):
    for i in range(2 ** n):
        s = "".join("AB"[(i >> j) & 1] for j in range(n))
        assert decodifica(codifica(s, prob), n, prob) == s
lunghezze, ideali = [], []
for _ in range(300):
    s = "".join(rng.choices("AB", weights=[99, 1], k=1000))
    bits = codifica(s, prob)
    assert decodifica(bits, 1000, prob) == s
    _, width, _ = intervallo(s, prob)
    lunghezze.append(len(bits))
    ideali.append(ceil(-log2(float(width))) if width > 0 else None)
H = 0.99 * log2(1 / 0.99) + 0.01 * log2(1 / 0.01)
media = sum(lunghezze) / len(lunghezze)
print("andata e ritorno esatti su 2046 sequenze corte e 300 da 1000 simboli")
errore = (sum((x - media) ** 2 for x in lunghezze) / (len(lunghezze) - 1)) ** 0.5 / len(lunghezze) ** 0.5
print(f"lunghezza media del codice aritmetico: {media:.1f} +/- {errore:.1f} bit per 1000 simboli (N·H = {1000 * H:.1f})")
print(f"scarto massimo da ceil(-log2 P(sequenza)): {max(abs(a - b) for a, b in zip(lunghezze, ideali))} bit")
print("Huffman simbolo per simbolo: 1000 bit per 1000 simboli")

# Output:
# andata e ritorno esatti su 2046 sequenze corte e 300 da 1000 simboli
# lunghezza media del codice aritmetico: 81.5 +/- 1.3 bit per 1000 simboli (N·H = 80.8)
# scarto massimo da ceil(-log2 P(sequenza)): 0 bit
# Huffman simbolo per simbolo: 1000 bit per 1000 simboli

Tutte le sequenze tornano identiche, e la lunghezza del codice coincide sempre con ⌈−log₂ P(sequenza)⌉ (scarto massimo 0 bit).

In media sono 81,5 ± 1,3 bit per 1000 simboli, contro N H = 80,8 e i 1000 bit di Huffman per simbolo: lo scarto di 0,7 bit è dentro l'errore e dello stesso ordine dell'arrotondamento per eccesso (al più 1 bit). Il listato usa razionali a precisione illimitata, che crescono con il messaggio: i codificatori reali usano aritmetica intera a precisione finita, e quanto perdano non l'ho misurato. Come per Huffman, il decodificatore deve conoscere lo stesso modello probabilistico.

I compressori reali, su dati costruiti qui

Il confronto con zlib, bz2 e lzma della libreria standard usa due sorgenti costruite da me, di 200.000 simboli, un byte per simbolo. La prima ha simboli indipendenti con probabilità (1/2, 1/4, 1/8, 1/8): Huffman per simbolo raggiunge l'entropia. La seconda è una catena di Markov a due stati che resta nello stato corrente con probabilità 0,95: l'entropia di un simbolo preso da solo è 1 bit, ma l'entropia condizionata al simbolo precedente, H(X_n|X_n−1) = 0,95 log₂(1/0,95) + 0,05 log₂(1/0,05) = 0,2864 bit, è molto più piccola (è la quantità dell'articolo 01 applicata a una sorgente con memoria).

Il listato stampa anche una riga di zlib limitato alla sola codifica di Huffman (Z_HUFFMAN_ONLY, senza corrispondenze).

import bz2
import lzma
import zlib

import numpy as np

rng = np.random.default_rng(3)
N = 200_000


def solo_huffman(d):
    c = zlib.compressobj(9, zlib.DEFLATED, 15, 9, zlib.Z_HUFFMAN_ONLY)   # niente corrispondenze LZ77
    return c.compress(d) + c.flush()


def bit_per_simbolo(dati):
    return {nome: 8 * len(f(dati)) / N for nome, f in (
        ("zlib -9", lambda d: zlib.compress(d, 9)),
        ("zlib solo Huffman", solo_huffman),
        ("bz2 -9", lambda d: bz2.compress(d, 9)),
        ("lzma -9", lambda d: lzma.compress(d, preset=9)),
    )}


# Sorgente 1: simboli indipendenti, p = (1/2, 1/4, 1/8, 1/8), un byte per simbolo
x = rng.choice(4, size=N, p=[0.5, 0.25, 0.125, 0.125])
dati = (x + 65).astype(np.uint8).tobytes()                                    # 'A'..'D'
freq = np.bincount(x) / N
H1 = float(-(freq * np.log2(freq)).sum())
L1 = float(lunghezza_media(dict(enumerate(np.bincount(x).tolist())), huffman(dict(enumerate(np.bincount(x).tolist())))))
print(f"sorgente 1 (i.i.d.): entropia empirica {H1:.4f} bit/simbolo; Huffman: {L1:.4f}")
for nome, v in bit_per_simbolo(dati).items():
    print(f"   {nome}: {v:.4f} bit/simbolo")
# DEFLATE aggiunge ai simboli un codice di fine blocco (circa uno ogni 32768 simboli): lo aggiungo ai conti
con_eob = dict(enumerate(np.bincount(x).tolist())) | {"EOB": 6}
parole = huffman(con_eob)
print(f"   Huffman sugli stessi conteggi + 6 fine-blocco: {float(lunghezza_media({k: v for k, v in con_eob.items() if k != 'EOB'}, parole)):.4f} bit/simbolo")

# Sorgente 2: catena di Markov a due stati, resta nello stato con probabilità 0.95
resta = rng.random(N) < 0.95
stato = np.zeros(N, dtype=np.int64)
for i in range(1, N):
    stato[i] = stato[i - 1] if resta[i] else 1 - stato[i - 1]
dati = (stato + 48).astype(np.uint8).tobytes()                                # '0' e '1'
marginale = np.bincount(stato) / N
H_marg = float(-(marginale * np.log2(marginale)).sum())
tasso = 0.95 * np.log2(1 / 0.95) + 0.05 * np.log2(1 / 0.05)   # H(X_n | X_{n-1})
print(f"sorgente 2 (Markov): entropia marginale {H_marg:.4f}; "
      f"entropia condizionata H(X_n|X_(n-1)) = {tasso:.4f} bit/simbolo")
conti = dict(enumerate(np.bincount(stato).tolist()))
print(f"   Huffman simbolo per simbolo: {float(lunghezza_media(conti, huffman(conti))):.4f} bit/simbolo")
con_eob = conti | {"EOB": 6}
print(f"   Huffman sugli stessi conteggi + 6 fine-blocco: {float(lunghezza_media(conti, huffman(con_eob))):.4f} bit/simbolo")
for nome, v in bit_per_simbolo(dati).items():
    print(f"   {nome}: {v:.4f} bit/simbolo")

# Output:
# sorgente 1 (i.i.d.): entropia empirica 1.7510 bit/simbolo; Huffman: 1.7510
#    zlib -9: 2.1342 bit/simbolo
#    zlib solo Huffman: 1.8789 bit/simbolo
#    bz2 -9: 2.0758 bit/simbolo
#    lzma -9: 2.0128 bit/simbolo
#    Huffman sugli stessi conteggi + 6 fine-blocco: 1.8759 bit/simbolo
# sorgente 2 (Markov): entropia marginale 1.0000; entropia condizionata H(X_n|X_(n-1)) = 0.2864 bit/simbolo
#    Huffman simbolo per simbolo: 1.0000 bit/simbolo
#    Huffman sugli stessi conteggi + 6 fine-blocco: 1.4986 bit/simbolo
#    zlib -9: 0.3930 bit/simbolo
#    zlib solo Huffman: 1.4930 bit/simbolo
#    bz2 -9: 0.3496 bit/simbolo
#    lzma -9: 0.4083 bit/simbolo

Sulla sorgente senza memoria nessun compressore generico scende sotto l'entropia, e nessuno la raggiunge: lzma 2,01, bz2 2,08, zlib 2,13 bit per simbolo contro 1,751.

Interessante è zlib: con le sole corrispondenze LZ77 disattivate scende a 1,879, meglio che con tutte attive, perché su un alfabeto di quattro lettere le ripetizioni sono casuali e sfruttarle costa più di quanto fa risparmiare. Perché 1,879 e non 1,751? Ipotizzo il codice di fine blocco di DEFLATE (circa uno ogni 32.768 simboli, stima mia che non ho verificato sul sorgente di zlib): è un simbolo raro in più, e Huffman allunga di un bit una delle parole.

Con 6 fine-blocco il calcolo dà 1,876, e per la sorgente di Markov 1,4986 contro un misurato di 1,4930: con due soli simboli, uno in più costa mezzo bit.

Sulla sorgente con memoria succede il contrario: Huffman per simbolo non può fare meglio di 1 bit, mentre i compressori arrivano a 0,35-0,41 bit per simbolo (bz2 0,3496, zlib 0,3930, lzma 0,4083), poco più di un terzo, pur non conoscendo il modello.

Non raggiungono i 0,2864 bit dell'entropia condizionata. È la conferma sperimentale che l'entropia di un simbolo preso da solo è la misura giusta solo per sorgenti senza memoria, e che un compressore che sfrutta le ripetizioni sfrutta in realtà la dipendenza fra simboli.

Quanto a DEFLATE, la specifica (RFC 1951) lo definisce come la combinazione dell'algoritmo LZ77 e della codifica di Huffman: il formato di zlib e gzip contiene quindi Huffman come ultimo stadio, ma non è «solo Huffman».

Errori tipici

• Credere che H sia raggiungibile da ogni codice per simbolo. Lo è solo se le probabilità sono potenze di 1/2; altrimenti la lunghezza media sta fra H e H + 1.
• Usare Kraft per concludere che un codice funziona. {0, 01, 10} ha somma 1 e non è univocamente decodificabile.
• Confondere il limite per simbolo con il teorema di Shannon. Il primo vale per ogni sorgente e dà H + 1; il secondo riguarda blocchi, con N che tende all'infinito.
• Dare per scontato che Huffman sia il meglio possibile. È il meglio fra i codici per simbolo: con (0,99, 0,01) fa 1 bit contro 0,0808.
• Aspettarsi che un compressore generico si avvicini all'entropia del simbolo. Sulla sorgente senza memoria resta sopra di 0,26-0,38 bit; su quella con memoria scende molto sotto.

I cinque errori nascono da una stessa confusione: parlare di «compressione» senza dire rispetto a quale modello, a quale classe di codici e con quale lunghezza di blocco. Ogni confronto con l'entropia ha senso solo se si dichiarano questi tre elementi.

Che cosa sai fare adesso

Sai verificare se un insieme di lunghezze è realizzabile da un codice prefisso e sai che la somma di Kraft non dice nulla su un codice dato. Sai dimostrare che nessun codice per simbolo sta sotto H e che esiste uno sotto H + 1, costruire il codice di Huffman, provarlo in andata e ritorno e confrontarlo con la ricerca esaustiva.

Sai anche dove si ferma: per sorgenti molto concentrate servono blocchi o la codifica aritmetica, e il teorema di Shannon dice dove è il confine, non come arrivarci. Il passo seguente, nell'articolo 04, sposta la questione dal comprimere al trasmettere su un canale che sbaglia, e nell'articolo 05 si costruiscono codici correttori. La compressione con perdita è nell'articolo 06, e l'indice è nella pagina della serie.