FiniteStateChannel: use Mealy-machine channels with dit channel codes

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`

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 <epsilon_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:

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

docs/references.bib

@@ -1612,4 +1612,14 @@ @book{Howie1995
   title     = {Fundamentals of Semigroup Theory},
   publisher = {Oxford University Press},
   year      = {1995},
-}
+}
+@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},
+}

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",

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,

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",

sofic/generators/capacity/_channel.py

@@ -98,15 +98,23 @@ 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])
     else:
         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:

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)