From 21cd061c6f6d1aca200c9a9a91ab1ea9337c68b8 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Mon, 26 Oct 2026 11:57:36 -0600 Subject: [PATCH] Add FiniteStateChannel so Mealy-machine channels attach to dit codes MealyMachine.__dit_channel__ wraps a stochastic transducer as a FiniteStateChannel implementing dit.coding's channel protocol: forward-algorithm likelihoods, state-path sampling, BCJR symbol likelihoods, and a no-feedback capacity bracket. Factor the initial-state law out of channel_tensors. Co-authored-by: Cursor --- CHANGELOG.md | 6 + docs/generators/channel_coding.rst | 67 +++++++++++ docs/generators/generators.rst | 1 + docs/references.bib | 12 +- sofic/__init__.py | 2 + sofic/automata/transducers.py | 6 + sofic/generators/__init__.py | 2 + sofic/generators/capacity/_channel.py | 12 +- sofic/generators/channel_coding.py | 153 ++++++++++++++++++++++++++ tests/test_channel_coding.py | 150 +++++++++++++++++++++++++ 10 files changed, 408 insertions(+), 3 deletions(-) create mode 100644 docs/generators/channel_coding.rst create mode 100644 sofic/generators/channel_coding.py create mode 100644 tests/test_channel_coding.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 652a406..a133805 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -135,6 +135,12 @@ brute-force reference or a closed-form value and has a regression test. ### New +- `FiniteStateChannel`: a stochastic `MealyMachine` (or ε-transducer) as a + channel for `dit.coding` channel codes. `MealyMachine.__dit_channel__` returns + it, so `hamming(3).probability_of_error(trapdoor())` works directly. The + likelihood comes from the forward algorithm, the soft decoders get BCJR symbol + likelihoods, and `capacity_gap` returns a `(lower, upper)` bracket from the + no-feedback capacity bounds. Needs a dit with `dit.coding.as_channel`. - `sofic.algebra`: algebraic automata theory for deterministic presentations. `TransformationSemigroup` (lazy closure, shortest words, sink completion of partial maps), Green's relations and maximal subgroups, `syntactic_monoid` diff --git a/docs/generators/channel_coding.rst b/docs/generators/channel_coding.rst new file mode 100644 index 0000000..a1909e7 --- /dev/null +++ b/docs/generators/channel_coding.rst @@ -0,0 +1,67 @@ +.. channel_coding.rst +.. py:currentmodule:: sofic.generators + +**************************************** +Channel Coding over Channels with Memory +**************************************** + +The channel codes in :mod:`dit.coding` (linear block, LDPC, polar and +convolutional codes) evaluate themselves against any object implementing dit's +channel protocol. A stochastic :class:`~sofic.automata.transducers.MealyMachine` +-- including every :doc:`ε-transducer ` -- with transitions +:math:`T(y, s' \mid s, x)` is a finite-state channel, and its +``__dit_channel__`` method wraps it as a :class:`FiniteStateChannel`. So the +machine can be handed to a dit code as is: + +.. code-block:: python + + from dit.coding import hamming + from sofic.examples.processes import trapdoor + + code = hamming(3) + code.decode((0, 1, 1, 0, 1, 0, 0), channel=trapdoor()) + code.probability_of_error(trapdoor(), method="exact") + code.capacity_gap(trapdoor()) # (lower, upper) + +What the wrapper provides +========================= + +- **Likelihoods.** :meth:`FiniteStateChannel.log_likelihood` is the forward + algorithm :math:`\log \pi^\top \prod_t T^{(x_t, y_t)} \mathbf{1}`, exact for + any finite-state channel, unifilar or not. Maximum-likelihood decoding and the + exact block-error probability are therefore exact over a channel with memory. +- **Sampling.** :meth:`FiniteStateChannel.sample` draws the state path from the + initial law, for Monte Carlo error estimates. +- **Soft information.** :meth:`FiniteStateChannel.symbol_log_likelihoods` runs + the forward-backward (BCJR) recursions :cite:`Bahl1974` with the inputs at + the other positions uniform and i.i.d., giving :math:`\log p(y^n \mid x_i)`. + Belief propagation, successive cancellation and soft Viterbi decode from + these per-position metrics. This is separate detection and decoding: it + ignores the correlations the channel induces between positions, so those + decoders are approximate over a channel with memory. +- **Capacity.** :meth:`FiniteStateChannel.capacity` returns the + ``(lower, upper)`` bracket of :func:`sofic.generators.capacity.capacity`, the + capacity *without* feedback (the codes do not use feedback). It requires a + strongly connected unifilar channel; see :doc:`channel_capacity`. + +Polar codes are designed from a channel's Bhattacharyya parameter, which a +channel with memory does not have. Design the polar code on a memoryless +channel and pass the finite-state channel to ``decode``. + +Symbols +======= + +dit's binary codes use the bits ``0`` and ``1``, while the example channels are +labeled with the strings ``"0"`` and ``"1"``. By default the wrapper maps digit +strings to integers; pass ``input_map`` and ``output_map`` to +:class:`FiniteStateChannel` for any other relabeling. + +A one-state machine is a memoryless channel: ``binary_channel(p, p)`` gives the +same likelihoods, soft information and error probabilities as dit's +``binary_symmetric_channel(p)``. + +API +=== + +.. autoclass:: FiniteStateChannel + :members: diff --git a/docs/generators/generators.rst b/docs/generators/generators.rst index 13bbdfb..54129b3 100644 --- a/docs/generators/generators.rst +++ b/docs/generators/generators.rst @@ -30,6 +30,7 @@ Computational mechanics epsilon_machine epsilon_transducer channel_capacity + channel_coding block_convergence bidirectional_epsilon_machine mixed_state_presentation diff --git a/docs/references.bib b/docs/references.bib index 049780c..006cd68 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1612,4 +1612,14 @@ @book{Howie1995 title = {Fundamentals of Semigroup Theory}, publisher = {Oxford University Press}, year = {1995}, -} \ No newline at end of file +} +@article{Bahl1974, + author = {Bahl, Lalit R. and Cocke, John and Jelinek, Frederick and Raviv, Josef}, + title = {Optimal Decoding of Linear Codes for Minimizing Symbol Error Rate}, + journal = {IEEE Transactions on Information Theory}, + volume = {20}, + number = {2}, + pages = {284--287}, + year = {1974}, + doi = {10.1109/TIT.1974.1055186}, +} diff --git a/sofic/__init__.py b/sofic/__init__.py index 2443353..6b86cf4 100644 --- a/sofic/__init__.py +++ b/sofic/__init__.py @@ -60,6 +60,7 @@ BlockEntropyEstimates, EpsilonMachine, EpsilonTransducer, + FiniteStateChannel, FunctionalGenerativeModel, GacsKornerGenerativeModel, HiddenMarkovModel, @@ -120,6 +121,7 @@ "EPSILON", "EpsilonMachine", "EpsilonTransducer", + "FiniteStateChannel", "FunctionalGenerativeModel", "GacsKornerGenerativeModel", "HiddenMarkovModel", diff --git a/sofic/automata/transducers.py b/sofic/automata/transducers.py index 817a788..911de88 100644 --- a/sofic/automata/transducers.py +++ b/sofic/automata/transducers.py @@ -219,6 +219,12 @@ def to_sliding_block_code(self, *, max_window: int = 4) -> Any: return TextileSystem.from_transducer(self).induced_code(max_window=max_window) + def __dit_channel__(self) -> Any: + """Return this stochastic transducer as a channel for :mod:`dit.coding` codes.""" + from sofic.generators.channel_coding import FiniteStateChannel + + return FiniteStateChannel(self) + def labeled_transition_matrices( self, node_ordering: Sequence[Hashable] | None = None, diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py index db8de0c..7eb8cac 100644 --- a/sofic/generators/__init__.py +++ b/sofic/generators/__init__.py @@ -4,6 +4,7 @@ from sofic.generators.bidirectional_epsilon_machine import BidirectionalEpsilonMachine from sofic.generators.block_convergence import BlockConvergenceDiagram, BlockConvergenceEstimates from sofic.generators.block_entropy import BlockEntropyDiagram, BlockEntropyEstimates +from sofic.generators.channel_coding import FiniteStateChannel from sofic.generators.channel_measures import ( channel_statistical_complexity, driven_entropy_rate, @@ -75,6 +76,7 @@ "BlockConvergenceEstimates", "EpsilonMachine", "EpsilonTransducer", + "FiniteStateChannel", "LumpabilityError", "channel_statistical_complexity", "driven_entropy_rate", diff --git a/sofic/generators/capacity/_channel.py b/sofic/generators/capacity/_channel.py index 4dd4e4c..d4f1bb8 100644 --- a/sofic/generators/capacity/_channel.py +++ b/sofic/generators/capacity/_channel.py @@ -98,6 +98,15 @@ def channel_tensors(channel: MealyMachine) -> ChannelTensors: if not allowed.any(axis=1).all(): raise StochasticValidationError("every state needs at least one admissible input") + return ChannelTensors(states, inputs, outputs, P, F, allowed, initial_law(channel, states)) + + +def initial_law(channel: MealyMachine, states: tuple[Hashable, ...] | list[Hashable]) -> np.ndarray: + """Return the initial state law of ``channel`` over ``states``. + + An explicit ``initial_distribution`` wins; otherwise the law is uniform on + ``initial_states``, or on every state when there are none. + """ distribution = getattr(channel, "initial_distribution", None) or {} if distribution: initial = np.array([float(distribution.get(state, 0.0)) for state in states]) @@ -105,8 +114,7 @@ def channel_tensors(channel: MealyMachine) -> ChannelTensors: initial = np.array([1.0 if state in channel.initial_states else 0.0 for state in states]) if initial.sum() == 0.0: initial = np.ones(len(states)) - initial = initial / initial.sum() - return ChannelTensors(states, inputs, outputs, P, F, allowed, initial) + return initial / initial.sum() def require_strongly_connected(tensors: ChannelTensors) -> None: diff --git a/sofic/generators/channel_coding.py b/sofic/generators/channel_coding.py new file mode 100644 index 0000000..b24d255 --- /dev/null +++ b/sofic/generators/channel_coding.py @@ -0,0 +1,153 @@ +"""Finite-state channels as channels for :mod:`dit.coding` channel codes. + +A stochastic :class:`~sofic.automata.transducers.MealyMachine` with +transitions ``T(y, s' | s, x)`` is a finite-state channel. Wrapped as a +:class:`FiniteStateChannel` it implements the channel protocol of +:func:`dit.coding.as_channel`, so a dit code can be decoded and evaluated over +it. ``MealyMachine.__dit_channel__`` returns this wrapper, so the machine can be +passed to dit directly. +""" + +from __future__ import annotations + +from collections.abc import Hashable, Mapping, Sequence +from typing import TYPE_CHECKING, Any + +import numpy as np + +from sofic.generators.capacity._channel import initial_law + +if TYPE_CHECKING: + from sofic.automata.transducers import MealyMachine + +__all__ = ["FiniteStateChannel"] + + +def _sorted(values: Any) -> tuple[Any, ...]: + try: + return tuple(sorted(values)) + except TypeError: + return tuple(sorted(values, key=repr)) + + +def _default_map(symbols: Sequence[Any]) -> dict[Any, Any]: + """Map digit strings to ``int`` (``"0" -> 0``), leaving other symbols alone.""" + if all(isinstance(s, str) and s.isdigit() for s in symbols): + mapping = {s: int(s) for s in symbols} + if len(set(mapping.values())) == len(mapping): + return mapping + return {s: s for s in symbols} + + +class FiniteStateChannel: + """A finite-state channel ``T(y, s' | s, x)`` in the :mod:`dit.coding` channel protocol. + + Parameters + ---------- + machine : MealyMachine + A stochastic Mealy machine whose transitions read ``x`` and emit ``y`` + with probability ``T(y, s' | s, x)``. Unifilarity is not required, + except by :meth:`capacity`. + input_map, output_map : Mapping, optional + Maps from the machine's symbols to the symbols dit sees. By default, + digit strings become integers, so ``"0"`` and ``"1"`` become the bits + ``0`` and ``1`` the binary codes use. + + Notes + ----- + :meth:`log_likelihood` is the forward algorithm, exact for any + finite-state channel. :meth:`symbol_log_likelihoods` runs the + forward-backward (BCJR) recursions :cite:`Bahl1974` with uniform i.i.d. + inputs at the other positions, giving ``log p(y^n | x_i)``; the per-symbol + decoders in dit use these as their soft information. + """ + + def __init__( + self, + machine: MealyMachine, + input_map: Mapping[Any, Any] | None = None, + output_map: Mapping[Any, Any] | None = None, + ) -> None: + matrices, states = machine.labeled_transition_matrices() + raw_inputs = _sorted({x for x, _ in matrices}) + raw_outputs = _sorted({y for _, y in matrices}) + input_map = dict(input_map) if input_map is not None else _default_map(raw_inputs) + output_map = dict(output_map) if output_map is not None else _default_map(raw_outputs) + + self.machine = machine + self.states: tuple[Hashable, ...] = tuple(states) + self.inputs = _sorted({input_map[x] for x in raw_inputs}) + self.outputs = _sorted({output_map[y] for y in raw_outputs}) + self._in_index = {x: i for i, x in enumerate(self.inputs)} + self._out_index = {y: j for j, y in enumerate(self.outputs)} + + n = len(self.states) + self.T = np.zeros((len(self.inputs), len(self.outputs), n, n)) + for (x, y), matrix in matrices.items(): + self.T[self._in_index[input_map[x]], self._out_index[output_map[y]]] += matrix + self.initial = initial_law(machine, self.states) + + def log_likelihood(self, codeword: Sequence[Any], received: Sequence[Any]) -> float: + """Return ``log p(received | codeword)`` by the scaled forward algorithm.""" + alpha = self.initial + total = 0.0 + for x, y in zip(codeword, received, strict=True): + alpha = alpha @ self.T[self._in_index[x], self._out_index[y]] + scale = alpha.sum() + if scale <= 0.0: + return -np.inf + total += np.log(scale) + alpha = alpha / scale + return float(total) + + def sample(self, codeword: Sequence[Any], prng: np.random.Generator) -> tuple[Any, ...]: + """Pass ``codeword`` through the channel once, starting from the initial law.""" + n = len(self.states) + state = prng.choice(n, p=self.initial) + received = [] + for x in codeword: + weights = self.T[self._in_index[x], :, state, :].ravel() + total = weights.sum() + if total <= 0.0: + raise ValueError(f"input {x!r} is not admissible in state {self.states[state]!r}") + y, state = divmod(prng.choice(weights.size, p=weights / total), n) + received.append(self.outputs[y]) + return tuple(received) + + def symbol_log_likelihoods(self, received: Sequence[Any]) -> np.ndarray: + """Return the ``(n, len(inputs))`` array of ``log p(y^n | x_i = a)``, up to row constants. + + The inputs at the other positions are taken to be uniform and i.i.d. + """ + cols = [self._out_index[y] for y in received] + n = len(cols) + M = [self.T[:, j].mean(axis=0) for j in cols] + + alphas = [self.initial] + for t in range(n): + alpha = alphas[-1] @ M[t] + scale = alpha.sum() + alphas.append(alpha / scale if scale > 0.0 else alpha) + betas = [np.ones(len(self.states))] + for t in reversed(range(n)): + beta = M[t] @ betas[-1] + scale = beta.sum() + betas.append(beta / scale if scale > 0.0 else beta) + betas.reverse() + + S = np.empty((n, len(self.inputs))) + with np.errstate(divide="ignore"): + for t, j in enumerate(cols): + S[t] = np.log(np.einsum("s,ast,t->a", alphas[t], self.T[:, j], betas[t + 1])) + return S + + def capacity(self, **kwargs: Any) -> tuple[float, float]: + """Bracket the capacity without feedback, as ``(lower, upper)`` bits per use. + + Delegates to :func:`sofic.generators.capacity.capacity`, which needs a + strongly connected unifilar channel; keyword arguments are passed on. + """ + from sofic.generators.capacity import capacity + + bounds = capacity(self.machine, **kwargs) + return (bounds.lower, bounds.upper) diff --git a/tests/test_channel_coding.py b/tests/test_channel_coding.py new file mode 100644 index 0000000..5233bf2 --- /dev/null +++ b/tests/test_channel_coding.py @@ -0,0 +1,150 @@ +"""Finite-state channels attached to dit channel codes.""" + +import itertools +import math + +import numpy as np +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.examples.processes import binary_channel, ising_channel, trapdoor +from sofic.generators import FiniteStateChannel + +dit_coding = pytest.importorskip("dit.coding") +if not hasattr(dit_coding, "as_channel"): + pytest.skip("dit.coding has no channel protocol", allow_module_level=True) + +from dit.coding import as_channel, binary_symmetric_channel, hamming, repetition # noqa: E402 +from dit.coding._channel import symbol_llrs # noqa: E402 + +CHANNELS = [trapdoor, ising_channel] + + +def path_sum(machine, codeword, received): + """``p(received | codeword)`` summed over state paths, from the raw machine.""" + matrices, states = machine.labeled_transition_matrices() + from sofic.generators.capacity._channel import initial_law + + pi = initial_law(machine, states) + total = 0.0 + for path in itertools.product(range(len(states)), repeat=len(codeword) + 1): + prob = pi[path[0]] + for t, (x, y) in enumerate(zip(codeword, received, strict=True)): + matrix = matrices.get((str(x), str(y))) + prob *= 0.0 if matrix is None else matrix[path[t], path[t + 1]] + total += prob + return total + + +words = st.integers(1, 5).flatmap( + lambda n: st.tuples( + st.lists(st.sampled_from((0, 1)), min_size=n, max_size=n), + st.lists(st.sampled_from((0, 1)), min_size=n, max_size=n), + ) +) + + +def test_mealy_machine_is_a_dit_channel(): + channel = as_channel(trapdoor()) + assert isinstance(channel, FiniteStateChannel) + assert channel.inputs == (0, 1) + assert channel.outputs == (0, 1) + + +@pytest.mark.parametrize("factory", CHANNELS) +@settings(max_examples=40, deadline=None) +@given(pair=words) +def test_likelihood_matches_path_sum(factory, pair): + machine = factory() + codeword, received = pair + expected = path_sum(machine, codeword, received) + value = FiniteStateChannel(machine).log_likelihood(codeword, received) + if expected == 0.0: + assert value == -math.inf + else: + assert math.exp(value) == pytest.approx(expected) + + +@pytest.mark.parametrize("factory", CHANNELS) +@settings(max_examples=20, deadline=None) +@given(received=st.lists(st.sampled_from((0, 1)), min_size=1, max_size=5)) +def test_symbol_likelihoods_match_marginal(factory, received): + """BCJR symbol likelihoods equal ``p(y^n | x_i)`` under uniform inputs, up to row constants.""" + machine = factory() + n = len(received) + S = FiniteStateChannel(machine).symbol_log_likelihoods(received) + for i in range(n): + marginal = [0.0, 0.0] + for codeword in itertools.product((0, 1), repeat=n): + marginal[codeword[i]] += path_sum(machine, codeword, received) + with np.errstate(divide="ignore"): + expected = np.log(marginal) + finite = np.isfinite(expected) + assert np.array_equal(np.isfinite(S[i]), finite) + if finite.all(): + assert S[i, 0] - S[i, 1] == pytest.approx(expected[0] - expected[1]) + + +@settings(max_examples=25, deadline=None) +@given( + p=st.floats(0.01, 0.45), + received=st.lists(st.sampled_from((0, 1)), min_size=7, max_size=7), + codeword=st.lists(st.sampled_from((0, 1)), min_size=7, max_size=7), +) +def test_one_state_channel_matches_bsc(p, received, codeword): + """A one-state binary channel is the binary symmetric channel.""" + stateful = as_channel(binary_channel(p, p)) + bsc = as_channel(binary_symmetric_channel(p)) + assert stateful.log_likelihood(codeword, received) == pytest.approx(bsc.log_likelihood(codeword, received)) + assert symbol_llrs(stateful, received) == pytest.approx(symbol_llrs(bsc, received)) + + +@pytest.mark.parametrize("p", [0.05, 0.2]) +def test_one_state_error_probability_matches_bsc(p): + code = hamming(3) + expected = code.probability_of_error(binary_symmetric_channel(p), method="exact") + assert code.probability_of_error(binary_channel(p, p), method="exact") == pytest.approx(expected) + + +@pytest.mark.parametrize("factory", CHANNELS) +def test_sample_frequencies_match_likelihood(factory): + channel = FiniteStateChannel(factory()) + codeword = (0, 1, 1) + prng = np.random.default_rng(0) + draws = [channel.sample(codeword, prng) for _ in range(4000)] + for received in itertools.product((0, 1), repeat=3): + frequency = sum(d == received for d in draws) / len(draws) + expected = math.exp(channel.log_likelihood(codeword, received)) + assert frequency == pytest.approx(expected, abs=0.03) + + +@pytest.mark.parametrize("factory", CHANNELS) +def test_exact_error_probability_over_memory(factory): + """The exact block-error probability is the direct sum over messages and received words.""" + machine = factory() + code = repetition(3) + expected = 0.0 + for message in [(0,), (1,)]: + for received in itertools.product((0, 1), repeat=3): + if tuple(code.decode(received, channel=machine)) != message: + expected += path_sum(machine, code.encode(message), received) / 2 + assert code.probability_of_error(machine, method="exact") == pytest.approx(expected) + + +def test_capacity_gap_is_a_bracket(): + p = 0.1 + code = hamming(3) + lower, upper = code.capacity_gap(binary_channel(p, p)) + bsc_gap = code.capacity_gap(binary_symmetric_channel(p)) + assert lower <= bsc_gap + 1e-6 + assert bsc_gap <= upper + 1e-6 + + +def test_symbol_maps(): + channel = FiniteStateChannel(trapdoor(), input_map={"0": "a", "1": "b"}) + assert channel.inputs == ("a", "b") + assert channel.outputs == (0, 1) + # Ball 0 left (prob 1/2): emit 0 surely, then emit the old 0 with prob 1/2. + # Ball 1 left: after inserting 0 and emitting 0, two 1s remain, so 0 is impossible. + assert channel.log_likelihood(("a", "b"), (0, 0)) == pytest.approx(math.log(0.25))