Generator search results: reproduction error and generator structure

CHANGELOG.md

@@ -7,6 +7,12 @@ brute-force reference or a closed-form value and has a regression test.
 
 ### Generator verification
 
+- `GenerativeModelResult` records the returned generator's word-process deviation from the source
+  (`reproduction_error`, from the exact difference-automaton check) and its `structure`: state count,
+  weights, whether the state is a function of `(S+, S-)` (`functional`), and the largest violation of the
+  lifting conditions (`lifting_residual`; `lifted` when it vanishes). `certified` now also requires
+  `reproduction_error <= 1e-6`, so a closed bracket is never reported on an approximate generator.
+
 - The one-step generator equations (edge law, forward stationarity, summed future-side consistency) are
   necessary but not sufficient: a kernel with kernel-LP violation 0 can emit the wrong word process,
   because nothing forces a state's belief `R_w` to be the generated future's actual conditional. On survey

docs/generators/generative_models.rst

@@ -228,9 +228,14 @@ runs the same search when the bidirectional machine is already known.
    result.cg_exceeds_g                            # True when no G optimum generates
    result.separating_variables_generate           # some G optimum passed kernel LP
    result.model                                   # the generator, verified against the source
-
-The searches are local; ``certified`` reports when the bracket closes (for
-example the golden mean, where the support bound meets :math:`C_\mu`).
+   result.reproduction_error                      # its word-process deviation from the source
+   result.structure                               # states, weights, functional, lifted
+
+The searches are local; ``certified`` reports when the bracket closes and the generator reproduces the
+source to the same tolerance (for example the golden mean, where the support bound meets :math:`C_\mu`).
+``structure["lifted"]`` says whether the generator's channel satisfies the lifting conditions, so that
+drawing its states independently at both ends of each bidirectional edge also generates the process;
+coupled generators such as the interleaved product flag do not.
 
 Wyner and functional generators
 -------------------------------

sofic/generators/generative_search.py

@@ -71,6 +71,7 @@
 _MAX_RECTANGLES = 7
 _FLOW_THRESHOLDS = ((0.6,), (0.55, 0.65))
 _CERT_TOL = 1e-6
+_LIFTED_ATOL = 1e-8
 _LINE_SPLIT_FLOORS = (1e-2, 1e-3, 1e-5)
 # Bracket width (bits) above which the costly targeted line-split search runs.
 _TARGETED_GAP = 1e-2
@@ -106,7 +107,11 @@ class GenerativeModelResult:
     generator state, so ``cg_exceeds_g=True`` (no optimum found passes) is evidence that :math:`C_g` exceeds
     :math:`\\op{G}{S^+ : S^-}`, though only over the optima found. A passing kernel LP is only a candidate:
     ``cg_exceeds_g=False`` does not prove :math:`C_g = \\op{G}{S^+ : S^-}`. ``model`` and ``upper`` always come
-    from a candidate whose word process was checked against the source (``_reproduction_error``).
+    from a candidate whose word process was checked against the source (``_reproduction_error``);
+    ``reproduction_error`` is that deviation, and ``certified`` requires it to be within the certification
+    tolerance. ``structure`` describes the generator: its state count and weights, whether its state is a
+    function of :math:`(S^+, S^-)` (``functional``), and the largest violation of the lifting conditions
+    (``lifting_residual``, bits; ``lifted`` when it vanishes, so the lift of its channel generates too).
     """
 
     model: CoupledGenerativeModel
@@ -124,11 +129,14 @@ class GenerativeModelResult:
     lower_certified: bool = True
     objective: str = "entropy"
     budget_exhausted: bool = False
+    reproduction_error: float | None = None
+    structure: dict[str, Any] | None = None
 
     @property
     def certified(self) -> bool:
-        """Whether the upper and lower bounds meet (to 1e-6 bits)."""
-        return self.upper - self.lower <= _CERT_TOL
+        """Whether the bounds meet (to 1e-6 bits) and the model's word process is the source's (to 1e-6)."""
+        reproduced = self.reproduction_error is None or self.reproduction_error <= _CERT_TOL
+        return self.upper - self.lower <= _CERT_TOL and reproduced
 
     @property
     def gap(self) -> float:
@@ -157,6 +165,26 @@ def _reliability_summary(result: GenerativeModelResult) -> dict[str, Any]:
         "separating_variables_generate": result.separating_variables_generate,
         "persistent_crypticity_theorem_applies": result.persistent_crypticity_theorem_applies,
         "budget_exhausted": result.budget_exhausted,
+        "reproduction_error": result.reproduction_error,
+        "structure": result.structure,
+    }
+
+
+def _generator_structure(problem: _GeneratorProblem, seed: tuple[np.ndarray, ...], cutoff: float) -> dict[str, Any]:
+    """State weights, functionality and lifting residual of a generator ``(pi, Q, R, K)``."""
+    pi, Q, R = (np.asarray(part, dtype=float) for part in seed[:3])
+    active = pi > cutoff
+    lifted = problem.lifted
+    rows_plus, rows_minus = lifted.plus.argmax(axis=1), lifted.minus.argmax(axis=1)
+    channel = (pi[active] * Q[active][:, rows_plus].T * R[active][:, rows_minus].T) / lifted.p[:, None]
+    channel /= channel.sum(axis=1, keepdims=True)
+    residual = float(lifted.markov_residual(channel))
+    return {
+        "n_states": int(active.sum()),
+        "weights": sorted((float(w) for w in pi[active] / pi[active].sum()), reverse=True),
+        "functional": bool(((channel > cutoff).sum(axis=1) == 1).all()),
+        "lifting_residual": residual,
+        "lifted": residual <= _LIFTED_ATOL,
     }
 
 
@@ -301,7 +329,8 @@ def _try_best_candidate(
         pi, Q, R, K = seed
         result = CoupledSearchResult(pi, Q, R, K, _entropy(pi), problem.residual(pi, Q, R, K), name)
         model = _model_from_kernel(problem, result, 1e-9)
-        if _reproduction_error(model, reference, max_length=length) > reproduction_atol:
+        error = float(_reproduction_error(model, reference, max_length=length))
+        if error > reproduction_atol:
             continue
         upper = float(value) if objective == "wyner" else float(model.state_entropy())
         return GenerativeModelResult(
@@ -319,6 +348,8 @@ def _try_best_candidate(
             cg_exceeds_g=cg_exceeds_g,
             lower_certified=lower_certified,
             objective=objective,
+            reproduction_error=error,
+            structure=_generator_structure(problem, seed, 1e-9),
         )
     return None
 

tests/test_coupled_generators.py

@@ -127,6 +127,21 @@ def test_memory_lift_of_the_flag_is_the_minimal_interleaved_generator() -> None:
     assert float(-(pi * np.log2(pi)).sum()) == pytest.approx(2 * _h(4 / 9), abs=1e-12)
 
 
+def test_generator_structure_tells_lifted_from_coupled_generators() -> None:
+    from sofic.generators.generative_search import _generator_structure
+
+    problem = cg._GeneratorProblem(_interleaved(0.1))
+    machine = _generator_structure(problem, problem.epsilon_seed(), 1e-9)
+    assert machine["n_states"] == 4
+    assert machine["weights"] == pytest.approx([0.25] * 4)
+    assert machine["functional"] and machine["lifted"]
+
+    coupled = _generator_structure(problem, cg.memory_lift(problem, _flag_channel(problem, 0.1)), 1e-9)
+    assert not coupled["functional"]
+    assert not coupled["lifted"]
+    assert coupled["lifting_residual"] == pytest.approx(0.019, abs=1e-3)
+
+
 def test_generator_kernel_decides_generator_states() -> None:
     rll = cg._GeneratorProblem(_rll())
     assert cg.generator_kernel(rll, np.eye(len(rll.lifted.pairs)))[0] == pytest.approx(0.0, abs=1e-10)

tests/test_generative_model.py

@@ -27,6 +27,31 @@ def test_generative_model_from_forward_epsilon_machine_certifies_golden_mean() -
     target = float(-(1 / 3) * np.log2(1 / 3) - (2 / 3) * np.log2(2 / 3))
     assert abs(result.upper - target) < 1e-5
     assert "universal coupled" not in result.upper_source
+    assert result.reproduction_error is not None and result.reproduction_error < 1e-9
+    assert result.structure is not None
+    assert result.structure["n_states"] == 2
+    assert result.structure["functional"] and result.structure["lifted"]
+    assert result.stages[-1]["structure"] == result.structure
+
+
+def test_certified_requires_the_model_to_reproduce_the_source() -> None:
+    from sofic.generators.generative_search import GenerativeModelResult
+
+    def result(error):
+        return GenerativeModelResult(
+            model=None,
+            upper=1.0,
+            lower=1.0,
+            upper_source="test",
+            lower_source="test",
+            exact_common_information_estimate=None,
+            separating_variables_generate=False,
+            reproduction_error=error,
+        )
+
+    assert result(None).certified
+    assert result(1e-9).certified
+    assert not result(1e-4).certified
 
 
 def test_generative_model_accepts_bidirectional_input() -> None:

tests/test_generative_model_invariants.py

@@ -0,0 +1,83 @@
+"""Invariants of the generator search on random ε-machines, and known values of C_g."""
+
+from __future__ import annotations
+
+import itertools
+
+import numpy as np
+import pytest
+from hypothesis import assume, given, settings
+from hypothesis import strategies as st
+
+from sofic.examples.epsilon_machines import from_symbol_matrices
+from sofic.generators.bidirectional_epsilon_machine import BidirectionalEpsilonMachine
+from sofic.generators.epsilon_machine import EpsilonMachine
+from sofic.generators.generative_search import generative_model_from_bidirectional
+
+_TOL = 1e-7
+
+
+def _topologies(n: int) -> list[tuple[tuple[int | None, ...], ...]]:
+    """Binary topologies on ``n`` states: a target or ``None`` per (state, symbol), strongly connected."""
+    out = []
+    for targets in itertools.product([None, *range(n)], repeat=2 * n):
+        rows = tuple(tuple(targets[2 * s + x] for x in range(2)) for s in range(n))
+        if any(all(t is None for t in row) for row in rows):
+            continue
+        reach = [{s} for s in range(n)]
+        for _ in range(n):
+            reach = [r | {t for u in r for t in rows[u] if t is not None} for r in reach]
+        if all(r == set(range(n)) for r in reach):
+            out.append(rows)
+    return out
+
+
+@st.composite
+def _epsilon_machines(draw):
+    rows = draw(st.sampled_from(_topologies(2) + _topologies(3)))
+    n = len(rows)
+    mats = {x: np.zeros((n, n)) for x in (0, 1)}
+    for s, row in enumerate(rows):
+        live = [x for x in (0, 1) if row[x] is not None]
+        p = draw(st.floats(min_value=0.1, max_value=0.9)) if len(live) == 2 else 1.0
+        for x, prob in zip(live, (p, 1 - p) if len(live) == 2 else (1.0,), strict=True):
+            mats[x][s, row[x]] += prob
+    machine = from_symbol_matrices(tuple(range(n)), (0, 1), mats)
+    return machine if isinstance(machine, EpsilonMachine) else EpsilonMachine.from_hmm(machine)
+
+
+@pytest.mark.slow
+@given(forward=_epsilon_machines())
+@settings(max_examples=12, deadline=None, derandomize=True)
+def test_generator_search_bracket_is_consistent(forward):
+    try:
+        bidir = BidirectionalEpsilonMachine.from_forward(forward)
+    except Exception:  # noqa: BLE001
+        assume(False)
+    result = generative_model_from_bidirectional(bidir, time_budget=20, rng=np.random.default_rng(0))
+    excess = float(bidir.excess_entropy())
+    c_mu = min(
+        float(bidir.forward_epsilon_machine().statistical_complexity()),
+        float(bidir.reverse_epsilon_machine().statistical_complexity()),
+    )
+    assert excess - _TOL <= result.lower <= result.upper + _TOL
+    assert result.upper <= c_mu + _TOL
+    assert result.reproduction_error is not None and result.reproduction_error <= 1e-5
+    assert float(result.model.state_entropy()) == pytest.approx(result.upper, abs=_TOL)
+    assert result.structure["n_states"] == len(list(result.model.states()))
+    assert sum(result.structure["weights"]) == pytest.approx(1.0)
+
+
+@pytest.mark.slow
+def test_run_length_limited_generative_complexity_is_three_halves():
+    """(2, inf) RLL: C_g = 3/2 = C_mu > G (computer-assisted proof in the cpfci manuscript)."""
+    machine = from_symbol_matrices(
+        ("A", "B", "C"),
+        (0, 1),
+        {0: np.array([[0.5, 0, 0], [0, 0, 1.0], [1.0, 0, 0]]), 1: np.array([[0, 0.5, 0], [0, 0, 0], [0, 0, 0]])},
+    )
+    bidir = BidirectionalEpsilonMachine.from_forward(machine)
+    result = generative_model_from_bidirectional(bidir, rng=np.random.default_rng(0))
+    assert result.upper == pytest.approx(1.5, abs=1e-6)
+    assert result.lower == pytest.approx(1.5, abs=1e-6)
+    assert result.certified