Coding

The dit.coding module builds both source codes (which compress a source Distribution) and channel codes (which add redundancy to protect against channel noise), and exposes their code-theoretic properties.

Source coding

A source code maps the outcomes of a source to codewords over a radix-ary alphabet (binary by default), with the goal of minimizing the expected number of code symbols per source symbol – the rate. By the source coding theorem [CT06], no uniquely decodable code can have a rate below the source entropy \(\H{X}\).

Base classes

All source codes derive from SourceCoding, which defines the encode(), decode(), and rate() interface and provides the comparisons to the source entropy that are common to every source code:

source_entropy() – \(\H{X}\) in radix-ary digits,
redundancy() – rate - source_entropy,
efficiency() – source_entropy / rate.

Channel codes derive from the companion ChannelCoding base class, described under Channel coding below.

Symbol codes

A SymbolCode assigns one codeword to each source outcome. The classic constructions are available:

shannon() – lengths \(\lceil \log_D 1/p \rceil\), codewords from the cumulative distribution,
fano() – top-down balanced partitioning,
shannon_fano_elias() – lengths \(\lceil \log_D 1/p \rceil + 1\), codewords from the midpoint cumulative distribution,
huffman() – the optimal symbol code,
length_limited_huffman() – the optimal code subject to a maximum codeword length (via package-merge),
golomb() / rice() – optimal codes for geometric integer sources.

Beyond the rate-based properties, a SymbolCode reports the Kraft sum (kraft_sum()), whether it is complete (is_complete()), prefix-free (is_prefix_free()), uniquely decodable (is_uniquely_decodable(), via Sardinas-Patterson), and optimal (is_optimal()).

Example

In [1]: from dit.coding import huffman, shannon

In [2]: d = dit.Distribution(['a', 'b', 'c', 'd', 'e'], [0.4, 0.2, 0.2, 0.1, 0.1])

In [3]: code = huffman(d)

In [4]: code.average_length(), code.source_entropy()
Out[4]: (2.2, 2.1219280948873624)

In [5]: code.is_optimal(), code.is_prefix_free(), code.is_complete()
Out[5]: (True, True, True)

In [6]: seq = list(d.outcomes)

In [7]: code.decode(code.encode(seq)) == seq
Out[7]: True

Huffman is optimal, so its average length is no larger than the Shannon code’s:

In [8]: shannon(d).average_length()
Out[8]: 2.8000000000000003

Universal integer codes

The universal codes encode the positive integers without reference to a source distribution, yet stay within a constant factor of optimal across a wide range of sources: unary(), elias_gamma(), elias_delta(), elias_omega(), and fibonacci(). The universal_code() helper wraps a chosen family into a SymbolCode over an integer-valued source.

Tunstall codes

Where a symbol code maps each symbol to a variable-length codeword, a tunstall() code parses the source into variable-length words and maps each word to a fixed-length codeword. It is the variable-to-fixed dual of Huffman and is realized by TunstallCode.

Polar source coding

Source polarization [Arikan10] applies the Arikan transform on the source side. For \(N = 2^m\) i.i.d. copies of a binary source \(X\) (optionally with side information \(Y\)), the transform \(U^N = X^N G_N\) produces synthesized coordinates whose conditional entropies \(\H{U_i \mid U^{i-1}, Y^N}\) polarize toward \(0\) (deterministic given the past) or \(1\) (uniform given the past) as \(N\) grows. The coordinates that stay near one – the high-entropy set – are exactly what a lossless source code must store; the rest are recovered by sequential decisions and the inverse transform. With side information this is the polar route to Slepian-Wolf coding.

The finite-block utilities are exact (no density evolution), so they are limited to small \(N\):

source_bhattacharyya() – the source Bhattacharyya parameter \(Z(X \mid Y) = 2 \sum_y \sqrt{p(0, y) p(1, y)}\), small when \(X\) is nearly determined by \(Y\) and near one when \(X\) is nearly uniform,
source_polarization_profile() – the per-coordinate conditional entropy and source Bhattacharyya (and an optional Goela-style max_correlation_with_past diagnostic [Goe14]),
source_high_entropy_set() – the coordinates a code keeps; by default the lossless set (every coordinate whose conditional entropy exceeds a tolerance).

The polar_source() constructor builds a PolarSourceCode – an exact finite-block source code (binary source, power-of-two block length, optional decoder side information). Its encode() returns the high-entropy coordinates and decode() fills the low-entropy coordinates by exact maximum a posteriori decisions before inverting the transform. Because the joint table is enumerated exactly, a max_states guard prevents accidental exponential blowups.

In [9]: from dit.coding import polar_source, source_polarization_profile

In [10]: dsbs = dit.Distribution(['00', '01', '10', '11'], [0.45, 0.05, 0.05, 0.45])

In [11]: [round(row['entropy'], 3) for row in source_polarization_profile(dsbs, 4, rv=0, crvs=[1])]
Out[11]: [0.875, 0.485, 0.438, 0.078]

In [12]: code = polar_source(dsbs, 8, rv=0, crvs=[1])

In [13]: code.rate(), code.high_entropy_set
Out[13]: (1.0, [0, 1, 2, 3, 4, 5, 6, 7])

