Feedback over finite-state channels, and the zero-error trapdoor code

CHANGELOG.md

@@ -141,6 +141,16 @@ brute-force reference or a closed-form value and has a regression test.
   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`.
+- `FiniteStateChannel.feedback_step` / `feedback_initial` step the channel one
+  use at a time (the belief is the normalized forward vector) and
+  `feedback_capacity` brackets the capacity with feedback, so dit's feedback
+  codes (`erasure_arq`, `stop_and_wait_arq`, `horstein`) run over a channel
+  with memory. Horstein's scheme on `trapdoor()` reaches rates above the
+  no-feedback upper bound (`0.584`). Needs a dit with `dit.coding.horstein`.
+- `sofic.generators.trapdoor_code.trapdoor_code(k)`: the zero-error feedback
+  code of Permuter et al. for the trapdoor channel (actions relative to the
+  known state, enumerative coding into words with no `11`, backward decoding
+  from output differences), with rate approaching `log2(phi)`.
 - `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

@@ -44,10 +44,65 @@ What the wrapper provides
   capacity *without* feedback (the codes do not use feedback). It requires a
   strongly connected unifilar channel; see :doc:`channel_capacity`.
 
+- **Feedback.** :meth:`FiniteStateChannel.feedback_step` steps the channel one
+  use at a time for dit's feedback codes, with the normalized forward vector
+  (the law of the state given the inputs sent and the outputs fed back) as the
+  transmitter's belief, and :meth:`FiniteStateChannel.feedback_capacity`
+  returns the bracket of :func:`sofic.generators.capacity.feedback_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``.
 
+Codes with feedback
+===================
+
+Feedback cannot raise the capacity of a memoryless channel, but it can raise
+the capacity of a channel with memory. The trapdoor channel's capacity without
+feedback lies in about :math:`[0.572, 0.584]` :cite:`Huleihel2021`, while with
+feedback it is :math:`\log_2 \varphi \approx 0.694` :cite:`Permuter2008`.
+Horstein's bisection scheme :cite:`Horstein1963` (``dit.coding.horstein``)
+tracks the joint posterior over the message and the channel state, so it runs
+over the trapdoor channel unchanged. Stopping once the posterior reaches
+:math:`1 - \delta`, it already beats every code without feedback:
+
+.. code-block:: python
+
+   import numpy as np
+   from dit.coding import horstein
+   from sofic.examples.processes import trapdoor
+
+   code = horstein(10, delta=1e-3)
+   pe, uses = code.performance(trapdoor(), method="montecarlo", samples=300,
+                               prng=np.random.default_rng(1))
+   10 / uses                      # about 0.62 bits per use, with no errors seen
+   code.capacity_gap(trapdoor())  # gap to the feedback capacity, as a bracket
+
+With ``k = 8, 10, 12`` message bits the rate is about ``0.60``, ``0.62`` and
+``0.63``, climbing toward :math:`\log_2 \varphi`. The scheme is not tuned to
+the trapdoor channel: it sends the median split regardless of the state.
+
+``sofic.generators.trapdoor_code.trapdoor_code`` is the scheme of
+:cite:`Permuter2008`, built for the trapdoor channel. The transmitter knows
+the channel's ball from the feedback and chooses whether to send a matching or
+an opposite ball; messages become words of such choices with no two
+"opposite" choices in a row, of which there are about :math:`\varphi^N`. The
+receiver reads the choices backward from whether consecutive outputs differ,
+so the code never errs. Each transmission first flushes the unknown initial
+state, 3.5 uses on average:
+
+.. code-block:: python
+
+   from sofic.generators.trapdoor_code import trapdoor_code
+
+   code = trapdoor_code(4)
+   code.performance(trapdoor(), method="exact", max_uses=70)  # about (0.0, 10.5): N = 7, plus flushing
+   code = trapdoor_code(64)
+   64 / code.expected_uses(trapdoor(), method="montecarlo", max_uses=400)  # about 0.663
+
+The module subclasses ``dit.coding.FeedbackCode``, so it needs a dit with
+feedback codes and is not imported by :mod:`sofic.generators`.
+
 Symbols
 =======
 
@@ -64,4 +119,4 @@ API
 ===
 
 .. autoclass:: FiniteStateChannel
-   :members:
+   :members:

docs/references.bib

@@ -1623,3 +1623,13 @@ @article{Bahl1974
   year    = {1974},
   doi     = {10.1109/TIT.1974.1055186},
 }
+@article{Horstein1963,
+  author  = {Horstein, Michael},
+  title   = {Sequential Transmission Using Noiseless Feedback},
+  journal = {IEEE Transactions on Information Theory},
+  volume  = {9},
+  number  = {3},
+  pages   = {136--143},
+  year    = {1963},
+  doi     = {10.1109/TIT.1963.1057832},
+}

sofic/generators/channel_coding.py

@@ -151,3 +151,35 @@ def capacity(self, **kwargs: Any) -> tuple[float, float]:
 
         bounds = capacity(self.machine, **kwargs)
         return (bounds.lower, bounds.upper)
+
+    def feedback_initial(self) -> np.ndarray:
+        """Return the transmitter's initial belief about the state: the initial law."""
+        return self.initial
+
+    def feedback_step(self, belief: np.ndarray, x: Any) -> list[tuple[Any, float, np.ndarray]]:
+        """Return every ``(y, p(y | x, belief), next_belief)`` of one channel use.
+
+        ``belief`` is the law of the current state given the inputs sent and
+        the outputs fed back so far; the next belief is the normalized forward
+        vector :math:`\\alpha T^{(x, y)} / \\alpha T^{(x, y)} \\mathbf{1}`.
+        """
+        T = self.T[self._in_index[x]]
+        outcomes = []
+        for j, y in enumerate(self.outputs):
+            alpha = belief @ T[j]
+            p = float(alpha.sum())
+            if p > 0.0:
+                outcomes.append((y, p, alpha / p))
+        return outcomes
+
+    def feedback_capacity(self, **kwargs: Any) -> tuple[float, float]:
+        """Bracket the capacity with feedback, as ``(lower, upper)`` bits per use.
+
+        Delegates to :func:`sofic.generators.capacity.feedback_capacity`, which
+        needs a strongly connected unifilar channel; keyword arguments are
+        passed on.
+        """
+        from sofic.generators.capacity import feedback_capacity
+
+        bounds = feedback_capacity(self.machine, **kwargs)
+        return (bounds.lower, bounds.upper)

sofic/generators/trapdoor_code.py

@@ -0,0 +1,143 @@
+"""A zero-error feedback code for the trapdoor channel.
+
+Permuter, Cuff, Van Roy and Weissman :cite:`Permuter2008` give a simple scheme
+that communicates through the trapdoor channel with feedback at rates
+approaching its feedback capacity :math:`\\log_2 \\varphi`, with no errors at
+all. This module implements it as a :class:`dit.coding.FeedbackCode`, so it
+needs a dit with feedback codes and is not imported by :mod:`sofic.generators`.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+from typing import Any
+
+from dit.coding.feedback import FeedbackCode  # ty: ignore[unresolved-import]
+from dit.exceptions import ditException
+
+__all__ = ["TrapdoorCode", "trapdoor_code"]
+
+
+def _counts(length: int) -> list[int]:
+    """``counts[r]`` is the number of binary words of length ``r`` with no ``11``."""
+    counts = [1, 2]
+    while len(counts) <= length:
+        counts.append(counts[-1] + counts[-2])
+    return counts
+
+
+def rll_encode(index: int, length: int) -> tuple[int, ...]:
+    """The `index`-th (lexicographic) binary word of `length` with no two adjacent ``1`` s."""
+    counts = _counts(length)
+    if not 0 <= index < counts[length]:
+        raise ValueError(f"index must lie in [0, {counts[length]})")
+    word, previous = [], 0
+    for i in range(length):
+        remaining = length - i - 1
+        if previous == 0 and index >= counts[remaining]:
+            index -= counts[remaining]
+            previous = 1
+        else:
+            previous = 0
+        word.append(previous)
+    return tuple(word)
+
+
+def rll_decode(word: Sequence[int]) -> int:
+    """The lexicographic index of a binary word with no two adjacent ``1`` s; inverse of :func:`rll_encode`."""
+    counts = _counts(len(word))
+    index = 0
+    for i, bit in enumerate(word):
+        if bit:
+            index += counts[len(word) - i - 1]
+    return index
+
+
+class TrapdoorCode(FeedbackCode):
+    """The capacity-achieving zero-error feedback code for the trapdoor channel :cite:`Permuter2008`.
+
+    The transmitter works in *actions*: ``0`` sends a ball of the same type as
+    the one in the channel, ``1`` the opposite type. It knows the channel's
+    state from the feedback. A message is mapped by enumerative coding to an
+    action word of length ``N`` that ends in ``0`` and has no two ``1`` s in a
+    row; there are about :math:`\\varphi^N` such words, so the rate approaches
+    :math:`\\log_2 \\varphi`. The receiver recovers the actions backward from
+    the last one, using only whether consecutive outputs differ.
+
+    The initial state is unknown, so each transmission first *flushes* the
+    channel with the inputs ``0, 1, 0, 1, ...`` until an output differs from
+    its input; from then on both ends know the state. Flushing takes 3.5 uses
+    on average from a uniform state, and is counted in the rate.
+
+    Parameters
+    ----------
+    k : int
+        The number of message bits.
+    """
+
+    def __init__(self, k: int) -> None:
+        super().__init__(k)
+        length = 0
+        while _counts(length)[length] < 2**k:
+            length += 1
+        self.length = length + 1
+
+    def __repr__(self) -> str:
+        return f"TrapdoorCode(k={self.k}, N={self.length})"
+
+    def start(self, channel: Any) -> tuple[int, ...]:
+        """The outputs received so far: none."""
+        if tuple(channel.inputs) != (0, 1) or tuple(channel.outputs) != (0, 1):
+            raise ditException("The trapdoor code needs a channel with inputs and outputs {0, 1}.")
+        return ()
+
+    @staticmethod
+    def _flushed(received: tuple[int, ...]) -> int | None:
+        """The number of flushing uses, or ``None`` if the channel is still being flushed."""
+        for t, y in enumerate(received):
+            if y != t % 2:
+                return t + 1
+        return None
+
+    def _actions(self, message: Sequence[int]) -> tuple[int, ...]:
+        index = int("".join(map(str, message)), 2)
+        return (*rll_encode(index, self.length - 1), 0)
+
+    def next_input(self, state: tuple[int, ...], message: Sequence[int]) -> int | None:
+        """The next flushing input, or the next action applied to the current state."""
+        flushed = self._flushed(state)
+        if flushed is None:
+            return len(state) % 2
+        data = state[flushed:]
+        if len(data) == self.length:
+            return None
+        actions = self._actions(message)
+        ball = (flushed - 1) % 2
+        for action, y in zip(actions, data, strict=False):
+            x = action ^ ball
+            ball ^= x ^ y
+        return actions[len(data)] ^ ball
+
+    def update(self, state: tuple[int, ...], y: int) -> tuple[int, ...]:
+        """Append the output `y`."""
+        return (*state, y)
+
+    def decode(self, state: tuple[int, ...]) -> tuple[int, ...]:
+        """Recover the actions backward from the last, then invert the enumerative code."""
+        flushed = self._flushed(state)
+        data = state[flushed:] if flushed is not None else ()
+        if len(data) < self.length:
+            return (0,) * self.k
+        actions = [0] * self.length
+        for t in range(self.length - 2, -1, -1):
+            differs = data[t + 1] ^ data[t]
+            actions[t] = int(differs and not actions[t + 1])
+        index = rll_decode(actions[:-1])
+        if index >= 2**self.k:
+            return (0,) * self.k
+        return tuple(int(b) for b in format(index, f"0{self.k}b"))
+
+
+def trapdoor_code(k: int) -> TrapdoorCode:
+    """A :class:`TrapdoorCode` for `k` message bits."""
+    return TrapdoorCode(k)