ε-machines of asymptotically synchronizing processes

CHANGELOG.md

@@ -5,6 +5,19 @@
 A second correctness review: every fix below was reproduced against a
 brute-force reference or a closed-form value and has a regression test.
 
+### ε-machine construction
+
+- Processes whose hidden state no finite word reveals (e.g. a period-2 process whose phase stays
+  uncertain) synchronize only asymptotically. Their mixed-state presentation either closed on hundreds of
+  beliefs within round-off of the simplex boundary (rounding acts as a reflecting barrier for the
+  posterior's random walk) or exceeded `max_states`, so `EpsilonMachine.from_hmm` returned a huge
+  machine or raised `MixedStateExplosionError` although the ε-machine is finite. It now seeds the
+  presentation from a limit belief (the heaviest near-boundary belief, or the belief after filtering a
+  long sample) with sub-`1e-8` entries set to zero, and accepts the result only if it generates the
+  source process exactly (difference-automaton test). On the 138 minimal binary 3-state topologies,
+  9 reverse-machine "explosions" are now 3-state reverse ε-machines, and finite/infinite now agrees
+  with `reverse_epsilon_machine_is_finite` on every topology.
+
 ### Generator verification
 
 - The one-step generator equations (edge law, forward stationarity, summed future-side consistency) are

sofic/generators/epsilon_construction.py

@@ -14,11 +14,11 @@
 
 import numpy as np
 
-from sofic.exceptions import StochasticValidationError
+from sofic.exceptions import MixedStateExplosionError, StochasticValidationError
 from sofic.generators.base import HiddenMarkovModel
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.mealy import MealyHMM
-from sofic.generators.mixed_state import MixedStatePresentation
+from sofic.generators.mixed_state import MixedState, MixedStatePresentation
 from sofic.generators.prob import (
     array_sum,
     as_prob,
@@ -34,19 +34,102 @@
 
 TransitionSignature = tuple[tuple[Any, int, Any], ...]
 
+_SUPPORT_ATOL = 1e-8
+_LIMIT_VERIFY_ATOL = 1e-9
+
 
 def build_epsilon_machine(hmm: HiddenMarkovModel, *, max_states: int = 10_000) -> EpsilonMachine:
     hmm = hmm.to_mealy()
-    presentation = _unifilar_presentation(hmm, max_states=max_states)
+    try:
+        presentation = _unifilar_presentation(hmm, max_states=max_states)
+    except MixedStateExplosionError:
+        limit = _limit_machine(hmm, _filtered_belief(hmm), max_states=max(64, max_states // 10))
+        if limit is None:
+            raise
+        return limit
+    machine = _minimize(presentation, constraints=getattr(hmm, "symbol_constraints", None))
+    seed = _boundary_belief(presentation)
+    if seed is None:
+        return machine
+    limit = _limit_machine(hmm, seed, max_states=len(list(presentation.states())) - 1)
+    if limit is None or len(list(limit.states())) >= len(list(machine.states())):
+        return machine
+    return limit
+
+
+def _minimize(presentation: MealyHMM, *, constraints: Any = None) -> EpsilonMachine:
     stationary = presentation.stationary_distribution()
     if is_zero(array_sum(stationary)):
         raise StochasticValidationError("generator must have a stationary distribution")
 
-    constraints = getattr(presentation, "symbol_constraints", None) or getattr(hmm, "symbol_constraints", None)
+    constraints = getattr(presentation, "symbol_constraints", None) or constraints
     partitions = _refine_probabilistic_partitions(presentation, constraints=constraints)
     return _quotient_machine(presentation, partitions, stationary, constraints=constraints)
 
 
+def _boundary_belief(presentation: MealyHMM) -> np.ndarray | None:
+    """The heaviest recurrent belief, if the recurrent beliefs sit within round-off of the simplex boundary.
+
+    When no finite word reveals the hidden state (e.g. the phase of a period-2 process), beliefs
+    converge to the causal states without reaching them. Rounding then closes the mixed-state
+    presentation on a reflecting boundary: a large recurrent class of beliefs with positive entries
+    below ``_SUPPORT_ATOL``, carrying geometrically decaying mass. Returns ``None`` for any other
+    presentation.
+    """
+    states = [state for state in presentation.states() if isinstance(state, MixedState)]
+    if not states or any(state.is_symbolic() for state in states):
+        return None
+    beliefs = np.array([state.as_array() for state in states])
+    if not ((beliefs > 0.0) & (beliefs < _SUPPORT_ATOL)).any():
+        return None
+    stationary = np.asarray(presentation.stationary_distribution(), dtype=float)
+    mass = dict(zip(presentation.reindex().states, stationary, strict=True))
+    return max(states, key=lambda state: mass[state]).as_array()
+
+
+def _filtered_belief(hmm: MealyHMM, *, length: int = 4096) -> np.ndarray | None:
+    """Belief over ``hmm``'s states after Bayes-filtering a sampled word of ``length`` symbols."""
+    from sofic.generators.matrices import symbol_matrices
+
+    if has_symbolic(hmm.stationary_distribution().ravel()):
+        return None
+    matrices = {symbol: np.asarray(matrix, dtype=float) for symbol, matrix in symbol_matrices(hmm).items()}
+    symbols, _ = hmm.sample(length, rng=np.random.default_rng(0))
+    belief = np.asarray(hmm.stationary_distribution(), dtype=float)
+    for symbol in symbols:
+        belief = belief @ matrices[symbol]
+        belief /= belief.sum()
+    return belief
+
+
+def _limit_machine(hmm: MealyHMM, seed: np.ndarray | None, *, max_states: int) -> EpsilonMachine | None:
+    """The ε-machine generated from a limit belief, if it is finite and generates ``hmm``'s process.
+
+    ``seed`` approximates a causal state of a process that synchronizes only asymptotically; entries
+    below ``_SUPPORT_ATOL`` are set to zero, and the mixed-state presentation it generates under exact
+    Bayes updates (at most ``max_states`` beliefs) is minimized. The result is accepted only if it
+    generates the same process as ``hmm``
+    (:func:`~sofic.generators.minimal_generative_model._equivalence_error`): a finite unifilar
+    presentation of the process minimizes to its ε-machine.
+    """
+    if seed is None:
+        return None
+    seed = np.where(seed < _SUPPORT_ATOL, 0.0, seed)
+    from sofic.generators.minimal_generative_model import _equivalence_error
+    from sofic.generators.mixed_state_construction import build_mixed_state_presentation
+
+    try:
+        limit = build_mixed_state_presentation(
+            hmm, initial_mixed_state=MixedState.from_vector(seed), max_states=max(1, max_states)
+        )
+    except MixedStateExplosionError:
+        return None
+    candidate = _minimize(_unifilar_presentation(limit))
+    if _equivalence_error(candidate, hmm) > _LIMIT_VERIFY_ATOL:
+        return None
+    return candidate
+
+
 def _unifilar_presentation(hmm: MealyHMM, *, max_states: int = 10_000) -> MealyHMM:
     """Return a row-unifilar presentation generating the same process."""
     if not hmm.is_unifilar():

tests/test_asymptotic_synchronization.py

@@ -0,0 +1,34 @@
+"""ε-machines of processes whose hidden state no finite word reveals."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from sofic.examples.epsilon_machines import from_symbol_matrices
+from sofic.generators.epsilon_machine import EpsilonMachine
+from sofic.generators.minimal_generative_model import _equivalence_error
+from sofic.generators.reversal import time_reverse_stochastic
+
+
+def _paired_bits(e: float):
+    """Period 2: a fair bit, then a copy of it flipped with probability ``e``; the phase is never certain."""
+    mats = {
+        0: np.array([[0, 0.5, 0], [1 - e, 0, 0], [e, 0, 0]]),
+        1: np.array([[0, 0, 0.5], [e, 0, 0], [1 - e, 0, 0]]),
+    }
+    return from_symbol_matrices(("A", "B0", "B1"), (0, 1), mats)
+
+
+# 0.4: rounding closes the mixed-state presentation on hundreds of near-limit beliefs;
+# 0.1: the mixed-state presentation exceeds ``max_states`` first.
+@pytest.mark.parametrize("e", [0.4, 0.1])
+def test_reverse_epsilon_machine_of_asymptotically_synchronizing_process(e):
+    forward = _paired_bits(e)
+    reversed_hmm = time_reverse_stochastic(forward)
+
+    reverse = EpsilonMachine.from_hmm(reversed_hmm)
+
+    assert len(list(reverse.states())) == 3
+    assert float(reverse.statistical_complexity()) == pytest.approx(1.5)
+    assert _equivalence_error(reverse, reversed_hmm) < 1e-12