Channel coding

A channel code adds redundancy to a message so that it can be recovered after transmission over a noisy channel. All channel codes in dit.coding are binary (over \(\mathrm{GF}(2)\)). The usual channel is a conditional Distribution \(p(Y \mid X)\) – a discrete memoryless channel applied independently to each transmitted symbol. The Channels page catalogs ready-made channels, including binary_symmetric_channel() and binary_erasure_channel() (re-exported from dit.coding for convenience). Channels with memory can be attached too; see Channels with memory.

Every channel code derives from ChannelCoding, which provides the information-theoretic evaluation common to all of them:

capacity_gap() – channel_capacity(channel) - rate, the gap to the channel coding theorem’s limit [CT06] (capacity is computed with the Blahut–Arimoto solver described on Channels; a channel that only brackets its capacity yields a (lower, upper) bracket on the gap),
probability_of_error() – the block-error probability under the code’s own decoder, computed exactly for small codes (enumerating every message and received word) or by Monte Carlo otherwise.

Linear block codes

A LinearCode is specified by a generator matrix G; the parity-check matrix H is derived as a basis for its null space. It supports hard-decision syndrome decoding (and maximum-likelihood decoding when a channel is given), and reports the minimum_distance(), weight_enumerator(), and error_correcting_capability(). The classical families are available as constructors: repetition(), parity_check(), hamming(), reed_muller(), and golay().

Modern codes

Three modern codes specialize the decoder:

ldpc() / gallager() build an LDPCCode decoded by sum-product belief propagation over its Tanner graph [CT06],
polar() builds a PolarCode, freezing the least reliable synthesized bit-channels (by Bhattacharyya parameter) and decoding by successive cancellation,
convolutional() builds a ConvolutionalCode from generator polynomials, decoded with the Viterbi algorithm (soft-decision when a channel is given).

Channel coding example

The Hamming [7, 4, 3] code corrects any single error; over a binary symmetric channel its exact block-error probability matches the closed form \(1 - (1-p)^7 - 7p(1-p)^6\):

In [14]: from dit.coding import hamming, binary_symmetric_channel

In [15]: code = hamming(3)

In [16]: code.length, code.dimension, code.minimum_distance()
Out[16]: (7, 4, 3)

In [17]: bsc = binary_symmetric_channel(0.05)

In [18]: code.capacity_gap(bsc)
Out[18]: 0.14217447145547246

In [19]: code.probability_of_error(bsc, method='exact')
Out[19]: 0.044380542187499934

In [20]: 1 - 0.95**7 - 7 * 0.05 * 0.95**6
Out[20]: 0.044380542187500316

Channels with memory

Codes do not require a Distribution: every channel-aware method first calls as_channel(), which wraps a conditional distribution as a MemorylessChannel and accepts any other object implementing a small duck-typed protocol:

inputs and outputs – the alphabets, as tuples,
log_likelihood(codeword, received) – \(\log p(y^n \mid x^n)\), used by maximum-likelihood decoding and by the exact block-error probability,
sample(codeword, prng) – one received word, used by Monte Carlo,
symbol_log_likelihoods(received) – an (n, len(inputs)) array of \(\log p(y^n \mid x_i = a)\) up to a per-row constant, from which belief propagation, successive cancellation, and soft Viterbi form their per-position metrics,
optionally capacity(**kwargs) (a float or a (lower, upper) bracket) and bhattacharyya() (needed only to design a polar code),
optionally feedback_initial(), feedback_step(belief, x) and feedback_capacity(**kwargs), used only by the feedback codes.

