From 65a6b8a76894829d2a2f30dc1f6932febf30fba6 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Wed, 7 Oct 2026 07:45:03 -0600 Subject: [PATCH 1/3] Port automata and shift ideas to processes: factor codes, HMM state splitting, minimal quasi-realization, bisimulation, active epsilon-machine learning, support decisions, conditioning, omega-regular probability, Renyi spectrum, experimental canonical residual HMM Co-authored-by: Cursor --- CHANGELOG.md | 49 ++ docs/generators/bisimulation.rst | 69 +++ docs/generators/canonical_residual_hmm.rst | 81 +++ docs/generators/conditioning.rst | 95 ++++ docs/generators/factor_codes.rst | 99 ++++ docs/generators/generators.rst | 9 + docs/generators/minimal_quasi_realization.rst | 69 +++ docs/generators/omega_probability.rst | 84 +++ docs/generators/renyi.rst | 108 ++++ docs/generators/state_splitting.rst | 106 ++++ docs/generators/support.rst | 77 +++ docs/inference/active_epsilon_learning.rst | 123 +++++ docs/inference/inference.rst | 1 + docs/references.bib | 83 +++ sofic/automata/dfa.py | 7 +- sofic/generators/__init__.py | 34 +- sofic/generators/base.py | 45 ++ sofic/generators/canonical_residual.py | 183 +++++++ sofic/generators/conditioning.py | 161 ++++++ sofic/generators/conversions.py | 9 +- sofic/generators/epsilon_construction.py | 39 +- sofic/generators/epsilon_machine.py | 6 + sofic/generators/factor_codes.py | 170 ++++++ sofic/generators/lumping.py | 213 +++++++- sofic/generators/markov.py | 8 +- sofic/generators/mealy.py | 8 +- sofic/generators/measures.py | 23 +- sofic/generators/minimal_quasi_realization.py | 237 +++++++++ sofic/generators/moore.py | 8 +- sofic/generators/omega_probability.py | 156 ++++++ sofic/generators/process_equivalence.py | 2 +- sofic/generators/relative_entropy_rate.py | 48 +- sofic/generators/renyi.py | 246 +++++++++ sofic/generators/state_splitting.py | 337 ++++++++++++ sofic/generators/support.py | 91 ++++ sofic/inference/__init__.py | 8 +- sofic/inference/active.py | 494 ++++++++++++++++++ sofic/shifts/sliding_block_code.py | 19 +- tests/test_active_epsilon_learning.py | 270 ++++++++++ tests/test_bisimulation.py | 212 ++++++++ tests/test_canonical_residual_hmm.py | 170 ++++++ tests/test_conditioning.py | 158 ++++++ tests/test_conversions.py | 20 +- tests/test_dfa.py | 7 + tests/test_epsilon_machine.py | 26 + tests/test_factor_codes.py | 179 +++++++ tests/test_hmm_state_splitting.py | 251 +++++++++ tests/test_mealy_hmm.py | 13 + tests/test_minimal_quasi_realization.py | 186 +++++++ tests/test_omega_probability.py | 172 ++++++ tests/test_process_equivalence.py | 8 + tests/test_renyi.py | 155 ++++++ tests/test_support.py | 125 +++++ 53 files changed, 5481 insertions(+), 76 deletions(-) create mode 100644 docs/generators/bisimulation.rst create mode 100644 docs/generators/canonical_residual_hmm.rst create mode 100644 docs/generators/conditioning.rst create mode 100644 docs/generators/factor_codes.rst create mode 100644 docs/generators/minimal_quasi_realization.rst create mode 100644 docs/generators/omega_probability.rst create mode 100644 docs/generators/renyi.rst create mode 100644 docs/generators/state_splitting.rst create mode 100644 docs/generators/support.rst create mode 100644 docs/inference/active_epsilon_learning.rst create mode 100644 sofic/generators/canonical_residual.py create mode 100644 sofic/generators/conditioning.py create mode 100644 sofic/generators/factor_codes.py create mode 100644 sofic/generators/minimal_quasi_realization.py create mode 100644 sofic/generators/omega_probability.py create mode 100644 sofic/generators/renyi.py create mode 100644 sofic/generators/state_splitting.py create mode 100644 sofic/generators/support.py create mode 100644 sofic/inference/active.py create mode 100644 tests/test_active_epsilon_learning.py create mode 100644 tests/test_bisimulation.py create mode 100644 tests/test_canonical_residual_hmm.py create mode 100644 tests/test_conditioning.py create mode 100644 tests/test_factor_codes.py create mode 100644 tests/test_hmm_state_splitting.py create mode 100644 tests/test_minimal_quasi_realization.py create mode 100644 tests/test_omega_probability.py create mode 100644 tests/test_process_equivalence.py create mode 100644 tests/test_renyi.py create mode 100644 tests/test_support.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 78d5b4b..22a5c6b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -112,6 +112,19 @@ brute-force reference or a closed-form value and has a regression test. paper; it is Chaos 28, 033116 (2018). - `automaton_to_regex` was ambiguous for multi-character symbols (`ab` vs `a`, `b`); such symbols are now quoted (`'ab'`), and `ε`, `∅` and `'` are escaped. +- `HiddenMarkovModel.to_support_dfa` accepted only terminal recurrent subsets, + so it rejected words of positive probability (e.g. `1` for the even process); + every non-empty subset now accepts. +- `entropy_rate` and `state_entropy` no longer raise `TypeError` when state + labels mix types (e.g. `0` and `(1, 0)` after a split). +- `is_equal_process` works on sympy-valued HMMs (object-dtype matrices broke + `matrix_rank`). +- `DFA.add_transition` adds a missing source state, as `NFA.add_transition` + does. +- `EpsilonMachine.from_hmm` split causal states that differed only by + floating-point round-off from belief updates (a fair-coin HMM gave 33 states + and C_mu ≈ 2.46); numeric transition probabilities are now compared within the + mixed-state tolerance. ### New @@ -168,6 +181,42 @@ brute-force reference or a closed-form value and has a regression test. Sierpinski (log₂ 3) examples. - `EpsilonMachine.from_hmm(hmm, max_states=...)` caps the mixed-state enumeration (it silently ignored keyword arguments). +- `minimal_quasi_realization` and `process_rank` (plus + `HiddenMarkovModel.process_rank`): Schützenberger/Fliess minimization of an + HMM's linear representation to a minimal-dimension `QuasiRealization` whose + dimension is the process rank; exact for sympy input. +- `bisimulation_partition` and `coarsest_lumping`: the coarsest strongly + lumpable partition (Larsen–Skou probabilistic bisimulation) by splitter-based + refinement; `lump(model)` and the `lump` methods accept `partition=None` to + mean the coarsest partition. +- `SlidingBlockCode.apply_to_process` and `factor_codes.image_process`: the + image process of a sliding block code; `HiddenMarkovModel.higher_block(k)` for + k-block presentations. Excess entropy and statistical complexity are not + conjugacy invariants (`E(β_k X) = E(X) + (k − 1) h_μ`); only `h_μ` is. +- `sofic.generators.state_splitting`: process-preserving out- and in-splitting + (`split_state`) and `amalgamate`. +- `sofic.generators.support`: `support_nfa`, `support_includes`, + `is_absolutely_continuous`, `support_equal` (stationary support; + finite-cylinder absolute continuity), with matching `HiddenMarkovModel` + methods. +- `condition_on_language`: conditions an HMM on a prefix-closed regular + constraint by the Doob h-transform; the uniform i.i.d. process conditioned on + an irreducible SFT is its Parry measure. +- `omega_probability` for deterministic Büchi properties of HMM output (bottom + strongly connected components of the product chain; Baier & Katoen 2008) and + `regular_language_probability` for exact `P(X_{0:n} ∈ L)`. +- `sofic.inference.learn_epsilon_machine_active` and `ProcessOracle`: an + L*-style active learner for ε-machines from probability and equivalence + queries; exact for finite, exactly synchronizable ε-machines, with `history` + seeding for infinite transients. +- `renyi_entropy_rate`, `pressure` and `rate_function` + (`sofic.generators.renyi`): Rényi entropy rates for α ∈ [0, ∞] (topological, + Shannon and min-entropy rates as special cases) and the large-deviation rate + function of −(1/n) log₂ P(X_{0:n}) by Legendre transform of the pressure. +- Experimental `canonical_residual_hmm` (`sofic.generators.canonical_residual`, + requires `experimental=True`): a generator whose states are the extreme future + morphs of a finite ε-machine, the process analogue of the canonical RFSA; it + can be strictly smaller than the ε-machine. ## 0.3.0 diff --git a/docs/generators/bisimulation.rst b/docs/generators/bisimulation.rst new file mode 100644 index 0000000..d3c9dd0 --- /dev/null +++ b/docs/generators/bisimulation.rst @@ -0,0 +1,69 @@ +.. bisimulation.rst + +************************** +Probabilistic Bisimulation +************************** + +:doc:`Lumping ` aggregates states along a *given* partition. The +coarsest partition that can be lumped is found automatically by computing +**probabilistic bisimulation** :cite:`LarsenSkou1991`: an equivalence relation +on states under which equivalent states send equal probability mass, for every +emitted symbol, into every equivalence class. For a +:class:`~sofic.generators.mealy.MealyHMM` with joint edge law +:math:`P(t, o \mid s)`, states :math:`s \sim s'` are bisimilar when + +.. math:: + + \sum_{t \in B} P(t, o \mid s) = \sum_{t \in B} P(t, o \mid s') + \qquad \text{for every class } B \text{ and every symbol } o. + +This is exactly the strong lumpability condition of Kemeny & Snell +:cite:`KemenySnell1976`, so the bisimulation classes form the coarsest strongly +lumpable partition: every lumpable partition refines it, and the lumped model +generates the same observed process. A :class:`~sofic.generators.moore.MooreHMM` +additionally requires equal state emission laws. A +:class:`~sofic.generators.markov.MarkovChain` has no edge labels, so its coarsest +lumpable partition is a single block unless an ``initial`` partition (for +example, the level sets of an observation function of the state) is supplied +to be refined. + +:func:`~sofic.generators.lumping.bisimulation_partition` computes the classes +by splitter-driven partition refinement in the style of Hopcroft and +Paige--Tarjan: each splitter set divides every block by the per-symbol mass its +states send into the splitter, and once a block that has already served as a +splitter breaks up, all of its pieces but the largest are queued as new +splitters. Probabilities that are all exact (sympy expressions, +:class:`~fractions.Fraction`, or integers) are compared exactly; floating-point +masses are grouped when they lie within ``tol`` of their sorted neighbors. +:func:`~sofic.generators.lumping.coarsest_lumping` returns the quotient model, +and :func:`~sofic.generators.lumping.lump` with ``partition=None`` uses the +coarsest partition. + +For a unifilar presentation whose states are all recurrent, bisimulation +minimization coincides with ε-machine minimization, so the quotient has as many +states as :meth:`EpsilonMachine.from_hmm +`; an ε-machine is +already bisimulation-minimal. For non-unifilar presentations the quotient is +generally not unifilar and can be far smaller than the ε-machine. + +.. ipython:: + + In [1]: from sofic.generators.mealy import MealyHMM + + In [2]: from sofic.generators.lumping import bisimulation_partition, coarsest_lumping + + In [3]: hmm = MealyHMM(initial_distribution={0: 1.0}) + + In [4]: _ = hmm.add_transition(0, 0, "a", 0.5), hmm.add_transition(0, 1, "b", 0.25), hmm.add_transition(0, 2, "b", 0.25) + + In [5]: _ = hmm.add_transition(1, 0, "a", 1.0), hmm.add_transition(2, 0, "a", 1.0) + + In [6]: bisimulation_partition(hmm) + + In [7]: sorted(coarsest_lumping(hmm).states(), key=str) + +API +=== + +.. autofunction:: sofic.generators.lumping.bisimulation_partition +.. autofunction:: sofic.generators.lumping.coarsest_lumping diff --git a/docs/generators/canonical_residual_hmm.rst b/docs/generators/canonical_residual_hmm.rst new file mode 100644 index 0000000..5beafc2 --- /dev/null +++ b/docs/generators/canonical_residual_hmm.rst @@ -0,0 +1,81 @@ +.. canonical_residual_hmm.rst +.. py:module:: sofic.generators.canonical_residual + +********************** +Canonical Residual HMM +********************** + +.. warning:: + + Experimental. :func:`canonical_residual_hmm` raises :class:`RuntimeError` + unless called with ``experimental=True``. The construction has no published + reference. Its correctness is checked against brute-force word + distributions in the test suite. + +The canonical residual finite-state automaton (RFSA) of a regular language has +the *prime* residual languages as its states. A residual is prime when it is not +the union of other residuals :cite:`Denis2002`; see :cite:`MaarandTamm2022` for +its dual, the átomaton. Its stochastic analogue for a process with a finite +ε-machine replaces residual languages by future morphs +:math:`f_s(w) = P(w \mid s)` of the causal states, and unions by nonnegative +combinations. Every morph has :math:`f_s(\lambda) = 1`, so nonnegative +combinations of morphs are convex combinations. + +Construction +============ + +1. Each causal state :math:`s` of positive stationary probability is + represented by its probabilities :math:`P(w_j \mid s)` on a basis of test + words. The test words are found breadth first: keep :math:`w` when + :math:`T^{(w)} \mathbf{1}` is linearly independent of the vectors already + kept. Linear relations among these coordinates are then exactly linear + relations among the morphs. +2. A morph is *extreme* when it is not a nonnegative combination of the other + morphs. This is decided by a linear program + (:func:`scipy.optimize.linprog`). The extreme morphs form the state set + :math:`R`. +3. Every morph is written as :math:`f_t = \sum_{r \in R} c_{t r} f_r` with + :math:`c \ge 0`, again by a linear program. +4. Since :math:`f_r(x w) = \sum_t T^{(x)}_{r t} f_t(w)`, the generator on + :math:`R` with joint transition matrices :math:`(T^{(x)} C)_{R,\cdot}` and + initial law :math:`\pi C` emits the future :math:`f_r` from state + :math:`r`. This follows by induction on word length. The generator + therefore produces the same process. + +The result is a :class:`~sofic.generators.mealy.MealyHMM`. It is in general not +unifilar, and it has at most as many states as the ε-machine. When the morphs +are affinely independent, every morph is extreme and the result is the +ε-machine itself. This holds for every two-state ε-machine, such as the golden +mean and even processes. + +A strictly smaller example +========================== + +Take two hidden states ``a`` and ``b``. Symbols ``2`` and ``3`` lead to ``a`` and +``b`` from either state. Symbol ``0`` is emitted with the same probability from +both states and keeps the state. Symbol ``1`` is emitted only by ``a`` and moves +to ``a`` or ``b`` with probabilities in the ratio ``split : 1 - split``. The +recurrent beliefs are ``a``, ``b`` and the mixture ``(split, 1 - split)``, so +the ε-machine has three causal states. The third morph is a convex combination +of the other two, and the canonical residual HMM has two states. With the +defaults used in the tests (``split = 1/2``) the ε-machine has statistical +complexity 1.55 bits and the canonical residual HMM has state entropy 0.99 +bits. + +Relation to generative complexity +================================= + +Löhr and Ay :cite:`LohrAy2009` distinguish prescient models, whose minimal +member is the ε-machine, from generative HMMs, which can be much smaller. Their +Example 3.6 is a two-state HMM with infinitely many causal states. For that +example :func:`canonical_residual_hmm` raises +:class:`~sofic.exceptions.MixedStateExplosionError`, because it needs a finite +ε-machine. The canonical residual HMM is a generator, so its state count and +its state entropy are upper bounds on the minimal generator size and the +minimal generator state entropy. Nothing else is computed. In particular, no +minimality among generators is claimed. + +API +=== + +.. autofunction:: canonical_residual_hmm diff --git a/docs/generators/conditioning.rst b/docs/generators/conditioning.rst new file mode 100644 index 0000000..fdd2d9b --- /dev/null +++ b/docs/generators/conditioning.rst @@ -0,0 +1,95 @@ +.. conditioning.rst +.. py:module:: sofic.generators.conditioning + +************************************** +Conditioning on a regular constraint +************************************** + +:func:`condition_on_language` conditions a stationary HMM on its output staying +inside the language of a DFA forever: the limit, as :math:`n \to \infty`, of +the process conditioned on every prefix of :math:`X_{0:n}` being accepted. + +Construction +============ + +With symbol matrices :math:`T^{(x)}` and DFA transition function +:math:`\delta`, the product of HMM states with *accepting* DFA states carries +the substochastic matrices + +.. math:: + + M^{(x)}_{(i, q), (j, \delta(q, x))} = T^{(x)}_{ij} + \quad \text{whenever } \delta(q, x) \text{ is accepting.} + +The DFA is completed first, so a missing transition is a move to a rejecting +trap. Only states reached from the stationary-positive HMM states paired with +the DFA's initial state are kept. + +Surviving paths of length :math:`n` have mass of order :math:`\lambda^n`, where +:math:`\lambda` is the largest Perron root among the strongly connected +components of :math:`M = \sum_x M^{(x)}`, and they spend all but a bounded +number of steps in the component attaining it. That component must be unique; +otherwise the limit depends on how the constraint is entered and a +``ValueError`` is raised. On it, with right and left Perron eigenvectors +:math:`h` and :math:`u`, the conditioned process is the Doob +:math:`h`-transform + +.. math:: + + P'(x, j \mid i) = \frac{M^{(x)}_{ij} h_j}{\lambda h_i}, + \qquad \pi'_i \propto u_i h_i, + +started in its stationary law :math:`\pi'`. For a 0-1 matrix this is Parry's +measure of maximal entropy :cite:`Parry1964` (:cite:`LindMarcus1995`, §13.3). +States of the result are pairs ``(hmm state, dfa state)``. + +Which constraints +================= + +The conditioning event is that the DFA run never leaves its accepting states, +i.e. every prefix is accepted. For a prefix-closed language this is just +:math:`X_{0:n} \in L`. For a factorial language, such as the language of a +shift of finite type or sofic shift, the event is also shift invariant, and the +result is the limit of the law of windows :math:`X_{k:k+m}` given +:math:`X_{0:n} \in L` with both :math:`k` and :math:`n - k` large (in the +Cesàro sense when the component is periodic). A language that is not +prefix-closed is effectively replaced by its largest prefix-closed sublanguage. + +Example: the Parry measure +========================== + +Conditioning the fair coin on the golden mean shift (forbid ``11``) gives the +golden mean process at :math:`p = 1/\varphi`, its measure of maximal entropy, +whose entropy rate is the topological entropy :math:`\log_2 \varphi`. The DFA of +the shift's language is obtained by determinizing its presentation with every +state initial and accepting. + +.. ipython:: + + In [1]: from sofic.automata.nfa import NFA + + In [2]: from sofic.examples import bernoulli + + In [3]: from sofic.generators.conditioning import condition_on_language + + In [4]: from sofic.shifts.sft import ShiftOfFiniteType + + In [5]: shift = ShiftOfFiniteType.from_forbidden_words({("1", "1")}, frozenset("01")) + + In [6]: states = frozenset(shift.states()) + + In [7]: dfa = NFA(graph=shift.graph.copy(), input_alphabet=shift.symbol_alphabet, + ...: initial_states=states, accepting_states=states).determinize() + + @doctest float + In [8]: condition_on_language(bernoulli(), dfa).entropy_rate() + Out[8]: 0.6942419136306173 + + @doctest float + In [9]: shift.topological_entropy() + Out[9]: 0.6942419136306162 + +API +=== + +.. autofunction:: condition_on_language diff --git a/docs/generators/factor_codes.rst b/docs/generators/factor_codes.rst new file mode 100644 index 0000000..262e0a0 --- /dev/null +++ b/docs/generators/factor_codes.rst @@ -0,0 +1,99 @@ +.. factor_codes.rst +.. py:module:: sofic.generators.factor_codes + +******************************* +Factor codes on processes +******************************* + +A sliding block code :math:`\Phi` with memory :math:`m`, anticipation :math:`a` +and window :math:`w = m + 1 + a` (:class:`~sofic.shifts.sliding_block_code.SlidingBlockCode`) +maps a process :math:`X` to its *image process* :math:`Y = \Phi(X)` +:cite:`LindMarcus1995` (§1.5): + +.. math:: + + Y_t = \Phi(X_{t:t+w}), + +which is :math:`(\Phi X)_{t+m}` in the two-sided indexing of Lind & Marcus; for a +stationary :math:`X` the shift by :math:`m` is immaterial. The law of +:math:`Y_{0:n}` is the pushforward of the law of :math:`X_{0:n+w-1}`: + +.. math:: + + \Pr(Y_{0:n} = y_{0:n}) = \sum_{x_{0:n+w-1} \,:\, \Phi(x_{0:n+w-1}) = y_{0:n}} \Pr(X_{0:n+w-1} = x_{0:n+w-1}). + +:meth:`SlidingBlockCode.apply_to_process ` +(equivalently :func:`image_process`) builds a Mealy HMM for :math:`Y` by composing +the code's sliding-window transducer with the generator +(:func:`~sofic.automata.transducer_operations.transduce_generator`). States are +pairs :math:`(q, x_{t-w+1:t})` of a generator state and the last :math:`w - 1` +symbols; the initial law is that of :math:`(Q_{w-1}, X_{0:w-1})`, so the composed +transducer starts synchronized with the source and the first output is +:math:`Y_0 = \Phi(X_{0:w})`. A ``ValueError`` is raised if the generator can +emit a :math:`w`-block missing from the code's ``block_map``. + +Invariants +========== + +A stationary coding never increases the entropy rate, :math:`h_\mu(\Phi(X)) \le +h_\mu(X)` :cite:`Gray1990`, so a conjugacy (an invertible code) preserves it; the +topological analogue is :cite:`LindMarcus1995` (§4.1). Neither the excess entropy +nor the statistical complexity is a conjugacy invariant. The cleanest witness is +the higher block code. + +Higher block presentations +========================== + +The higher block code :math:`\beta_k` :cite:`LindMarcus1995` (§1.4) recodes +:math:`X` to overlapping :math:`k`-tuples, + +.. math:: + + Y_t = X_{t:t+k} = (X_t, \ldots, X_{t+k-1}), + +so :math:`\Pr(Y_{0:n} = y_{0:n}) = \Pr(X_{0:n+k-1} = x_{0:n+k-1})` when the tuples +overlap consistently and spell :math:`x_{0:n+k-1}`, and zero otherwise. +:meth:`HiddenMarkovModel.higher_block ` +returns this presentation; it is unifilar whenever the source is. The 1-block +code :math:`y \mapsto y[0]` inverts :math:`\beta_k`, so :math:`\beta_k` is a +conjugacy and :math:`h_\mu(Y) = h_\mu(X)`. However, since +:math:`H[Y_{0:L}] = H[X_{0:L+k-1}]` and +:math:`\mathbf{E} = \lim_L (H[X_{0:L}] - L\, h_\mu)` :cite:`CrutchfieldFeldman2003`, + +.. math:: + + \mathbf{E}(Y) = \mathbf{E}(X) + (k - 1)\, h_\mu(X), + +and, because a :math:`Y`-past fixes both the causal state of :math:`X` and the +last :math:`k - 1` symbols (which the next tuple repeats), + +.. math:: + + C_\mu(Y) = H[\mathcal{S}_0, X_{-(k-1):0}] = C_\mu(X) + H[X_{-(k-1):0} \mid \mathcal{S}_0], + +with :math:`\mathcal{S}_0` the causal state of :math:`X` after :math:`X_{-1}`. +Both identities are checked in ``tests/test_factor_codes.py``. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean + + In [2]: from sofic.shifts.sliding_block_code import SlidingBlockCode + + In [3]: machine = golden_mean(0.5) + + In [4]: pairs = machine.higher_block(2) + + In [5]: sorted(pairs.observation_alphabet) + + In [6]: float(pairs.entropy_rate()), float(machine.entropy_rate()) + + In [7]: first = SlidingBlockCode({(y,): y[0] for y in pairs.observation_alphabet}) + + In [8]: first.apply_to_process(pairs).is_equal_process(machine) + +API +=== + +.. autofunction:: image_process +.. autofunction:: higher_block diff --git a/docs/generators/generators.rst b/docs/generators/generators.rst index 49c24ae..4c7aa7b 100644 --- a/docs/generators/generators.rst +++ b/docs/generators/generators.rst @@ -36,6 +36,11 @@ Computational mechanics edge_machine information_anatomy relative_entropy_rate + support + conditioning + omega_probability + renyi + canonical_residual_hmm correlations predictive_rate_distortion symbolic_hmm @@ -53,6 +58,9 @@ Constructions and conversions constructions conversions lumping + bisimulation + factor_codes + state_splitting Advanced ======== @@ -62,5 +70,6 @@ Advanced nmachine quasi_realization + minimal_quasi_realization synchronization topological_epsilon_enumeration diff --git a/docs/generators/minimal_quasi_realization.rst b/docs/generators/minimal_quasi_realization.rst new file mode 100644 index 0000000..4145bf4 --- /dev/null +++ b/docs/generators/minimal_quasi_realization.rst @@ -0,0 +1,69 @@ +.. minimal_quasi_realization.rst +.. py:module:: sofic.generators.minimal_quasi_realization + +************************* +Minimal Quasi-Realization +************************* + +A hidden Markov model with symbol-labeled matrices :math:`T^{(x)}`, start vector +:math:`\pi`, and final vector :math:`\mathbf{1}` is a *linear representation* of +its word probabilities, + +.. math:: + + P(x_0 x_1 \cdots x_{L-1}) = \pi \, T^{(x_0)} T^{(x_1)} \cdots T^{(x_{L-1})} \, \mathbf{1}. + +Many representations, of different dimensions, define the same process. +Schützenberger's minimization of weighted automata :cite:`Schutzenberger1961` +finds the smallest: only the part of state space that is both *reachable*, +spanned by the row vectors :math:`\pi T^{(w)}`, and *observable*, spanned by the +column vectors :math:`T^{(w)} \mathbf{1}`, matters. Its dimension is the +**process rank**, which equals the rank of the Hankel matrix + +.. math:: + + H_{u, v} = P(uv) + +over all pairs of words :math:`u, v`, and is the dimension of Upper's +generalized-state space :cite:`Upper1997`. It is a lower bound on the number of +states of every HMM presentation of the process, including the ε-machine, and +can be strictly smaller than all of them: the minimal representation is in +general only a quasi-realization with signed entries. + +:func:`minimal_quasi_realization` grows bases :math:`R` (rows) and :math:`O` +(columns) of the reachable and observable spaces breadth-first over words, +factors the finite Hankel core :math:`M = R O = C F` with +:math:`k = \operatorname{rank} M`, and returns the projected representation + +.. math:: + + D^{(x)} = C^{+} R \, T^{(x)} \, O F^{+} + +as a :class:`~sofic.generators.quasi_realization.QuasiRealization` of +dimension :math:`k`. A final change of basis makes the final vector ``tau`` the +all-ones vector, so ``pi`` sums to one. Models whose probabilities are sympy +expressions are reduced exactly and keep object-dtype sympy entries; numeric +models use orthonormal Krylov bases and drop singular values of :math:`M` +below ``tol`` times the largest. :func:`process_rank` (also +:meth:`HiddenMarkovModel.process_rank +`) returns :math:`k`. + +Word indices are 0-based: a length-:math:`L` word is :math:`x_{0:L} = x_0 \cdots x_{L-1}`. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean, iid + + In [2]: from sofic.generators.minimal_quasi_realization import minimal_quasi_realization + + In [3]: iid(3).process_rank() + + In [4]: qr = minimal_quasi_realization(golden_mean(0.5)) + + In [5]: qr.pi.size, qr.pi.sum() + +API +=== + +.. autofunction:: minimal_quasi_realization +.. autofunction:: process_rank diff --git a/docs/generators/omega_probability.rst b/docs/generators/omega_probability.rst new file mode 100644 index 0000000..56f211e --- /dev/null +++ b/docs/generators/omega_probability.rst @@ -0,0 +1,84 @@ +.. omega_probability.rst +.. py:module:: sofic.generators.omega_probability + +***************************************** +Probabilities of regular properties +***************************************** + +A deterministic automaton reading the output of an HMM is in a state that is a +function of the symbols read so far, so the pair (HMM state, automaton state) is +again a finite Markov chain. Both functions below work on this product chain +:cite:`BaierKatoen2008`. The HMM starts in its stationary law unless a +``start`` (a state, a state-to-mass mapping, or a vector) is given. + +Finite words +============ + +:func:`regular_language_probability` returns +:math:`P(X_{0:n} \in L)` for the language :math:`L` of a DFA exactly, by +pushing one forward vector per DFA state through :math:`T^{(x)}` along +:math:`\delta(q, x)` for :math:`n` steps. For the fair coin and the golden mean +language, half of the 16 words of length 4 avoid ``11``: + +.. ipython:: + + In [1]: from sofic.automata.nfa import NFA + + In [2]: from sofic.examples import bernoulli, golden_mean + + In [3]: from sofic.generators.omega_probability import omega_probability, regular_language_probability + + In [4]: from sofic.shifts.sft import ShiftOfFiniteType + + In [5]: shift = ShiftOfFiniteType.from_forbidden_words({("1", "1")}, frozenset("01")) + + In [6]: states = frozenset(shift.states()) + + In [7]: dfa = NFA(graph=shift.graph.copy(), input_alphabet=shift.symbol_alphabet, + ...: initial_states=states, accepting_states=states).determinize() + + @doctest float + In [8]: regular_language_probability(bernoulli(), dfa, 4) + Out[8]: 0.5 + +Infinite words +============== + +:func:`omega_probability` returns the probability that the infinite output +:math:`X_0 X_1 \cdots` is accepted by a *deterministic* Büchi automaton (DBA). +A finite Markov chain almost surely enters a bottom strongly connected component +of the product and then visits each of its states infinitely often, so the run +is accepted iff that component contains an accepting automaton state. The +answer is the mass absorbed into such components, a linear solve on the +transient states :cite:`BaierKatoen2008`. A missing DBA transition moves to an +absorbing rejecting state. + +DBAs recognize a strict subclass of the ω-regular languages. Safety properties +("never ``11``"), reachability ("eventually ``1``"), and recurrence +("infinitely many ``1``") are DBA-expressible; persistence ("eventually always +``0``", i.e. finitely many ``1``) is not, but its probability is one minus that +of its complement "infinitely many ``1``". A non-deterministic Büchi automaton +raises :class:`~sofic.exceptions.NonDeterministicError`. + +.. ipython:: + + In [9]: from sofic.automata.buchi import BuchiAutomaton + + In [10]: ones = BuchiAutomaton(input_alphabet=frozenset("01"), initial_states=frozenset({"a"}), + ....: accepting_states=frozenset({"b"})) + + In [11]: ones.graph.add_state("a"); ones.graph.add_state("b") + + In [12]: for source, target, symbol in [("a", "a", "0"), ("a", "b", "1"), ("b", "a", "0"), ("b", "b", "1")]: + ....: ones.add_transition(source, target, symbol) + + @doctest float + In [13]: omega_probability(golden_mean(0.4), ones) + Out[13]: 1.0 + +API +=== + +.. autofunction:: regular_language_probability + +.. autofunction:: omega_probability diff --git a/docs/generators/renyi.rst b/docs/generators/renyi.rst new file mode 100644 index 0000000..307c4df --- /dev/null +++ b/docs/generators/renyi.rst @@ -0,0 +1,108 @@ +.. renyi.rst +.. py:module:: sofic.generators.renyi + +********************************************** +Rényi Entropy Rates and Thermodynamic Spectrum +********************************************** + +The Rényi entropy rate of order :math:`\alpha \in [0, \infty]` of a stationary +process, in bits per symbol, is + +.. math:: + + h_\alpha = \lim_{n \to \infty} \tfrac{1}{n} H_\alpha[X_{0:n}], + \qquad H_\alpha[X_{0:n}] = \tfrac{1}{1 - \alpha} \log_2 \sum_w P(w)^\alpha, + +where :math:`w` ranges over the words of length :math:`n`. The pressure is +:math:`\mathcal{P}(\beta) = \lim_n \tfrac{1}{n} \log_2 \sum_w P(w)^\beta`, so +:math:`h_\alpha = \mathcal{P}(\alpha) / (1 - \alpha)`. + +Transfer matrix formula +======================= + +On a unifilar presentation with labeled matrices :math:`T^{(x)}`, restricted to +the states of positive stationary probability, + +.. math:: + + \mathcal{P}(\beta) = \log_2 \rho(A_\beta), \qquad + A_\beta = \sum_x \big(T^{(x)}\big)^{\circ \beta}, \qquad \beta \ge 0, + +where :math:`\circ \beta` raises the positive entries to the power +:math:`\beta` and :math:`\rho` is the spectral radius. For Markov chains this is +:cite:`Rached2001`. For a unifilar presentation, each start state :math:`s` +emits a word :math:`w` along at most one path, with probability +:math:`p(w \mid s)`, and :math:`P(w) = \sum_s \pi_s p(w \mid s)` over the +:math:`N` support states. For :math:`\beta \ge 0`, + +.. math:: + + N^{-1} \sum_s \big(\pi_s p(w \mid s)\big)^\beta \le P(w)^\beta + \le N^\beta \sum_s \big(\pi_s p(w \mid s)\big)^\beta . + +Summing over :math:`w` gives :math:`(\pi^{\circ \beta})^\top A_\beta^n \mathbf{1}` +up to these constant factors. Both outer vectors are positive on the support, so +the growth rate is :math:`\rho(A_\beta)`. Negative :math:`\beta` is not +supported, because then each word is weighted by its *least* likely start +state. + +Special orders: + +* :math:`h_0` is the topological entropy of the process support: :math:`A_0` is + the support adjacency matrix. +* :math:`h_1 = h_\mu`, the Shannon entropy rate, is the limit of the formula at + :math:`\alpha \to 1`. It is computed by + :meth:`~sofic.generators.base.HiddenMarkovModel.entropy_rate`. +* :math:`h_\infty = -\lim_n \tfrac{1}{n} \log_2 \max_w P(w)` is the min-entropy + rate. It equals minus the maximum mean :math:`\log_2` edge probability over + cycles of the support graph, computed with Karp's algorithm. + +:math:`h_\alpha` is non-increasing in :math:`\alpha`. + +Non-unifilar generators are first converted with +:meth:`EpsilonMachine.from_hmm(model, max_states=...) +`. If the mixed-state +presentation does not close within ``max_states`` this raises +:class:`~sofic.exceptions.MixedStateExplosionError`. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean + + In [2]: from sofic.generators.renyi import renyi_entropy_rate + + @doctest float + In [3]: renyi_entropy_rate(golden_mean(0.5), 0.0) + Out[3]: 0.6942419136306174 + + @doctest float + In [4]: renyi_entropy_rate(golden_mean(0.5), float("inf")) + Out[4]: 0.5 + +Large deviations +================ + +Let :math:`U_n = -\tfrac{1}{n} \log_2 P(X_{0:n})`. Its scaled cumulant generating +function in base 2 is :math:`\lim_n \tfrac{1}{n} \log_2 \mathbb{E}[2^{n t U_n}] = \mathcal{P}(1 - t)`. +By the Gärtner-Ellis theorem :cite:`DemboZeitouni1998`, the rate function of +:math:`U_n` is the Legendre transform + +.. math:: I(u) = \sup_\beta \big[(1 - \beta) u - \mathcal{P}(\beta)\big]. + +:func:`rate_function` takes the supremum over :math:`\beta \in [0, \beta_{\max}]` +on a grid, refined with :func:`scipy.optimize.minimize_scalar`. That gives +:math:`I` on :math:`[h_\infty, u_0]`, where :math:`u_0 = -\mathcal{P}'(0)`. +:math:`I` is non-negative and vanishes at :math:`u = h_\mu`. It is ``inf`` below +:math:`h_\infty`, and ``nan`` above :math:`u_0`, where :math:`\beta < 0` would +be needed. Near :math:`h_\infty` the optimal :math:`\beta` diverges, so values +there are lower bounds that tighten as ``beta_max`` grows. For an i.i.d. source +the result is Cramér's rate function. + +API +=== + +.. autofunction:: renyi_entropy_rate + +.. autofunction:: pressure + +.. autofunction:: rate_function diff --git a/docs/generators/state_splitting.rst b/docs/generators/state_splitting.rst new file mode 100644 index 0000000..a1ff41d --- /dev/null +++ b/docs/generators/state_splitting.rst @@ -0,0 +1,106 @@ +.. state_splitting.rst +.. py:module:: sofic.generators.state_splitting + +******************************************* +State splitting and amalgamation of HMMs +******************************************* + +State splitting :cite:`LindMarcus1995` (§2.4; for edge shifts see +:mod:`sofic.shifts.state_splitting`) replaces a state by several copies, +partitioning either its out-edges or its in-edges. For a Mealy HMM the edge +probabilities must be redistributed so that the generated *process* is +unchanged; such splittings and their inverses are among the equivalence moves +between HMM presentations studied by :cite:`Upper1997`. Edges are named +``(source, target, key)`` (:func:`out_edges`, :func:`in_edges`), and splitting +state :math:`s` into :math:`m` copies names them ``(s, 0), ..., (s, m - 1)``. + +Write :math:`p(e)` for the joint probability :math:`\Pr(\text{target}, \text{symbol} +\mid \text{source})` of edge :math:`e` and :math:`\pi` for the initial distribution. + +Out-splitting +============= + +Partition the out-edges of :math:`s` into :math:`E_0, \ldots, E_{m-1}` with masses +:math:`w_i = \sum_{e \in E_i} p(e)`. In the split model + +* copy :math:`(s, i)` keeps the edges of :math:`E_i` with probabilities + :math:`p(e) / w_i`; +* every in-edge :math:`e` of :math:`s` becomes :math:`m` edges, one into each copy, + with probabilities :math:`p(e)\, w_i`; +* :math:`\pi(s, i) = \pi(s)\, w_i`. + +By induction on :math:`t`, the forward vectors satisfy +:math:`\alpha_t(s, i) = w_i\, \alpha_t(s)`: inflow into a copy is the inflow into +:math:`s` scaled by :math:`w_i`, and outflow from :math:`(s, i)` is +:math:`w_i \alpha_t(s) \cdot p(e) / w_i = \alpha_t(s)\, p(e)` along each +:math:`e \in E_i`. Summing over copies, every word :math:`x_{0:n}` has the same +probability. If :math:`\pi` is stationary, so is the split :math:`\pi`. +Out-splitting generally breaks unifilarity, since an in-edge becomes :math:`m` +parallel edges with the same symbol. + +In-splitting +============ + +Partition the in-edges of :math:`s` into :math:`F_0, \ldots, F_{m-1}`. Edge +:math:`e \in F_i` enters only copy :math:`(s, i)`, and every copy keeps all +out-edges of :math:`s` with their original probabilities. The copies have +identical futures, so the backward vectors agree, +:math:`\beta_t(s, i) = \beta_t(s)`, and :math:`\pi(s)` may be shared among the +copies in any proportion. It is shared in proportion to the stationary inflow +:math:`\sum_{e \in F_i} \mu(\operatorname{source}(e))\, p(e)` (uniformly when +:math:`s` has none), which keeps a stationary :math:`\pi` stationary. In-splitting +preserves unifilarity. + +Amalgamation +============ + +:func:`amalgamate` merges the blocks of a state partition when one of two exact +rules applies: + +* **identical futures** -- the partition is strongly lumpable and is delegated to + :func:`~sofic.generators.lumping.lump` :cite:`KemenySnell1976`; this inverts + in-splitting; +* **proportional inflow** -- for each block :math:`B` there is a weight vector + :math:`w` on :math:`B` such that, for every source state and symbol, the edge + mass into :math:`B` and the initial mass on :math:`B` are multiples of + :math:`w`. The forward vector restricted to :math:`B` is then always a multiple + of :math:`w`, so :math:`B` merges into a single state whose out-edges are the + :math:`w`-mixture of its members' out-edges; this inverts out-splitting. + +The whole partition is tried under each rule, then blocks are merged one at a +time, so a model that was both in- and out-split amalgamates in one call. A +block with neither inflow nor initial mass is never visited and is merged with +uniform weights: the process is unchanged, but its original out-split weights +cannot be recovered. + +Splits and amalgamations leave every process quantity unchanged -- word +distributions, :math:`h_\mu`, :math:`\mathbf{E}`, and the ε-machine (hence +:math:`C_\mu`) -- because they change only the presentation, not the process. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean + + In [2]: from sofic.generators.state_splitting import amalgamate, out_edges, split_state + + In [3]: machine = golden_mean(0.5) + + In [4]: edges = out_edges(machine, "A") + + In [5]: split = split_state(machine, "A", [edges[:1], edges[1:]]) + + In [6]: sorted(split.states(), key=repr), split.is_equal_process(machine) + + In [7]: copies = frozenset({("A", 0), ("A", 1)}) + + In [8]: merged = amalgamate(split, [copies, {"B"}], labels={copies: "A"}) + + In [9]: sorted(merged.states()), merged.initial_distribution + +API +=== + +.. autofunction:: split_state +.. autofunction:: amalgamate +.. autofunction:: out_edges +.. autofunction:: in_edges diff --git a/docs/generators/support.rst b/docs/generators/support.rst new file mode 100644 index 0000000..54675e1 --- /dev/null +++ b/docs/generators/support.rst @@ -0,0 +1,77 @@ +.. support.rst +.. py:module:: sofic.generators.support + +******************* +Process supports +******************* + +The *stationary support* of an HMM is the set of finite words +:math:`w = x_0 \cdots x_{n-1}` with :math:`P(X_{0:n} = w) > 0` under the +stationary law. It is a factorial regular language, recognized by the NFA +:func:`support_nfa`: the states of positive stationary probability, all initial +and accepting, joined by the positive-probability transitions. + +:meth:`~sofic.generators.base.HiddenMarkovModel.to_support_nfa` is a different +automaton: it starts from *every* state, transient ones included, so it also +accepts words that only a transient state can emit. Those words have stationary +probability zero and are rejected by :func:`support_nfa`. + +Support decisions +================= + +Comparisons of supports are inclusions of regular languages, decided by +antichain inclusion without determinizing :cite:`DeWulf2006`. + +:func:`support_includes` ``(p, q)`` + Every word of positive :math:`P` probability has positive :math:`Q` + probability. + +:func:`is_absolutely_continuous` ``(p, q)`` + :math:`P_{0:n} \ll Q_{0:n}` for every block length :math:`n`, i.e. + absolute continuity on the finite-dimensional cylinders. This is the same + decision as :func:`support_includes`, and it is exactly the condition for the + block divergences :math:`D(P_{0:n} \| Q_{0:n})` and the + :doc:`relative entropy rate ` to be finite. It is not + absolute continuity of the laws on infinite sequences: distinct ergodic + stationary measures are mutually singular there. + +:func:`support_equal` ``(p, q)`` + Inclusion both ways: the two processes allow the same finite words, whatever + their probabilities. + +The golden mean process forbids ``11``, so it is absolutely continuous with +respect to the fair coin but not the other way round. + +.. ipython:: + + In [1]: from sofic.examples import bernoulli, golden_mean + + In [2]: from sofic.generators.support import is_absolutely_continuous, support_equal, support_includes + + @doctest + In [3]: support_includes(golden_mean(0.5), bernoulli()) + Out[3]: True + + @doctest + In [4]: is_absolutely_continuous(bernoulli(), golden_mean(0.5)) + Out[4]: False + + @doctest + In [5]: support_equal(golden_mean(0.3), golden_mean(0.8)) + Out[5]: True + +The same decisions are available as +:meth:`~sofic.generators.base.HiddenMarkovModel.support_includes`, +:meth:`~sofic.generators.base.HiddenMarkovModel.is_absolutely_continuous`, and +:meth:`~sofic.generators.base.HiddenMarkovModel.support_equal`. + +API +=== + +.. autofunction:: support_nfa + +.. autofunction:: support_includes + +.. autofunction:: is_absolutely_continuous + +.. autofunction:: support_equal diff --git a/docs/inference/active_epsilon_learning.rst b/docs/inference/active_epsilon_learning.rst new file mode 100644 index 0000000..32f12c8 --- /dev/null +++ b/docs/inference/active_epsilon_learning.rst @@ -0,0 +1,123 @@ +.. active_epsilon_learning.rst +.. py:module:: sofic.inference.active + +**************************************** +Active learning of ε-machines +**************************************** + +Sample-based reconstruction (:doc:`cssr`, :doc:`spectral`) estimates an +ε-machine :cite:`Shalizi2001` from data. Active learning instead asks a +*teacher* questions, in the manner of Angluin's L\* for regular languages +:cite:`Angluin1987` (see :func:`sofic.automata.learning.active.learn_dfa_lstar`). +The teacher answers + +* **probability queries** ``oracle.word_probability(w)``: the stationary + probability :math:`P(w)` of a finite word, and optionally +* **equivalence queries** ``oracle.equivalent(machine)``: ``None`` when + ``machine`` generates the target process, otherwise a counterexample word + whose probability differs. + +:class:`ProcessOracle` answers both exactly from any finite HMM: probabilities +from the symbol-labeled matrices and the stationary state law (in sympy +arithmetic for symbolic models), equivalence with +:func:`~sofic.generators.process_equivalence.is_equal_process`. + +.. code-block:: python + + from sofic.examples import nemo_process + from sofic.inference.active import ProcessOracle, learn_epsilon_machine_active + + target = nemo_process(0.3, 0.6) + oracle = ProcessOracle(target) + machine = learn_epsilon_machine_active(oracle, ("0", "1")) + len(list(machine.states())) # 3 + machine.is_equal_process(target) # True + oracle.probability_queries, oracle.equivalence_queries # (37, 2) + +The observation table +===================== + +Rows are histories :math:`h` (prefixes), columns are test suffixes :math:`s`, and +the entry is the conditional future probability :math:`P(s \mid h) = P(hs)/P(h)`. +A row is the finite-history *mixed state* :cite:`Ellison2009` seen through the +tests. Tests start as the single symbols, so every row carries its next-symbol +distribution. Rows match when their entries differ by at most ``tol``; answers +that are ``int``, :class:`~fractions.Fraction` or sympy are compared exactly. + +The learner repeats: + +1. **Close** the table: every positive-probability one-symbol extension + :math:`hx` of a state row must match a state row; otherwise :math:`hx` + becomes a state row. State rows stay pairwise distinct, and adding tests + only refines rows, so the table is also *consistent*. +2. **Hypothesize**: one state per state row, the transition from :math:`h` on + :math:`x` going to the row matching :math:`hx` with probability + :math:`P(x \mid h)`. The hypothesis is unifilar by construction and starts at + the row of the seed ``history`` (empty by default); when seeded, the oracle is + asked ``equivalent(machine, history=history)`` and compares with the process + conditioned on the seed. +3. **Query** equivalence. If the oracle has no ``equivalent``, every word of + length at most ``max_counterexample_length`` is compared instead. +4. **Refine**: every suffix of the counterexample becomes a test. + +The returned ε-machine is the hypothesis restricted to its closed class, started +at its stationary distribution; the transient rows (histories that have not yet +synchronized) are dropped. + +Guarantees +========== + +Suppose the answers and the equivalence oracle are exact, and let :math:`N` be +the number of distinct predictive distributions :math:`P(\cdot \mid h)` over the +positive-probability histories :math:`h` extending the seed ``history``. A closed +table whose tests include every suffix of a word :math:`w = x_0 \cdots x_{L-1}` +predicts :math:`P(w)` correctly, by the chain rule +:math:`P(w) = \prod_{t=0}^{L-1} P(x_t \mid x_{0:t})` and the row match of each +extension on the suffix :math:`x_{t+1:L}`. So each counterexample adds at least +one state row, and when :math:`N` is finite the learner stops after at most +:math:`N` equivalence queries with exactly :math:`N` hypothesis states. Rows are +needed for the :math:`N` state rows and their one-symbol extensions, answers are +cached, and so with an equivalence oracle the number of probability queries is +at most +:math:`2 N (|\mathcal{A}| + 1)(|\mathcal{A}| + \sum_k |w_k|)` for counterexamples +:math:`w_k` -- polynomial in :math:`N`, the alphabet size and the counterexample +lengths. + +If the target is a stationary ergodic process whose ε-machine is finite and +exactly synchronizable :cite:`Travers2010`, every predictive state reaches a +causal state, so the closed class of the final hypothesis *is* the ε-machine: +the same number of states, the same process, and the same :math:`C_\mu` and +:math:`h_\mu`. + +:math:`N` is finite when + +* the mixed-state presentation from the stationary law is finite (golden mean, + even, Nemo, butterfly, restricted golden mean, ...), or +* the seed ``history`` is a synchronizing word: every row is then a causal + state, :math:`N` is the number of causal states, and the hypothesis is the + ε-machine itself. + +A finite ε-machine that is not exactly synchronizable has infinitely many mixed +states, and so can an exactly synchronizable one along a non-synchronizing +transient (for example a symbol that fixes two states with different +probabilities). The table then never closes, and the learner raises +:class:`~sofic.exceptions.MixedStateExplosionError` at ``max_states``; seed a +synchronizing ``history`` in that case. More than ``max_rounds`` equivalence +queries raises :class:`RuntimeError`. + +With float answers the result is approximate. Matching within ``tol`` is not +transitive; converging transient mixed states are merged once they agree to +``tol``, which keeps the table finite at the cost of extra rows; and a +counterexample that separates no rows at ``tol`` ends learning with the current +hypothesis, which then merges causal states whose predictions differ by less than +``tol``. + +API +=== + +.. autofunction:: learn_epsilon_machine_active + +.. autoclass:: ProcessOracle + :members: word_probability, equivalent + +.. autoclass:: ProbabilityOracle diff --git a/docs/inference/inference.rst b/docs/inference/inference.rst index b0105a0..c56c795 100644 --- a/docs/inference/inference.rst +++ b/docs/inference/inference.rst @@ -48,6 +48,7 @@ The historical names ``InferMC`` and ``InferEM`` are retained as aliases for markov epsilon spectral + active_epsilon_learning model_selection diagnostics hdp_hmm diff --git a/docs/references.bib b/docs/references.bib index 1456962..4b7b3ca 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1300,3 +1300,86 @@ @incollection{blackwell1957entropy address = {Prague}, year = {1957}, } + +@phdthesis{Upper1997, + author = {Upper, Daniel R.}, + title = {Theory and Algorithms for Hidden {Markov} Models and Generalized Hidden {Markov} Models}, + school = {University of California, Berkeley}, + year = {1997}, +} + +@article{Schutzenberger1961, + author = {Sch{\"u}tzenberger, M. P.}, + title = {On the Definition of a Family of Automata}, + journal = {Information and Control}, + volume = {4}, + number = {2--3}, + pages = {245--270}, + year = {1961}, + doi = {10.1016/S0019-9958(61)80020-X}, +} + +@article{LarsenSkou1991, + author = {Larsen, Kim G. and Skou, Arne}, + title = {Bisimulation through Probabilistic Testing}, + journal = {Information and Computation}, + volume = {94}, + number = {1}, + pages = {1--28}, + year = {1991}, + doi = {10.1016/0890-5401(91)90030-6}, +} + +@article{Shalizi2001, + author = {Shalizi, Cosma Rohilla and Crutchfield, James P.}, + title = {Computational Mechanics: Pattern and Prediction, Structure and Simplicity}, + journal = {Journal of Statistical Physics}, + volume = {104}, + number = {3/4}, + pages = {817--879}, + year = {2001}, + doi = {10.1023/A:1010388907793}, +} + +@article{CrutchfieldFeldman2003, + author = {Crutchfield, James P. and Feldman, David P.}, + title = {Regularities Unseen, Randomness Observed: Levels of Entropy Convergence}, + journal = {Chaos: An Interdisciplinary Journal of Nonlinear Science}, + volume = {13}, + number = {1}, + pages = {25--54}, + year = {2003}, + doi = {10.1063/1.1530990}, +} + +@article{LohrAy2009, + author = {L{\"o}hr, Wolfgang and Ay, Nihat}, + title = {On the Generative Nature of Prediction}, + journal = {Advances in Complex Systems}, + volume = {12}, + number = {2}, + pages = {169--194}, + year = {2009}, + doi = {10.1142/S0219525909002143}, +} + +@article{Rached2001, + author = {Rached, Ziad and Alajaji, Fady and Campbell, L. Lorne}, + title = {R{\'e}nyi's Divergence and Entropy Rates for Finite Alphabet {Markov} Sources}, + journal = {IEEE Transactions on Information Theory}, + volume = {47}, + number = {4}, + pages = {1553--1561}, + year = {2001}, + doi = {10.1109/18.923736}, +} + +@book{DemboZeitouni1998, + author = {Dembo, Amir and Zeitouni, Ofer}, + title = {Large Deviations Techniques and Applications}, + edition = {2}, + publisher = {Springer}, + address = {New York}, + year = {1998}, + doi = {10.1007/978-1-4612-5320-4}, +} diff --git a/sofic/automata/dfa.py b/sofic/automata/dfa.py index 382b2c1..3619f5b 100644 --- a/sofic/automata/dfa.py +++ b/sofic/automata/dfa.py @@ -44,9 +44,10 @@ def add_transition(self, source: Hashable, target: Hashable, symbol: Any, **attr if symbol is EPSILON: raise NonDeterministicError("DFA does not allow epsilon transitions") key = (source, symbol) - for transition in self.graph.out_transitions(source): - if transition.data.get(ATTR_SYMBOL) == symbol and transition.target != target: - raise NonDeterministicError(f"non-deterministic transition on {key}") + if self.graph.has_state(source): + for transition in self.graph.out_transitions(source): + if transition.data.get(ATTR_SYMBOL) == symbol and transition.target != target: + raise NonDeterministicError(f"non-deterministic transition on {key}") return self.graph.add_transition(source, target, **{ATTR_SYMBOL: symbol, **attrs}) def recognizes(self, word: Sequence[Any]) -> bool: diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py index db8de0c..810786e 100644 --- a/sofic/generators/__init__.py +++ b/sofic/generators/__init__.py @@ -9,6 +9,7 @@ driven_entropy_rate, ) from sofic.generators.complexity_dimension import StatisticalComplexityDimension, statistical_complexity_dimension +from sofic.generators.conditioning import condition_on_language from sofic.generators.correlations import autocorrelation, mutual_information_function, power_spectrum from sofic.generators.directional_flow import ( directed_information, @@ -20,7 +21,15 @@ ) from sofic.generators.epsilon_machine import EpsilonMachine from sofic.generators.epsilon_transducer import EpsilonTransducer -from sofic.generators.lumping import LumpabilityError, is_lumpable, lump, normalize_partition +from sofic.generators.factor_codes import higher_block, image_process +from sofic.generators.lumping import ( + LumpabilityError, + bisimulation_partition, + coarsest_lumping, + is_lumpable, + lump, + normalize_partition, +) from sofic.generators.markov import MarkovChain from sofic.generators.mealy import MealyHMM from sofic.generators.measures import EntropyRateEstimate, entropy_rate_blackwell, entropy_rate_bounds @@ -34,9 +43,11 @@ minimal_generative_model, wyner_generative_model, ) +from sofic.generators.minimal_quasi_realization import minimal_quasi_realization, process_rank from sofic.generators.mixed_state import MixedState, MixedStatePresentation from sofic.generators.moore import MooreHMM from sofic.generators.nmachine import NMachine +from sofic.generators.omega_probability import omega_probability, regular_language_probability from sofic.generators.pfa import ProbabilisticFiniteAutomaton from sofic.generators.predictive_rd import PredictiveRateDistortionCurve, predictive_rate_distortion from sofic.generators.quasi_realization import QuasiRealization @@ -45,7 +56,10 @@ relative_entropy_rate, relative_entropy_rate_bounds, ) +from sofic.generators.renyi import pressure, rate_function, renyi_entropy_rate from sofic.generators.stack_hmm import HiddenMarkovStackModel +from sofic.generators.state_splitting import amalgamate, split_state +from sofic.generators.support import is_absolutely_continuous, support_equal, support_includes, support_nfa from sofic.generators.topological_epsilon_enumeration import ( count_topological_epsilon_machines, epsilon_machine_to_idfa_string, @@ -136,4 +150,22 @@ "entropy_rate_blackwell", "entropy_rate_bounds", "statistical_complexity_dimension", + "bisimulation_partition", + "coarsest_lumping", + "minimal_quasi_realization", + "process_rank", + "amalgamate", + "condition_on_language", + "higher_block", + "image_process", + "is_absolutely_continuous", + "omega_probability", + "regular_language_probability", + "split_state", + "support_equal", + "support_includes", + "support_nfa", + "pressure", + "rate_function", + "renyi_entropy_rate", ] diff --git a/sofic/generators/base.py b/sofic/generators/base.py index b0d7fcc..862effa 100644 --- a/sofic/generators/base.py +++ b/sofic/generators/base.py @@ -279,6 +279,15 @@ def relative_entropy_rate(self, other: HiddenMarkovModel) -> float: return relative_entropy_rate(self, other) + def higher_block(self, k: int) -> MealyHMM: + """Return the ``k``-block presentation, emitting ``Y_t = X_{t:t+k}`` :cite:`LindMarcus1995`. + + See :func:`~sofic.generators.factor_codes.higher_block`. + """ + from sofic.generators.factor_codes import higher_block + + return higher_block(self, k) + def to_sofic_shift(self) -> SoficShift: """Strip probabilities and return a sofic shift with the same support.""" from sofic.generators.conversions import hmm_to_sofic_shift @@ -297,6 +306,42 @@ def to_support_dfa(self) -> DFA: return hmm_to_support_dfa(self) + def process_rank(self, *, tol: float = 1e-9) -> int: + """Return the process rank (Hankel rank of word probabilities) :cite:`Upper1997`. + + See :func:`~sofic.generators.minimal_quasi_realization.process_rank`. + """ + from sofic.generators.minimal_quasi_realization import process_rank + + return process_rank(self, tol=tol) + + def support_includes(self, other: HiddenMarkovModel) -> bool: + """Return whether every positive-probability word of ``self`` has positive ``other`` probability. + + See :func:`~sofic.generators.support.support_includes` :cite:`DeWulf2006`. + """ + from sofic.generators.support import support_includes + + return support_includes(self, other) + + def is_absolutely_continuous(self, other: HiddenMarkovModel) -> bool: + """Return whether ``self`` is absolutely continuous w.r.t. ``other`` on every finite block. + + See :func:`~sofic.generators.support.is_absolutely_continuous`. + """ + from sofic.generators.support import is_absolutely_continuous + + return is_absolutely_continuous(self, other) + + def support_equal(self, other: HiddenMarkovModel) -> bool: + """Return whether ``self`` and ``other`` have the same stationary finite-word support. + + See :func:`~sofic.generators.support.support_equal`. + """ + from sofic.generators.support import support_equal + + return support_equal(self, other) + class QuasiStochasticModel(StateMachine): """Generator allowing signed quasiprobabilities on internal weights.""" diff --git a/sofic/generators/canonical_residual.py b/sofic/generators/canonical_residual.py new file mode 100644 index 0000000..0f960cf --- /dev/null +++ b/sofic/generators/canonical_residual.py @@ -0,0 +1,183 @@ +"""Canonical residual HMM: the process analogue of the canonical RFSA. + +The canonical residual finite-state automaton of a regular language has the +*prime* residuals as states, those that are not a union of other residuals +:cite:`Denis2002` (see also :cite:`MaarandTamm2022`). For a stationary process +with a finite ε-machine, the residual of a causal state ``s`` is its future +morph :math:`f_s(w) = P(w \\mid s)`. Replacing "union" by "nonnegative +combination", the prime morphs are the extreme rays of the cone spanned by the +morphs. Every morph satisfies :math:`f_s(\\lambda) = 1` on the empty word, so a +nonnegative combination of morphs is a convex one. + +Writing each morph as :math:`f_t = \\sum_{r \\in R} c_{t r} f_r` over the extreme +set :math:`R`, the identity :math:`f_r(x w) = \\sum_t T^{(x)}_{r t} f_t(w)` gives +a generator on :math:`R` with joint transition matrices +:math:`T^{(x)} C` restricted to the rows in :math:`R`, and initial law +:math:`\\pi C`. By induction on word length state ``r`` of this generator emits +the futures :math:`f_r`, so it generates the same process with at most as many +states as the ε-machine. +""" + +from __future__ import annotations + +from collections.abc import Hashable +from typing import TYPE_CHECKING, Any + +import numpy as np +from scipy.optimize import linprog + +from sofic.generators.mealy import MealyHMM + +if TYPE_CHECKING: + from sofic.generators.base import HiddenMarkovModel + +__all__ = ["canonical_residual_hmm"] + + +def _unifilar_support( + model: HiddenMarkovModel, max_states: int +) -> tuple[list[Hashable], np.ndarray, dict[Any, np.ndarray]]: + from sofic.generators.epsilon_machine import EpsilonMachine + from sofic.generators.matrices import symbol_matrices + from sofic.generators.renyi import _stationary_support + + mealy = model.to_mealy() + if not mealy.is_unifilar(): + mealy = EpsilonMachine.from_hmm(mealy, max_states=max_states) + states = list(mealy.reindex().states) + pi = np.asarray(mealy.stationary_distribution(), dtype=float) + raw = {symbol: np.asarray(m, dtype=float) for symbol, m in symbol_matrices(mealy).items()} + keep = _stationary_support(pi, list(raw.values())) + matrices = {symbol: m[np.ix_(keep, keep)] for symbol, m in raw.items()} + return [s for s, k in zip(states, keep, strict=True) if k], pi[keep] / pi[keep].sum(), matrices + + +def _morph_coordinates(matrices: dict[Any, np.ndarray], size: int, tol: float) -> np.ndarray: + """Return the state-by-test-word matrix ``F[s, j] = P(w_j | s)`` on a basis of test words. + + Test words are found breadth first from the empty word, keeping ``w`` when + :math:`T^{(w)} \\mathbf{1}` is independent of those kept so far. Their span + is the span of all :math:`T^{(w)} \\mathbf{1}`, so linear relations among + rows of ``F`` are exactly linear relations among the morphs. + """ + basis: list[np.ndarray] = [] + queue = [np.ones(size)] + while queue and len(basis) < size: + vector = queue.pop(0) + candidate = np.column_stack([*basis, vector]) + if np.linalg.matrix_rank(candidate, tol=tol * max(1.0, np.abs(candidate).max())) > len(basis): + basis.append(vector) + queue.extend(m @ vector for m in matrices.values()) + return np.column_stack(basis) + + +def _conic_residual(target: np.ndarray, generators: np.ndarray) -> tuple[float, np.ndarray]: + """Return ``min ||generators @ c - target||_1`` over ``c >= 0`` and a minimizer.""" + rows, cols = generators.shape + if cols == 0: + return float(np.abs(target).sum()), np.zeros(0) + identity = np.eye(rows) + result = linprog( + np.concatenate([np.zeros(cols), np.ones(2 * rows)]), + A_eq=np.hstack([generators, identity, -identity]), + b_eq=target, + bounds=(0, None), + method="highs", + ) + if not result.success: + raise RuntimeError(f"linear program failed: {result.message}") + return float(result.fun), result.x[:cols] + + +def canonical_residual_hmm( + machine: HiddenMarkovModel, + *, + experimental: bool = False, + tol: float = 1e-7, + max_states: int = 10_000, +) -> MealyHMM: + """Return the canonical residual HMM of a process with a finite ε-machine. + + .. warning:: + Experimental, and must be enabled with ``experimental=True``. The + construction is a process analogue of the canonical RFSA + :cite:`Denis2002,MaarandTamm2022` with no published reference; it is + verified against brute-force word distributions in the tests only. + + The states are the causal states whose future morphs are extreme rays of + the cone spanned by all morphs (see the module docstring). Morphs are + compared through their probabilities on a basis of test words, and + extremality and the combination coefficients are decided by + :func:`scipy.optimize.linprog` with residual tolerance ``tol`` (floating + point; no exact arithmetic). The default sits at the HiGHS feasibility + tolerance; morphs closer than ``tol`` are treated as equal, so ε-machines + whose states accumulate numerically (finite approximations of infinite + ones) lose those states. Causal states with identical morphs are + merged, so a unifilar but non-minimal presentation gives the same result + as its ε-machine. The result is a possibly non-unifilar + :class:`~sofic.generators.mealy.MealyHMM` on the extreme states, keeping + their labels, that generates the same process with at most as many states + as the ε-machine. When the morphs are affinely independent (for example + when the ε-machine has two states) it is the ε-machine itself. + + Since it is a generator of the process, its state entropy is an upper + bound on the minimal state entropy of any generating HMM, and its state + count an upper bound on the minimal generator size studied by Löhr and Ay + :cite:`LohrAy2009`, which can lie strictly below the ε-machine. No lower + bound or minimality claim is made. + + Non-unifilar input is converted with + :meth:`EpsilonMachine.from_hmm(machine, max_states=max_states) + `, which raises + :class:`~sofic.exceptions.MixedStateExplosionError` when the ε-machine is + not finite. + + Raises + ------ + RuntimeError + If ``experimental`` is not ``True``. + """ + if experimental is not True: + raise RuntimeError("canonical_residual_hmm is experimental; pass experimental=True to use it") + states, pi, matrices = _unifilar_support(machine, max_states) + size = len(states) + morphs = _morph_coordinates(matrices, size, tol) + scale = max(1.0, float(np.abs(morphs).max())) + + representatives: list[int] = [] + for i in range(size): + if not any(np.abs(morphs[i] - morphs[j]).sum() <= tol * scale for j in representatives): + representatives.append(i) + extreme = [ + i + for i in representatives + if _conic_residual(morphs[i], morphs[[j for j in representatives if j != i]].T)[0] > tol * scale + ] + + coefficients = np.zeros((size, len(extreme))) + for t in range(size): + if t in extreme: + coefficients[t, extreme.index(t)] = 1.0 + continue + _residual, weights = _conic_residual(morphs[t], morphs[extreme].T) + coefficients[t] = weights / weights.sum() + + initial = pi @ coefficients + hmm = MealyHMM( + initial_distribution={states[r]: float(mass) for r, mass in zip(extreme, initial / initial.sum(), strict=True)}, + observation_alphabet=frozenset(matrices), + ) + for r in extreme: + hmm.graph.add_state(states[r]) + joint = {symbol: (m @ coefficients)[extreme] for symbol, m in matrices.items()} + for a, r in enumerate(extreme): + edges = [ + (symbol, b, block[a, b]) + for symbol, block in joint.items() + for b in range(len(extreme)) + if block[a, b] > tol + ] + total = sum(prob for *_, prob in edges) + for symbol, b, prob in edges: + hmm.add_transition(states[r], states[extreme[b]], symbol, float(prob / total)) + return hmm diff --git a/sofic/generators/conditioning.py b/sofic/generators/conditioning.py new file mode 100644 index 0000000..2a2775d --- /dev/null +++ b/sofic/generators/conditioning.py @@ -0,0 +1,161 @@ +"""Conditioning a stationary process on a regular constraint. + +Given an HMM with symbol matrices ``T^(x)`` and a DFA, the words whose every +prefix the DFA accepts are those whose run stays in accepting DFA states. On +the product of HMM states and accepting DFA states this is the substochastic +matrix + +.. math:: + + M^{(x)}_{(i, q), (j, \\delta(q, x))} = T^{(x)}_{ij} + \\quad \\text{whenever } \\delta(q, x) \\text{ is accepting.} + +Conditioning on surviving ``n`` steps and letting ``n -> infinity`` gives +the Doob ``h``-transform of ``M`` by its right Perron eigenvector ``h``, + +.. math:: + + P'(x, j \\mid i) = \\frac{M^{(x)}_{ij} h_j}{\\lambda h_i}, + +the construction of Parry's measure of maximal entropy :cite:`Parry1964` +(:cite:`LindMarcus1995`, §13.3) applied to a weighted rather than a 0-1 matrix. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Hashable +from typing import TYPE_CHECKING, Any + +import networkx as nx +import numpy as np + +from sofic.exceptions import NonDeterministicError +from sofic.graph import ATTR_SYMBOL + +if TYPE_CHECKING: + from sofic.automata.dfa import DFA + from sofic.generators.base import HiddenMarkovModel + from sofic.generators.mealy import MealyHMM + +_TOL = 1e-12 +_RTOL = 1e-9 + + +def condition_on_language(hmm: HiddenMarkovModel, dfa: DFA) -> MealyHMM: + """Return the stationary process ``hmm`` conditioned to stay in ``dfa``'s language forever. + + The conditioning event for length ``n`` is that every prefix of + ``X_{0:n}`` is accepted by ``dfa``, i.e. its run never leaves the accepting + states. For a prefix-closed language this is just ``X_{0:n} ∈ L``, and for a + factorial one (the language of a shift, such as an SFT) it is moreover + shift invariant, so the result is the ``n -> infinity`` limit of the law of + windows ``X_{k:k+m}`` given ``X_{0:n} ∈ L`` with ``k`` and ``n - k`` both + large. A non-prefix-closed ``dfa`` is effectively replaced by the largest + prefix-closed sublanguage of its language. + + ``dfa`` is completed over the union of both alphabets (the trap state is + rejecting), and its product with the stationary-positive states of ``hmm`` + is restricted to accepting DFA states. Surviving paths of length ``n`` + concentrate on the strongly connected component of the restricted matrix + with the largest Perron root ``λ``; this component is required to be + unique. On it the result is the Doob ``h``-transform + + .. math:: + + P'(x, (j, q') \\mid (i, q)) = \\frac{T^{(x)}_{ij} h_{(j, q')}}{\\lambda h_{(i, q)}}, + + with ``h`` the right Perron eigenvector, started in its stationary law + ``u_s h_s`` (``u`` the left Perron eigenvector). This is the Parry + construction :cite:`Parry1964` (:cite:`LindMarcus1995`, §13.3): conditioning + the uniform i.i.d. process on an irreducible SFT returns its measure of + maximal entropy. When the component is periodic the window limit holds in + the Cesàro sense. + + States of the returned HMM are pairs ``(hmm state, dfa state)``; its + observation alphabet is that of ``hmm``. + + Raises + ------ + NonDeterministicError + If ``dfa`` is not deterministic. + ValueError + If no infinite word survives the constraint, or several components + share the maximal Perron root, so the limit depends on how the + constraint is entered. + """ + from sofic.automata.algorithms import complete + from sofic.generators.matrices import emission_tensors + from sofic.generators.mealy import MealyHMM + from sofic.shifts.parry_construction import _perron_pair + + if not dfa.is_deterministic(): + raise NonDeterministicError("condition_on_language needs a deterministic automaton") + mealy = hmm.to_mealy() + states = tuple(mealy.reindex().states) + pi, raw = emission_tensors(mealy, policy="stationary") + matrices = {symbol: np.asarray(matrix, dtype=float) for symbol, matrix in raw.items()} + pi = np.asarray(pi, dtype=float) + dfa = complete(dfa, frozenset(matrices) | dfa.input_alphabet) + delta = {(t.source, t.data[ATTR_SYMBOL]): t.target for t in dfa.transitions()} + accepting = dfa.accepting_states + (start,) = dfa.initial_states + if start not in accepting: + raise ValueError("the empty word is rejected, so no word survives the constraint") + + nodes: dict[tuple[int, Hashable], int] = {} + frontier: deque[tuple[int, Hashable]] = deque() + for i in np.flatnonzero(pi > _TOL): + nodes[(int(i), start)] = len(nodes) + frontier.append((int(i), start)) + edges: list[tuple[int, int, Any, float]] = [] + while frontier: + i, q = frontier.popleft() + for symbol, matrix in matrices.items(): + target = delta[(q, symbol)] + if target not in accepting: + continue + for j in np.flatnonzero(matrix[i] > _TOL): + node = (int(j), target) + if node not in nodes: + nodes[node] = len(nodes) + frontier.append(node) + edges.append((nodes[(i, q)], nodes[node], symbol, float(matrix[i, j]))) + + chain = np.zeros((len(nodes), len(nodes))) + for source, target, _symbol, prob in edges: + chain[source, target] += prob + graph = nx.DiGraph() + graph.add_nodes_from(range(len(nodes))) + graph.add_edges_from((source, target) for source, target, _symbol, _prob in edges) + radii = [] + for members in nx.strongly_connected_components(graph): + component = sorted(members) + block = chain[np.ix_(component, component)] + radii.append((float(np.max(np.abs(np.linalg.eigvals(block)))), component)) + lam = max((radius for radius, _ in radii), default=0.0) + if lam <= _TOL: + raise ValueError("no infinite word survives the constraint") + maximal = [component for radius, component in radii if radius >= lam * (1.0 - _RTOL)] + if len(maximal) > 1: + raise ValueError( + f"{len(maximal)} components share the maximal growth rate; the conditioned limit is not unique" + ) + (component,) = maximal + position = {node: k for k, node in enumerate(component)} + lam, left, right = _perron_pair(chain[np.ix_(component, component)]) + stationary = left * right + stationary /= stationary.sum() + + labels = {index: (states[i], q) for (i, q), index in nodes.items()} + result = MealyHMM( + initial_distribution={labels[node]: float(stationary[position[node]]) for node in component}, + observation_alphabet=frozenset(mealy.observation_alphabet), + ) + for node in component: + result.graph.add_state(labels[node]) + for source, target, symbol, prob in edges: + if source in position and target in position: + weight = prob * right[position[target]] / (lam * right[position[source]]) + result.add_transition(labels[source], labels[target], symbol, weight) + return result diff --git a/sofic/generators/conversions.py b/sofic/generators/conversions.py index b5abf5b..cd1f6dc 100644 --- a/sofic/generators/conversions.py +++ b/sofic/generators/conversions.py @@ -90,7 +90,12 @@ def hmm_to_support_nfa(hmm: HiddenMarkovModel) -> NFA: def hmm_to_support_dfa(hmm: HiddenMarkovModel) -> DFA: - """Determinize the HMM support NFA from the all-states subset.""" + """Determinize the HMM support NFA from the all-states subset. + + The language is the set of finite words some state can emit, which is + factorial, so every non-empty reachable subset accepts (the empty word + included). + """ support = _mealy_support(hmm) states = frozenset(support.states()) nfa = hmm_to_support_nfa(support) @@ -100,7 +105,7 @@ def hmm_to_support_dfa(hmm: HiddenMarkovModel) -> DFA: empty_subset = frozenset() if dfa.graph.has_state(empty_subset): dfa.graph.nx.remove_node(empty_subset) - dfa.accepting_states = dfa.graph.terminal_recurrent_states() + dfa.accepting_states = frozenset(dfa.states()) return dfa diff --git a/sofic/generators/epsilon_construction.py b/sofic/generators/epsilon_construction.py index 2163bb0..8d22c0a 100644 --- a/sofic/generators/epsilon_construction.py +++ b/sofic/generators/epsilon_construction.py @@ -9,7 +9,6 @@ from __future__ import annotations -from collections import defaultdict from typing import Any import numpy as np @@ -26,6 +25,7 @@ has_symbolic, is_positive_mass, is_zero, + probs_equal, simplify_prob, sum_probs, ) @@ -92,11 +92,36 @@ def _split_block( *, constraints: Any = None, ) -> list[set[Any]]: - """Split a block when states disagree on labeled successor blocks.""" - futures: dict[TransitionSignature, set[Any]] = defaultdict(set) - for state in block: - futures[_transition_signature(hmm, state, state_to_block, constraints=constraints)].add(state) - return list(futures.values()) + """Split a block when states disagree on labeled successor blocks. + + Numeric probabilities are compared within the tolerance the mixed-state + construction uses to merge beliefs; float round-off from belief updates + would otherwise split states with identical futures. + """ + groups: list[tuple[TransitionSignature, set[Any]]] = [] + for state in sorted(block, key=repr): + signature = _transition_signature(hmm, state, state_to_block, constraints=constraints) + for representative, members in groups: + if _signatures_match(representative, signature): + members.add(state) + break + else: + groups.append((signature, {state})) + return [members for _, members in groups] + + +def _signatures_match(left: TransitionSignature, right: TransitionSignature) -> bool: + if len(left) != len(right): + return False + for (symbol_a, block_a, prob_a), (symbol_b, block_b, prob_b) in zip(left, right, strict=True): + if symbol_a != symbol_b or block_a != block_b: + return False + if isinstance(prob_a, float) and isinstance(prob_b, float): + if not probs_equal(prob_a, prob_b): + return False + elif prob_a != prob_b: + return False + return True def _transition_signature( @@ -119,7 +144,7 @@ def _transition_signature( canonical_prob_key(prob, constraints), ) ) - return tuple(sorted(triples, key=repr)) + return tuple(sorted(triples, key=lambda triple: (repr(triple[0]), triple[1]))) def _quotient_machine( diff --git a/sofic/generators/epsilon_machine.py b/sofic/generators/epsilon_machine.py index 24efed7..237e267 100644 --- a/sofic/generators/epsilon_machine.py +++ b/sofic/generators/epsilon_machine.py @@ -47,6 +47,12 @@ def from_hmm(cls, hmm: MealyHMM | MooreHMM, *, max_states: int = 10_000) -> Epsi ``max_states`` caps the mixed-state enumeration used for non-unifilar input; exceeding it raises :class:`~sofic.exceptions.MixedStateExplosionError`. + + With floating-point probabilities, beliefs within the mixed-state + tolerance (``rtol=1e-9``, ``atol=1e-12``) are merged. A belief set that + converges to a limit point -- an infinite ε-machine -- can therefore + close at a finite approximation instead of raising. Use exact + (sympy) probabilities when the answer must be exact. """ from sofic.generators.epsilon_construction import build_epsilon_machine diff --git a/sofic/generators/factor_codes.py b/sofic/generators/factor_codes.py new file mode 100644 index 0000000..ce8ef90 --- /dev/null +++ b/sofic/generators/factor_codes.py @@ -0,0 +1,170 @@ +"""Sliding block codes applied to stochastic processes. + +A sliding block code ``Phi`` with memory ``m``, anticipation ``a`` and window +``w = m + 1 + a`` maps a process ``X`` to the *image process* ``Y = Phi(X)``, + +.. math:: + + Y_t = \\Phi(X_{t:t+w}), + +i.e. ``(Phi X)_{t+m}`` in the two-sided indexing of :cite:`LindMarcus1995` +(§1.5). For a stationary ``X`` the index shift by ``m`` is immaterial, and the +law of ``Y_{0:n}`` is the pushforward of the law of ``X_{0:n+w-1}`` under the +word map ``x_{0:n+w-1} -> Phi(x_{0:n+w-1})``. A stationary coding cannot +increase the entropy rate, so a conjugacy (an invertible code) preserves it +:cite:`Gray1990` (the topological analogue is :cite:`LindMarcus1995`, §4.1); +neither the excess entropy nor the statistical complexity is a conjugacy +invariant (see :func:`higher_block`). +""" + +from __future__ import annotations + +from collections import defaultdict +from collections.abc import Hashable +from typing import TYPE_CHECKING, Any + +from sofic.automata.transducer_operations import transduce_generator +from sofic.generators.base import HiddenMarkovModel +from sofic.generators.mealy import MealyHMM +from sofic.graph import ATTR_EMISSION, ATTR_PROB + +if TYPE_CHECKING: + from sofic.shifts.sliding_block_code import SlidingBlockCode + + +def image_process(hmm: HiddenMarkovModel, code: SlidingBlockCode) -> MealyHMM: + """Return a Mealy HMM generating the image process ``Phi(X)`` of ``hmm``. + + The code's sliding-window transducer (:meth:`SlidingBlockCode.to_transducer + `, whose + state is the last ``w - 1`` input symbols) is composed with ``hmm`` via + :func:`~sofic.automata.transducer_operations.transduce_generator`. The + composed state ``(q, x_{t-w+1:t})`` must start *synchronized*: its initial + law is that of ``(Q_{w-1}, X_{0:w-1})`` under ``hmm``, so the first emitted + symbol is ``Y_0 = Phi(X_{0:w})`` and the image of ``X_{0:n+w-1}`` is + ``Y_{0:n}`` exactly :cite:`LindMarcus1995` (§1.5, §6.1). Only states + reachable from that initial law are kept. + + Parameters + ---------- + hmm + Generator of the source process ``X``; its initial distribution is used + as given (stationary input gives the stationary image). + code + The sliding block code ``Phi``. + + Returns + ------- + MealyHMM + States are ``(hmm_state, context)`` pairs with ``context`` a tuple of + the last ``w - 1`` source symbols; the observation alphabet is + ``code.output_alphabet``. + + Raises + ------ + ValueError + If ``hmm`` can emit a ``w``-block that is missing from + ``code.block_map``. + """ + gen = hmm.to_mealy() + composed = transduce_generator(code.to_transducer(), gen, complete=False, normalize=False) + + initial: dict[tuple[Hashable, tuple[Any, ...]], Any] = { + (state, ()): mass for state, mass in gen.initial_distribution.items() if mass + } + for _ in range(code.window - 1): + advanced: dict[tuple[Hashable, tuple[Any, ...]], Any] = defaultdict(int) + for (state, context), mass in initial.items(): + for edge in gen.graph.out_transitions(state): + advanced[(edge.target, (*context, edge.data.get(ATTR_EMISSION)))] += mass * edge.data.get(ATTR_PROB) + initial = {key: mass for key, mass in advanced.items() if mass} + + for key in initial: + if not composed.graph.has_state(key): + _state, context = key + raise ValueError(f"context {context!r} uses symbols outside the code's input alphabet") + reachable = composed.graph.forward_reachable(frozenset(initial)) + for state, context in reachable: + for edge in gen.graph.out_transitions(state): + block = (*context, edge.data.get(ATTR_EMISSION)) + if block not in code.block_map: + raise ValueError(f"block {block!r} of the process is not in block_map") + + image = MealyHMM(initial_distribution=initial, observation_alphabet=frozenset(code.output_alphabet)) + for state in composed.states(): + if state in reachable: + image.graph.add_state(state) + for transition in composed.transitions(): + if transition.source in reachable: + image.add_transition( + transition.source, + transition.target, + transition.data[ATTR_EMISSION], + transition.data[ATTR_PROB], + ) + return image + + +def higher_block(hmm: HiddenMarkovModel, k: int) -> MealyHMM: + """Return the ``k``-block presentation of ``hmm``'s process. + + The higher block code ``beta_k`` :cite:`LindMarcus1995` (§1.4, Example + 1.5.10) recodes ``X`` to the process of overlapping ``k``-tuples + + .. math:: + + Y_t = X_{t:t+k} = (X_t, \\ldots, X_{t+k-1}), + + so ``P(Y_{0:n} = y_{0:n}) = P(X_{0:n+k-1} = x_{0:n+k-1})`` when the tuples + overlap consistently (``y_{t+1}[:-1] == y_t[1:]``) and ``x_{0:n+k-1}`` is + the word they spell, and is zero otherwise. ``beta_k`` is a conjugacy: the + 1-block code ``y -> y[0]`` inverts it. Hence the entropy rate is unchanged, + but the excess entropy and statistical complexity grow: + + .. math:: + + \\mathbf{E}(Y) = \\mathbf{E}(X) + (k - 1)\\, h_\\mu(X), \\qquad + C_\\mu(Y) = H[\\mathcal{S}_0, X_{-(k-1):0}] + = C_\\mu(X) + H[X_{-(k-1):0} \\mid \\mathcal{S}_0], + + with ``S_0`` the causal state of ``X`` after ``X_{-1}``. The first follows + from ``H[Y_{0:L}] = H[X_{0:L+k-1}]`` and ``E = lim (H[X_{0:L}] - L h_mu)`` + :cite:`CrutchfieldFeldman2003`; the second because a ``Y``-past fixes both the + last ``k - 1`` symbols (read again by the next tuple) and the causal state + of ``X``. + + Parameters + ---------- + hmm + Generator of ``X``. + k + Block length, ``k >= 1`` (``k = 1`` recodes symbols to 1-tuples). + + Returns + ------- + MealyHMM + Unifilar whenever ``hmm`` is. + """ + from sofic.shifts.sliding_block_code import SlidingBlockCode + + if k < 1: + raise ValueError("k must be at least 1") + gen = hmm.to_mealy() + frontier: set[tuple[Hashable, tuple[Any, ...]]] = {(state, ()) for state in gen.states()} + for _ in range(k): + frontier = { + (edge.target, (*word, edge.data.get(ATTR_EMISSION))) + for state, word in frontier + for edge in gen.graph.out_transitions(state) + } + blocks = {word for _state, word in frontier} + code = SlidingBlockCode( + {block: block for block in blocks}, + memory=0, + anticipation=k - 1, + input_alphabet=gen.observation_alphabet, + ) + return image_process(gen, code) + + +__all__ = ["higher_block", "image_process"] diff --git a/sofic/generators/lumping.py b/sofic/generators/lumping.py index 766cf1b..2ed62ea 100644 --- a/sofic/generators/lumping.py +++ b/sofic/generators/lumping.py @@ -18,16 +18,24 @@ The public entry points are :func:`is_lumpable` (predicate) and :func:`lump` (constructor). ``lump`` raises :class:`~sofic.exceptions.LumpabilityError` when the partition is not strongly lumpable unless ``check=False``. + +The coarsest strongly lumpable partition is the probabilistic bisimulation +equivalence of Larsen & Skou (:cite:`LarsenSkou1991`); :func:`bisimulation_partition` +computes it by partition refinement and :func:`coarsest_lumping` builds the +corresponding quotient. """ from __future__ import annotations +from collections import deque from collections.abc import Callable, Hashable, Iterable, Mapping +from fractions import Fraction from typing import TYPE_CHECKING, Any, overload import numpy as np from sofic.exceptions import LumpabilityError +from sofic.generators.prob import _sympy, canonical_prob_key, is_symbolic from sofic.graph import ATTR_EMISSION, ATTR_EMISSION_DIST, ATTR_PROB, TransitionGraph if TYPE_CHECKING: @@ -345,7 +353,7 @@ def is_lumpable(model: StateMachine, partition: PartitionLike, *, rtol: float = @overload def lump( model: MarkovChain, - partition: PartitionLike, + partition: PartitionLike | None = ..., *, check: bool = ..., labels: LabelsLike | None = ..., @@ -357,7 +365,7 @@ def lump( @overload def lump( model: MooreHMM, - partition: PartitionLike, + partition: PartitionLike | None = ..., *, check: bool = ..., labels: LabelsLike | None = ..., @@ -369,7 +377,7 @@ def lump( @overload def lump( model: MealyHMM, - partition: PartitionLike, + partition: PartitionLike | None = ..., *, check: bool = ..., labels: LabelsLike | None = ..., @@ -381,7 +389,7 @@ def lump( @overload def lump( model: StateMachine, - partition: PartitionLike, + partition: PartitionLike | None = ..., *, check: bool = ..., labels: LabelsLike | None = ..., @@ -392,7 +400,7 @@ def lump( def lump( model: StateMachine, - partition: PartitionLike, + partition: PartitionLike | None = None, *, check: bool = True, labels: LabelsLike | None = None, @@ -417,7 +425,8 @@ def lump( The generator to lump. partition Blocks (iterable of iterables) or a state-to-block mapping; must cover - every state exactly. + every state exactly. ``None`` (default) uses the coarsest strongly + lumpable partition, :func:`bisimulation_partition` with ``tol=atol``. check When ``True`` (default), raise :class:`~sofic.exceptions.LumpabilityError` if ``partition`` is not strongly lumpable. When ``False``, build the model @@ -441,7 +450,199 @@ def lump( If ``partition`` is invalid or the resolved block labels collide. """ checker, builder = _dispatch(model) + if partition is None: + partition = bisimulation_partition(model, tol=atol) blocks = normalize_partition(model, partition) if check and not checker(model, blocks, rtol=rtol, atol=atol): raise LumpabilityError("partition is not strongly lumpable for this model") return builder(model, blocks, labels) + + +def bisimulation_partition( + model: StateMachine, *, initial: PartitionLike | None = None, tol: float = 1e-10 +) -> list[frozenset[Hashable]]: + """Return the coarsest strongly lumpable partition of ``model``'s states. + + Two states are probabilistically bisimilar (Larsen & Skou, + :cite:`LarsenSkou1991`) when, for every emitted symbol, they send equal + probability mass into every block of the bisimulation; the bisimulation + classes are therefore exactly the coarsest partition satisfying the strong + lumpability condition of :func:`is_lumpable` (:cite:`KemenySnell1976`), and + every lumpable partition refines it. A Moore HMM's initial partition groups + states with equal emission laws. A + :class:`~sofic.generators.markov.MarkovChain` carries no edge labels, so its + coarsest lumpable partition is the single block unless ``initial`` (for + example, the level sets of an observation function of the state) is given. + + The partition is computed by splitter-driven refinement in the style of + Hopcroft and Paige--Tarjan: each splitter set ``C`` divides every block by + the per-symbol mass its states send into ``C``, and when a block that has + already served as a splitter breaks up, all of its pieces except the largest + become splitters (the largest piece's masses are then determined by + subtraction). + + Probabilities that are all exact (sympy expressions, :class:`~fractions.Fraction` + or integers) are compared exactly; otherwise masses within ``tol`` of their + sorted neighbors are grouped together. + + Parameters + ---------- + model + A :class:`~sofic.generators.markov.MarkovChain`, + :class:`~sofic.generators.mealy.MealyHMM` (including + :class:`~sofic.generators.epsilon_machine.EpsilonMachine`), or + :class:`~sofic.generators.moore.MooreHMM`. + initial + Optional partition to refine; the result is the coarsest strongly + lumpable partition finer than it. Defaults to the single block. + tol + Absolute tolerance for numeric masses; ignored for exact input. + + Returns + ------- + list of frozenset + Blocks ordered as by :func:`normalize_partition`. + + Raises + ------ + TypeError + If ``model`` is not a supported generator type. + ValueError + If ``initial`` is not a valid partition of ``model``'s states. + """ + from sofic.generators.moore import MooreHMM + + _dispatch(model) + order = {state: index for index, state in enumerate(model.states())} + exact = _exact_converter(model) + + def split(block: frozenset[Hashable], value: Callable[[Hashable], Any]) -> list[frozenset[Hashable]]: + return _split_block(block, value, order, exact=exact, tol=tol) + + blocks = normalize_partition(model, [order] if initial is None else initial) + if isinstance(model, MooreHMM): + emissions = {state: _emission_dist(model, state) for state in order} + symbols = sorted({symbol for law in emissions.values() for symbol in law}, key=repr) + for symbol in symbols: + blocks = [piece for block in blocks for piece in split(block, lambda s, x=symbol: emissions[s].get(x, 0))] + + worklist: deque[frozenset[Hashable]] = deque([frozenset(order), *blocks]) + pending = set(worklist) + while worklist: + splitter = worklist.popleft() + pending.discard(splitter) + masses = _masses_into(model, splitter, exact=exact) + labels = sorted({label for row in masses.values() for label in row}, key=repr) + refined: list[frozenset[Hashable]] = [] + for block in blocks: + pieces = [block] + for label in labels: + pieces = [ + piece + for part in pieces + for piece in split(part, lambda s, x=label, m=masses: m.get(s, {}).get(x, 0)) + ] + if len(pieces) > 1: + if block in pending: + pending.discard(block) + worklist.remove(block) + new_splitters = pieces + else: + largest = max(pieces, key=len) + new_splitters = [piece for piece in pieces if piece is not largest] + for piece in new_splitters: + if piece not in pending: + pending.add(piece) + worklist.append(piece) + refined.extend(pieces) + blocks = refined + return normalize_partition(model, blocks) + + +def coarsest_lumping( + model: StateMachine, + *, + initial: PartitionLike | None = None, + tol: float = 1e-10, + labels: LabelsLike | None = None, +) -> StateMachine: + """Return the quotient of ``model`` by its :func:`bisimulation_partition`. + + The result is the smallest strongly lumped model (Larsen & Skou, + :cite:`LarsenSkou1991`; Kemeny & Snell, :cite:`KemenySnell1976`) and, for + hidden Markov models, generates the same observed process as ``model`` from + the summed initial law. ``initial`` and ``tol`` are as in + :func:`bisimulation_partition`; ``labels`` is as in :func:`lump`. + """ + return lump(model, bisimulation_partition(model, initial=initial, tol=tol), check=False, labels=labels) + + +def _exact_converter(model: StateMachine) -> Callable[[Any], Any] | None: + """Return the exact scalar type for ``model``'s probabilities, or ``None`` when any is inexact.""" + values: list[Any] = [transition.data.get(ATTR_PROB, 0) for transition in model.graph.transitions()] + for state in model.states(): + values.extend(_emission_dist(model, state).values()) + if any( + isinstance(value, bool) or not (isinstance(value, (int, Fraction)) or is_symbolic(value)) for value in values + ): + return None + if any(is_symbolic(value) for value in values): + return _sympy().sympify + if any(isinstance(value, Fraction) for value in values): + return Fraction + return None + + +def _masses_into( + model: StateMachine, splitter: frozenset[Hashable], *, exact: Callable[[Any], Any] | None +) -> dict[Hashable, dict[Any, Any]]: + """Return ``state -> {symbol: mass into splitter}`` (``symbol`` is ``None`` without edge emissions).""" + from sofic.generators.mealy import MealyHMM + + by_symbol = isinstance(model, MealyHMM) + masses: dict[Hashable, dict[Any, Any]] = {} + for state in model.states(): + for transition in model.graph.out_transitions(state): + if transition.target not in splitter: + continue + label = transition.data.get(ATTR_EMISSION) if by_symbol else None + prob = transition.data.get(ATTR_PROB, 0) + row = masses.setdefault(state, {}) + if exact is not None: + row[label] = row.get(label, exact(0)) + exact(prob) + else: + row[label] = row.get(label, 0.0) + float(prob) + return masses + + +def _split_block( + block: frozenset[Hashable], + value: Callable[[Hashable], Any], + order: Mapping[Hashable, int], + *, + exact: Callable[[Any], Any] | None, + tol: float, +) -> list[frozenset[Hashable]]: + members = sorted(block, key=order.__getitem__) + if len(members) == 1: + return [block] + groups: list[list[Hashable]] = [] + if exact is not None: + keyed: dict[Any, list[Hashable]] = {} + for state in members: + mass = exact(value(state)) + keyed.setdefault(canonical_prob_key(mass) if is_symbolic(mass) else mass, []).append(state) + groups = list(keyed.values()) + else: + ranked = sorted(members, key=lambda state: float(value(state))) + previous: float | None = None + for state in ranked: + current = float(value(state)) + if previous is None or current - previous > tol: + groups.append([]) + groups[-1].append(state) + previous = current + if len(groups) == 1: + return [block] + groups.sort(key=lambda group: min(order[state] for state in group)) + return [frozenset(group) for group in groups] diff --git a/sofic/generators/markov.py b/sofic/generators/markov.py index 60a2651..f930349 100644 --- a/sofic/generators/markov.py +++ b/sofic/generators/markov.py @@ -36,14 +36,18 @@ def is_lumpable(self, partition: PartitionLike, *, rtol: float = 1e-8, atol: flo def lump( self, - partition: PartitionLike, + partition: PartitionLike | None = None, *, check: bool = True, labels: LabelsLike | None = None, rtol: float = 1e-8, atol: float = 1e-10, ) -> MarkovChain: - """Aggregate states into blocks, returning the lumped chain.""" + """Aggregate states into blocks, returning the lumped chain. + + ``partition=None`` uses the coarsest lumpable partition + (:func:`~sofic.generators.lumping.bisimulation_partition`). + """ from sofic.generators.lumping import lump return lump(self, partition, check=check, labels=labels, rtol=rtol, atol=atol) diff --git a/sofic/generators/mealy.py b/sofic/generators/mealy.py index 2777b56..c5700db 100644 --- a/sofic/generators/mealy.py +++ b/sofic/generators/mealy.py @@ -95,14 +95,18 @@ def is_lumpable(self, partition: PartitionLike, *, rtol: float = 1e-8, atol: flo def lump( self, - partition: PartitionLike, + partition: PartitionLike | None = None, *, check: bool = True, labels: LabelsLike | None = None, rtol: float = 1e-8, atol: float = 1e-10, ) -> MealyHMM: - """Aggregate states into blocks, returning the lumped Mealy HMM.""" + """Aggregate states into blocks, returning the lumped Mealy HMM. + + ``partition=None`` uses the coarsest lumpable partition + (:func:`~sofic.generators.lumping.bisimulation_partition`). + """ from sofic.generators.lumping import lump return lump(self, partition, check=check, labels=labels, rtol=rtol, atol=atol) diff --git a/sofic/generators/measures.py b/sofic/generators/measures.py index 5b83468..396e6ff 100644 --- a/sofic/generators/measures.py +++ b/sofic/generators/measures.py @@ -3,7 +3,7 @@ from __future__ import annotations import warnings -from collections.abc import Iterator +from collections.abc import Iterable, Iterator from typing import Any, Literal, NamedTuple import numpy as np @@ -56,6 +56,21 @@ def dit_state_label(state: Any) -> Any: return state +def dit_state_labels(states: Iterable[Any]) -> dict[Any, Any]: + """Map each state to a dit-safe label that dit can sort against the others. + + Per-state labels from :func:`dit_state_label` are kept when they share a + type. Mixed types (e.g. ``0`` beside an encoded tuple after a state split) + cannot be ordered, so every label becomes its ``repr``, which keeps them + distinct. + """ + states = list(states) + labels = [dit_state_label(state) for state in states] + if len({type(label) for label in labels}) > 1: + labels = [repr(label) for label in labels] + return dict(zip(states, labels, strict=True)) + + def state_distribution(model: StochasticModel) -> Any: """Return the stationary state law as a ``dit.Distribution``. @@ -65,7 +80,8 @@ def state_distribution(model: StochasticModel) -> Any: dit = _require_dit() idx = model.reindex() pi = model.stationary_distribution() - outcomes = [(dit_state_label(idx.state(i)),) for i in range(len(idx))] + labels = dit_state_labels(idx.states) + outcomes = [(labels[idx.state(i)],) for i in range(len(idx))] from sofic.generators.prob import as_prob, has_symbolic, simplify_prob probs = [as_prob(pi[i]) for i in range(len(idx))] @@ -203,9 +219,10 @@ def entropy_rate_hmm(hmm: HiddenMarkovModel) -> Any: # states like ("A", "0", "A")) do not break dit.Distribution. outcomes: list[tuple[Any, Any]] = [] probs: list[Any] = [] + labels = dit_state_labels(idx.states) for state in idx.states: i = idx.index(state) - label = dit_state_label(state) + label = labels[state] for symbol, matrix in joint.items(): row_mass = as_prob(pi[i]) * array_sum(matrix[i]) row_mass = simplify_prob(row_mass) if symbolic or is_symbolic(row_mass) else float(row_mass) diff --git a/sofic/generators/minimal_quasi_realization.py b/sofic/generators/minimal_quasi_realization.py new file mode 100644 index 0000000..9cb4131 --- /dev/null +++ b/sofic/generators/minimal_quasi_realization.py @@ -0,0 +1,237 @@ +"""Minimal linear (quasi-)realizations of finite-state processes. + +A hidden Markov model with symbol-labeled matrices :math:`T^{(x)}`, start vector +:math:`\\pi`, and final vector :math:`\\mathbf{1}` is a *linear representation* of +its word function :math:`f(w) = \\pi T^{(w)} \\mathbf{1}`. Schützenberger's +minimization of weighted automata :cite:`Schutzenberger1961` (equivalently, +Fliess' Hankel-rank theorem) reduces any such representation to one of minimal +dimension: the dimension of the subspace that is both reachable (spanned by +:math:`\\pi T^{(w)}`) and observable (spanned by :math:`T^{(w)} \\mathbf{1}`). That +dimension is the *process rank*, the rank of the Hankel matrix +:math:`H_{u, v} = f(uv)`, and the finite dimension of Upper's generalized-state +space :cite:`Upper1997`. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Callable, Mapping +from typing import Any + +import numpy as np + +from sofic.generators.base import HiddenMarkovModel +from sofic.generators.matrices import emission_tensors +from sofic.generators.prob import _sympy +from sofic.generators.quasi_realization import QuasiRealization + +_DEFAULT_TOL = 1e-9 + + +def minimal_quasi_realization( + model: HiddenMarkovModel | QuasiRealization, *, tol: float = _DEFAULT_TOL +) -> QuasiRealization: + """Return a minimal-dimension quasi-realization of ``model``'s word process. + + The representation :math:`(\\pi, T^{(x)}, \\tau)` of ``model`` (for an HMM the + model start vector, its joint symbol matrices, and :math:`\\tau = \\mathbf{1}`) + is reduced by Schützenberger/Fliess minimization :cite:`Schutzenberger1961`: + bases :math:`R` of the reachable row space :math:`\\{\\pi T^{(w)}\\}` and + :math:`O` of the observable column space :math:`\\{T^{(w)} \\tau\\}` are grown + breadth-first, the finite Hankel core :math:`M = R O` is factored as + :math:`M = C F` with :math:`k = \\operatorname{rank} M`, and the projected maps + + .. math:: D^{(x)} = C^{+} R \\, T^{(x)} \\, O F^{+} + + reproduce every word probability with :math:`k` equal to the process rank + :cite:`Upper1997`. A final change of basis makes ``tau`` the all-ones vector, + so ``pi`` sums to :math:`f(\\varepsilon)` (one for a normalized process). + + Symbolic (sympy) input is reduced exactly and returns object-dtype arrays of + sympy expressions; numeric input uses orthonormal Krylov bases and truncates + the singular values of :math:`M` below ``tol`` times the largest. + + Parameters + ---------- + model + A hidden Markov model (any presentation with ``to_mealy()``) or an + existing :class:`~sofic.generators.quasi_realization.QuasiRealization`. + tol + Relative numeric rank tolerance; ignored for exact input. + """ + pi, matrices, tau = _representation(model) + symbols = list(matrices) + if _is_exact(pi, tau, *matrices.values()): + start, maps, final = _minimize_exact(pi, [matrices[x] for x in symbols], tau) + else: + start, maps, final = _minimize_numeric(pi, [matrices[x] for x in symbols], tau, tol=tol) + return _quasi_realization(start, dict(zip(symbols, maps, strict=True)), final) + + +def process_rank(model: HiddenMarkovModel | QuasiRealization, *, tol: float = _DEFAULT_TOL) -> int: + """Return the process rank: the dimension of :func:`minimal_quasi_realization`. + + This equals the rank of the Hankel matrix :math:`H_{u, v} = P(uv)` of word + probabilities :cite:`Upper1997`, and lower-bounds the state count of every + HMM presentation of the process, including its ε-machine. + """ + return int(minimal_quasi_realization(model, tol=tol).pi.size) + + +def _representation( + model: HiddenMarkovModel | QuasiRealization, +) -> tuple[np.ndarray, dict[Any, np.ndarray], np.ndarray]: + if isinstance(model, QuasiRealization): + return model.pi, dict(model.symbol_maps), model.tau + pi, matrices = emission_tensors(model) + symbolic = pi.dtype == object or any(matrix.dtype == object for matrix in matrices.values()) + tau = np.ones(len(pi), dtype=object if symbolic else float) + if symbolic: + sp = _sympy() + tau.fill(sp.Integer(1)) + return pi, matrices, tau + + +def _is_exact(*arrays: np.ndarray) -> bool: + return any(np.asarray(array).dtype == object for array in arrays) + + +def _minimize_numeric( + pi: np.ndarray, matrices: list[np.ndarray], tau: np.ndarray, *, tol: float +) -> tuple[np.ndarray, list[np.ndarray], np.ndarray]: + pi = np.asarray(pi, dtype=float) + tau = np.asarray(tau, dtype=float) + matrices = [np.asarray(matrix, dtype=float) for matrix in matrices] + n = pi.size + + reachable = _orthonormal_krylov(pi, [lambda v, m=m: v @ m for m in matrices], tol=tol, n=n) + observable = _orthonormal_krylov(tau, [lambda v, m=m: m @ v for m in matrices], tol=tol, n=n).T + core = reachable @ observable + if core.size == 0: + return _empty_numeric(len(matrices)) + u, s, vt = np.linalg.svd(core) + k = int(np.sum(s > tol * s[0])) if s.size and s[0] > 0.0 else 0 + if k == 0: + return _empty_numeric(len(matrices)) + root = np.sqrt(s[:k]) + left = (u[:, :k] / root).T @ reachable + right = observable @ (vt[:k].T / root) + + start = pi @ right + maps = [left @ matrix @ right for matrix in matrices] + final = left @ tau + pivot = int(np.argmax(np.abs(final))) + shift, unshift = _unit_final_basis(final, pivot, identity=np.eye(k)) + return start @ unshift, [shift @ matrix @ unshift for matrix in maps], np.ones(k) + + +def _orthonormal_krylov( + seed: np.ndarray, steps: list[Callable[[np.ndarray], np.ndarray]], *, tol: float, n: int +) -> np.ndarray: + basis: list[np.ndarray] = [] + queue: deque[np.ndarray] = deque([seed]) + while queue and len(basis) < n: + vector = queue.popleft() + scale = float(np.linalg.norm(vector)) + residual = vector + for _ in range(2): + for q in basis: + residual = residual - (q @ residual) * q + norm = float(np.linalg.norm(residual)) + if scale == 0.0 or norm <= tol * max(1.0, scale): + continue + direction = residual / norm + basis.append(direction) + queue.extend(step(direction) for step in steps) + return np.asarray(basis, dtype=float).reshape(len(basis), n) + + +def _empty_numeric(num_symbols: int) -> tuple[np.ndarray, list[np.ndarray], np.ndarray]: + return np.zeros(0), [np.zeros((0, 0)) for _ in range(num_symbols)], np.zeros(0) + + +def _minimize_exact( + pi: np.ndarray, matrices: list[np.ndarray], tau: np.ndarray +) -> tuple[np.ndarray, list[np.ndarray], np.ndarray]: + sp = _sympy() + n = len(pi) + start_row = sp.Matrix([list(pi)]).applyfunc(sp.sympify) + final_col = sp.Matrix(list(tau)).applyfunc(sp.sympify) + maps = [sp.Matrix(np.asarray(matrix, dtype=object).tolist()).applyfunc(sp.sympify) for matrix in matrices] + + reachable = _exact_krylov(start_row, [lambda v, m=m: v * m for m in maps], n=n) + observable = _exact_krylov(final_col.T, [lambda v, m=m: (m * v.T).T for m in maps], n=n).T + core = (reachable * observable).applyfunc(sp.simplify) + if core.rank(simplify=True) == 0: + return _empty_exact(len(maps)) + c, f = core.rank_decomposition(simplify=True) + k = c.shape[1] + left = (c.T * c).inv() * c.T * reachable + right = observable * f.T * (f * f.T).inv() + + start = (start_row * right).applyfunc(sp.simplify) + reduced = [(left * m * right).applyfunc(sp.simplify) for m in maps] + final = (left * final_col).applyfunc(sp.simplify) + pivot = next(i for i in range(k) if sp.simplify(final[i]) != 0) + shift, unshift = _unit_final_basis(final, pivot, identity=sp.eye(k)) + return ( + _to_object_vector(start * unshift), + [_to_object(shift * m * unshift) for m in reduced], + _to_object_vector(sp.ones(k, 1)), + ) + + +def _exact_krylov(seed: Any, steps: list[Callable[[Any], Any]], *, n: int) -> Any: + sp = _sympy() + rows: list[Any] = [] + queue: deque[Any] = deque([seed]) + while queue and len(rows) < n: + row = queue.popleft().applyfunc(sp.simplify) + if sp.Matrix.vstack(*rows, row).rank(simplify=True) <= len(rows): + continue + rows.append(row) + queue.extend(step(row) for step in steps) + return sp.Matrix.vstack(*rows) if rows else sp.zeros(0, n) + + +def _empty_exact(num_symbols: int) -> tuple[np.ndarray, list[np.ndarray], np.ndarray]: + return ( + np.zeros(0, dtype=object), + [np.zeros((0, 0), dtype=object) for _ in range(num_symbols)], + np.zeros(0, dtype=object), + ) + + +def _unit_final_basis(final: Any, pivot: int, *, identity: Any) -> tuple[Any, Any]: + """Return ``(S, S^{-1})`` with ``S @ final`` the all-ones vector. + + ``S = I + (1 - final) e_pivot^T / final[pivot]`` has determinant + ``1 / final[pivot]`` and inverse ``I - (1 - final) e_pivot^T``. + """ + k = identity.shape[0] + deficit = [1 - final[i] for i in range(k)] + shift = identity.copy() + unshift = identity.copy() + for i in range(k): + shift[i, pivot] += deficit[i] / final[pivot] + unshift[i, pivot] -= deficit[i] + return shift, unshift + + +def _to_object(matrix: Any) -> np.ndarray: + sp = _sympy() + return np.asarray(matrix.applyfunc(sp.simplify).tolist(), dtype=object).reshape(matrix.shape) + + +def _to_object_vector(matrix: Any) -> np.ndarray: + return _to_object(matrix).ravel() + + +def _quasi_realization(pi: np.ndarray, symbol_maps: Mapping[Any, np.ndarray], tau: np.ndarray) -> QuasiRealization: + if pi.dtype != object: + return QuasiRealization(pi=pi, tau=tau, symbol_maps=dict(symbol_maps)) + realization = QuasiRealization(pi=np.zeros(0), tau=np.zeros(0), symbol_maps={}) + realization.pi = pi + realization.tau = tau + realization.symbol_maps = dict(symbol_maps) + return realization diff --git a/sofic/generators/moore.py b/sofic/generators/moore.py index eeb5af4..007f4a3 100644 --- a/sofic/generators/moore.py +++ b/sofic/generators/moore.py @@ -44,14 +44,18 @@ def is_lumpable(self, partition: PartitionLike, *, rtol: float = 1e-8, atol: flo def lump( self, - partition: PartitionLike, + partition: PartitionLike | None = None, *, check: bool = True, labels: LabelsLike | None = None, rtol: float = 1e-8, atol: float = 1e-10, ) -> MooreHMM: - """Aggregate states into blocks, returning the lumped Moore HMM.""" + """Aggregate states into blocks, returning the lumped Moore HMM. + + ``partition=None`` uses the coarsest lumpable partition + (:func:`~sofic.generators.lumping.bisimulation_partition`). + """ from sofic.generators.lumping import lump return lump(self, partition, check=check, labels=labels, rtol=rtol, atol=atol) diff --git a/sofic/generators/omega_probability.py b/sofic/generators/omega_probability.py new file mode 100644 index 0000000..7c12c33 --- /dev/null +++ b/sofic/generators/omega_probability.py @@ -0,0 +1,156 @@ +"""Probabilities of regular and ω-regular properties of HMM output. + +Both quantities are computed on the product of the HMM with a deterministic +automaton: the automaton state after reading ``X_{0:t}`` is a function of the +output, so the pair (HMM state, automaton state) is again a finite Markov chain +:cite:`BaierKatoen2008`. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Hashable +from typing import TYPE_CHECKING, Any + +import networkx as nx +import numpy as np + +from sofic.exceptions import NonDeterministicError +from sofic.graph import ATTR_SYMBOL + +if TYPE_CHECKING: + from sofic.automata.base import LabeledAutomaton + from sofic.automata.buchi import BuchiAutomaton + from sofic.automata.dfa import DFA + from sofic.generators.base import HiddenMarkovModel + from sofic.generators.matrices import StartSpec + +_TOL = 1e-12 +_DEAD = object() + + +def _tensors(hmm: HiddenMarkovModel, start: StartSpec) -> tuple[np.ndarray, dict[Any, np.ndarray]]: + from sofic.generators.matrices import start_vector, symbol_matrices + + mealy = hmm.to_mealy() + pi = np.asarray(start_vector(mealy, start, policy="stationary"), dtype=float) + matrices = {symbol: np.asarray(matrix, dtype=float) for symbol, matrix in symbol_matrices(mealy).items()} + return pi, matrices + + +def _transition_table(aut: LabeledAutomaton) -> tuple[Hashable, dict[tuple[Hashable, Any], Hashable]]: + if not aut.is_deterministic(): + raise NonDeterministicError(f"{type(aut).__name__} must be deterministic") + (initial,) = aut.initial_states + return initial, {(t.source, t.data[ATTR_SYMBOL]): t.target for t in aut.transitions()} + + +def omega_probability(hmm: HiddenMarkovModel, buchi: BuchiAutomaton, *, start: StartSpec = None) -> float: + """Return the probability that the output ``X_0 X_1 ...`` of ``hmm`` is accepted by ``buchi``. + + ``buchi`` must be a deterministic Büchi automaton (one initial state, no + ε-moves, at most one successor per state and symbol); a missing transition + rejects. The HMM starts in its stationary law, or in ``start`` (a state, a + state-to-mass mapping, or a vector in ``hmm.to_mealy().reindex()`` order). + + The product chain on (HMM state, automaton state), with an absorbing + rejecting state for missing transitions, almost surely enters a bottom + strongly connected component and then visits each of its states infinitely + often. The run is therefore accepted iff that component contains an + accepting automaton state, and the probability is the mass absorbed into + such components :cite:`BaierKatoen2008`. + + Deterministic Büchi automata cannot express every ω-regular property + (``finitely many 1s`` is the standard counterexample); its probability is + one minus that of the DBA ``infinitely many 1s``. + + Raises + ------ + NonDeterministicError + If ``buchi`` is not deterministic. + """ + pi, matrices = _tensors(hmm, start) + initial, delta = _transition_table(buchi) + + nodes: dict[Any, int] = {_DEAD: 0} + frontier: deque[tuple[int, Hashable]] = deque() + for i in np.flatnonzero(pi > _TOL): + nodes[(int(i), initial)] = len(nodes) + frontier.append((int(i), initial)) + edges: list[tuple[int, int, float]] = [(0, 0, 1.0)] + while frontier: + i, q = frontier.popleft() + for symbol, matrix in matrices.items(): + target = delta.get((q, symbol)) + for j in np.flatnonzero(matrix[i] > _TOL): + node = _DEAD if target is None else (int(j), target) + if node not in nodes: + nodes[node] = len(nodes) + frontier.append(node) + edges.append((nodes[(i, q)], nodes[node], float(matrix[i, j]))) + + size = len(nodes) + chain = np.zeros((size, size)) + for source, target, prob in edges: + chain[source, target] += prob + graph = nx.DiGraph() + graph.add_nodes_from(range(size)) + graph.add_edges_from((source, target) for source, target, _ in edges) + condensation = nx.condensation(graph) + accepting = {index for node, index in nodes.items() if node is not _DEAD and node[1] in buchi.accepting_states} + good: set[int] = set() + for component in condensation.nodes: + members = condensation.nodes[component]["members"] + if condensation.out_degree(component) == 0 and members & accepting: + good |= members + if not good: + return 0.0 + + reaching = set(good) + for node in good: + reaching |= nx.ancestors(graph, node) + transient = sorted(reaching - good) + absorbed = np.zeros(size) + absorbed[sorted(good)] = 1.0 + if transient: + sub = chain[np.ix_(transient, transient)] + into = chain[np.ix_(transient, sorted(good))].sum(axis=1) + absorbed[transient] = np.linalg.solve(np.eye(len(transient)) - sub, into) + + total = sum(pi[i] * absorbed[nodes[(int(i), initial)]] for i in np.flatnonzero(pi > _TOL)) + return float(min(max(total, 0.0), 1.0)) + + +def regular_language_probability(hmm: HiddenMarkovModel, dfa: DFA, n: int, *, start: StartSpec = None) -> float: + """Return ``P(X_{0:n} ∈ L(dfa))`` exactly. + + Dynamic programming over the product of HMM states and DFA states: the + forward vector of each DFA state is pushed through ``T^(x)`` along + ``δ(q, x)`` for ``n`` steps, and the mass in accepting DFA states is summed + :cite:`HopcroftUllman1979,Rabiner1989`. Cost is ``O(n |Q| |A| |S|^2)``. + ``start`` is as in :func:`omega_probability`, defaulting to the stationary + law. A missing DFA transition rejects. + + Raises + ------ + NonDeterministicError + If ``dfa`` is not deterministic. + ValueError + If ``n`` is negative. + """ + if n < 0: + raise ValueError("n must be nonnegative") + pi, matrices = _tensors(hmm, start) + initial, delta = _transition_table(dfa) + alpha: dict[Hashable, np.ndarray] = {initial: pi} + for _ in range(n): + following: dict[Hashable, np.ndarray] = {} + for q, vector in alpha.items(): + for symbol, matrix in matrices.items(): + target = delta.get((q, symbol)) + if target is None: + continue + pushed = vector @ matrix + following[target] = following[target] + pushed if target in following else pushed + alpha = following + return float(sum(vector.sum() for q, vector in alpha.items() if q in dfa.accepting_states)) diff --git a/sofic/generators/process_equivalence.py b/sofic/generators/process_equivalence.py index bd5881f..d207c94 100644 --- a/sofic/generators/process_equivalence.py +++ b/sofic/generators/process_equivalence.py @@ -78,7 +78,7 @@ def from_hmm( mealy = hmm.to_mealy() return cls( alphabet=tuple(sorted(mealy.observation_alphabet, key=repr)), - matrices=symbol_matrices(mealy), + matrices={symbol: np.asarray(matrix, dtype=float) for symbol, matrix in symbol_matrices(mealy).items()}, start=np.asarray(start_vector(mealy, start), dtype=float), ) diff --git a/sofic/generators/relative_entropy_rate.py b/sofic/generators/relative_entropy_rate.py index b034fac..7bb2a3a 100644 --- a/sofic/generators/relative_entropy_rate.py +++ b/sofic/generators/relative_entropy_rate.py @@ -22,10 +22,9 @@ import networkx as nx import numpy as np -from sofic.graph import ATTR_SYMBOL +from sofic.generators.support import support_includes, support_nfa if TYPE_CHECKING: - from sofic.automata.nfa import NFA from sofic.generators.base import HiddenMarkovModel _TOL = 1e-12 @@ -71,41 +70,8 @@ def _aligned(p: _Presentation, q: _Presentation) -> tuple[list[Any], dict[Any, n return symbols, {symbol: q.matrices.get(symbol, zeros) for symbol in symbols} -def _support_nfa(model: HiddenMarkovModel) -> NFA: - """Return an NFA accepting exactly the words of positive stationary probability. - - States are the generator states of positive stationary probability, all of - them initial and accepting; edges are the positive-probability transitions. - The stationary support is closed under such transitions, so a word is - accepted iff some stationary-positive state can emit it. - """ - from sofic.automata.nfa import NFA - - presentation = _presentation(model) - live = [i for i, mass in enumerate(presentation.pi) if mass > _TOL] - states = frozenset(presentation.states[i] for i in live) - nfa = NFA( - input_alphabet=frozenset(presentation.matrices), - initial_states=states, - accepting_states=states, - ) - for state in states: - nfa.graph.add_state(state) - for symbol, matrix in presentation.matrices.items(): - for i in live: - for j in np.flatnonzero(matrix[i] > _TOL): - nfa.graph.add_transition(presentation.states[i], presentation.states[j], **{ATTR_SYMBOL: symbol}) - return nfa - - -def _support_includes(p: HiddenMarkovModel, q: HiddenMarkovModel) -> bool: - """Return whether every word with positive ``P`` probability has positive ``Q`` probability. - - This is absolute continuity ``P_{0:n} << Q_{0:n}`` for every ``n``, decided - by antichain inclusion of the stationary support automata - :cite:`DeWulf2006`. - """ - return _support_nfa(q).includes(_support_nfa(p)) +_support_nfa = support_nfa +_support_includes = support_includes def _unifilar_successors(matrices: dict[Any, np.ndarray]) -> dict[Any, dict[int, int]] | None: @@ -126,8 +92,8 @@ def relative_entropy_rate(p: HiddenMarkovModel, q: HiddenMarkovModel) -> float: """Return the relative entropy rate ``D(P || Q)`` in bits per symbol, exactly. The result is :data:`math.inf` when some word of positive ``P`` probability - has ``Q`` probability zero, decided by inclusion of the support automata - (:meth:`~sofic.automata.base.LabeledAutomaton.includes`). + has ``Q`` probability zero, decided by + :func:`~sofic.generators.support.support_includes`. Otherwise ``q`` must be unifilar. Writing ``h_μ(P)`` for the entropy rate of ``p``, @@ -159,7 +125,7 @@ def relative_entropy_rate(p: HiddenMarkovModel, q: HiddenMarkovModel) -> float: If ``q`` is not unifilar, or is not exactly synchronized by ``P``-typical pasts. Use :func:`relative_entropy_rate_bounds` instead. """ - if not _support_includes(p, q): + if not support_includes(p, q): return math.inf pp, qp = _presentation(p), _presentation(q) symbols, q_matrices = _aligned(pp, qp) @@ -280,7 +246,7 @@ def relative_entropy_rate_bounds( if block_length < 1: raise ValueError("block_length must be positive") n = block_length - if not _support_includes(p, q): + if not support_includes(p, q): return RelativeEntropyRateBounds(math.inf, math.inf, n, math.inf, math.inf) pp, qp = _presentation(p), _presentation(q) symbols, q_matrices = _aligned(pp, qp) diff --git a/sofic/generators/renyi.py b/sofic/generators/renyi.py new file mode 100644 index 0000000..9f084e4 --- /dev/null +++ b/sofic/generators/renyi.py @@ -0,0 +1,246 @@ +"""Rényi entropy rates and the thermodynamic spectrum of a stationary process. + +For a stationary process with length-``n`` block law :math:`P`, the Rényi +entropy rate of order :math:`\\alpha \\ge 0` and the pressure are + +.. math:: + + h_\\alpha = \\lim_{n \\to \\infty} \\tfrac{1}{n} H_\\alpha[X_{0:n}], \\qquad + \\mathcal{P}(\\beta) = \\lim_{n \\to \\infty} \\tfrac{1}{n} \\log_2 \\sum_{w} P(w)^\\beta, + +so that :math:`h_\\alpha = \\mathcal{P}(\\alpha) / (1 - \\alpha)`. On a unifilar +presentation with labeled matrices :math:`T^{(x)}`, restricted to the +stationary support, + +.. math:: \\mathcal{P}(\\beta) = \\log_2 \\rho\\Big(\\sum_x \\big(T^{(x)}\\big)^{\\circ \\beta}\\Big) + +for :math:`\\beta \\ge 0`, where :math:`\\circ \\beta` is the entrywise power of +the positive entries and :math:`\\rho` the spectral radius. For Markov chains +this is :cite:`Rached2001`. For a unifilar presentation each start state ``s`` +emits ``w`` along at most one path, of probability :math:`p(w \\mid s)`, and +:math:`P(w) = \\sum_s \\pi_s p(w \\mid s)`. For :math:`\\beta \\ge 0` and +:math:`N` support states, + +.. math:: + + N^{-1} \\textstyle\\sum_s (\\pi_s p(w \\mid s))^\\beta \\le P(w)^\\beta + \\le N^\\beta \\textstyle\\sum_s (\\pi_s p(w \\mid s))^\\beta, + +and :math:`\\sum_w \\sum_s (\\pi_s p(w \\mid s))^\\beta = (\\pi^{\\circ\\beta})^\\top A_\\beta^n \\mathbf{1}` +with :math:`A_\\beta = \\sum_x (T^{(x)})^{\\circ \\beta}`. Both vectors are +positive on the support, so the growth rate is :math:`\\rho(A_\\beta)`. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from typing import TYPE_CHECKING + +import numpy as np +from scipy.optimize import minimize_scalar + +if TYPE_CHECKING: + from sofic.generators.base import HiddenMarkovModel + from sofic.generators.mealy import MealyHMM + +__all__ = ["pressure", "rate_function", "renyi_entropy_rate"] + +_TOL = 1e-12 +_BETA_MAX = 64.0 + + +def _unifilar_mealy(model: HiddenMarkovModel, max_states: int) -> MealyHMM: + from sofic.generators.epsilon_machine import EpsilonMachine + + mealy = model.to_mealy() + if not mealy.is_unifilar(): + mealy = EpsilonMachine.from_hmm(mealy, max_states=max_states) + return mealy + + +def _stationary_support(pi: np.ndarray, matrices: Sequence[np.ndarray]) -> np.ndarray: + """Return the mask of states reachable from those of positive stationary mass. + + Closing under transitions keeps recurrent states whose computed stationary + mass underflows the threshold. + """ + adjacency = sum((np.asarray(m, dtype=float) > 0.0) for m in matrices) + keep = pi > _TOL + frontier = keep.copy() + while frontier.any(): + reached = np.asarray(adjacency[frontier].sum(axis=0) > 0) & ~keep + keep |= reached + frontier = reached + return keep + + +class _Spectrum: + """Labeled matrices of a unifilar presentation on its stationary support.""" + + def __init__(self, model: HiddenMarkovModel, *, max_states: int) -> None: + from sofic.generators.matrices import symbol_matrices + + self.mealy = _unifilar_mealy(model, max_states) + pi = np.asarray(self.mealy.stationary_distribution(), dtype=float) + raw = list(symbol_matrices(self.mealy).values()) + keep = _stationary_support(pi, raw) + matrices = [np.asarray(m, dtype=float)[np.ix_(keep, keep)] for m in raw] + self.matrices = [m for m in matrices if np.any(m > 0.0)] + self.size = int(keep.sum()) + self.mu = self._max_cycle_mean() + + def _max_cycle_mean(self) -> float: + """Karp's maximum mean cycle weight of the ``log2`` edge probabilities.""" + n = self.size + weights = np.full((n, n), -np.inf) + for m in self.matrices: + with np.errstate(divide="ignore"): + weights = np.maximum(weights, np.where(m > 0.0, np.log2(np.where(m > 0.0, m, 1.0)), -np.inf)) + walks = np.zeros((n + 1, n)) + for k in range(1, n + 1): + walks[k] = np.max(walks[k - 1][:, None] + weights, axis=0) + best = -np.inf + for v in range(n): + if not np.isfinite(walks[n, v]): + continue + finite = np.isfinite(walks[:n, v]) + ks = np.arange(n)[finite] + best = max(best, float(np.min((walks[n, v] - walks[ks, v]) / (n - ks)))) + return best + + def pressure(self, beta: float) -> float: + scale = beta * self.mu + total = np.zeros((self.size, self.size)) + for m in self.matrices: + positive = m > 0.0 + logs = np.log2(np.where(positive, m, 1.0)) + total += np.where(positive, np.exp2(beta * logs - scale), 0.0) + return scale + float(np.log2(np.max(np.abs(np.linalg.eigvals(total))))) + + def pressure_slope_at_zero(self, step: float = 1e-5) -> float: + return (-3.0 * self.pressure(0.0) + 4.0 * self.pressure(step) - self.pressure(2.0 * step)) / (2.0 * step) + + +def _check_order(value: float, name: str) -> float: + value = float(value) + if np.isnan(value) or value < 0.0: + raise ValueError(f"{name} must be non-negative, got {value}") + return value + + +def renyi_entropy_rate(model: HiddenMarkovModel, alpha: float, *, max_states: int = 10_000) -> float: + """Return the Rényi entropy rate :math:`h_\\alpha` in bits per symbol. + + :math:`h_\\alpha = \\lim_n \\tfrac{1}{n} H_\\alpha[X_{0:n}] + = \\tfrac{1}{1 - \\alpha} \\log_2 \\rho\\big(\\sum_x (T^{(x)})^{\\circ \\alpha}\\big)` + on a unifilar presentation (:cite:`Rached2001` for Markov chains; see the + module docstring for unifilar HMMs). ``alpha`` ranges over + :math:`[0, \\infty]`: + + * ``alpha = 0`` is the topological entropy of the process support, + :math:`\\log_2 \\rho` of the support adjacency matrix; + * ``alpha = 1`` is the Shannon entropy rate :math:`h_\\mu` + (:meth:`~sofic.generators.base.HiddenMarkovModel.entropy_rate`); + * ``alpha = inf`` is the min-entropy rate + :math:`-\\lim_n \\tfrac{1}{n} \\log_2 \\max_w P(w)`, minus the maximum mean + :math:`\\log_2` edge probability over cycles (Karp's algorithm). + + Non-unifilar input is first converted with + :meth:`EpsilonMachine.from_hmm(model, max_states=max_states) + `, which raises + :class:`~sofic.exceptions.MixedStateExplosionError` when the mixed states + do not close. + + Examples + -------- + >>> from sofic.examples import bernoulli + >>> round(renyi_entropy_rate(bernoulli(0.25), 2.0), 6) + 0.678072 + """ + alpha = _check_order(alpha, "alpha") + spectrum = _Spectrum(model, max_states=max_states) + if np.isinf(alpha): + return -spectrum.mu + if alpha == 1.0: + return float(spectrum.mealy.entropy_rate()) + return spectrum.pressure(alpha) / (1.0 - alpha) + + +def pressure(model: HiddenMarkovModel, beta: float, *, max_states: int = 10_000) -> float: + """Return the pressure :math:`\\mathcal{P}(\\beta) = \\lim_n \\tfrac{1}{n} \\log_2 \\sum_w P(w)^\\beta`. + + Computed as :math:`\\log_2 \\rho\\big(\\sum_x (T^{(x)})^{\\circ \\beta}\\big)` on a + unifilar presentation for finite :math:`\\beta \\ge 0` (see the module + docstring). :math:`\\mathcal{P}` is convex and non-increasing, with + :math:`\\mathcal{P}(0) = h_0`, :math:`\\mathcal{P}(1) = 0` and + :math:`\\mathcal{P}(\\beta) = (1 - \\beta) h_\\beta`. Negative :math:`\\beta` + is rejected: there the matrix formula weights each word by its least + likely start state and need not equal the pressure. Non-unifilar input is + handled as in :func:`renyi_entropy_rate`. + + Examples + -------- + >>> from sofic.examples import golden_mean + >>> abs(pressure(golden_mean(0.5), 1.0)) < 1e-12 + True + """ + beta = _check_order(beta, "beta") + if np.isinf(beta): + raise ValueError("beta must be finite") + return _Spectrum(model, max_states=max_states).pressure(beta) + + +def rate_function( + model: HiddenMarkovModel, + values: float | Sequence[float] | np.ndarray, + *, + beta_max: float = _BETA_MAX, + max_states: int = 10_000, +) -> np.ndarray: + """Return the large-deviation rate function of :math:`-\\tfrac{1}{n} \\log_2 P(X_{0:n})`. + + The scaled cumulant generating function of + :math:`U_n = -\\tfrac{1}{n} \\log_2 P(X_{0:n})` in base 2 is + :math:`\\lim_n \\tfrac{1}{n} \\log_2 \\mathbb{E}[2^{n t U_n}] = \\mathcal{P}(1 - t)`, + so by the Gärtner-Ellis theorem :cite:`DemboZeitouni1998` (§2.3) the rate + function, in bits per symbol, is the Legendre transform + + .. math:: I(u) = \\sup_{\\beta} \\big[(1 - \\beta) u - \\mathcal{P}(\\beta)\\big]. + + Only :math:`\\beta \\in [0, \\infty)` is available (see :func:`pressure`), so + ``I`` is computed on :math:`[h_\\infty, u_0]` with + :math:`u_0 = -\\mathcal{P}'(0)`, the typical value of ``U_n`` under the + measure of maximal entropy on the support. The supremum is taken over + :math:`[0, \\beta_{\\max}]` by a grid search refined with + :func:`scipy.optimize.minimize_scalar`; near :math:`h_\\infty` the optimum + escapes to :math:`\\beta \\to \\infty`, so values there are lower bounds that + tighten with ``beta_max``. ``I(h_\\mu) = 0`` and ``I \\ge 0``. Values below + :math:`h_\\infty` give ``inf`` (no word is that likely) and values above + :math:`u_0` give ``nan``. Returns an array shaped like ``values``. + """ + beta_max = _check_order(beta_max, "beta_max") + spectrum = _Spectrum(model, max_states=max_states) + h_inf = -spectrum.mu + u_0 = -spectrum.pressure_slope_at_zero() + grid = np.concatenate([[0.0], np.geomspace(1e-3, beta_max, 96)]) if beta_max > 0 else np.zeros(1) + pressures = np.array([spectrum.pressure(beta) for beta in grid]) + + def evaluate(u: float) -> float: + if u < h_inf - 1e-9: + return np.inf + if u > u_0 + 1e-6: + return np.nan + objective = (1.0 - grid) * u - pressures + best = int(np.argmax(objective)) + if 0 < best < grid.size - 1: + result = minimize_scalar( + lambda beta: -((1.0 - beta) * u - spectrum.pressure(beta)), + bounds=(grid[best - 1], grid[best + 1]), + method="bounded", + options={"xatol": 1e-10}, + ) + return max(float(-result.fun), float(objective[best]), 0.0) + return max(float(objective[best]), 0.0) + + flat = np.asarray(values, dtype=float) + return np.vectorize(evaluate, otypes=[float])(flat) diff --git a/sofic/generators/state_splitting.py b/sofic/generators/state_splitting.py new file mode 100644 index 0000000..d36685c --- /dev/null +++ b/sofic/generators/state_splitting.py @@ -0,0 +1,337 @@ +"""Process-preserving state splitting and amalgamation of hidden Markov models. + +Lind & Marcus state splitting (:cite:`LindMarcus1995`, §2.4; implemented for +edge shifts in :mod:`sofic.shifts.state_splitting`) is lifted to Mealy HMMs so +that the generated *process* is unchanged, following the equivalence moves of +:cite:`Upper1997`. Edges are named ``(source, target, key)`` as in +:class:`~sofic.graph.Transition`; :func:`out_edges` and :func:`in_edges` list +them. Splitting state ``s`` into ``m`` copies names them ``(s, 0), ..., +(s, m - 1)``. + +Write ``p(e)`` for the joint probability ``P(target, symbol | source)`` of edge +``e`` and ``pi`` for the initial distribution. + +*Out-splitting* partitions the out-edges of ``s`` into ``E_0, ..., E_{m-1}`` +with masses ``w_i = sum_{e in E_i} p(e)``. Copy ``(s, i)`` keeps the edges of +``E_i`` renormalized to ``p(e) / w_i``; every in-edge ``e`` of ``s`` becomes +``m`` edges into the copies with probabilities ``p(e) w_i``; and +``pi(s, i) = pi(s) w_i``. By induction on ``t`` the forward vector of the split +model is ``alpha_t(s, i) = w_i alpha_t(s)``, so every word has the same +probability. If ``pi`` is stationary, so is the split ``pi``. + +*In-splitting* partitions the in-edges of ``s`` into ``F_0, ..., F_{m-1}``. +Edge ``e in F_i`` enters copy ``(s, i)`` only, and every copy keeps all +out-edges of ``s`` with their probabilities. Copies have identical futures, so +the backward vectors agree, ``beta_t(s, i) = beta_t(s)``, and the initial mass +``pi(s)`` may be shared out arbitrarily; it is split in proportion to the +stationary inflow ``sum_{e in F_i} mu(source(e)) p(e)`` (uniformly if ``s`` +has no stationary inflow), which keeps a stationary ``pi`` stationary. + +:func:`amalgamate` inverts both: in-split copies have identical futures and +merge by Kemeny-Snell strong lumping (:func:`~sofic.generators.lumping.lump`, +:cite:`KemenySnell1976`), while out-split copies receive proportional inflow +and merge into the inflow-weighted mixture of their out-edges. +""" + +from __future__ import annotations + +from collections import defaultdict +from collections.abc import Collection, Hashable, Mapping, Sequence +from typing import Any, Literal + +import numpy as np + +from sofic.exceptions import LumpabilityError +from sofic.generators.base import HiddenMarkovModel +from sofic.generators.lumping import ( + LabelsLike, + PartitionLike, + _block_of, + _resolve_labels, + is_lumpable, + lump, + normalize_partition, +) +from sofic.generators.mealy import MealyHMM +from sofic.graph import ATTR_EMISSION, ATTR_PROB, Transition, TransitionGraph + +Edge = tuple[Hashable, Hashable, int] + + +def out_edges(hmm: HiddenMarkovModel, state: Hashable) -> list[Edge]: + """Return the ``(source, target, key)`` edges leaving ``state`` :cite:`LindMarcus1995`.""" + return [(t.source, t.target, t.key) for t in hmm.to_mealy().transitions() if t.source == state] + + +def in_edges(hmm: HiddenMarkovModel, state: Hashable) -> list[Edge]: + """Return the ``(source, target, key)`` edges entering ``state`` :cite:`LindMarcus1995`.""" + return [(t.source, t.target, t.key) for t in hmm.to_mealy().transitions() if t.target == state] + + +def split_state( + hmm: HiddenMarkovModel, + state: Hashable, + partition: Sequence[Collection[Edge]], + *, + kind: Literal["out", "in"] = "out", +) -> MealyHMM: + """Split ``state`` into copies without changing the generated process. + + The process-preserving lift of the out-/in-splittings of + :cite:`LindMarcus1995` (Definitions 2.4.3 and 2.4.7) described in the + module docstring; such splittings are among the HMM equivalence moves of + :cite:`Upper1997`. Probabilities are combined with exact arithmetic, so + sympy-valued models stay exact. + + Parameters + ---------- + hmm + The generator; non-Mealy presentations are converted with ``to_mealy``. + state + The state to split into ``(state, 0), ..., (state, m - 1)``. + partition + ``m`` nonempty, disjoint parts covering :func:`out_edges` (``kind="out"``) + or :func:`in_edges` (``kind="in"``) of ``state``. + kind + ``"out"`` partitions out-edges (copies receive proportionally split + inflow); ``"in"`` partitions in-edges (copies share every out-edge). + + Returns + ------- + MealyHMM + A presentation of the same process. In-splitting preserves + unifilarity; out-splitting generally does not, since one in-edge + becomes ``m`` parallel edges with the same symbol. + + Raises + ------ + ValueError + If ``state`` is unknown, the parts are not a partition of the relevant + edges, ``kind`` is invalid, or a copy name collides with a state. + """ + if kind not in ("out", "in"): + raise ValueError(f"kind must be 'out' or 'in', not {kind!r}") + gen = hmm.to_mealy() + if not gen.graph.has_state(state): + raise ValueError(f"unknown state {state!r}") + edges = out_edges(gen, state) if kind == "out" else in_edges(gen, state) + part_of = _part_index(state, edges, partition) + copies = [(state, index) for index in range(len(partition))] + if any(gen.graph.has_state(copy) for copy in copies): + raise ValueError("split state names collide with existing states") + + if kind == "out": + shares: list[Any] = [0] * len(copies) + for t in gen.transitions(): + if t.source == state: + shares[part_of[(t.source, t.target, t.key)]] += t.data[ATTR_PROB] + else: + shares = _stationary_inflow(gen, state, part_of, len(copies)) + + def source_copies(t: Transition) -> list[tuple[Hashable, Any]]: + if t.source != state: + return [(t.source, 1)] + if kind == "out": + index = part_of[(t.source, t.target, t.key)] + return [(copies[index], 1 / shares[index])] + return [(copy, 1) for copy in copies] + + def target_copies(t: Transition) -> list[tuple[Hashable, Any]]: + if t.target != state: + return [(t.target, 1)] + if kind == "out": + return list(zip(copies, shares, strict=True)) + return [(copies[part_of[(t.source, t.target, t.key)]], 1)] + + result = MealyHMM(observation_alphabet=gen.observation_alphabet) + for name in gen.states(): + for copy in copies if name == state else [name]: + result.graph.add_state(copy) + for t in gen.transitions(): + targets = target_copies(t) + for source, source_factor in source_copies(t): + for target, target_factor in targets: + prob = t.data[ATTR_PROB] * source_factor * target_factor + result.add_transition(source, target, t.data.get(ATTR_EMISSION), prob) + for name, mass in gen.initial_distribution.items(): + if name != state: + result.initial_distribution[name] = mass + continue + for copy, share in zip(copies, shares, strict=True): + result.initial_distribution[copy] = mass * share + return result + + +def amalgamate( + hmm: HiddenMarkovModel, + partition: PartitionLike, + *, + labels: LabelsLike | None = None, + rtol: float = 1e-8, + atol: float = 1e-10, +) -> MealyHMM: + """Merge the blocks of ``partition`` without changing the generated process. + + Inverts :func:`split_state` (the amalgamations of :cite:`LindMarcus1995`, + Definition 2.4.9, lifted to processes as in :cite:`Upper1997`). Two exact + merges are recognized: + + * *same futures* -- the partition is strongly lumpable + (:func:`~sofic.generators.lumping.is_lumpable`), as for in-split copies, + and is delegated to :func:`~sofic.generators.lumping.lump` + :cite:`KemenySnell1976`; + * *proportional inflow* -- for each block ``B`` there is a weight vector + ``w`` on ``B`` such that, for every source state and symbol, the edge + mass into ``B`` and the initial mass on ``B`` are multiples of ``w``, as + for out-split copies. The forward vector on ``B`` is then always a + multiple of ``w``, so ``B`` merges exactly into one state whose + out-edges are the ``w``-mixture of its members' out-edges. A block with + neither inflow nor initial mass is never visited and is merged with + uniform ``w``; the process is unchanged, but the out-split weights of + such a block are unrecoverable, so the original rows need not return. + + The whole partition is tried with each rule; failing that, blocks are + merged one at a time with whichever rule applies. + + Parameters + ---------- + hmm + The generator; non-Mealy presentations are converted with ``to_mealy``. + partition + Blocks (iterable of iterables) or a state-to-block mapping covering + every state, as in :func:`~sofic.generators.lumping.normalize_partition`. + labels + Merged-state labels, as in :func:`~sofic.generators.lumping.lump`. + rtol, atol + Tolerances forwarded to :func:`numpy.isclose`. + + Returns + ------- + MealyHMM + A presentation of the same process with one state per block. + + Raises + ------ + LumpabilityError + If some block can be merged by neither rule. + """ + gen = hmm.to_mealy() + blocks = normalize_partition(gen, partition) + resolved = _resolve_labels(blocks, labels) + if is_lumpable(gen, blocks, rtol=rtol, atol=atol): + return lump(gen, blocks, labels=labels) + weights = _inflow_weights(gen, blocks, rtol=rtol, atol=atol) + if weights is not None: + return _merge_by_inflow(gen, blocks, resolved, weights) + + result = gen + for block, label in zip(blocks, resolved, strict=True): + if len(block) == 1 and next(iter(block)) == label: + continue + step = [block, *({state} for state in result.states() if state not in block)] + step_labels = {block: label} + if is_lumpable(result, step, rtol=rtol, atol=atol): + result = lump(result, step, labels=step_labels) + continue + step_blocks = normalize_partition(result, step) + step_weights = _inflow_weights(result, step_blocks, rtol=rtol, atol=atol) + if step_weights is None: + raise LumpabilityError( + f"block {sorted(map(str, block))} has neither identical futures nor proportional inflow" + ) + result = _merge_by_inflow(result, step_blocks, _resolve_labels(step_blocks, step_labels), step_weights) + return result + + +def _part_index(state: Hashable, edges: list[Edge], partition: Sequence[Collection[Edge]]) -> dict[Edge, int]: + expected = set(edges) + part_of: dict[Edge, int] = {} + for index, part in enumerate(partition): + part = set(part) + if not part: + raise ValueError(f"partition of {state!r} has an empty part") + if not part <= expected or part & set(part_of): + raise ValueError(f"parts for {state!r} must be disjoint edges of that state") + part_of.update(dict.fromkeys(part, index)) + if not partition or set(part_of) != expected: + raise ValueError(f"parts for {state!r} do not cover all of its edges") + return part_of + + +def _stationary_inflow(gen: MealyHMM, state: Hashable, part_of: Mapping[Edge, int], m: int) -> list[Any]: + mu = dict(zip(gen.states(), gen.stationary_distribution().tolist(), strict=True)) + inflow: list[Any] = [0] * m + for t in gen.transitions(): + if t.target == state: + inflow[part_of[(t.source, t.target, t.key)]] += mu[t.source] * t.data[ATTR_PROB] + total = sum(inflow) + if not np.isclose(float(total), 0.0): + return [mass / total for mass in inflow] + return [1 / m] * m + + +def _inflow_weights( + gen: MealyHMM, blocks: list[frozenset[Hashable]], *, rtol: float, atol: float +) -> list[dict[Hashable, float]] | None: + """Return per-block inflow weights ``w`` when every block has proportional inflow.""" + block_of = _block_of(blocks) + columns: dict[tuple[int, Hashable, Any], dict[Hashable, float]] = defaultdict(lambda: defaultdict(float)) + for t in gen.transitions(): + key = (block_of[t.target], t.source, t.data.get(ATTR_EMISSION)) + columns[key][t.target] += float(t.data[ATTR_PROB]) + vectors: list[list[dict[Hashable, float]]] = [[] for _ in blocks] + for (index, _source, _symbol), column in columns.items(): + vectors[index].append(column) + initial = defaultdict(dict) + for state, mass in gen.initial_distribution.items(): + initial[block_of[state]][state] = float(mass) + for index, column in initial.items(): + vectors[index].append(column) + + weights: list[dict[Hashable, float]] = [] + for index, block in enumerate(blocks): + reference: np.ndarray | None = None + order = sorted(block, key=repr) + for column in vectors[index]: + vector = np.array([column.get(state, 0.0) for state in order]) + total = vector.sum() + if total <= atol: + continue + vector = vector / total + if reference is None: + reference = vector + elif not np.allclose(vector, reference, rtol=rtol, atol=atol): + return None + if reference is None: + reference = np.full(len(order), 1.0 / len(order)) + weights.append(dict(zip(order, reference.tolist(), strict=True))) + return weights + + +def _merge_by_inflow( + gen: MealyHMM, + blocks: list[frozenset[Hashable]], + labels: list[Hashable], + weights: list[dict[Hashable, float]], +) -> MealyHMM: + block_of = _block_of(blocks) + graph = TransitionGraph() + for label in labels: + graph.add_state(label) + result = MealyHMM(graph=graph, observation_alphabet=gen.observation_alphabet) + merged: dict[tuple[int, int, Any], float] = defaultdict(float) + for t in gen.transitions(): + source = block_of[t.source] + key = (source, block_of[t.target], t.data.get(ATTR_EMISSION)) + merged[key] += weights[source][t.source] * float(t.data[ATTR_PROB]) + for (source, target, symbol), prob in merged.items(): + if prob > 0.0: + result.add_transition(labels[source], labels[target], symbol, prob) + for state, mass in gen.initial_distribution.items(): + label = labels[block_of[state]] + result.initial_distribution[label] = result.initial_distribution.get(label, 0.0) + float(mass) + result.validate() + return result + + +__all__ = ["amalgamate", "in_edges", "out_edges", "split_state"] diff --git a/sofic/generators/support.py b/sofic/generators/support.py new file mode 100644 index 0000000..47904f9 --- /dev/null +++ b/sofic/generators/support.py @@ -0,0 +1,91 @@ +"""Support languages of stationary processes and the decisions built on them. + +The *stationary support* of an HMM is the set of finite words of positive +probability under its stationary law. It is a factorial regular language, so +comparisons between supports reduce to inclusion of finite automata, decided +here by antichain inclusion :cite:`DeWulf2006`. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np + +from sofic.graph import ATTR_SYMBOL + +if TYPE_CHECKING: + from sofic.automata.nfa import NFA + from sofic.generators.base import HiddenMarkovModel + +_TOL = 1e-12 + + +def support_nfa(model: HiddenMarkovModel) -> NFA: + """Return an NFA accepting exactly the words of positive stationary probability. + + States are the generator states of positive stationary probability, all of + them initial and accepting; edges are the positive-probability transitions. + The stationary support is closed under such transitions, so a word is + accepted iff some stationary-positive state can emit it. + + This differs from :meth:`~sofic.generators.base.HiddenMarkovModel.to_support_nfa`, + which starts from *every* state, transient ones included, and so accepts + words that only a transient state can emit (words of stationary + probability zero). + """ + from sofic.automata.nfa import NFA + from sofic.generators.matrices import symbol_matrices + + mealy = model.to_mealy() + states = tuple(mealy.reindex().states) + pi = np.asarray(mealy.stationary_distribution(), dtype=float) + matrices = {symbol: np.asarray(matrix, dtype=float) for symbol, matrix in symbol_matrices(mealy).items()} + live = [i for i, mass in enumerate(pi) if mass > _TOL] + live_states = frozenset(states[i] for i in live) + nfa = NFA( + input_alphabet=frozenset(matrices), + initial_states=live_states, + accepting_states=live_states, + ) + for state in live_states: + nfa.graph.add_state(state) + for symbol, matrix in matrices.items(): + for i in live: + for j in np.flatnonzero(matrix[i] > _TOL): + nfa.graph.add_transition(states[i], states[j], **{ATTR_SYMBOL: symbol}) + return nfa + + +def support_includes(p: HiddenMarkovModel, q: HiddenMarkovModel) -> bool: + """Return whether every word with positive ``P`` probability has positive ``Q`` probability. + + Both processes are taken in their stationary laws. Decided by antichain + inclusion of the :func:`support_nfa` automata :cite:`DeWulf2006`, without + determinizing. + """ + return support_nfa(q).includes(support_nfa(p)) + + +def is_absolutely_continuous(p: HiddenMarkovModel, q: HiddenMarkovModel) -> bool: + """Return whether ``P_{0:n} << Q_{0:n}`` for every block length ``n``. + + This is absolute continuity on the finite-dimensional cylinders, the + condition under which every block divergence ``D(P_{0:n} || Q_{0:n})`` and + the relative entropy rate can be finite; it is the same decision as + :func:`support_includes`. It is *not* absolute continuity of the laws on + infinite sequences: two distinct ergodic stationary measures are mutually + singular there, even when their finite-word supports coincide. + """ + return support_includes(p, q) + + +def support_equal(p: HiddenMarkovModel, q: HiddenMarkovModel) -> bool: + """Return whether ``P`` and ``Q`` give positive probability to the same finite words. + + Equivalently, ``P_{0:n}`` and ``Q_{0:n}`` are mutually absolutely continuous + for every ``n``: inclusion of the :func:`support_nfa` languages both ways + :cite:`DeWulf2006`. + """ + left, right = support_nfa(p), support_nfa(q) + return right.includes(left) and left.includes(right) diff --git a/sofic/inference/__init__.py b/sofic/inference/__init__.py index aa0d8d2..47c2b7d 100644 --- a/sofic/inference/__init__.py +++ b/sofic/inference/__init__.py @@ -1,13 +1,14 @@ """Inference algorithms for stochastic generators. Learners are named ``learn__``; the subpackages and modules -(:mod:`~sofic.inference.bayesian`, :mod:`~sofic.inference.cssr`, +(:mod:`~sofic.inference.active`, :mod:`~sofic.inference.bayesian`, :mod:`~sofic.inference.cssr`, :mod:`~sofic.inference.diagnostics`, :mod:`~sofic.inference.hmm`, :mod:`~sofic.inference.model_selection`, :mod:`~sofic.inference.spectral`) are bound under their module names. """ -from sofic.inference import bayesian, cssr, diagnostics, hmm, model_selection, spectral +from sofic.inference import active, bayesian, cssr, diagnostics, hmm, model_selection, spectral +from sofic.inference.active import ProcessOracle, learn_epsilon_machine_active from sofic.inference.cssr import ( learn_epsilon_machine_cssr, learn_epsilon_machine_subtree, @@ -63,6 +64,9 @@ ) __all__ = [ + "active", + "ProcessOracle", + "learn_epsilon_machine_active", "bayesian", "cssr", "diagnostics", diff --git a/sofic/inference/active.py b/sofic/inference/active.py new file mode 100644 index 0000000..47a9fa7 --- /dev/null +++ b/sofic/inference/active.py @@ -0,0 +1,494 @@ +r"""Active (query) learning of ε-machines: an L\*-style learner over predictive rows. + +Angluin's L\* :cite:`Angluin1987` learns a minimal DFA from membership and +equivalence queries. Its stochastic analogue here learns an ε-machine +:cite:`Shalizi2001` from a teacher answering + +* **probability queries** ``oracle.word_probability(w)``: the stationary + probability :math:`P(w)` of a finite word, and +* optionally **equivalence queries** ``oracle.equivalent(machine)``: ``None`` when + ``machine`` generates the target process, else a counterexample word whose + probability differs. + +Rows of the observation table are histories :math:`h` (prefixes) and a row's +content is the vector of conditional future probabilities :math:`P(s \mid h)` over +the test suffixes :math:`s`. A row is a finite-history *mixed state* +:cite:`Ellison2009`; distinct rows are distinct predictive distributions, so the +hypothesis is unifilar and minimal by construction. :class:`ProcessOracle` +answers both kinds of query exactly from any finite HMM. +""" + +from __future__ import annotations + +from collections.abc import Iterable, Sequence +from fractions import Fraction +from typing import Any, Protocol, cast, runtime_checkable + +import numpy as np + +from sofic.exceptions import MixedStateExplosionError, StochasticValidationError +from sofic.generators.base import HiddenMarkovModel +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.matrices import emission_tensors +from sofic.generators.mealy import MealyHMM +from sofic.generators.prob import as_prob, has_symbolic, is_symbolic, probs_equal, simplify_prob, zeros +from sofic.generators.process_equivalence import ( + _DEFAULT_ATOL, + _DEFAULT_RTOL, + _HistoryFutureWordList, + _sorted_words, + is_equal_process, +) +from sofic.generators.stationary import stationary_distribution_from_transition +from sofic.graph import ATTR_EMISSION, ATTR_PROB, TransitionGraph +from sofic.inference.cssr.process import _recurrent_states + +__all__ = [ + "ProbabilityOracle", + "ProcessOracle", + "learn_epsilon_machine_active", +] + +Word = tuple[Any, ...] + +_MASS_ATOL = 1e-12 + + +@runtime_checkable +class ProbabilityOracle(Protocol): + """Answers the stationary probability of a finite word.""" + + def word_probability(self, word: Sequence[Any]) -> Any: ... + + +class ProcessOracle: + r"""Probability and equivalence oracle answering exactly from a finite HMM. + + ``word_probability(w)`` is :math:`\pi T^{(w_0)} \cdots T^{(w_{L-1})} \mathbf{1}` + with :math:`\\pi` the stationary state law, in exact sympy arithmetic when the + HMM's probabilities are symbolic. ``equivalent(machine, history=h)`` compares + ``machine`` (from its initial distribution) with the target process conditioned + on the history ``h`` using + :func:`~sofic.generators.process_equivalence.is_equal_process`, and returns + ``None`` or a shortest word among the history/future basis words of the two + presentations whose probability differs. + + ``probability_queries`` and ``equivalence_queries`` count the calls. + """ + + def __init__(self, hmm: HiddenMarkovModel, *, rtol: float | None = None, atol: float | None = None) -> None: + self._hmm = hmm + self._start, self._matrices = emission_tensors(hmm, policy="stationary") + self._rtol = _DEFAULT_RTOL if rtol is None else float(rtol) + self._atol = _DEFAULT_ATOL if atol is None else float(atol) + self.probability_queries = 0 + self.equivalence_queries = 0 + + def _mass(self, word: Sequence[Any]) -> np.ndarray: + mass = self._start + n = len(mass) + for symbol in word: + matrix = self._matrices.get(symbol) + if matrix is None: + return zeros((n,)) + mass = mass @ matrix + return mass + + def word_probability(self, word: Sequence[Any]) -> Any: + self.probability_queries += 1 + mass = self._mass(tuple(word)) + if mass.dtype == object: + return simplify_prob(sum(mass.tolist())) + return float(mass.sum()) + + def equivalent(self, machine: HiddenMarkovModel, history: Sequence[Any] = ()) -> Word | None: + self.equivalence_queries += 1 + mass = np.asarray(self._mass(tuple(history)), dtype=float) + total = float(mass.sum()) + if total <= 0.0: + raise ValueError(f"history {tuple(history)!r} has zero probability") + target = _float_presentation(self._hmm) + start = dict(zip(target.reindex().states, mass / total, strict=True)) + hypothesis = _float_presentation(machine) + if is_equal_process(hypothesis, target, start2=start, rtol=self._rtol, atol=self._atol): + return None + return _process_counterexample(hypothesis, target, start, rtol=self._rtol, atol=self._atol) + + +def _float_presentation(hmm: HiddenMarkovModel) -> MealyHMM: + """``hmm`` as a Mealy presentation with float probabilities and the same state names. + + :func:`~sofic.generators.process_equivalence.is_equal_process` needs float + matrices: its rank tests fail on sympy-valued ones. + """ + mealy = hmm.to_mealy() + transitions = list(mealy.transitions()) + if not has_symbolic(t.data.get(ATTR_PROB, 0.0) for t in transitions) and not has_symbolic( + mealy.initial_distribution.values() + ): + return mealy + graph = TransitionGraph() + for state in mealy.reindex().states: + graph.add_state(state) + for transition in transitions: + data = dict(transition.data) + data[ATTR_PROB] = float(data[ATTR_PROB]) + graph.add_transition(transition.source, transition.target, **data) + return MealyHMM( + graph=graph, + initial_distribution={state: float(mass) for state, mass in mealy.initial_distribution.items()}, + observation_alphabet=mealy.observation_alphabet, + ) + + +def _process_counterexample( + machine: HiddenMarkovModel, + target: HiddenMarkovModel, + start: dict[Any, float], + *, + rtol: float, + atol: float, +) -> Word | None: + """Shortest word ``h + f`` over both presentations' basis words whose probability differs. + + :func:`~sofic.generators.process_equivalence.is_equal_process` compares future + probabilities and conditionals ``P(f | h)`` over these words, so a difference it + detects is a difference in ``P(h)`` or ``P(h + f)``. + """ + left = _HistoryFutureWordList.from_hmm(machine) + right = _HistoryFutureWordList.from_hmm(target, start=start) + histories = set(left.history_word_list()) | set(right.history_word_list()) + futures = set(left.future_word_list()) | set(right.future_word_list()) + futures |= {(symbol, *word) for word in list(futures) for symbol in left.alphabet} + best, best_gap = None, 0.0 + for word in _sorted_words({h + f for h in histories for f in futures}): + p = float(left.start @ left.future_vector(word)) + q = float(right.start @ right.future_vector(word)) + if not np.isclose(p, q, rtol=rtol, atol=atol): + return word + if abs(p - q) > best_gap: + best, best_gap = word, abs(p - q) + return best + + +def _is_exact(value: Any) -> bool: + return is_symbolic(value) or isinstance(value, (int, Fraction)) + + +def _is_zero(value: Any, *, atol: float = 0.0) -> bool: + if is_symbolic(value): + return probs_equal(value, 0) + if isinstance(value, (int, Fraction)): + return value == 0 + return float(value) <= atol + + +def _divide(numerator: Any, denominator: Any) -> Any: + if is_symbolic(numerator) or is_symbolic(denominator): + import sympy as sp + + return simplify_prob(sp.sympify(numerator) / sp.sympify(denominator)) + if _is_exact(numerator) and _is_exact(denominator): + return Fraction(numerator) / Fraction(denominator) + return float(numerator) / float(denominator) + + +def _values_match(left: Any, right: Any, tol: float) -> bool: + if _is_exact(left) and _is_exact(right): + if is_symbolic(left) or is_symbolic(right): + return probs_equal(left, right) + return left == right + return abs(float(left) - float(right)) <= tol + + +def _machine_prob(value: Any) -> Any: + if isinstance(value, Fraction): + import sympy as sp + + return sp.Rational(value.numerator, value.denominator) + return as_prob(value) + + +class _PredictiveTable: + """Observation table of histories (``seed + access``) by test suffixes. + + State rows are kept pairwise distinct, so the table is consistent by + construction (the Maler-Pnueli treatment of counterexamples). + """ + + def __init__(self, oracle: ProbabilityOracle, alphabet: Sequence[Any], seed: Word, tol: float) -> None: + self._oracle = oracle + self._alphabet = tuple(alphabet) + self._seed = seed + self._tol = tol + self._probabilities: dict[Word, Any] = {} + self.states: list[Word] = [()] + self.tests: list[Word] = [(symbol,) for symbol in self._alphabet] + + @property + def exact(self) -> bool: + return all(_is_exact(value) for value in self._probabilities.values()) + + def probability(self, word: Word) -> Any: + word = self._seed + word + value = self._probabilities.get(word) + if value is None: + value = self._oracle.word_probability(word) + self._probabilities[word] = value + return value + + def conditional(self, history: Word, suffix: Word) -> Any: + return _divide(self.probability(history + suffix), self.probability(history)) + + def row(self, history: Word) -> tuple[Any, ...]: + return tuple(self.conditional(history, test) for test in self.tests) + + def match(self, history: Word) -> int | None: + target = self.row(history) + for index, state in enumerate(self.states): + if all(_values_match(a, b, self._tol) for a, b in zip(self.row(state), target, strict=True)): + return index + return None + + def successors(self, history: Word) -> list[tuple[Any, Any]]: + """``(symbol, P(symbol | history))`` for each positive-probability next symbol.""" + out = [] + for symbol in self._alphabet: + if _is_zero(self.probability(history + (symbol,))): + continue + prob = self.conditional(history, (symbol,)) + if not _is_zero(prob, atol=_MASS_ATOL): + out.append((symbol, prob)) + return out + + def add_tests(self, counterexample: Word) -> None: + for index in range(len(counterexample)): + suffix = tuple(counterexample[index:]) + if suffix not in self.tests: + self.tests.append(suffix) + + def close(self, max_states: int) -> None: + """Add one-symbol extensions until every extension's row matches a state row.""" + index = 0 + while index < len(self.states): + history = self.states[index] + for symbol, _prob in self.successors(history): + extension = history + (symbol,) + if self.match(extension) is None: + self.states.append(extension) + if len(self.states) > max_states: + raise MixedStateExplosionError( + f"the predictive rows did not close within max_states={max_states}; the process " + "may not be exactly synchronizable from this history, or it has infinitely many " + "transient mixed states (seed a synchronizing `history`)" + ) + index += 1 + + def transitions(self) -> dict[int, dict[Any, tuple[int, Any]]]: + """``delta[state][symbol] = (target, probability)`` of the closed table.""" + delta: dict[int, dict[Any, tuple[int, Any]]] = {} + for index, history in enumerate(self.states): + delta[index] = {} + for symbol, prob in self.successors(history): + target = self.match(history + (symbol,)) + if target is None: # pragma: no cover - guaranteed by close() + raise RuntimeError("observation table is not closed") + delta[index][symbol] = (target, prob) + return delta + + def presentation(self, delta: dict[int, dict[Any, tuple[int, Any]]]) -> MealyHMM: + """The unifilar presentation of all state rows, started at the seed history's row.""" + return _build_machine(MealyHMM, sorted(delta), delta, {0: 1}, self._alphabet) + + +def _build_machine( + machine_type: type[MealyHMM], + keep: Sequence[int], + delta: dict[int, dict[Any, tuple[int, Any]]], + initial: dict[int, Any], + alphabet: Sequence[Any], +) -> MealyHMM: + labels = {state: f"s{rank}" for rank, state in enumerate(keep)} + graph = TransitionGraph() + for state in keep: + graph.add_state(labels[state]) + for state in keep: + for symbol, (target, prob) in delta[state].items(): + graph.add_transition( + labels[state], labels[target], **{ATTR_PROB: _machine_prob(prob), ATTR_EMISSION: symbol} + ) + machine = machine_type( + graph=graph, + initial_distribution={labels[state]: _machine_prob(mass) for state, mass in initial.items()}, + observation_alphabet=frozenset(alphabet), + ) + machine.validate() + return machine + + +def _recurrent_epsilon_machine(delta: dict[int, dict[Any, tuple[int, Any]]], alphabet: Sequence[Any]) -> EpsilonMachine: + """The ε-machine on the presentation's unique closed class, with its stationary law.""" + classes = _recurrent_states( + {state: {symbol: (target,) for symbol, (target, _) in edges.items()} for state, edges in delta.items()} + ) + if len(classes) != 1: + raise StochasticValidationError( + f"the learned presentation has {len(classes)} closed classes; the process is not ergodic" + ) + keep = sorted(classes[0]) + position = {state: i for i, state in enumerate(keep)} + probs = [_machine_prob(prob) for state in keep for _, prob in delta[state].values()] + transition = zeros((len(keep), len(keep)), symbolic=has_symbolic(probs)) + for state in keep: + for target, prob in delta[state].values(): + transition[position[state], position[target]] += _machine_prob(prob) + pi = stationary_distribution_from_transition(transition) + initial = {state: pi[position[state]] for state in keep} + return cast(EpsilonMachine, _build_machine(EpsilonMachine, keep, delta, initial, alphabet)) + + +def _presentation_probability(delta: dict[int, dict[Any, tuple[int, Any]]], word: Word) -> Any: + state, prob = 0, 1 + for symbol in word: + edge = delta[state].get(symbol) + if edge is None: + return 0 + state = edge[0] + prob = prob * edge[1] + return prob + + +def _bounded_counterexample( + table: _PredictiveTable, + delta: dict[int, dict[Any, tuple[int, Any]]], + alphabet: Sequence[Any], + *, + max_length: int, + tol: float, +) -> Word | None: + """Shortest word up to ``max_length`` where ``P(w | seed)`` and the presentation disagree.""" + frontier: list[Word] = [()] + for _ in range(max_length): + extended = [] + for word in frontier: + for symbol in alphabet: + candidate = word + (symbol,) + target = table.conditional((), candidate) + hypothesis = _presentation_probability(delta, candidate) + if not _values_match(target, hypothesis, tol): + return candidate + if not _is_zero(target): + extended.append(candidate) + frontier = extended + return None + + +def learn_epsilon_machine_active( + oracle: ProbabilityOracle, + alphabet: Iterable[Any], + *, + tol: float = 1e-9, + max_states: int = 1000, + max_rounds: int = 100, + history: Sequence[Any] = (), + max_counterexample_length: int = 8, +) -> EpsilonMachine: + r"""Learn an ε-machine from probability (and equivalence) queries, L\*-style :cite:`Angluin1987`. + + The observation table has histories ``history + u`` as rows and test suffixes + ``s`` as columns, with entries :math:`P(s \mid \text{history}\, u)` from + ``oracle.word_probability``. Tests start as the single symbols. The table is + kept *closed* -- every positive-probability one-symbol extension of a state row + matches some state row -- and *consistent*, since state rows are pairwise + distinct and adding tests only refines them. The hypothesis has one state per + distinct row and the transition :math:`P(x \mid h)` from row :math:`h` on + :math:`x`, so it is unifilar by construction; started at the row of + ``history`` it is a candidate for the conditional process + :math:`P(\cdot \mid \text{history})`. It is checked with + ``oracle.equivalent(hypothesis)`` (``oracle.equivalent(hypothesis, + history=history)`` when seeded), or, if the oracle has no ``equivalent``, + against every word of length at most ``max_counterexample_length``. Every + suffix of a counterexample becomes a test, and the loop repeats. The returned + ε-machine :cite:`Shalizi2001` is the hypothesis restricted to its closed class, + started at its stationary distribution. + + Parameters + ---------- + oracle + Has ``word_probability(word)`` (stationary probability) and optionally + ``equivalent(machine[, history])``; see :class:`ProcessOracle`. + alphabet + The process alphabet. + tol + Rows match when entries differ by at most ``tol``. Exact answers + (``int``, :class:`~fractions.Fraction`, sympy) are compared exactly. + max_states + Cap on table rows; exceeding it raises + :class:`~sofic.exceptions.MixedStateExplosionError`. + max_rounds + Cap on equivalence queries; exceeding it raises :class:`RuntimeError`. + history + Seed history prepended to every row (default empty). + max_counterexample_length + Word length searched when the oracle has no ``equivalent``. + + Notes + ----- + With exact answers and an exact equivalence oracle, let :math:`N` be the number + of distinct predictive distributions :math:`P(\cdot \mid \text{history}\, u)` over + positive-probability ``u``. If :math:`N` is finite, the learner stops after at + most :math:`N` equivalence queries (a counterexample whose suffixes are all tests + is predicted correctly from a closed table by the chain rule, so each one adds a + state) with exactly :math:`N` hypothesis states, and every table entry costs one + or two probability queries. For a stationary ergodic process whose ε-machine is + finite and exactly synchronizable :cite:`Travers2010`, every predictive state + reaches a causal state, so the closed class is the ε-machine: same number of + states, same process, same :math:`C_\mu` and :math:`h_\mu`. + + :math:`N` is finite when ``history`` is a synchronizing word (then all rows are + causal states) or when the mixed-state presentation :cite:`Ellison2009` from the + stationary law is finite. A finite ε-machine that is not exactly synchronizable + has infinitely many mixed states, as can an exactly synchronizable one along a + non-synchronizing transient; those exceed ``max_states`` unless seeded with a + synchronizing ``history``. With float answers the result is approximate: row + matching within ``tol`` is not transitive, and a counterexample that separates + no rows at ``tol`` ends learning with the current hypothesis. + + Examples + -------- + >>> from sofic.examples import even_process + >>> from sofic.inference.active import ProcessOracle, learn_epsilon_machine_active + >>> machine = learn_epsilon_machine_active(ProcessOracle(even_process(0.5)), ("0", "1")) + >>> len(list(machine.states())) + 2 + """ + symbols = tuple(sorted(set(alphabet), key=repr)) + seed = tuple(history) + table = _PredictiveTable(oracle, symbols, seed, float(tol)) + if _is_zero(table.probability(())): + raise ValueError(f"history {seed!r} has zero probability") + equivalent = getattr(oracle, "equivalent", None) + + table.close(max_states) + for _ in range(max_rounds): + delta = table.transitions() + if equivalent is None: + counterexample = _bounded_counterexample( + table, delta, symbols, max_length=max_counterexample_length, tol=float(tol) + ) + else: + hypothesis = table.presentation(delta) + counterexample = equivalent(hypothesis, history=seed) if seed else equivalent(hypothesis) + if counterexample is None: + return _recurrent_epsilon_machine(delta, symbols) + size = len(table.states) + table.add_tests(tuple(counterexample)) + table.close(max_states) + if len(table.states) == size: + if table.exact: + raise RuntimeError( + f"counterexample {tuple(counterexample)!r} separates no rows; the oracles are inconsistent" + ) + return _recurrent_epsilon_machine(delta, symbols) + + raise RuntimeError("active ε-machine learning did not converge within max_rounds; check the equivalence oracle") diff --git a/sofic/shifts/sliding_block_code.py b/sofic/shifts/sliding_block_code.py index 577e91c..62100bf 100644 --- a/sofic/shifts/sliding_block_code.py +++ b/sofic/shifts/sliding_block_code.py @@ -12,11 +12,15 @@ from collections.abc import Iterable, Mapping, Sequence from itertools import product -from typing import Any +from typing import TYPE_CHECKING, Any from sofic.graph import ATTR_OUTPUT, ATTR_SYMBOL, TransitionGraph from sofic.shifts.sofic import SoficShift +if TYPE_CHECKING: + from sofic.generators.base import HiddenMarkovModel + from sofic.generators.mealy import MealyHMM + class SlidingBlockCode: """A block map ``Phi`` with given memory and anticipation. @@ -113,6 +117,19 @@ def apply(self, shift: Any) -> SoficShift: image.symbol_alphabet = frozenset(used_outputs) return image.trim_transient() + def apply_to_process(self, hmm: HiddenMarkovModel) -> MealyHMM: + """Return an HMM generating the image process ``Phi(X)`` of ``hmm``. + + The image is ``Y_t = Phi(X_{t:t+window})`` :cite:`LindMarcus1995` + (§1.5): its length-``n`` word distribution is the pushforward of + ``hmm``'s length-``(n + window - 1)`` word distribution, and its entropy + rate is at most that of ``hmm`` :cite:`Gray1990`. See + :func:`~sofic.generators.factor_codes.image_process`. + """ + from sofic.generators.factor_codes import image_process + + return image_process(hmm, self) + def is_right_resolving(self, shift: Any | None = None) -> bool: """Return whether the induced image presentation is right-resolving.""" source = shift if shift is not None else full_shift(self.input_alphabet) diff --git a/tests/test_active_epsilon_learning.py b/tests/test_active_epsilon_learning.py new file mode 100644 index 0000000..a0cead3 --- /dev/null +++ b/tests/test_active_epsilon_learning.py @@ -0,0 +1,270 @@ +"""Tests for L*-style active learning of ε-machines.""" + +from __future__ import annotations + +from fractions import Fraction + +import numpy as np +import pytest +import sympy as sp +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.examples import butterfly_process, even_process, golden_mean, nemo_process, restricted_golden_mean +from sofic.examples.epsilon_machines import from_symbol_matrices +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.topological_epsilon_enumeration import iter_topological_epsilon_machines +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.inference.active import ProbabilityOracle, ProcessOracle, learn_epsilon_machine_active +from sofic.testing.strategies import epsilon_machines + +R = sp.Rational + + +def _n_states(machine): + return len(list(machine.states())) + + +def _alphabet(machine): + return sorted(machine.observation_alphabet, key=repr) + + +def _assert_recovers(learned, target): + assert isinstance(learned, EpsilonMachine) + assert _n_states(learned) == _n_states(target) + assert learned.is_equal_process(target) + assert learned.statistical_complexity() == pytest.approx(target.statistical_complexity(), abs=1e-9) + assert learned.entropy_rate() == pytest.approx(target.entropy_rate(), abs=1e-9) + + +def _with_probabilities(machine, rng): + """A new ε-machine on ``machine``'s topology with Dirichlet-random edge probabilities.""" + states = list(machine.states()) + symbols = _alphabet(machine) + index = {state: i for i, state in enumerate(states)} + matrices = {symbol: np.zeros((len(states), len(states))) for symbol in symbols} + for state in states: + out = list(machine.graph.out_transitions(state)) + for transition, prob in zip(out, rng.dirichlet(np.full(len(out), 2.0)), strict=True): + matrices[transition.data[ATTR_EMISSION]][index[state], index[transition.target]] = prob + return from_symbol_matrices(states, symbols, matrices) + + +def _synchronizing_word(machine): + """Shortest word whose subset-construction image is one state.""" + symbols = _alphabet(machine) + start = frozenset(machine.states()) + seen, queue = {start}, [(start, ())] + for subset, word in queue: + if len(subset) == 1: + return word + for symbol in symbols: + image = frozenset( + t.target for s in subset for t in machine.graph.out_transitions(s) if t.data[ATTR_EMISSION] == symbol + ) + if image and image not in seen: + seen.add(image) + queue.append((image, (*word, symbol))) + raise AssertionError("no synchronizing word") + + +def _exact(states, symbols, matrices): + return from_symbol_matrices(states, symbols, {k: np.array(v, dtype=object) for k, v in matrices.items()}) + + +def _exact_even(): + return _exact(("A", "B"), ("0", "1"), {"0": [[R(1, 3), 0], [0, 0]], "1": [[0, R(2, 3)], [1, 0]]}) + + +def _infinite_transients(exact=True): + """``1`` fixes both states with different probabilities, ``0`` swaps them, ``2`` synchronizes to ``A``. + + Along ``1^n`` the mixed states never repeat, so the mixed-state presentation is infinite. + """ + q = R(1, 4) if exact else 0.25 + a, b = (R(1, 2), R(1, 4)) if exact else (0.5, 0.25) + matrices = {"0": [[0, q], [q, 0]], "1": [[a, 0], [0, b]], "2": [[1 - a - q, 0], [1 - b - q, 0]]} + if exact: + return _exact(("A", "B"), ("0", "1", "2"), matrices) + return from_symbol_matrices(("A", "B"), ("0", "1", "2"), {k: np.array(v, dtype=float) for k, v in matrices.items()}) + + +class _ProbabilityOnly: + def __init__(self, target): + self._oracle = ProcessOracle(target) + + def word_probability(self, word): + return self._oracle.word_probability(word) + + +class _FractionOracle: + def __init__(self, target): + self._oracle = ProcessOracle(target) + + def word_probability(self, word): + value = self._oracle.word_probability(word) + return Fraction(int(value.p), int(value.q)) + + def equivalent(self, machine): + return self._oracle.equivalent(machine) + + +# --- named processes ------------------------------------------------------------ + + +@pytest.mark.parametrize( + "target", + [golden_mean(0.3), even_process(0.4), nemo_process(0.3, 0.6), butterfly_process(), restricted_golden_mean(3)], + ids=["golden_mean", "even", "nemo", "butterfly", "rgm3"], +) +def test_recovers_named_processes(target): + oracle = ProcessOracle(target) + learned = learn_epsilon_machine_active(oracle, _alphabet(target)) + _assert_recovers(learned, target) + assert oracle.equivalence_queries <= _n_states(target) + 2 + + +def test_learned_machine_is_unifilar_with_stationary_start(): + learned = learn_epsilon_machine_active(ProcessOracle(nemo_process(0.3, 0.6)), ("0", "1")) + learned.validate() + pi = learned.stationary_distribution() + start = [learned.initial_distribution[state] for state in learned.reindex().states] + assert np.allclose(start, pi) + + +def test_presentation_passed_to_oracle_starts_at_seed_row(): + seen = [] + + class Recording(ProcessOracle): + def equivalent(self, machine, history=()): + seen.append(machine) + return super().equivalent(machine, history) + + target = even_process(0.5) + learn_epsilon_machine_active(Recording(target), ("0", "1")) + presentation = seen[-1] + assert presentation.initial_distribution == {"s0": 1} + assert presentation.is_equal_process(target) + assert _n_states(presentation) == 4 + + +# --- random ε-machines ------------------------------------------------------------ + + +@pytest.mark.parametrize(("k", "n"), [(2, 1), (2, 2), (2, 3), (3, 2)]) +def test_recovers_random_enumerated_epsilon_machines(k, n): + rng = np.random.default_rng(1000 * k + n) + symbols = tuple("abc"[:k]) + for count, topology in enumerate(iter_topological_epsilon_machines(k, n, alphabet=symbols)): + if count >= 25: + break + target = _with_probabilities(topology, rng) + oracle = ProcessOracle(target) + learned = learn_epsilon_machine_active(oracle, symbols, history=_synchronizing_word(target)) + _assert_recovers(learned, target) + assert oracle.equivalence_queries <= n + 1 + + +@settings(max_examples=25, deadline=None) +@given(epsilon_machines(max_states=3), st.integers(min_value=0, max_value=2**32 - 1)) +def test_recovers_random_strategy_epsilon_machines(topology, seed): + target = _with_probabilities(topology, np.random.default_rng(seed)) + learned = learn_epsilon_machine_active( + ProcessOracle(target), _alphabet(target), history=_synchronizing_word(target) + ) + _assert_recovers(learned, target) + + +# --- exact arithmetic ------------------------------------------------------------ + + +def test_exact_sympy_oracle_gives_exact_machine(): + target = _exact_even() + learned = learn_epsilon_machine_active(ProcessOracle(target), ("0", "1")) + probs = sorted((t.data[ATTR_EMISSION], t.data[ATTR_PROB]) for t in learned.transitions()) + assert probs == [("0", R(1, 3)), ("1", R(2, 3)), ("1", 1)] + assert sorted(learned.initial_distribution.values()) == [R(2, 5), R(3, 5)] + assert sp.simplify(learned.statistical_complexity() - target.statistical_complexity()) == 0 + assert ProcessOracle(target).equivalent(learned) is None + + +def test_fraction_oracle_gives_exact_machine(): + learned = learn_epsilon_machine_active(_FractionOracle(_exact_even()), ("0", "1")) + assert sorted(t.data[ATTR_PROB] for t in learned.transitions()) == [R(1, 3), R(2, 3), 1] + + +def test_infinite_transients_explode_unless_seeded(): + target = _infinite_transients() + with pytest.raises(MixedStateExplosionError): + learn_epsilon_machine_active(ProcessOracle(target), ("0", "1", "2"), max_states=30) + oracle = ProcessOracle(target) + learned = learn_epsilon_machine_active(oracle, ("0", "1", "2"), history=("2",)) + assert _n_states(learned) == 2 + assert ProcessOracle(target).equivalent(learned) is None + assert oracle.equivalence_queries == 1 + + +def test_float_tolerance_truncates_converging_transients(): + target = _infinite_transients(exact=False) + learned = learn_epsilon_machine_active(ProcessOracle(target), ("0", "1", "2")) + _assert_recovers(learned, target) + + +# --- queries, termination, tolerance ------------------------------------------------ + + +def test_query_counts_are_polynomial(): + counts = [] + for k in range(1, 9): + oracle = ProcessOracle(restricted_golden_mean(k)) + learned = learn_epsilon_machine_active(oracle, ("0", "1")) + n = k + 1 + assert _n_states(learned) == n + assert oracle.equivalence_queries <= 2 + assert oracle.probability_queries <= 4 * n**2 + counts.append(oracle.probability_queries) + assert counts == sorted(counts) + + +def test_probability_only_oracle_uses_bounded_search(): + target = nemo_process(0.3, 0.6) + learned = learn_epsilon_machine_active(_ProbabilityOnly(target), ("0", "1")) + _assert_recovers(learned, target) + assert isinstance(_ProbabilityOnly(target), ProbabilityOracle) + + +def test_max_rounds_is_enforced(): + with pytest.raises(RuntimeError, match="max_rounds"): + learn_epsilon_machine_active(ProcessOracle(nemo_process(0.3, 0.6)), ("0", "1"), max_rounds=1) + + +def test_max_states_is_enforced(): + with pytest.raises(MixedStateExplosionError): + learn_epsilon_machine_active(ProcessOracle(butterfly_process()), _alphabet(butterfly_process()), max_states=3) + + +def test_zero_probability_history_is_rejected(): + with pytest.raises(ValueError, match="zero probability"): + learn_epsilon_machine_active(ProcessOracle(golden_mean(0.5)), ("0", "1"), history=("1", "1")) + + +def _near_iid(delta): + return from_symbol_matrices( + ("A", "B"), + ("0", "1"), + {"0": np.array([[0.5, 0.0], [0.5 + delta, 0.0]]), "1": np.array([[0.0, 0.5], [0.0, 0.5 - delta]])}, + ) + + +def test_tight_tolerance_resolves_close_states(): + target = _near_iid(0.01) + _assert_recovers(learn_epsilon_machine_active(ProcessOracle(target), ("0", "1"), tol=1e-9), target) + + +def test_loose_tolerance_merges_close_states(): + target = _near_iid(0.01) + learned = learn_epsilon_machine_active(ProcessOracle(target), ("0", "1"), tol=0.05) + assert _n_states(learned) == 1 + assert not learned.is_equal_process(target) + assert learned.entropy_rate() == pytest.approx(target.entropy_rate(), abs=1e-3) diff --git a/tests/test_bisimulation.py b/tests/test_bisimulation.py new file mode 100644 index 0000000..23692df --- /dev/null +++ b/tests/test_bisimulation.py @@ -0,0 +1,212 @@ +"""Tests for the coarsest probabilistic bisimulation (coarsest strong lumping).""" + +from fractions import Fraction + +import pytest +import sympy as sp +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.lumping import bisimulation_partition, coarsest_lumping, is_lumpable, lump +from sofic.generators.markov import MarkovChain +from sofic.generators.mealy import MealyHMM +from sofic.generators.moore import MooreHMM +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.testing.strategies import epsilon_machines, mealy_hmms +from tests import oracles +from tests.test_lumping import _moore, _symmetric_chain + +ALPHABET = (0, 1) + + +def _golden_mean() -> MealyHMM: + hmm = MealyHMM(initial_distribution={0: 2 / 3, 1: 1 / 3}) + hmm.add_transition(0, 0, 0, 0.5) + hmm.add_transition(0, 1, 1, 0.5) + hmm.add_transition(1, 0, 0, 1.0) + return hmm + + +def _set_partitions(items: list) -> list[list[list]]: + if not items: + return [[]] + first, rest = items[0], items[1:] + partitions = [] + for partition in _set_partitions(rest): + partitions.append([[first], *partition]) + for i in range(len(partition)): + partitions.append([*partition[:i], [first, *partition[i]], *partition[i + 1 :]]) + return partitions + + +def _refines(fine: list, coarse: list[frozenset]) -> bool: + return all(any(set(block) <= other for other in coarse) for block in fine) + + +def _clone_state(hmm: MealyHMM, state: int, route: list[float]) -> MealyHMM: + """Add a copy of ``state``; the ``i``-th edge into ``state`` sends ``route[i]`` of its mass to the original.""" + clone = max(hmm.states()) + 1 + out = MealyHMM(observation_alphabet=hmm.observation_alphabet) + for s in [*hmm.states(), clone]: + out.graph.add_state(s) + incoming = 0 + for transition in hmm.transitions(): + symbol, prob = transition.data[ATTR_EMISSION], transition.data[ATTR_PROB] + sources = [transition.source, clone] if transition.source == state else [transition.source] + for source in sources: + if transition.target != state: + out.add_transition(source, transition.target, symbol, prob) + continue + share = route[incoming % len(route)] + if share > 0: + out.add_transition(source, state, symbol, share * prob) + if share < 1: + out.add_transition(source, clone, symbol, (1 - share) * prob) + incoming += transition.target == state + initial = dict(hmm.initial_distribution) + mass = initial.pop(state, 0.0) + out.initial_distribution = {**initial, state: mass / 2, clone: mass / 2} + return out + + +@st.composite +def redundant_hmms(draw, *, max_states=4, unifilar_routes=False): + hmm = draw(mealy_hmms(alphabet=ALPHABET, max_states=max_states)) + state = draw(st.sampled_from(sorted(hmm.states()))) + shares = st.sampled_from([0.0, 1.0]) if unifilar_routes else st.sampled_from([0.0, 0.25, 0.5, 1.0]) + return _clone_state(hmm, state, draw(st.lists(shares, min_size=1, max_size=4))) + + +any_small_hmm = st.one_of(mealy_hmms(alphabet=ALPHABET, max_states=5), redundant_hmms()) + + +@settings(max_examples=80, deadline=None) +@given(any_small_hmm) +def test_partition_is_lumpable(hmm): + assert is_lumpable(hmm, bisimulation_partition(hmm)) + + +@settings(max_examples=60, deadline=None) +@given(any_small_hmm) +def test_partition_is_coarsest(hmm): + blocks = bisimulation_partition(hmm) + for candidate in _set_partitions(sorted(hmm.states())): + if is_lumpable(hmm, candidate): + assert _refines(candidate, blocks) + + +@settings(max_examples=60, deadline=None) +@given(any_small_hmm) +def test_lumped_model_generates_same_process(hmm): + lumped = coarsest_lumping(hmm) + assert len(list(lumped.states())) == len(bisimulation_partition(hmm)) + for word in oracles.words(ALPHABET, 4): + assert oracles.word_probability(lumped, word) == pytest.approx(oracles.word_probability(hmm, word), abs=1e-12) + + +@settings(max_examples=60, deadline=None) +@given(epsilon_machines(alphabet=ALPHABET, max_states=4)) +def test_epsilon_machine_is_bisimulation_minimal(machine): + assert all(len(block) == 1 for block in bisimulation_partition(machine)) + + +@settings(max_examples=60, deadline=None) +@given(mealy_hmms(alphabet=ALPHABET, max_states=4), st.data()) +def test_duplicated_state_collapses(hmm, data): + state = data.draw(st.sampled_from(sorted(hmm.states()))) + redundant = _clone_state(hmm, state, [0.5]) + clone = max(redundant.states()) + blocks = bisimulation_partition(redundant) + assert any({state, clone} <= block for block in blocks) + assert len(blocks) == len(bisimulation_partition(hmm)) + + +@settings(max_examples=60, deadline=None) +@given(epsilon_machines(alphabet=ALPHABET, max_states=4), st.data()) +def test_unifilar_quotient_matches_epsilon_machine(machine, data): + state = data.draw(st.sampled_from(sorted(machine.states()))) + routes = data.draw(st.lists(st.sampled_from([0.0, 1.0]), min_size=1, max_size=4)) + redundant = _clone_state(machine, state, routes) + assert redundant.is_unifilar() + blocks = bisimulation_partition(redundant) + assert len(blocks) == len(list(EpsilonMachine.from_hmm(redundant).states())) + assert len(blocks) == len(list(machine.states())) + + +def test_lump_without_partition_uses_coarsest(): + hmm = _clone_state(_golden_mean(), 1, [0.5]) + assert bisimulation_partition(hmm) == [frozenset({0}), frozenset({1, 2})] + assert sorted(lump(hmm).states(), key=str) == [0, "1+2"] + + +def test_markov_chain_without_labels_lumps_to_one_block(): + assert bisimulation_partition(_symmetric_chain()) == [frozenset({"A", "B", "C"})] + + +def test_markov_chain_refines_initial_partition(): + chain = _symmetric_chain() + blocks = bisimulation_partition(chain, initial={"A": 0, "B": 1, "C": 1}) + assert blocks == [frozenset({"A"}), frozenset({"B", "C"})] + assert bisimulation_partition(chain, initial=[{"A"}, {"B"}, {"C"}]) == [ + frozenset({"A"}), + frozenset({"B"}), + frozenset({"C"}), + ] + + +def test_moore_partition_respects_emissions(): + assert bisimulation_partition(_moore()) == [frozenset({"A"}), frozenset({"B", "C"})] + assert isinstance(coarsest_lumping(_moore()), MooreHMM) + + +def test_markov_chain_quotient_type(): + lumped = coarsest_lumping(_symmetric_chain(), initial=[{"A"}, {"B", "C"}]) + assert isinstance(lumped, MarkovChain) + assert sorted(lumped.states()) == ["A", "B+C"] + + +def test_symbolic_probabilities_compare_exactly(): + a = sp.Symbol("a", positive=True) + hmm = MealyHMM(initial_distribution={0: sp.Integer(1)}) + hmm.add_transition(0, 1, 0, a / 2) + hmm.add_transition(0, 2, 0, a / 2) + hmm.add_transition(0, 0, 1, 1 - a) + hmm.add_transition(1, 0, 1, sp.Integer(1)) + hmm.add_transition(2, 0, 1, sp.Integer(1)) + assert bisimulation_partition(hmm) == [frozenset({0}), frozenset({1, 2})] + + +def test_fraction_probabilities_compare_exactly(): + chain = MarkovChain(initial_distribution={0: 1.0}) + for state in range(3): + chain.graph.add_state(state) + tiny = Fraction(1, 10**20) + chain.graph.add_transition(0, 1, **{ATTR_PROB: Fraction(1, 2)}) + chain.graph.add_transition(0, 2, **{ATTR_PROB: Fraction(1, 2)}) + chain.graph.add_transition(1, 0, **{ATTR_PROB: Fraction(1, 2) + tiny}) + chain.graph.add_transition(1, 1, **{ATTR_PROB: Fraction(1, 2) - tiny}) + chain.graph.add_transition(2, 0, **{ATTR_PROB: Fraction(1, 2)}) + chain.graph.add_transition(2, 2, **{ATTR_PROB: Fraction(1, 2)}) + blocks = bisimulation_partition(chain, initial=[{0}, {1, 2}]) + assert blocks == [frozenset({0}), frozenset({1}), frozenset({2})] + + +def test_float_tolerance_merges_rounding_noise(): + chain = MarkovChain(initial_distribution={0: 1.0}) + for state in range(3): + chain.graph.add_state(state) + chain.add_transition(0, 1, 0.5) + chain.add_transition(0, 2, 0.5) + chain.add_transition(1, 0, 0.1 + 0.2) + chain.add_transition(1, 1, 0.7) + chain.add_transition(2, 0, 0.3) + chain.add_transition(2, 2, 0.7) + assert bisimulation_partition(chain, initial=[{0}, {1, 2}]) == [frozenset({0}), frozenset({1, 2})] + + +def test_unsupported_model_raises(): + from sofic.automata.dfa import DFA + + with pytest.raises(TypeError): + bisimulation_partition(DFA()) diff --git a/tests/test_canonical_residual_hmm.py b/tests/test_canonical_residual_hmm.py new file mode 100644 index 0000000..2f9db58 --- /dev/null +++ b/tests/test_canonical_residual_hmm.py @@ -0,0 +1,170 @@ +"""Canonical residual HMM: extreme future morphs of a finite ε-machine.""" + +from __future__ import annotations + +import numpy as np +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.examples import bernoulli, even_process, golden_mean +from sofic.examples.epsilon_machines import from_symbol_matrices +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.canonical_residual import canonical_residual_hmm +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.matrices import symbol_matrices +from sofic.generators.mealy import MealyHMM +from sofic.testing.strategies import epsilon_machines + +MAX_LENGTH = 6 + + +def _same_process(left, right, max_length: int = MAX_LENGTH) -> None: + for n in range(1, max_length + 1): + a, b = left.words_of_length(n), right.words_of_length(n) + for word in set(a) | set(b): + assert a.get(word, 0.0) == pytest.approx(b.get(word, 0.0), abs=1e-9), word + + +def _hmm(edges: list[tuple[str, str, str, float]]) -> MealyHMM: + hmm = MealyHMM(observation_alphabet=frozenset(symbol for _s, _t, symbol, _p in edges)) + for source, target, symbol, prob in edges: + hmm.graph.add_state(source) + hmm.graph.add_state(target) + hmm.add_transition(source, target, symbol, prob) + pi = hmm.stationary_distribution() + hmm.initial_distribution = dict(zip(hmm.reindex().states, map(float, pi), strict=True)) + hmm.validate() + return hmm + + +def _mixing_hmm(q: float = 0.3, r: float = 0.4, split: float = 0.5, s: float = 0.2, t: float = 0.5) -> MealyHMM: + """Two-state generator whose ε-machine has a third causal state, a mixture of the other two. + + ``2`` and ``3`` synchronize to ``a`` and ``b``; ``0`` is emitted with + probability ``q`` from both states and keeps the belief; ``1`` is emitted + only by ``a`` and moves to ``a`` or ``b`` in proportion ``split``. The + recurrent beliefs are ``a``, ``b`` and ``(split, 1 - split)``. + """ + return _hmm( + [ + ("a", "a", "0", q), + ("a", "a", "1", r * split), + ("a", "b", "1", r * (1 - split)), + ("a", "a", "2", (1 - q - r) * s), + ("a", "b", "3", (1 - q - r) * (1 - s)), + ("b", "b", "0", q), + ("b", "a", "2", (1 - q) * t), + ("b", "b", "3", (1 - q) * (1 - t)), + ] + ) + + +def _lohr_ay(p: float = 0.1) -> MealyHMM: + """Example 3.6 of :cite:`LohrAy2009`: a two-state generator with infinitely many causal states.""" + return _hmm( + [ + ("0", "0", "0", 1 - 2 * p), + ("0", "0", "1", p), + ("0", "1", "1", p), + ("1", "1", "1", 1 - 2 * p), + ("1", "0", "0", p), + ("1", "1", "0", p), + ] + ) + + +@st.composite +def weighted_machines(draw, max_states: int = 3) -> EpsilonMachine: + """Topological ε-machine skeleton with random positive edge probabilities.""" + skeleton = draw(epsilon_machines(max_states=max_states)) + states = list(skeleton.reindex().states) + matrices = {symbol: np.array(m, dtype=float) for symbol, m in symbol_matrices(skeleton).items()} + for i in range(len(states)): + edges = [(symbol, j) for symbol, m in matrices.items() for j in np.flatnonzero(m[i])] + weights = np.array([draw(st.floats(0.1, 1.0)) for _ in edges]) + for (symbol, j), weight in zip(edges, weights / weights.sum(), strict=True): + matrices[symbol][i, j] = weight + return from_symbol_matrices(states, tuple(matrices), matrices) + + +def test_requires_experimental_flag() -> None: + with pytest.raises(RuntimeError, match="experimental"): + canonical_residual_hmm(golden_mean(0.5)) + + +@pytest.mark.parametrize("model", [golden_mean(0.4), even_process(0.3), bernoulli(0.3)], ids=["gm", "even", "iid"]) +def test_affinely_independent_morphs_return_the_epsilon_machine(model: EpsilonMachine) -> None: + residual = canonical_residual_hmm(model, experimental=True) + assert set(residual.states()) == set(model.states()) + expected = {(t.source, t.target, t.data["emission"]): t.data["prob"] for t in model.transitions()} + actual = {(t.source, t.target, t.data["emission"]): t.data["prob"] for t in residual.transitions()} + assert actual.keys() == expected.keys() + for key, prob in expected.items(): + assert actual[key] == pytest.approx(prob, abs=1e-9) + _same_process(residual, model) + + +def test_strictly_smaller_than_epsilon_machine() -> None: + hmm = _mixing_hmm() + eps = EpsilonMachine.from_hmm(hmm) + assert len(list(eps.states())) == 3 + residual = canonical_residual_hmm(eps, experimental=True) + residual.validate() + assert len(list(residual.states())) == 2 + assert not residual.is_unifilar() + _same_process(residual, eps) + _same_process(residual, hmm) + + +def test_non_unifilar_input_is_converted() -> None: + hmm = _mixing_hmm() + assert len(list(canonical_residual_hmm(hmm, experimental=True).states())) == 2 + + +def test_infinite_epsilon_machine_raises() -> None: + with pytest.raises(MixedStateExplosionError): + canonical_residual_hmm(_lohr_ay(), experimental=True, max_states=200) + + +def test_duplicate_states_are_merged() -> None: + model = from_symbol_matrices( + ("A", "B"), + ("0", "1"), + {"0": np.array([[0.3, 0.0], [0.0, 0.3]]), "1": np.array([[0.0, 0.7], [0.7, 0.0]])}, + ) + residual = canonical_residual_hmm(model, experimental=True) + assert len(list(residual.states())) == 1 + _same_process(residual, bernoulli(0.7)) + + +@settings(max_examples=40) +@given(weighted_machines()) +def test_random_epsilon_machines(model: EpsilonMachine) -> None: + residual = canonical_residual_hmm(model, experimental=True) + residual.validate() + size = len(list(model.states())) + assert len(list(residual.states())) <= size + _same_process(residual, model, max_length=5) + futures = np.array( + [ + [model.word_probability(word, start=state) for n in range(size) for word in model.words_of_length(n)] + for state in model.reindex().states + ] + ) + if np.linalg.matrix_rank(futures, tol=1e-6) == size: + assert set(residual.states()) == set(model.states()) + + +unit = st.floats(0.1, 0.9) + + +@settings(max_examples=25) +@given(st.floats(0.05, 0.45), st.floats(0.05, 0.45), unit, unit, unit) +def test_mixed_causal_state_is_dropped(q: float, r: float, split: float, s: float, t: float) -> None: + hmm = _mixing_hmm(q, r, split, s, t) + eps = EpsilonMachine.from_hmm(hmm) + assert len(list(eps.states())) == 3 + residual = canonical_residual_hmm(eps, experimental=True) + assert len(list(residual.states())) == 2 + _same_process(residual, hmm, max_length=4) diff --git a/tests/test_conditioning.py b/tests/test_conditioning.py new file mode 100644 index 0000000..7aff509 --- /dev/null +++ b/tests/test_conditioning.py @@ -0,0 +1,158 @@ +"""Conditioning a stationary process on a regular constraint (Doob h-transform).""" + +from __future__ import annotations + +import math + +import pytest + +from sofic.automata.dfa import DFA +from sofic.automata.nfa import NFA +from sofic.examples import bernoulli, binary_markov_chain, golden_mean +from sofic.exceptions import NonDeterministicError +from sofic.generators.conditioning import condition_on_language +from sofic.generators.mealy import MealyHMM +from sofic.shifts.sft import ShiftOfFiniteType +from sofic.shifts.topological_anatomy import parry_measure_sofic + +PHI = (1 + math.sqrt(5)) / 2 + + +def _uniform(symbols: str) -> MealyHMM: + hmm = MealyHMM(observation_alphabet=frozenset(symbols), initial_distribution={"A": 1.0}) + for symbol in symbols: + hmm.add_transition("A", "A", symbol, 1 / len(symbols)) + return hmm + + +def _sft(forbidden: set[str], symbols: str) -> ShiftOfFiniteType: + return ShiftOfFiniteType.from_forbidden_words({tuple(word) for word in forbidden}, frozenset(symbols)) + + +def _language_dfa(shift: ShiftOfFiniteType) -> DFA: + states = frozenset(shift.states()) + nfa = NFA( + graph=shift.graph.copy(), + input_alphabet=shift.symbol_alphabet, + initial_states=states, + accepting_states=states, + ) + return nfa.determinize() + + +def _dfa(edges: list[tuple[str, str, str]], accepting: str, initial: str = "q") -> DFA: + dfa = DFA( + input_alphabet=frozenset("01"), initial_states=frozenset({initial}), accepting_states=frozenset(accepting) + ) + dfa.graph.add_state(initial) + for source, target, _symbol in edges: + dfa.graph.add_state(source) + dfa.graph.add_state(target) + for source, target, symbol in edges: + dfa.add_transition(source, target, symbol) + return dfa + + +def _stationary_words(model, length: int) -> dict[tuple[str, ...], float]: + return model.word_probabilities(length, start=model.to_mealy().stationary_distribution()) + + +def _window_law(hmm, dfa: DFA, n: int, offset: int, width: int) -> dict[tuple[str, ...], float]: + """Brute-force law of ``X_{offset:offset+width}`` given ``X_{0:n} ∈ L``.""" + law: dict[tuple[str, ...], float] = {} + for word, prob in _stationary_words(hmm, n).items(): + if dfa.recognizes(word): + window = word[offset : offset + width] + law[window] = law.get(window, 0.0) + prob + total = sum(law.values()) + return {window: prob / total for window, prob in law.items()} + + +def _distance(a: dict, b: dict) -> float: + return max(abs(a.get(word, 0.0) - b.get(word, 0.0)) for word in set(a) | set(b)) + + +def test_fair_coin_conditioned_on_golden_mean_is_the_parry_measure(): + shift = _sft({"11"}, "01") + conditioned = condition_on_language(bernoulli(), _language_dfa(shift)) + assert conditioned.is_equal_process(golden_mean(1 / PHI), rtol=1e-7, atol=1e-10) + assert conditioned.entropy_rate() == pytest.approx(shift.topological_entropy(), abs=1e-10) + assert conditioned.entropy_rate() == pytest.approx(math.log2(PHI), abs=1e-10) + + +@pytest.mark.parametrize( + ("forbidden", "symbols"), + [({"11"}, "01"), ({"111"}, "01"), ({"11", "20"}, "012"), ({"00", "12", "21"}, "012")], +) +def test_uniform_iid_conditioned_on_sft_is_its_parry_measure(forbidden, symbols): + shift = _sft(forbidden, symbols) + conditioned = condition_on_language(_uniform(symbols), _language_dfa(shift)) + parry = parry_measure_sofic(shift) + for length in range(1, 6): + assert _distance(_stationary_words(conditioned, length), _stationary_words(parry, length)) < 1e-10 + assert conditioned.entropy_rate() == pytest.approx(shift.topological_entropy(), abs=1e-10) + + +@pytest.mark.parametrize("model", [binary_markov_chain(0.4, 0.3), binary_markov_chain(0.7, 0.2), golden_mean(0.3)]) +def test_conditioned_markov_chain_is_the_limit_of_conditioned_windows(model): + dfa = _language_dfa(_sft({"00"}, "01")) + conditioned = condition_on_language(model, dfa) + width = 2 + limit = _stationary_words(conditioned, width) + errors = [_distance(_window_law(model, dfa, n, (n - width) // 2, width), limit) for n in range(4, 13, 2)] + assert all(later <= earlier + 1e-12 for earlier, later in zip(errors, errors[1:], strict=False)) + assert errors[-1] < 1e-4 + + +def test_conditioning_on_the_full_language_changes_nothing(): + universal = _dfa([("q", "q", "0"), ("q", "q", "1")], "q") + model = binary_markov_chain(0.4, 0.3) + conditioned = condition_on_language(model, universal) + for length in range(1, 5): + assert _distance(_stationary_words(conditioned, length), _stationary_words(model, length)) < 1e-10 + + +def test_incomplete_dfa_rejects_missing_transitions(): + no_ones = _dfa([("q", "q", "0")], "q") + conditioned = condition_on_language(bernoulli(0.3), no_ones) + assert _stationary_words(conditioned, 3) == pytest.approx({("0", "0", "0"): 1.0}) + assert conditioned.entropy_rate() == pytest.approx(0.0, abs=1e-12) + + +_SPLIT = [("s", "a", "0"), ("s", "b", "1"), ("a", "a", "0"), ("b", "b", "1")] + + +def test_two_maximal_components_are_rejected(): + split = _dfa(_SPLIT, "sab", initial="s") + with pytest.raises(ValueError, match="maximal growth"): + condition_on_language(bernoulli(), split) + + +def test_dominant_component_wins(): + split = _dfa(_SPLIT, "sab", initial="s") + conditioned = condition_on_language(bernoulli(0.8), split) + assert _stationary_words(conditioned, 2) == pytest.approx({("1", "1"): 1.0}) + + +def test_empty_constraints_are_rejected(): + rejecting = _dfa([("q", "q", "0")], "") + with pytest.raises(ValueError, match="empty word"): + condition_on_language(bernoulli(), rejecting) + finite = _dfa([("q", "r", "0")], "qr") + with pytest.raises(ValueError, match="no infinite word"): + condition_on_language(bernoulli(), finite) + + +def test_nondeterministic_constraint_is_rejected(): + nfa = NFA(input_alphabet=frozenset("01"), initial_states=frozenset({"q", "r"}), accepting_states=frozenset({"q"})) + nfa.add_transition("q", "q", "0") + with pytest.raises(NonDeterministicError): + condition_on_language(bernoulli(), nfa) + + +def test_a_process_already_in_the_constraint_is_unchanged(): + model = golden_mean(0.3) + conditioned = condition_on_language(model, _language_dfa(_sft({"11"}, "01"))) + assert conditioned.is_equal_process(model, rtol=1e-7, atol=1e-10) + assert conditioned.observation_alphabet == frozenset("01") + assert {state[0] for state in conditioned.states()} <= set(model.states()) diff --git a/tests/test_conversions.py b/tests/test_conversions.py index c401650..5bb82e8 100644 --- a/tests/test_conversions.py +++ b/tests/test_conversions.py @@ -83,15 +83,15 @@ def test_hmm_to_support_nfa_uses_fresh_epsilon_start(): nfa.validate() -def test_hmm_to_support_dfa_starts_from_all_states_and_accepts_recurrent_subsets(): +def test_hmm_to_support_dfa_starts_from_all_states_and_accepts_every_nonempty_subset(): hmm = _golden_mean_support_hmm() dfa = hmm.to_support_dfa() assert isinstance(dfa, DFA) assert dfa.initial_states == frozenset({frozenset({"A", "B"})}) assert frozenset() not in set(dfa.states()) - assert dfa.accepting_states == frozenset({frozenset({"A"}), frozenset({"B"})}) - assert not dfa.recognizes(()) + assert dfa.accepting_states == frozenset(dfa.states()) + assert dfa.recognizes(()) assert dfa.recognizes((0,)) assert dfa.recognizes((1,)) assert dfa.recognizes((1, 0, 1)) @@ -135,3 +135,17 @@ def to_mealy(self) -> MealyHMM: wrapped = WrappedHMM(_golden_mean_support_hmm()) assert wrapped.to_support_nfa().recognizes((1, 0, 1)) assert not wrapped.to_support_nfa().recognizes((1, 1)) + + +def test_hmm_to_support_dfa_accepts_every_word_of_positive_probability(): + """Regression: only terminal recurrent subsets accepted, rejecting e.g. ``1`` for the even process.""" + from itertools import product + + from sofic.examples import even_process, nemo_process + + for model in (even_process(), nemo_process()): + dfa = model.to_support_dfa() + alphabet = sorted(model.observation_alphabet) + for length in range(6): + for word in product(alphabet, repeat=length): + assert dfa.recognizes(word) == (model.word_probability(word) > 0) diff --git a/tests/test_dfa.py b/tests/test_dfa.py index 9f136e2..dc1a6cd 100644 --- a/tests/test_dfa.py +++ b/tests/test_dfa.py @@ -120,3 +120,10 @@ def test_standard_regular_operations_methods(): assert star.recognizes(()) assert star.recognizes(("a",)) assert not star.recognizes(("b",)) + + +def test_add_transition_creates_missing_source_state(): + """Regression: ``DFA.add_transition`` raised ``NetworkXError`` for a new source, unlike ``NFA``.""" + dfa = DFA(input_alphabet=frozenset("ab"), initial_states=frozenset({0}), accepting_states=frozenset({1})) + dfa.add_transition(0, 1, "a") + assert dfa.recognizes("a") diff --git a/tests/test_epsilon_machine.py b/tests/test_epsilon_machine.py index ce23850..b4b43a3 100644 --- a/tests/test_epsilon_machine.py +++ b/tests/test_epsilon_machine.py @@ -270,3 +270,29 @@ def test_from_hmm_honors_max_states(): with pytest.raises(MixedStateExplosionError): EpsilonMachine.from_hmm(sns(), max_states=3) + + +def test_from_hmm_merges_states_split_only_by_round_off(): + """Regression: belief-update round-off (0.49999999999999994 vs 0.5) split causal states. + + Every state of this non-unifilar HMM emits a fair coin, so the process is i.i.d. + and its ε-machine has one state; exact float comparison gave 33. + """ + import pytest + + from sofic.generators.mealy import MealyHMM + + hmm = MealyHMM(initial_distribution={"A": 1.0}, observation_alphabet=frozenset({0, 1})) + for state in ("A", "B", "C"): + hmm.graph.add_state(state) + hmm.add_transition("A", "B", 0, 0.5) + hmm.add_transition("A", "C", 1, 0.5) + hmm.add_transition("B", "A", 0, 0.5) + hmm.add_transition("B", "B", 1, 0.25) + hmm.add_transition("B", "C", 1, 0.25) + hmm.add_transition("C", "A", 0, 0.5) + hmm.add_transition("C", "B", 1, 0.5) + + machine = EpsilonMachine.from_hmm(hmm) + assert len(list(machine.states())) == 1 + assert machine.statistical_complexity() == pytest.approx(0.0, abs=1e-12) diff --git a/tests/test_factor_codes.py b/tests/test_factor_codes.py new file mode 100644 index 0000000..1697c5b --- /dev/null +++ b/tests/test_factor_codes.py @@ -0,0 +1,179 @@ +"""Sliding block codes on processes: image processes and higher block presentations. + +Word distributions are checked against brute-force pushforwards of the +path-enumeration oracle in :mod:`tests.oracles`. +""" + +from __future__ import annotations + +import math +from collections import defaultdict +from itertools import product +from typing import Any + +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.examples import golden_mean +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.factor_codes import higher_block, image_process +from sofic.generators.measures import entropy_rate_bounds +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.shifts.sliding_block_code import SlidingBlockCode +from sofic.testing.strategies import epsilon_machines, mealy_hmms +from tests import oracles + +ALPHABET = ("0", "1") +OUTPUTS = ("a", "b", "c") +TOL = 1e-9 + +small_hmms = mealy_hmms(alphabet=ALPHABET, max_states=3) +small_eps = epsilon_machines(alphabet=ALPHABET, max_states=3) +reversible_eps = small_eps.filter(lambda m: m.reverse_is_finite()) +heavy = settings(max_examples=max(10, settings.default.max_examples // 4)) + + +@st.composite +def codes(draw: Any, max_window: int = 3) -> SlidingBlockCode: + window = draw(st.integers(1, max_window)) + memory = draw(st.integers(0, window - 1)) + blocks = list(product(ALPHABET, repeat=window)) + outputs = draw(st.lists(st.sampled_from(OUTPUTS), min_size=len(blocks), max_size=len(blocks))) + return SlidingBlockCode(dict(zip(blocks, outputs, strict=True)), memory=memory, anticipation=window - 1 - memory) + + +def pushforward(hmm: Any, code: SlidingBlockCode, n: int) -> dict[tuple[Any, ...], float]: + """Law of ``Phi(X_{0:n+w-1})`` by brute-force enumeration of source words.""" + image: dict[tuple[Any, ...], float] = defaultdict(float) + for word, prob in oracles.word_distribution(hmm, n + code.window - 1, ALPHABET).items(): + image[code.apply_word(word)] += prob + return image + + +def entropy(distribution: dict[Any, float]) -> float: + return -sum(p * math.log2(p) for p in distribution.values() if p > 0.0) + + +def state_and_context_entropy(machine: EpsilonMachine, length: int) -> float: + """``H[S_0, X_{-length:0}]`` for a stationary unifilar machine.""" + frontier = {(state, ()): float(mass) for state, mass in machine.initial_distribution.items()} + for _ in range(length): + advanced: dict[tuple[Any, tuple[Any, ...]], float] = defaultdict(float) + for (state, word), mass in frontier.items(): + for edge in machine.graph.out_transitions(state): + advanced[(edge.target, (*word, edge.data[ATTR_EMISSION]))] += mass * float(edge.data[ATTR_PROB]) + frontier = advanced + return entropy(frontier) + + +# --------------------------------------------------------------------------- image processes + + +@given(small_hmms, codes(), st.integers(0, 2)) +def test_image_word_distribution_is_pushforward(hmm, code, n): + image = code.apply_to_process(hmm) + expected = pushforward(hmm, code, n) + for word in product(OUTPUTS, repeat=n): + assert image.word_probability(word) == pytest.approx(expected.get(word, 0.0), abs=TOL) + + +@given(small_hmms, codes(), st.integers(1, 3)) +def test_image_block_entropy_bounded_by_source_window(hmm, code, n): + image = code.apply_to_process(hmm) + source = oracles.word_distribution(hmm, n + code.window - 1, ALPHABET) + assert oracles.block_entropy(image.words_of_length(n)) <= oracles.block_entropy(source) + TOL + + +@heavy +@given(small_eps, codes(max_window=2)) +def test_factor_code_does_not_increase_entropy_rate(machine, code): + image = code.apply_to_process(machine) + h_source = float(machine.entropy_rate()) + lower, _upper = entropy_rate_bounds(image, 8) + assert lower <= h_source + TOL + try: + h_image = float(EpsilonMachine.from_hmm(image, max_states=500).entropy_rate()) + except MixedStateExplosionError: + return + assert h_image <= h_source + TOL + + +def test_image_matches_apply_word_on_golden_mean(): + xor = SlidingBlockCode({(a, b): str(int(a) ^ int(b)) for a, b in product(ALPHABET, repeat=2)}, memory=1) + image = xor.apply_to_process(golden_mean(0.5)) + assert image.words_of_length(2) == pytest.approx( + {("0", "0"): 1 / 6, ("0", "1"): 1 / 6, ("1", "0"): 1 / 6, ("1", "1"): 1 / 2} + ) + + +def test_constant_code_collapses_entropy_rate(): + machine = golden_mean(0.5) + constant = SlidingBlockCode({(a,): "*" for a in ALPHABET}) + image = constant.apply_to_process(machine) + assert image.observation_alphabet == frozenset({"*"}) + assert float(image.entropy_rate()) == pytest.approx(0.0, abs=TOL) + + +def test_image_is_stationary_for_stationary_source(): + code = SlidingBlockCode({word: word.count("1") for word in product(ALPHABET, repeat=3)}, memory=1, anticipation=1) + image = image_process(golden_mean(0.5), code) + assert image.is_stationary() + + +def test_missing_block_raises(): + partial = SlidingBlockCode({("0",): "a"}, input_alphabet=ALPHABET) + with pytest.raises(ValueError, match="not in block_map"): + partial.apply_to_process(golden_mean(0.5)) + + +# --------------------------------------------------------------------------- higher block presentations + + +@given(small_hmms, st.integers(1, 3), st.integers(0, 2)) +def test_higher_block_words_are_recoded_source_words(hmm, k, n): + recoded = hmm.higher_block(k) + source = oracles.word_distribution(hmm, n + k - 1, ALPHABET) if n else {(): 1.0} + expected = {tuple(word[t : t + k] for t in range(n)): prob for word, prob in source.items() if prob > 0.0} + observed = recoded.words_of_length(n) + assert set(observed) == set(expected) + for word, prob in expected.items(): + assert observed[word] == pytest.approx(prob, abs=TOL) + + +@given(small_hmms, st.integers(1, 3)) +def test_one_block_code_inverts_higher_block(hmm, k): + recoded = higher_block(hmm, k) + first = SlidingBlockCode( + {(block,): block[0] for block in recoded.observation_alphabet}, output_alphabet=hmm.observation_alphabet + ) + assert first.apply_to_process(recoded).is_equal_process(hmm) + + +@given(small_eps, st.integers(1, 3)) +def test_higher_block_preserves_unifilarity_and_entropy_rate(machine, k): + recoded = machine.higher_block(k) + assert recoded.is_unifilar() + assert float(recoded.entropy_rate()) == pytest.approx(float(machine.entropy_rate()), abs=1e-9) + + +@heavy +@given(reversible_eps, st.integers(1, 3)) +def test_higher_block_shifts_excess_entropy_by_k_minus_one_entropy_rates(machine, k): + recoded = EpsilonMachine.from_hmm(machine.higher_block(k)) + h = float(machine.entropy_rate()) + assert float(recoded.excess_entropy()) == pytest.approx(float(machine.excess_entropy()) + (k - 1) * h, abs=1e-6) + + +@given(small_eps, st.integers(1, 3)) +def test_higher_block_statistical_complexity_is_state_and_context_entropy(machine, k): + recoded = EpsilonMachine.from_hmm(machine.higher_block(k)) + expected = state_and_context_entropy(machine, k - 1) + assert float(recoded.statistical_complexity()) == pytest.approx(expected, abs=1e-6) + assert float(recoded.statistical_complexity()) >= float(machine.statistical_complexity()) - 1e-9 + + +def test_higher_block_rejects_nonpositive_k(): + with pytest.raises(ValueError, match="at least 1"): + golden_mean(0.5).higher_block(0) diff --git a/tests/test_hmm_state_splitting.py b/tests/test_hmm_state_splitting.py new file mode 100644 index 0000000..2c9d3c2 --- /dev/null +++ b/tests/test_hmm_state_splitting.py @@ -0,0 +1,251 @@ +"""Process-preserving state splitting and amalgamation of Mealy HMMs. + +Metamorphic checks: a split generates the same word distributions (brute force +via :mod:`tests.oracles`), entropy rate, excess entropy and ε-machine, and +amalgamating the copies recovers the original presentation. +""" + +from __future__ import annotations + +from collections import defaultdict +from typing import Any + +import pytest +import sympy +from hypothesis import assume, given, settings +from hypothesis import strategies as st + +from sofic.examples import golden_mean +from sofic.exceptions import LumpabilityError, MixedStateExplosionError +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.mealy import MealyHMM +from sofic.generators.measures import entropy_rate_bounds +from sofic.generators.state_splitting import amalgamate, in_edges, out_edges, split_state +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(alphabet=ALPHABET, max_states=3) +small_eps = epsilon_machines(alphabet=ALPHABET, max_states=3) +reversible_eps = small_eps.filter(lambda m: m.reverse_is_finite()) +heavy = settings(max_examples=max(10, settings.default.max_examples // 4)) + + +@st.composite +def splits(draw: Any, models: Any, kind: str) -> tuple[MealyHMM, Any, list[list[Any]]]: + """Draw ``(hmm, state, partition)`` with at least two parts.""" + hmm = draw(models) + state = draw(st.sampled_from(sorted(hmm.states(), key=repr))) + edges = out_edges(hmm, state) if kind == "out" else in_edges(hmm, state) + assume(len(edges) >= 2) + parts = draw(st.integers(2, len(edges))) + labels = list(range(parts)) + draw( + st.lists(st.integers(0, parts - 1), min_size=len(edges) - parts, max_size=len(edges) - parts) + ) + order = draw(st.permutations(range(len(edges)))) + partition: list[list[Any]] = [[] for _ in range(parts)] + for index, label in zip(order, labels, strict=True): + partition[label].append(edges[index]) + return hmm, state, partition + + +def edge_table(hmm: MealyHMM) -> dict[tuple[Any, Any, Any], float]: + table: dict[tuple[Any, Any, Any], float] = defaultdict(float) + for t in hmm.transitions(): + table[(t.source, t.target, t.data[ATTR_EMISSION])] += float(t.data[ATTR_PROB]) + return dict(table) + + +def assert_same_presentation(left: MealyHMM, right: MealyHMM) -> None: + assert set(left.states()) == set(right.states()) + left_edges, right_edges = edge_table(left), edge_table(right) + assert set(left_edges) == set(right_edges) + for key, prob in left_edges.items(): + assert right_edges[key] == pytest.approx(prob, abs=TOL) + states = set(left.initial_distribution) | set(right.initial_distribution) + for state in states: + assert float(left.initial_distribution.get(state, 0.0)) == pytest.approx( + float(right.initial_distribution.get(state, 0.0)), abs=TOL + ) + + +def assert_same_words(left: MealyHMM, right: MealyHMM, max_length: int = 3) -> None: + for n in range(max_length + 1): + expected = oracles.word_distribution(left, n, ALPHABET) + observed = oracles.word_distribution(right, n, ALPHABET) + for word, prob in expected.items(): + assert observed[word] == pytest.approx(prob, abs=TOL) + + +def copies_of(state: Any, partition: list[list[Any]]) -> frozenset[Any]: + return frozenset((state, index) for index in range(len(partition))) + + +def has_inflow(hmm: MealyHMM, state: Any) -> bool: + """Whether out-split weights of ``state`` are recorded in its copies' inflow.""" + return bool(in_edges(hmm, state)) or float(hmm.initial_distribution.get(state, 0.0)) > 0.0 + + +def rebuild_epsilon_machine(split: MealyHMM) -> EpsilonMachine: + """ε-machine of a split; out-splits are non-unifilar and may have infinitely many mixed states.""" + try: + return EpsilonMachine.from_hmm(split, max_states=200) + except MixedStateExplosionError: + assume(False) + raise + + +def merge_back(split: MealyHMM, state: Any, partition: list[list[Any]]) -> MealyHMM: + block = copies_of(state, partition) + blocks = [block, *({other} for other in split.states() if other not in block)] + return amalgamate(split, blocks, labels={block: state}) + + +# --------------------------------------------------------------------------- word distributions + + +@pytest.mark.parametrize("kind", ["out", "in"]) +@given(data=st.data()) +def test_split_preserves_word_distributions(kind, data): + hmm, state, partition = data.draw(splits(small_hmms, kind)) + split = split_state(hmm, state, partition, kind=kind) + split.validate() + assert_same_words(hmm, split) + assert split.is_equal_process(hmm) + + +@pytest.mark.parametrize("kind", ["out", "in"]) +@given(data=st.data()) +def test_amalgamating_a_split_recovers_the_original(kind, data): + hmm, state, partition = data.draw(splits(small_hmms, kind)) + assume(kind == "in" or has_inflow(hmm, state)) + split = split_state(hmm, state, partition, kind=kind) + assert_same_presentation(merge_back(split, state, partition), hmm) + + +@given(data=st.data()) +def test_amalgamating_an_unreachable_out_split_preserves_the_process(data): + hmm, state, partition = data.draw(splits(small_hmms, "out")) + split = split_state(hmm, state, partition) + assert merge_back(split, state, partition).is_equal_process(hmm) + + +@given(data=st.data()) +def test_in_then_out_split_amalgamates_in_one_call(data): + hmm, first, in_partition = data.draw(splits(small_hmms, "in")) + split = split_state(hmm, first, in_partition, kind="in") + second = data.draw(st.sampled_from(sorted(set(hmm.states()) - {first}, key=repr) or [first])) + assume(second != first and len(out_edges(split, second)) >= 2 and has_inflow(split, second)) + edges = out_edges(split, second) + out_partition = [edges[:1], edges[1:]] + twice = split_state(split, second, out_partition, kind="out") + assert_same_words(hmm, twice, max_length=2) + + in_block, out_block = copies_of(first, in_partition), copies_of(second, out_partition) + rest = [{state} for state in twice.states() if state not in in_block | out_block] + merged = amalgamate(twice, [in_block, out_block, *rest], labels={in_block: first, out_block: second}) + assert_same_presentation(merged, hmm) + + +# --------------------------------------------------------------------------- information measures + + +@pytest.mark.parametrize("kind", ["out", "in"]) +@heavy +@given(data=st.data()) +def test_split_preserves_entropy_rate_and_epsilon_machine(kind, data): + machine, state, partition = data.draw(splits(small_eps, kind)) + split = split_state(machine, state, partition, kind=kind) + rebuilt = rebuild_epsilon_machine(split) + assert rebuilt.is_equal_process(machine) + assert len(list(rebuilt.states())) == len(list(machine.states())) + assert float(rebuilt.statistical_complexity()) == pytest.approx(float(machine.statistical_complexity()), abs=1e-9) + assert float(rebuilt.entropy_rate()) == pytest.approx(float(machine.entropy_rate()), abs=1e-9) + lower, upper = entropy_rate_bounds(split, 6) + assert lower - TOL <= float(machine.entropy_rate()) <= upper + TOL + + +@pytest.mark.parametrize("kind", ["out", "in"]) +@heavy +@given(data=st.data()) +def test_split_preserves_excess_entropy(kind, data): + machine, state, partition = data.draw(splits(reversible_eps, kind)) + rebuilt = rebuild_epsilon_machine(split_state(machine, state, partition, kind=kind)) + assert float(rebuilt.excess_entropy()) == pytest.approx(float(machine.excess_entropy()), abs=1e-6) + + +@pytest.mark.parametrize("kind", ["out", "in"]) +@given(data=st.data()) +def test_stationary_initial_stays_stationary(kind, data): + machine, state, partition = data.draw(splits(small_eps, kind)) + assert split_state(machine, state, partition, kind=kind).is_stationary() + + +def test_in_split_preserves_unifilarity(): + machine = golden_mean(0.5) + edges = in_edges(machine, "A") + split = split_state(machine, "A", [edges[:1], edges[1:]], kind="in") + assert split.is_unifilar() + assert split.is_equal_process(machine) + + +def test_out_split_is_exact_for_symbolic_probabilities(): + p = sympy.Symbol("p", positive=True) + hmm = MealyHMM(initial_distribution={"A": 1 / (1 + p), "B": p / (1 + p)}) + hmm.graph.add_state("A") + hmm.graph.add_state("B") + hmm.add_transition("A", "A", "0", 1 - p) + hmm.add_transition("A", "B", "1", p) + hmm.add_transition("B", "A", "0", 1) + edges = out_edges(hmm, "A") + split = split_state(hmm, "A", [edges[:1], edges[1:]]) + for state in split.states(): + row = sum(t.data[ATTR_PROB] for t in split.graph.out_transitions(state)) + assert sympy.simplify(row - 1) == 0 + assert sympy.simplify(sum(split.initial_distribution.values()) - 1) == 0 + + value = sympy.Rational(1, 3) + numeric_split, numeric_hmm = substitute(split, p, value), substitute(hmm, p, value) + assert numeric_split.is_equal_process(numeric_hmm) + + +def substitute(hmm: MealyHMM, symbol: sympy.Symbol, value: Any) -> MealyHMM: + numeric = MealyHMM( + observation_alphabet=hmm.observation_alphabet, + initial_distribution={ + s: float(sympy.sympify(m).subs(symbol, value)) for s, m in hmm.initial_distribution.items() + }, + ) + for state in hmm.states(): + numeric.graph.add_state(state) + for t in hmm.transitions(): + prob = float(sympy.sympify(t.data[ATTR_PROB]).subs(symbol, value)) + numeric.add_transition(t.source, t.target, t.data[ATTR_EMISSION], prob) + return numeric + + +# --------------------------------------------------------------------------- errors + + +def test_split_rejects_bad_partitions(): + machine = golden_mean(0.5) + edges = out_edges(machine, "A") + with pytest.raises(ValueError, match="unknown state"): + split_state(machine, "Z", [edges]) + with pytest.raises(ValueError, match="empty part"): + split_state(machine, "A", [edges, []]) + with pytest.raises(ValueError, match="do not cover"): + split_state(machine, "A", [edges[:1]]) + with pytest.raises(ValueError, match="disjoint"): + split_state(machine, "A", [edges, edges[:1]]) + with pytest.raises(ValueError, match="kind"): + split_state(machine, "A", [edges], kind="sideways") + + +def test_amalgamate_rejects_inequivalent_states(): + machine = golden_mean(0.5) + with pytest.raises(LumpabilityError): + amalgamate(machine, [{"A", "B"}]) diff --git a/tests/test_mealy_hmm.py b/tests/test_mealy_hmm.py index 3d61150..266f17d 100644 --- a/tests/test_mealy_hmm.py +++ b/tests/test_mealy_hmm.py @@ -102,3 +102,16 @@ def test_mealy_without_alphabet_or_initial_distribution(): assert len(symbols) == 50 assert hmm.word_probability("ab") == pytest.approx(2 / 3 * 0.1) assert hmm.entropy_rate() == pytest.approx(0.5533064, abs=1e-6) + + +def test_entropy_with_mixed_int_and_tuple_state_labels(): + """Regression: dit could not sort labels mixing ``0`` and an encoded tuple.""" + hmm = MealyHMM(initial_distribution={0: 1.0}, observation_alphabet=frozenset({"0", "1"})) + for state in (0, (1, 0)): + hmm.graph.add_state(state) + hmm.add_transition(0, (1, 0), "0", 0.5) + hmm.add_transition(0, 0, "1", 0.5) + hmm.add_transition((1, 0), 0, "1", 1.0) + + assert hmm.entropy_rate() == pytest.approx(2 / 3) + assert hmm.state_entropy() == pytest.approx(0.9182958340544896) diff --git a/tests/test_minimal_quasi_realization.py b/tests/test_minimal_quasi_realization.py new file mode 100644 index 0000000..039fee6 --- /dev/null +++ b/tests/test_minimal_quasi_realization.py @@ -0,0 +1,186 @@ +"""Tests for Schützenberger/Fliess minimal quasi-realizations and the process rank.""" + +from functools import cache + +import numpy as np +import pytest +import sympy as sp +from hypothesis import assume, given, settings +from hypothesis import strategies as st + +from sofic.examples import even_process, golden_mean, iid +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.mealy import MealyHMM +from sofic.generators.minimal_quasi_realization import minimal_quasi_realization, process_rank +from sofic.generators.quasi_realization import QuasiRealization +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.testing.strategies import epsilon_machines, mealy_hmms +from tests import oracles + +ALPHABET = (0, 1) + + +def _duplicate_state(hmm: MealyHMM, state: int, weight: float) -> MealyHMM: + """Split ``state`` into itself and a clone, sending ``weight`` of its inflow to the original.""" + clone = max(hmm.states()) + 1 + out = MealyHMM(observation_alphabet=hmm.observation_alphabet) + for s in [*hmm.states(), clone]: + out.graph.add_state(s) + for transition in hmm.transitions(): + symbol, prob = transition.data[ATTR_EMISSION], transition.data[ATTR_PROB] + sources = [transition.source, clone] if transition.source == state else [transition.source] + for source in sources: + if transition.target == state: + out.add_transition(source, state, symbol, weight * prob) + out.add_transition(source, clone, symbol, (1 - weight) * prob) + else: + out.add_transition(source, transition.target, symbol, prob) + initial = dict(hmm.initial_distribution) + mass = initial.pop(state, 0.0) + out.initial_distribution = {**initial, state: weight * mass, clone: (1 - weight) * mass} + return out + + +def _words(max_length: int) -> list[tuple[int, ...]]: + return oracles.words(ALPHABET, max_length) + + +def _hankel_rank(hmm: MealyHMM, length: int) -> int: + probability = cache(lambda word: oracles.word_probability(hmm, word)) + prefixes = _words(length) + hankel = np.array([[probability(u + v) for v in prefixes] for u in prefixes]) + singular = np.linalg.svd(hankel, compute_uv=False) + return int(np.sum(singular > 1e-9 * singular[0])) + + +@settings(max_examples=60, deadline=None) +@given(mealy_hmms(alphabet=ALPHABET, max_states=3)) +def test_reproduces_every_word_probability(hmm): + realization = minimal_quasi_realization(hmm) + for word in _words(5): + assert realization.word_probability(word) == pytest.approx(oracles.word_probability(hmm, word), abs=1e-9) + + +@settings(max_examples=60, deadline=None) +@given(mealy_hmms(alphabet=ALPHABET, max_states=3)) +def test_dimension_equals_brute_force_hankel_rank(hmm): + n = len(list(hmm.states())) + assert process_rank(hmm) == _hankel_rank(hmm, n) + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(alphabet=ALPHABET, max_states=3), st.data()) +def test_dimension_bounded_by_every_presentation(hmm, data): + rank = process_rank(hmm) + assert 1 <= rank <= len(list(hmm.states())) + + state = data.draw(st.sampled_from(sorted(hmm.states()))) + weight = data.draw(st.sampled_from([0.25, 0.5, 0.75])) + redundant = _duplicate_state(hmm, state, weight) + assert process_rank(redundant) == rank + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(alphabet=ALPHABET, max_states=3)) +def test_stationary_dimension_bounded_by_epsilon_machine_of_hmm(hmm): + assume(hmm.is_irreducible()) + try: + machine = EpsilonMachine.from_hmm(hmm, max_states=200) + except MixedStateExplosionError: + assume(False) + stationary = hmm.copy() + stationary.initial_distribution = dict(zip(hmm.reindex().states, hmm.stationary_distribution(), strict=True)) + assert process_rank(stationary) <= len(list(machine.states())) + + +@settings(max_examples=40, deadline=None) +@given(epsilon_machines(alphabet=ALPHABET, max_states=4)) +def test_dimension_bounded_by_epsilon_machine(machine): + assert process_rank(machine) <= len(list(machine.states())) + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(alphabet=ALPHABET, max_states=3)) +def test_idempotent(hmm): + once = minimal_quasi_realization(hmm) + twice = minimal_quasi_realization(once) + assert twice.pi.size == once.pi.size + for word in _words(4): + assert twice.word_probability(word) == pytest.approx(once.word_probability(word), abs=1e-9) + + +def test_iid_has_rank_one(): + assert process_rank(iid(3)) == 1 + + +def test_golden_mean_has_rank_two(): + assert golden_mean(0.5).process_rank() == 2 + + +def test_even_process_has_rank_two(): + assert process_rank(even_process()) == 2 + + +def test_redundant_states_reduce(): + hmm = MealyHMM(observation_alphabet=frozenset(ALPHABET), initial_distribution={0: 0.5, 1: 0.25, 2: 0.25}) + hmm.add_transition(0, 1, 0, 0.25) + hmm.add_transition(0, 2, 0, 0.25) + hmm.add_transition(0, 0, 1, 0.5) + for state in (1, 2): + hmm.add_transition(state, 0, 1, 1.0) + assert len(list(hmm.states())) == 3 + realization = minimal_quasi_realization(hmm) + assert realization.pi.size == 2 + assert realization.pi.sum() == pytest.approx(1.0) + np.testing.assert_allclose(realization.tau, np.ones(2)) + + +def test_returns_quasi_realization_with_unit_final_vector(): + realization = minimal_quasi_realization(golden_mean(0.5)) + assert isinstance(realization, QuasiRealization) + realization.validate_quasistochastic() + np.testing.assert_allclose(realization.tau, np.ones(2)) + + +def test_exact_rational_input_stays_exact(): + third = sp.Rational(1, 3) + hmm = MealyHMM(initial_distribution={0: sp.Rational(3, 4), 1: sp.Rational(1, 4)}) + hmm.add_transition(0, 0, 0, third) + hmm.add_transition(0, 1, 1, 1 - third) + hmm.add_transition(1, 0, 0, sp.Integer(1)) + realization = minimal_quasi_realization(hmm) + assert realization.pi.dtype == object + assert realization.pi.size == 2 + for word in _words(4): + result = realization.pi + for symbol in word: + result = result @ realization.symbol_maps[symbol] + assert sp.simplify(result @ realization.tau - _exact_word_probability(hmm, word)) == 0 + + +def test_symbolic_parameter(): + a = sp.Symbol("a", positive=True) + hmm = MealyHMM(initial_distribution={0: sp.Integer(1)}) + hmm.add_transition(0, 0, 0, a) + hmm.add_transition(0, 1, 1, 1 - a) + hmm.add_transition(1, 0, 0, sp.Integer(1)) + realization = minimal_quasi_realization(hmm) + assert realization.pi.size == 2 + word = (0, 1, 0) + result = realization.pi + for symbol in word: + result = result @ realization.symbol_maps[symbol] + assert sp.simplify(result @ realization.tau - a * (1 - a)) == 0 + + +def _exact_word_probability(hmm: MealyHMM, word: tuple[int, ...]) -> sp.Expr: + states = sorted(hmm.states()) + vector = {s: sp.sympify(hmm.initial_distribution.get(s, 0)) for s in states} + for symbol in word: + nxt = dict.fromkeys(states, sp.Integer(0)) + for transition in hmm.transitions(): + if transition.data[ATTR_EMISSION] == symbol: + nxt[transition.target] += vector[transition.source] * transition.data[ATTR_PROB] + vector = nxt + return sp.simplify(sum(vector.values())) diff --git a/tests/test_omega_probability.py b/tests/test_omega_probability.py new file mode 100644 index 0000000..e62f0ae --- /dev/null +++ b/tests/test_omega_probability.py @@ -0,0 +1,172 @@ +"""Probabilities of regular and deterministic-Büchi properties of HMM output.""" + +from __future__ import annotations + +from itertools import product + +import numpy as np +import pytest + +from sofic.automata.buchi import BuchiAutomaton +from sofic.automata.dfa import DFA +from sofic.examples import bernoulli, binary_markov_chain, even_process, golden_mean +from sofic.exceptions import NonDeterministicError +from sofic.generators.mealy import MealyHMM +from sofic.generators.omega_probability import omega_probability, regular_language_probability + + +def _automaton(cls, edges, accepting: str, initial: str = "a"): + aut = cls( + input_alphabet=frozenset("01"), initial_states=frozenset({initial}), accepting_states=frozenset(accepting) + ) + aut.graph.add_state(initial) + for source, target, _symbol in edges: + aut.graph.add_state(source) + aut.graph.add_state(target) + for source, target, symbol in edges: + aut.add_transition(source, target, symbol) + return aut + + +# State b is entered exactly on reading a 1. +INFINITELY_MANY_ONES = [("a", "a", "0"), ("a", "b", "1"), ("b", "a", "0"), ("b", "b", "1")] +# b is absorbing once a 1 (resp. 11) has been read; a, c never saw one. +EVENTUALLY_ONE = [("a", "a", "0"), ("a", "b", "1"), ("b", "b", "0"), ("b", "b", "1")] +ALWAYS_ZERO = [("a", "a", "0")] +EVENTUALLY_ONE_ONE = [ + ("a", "a", "0"), + ("a", "c", "1"), + ("c", "a", "0"), + ("c", "b", "1"), + ("b", "b", "0"), + ("b", "b", "1"), +] +NEVER_ONE_ONE = [("a", "a", "0"), ("a", "c", "1"), ("c", "a", "0")] + + +def _gf1() -> BuchiAutomaton: + return _automaton(BuchiAutomaton, INFINITELY_MANY_ONES, "b") + + +def _split(mass_to_coin: float = 0.3) -> MealyHMM: + """From ``S`` emit ``0`` and move to a fair coin ``B`` or to the all-zeros state ``Z``.""" + hmm = MealyHMM(observation_alphabet=frozenset("01"), initial_distribution={"S": 1.0}) + hmm.add_transition("S", "B", "0", mass_to_coin) + hmm.add_transition("S", "Z", "0", 1 - mass_to_coin) + hmm.add_transition("B", "B", "0", 0.5) + hmm.add_transition("B", "B", "1", 0.5) + hmm.add_transition("Z", "Z", "0", 1.0) + return hmm + + +def _stops_emitting_ones(stop: float = 0.25) -> MealyHMM: + """A fair coin that switches to emitting ``0`` forever with probability ``stop`` per step.""" + hmm = MealyHMM(observation_alphabet=frozenset("01"), initial_distribution={"C": 1.0}) + hmm.add_transition("C", "C", "0", 0.5 * (1 - stop)) + hmm.add_transition("C", "C", "1", 0.5 * (1 - stop)) + hmm.add_transition("C", "Z", "0", stop) + hmm.add_transition("Z", "Z", "0", 1.0) + return hmm + + +@pytest.mark.parametrize("p", [0.01, 0.3, 0.5, 0.99]) +def test_infinitely_many_ones_is_almost_sure_for_iid(p): + assert omega_probability(bernoulli(p), _gf1()) == pytest.approx(1.0) + + +@pytest.mark.parametrize("model", [golden_mean(0.4), even_process(0.3), binary_markov_chain(0.2, 0.7)]) +def test_infinitely_many_ones_is_almost_sure_for_ergodic_processes(model): + assert omega_probability(model, _gf1()) == pytest.approx(1.0) + + +def test_infinitely_many_ones_fails_once_ones_stop(): + model = _stops_emitting_ones() + assert omega_probability(model, _gf1(), start="C") == pytest.approx(0.0) + assert omega_probability(model, _gf1(), start={"C": 1.0}) == pytest.approx(0.0) + # 'Eventually always 0' (finitely many 1s) is not DBA-expressible; it is the complement. + assert 1 - omega_probability(_split(0.3), _gf1(), start="S") == pytest.approx(0.7) + assert omega_probability(_split(0.3), _gf1(), start="B") == pytest.approx(1.0) + + +def test_eventually_one_is_the_mass_that_reaches_the_coin(): + eventually = _automaton(BuchiAutomaton, EVENTUALLY_ONE, "b") + assert omega_probability(bernoulli(0.5), eventually) == pytest.approx(1.0) + assert omega_probability(_split(0.0001), eventually, start="S") == pytest.approx(0.0001) + assert omega_probability(_split(0.3), eventually, start="Z") == pytest.approx(0.0) + + +@pytest.mark.parametrize( + ("prop", "complement"), + [ + ((EVENTUALLY_ONE, "b"), (ALWAYS_ZERO, "a")), + ((EVENTUALLY_ONE_ONE, "b"), (NEVER_ONE_ONE, "ac")), + ], +) +@pytest.mark.parametrize("mass", [0.0, 0.3, 0.8]) +def test_dba_and_complement_dba_sum_to_one(prop, complement, mass): + model = _split(mass) + p = omega_probability(model, _automaton(BuchiAutomaton, *prop), start="S") + q = omega_probability(model, _automaton(BuchiAutomaton, *complement), start="S") + assert p + q == pytest.approx(1.0) + assert p == pytest.approx(mass) + + +def test_never_one_one_has_stationary_probability_zero_unless_forbidden(): + never = _automaton(BuchiAutomaton, NEVER_ONE_ONE, "ac") + assert omega_probability(golden_mean(0.4), never) == pytest.approx(1.0) + assert omega_probability(bernoulli(0.3), never) == pytest.approx(0.0) + + +def test_safety_probability_is_the_limit_of_prefix_probabilities(): + never = _automaton(DFA, NEVER_ONE_ONE, "ac") + model = _stops_emitting_ones(0.4) + omega = omega_probability(model, _automaton(BuchiAutomaton, NEVER_ONE_ONE, "ac"), start="C") + finite = [regular_language_probability(model, never, n, start="C") for n in (10, 40, 160)] + assert finite[0] >= finite[1] >= finite[2] >= omega - 1e-12 + assert finite[-1] == pytest.approx(omega, abs=1e-9) + assert 0.0 < omega < 1.0 + + +def test_monte_carlo_agrees_with_omega_probability(): + model = _split(0.3) + exact = omega_probability(model, _gf1(), start="S") + rng = np.random.default_rng(20261007) + runs = 1000 + hits = sum("1" in model.sample(40, rng)[0][-15:] for _ in range(runs)) + assert hits / runs == pytest.approx(exact, abs=4 * np.sqrt(exact * (1 - exact) / runs)) + + +def test_nondeterministic_buchi_is_rejected(): + nba = _automaton(BuchiAutomaton, [("a", "a", "0"), ("a", "a", "1"), ("a", "b", "0"), ("b", "b", "0")], "b") + with pytest.raises(NonDeterministicError): + omega_probability(bernoulli(), nba) + + +def _brute_language_probability(model, accepts, n: int, start=None) -> float: + if start is None: + start = model.to_mealy().stationary_distribution() + return sum(model.word_probability(word, start=start) for word in product("01", repeat=n) if accepts(word)) + + +@pytest.mark.parametrize("model", [bernoulli(0.3), golden_mean(0.4), even_process(0.6), binary_markov_chain(0.2, 0.7)]) +@pytest.mark.parametrize("n", [0, 1, 2, 5, 8]) +def test_regular_language_probability_matches_word_enumeration(model, n): + even_ones = _automaton(DFA, [("a", "a", "0"), ("a", "b", "1"), ("b", "b", "0"), ("b", "a", "1")], "a") + no_one_one = _automaton(DFA, NEVER_ONE_ONE, "ac") + for dfa in (even_ones, no_one_one): + expected = _brute_language_probability(model, dfa.recognizes, n) + assert regular_language_probability(model, dfa, n) == pytest.approx(expected, abs=1e-12) + + +def test_regular_language_probability_with_explicit_start(): + model = _split(0.3) + contains_one = _automaton(DFA, EVENTUALLY_ONE, "b") + for n in range(6): + expected = _brute_language_probability(model, contains_one.recognizes, n, start={"S": 1.0}) + assert regular_language_probability(model, contains_one, n, start="S") == pytest.approx(expected) + + +def test_regular_language_probability_validates_input(): + dfa = _automaton(DFA, ALWAYS_ZERO, "a") + with pytest.raises(ValueError, match="nonnegative"): + regular_language_probability(bernoulli(), dfa, -1) diff --git a/tests/test_process_equivalence.py b/tests/test_process_equivalence.py new file mode 100644 index 0000000..1716e76 --- /dev/null +++ b/tests/test_process_equivalence.py @@ -0,0 +1,8 @@ +def test_is_equal_process_accepts_sympy_probabilities(): + """Regression: object-dtype symbol matrices broke ``matrix_rank``.""" + import sympy as sp + + from sofic.examples import golden_mean + + assert golden_mean(sp.Rational(1, 2)).is_equal_process(golden_mean(0.5)) + assert not golden_mean(sp.Rational(1, 3)).is_equal_process(golden_mean(0.5)) diff --git a/tests/test_renyi.py b/tests/test_renyi.py new file mode 100644 index 0000000..0140cc9 --- /dev/null +++ b/tests/test_renyi.py @@ -0,0 +1,155 @@ +"""Rényi entropy rates, pressure and the large-deviation rate function.""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.examples import bernoulli, even_process, golden_mean +from sofic.examples.epsilon_machines import from_symbol_matrices +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.matrices import symbol_matrices +from sofic.generators.mealy import MealyHMM +from sofic.generators.renyi import pressure, rate_function, renyi_entropy_rate +from sofic.testing.strategies import epsilon_machines + +TOL = 1e-9 + + +@st.composite +def weighted_machines(draw, max_states: int = 3) -> EpsilonMachine: + """Topological ε-machine skeleton with random positive edge probabilities.""" + skeleton = draw(epsilon_machines(max_states=max_states)) + states = list(skeleton.reindex().states) + matrices = {symbol: np.array(m, dtype=float) for symbol, m in symbol_matrices(skeleton).items()} + for i in range(len(states)): + edges = [(symbol, j) for symbol, m in matrices.items() for j in np.flatnonzero(m[i])] + weights = np.array([draw(st.floats(0.1, 1.0)) for _ in edges]) + for (symbol, j), weight in zip(edges, weights / weights.sum(), strict=True): + matrices[symbol][i, j] = weight + return from_symbol_matrices(states, tuple(matrices), matrices) + + +def _block_power_sum(model, n: int, alpha: float) -> float: + probs = np.array([p for p in model.words_of_length(n).values() if p > 0.0]) + return float(np.log2(np.sum(probs**alpha))) + + +def _support_entropy(model) -> float: + adjacency = sum((np.asarray(m, dtype=float) > 0).astype(float) for m in symbol_matrices(model).values()) + return float(np.log2(np.max(np.abs(np.linalg.eigvals(adjacency))))) + + +def _noisy_hmm() -> MealyHMM: + hmm = MealyHMM(observation_alphabet=frozenset({"0", "1"}), initial_distribution={"A": 0.5, "B": 0.5}) + for state in ("A", "B"): + hmm.graph.add_state(state) + hmm.add_transition("A", "A", "0", 0.5) + hmm.add_transition("A", "B", "1", 0.5) + hmm.add_transition("B", "A", "1", 0.5) + hmm.add_transition("B", "B", "1", 0.5) + return hmm + + +@pytest.mark.parametrize("p", [0.1, 0.25, 0.5]) +@pytest.mark.parametrize("alpha", [0.0, 0.5, 2.0, 7.0]) +def test_iid_closed_form(p: float, alpha: float) -> None: + expected = math.log2((1 - p) ** alpha + p**alpha) / (1 - alpha) + assert renyi_entropy_rate(bernoulli(p), alpha) == pytest.approx(expected, abs=TOL) + assert pressure(bernoulli(p), alpha) == pytest.approx((1 - alpha) * expected, abs=TOL) + + +def test_iid_min_entropy() -> None: + assert renyi_entropy_rate(bernoulli(0.2), np.inf) == pytest.approx(-math.log2(0.8), abs=TOL) + + +def test_golden_mean_limits() -> None: + model = golden_mean(0.5) + assert renyi_entropy_rate(model, 0.0) == pytest.approx(math.log2((1 + math.sqrt(5)) / 2), abs=TOL) + assert renyi_entropy_rate(model, 1.0) == pytest.approx(2 / 3, abs=TOL) + assert renyi_entropy_rate(model, np.inf) == pytest.approx(0.5, abs=TOL) + + +def test_non_unifilar_input_is_converted() -> None: + hmm = _noisy_hmm() + eps = EpsilonMachine.from_hmm(hmm) + for alpha in (0.0, 0.5, 1.0, 2.0, np.inf): + assert renyi_entropy_rate(hmm, alpha) == pytest.approx(renyi_entropy_rate(eps, alpha), abs=1e-9) + + +@pytest.mark.parametrize("bad", [-0.5, float("nan")]) +def test_rejects_negative_orders(bad: float) -> None: + with pytest.raises(ValueError, match="non-negative"): + renyi_entropy_rate(golden_mean(0.5), bad) + with pytest.raises(ValueError, match="non-negative"): + pressure(golden_mean(0.5), bad) + + +def test_brute_force_even_process() -> None: + model = even_process(0.5) + for alpha in (0.5, 2.0): + estimate = (_block_power_sum(model, 20, alpha) - _block_power_sum(model, 10, alpha)) / (10 * (1 - alpha)) + assert estimate == pytest.approx(renyi_entropy_rate(model, alpha), abs=0.01) + + +@pytest.mark.parametrize("p", [0.2, 0.4]) +def test_iid_rate_function_is_cramer(p: float) -> None: + low, high = -math.log2(1 - p), -math.log2(p) + values = np.linspace(low + 0.05, (low + high) / 2, 5) + q = (values - low) / (high - low) + expected = q * np.log2(q / p) + (1 - q) * np.log2((1 - q) / (1 - p)) + np.testing.assert_allclose(rate_function(bernoulli(p), values), expected, atol=1e-7) + + +def test_rate_function_domain() -> None: + model = golden_mean(0.5) + below, above = rate_function(model, [0.4, 5.0]) + assert below == np.inf + assert np.isnan(above) + assert rate_function(model, 2 / 3) == pytest.approx(0.0, abs=1e-9) + + +@settings(max_examples=40) +@given(weighted_machines(), st.floats(0.05, 5.0)) +def test_properties(model: EpsilonMachine, alpha: float) -> None: + h_mu = float(model.entropy_rate()) + assert renyi_entropy_rate(model, 1.0) == pytest.approx(h_mu, abs=TOL) + assert renyi_entropy_rate(model, 1.0 + 1e-6) == pytest.approx(h_mu, abs=1e-5) + assert renyi_entropy_rate(model, 0.0) == pytest.approx(_support_entropy(model), abs=TOL) + rates = [renyi_entropy_rate(model, a) for a in (0.0, alpha, alpha + 0.5, 2 * alpha + 1.0, np.inf)] + assert all(a >= b - 1e-9 for a, b in zip(rates, rates[1:], strict=False)) + + +@settings(max_examples=20) +@given(weighted_machines(), st.sampled_from([0.5, 2.0, 3.0]), st.integers(1, 10)) +def test_brute_force_sandwich(model: EpsilonMachine, alpha: float, n: int) -> None: + pi = np.asarray(model.stationary_distribution(), dtype=float) + keep = pi > 1e-12 + powered = sum( + np.where(m > 0, m, 1.0) ** alpha * (m > 0) + for m in (np.asarray(m, dtype=float)[np.ix_(keep, keep)] for m in symbol_matrices(model).values()) + ) + middle = math.log2(pi[keep] ** alpha @ np.linalg.matrix_power(powered, n) @ np.ones(keep.sum())) + size = math.log2(keep.sum()) + brute = _block_power_sum(model, n, alpha) + assert middle - size - 1e-9 <= brute <= middle + alpha * size + 1e-9 + assert renyi_entropy_rate(model, alpha) == pytest.approx( + math.log2(np.max(np.abs(np.linalg.eigvals(powered)))) / (1 - alpha), abs=1e-9 + ) + + +@settings(max_examples=20) +@given(weighted_machines()) +def test_rate_function_nonnegative_with_zero_at_entropy_rate(model: EpsilonMachine) -> None: + h_mu = float(model.entropy_rate()) + h_inf = renyi_entropy_rate(model, np.inf) + h_0 = renyi_entropy_rate(model, 0.0) + values = np.linspace(h_inf, h_0, 9) + rates = rate_function(model, values) + finite = rates[np.isfinite(rates)] + assert np.all(finite >= 0.0) + assert rate_function(model, h_mu) == pytest.approx(0.0, abs=1e-7) diff --git a/tests/test_support.py b/tests/test_support.py new file mode 100644 index 0000000..12cda5e --- /dev/null +++ b/tests/test_support.py @@ -0,0 +1,125 @@ +"""Stationary support languages and the decisions built on them.""" + +from __future__ import annotations + +import math +from itertools import product + +import pytest +from hypothesis import given, settings + +from sofic.examples import bernoulli, binary_markov_chain, even_process, golden_mean +from sofic.generators.mealy import MealyHMM +from sofic.generators.relative_entropy_rate import relative_entropy_rate +from sofic.generators.support import is_absolutely_continuous, support_equal, support_includes, support_nfa +from sofic.testing.strategies import mealy_hmms + +POSITIVE = 1e-12 + + +def _transient_prefix() -> MealyHMM: + """A transient state emits one ``1`` and then the chain emits ``0`` forever.""" + hmm = MealyHMM(observation_alphabet=frozenset({"0", "1"}), initial_distribution={"T": 1.0}) + hmm.add_transition("T", "A", "1", 1.0) + hmm.add_transition("A", "A", "0", 1.0) + return hmm + + +def _period_two() -> MealyHMM: + hmm = MealyHMM(observation_alphabet=frozenset({"0", "1"}), initial_distribution={"A": 0.5, "B": 0.5}) + hmm.add_transition("A", "B", "0", 1.0) + hmm.add_transition("B", "A", "1", 1.0) + return hmm + + +MODELS = { + "iid": bernoulli(0.3), + "golden": golden_mean(0.4), + "golden_parry": golden_mean(2 / (1 + math.sqrt(5))), + "even": even_process(0.6), + "markov": binary_markov_chain(0.2, 0.7), + "period_two": _period_two(), + "transient": _transient_prefix(), +} + + +def _positive_words(model, max_length: int) -> set[tuple[str, ...]]: + pi = model.to_mealy().stationary_distribution() + alphabet = sorted(model.observation_alphabet) + return { + word + for length in range(max_length + 1) + for word in product(alphabet, repeat=length) + if model.word_probability(word, start=pi) > POSITIVE + } + + +def _brute_includes(p, q, max_length: int) -> bool: + return _positive_words(p, max_length) <= _positive_words(q, max_length) + + +@pytest.mark.parametrize("name", sorted(MODELS)) +def test_support_nfa_accepts_exactly_the_positive_words(name): + model = MODELS[name] + nfa = support_nfa(model) + positive = _positive_words(model, 6) + for length in range(7): + for word in product(("0", "1"), repeat=length): + assert nfa.recognizes(word) == (word in positive) + + +def test_support_nfa_drops_words_only_transient_states_emit(): + model = _transient_prefix() + assert model.to_support_nfa().recognizes(("1",)) + assert not support_nfa(model).recognizes(("1",)) + assert support_nfa(model).recognizes(("0", "0", "0")) + + +@pytest.mark.parametrize("p", sorted(MODELS)) +@pytest.mark.parametrize("q", sorted(MODELS)) +def test_support_includes_matches_word_enumeration(p, q): + expected = _brute_includes(MODELS[p], MODELS[q], 8) + assert support_includes(MODELS[p], MODELS[q]) == expected + assert is_absolutely_continuous(MODELS[p], MODELS[q]) == expected + + +@pytest.mark.filterwarnings("ignore:stationary distribution is not unique") +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(max_states=2), mealy_hmms(max_states=2)) +def test_support_includes_matches_enumeration_on_random_hmms(p, q): + # A shortest word in L(p) \ L(q) has length below |p states| * 2^|q states| = 8. + assert support_includes(p, q) == _brute_includes(p, q, 8) + + +def test_golden_mean_is_below_the_full_shift(): + assert support_includes(golden_mean(0.5), bernoulli()) + assert not support_includes(bernoulli(), golden_mean(0.5)) + assert golden_mean(0.5).is_absolutely_continuous(bernoulli()) + assert not bernoulli().support_includes(golden_mean(0.5)) + + +def test_support_equal_ignores_probabilities(): + assert support_equal(golden_mean(0.3), golden_mean(0.8)) + assert support_equal(bernoulli(0.1), binary_markov_chain(0.2, 0.7)) + assert not support_equal(golden_mean(0.5), bernoulli()) + assert golden_mean(0.3).support_equal(golden_mean(0.6)) + + +@pytest.mark.parametrize("name", sorted(MODELS)) +def test_support_includes_is_reflexive(name): + assert support_includes(MODELS[name], MODELS[name]) + assert support_equal(MODELS[name], MODELS[name]) + + +def test_support_includes_is_transitive(): + names = sorted(MODELS) + included = {(a, b): support_includes(MODELS[a], MODELS[b]) for a in names for b in names} + for a, b, c in product(names, repeat=3): + if included[a, b] and included[b, c]: + assert included[a, c] + + +@pytest.mark.parametrize("p", ["iid", "golden", "even", "markov", "period_two"]) +@pytest.mark.parametrize("q", ["iid", "golden", "golden_parry", "markov"]) +def test_finite_relative_entropy_rate_iff_support_included(p, q): + assert math.isfinite(relative_entropy_rate(MODELS[p], MODELS[q])) == support_includes(MODELS[p], MODELS[q]) From 8f84c34275a49351bc59d3d15b2afdfa84ca21af Mon Sep 17 00:00:00 2001 From: Ryan James Date: Wed, 7 Oct 2026 07:46:17 -0600 Subject: [PATCH 2/3] CHANGELOG: finalize 0.4.0 (tag left to the maintainer) Co-authored-by: Cursor --- CHANGELOG.md | 144 ++++++++++++++++++++++++++++++--------------------- 1 file changed, 86 insertions(+), 58 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 22a5c6b..b3dd0bd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,9 +1,11 @@ # Changelog -## Unreleased (0.4.0) +## 0.4.0 (unreleased; tag pending) -A second correctness review: every fix below was reproduced against a -brute-force reference or a closed-form value and has a regression test. +A second correctness review, a property and metamorphic test suite, the +literature gaps it surfaced, and ports of automata and shift constructions to +processes. Every fix below was reproduced against a brute-force reference or a +closed-form value and has a regression test. ### Breaking changes @@ -45,6 +47,17 @@ brute-force reference or a closed-form value and has a regression test. | `period7()` | `period8()` (they used the same word) | | `tetris_tgm()` | `tetris_history()` | +- `HiddenMarkovModel.to_support_dfa()` accepts every word of positive + probability, the empty word included (it accepted only words ending in a + terminal recurrent subset). +- `HiddenMarkovModel.entropy_rate()` no longer raises `NotImplementedError` for + non-unifilar presentations; see `method=` below. +- `EpsilonMachine.from_hmm` takes a keyword-only `max_states` and rejects other + keyword arguments (they were silently ignored). Numeric non-unifilar input can + now yield fewer causal states, because round-off no longer splits them. +- `automaton_to_regex` quotes multi-character symbols and escapes `ε`, `∅` and + `'`. + ### Fixes - **Shifts:** `SoficShift.topological_entropy` counts words on a @@ -106,38 +119,42 @@ brute-force reference or a closed-form value and has a regression test. attributes; Graphviz/TikZ node names are unique; TikZ escapes `^`, `~`, `\` for math mode and handles newlines in labels; YAML uses libyaml's C loader and dumper when available. -- `entropy_rate_hmm` and `state_distribution` no longer crash on presentations - whose states are `MixedState` objects. -- The `riechers2017spectral2` citation pointed at an unrelated Phys. Rev. E - paper; it is Chaos 28, 033116 (2018). -- `automaton_to_regex` was ambiguous for multi-character symbols (`ab` vs `a`, - `b`); such symbols are now quoted (`'ab'`), and `ε`, `∅` and `'` are escaped. -- `HiddenMarkovModel.to_support_dfa` accepted only terminal recurrent subsets, - so it rejected words of positive probability (e.g. `1` for the even process); - every non-empty subset now accepts. -- `entropy_rate` and `state_entropy` no longer raise `TypeError` when state - labels mix types (e.g. `0` and `(1, 0)` after a split). -- `is_equal_process` works on sympy-valued HMMs (object-dtype matrices broke - `matrix_rank`). -- `DFA.add_transition` adds a missing source state, as `NFA.add_transition` - does. -- `EpsilonMachine.from_hmm` split causal states that differed only by - floating-point round-off from belief updates (a fair-coin HMM gave 33 states - and C_mu ≈ 2.46); numeric transition probabilities are now compared within the - mixed-state tolerance. +- **More generators and automata:** `EpsilonMachine.from_hmm` no longer splits causal + states that differ only by belief-update round-off (a fair-coin HMM gave 33 + states and C_mu ≈ 2.46); `to_support_dfa` no longer rejects words of positive + probability (e.g. `1` for the even process); `entropy_rate` and + `state_entropy` handle mixed state-label types (e.g. `0` and `(1, 0)` after a + split) and `MixedState` labels; `is_equal_process` works on sympy-valued + HMMs; `DFA.add_transition` adds a missing source state, as + `NFA.add_transition` does; `automaton_to_regex` was ambiguous for + multi-character symbols (`ab` vs `a`, `b`). +- **Docs:** the `riechers2017spectral2` citation pointed at an unrelated + Phys. Rev. E paper; it is Chaos 28, 033116 (2018). ### New +**Testing** + - 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. + +**Shifts** + +- `periodic_points(n)` and `zeta_function()` for topological Markov chains, + shifts of finite type and sofic shifts (Manning's signed-subset formula for + sofic shifts; Lind & Marcus §6.4); counts are exact integers, and + `zeta_function` needs the `symbolic` extra. +- In- and out-state splitting and amalgamation of topological Markov chains with + division and edge matrices (`A = DE`, `A' = ED`; Lind & Marcus §2.4); + `bowen_franks_group` (with the sign of `det(I − A)`) and + `jordan_form_away_from_zero` as conjugacy invariants (§7.4). + +**Automata** + - `parse_regex`, `regex_to_nfa` and `NFA.from_regex` build an NFA from a regular expression by Thompson's construction (Thompson 1968); malformed input raises `sofic.exceptions.RegexSyntaxError`. @@ -149,18 +166,34 @@ brute-force reference or a closed-form value and has a regression test. overrides them with ω-semantics: `is_empty` and the new `accepted_lasso` decide Büchi emptiness; universality and inclusion raise `NotImplementedError`. + +**Entropy rates and divergences** + +- `HiddenMarkovModel.entropy_rate(method="auto"|"exact"|"bounds"|"blackwell")` + handles non-unifilar presentations: exact when the mixed states close, + otherwise the converged Cover & Thomas (Thm 4.5.1) bounds, with a warning if + they have not met `tol`; `auto` never returns a stochastic estimate. New + `entropy_rate_bounds`, `entropy_rate_blackwell` (batch-means standard error) + and `mixed_state_walk`. +- `statistical_complexity_dimension` and `ifs_lyapunov_dimension` (Jurgens & + Crutchfield 2021): the Lyapunov dimension of the Blackwell measure from the + mixed-state random walk; reproduces the paper's Cantor (log 2 / log 3) and + Sierpinski (log₂ 3) examples. - `sofic.generators.relative_entropy_rate`: `relative_entropy_rate(p, q)` in bits (exact for unifilar `q`, `inf` when P is not absolutely continuous with respect to Q on finite words), `relative_entropy_rate_bounds(p, q, n)` for general HMMs, and `HiddenMarkovModel.relative_entropy_rate`. -- `periodic_points(n)` and `zeta_function()` for topological Markov chains, - shifts of finite type and sofic shifts (Manning's signed-subset formula for - sofic shifts; Lind & Marcus §6.4); counts are exact integers, and - `zeta_function` needs the `symbolic` extra. -- In- and out-state splitting and amalgamation of topological Markov chains with - division and edge matrices (`A = DE`, `A' = ED`; Lind & Marcus §2.4); - `bowen_franks_group` (with the sign of `det(I − A)`) and - `jordan_form_away_from_zero` as conjugacy invariants (§7.4). +- `renyi_entropy_rate`, `pressure` and `rate_function` + (`sofic.generators.renyi`): Rényi entropy rates for α ∈ [0, ∞] (topological, + Shannon and min-entropy rates as special cases) and the large-deviation rate + function of −(1/n) log₂ P(X_{0:n}) by Legendre transform of the pressure. +- `sofic.generators.support`: `support_nfa`, `support_includes`, + `is_absolutely_continuous`, `support_equal` (stationary support; + finite-cylinder absolute continuity), with matching `HiddenMarkovModel` + methods. + +**Spectral and predictive structure** + - `sofic.generators.correlations`: closed-form `autocorrelation`, `power_spectrum` (continuous part; delta peaks at unit-circle eigenvalues excluded) and `mutual_information_function` for any finite HMM (Riechers & @@ -169,18 +202,9 @@ brute-force reference or a closed-form value and has a regression test. information bottleneck / predictive rate-distortion curve of a finite ε-machine, by an annealed, seeded Blahut–Arimoto iteration over exact future morphs (Still et al. 2010; Marzen & Crutchfield 2016; Tishby et al. 2000). -- `HiddenMarkovModel.entropy_rate(method="auto"|"exact"|"bounds"|"blackwell")` - handles non-unifilar presentations: exact when the mixed states close, - otherwise the converged Cover & Thomas (Thm 4.5.1) bounds, with a warning if - they have not met `tol`; `auto` never returns a stochastic estimate. New - `entropy_rate_bounds`, `entropy_rate_blackwell` (batch-means standard error) - and `mixed_state_walk`. -- `statistical_complexity_dimension` and `ifs_lyapunov_dimension` (Jurgens & - Crutchfield 2021): the Lyapunov dimension of the Blackwell measure from the - mixed-state random walk; reproduces the paper's Cantor (log 2 / log 3) and - Sierpinski (log₂ 3) examples. -- `EpsilonMachine.from_hmm(hmm, max_states=...)` caps the mixed-state - enumeration (it silently ignored keyword arguments). + +**Presentations and constructions** + - `minimal_quasi_realization` and `process_rank` (plus `HiddenMarkovModel.process_rank`): Schützenberger/Fliess minimization of an HMM's linear representation to a minimal-dimension `QuasiRealization` whose @@ -195,29 +219,33 @@ brute-force reference or a closed-form value and has a regression test. conjugacy invariants (`E(β_k X) = E(X) + (k − 1) h_μ`); only `h_μ` is. - `sofic.generators.state_splitting`: process-preserving out- and in-splitting (`split_state`) and `amalgamate`. -- `sofic.generators.support`: `support_nfa`, `support_includes`, - `is_absolutely_continuous`, `support_equal` (stationary support; - finite-cylinder absolute continuity), with matching `HiddenMarkovModel` - methods. - `condition_on_language`: conditions an HMM on a prefix-closed regular constraint by the Doob h-transform; the uniform i.i.d. process conditioned on an irreducible SFT is its Parry measure. - `omega_probability` for deterministic Büchi properties of HMM output (bottom strongly connected components of the product chain; Baier & Katoen 2008) and `regular_language_probability` for exact `P(X_{0:n} ∈ L)`. -- `sofic.inference.learn_epsilon_machine_active` and `ProcessOracle`: an - L*-style active learner for ε-machines from probability and equivalence - queries; exact for finite, exactly synchronizable ε-machines, with `history` - seeding for infinite transients. -- `renyi_entropy_rate`, `pressure` and `rate_function` - (`sofic.generators.renyi`): Rényi entropy rates for α ∈ [0, ∞] (topological, - Shannon and min-entropy rates as special cases) and the large-deviation rate - function of −(1/n) log₂ P(X_{0:n}) by Legendre transform of the pressure. +- `EpsilonMachine.from_hmm(hmm, max_states=...)` caps the mixed-state + enumeration (it silently ignored keyword arguments). - Experimental `canonical_residual_hmm` (`sofic.generators.canonical_residual`, requires `experimental=True`): a generator whose states are the extreme future morphs of a finite ε-machine, the process analogue of the canonical RFSA; it can be strictly smaller than the ε-machine. +**Time reversal** + +- 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. + +**Inference** + +- `sofic.inference.learn_epsilon_machine_active` and `ProcessOracle`: an + L*-style active learner for ε-machines from probability and equivalence + queries; exact for finite, exactly synchronizable ε-machines, with `history` + seeding for infinite transients. + ## 0.3.0 A correctness review and breaking refactor. Old names are **not** kept as From 38f22e74e3e01fb501f784f9eb4a5181abfb74b6 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Wed, 7 Oct 2026 10:10:37 -0600 Subject: [PATCH 3/3] Tests: skip sympy-dependent cases when sympy is not installed (CI installs only the test extra) Co-authored-by: Cursor --- tests/test_active_epsilon_learning.py | 7 +++++-- tests/test_bisimulation.py | 2 +- tests/test_hmm_state_splitting.py | 5 +++-- tests/test_minimal_quasi_realization.py | 7 +++++-- tests/test_process_equivalence.py | 4 +++- 5 files changed, 17 insertions(+), 8 deletions(-) diff --git a/tests/test_active_epsilon_learning.py b/tests/test_active_epsilon_learning.py index a0cead3..cb42182 100644 --- a/tests/test_active_epsilon_learning.py +++ b/tests/test_active_epsilon_learning.py @@ -6,7 +6,6 @@ import numpy as np import pytest -import sympy as sp from hypothesis import given, settings from hypothesis import strategies as st @@ -19,7 +18,10 @@ from sofic.inference.active import ProbabilityOracle, ProcessOracle, learn_epsilon_machine_active from sofic.testing.strategies import epsilon_machines -R = sp.Rational + +def R(numerator, denominator): + """Exact rational; skips the calling test when sympy is not installed.""" + return pytest.importorskip("sympy").Rational(numerator, denominator) def _n_states(machine): @@ -180,6 +182,7 @@ def test_recovers_random_strategy_epsilon_machines(topology, seed): def test_exact_sympy_oracle_gives_exact_machine(): + sp = pytest.importorskip("sympy") target = _exact_even() learned = learn_epsilon_machine_active(ProcessOracle(target), ("0", "1")) probs = sorted((t.data[ATTR_EMISSION], t.data[ATTR_PROB]) for t in learned.transitions()) diff --git a/tests/test_bisimulation.py b/tests/test_bisimulation.py index 23692df..ff89891 100644 --- a/tests/test_bisimulation.py +++ b/tests/test_bisimulation.py @@ -3,7 +3,6 @@ from fractions import Fraction import pytest -import sympy as sp from hypothesis import given, settings from hypothesis import strategies as st @@ -167,6 +166,7 @@ def test_markov_chain_quotient_type(): def test_symbolic_probabilities_compare_exactly(): + sp = pytest.importorskip("sympy") a = sp.Symbol("a", positive=True) hmm = MealyHMM(initial_distribution={0: sp.Integer(1)}) hmm.add_transition(0, 1, 0, a / 2) diff --git a/tests/test_hmm_state_splitting.py b/tests/test_hmm_state_splitting.py index 2c9d3c2..01b7ded 100644 --- a/tests/test_hmm_state_splitting.py +++ b/tests/test_hmm_state_splitting.py @@ -11,7 +11,6 @@ from typing import Any import pytest -import sympy from hypothesis import assume, given, settings from hypothesis import strategies as st @@ -193,6 +192,7 @@ def test_in_split_preserves_unifilarity(): def test_out_split_is_exact_for_symbolic_probabilities(): + sympy = pytest.importorskip("sympy") p = sympy.Symbol("p", positive=True) hmm = MealyHMM(initial_distribution={"A": 1 / (1 + p), "B": p / (1 + p)}) hmm.graph.add_state("A") @@ -212,7 +212,8 @@ def test_out_split_is_exact_for_symbolic_probabilities(): assert numeric_split.is_equal_process(numeric_hmm) -def substitute(hmm: MealyHMM, symbol: sympy.Symbol, value: Any) -> MealyHMM: +def substitute(hmm: MealyHMM, symbol: Any, value: Any) -> MealyHMM: + sympy = pytest.importorskip("sympy") numeric = MealyHMM( observation_alphabet=hmm.observation_alphabet, initial_distribution={ diff --git a/tests/test_minimal_quasi_realization.py b/tests/test_minimal_quasi_realization.py index 039fee6..a17d2ce 100644 --- a/tests/test_minimal_quasi_realization.py +++ b/tests/test_minimal_quasi_realization.py @@ -1,10 +1,10 @@ """Tests for Schützenberger/Fliess minimal quasi-realizations and the process rank.""" from functools import cache +from typing import Any import numpy as np import pytest -import sympy as sp from hypothesis import assume, given, settings from hypothesis import strategies as st @@ -144,6 +144,7 @@ def test_returns_quasi_realization_with_unit_final_vector(): def test_exact_rational_input_stays_exact(): + sp = pytest.importorskip("sympy") third = sp.Rational(1, 3) hmm = MealyHMM(initial_distribution={0: sp.Rational(3, 4), 1: sp.Rational(1, 4)}) hmm.add_transition(0, 0, 0, third) @@ -160,6 +161,7 @@ def test_exact_rational_input_stays_exact(): def test_symbolic_parameter(): + sp = pytest.importorskip("sympy") a = sp.Symbol("a", positive=True) hmm = MealyHMM(initial_distribution={0: sp.Integer(1)}) hmm.add_transition(0, 0, 0, a) @@ -174,7 +176,8 @@ def test_symbolic_parameter(): assert sp.simplify(result @ realization.tau - a * (1 - a)) == 0 -def _exact_word_probability(hmm: MealyHMM, word: tuple[int, ...]) -> sp.Expr: +def _exact_word_probability(hmm: MealyHMM, word: tuple[int, ...]) -> Any: + sp = pytest.importorskip("sympy") states = sorted(hmm.states()) vector = {s: sp.sympify(hmm.initial_distribution.get(s, 0)) for s in states} for symbol in word: diff --git a/tests/test_process_equivalence.py b/tests/test_process_equivalence.py index 1716e76..b5e2715 100644 --- a/tests/test_process_equivalence.py +++ b/tests/test_process_equivalence.py @@ -1,6 +1,8 @@ def test_is_equal_process_accepts_sympy_probabilities(): """Regression: object-dtype symbol matrices broke ``matrix_rank``.""" - import sympy as sp + import pytest + + sp = pytest.importorskip("sympy") from sofic.examples import golden_mean