diff --git a/CHANGELOG.md b/CHANGELOG.md index ce47d10..3db7e9a 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/generators/generative_models.rst b/docs/generators/generative_models.rst index 44df4d3..d42637e 100644 --- a/docs/generators/generative_models.rst +++ b/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 ------------------------------- diff --git a/sofic/generators/generative_search.py b/sofic/generators/generative_search.py index 4534e84..def710d 100644 --- a/sofic/generators/generative_search.py +++ b/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 diff --git a/tests/test_coupled_generators.py b/tests/test_coupled_generators.py index d8a2aaf..b56f74d 100644 --- a/tests/test_coupled_generators.py +++ b/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) diff --git a/tests/test_generative_model.py b/tests/test_generative_model.py index f8b1029..62c797a 100644 --- a/tests/test_generative_model.py +++ b/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: diff --git a/tests/test_generative_model_invariants.py b/tests/test_generative_model_invariants.py new file mode 100644 index 0000000..73e2d47 --- /dev/null +++ b/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