A finite-state channel \(T(y, s' \mid s, x)\) computes log_likelihood with the forward algorithm and symbol_log_likelihoods with the forward-backward (BCJR) recursions [BCJR74] under uniform i.i.d. inputs – the standard separate detect-then-decode approach. Maximum-likelihood decoding and the exact error probability remain exact over such a channel; the per-symbol decoders are approximate because they ignore the correlations the channel induces between positions.

An object may also define __dit_channel__() returning a protocol object, so another library can make its own channel types attachable without dit knowing about them. For example, sofic finite-state channels (Mealy machines such as the trapdoor channel) work directly:

from dit.coding import hamming
from sofic.examples import trapdoor

code = hamming(3)
code.probability_of_error(trapdoor(), method='exact')
code.capacity_gap(trapdoor())   # a (lower, upper) bracket

Codes with feedback

With a noiseless feedback link the transmitter sees each output before choosing the next input, so a code is a strategy \(x_t = f(m, y_{0:t})\) and the number of channel uses \(N\) may be random. Feedback does not increase the capacity of a memoryless channel [CT06], but it makes capacity-achieving schemes much simpler and improves the error exponent [Bur76]. On channels with memory it can raise the capacity itself; the trapdoor channel goes from about \(0.58\) bits without feedback to \(\log_2 \varphi \approx 0.694\) with it [PCVRW08].

A FeedbackCode keeps a state that both ends know, updated from each output alone. The channel is stepped one use at a time through feedback_initial() and feedback_step(belief, x), which returns every (y, p(y | x, belief), next_belief); MemorylessChannel provides both, and a finite-state channel uses the normalized forward vector as its belief. The exact block-error probability and \(\mathbb{E}[N]\) come from walking the tree of outputs, merging branches that reach the same state, and the rate is \(k / \mathbb{E}[N]\).

erasure_arq() resends each bit until it arrives unerased. On the binary erasure channel it is zero-error with \(\mathbb{E}[N] = k / (1 - \epsilon)\), so it achieves capacity exactly.
stop_and_wait_arq() resends a block code’s codeword until the received word is a codeword; it fails only when the error pattern is itself a codeword.
horstein() is Horstein’s bisection scheme [Hor63]: send whether the message lies above the receiver’s posterior median, stopping after n uses or once the posterior reaches 1 - delta. Posterior matching [SF11] generalizes it, and the Schalkwijk–Kailath scheme [SK66] is its Gaussian-channel analogue (not implemented, since dit channels are discrete).

In [21]: from dit.coding import erasure_arq, horstein, binary_erasure_channel, binary_symmetric_channel

In [22]: erasure_arq(3).performance(binary_erasure_channel(0.3), method='exact')
Out[22]: (1.9231633866980787e-30, 4.285714285714281)

In [23]: 3 / 0.7
Out[23]: 4.285714285714286

In [24]: bsc = binary_symmetric_channel(0.1)

In [25]: [horstein(3, n=n).probability_of_error(bsc, method='exact') for n in (4, 8, 12)]
Out[25]: [0.2709999999999998, 0.07529265999999983, 0.017922541167996155]

APIs

class SourceCoding(dist=None, radix=2)[source]

Abstract base class for source codes.

A source code maps the outcomes of a Distribution to codewords over a radix-ary code alphabet. Concrete subclasses define how outcomes are encoded and decoded, and what the code’s rate is; this base class provides the comparisons to the source entropy that are common to all source codes.

Parameters:

dist (Distribution, None) – The source distribution the code is built for. Required for any of the rate-based properties.
radix (int) – The size of the code alphabet. Default is 2 (binary).

abstractmethod encode(source)[source]: Encode a sequence of source outcomes into a sequence of code symbols.

abstractmethod decode(encoded)[source]: Decode a sequence of code symbols back into source outcomes.

abstractmethod rate()[source]: The expected number of code symbols emitted per source symbol.

source_entropy()[source]

The entropy of the source, in units of radix-ary digits.

This is the fundamental lower bound on the rate of any uniquely decodable source code (the source coding theorem; Cover & Thomas, Ch. 5).

Returns:: H – The source entropy, H[X] / log2(radix).
Return type:: float

redundancy()[source]

The rate in excess of the source entropy, rate - source_entropy.

Returns:: r – The redundancy, in radix-ary digits per source symbol. Always non-negative for a uniquely decodable code.
Return type:: float

efficiency()[source]

The ratio of the source entropy to the rate, in [0, 1].

Returns:: e – The coding efficiency, source_entropy / rate.
Return type:: float

class ChannelCoding(channel=None, radix=2)[source]

Abstract base class for channel codes.

A channel code adds redundancy to a message so that it can be recovered after transmission across a noisy channel. This class is scaffolding for future channel codes (repetition, polar, LDPC, …); no concrete subclass is provided yet.

Concrete subclasses are expected to implement encode() (message -> codeword), decode() (received word -> message), and rate() (message symbols per channel use), and may additionally expose channel-code properties such as the minimum distance, the block error probability for a given channel, and the gap to capacity (via dit.algorithms.channel_capacity()).

Parameters:

channel (Distribution, object, None) – The channel, as a conditional distribution p(output | input) or any object implementing the channel protocol (see as_channel()).
radix (int) – The size of the code alphabet. Default is 2 (binary).

abstractmethod encode(message)[source]: Encode a message into a channel codeword.

abstractmethod decode(received, channel=None)[source]

Decode a received word back into a message.

Soft-decision decoders use channel (a conditional p(Y|X) distribution, or any object implementing the channel protocol) to form log-likelihoods; hard-decision decoders ignore it.

abstractmethod rate()[source]: The number of message symbols transmitted per channel use.

capacity_gap(channel, **kwargs)[source]

The gap between the channel capacity and the code rate.

Parameters:

channel (Distribution or object) – The channel: a conditional distribution p(Y|X), or any object implementing the channel protocol (see as_channel()) and providing capacity().
kwargs (dict) – Passed to the channel’s capacity().

Returns:

gap – capacity(channel) - rate, in bits per channel use. Positive when the code operates below capacity. When the channel only brackets its capacity, a (lower, upper) bracket on the gap.

Return type:

float or tuple

probability_of_error(channel, method='auto', samples=10000, prng=None)[source]

The block error probability of the code over a channel.

The code’s own decode() is used (passing channel so soft-decision decoders engage). With method='exact' every message and every received word is enumerated – feasible only for small codes – while method='montecarlo' estimates the probability by sampling.

Parameters:

channel (Distribution or object) – The channel: a conditional distribution p(Y|X), or any object implementing the channel protocol (see as_channel()), such as a finite-state channel with memory.
method (str) – 'exact', 'montecarlo', or 'auto' (exact when small).
samples (int) – The number of Monte Carlo samples.
prng (numpy.random.Generator, None) – The random number generator for Monte Carlo sampling.

Returns:

pe – The probability that the decoded message differs from the sent one.

Return type:

float

class SymbolCode(codebook, dist=None, radix=2)[source]

A source code that maps each outcome to a single codeword.

Parameters:

codebook (dict) – A mapping from source outcomes to codeword strings (e.g. '010').
dist (Distribution, None) – The source distribution. Required for the rate-based properties (rate(), average_length(), redundancy(), …).
radix (int) – The size of the code alphabet. Default is 2 (binary).

encode(source)[source]

Encode a sequence of source outcomes into a string of code symbols.

Parameters:: source (iterable) – A sequence of source outcomes.
Returns:: encoded – The concatenated codewords.
Return type:: str

decode(encoded)[source]

Decode a string of code symbols back into source outcomes.

Parameters:: encoded (str) – A concatenation of codewords.
Returns:: source – The decoded sequence of source outcomes.
Return type:: list

average_length()[source]

The expected codeword length, sum_x p(x) * len(codeword(x)).

Returns:: L
Return type:: float

rate()[source]

The expected number of code symbols per source symbol.

For a symbol code this is exactly the average codeword length.

length_variance()[source]

The variance of the codeword length under the source distribution.

Among optimal codes (which share the minimal average length) one often prefers the one of least length variance.

Returns:: var
Return type:: float

kraft_sum()[source]

The Kraft sum sum_x radix ** -len(codeword(x)).

By the Kraft-McMillan inequality this is <= 1 for any uniquely decodable code, with equality iff the code is complete.

Returns:: K
Return type:: float

is_complete()[source]

Whether the code is complete (its Kraft sum equals 1).

Returns:: complete
Return type:: bool

is_prefix_free()[source]

Whether no codeword is a prefix of another (the code is instantaneous).

Returns:: prefix_free
Return type:: bool

is_uniquely_decodable()[source]

Whether the code is uniquely decodable, via the Sardinas-Patterson test.

Returns:: unique
Return type:: bool

is_optimal()[source]

Whether the code achieves the minimal average length (Huffman optimal).

Returns:: optimal
Return type:: bool

class TunstallCode(word_to_code, word_probs, code_length, dist=None, radix=2)[source]

A variable-to-fixed-length source code.

Parameters:

word_to_code (dict) – A mapping from source words (tuples of outcomes) to fixed-length codeword strings.
word_probs (dict) – A mapping from source words to their probabilities (used for the rate).
code_length (int) – The fixed codeword length, in code symbols.
dist (Distribution, None) – The source distribution.
radix (int) – The size of the code alphabet. Default is 2.

expected_word_length()[source]

The expected number of source symbols per dictionary word.

Returns:: L
Return type:: float

rate()[source]

The expected number of code symbols per source symbol.

Each word maps to code_length code symbols and spans expected_word_length source symbols on average.

encode(source)[source]

Encode a sequence of source outcomes by greedily parsing it into words.

Parameters:: source (iterable) – A sequence of source outcomes whose symbols form complete words.
Returns:: encoded
Return type:: str

decode(encoded)[source]

Decode a string of code symbols back into source outcomes.

Parameters:: encoded (str) – A concatenation of fixed-length codewords.
Returns:: source
Return type:: list

class PolarSourceCode(dist, block_length, rv=0, crvs=None, rank_by='entropy', rate=None, size=None, tol=1e-09, max_states=65536)[source]

An exact finite-block polar source code for a small binary source.

The code applies the Arikan transform U^N = X^N G_N to a block of N source bits and stores only the coordinates in the high-entropy set (those that stay nearly uniform given the past and any side information). Decoding fills in the remaining low-entropy coordinates by sequential maximum a posteriori (MAP) decisions against the exact p(U^N, Y^N) table and then inverts the transform.

Because the code enumerates the joint distribution exactly (no density evolution or list decoding), it is intended for small blocks: the number of enumerated states is |copy support|^N, guarded by max_states.

Parameters:

dist (Distribution) – A single-copy source (X) or source-with-side-information (X, Y, ...) distribution. X (indexed by rv) must be binary.
block_length (int) – The block length N; must be a power of two.
rv (int) – The index of the binary source variable. Default is 0.
crvs (list, None) – Indices of the side-information variables available at the decoder. If None, the code is a plain (side-information-free) polar source code.
rank_by (str) – How to rank coordinates for the high-entropy set: "entropy" (default) or "bhattacharyya".
rate (float, None) – A fixed target rate in [0, 1] (keep round(rate * N) coordinates). Mutually exclusive with size.
size (int, None) – The exact high-entropy set size. By default the set is chosen losslessly – every coordinate whose conditional entropy exceeds tol is kept, so the code reconstructs each block exactly. Fixing size (or rate) below the lossless size yields a fixed-rate code that may not be exact.
tol (float) – The conditional-entropy threshold used by the lossless default.
max_states (int) – A guard on the number of enumerated joint outcomes.

property message_length: The number of stored (high-entropy) coordinates per block.

encode(x_block)[source]

Encode a block of N source bits into the high-entropy coordinates.

Parameters:: x_block (sequence) – A length-N block of source outcomes (0/1 or the source’s own binary alphabet symbols).
Returns:: encoded – The transform bits at the high-entropy indices.
Return type:: list of int

decode(encoded, side_information=None)[source]

Decode a block from its high-entropy coordinates (and any side info).

The high-entropy coordinates are filled from encoded; the low-entropy coordinates are recovered one at a time by an exact MAP decision against p(U_i | U^{i-1}, Y^N). The inverse Arikan transform then returns the source bits.

Parameters:

encoded (sequence of int) – The stored high-entropy bits, in ascending index order.
side_information (sequence, None) – The decoder side information Y^N: one symbol per copy when the code was built with crvs (flattened copy-major when multiple side variables). Required iff the code has side information.

Returns:

x_block – The recovered length-N block of source bits.

Return type:

list of int

rate()[source]: The code rate, |high_entropy_set| / block_length.

polar_source(dist, block_length, rv=0, crvs=None, rank_by='entropy', rate=None, size=None, tol=1e-09, max_states=65536)[source]

Build an exact finite-block polar source code.

Parameters:

dist (Distribution) – A single-copy binary source (X) or (X, Y, ...) distribution.
block_length (int) – The block length N; must be a power of two.
rv (int) – The index of the binary source variable. Default is 0.
crvs (list, None) – Indices of the side-information variables available at the decoder.
rank_by (str) – "entropy" (default) or "bhattacharyya".
rate (float, None) – A fixed target rate in [0, 1]. Mutually exclusive with size.
size (int, None) – The high-entropy set size. Defaults to the lossless set (see PolarSourceCode).
tol (float) – The conditional-entropy threshold used by the lossless default.
max_states (int) – A guard on the number of enumerated joint outcomes.

Returns:

code

Return type:

PolarSourceCode

source_bhattacharyya(dist, rv=0, crvs=None)[source]

The source Bhattacharyya parameter of a binary variable.

For a binary X with side information Y,

\[Z(X \mid Y) = 2 \sum_y \sqrt{p(0, y)\,p(1, y)},\]

which reduces to 2 sqrt(p(0) p(1)) when there is no side information. It satisfies 0 <= Z <= 1, with Z near 1 when X is nearly uniform given Y and near 0 when X is nearly determined by Y.

Parameters:

dist (Distribution) – The source distribution.
rv (int) – The index of the binary source variable. Default is 0.
crvs (list, None) – Indices of the side-information variables. If None, Z(X) is computed.

Returns:

Z

Return type:

float

source_polarization_profile(dist, block_length, rv=0, crvs=None, metrics=('entropy', 'bhattacharyya'))[source]

The exact source-polarization profile of N = block_length i.i.d. copies.

Applies the Arikan transform U^N = X^N G_N and reports, for each synthesized coordinate i, the requested metrics conditioned on the past U^{i-1} and the side information Y^N:

"entropy" – the conditional entropy H(U_i | U^{i-1}, Y^N),
"bhattacharyya" – the conditional source Bhattacharyya parameter Z(U_i | U^{i-1}, Y^N).

The conditional entropies sum to N * H(X | Y) (entropy conservation), and as N grows they polarize toward 0 or 1.

A third, optional Goela-style diagnostic reuses dit.divergences.maximum_correlation():

"max_correlation_with_past" – the maximal correlation rho_m(U_i ; U^{i-1}) between each coordinate and its past (0.0 at i = 0, and ignoring side information). Goela et al. (2014) show these also polarize.

Parameters:

dist (Distribution) – The source distribution for a single copy (X) or (X, Y, ...).
block_length (int) – The number of i.i.d. copies N; must be a power of two.
rv (int) – The index of the binary source variable. Default is 0.
crvs (list, None) – Indices of the side-information variables. If None, no side information.
metrics (tuple of str) – Which metrics to compute per coordinate.

Returns:

profile – One dict per coordinate i (in [0, N)), keyed "index" plus the requested metric names.

Return type:

list of dict

source_high_entropy_set(dist, block_length, rate=None, size=None, rv=0, crvs=None, rank_by='entropy', tol=1e-09)[source]

The high-entropy indices selected by a polar source code.

Ranks the N synthesized coordinates by their conditional metric (largest conditional entropy, or largest source Bhattacharyya, given the past and any side information) and returns the top indices – those a lossless polar source code must transmit. At most one of rate or size may be given:

size – keep exactly this many coordinates,
rate – keep round(rate * N) coordinates,
neither – keep every coordinate whose conditional entropy exceeds tol, i.e. the lossless set (the dropped coordinates are deterministic given the past and side information, so the code reconstructs the block exactly).

Parameters:

dist (Distribution) – The source distribution.
block_length (int) – The number of i.i.d. copies N; must be a power of two.
rate (float, None) – The target fraction of coordinates to keep, in [0, 1]. The set size is round(rate * N).
size (int, None) – The exact number of coordinates to keep.
rv (int) – The index of the binary source variable. Default is 0.
crvs (list, None) – Indices of the side-information variables.
rank_by (str) – "entropy" (default) or "bhattacharyya".
tol (float) – The conditional-entropy threshold for the lossless default; coordinates at or below tol are treated as deterministic and dropped.

Returns:

indices – The selected coordinate indices, sorted ascending.

Return type:

list of int

shannon(dist, radix=2)[source]

Build a Shannon code.

Codeword lengths are ceil(log_radix(1/p)) and codewords are read off the base-radix expansion of the cumulative distribution (outcomes sorted by decreasing probability).

Parameters:

dist (Distribution) – The source distribution.
radix (int) – The size of the code alphabet. Default is 2.

Returns:

code

Return type:

SymbolCode

fano(dist, radix=2)[source]

Build a Fano code (the Shannon-Fano top-down splitting method).

Outcomes are sorted by probability and recursively partitioned into two groups of as-equal-as-possible total probability; one bit distinguishes the groups at each level.

Parameters:

dist (Distribution) – The source distribution.
radix (int) – The size of the code alphabet. Only radix=2 is currently supported.

Returns:

code

Return type:

SymbolCode

shannon_fano_elias(dist, radix=2)[source]

Build a Shannon-Fano-Elias code.

Codeword lengths are ceil(log_radix(1/p)) + 1 and codewords are read off the base-radix expansion of the midpoint cumulative distribution F_bar(x) = sum_{y<x} p(y) + p(x)/2. The result is prefix-free.

Parameters:

dist (Distribution) – The source distribution.
radix (int) – The size of the code alphabet. Default is 2.

Returns:

code

Return type:

SymbolCode

huffman(dist, radix=2)[source]

Build a Huffman code, the optimal symbol code for the source.

For radix > 2 the alphabet is padded with zero-probability dummy symbols so that every internal node has exactly radix children.

Parameters:

dist (Distribution) – The source distribution.
radix (int) – The size of the code alphabet. Default is 2.

Returns:

code

Return type:

SymbolCode

length_limited_huffman(dist, max_length, radix=2)[source]

Build an optimal prefix code whose codewords are at most max_length long.

Uses the package-merge algorithm (Larmore & Hirschberg, 1990).

Parameters:

dist (Distribution) – The source distribution.
max_length (int) – The maximum allowed codeword length.
radix (int) – The size of the code alphabet. Only radix=2 is currently supported.

Returns:

code

Return type:

SymbolCode

golomb(dist, m=None, radix=2)[source]

Build a Golomb code over non-negative-integer outcomes.

Each codeword is the unary-coded quotient n // m followed by the truncated-binary remainder n % m. Golomb codes are optimal prefix codes for geometrically distributed sources.

Parameters:

dist (Distribution) – The source distribution; its outcomes must be non-negative integers.
m (int, None) – The Golomb parameter. If None, the parameter optimal for the (assumed geometric) source is chosen from its mean.
radix (int) – The size of the code alphabet. Only radix=2 is currently supported.

Returns:

code

Return type:

SymbolCode

rice(dist, k=None)[source]

Build a Rice code, the Golomb code with parameter m = 2 ** k.

Parameters:

dist (Distribution) – The source distribution; its outcomes must be non-negative integers.
k (int, None) – The Rice parameter. If None, it is chosen from the optimal Golomb parameter for the source.

Returns:

code

Return type:

SymbolCode

tunstall(dist, code_length, radix=2)[source]

Build a Tunstall code for a memoryless source.

Parameters:

dist (Distribution) – The source distribution over single symbols (assumed i.i.d.).
code_length (int) – The fixed codeword length, in code symbols. The dictionary holds up to radix ** code_length words.
radix (int) – The size of the code alphabet. Default is 2.

Returns:

code

Return type:

TunstallCode

universal_code(dist, kind='gamma')[source]

Build a SymbolCode over integer outcomes using a universal code.

Parameters:

dist (Distribution) – The source distribution; its outcomes must be positive integers.
kind (str) – One of 'unary', 'gamma', 'delta', 'omega', or 'fibonacci'.

Returns:

code

Return type:

SymbolCode

unary(n)[source]: The unary codeword for n >= 1: n - 1 ones followed by a zero.

elias_gamma(n)[source]

The Elias gamma codeword for n >= 1.

The codeword is floor(log2(n)) zeros followed by the binary representation of n.

elias_delta(n)[source]

The Elias delta codeword for n >= 1.

The codeword is the gamma code of the bit-length of n followed by the bits of n below its leading one.

elias_omega(n)[source]: The Elias omega (recursive) codeword for n >= 1.

fibonacci(n)[source]: The Fibonacci codeword for n >= 1 (Zeckendorf representation + a 1).

as_channel(channel)[source]

Coerce channel to an object implementing the channel protocol.

Parameters:: channel (Distribution or object) – A conditional distribution p(Y | X); an object with a __dit_channel__() method; or an object already implementing the channel protocol described in dit.coding._channel.
Returns:: channel – An object implementing the channel protocol.
Return type:: object

class MemorylessChannel(dist)[source]

A discrete memoryless channel, applying p(Y | X) independently to each symbol.

Parameters:: dist (Distribution) – A conditional distribution p(Y | X).

log_likelihood(codeword, received)[source]: The log-probability log p(received | codeword).

sample(codeword, prng)[source]: Pass codeword through the channel once.

symbol_log_likelihoods(received)[source]: The (n, len(inputs)) array of log p(y_i | x).

feedback_initial()[source]: The transmitter’s initial belief about the channel; None, since there is no state.

feedback_step(belief, x)[source]: The outcomes (y, p(y | x), None) of one channel use with input x.

capacity()[source]: The channel capacity, in bits per channel use.

bhattacharyya()[source]: The Bhattacharyya parameter of a binary-input channel.

class LinearCode(G, channel=None)[source]

A linear block code over GF(2).

Parameters:

G (array-like) – A k x n binary generator matrix of full row rank.
channel (Distribution, None) – A default channel, as a conditional distribution p(Y|X).

property length: The block length n.

property dimension: The number of information bits k.

property message_length: The number of message bits per codeword (k).

rate()[source]: The code rate k / n.

encode(message)[source]: Encode a length-k message into a length-n codeword.

decode(received, channel=None)[source]

Decode a received word.

With no channel, hard-decision syndrome decoding is used. With a channel, maximum-likelihood decoding over the codebook is used.

codewords()[source]: The list of all codewords.

weight_enumerator()[source]: A Counter mapping each codeword weight to its multiplicity.

minimum_distance()[source]: The minimum distance: the least weight among nonzero codewords.

error_correcting_capability()[source]: The number of errors the code can correct, floor((d - 1) / 2).

property generator_matrix: The generator matrix G.

property parity_check_matrix: The parity-check matrix H.

class LDPCCode(H, channel=None, max_iterations=50)[source]

A linear code decoded by belief propagation on a sparse parity-check matrix.

Parameters:

H (array-like) – An m x n binary parity-check matrix (typically sparse).
channel (Distribution, None) – A default channel.
max_iterations (int) – The maximum number of belief-propagation iterations.

decode(received, channel=None)[source]: Decode by belief propagation when a channel is given, else hard decoding.

class PolarCode(n, k, channel)[source]

A polar code of length n = 2^m with k information bits.

Parameters:

n (int) – The block length; must be a power of two.
k (int) – The number of information bits.
channel (Distribution or object) – The channel used both to select the frozen set and to decode. It must provide a Bhattacharyya parameter; to decode over a channel that does not (e.g. one with memory), design on a memoryless channel and pass the other channel to decode().

property message_length: The number of message bits per codeword (k).

decode(received, channel=None)[source]: Decode by successive cancellation using channel LLRs.

class ConvolutionalCode(generators, message_length, channel=None)[source]

A rate-1/n convolutional code.

Parameters:

generators (sequence of int) – The generator polynomials, in octal (e.g. (0o7, 0o5)).
message_length (int) – The number of information bits per encoded block.
channel (Distribution, None) – A default channel.

rate()[source]: The code rate 1 / n_outputs.

encode(message)[source]: Encode a message, appending K - 1 zeros to terminate the trellis.

decode(received, channel=None)[source]: Decode by the Viterbi algorithm (soft when a channel is given).

repetition(n, channel=None)[source]

The [n, 1, n] repetition code.

Parameters:

n (int) – The block length.
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

LinearCode

parity_check(k, channel=None)[source]

The [k+1, k, 2] single-parity-check code.

Parameters:

k (int) – The number of information bits.
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

LinearCode

hamming(r, channel=None)[source]

The [2^r - 1, 2^r - 1 - r, 3] Hamming code.

Parameters:

r (int) – The number of parity bits (r >= 2).
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

LinearCode

reed_muller(r, m, channel=None)[source]

The Reed-Muller code RM(r, m).

Parameters:

r (int) – The order (0 <= r <= m).
m (int) – The number of variables; the block length is 2^m.
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

LinearCode

golay(extended=False, channel=None)[source]

The binary Golay code.

Parameters:

extended (bool) – If True, return the extended [24, 12, 8] code; otherwise the perfect [23, 12, 7] code.
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

LinearCode

ldpc(H, channel=None, max_iterations=50)[source]

Build an LDPC code from a parity-check matrix.

Parameters:

H (array-like) – An m x n binary parity-check matrix.
channel (Distribution, None) – A default channel.
max_iterations (int) – The maximum number of belief-propagation iterations.

Returns:

code

Return type:

LDPCCode

gallager(n, wc, wr, prng=None, channel=None)[source]

Build a regular Gallager LDPC code.

Parameters:

n (int) – The block length; must be divisible by wr.
wc (int) – The column weight (variable degree).
wr (int) – The row weight (check degree); wc < wr.
prng (numpy.random.Generator, None) – The random number generator used to permute the sub-bands.
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

LDPCCode

polar(n, k, channel)[source]

Build a polar code.

Parameters:

n (int) – The block length; must be a power of two.
k (int) – The number of information bits.
channel (Distribution) – The channel used to select the frozen set and to decode.

Returns:

code

Return type:

PolarCode

convolutional(generators, message_length, channel=None)[source]

Build a rate-1/n convolutional code.

Parameters:

generators (sequence of int) – The generator polynomials in octal (e.g. (0o7, 0o5)).
message_length (int) – The number of information bits per encoded block.
channel (Distribution, None) – A default channel.

Returns:

code

Return type:

ConvolutionalCode

class FeedbackCode(k)[source]

Abstract base class for channel codes with noiseless feedback.

Subclasses implement start(), next_input(), update() and decode(). States must be treated as immutable: update() returns a new state, so that probability_of_error() can branch on every channel output.

Parameters:: k (int) – The number of message bits.

abstractmethod start(channel)[source]: The shared state before the first channel use.

abstractmethod next_input(state, message)[source]: The next channel input for message, or None to stop transmitting.

abstractmethod update(state, y)[source]: The shared state after the channel output y is received (and fed back).

abstractmethod decode(state)[source]: The receiver’s estimate of the message, as a tuple of bits.

transmit(message, channel, prng=None, max_uses=64)[source]

Send message across channel once.

Parameters:

message (sequence of int) – The k message bits.
channel (Distribution or object) – The channel, implementing feedback_initial and feedback_step.
prng (numpy.random.Generator, None) – The random number generator.
max_uses (int) – Stop (and decode) after this many channel uses regardless.

Returns:

decoded (tuple) – The receiver’s estimate of the message.
uses (int) – The number of channel uses.

performance(channel, method='auto', samples=10000, prng=None, max_uses=64)[source]

The block error probability and the expected number of channel uses.

Parameters:

channel (Distribution or object) – The channel, implementing feedback_initial and feedback_step.
method (str) – 'exact' enumerates every message and every output sequence; 'montecarlo' samples; 'auto' is exact unless the output tree exceeds about a million branches.
samples (int) – The number of Monte Carlo samples.
prng (numpy.random.Generator, None) – The random number generator for Monte Carlo sampling.
max_uses (int) – Transmission is cut off (and decoded) after this many channel uses, so the exact values are those of the truncated code.

Returns:

pe (float) – The probability that the decoded message differs from the sent one.
uses (float) – The expected number of channel uses.

probability_of_error(channel, **kwargs)[source]: The block error probability over channel; keyword arguments as in performance().

expected_uses(channel, **kwargs)[source]: The expected number of channel uses; keyword arguments as in performance().

rate(channel, **kwargs)[source]

The rate k / E[N] in bits per channel use, where N is the number of uses.

The rate of a variable-length code depends on the channel, so unlike ChannelCoding.rate it takes one. Keyword arguments are as in performance().

capacity_gap(channel, capacity_kwargs=None, **kwargs)[source]

The feedback capacity of channel minus the code’s rate on it.

The channel’s feedback_capacity() is used when it has one, and otherwise its capacity(), which equals the feedback capacity for a memoryless channel [CT06].

Parameters:

channel (Distribution or object) – The channel.
capacity_kwargs (dict, None) – Passed to the channel’s capacity method.
kwargs (dict) – Passed to rate().

Returns:

gap – capacity - rate, or a (lower, upper) bracket when the channel only brackets its capacity.

Return type:

float or tuple

class ErasureARQ(k, erasure=2)[source]

Send each bit until it arrives unerased.

On the binary erasure channel this is zero-error and uses k / (1 - epsilon) channel uses on average, so it achieves the capacity 1 - epsilon exactly [CT06]. Any output other than the erasure symbol is taken at face value, so on a channel that also flips bits it makes errors.

Parameters:

k (int) – The number of message bits.
erasure (object) – The channel’s erasure symbol.

class StopAndWaitARQ(code, max_rounds=None)[source]

Stop-and-wait automatic repeat request over a linear block code.

The codeword is sent; if the received word is not a codeword (a nonzero syndrome, or any non-binary output such as an erasure) the receiver asks for it again. An error goes undetected only when the error pattern is itself a nonzero codeword. After max_rounds failed attempts the last received word is syndrome-decoded.

Parameters:

code (LinearCode) – The block code used for error detection.
max_rounds (int, None) – The maximum number of transmissions of the block; unlimited by default.

class HorsteinCode(k, n=None, delta=None)[source]

Horstein’s bisection scheme over the 2^k messages.

Both ends keep the receiver’s posterior over the messages, in a fixed order. At each use the transmitter sends 1 if its message lies above the posterior median and 0 otherwise, so every output roughly halves the receiver’s uncertainty [Hor63]. Posterior matching [SF11] generalizes the idea; Schalkwijk and Kailath’s scheme is its Gaussian-channel counterpart [SK66].

The posterior is updated through the channel’s feedback_step, as a joint belief over the message and the channel’s hidden state, so the same code runs on finite-state channels. Inputs must be the bits {0, 1}.

Parameters:

k (int) – The number of message bits.
n (int, None) – Stop after n channel uses.
delta (float, None) – Stop once the largest posterior probability reaches 1 - delta. At least one of n and delta must be given; with both, whichever comes first.

erasure_arq(k, erasure=2)[source]: An ErasureARQ code for k bits: resend each bit until it is not erased.

stop_and_wait_arq(code, max_rounds=None)[source]: A StopAndWaitARQ code: resend code’s block until it is received as a codeword.

horstein(k, n=None, delta=None)[source]: A HorsteinCode for k bits, stopping after n uses or at posterior 1 - delta.