From ebc140257386a5936b4ccbdce0ed53177afd4053 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Mon, 26 Oct 2026 12:40:18 -0600 Subject: [PATCH 1/4] Step FiniteStateChannel one use at a time for dit's feedback codes feedback_initial / feedback_step use the normalized forward vector as the transmitter's belief, and feedback_capacity brackets the capacity with feedback. Horstein's scheme on the trapdoor channel beats the no-feedback upper bound. Co-authored-by: Cursor --- CHANGELOG.md | 6 +++ docs/generators/channel_coding.rst | 34 +++++++++++++ docs/references.bib | 10 ++++ sofic/generators/channel_coding.py | 32 +++++++++++++ tests/test_channel_feedback.py | 76 ++++++++++++++++++++++++++++++ 5 files changed, 158 insertions(+) create mode 100644 tests/test_channel_feedback.py diff --git a/CHANGELOG.md b/CHANGELOG.md index a133805..f9f7aed 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -141,6 +141,12 @@ 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.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 index a1909e7..f2fccfc 100644 --- a/docs/generators/channel_coding.rst +++ b/docs/generators/channel_coding.rst @@ -44,10 +44,44 @@ 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. + Symbols ======= diff --git a/docs/references.bib b/docs/references.bib index 006cd68..a05dad1 100644 --- a/docs/references.bib +++ b/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}, +} diff --git a/sofic/generators/channel_coding.py b/sofic/generators/channel_coding.py index b24d255..04ca153 100644 --- a/sofic/generators/channel_coding.py +++ b/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) diff --git a/tests/test_channel_feedback.py b/tests/test_channel_feedback.py new file mode 100644 index 0000000..1864dd9 --- /dev/null +++ b/tests/test_channel_feedback.py @@ -0,0 +1,76 @@ +"""Finite-state channels driven one use at a time by dit's feedback codes.""" + +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, "horstein"): + pytest.skip("dit.coding has no feedback codes", allow_module_level=True) + +from dit.coding import binary_symmetric_channel, horstein # noqa: E402 + +CHANNELS = [trapdoor, ising_channel] +LOG2_PHI = math.log2((1 + math.sqrt(5)) / 2) +NO_FEEDBACK_UPPER = 0.584 # Huleihel et al. (2021) + +words = st.integers(1, 6).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), + ) +) + + +@pytest.mark.parametrize("factory", CHANNELS) +@settings(max_examples=40, deadline=None) +@given(pair=words) +def test_feedback_steps_reproduce_log_likelihood(factory, pair): + channel = FiniteStateChannel(factory()) + codeword, received = pair + belief, total = channel.feedback_initial(), 0.0 + for x, y in zip(codeword, received, strict=True): + outcomes = {out: (p, b) for out, p, b in channel.feedback_step(belief, x)} + assert sum(p for p, _ in outcomes.values()) == pytest.approx(1.0) + if y not in outcomes: + total = -math.inf + break + p, belief = outcomes[y] + assert belief.sum() == pytest.approx(1.0) + total += math.log(p) + assert total == pytest.approx(channel.log_likelihood(codeword, received)) + + +def test_one_state_channel_matches_bsc_under_feedback(): + code = horstein(3, n=7) + expected = code.probability_of_error(binary_symmetric_channel(0.1), method="exact") + assert code.probability_of_error(binary_channel(0.1, 0.1), method="exact") == pytest.approx(expected) + + +def test_horstein_exact_matches_montecarlo_on_trapdoor(): + code = horstein(2, n=6) + exact = code.probability_of_error(trapdoor(), method="exact") + sampled = code.probability_of_error(trapdoor(), method="montecarlo", samples=6000, prng=np.random.default_rng(0)) + assert sampled == pytest.approx(exact, abs=0.015) + + +def test_feedback_capacity_brackets_log2_phi(): + lower, upper = FiniteStateChannel(trapdoor()).feedback_capacity() + assert lower - 1e-6 <= LOG2_PHI <= upper + 1e-6 + + +def test_horstein_beats_the_no_feedback_capacity_of_the_trapdoor(): + """Feedback lifts the trapdoor's rate above every no-feedback code, but not above log2(phi).""" + code = horstein(8, delta=1e-3) + pe, uses = code.performance(trapdoor(), method="montecarlo", samples=400, prng=np.random.default_rng(1)) + rate = 8 / uses + assert pe <= 0.01 + assert NO_FEEDBACK_UPPER < rate < LOG2_PHI + lower, upper = code.capacity_gap(trapdoor(), method="montecarlo", samples=100, prng=np.random.default_rng(2)) + assert 0 < lower <= upper From e58bd657717afe4b0abc7defa44f51bea57dda1f Mon Sep 17 00:00:00 2001 From: Ryan James Date: Mon, 26 Oct 2026 12:42:46 -0600 Subject: [PATCH 2/4] Add the zero-error trapdoor feedback code of Permuter et al. Actions relative to the known ball, enumerative coding into words with no two 1s, and backward decoding from output differences; the rate approaches log2(phi), the trapdoor's feedback capacity. Co-authored-by: Cursor --- CHANGELOG.md | 4 + docs/generators/channel_coding.rst | 24 +++++ sofic/generators/trapdoor_code.py | 143 +++++++++++++++++++++++++++++ tests/test_channel_feedback.py | 42 +++++++++ 4 files changed, 213 insertions(+) create mode 100644 sofic/generators/trapdoor_code.py diff --git a/CHANGELOG.md b/CHANGELOG.md index f9f7aed..556b653 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -147,6 +147,10 @@ brute-force reference or a closed-form value and has a regression test. 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` diff --git a/docs/generators/channel_coding.rst b/docs/generators/channel_coding.rst index f2fccfc..8aa945c 100644 --- a/docs/generators/channel_coding.rst +++ b/docs/generators/channel_coding.rst @@ -82,6 +82,27 @@ 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. +:func:`~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 ======= @@ -99,3 +120,6 @@ API .. autoclass:: FiniteStateChannel :members: + +.. automodule:: sofic.generators.trapdoor_code + :members: TrapdoorCode, trapdoor_code, rll_encode, rll_decode diff --git a/sofic/generators/trapdoor_code.py b/sofic/generators/trapdoor_code.py new file mode 100644 index 0000000..8bfa0b3 --- /dev/null +++ b/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 +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) diff --git a/tests/test_channel_feedback.py b/tests/test_channel_feedback.py index 1864dd9..b269be1 100644 --- a/tests/test_channel_feedback.py +++ b/tests/test_channel_feedback.py @@ -16,6 +16,8 @@ from dit.coding import binary_symmetric_channel, horstein # noqa: E402 +from sofic.generators.trapdoor_code import rll_decode, rll_encode, trapdoor_code # noqa: E402 + CHANNELS = [trapdoor, ising_channel] LOG2_PHI = math.log2((1 + math.sqrt(5)) / 2) NO_FEEDBACK_UPPER = 0.584 # Huleihel et al. (2021) @@ -74,3 +76,43 @@ def test_horstein_beats_the_no_feedback_capacity_of_the_trapdoor(): assert NO_FEEDBACK_UPPER < rate < LOG2_PHI lower, upper = code.capacity_gap(trapdoor(), method="montecarlo", samples=100, prng=np.random.default_rng(2)) assert 0 < lower <= upper + + +# ── the zero-error trapdoor code ──────────────────────────────────────── + + +@settings(max_examples=30, deadline=None) +@given(st.integers(0, 12).flatmap(lambda n: st.tuples(st.just(n), st.integers(0, 10**6)))) +def test_rll_enumeration_is_a_bijection(pair): + length, seed = pair + words = {rll_encode(i, length) for i in range(_rll_count(length))} + assert len(words) == _rll_count(length) + assert all("11" not in "".join(map(str, w)) for w in words) + index = seed % _rll_count(length) + assert rll_decode(rll_encode(index, length)) == index + + +def _rll_count(length): + a, b = 1, 2 + for _ in range(length): + a, b = b, a + b + return a + + +@pytest.mark.parametrize("k", [1, 2, 3, 4]) +def test_trapdoor_code_is_zero_error(k): + """Exact over every message and output: no errors, and E[N] = N + 3.5 flushing uses.""" + code = trapdoor_code(k) + pe, uses = code.performance(trapdoor(), method="exact", max_uses=70) + assert pe == pytest.approx(0.0, abs=1e-6) + assert uses == pytest.approx(code.length + 3.5, abs=1e-6) + + +def test_trapdoor_code_approaches_the_feedback_capacity(): + code = trapdoor_code(32) + pe, uses = code.performance( + trapdoor(), method="montecarlo", samples=200, prng=np.random.default_rng(3), max_uses=300 + ) + assert pe == 0.0 + assert NO_FEEDBACK_UPPER < 32 / uses < LOG2_PHI + assert 32 / (code.length - 1) == pytest.approx(LOG2_PHI, abs=0.03) From a026ac18cb4254677dd005fc4f87c55718f7292c Mon Sep 17 00:00:00 2001 From: Ryan James Date: Mon, 26 Oct 2026 13:14:24 -0600 Subject: [PATCH 3/4] Silence ty on dit feedback-code imports not yet in a released dit Co-authored-by: Cursor --- sofic/generators/trapdoor_code.py | 2 +- tests/test_channel_feedback.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/sofic/generators/trapdoor_code.py b/sofic/generators/trapdoor_code.py index 8bfa0b3..98f88e6 100644 --- a/sofic/generators/trapdoor_code.py +++ b/sofic/generators/trapdoor_code.py @@ -12,7 +12,7 @@ from collections.abc import Sequence from typing import Any -from dit.coding.feedback import FeedbackCode +from dit.coding.feedback import FeedbackCode # ty: ignore[unresolved-import] from dit.exceptions import ditException __all__ = ["TrapdoorCode", "trapdoor_code"] diff --git a/tests/test_channel_feedback.py b/tests/test_channel_feedback.py index b269be1..fcc7ee0 100644 --- a/tests/test_channel_feedback.py +++ b/tests/test_channel_feedback.py @@ -14,7 +14,7 @@ if not hasattr(dit_coding, "horstein"): pytest.skip("dit.coding has no feedback codes", allow_module_level=True) -from dit.coding import binary_symmetric_channel, horstein # noqa: E402 +from dit.coding import binary_symmetric_channel, horstein # noqa: E402 # ty: ignore[unresolved-import] from sofic.generators.trapdoor_code import rll_decode, rll_encode, trapdoor_code # noqa: E402 From 9a7de821bc8312b9869a84bafd033b5300d034eb Mon Sep 17 00:00:00 2001 From: Ryan James Date: Mon, 26 Oct 2026 13:26:01 -0600 Subject: [PATCH 4/4] Docs: drop trapdoor_code autodoc until dit ships feedback codes Co-authored-by: Cursor --- docs/generators/channel_coding.rst | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/docs/generators/channel_coding.rst b/docs/generators/channel_coding.rst index 8aa945c..6899614 100644 --- a/docs/generators/channel_coding.rst +++ b/docs/generators/channel_coding.rst @@ -82,7 +82,7 @@ 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. -:func:`~sofic.generators.trapdoor_code.trapdoor_code` is the scheme of +``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 @@ -119,7 +119,4 @@ API === .. autoclass:: FiniteStateChannel - :members: - -.. automodule:: sofic.generators.trapdoor_code - :members: TrapdoorCode, trapdoor_code, rll_encode, rll_decode + :members: \ No newline at end of file