Automata and shift ideas ported to processes (4/4)

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,25 +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.
+- **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`.
@@ -136,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 &
@@ -156,18 +202,49 @@ 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.
+
+**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
+  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`.
+- `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)`.
 - `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
 

docs/generators/bisimulation.rst

@@ -0,0 +1,69 @@
+.. bisimulation.rst
+
+**************************
+Probabilistic Bisimulation
+**************************
+
+:doc:`Lumping <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
+<sofic.generators.epsilon_machine.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

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