From 9e03ace91e4750c40c548a9eaa14291e06491138 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Wed, 7 Oct 2026 01:48:04 -0600 Subject: [PATCH 1/3] Property and metamorphic test suite; fixes it found MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds tests/test_properties_*.py: shifts, generators, inference, automata (DFA/NFA, Büchi, Wheeler, VPA/NWA, transducers, learners), serialization, viz and examples, checked against brute-force oracles and invariants. Fixes found by the suite: Fischer cover on reducible shifts, Parry measure transient states, TikZ newlines, from_hmm beyond 26 causal states, word_probability below 1e-15, backward(normalize=True) final row, complete() reusing an accepting trap after complement. reverse_is_finite is documented as mixed-state finiteness; experimental reverse_epsilon_machine_is_finite decides the recurrent question, and reversal warns on numeric truncation. YAML uses libyaml's C loader/dumper when available. Co-authored-by: Cursor --- CHANGELOG.md | 27 +- docs/generators/epsilon_machine.rst | 23 +- sofic/automata/algorithms.py | 31 +- .../generators/bidirectional_construction.py | 3 +- sofic/generators/epsilon_construction.py | 3 +- sofic/generators/epsilon_machine.py | 41 +- sofic/generators/reversal.py | 91 ++- sofic/generators/synchronization.py | 12 +- sofic/generators/words.py | 4 +- sofic/inference/hmm/filtering.py | 2 +- sofic/serialization.py | 8 +- sofic/shifts/cover_construction.py | 36 ++ sofic/shifts/topological_anatomy.py | 5 +- sofic/viz/_tikz_format.py | 5 +- tests/test_epsilon_machine.py | 59 ++ tests/test_properties_automata.py | 270 ++++++++ tests/test_properties_automata_learning.py | 164 +++++ .../test_properties_automata_omega_wheeler.py | 146 +++++ tests/test_properties_automata_transducers.py | 319 ++++++++++ tests/test_properties_automata_visibly.py | 328 ++++++++++ tests/test_properties_examples.py | 221 +++++++ tests/test_properties_generators.py | 452 +++++++++++++ tests/test_properties_inference.py | 315 ++++++++++ tests/test_properties_serialization.py | 252 ++++++++ tests/test_properties_shifts.py | 592 ++++++++++++++++++ tests/test_properties_viz.py | 185 ++++++ 26 files changed, 3548 insertions(+), 46 deletions(-) create mode 100644 tests/test_properties_automata.py create mode 100644 tests/test_properties_automata_learning.py create mode 100644 tests/test_properties_automata_omega_wheeler.py create mode 100644 tests/test_properties_automata_transducers.py create mode 100644 tests/test_properties_automata_visibly.py create mode 100644 tests/test_properties_examples.py create mode 100644 tests/test_properties_generators.py create mode 100644 tests/test_properties_inference.py create mode 100644 tests/test_properties_serialization.py create mode 100644 tests/test_properties_shifts.py create mode 100644 tests/test_properties_viz.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 730fbb2..9bcecc4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -54,13 +54,24 @@ brute-force reference or a closed-form value and has a regression test. `reverse()` returns the mirror-image shift; `from_adjacency` assigns symbols in sorted order and honors subclasses; the synchronization methods on `SoficShift` (`markov_order`, `reset_threshold`, ...) no longer raise - `AttributeError`. + `AttributeError`; the Fischer cover raises on reducible shifts instead of + returning a single terminal component; the Parry measure drops transient + states. - **Generators:** `cryptic_order` is 0 for zero-crypticity processes (even, periodic) via a new algorithm checked against brute force; `log_word_probability` no longer underflows; channel statistical complexity uses the driven process's occupation law (it was hash-seed dependent); mixed-state enumeration of infinite belief sets fails in under a second; - numeric stationary distributions warn when not unique. + numeric stationary distributions warn when not unique; + `EpsilonMachine.from_hmm` handles more than 26 causal states; + `word_probability` no longer returns 0 below 1e-15 (and + `conditional_word_probability` no longer raises `ZeroDivisionError` there); + `backward(normalize=True)` normalizes the final row; `reverse_is_finite` is + documented as deciding finiteness of the reverse mixed states, a sufficient + but not necessary condition for a finite reverse ε-machine (alternating + biased coins has infinitely many reverse beliefs but a two-state reverse + machine); `from_time_reversed` and `to_bidirectional` warn when numeric belief + merging returns a finite truncation of an infinite reverse machine. - **Inference:** Bayesian posterior machines start where the data start (they started in its final state, giving −∞ likelihoods and infinite BIC); WAIC scores from the stationary distribution; `learn_stack_hmm_mle` and the stack @@ -75,6 +86,8 @@ brute-force reference or a closed-form value and has a regression test. `wheeler_index`, word counting, empty marks, `minimize_wheeler` on dead states, and `colex_width`; modular VPA `minimize` on valid well-matched VPAs; intersection, difference and complement with empty-language operands; + `complete()` no longer reuses an accepting trap left by a previous + complement (double complement could return the wrong language); `right_quotient` by the empty word; `sofic.automata.languages.residuals` callable again; IDFA enumeration includes missing transitions before the first flag (`count_accessible_idfa(2, 2)` is 45), which changes some canonical @@ -88,13 +101,21 @@ brute-force reference or a closed-form value and has a regression test. - **Serialization and viz:** YAML stores observation, output and stack alphabets and supports sympy and `Fraction` values; `copy()` deep-copies attributes; Graphviz/TikZ node names are unique; TikZ escapes `^`, `~`, `\` - for math mode. + for math mode and handles newlines in labels; YAML uses libyaml's C loader + and dumper when available. ### New - Hypothesis strategies in `sofic.testing` for NFAs, Büchi automata and lassos, Wheeler NFAs, Markov chains, Mealy HMMs, sofic shifts, SFTs, VPAs, NWAs, and Mealy transducers; `ci` and `nightly` Hypothesis profiles (`HYPOTHESIS_PROFILE`). +- Experimental `EpsilonMachine.reverse_epsilon_machine_is_finite()` (and + `sofic.generators.reversal.reverse_epsilon_machine_is_finite`): decides whether + the reverse ε-machine has finitely many *recurrent* causal states, which + `reverse_is_finite` only bounds from one side. +- A property and metamorphic test suite (`tests/test_properties_*.py`) checking + shifts, generators, inference, automata, serialization, viz and examples + against brute-force oracles and invariants. ## 0.3.0 diff --git a/docs/generators/epsilon_machine.rst b/docs/generators/epsilon_machine.rst index e6e3264..ff9c5c3 100644 --- a/docs/generators/epsilon_machine.rst +++ b/docs/generators/epsilon_machine.rst @@ -75,12 +75,21 @@ infinite — that is the entropy of the retrodictive stationary distribution, an it converges when those weights decay geometrically, as in the example above. The obstruction here is cardinality, not divergence. -:meth:`reverse_is_finite` decides this in polynomial time without enumerating -beliefs. Since the retrodictive causal states are the normalized vectors -:math:`(\Pr(x \mid s))_s`, the reverse machine is finite exactly when every -cycle of the pair graph on :math:`S \times S`, weighted by +:meth:`reverse_is_finite` gives a polynomial-time sufficient condition without +enumerating beliefs. The reverse mixed states are the normalized vectors +:math:`(\Pr(x \mid s))_s`, and they are finite exactly when every cycle of the +pair graph on :math:`S \times S`, weighted by :math:`p(a \mid s) / p(a \mid s')`, has weight one — the twins property -:cite:`Mohri2009` :cite:`AllauzenMohri2003`. +:cite:`Mohri2009` :cite:`AllauzenMohri2003`. Finitely many mixed states implies +a finite reverse machine, but not conversely: for alternating biased coins the +mixed states are infinite, yet they converge to two recurrent causal states. + +:meth:`reverse_epsilon_machine_is_finite` (experimental) decides the recurrent +question instead: the reverse machine is infinite exactly when a pair-graph +component with a non-unit cycle can reach the diagonal, so that a grown ratio +is frozen by the two forward paths merging. :meth:`from_time_reversed` warns +when it says infinite but numeric belief merging closed the construction +anyway. This is not a structural property: two machines with identical transition structure can differ, because a cycle weight may equal one only by algebraic @@ -93,10 +102,12 @@ API === .. autoclass:: EpsilonMachine - :members: from_hmm, from_sequence, from_time_reversed, to_bidirectional, block_entropy_diagram, block_entropy_estimates, plot_block_entropy_diagram, block_convergence_diagram, block_convergence_estimates, plot_block_convergence_diagram, caekl_block_information, caekl_rate, caekl_intercept, caekl_rate_converged, approximate_entropy_rate, approximate_excess_entropy, approximate_information_anatomy, statistical_complexity, bidirectional_statistical_complexity, excess_entropy, predicted_information, bound_information, ephemeral_information, information_anatomy, caekl_causal_information, crypticity, bidirectional_crypticity, causal_irreversibility, stored_information_decomposition, transient_information, oracular_information, gauge_information, predictability_gain, structural_information, thermodynamic_depth, spectral_complexity, markov_order, is_markov, cryptic_order, reset_threshold, synchronizing_word, is_exactly_synchronizable, is_asymptotically_synchronizable, is_definite, reverse_is_finite + :members: from_hmm, from_sequence, from_time_reversed, to_bidirectional, block_entropy_diagram, block_entropy_estimates, plot_block_entropy_diagram, block_convergence_diagram, block_convergence_estimates, plot_block_convergence_diagram, caekl_block_information, caekl_rate, caekl_intercept, caekl_rate_converged, approximate_entropy_rate, approximate_excess_entropy, approximate_information_anatomy, statistical_complexity, bidirectional_statistical_complexity, excess_entropy, predicted_information, bound_information, ephemeral_information, information_anatomy, caekl_causal_information, crypticity, bidirectional_crypticity, causal_irreversibility, stored_information_decomposition, transient_information, oracular_information, gauge_information, predictability_gain, structural_information, thermodynamic_depth, spectral_complexity, markov_order, is_markov, cryptic_order, reset_threshold, synchronizing_word, is_exactly_synchronizable, is_asymptotically_synchronizable, is_definite, reverse_is_finite, reverse_epsilon_machine_is_finite .. autofunction:: sofic.generators.reversal.reverse_is_finite +.. autofunction:: sofic.generators.reversal.reverse_epsilon_machine_is_finite + .. autoclass:: sofic.generators.block_entropy.BlockEntropyDiagram :members: plot, transient_information .. autoclass:: sofic.generators.block_entropy.BlockEntropyEstimates diff --git a/sofic/automata/algorithms.py b/sofic/automata/algorithms.py index ba49bb4..40d066a 100644 --- a/sofic/automata/algorithms.py +++ b/sofic/automata/algorithms.py @@ -37,10 +37,22 @@ def trim(aut: L) -> L: # noqa: UP047 - keep Python 3.11 compatibility. def complete(dfa: DFA, alphabet: frozenset[Any] | None = None) -> DFA: """Add a trap state so every state has one outgoing transition per symbol.""" + return _complete_with_trap(dfa, alphabet)[0] + + +def _complete_with_trap(dfa: DFA, alphabet: frozenset[Any] | None = None) -> tuple[DFA, Hashable | None]: + """Complete ``dfa``; also return the trap state if it was newly added. + + An existing ``_TRAP`` state is reused only while it is still a rejecting + sink: a complemented DFA keeps its trap but makes it accepting. + """ symbols = alphabet if alphabet is not None else _effective_alphabet(dfa) result = dfa.copy() - trap = _TRAP - if trap not in result.graph.nx: + trap: Hashable = _TRAP + if trap in result.graph.nx and not _is_rejecting_sink(result, trap): + trap = Sentinel("trap") + added = trap not in result.graph.nx + if added: result.graph.add_state(trap) for state in list(result.states()): @@ -51,11 +63,18 @@ def complete(dfa: DFA, alphabet: frozenset[Any] | None = None) -> DFA: result.add_transition(state, trap, symbol) for symbol in symbols: - result.add_transition(trap, trap, symbol) + if not result.delta(trap, symbol): + result.add_transition(trap, trap, symbol) if symbols: result.input_alphabet = frozenset(symbols) | result.input_alphabet - return result + return result, trap if added else None + + +def _is_rejecting_sink(dfa: DFA, state: Hashable) -> bool: + return state not in dfa.accepting_states and all( + transition.target == state for transition in dfa.graph.out_transitions(state) + ) def reverse(aut: NFA | DFA) -> NFA: @@ -160,13 +179,13 @@ def nerode_partition( -- need that distinction. """ symbols = alphabet if alphabet is not None else _effective_alphabet(dfa) - work = complete(dfa, symbols) + work, trap = _complete_with_trap(dfa, symbols) states = sorted(work.states(), key=repr) if not states: return () partition = _moore_refine(work, states, symbols) - blocks = (frozenset(block) - {_TRAP} for block in partition) + blocks = (frozenset(block) - {trap} for block in partition) return tuple(block for block in blocks if block) diff --git a/sofic/generators/bidirectional_construction.py b/sofic/generators/bidirectional_construction.py index cd25da4..682e37f 100644 --- a/sofic/generators/bidirectional_construction.py +++ b/sofic/generators/bidirectional_construction.py @@ -23,7 +23,7 @@ sum_probs, zeros, ) -from sofic.generators.reversal import time_reverse_stochastic +from sofic.generators.reversal import time_reverse_stochastic, warn_if_reverse_truncated from sofic.generators.stationary import ( stationary_distribution_from_transition, stationary_distribution_hmm, @@ -539,6 +539,7 @@ def infer_reverse_epsilon_machine(forward: EpsilonMachine) -> EpsilonMachine: """ rev_hmm = time_reverse_stochastic(forward) reverse = EpsilonMachine.from_hmm(rev_hmm) + warn_if_reverse_truncated(forward) return _relabel_collision_free(reverse, forward) diff --git a/sofic/generators/epsilon_construction.py b/sofic/generators/epsilon_construction.py index 490e444..4a80e77 100644 --- a/sofic/generators/epsilon_construction.py +++ b/sofic/generators/epsilon_construction.py @@ -135,7 +135,8 @@ def _quotient_machine( state_map[state] = index graph = TransitionGraph() - labels = sequential_labels(len(partitions)) + count = len(partitions) + labels = sequential_labels(count) if count <= 26 else tuple(range(count)) for label in labels: graph.add_state(label) label_for_index = {index: labels[index] for index in range(len(partitions))} diff --git a/sofic/generators/epsilon_machine.py b/sofic/generators/epsilon_machine.py index 55536b5..5a5055b 100644 --- a/sofic/generators/epsilon_machine.py +++ b/sofic/generators/epsilon_machine.py @@ -526,20 +526,34 @@ def is_asymptotically_synchronizable(self) -> bool: return is_asymptotically_synchronizable_from_graph(graph_from_epsilon_machine(self)) def reverse_is_finite(self) -> bool: - """Return whether the reverse ε-machine has finitely many causal states. + """Return whether the reverse mixed-state presentation is finite. Decided in polynomial time by the twins property rather than by enumerating beliefs; see - :func:`~sofic.generators.reversal.reverse_is_finite`. Use this before - :meth:`from_time_reversed` or :meth:`causal_irreversibility`, both of - which raise - :class:`~sofic.exceptions.MixedStateExplosionError` when it returns - ``False``. + :func:`~sofic.generators.reversal.reverse_is_finite`. ``True`` + guarantees that :meth:`from_time_reversed` and + :meth:`causal_irreversibility` succeed. ``False`` means infinitely many + reverse beliefs, which may still converge to a finite reverse + ε-machine; those methods raise + :class:`~sofic.exceptions.MixedStateExplosionError` only when the + recurrent part does not close. """ from sofic.generators.reversal import reverse_is_finite return reverse_is_finite(self) + def reverse_epsilon_machine_is_finite(self) -> bool: + """Return whether the reverse ε-machine has finitely many recurrent causal states. + + Experimental; see + :func:`~sofic.generators.reversal.reverse_epsilon_machine_is_finite`. + Unlike :meth:`reverse_is_finite` it ignores transient beliefs, so it is + ``True`` for alternating biased coins. + """ + from sofic.generators.reversal import reverse_epsilon_machine_is_finite + + return reverse_epsilon_machine_is_finite(self) + def is_definite(self) -> bool: """Return whether the ε-machine is a definite automaton (finite Markov order). @@ -562,17 +576,26 @@ def from_time_reversed(cls, forward: EpsilonMachine) -> EpsilonMachine: ε-machine does not imply a finite reverse one -- the "explosive irreversibility" of :cite:`Ellison2011` -- and when the reverse is infinite there is no presentation to return. + + Warns + ----- + RuntimeWarning + If the belief set closed only because numeric beliefs merge within + tolerance while + :func:`~sofic.generators.reversal.reverse_epsilon_machine_is_finite` + says the reverse machine is infinite; the result is then a finite + approximation. """ from sofic.exceptions import ( MixedStateExplosionError, StochasticValidationError, UnifilarityError, ) - from sofic.generators.reversal import time_reverse_stochastic + from sofic.generators.reversal import time_reverse_stochastic, warn_if_reverse_truncated rev_hmm = time_reverse_stochastic(forward) try: - return cls.from_hmm(rev_hmm) + reverse = cls.from_hmm(rev_hmm) except MixedStateExplosionError: # Row-normalizing would return a non-unifilar machine whose state # entropy coincides with the forward C_mu, reporting Delta C_mu = 0 @@ -580,6 +603,8 @@ def from_time_reversed(cls, forward: EpsilonMachine) -> EpsilonMachine: raise except (StochasticValidationError, UnifilarityError): return _row_normalized_presentation(rev_hmm) + warn_if_reverse_truncated(forward) + return reverse def _row_normalized_presentation(hmm: MealyHMM) -> EpsilonMachine: diff --git a/sofic/generators/reversal.py b/sofic/generators/reversal.py index d8fb301..e867cd2 100644 --- a/sofic/generators/reversal.py +++ b/sofic/generators/reversal.py @@ -34,22 +34,29 @@ def is_markov_like(model: StateMachine) -> bool: def reverse_is_finite(model: StateMachine, *, rtol: float = 1e-9) -> bool: - """Decide whether the reverse ε-machine of a unifilar ``model`` has finitely many states. + """Decide whether the reverse mixed-state presentation of a unifilar ``model`` is finite. A finite forward ε-machine does not imply a finite reverse one (:cite:`Ellison2011`). Because the - seed belief is uniform and ``model`` is unifilar, the retrodictive causal - states are exactly the normalized vectors :math:`(\\Pr(x \\mid s))_{s}` over - all words :math:`x`, so the reverse machine is finite iff every - log-likelihood ratio :math:`\\log \\Pr(x \\mid s) - \\log \\Pr(x \\mid s')` - takes finitely many values. + seed belief is uniform and ``model`` is unifilar, the reverse mixed states + -- the retrodictive beliefs after finitely many symbols, transient and + recurrent -- are exactly the normalized vectors + :math:`(\\Pr(x \\mid s))_{s}` over all words :math:`x`, so they are finite + iff every log-likelihood ratio + :math:`\\log \\Pr(x \\mid s) - \\log \\Pr(x \\mid s')` takes finitely many + values. + + ``True`` implies a finite reverse ε-machine. ``False`` does not imply an + infinite one: infinitely many transient beliefs may converge to finitely + many recurrent causal states. Alternating biased coins is an example -- its + likelihood ratio diverges, yet its reverse ε-machine has two states. Reading a symbol ``a`` moves the pair ``(s, s')`` to ``(delta(s, a), delta(s', a))`` and multiplies the ratio by ``p(a|s) / p(a|s')``. Placing that weight on the pair graph over ``states x states`` gives - the reverse ε-machine is finite + the reverse mixed states are finite iff every directed cycle of the pair graph has weight one. A cycle of weight :math:`\\gamma \\neq 1` traversed :math:`k` times yields @@ -82,6 +89,65 @@ def reverse_is_finite(model: StateMachine, *, rtol: float = 1e-9) -> bool: transition structure can differ, since a cycle weight can equal one by algebraic coincidence. """ + _pair_graph, inconsistent = _inconsistent_pair_components(model, rtol=rtol) + return not inconsistent + + +def reverse_epsilon_machine_is_finite(model: StateMachine, *, rtol: float = 1e-9) -> bool: + """Decide whether the reverse ε-machine of a unifilar ``model`` has finitely many recurrent states. + + .. warning:: + Experimental. The criterion below has a proof sketch and agrees with + sampled posterior limits on every strongly connected 3-state binary + machine tested, but it has not been published. + + :func:`reverse_is_finite` tests the stronger property that *every* reverse + belief after finitely many symbols is one of finitely many. A pair of + forward states ``(s, s')`` whose paths never merge has a log-likelihood + ratio that diverges almost surely, sending the posterior to a face of the + simplex; those infinitely many transient beliefs leave no recurrent trace. + Infinitely many recurrent states need a ratio that grows around a cycle and + is then *frozen* by the two paths merging, so + + the reverse ε-machine is finite + iff no strongly connected component of the pair graph that contains + a cycle of weight other than one can reach the diagonal. + + If such a component exists, traversing its cycle ``k`` times and then + merging happens with positive probability and freezes the ratio at + :math:`c \\gamma^k`, a distinct recurrent belief for each ``k``. If none + exists, merging pairs only cross components with consistent potentials, so + their frozen ratios take finitely many values. :cite:`Ellison2011`'s + explosive example and the single-atom binary witness fail it; alternating + biased coins passes it although :func:`reverse_is_finite` is ``False``. + """ + import networkx as nx + + pair_graph, inconsistent = _inconsistent_pair_components(model, rtol=rtol) + if not inconsistent: + return True + diagonal = [node for node in pair_graph.nodes if node[0] == node[1]] + can_merge = set(diagonal) + for node in diagonal: + can_merge |= nx.ancestors(pair_graph, node) + return not any(component & can_merge for component in inconsistent) + + +def warn_if_reverse_truncated(forward: StateMachine, *, stacklevel: int = 3) -> None: + """Warn when a numerically built reverse ε-machine must be a finite truncation.""" + import warnings + + if not reverse_epsilon_machine_is_finite(forward): + warnings.warn( + "the reverse ε-machine of this process is infinite; the returned machine merges " + "beliefs within numerical tolerance and is a finite approximation", + RuntimeWarning, + stacklevel=stacklevel, + ) + + +def _inconsistent_pair_components(model: StateMachine, *, rtol: float) -> tuple[Any, list[set[Any]]]: + """Return the pair graph and its strongly connected components with a non-unit cycle.""" import math import networkx as nx @@ -126,6 +192,7 @@ def agrees(left: Any, right: Any) -> bool: return is_zero(simplify_prob(left - right)) return abs(left - right) <= rtol + inconsistent: list[set[Any]] = [] for component in nx.strongly_connected_components(pair_graph): if len(component) == 1: node = next(iter(component)) @@ -141,10 +208,12 @@ def agrees(left: Any, right: Any) -> bool: if successor not in potential: potential[successor] = combine(potential[current], data["weight"]) stack.append(successor) - for source, target, data in sub.edges(data=True): - if not agrees(combine(potential[source], data["weight"]), potential[target]): - return False - return True + if any( + not agrees(combine(potential[source], data["weight"]), potential[target]) + for source, target, data in sub.edges(data=True) + ): + inconsistent.append(set(component)) + return pair_graph, inconsistent def time_reverse_stochastic(model: S) -> S: # noqa: UP047 - keep Python 3.11 compatibility. diff --git a/sofic/generators/synchronization.py b/sofic/generators/synchronization.py index a722493..2d04c5e 100644 --- a/sofic/generators/synchronization.py +++ b/sofic/generators/synchronization.py @@ -222,7 +222,14 @@ def _bellman_ford_longest_transient_path(pa: PowerAutomaton) -> float: def markov_order_from_graph(graph: TopologicalUnifilarGraph) -> int | float: - """Longest prefix-free synchronizing word length (Markov order ``R``).""" + """Longest prefix-free synchronizing word length (Markov order ``R``). + + Power-automaton subsets that can never synchronize are skipped, so the + result is the Markov order only for a minimal presentation (an + ε-machine, or a minimal right-resolving cover). On a non-minimal + presentation, equivalent states that no word separates can make it finite + even though some infinite words never fix the presentation state. + """ if not graph.states: return 0 pa = power_automaton(graph) @@ -377,7 +384,8 @@ def is_definite_from_graph(graph: TopologicalUnifilarGraph) -> bool: state is fixed by the last ``R`` symbols regardless of the start state :cite:`Cerny1964`; this coincides with finite Markov order ``R`` for a right-resolving epsilon-machine :cite:`James2010`. Definiteness implies - exact synchronizability, but not conversely. + exact synchronizability, but not conversely. Like + :func:`markov_order_from_graph`, it assumes a minimal presentation. """ return markov_order_from_graph(graph) != math.inf diff --git a/sofic/generators/words.py b/sofic/generators/words.py index f9ba213..4e0f14b 100644 --- a/sofic/generators/words.py +++ b/sofic/generators/words.py @@ -104,7 +104,7 @@ def _hmm_word_probability( for symbol in word: matrix = joint.get(symbol, zero) mass = mass @ matrix - if not np.any(np.abs(mass) > _TOL): + if not np.any(mass): return 0.0 return float(mass.sum()) @@ -174,7 +174,7 @@ def _hmm_conditional_word_probability( ) -> float: """Return ``P(word | condition)`` from the requested start distribution.""" condition_probability = _hmm_word_probability(hmm, condition, start=start) - if condition_probability <= _TOL: + if condition_probability <= 0.0: raise ZeroDivisionError("condition has zero probability") joint_word = tuple(condition) + tuple(word) return _hmm_word_probability(hmm, joint_word, start=start) / condition_probability diff --git a/sofic/inference/hmm/filtering.py b/sofic/inference/hmm/filtering.py index d30df86..e036af0 100644 --- a/sofic/inference/hmm/filtering.py +++ b/sofic/inference/hmm/filtering.py @@ -89,7 +89,7 @@ def backward(hmm: HiddenMarkovModel, observations: Sequence[Any], *, normalize: n = next(iter(joint.values())).shape[0] if joint else len(_as_mealy_hmm(hmm).reindex()) obs = list(observations) beta = np.zeros((len(obs) + 1, n), dtype=float) - beta[len(obs)] = 1.0 + beta[len(obs)] = 1.0 / n if normalize and n else 1.0 for t in range(len(obs) - 1, -1, -1): matrix = joint.get(obs[t]) if matrix is None: diff --git a/sofic/serialization.py b/sofic/serialization.py index c91e429..04e950e 100644 --- a/sofic/serialization.py +++ b/sofic/serialization.py @@ -30,6 +30,10 @@ VERSION = 1 _TYPE_KEY = "__sofic_type__" +# libyaml's C implementations are ~10x faster on large models; same safe subset. +_SafeLoader = getattr(yaml, "CSafeLoader", yaml.SafeLoader) +_SafeDumper = getattr(yaml, "CSafeDumper", yaml.SafeDumper) + @dataclass(frozen=True, slots=True) class _ModelSpec: @@ -40,12 +44,12 @@ class _ModelSpec: def model_to_yaml(model: StateMachine) -> str: """Return a YAML representation of ``model``.""" - return yaml.safe_dump(model_to_dict(model), sort_keys=False) + return yaml.dump(model_to_dict(model), Dumper=_SafeDumper, sort_keys=False) def model_from_yaml(text: str, *, validate: bool = True) -> StateMachine: """Reconstruct a sofic model from YAML text.""" - loaded = yaml.safe_load(text) + loaded = yaml.load(text, Loader=_SafeLoader) if not isinstance(loaded, dict): raise TypeError("sofic model YAML must load to a mapping") return model_from_dict(loaded, validate=validate) diff --git a/sofic/shifts/cover_construction.py b/sofic/shifts/cover_construction.py index ddbf5f1..ceec54e 100644 --- a/sofic/shifts/cover_construction.py +++ b/sofic/shifts/cover_construction.py @@ -142,10 +142,46 @@ def right_fischer_cover(shift: SoficShift) -> RightFischerCover: members = set(condensation.nodes[terminal[0]]["members"]) if not any(classes[s] in members and delta[s] for s in delta): return RightFischerCover(symbol_alphabet=shift.symbol_alphabet) + if not _presents_every_word(frozenset(trimmed.states()), delta, classes, members): + raise SoficValidationError("the Fischer cover is defined for irreducible sofic shifts; this shift is reducible") vertices = [subset for subset in delta if classes[subset] in members] return _quotient_shift(RightFischerCover, shift, vertices, delta, classes) +def _presents_every_word( + start: Subset, + delta: dict[Subset, dict[Any, Subset]], + classes: dict[Subset, int], + members: set[int], +) -> bool: + """Whether the terminal component reads every word the whole presentation reads. + + A unique terminal component does not make the shift irreducible (``0^inf``, + ``1^inf`` and ``0...01...1`` have one), but the shift is irreducible exactly + when that component presents all of it. Both sides are deterministic, so + walk their product from (all vertices, all component classes). + """ + moves: dict[int, dict[Any, int]] = defaultdict(dict) + for subset, edges in delta.items(): + if classes[subset] in members: + for symbol, target in edges.items(): + moves[classes[subset]][symbol] = classes[target] + initial = (start, frozenset(members)) + seen = {initial} + queue = deque([initial]) + while queue: + subset, component = queue.popleft() + for symbol, target in delta[subset].items(): + following = frozenset(moves[c][symbol] for c in component if symbol in moves[c]) + if not following: + return False + pair = (target, following) + if pair not in seen: + seen.add(pair) + queue.append(pair) + return True + + def left_fischer_cover(shift: SoficShift) -> LeftFischerCover: """Return the minimal left-resolving presentation (mirror of the right Fischer cover).""" return _mirror(LeftFischerCover, right_fischer_cover(shift.reverse())) diff --git a/sofic/shifts/topological_anatomy.py b/sofic/shifts/topological_anatomy.py index a2f07a9..13ace73 100644 --- a/sofic/shifts/topological_anatomy.py +++ b/sofic/shifts/topological_anatomy.py @@ -107,7 +107,8 @@ def _right_resolving(shift: SoficShift) -> SoficShift: def parry_measure_sofic(shift: SoficShift) -> MealyHMM: """Return the measure of maximal entropy (Parry measure) as a labeled MealyHMM. - Puts ``shift`` in right-resolving form (:func:`_right_resolving`) and builds + Drops states off every bi-infinite path, puts ``shift`` in right-resolving + form (:func:`_right_resolving`), and builds the Parry chain ``P[i, j] = A[i, j] v_j / (lambda v_i)`` with stationary ``pi ~ u * v`` from the Perron data of the adjacency matrix ``A`` (Parry 1964), preserving each edge's emitted symbol. Its entropy rate is @@ -116,7 +117,7 @@ def parry_measure_sofic(shift: SoficShift) -> MealyHMM: """ from sofic.shifts.parry_construction import parry_measure - return parry_measure(_right_resolving(shift)) + return parry_measure(_right_resolving(shift.trim_transient())) def topological_anatomy(shift: SoficShift) -> dict[str, float]: diff --git a/sofic/viz/_tikz_format.py b/sofic/viz/_tikz_format.py index 2eae134..6b819d7 100644 --- a/sofic/viz/_tikz_format.py +++ b/sofic/viz/_tikz_format.py @@ -23,6 +23,9 @@ "_": r"\_", "^": r"\mbox{\^{}}", "~": r"\mbox{\textasciitilde}", + # A blank line is ``\par``, which breaks node text and math mode alike. + "\n": " ", + "\r": " ", } @@ -38,7 +41,7 @@ def _latex_arg(text: str) -> str: """Escape content for a LaTeX macro argument (no math mode).""" if text == "": return "{}" - if any(char in text for char in "{}\\#%&_^~$"): + if any(char in _LATEX_SPECIAL for char in text): return f"{{{latex_escape(text)}}}" return text diff --git a/tests/test_epsilon_machine.py b/tests/test_epsilon_machine.py index 864f236..c9f4a31 100644 --- a/tests/test_epsilon_machine.py +++ b/tests/test_epsilon_machine.py @@ -201,3 +201,62 @@ def test_from_time_reversed_propagates_explosion_instead_of_degrading(): pytest.raises(MixedStateExplosionError), ): EpsilonMachine.from_time_reversed(forward) + + +def test_reverse_is_finite_false_does_not_mean_infinite_reverse_machine(): + """Infinitely many transient reverse beliefs can converge to finitely many causal states. + + Alternating biased coins has a likelihood ratio that diverges around a + pair-graph cycle, so the twins test fails, yet its reverse ε-machine has two + states. + """ + from sofic.examples import alternating_biased_coins + + forward = alternating_biased_coins() + assert not forward.reverse_is_finite() + assert forward.reverse_epsilon_machine_is_finite() + reverse = EpsilonMachine.from_time_reversed(forward) + assert len(list(reverse.states())) == 2 + + +def test_reverse_epsilon_machine_is_finite_on_known_cases(): + from sofic.examples import even_process, golden_mean + + assert not _explosive_forward().reverse_epsilon_machine_is_finite() + assert golden_mean(0.5).reverse_epsilon_machine_is_finite() + assert even_process(0.5).reverse_epsilon_machine_is_finite() + assert irreversible_two_state().reverse_epsilon_machine_is_finite() + + +def test_reverse_is_finite_implies_reverse_epsilon_machine_is_finite(): + from sofic.generators.topological_epsilon_enumeration import iter_topological_epsilon_machines + + for machine in iter_topological_epsilon_machines(2, 3, alphabet=("0", "1")): + if machine.reverse_is_finite(): + assert machine.reverse_epsilon_machine_is_finite() + + +def test_numerically_truncated_reverse_machine_warns(): + """Numeric beliefs merging within tolerance must not pass for a finite reverse machine. + + The pair ``(0, 1)`` gains a factor of two per ``1`` before the paths merge, so + the frozen posterior ratios are ``2^k`` -- infinitely many recurrent states -- + yet tolerance merging closes the belief set. + """ + import pytest + + spec = { + "A": [("1", "B", 1.0)], + "B": [("0", "A", 0.5), ("1", "C", 0.5)], + "C": [("0", "B", 0.5), ("1", "C", 0.5)], + } + forward = EpsilonMachine(initial_distribution={"A": 1.0}, observation_alphabet=frozenset({"0", "1"})) + for state in spec: + forward.graph.add_state(state) + for state, edges in spec.items(): + for symbol, target, prob in edges: + forward.graph.add_transition(state, target, **{ATTR_PROB: prob, ATTR_EMISSION: symbol}) + + assert not forward.reverse_epsilon_machine_is_finite() + with pytest.warns(RuntimeWarning, match="finite approximation"): + EpsilonMachine.from_time_reversed(forward) diff --git a/tests/test_properties_automata.py b/tests/test_properties_automata.py new file mode 100644 index 0000000..db0b560 --- /dev/null +++ b/tests/test_properties_automata.py @@ -0,0 +1,270 @@ +"""Property-based and metamorphic tests for finite-automaton constructions. + +Every language is compared against :mod:`tests.oracles` on all words up to a +small length; the automata are tiny enough that this is exhaustive in practice. +""" + +from __future__ import annotations + +import re +from collections.abc import Callable, Hashable +from typing import Any + +from hypothesis import given +from hypothesis import strategies as st + +from sofic.automata import ( + DFA, + NFA, + CanonicalRFSA, + automaton_to_regex, + complete, + determinize, + equivalent, + minimize, +) +from sofic.automata.languages import atoms, left_quotient, right_quotient +from sofic.graph import ATTR_SYMBOL, EPSILON +from sofic.testing import dfas, nfas +from tests import oracles + +AB = ("0", "1") +SIGMA = frozenset(AB) +N = 5 + +word = st.lists(st.sampled_from(AB), max_size=3).map(tuple) + + +def lang(aut: Any, n: int = N, alphabet: tuple[Any, ...] = AB) -> frozenset[tuple[Any, ...]]: + return frozenset(w for w in oracles.words(alphabet, n) if aut.recognizes(w)) + + +def truth(aut: Any, n: int = N, alphabet: tuple[Any, ...] = AB) -> frozenset[tuple[Any, ...]]: + return oracles.language(lambda w: oracles.nfa_accepts(aut, w), alphabet, n) + + +def relabel( + aut: Any, + state_map: Callable[[Hashable], Hashable], + symbol_map: Callable[[Any], Any], +) -> Any: + result = type(aut)( + input_alphabet=frozenset(symbol_map(a) for a in aut.input_alphabet), + initial_states=frozenset(state_map(s) for s in aut.initial_states), + accepting_states=frozenset(state_map(s) for s in aut.accepting_states), + ) + for state in aut.states(): + result.graph.add_state(state_map(state)) + for t in aut.transitions(): + symbol = t.data.get(ATTR_SYMBOL) + mapped = symbol if symbol is EPSILON else symbol_map(symbol) + result.graph.add_transition(state_map(t.source), state_map(t.target), **{ATTR_SYMBOL: mapped}) + return result + + +def epsilon_nfa() -> NFA: + nfa = NFA(input_alphabet=SIGMA, initial_states=frozenset({"e"}), accepting_states=frozenset({"e"})) + nfa.graph.add_state("e") + return nfa + + +def has_empty_residual(accepts: Callable[[tuple[Any, ...]], bool], k: int) -> bool: + suffixes = oracles.words(AB, k) + return any(not any(accepts(u + v) for v in suffixes) for u in oracles.words(AB, k)) + + +def star_member(language: frozenset[tuple[Any, ...]], w: tuple[Any, ...]) -> bool: + reach = [True] + [False] * len(w) + for end in range(1, len(w) + 1): + reach[end] = any(reach[start] and w[start:end] in language for start in range(end)) + return reach[len(w)] + + +# --------------------------------------------------------------------------- Boolean / regular operations + + +@given(nfas()) +def test_recognizes_matches_oracle(a: NFA) -> None: + assert lang(a) == truth(a) + + +@given(nfas(), nfas()) +def test_boolean_operations_match_brute_force(a: NFA, b: NFA) -> None: + la, lb = truth(a), truth(b) + everything = frozenset(oracles.words(AB, N)) + assert lang(a.union(b)) == la | lb + assert lang(a.intersection(b)) == la & lb + assert lang(a.complement(SIGMA)) == everything - la + assert lang(a.difference(b, SIGMA)) == la - lb + + +@given(nfas(), nfas()) +def test_concat_matches_brute_force(a: NFA, b: NFA) -> None: + la, lb = truth(a), truth(b) + expected = frozenset(w for w in oracles.words(AB, N) if any(w[:i] in la and w[i:] in lb for i in range(len(w) + 1))) + assert lang(a.concat(b)) == expected + + +@given(nfas()) +def test_kleene_star_matches_brute_force(a: NFA) -> None: + la = truth(a) + assert lang(a.kleene_star()) == frozenset(w for w in oracles.words(AB, N) if star_member(la, w)) + + +@given(nfas()) +def test_reverse_matches_brute_force(a: NFA) -> None: + assert lang(a.reverse()) == frozenset(w[::-1] for w in truth(a)) + + +@given(nfas(), word) +def test_quotients_match_brute_force(a: NFA, u: tuple[Any, ...]) -> None: + left, right = left_quotient(u, a), right_quotient(a, u) + for w in oracles.words(AB, 4): + expected_left = oracles.nfa_accepts(a, u + w) + expected_right = oracles.nfa_accepts(a, w + u) + assert (w in left) == expected_left + assert (w in right) == expected_right + + +# --------------------------------------------------------------------------- determinize / minimize / equivalent + + +@given(nfas()) +def test_determinize_preserves_language(a: NFA) -> None: + dfa = determinize(a) + assert dfa.is_deterministic() + assert lang(dfa) == truth(a) + + +@given(nfas()) +def test_minimization_algorithms_agree(a: NFA) -> None: + results = [minimize(a, algorithm=algorithm) for algorithm in ("hopcroft", "moore", "brzozowski")] + sizes = {len(list(m.states())) for m in results} + assert len(sizes) == 1 + expected = truth(a) + for m in results: + assert lang(m) == expected + assert equivalent(m, results[0]) + + +@given(nfas(max_states=2) | dfas()) +def test_minimal_size_is_brute_force_nerode_count(a: Any) -> None: + def accepts(w: tuple[Any, ...]) -> bool: + return oracles.nfa_accepts(a, w) + + k = 3 + expected = oracles.nerode_class_count(accepts, AB, k) - has_empty_residual(accepts, k) + assert len(list(minimize(a).states())) == expected + + +@given(nfas()) +def test_minimize_is_idempotent(a: NFA) -> None: + once = minimize(a) + twice = minimize(once) + assert len(list(twice.states())) == len(list(once.states())) + assert equivalent(once, twice) + + +@given(nfas(max_states=2), nfas(max_states=2)) +def test_equivalent_matches_brute_force(a: NFA, b: NFA) -> None: + # Complete DFAs here have <= 4 states, so a shortest distinguishing word has length <= 6. + assert equivalent(a, b) == (truth(a, 6) == truth(b, 6)) + + +@given(nfas()) +def test_equivalent_to_own_constructions(a: NFA) -> None: + assert equivalent(a, determinize(a)) + assert equivalent(a, minimize(a, algorithm="brzozowski")) + assert equivalent(a, a.union(a)) + + +# --------------------------------------------------------------------------- regex, RFSA, atoms + + +@given(nfas()) +def test_regex_round_trip_via_python_re(a: NFA) -> None: + pattern = re.compile(automaton_to_regex(a)) + for w in oracles.words(AB, N): + assert (pattern.fullmatch("".join(w)) is not None) == oracles.nfa_accepts(a, w) + + +@given(nfas()) +def test_canonical_rfsa_is_no_larger_than_minimal_dfa(a: NFA) -> None: + rfsa = CanonicalRFSA.from_language(a) + assert len(list(rfsa.states())) <= len(list(minimize(a).states())) + assert lang(rfsa) == truth(a) + + +@given(nfas()) +def test_atoms_partition_sigma_star(a: NFA) -> None: + pieces = list(atoms(a)) + for w in oracles.words(AB, N): + assert sum(w in atom for atom in pieces) == 1 + + +@given(nfas()) +def test_words_in_one_atom_share_quotient_signature(a: NFA) -> None: + prefixes = oracles.words(AB, 2) + for atom in atoms(a): + members = [w for w in oracles.words(AB, 4) if w in atom] + signatures = {tuple(oracles.nfa_accepts(a, p + w) for p in prefixes) for w in members} + assert len(signatures) <= 1 + + +# --------------------------------------------------------------------------- metamorphic + + +@given(nfas()) +def test_relabeling_states_and_symbols_maps_language(a: NFA) -> None: + to_letter = {"0": "x", "1": "y"} + renamed = relabel(a, lambda s: ("q", s), to_letter.__getitem__) + expected = frozenset(tuple(to_letter[c] for c in w) for w in truth(a)) + assert lang(renamed, alphabet=("x", "y")) == expected + assert len(list(minimize(renamed).states())) == len(list(minimize(a).states())) + + +@given(nfas(), st.lists(st.tuples(st.integers(0, 2), st.sampled_from(AB)), max_size=4), st.booleans()) +def test_adding_unreachable_states_preserves_language(a: NFA, edges: list[tuple[int, str]], accepting: bool) -> None: + padded = relabel(a, lambda s: s, lambda c: c) + states = list(a.states()) + padded.graph.add_state("junk") + if accepting: + padded.accepting_states = padded.accepting_states | {"junk"} + for index, symbol in edges: + padded.graph.add_transition("junk", states[index % len(states)], **{ATTR_SYMBOL: symbol}) + assert lang(padded) == truth(a) + assert len(list(minimize(padded).states())) == len(list(minimize(a).states())) + + +@given(nfas()) +def test_double_reversal_minimizes_to_same_dfa(a: NFA) -> None: + left = minimize(determinize(a.reverse().reverse())) + right = minimize(a) + assert len(list(left.states())) == len(list(right.states())) + assert equivalent(left, right) + + +@given(nfas()) +def test_concat_with_epsilon_is_identity(a: NFA) -> None: + expected = truth(a) + assert lang(a.concat(epsilon_nfa())) == expected + assert lang(epsilon_nfa().concat(a)) == expected + + +@given(nfas() | dfas()) +def test_double_complement_is_identity(a: NFA | DFA) -> None: + assert lang(a.complement(SIGMA).complement(SIGMA)) == truth(a) + + +def test_double_complement_regression_accepting_trap() -> None: + # The first complement's trap state is accepting; re-completing must not route missing edges into it. + a = NFA(input_alphabet=SIGMA, initial_states=frozenset({1}), accepting_states=frozenset({0})) + a.graph.add_state(0) + a.graph.add_state(1) + a.add_transition(0, 0, "0") + a.add_transition(0, 0, "1") + a.add_transition(1, 0, "0") + once = a.complement(SIGMA) + assert lang(once) == frozenset(oracles.words(AB, N)) - truth(a) + assert lang(once.complement(SIGMA)) == truth(a) + assert lang(complete(once, SIGMA)) == lang(once) diff --git a/tests/test_properties_automata_learning.py b/tests/test_properties_automata_learning.py new file mode 100644 index 0000000..d79abe3 --- /dev/null +++ b/tests/test_properties_automata_learning.py @@ -0,0 +1,164 @@ +"""Property-based and metamorphic tests for automaton learners and observation tables.""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Callable, Hashable +from itertools import combinations +from typing import Any + +from hypothesis import given +from hypothesis import strategies as st + +from sofic.automata import ( + DFA, + AutomatonEquivalenceOracle, + CanonicalRFSA, + LanguageMembershipOracle, + ObservationTable, + equivalent, + learn_dfa_edsm, + learn_dfa_lstar, + learn_dfa_rpni, + learn_dfa_ttt, + learn_rfsa_nlstar, +) +from sofic.testing import dfas +from tests import oracles +from tests.test_properties_automata import relabel + +AB = ("0", "1") +PRESERVING = {"0": "a", "1": "b"} +REVERSING = {"0": "b", "1": "a"} + +Word = tuple[Any, ...] + + +def lang(aut: Any, alphabet: tuple[Any, ...] = AB, n: int = 5) -> frozenset[Word]: + return frozenset(w for w in oracles.words(alphabet, n) if aut.recognizes(w)) + + +def mapped(words: frozenset[Word], symbol_map: dict[str, str]) -> frozenset[Word]: + return frozenset(tuple(symbol_map[c] for c in w) for w in words) + + +def size(aut: Any) -> int: + return len(list(aut.states())) + + +def nerode(dfa: DFA) -> int: + return oracles.nerode_class_count(lambda w: oracles.nfa_accepts(dfa, w), AB, 3) + + +def active(learner: Callable[..., Any], target: Any, alphabet: tuple[str, ...]) -> Any: + return learner(alphabet, LanguageMembershipOracle(target), AutomatonEquivalenceOracle(target, alphabet)) + + +# --------------------------------------------------------------------------- active learners + + +@given(dfas(), st.sampled_from([learn_dfa_lstar, learn_dfa_ttt])) +def test_active_dfa_learners_find_the_minimal_dfa(target: DFA, learner: Callable[..., DFA]) -> None: + learned = active(learner, target, AB) + assert equivalent(learned, target) + assert size(learned) <= nerode(target) + + +@given(dfas(), st.sampled_from([learn_dfa_lstar, learn_dfa_ttt]), st.sampled_from([PRESERVING, REVERSING])) +def test_active_dfa_learners_are_relabeling_invariant( + target: DFA, learner: Callable[..., DFA], symbol_map: dict[str, str] +) -> None: + renamed = relabel(target, lambda s: ("r", s), symbol_map.__getitem__) + original = active(learner, target, AB) + learned = active(learner, renamed, ("a", "b")) + assert size(learned) == size(original) + assert lang(learned, ("a", "b")) == mapped(lang(original), symbol_map) + + +@given(dfas(), st.sampled_from([PRESERVING, REVERSING])) +def test_nlstar_learns_canonical_rfsa_invariantly(target: DFA, symbol_map: dict[str, str]) -> None: + learned = active(learn_rfsa_nlstar, target, AB) + assert lang(learned) == lang(target) + assert size(learned) == size(CanonicalRFSA.from_language(target)) + renamed = active(learn_rfsa_nlstar, relabel(target, lambda s: ("r", s), symbol_map.__getitem__), ("a", "b")) + assert size(renamed) == size(learned) + assert lang(renamed, ("a", "b")) == mapped(lang(learned), symbol_map) + + +# --------------------------------------------------------------------------- passive learners + + +@given( + dfas(), + st.sets(st.sampled_from(oracles.words(AB, 3)), min_size=1), + st.sampled_from([learn_dfa_rpni, learn_dfa_edsm]), +) +def test_passive_learners_are_consistent_and_relabeling_invariant( + target: DFA, sample: set[Word], learner: Callable[..., DFA] +) -> None: + ordered = sorted(sample, key=lambda w: (len(w), w)) + positive = [w for w in ordered if target.recognizes(w)] + negative = [w for w in ordered if not target.recognizes(w)] + learned = learner(positive, negative) + assert all(learned.recognizes(w) for w in positive) + assert not any(learned.recognizes(w) for w in negative) + + def rename(words: list[Word]) -> list[Word]: + return [tuple(PRESERVING[c] for c in w) for w in words] + + renamed = learner(rename(positive), rename(negative)) + assert size(renamed) == size(learned) + assert lang(renamed, ("a", "b"), 4) == mapped(lang(learned, AB, 4), PRESERVING) + + +# --------------------------------------------------------------------------- observation tables + + +def union_witness(target: DFA, state: Hashable, others: frozenset[Hashable]) -> Word | None: + """Shortest word in the right language of ``state`` but of none of ``others`` (complete DFA).""" + start = (state, others) + queue, seen = deque([(start, ())]), {start} + while queue: + (p, qs), w = queue.popleft() + if p in target.accepting_states and not qs & target.accepting_states: + return w + for a in AB: + nxt = (next(iter(target.delta(p, a))), frozenset(next(iter(target.delta(q, a))) for q in qs)) + if nxt not in seen: + seen.add(nxt) + queue.append((nxt, (*w, a))) + return None + + +def full_table(target: DFA, n: int = 3) -> ObservationTable: + # Extraction reads transitions u -> ua only among access words, so access must hold every word < n and + # its one-symbol extensions. Experiments must witness every residual non-inclusion L_p not in the union + # of L_q (q in Q), or the table is closed and consistent but not RFSA-closed and RFSA-consistent. + access = frozenset(oracles.words(AB, n)) + states = list(target.states()) + witnesses = { + union_witness(target, p, frozenset(others)) + for p in states + for k in range(len(states)) + for others in combinations([q for q in states if q != p], k) + } + experiments = frozenset(oracles.words(AB, 2)) | {w[i:] for w in witnesses if w is not None for i in range(len(w))} + membership = { + u + a + e: oracles.nfa_accepts(target, u + a + e) + for u in access + for a in ((), *((s,) for s in AB)) + for e in experiments + } + return ObservationTable(access_words=access, experiments=experiments, membership=membership) + + +@given(dfas()) +def test_observation_table_extractions(target: DFA) -> None: + table = full_table(target) + dfa = table.to_minimal_dfa() + assert lang(dfa) == lang(target) + assert size(dfa) == nerode(target) + rfsa = table.to_canonical_rfsa() + assert lang(rfsa) == lang(target) + assert size(rfsa) == size(CanonicalRFSA.from_language(target)) + assert lang(table.to_atomaton()) == lang(target) diff --git a/tests/test_properties_automata_omega_wheeler.py b/tests/test_properties_automata_omega_wheeler.py new file mode 100644 index 0000000..f643550 --- /dev/null +++ b/tests/test_properties_automata_omega_wheeler.py @@ -0,0 +1,146 @@ +"""Property-based tests for Büchi automata and Wheeler automata / indexes.""" + +from __future__ import annotations + +from collections.abc import Hashable, Sequence +from itertools import permutations +from typing import Any + +from hypothesis import given +from hypothesis import strategies as st + +from sofic.automata import ( + NFA, + BuchiAutomaton, + determinize_wheeler, + is_wheeler, + minimize, + minimize_wheeler, + wheeler_index, + wheeler_order, +) +from sofic.graph import ATTR_SYMBOL +from sofic.testing import buchi_automata, lassos, nfas, wheeler_nfas +from tests import oracles + +AB = ("0", "1") +N = 5 + + +def lang(aut: Any, n: int = N) -> frozenset[tuple[Any, ...]]: + return frozenset(w for w in oracles.words(AB, n) if aut.recognizes(w)) + + +def truth(aut: Any, n: int = N) -> frozenset[tuple[Any, ...]]: + return oracles.language(lambda w: oracles.nfa_accepts(aut, w), AB, n) + + +def edges(aut: Any) -> list[tuple[Hashable, Any, Hashable]]: + return [(t.source, t.data.get(ATTR_SYMBOL), t.target) for t in aut.transitions()] + + +def satisfies_wheeler_axioms(aut: Any, order: Sequence[Hashable]) -> bool: + """Gagie-Manzini-Sirén axioms, with labels ordered by ``repr``.""" + rank = {state: index for index, state in enumerate(order)} + arcs = edges(aut) + entered = {target for _source, _symbol, target in arcs} + sources = [rank[s] for s in rank if s not in entered] + if sources and any(rank[t] < max(sources) for t in entered): + return False + for u, a, u2 in arcs: + for v, b, v2 in arcs: + if repr(a) < repr(b) and not rank[u2] < rank[v2]: + return False + if a == b and rank[u] < rank[v] and not rank[u2] <= rank[v2]: + return False + return True + + +def admits_wheeler_order(aut: Any) -> bool: + return any(satisfies_wheeler_axioms(aut, order) for order in permutations(aut.states())) + + +# --------------------------------------------------------------------------- Büchi + + +@given(buchi_automata() | buchi_automata(allow_epsilon=True), lassos()) +def test_accepts_lasso_matches_oracle(ba: BuchiAutomaton, lasso: tuple[tuple[Any, ...], tuple[Any, ...]]) -> None: + prefix, loop = lasso + assert ba.accepts_lasso(prefix, loop) == oracles.buchi_accepts_lasso(ba, prefix, loop) + + +@given(buchi_automata(), lassos(), st.integers(0, 3)) +def test_accepts_lasso_invariant_under_unrolling_and_rotation( + ba: BuchiAutomaton, lasso: tuple[tuple[Any, ...], tuple[Any, ...]], shift: int +) -> None: + prefix, loop = lasso + expected = ba.accepts_lasso(prefix, loop) + k = shift % len(loop) + assert ba.accepts_lasso(prefix + loop, loop) == expected + assert ba.accepts_lasso(prefix, loop + loop) == expected + assert ba.accepts_lasso(prefix + loop[:k], loop[k:] + loop[:k]) == expected + + +# --------------------------------------------------------------------------- Wheeler order + + +@given(nfas(allow_epsilon=False, max_states=4) | wheeler_nfas()) +def test_wheeler_order_is_sound_and_complete(a: NFA) -> None: + order = wheeler_order(a) + exists = admits_wheeler_order(a) + assert (order is not None) == exists + assert is_wheeler(a) == exists + if order is not None: + assert sorted(order.states, key=repr) == sorted(a.states(), key=repr) + assert satisfies_wheeler_axioms(a, order.states) + + +@given(wheeler_nfas()) +def test_wheeler_strategy_generates_wheeler_nfas(a: NFA) -> None: + assert wheeler_order(a) is not None + + +# --------------------------------------------------------------------------- Wheeler determinization / minimization + + +@given(wheeler_nfas()) +def test_determinize_wheeler_preserves_language_and_wheelerness(a: NFA) -> None: + dfa = determinize_wheeler(a) + assert dfa.is_deterministic() + assert lang(dfa) == truth(a) + assert is_wheeler(dfa) + + +@given(wheeler_nfas()) +def test_minimize_wheeler_preserves_language_and_is_idempotent(a: NFA) -> None: + dfa = determinize_wheeler(a) + smallest = minimize_wheeler(dfa) + assert lang(smallest) == truth(a) + assert is_wheeler(smallest) + assert len(list(minimize(a).states())) <= len(list(smallest.states())) <= len(list(dfa.states())) + assert len(list(minimize_wheeler(smallest).states())) == len(list(smallest.states())) + + +# --------------------------------------------------------------------------- Wheeler index + + +@given(wheeler_nfas()) +def test_wheeler_index_membership_and_reachability(a: NFA) -> None: + index = wheeler_index(a) + for w in oracles.words(AB, N): + assert index.contains(w) == oracles.nfa_accepts(a, w) + assert frozenset(index.states_reached(w)) == oracles.nfa_reachable(a, w) + assert index.count_states(w) == len(oracles.nfa_reachable(a, w)) + + +@given(wheeler_nfas(), st.integers(0, N)) +def test_wheeler_index_counts_ranks_and_enumerates(a: NFA, n: int) -> None: + index = wheeler_index(a) + expected = frozenset(w for w in truth(a, n) if len(w) == n) + assert index.count_words(n) == len(expected) + listed = list(index.words_of_length(n)) + assert frozenset(listed) == expected + assert len(listed) == len(expected) + for i, w in enumerate(listed): + assert index.rank_word(w) == i + assert index.unrank_word(i, n) == w diff --git a/tests/test_properties_automata_transducers.py b/tests/test_properties_automata_transducers.py new file mode 100644 index 0000000..3977e5b --- /dev/null +++ b/tests/test_properties_automata_transducers.py @@ -0,0 +1,319 @@ +"""Property-based tests for transducers, transducer/generator operations, and unifilar automata.""" + +from __future__ import annotations + +import math +from collections.abc import Hashable +from itertools import product +from typing import Any + +import pytest +from hypothesis import assume, given +from hypothesis import strategies as st + +from sofic.automata import ( + MealyMachine, + SubsequentialTransducer, + UnifilarAutomaton, + WeightedFiniteStateTransducer, + compose_transducer_generator, + compose_transducers, + generator_product, + transduce_generator, + transducer_product, +) +from sofic.generators.mealy import MealyHMM +from sofic.graph import ATTR_OUTPUT, ATTR_PROB, ATTR_SYMBOL, EPSILON +from sofic.testing import dfas, mealy_hmms, mealy_transducers +from tests import oracles + +AB = ("0", "1") +OUT = ("a", "b") +FINAL = ("x", "y") + +Word = tuple[Any, ...] + + +@st.composite +def stochastic_transducers( + draw: Any, inputs: tuple[str, ...] = AB, outputs: tuple[str, ...] = OUT, max_states: int = 2 +) -> MealyMachine: + """Complete, epsilon-free Mealy transducers with ``P(. | state, input)`` summing to one.""" + n = draw(st.integers(1, max_states)) + machine = MealyMachine( + input_alphabet=frozenset(inputs), + output_alphabet=frozenset(outputs), + initial_states=draw(st.frozensets(st.integers(0, n - 1), min_size=1)), + ) + for state in range(n): + machine.graph.add_state(state) + edge = st.tuples(st.integers(0, n - 1), st.sampled_from(outputs)) + for state in range(n): + for symbol in inputs: + chosen = draw(st.lists(edge, min_size=1, max_size=2, unique=True)) + weights = draw(st.lists(st.integers(1, 5), min_size=len(chosen), max_size=len(chosen))) + for (target, output), weight in zip(chosen, weights, strict=True): + machine.add_transition(state, target, symbol, output, prob=weight / sum(weights)) + return machine + + +def conditional(machine: MealyMachine, x: Word, y: Word) -> float: + """``P(y | x)`` by path enumeration, starting uniformly over the initial states.""" + + def paths(state: Hashable, position: int) -> float: + if position == len(x): + return 1.0 + total = 0.0 + for t in machine.graph.out_transitions(state): + if t.data.get(ATTR_SYMBOL) == x[position] and t.data.get(ATTR_OUTPUT) == y[position]: + total += float(t.data.get(ATTR_PROB, 1.0)) * paths(t.target, position + 1) + return total + + starts = tuple(machine.initial_states) + return sum(paths(state, 0) for state in starts) / len(starts) + + +def path_weights(machine: MealyMachine, x: Word, y: Word) -> list[float]: + """Probability of every complete path reading ``x`` and writing ``y``.""" + found: list[float] = [] + + def walk(state: Hashable, position: int, weight: float) -> None: + if position == len(x): + found.append(weight) + return + for t in machine.graph.out_transitions(state): + if t.data.get(ATTR_SYMBOL) == x[position] and t.data.get(ATTR_OUTPUT) == y[position]: + walk(t.target, position + 1, weight * float(t.data.get(ATTR_PROB, 1.0))) + + for state in machine.initial_states: + walk(state, 0, 1.0) + return found + + +def close(left: float, right: float) -> bool: + return math.isclose(left, right, rel_tol=1e-9, abs_tol=1e-12) + + +# --------------------------------------------------------------------------- relations + + +@given(mealy_transducers(), st.lists(st.sampled_from(AB), max_size=5).map(tuple)) +def test_transduce_matches_oracle(machine: MealyMachine, x: Word) -> None: + assert machine.transduce(x) == oracles.transducer_outputs(machine, x) + + +@given( + mealy_transducers(), + mealy_transducers(input_alphabet=OUT, output_alphabet=FINAL), +) +def test_compose_transducers_is_relational_composition(first: MealyMachine, second: MealyMachine) -> None: + inputs = oracles.words(AB, 4) + composed = compose_transducers((first, second), complete=False) + expected = oracles.compose_relations( + oracles.transducer_relation(first, inputs), lambda y: oracles.transducer_outputs(second, y) + ) + assert oracles.transducer_relation(composed, inputs) == expected + + +@given( + mealy_transducers(allow_epsilon_output=False), + mealy_transducers(input_alphabet=("p", "q"), output_alphabet=FINAL, allow_epsilon_output=False), +) +def test_transducer_product_is_componentwise(first: MealyMachine, second: MealyMachine) -> None: + joint = transducer_product((first, second)) + for n in range(4): + for x1, x2 in product(oracles.words(AB, n, min_length=n), oracles.words(("p", "q"), n, min_length=n)): + expected = { + tuple(zip(y1, y2, strict=True)) + for y1 in oracles.transducer_outputs(first, x1) + for y2 in oracles.transducer_outputs(second, x2) + } + assert joint.transduce(tuple(zip(x1, x2, strict=True))) == expected + + +@given(stochastic_transducers(), stochastic_transducers(inputs=OUT, outputs=FINAL)) +def test_compose_transducers_multiplies_conditionals(first: MealyMachine, second: MealyMachine) -> None: + composed = compose_transducers((first, second)) + for n in range(4): + middles = oracles.words(OUT, n, min_length=n) + for x in oracles.words(AB, n, min_length=n): + for z in oracles.words(FINAL, n, min_length=n): + expected = sum(conditional(first, x, y) * conditional(second, y, z) for y in middles) + assert close(conditional(composed, x, z), expected) + + +# --------------------------------------------------------------------------- generator pushforwards + + +@given(stochastic_transducers(), mealy_hmms(max_states=2), st.integers(0, 3)) +def test_transduce_generator_is_pushforward(machine: MealyMachine, generator: MealyHMM, n: int) -> None: + output = transduce_generator(machine, generator) + inputs = oracles.words(AB, n, min_length=n) + for y in oracles.words(OUT, n, min_length=n): + expected = sum(oracles.word_probability(generator, x) * conditional(machine, x, y) for x in inputs) + assert close(oracles.word_probability(output, y), expected) + + +@given(stochastic_transducers(), mealy_hmms(max_states=2), st.integers(0, 3)) +def test_compose_transducer_generator_is_joint_law(machine: MealyMachine, generator: MealyHMM, n: int) -> None: + joint = compose_transducer_generator(machine, generator) + for x in oracles.words(AB, n, min_length=n): + px = oracles.word_probability(generator, x) + for y in oracles.words(OUT, n, min_length=n): + expected = px * conditional(machine, x, y) + assert close(oracles.word_probability(joint, tuple(zip(x, y, strict=True))), expected) + + +@given(mealy_hmms(max_states=2), mealy_hmms(alphabet=OUT, max_states=2), st.integers(0, 3)) +def test_generator_product_is_independent_coupling(left: MealyHMM, right: MealyHMM, n: int) -> None: + joint = generator_product((left, right)) + for a in oracles.words(AB, n, min_length=n): + pa = oracles.word_probability(left, a) + for b in oracles.words(OUT, n, min_length=n): + expected = pa * oracles.word_probability(right, b) + assert close(oracles.word_probability(joint, tuple(zip(a, b, strict=True))), expected) + + +@given( + stochastic_transducers(max_states=1), + stochastic_transducers(inputs=OUT, outputs=FINAL, max_states=2), + mealy_hmms(max_states=2), +) +def test_transducing_twice_equals_transducing_by_composition( + first: MealyMachine, second: MealyMachine, generator: MealyHMM +) -> None: + staged = transduce_generator(second, transduce_generator(first, generator)) + direct = transduce_generator(compose_transducers((first, second)), generator) + for z in oracles.words(FINAL, 2, min_length=2): + assert close(oracles.word_probability(staged, z), oracles.word_probability(direct, z)) + + +# --------------------------------------------------------------------------- subsequential and weighted + + +@st.composite +def subsequential_transducers(draw: Any) -> SubsequentialTransducer: + n = draw(st.integers(1, 3)) + output_word = st.lists(st.sampled_from(OUT), max_size=2).map(tuple) + machine = SubsequentialTransducer( + input_alphabet=frozenset(AB), + output_alphabet=frozenset(OUT), + initial_states=frozenset({0}), + final_output=draw(st.dictionaries(st.integers(0, n - 1), output_word)), + ) + for state in range(n): + machine.graph.add_state(state) + for state in range(n): + for symbol in AB: + if draw(st.booleans()): + target = draw(st.integers(0, n - 1)) + machine.add_transition(state, target, symbol, draw(st.sampled_from((*OUT, EPSILON)))) + return machine + + +@given(subsequential_transducers(), st.lists(st.sampled_from(AB), max_size=5).map(tuple)) +def test_subsequential_transduce_appends_final_output(machine: SubsequentialTransducer, x: Word) -> None: + machine.validate() + assert machine.is_subsequential() + state: Hashable | None = 0 + out: Word = () + for symbol in x: + step = [t for t in machine.graph.out_transitions(state) if t.data.get(ATTR_SYMBOL) == symbol] + if not step: + state = None + break + emitted = step[0].data.get(ATTR_OUTPUT) + out += () if emitted is EPSILON else (emitted,) + state = step[0].target + expected = set() if state is None else {out + machine.final_output.get(state, ())} + assert machine.transduce(x) == expected + + +@given(stochastic_transducers(), st.integers(0, 4), st.data()) +def test_wfst_weights_match_path_sums(machine: MealyMachine, n: int, data: Any) -> None: + x = data.draw(st.lists(st.sampled_from(AB), min_size=n, max_size=n).map(tuple)) + y = data.draw(st.lists(st.sampled_from(OUT), min_size=n, max_size=n).map(tuple)) + found = path_weights(machine, x, y) + probability = WeightedFiniteStateTransducer.from_transducer(machine) + tropical = WeightedFiniteStateTransducer.from_transducer(machine, semiring="tropical") + assert close(probability.weight(x, y), sum(found)) + expected_cost = -math.log(max(found)) if found else math.inf + assert tropical.weight(x, y) == pytest.approx(expected_cost) + + +# --------------------------------------------------------------------------- unifilar automata + + +@st.composite +def unifilar_automata(draw: Any) -> UnifilarAutomaton: + dfa = draw(dfas()) + result = UnifilarAutomaton( + input_alphabet=frozenset(AB), + initial_states=frozenset(dfa.initial_states), + accepting_states=frozenset(dfa.accepting_states), + ) + for state in dfa.states(): + result.graph.add_state(state) + for t in dfa.transitions(): + if draw(st.integers(0, 3)): + result.graph.add_transition(t.source, t.target, **{ATTR_SYMBOL: t.data.get(ATTR_SYMBOL)}) + return result + + +def delta_set(aut: UnifilarAutomaton, states: frozenset[Hashable], w: Word) -> frozenset[Hashable]: + for symbol in w: + states = oracles._symbol_step(aut, states, symbol) + return states + + +def brute_reset_threshold(aut: UnifilarAutomaton, limit: int = 8) -> float: + full = frozenset(aut.states()) + for length in range(limit + 1): + if any(len(delta_set(aut, full, w)) == 1 for w in oracles.words(AB, length, min_length=length)): + return length + return math.inf + + +def brute_markov_order(aut: UnifilarAutomaton, limit: int = 4) -> float: + full = frozenset(aut.states()) + for length in range(limit + 1): + images = [delta_set(aut, full, w) for w in oracles.words(AB, length, min_length=length)] + if all(len(image) <= 1 for image in images): + return length + return math.inf + + +@given(unifilar_automata()) +def test_unifilar_recognizes_matches_oracle(aut: UnifilarAutomaton) -> None: + aut.validate() + assert aut.is_unifilar() + for w in oracles.words(AB, 5): + assert aut.recognizes(w) == oracles.nfa_accepts(aut, w) + + +@given(unifilar_automata()) +def test_synchronizing_word_is_shortest_reset_word(aut: UnifilarAutomaton) -> None: + threshold = aut.reset_threshold() + assert threshold == brute_reset_threshold(aut) + assert aut.is_exactly_synchronizable() == (threshold != math.inf) + word = aut.synchronizing_word() + if word is None: + assert threshold == math.inf + else: + assert len(word) == threshold + assert len(delta_set(aut, frozenset(aut.states()), tuple(word))) == 1 + + +def futures(aut: UnifilarAutomaton, state: Hashable, n: int = 3) -> frozenset[Word]: + return frozenset(w for w in oracles.words(AB, n) if delta_set(aut, frozenset({state}), w)) + + +@given(unifilar_automata()) +def test_markov_order_matches_brute_force(aut: UnifilarAutomaton) -> None: + # markov_order ignores subsets that can never synchronize, which only exist + # when two states share their futures, so require a minimal presentation. + assume(all(any(True for _ in aut.graph.out_transitions(state)) for state in aut.states())) + assume(len({futures(aut, state) for state in aut.states()}) == len(list(aut.states()))) + order = aut.markov_order() + assert order == brute_markov_order(aut) + assert aut.is_definite() == (order != math.inf) diff --git a/tests/test_properties_automata_visibly.py b/tests/test_properties_automata_visibly.py new file mode 100644 index 0000000..032bf27 --- /dev/null +++ b/tests/test_properties_automata_visibly.py @@ -0,0 +1,328 @@ +"""Property-based tests for visibly pushdown automata, nested word automata, and PAPNI encoding.""" + +from __future__ import annotations + +from typing import Any + +import pytest +from hypothesis import assume, given +from hypothesis import strategies as st + +from sofic.automata import ( + DyckAlphabet, + NestedWordAutomaton, + VisiblyPushdownAutomaton, + encode_dyck_samples, + encode_dyck_word, + is_well_matched, +) +from sofic.automata.vpa import to_multiple_entry, to_single_entry +from sofic.exceptions import NonWellMatchedLanguageError +from sofic.testing import nwas, vpas +from tests import oracles + +VISIBLE = ("c", "r", "i") +N = 5 +WORDS = oracles.words(VISIBLE, N) +EVERYTHING = frozenset(WORDS) + +Word = tuple[Any, ...] + + +def vpa_truth(vpa: VisiblyPushdownAutomaton) -> frozenset[Word]: + return frozenset(w for w in WORDS if oracles.vpa_accepts(vpa, w)) + + +def vpa_lang(vpa: VisiblyPushdownAutomaton) -> frozenset[Word]: + return frozenset(w for w in WORDS if vpa.recognizes(w)) + + +def nwa_truth(nwa: NestedWordAutomaton) -> frozenset[Word]: + return frozenset(w for w in WORDS if oracles.nwa_accepts(nwa, w)) + + +def nwa_lang(nwa: NestedWordAutomaton) -> frozenset[Word]: + return frozenset(w for w in WORDS if nwa.recognizes_visible(w)) + + +def concat_of(left: frozenset[Word], right: frozenset[Word]) -> frozenset[Word]: + return frozenset(w for w in WORDS if any(w[:i] in left and w[i:] in right for i in range(len(w) + 1))) + + +def star_of(language: frozenset[Word]) -> frozenset[Word]: + def member(w: Word) -> bool: + reach = [True] + [False] * len(w) + for end in range(1, len(w) + 1): + reach[end] = any(reach[start] and w[start:end] in language for start in range(end)) + return reach[-1] + + return frozenset(w for w in WORDS if member(w)) + + +def well_matched(w: Word) -> bool: + depth = 0 + for symbol in w: + depth += {"c": 1, "r": -1}.get(symbol, 0) + if depth < 0: + return False + return depth == 0 + + +def well_matched_vpa() -> VisiblyPushdownAutomaton: + """VPA of all well-matched words: state ``top`` means the stack is empty.""" + vpa = VisiblyPushdownAutomaton( + call_alphabet=frozenset({"c"}), + return_alphabet=frozenset({"r"}), + internal_alphabet=frozenset({"i"}), + stack_alphabet=frozenset({"d0", "d1"}), + initial_state="top", + accepting_states=frozenset({"top"}), + ) + for state in ("top", "deep"): + vpa.graph.add_state(state) + vpa.add_internal_transition(state, state, "i") + vpa.add_call_transition("top", "deep", "c", "d0") + vpa.add_call_transition("deep", "deep", "c", "d1") + vpa.add_return_transition("deep", "deep", "r", "d1") + vpa.add_return_transition("deep", "top", "r", "d0") + return vpa + + +# --------------------------------------------------------------------------- VPA + + +def test_well_matched_vpa_fixture() -> None: + assert vpa_truth(well_matched_vpa()) == frozenset(w for w in WORDS if well_matched(w)) + + +@given(vpas()) +def test_vpa_recognizes_matches_oracle(a: VisiblyPushdownAutomaton) -> None: + assert vpa_lang(a) == vpa_truth(a) + + +@given(vpas(), vpas()) +def test_vpa_boolean_operations_match_brute_force(a: VisiblyPushdownAutomaton, b: VisiblyPushdownAutomaton) -> None: + la, lb = vpa_truth(a), vpa_truth(b) + assert vpa_lang(a.union(b)) == la | lb + assert vpa_lang(a.intersection(b)) == la & lb + assert vpa_lang(a.complement()) == EVERYTHING - la + assert vpa_lang(a.difference(b)) == la - lb + + +@given(vpas(max_states=2), vpas(max_states=2)) +def test_vpa_concat_matches_brute_force(a: VisiblyPushdownAutomaton, b: VisiblyPushdownAutomaton) -> None: + assert vpa_lang(a.concat(b)) == concat_of(vpa_truth(a), vpa_truth(b)) + + +@given(vpas(max_states=2)) +def test_vpa_kleene_star_matches_brute_force(a: VisiblyPushdownAutomaton) -> None: + assert vpa_lang(a.kleene_star()) == star_of(vpa_truth(a)) + + +@given(vpas()) +def test_vpa_determinize_preserves_language(a: VisiblyPushdownAutomaton) -> None: + det = a.determinize() + det.validate() + assert vpa_lang(det) == vpa_truth(a) + assert det.equivalent(a) + + +@given(vpas()) +def test_vpa_emptiness_and_witness(a: VisiblyPushdownAutomaton) -> None: + accepted = vpa_truth(a) + witness = a.accepted_word() + assert (witness is None) == a.is_empty() + if witness is not None: + assert oracles.vpa_accepts(a, witness) + else: + assert not accepted + if accepted: + assert not a.is_empty() + + +@given(vpas()) +def test_vpa_universality(a: VisiblyPushdownAutomaton) -> None: + accepted = vpa_truth(a) + if a.is_universal(): + assert accepted == EVERYTHING + if accepted != EVERYTHING: + assert not a.is_universal() + + +@given(vpas(), vpas()) +def test_vpa_inclusion_and_equivalence(a: VisiblyPushdownAutomaton, b: VisiblyPushdownAutomaton) -> None: + la, lb = vpa_truth(a), vpa_truth(b) + if a.includes(b): + assert lb <= la + if not lb <= la: + assert not a.includes(b) + if a.equivalent(b): + assert la == lb + if la != lb: + assert not a.equivalent(b) + assert a.includes(a.intersection(b)) + assert a.union(b).includes(b) + + +@given(vpas()) +def test_vpa_unmatched_word_detection(a: VisiblyPushdownAutomaton) -> None: + accepted = vpa_truth(a) + if any(not well_matched(w) for w in accepted): + assert a.has_unmatched_word() + if not a.has_unmatched_word(): + assert all(well_matched(w) for w in accepted) + + +@given(vpas()) +def test_modular_conversions_on_well_matched_languages(a: VisiblyPushdownAutomaton) -> None: + restricted = a.intersection(well_matched_vpa()) + assert not restricted.has_unmatched_word() + expected = vpa_truth(a) & vpa_truth(well_matched_vpa()) + single = to_single_entry(restricted) + multiple = to_multiple_entry(restricted) + assert vpa_lang(single) == expected + assert vpa_lang(multiple) == expected + + +@given(vpas()) +def test_modular_conversion_rejects_unmatched_languages(a: VisiblyPushdownAutomaton) -> None: + assume(a.has_unmatched_word()) + with pytest.raises(NonWellMatchedLanguageError): + to_single_entry(a) + with pytest.raises(NonWellMatchedLanguageError): + to_multiple_entry(a) + + +# --------------------------------------------------------------------------- NWA + + +@given(nwas()) +def test_nwa_recognizes_matches_oracle(a: NestedWordAutomaton) -> None: + assert nwa_lang(a) == nwa_truth(a) + + +@given(nwas()) +def test_nwa_vpa_round_trips(a: NestedWordAutomaton) -> None: + expected = nwa_truth(a) + assert vpa_lang(a.to_vpa(tag_symbols=False)) == expected + assert nwa_lang(NestedWordAutomaton.from_vpa(a.to_vpa(tag_symbols=False))) == expected + + +@given(vpas()) +def test_nwa_from_vpa_preserves_language(a: VisiblyPushdownAutomaton) -> None: + assert nwa_lang(NestedWordAutomaton.from_vpa(a)) == vpa_truth(a) + + +@given(nwas(), nwas()) +def test_nwa_boolean_operations_match_brute_force(a: NestedWordAutomaton, b: NestedWordAutomaton) -> None: + la, lb = nwa_truth(a), nwa_truth(b) + assert nwa_lang(a.union(b)) == la | lb + assert nwa_lang(a.intersection(b)) == la & lb + assert nwa_lang(a.complement()) == EVERYTHING - la + assert nwa_lang(a.difference(b)) == la - lb + + +@given(nwas(max_states=2), nwas(max_states=2)) +def test_nwa_concat_and_star_match_brute_force(a: NestedWordAutomaton, b: NestedWordAutomaton) -> None: + la, lb = nwa_truth(a), nwa_truth(b) + assert nwa_lang(a.concat(b)) == concat_of(la, lb) + assert nwa_lang(a.kleene_star()) == star_of(la) + + +@given(nwas(), nwas()) +def test_nwa_decision_procedures(a: NestedWordAutomaton, b: NestedWordAutomaton) -> None: + la, lb = nwa_truth(a), nwa_truth(b) + if la: + assert not a.is_empty() + if a.is_empty(): + assert not la + if la != EVERYTHING: + assert not a.is_universal() + if a.is_universal(): + assert la == EVERYTHING + if not lb <= la: + assert not a.includes(b) + if la != lb: + assert not a.equivalent(b) + assert a.equivalent(a.union(a.intersection(b))) + assert a.union(b).includes(a) + + +# --------------------------------------------------------------------------- PAPNI encoding + +DYCK = DyckAlphabet( + call_alphabet=frozenset({"(", "["}), + return_alphabet=frozenset({")", "]"}), + internal_alphabet=frozenset({"i"}), +) +DYCK_SYMBOLS = ("(", "[", ")", "]", "i") + + +@st.composite +def dyck_words(draw: Any) -> Word: + tokens = draw(st.lists(st.sampled_from(("open", "close", "internal")), max_size=8)) + word: list[str] = [] + depth = 0 + for token in tokens: + if token == "open": + word.append(draw(st.sampled_from(("(", "[")))) + depth += 1 + elif token == "close" and depth: + word.append(draw(st.sampled_from((")", "]")))) + depth -= 1 + elif token == "internal": + word.append("i") + word.extend(draw(st.sampled_from((")", "]"))) for _ in range(depth)) + return tuple(word) + + +def counter_well_matched(w: Word) -> bool: + depth = 0 + for symbol in w: + if symbol in DYCK.call_alphabet: + depth += 1 + elif symbol in DYCK.return_alphabet: + depth -= 1 + if depth < 0: + return False + elif symbol not in DYCK.internal_alphabet: + return False + return depth == 0 + + +any_word = st.lists(st.sampled_from(DYCK_SYMBOLS), max_size=6).map(tuple) + + +@given(any_word | dyck_words()) +def test_is_well_matched_matches_counter(w: Word) -> None: + assert is_well_matched(w, DYCK) == counter_well_matched(w) + + +@given(dyck_words()) +def test_encode_dyck_word_pairs_returns_with_matching_calls(w: Word) -> None: + assert is_well_matched(w, DYCK) + encoded = encode_dyck_word(w, DYCK) + assert len(encoded) == len(w) + stack: list[str] = [] + for symbol, code in zip(w, encoded, strict=True): + if symbol in DYCK.return_alphabet: + assert code == (symbol, stack.pop()) + else: + assert code == symbol + if symbol in DYCK.call_alphabet: + stack.append(symbol) + decoded = tuple(code[0] if isinstance(code, tuple) else code for code in encoded) + assert decoded == w + + +@given(st.lists(any_word | dyck_words(), max_size=6)) +def test_encode_dyck_samples_filters_and_encodes(samples: list[Word]) -> None: + kept = [encode_dyck_word(w, DYCK) for w in samples if counter_well_matched(w)] + assert encode_dyck_samples(samples, DYCK) == kept + if all(counter_well_matched(w) for w in samples): + assert encode_dyck_samples(samples, DYCK, drop_non_well_matched=False) == kept + else: + with pytest.raises(ValueError): + encode_dyck_samples(samples, DYCK, drop_non_well_matched=False) + with pytest.raises(ValueError): + encode_dyck_word(next(w for w in samples if not counter_well_matched(w)), DYCK) diff --git a/tests/test_properties_examples.py b/tests/test_properties_examples.py new file mode 100644 index 0000000..f7d779f --- /dev/null +++ b/tests/test_properties_examples.py @@ -0,0 +1,221 @@ +"""Every example factory builds, validates, and matches its documented numbers.""" + +from __future__ import annotations + +import inspect +import math +from collections import defaultdict +from typing import Any + +import numpy as np +import pytest +from hypothesis import given +from hypothesis import strategies as st + +import sofic.examples as examples +from sofic.base import StateMachine +from sofic.examples import processes +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.graph import ATTR_EMISSION, ATTR_EMISSION_DIST, ATTR_OUTPUT, ATTR_PROB, ATTR_SYMBOL, EPSILON + +ALPHABET_FIELDS = ( + "symbol_alphabet", + "observation_alphabet", + "input_alphabet", + "output_alphabet", + "call_alphabet", + "return_alphabet", + "internal_alphabet", +) +MAX_ORACLE_STATES = 200 + +PARAMETRIZED_FACTORIES = { + "afc": lambda: examples.afc(3), + "afc2": lambda: examples.afc2(3), + "bmc_em": lambda: examples.bmc_em(0.3, 0.6), + "bmc_gen": lambda: examples.bmc_gen(0.3, 0.6), + "bmc_lohr": lambda: examples.bmc_lohr(0.3), + "bmc_param": lambda: examples.bmc_param(0.3, 0.6, 0.2, 0.7), + "cyclic_branching": lambda: examples.cyclic_branching(3, 2), + "iid": lambda: examples.iid(3), + "lollipop": lambda: examples.lollipop(2, 3), + "markov_skeleton": lambda: examples.markov_skeleton(2, 2), + "multiple_n": lambda: examples.multiple_n(3), + "period": lambda: examples.period(3), + "periodic": lambda: examples.periodic("0011"), + "rk_gm": lambda: examples.rk_gm(2, 2), + "rn_gm": lambda: examples.rn_gm(2, 2), + "stretched_gm": lambda: examples.stretched_gm(2), + "uniform_mealymc": lambda: examples.uniform_mealymc(1, 2), +} + + +def _required(factory: Any) -> list[str]: + return [ + p.name + for p in inspect.signature(factory).parameters.values() + if p.default is p.empty and p.kind not in (p.VAR_POSITIONAL, p.VAR_KEYWORD) + ] + + +def _model_factories() -> dict[str, Any]: + found: dict[str, Any] = {} + modules = (examples, processes) + candidates = [(name, getattr(module, name)) for module in modules for name in sorted(set(module.__all__))] + for name, factory in candidates: + if not inspect.isfunction(factory) or inspect.signature(factory).return_annotation in ("float", "dict"): + continue + if not _required(factory): + found[name] = factory + found.update(PARAMETRIZED_FACTORIES) + return found + + +FACTORIES = _model_factories() + + +def test_every_required_argument_factory_is_curated(): + names = {name for module in (examples, processes) for name in module.__all__} + needs_args = { + name + for name in names + if inspect.isfunction(getattr(examples, name, None) or getattr(processes, name)) + and _required(getattr(examples, name, None) or getattr(processes, name)) + and inspect.signature(getattr(examples, name, None) or getattr(processes, name)).return_annotation + not in ("float", "dict") + } + missing = needs_args - set(PARAMETRIZED_FACTORIES) + assert not {name for name in missing if not name.startswith(("tent_map", "logic_machine", "uniform_mealyhmm"))} + + +def _emitted_symbols(model: StateMachine) -> set[Any]: + symbols: set[Any] = set() + for transition in model.transitions(): + for key in (ATTR_SYMBOL, ATTR_EMISSION, ATTR_OUTPUT): + if key in transition.data and transition.data[key] is not EPSILON: + symbols.add(transition.data[key]) + for state in model.states(): + symbols.update(model.graph.state_attrs(state).get(ATTR_EMISSION_DIST, {}) or {}) + return symbols + + +@pytest.mark.parametrize("name", sorted(FACTORIES)) +def test_factory_builds_validates_and_uses_string_symbols(name): + model = FACTORIES[name]() + if not isinstance(model, StateMachine): + pytest.skip(f"{name} returns {type(model).__name__}") + model.validate() + assert list(model.states()), name + emitted = _emitted_symbols(model) + assert all(isinstance(symbol, str) for symbol in emitted), sorted(map(repr, emitted)) + for field in ALPHABET_FIELDS: + alphabet = getattr(model, field, None) + if alphabet: + assert all(isinstance(symbol, str) for symbol in alphabet), (field, sorted(map(repr, alphabet))) + + +# --------------------------------------------------------------------------- raw-graph oracle + + +def _entropy(probs: Any) -> float: + return -sum(p * math.log2(p) for p in probs if p > 1e-15) + + +def unifilar_oracle(machine: Any) -> tuple[float, float]: + """``(h_mu, C_mu)`` from the raw edge list: stationary ``pi`` by eigenvector, then ``sum pi H(row)`` and ``H(pi)``.""" + states = list(machine.states()) + index = {state: i for i, state in enumerate(states)} + matrix = np.zeros((len(states), len(states))) + rows: dict[Any, list[float]] = defaultdict(list) + for transition in machine.transitions(): + prob = float(transition.data[ATTR_PROB]) + matrix[index[transition.source], index[transition.target]] += prob + rows[transition.source].append(prob) + eigenvalues, vectors = np.linalg.eig(matrix.T) + pi = np.real(vectors[:, np.argmin(np.abs(eigenvalues - 1.0))]) + pi = pi / pi.sum() + h = sum(pi[index[state]] * _entropy(rows[state]) for state in states) + return h, _entropy(pi) + + +EPSILON_MACHINES = sorted( + name for name in FACTORIES if FACTORIES[name].__annotations__.get("return") in ("EpsilonMachine", EpsilonMachine) +) + + +@pytest.mark.parametrize("name", EPSILON_MACHINES) +def test_epsilon_machine_examples_match_raw_graph_oracle(name): + machine = FACTORIES[name]() + if not isinstance(machine, EpsilonMachine) or len(list(machine.states())) > MAX_ORACLE_STATES: + pytest.skip("not a small epsilon-machine") + if not machine.is_irreducible(): + pytest.skip("stationary distribution is not unique") + h, c = unifilar_oracle(machine) + assert machine.entropy_rate() == pytest.approx(h, abs=1e-8) + assert machine.statistical_complexity() == pytest.approx(c, abs=1e-8) + + +# --------------------------------------------------------------------------- documented numbers + +PHI = (1 + math.sqrt(5)) / 2 +H_THIRD = _entropy([1 / 3, 2 / 3]) + +DOCUMENTED = { + "golden_mean()": (examples.golden_mean, {"entropy_rate": 2 / 3, "statistical_complexity": H_THIRD}), + "golden_mean(1/phi)": ( + lambda: examples.golden_mean(1 / PHI), + {"entropy_rate": math.log2(PHI)}, + ), + "even_process()": (examples.even_process, {"entropy_rate": 2 / 3, "statistical_complexity": H_THIRD}), + "butterfly_process()": ( + examples.butterfly_process, + {"entropy_rate": 1.0, "statistical_complexity": math.log2(5), "crypticity": 0.3}, + ), + "bernoulli()": (examples.bernoulli, {"entropy_rate": 1.0, "statistical_complexity": 0.0}), +} + + +@pytest.mark.parametrize("label", sorted(DOCUMENTED)) +def test_documented_numbers(label): + factory, expected = DOCUMENTED[label] + machine = factory() + for measure, value in expected.items(): + assert getattr(machine, measure)() == pytest.approx(value, abs=1e-9), measure + + +@given(st.floats(0.01, 0.99)) +def test_bernoulli_entropy_rate_is_binary_entropy(p): + machine = examples.bernoulli(p) + assert machine.entropy_rate() == pytest.approx(_entropy([p, 1 - p]), abs=1e-9) + assert machine.statistical_complexity() == pytest.approx(0.0, abs=1e-9) + + +@given(st.floats(0.05, 0.95)) +def test_golden_mean_family_matches_oracle(p): + machine = examples.golden_mean(p) + h, c = unifilar_oracle(machine) + assert machine.entropy_rate() == pytest.approx(h, abs=1e-9) + assert machine.statistical_complexity() == pytest.approx(c, abs=1e-9) + assert machine.entropy_rate() == pytest.approx(_entropy([p, 1 - p]) / (2 - p), abs=1e-9) + + +@given(st.floats(0.05, 0.95), st.floats(0.05, 0.95)) +def test_nemo_process_matches_oracle(p, q): + machine = examples.nemo_process(p, q) + h, c = unifilar_oracle(machine) + assert machine.entropy_rate() == pytest.approx(h, abs=1e-8) + assert machine.statistical_complexity() == pytest.approx(c, abs=1e-8) + + +@pytest.mark.parametrize("k", range(1, 7)) +def test_period_k_is_deterministic_with_log_k_memory(k): + machine = EpsilonMachine.from_hmm(examples.period(k)) + assert len(list(machine.states())) == k + assert machine.entropy_rate() == pytest.approx(0.0, abs=1e-9) + assert machine.statistical_complexity() == pytest.approx(math.log2(k), abs=1e-9) + + +def test_tent_map_forward_matches_documented_information(): + machine = examples.tent_map_misiurewicz_forward() + expected = examples.tent_map_misiurewicz_information_expected() + assert machine.entropy_rate() == pytest.approx(float(expected["entropy_rate"]), abs=1e-8) diff --git a/tests/test_properties_generators.py b/tests/test_properties_generators.py new file mode 100644 index 0000000..74325a2 --- /dev/null +++ b/tests/test_properties_generators.py @@ -0,0 +1,452 @@ +"""Property-based and metamorphic tests for stochastic generators. + +Word distributions, block entropies, and process equivalence are checked +against the brute-force path-enumeration oracles in :mod:`tests.oracles`. +""" + +from __future__ import annotations + +import math +from collections.abc import Callable, Hashable +from itertools import product +from typing import Any + +import numpy as np +import pytest +from hypothesis import assume, given, note, settings +from hypothesis import strategies as st + +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.lumping import is_lumpable, lump +from sofic.generators.mealy import MealyHMM +from sofic.generators.mixed_state_construction import build_mixed_state_presentation +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.testing.strategies import epsilon_machines, mealy_hmms +from tests import oracles + +ALPHABET = ("0", "1") +TOL = 1e-9 + +small_hmms = mealy_hmms(max_states=3) +small_eps = epsilon_machines(max_states=3) +# Bidirectional quantities (E, anatomy, reversal) need a finite reverse machine. +reversible_eps = small_eps.filter(lambda m: m.reverse_is_finite()) +# Bidirectional constructions cost ~50 ms each; scale their example budget with the profile. +heavy = settings(max_examples=max(10, settings.default.max_examples // 4)) + + +# --------------------------------------------------------------------------- helpers + + +def rebuild( + hmm: MealyHMM, + *, + state_map: Callable[[Hashable], Hashable] = lambda s: s, + symbol_map: Callable[[Any], Any] = lambda x: x, + cls: type[MealyHMM] | None = None, + initial: dict[Hashable, float] | None = None, +) -> MealyHMM: + """Copy ``hmm`` with states and symbols renamed through the given maps.""" + cls = type(hmm) if cls is None else cls + source_initial = hmm.initial_distribution if initial is None else initial + out = cls( + observation_alphabet=frozenset(symbol_map(x) for x in hmm.observation_alphabet), + initial_distribution={state_map(s): float(p) for s, p in source_initial.items()}, + ) + for state in hmm.states(): + out.graph.add_state(state_map(state)) + for tr in hmm.transitions(): + out.add_transition( + state_map(tr.source), state_map(tr.target), symbol_map(tr.data[ATTR_EMISSION]), float(tr.data[ATTR_PROB]) + ) + return out + + +def with_stationary_start(hmm: MealyHMM) -> MealyHMM: + pi = hmm.stationary_distribution() + states = hmm.reindex().states + return rebuild(hmm, initial={s: float(p) for s, p in zip(states, pi, strict=True) if p > 0}) + + +def add_unreachable_state(hmm: MealyHMM) -> MealyHMM: + out = rebuild(hmm) + out.graph.add_state("dead") + out.add_transition("dead", "dead", "0", 0.5) + out.add_transition("dead", next(iter(hmm.states())), "1", 0.5) + return out + + +def add_unreachable_copy(hmm: MealyHMM) -> MealyHMM: + out = rebuild(hmm) + for state in hmm.states(): + out.graph.add_state(("copy", state)) + for tr in hmm.transitions(): + out.add_transition(("copy", tr.source), ("copy", tr.target), tr.data[ATTR_EMISSION], float(tr.data[ATTR_PROB])) + return out + + +def split_states(hmm: MealyHMM, weights: list[float]) -> MealyHMM: + """Split every state ``s`` into ``(s, 0)``/``(s, 1)``; ``{(s,0),(s,1)}`` is lumpable.""" + states = list(hmm.states()) + w = {s: weights[i % len(weights)] for i, s in enumerate(states)} + out = MealyHMM( + observation_alphabet=hmm.observation_alphabet, + initial_distribution={ + (s, k): float(p) * (w[s] if k == 0 else 1 - w[s]) + for s, p in hmm.initial_distribution.items() + for k in (0, 1) + }, + ) + for s in states: + for k in (0, 1): + out.graph.add_state((s, k)) + for tr in hmm.transitions(): + p = float(tr.data[ATTR_PROB]) + wt = w[tr.target] + for k in (0, 1): + out.add_transition((tr.source, k), (tr.target, 0), tr.data[ATTR_EMISSION], p * wt) + out.add_transition((tr.source, k), (tr.target, 1), tr.data[ATTR_EMISSION], p * (1 - wt)) + return out + + +def word_dist(hmm: MealyHMM, length: int) -> dict[tuple[Any, ...], float]: + return hmm.word_probabilities(length, sparse=False) + + +def assert_same_dist(left: dict, right: dict, *, atol: float = TOL) -> None: + for word in set(left) | set(right): + assert left.get(word, 0.0) == pytest.approx(right.get(word, 0.0), abs=atol), word + + +def oracle_agrees(a: MealyHMM, b: MealyHMM, *, atol: float = 1e-9) -> bool: + """Brute-force process equality: words up to length ``n_a + n_b - 1`` suffice.""" + depth = len(list(a.states())) + len(list(b.states())) - 1 + for length in range(depth + 1): + da, db = oracles.word_distribution(a, length, ALPHABET), oracles.word_distribution(b, length, ALPHABET) + if any(abs(da[w] - db[w]) > atol for w in da): + return False + return True + + +def perturb(hmm: MealyHMM, data: st.DataObject) -> MealyHMM: + """Move probability mass between two out-edges of a state with at least two.""" + branching = [s for s in hmm.states() if len(list(hmm.graph.out_transitions(s))) >= 2] + assume(branching) + state = data.draw(st.sampled_from(sorted(branching))) + edges = list(hmm.graph.out_transitions(state)) + i, j = data.draw(st.permutations(range(len(edges))))[:2] + eps = float(edges[i].data[ATTR_PROB]) / 2 + out = MealyHMM(observation_alphabet=hmm.observation_alphabet, initial_distribution=dict(hmm.initial_distribution)) + for s in hmm.states(): + out.graph.add_state(s) + for tr in hmm.transitions(): + p = float(tr.data[ATTR_PROB]) + if tr.source == state: + k = edges.index(tr) + p += -eps if k == i else eps if k == j else 0.0 + out.add_transition(tr.source, tr.target, tr.data[ATTR_EMISSION], p) + return out + + +def block_entropies(hmm: MealyHMM, max_length: int) -> list[float]: + return [oracles.block_entropy(oracles.word_distribution(hmm, n, ALPHABET)) for n in range(max_length + 1)] + + +# --------------------------------------------------------------------------- word distributions + + +@given(small_hmms, st.integers(0, 4)) +def test_word_probabilities_sum_to_one_and_match_oracle(hmm, length): + dist = word_dist(hmm, length) + assert sum(dist.values()) == pytest.approx(1.0, abs=TOL) + assert_same_dist(dist, oracles.word_distribution(hmm, length)) + for word, prob in dist.items(): + assert hmm.word_probability(word) == pytest.approx(prob, abs=TOL) + + +@given(small_hmms, st.integers(1, 4)) +def test_words_of_length_is_support_of_word_distribution(hmm, length): + sparse = hmm.words_of_length(length) + dense = oracles.word_distribution(hmm, length) + assert set(sparse) == {w for w, p in dense.items() if p > 0} + assert_same_dist(sparse, dense) + + +@given(small_hmms, st.lists(st.sampled_from(ALPHABET), max_size=6)) +def test_log_word_probability_is_log2_of_word_probability(hmm, word): + prob = hmm.word_probability(word) + logp = hmm.log_word_probability(word) + if prob > 0: + assert logp == pytest.approx(math.log2(prob), abs=1e-9) + else: + assert logp == -math.inf + + +@heavy +@given(small_eps, st.integers(0, 2**32 - 1)) +def test_log_word_probability_does_not_underflow(machine, seed): + word, _ = machine.sample(2000, np.random.default_rng(seed)) + logp = machine.log_word_probability(word) + assert np.isfinite(logp) + assert logp == pytest.approx(machine.log_likelihood(word), rel=1e-9, abs=1e-6) + assert machine.word_probability(word) == pytest.approx(2.0**logp, rel=1e-6, abs=1e-300) + + +def test_tiny_word_probabilities_are_not_truncated(): + """Regression: words with probability below 1e-15 were reported as impossible.""" + coin = MealyHMM(observation_alphabet=frozenset(ALPHABET), initial_distribution={0: 1.0}) + coin.graph.add_state(0) + for symbol in ALPHABET: + coin.add_transition(0, 0, symbol, 0.5) + word = ["0"] * 60 + assert coin.word_probability(word) == 2.0**-60 + assert coin.conditional_word_probability(["1"], word) == 0.5 + + +@given( + small_hmms, + st.lists(st.sampled_from(ALPHABET), max_size=3), + st.lists(st.sampled_from(ALPHABET), max_size=3), +) +def test_conditional_word_probability_is_joint_over_marginal(hmm, word, condition): + marginal = oracles.word_probability(hmm, condition) + assume(marginal > 1e-12) + joint = oracles.word_probability(hmm, [*condition, *word]) + assert hmm.conditional_word_probability(word, condition) == pytest.approx(joint / marginal, abs=1e-9) + + +@given(small_hmms, st.lists(st.sampled_from(ALPHABET), max_size=3)) +def test_conditional_words_of_next_symbol_sum_to_one(hmm, condition): + assume(oracles.word_probability(hmm, condition) > 1e-12) + total = sum(hmm.conditional_word_probability([x], condition) for x in ALPHABET) + assert total == pytest.approx(1.0, abs=1e-9) + + +# --------------------------------------------------------------------------- block entropy / measures + + +@heavy +@given(reversible_eps) +def test_block_entropy_diagram_matches_oracle(machine): + diagram = machine.block_entropy_diagram(4) + assert diagram.block_entropy[0] == 0.0 + np.testing.assert_allclose(diagram.block_entropy, block_entropies(machine, 4), atol=1e-9) + + +@given(small_eps) +def test_entropy_rate_bounds_and_concavity(machine): + h = machine.entropy_rate() + entropies = block_entropies(machine, 5) + assert entropies[0] == pytest.approx(0.0, abs=1e-12) + deltas = np.diff(entropies) + for length in range(1, 6): + assert h <= entropies[length] / length + 1e-9 + assert h <= deltas[length - 1] + 1e-9 + assert np.all(np.diff(deltas) <= 1e-9) + + +@heavy +@given(reversible_eps) +def test_excess_entropy_bounded_by_statistical_complexity(machine): + e, c = machine.excess_entropy(), machine.statistical_complexity() + assert -1e-9 <= e <= c + 1e-9 + assert machine.crypticity() == pytest.approx(c - e, abs=1e-9) + assert machine.crypticity() >= -1e-9 + + +@heavy +@given(reversible_eps) +def test_excess_entropy_lower_bounds_block_entropy_excess(machine): + """``H(L) - L h`` increases to ``E`` from below.""" + h, e = machine.entropy_rate(), machine.excess_entropy() + for length, entropy in enumerate(block_entropies(machine, 5)): + assert entropy - length * h <= e + 1e-9 + + +@heavy +@given(reversible_eps) +def test_information_anatomy_identities(machine): + anatomy = machine.information_anatomy() + h, b, r, rho = (anatomy[k] for k in ("entropy_rate", "bound_mu", "ephemeral_mu", "rho_mu")) + h1 = oracles.block_entropy(oracles.word_distribution(machine, 1, ALPHABET)) + assert h == pytest.approx(machine.entropy_rate(), abs=1e-9) + assert h == pytest.approx(b + r, abs=1e-9) + assert rho == pytest.approx(h1 - h, abs=1e-9) + assert b == pytest.approx(anatomy["bound_structural"] + anatomy["bound_gauge"], abs=1e-9) + assert r == pytest.approx(anatomy["ephemeral_structural"] + anatomy["ephemeral_gauge"], abs=1e-9) + assert min(b, r, rho) >= -1e-9 + assert anatomy["excess_entropy"] == pytest.approx(machine.excess_entropy(), abs=1e-9) + assert anatomy["crypticity"] == pytest.approx(machine.statistical_complexity() - machine.excess_entropy(), abs=1e-9) + + +@given(small_eps) +def test_cryptic_order_at_most_markov_order(machine): + markov, cryptic = machine.markov_order(), machine.cryptic_order() + if math.isfinite(markov): + assert cryptic <= markov + if cryptic == 0: + assert machine.crypticity() == pytest.approx(0.0, abs=1e-9) + + +# --------------------------------------------------------------------------- constructions + + +@given(small_hmms, st.integers(1, 4)) +def test_mixed_state_presentation_reproduces_word_distribution(hmm, length): + try: + msp = build_mixed_state_presentation(hmm, initial_mixed_state=hmm.initial_distribution, max_states=100) + except MixedStateExplosionError: + assume(False) + assert msp.is_unifilar() + assert_same_dist(word_dist(msp, length), oracles.word_distribution(hmm, length), atol=1e-8) + + +@heavy +@given(small_hmms) +def test_epsilon_machine_from_hmm_reproduces_stationary_process(hmm): + assume(hmm.is_irreducible()) + stationary = with_stationary_start(hmm) + if not hmm.is_unifilar(): + try: + build_mixed_state_presentation(stationary, max_states=100) + except MixedStateExplosionError: + assume(False) + machine = EpsilonMachine.from_hmm(stationary) + note(repr(machine)) + for length in range(1, 5): + assert_same_dist(word_dist(machine, length), oracles.word_distribution(stationary, length), atol=1e-7) + + +def test_epsilon_machine_from_hmm_beyond_26_causal_states(): + """Regression: the causal-state quotient crashed on more than 26 states (A-Z labels).""" + rng = np.random.default_rng(0) + histories = list(product(ALPHABET, repeat=5)) + chain = MealyHMM(observation_alphabet=frozenset(ALPHABET), initial_distribution=dict.fromkeys(histories, 1 / 32)) + for history in histories: + chain.graph.add_state(history) + for history in histories: + p = float(rng.uniform(0.2, 0.8)) + for symbol, prob in zip(ALPHABET, (p, 1 - p), strict=True): + chain.add_transition(history, history[1:] + (symbol,), symbol, prob) + machine = EpsilonMachine.from_hmm(chain) + assert len(list(machine.states())) == 32 + assert machine.markov_order() == 5 + assert machine.is_equal_process(chain) + + +@given(small_eps) +def test_epsilon_machine_from_hmm_is_idempotent_on_minimal_machines(machine): + rebuilt = EpsilonMachine.from_hmm(machine) + assert len(list(rebuilt.states())) <= len(list(machine.states())) + assert rebuilt.is_equal_process(machine) + assert rebuilt.statistical_complexity() <= machine.statistical_complexity() + 1e-9 + + +@heavy +@given(reversible_eps) +def test_time_reversal_reverses_words_and_is_an_involution(machine): + reverse = machine.reverse() + for length in range(1, 5): + forward = word_dist(machine, length) + backward = word_dist(reverse, length) + assert_same_dist(backward, {w[::-1]: p for w, p in forward.items()}, atol=1e-9) + assert reverse.reverse().is_equal_process(machine) + assert reverse.entropy_rate() == pytest.approx(machine.entropy_rate(), abs=1e-9) + assert reverse.excess_entropy() == pytest.approx(machine.excess_entropy(), abs=1e-8) + + +@given(small_hmms, st.lists(st.sampled_from([0.25, 0.5, 0.75]), min_size=1, max_size=3)) +def test_lumping_split_states_recovers_process(hmm, weights): + split = split_states(hmm, weights) + partition = [{(s, 0), (s, 1)} for s in hmm.states()] + assert is_lumpable(split, partition) + lumped = lump(split, partition) + assert len(list(lumped.states())) == len(list(hmm.states())) + for length in range(4): + expected = oracles.word_distribution(hmm, length) + assert_same_dist(word_dist(split, length), expected) + assert_same_dist(word_dist(lumped, length), expected) + + +# --------------------------------------------------------------------------- process equivalence + + +@given(small_hmms) +def test_is_equal_process_reflexive(hmm): + assert hmm.is_equal_process(hmm) + assert hmm.is_equal_process(rebuild(hmm)) + + +@given(small_hmms, small_hmms) +def test_is_equal_process_symmetric_and_matches_oracle(a, b): + forward = a.is_equal_process(b) + assert forward == b.is_equal_process(a) + assert forward == oracle_agrees(a, b) + + +@given(small_hmms, st.data()) +def test_is_equal_process_detects_perturbations(hmm, data): + perturbed = perturb(hmm, data) + assert hmm.is_equal_process(perturbed) == oracle_agrees(hmm, perturbed) + + +# --------------------------------------------------------------------------- metamorphic + + +@given(small_hmms, st.data()) +def test_state_relabeling_preserves_process(hmm, data): + states = list(hmm.states()) + names = data.draw(st.permutations([f"s{i}" for i in range(len(states))])) + mapping = dict(zip(states, names, strict=True)) + relabeled = rebuild(hmm, state_map=mapping.__getitem__) + word = data.draw(st.lists(st.sampled_from(ALPHABET), max_size=6)) + for length in range(4): + assert_same_dist(word_dist(relabeled, length), word_dist(hmm, length)) + assert relabeled.log_likelihood(word) == pytest.approx(hmm.log_likelihood(word), abs=1e-9) or ( + hmm.log_likelihood(word) == relabeled.log_likelihood(word) == -math.inf + ) + assert relabeled.is_equal_process(hmm) + + +@heavy +@given(reversible_eps, st.data()) +def test_state_relabeling_preserves_measures(machine, data): + states = list(machine.states()) + names = data.draw(st.permutations(range(len(states)))) + mapping = dict(zip(states, names, strict=True)) + relabeled = rebuild(machine, state_map=mapping.__getitem__) + assert relabeled.entropy_rate() == pytest.approx(machine.entropy_rate(), abs=1e-9) + assert relabeled.statistical_complexity() == pytest.approx(machine.statistical_complexity(), abs=1e-9) + assert relabeled.excess_entropy() == pytest.approx(machine.excess_entropy(), abs=1e-9) + + +@heavy +@given(reversible_eps, st.sampled_from([{"0": "1", "1": "0"}, {"0": "a", "1": "b"}])) +def test_symbol_relabeling_preserves_measures_and_maps_words(machine, mapping): + relabeled = rebuild(machine, symbol_map=mapping.__getitem__) + assert relabeled.entropy_rate() == pytest.approx(machine.entropy_rate(), abs=1e-9) + assert relabeled.statistical_complexity() == pytest.approx(machine.statistical_complexity(), abs=1e-9) + assert relabeled.excess_entropy() == pytest.approx(machine.excess_entropy(), abs=1e-9) + assert relabeled.crypticity() == pytest.approx(machine.crypticity(), abs=1e-9) + for length in range(4): + mapped = {tuple(mapping[x] for x in w): p for w, p in word_dist(machine, length).items()} + assert_same_dist(word_dist(relabeled, length), mapped) + + +@given(small_hmms, st.sampled_from([add_unreachable_state, add_unreachable_copy])) +def test_unreachable_states_leave_process_unchanged(hmm, extend): + extended = extend(hmm) + for length in range(4): + assert_same_dist(word_dist(extended, length), word_dist(hmm, length)) + assert extended.is_equal_process(hmm) + assert hmm.is_equal_process(extended) + + +@heavy +@given(reversible_eps) +def test_epsilon_machine_ignores_unreachable_copy(machine): + extended = rebuild(add_unreachable_copy(machine), cls=MealyHMM) + rebuilt = EpsilonMachine.from_hmm(extended) + assert rebuilt.is_equal_process(machine) + assert rebuilt.statistical_complexity() == pytest.approx(machine.statistical_complexity(), abs=1e-9) diff --git a/tests/test_properties_inference.py b/tests/test_properties_inference.py new file mode 100644 index 0000000..181bde4 --- /dev/null +++ b/tests/test_properties_inference.py @@ -0,0 +1,315 @@ +"""Property-based tests for HMM filtering, decoding, EM, and model selection. + +Every message, posterior, and likelihood is checked against the brute-force +path-enumeration oracles in :mod:`tests.oracles` on random Mealy HMMs. +""" + +from __future__ import annotations + +import math +from collections.abc import Hashable +from typing import Any + +import numpy as np +import pytest +from hypothesis import assume, given, settings +from hypothesis import strategies as st + +from sofic.generators.mealy import MealyHMM +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.inference import hmm as inf +from sofic.inference import learn_epsilon_machine_cssr, learn_epsilon_machine_spectral, learn_epsilon_machine_subtree +from sofic.inference.model_selection import count_free_parameters, score_model +from sofic.testing.strategies import epsilon_machines, mealy_hmms +from tests import oracles + +ALPHABET = ("0", "1") +LN2 = math.log(2.0) + +small_hmms = mealy_hmms(max_states=3) +observations = st.lists(st.sampled_from(ALPHABET), max_size=5) + +# Learners need long samples (~0.5 s each); keep their example budget modest. +learner_settings = settings(max_examples=max(4, settings.default.max_examples // 20)) +markov_eps = epsilon_machines(max_states=3).filter(lambda m: m.markov_order() <= 3) +SWAP = {"0": "1", "1": "0"} +LEARNERS = { + "cssr": (lambda x, alphabet: learn_epsilon_machine_cssr(x, alphabet=alphabet, max_history=4), 1e-12), + "subtree": (lambda x, alphabet: learn_epsilon_machine_subtree(x, alphabet=alphabet, max_history=4), 1e-12), + "spectral": (lambda x, alphabet: learn_epsilon_machine_spectral([x], alphabet=alphabet), 1e-5), +} + + +def order(hmm: MealyHMM) -> list[Hashable]: + return list(hmm.reindex().states) + + +def edge_table(hmm: MealyHMM) -> dict[tuple[Hashable, Any, Hashable], float]: + return {(tr.source, tr.data[ATTR_EMISSION], tr.target): float(tr.data[ATTR_PROB]) for tr in hmm.transitions()} + + +def with_edges(hmm: MealyHMM, edges: dict[tuple[Hashable, Any, Hashable], float]) -> MealyHMM: + out = MealyHMM(observation_alphabet=hmm.observation_alphabet, initial_distribution=dict(hmm.initial_distribution)) + for state in hmm.states(): + out.graph.add_state(state) + for (source, symbol, target), prob in edges.items(): + out.add_transition(source, target, symbol, prob) + return out + + +def ln_likelihood(hmm: MealyHMM, obs: list[Any]) -> float: + return math.log(oracles.word_probability(hmm, obs)) + + +def possible(hmm: MealyHMM, obs: list[Any], *, floor: float = 0.0) -> bool: + return oracles.word_probability(hmm, obs) > floor + + +# --------------------------------------------------------------------------- filtering and smoothing + + +@given(small_hmms, observations) +def test_forward_backward_match_oracle(hmm, obs): + states = order(hmm) + np.testing.assert_allclose(inf.forward(hmm, obs), oracles.forward(hmm, obs, states=states), atol=1e-12) + np.testing.assert_allclose(inf.backward(hmm, obs), oracles.backward(hmm, obs, states=states), atol=1e-12) + + +@given(small_hmms, observations) +def test_forward_backward_product_is_constant_likelihood(hmm, obs): + alpha, beta = inf.forward(hmm, obs), inf.backward(hmm, obs) + np.testing.assert_allclose((alpha * beta).sum(axis=1), oracles.word_probability(hmm, obs), atol=1e-12) + + +@given(small_hmms, observations) +def test_normalized_messages_are_rescaled_raw_messages(hmm, obs): + assume(possible(hmm, obs)) + alpha = inf.forward(hmm, obs) + np.testing.assert_allclose(inf.forward(hmm, obs, normalize=True), alpha / alpha.sum(axis=1, keepdims=True)) + beta = inf.backward(hmm, obs) + np.testing.assert_allclose(inf.backward(hmm, obs, normalize=True), beta / beta.sum(axis=1, keepdims=True)) + + +def test_normalized_backward_final_row_sums_to_one(): + """Regression: ``backward(normalize=True)`` left ``beta[n]`` as all ones.""" + hmm = MealyHMM(observation_alphabet=frozenset(ALPHABET), initial_distribution={0: 1.0}) + for state in (0, 1, 2): + hmm.graph.add_state(state) + hmm.add_transition(state, (state + 1) % 3, "0", 0.5) + hmm.add_transition(state, state, "1", 0.5) + for obs in ([], ["0"], ["0", "1"]): + np.testing.assert_allclose(inf.backward(hmm, obs, normalize=True).sum(axis=1), 1.0) + + +@given(small_hmms, observations) +def test_smooth_matches_oracle_and_rows_sum_to_one(hmm, obs): + gamma = inf.smooth(hmm, obs) + expected = oracles.gamma(hmm, obs, states=order(hmm)) + np.testing.assert_allclose(gamma, expected, atol=1e-10) + if possible(hmm, obs): + np.testing.assert_allclose(gamma.sum(axis=1), 1.0, atol=1e-12) + assert np.all(gamma >= 0.0) + else: + assert not gamma.any() + + +@given(small_hmms, observations) +def test_two_slice_marginals_match_oracle_and_marginalize_to_gamma(hmm, obs): + xi = inf.two_slice_marginals(hmm, obs) + states = order(hmm) + assert xi.shape == (len(obs), len(states), len(states)) + np.testing.assert_allclose(xi, oracles.xi(hmm, obs, states=states), atol=1e-10) + gamma = inf.smooth(hmm, obs) + np.testing.assert_allclose(xi.sum(axis=2), gamma[:-1], atol=1e-10) + np.testing.assert_allclose(xi.sum(axis=1), gamma[1:], atol=1e-10) + + +@given(small_hmms, observations) +def test_log_likelihood_is_log2_oracle_probability(hmm, obs): + prob = oracles.word_probability(hmm, obs) + ll = inf.log_likelihood(hmm, obs) + if prob > 0.0: + assert ll == pytest.approx(math.log2(prob), abs=1e-9) + else: + assert ll == -math.inf + + +@given(small_hmms, st.lists(observations, min_size=1, max_size=3)) +def test_log_likelihood_of_iid_sequences_is_additive(hmm, seqs): + assume(all(possible(hmm, seq) for seq in seqs)) + scores = score_model(hmm, seqs) + assert scores.log_likelihood == pytest.approx(sum(inf.log_likelihood(hmm, seq) for seq in seqs), abs=1e-9) + + +# --------------------------------------------------------------------------- decoding + + +@given(small_hmms, observations) +def test_viterbi_achieves_oracle_maximum(hmm, obs): + path = inf.viterbi(hmm, obs) + best, argmaxes = oracles.viterbi_paths(hmm, obs) + if best == 0.0: + assert path == [] + return + assert len(path) == len(obs) + 1 + assert oracles.path_probability(hmm, path, obs) == pytest.approx(best, rel=1e-9) + assert tuple(path) in argmaxes + + +# --------------------------------------------------------------------------- score and information + + +@given(small_hmms, st.lists(st.sampled_from(ALPHABET), min_size=1, max_size=5)) +def test_score_matches_finite_differences(hmm, obs): + assume(possible(hmm, obs, floor=1e-8)) + edges = edge_table(hmm) + gradient = inf.score(hmm, obs) + assert set(gradient) == set(edges) + step = 1e-6 + for key, prob in edges.items(): + up = with_edges(hmm, {**edges, key: prob + step}) + down = with_edges(hmm, {**edges, key: prob - step}) + numeric = (ln_likelihood(up, obs) - ln_likelihood(down, obs)) / (2 * step) + assert gradient[key] == pytest.approx(numeric, rel=1e-5, abs=1e-5), key + + +@given(small_hmms, st.lists(st.sampled_from(ALPHABET), min_size=1, max_size=4)) +def test_observed_information_matches_finite_difference_hessian(hmm, obs): + assume(possible(hmm, obs, floor=1e-6)) + edges = edge_table(hmm) + labels = inf.free_parameter_labels(hmm) + assume(labels) + references = {} + for source in {label[0] for label in labels}: + out = sorted((k for k in edges if k[0] == source), key=lambda k: (str(k[1]), k[2])) + references[source] = out[-1] + + def shifted(deltas: dict[int, float]) -> float: + table = dict(edges) + for index, delta in deltas.items(): + label = labels[index] + table[label] += delta + table[references[label[0]]] -= delta + return ln_likelihood(with_edges(hmm, table), obs) + + h = 1e-4 + size = len(labels) + numeric = np.zeros((size, size)) + for a in range(size): + for b in range(size): + if a == b: + numeric[a, a] = (shifted({a: h}) - 2 * shifted({}) + shifted({a: -h})) / h**2 + else: + numeric[a, b] = ( + shifted({a: h, b: h}) - shifted({a: h, b: -h}) - shifted({a: -h, b: h}) + shifted({a: -h, b: -h}) + ) / (4 * h**2) + np.testing.assert_allclose(inf.observed_information(hmm, obs), -numeric, rtol=1e-3, atol=1e-3) + + +# --------------------------------------------------------------------------- EM + + +@given(small_hmms, st.lists(st.lists(st.sampled_from(ALPHABET), min_size=1, max_size=8), min_size=1, max_size=3)) +def test_baum_welch_trace_is_non_decreasing(hmm, seqs): + assume(all(possible(hmm, seq) for seq in seqs)) + fitted, trace = inf.baum_welch(hmm, seqs, max_iter=15, tol=0.0) + assert trace + assert trace[0] == pytest.approx(sum(inf.log_likelihood(hmm, seq) for seq in seqs), abs=1e-9) + assert np.all(np.diff(trace) >= -1e-9) + final = sum(inf.log_likelihood(fitted, seq) for seq in seqs) + assert final >= trace[-1] - 1e-9 + assert set(edge_table(fitted)) <= set(edge_table(hmm)) + fitted.validate() + + +# --------------------------------------------------------------------------- model selection + + +@given(small_hmms, st.lists(st.sampled_from(ALPHABET), min_size=1, max_size=6), st.booleans()) +def test_information_criteria_formulas(hmm, obs, include_initial): + scores = score_model(hmm, obs, include_initial=include_initial) + k = count_free_parameters(hmm, include_initial=include_initial) + n = len(obs) + assert (scores.num_parameters, scores.num_observations) == (k, n) + assert k == sum(max(0, len(list(hmm.graph.out_transitions(s))) - 1) for s in hmm.states()) + ( + max(0, sum(1 for p in hmm.initial_distribution.values() if p > 0) - 1) if include_initial else 0 + ) + if not possible(hmm, obs): + assert scores.log_likelihood == -math.inf + assert scores.aic == scores.bic == scores.mdl == math.inf + return + ln_l = LN2 * inf.log_likelihood(hmm, obs) + assert scores.log_likelihood == pytest.approx(ln_l / LN2, abs=1e-9) + assert scores.aic == pytest.approx(2 * k - 2 * ln_l, abs=1e-9) + assert scores.bic == pytest.approx(k * math.log(n) - 2 * ln_l, abs=1e-9) + assert scores.mdl == pytest.approx(0.5 * k * math.log2(n) - ln_l / LN2, abs=1e-9) + if n - k - 1 > 0: + assert scores.aicc == pytest.approx(scores.aic + 2 * k * (k + 1) / (n - k - 1), abs=1e-9) + else: + assert scores.aicc == math.inf + for criterion in ("aic", "aicc", "bic", "mdl"): + assert scores.value(criterion) == getattr(scores, criterion) + + +# --------------------------------------------------------------------------- metamorphic + + +@given(small_hmms, observations, st.data()) +def test_state_relabeling_permutes_posteriors(hmm, obs, data): + states = order(hmm) + names = data.draw(st.permutations([f"q{i}" for i in range(len(states))])) + mapping = dict(zip(states, names, strict=True)) + relabeled = MealyHMM( + observation_alphabet=hmm.observation_alphabet, + initial_distribution={mapping[s]: p for s, p in hmm.initial_distribution.items()}, + ) + for state in states: + relabeled.graph.add_state(mapping[state]) + for (source, symbol, target), prob in edge_table(hmm).items(): + relabeled.add_transition(mapping[source], mapping[target], symbol, prob) + ll, ll2 = inf.log_likelihood(hmm, obs), inf.log_likelihood(relabeled, obs) + assert ll == ll2 == -math.inf or ll == pytest.approx(ll2, abs=1e-12) + position = [order(relabeled).index(mapping[s]) for s in states] + np.testing.assert_allclose(inf.smooth(relabeled, obs)[:, position], inf.smooth(hmm, obs), atol=1e-12) + path = inf.viterbi(relabeled, obs) + best, _ = oracles.viterbi_paths(hmm, obs) + if best > 0.0: + inverse = {v: k for k, v in mapping.items()} + assert oracles.path_probability(hmm, [inverse[s] for s in path], obs) == pytest.approx(best, rel=1e-9) + + +# --------------------------------------------------------------------------- learners + + +def words4(machine: MealyHMM, rename: dict[str, str] | None = None) -> dict[tuple[str, ...], float]: + dist = machine.word_probabilities(4, sparse=False) + return dist if rename is None else {tuple(rename[x] for x in w): p for w, p in dist.items()} + + +def total_variation(left: dict, right: dict) -> float: + return 0.5 * sum(abs(left.get(w, 0.0) - right.get(w, 0.0)) for w in set(left) | set(right)) + + +@pytest.mark.parametrize("name", sorted(LEARNERS)) +@learner_settings +@given(machine=markov_eps) +def test_learners_recover_small_markov_machines(name, machine): + learn, _ = LEARNERS[name] + sample, _ = machine.sample(20_000, np.random.default_rng(0)) + learned = learn(sample, ALPHABET) + assert total_variation(words4(learned), words4(machine)) < 0.05 + assert len(list(learned.states())) == len(list(machine.states())) + + +@pytest.mark.parametrize("name", sorted(LEARNERS)) +@learner_settings +@given(machine=markov_eps) +def test_learners_are_equivariant_under_symbol_relabeling(name, machine): + learn, atol = LEARNERS[name] + sample, _ = machine.sample(20_000, np.random.default_rng(1)) + learned = learn(sample, ALPHABET) + relabeled = learn([SWAP[x] for x in sample], ALPHABET[::-1]) + assert len(list(relabeled.states())) == len(list(learned.states())) + assert total_variation(words4(relabeled, SWAP), words4(learned)) < atol + assert relabeled.entropy_rate() == pytest.approx(learned.entropy_rate(), abs=atol) + assert relabeled.statistical_complexity() == pytest.approx(learned.statistical_complexity(), abs=atol) diff --git a/tests/test_properties_serialization.py b/tests/test_properties_serialization.py new file mode 100644 index 0000000..db5734f --- /dev/null +++ b/tests/test_properties_serialization.py @@ -0,0 +1,252 @@ +"""Property-based round-trip tests for YAML serialization and ``copy()``.""" + +from __future__ import annotations + +import inspect +from collections import Counter +from collections.abc import Hashable, Mapping +from typing import Any + +import networkx as nx +import numpy as np +import pytest +from hypothesis import assume, given +from hypothesis import strategies as st + +import sofic.examples as examples +from sofic.base import StateMachine +from sofic.exceptions import SoficValidationError +from sofic.graph import TransitionGraph +from sofic.serialization import _registry_by_type, model_from_dict, model_from_yaml, model_to_dict, model_to_yaml +from sofic.shifts import ( + LeftFischerCover, + LeftKriegerCover, + RightFischerCover, + RightKriegerCover, + TopologicalMarkovChain, +) +from sofic.testing import ( + buchi_automata, + dfas, + epsilon_machines, + markov_chains, + mealy_hmms, + mealy_transducers, + nfas, + nwas, + sfts, + sofic_shifts, + vpas, + wheeler_nfas, +) + +STATE_FIELDS = frozenset({"initial_states", "accepting_states", "initial_distribution", "initial_state"}) + + +def _tagged(value: Any) -> Any: + """Type-tagged structural key so ``1``, ``1.0``, ``True`` and ``"1"`` stay distinct.""" + if isinstance(value, tuple): + return ("tuple", tuple(_tagged(item) for item in value)) + if isinstance(value, frozenset | set): + return (type(value).__name__, frozenset(_tagged(item) for item in value)) + if isinstance(value, Mapping): + return ("dict", frozenset((_tagged(k), _tagged(v)) for k, v in value.items())) + if isinstance(value, list): + return ("list", tuple(_tagged(item) for item in value)) + if isinstance(value, np.ndarray): + return ("ndarray", str(value.dtype), value.shape, _tagged(value.tolist())) + if isinstance(value, StateMachine): + return ("model", repr(model_to_dict(value))) + return (type(value).__name__, value if isinstance(value, Hashable) else repr(value)) + + +def assert_same_model(restored: StateMachine, original: StateMachine) -> None: + assert type(restored) is type(original) + nodes = lambda m: {_tagged(s): _tagged(dict(a)) for s, a in m.graph.nx.nodes(data=True)} # noqa: E731 + edges = lambda m: Counter( # noqa: E731 + (_tagged(s), _tagged(t), k, _tagged(dict(d))) for s, t, k, d in m.graph.nx.edges(keys=True, data=True) + ) + assert nodes(restored) == nodes(original) + assert edges(restored) == edges(original) + for field in _registry_by_type()[type(original)].fields: + assert _tagged(getattr(restored, field)) == _tagged(getattr(original, field)), field + + +LARGE_MODEL_STATES = 500 + + +def round_trip(model: StateMachine) -> StateMachine: + """YAML round trip; very large models skip the (slow, pure-Python) YAML text layer.""" + if len(model.graph.nx) > LARGE_MODEL_STATES: + restored = model_from_dict(model_to_dict(model)) + else: + restored = model_from_yaml(model_to_yaml(model)) + assert_same_model(restored, model) + assert model_to_dict(restored) == model_to_dict(model) + return restored + + +def relabel_states(model: StateMachine, mapping: dict[Hashable, Hashable]) -> StateMachine: + """Copy of ``model`` with every state renamed, including state-valued metadata.""" + result = model.copy() + result.graph = TransitionGraph(nx.relabel_nodes(model.graph.nx, mapping, copy=True)) + for field in STATE_FIELDS & set(vars(result)): + value = getattr(result, field) + if isinstance(value, frozenset): + setattr(result, field, frozenset(mapping[s] for s in value)) + elif isinstance(value, dict): + setattr(result, field, {mapping[s]: p for s, p in value.items()}) + elif value is not None: + setattr(result, field, mapping[value]) + return result + + +# --------------------------------------------------------------------------- strategy-generated models + + +ADVERSARIAL_ALPHABETS = [("0", "1"), (1, "1"), ("yes", "null"), ("~", "1.0"), ("a b", 'q"')] + +MODEL_STRATEGIES = { + "dfa": st.sampled_from(ADVERSARIAL_ALPHABETS).flatmap(lambda a: dfas(alphabet=a)), + "nfa": st.sampled_from(ADVERSARIAL_ALPHABETS).flatmap(lambda a: nfas(alphabet=a)), + "buchi": buchi_automata(allow_epsilon=True), + "wheeler": wheeler_nfas(), + "markov": markov_chains(), + "mealy_hmm": st.sampled_from(ADVERSARIAL_ALPHABETS).flatmap(lambda a: mealy_hmms(alphabet=a)), + "epsilon_machine": epsilon_machines(), + "sofic": st.sampled_from(ADVERSARIAL_ALPHABETS).flatmap(lambda a: sofic_shifts(alphabet=a)), + "sft": sfts(), + "vpa": vpas(), + "nwa": nwas(), + "transducer": mealy_transducers(), +} + +any_model = st.one_of(*MODEL_STRATEGIES.values()) + +state_labels = st.one_of( + st.integers(-(2**70), 2**70), + st.text(max_size=6), + st.sampled_from(["yes", "no", "null", "~", "1", "1.0", "0x1F", "1_000", "true", "", " ", "- a", "&x", "*y"]), + st.tuples(st.integers(0, 3), st.text(max_size=2)), + st.frozensets(st.integers(0, 3), max_size=3), +) + + +@pytest.mark.parametrize("name", sorted(MODEL_STRATEGIES)) +@given(data=st.data()) +def test_strategy_models_round_trip(name, data): + round_trip(data.draw(MODEL_STRATEGIES[name])) + + +@given(any_model, st.data()) +def test_round_trip_survives_adversarial_state_labels(model, data): + states = list(model.states()) + labels = data.draw(st.lists(state_labels, min_size=len(states), max_size=len(states), unique_by=_tagged)) + assume(len(set(labels)) == len(labels)) + relabeled = relabel_states(model, dict(zip(states, labels, strict=True))) + round_trip(relabeled) + + +@given(sofic_shifts(max_states=3)) +def test_cover_models_round_trip(shift): + round_trip(RightKriegerCover.from_presentation(shift)) + round_trip(LeftKriegerCover.from_presentation(shift)) + try: + round_trip(RightFischerCover.from_presentation(shift)) + round_trip(LeftFischerCover.from_presentation(shift)) + except SoficValidationError: + pass + + +@given(st.integers(1, 3).flatmap(lambda n: st.lists(st.integers(0, 2), min_size=n * n, max_size=n * n))) +def test_topological_markov_chain_round_trips_with_multiplicity(entries): + n = int(len(entries) ** 0.5) + round_trip(TopologicalMarkovChain.from_adjacency(np.array(entries).reshape(n, n))) + + +# --------------------------------------------------------------------------- copy() + + +def _mutate(model: StateMachine) -> None: + model.graph.add_state("__fresh__", note=["mutable"]) + for _source, _target, _key, data in model.graph.nx.edges(keys=True, data=True): + data["mutated"] = True + break + for name, value in vars(model).items(): + if isinstance(value, dict): + value["__fresh__"] = 0.0 + elif isinstance(value, set | list): + getattr(model, name).clear() + + +@given(any_model) +def test_copy_equals_original_and_is_independent(model): + snapshot = model_to_yaml(model) + clone = model.copy() + assert_same_model(clone, model) + _mutate(clone) + assert model_to_yaml(model) == snapshot + assert "__fresh__" not in set(model.states()) + + +# --------------------------------------------------------------------------- example factories + + +def _required(factory: Any) -> list[str]: + return [ + p.name + for p in inspect.signature(factory).parameters.values() + if p.default is p.empty and p.kind not in (p.VAR_POSITIONAL, p.VAR_KEYWORD) + ] + + +PARAMETRIZED_FACTORIES = { + "afc": lambda: examples.afc(3), + "afc2": lambda: examples.afc2(3), + "bmc_em": lambda: examples.bmc_em(0.3, 0.6), + "bmc_gen": lambda: examples.bmc_gen(0.3, 0.6), + "bmc_lohr": lambda: examples.bmc_lohr(0.3), + "cyclic_branching": lambda: examples.cyclic_branching(3, 2), + "iid": lambda: examples.iid(3), + "lollipop": lambda: examples.lollipop(2, 3), + "markov_skeleton": lambda: examples.markov_skeleton(2, 2), + "multiple_n": lambda: examples.multiple_n(3), + "period": lambda: examples.period(3), + "periodic": lambda: examples.periodic("0011"), + "rk_gm": lambda: examples.rk_gm(2, 2), + "rn_gm": lambda: examples.rn_gm(2, 2), + "stretched_gm": lambda: examples.stretched_gm(2), + "uniform_mealymc": lambda: examples.uniform_mealymc(1, 2), +} + + +def example_factories() -> dict[str, Any]: + """Every model factory in :mod:`sofic.examples`: defaults where possible, curated arguments otherwise.""" + found: dict[str, Any] = {} + candidates = [(name, getattr(examples, name)) for name in sorted(set(examples.__all__))] + for list_name in ("process_list", "transducer_list", "epsilon_transducer_list"): + candidates += [(factory.__name__, factory) for factory in getattr(examples, list_name)] + for name, factory in candidates: + if not inspect.isfunction(factory) or inspect.signature(factory).return_annotation in ("float", "dict"): + continue + if not _required(factory): + found[name] = factory + found.update(PARAMETRIZED_FACTORIES) + return found + + +EXAMPLE_FACTORIES = example_factories() + + +@pytest.mark.parametrize("name", sorted(EXAMPLE_FACTORIES)) +def test_example_factory_round_trips_and_copies(name): + model = EXAMPLE_FACTORIES[name]() + if not isinstance(model, StateMachine): + pytest.skip(f"{name} returns {type(model).__name__}, not a model") + assert type(model) in _registry_by_type(), f"{name}: {type(model).__qualname__} is not registered" + round_trip(model) + snapshot = model_to_dict(model) + clone = model.copy() + assert_same_model(clone, model) + _mutate(clone) + assert model_to_dict(model) == snapshot diff --git a/tests/test_properties_shifts.py b/tests/test_properties_shifts.py new file mode 100644 index 0000000..461663a --- /dev/null +++ b/tests/test_properties_shifts.py @@ -0,0 +1,592 @@ +"""Property-based and metamorphic tests for symbolic shifts. + +Every oracle here touches only the raw labeled graph (networkx reachability, +explicit path and word enumeration), never sofic's own trimming, subset +construction, or cover code. +""" + +from __future__ import annotations + +import math +from collections import defaultdict +from collections.abc import Hashable +from itertools import product +from typing import Any + +import networkx as nx +import numpy as np +import pytest +from hypothesis import assume, given, settings +from hypothesis import strategies as st + +from sofic.exceptions import SoficValidationError +from sofic.graph import ATTR_EMISSION, ATTR_PROB, ATTR_SYMBOL +from sofic.shifts import ( + LeftFischerCover, + LeftKriegerCover, + RightFischerCover, + RightKriegerCover, + ShiftOfFiniteType, + SlidingBlockCode, + SoficDyckShift, + SoficShift, + TopologicalMarkovChain, + higher_block_presentation, + right_resolve, +) +from sofic.testing import sfts, sofic_shifts +from tests import oracles + +ALPHABET = ("0", "1") +TOL = 1e-9 +_FEW = settings(max_examples=max(10, settings.default.max_examples // 4)) + + +# --------------------------------------------------------------------------- graph oracles + + +def _edges(shift: Any) -> list[tuple[Hashable, Any, Hashable]]: + return [(t.source, t.data.get(ATTR_SYMBOL), t.target) for t in shift.transitions()] + + +def _digraph(shift: Any) -> nx.DiGraph: + graph = nx.DiGraph() + graph.add_nodes_from(shift.states()) + graph.add_edges_from((source, target) for source, _symbol, target in _edges(shift)) + return graph + + +def _on_cycle(graph: nx.DiGraph) -> set[Hashable]: + cyclic: set[Hashable] = set() + for component in nx.strongly_connected_components(graph): + node = next(iter(component)) + if len(component) > 1 or graph.has_edge(node, node): + cyclic |= component + return cyclic + + +def _bi_infinite_states(shift: Any) -> tuple[set[Hashable], set[Hashable]]: + """``(states with an infinite past, states with an infinite future)``.""" + graph = _digraph(shift) + cyclic = _on_cycle(graph) + past = set(cyclic).union(*(nx.descendants(graph, node) for node in cyclic)) + future = set(cyclic).union(*(nx.ancestors(graph, node) for node in cyclic)) + return past, future + + +def brute_language(shift: Any, n: int) -> frozenset[tuple[Any, ...]]: + """Labels of length-``n`` paths that extend to bi-infinite paths.""" + past, future = _bi_infinite_states(shift) + if not past & future: + return frozenset() + if n == 0: + return frozenset({()}) + out: dict[Hashable, list[tuple[Any, Hashable]]] = defaultdict(list) + for source, symbol, target in _edges(shift): + out[source].append((symbol, target)) + words: set[tuple[Any, ...]] = set() + frontier = [(state, ()) for state in past] + for _ in range(n): + frontier = [(target, word + (symbol,)) for state, word in frontier for symbol, target in out[state]] + words.update(word for state, word in frontier if state in future) + return frozenset(words) + + +def language(shift: Any, n: int) -> frozenset[tuple[Any, ...]]: + words = list(shift.factor_language(n)) + assert len(words) == len(set(words)), "factor_language yielded duplicates" + return frozenset(words) + + +def is_right_resolving(shift: Any) -> bool: + seen = set() + for source, symbol, _target in _edges(shift): + if (source, symbol) in seen: + return False + seen.add((source, symbol)) + return True + + +def is_left_resolving(shift: Any) -> bool: + seen = set() + for _source, symbol, target in _edges(shift): + if (target, symbol) in seen: + return False + seen.add((target, symbol)) + return True + + +def log2_spectral_radius(shift: Any) -> float: + states = list(shift.states()) + if not states: + return 0.0 + index = {state: i for i, state in enumerate(states)} + matrix = np.zeros((len(states), len(states))) + for source, _symbol, target in _edges(shift): + matrix[index[source], index[target]] += 1.0 + radius = float(np.max(np.abs(np.linalg.eigvals(matrix)))) + return math.log2(radius) if radius > TOL else 0.0 + + +def is_irreducible_presentation(shift: Any) -> bool: + """Whether the bi-infinite part of the presentation is one strongly connected component.""" + past, future = _bi_infinite_states(shift) + core = past & future + return bool(core) and nx.is_strongly_connected(_digraph(shift).subgraph(core)) + + +def relabel(shift: SoficShift, state_map: dict, symbol_map: dict) -> SoficShift: + result = SoficShift(symbol_alphabet=frozenset(symbol_map[a] for a in shift.symbol_alphabet)) + for state in shift.states(): + result.graph.add_state(state_map[state]) + for source, symbol, target in _edges(shift): + result.add_transition(state_map[source], state_map[target], symbol_map[symbol]) + return result + + +def follower_signature(shift: Any, state: Hashable, depth: int) -> frozenset[tuple[Any, ...]]: + out: dict[Hashable, list[tuple[Any, Hashable]]] = defaultdict(list) + for source, symbol, target in _edges(shift): + out[source].append((symbol, target)) + words: set[tuple[Any, ...]] = set() + frontier = [(state, ())] + for _ in range(depth): + frontier = [(target, word + (symbol,)) for current, word in frontier for symbol, target in out[current]] + words.update(word for _state, word in frontier) + return frozenset(words) + + +irreducible_shifts = sofic_shifts(max_states=4).filter(is_irreducible_presentation) + + +# --------------------------------------------------------------------------- factor language + + +@given(sofic_shifts(max_states=4), st.integers(0, 5)) +def test_factor_language_matches_bi_extendable_paths(shift, n): + assert language(shift, n) == brute_language(shift, n) + + +@given(sofic_shifts(max_states=4)) +def test_factor_language_is_factorial_and_extendable(shift): + for n in range(1, 5): + longer, shorter = language(shift, n + 1), language(shift, n) + assert {word[1:] for word in longer} == shorter + assert {word[:-1] for word in longer} == shorter + + +# --------------------------------------------------------------------------- SFTs from forbidden words + + +def _contains_forbidden(word: tuple[Any, ...], forbidden: frozenset[tuple[Any, ...]]) -> bool: + return any(word[i : i + len(f)] == f for f in forbidden for i in range(len(word) - len(f) + 1)) + + +def sft_language_oracle(forbidden: frozenset, alphabet: tuple, max_n: int) -> dict[int, frozenset]: + """Words of length ``<= max_n`` that sit inside a long ``F``-free word with ``K`` symbols on each side. + + With memory ``M = max|f| - 1``, ``K = |A|^M + M`` symbols on each side force + a repeated ``M``-block, so the extension pumps to a bi-infinite point. + """ + if not forbidden: + return {n: frozenset(product(alphabet, repeat=n)) for n in range(max_n + 1)} + memory = max(len(word) for word in forbidden) - 1 + pad = len(alphabet) ** memory + memory + allowed = [()] + for _ in range(max_n + 2 * pad): + allowed = [ + word + (a,) + for word in allowed + for a in alphabet + if not _contains_forbidden((word + (a,))[-memory - 1 :], forbidden) + ] + return ( + {n: frozenset(word[pad : pad + n] for word in allowed) for n in range(max_n + 1)} + if allowed + else {n: frozenset() for n in range(max_n + 1)} + ) + + +@given(sfts(max_word_length=3)) +def test_sft_language_equals_forbidden_word_oracle(sft): + forbidden = frozenset(sft.forbidden_words()) + expected = sft_language_oracle(forbidden, ALPHABET, 4) + for n in range(1, 5): + words = language(sft, n) + assert not any(_contains_forbidden(word, forbidden) for word in words) + assert words == expected[n] + + +@given(sfts(max_word_length=3)) +def test_sft_inferred_forbidden_words_regenerate_the_shift(sft): + forbidden = sft.forbidden_words(max_length=3, minimal=True) + rebuilt = ShiftOfFiniteType.from_forbidden_words(set(forbidden), frozenset(ALPHABET)) + for n in range(5): + assert language(rebuilt, n) == language(sft, n) + + +# --------------------------------------------------------------------------- entropy across presentations + + +def _presentations(shift: SoficShift) -> dict[str, SoficShift]: + found = { + "right_resolve": right_resolve(shift), + "right_krieger": RightKriegerCover.from_presentation(shift), + "left_krieger": LeftKriegerCover.from_presentation(shift), + "higher_block_1": higher_block_presentation(shift, 1), + "higher_block_2": higher_block_presentation(shift, 2), + } + try: + found["right_fischer"] = RightFischerCover.from_presentation(shift) + found["left_fischer"] = LeftFischerCover.from_presentation(shift) + except SoficValidationError: + pass + return found + + +@given(sofic_shifts(max_states=3)) +def test_topological_entropy_is_presentation_invariant(shift): + h = shift.topological_entropy() + for name, presentation in _presentations(shift).items(): + assert presentation.topological_entropy() == pytest.approx(h, abs=1e-7), name + for n in range(4): + assert language(presentation, n) == language(shift, n), name + + +@given(sofic_shifts(max_states=3)) +def test_topological_entropy_is_log_spectral_radius_of_right_resolving_cover(shift): + h = shift.topological_entropy() + for cover in (RightKriegerCover.from_presentation(shift), right_resolve(shift)): + assert is_right_resolving(cover) + assert log2_spectral_radius(cover) == pytest.approx(h, abs=1e-7) + + +@given(sofic_shifts(max_states=3)) +def test_topological_entropy_is_below_every_word_count_rate(shift): + """``|L_{m+n}| <= |L_m| |L_n|``, so ``h = inf_n log|L_n| / n``.""" + h = shift.topological_entropy() + for n in range(1, 8): + assert n * h <= math.log2(len(brute_language(shift, n))) + 1e-7 + + +def test_right_fischer_cover_rejects_reducible_shift_with_one_terminal_component(): + shift = SoficShift(symbol_alphabet=frozenset(ALPHABET)) + shift.add_transition("A", "A", "0") + shift.add_transition("A", "B", "0") + shift.add_transition("A", "B", "1") + shift.add_transition("B", "B", "1") + with pytest.raises(SoficValidationError, match="reducible"): + RightFischerCover.from_presentation(shift) + with pytest.raises(SoficValidationError, match="reducible"): + LeftFischerCover.from_presentation(shift) + assert RightKriegerCover.from_presentation(shift).topological_entropy() == pytest.approx(0.0) + + +def _component_shifts(shift: SoficShift) -> list[SoficShift]: + graph = _digraph(shift) + components = [] + for component in nx.strongly_connected_components(graph): + node = next(iter(component)) + if len(component) == 1 and not graph.has_edge(node, node): + continue + sub = SoficShift(symbol_alphabet=shift.symbol_alphabet) + for source, symbol, target in _edges(shift): + if source in component and target in component: + sub.add_transition(source, target, symbol) + components.append(sub) + return components + + +@given(sofic_shifts(max_states=3)) +def test_right_fischer_cover_exists_iff_one_component_presents_the_shift(shift): + """An irreducible shift is presented by one irreducible component of any presentation. + + The left half of a left-transitive point contains every word and eventually + stays inside one component (Lind & Marcus, ch. 3-4). + """ + try: + RightFischerCover.from_presentation(shift) + irreducible = True + except SoficValidationError: + irreducible = False + lengths = range(1, 7) + full = {n: brute_language(shift, n) for n in lengths} + presenting = [all(brute_language(c, n) == full[n] for n in lengths) for c in _component_shifts(shift)] + assert any(presenting) == irreducible + + +# --------------------------------------------------------------------------- Fischer and Krieger covers + + +@given(irreducible_shifts) +def test_right_fischer_cover_is_minimal_irreducible_right_resolving(shift): + cover = RightFischerCover.from_presentation(shift) + assert isinstance(cover, RightFischerCover) + assert is_right_resolving(cover) + assert list(cover.states()) and nx.is_strongly_connected(_digraph(cover)) + for n in range(6): + assert language(cover, n) == language(shift, n) + states = list(cover.states()) + depth = len(states) + signatures = {follower_signature(cover, state, depth) for state in states} + assert len(signatures) == len(states), "Fischer cover states must be follower-separated" + for other in (RightKriegerCover.from_presentation(shift), right_resolve(shift)): + assert len(states) <= len(list(other.trim_transient().states())) + + +@given(irreducible_shifts) +def test_left_fischer_cover_mirrors_right_fischer_cover_of_reverse(shift): + left = LeftFischerCover.from_presentation(shift) + mirrored = RightFischerCover.from_presentation(shift.reverse()) + assert is_left_resolving(left) + assert len(list(left.states())) == len(list(mirrored.states())) + for n in range(6): + assert language(left, n) == language(shift, n) + assert language(left.reverse(), n) == language(mirrored, n) + + +@given(sofic_shifts(max_states=4)) +def test_krieger_covers_resolve_and_preserve_language(shift): + right = RightKriegerCover.from_presentation(shift) + left = LeftKriegerCover.from_presentation(shift) + mirrored = RightKriegerCover.from_presentation(shift.reverse()) + assert is_right_resolving(right) + assert is_left_resolving(left) + assert len(list(left.states())) == len(list(mirrored.states())) + for n in range(6): + assert language(right, n) == language(shift, n) + assert language(left, n) == language(shift, n) + assert language(left.reverse(), n) == language(mirrored, n) + assert language(shift.reverse(), n) == {word[::-1] for word in language(shift, n)} + + +# --------------------------------------------------------------------------- Parry measure + + +def _unifilar_entropy_rate(hmm: Any) -> float: + """``sum_i pi_i H(row_i)`` from raw edges; checks ``pi`` is stationary on the way.""" + pi = {state: float(hmm.initial_distribution.get(state, 0.0)) for state in hmm.states()} + rows: dict[Hashable, list[float]] = defaultdict(list) + flow: dict[Hashable, float] = defaultdict(float) + seen = set() + for transition in hmm.transitions(): + prob = float(transition.data[ATTR_PROB]) + key = (transition.source, transition.data[ATTR_EMISSION]) + assert key not in seen, "Parry measure must be unifilar" + seen.add(key) + rows[transition.source].append(prob) + flow[transition.target] += pi[transition.source] * prob + for state in pi: + assert sum(rows[state]) == pytest.approx(1.0, abs=1e-9) + assert flow[state] == pytest.approx(pi[state], abs=1e-9) + return -sum(pi[s] * sum(p * math.log2(p) for p in rows[s] if p > 0) for s in pi) + + +@st.composite +def irreducible_adjacency(draw: Any) -> np.ndarray: + n = draw(st.integers(1, 4)) + matrix = np.array(draw(st.lists(st.integers(0, 2), min_size=n * n, max_size=n * n)), dtype=float).reshape(n, n) + graph = nx.DiGraph([(i, j) for i in range(n) for j in range(n) if matrix[i, j] > 0]) + graph.add_nodes_from(range(n)) + assume(matrix.any() and nx.is_strongly_connected(graph)) + return matrix + + +@given(irreducible_adjacency()) +def test_parry_measure_entropy_rate_equals_h_top_for_tmc(matrix): + tmc = TopologicalMarkovChain.from_adjacency(matrix) + expected = math.log2(float(np.max(np.abs(np.linalg.eigvals(matrix))))) + assert tmc.topological_entropy() == pytest.approx(expected, abs=1e-9) + assert _unifilar_entropy_rate(tmc.parry_measure()) == pytest.approx(expected, abs=1e-7) + + +@_FEW +@given(sfts(max_word_length=3).filter(is_irreducible_presentation)) +def test_parry_measure_entropy_rate_equals_h_top_for_sft(sft): + shift = SoficShift(graph=sft.graph.copy(), symbol_alphabet=sft.symbol_alphabet) + parry = shift.parry_measure() + parry.validate_stochastic() + assert _unifilar_entropy_rate(parry) == pytest.approx(sft.topological_entropy(), abs=1e-7) + + +def test_parry_measure_of_unifilar_presentation_drops_transient_states(): + """Golden mean core ``A <-> B`` plus a source ``S`` and a dead end ``D``.""" + shift = SoficShift(symbol_alphabet=frozenset(ALPHABET)) + shift.add_transition("A", "A", "0") + shift.add_transition("A", "B", "1") + shift.add_transition("B", "A", "0") + shift.add_transition("B", "D", "1") + shift.add_transition("S", "A", "1") + parry = shift.parry_measure() + parry.validate_stochastic() + assert set(parry.states()) == {"A", "B"} + assert _unifilar_entropy_rate(parry) == pytest.approx(math.log2((1 + math.sqrt(5)) / 2)) + + +@_FEW +@given(irreducible_shifts) +def test_parry_measure_of_sofic_shift_is_maximal_and_supported_on_language(shift): + parry = shift.parry_measure() + assert _unifilar_entropy_rate(parry) == pytest.approx(shift.topological_entropy(), abs=1e-7) + for n in range(1, 4): + support = {word for word, p in oracles.word_distribution(parry, n, ALPHABET).items() if p > 1e-12} + assert support == language(shift, n) + + +@_FEW +@given(irreducible_shifts) +def test_topological_anatomy_splits_h_top(shift): + anatomy = shift.topological_anatomy() + assert anatomy["h_top"] == pytest.approx(shift.topological_entropy(), abs=1e-7) + assert anatomy["h_top"] == pytest.approx(anatomy["b_top"] + anatomy["r_top"], abs=1e-7) + assert anatomy["b_top"] >= -1e-9 + assert anatomy["r_top"] >= -1e-9 + + +# --------------------------------------------------------------------------- sliding block codes + + +@st.composite +def block_codes(draw: Any, alphabet: tuple = ALPHABET, outputs: tuple = ("a", "b", "c")) -> SlidingBlockCode: + memory = draw(st.integers(0, 1)) + anticipation = draw(st.integers(0, 1)) + window = memory + anticipation + 1 + blocks = list(product(alphabet, repeat=window)) + images = draw(st.lists(st.sampled_from(outputs), min_size=len(blocks), max_size=len(blocks))) + return SlidingBlockCode( + dict(zip(blocks, images, strict=True)), + memory=memory, + anticipation=anticipation, + input_alphabet=alphabet, + ) + + +words_over_alphabet = st.lists(st.sampled_from(ALPHABET), max_size=8).map(tuple) + + +@given(block_codes(), words_over_alphabet) +def test_apply_word_is_pointwise_block_map(code, word): + image = code.apply_word(word) + assert len(image) == max(0, len(word) - code.window + 1) + for i, symbol in enumerate(image): + assert symbol == code.block_map[word[i : i + code.window]] + + +@given(block_codes(), block_codes(alphabet=("a", "b", "c"), outputs=ALPHABET), words_over_alphabet) +def test_compose_is_function_composition(first, second, word): + composed = first.compose(second) + assert composed.window == first.window + second.window - 1 + assert composed.apply_word(word) == second.apply_word(first.apply_word(word)) + + +@given(sofic_shifts(max_states=3), block_codes()) +def test_image_language_is_image_of_language(shift, code): + image = code.apply(shift) + for n in range(4): + expected = {code.apply_word(word) for word in language(shift, n + code.window - 1)} + assert language(image, n) == expected + assert image.topological_entropy() <= shift.topological_entropy() + 1e-7 + + +def _higher_block_code(window: int) -> SlidingBlockCode: + return SlidingBlockCode({block: "".join(block) for block in product(ALPHABET, repeat=window)}, memory=window - 1) + + +@given(sofic_shifts(max_states=3), st.integers(1, 3)) +def test_higher_block_conjugacy_preserves_entropy(shift, window): + image = _higher_block_code(window).apply(shift) + assert image.topological_entropy() == pytest.approx(shift.topological_entropy(), abs=1e-7) + for n in range(1, 4): + assert len(language(image, n)) == len(language(shift, n + window - 1)) + + +# --------------------------------------------------------------------------- Dyck shifts + + +@st.composite +def dyck_shifts(draw: Any) -> SoficDyckShift: + n = draw(st.integers(1, 3)) + shift = SoficDyckShift( + call_alphabet=frozenset({"(", "["}), + return_alphabet=frozenset({")", "]"}), + internal_alphabet=frozenset({"i"}), + ) + for state in range(n): + shift.graph.add_state(state) + state = st.integers(0, n - 1) + calls = [ + shift.add_call_transition(s, t, a) + for s, t, a in draw(st.lists(st.tuples(state, state, st.sampled_from("([")), min_size=1, max_size=3)) + ] + returns = [ + shift.add_return_transition(s, t, a) + for s, t, a in draw(st.lists(st.tuples(state, state, st.sampled_from(")]")), min_size=1, max_size=3)) + ] + for s, t in draw(st.lists(st.tuples(state, state), max_size=2)): + shift.add_internal_transition(s, t, "i") + for call, ret in draw(st.lists(st.tuples(st.sampled_from(calls), st.sampled_from(returns)), max_size=4)): + shift.add_matched_pair(call, ret) + shift.validate() + return shift + + +@given(dyck_shifts()) +def test_dyck_reverse_mirrors_language_and_is_an_involution(shift): + reversed_shift = shift.reverse() + reversed_shift.validate() + twice = reversed_shift.reverse() + twice.validate() + assert twice.call_alphabet == shift.call_alphabet + assert twice.return_alphabet == shift.return_alphabet + for n in range(6): + words = language(shift, n) + assert language(reversed_shift, n) == {word[::-1] for word in words} + assert language(twice, n) == words + + +@given(dyck_shifts(), st.lists(st.sampled_from(["(", "[", ")", "]", "i"]), max_size=5).map(tuple)) +def test_dyck_admissibility_agrees_with_factor_language(shift, word): + assert shift.is_admissible_word(word) == (word in language(shift, len(word))) + + +# --------------------------------------------------------------------------- metamorphic + + +_LABELS = st.one_of(st.integers(-5, 50), st.text(max_size=3), st.tuples(st.integers(0, 3), st.text(max_size=1))) + + +@given(sofic_shifts(max_states=4), st.data()) +def test_relabeling_states_and_symbols_preserves_entropy_and_language(shift, data): + states = list(shift.states()) + new_states = data.draw(st.lists(_LABELS, min_size=len(states), max_size=len(states), unique=True)) + new_symbols = data.draw(st.permutations(["x", "y"])) + symbol_map = dict(zip(ALPHABET, new_symbols, strict=True)) + relabeled = relabel(shift, dict(zip(states, new_states, strict=True)), symbol_map) + assert relabeled.topological_entropy() == pytest.approx(shift.topological_entropy(), abs=1e-9) + for n in range(5): + mapped = {tuple(symbol_map[a] for a in word) for word in language(shift, n)} + assert language(relabeled, n) == mapped + + +@given( + sofic_shifts(max_states=3), + st.lists(st.tuples(st.booleans(), st.integers(0, 2), st.sampled_from(ALPHABET)), min_size=1, max_size=4), +) +def test_adding_transient_states_changes_nothing(shift, extra): + """Sources (no in-edges) and sinks (no out-edges) carry no bi-infinite path.""" + states = list(shift.states()) + grown = shift.copy() + for i, (is_source, target, symbol) in enumerate(extra): + name = f"transient{i}" + grown.graph.add_state(name) + anchor = states[target % len(states)] + if is_source: + grown.add_transition(name, anchor, symbol) + else: + grown.add_transition(anchor, name, symbol) + assert grown.topological_entropy() == pytest.approx(shift.topological_entropy(), abs=1e-9) + assert len(list(RightKriegerCover.from_presentation(grown).states())) == len( + list(RightKriegerCover.from_presentation(shift).states()) + ) + for n in range(5): + assert language(grown, n) == language(shift, n) diff --git a/tests/test_properties_viz.py b/tests/test_properties_viz.py new file mode 100644 index 0000000..a2e13ff --- /dev/null +++ b/tests/test_properties_viz.py @@ -0,0 +1,185 @@ +"""Property-based tests for Graphviz / TikZ output on adversarially labeled models.""" + +from __future__ import annotations + +import json +import re +import shutil +from collections.abc import Hashable +from typing import Any + +import networkx as nx +import pytest +from hypothesis import assume, given, settings +from hypothesis import strategies as st + +from sofic.base import StateMachine +from sofic.graph import TransitionGraph +from sofic.testing import dfas, markov_chains, mealy_hmms, mealy_transducers, nfas, sofic_shifts, vpas +from sofic.viz._names import node_names + +# Each example spawns a Graphviz subprocess (~0.1 s), so the DOT test runs a quarter of the profile. +_FEW = settings(max_examples=max(10, settings.default.max_examples // 4)) +STATE_FIELDS = frozenset({"initial_states", "accepting_states", "initial_distribution", "initial_state"}) +NASTY = " \"'\\{}$%&#_^~<>,;=[]()|@*-+.:/éλ→\t\n" + +labels = st.one_of( + st.integers(-20, 20), + st.text(alphabet=NASTY + "ab01", max_size=5), + st.sampled_from(["node", "edge", "graph", "digraph", "subgraph", "strict", "__start__", "start", "1", "-1"]), + st.tuples(st.integers(0, 2), st.text(alphabet="a,{}", max_size=2)), + st.frozensets(st.integers(0, 2), max_size=2), +) +symbol_alphabets = st.sampled_from([("0", "1"), (0, "0"), ("a b", 'q"'), ("\\", "{}"), ("$x$", "%"), ("é", "→")]) + +models = st.one_of( + symbol_alphabets.flatmap(lambda a: sofic_shifts(alphabet=a)), + symbol_alphabets.flatmap(lambda a: dfas(alphabet=a)), + symbol_alphabets.flatmap(lambda a: nfas(alphabet=a)), + symbol_alphabets.flatmap(lambda a: mealy_hmms(alphabet=a)), + markov_chains(), + mealy_transducers(input_alphabet=("^", "_"), output_alphabet=("~", "#")), + vpas(), +) + + +def relabel_states(model: StateMachine, mapping: dict[Hashable, Hashable]) -> StateMachine: + result = model.copy() + result.graph = TransitionGraph(nx.relabel_nodes(model.graph.nx, mapping, copy=True)) + for field in STATE_FIELDS & set(vars(result)): + value = getattr(result, field) + if isinstance(value, frozenset): + setattr(result, field, frozenset(mapping[s] for s in value)) + elif isinstance(value, dict): + setattr(result, field, {mapping[s]: p for s, p in value.items()}) + elif value is not None: + setattr(result, field, mapping[value]) + return result + + +@st.composite +def adversarial_models(draw: Any) -> StateMachine: + model = draw(models) + states = list(model.states()) + new = draw(st.lists(labels, min_size=len(states), max_size=len(states), unique_by=repr)) + assume(len(set(new)) == len(new)) + return relabel_states(model, dict(zip(states, new, strict=True))) + + +# --------------------------------------------------------------------------- node ids + + +@given(st.lists(labels, max_size=8)) +def test_node_names_are_unique_identifiers(states): + names = node_names(states) + assert set(names) == set(states) + assert len(set(names.values())) == len(names) + assert all(re.fullmatch(r"[A-Za-z_][A-Za-z0-9_]*", name) for name in names.values()) + assert "__start__" not in names.values() + + +# --------------------------------------------------------------------------- Graphviz + + +def _balanced_dot(source: str) -> None: + depth, in_string, escaped = 0, False, False + for char in source: + if in_string: + if escaped: + escaped = False + elif char == "\\": + escaped = True + elif char == '"': + in_string = False + elif char == '"': + in_string = True + elif char in "{[": + depth += 1 + elif char in "}]": + depth -= 1 + assert depth >= 0 + assert not in_string and depth == 0 + + +@_FEW +@given(adversarial_models()) +def test_dot_output_parses_with_unique_node_ids(model): + graphviz = pytest.importorskip("graphviz") + source = model.to_graphviz().source + _balanced_dot(source) + if shutil.which("dot") is None: + return + parsed = json.loads(graphviz.Source(source).pipe(format="json")) + nodes = parsed.get("objects", []) + names = [node["name"] for node in nodes] + start = "__start__" in names + assert len(names) == len(set(names)) + assert len(names) == len(list(model.states())) + start + starts = sum(1 for edge in parsed.get("edges", []) if nodes[edge["tail"]]["name"] == "__start__") + assert len(parsed.get("edges", [])) - starts == len(list(model.transitions())) + + +# --------------------------------------------------------------------------- TikZ + + +def _strip_controls(text: str) -> str: + """Drop every ``\\x`` control symbol so escaped braces and dollars don't count.""" + return re.sub(r"\\.", "", text, flags=re.DOTALL) + + +def _balanced_braces(text: str) -> bool: + depth = 0 + for char in _strip_controls(text): + depth += {"{": 1, "}": -1}.get(char, 0) + if depth < 0: + return False + return depth == 0 + + +def _remove_mbox_groups(math: str) -> str: + out, i = [], 0 + while i < len(math): + if math.startswith(r"\mbox{", i): + depth, i = 1, i + len(r"\mbox{") + while depth: + if math[i] == "\\": + i += 2 + continue + depth += {"{": 1, "}": -1}.get(math[i], 0) + i += 1 + continue + out.append(math[i]) + i += 1 + return "".join(out) + + +def _math_segments(text: str) -> list[str]: + positions = [m.start() for m in re.finditer(r"(? Date: Wed, 7 Oct 2026 08:09:45 -0600 Subject: [PATCH 2/3] Make the bidirectional joint-class choice independent of float round-off Two closed joint classes matched the forward/reverse marginals and tied on the anatomy gap up to round-off; exact float comparison picked the wrong one on x86-64 (E = 1.0 instead of the brute-force 1.5), failing CI on Linux and Windows. Co-authored-by: Cursor --- CHANGELOG.md | 5 +++- .../generators/bidirectional_construction.py | 11 +++++--- tests/test_bidirectional_epsilon.py | 25 +++++++++++++++++++ 3 files changed, 36 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9bcecc4..7cf0e08 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -71,7 +71,10 @@ brute-force reference or a closed-form value and has a regression test. but not necessary condition for a finite reverse ε-machine (alternating biased coins has infinitely many reverse beliefs but a two-state reverse machine); `from_time_reversed` and `to_bidirectional` warn when numeric belief - merging returns a finite truncation of an infinite reverse machine. + merging returns a finite truncation of an infinite reverse machine; the bidirectional machine's choice + between joint classes that tie on the anatomy identity no longer depends on + floating-point round-off (it chose a wrong class on x86-64, giving E = 1.0 + instead of 1.5 for a 3-state machine). - **Inference:** Bayesian posterior machines start where the data start (they started in its final state, giving −∞ likelihoods and infinite BIC); WAIC scores from the stationary distribution; `learn_stack_hmm_mle` and the stack diff --git a/sofic/generators/bidirectional_construction.py b/sofic/generators/bidirectional_construction.py index 682e37f..cf5c395 100644 --- a/sofic/generators/bidirectional_construction.py +++ b/sofic/generators/bidirectional_construction.py @@ -389,6 +389,7 @@ def _joint_pi_minimum_support( *, tol: float = 1e-12, marginal_tol: float = 1e-6, + anatomy_tol: float = 1e-9, ) -> dict[tuple[Hashable, Hashable], Any]: """Pick the closed undirected component that is the true bidirectional class. @@ -431,10 +432,12 @@ class is the one whose forward/reverse marginals equal the forward/reverse matching = [candidate for candidate in candidates if candidate[0] == 0.0] else: matching = [candidate for candidate in candidates if candidate[0] <= marginal_tol] - if matching: - best = min(matching, key=lambda candidate: (candidate[1], candidate[2])) - else: - best = min(candidates, key=lambda candidate: (candidate[0], candidate[1], candidate[2])) + pool = matching or [min(candidates, key=lambda candidate: candidate[0])] + # Anatomy gaps that agree to round-off are ties; comparing them exactly made + # the choice depend on the platform's floating-point kernels. + best_gap = min(candidate[1] for candidate in pool) + tied = [candidate for candidate in pool if candidate[1] <= best_gap + anatomy_tol] + best = min(tied, key=lambda candidate: (candidate[2], sorted(map(repr, candidate[3])))) return best[3] diff --git a/tests/test_bidirectional_epsilon.py b/tests/test_bidirectional_epsilon.py index 8b0a098..50a164b 100644 --- a/tests/test_bidirectional_epsilon.py +++ b/tests/test_bidirectional_epsilon.py @@ -504,3 +504,28 @@ def test_ellison_fig15_single_weak_component(): bidir = ellison_fig15_bidirectional() graph = bidir.to_networkx() assert len(list(nx.weakly_connected_components(graph))) == 1 + + +def test_joint_class_choice_is_platform_independent(): + """Regression: two marginal-matching joint classes tied on the anatomy gap up to round-off. + + Exact float comparison picked the 4-pair class on x86-64 (E = 1.0) and the + diagonal class on arm64 (E = 1.5). Brute-force block entropies converge to + E = 1.5, so the crypticity of this machine is 0. + """ + import pytest + + from sofic.generators.epsilon_machine import EpsilonMachine + from sofic.graph import ATTR_EMISSION, ATTR_PROB + + machine = EpsilonMachine( + initial_distribution={0: 0.25, 1: 0.5, 2: 0.25}, observation_alphabet=frozenset({"0", "1"}) + ) + for state in (0, 1, 2): + machine.graph.add_state(state) + for source, symbol, target, prob in [(0, "1", 1, 1.0), (1, "0", 1, 0.5), (1, "1", 2, 0.5), (2, "0", 0, 1.0)]: + machine.graph.add_transition(source, target, **{ATTR_PROB: prob, ATTR_EMISSION: symbol}) + + assert machine.excess_entropy() == pytest.approx(1.5, abs=1e-9) + assert machine.crypticity() == pytest.approx(0.0, abs=1e-9) + assert len(machine.to_bidirectional().joint_distribution()) == 3 From 53325ed7aee184be8aabdbbc524e261852d06163 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Wed, 7 Oct 2026 10:38:27 -0600 Subject: [PATCH 3/3] CI: test with the symbolic extra; keep one job without it The test matrix installed only the test extra, so every sympy-dependent test was skipped. Install symbolic across the matrix, and keep a single Ubuntu job without it so optional-dependency guards stay covered. Co-authored-by: Cursor --- .github/workflows/build.yml | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml index e6f2daa..d6d4f6b 100644 --- a/.github/workflows/build.yml +++ b/.github/workflows/build.yml @@ -30,7 +30,7 @@ jobs: uses: ts-graphviz/setup-graphviz@v2 - name: Install package and dependencies - run: uv sync --extra test + run: uv sync --extra test --extra symbolic - name: Run tests run: uv run pytest --cov-report=xml @@ -43,6 +43,30 @@ jobs: files: coverage.xml fail_ci_if_error: false + test-without-symbolic: + # sympy is optional; this keeps the import guards and pytest.importorskip paths honest. + runs-on: ubuntu-latest + steps: + - name: Checkout repository + uses: actions/checkout@v6 + + - name: Set up uv + uses: astral-sh/setup-uv@v8.1.0 + with: + version: "latest" + + - name: Install Python + run: uv python install 3.13 + + - name: Install Graphviz + uses: ts-graphviz/setup-graphviz@v2 + + - name: Install package and dependencies + run: uv sync --extra test + + - name: Run tests + run: uv run pytest --no-cov + lint: runs-on: ubuntu-latest steps: