diff --git a/CHANGELOG.md b/CHANGELOG.md index 7cf0e08..78d5b4b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -106,6 +106,12 @@ 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. ### New @@ -119,6 +125,49 @@ brute-force reference or a closed-form value and has a regression test. - 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. +- `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`. +- `learn_dfa_k_testable`, the García–Vidal learner for strictly k-testable + languages: the smallest such language containing a set of positive samples. +- `LabeledAutomaton.is_empty`, `accepted_word` (a shortest witness), + `is_universal` and `includes` for NFAs and DFAs, by the antichain algorithms + of De Wulf et al. (2006) without full determinization. `BuchiAutomaton` + overrides them with ω-semantics: `is_empty` and the new `accepted_lasso` + decide Büchi emptiness; universality and inclusion raise + `NotImplementedError`. +- `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). +- `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 & + Crutchfield 2018). +- `predictive_rate_distortion` and `PredictiveRateDistortionCurve`: the causal + 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). ## 0.3.0 diff --git a/docs/automata/algorithms.rst b/docs/automata/algorithms.rst index bbabcdb..a48d78e 100644 --- a/docs/automata/algorithms.rst +++ b/docs/automata/algorithms.rst @@ -33,6 +33,39 @@ All three minimizers return the same trimmed (partial) minimal DFA; use In [6]: equivalent(dfa, trimmed, frozenset({0})) Out[6]: True +Regular expressions +=================== + +:func:`~sofic.automata.regex.automaton_to_regex` converts an automaton to a +regular expression by state elimination :cite:`Kleene1956,HopcroftUllman1979`, +and :func:`~sofic.automata.regex.regex_to_nfa` (also +:meth:`NFA.from_regex `) converts back with +Thompson's construction :cite:`Thompson1968`, producing an ε-NFA with one +initial and one accepting state. The two read and write the same dialect, which +is Python :mod:`re` syntax whenever every symbol is a single-character string: + +* a single character is a symbol; ``\c`` escapes any character ``c``; +* ``'ab'`` quotes a multi-character symbol (``\'`` and ``\\`` escape inside the + quotes), so ``'ab'`` is one symbol while ``ab`` is two; +* ``ε`` or ``(?:)`` is the empty word, ``∅`` or ``(?!)`` the empty language; +* ``|``, juxtaposition, and postfix ``*``, ``+``, ``?`` are alternation, + concatenation, and repetition; ``(...)`` and ``(?:...)`` group; +* ``. [ ] { } ^ $`` are reserved and must be escaped. + +Malformed input (unbalanced parentheses, dangling operators, empty +alternatives) raises :class:`~sofic.exceptions.RegexSyntaxError`, a +:class:`ValueError`. Passing ``alphabet`` resolves each parsed symbol to the +alphabet element with the same ``str``, so automata over non-string symbols +round-trip: + +.. code-block:: python + + from sofic.automata import NFA, automaton_to_regex, equivalent, regex_to_nfa + + nfa = NFA.from_regex("'ab'(a|b)*") + assert nfa.recognizes(("ab", "a", "b")) + assert equivalent(nfa, regex_to_nfa(automaton_to_regex(nfa), alphabet=nfa.input_alphabet)) + API === @@ -47,5 +80,10 @@ API .. autofunction:: nerode_partition .. autofunction:: equivalent .. autofunction:: sofic.automata.regex.automaton_to_regex +.. autofunction:: sofic.automata.regex.parse_regex +.. autofunction:: sofic.automata.regex.regex_to_nfa +.. autoclass:: sofic.automata.regex.Regex + :members: pattern +.. autoexception:: sofic.exceptions.RegexSyntaxError .. autodata:: MinimizationAlgorithm diff --git a/docs/automata/buchi.rst b/docs/automata/buchi.rst index fb85a39..935daa7 100644 --- a/docs/automata/buchi.rst +++ b/docs/automata/buchi.rst @@ -12,7 +12,7 @@ API === .. autoclass:: BuchiAutomaton - :members: accepts_lasso, accepts_omega + :members: accepts_lasso, accepts_omega, accepted_lasso, is_empty .. autofunction:: sofic.automata.buchi_simulation.accepts_lasso_buchi .. autofunction:: sofic.automata.buchi_simulation.accepts_omega_buchi diff --git a/docs/automata/dfa.rst b/docs/automata/dfa.rst index 70b4462..4ded56e 100644 --- a/docs/automata/dfa.rst +++ b/docs/automata/dfa.rst @@ -35,8 +35,21 @@ and no ε-transitions, following the standard finite-automata model In [8]: dfa.recognizes(()) Out[8]: False + @doctest + In [9]: dfa.accepted_word() + Out[9]: (0,) + + @doctest + In [10]: dfa.is_universal() + Out[10]: False + +The emptiness, shortest-witness, universality, and inclusion checks are shared +with :class:`~sofic.automata.nfa.NFA`; see that page for the antichain +algorithms :cite:`DeWulf2006`. + API === .. autoclass:: DFA - :members: add_transition, recognizes, union, intersection, complement, difference, concat, kleene_star, minimize, from_nfa + :members: add_transition, recognizes, union, intersection, complement, difference, concat, kleene_star, minimize, from_nfa, + is_empty, accepted_word, is_universal, includes diff --git a/docs/automata/learning.rst b/docs/automata/learning.rst index 3465398..a840e51 100644 --- a/docs/automata/learning.rst +++ b/docs/automata/learning.rst @@ -138,6 +138,29 @@ satisfiable ``k``. It requires the optional `python-sat .. autofunction:: sofic.automata.learning.dfasat.learn_dfa_sat +Passive learning from positive data (k-testable) +================================================ + +A language is *strictly k-testable* (k-TSS) when membership depends only on a +word's length-``k - 1`` prefix and suffix, the set of its length-``k`` factors, +and, for words shorter than ``k - 1``, the word itself. García & Vidal's k-TSSI +algorithm :cite:`GarciaVidal1990` collects those prefixes, suffixes, factors, +and short words from **positive samples only** and returns a DFA for the +smallest k-TSS language containing the sample. Its states are the sample +prefixes shorter than ``k - 1`` and the reachable ``(k - 1)``-windows. The +learner identifies every k-TSS language in the limit from positive data, and a +larger ``k`` gives a (weakly) smaller language: + +.. code-block:: python + + from sofic.automata import learn_dfa_k_testable + + dfa = learn_dfa_k_testable(["", "ab", "abab"], k=2) # (ab)* + assert dfa.recognizes("ababab") + assert not dfa.recognizes("aab") + +.. autofunction:: sofic.automata.learning.k_testable.learn_dfa_k_testable + Probabilistic passive learning (ALERGIA) ======================================== diff --git a/docs/automata/nfa.rst b/docs/automata/nfa.rst index 64983e6..4686d09 100644 --- a/docs/automata/nfa.rst +++ b/docs/automata/nfa.rst @@ -27,8 +27,33 @@ languages :cite:`RabinScott1959,HopcroftUllman1979`. In [6]: dfa = determinize(nfa) +Decision procedures +=================== + +``is_empty`` and ``accepted_word`` search the state graph breadth-first, so the +witness is a shortest accepted word (ε-moves cost nothing). ``is_universal`` +and ``includes`` use the forward antichain algorithms of De Wulf et al. +:cite:`DeWulf2006`, which explore the subset construction on the fly and keep +only ``⊆``-minimal subsets instead of determinizing. Universality is relative +to ``input_alphabet`` by default, matching :meth:`~NFA.complement`. + +.. ipython:: + + @doctest + In [7]: nfa.accepted_word() + Out[7]: (0,) + + @doctest + In [8]: nfa.is_universal() + Out[8]: False + + @doctest + In [9]: dfa.includes(nfa) + Out[9]: True + API === .. autoclass:: NFA - :members: add_transition, recognizes, union, intersection, complement, difference, concat, kleene_star, determinize, minimize + :members: add_transition, recognizes, union, intersection, complement, difference, concat, kleene_star, determinize, minimize, + is_empty, accepted_word, is_universal, includes diff --git a/docs/generators/correlations.rst b/docs/generators/correlations.rst new file mode 100644 index 0000000..559c9f6 --- /dev/null +++ b/docs/generators/correlations.rst @@ -0,0 +1,104 @@ +.. correlations.rst +.. py:module:: sofic.generators.correlations + +************************************** +Correlations and power spectra +************************************** + +Pairwise statistics of a process generated by any finite HMM presentation +follow in closed form from its symbol-labeled transition matrices +:math:`T^{(x)}`, :math:`T = \sum_x T^{(x)}`, and stationary state law +:math:`\pi` (Riechers & Crutchfield :cite:`riechers2017spectral1`, §V.A). For a +real observable :math:`f` of the symbols set :math:`\Omega = \sum_x f(x) T^{(x)}`. + +Autocorrelation +=============== + +.. math:: + + \gamma(0) = E[f(X_0)^2], \qquad + \gamma(\tau) = E[f(X_0) f(X_\tau)] + = \langle \pi | \Omega \, T^{\tau - 1} \, \Omega | \mathbf{1} \rangle + \quad (\tau \geq 1). + +``centered=True`` returns the autocovariance :math:`\gamma(\tau) - E[f(X_0)]^2`. +The default observable is ``float(symbol)``, so the binary alphabet ``"0"``, +``"1"`` maps to 0, 1; pass a mapping or callable for anything else. + +Power spectrum +============== + +The power spectrum :math:`P(\omega) = \lim_{N\to\infty} \frac{1}{N} +E\,|\sum_{t=0}^{N-1} f(X_t) e^{-i\omega t}|^2 = \sum_{\tau} \gamma(\tau) +e^{-i\omega\tau}` (Wiener–Khinchin) of a finite HMM is a continuous part plus +Dirac deltas. :func:`power_spectrum` returns the continuous part +(:cite:`riechers2017spectral1`, Eq. (5)): + +.. math:: + + P_c(\omega) = E[f(X_0)^2] + + 2\,\mathrm{Re}\,\langle \pi | \Omega \, (e^{i\omega} I - T)^{-1} \, \Omega | \mathbf{1} \rangle . + +The resolvent is singular whenever :math:`e^{i\omega}` is an eigenvalue of +:math:`T` — always at :math:`\omega = 0` — so the implementation projects out +the peripheral (unit-modulus) spectrum with its spectral projector :math:`Q` +and evaluates the equivalent, everywhere-finite form + +.. math:: + + P_c(\omega) = E[f(X_0)^2] - \langle \pi | \Omega (TQ)^{\#} \Omega | \mathbf{1} \rangle + + 2\,\mathrm{Re}\,\langle \pi | \Omega \, (e^{i\omega} I - T(I-Q))^{-1} (I - Q)\, \Omega | \mathbf{1} \rangle . + +.. note:: + + **Delta peaks are excluded.** The DC peak :math:`2\pi E[f]^2\delta(\omega)` + and, for processes with periodic components, peaks at the other roots of + unity in the spectrum of :math:`T` are not part of the returned array. In + particular a purely periodic process has :math:`P_c \equiv 0`, and + :math:`P_c` is unchanged by centering :math:`f`. At :math:`\omega = 0` the + returned value is the (finite) Green–Kubo limit of the continuous part. + +A flat :math:`P_c` does not imply an i.i.d. process: Riechers & Crutchfield +:cite:`riechers2019fraudulent` construct arbitrarily structured "fraudulent +white noise" processes whose pairwise statistics vanish. + +Mutual information function +=========================== + +:func:`mutual_information_function` returns :math:`I[X_0 : X_\tau]` in bits +for :math:`\tau = 0, \ldots, L` (entry 0 is :math:`H[X_0]`) from the exact lagged +joint law + +.. math:: + + P(X_0 = x, X_\tau = y) = \langle \pi | T^{(x)} \, T^{\tau-1} \, T^{(y)} | \mathbf{1} \rangle . + +Unlike :math:`\gamma`, it depends only on which symbols are distinct, not on +their numeric values, and it decays at roughly twice the rate set by the +subdominant eigenvalue of :math:`T` (for the golden mean process, +:math:`I[X_0 : X_{\tau+1}] / I[X_0 : X_\tau] \to \lambda_2^2`). + +Example +======= + +.. code-block:: python + + import numpy as np + from sofic.examples import golden_mean + from sofic.generators.correlations import ( + autocorrelation, + mutual_information_function, + power_spectrum, + ) + + process = golden_mean(0.5) + autocorrelation(process, 5, centered=True) # Var * (p - 1)**tau + power_spectrum(process, np.linspace(0, np.pi, 64)) + mutual_information_function(process, 10) + +API +=== + +.. autofunction:: autocorrelation +.. autofunction:: power_spectrum +.. autofunction:: mutual_information_function diff --git a/docs/generators/generators.rst b/docs/generators/generators.rst index 87b695b..49c24ae 100644 --- a/docs/generators/generators.rst +++ b/docs/generators/generators.rst @@ -32,8 +32,12 @@ Computational mechanics block_convergence bidirectional_epsilon_machine mixed_state_presentation + nonunifilar_entropy_rate edge_machine information_anatomy + relative_entropy_rate + correlations + predictive_rate_distortion symbolic_hmm generative_models directional_flow diff --git a/docs/generators/nonunifilar_entropy_rate.rst b/docs/generators/nonunifilar_entropy_rate.rst new file mode 100644 index 0000000..4b6308b --- /dev/null +++ b/docs/generators/nonunifilar_entropy_rate.rst @@ -0,0 +1,83 @@ +.. nonunifilar_entropy_rate.rst +.. py:module:: sofic.generators.measures + :no-index: + +************************************* +Entropy Rate of Non-Unifilar Models +************************************* + +For a unifilar presentation the entropy rate is the state-averaged +next-symbol entropy. For a non-unifilar hidden Markov model there is no closed +form :cite:`blackwell1957entropy`: the observer's belief about the hidden state +(the mixed state) generically ranges over an infinite set +:cite:`jurgens2021shannon`. :meth:`~sofic.generators.base.HiddenMarkovModel.entropy_rate` +therefore takes a ``method``: + +``"exact"`` + The unifilar closed form, applied to the mixed-state presentation + :cite:`Ellison2009` when the presentation is non-unifilar and its mixed + states close within ``max_states``. +``"bounds"`` + The Cover & Thomas sandwich (Thm. 4.5.1 :cite:`Cover2006`, written with + 0-based indices), + + .. math:: + + H[X_n \mid X_{0:n}, S_0] \;\le\; h_\mu \;\le\; H[X_n \mid X_{0:n}], + + with the lower bound nondecreasing and the upper bound nonincreasing in + :math:`n`, iterated until the gap is below ``tol``. Convergence is often + exponential :cite:`Travers2013`. +``"blackwell"`` + A Monte Carlo average of :math:`H[X \mid \eta_t]` along a sampled + mixed-state trajectory :cite:`blackwell1957entropy,jurgens2021shannon`. +``"auto"`` (default) + ``"exact"`` when the mixed states close, otherwise ``"bounds"`` (warning + if they have not met ``tol``). It never silently returns a stochastic + estimate. + +The simple nonunifilar source has countably many mixed states, so ``"auto"`` +falls back to the bounds: + +.. ipython:: + + In [1]: from sofic.examples import sns + + In [2]: from sofic.generators.measures import entropy_rate_bounds, entropy_rate_blackwell + + In [3]: entropy_rate_bounds(sns(), 10) + + In [4]: sns().entropy_rate() + + In [5]: entropy_rate_blackwell(sns(), n_samples=20_000, seed=0) + +Statistical complexity dimension +================================ + +When the mixed states do not close, the statistical complexity diverges. Its +divergence rate, the statistical complexity dimension :math:`d_\mu`, is the +information dimension of the Blackwell measure on the mixed states +:cite:`jurgens2021divergent`. The mixed-state dynamic is an iterated function +system on the simplex, and :math:`d_\mu` is bounded above by (and, under the +open set condition, equal to) its Lyapunov dimension. That dimension is the +Kaplan–Yorke formula with the entropy rate as the leading exponent. A process +with finitely many causal states has :math:`d_\mu = 0`. + +.. ipython:: + + In [6]: from sofic.generators.complexity_dimension import statistical_complexity_dimension + + In [7]: statistical_complexity_dimension(sns(), n_samples=5_000, seed=0).dimension + +API +=== + +.. autofunction:: sofic.generators.measures.entropy_rate +.. autofunction:: sofic.generators.measures.entropy_rate_bounds +.. autofunction:: sofic.generators.measures.entropy_rate_blackwell +.. autoclass:: sofic.generators.measures.EntropyRateEstimate +.. autofunction:: sofic.generators.measures.mixed_state_walk +.. autofunction:: sofic.generators.measures.batch_means_stderr +.. autofunction:: sofic.generators.complexity_dimension.statistical_complexity_dimension +.. autofunction:: sofic.generators.complexity_dimension.ifs_lyapunov_dimension +.. autoclass:: sofic.generators.complexity_dimension.StatisticalComplexityDimension diff --git a/docs/generators/predictive_rate_distortion.rst b/docs/generators/predictive_rate_distortion.rst new file mode 100644 index 0000000..eea3596 --- /dev/null +++ b/docs/generators/predictive_rate_distortion.rst @@ -0,0 +1,78 @@ +.. predictive_rate_distortion.rst +.. py:module:: sofic.generators.predictive_rd + +************************** +Predictive rate-distortion +************************** + +How much of a process's predictive information survives when the observer may +only store a limited amount of memory? Predictive rate-distortion answers this +with lossy causal states. Optimal causal inference (Still, Crutchfield & +Ellison :cite:`still2010optimal`) and predictive rate-distortion (Marzen & +Crutchfield :cite:`Marzen2016`) compress the past into features :math:`R` via +the information bottleneck (Tishby, Pereira & Bialek +:cite:`tishby2000information`). The causal state :math:`S` is a sufficient +statistic of the past, so it is enough to compress :math:`S`: + +.. math:: + + \min_{q(r \mid s)} \; I[S; R] - \beta \, I[R; X_{0:L}], + \qquad R - S - X_{0:L} . + +The rate :math:`I[S;R]` lies in :math:`[0, C_\mu]` and the retained predictive +information :math:`I[R; X_{0:L}]` in :math:`[0, I[S; X_{0:L}]]`, where +:math:`I[S; X_{0:L}] \uparrow \mathbf{E}` as :math:`L \to \infty`. The +distortion reported is the predictive information lost, +:math:`I[S; X_{0:L}] - I[R; X_{0:L}]`. + +Algorithm +========= + +* **Future morphs** :math:`P(X_{0:L} = w \mid S = s) = \langle \delta_s | T^{(w)} + | \mathbf{1} \rangle` are exact. Futures are built by prepending symbols + to the column vectors :math:`T^{(w)}\mathbf{1}`. Words whose columns are + proportional induce the same posterior over :math:`S`, so they are merged. + The merge is lossless for the bottleneck, and long futures stay cheap for + processes whose reverse-time mixed states are finite. +* **Default future length**: the smallest :math:`L` with :math:`\mathbf{E} - + I[S; X_{0:L}] \le` ``tol`` (``1e-6``), capped at ``max_future_length`` + (a ``RuntimeWarning`` is issued if the cap is reached first). +* **Self-consistent iteration** (Blahut–Arimoto, Cover & Thomas + :cite:`Cover2006`, §10.8; Tishby et al. :cite:`tishby2000information`): + + .. math:: + + q(r \mid s) \propto q(r) \, e^{-\beta D_{KL}[P(X_{0:L} \mid s) \,\|\, q(X_{0:L} \mid r)]}, + \quad q(r) = \sum_s \pi(s) q(r \mid s), + \quad q(w \mid r) = \sum_s P(w \mid s) \, q(s \mid r), + + with :math:`|R| = |S|`. +* **Annealing and restarts**: the :math:`\beta` grid is swept upward. Each + :math:`\beta` is started from a perturbation of the previous optimum, from + the one-to-one encoder :math:`R = S`, and from ``restarts`` random encoders + (seeded by ``seed``). The candidate with the lowest Lagrangian wins. For + :math:`\beta \le 1` the optimum is always the trivial :math:`R` (zero rate); + as :math:`\beta \to \infty` the rate tends to :math:`C_\mu`, because the + length-:math:`L` morphs of distinct causal states differ once + :math:`I[S; X_{0:L}] = \mathbf{E}`. + +Example +======= + +.. code-block:: python + + from sofic.examples import golden_mean + from sofic.generators.predictive_rd import predictive_rate_distortion + + curve = predictive_rate_distortion(golden_mean(0.5), 40) + curve.rate # I[S; R] in bits, rises from 0 to C_mu + curve.relevant_information # I[R; X_{0:L}], rises from 0 to E + curve.distortion # I[S; X_{0:L}] - I[R; X_{0:L}] + +API +=== + +.. autoclass:: PredictiveRateDistortionCurve + :members: + +.. autofunction:: predictive_rate_distortion diff --git a/docs/generators/relative_entropy_rate.rst b/docs/generators/relative_entropy_rate.rst new file mode 100644 index 0000000..f2b9b98 --- /dev/null +++ b/docs/generators/relative_entropy_rate.rst @@ -0,0 +1,76 @@ +.. relative_entropy_rate.rst +.. py:module:: sofic.generators.relative_entropy_rate + +********************* +Relative Entropy Rate +********************* + +The relative entropy rate between stationary processes :math:`P` and :math:`Q`, +in bits per symbol, is + +.. math:: + + D(P \| Q) = \lim_{n \to \infty} \tfrac{1}{n} D(P_{0:n} \| Q_{0:n}) + +:cite:`Gray1990,Cover2006`. Both processes are taken in their stationary laws. + +Absolute continuity +=================== + +The rate is infinite as soon as some word with positive :math:`P` probability +has :math:`Q` probability zero. The stationary support of an HMM is a regular +language, so this is checked by building both support automata and asking +:meth:`~sofic.automata.base.LabeledAutomaton.includes`, which uses antichain +inclusion :cite:`DeWulf2006`. For example, the golden mean process is +absolutely continuous with respect to the fair coin, but not the other way +round. + +Exact rate for unifilar references +================================== + +When :math:`Q` is unifilar (an ε-machine or other unifilar HMM), + +.. math:: + + D(P \| Q) = -h_\mu(P) - \lim_{n \to \infty} \mathbb{E}_P[\log_2 Q(X_n \mid X_{0:n})]. + +For a :math:`k`-step Markov :math:`Q` this is Lemma 3.10 of Gray +:cite:`Gray1990` (§3.5). In general the past determines the set of :math:`Q` +states it is consistent with. The limit is then an average over the long-run +law of the finite chain of pairs (:math:`P` state, set of :math:`Q` states). +This requires :math:`P` to be ergodic and :math:`Q` to be exactly synchronized +by typical :math:`P` pasts. + +.. ipython:: + + In [1]: from sofic.examples import bernoulli, golden_mean + + In [2]: from sofic.generators.relative_entropy_rate import relative_entropy_rate + + @doctest float + In [3]: relative_entropy_rate(golden_mean(0.5), bernoulli(0.5)) + Out[3]: 0.3333333333333335 + + @doctest + In [4]: relative_entropy_rate(bernoulli(0.5), golden_mean(0.5)) + Out[4]: inf + +Bounds for general references +============================= + +:func:`relative_entropy_rate_bounds` enumerates the words of length :math:`n`. +It brackets the cross-entropy rate between superadditive and subadditive block +quantities built from :math:`\max_s Q_s` and :math:`\min_s Q_s`. It brackets the +entropy rate of :math:`P` between the hidden-Markov bounds +:math:`H[X_{n-1} \mid X_{0:n-1}, S_0] \le h_\mu \le H[X_{n-1} \mid X_{0:n-1}]` +(:cite:`Cover2006`, Theorem 4.5.1). + +API +=== + +.. autofunction:: relative_entropy_rate + +.. autofunction:: relative_entropy_rate_bounds + +.. autoclass:: RelativeEntropyRateBounds + :members: diff --git a/docs/references.bib b/docs/references.bib index 7a1b6e3..1456962 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -520,12 +520,12 @@ @article{Feldman1998b @article{riechers2017spectral2, author = {Riechers, Paul M. and Crutchfield, James P.}, title = {Spectral Simplicity of Apparent Complexity, Part {II}: Exact Complexities and Complexity Spectra}, - journal = {Physical Review E}, - volume = {96}, + journal = {Chaos: An Interdisciplinary Journal of Nonlinear Science}, + volume = {28}, number = {3}, - pages = {032304}, - year = {2017}, - doi = {10.1103/PhysRevE.96.032304}, + pages = {033116}, + year = {2018}, + doi = {10.1063/1.4986248}, } @article{bialek2001predictability, @@ -1143,3 +1143,160 @@ @article{Burman1994 year = {1994}, doi = {10.2307/2336965}, } + +@article{Thompson1968, + author = {Thompson, Ken}, + title = {Programming Techniques: Regular Expression Search Algorithm}, + journal = {Communications of the ACM}, + volume = {11}, + number = {6}, + pages = {419--422}, + year = {1968}, + doi = {10.1145/363347.363387}, +} + +@article{GarciaVidal1990, + author = {Garc{\'\i}a, Pedro and Vidal, Enrique}, + title = {Inference of k-Testable Languages in the Strict Sense and Application to Syntactic Pattern Recognition}, + journal = {IEEE Transactions on Pattern Analysis and Machine Intelligence}, + volume = {12}, + number = {9}, + pages = {920--925}, + year = {1990}, + doi = {10.1109/34.57687}, +} + +@book{Gray1990, + author = {Gray, Robert M.}, + title = {Entropy and Information Theory}, + publisher = {Springer-Verlag}, + address = {New York}, + year = {1990}, + doi = {10.1007/978-1-4757-3982-4}, +} + +@book{Cover2006, + author = {Cover, Thomas M. and Thomas, Joy A.}, + title = {Elements of Information Theory}, + edition = {2nd}, + publisher = {Wiley-Interscience}, + address = {Hoboken, NJ}, + year = {2006}, + doi = {10.1002/047174882X}, +} + +@inproceedings{DeWulf2006, + author = {De Wulf, Martin and Doyen, Laurent and Henzinger, Thomas A. and Raskin, Jean-Fran{\c{c}}ois}, + title = {Antichains: A New Algorithm for Checking Universality of Finite Automata}, + booktitle = {Computer Aided Verification (CAV 2006)}, + series = {Lecture Notes in Computer Science}, + volume = {4144}, + pages = {17--30}, + publisher = {Springer}, + year = {2006}, + doi = {10.1007/11817963_5}, +} + +@book{BaierKatoen2008, + author = {Baier, Christel and Katoen, Joost-Pieter}, + title = {Principles of Model Checking}, + publisher = {MIT Press}, + address = {Cambridge, MA}, + year = {2008}, + isbn = {978-0-262-02649-9}, +} + +@article{riechers2017spectral1, + author = {Riechers, Paul M. and Crutchfield, James P.}, + title = {Spectral Simplicity of Apparent Complexity, Part {I}: The Nondiagonalizable Metadynamics of Prediction}, + journal = {Chaos: An Interdisciplinary Journal of Nonlinear Science}, + volume = {28}, + number = {3}, + pages = {033115}, + year = {2018}, + doi = {10.1063/1.4985199}, +} + +@article{riechers2019fraudulent, + author = {Riechers, Paul M. and Crutchfield, James P.}, + title = {Fraudulent White Noise: Flat Power Spectra Belie Arbitrarily Complex Processes}, + journal = {Physical Review Research}, + volume = {3}, + number = {1}, + pages = {013170}, + year = {2021}, + doi = {10.1103/PhysRevResearch.3.013170}, +} + +@article{Marzen2016, + author = {Marzen, Sarah E. and Crutchfield, James P.}, + title = {Predictive Rate-Distortion for Infinite-Order {Markov} Processes}, + journal = {Journal of Statistical Physics}, + volume = {163}, + number = {6}, + pages = {1312--1338}, + year = {2016}, + doi = {10.1007/s10955-016-1520-1}, +} + +@article{tishby2000information, + author = {Tishby, Naftali and Pereira, Fernando C. and Bialek, William}, + title = {The Information Bottleneck Method}, + journal = {arXiv preprint physics/0004057}, + year = {2000}, +} + +@article{still2010optimal, + author = {Still, Susanne and Crutchfield, James P. and Ellison, Christopher J.}, + title = {Optimal Causal Inference: Estimating Stored Information and Approximating Causal Architecture}, + journal = {Chaos: An Interdisciplinary Journal of Nonlinear Science}, + volume = {20}, + number = {3}, + pages = {037111}, + year = {2010}, + doi = {10.1063/1.3489885}, +} + +@article{jurgens2021shannon, + author = {Jurgens, Alexandra M. and Crutchfield, James P.}, + title = {Shannon Entropy Rate of Hidden {Markov} Processes}, + journal = {Journal of Statistical Physics}, + volume = {183}, + number = {2}, + pages = {32}, + year = {2021}, + doi = {10.1007/s10955-021-02769-3}, + eprint = {2008.12886}, + archivePrefix = {arXiv}, +} + +@article{jurgens2021divergent, + author = {Jurgens, Alexandra M. and Crutchfield, James P.}, + title = {Divergent Predictive States: The Statistical Complexity Dimension of Stationary, Ergodic Hidden {Markov} Processes}, + journal = {Chaos}, + volume = {31}, + number = {8}, + pages = {083114}, + year = {2021}, + doi = {10.1063/5.0050460}, + eprint = {2102.10487}, + archivePrefix = {arXiv}, +} + +@phdthesis{Travers2013, + author = {Travers, Nicholas F.}, + title = {Bounds on Convergence of Entropy Rate Approximations in Hidden {Markov} Processes}, + school = {University of California, Davis}, + year = {2013}, + url = {https://csc.ucdavis.edu/~cmg/papers/Travers.UCDPhDDissertation2013.pdf}, +} + +@incollection{blackwell1957entropy, + author = {Blackwell, David}, + title = {The Entropy of Functions of Finite-State {Markov} Chains}, + booktitle = {Transactions of the First Prague Conference on Information Theory, Statistical Decision Functions, Random Processes}, + pages = {13--20}, + publisher = {Publishing House of the Czechoslovak Academy of Sciences}, + address = {Prague}, + year = {1957}, +} diff --git a/docs/shifts/shifts.rst b/docs/shifts/shifts.rst index 03a0425..8279236 100644 --- a/docs/shifts/shifts.rst +++ b/docs/shifts/shifts.rst @@ -17,6 +17,8 @@ type, Sofic shifts, Dyck shifts, topological Markov chains, and covers sofic_dyck_shift markov_dyck_shift topological_markov_chain + zeta_function + state_splitting sliding_block_code product_alphabet_shift textile diff --git a/docs/shifts/state_splitting.rst b/docs/shifts/state_splitting.rst new file mode 100644 index 0000000..1ec8892 --- /dev/null +++ b/docs/shifts/state_splitting.rst @@ -0,0 +1,73 @@ +.. state_splitting.rst +.. py:module:: sofic.shifts.state_splitting + +**************************************** +State Splitting and Conjugacy Invariants +**************************************** + +Out- and in-splitting a graph replaces a state by copies that share its +out-edges (respectively in-edges) according to a partition, giving an edge shift +conjugate to the original :cite:`LindMarcus1995` (§2.4). Amalgamation is the +inverse operation. By the Decomposition Theorem every conjugacy of edge shifts +is a composition of splitting and amalgamation codes (§7.1). + +Edges of a :class:`~sofic.shifts.tmc.TopologicalMarkovChain` carry +multiplicities, so a single edge is named ``(source, target, key, copy)``; +:func:`out_edges` and :func:`in_edges` list them, and a partition is a mapping +``{state: [part, ...]}`` of those edges. The division and edge matrices ``D`` +and ``E`` of a splitting satisfy ``A = DE`` and ``A' = ED`` (Theorem 2.4.12). + +.. ipython:: + + In [1]: import numpy as np; from sofic.shifts import TopologicalMarkovChain, out_edges, out_split, out_split_matrices + + In [2]: full = TopologicalMarkovChain.from_adjacency(np.array([[2]])) + + In [3]: (state,) = full.states(); partition = {state: [[edge] for edge in out_edges(full, state)]} + + In [4]: d, e = out_split_matrices(full, partition) + + @doctest + In [5]: (e @ d).tolist() + Out[5]: [[1, 1], [1, 1]] + +The complete out-splitting of the full 2-shift is its 2-block presentation +(Example 2.4.5). + +Conjugacy invariants +==================== + +The Bowen-Franks group ``BF(A) = Z^r / Z^r (I - A)`` is read from the Smith form +of ``I - A`` (Definition 7.4.15, Theorem 7.4.17); the sign of ``det(I - A)`` +completes the Parry-Sullivan / Bowen-Franks flow-equivalence invariant (§13.6). +The Jordan form away from zero ``J^x(A)`` keeps the Jordan blocks of nonzero +eigenvalues (Definition 7.4.9, Theorem 7.4.10). Both, together with the zeta +function (:doc:`zeta_function`), are unchanged by state splitting. + +.. ipython:: + + In [6]: from sofic.shifts import bowen_franks_group, jordan_form_away_from_zero + + @doctest + In [7]: bowen_franks_group(np.array([[4, 1], [1, 0]])) + Out[7]: BowenFranksGroup(invariant_factors=(4,), det_sign=-1) + + @doctest + In [8]: jordan_form_away_from_zero(np.array([[3, 1, 1], [2, 2, 1], [1, 2, 2]])) + Out[8]: ((1, 2), (5, 1)) + +API +=== + +.. autofunction:: out_edges +.. autofunction:: in_edges +.. autofunction:: out_split +.. autofunction:: in_split +.. autofunction:: out_split_matrices +.. autofunction:: in_split_matrices +.. autofunction:: out_amalgamate +.. autofunction:: in_amalgamate + +.. autoclass:: sofic.shifts.invariants.BowenFranksGroup +.. autofunction:: sofic.shifts.invariants.bowen_franks_group +.. autofunction:: sofic.shifts.invariants.jordan_form_away_from_zero diff --git a/docs/shifts/zeta_function.rst b/docs/shifts/zeta_function.rst new file mode 100644 index 0000000..75b57be --- /dev/null +++ b/docs/shifts/zeta_function.rst @@ -0,0 +1,51 @@ +.. zeta_function.rst + +************************************* +Periodic Points and the Zeta Function +************************************* + +The number ``p_n`` of points ``x`` with ``sigma^n x = x`` and the zeta function +``zeta(t) = exp(sum_n p_n t^n / n)`` are conjugacy invariants +:cite:`LindMarcus1995` (§6.4). All counts are exact Python integers; the zeta +function is a sympy rational function (install ``sofic[symbolic]``). + +* A :class:`~sofic.shifts.tmc.TopologicalMarkovChain` is an edge shift, so + ``p_n = tr A^n`` with edge multiplicities and ``zeta(t) = 1 / det(I - tA)`` + (Theorem 6.4.6). +* A :class:`~sofic.shifts.sft.ShiftOfFiniteType` given by forbidden words of + length at most ``L`` is recoded as the edge shift on allowed ``L``-blocks, + and the same trace formula applies. +* A :class:`~sofic.shifts.sofic.SoficShift` (or a presentation-only SFT) uses + Manning's formula on a right-resolving presentation: with ``A_j`` the + ``j``-th signed subset matrix, ``p_n = sum_j (-1)^(j+1) tr A_j^n`` and + ``zeta(t) = prod_j det(I - tA_j)^((-1)^j)`` (Theorem 6.4.8). Counting cycles + of the presentation itself overcounts points that have several + presentations, such as ``0^infinity`` in the even shift. + +.. ipython:: + + In [1]: from sofic.shifts import SoficShift + + In [2]: even = SoficShift(symbol_alphabet=frozenset("ab")) + + In [3]: even.add_transition(0, 0, "a"); even.add_transition(0, 1, "b"); even.add_transition(1, 0, "b"); + + @doctest + In [4]: [even.periodic_points(n) for n in range(1, 7)] + Out[4]: [2, 2, 5, 6, 12, 17] + + In [5]: even.zeta_function() + +API +=== + +.. autofunction:: sofic.shifts.algorithms.periodic_points +.. autofunction:: sofic.shifts.algorithms.zeta_function +.. autofunction:: sofic.shifts.algorithms.signed_subset_matrices +.. autofunction:: sofic.shifts.algorithms.reciprocal_characteristic_polynomial +.. automethod:: sofic.shifts.tmc.TopologicalMarkovChain.periodic_points +.. automethod:: sofic.shifts.tmc.TopologicalMarkovChain.zeta_function +.. automethod:: sofic.shifts.sft.ShiftOfFiniteType.periodic_points +.. automethod:: sofic.shifts.sft.ShiftOfFiniteType.zeta_function +.. automethod:: sofic.shifts.sofic.SoficShift.periodic_points +.. automethod:: sofic.shifts.sofic.SoficShift.zeta_function diff --git a/sofic/automata/__init__.py b/sofic/automata/__init__.py index 1b7ba55..cbaa2d1 100644 --- a/sofic/automata/__init__.py +++ b/sofic/automata/__init__.py @@ -64,6 +64,7 @@ from sofic.automata.learning.alergia import learn_pfa_alergia from sofic.automata.learning.dfasat import learn_dfa_sat from sofic.automata.learning.edsm import learn_dfa_edsm +from sofic.automata.learning.k_testable import learn_dfa_k_testable from sofic.automata.learning.nlstar import learn_prime_atomaton_nlstar, learn_rfsa_from_language, learn_rfsa_nlstar from sofic.automata.learning.observation import ObservationTable from sofic.automata.learning.papni import ( @@ -77,7 +78,7 @@ from sofic.automata.learning.rpni import learn_dfa_rpni from sofic.automata.nfa import NFA from sofic.automata.nwa import NestedWord, NestedWordAutomaton -from sofic.automata.regex import automaton_to_regex +from sofic.automata.regex import Regex, RegexSyntaxError, automaton_to_regex, parse_regex, regex_to_nfa from sofic.automata.subsequential import ( SubsequentialTransducer, WeightedFiniteStateTransducer, @@ -153,6 +154,8 @@ "NestedWord", "NestedWordAutomaton", "ObservationTable", + "Regex", + "RegexSyntaxError", "RegularLanguage", "ResidualFiniteStateAutomaton", "SingleEntryVisiblyPushdownAutomaton", @@ -192,6 +195,7 @@ "MISSING_TRANSITION", "learn_dfa_edsm", "learn_dfa_from_language", + "learn_dfa_k_testable", "learn_dfa_lstar", "learn_dfa_rpni", "learn_dfa_sat", @@ -208,6 +212,8 @@ "minimize_wheeler", "icdfa_next_flags", "next_icdfa_empty_string", + "parse_regex", + "regex_to_nfa", "encode_dyck_word", "encode_dyck_samples", "rank_idfa_string", diff --git a/sofic/automata/base.py b/sofic/automata/base.py index 6c3a226..ae76d07 100644 --- a/sofic/automata/base.py +++ b/sofic/automata/base.py @@ -103,6 +103,45 @@ def kleene_star(self) -> NFA: return kleene_star_nfa(self) + def is_empty(self) -> bool: + """Return whether no finite word is accepted :cite:`HopcroftUllman1979`.""" + from sofic.automata.decisions import is_empty + + return is_empty(self) + + def accepted_word(self) -> tuple[Any, ...] | None: + """Return a shortest accepted word, or ``None`` when the language is empty. + + Breadth-first search layered by the number of symbols read, so ε-moves + cost nothing :cite:`HopcroftUllman1979`. + """ + from sofic.automata.decisions import accepted_word + + return accepted_word(self) + + def is_universal(self, alphabet: frozenset[Any] | None = None) -> bool: + """Return whether every finite word over ``alphabet`` is accepted. + + ``alphabet`` defaults to :attr:`input_alphabet`, the alphabet + :meth:`complement` uses, so ``is_universal()`` agrees with + ``complement().is_empty()``. Decided by the forward antichain algorithm + on the subset construction, without determinizing :cite:`DeWulf2006`. + """ + from sofic.automata.decisions import is_universal + + return is_universal(self, self.input_alphabet if alphabet is None else alphabet) + + def includes(self, other: LabeledAutomaton) -> bool: + """Return whether ``other``'s language is contained in this one. + + Decided by the forward antichain inclusion algorithm, which explores + pairs of an ``other`` state and a subset of this automaton's states + without determinizing :cite:`DeWulf2006`. + """ + from sofic.automata.decisions import includes + + return includes(self, other) + def is_deterministic(self) -> bool: """Return whether this automaton is DFA-deterministic.""" from sofic.properties import is_deterministic_automaton diff --git a/sofic/automata/buchi.py b/sofic/automata/buchi.py index 7cae49b..6558e2d 100644 --- a/sofic/automata/buchi.py +++ b/sofic/automata/buchi.py @@ -21,3 +21,22 @@ def accepts_omega(self, word: Sequence[Any]) -> bool: from sofic.automata.buchi_simulation import accepts_omega_buchi return accepts_omega_buchi(self, word) + + def accepted_lasso(self) -> tuple[tuple[Any, ...], tuple[Any, ...]] | None: + """Return ``(prefix, loop)`` with ``prefix loop^omega`` accepted, or ``None`` when the ω-language is empty.""" + from sofic.automata.buchi_simulation import accepted_lasso_buchi + + return accepted_lasso_buchi(self) + + def is_empty(self) -> bool: + """Return whether no infinite word is accepted (Büchi emptiness :cite:`BaierKatoen2008`).""" + return self.accepted_lasso() is None + + def accepted_word(self) -> tuple[Any, ...] | None: + raise NotImplementedError("a Büchi automaton accepts infinite words; use accepted_lasso()") + + def is_universal(self, alphabet: frozenset[Any] | None = None) -> bool: + raise NotImplementedError("Büchi universality needs Büchi complementation, which sofic does not implement") + + def includes(self, other: Any) -> bool: + raise NotImplementedError("Büchi inclusion needs Büchi complementation, which sofic does not implement") diff --git a/sofic/automata/buchi_simulation.py b/sofic/automata/buchi_simulation.py index a0f9c96..28dc5df 100644 --- a/sofic/automata/buchi_simulation.py +++ b/sofic/automata/buchi_simulation.py @@ -8,7 +8,7 @@ import networkx as nx from sofic.automata.buchi import BuchiAutomaton -from sofic.graph import EPSILON +from sofic.graph import ATTR_SYMBOL, EPSILON def accepts_lasso_buchi(ba: BuchiAutomaton, prefix: Sequence[Any], loop: Sequence[Any]) -> bool: @@ -38,6 +38,52 @@ def accepts_omega_buchi(ba: BuchiAutomaton, word: Sequence[Any]) -> bool: raise NotImplementedError("BuchiAutomaton.accepts_omega supports ultimately periodic inputs as (prefix, loop) only") +def accepted_lasso_buchi(ba: BuchiAutomaton) -> tuple[tuple[Any, ...], tuple[Any, ...]] | None: + """Return ``(prefix, loop)`` with ``prefix loop^omega`` accepted, or ``None`` when the ω-language is empty. + + The ω-language is non-empty iff some accepting state reachable from an + initial state lies in a strongly connected component with a + symbol-consuming edge :cite:`BaierKatoen2008`. + """ + graph = nx.MultiDiGraph() + graph.add_nodes_from(ba.states()) + for transition in ba.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is not None: + graph.add_edge(transition.source, transition.target, symbol=symbol) + + def word(path: list[Hashable]) -> list[Any]: + symbols: list[Any] = [] + for source, target in zip(path, path[1:], strict=False): + labels = [data["symbol"] for data in graph.get_edge_data(source, target).values()] + consuming = sorted((label for label in labels if label is not EPSILON), key=repr) + if consuming and EPSILON not in labels: + symbols.append(consuming[0]) + return symbols + + reachable: set[Hashable] = set() + for initial in ba.initial_states: + reachable |= {initial} | nx.descendants(graph, initial) + for component in sorted(nx.strongly_connected_components(graph.subgraph(reachable)), key=repr): + accepting = sorted((state for state in component if state in ba.accepting_states), key=repr) + if not accepting: + continue + sub = graph.subgraph(component) + consuming = [(u, v, d["symbol"]) for u, v, d in sub.edges(data=True) if d["symbol"] is not EPSILON] + if not consuming: + continue + target = accepting[0] + source = min( + ba.initial_states, + key=lambda s: (not nx.has_path(graph, s, target), repr(s)), + ) + prefix = word(nx.shortest_path(graph, source, target)) + u, v, symbol = sorted(consuming, key=repr)[0] + loop = [*word(nx.shortest_path(sub, target, u)), symbol, *word(nx.shortest_path(sub, v, target))] + return tuple(prefix), tuple(loop) + return None + + def _accepting_cycle_reachable(ba: BuchiAutomaton, start: set[Hashable], loop: Sequence[Any]) -> bool: """Search the product of ``ba`` with loop positions for a reachable accepting cycle. diff --git a/sofic/automata/decisions.py b/sofic/automata/decisions.py new file mode 100644 index 0000000..c163cac --- /dev/null +++ b/sofic/automata/decisions.py @@ -0,0 +1,180 @@ +"""Emptiness, universality, and inclusion for finite-word automata. + +Universality and inclusion use the forward antichain algorithms of De Wulf, +Doyen, Henzinger, and Raskin :cite:`DeWulf2006`: the subset construction is +explored on the fly, and a subset is pruned when a subset of it has already been +reached, because every word rejected from the larger subset is also rejected +from the smaller one. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Hashable, Iterable +from typing import TYPE_CHECKING, Any + +from sofic.graph import ATTR_SYMBOL, EPSILON + +if TYPE_CHECKING: + from sofic.automata.base import LabeledAutomaton + +Successors = dict[Hashable, dict[Any, frozenset[Hashable]]] + + +def _successors(aut: LabeledAutomaton) -> Successors: + table: dict[Hashable, dict[Any, set[Hashable]]] = {state: {} for state in aut.states()} + for transition in aut.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is None: + continue + table.setdefault(transition.source, {}).setdefault(symbol, set()).add(transition.target) + return {state: {symbol: frozenset(targets) for symbol, targets in row.items()} for state, row in table.items()} + + +def _closure(table: Successors, states: Iterable[Hashable]) -> frozenset[Hashable]: + closure = set(states) + stack = list(closure) + while stack: + for target in table.get(stack.pop(), {}).get(EPSILON, ()): + if target not in closure: + closure.add(target) + stack.append(target) + return frozenset(closure) + + +def _post(table: Successors, states: Iterable[Hashable], symbol: Any) -> frozenset[Hashable]: + targets: set[Hashable] = set() + for state in states: + targets.update(table.get(state, {}).get(symbol, ())) + return _closure(table, targets) + + +def _symbols(table: Successors, state: Hashable) -> list[Any]: + return sorted((symbol for symbol in table.get(state, {}) if symbol is not EPSILON), key=repr) + + +def _dominated(antichain: list[frozenset[Hashable]], candidate: frozenset[Hashable]) -> bool: + return any(member <= candidate for member in antichain) + + +def _insert(antichain: list[frozenset[Hashable]], candidate: frozenset[Hashable]) -> None: + antichain[:] = [member for member in antichain if not candidate <= member] + antichain.append(candidate) + + +def accepted_word(aut: LabeledAutomaton) -> tuple[Any, ...] | None: + """Return a shortest accepted word, or ``None`` when the language is empty. + + Breadth-first search over states, layered by the number of non-ε symbols + read :cite:`HopcroftUllman1979`; ties are broken by ``repr`` order of the + symbols. + """ + table = _successors(aut) + parent: dict[Hashable, tuple[Hashable, Any] | None] = dict.fromkeys(aut.initial_states) + + def close(layer: list[Hashable]) -> list[Hashable]: + stack = list(layer) + while stack: + state = stack.pop() + for target in table.get(state, {}).get(EPSILON, ()): + if target not in parent: + parent[target] = (state, EPSILON) + layer.append(target) + stack.append(target) + return layer + + layer = close(list(parent)) + while layer: + for state in layer: + if state in aut.accepting_states: + return _trace(parent, state) + following: list[Hashable] = [] + for state in layer: + for symbol in _symbols(table, state): + for target in table[state][symbol]: + if target not in parent: + parent[target] = (state, symbol) + following.append(target) + layer = close(following) + return None + + +def _trace(parent: dict[Hashable, tuple[Hashable, Any] | None], state: Hashable) -> tuple[Any, ...]: + word: list[Any] = [] + step = parent[state] + while step is not None: + source, symbol = step + if symbol is not EPSILON: + word.append(symbol) + step = parent[source] + return tuple(reversed(word)) + + +def is_empty(aut: LabeledAutomaton) -> bool: + """Return whether ``aut`` accepts no finite word :cite:`HopcroftUllman1979`.""" + return accepted_word(aut) is None + + +def is_universal(aut: LabeledAutomaton, alphabet: frozenset[Any]) -> bool: + """Return whether ``aut`` accepts every finite word over ``alphabet``. + + Forward antichain search on the subset construction :cite:`DeWulf2006`: + a reachable subset without an accepting state witnesses a rejected word, + and only ``⊆``-minimal subsets need to be explored. + """ + table = _successors(aut) + accepting = aut.accepting_states + symbols = sorted((symbol for symbol in alphabet if symbol is not EPSILON), key=repr) + start = _closure(table, aut.initial_states) + if not start & accepting: + return False + antichain = [start] + queue = deque([start]) + while queue: + current = queue.popleft() + for symbol in symbols: + successor = _post(table, current, symbol) + if not successor & accepting: + return False + if _dominated(antichain, successor): + continue + _insert(antichain, successor) + queue.append(successor) + return True + + +def includes(larger: LabeledAutomaton, smaller: LabeledAutomaton) -> bool: + """Return whether ``L(smaller) ⊆ L(larger)``. + + Forward antichain search over pairs ``(q, S)`` of a state ``q`` of + ``smaller`` and the subset ``S`` of ``larger`` reached on the same word + :cite:`DeWulf2006`. A pair with ``q`` accepting and ``S`` non-accepting + witnesses a word in the difference; ``(q, S)`` is pruned when some + ``(q, S')`` with ``S' ⊆ S`` has already been reached. + """ + big = _successors(larger) + small = _successors(smaller) + antichains: dict[Hashable, list[frozenset[Hashable]]] = {} + queue: deque[tuple[Hashable, frozenset[Hashable]]] = deque() + + def visit(state: Hashable, subset: frozenset[Hashable]) -> bool: + if state in smaller.accepting_states and not subset & larger.accepting_states: + return False + antichain = antichains.setdefault(state, []) + if not _dominated(antichain, subset): + _insert(antichain, subset) + queue.append((state, subset)) + return True + + start = _closure(big, larger.initial_states) + for state in _closure(small, smaller.initial_states): + if not visit(state, start): + return False + while queue: + state, subset = queue.popleft() + for symbol in _symbols(small, state): + successor = _post(big, subset, symbol) + for target in _post(small, (state,), symbol): + if not visit(target, successor): + return False + return True diff --git a/sofic/automata/learning/__init__.py b/sofic/automata/learning/__init__.py index 88d9c5d..52f4cc7 100644 --- a/sofic/automata/learning/__init__.py +++ b/sofic/automata/learning/__init__.py @@ -22,6 +22,7 @@ from sofic.automata.learning.alergia import learn_pfa_alergia from sofic.automata.learning.dfasat import learn_dfa_sat from sofic.automata.learning.edsm import learn_dfa_edsm +from sofic.automata.learning.k_testable import learn_dfa_k_testable from sofic.automata.learning.nlstar import learn_prime_atomaton_nlstar, learn_rfsa_from_language, learn_rfsa_nlstar from sofic.automata.learning.observation import ObservationTable from sofic.automata.learning.papni import ( @@ -52,6 +53,7 @@ "is_well_matched", "learn_dfa_edsm", "learn_dfa_from_language", + "learn_dfa_k_testable", "learn_dfa_lstar", "learn_dfa_rpni", "learn_dfa_sat", diff --git a/sofic/automata/learning/k_testable.py b/sofic/automata/learning/k_testable.py new file mode 100644 index 0000000..6c13d9a --- /dev/null +++ b/sofic/automata/learning/k_testable.py @@ -0,0 +1,80 @@ +"""Passive DFA learning of k-testable languages in the strict sense.""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Iterable, Sequence +from typing import Any + +from sofic.automata.dfa import DFA + +__all__ = ["learn_dfa_k_testable"] + + +def learn_dfa_k_testable( + samples: Sequence[Sequence[Any]], + k: int, + *, + alphabet: Iterable[Any] | None = None, +) -> DFA: + """Learn the smallest strictly k-testable language containing ``samples``. + + Implements the k-TSSI algorithm of García & Vidal :cite:`GarciaVidal1990`. + From the positive samples it collects the length-``k - 1`` prefixes ``I`` + and suffixes ``F``, the length-``k`` factors ``T``, and the sample words + shorter than ``k - 1``. A word ``w`` is accepted iff it is one of those + short sample words, or ``|w| >= k - 1``, its prefix is in ``I``, its suffix + is in ``F``, and every length-``k`` factor is in ``T``. The language grows + monotonically with the sample and identifies every strictly k-testable + language in the limit from positive data; it shrinks (weakly) as ``k`` grows. + + States are the sample prefixes shorter than ``k - 1`` together with the + reachable length-``k - 1`` windows, labelled by the tuples themselves; the + DFA is partial and trim on the accessible side. ``alphabet`` defaults to + the symbols occurring in ``samples`` and must contain all of them. + """ + if k < 1: + raise ValueError("k must be at least 1") + words = [tuple(word) for word in samples] + window = k - 1 + occurring = {symbol for word in words for symbol in word} + if alphabet is None: + symbols = tuple(sorted(occurring, key=repr)) + else: + symbols = tuple(sorted(set(alphabet), key=repr)) + if missing := occurring - set(symbols): + raise ValueError(f"sample symbols {sorted(missing, key=repr)!r} are not in the alphabet") + + heads = {word[:length] for word in words for length in range(min(len(word), window) + 1)} + suffixes = {word[len(word) - window :] for word in words if len(word) >= window} + factors = {word[index : index + k] for word in words for index in range(len(word) - k + 1)} + short = {word for word in words if len(word) < window} + + dfa = DFA(input_alphabet=frozenset(symbols), initial_states=frozenset({()})) + dfa.graph.add_state(()) + accepting: set[tuple[Any, ...]] = set() + queue: deque[tuple[Any, ...]] = deque([()]) + seen = {()} + while queue: + state = queue.popleft() + if state in short or (len(state) == window and state in suffixes): + accepting.add(state) + for symbol in symbols: + extended = (*state, symbol) + if len(state) < window: + if extended not in heads: + continue + target = extended + else: + if extended not in factors: + continue + target = extended[1:] + if target not in seen: + seen.add(target) + dfa.graph.add_state(target) + queue.append(target) + dfa.add_transition(state, target, symbol) + + dfa.accepting_states = frozenset(accepting) + dfa.validate() + return dfa diff --git a/sofic/automata/nfa.py b/sofic/automata/nfa.py index 020e102..77e8bc5 100644 --- a/sofic/automata/nfa.py +++ b/sofic/automata/nfa.py @@ -2,7 +2,7 @@ from __future__ import annotations -from collections.abc import Hashable, Sequence +from collections.abc import Hashable, Iterable, Sequence from typing import TYPE_CHECKING, Any from sofic.automata.base import LabeledAutomaton @@ -37,3 +37,10 @@ def minimize( from sofic.automata.algorithms import minimize return minimize(self, algorithm=algorithm, alphabet=alphabet) + + @classmethod + def from_regex(cls, text: str, *, alphabet: Iterable[Any] | None = None) -> NFA: + """Return a Thompson ε-NFA for ``text``; see :func:`~sofic.automata.regex.regex_to_nfa`.""" + from sofic.automata.regex import regex_to_nfa + + return regex_to_nfa(text, alphabet=alphabet) diff --git a/sofic/automata/regex.py b/sofic/automata/regex.py index 34311bc..c0575b9 100644 --- a/sofic/automata/regex.py +++ b/sofic/automata/regex.py @@ -1,52 +1,91 @@ -"""Automata to regular-expression conversion.""" +"""Conversions between automata and regular expressions. + +Regular expressions use a small dialect that coincides with Python's :mod:`re` +syntax whenever every symbol is a single-character string: + +* a single character is a literal symbol; ``\\c`` escapes any character ``c``; +* ``'...'`` quotes a multi-character symbol, with ``\\'`` and ``\\\\`` escapes + inside the quotes (``'ab'`` is one symbol, ``ab`` is two); +* ``ε`` or ``(?:)`` denotes the empty word and ``∅`` or ``(?!)`` the empty + language; +* ``|`` is alternation, juxtaposition is concatenation, and the postfix + operators ``*``, ``+``, ``?`` bind tightest; +* ``(...)`` and ``(?:...)`` group. + +The characters ``. [ ] { } ^ $`` are reserved (their :mod:`re` meaning is not +supported) and must be escaped to be used as symbols. +""" from __future__ import annotations import re +from collections.abc import Callable, Iterable from dataclasses import dataclass from typing import Any from sofic.automata.base import LabeledAutomaton +from sofic.automata.nfa import NFA +from sofic.exceptions import RegexSyntaxError from sofic.graph import ATTR_SYMBOL, EPSILON +__all__ = ["Regex", "RegexSyntaxError", "automaton_to_regex", "parse_regex", "regex_to_nfa"] + @dataclass(frozen=True, slots=True) -class _Regex: +class Regex: + """Regular-expression syntax tree. + + ``kind`` is one of ``"empty"``, ``"epsilon"``, ``"literal"`` (``parts`` is + the symbol), ``"union"``, ``"concat"``, ``"star"``, ``"plus"``, or + ``"optional"`` (``parts`` holds the operand nodes). + """ + kind: str parts: tuple[Any, ...] = () def pattern(self) -> str: + """Render this tree in the module's regular-expression dialect.""" if self.kind == "empty": return "(?!)" if self.kind == "epsilon": return "(?:)" if self.kind == "literal": - return re.escape(str(self.parts[0])) + return _render_symbol(self.parts[0]) if self.kind == "union": return "(?:" + "|".join(part.pattern() for part in self.parts) + ")" if self.kind == "concat": return "".join(_atom(part) for part in self.parts) - if self.kind == "star": - return _atom(self.parts[0]) + "*" + if self.kind in _POSTFIX_KINDS: + return _atom(self.parts[0]) + _POSTFIX_KINDS[self.kind] raise ValueError(f"unknown regex node {self.kind!r}") -_EMPTY = _Regex("empty") -_EPSILON = _Regex("epsilon") +_EMPTY = Regex("empty") +_EPSILON = Regex("epsilon") +_POSTFIX_KINDS = {"star": "*", "plus": "+", "optional": "?"} +_POSTFIX_OPERATORS = {operator: kind for kind, operator in _POSTFIX_KINDS.items()} +_QUOTE = "'" +_RESERVED_SYMBOLS = frozenset({"ε", "∅", _QUOTE}) +_UNSUPPORTED_METACHARACTERS = frozenset(".[]{}^$") def automaton_to_regex(automaton: LabeledAutomaton) -> str: - """Return a Python-regex-compatible expression for ``automaton``. + """Return a regular expression for ``automaton`` by state elimination. - The result uses ``(?!)`` for the empty language and ``(?:)`` for epsilon. - Symbols are rendered with ``str(symbol)`` and escaped for Python's - :mod:`re` engine. + The construction follows Kleene's theorem :cite:`Kleene1956` in the + state-elimination form of :cite:`HopcroftUllman1979`. The result uses + ``(?!)`` for the empty language and ``(?:)`` for epsilon. Symbols are + rendered with ``str(symbol)``: single characters are escaped as for + Python's :mod:`re` engine, so the output is a valid :mod:`re` pattern when + every symbol is a single-character string, and longer symbols are quoted + (``'ab'``) so the output stays unambiguous. :func:`parse_regex` and + :func:`regex_to_nfa` read the output back. """ start = object() final = object() states = sorted(automaton.states(), key=repr) remaining: list[Any] = [start, *states, final] - labels: dict[tuple[Any, Any], _Regex] = {} + labels: dict[tuple[Any, Any], Regex] = {} for state in automaton.initial_states: labels[(start, state)] = _union(labels.get((start, state), _EMPTY), _EPSILON) @@ -56,7 +95,7 @@ def automaton_to_regex(automaton: LabeledAutomaton) -> str: symbol = transition.data.get(ATTR_SYMBOL) if symbol is None: continue - label = _EPSILON if symbol is EPSILON else _Regex("literal", (symbol,)) + label = _EPSILON if symbol is EPSILON else Regex("literal", (symbol,)) key = (transition.source, transition.target) labels[key] = _union(labels.get(key, _EMPTY), label) @@ -79,8 +118,206 @@ def automaton_to_regex(automaton: LabeledAutomaton) -> str: return labels.get((start, final), _EMPTY).pattern() -def _union(*terms: _Regex) -> _Regex: - parts: list[_Regex] = [] +def parse_regex(text: str) -> Regex: + """Parse ``text`` in the module's regular-expression dialect. + + The tree is returned unsimplified. Raises :class:`RegexSyntaxError` on + unbalanced parentheses, operators without operands, empty alternatives, + unterminated quotes or escapes, and reserved metacharacters. + """ + return _Parser(text).parse() + + +def regex_to_nfa(text: str, *, alphabet: Iterable[Any] | None = None) -> NFA: + """Return an ε-NFA for ``text`` via Thompson's construction. + + Implements :cite:`Thompson1968`: every literal and ε becomes a two-state + fragment, and union, concatenation, and the postfix operators wire + fragments together with ε-transitions, so the NFA has one initial and one + accepting state and at most two states per syntax-tree node. + + Without ``alphabet``, symbols are the parsed strings and the input alphabet + is the set of symbols occurring in ``text``. With ``alphabet``, each parsed + string is resolved to the alphabet element equal to it, or else to the + unique element whose ``str`` matches it, so non-string symbols round-trip + through :func:`automaton_to_regex`; the NFA's input alphabet is + ``alphabet``. + """ + tree = parse_regex(text) + symbols = None if alphabet is None else frozenset(alphabet) + resolve = _symbol_resolver(symbols) + transitions: list[tuple[int, int, Any]] = [] + count = 0 + + def new_state() -> int: + nonlocal count + count += 1 + return count - 1 + + def build(node: Regex) -> tuple[int, int]: + if node.kind == "concat": + fragments = [build(part) for part in node.parts] + for (_, accept), (start, _) in zip(fragments, fragments[1:], strict=False): + transitions.append((accept, start, EPSILON)) + return fragments[0][0], fragments[-1][1] + start, accept = new_state(), new_state() + if node.kind == "epsilon": + transitions.append((start, accept, EPSILON)) + elif node.kind == "literal": + transitions.append((start, accept, resolve(node.parts[0]))) + elif node.kind == "union": + for part in node.parts: + inner_start, inner_accept = build(part) + transitions.append((start, inner_start, EPSILON)) + transitions.append((inner_accept, accept, EPSILON)) + elif node.kind in _POSTFIX_KINDS: + inner_start, inner_accept = build(node.parts[0]) + transitions.append((start, inner_start, EPSILON)) + transitions.append((inner_accept, accept, EPSILON)) + if node.kind != "optional": + transitions.append((inner_accept, inner_start, EPSILON)) + if node.kind != "plus": + transitions.append((start, accept, EPSILON)) + return start, accept + + start, accept = build(tree) + if symbols is None: + symbols = frozenset(symbol for _, _, symbol in transitions if symbol is not EPSILON) + nfa = NFA(input_alphabet=symbols, initial_states=frozenset({start}), accepting_states=frozenset({accept})) + for state in range(count): + nfa.graph.add_state(state) + for source, target, symbol in transitions: + nfa.add_transition(source, target, symbol) + return nfa + + +class _Parser: + def __init__(self, text: str) -> None: + self.text = text + self.position = 0 + + def parse(self) -> Regex: + tree = self._union() + if self.position < len(self.text): + self._fail("unbalanced parenthesis ')'") + return tree + + def _fail(self, message: str, position: int | None = None) -> None: + raise RegexSyntaxError(message, self.text, self.position if position is None else position) + + def _peek(self) -> str | None: + return self.text[self.position] if self.position < len(self.text) else None + + def _union(self) -> Regex: + alternatives = [self._concat()] + while self._peek() == "|": + self.position += 1 + alternatives.append(self._concat()) + return alternatives[0] if len(alternatives) == 1 else Regex("union", tuple(alternatives)) + + def _concat(self) -> Regex: + factors: list[Regex] = [] + while self._peek() not in {None, "|", ")"}: + factors.append(self._repeat()) + if not factors: + self._fail("expected an expression") + return factors[0] if len(factors) == 1 else Regex("concat", tuple(factors)) + + def _repeat(self) -> Regex: + if self._peek() in _POSTFIX_OPERATORS: + self._fail(f"nothing to repeat before {self._peek()!r}") + node = self._atom() + if (operator := self._peek()) in _POSTFIX_OPERATORS: + self.position += 1 + node = Regex(_POSTFIX_OPERATORS[operator], (node,)) + if self._peek() in _POSTFIX_OPERATORS: + self._fail("multiple repeat") + return node + + def _atom(self) -> Regex: + text = self.text + position = self.position + char = text[position] + if text.startswith("(?:)", position): + self.position += 4 + return _EPSILON + if text.startswith("(?!)", position): + self.position += 4 + return _EMPTY + if char == "(": + if text.startswith("(?:", position): + self.position += 3 + elif text.startswith("(?", position): + self._fail("unsupported group extension '(?'") + else: + self.position += 1 + node = self._union() + if self._peek() != ")": + self._fail("missing ')' for parenthesis", position) + self.position += 1 + return node + self.position += 1 + if char == "ε": + return _EPSILON + if char == "∅": + return _EMPTY + if char == "\\": + if self.position >= len(text): + self._fail("dangling backslash", position) + self.position += 1 + return Regex("literal", (text[self.position - 1],)) + if char == _QUOTE: + return Regex("literal", (self._quoted(position),)) + if char in _UNSUPPORTED_METACHARACTERS: + self._fail(f"unsupported metacharacter {char!r} (escape it as '\\{char}')", position) + return Regex("literal", (char,)) + + def _quoted(self, opening: int) -> str: + chars: list[str] = [] + text = self.text + while self.position < len(text): + char = text[self.position] + self.position += 1 + if char == _QUOTE: + return "".join(chars) + if char == "\\": + if self.position >= len(text): + break + char = text[self.position] + self.position += 1 + chars.append(char) + self._fail("unterminated quoted symbol", opening) + raise AssertionError # pragma: no cover + + +def _symbol_resolver(symbols: frozenset[Any] | None) -> Callable[[str], Any]: + if symbols is None: + return lambda token: token + by_text: dict[str, list[Any]] = {} + for symbol in symbols: + by_text.setdefault(str(symbol), []).append(symbol) + + def resolve(token: str) -> Any: + if token in symbols: + return token + matches = by_text.get(token, []) + if len(matches) != 1: + reason = "not in" if not matches else "ambiguous in" + raise ValueError(f"regex symbol {token!r} is {reason} the alphabet") + return matches[0] + + return resolve + + +def _render_symbol(symbol: Any) -> str: + text = str(symbol) + if len(text) == 1: + return "\\" + text if text in _RESERVED_SYMBOLS else re.escape(text) + return _QUOTE + text.replace("\\", "\\\\").replace(_QUOTE, "\\" + _QUOTE) + _QUOTE + + +def _union(*terms: Regex) -> Regex: + parts: list[Regex] = [] for term in terms: if term == _EMPTY: continue @@ -93,11 +330,11 @@ def _union(*terms: _Regex) -> _Regex: return _EMPTY if len(unique) == 1: return unique[0] - return _Regex("union", tuple(unique)) + return Regex("union", tuple(unique)) -def _concat(*terms: _Regex) -> _Regex: - parts: list[_Regex] = [] +def _concat(*terms: Regex) -> Regex: + parts: list[Regex] = [] for term in terms: if term == _EMPTY: return _EMPTY @@ -111,18 +348,18 @@ def _concat(*terms: _Regex) -> _Regex: return _EPSILON if len(parts) == 1: return parts[0] - return _Regex("concat", tuple(parts)) + return Regex("concat", tuple(parts)) -def _star(term: _Regex) -> _Regex: +def _star(term: Regex) -> Regex: if term in {_EMPTY, _EPSILON}: return _EPSILON if term.kind == "star": return term - return _Regex("star", (term,)) + return Regex("star", (term,)) -def _atom(term: _Regex) -> str: - if term.kind in {"literal", "epsilon", "empty", "star"}: +def _atom(term: Regex) -> str: + if term.kind in {"literal", "epsilon", "empty", *_POSTFIX_KINDS}: return term.pattern() return "(?:" + term.pattern() + ")" diff --git a/sofic/exceptions.py b/sofic/exceptions.py index e9d39ff..49259c1 100644 --- a/sofic/exceptions.py +++ b/sofic/exceptions.py @@ -47,3 +47,12 @@ class InfiniteTransductionError(SoficError): class NonWellMatchedLanguageError(SoficValidationError): """Raised when a construction defined for well-matched languages meets a pending call or return.""" + + +class RegexSyntaxError(SoficError, ValueError): + """Raised when a regular expression cannot be parsed.""" + + def __init__(self, message: str, text: str, position: int) -> None: + super().__init__(f"{message} at position {position} in {text!r}") + self.text = text + self.position = position diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py index 92ce391..db8de0c 100644 --- a/sofic/generators/__init__.py +++ b/sofic/generators/__init__.py @@ -8,6 +8,8 @@ channel_statistical_complexity, driven_entropy_rate, ) +from sofic.generators.complexity_dimension import StatisticalComplexityDimension, statistical_complexity_dimension +from sofic.generators.correlations import autocorrelation, mutual_information_function, power_spectrum from sofic.generators.directional_flow import ( directed_information, independent_pair_generator, @@ -21,6 +23,7 @@ from sofic.generators.lumping import LumpabilityError, 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 from sofic.generators.minimal_generative_model import ( FunctionalGenerativeModel, GacsKornerGenerativeModel, @@ -35,7 +38,13 @@ from sofic.generators.moore import MooreHMM from sofic.generators.nmachine import NMachine from sofic.generators.pfa import ProbabilisticFiniteAutomaton +from sofic.generators.predictive_rd import PredictiveRateDistortionCurve, predictive_rate_distortion from sofic.generators.quasi_realization import QuasiRealization +from sofic.generators.relative_entropy_rate import ( + RelativeEntropyRateBounds, + relative_entropy_rate, + relative_entropy_rate_bounds, +) from sofic.generators.stack_hmm import HiddenMarkovStackModel from sofic.generators.topological_epsilon_enumeration import ( count_topological_epsilon_machines, @@ -114,4 +123,17 @@ "wheeler_statistical_complexity", "word_cylinder_measure", "wyner_generative_model", + "PredictiveRateDistortionCurve", + "RelativeEntropyRateBounds", + "autocorrelation", + "mutual_information_function", + "power_spectrum", + "predictive_rate_distortion", + "relative_entropy_rate", + "relative_entropy_rate_bounds", + "EntropyRateEstimate", + "StatisticalComplexityDimension", + "entropy_rate_blackwell", + "entropy_rate_bounds", + "statistical_complexity_dimension", ] diff --git a/sofic/generators/base.py b/sofic/generators/base.py index bc28e17..b0d7fcc 100644 --- a/sofic/generators/base.py +++ b/sofic/generators/base.py @@ -14,6 +14,7 @@ from sofic.automata.dfa import DFA from sofic.automata.nfa import NFA from sofic.generators.mealy import MealyHMM + from sofic.generators.measures import EntropyRateMethod from sofic.generators.prob import SymbolConstraints from sofic.shifts.sofic import SoficShift @@ -187,10 +188,16 @@ def to_mealy(self) -> MealyHMM: """Return an equivalent Mealy-style presentation.""" raise NotImplementedError(f"{type(self).__name__} must implement to_mealy()") - def entropy_rate(self) -> float: - from sofic.generators.measures import entropy_rate_hmm + def entropy_rate(self, method: EntropyRateMethod = "auto", **kwargs: Any) -> Any: + """Shannon entropy rate in bits per symbol; see :func:`sofic.generators.measures.entropy_rate`. - return entropy_rate_hmm(self) + The default ``method="auto"`` is exact for unifilar presentations and for + non-unifilar ones whose mixed states close finitely, and otherwise uses + converged upper/lower bounds :cite:`Cover2006`. + """ + from sofic.generators.measures import entropy_rate + + return entropy_rate(self, method, **kwargs) def joint_block_distribution(self, block_length: int = 2) -> Any: from sofic.generators.measures import joint_block_distribution @@ -263,6 +270,15 @@ def is_equal_process( return is_equal_process(self, other, start1=start1, start2=start2, rtol=rtol, atol=atol) + def relative_entropy_rate(self, other: HiddenMarkovModel) -> float: + """Return ``D(self || other)`` in bits per symbol :cite:`Gray1990`. + + See :func:`~sofic.generators.relative_entropy_rate.relative_entropy_rate`. + """ + from sofic.generators.relative_entropy_rate import relative_entropy_rate + + return relative_entropy_rate(self, other) + 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 diff --git a/sofic/generators/complexity_dimension.py b/sofic/generators/complexity_dimension.py new file mode 100644 index 0000000..0957be4 --- /dev/null +++ b/sofic/generators/complexity_dimension.py @@ -0,0 +1,164 @@ +"""Statistical complexity dimension of hidden Markov processes. + +The mixed states of an :math:`N`-state HMM evolve on the +:math:`(N-1)`-simplex as a place-dependent iterated function system (IFS) +with maps :math:`f^{(x)}(\\eta) = \\eta T^{(x)} / \\eta T^{(x)} \\mathbf{1}` +chosen with probabilities :math:`\\eta T^{(x)} \\mathbf{1}`. When the mixed +states do not close finitely the statistical complexity diverges, and its +divergence rate is the information dimension :math:`d_\\mu` of the Blackwell +measure. Jurgens & Crutchfield bound (and, under the open set condition, +equate) it with the Lyapunov dimension of the IFS, in which the entropy rate +plays the role of the leading Lyapunov exponent :cite:`jurgens2021divergent`. +""" + +from __future__ import annotations + +from dataclasses import dataclass + +import numpy as np +from scipy.linalg import null_space + +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.base import StochasticModel +from sofic.generators.markov import MarkovChain +from sofic.generators.measures import ( + _numeric_symbol_tensors, + _row_entropies_bits, + batch_means_stderr, + entropy_rate_hmm, + mixed_state_walk, +) + + +@dataclass(frozen=True) +class StatisticalComplexityDimension: + """Statistical complexity dimension and the quantities it is computed from. + + ``entropy_rate`` and ``lyapunov_exponents`` are in bits (base-2 logarithms) + per symbol, exponents sorted in decreasing order. ``finite`` is true when + the presentation is unifilar or its mixed states close within the state + cap, in which case ``dimension`` is zero by definition and + ``entropy_rate`` is exact. + """ + + dimension: float + entropy_rate: float + entropy_rate_stderr: float + lyapunov_exponents: np.ndarray + finite: bool + + +def ifs_lyapunov_dimension(entropy_rate: float, exponents: np.ndarray) -> float: + r"""Lyapunov dimension of a mixed-state IFS (:cite:`jurgens2021divergent`, Eq. 13). + + With the (negative) exponents sorted as :math:`\lambda_1 \ge \lambda_2 \ge + \cdots` and :math:`k` the largest index with + :math:`h_\mu + \sum_{i \le k} \lambda_i > 0`, + + .. math:: d_\Gamma = k + \frac{h_\mu + \sum_{i \le k} \lambda_i}{|\lambda_{k+1}|}, + + saturating at the simplex dimension when :math:`h_\mu + \sum_i \lambda_i + \ge 0`. This is the Kaplan-Yorke formula with :math:`h_\mu` as the leading + exponent; for two-state HMMs it is :math:`-h_\mu / \lambda_1` (Eq. 15). + Both arguments must use the same logarithm base. + """ + ordered = np.sort(np.asarray(exponents, dtype=float))[::-1] + total = float(entropy_rate) + for k, value in enumerate(ordered): + if total + value <= 0.0: + return float(k + (total / abs(value) if total > 0.0 else 0.0)) + total += value + return float(len(ordered)) + + +def statistical_complexity_dimension( + model: StochasticModel, + *, + n_samples: int = 100_000, + burn_in: int = 1_000, + seed: int | np.random.Generator | None = None, + n_batches: int = 20, + max_states: int | None = 1_000, +) -> StatisticalComplexityDimension: + r"""Estimate the statistical complexity dimension :math:`d_\mu` :cite:`jurgens2021divergent`. + + Follows the algorithm of Jurgens & Crutchfield (Sec. V.B): sample a + mixed-state trajectory from the stationary state law (discarding + ``burn_in`` steps), estimate the entropy rate as the Blackwell time average + of :math:`H[X \mid \eta_t]` :cite:`jurgens2021shannon`, estimate the + :math:`N - 1` Lyapunov exponents of the IFS by QR (pull-back) iteration of + the Jacobians of :math:`f^{(x_t)}` restricted to the simplex tangent space + :math:`\{v : v \mathbf{1} = 0\}`, and combine them with + :func:`ifs_lyapunov_dimension`. Constant (synchronizing) maps contribute + exponents of :math:`-\infty`. + + The Lyapunov dimension equals the information dimension of the Blackwell + measure when the IFS satisfies the open set condition and is otherwise an + upper bound (Eq. 14). It describes the mixed states of *this* + presentation, which are the causal states unless distinct mixed states + predict identically (the paper's Appendix B). + + A process with finitely many causal states has :math:`d_\mu = 0` by + definition (Sec. V). Unifilar presentations are therefore reported as + finite, as are presentations whose mixed-state construction closes within + ``max_states`` beliefs (``None`` skips that check). + """ + pi, tensors = _numeric_symbol_tensors(model) + rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed) + beliefs, symbols, probabilities = mixed_state_walk(pi, tensors, n_samples=n_samples, burn_in=burn_in, rng=rng) + entropies = _row_entropies_bits(probabilities) + rate = float(entropies.mean()) + stderr = batch_means_stderr(entropies, n_batches) + exponents = _lyapunov_spectrum(tensors, beliefs, symbols) + + exact_rate = _finite_entropy_rate(model, max_states) + if exact_rate is not None: + return StatisticalComplexityDimension(0.0, exact_rate, 0.0, exponents, True) + return StatisticalComplexityDimension(ifs_lyapunov_dimension(rate, exponents), rate, stderr, exponents, False) + + +def _finite_entropy_rate(model: StochasticModel, max_states: int | None) -> float | None: + """Return the exact entropy rate when the process has finitely many causal states, else ``None``.""" + if isinstance(model, MarkovChain): + return float(model.entropy_rate()) + is_unifilar = getattr(model, "is_unifilar", None) + if is_unifilar is not None and is_unifilar(): + return float(entropy_rate_hmm(model)) + if max_states is None: + return None + from sofic.generators.mixed_state_construction import build_mixed_state_presentation + + try: + presentation = build_mixed_state_presentation(model.to_mealy(), max_states=max_states) + except MixedStateExplosionError: + return None + return float(entropy_rate_hmm(presentation.to_recurrent())) + + +def _lyapunov_spectrum( + tensors: np.ndarray, + beliefs: np.ndarray, + symbols: np.ndarray, +) -> np.ndarray: + """QR-iterated Lyapunov exponents (bits per step) of the projective maps along a mixed-state path. + + For row-vector beliefs the Jacobian of :math:`f^{(x)}` at :math:`\\eta` acts as + :math:`v \\mapsto v (T^{(x)} - T^{(x)}\\mathbf{1} \\, f^{(x)}(\\eta)) / \\eta T^{(x)} \\mathbf{1}` + and preserves the tangent space spanned by the columns of ``basis``. + """ + n_states = tensors.shape[1] + if n_states < 2: + return np.zeros(0) + basis = null_space(np.ones((1, n_states))) + frame = np.eye(n_states - 1) + log_growth = np.zeros(n_states - 1) + with np.errstate(divide="ignore"): + for eta, x in zip(beliefs, symbols, strict=True): + matrix = tensors[x] + mass = eta @ matrix + p = float(mass.sum()) + jacobian = (matrix - np.outer(matrix.sum(axis=1), mass / p)) / p + reduced = basis.T @ jacobian.T @ basis + frame, upper = np.linalg.qr(reduced @ frame) + log_growth += np.log2(np.abs(np.diag(upper))) + return np.sort(log_growth / len(symbols))[::-1] diff --git a/sofic/generators/correlations.py b/sofic/generators/correlations.py new file mode 100644 index 0000000..2b9b5d8 --- /dev/null +++ b/sofic/generators/correlations.py @@ -0,0 +1,240 @@ +"""Pairwise correlations of hidden Markov processes in closed form. + +Every quantity here is a function of the stationary state law :math:`\\pi` and +the symbol-labeled transition matrices :math:`T^{(x)}` (see +:mod:`sofic.generators.matrices`), with :math:`T = \\sum_x T^{(x)}`. For a real +observable :math:`f` of the emitted symbol let :math:`\\Omega = \\sum_x f(x) +T^{(x)}`. The autocorrelation is :math:`\\gamma(\\tau) = \\langle \\pi | \\Omega +T^{|\\tau| - 1} \\Omega | \\mathbf{1} \\rangle` for :math:`\\tau \\neq 0` and the +continuous part of the power spectrum is a resolvent of :math:`T` (Riechers & +Crutchfield :cite:`riechers2017spectral1`, Eqs. (3) and (5)). The lagged joint +symbol law :math:`P(X_0 = x, X_\\tau = y) = \\langle \\pi | T^{(x)} T^{\\tau - 1} +T^{(y)} | \\mathbf{1} \\rangle` gives the mutual information function. Pairwise +statistics can be flat (white) even for highly structured processes (Riechers & +Crutchfield :cite:`riechers2019fraudulent`). +""" + +from __future__ import annotations + +from collections.abc import Callable, Hashable, Mapping +from typing import Any + +import numpy as np +import scipy.linalg + +from sofic.generators.matrices import emission_tensors + +Observable = Mapping[Hashable, float] | Callable[[Hashable], float] + +_PERIPHERAL_TOL = 1e-8 + + +def _process_tensors(hmm: Any) -> tuple[np.ndarray, list[Hashable], np.ndarray]: + """Return the stationary law, the symbol order, and the stacked ``T^(x)`` of shape (k, n, n).""" + pi, joint = emission_tensors(hmm, policy="stationary") + symbols = list(joint) + n = len(pi) + stacked = np.array([np.asarray(joint[symbol], dtype=float) for symbol in symbols]).reshape(len(symbols), n, n) + return np.asarray(pi, dtype=float), symbols, stacked + + +def _observable_values(symbols: list[Hashable], observable: Observable | None) -> np.ndarray: + if observable is None: + values = [] + for symbol in symbols: + try: + values.append(float(symbol)) + except (TypeError, ValueError) as exc: + raise ValueError(f"symbol {symbol!r} is not numeric; pass an explicit observable") from exc + return np.asarray(values, dtype=float) + if isinstance(observable, Mapping): + missing = [symbol for symbol in symbols if symbol not in observable] + if missing: + raise ValueError(f"observable has no value for symbols {missing!r}") + return np.asarray([float(observable[symbol]) for symbol in symbols], dtype=float) + return np.asarray([float(observable(symbol)) for symbol in symbols], dtype=float) + + +def _correlation_operators( + hmm: Any, observable: Observable | None +) -> tuple[np.ndarray, np.ndarray, np.ndarray, float, float]: + """Return ``T``, ````, ``E[f(X_0)]`` and ``E[f(X_0)^2]``.""" + pi, symbols, stacked = _process_tensors(hmm) + values = _observable_values(symbols, observable) + ones = np.ones(len(pi)) + transition = stacked.sum(axis=0) + omega = np.tensordot(values, stacked, axes=1) + symbol_mass = np.einsum("i,kij,j->k", pi, stacked, ones) + return transition, pi @ omega, omega @ ones, float(values @ symbol_mass), float(values**2 @ symbol_mass) + + +def autocorrelation( + hmm: Any, + max_lag: int, + *, + observable: Observable | None = None, + centered: bool = False, +) -> np.ndarray: + """Return the autocorrelation ``γ(τ) = E[f(X_0) f(X_τ)]`` for ``τ = 0..max_lag``. + + Computed exactly as :math:`\\gamma(0) = E[f(X_0)^2]` and + :math:`\\gamma(\\tau) = \\langle \\pi | \\Omega \\, T^{\\tau-1} \\, \\Omega | + \\mathbf{1} \\rangle` for :math:`\\tau \\geq 1`, where :math:`\\Omega = \\sum_x + f(x) T^{(x)}` (Riechers & Crutchfield :cite:`riechers2017spectral1`, Eq. (3)), + under the stationary state law. + + Parameters + ---------- + hmm + Hidden Markov model (anything with ``to_mealy()``) generating the process. + max_lag + Largest lag ``τ`` returned. + observable + Real value ``f(x)`` of each symbol, as a mapping or callable. Defaults to + ``float(symbol)``, so string symbols ``"0"``, ``"1"`` map to 0, 1. + centered + Return the autocovariance ``γ(τ) − E[f(X_0)]²`` instead. + + Returns + ------- + numpy.ndarray + Array of length ``max_lag + 1`` indexed by ``τ``. + """ + if max_lag < 0: + raise ValueError("max_lag must be nonnegative") + transition, row, col, mean, second_moment = _correlation_operators(hmm, observable) + gamma = np.empty(max_lag + 1) + gamma[0] = second_moment + for lag in range(1, max_lag + 1): + gamma[lag] = row @ col + row = row @ transition + if centered: + gamma -= mean**2 + return gamma + + +def _peripheral_parts(transition: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Return ``(Q, X)``: the spectral projector onto ``|λ| = 1`` and ``X = (TQ)^#``. + + ``Q`` projects onto the invariant subspace of the unit-modulus eigenvalues + along the complementary invariant subspace; ``X`` inverts ``T`` on the range + of ``Q`` and vanishes on its kernel. + """ + + def peripheral(z: complex) -> bool: + return abs(z) > 1.0 - _PERIPHERAL_TOL + + _, right, dim = scipy.linalg.schur(transition.astype(complex), output="complex", sort=peripheral) + _, left, left_dim = scipy.linalg.schur(transition.T.astype(complex), output="complex", sort=peripheral) + if dim != left_dim: + raise np.linalg.LinAlgError("inconsistent peripheral spectrum") + u = right[:, :dim] + v_h = left[:, :dim].T + projector = u @ np.linalg.solve(v_h @ u, v_h) + inverse = u @ np.linalg.solve(v_h @ transition @ u, v_h) + return projector, inverse + + +def power_spectrum( + hmm: Any, + omegas: Any, + *, + observable: Observable | None = None, +) -> np.ndarray: + """Return the continuous part ``P_c(ω)`` of the power spectrum at angular frequencies ``omegas``. + + The power spectrum is :math:`P(\\omega) = \\lim_{N \\to \\infty} \\frac{1}{N} + E \\left| \\sum_{t=0}^{N-1} f(X_t) e^{-i \\omega t} \\right|^2 + = \\sum_{\\tau \\in \\mathbb{Z}} \\gamma(\\tau) e^{-i \\omega \\tau}`. + For a finite HMM it splits into a continuous part plus Dirac deltas at the + arguments of the eigenvalues of :math:`T` on the unit circle. The continuous + part is (Riechers & Crutchfield :cite:`riechers2017spectral1`, Eq. (5)) + + .. math:: P_c(\\omega) = E[f(X_0)^2] + 2 \\, \\mathrm{Re} \\, \\langle \\pi | + \\Omega (e^{i \\omega} I - T)^{-1} \\Omega | \\mathbf{1} \\rangle . + + The resolvent is singular at the delta frequencies (always at + :math:`\\omega = 0`, from :math:`\\lambda = 1`), so it is evaluated with the + peripheral spectrum projected out. With :math:`Q` the spectral projector + onto the unit-modulus eigenvalues of :math:`T`, + + .. math:: P_c(\\omega) = E[f(X_0)^2] - \\langle \\pi | \\Omega (TQ)^{\\#} + \\Omega | \\mathbf{1} \\rangle + 2 \\, \\mathrm{Re} \\, \\langle \\pi | + \\Omega (e^{i \\omega} I - T(I - Q))^{-1} (I - Q) \\Omega | + \\mathbf{1} \\rangle , + + which equals the formula above away from the delta frequencies and is its + continuous extension onto them (``P_c(0)`` is the Green--Kubo limit). The + middle term is the periodic part of :math:`\\gamma` at lag zero; it contains + :math:`E[f]^2`, so ``P_c`` is the same for ``f`` and ``f − E[f]``. + + The delta peaks -- the DC term :math:`2\\pi E[f]^2 \\delta(\\omega)` and, for + processes with periodic components, peaks at the other roots of unity in + the spectrum of :math:`T` -- are **excluded**. A purely periodic process + has ``P_c ≡ 0``. + + Parameters + ---------- + hmm + Hidden Markov model generating the process. + omegas + Angular frequencies (scalar or array-like, radians per time step). + observable + Real value of each symbol; see :func:`autocorrelation`. + + Returns + ------- + numpy.ndarray + Real array with the shape of ``omegas``. + """ + transition, row, col, _mean, second_moment = _correlation_operators(hmm, observable) + frequencies = np.asarray(omegas, dtype=float) + projector, inverse = _peripheral_parts(transition) + eye = np.eye(len(row)) + damped = transition @ (eye - projector) + damped_col = (eye - projector) @ col + baseline = second_moment - float(np.real(row @ inverse @ col)) + spectrum = np.empty(frequencies.size) + for k, omega in enumerate(frequencies.ravel()): + resolvent_col = np.linalg.solve(np.exp(1j * omega) * eye - damped, damped_col) + spectrum[k] = baseline + 2.0 * float(np.real(row @ resolvent_col)) + return spectrum.reshape(frequencies.shape) + + +def _mutual_information(joint: np.ndarray) -> float: + """Mutual information (bits) between the row and column variables of ``joint``.""" + joint = np.asarray(joint, dtype=float) + rows = joint.sum(axis=1, keepdims=True) + cols = joint.sum(axis=0, keepdims=True) + mask = joint > 0.0 + ratio = joint[mask] / (rows @ cols)[mask] + return max(float(np.sum(joint[mask] * np.log2(ratio))), 0.0) + + +def mutual_information_function(hmm: Any, max_lag: int) -> np.ndarray: + """Return the mutual information function ``I[X_0 : X_τ]`` (bits) for ``τ = 0..max_lag``. + + Uses the exact lagged joint law :math:`P(X_0 = x, X_\\tau = y) = \\langle \\pi | + T^{(x)} T^{\\tau-1} T^{(y)} | \\mathbf{1} \\rangle` under the stationary state + law (cf. Riechers & Crutchfield :cite:`riechers2017spectral1`, §V.A). Entry + ``τ = 0`` is ``I[X_0 : X_0] = H[X_0]``. Unlike :func:`autocorrelation`, the + result depends only on which symbols are distinct, not on their values. + + Returns + ------- + numpy.ndarray + Array of length ``max_lag + 1`` indexed by ``τ``. + """ + if max_lag < 0: + raise ValueError("max_lag must be nonnegative") + pi, _symbols, stacked = _process_tensors(hmm) + ones = np.ones(len(pi)) + transition = stacked.sum(axis=0) + heads = np.einsum("i,kij->kj", pi, stacked) + tails = stacked @ ones + info = np.empty(max_lag + 1) + info[0] = _mutual_information(np.diag(heads @ ones)) + for lag in range(1, max_lag + 1): + info[lag] = _mutual_information(heads @ tails.T) + heads = heads @ transition + return info diff --git a/sofic/generators/epsilon_construction.py b/sofic/generators/epsilon_construction.py index 4a80e77..2163bb0 100644 --- a/sofic/generators/epsilon_construction.py +++ b/sofic/generators/epsilon_construction.py @@ -35,9 +35,9 @@ TransitionSignature = tuple[tuple[Any, int, Any], ...] -def build_epsilon_machine(hmm: HiddenMarkovModel) -> EpsilonMachine: +def build_epsilon_machine(hmm: HiddenMarkovModel, *, max_states: int = 10_000) -> EpsilonMachine: hmm = hmm.to_mealy() - presentation = _unifilar_presentation(hmm) + presentation = _unifilar_presentation(hmm, max_states=max_states) stationary = presentation.stationary_distribution() if is_zero(array_sum(stationary)): raise StochasticValidationError("generator must have a stationary distribution") @@ -47,10 +47,10 @@ def build_epsilon_machine(hmm: HiddenMarkovModel) -> EpsilonMachine: return _quotient_machine(presentation, partitions, stationary, constraints=constraints) -def _unifilar_presentation(hmm: MealyHMM) -> MealyHMM: +def _unifilar_presentation(hmm: MealyHMM, *, max_states: int = 10_000) -> MealyHMM: """Return a row-unifilar presentation generating the same process.""" if not hmm.is_unifilar(): - hmm = hmm.mixed_state_presentation() + hmm = hmm.mixed_state_presentation(max_states=max_states) if isinstance(hmm, MixedStatePresentation) and hmm.recurrent_states: return hmm.to_recurrent() return hmm diff --git a/sofic/generators/epsilon_machine.py b/sofic/generators/epsilon_machine.py index 5a5055b..24efed7 100644 --- a/sofic/generators/epsilon_machine.py +++ b/sofic/generators/epsilon_machine.py @@ -41,10 +41,16 @@ def validate(self) -> None: self._check_unifilar() @classmethod - def from_hmm(cls, hmm: MealyHMM | MooreHMM, **kwargs: Any) -> EpsilonMachine: + def from_hmm(cls, hmm: MealyHMM | MooreHMM, *, max_states: int = 10_000) -> EpsilonMachine: + """Build the ε-machine of ``hmm``'s process. + + ``max_states`` caps the mixed-state enumeration used for non-unifilar + input; exceeding it raises + :class:`~sofic.exceptions.MixedStateExplosionError`. + """ from sofic.generators.epsilon_construction import build_epsilon_machine - return build_epsilon_machine(hmm) + return build_epsilon_machine(hmm, max_states=max_states) @classmethod def from_sequence( diff --git a/sofic/generators/mealy.py b/sofic/generators/mealy.py index 82b7ef3..2777b56 100644 --- a/sofic/generators/mealy.py +++ b/sofic/generators/mealy.py @@ -147,11 +147,14 @@ def mixed_state_presentation( self, *, initial_mixed_state: MixedState | Mapping[Hashable, float] | Sequence[float] | None = None, + max_states: int = 10_000, ) -> MixedStatePresentation: """Build the mixed-state presentation (observer belief dynamics).""" from sofic.generators.mixed_state import MixedStatePresentation - return MixedStatePresentation.from_presentation(self, initial_mixed_state=initial_mixed_state) + return MixedStatePresentation.from_presentation( + self, initial_mixed_state=initial_mixed_state, max_states=max_states + ) def to_edge_machine(self, iterations: int = 1, style: int = 0) -> MealyHMM: from sofic.generators.edge_machine import hmm_to_edge_machine diff --git a/sofic/generators/measures.py b/sofic/generators/measures.py index 0bf1013..5b83468 100644 --- a/sofic/generators/measures.py +++ b/sofic/generators/measures.py @@ -2,14 +2,19 @@ from __future__ import annotations -from typing import Any +import warnings +from collections.abc import Iterator +from typing import Any, Literal, NamedTuple import numpy as np +from sofic.exceptions import MixedStateExplosionError from sofic.generators.base import HiddenMarkovModel, QuasiStochasticModel, StochasticModel from sofic.generators.markov import MarkovChain from sofic.graph import ATTR_EMISSION, ATTR_PROB +EntropyRateMethod = Literal["auto", "exact", "bounds", "blackwell"] + def require_dit(feature: str = "entropy measures") -> Any: """Import and return the :mod:`dit` package, or raise a helpful error. @@ -36,8 +41,14 @@ def dit_state_label(state: Any) -> Any: ``("A", "0", "A")``) is misread as multi-dimensional coordinate data and raises ``MissingDimensionsError``. Tuple states are encoded to a lossless string via :func:`sofic.generators.edge_machine.edge_state_label` (invert - with ``parse_edge_state_label``); scalar states pass through unchanged. + with ``parse_edge_state_label``); mixed states, which dit cannot sort, are + labeled by the ``repr`` of their belief; other scalar states pass through + unchanged. """ + from sofic.generators.mixed_state import MixedState + + if isinstance(state, MixedState): + return repr(state.belief) if isinstance(state, tuple): from sofic.generators.edge_machine import edge_state_label @@ -234,6 +245,285 @@ def entropy_rate_markov(chain: MarkovChain) -> Any: return _entropy_rate_from_transitions(chain, pi, idx) +class EntropyRateEstimate(NamedTuple): + """Monte Carlo entropy-rate estimate in bits per symbol and its batch-means standard error.""" + + estimate: float + stderr: float + + +def entropy_rate( + model: StochasticModel, + method: EntropyRateMethod = "auto", + *, + tol: float = 1e-6, + max_length: int = 30, + max_words: int = 2**20, + max_states: int = 1_000, + n_samples: int = 100_000, + burn_in: int = 1_000, + seed: int | np.random.Generator | None = None, +) -> Any: + r"""Shannon entropy rate :math:`h_\mu` in bits per symbol of any finite HMM presentation. + + ``method`` selects the algorithm: + + ``"exact"`` + Unifilar presentations use the closed form over their states + (:func:`entropy_rate_hmm`); visible Markov chains use + :func:`entropy_rate_markov`. Non-unifilar presentations are unifilarized + by mixed-state construction :cite:`Ellison2009` with at most + ``max_states`` beliefs, and the closed form is applied to its recurrent + component. Raises :class:`~sofic.exceptions.MixedStateExplosionError` + when the mixed states do not close, which is the generic case + :cite:`blackwell1957entropy,jurgens2021shannon`. + ``"bounds"`` + The midpoint of the sandwich + :math:`H[X_n \mid X_{0:n}, S_0] \le h_\mu \le H[X_n \mid X_{0:n}]` + (see :func:`entropy_rate_bounds`), with ``n`` increased until the gap is + at most ``tol``. A :class:`RuntimeWarning` is issued when the gap still + exceeds ``tol`` at ``n = max_length`` or once more than ``max_words`` + positive-probability words are enumerated. + ``"blackwell"`` + The stochastic mixed-state estimate of :func:`entropy_rate_blackwell` + (``n_samples``, ``burn_in``, ``seed``); call that function directly to + also obtain its standard error. + ``"auto"`` + ``"exact"`` when it succeeds, otherwise ``"bounds"``. The result is + deterministic: the stochastic estimator is never chosen implicitly. + + Exact results are sympy expressions for symbolic presentations; bounds and + Blackwell estimates are numeric only. + """ + if method not in ("auto", "exact", "bounds", "blackwell"): + raise ValueError(f"unknown entropy-rate method {method!r}; expected 'auto', 'exact', 'bounds', or 'blackwell'") + if method == "blackwell": + return entropy_rate_blackwell(model, n_samples=n_samples, burn_in=burn_in, seed=seed).estimate + if method in ("auto", "exact"): + try: + return _entropy_rate_exact(model, max_states=max_states) + except MixedStateExplosionError: + if method == "exact": + raise + return _entropy_rate_from_bounds(model, tol=tol, max_length=max_length, max_words=max_words) + + +def _entropy_rate_exact(model: StochasticModel, *, max_states: int) -> Any: + if isinstance(model, MarkovChain): + return entropy_rate_markov(model) + if not isinstance(model, HiddenMarkovModel): + raise TypeError(f"entropy rate requires a HiddenMarkovModel or MarkovChain, not {type(model)!r}") + is_unifilar = getattr(model, "is_unifilar", None) + if is_unifilar is not None and is_unifilar(): + return entropy_rate_hmm(model) + from sofic.generators.mixed_state_construction import build_mixed_state_presentation + + presentation = build_mixed_state_presentation(model.to_mealy(), max_states=max_states) + return entropy_rate_hmm(presentation.to_recurrent()) + + +def _entropy_rate_from_bounds(model: StochasticModel, *, tol: float, max_length: int, max_words: int) -> float: + if max_length < 0: + raise ValueError("max_length must be nonnegative") + pi, tensors = _numeric_symbol_tensors(model) + for n, (lower, upper, n_words) in enumerate(_iter_entropy_rate_bounds(pi, tensors)): + if upper - lower <= tol: + return 0.5 * (lower + upper) + if n >= max_length or n_words > max_words: + break + warnings.warn( + f"entropy-rate bounds did not converge to tol={tol:g} by n={n}: [{lower:.10g}, {upper:.10g}]; " + "returning the midpoint", + RuntimeWarning, + stacklevel=3, + ) + return 0.5 * (lower + upper) + + +def entropy_rate_bounds(model: StochasticModel, n: int) -> tuple[float, float]: + r"""Return ``(lower, upper)`` bounds on the entropy rate in bits per symbol. + + For a stationary HMM with hidden states :math:`S_t` and symbols :math:`X_t`, + + .. math:: + + H[X_n \mid X_{0:n}, S_0] \;\le\; h_\mu \;\le\; H[X_n \mid X_{0:n}], + + the lower bound is nondecreasing and the upper bound nonincreasing in + ``n``, and both converge to :math:`h_\mu` (Cover & Thomas, Thm. 4.5.1, there + stated for functions of Markov chains with 1-based indices + :cite:`Cover2006`). For edge-emitting presentations :math:`S_0` is the + state *before* :math:`X_0` is emitted; by the Markov property it screens + the future off from the infinite past. The upper bound is the block-entropy + difference :math:`H[X_{0:n+1}] - H[X_{0:n}]`, which for many HMMs converges + exponentially fast :cite:`Travers2013`. + + Both bounds are exact: every positive-probability word of length ``n + 1`` + is enumerated from the stationary state law, so the cost grows with the + number of such words (about :math:`2^{n h_\mu}`). + """ + if n < 0: + raise ValueError("n must be nonnegative") + pi, tensors = _numeric_symbol_tensors(model) + for index, (lower, upper, _n_words) in enumerate(_iter_entropy_rate_bounds(pi, tensors)): + if index == n: + return lower, upper + raise AssertionError("unreachable") + + +def _iter_entropy_rate_bounds(pi: np.ndarray, tensors: np.ndarray) -> Iterator[tuple[float, float, int]]: + """Yield ``(H[X_n | X_{0:n}, S_0], H[X_n | X_{0:n}], word count)`` for ``n = 0, 1, ...``. + + Rows of ``backward`` are :math:`T^{(w)} \\mathbf{1}` over positive-probability + words ``w``; prepending a symbol ``x`` maps a row ``b`` to :math:`T^{(x)} b`. + """ + n_states = len(pi) + backward = np.ones((1, n_states)) + block_entropy = 0.0 + state_block_entropy = _entropy_bits(pi) + while True: + backward = np.einsum("xij,wj->xwi", tensors, backward).reshape(-1, n_states) + joint = backward * pi + words = joint.sum(axis=1) + keep = words > 0.0 + backward, joint, words = backward[keep], joint[keep], words[keep] + next_block = _entropy_bits(words) + next_state_block = _entropy_bits(joint.ravel()) + yield next_state_block - state_block_entropy, next_block - block_entropy, len(words) + block_entropy, state_block_entropy = next_block, next_state_block + + +def entropy_rate_blackwell( + model: StochasticModel, + *, + n_samples: int = 100_000, + burn_in: int = 1_000, + seed: int | np.random.Generator | None = None, + n_batches: int = 20, +) -> EntropyRateEstimate: + r"""Estimate the entropy rate by a random walk on mixed states. + + Blackwell's formula :cite:`blackwell1957entropy` writes the entropy rate as + the average next-symbol uncertainty over the Blackwell measure + :math:`\mu_B` on beliefs (mixed states) :math:`\eta`, + + .. math:: + + h_\mu = \int H[X \mid \eta] \, d\mu_B(\eta), \qquad + \Pr(x \mid \eta) = \eta T^{(x)} \mathbf{1}, \qquad + \eta \mapsto \frac{\eta T^{(x)}}{\eta T^{(x)} \mathbf{1}} . + + Following Jurgens & Crutchfield (:cite:`jurgens2021shannon`, Eqs. 16-17), + the integral is replaced by a time average along a sampled mixed-state + trajectory started from the stationary state law. The first ``burn_in`` + steps are discarded; the standard error comes from ``n_batches`` + non-overlapping batch means, which accounts for autocorrelation. + + ``seed`` is an integer seed or a :class:`numpy.random.Generator`. + """ + if n_batches < 2: + raise ValueError("n_batches must be at least 2") + if n_samples < n_batches: + raise ValueError("n_samples must be at least n_batches") + pi, tensors = _numeric_symbol_tensors(model) + rng = seed if isinstance(seed, np.random.Generator) else np.random.default_rng(seed) + _beliefs, _symbols, probabilities = mixed_state_walk(pi, tensors, n_samples=n_samples, burn_in=burn_in, rng=rng) + entropies = _row_entropies_bits(probabilities) + return EntropyRateEstimate(float(entropies.mean()), batch_means_stderr(entropies, n_batches)) + + +def mixed_state_walk( + pi: np.ndarray, + tensors: np.ndarray, + *, + n_samples: int, + burn_in: int, + rng: np.random.Generator, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Sample the mixed-state process started from belief ``pi``. + + ``tensors`` stacks the symbol-labeled matrices :math:`T^{(x)}` along axis 0. + After discarding ``burn_in`` steps, returns the beliefs ``eta_t`` (one row + per step), the symbol indices ``x_t`` emitted from them, and the + next-symbol distributions ``Pr(x | eta_t)``. + """ + if burn_in < 0: + raise ValueError("burn_in must be nonnegative") + if n_samples < 1: + raise ValueError("n_samples must be positive") + n_symbols, n_states, _ = tensors.shape + beliefs = np.empty((n_samples, n_states)) + symbols = np.empty(n_samples, dtype=np.intp) + probabilities = np.empty((n_samples, n_symbols)) + uniforms = rng.random(burn_in + n_samples) + eta = np.asarray(pi, dtype=float) / float(np.sum(pi)) + for t, u in enumerate(uniforms): + masses = eta @ tensors + p = masses.sum(axis=1) + total = float(p.sum()) + if total <= 0.0: + raise ValueError("mixed-state walk reached a belief with no outgoing probability") + p /= total + x = min(int(np.searchsorted(np.cumsum(p), u, side="right")), n_symbols - 1) + if t >= burn_in: + row = t - burn_in + beliefs[row] = eta + symbols[row] = x + probabilities[row] = p + eta = masses[x] / float(masses[x].sum()) + return beliefs, symbols, probabilities + + +def batch_means_stderr(values: np.ndarray, n_batches: int) -> float: + """Standard error of the mean of an autocorrelated stationary series by non-overlapping batch means.""" + size = len(values) // n_batches + means = np.asarray(values[: size * n_batches], dtype=float).reshape(n_batches, size).mean(axis=1) + return float(means.std(ddof=1) / np.sqrt(n_batches)) + + +def _numeric_symbol_tensors(model: StochasticModel) -> tuple[np.ndarray, np.ndarray]: + """Return the stationary state law and the stacked ``T^(x)`` matrices as float arrays. + + A visible Markov chain emits its next state, so ``T^(j)`` keeps column ``j`` + of the transition matrix. + """ + if isinstance(model, MarkovChain): + from sofic.properties import transition_matrix + + idx = model.reindex() + transition, _states = transition_matrix(model, attr=ATTR_PROB, states=idx.states) + transition = np.asarray(transition, dtype=float) + pi = np.asarray(model.stationary_distribution(), dtype=float) + n_states = len(pi) + tensors = np.zeros((n_states, n_states, n_states)) + for j in range(n_states): + tensors[j, :, j] = transition[:, j] + return pi / pi.sum(), tensors + + from sofic.generators.matrices import emission_tensors + + pi, joint = emission_tensors(model, policy="stationary") + if not joint: + raise ValueError("model emits no symbols") + try: + pi_f = np.asarray(pi, dtype=float) + tensors = np.stack([np.asarray(matrix, dtype=float) for matrix in joint.values()]) + except TypeError as exc: + raise NotImplementedError("entropy-rate bounds and estimators need numeric probabilities") from exc + return pi_f / pi_f.sum(), tensors + + +def _entropy_bits(probabilities: np.ndarray) -> float: + positive = probabilities[probabilities > 0.0] + return float(-np.sum(positive * np.log2(positive))) + + +def _row_entropies_bits(probabilities: np.ndarray) -> np.ndarray: + with np.errstate(divide="ignore", invalid="ignore"): + terms = np.where(probabilities > 0.0, probabilities * np.log2(probabilities), 0.0) + return -terms.sum(axis=1) + + def collision_entropy(quasi_model: QuasiStochasticModel) -> float: """Second Renyi entropy rate (bits) from quasi transition matrices.""" matrices = quasi_model.symbol_matrices() diff --git a/sofic/generators/mixed_state.py b/sofic/generators/mixed_state.py index 79cb43a..c65805f 100644 --- a/sofic/generators/mixed_state.py +++ b/sofic/generators/mixed_state.py @@ -137,10 +137,11 @@ def from_presentation( hmm: MealyHMM, *, initial_mixed_state: MixedState | Mapping[Hashable, float] | Sequence[float] | None = None, + max_states: int = 10_000, ) -> MixedStatePresentation: from sofic.generators.mixed_state_construction import build_mixed_state_presentation - return build_mixed_state_presentation(hmm, initial_mixed_state=initial_mixed_state) + return build_mixed_state_presentation(hmm, initial_mixed_state=initial_mixed_state, max_states=max_states) def to_recurrent(self) -> MealyHMM: """Return the recurrent component with stationary initial weights. diff --git a/sofic/generators/predictive_rd.py b/sofic/generators/predictive_rd.py new file mode 100644 index 0000000..ccffa6b --- /dev/null +++ b/sofic/generators/predictive_rd.py @@ -0,0 +1,260 @@ +"""Predictive rate-distortion (causal information bottleneck) for finite ε-machines. + +Lossy predictive features :math:`R` compress the causal state :math:`S` while +keeping information about the length-:math:`L` future :math:`X_{0:L}`. They +solve the information bottleneck (Tishby, Pereira & Bialek +:cite:`tishby2000information`) + +.. math:: \\min_{q(r \\mid s)} \\; I[S; R] - \\beta \\, I[R; X_{0:L}], + +with the Markov chain :math:`R - S - X_{0:L}`. Because the causal state is a +sufficient statistic of the past, compressing :math:`S` is equivalent to +compressing the whole past (optimal causal inference, Still, Crutchfield & +Ellison :cite:`still2010optimal`; predictive rate-distortion, Marzen & +Crutchfield :cite:`Marzen2016`). The self-consistent equations are iterated +Blahut--Arimoto style (Cover & Thomas :cite:`Cover2006`, §10.8). +""" + +from __future__ import annotations + +import warnings +from dataclasses import dataclass +from typing import Any + +import numpy as np + +from sofic.generators.correlations import _mutual_information +from sofic.generators.matrices import emission_tensors + +_DIRECTION_DECIMALS = 12 + + +@dataclass(frozen=True) +class PredictiveRateDistortionCurve: + """Points of a predictive rate-distortion (causal information bottleneck) curve. + + Attributes + ---------- + beta + Trade-off parameters, ascending. + rate + Coding rate ``I[S; R]`` in bits at each ``beta``. + relevant_information + Retained predictive information ``I[R; X_{0:L}]`` in bits. + distortion + Lost predictive information ``I[S; X_{0:L}] − I[R; X_{0:L}]`` in bits. + future_length + Future length ``L`` used for the morphs ``P(X_{0:L} | S)``. + predictive_information + ``I[S; X_{0:L}]`` in bits, the supremum of ``relevant_information`` + (at most the excess entropy ``E``). + statistical_complexity + ``C_μ = H[S]`` in bits, the supremum of ``rate``. + """ + + beta: np.ndarray + rate: np.ndarray + relevant_information: np.ndarray + distortion: np.ndarray + future_length: int + predictive_information: float + statistical_complexity: float + + +def _future_classes(transition_stack: np.ndarray, columns: np.ndarray) -> np.ndarray: + """Prepend every symbol to each future class and merge classes with proportional morphs. + + ``columns[:, k]`` is ``P(w | S = s)`` summed over the words ``w`` in class ``k``. + Words with proportional columns induce the same posterior over ``S``, so + merging them leaves every quantity of the bottleneck ``R - S - X_{0:L}`` + unchanged; prepending is linear, so merged classes extend consistently. + """ + extended = np.concatenate([matrix @ columns for matrix in transition_stack], axis=1) + totals = extended.sum(axis=0) + extended = extended[:, totals > 0.0] + totals = totals[totals > 0.0] + keys = np.round(extended / totals, _DIRECTION_DECIMALS) + _unique, inverse = np.unique(keys, axis=1, return_inverse=True) + merged = np.zeros((extended.shape[0], int(inverse.max()) + 1)) + np.add.at(merged.T, inverse.ravel(), extended.T) + return merged + + +def _bottleneck_inputs( + machine: Any, + future_length: int | None, + *, + tol: float, + max_future_length: int, +) -> tuple[np.ndarray, np.ndarray, int]: + """Return ``p(s)``, the morph matrix ``P(y | s)`` over merged futures, and ``L``.""" + pi, joint = emission_tensors(machine, policy="stationary") + pi = np.asarray(pi, dtype=float) + stacked = np.array([np.asarray(matrix, dtype=float) for matrix in joint.values()]) + columns = np.ones((len(pi), 1)) + if future_length is not None: + if future_length < 0: + raise ValueError("future_length must be nonnegative") + for _ in range(future_length): + columns = _future_classes(stacked, columns) + return pi, columns, future_length + + excess = float(machine.excess_entropy()) + for length in range(max_future_length + 1): + if excess - _mutual_information(pi[:, None] * columns) <= tol: + return pi, columns, length + if length < max_future_length: + columns = _future_classes(stacked, columns) + warnings.warn( + f"I[S; X_0:L] did not reach E within tol={tol} by L={max_future_length}", + RuntimeWarning, + stacklevel=3, + ) + return pi, columns, max_future_length + + +def _iterate( + p_s: np.ndarray, + morphs: np.ndarray, + log_morphs: np.ndarray, + beta: float, + encoder: np.ndarray, + *, + max_iter: int, + atol: float, +) -> np.ndarray: + """Iterate the information-bottleneck self-consistent equations from ``encoder``.""" + negentropy = np.sum(morphs * log_morphs, axis=1, keepdims=True) + for _ in range(max_iter): + q_r = p_s @ encoder + joint_ry = (p_s[:, None] * encoder).T @ morphs + with np.errstate(divide="ignore", invalid="ignore"): + decoder = np.where(q_r[:, None] > 0.0, joint_ry / q_r[:, None], 0.0) + log_q_r = np.where(q_r > 0.0, np.log(np.maximum(q_r, 1e-300)), -np.inf) + divergence = negentropy - morphs @ np.log(np.maximum(decoder, 1e-300)).T + logits = log_q_r[None, :] - beta * divergence + logits -= logits.max(axis=1, keepdims=True) + updated = np.exp(logits) + updated /= updated.sum(axis=1, keepdims=True) + if np.max(np.abs(updated - encoder)) < atol: + return updated + encoder = updated + return encoder + + +def _information(p_s: np.ndarray, morphs: np.ndarray, encoder: np.ndarray) -> tuple[float, float]: + joint_sr = p_s[:, None] * encoder + return _mutual_information(joint_sr), _mutual_information(joint_sr.T @ morphs) + + +def predictive_rate_distortion( + machine: Any, + betas: Any = 50, + *, + beta_range: tuple[float, float] = (0.5, 500.0), + future_length: int | None = None, + tol: float = 1e-6, + max_future_length: int = 64, + restarts: int = 4, + max_iter: int = 2000, + seed: int | None = 0, +) -> PredictiveRateDistortionCurve: + """Return the predictive rate-distortion curve of a finite ε-machine. + + For each trade-off ``β`` the causal states ``S`` are compressed into + features ``R`` (with ``|R| = |S|``, which suffices) by minimizing + :math:`I[S; R] - \\beta \\, I[R; X_{0:L}]` (Tishby et al. + :cite:`tishby2000information`; Still et al. :cite:`still2010optimal`; + Marzen & Crutchfield :cite:`Marzen2016`). The future morphs + :math:`P(X_{0:L} = w \\mid S = s) = \\langle \\delta_s | T^{(w)} | \\mathbf{1} + \\rangle` are exact; futures with proportional morph columns are merged, + which is lossless for the bottleneck and keeps long futures cheap. + + The encoder is found by iterating the self-consistent equations + :math:`q(r \\mid s) \\propto q(r) \\exp(-\\beta D_{KL}[P(X_{0:L} \\mid s) \\,\\|\\, + q(X_{0:L} \\mid r)])` (Blahut--Arimoto, Cover & Thomas :cite:`Cover2006`, + §10.8). The ``β`` grid is annealed in ascending order: each ``β`` starts + from a perturbation of the previous solution, the one-to-one encoder + ``R = S`` (near the bifurcation from the collapsed ``R`` the nontrivial + branch has a small basin that random starts tend to miss), and + ``restarts`` random encoders; the candidate with the lowest Lagrangian is + kept. + + Parameters + ---------- + machine + A finite :class:`~sofic.generators.epsilon_machine.EpsilonMachine`. + Any other HMM is converted with ``EpsilonMachine.from_hmm`` first. + betas + Trade-off values, or an integer number of geometrically spaced values + over ``beta_range``. For ``β ≤ 1`` the optimum is the trivial ``R`` + (zero rate). + beta_range + ``(low, high)`` used when ``betas`` is an integer. + future_length + Future length ``L``. By default the smallest ``L`` with + ``E − I[S; X_{0:L}] ≤ tol`` is used, where ``E`` is the excess entropy; + ``I[S; X_{0:L}]`` increases to ``I[S; X_{0:∞}] = E``. + tol + Tolerance for choosing the default ``future_length``. + max_future_length + Cap on the default ``future_length`` (a ``RuntimeWarning`` is issued if + it is reached before the tolerance). + restarts + Random restarts per ``β``. + max_iter + Maximum self-consistent iterations per start. + seed + Seed for the random restarts. + + Returns + ------- + PredictiveRateDistortionCurve + """ + from sofic.generators.epsilon_machine import EpsilonMachine + + if not isinstance(machine, EpsilonMachine): + machine = EpsilonMachine.from_hmm(machine) + if isinstance(betas, int | np.integer): + beta_grid = np.geomspace(beta_range[0], beta_range[1], int(betas)) + else: + beta_grid = np.sort(np.asarray(betas, dtype=float).ravel()) + if np.any(beta_grid < 0.0): + raise ValueError("betas must be nonnegative") + + p_s, morphs_joint, length = _bottleneck_inputs(machine, future_length, tol=tol, max_future_length=max_future_length) + keep = p_s > 0.0 + p_s = p_s[keep] + morphs = morphs_joint[keep] + with np.errstate(divide="ignore"): + log_morphs = np.where(morphs > 0.0, np.log(np.maximum(morphs, 1e-300)), 0.0) + total = _mutual_information(p_s[:, None] * morphs) + complexity = _mutual_information(np.diag(p_s)) + + rng = np.random.default_rng(seed) + n = len(p_s) + encoder = np.full((n, n), 1.0 / n) + rates = np.empty(len(beta_grid)) + relevant = np.empty(len(beta_grid)) + for k, beta in enumerate(beta_grid): + warm = 0.9 * encoder + 0.1 * rng.dirichlet(np.ones(n), size=n) + starts = [warm, np.eye(n), *(rng.dirichlet(np.ones(n), size=n) for _ in range(restarts))] + best = None + for start in starts: + candidate = _iterate(p_s, morphs, log_morphs, beta, start, max_iter=max_iter, atol=1e-12) + rate, info = _information(p_s, morphs, candidate) + lagrangian = rate - beta * info + if best is None or lagrangian < best[0] - 1e-12: + best = (lagrangian, rate, info, candidate) + _, rates[k], relevant[k], encoder = best + + relevant = np.minimum(relevant, total) + return PredictiveRateDistortionCurve( + beta=beta_grid, + rate=rates, + relevant_information=relevant, + distortion=total - relevant, + future_length=length, + predictive_information=total, + statistical_complexity=complexity, + ) diff --git a/sofic/generators/relative_entropy_rate.py b/sofic/generators/relative_entropy_rate.py new file mode 100644 index 0000000..b034fac --- /dev/null +++ b/sofic/generators/relative_entropy_rate.py @@ -0,0 +1,350 @@ +"""Relative entropy (Kullback-Leibler divergence) rate between stationary processes. + +For stationary processes ``P`` and ``Q`` over a finite alphabet the relative +entropy rate, in bits per symbol, is + +.. math:: + + D(P \\| Q) = \\lim_{n \\to \\infty} \\tfrac{1}{n} D(P_{0:n} \\| Q_{0:n}), + +where ``P_{0:n}`` is the law of the length-``n`` block :cite:`Gray1990,Cover2006`. +Both processes are taken in their stationary laws (each generator's +:meth:`~sofic.generators.base.StochasticModel.stationary_distribution`). +""" + +from __future__ import annotations + +import math +from collections.abc import Hashable +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import networkx as nx +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 + + +@dataclass(frozen=True) +class RelativeEntropyRateBounds: + """Bounds on a relative entropy rate in bits per symbol. + + ``lower <= D(P || Q) <= upper``. ``block_estimate`` is + ``D(P_{0:n} || Q_{0:n}) / n`` and ``conditional_estimate`` is + ``D(P_{0:n} || Q_{0:n}) - D(P_{0:n-1} || Q_{0:n-1})``, the expected + divergence of the next-symbol predictions given the length ``n - 1`` past. + """ + + lower: float + upper: float + block_length: int + block_estimate: float + conditional_estimate: float + + +@dataclass(frozen=True) +class _Presentation: + states: tuple[Hashable, ...] + pi: np.ndarray + matrices: dict[Any, np.ndarray] + + +def _presentation(model: HiddenMarkovModel) -> _Presentation: + 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()} + return _Presentation(states, pi, matrices) + + +def _aligned(p: _Presentation, q: _Presentation) -> tuple[list[Any], dict[Any, np.ndarray]]: + symbols = sorted(set(p.matrices) | set(q.matrices), key=repr) + zeros = np.zeros((len(q.states), len(q.states))) + 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)) + + +def _unifilar_successors(matrices: dict[Any, np.ndarray]) -> dict[Any, dict[int, int]] | None: + successors: dict[Any, dict[int, int]] = {} + for symbol, matrix in matrices.items(): + row_targets: dict[int, int] = {} + for i, row in enumerate(matrix): + targets = np.flatnonzero(row > _TOL) + if len(targets) > 1: + return None + if len(targets) == 1: + row_targets[i] = int(targets[0]) + successors[symbol] = row_targets + return successors + + +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`). + + Otherwise ``q`` must be unifilar. Writing ``h_μ(P)`` for the entropy rate of + ``p``, + + .. math:: + + D(P \\| Q) = -h_\\mu(P) - \\lim_{n \\to \\infty} + \\mathbb{E}_P\\left[\\log_2 Q(X_n \\mid X_{0:n})\\right], + + which for a ``k``-step Markov ``Q`` is Gray's formula + :math:`-h_\\mu(P) - \\mathbb{E}_P[\\log_2 Q(X_k \\mid X_{0:k})]` + (:cite:`Gray1990`, Lemma 3.10, stated there in nats). For a unifilar ``q`` + the past determines the set of ``q`` states consistent with it. The + expectation is computed on the finite Markov chain of pairs + ``(p state, set of q states)`` started from ``p``'s stationary law and the + stationary support of ``q``: its Cesàro-limit law is built from the + absorption probabilities into, and stationary laws of, its closed classes. + In each closed class the ``q`` set must be a single state (``q`` is exactly + synchronized by ``P``-typical pasts), and ``Q(x | past)`` is then that + state's emission probability. + + ``h_μ(P)`` comes from :meth:`~sofic.generators.base.HiddenMarkovModel.entropy_rate` + (exact for unifilar ``p``). ``p`` should have a unique stationary law + (ergodic). + + Raises + ------ + NotImplementedError + 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): + return math.inf + pp, qp = _presentation(p), _presentation(q) + symbols, q_matrices = _aligned(pp, qp) + successors = _unifilar_successors(q_matrices) + if successors is None: + raise NotImplementedError("exact relative entropy rate needs a unifilar q; use relative_entropy_rate_bounds") + + start_set = frozenset(int(i) for i in np.flatnonzero(qp.pi > _TOL)) + nodes: dict[tuple[int, frozenset[int]], int] = {} + edges: list[tuple[int, int, float]] = [] + frontier = [(int(i), start_set) for i in np.flatnonzero(pp.pi > _TOL)] + for node in frontier: + nodes[node] = len(nodes) + while frontier: + state, subset = frontier.pop() + for symbol in symbols: + matrix = pp.matrices.get(symbol) + if matrix is None: + continue + following = frozenset(successors[symbol][r] for r in subset if r in successors[symbol]) + for target in np.flatnonzero(matrix[state] > _TOL): + if not following: + return math.inf + node = (int(target), following) + if node not in nodes: + nodes[node] = len(nodes) + frontier.append(node) + edges.append((nodes[(state, subset)], nodes[node], float(matrix[state, target]))) + + size = len(nodes) + chain = np.zeros((size, size)) + for source, target, prob in edges: + chain[source, target] += prob + initial = np.zeros(size) + for i in np.flatnonzero(pp.pi > _TOL): + initial[nodes[(int(i), start_set)]] = pp.pi[i] + limit = _cesaro_limit(chain, initial / initial.sum()) + + labels = {index: node for node, index in nodes.items()} + emit_p = {symbol: pp.matrices[symbol].sum(axis=1) for symbol in pp.matrices} + emit_q = {symbol: q_matrices[symbol].sum(axis=1) for symbol in symbols} + cross = 0.0 + for index in np.flatnonzero(limit > _TOL): + state, subset = labels[int(index)] + if len(subset) != 1: + raise NotImplementedError( + "q is not exactly synchronized by typical pasts of p; use relative_entropy_rate_bounds" + ) + (r,) = subset + for symbol, probs in emit_p.items(): + if probs[state] <= _TOL: + continue + if emit_q[symbol][r] <= _TOL: + return math.inf + cross -= limit[index] * probs[state] * math.log2(emit_q[symbol][r]) + return max(cross - float(p.entropy_rate()), 0.0) + + +def _cesaro_limit(chain: np.ndarray, initial: np.ndarray) -> np.ndarray: + graph = nx.DiGraph() + graph.add_nodes_from(range(len(chain))) + graph.add_edges_from(zip(*np.nonzero(chain), strict=True)) + condensation = nx.condensation(graph) + closed = [sorted(condensation.nodes[c]["members"]) for c in condensation.nodes if condensation.out_degree(c) == 0] + recurrent = {state for members in closed for state in members} + transient = [state for state in range(len(chain)) if state not in recurrent] + + if transient: + sub = chain[np.ix_(transient, transient)] + visits = np.linalg.solve((np.eye(len(transient)) - sub).T, initial[transient]) + limit = np.zeros(len(chain)) + for members in closed: + mass = float(initial[members].sum()) + if transient: + mass += float(visits @ chain[np.ix_(transient, members)].sum(axis=1)) + block = chain[np.ix_(members, members)] + system = np.vstack([block.T - np.eye(len(members)), np.ones(len(members))]) + rhs = np.zeros(len(members) + 1) + rhs[-1] = 1.0 + stationary = np.linalg.lstsq(system, rhs, rcond=None)[0] + limit[members] = mass * stationary + return limit + + +def relative_entropy_rate_bounds( + p: HiddenMarkovModel, q: HiddenMarkovModel, block_length: int +) -> RelativeEntropyRateBounds: + """Return bounds on ``D(P || Q)`` from exact enumeration of length-``n`` words. + + Valid for any finite HMMs ``p`` and ``q`` (``q`` need not be unifilar). With + ``Q_s`` the law of ``q`` started in state ``s`` and ``s`` ranging over the + stationary support of ``q``, let + + .. math:: + + L_n = \\mathbb{E}_P[-\\log_2 \\max_s Q_s(X_{0:n})], \\qquad + U_n = \\mathbb{E}_P[-\\log_2 \\min_s Q_s(X_{0:n})]. + + Because ``Q(uv)`` lies between ``Q(u) \\min_s Q_s(v)`` and + ``Q(u) \\max_s Q_s(v)``, ``L_n`` is superadditive and ``U_n`` subadditive, + so the cross-entropy rate ``c = \\lim -\\tfrac{1}{n}\\mathbb{E}_P \\log_2 Q(X_{0:n})`` + satisfies ``L_n / n \\le c \\le U_n / n``, and ``L_n / n`` increases to + ``c`` when ``q``'s stationary law has full support on its recurrent states. + ``U_n`` is infinite whenever some recurrent ``q`` state cannot emit a + ``P``-typical word; it is finite, and ``U_n / n`` decreases to ``c``, when + every ``q`` state can emit every word. + + The entropy rate is bracketed by the hidden-Markov bounds + ``H[X_{n-1} | X_{0:n-1}, S_0] \\le h_μ(P) \\le H[X_{n-1} | X_{0:n-1}]`` + (:cite:`Cover2006`, Theorem 4.5.1), which converge from both sides. So + ``lower = max(0, L_n/n - H[X_{n-1} | X_{0:n-1}])`` and + ``upper = U_n/n - H[X_{n-1} | X_{0:n-1}, S_0]``. + + Both bounds are :data:`math.inf` when ``P`` is not absolutely continuous + with respect to ``Q`` on finite words, in which case the rate is infinite. + Enumeration costs ``O(|A|^n)`` words. + """ + if block_length < 1: + raise ValueError("block_length must be positive") + n = block_length + 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) + q_live = np.flatnonzero(qp.pi > _TOL) + q_pi = qp.pi[q_live] + q_sub = {symbol: matrix[np.ix_(q_live, q_live)] for symbol, matrix in q_matrices.items()} + + block = np.zeros(n + 1) + cross = np.zeros(n + 1) + low = np.zeros(n + 1) + high = np.zeros(n + 1) + + def walk(alpha: np.ndarray, product: np.ndarray, length: int) -> None: + prob = float(alpha.sum()) + if length: + from_states = product.sum(axis=1) + block[length] -= prob * math.log2(prob) + cross[length] -= prob * math.log2(float(q_pi @ from_states)) + low[length] -= prob * math.log2(float(from_states.max())) + smallest = float(from_states.min()) + high[length] = math.inf if smallest <= 0.0 else high[length] - prob * math.log2(smallest) + if length == n: + return + for symbol in symbols: + matrix = pp.matrices.get(symbol) + if matrix is None: + continue + following = alpha @ matrix + if following.sum() > _TOL: + walk(following, product @ q_sub[symbol], length + 1) + + walk(pp.pi, np.eye(len(q_live)), 0) + + conditioned = np.zeros(n + 1) + for i in np.flatnonzero(pp.pi > _TOL): + start = np.zeros(len(pp.pi)) + start[i] = 1.0 + conditioned += pp.pi[i] * _block_entropies(start, pp.matrices, n) + + divergence = cross - block + h_upper = block[n] - block[n - 1] + h_lower = conditioned[n] - conditioned[n - 1] + return RelativeEntropyRateBounds( + lower=max(0.0, float(low[n] / n - h_upper)), + upper=float(high[n] / n - h_lower), + block_length=n, + block_estimate=float(divergence[n] / n), + conditional_estimate=float(divergence[n] - divergence[n - 1]), + ) + + +def _block_entropies(start: np.ndarray, matrices: dict[Any, np.ndarray], n: int) -> np.ndarray: + entropies = np.zeros(n + 1) + + def walk(alpha: np.ndarray, length: int) -> None: + if length: + prob = float(alpha.sum()) + entropies[length] -= prob * math.log2(prob) + if length == n: + return + for matrix in matrices.values(): + following = alpha @ matrix + if following.sum() > _TOL: + walk(following, length + 1) + + walk(start, 0) + return entropies diff --git a/sofic/shifts/__init__.py b/sofic/shifts/__init__.py index 1012474..4b9eb1c 100644 --- a/sofic/shifts/__init__.py +++ b/sofic/shifts/__init__.py @@ -1,5 +1,6 @@ """Symbolic shifts and subshifts.""" +from sofic.shifts.algorithms import periodic_points, zeta_function from sofic.shifts.base import SymbolicModel from sofic.shifts.covers import ( LeftFischerCover, @@ -16,12 +17,23 @@ iter_sofic_dyck_topologies, shift_to_dyck_graph_string, ) +from sofic.shifts.invariants import BowenFranksGroup, bowen_franks_group, jordan_form_away_from_zero from sofic.shifts.markov_dyck import MarkovDyckShift from sofic.shifts.product_alphabet_shift import ProductAlphabetShift from sofic.shifts.sft import ShiftOfFiniteType from sofic.shifts.sliding_block_code import SlidingBlockCode, full_shift from sofic.shifts.sofic import SoficShift from sofic.shifts.sofic_dyck import SoficDyckShift +from sofic.shifts.state_splitting import ( + in_amalgamate, + in_edges, + in_split, + in_split_matrices, + out_amalgamate, + out_edges, + out_split, + out_split_matrices, +) from sofic.shifts.textile import TextileSystem from sofic.shifts.tmc import TopologicalMarkovChain from sofic.shifts.wheeler import ( @@ -34,6 +46,19 @@ ) __all__ = [ + "BowenFranksGroup", + "bowen_franks_group", + "in_amalgamate", + "in_edges", + "in_split", + "in_split_matrices", + "jordan_form_away_from_zero", + "out_amalgamate", + "out_edges", + "out_split", + "out_split_matrices", + "periodic_points", + "zeta_function", "DyckGraphString", "count_dyck_graph_strings", "dyck_graph_string_to_shift", diff --git a/sofic/shifts/algorithms.py b/sofic/shifts/algorithms.py index d85cb88..9875ec7 100644 --- a/sofic/shifts/algorithms.py +++ b/sofic/shifts/algorithms.py @@ -2,10 +2,12 @@ from __future__ import annotations -from collections import deque +from collections import defaultdict, deque from collections.abc import Hashable, Iterator +from itertools import combinations from typing import Any +import networkx as nx import numpy as np from sofic.graph import ATTR_MULTIPLICITY, ATTR_SYMBOL @@ -100,3 +102,235 @@ def sofic_topological_entropy(model: SymbolicModel) -> float: matrix, _ = adjacency_matrix(right_resolve(model), multiplicity=False) return topological_entropy_from_matrix(matrix) + + +def integer_adjacency_matrix(model: SymbolicModel) -> tuple[np.ndarray, tuple[Hashable, ...]]: + """Return the exact adjacency matrix (Python ``int`` entries) and its state order. + + Each transition counts its ``multiplicity`` attribute, as for an edge shift. + The matrix has ``object`` dtype so powers never overflow. + """ + states = tuple(model.states()) + index = {state: i for i, state in enumerate(states)} + matrix = np.zeros((len(states), len(states)), dtype=object) + for transition in model.transitions(): + matrix[index[transition.source], index[transition.target]] += edge_multiplicity(transition.data) + return matrix, states + + +def edge_multiplicity(data: dict[str, Any]) -> int: + """Return the integer ``multiplicity`` attribute of a transition (default 1).""" + multiplicity = data.get(ATTR_MULTIPLICITY, 1) + if int(multiplicity) != multiplicity: + raise ValueError(f"edge multiplicity must be an integer, got {multiplicity!r}") + return int(multiplicity) + + +def periodic_points(shift: SymbolicModel, n: int) -> int: + """Return ``p_n``, the number of points ``x`` with ``sigma^n x = x``. + + * :class:`~sofic.shifts.tmc.TopologicalMarkovChain`: ``tr A^n`` for the + adjacency matrix with multiplicities (edge-shift semantics) + :cite:`LindMarcus1995` (Proposition 2.2.12). + * :class:`~sofic.shifts.sft.ShiftOfFiniteType` given by forbidden words of + length at most ``L``: ``tr A^n`` for the edge graph whose vertices are + the allowed ``(L - 1)``-blocks and whose edges are the allowed + ``L``-blocks, an edge shift conjugate to the SFT + :cite:`LindMarcus1995` (Theorem 2.3.2). + * Any other labeled presentation (sofic shifts, presentation-only SFTs): + Manning's formula ``p_n = sum_j (-1)^(j+1) tr A_j^n`` over the signed + subset matrices ``A_j`` of a right-resolving presentation + :cite:`LindMarcus1995` (Theorem 6.4.8). Counting cycles of the + presentation itself would overcount points with several presentations. + + Arithmetic is exact (Python integers). + + Parameters + ---------- + shift : SymbolicModel + A TMC, SFT, or sofic-shift presentation. + n : int + Period, ``n >= 1``. + + Returns + ------- + int + The number of points of period ``n`` (not necessarily least period). + """ + if n < 1: + raise ValueError("n must be a positive integer") + return sum(sign * _trace_power(matrix, n) for matrix, sign in _periodic_point_matrices(shift)) + + +def zeta_function(shift: SymbolicModel, t: Any = None) -> Any: + """Return the Artin-Mazur zeta function ``exp(sum_n p_n t^n / n)`` as a sympy rational function. + + For an edge shift ``X_A`` this is ``1 / det(I - tA)`` + :cite:`LindMarcus1995` (Theorem 6.4.6); for a sofic shift with signed + subset matrices ``A_j`` it is ``prod_j det(I - tA_j)^((-1)^j)`` + :cite:`LindMarcus1995` (Theorem 6.4.8). The result is reduced to lowest + terms with numerator and denominator normalized to constant term ``1``. + See :func:`periodic_points` for how each shift class is presented. + Requires sympy (``sofic[symbolic]``). + + Parameters + ---------- + shift : SymbolicModel + A TMC, SFT, or sofic-shift presentation. + t : sympy.Symbol, optional + Variable of the rational function; defaults to ``Symbol("t")``. + + Returns + ------- + sympy.Expr + ``numerator / denominator`` in ``t``. + """ + import sympy as sp + + t = sp.Symbol("t") if t is None else t + numerator = sp.Poly(1, t) + denominator = sp.Poly(1, t) + for matrix, sign in _periodic_point_matrices(shift): + coefficients = reciprocal_characteristic_polynomial(matrix) + factor = sp.Poly(list(reversed(coefficients)), t) + if sign > 0: + denominator *= factor + else: + numerator *= factor + common = sp.gcd(numerator, denominator) + numerator = sp.quo(numerator, common) + denominator = sp.quo(denominator, common) + if denominator.eval(0) < 0: + numerator, denominator = -numerator, -denominator + return numerator.as_expr() / denominator.as_expr() + + +def reciprocal_characteristic_polynomial(matrix: np.ndarray) -> list[int]: + """Return the coefficients ``[1, c_1, ..., c_r]`` of ``det(I - tA)`` in increasing powers of ``t``. + + Exact Faddeev-LeVerrier recursion on an integer matrix; ``det(I - tA)`` is + ``t^r chi_A(1/t)`` :cite:`LindMarcus1995` (proof of Theorem 6.4.6). + """ + a = np.asarray(matrix, dtype=object) + size = a.shape[0] + identity = np.eye(size, dtype=int).astype(object) + coefficients = [1] + m = np.zeros((size, size), dtype=object) + for k in range(1, size + 1): + m = a.dot(m) + coefficients[-1] * identity + quotient, remainder = divmod(-int(np.trace(a.dot(m))), k) + if remainder: + raise ValueError("matrix must have integer entries") + coefficients.append(quotient) + return coefficients + + +def _trace_power(matrix: np.ndarray, n: int) -> int: + if matrix.shape[0] == 0: + return 0 + return int(np.trace(np.linalg.matrix_power(matrix, n))) + + +def _periodic_point_matrices(shift: SymbolicModel) -> list[tuple[np.ndarray, int]]: + """Integer matrices ``M`` with signs ``s`` such that ``p_n = sum s * tr M^n``.""" + from sofic.shifts.sft import ShiftOfFiniteType + from sofic.shifts.sofic_dyck import SoficDyckShift + from sofic.shifts.tmc import TopologicalMarkovChain + + if isinstance(shift, SoficDyckShift): + raise TypeError("Dyck shifts are not sofic; periodic points are not supported") + if isinstance(shift, TopologicalMarkovChain): + return [(integer_adjacency_matrix(shift)[0], 1)] + if isinstance(shift, ShiftOfFiniteType) and shift._has_forbidden_word_spec: + return [(_forbidden_block_matrix(shift), 1)] + return signed_subset_matrices(shift) + + +def _forbidden_block_matrix(shift: Any) -> np.ndarray: + """Adjacency of the edge graph with allowed ``(L-1)``-blocks as vertices and ``L``-blocks as edges.""" + forbidden = shift.forbidden_words() + if () in forbidden: + return np.zeros((0, 0), dtype=object) + length = max((len(word) for word in forbidden), default=1) + + def allowed(word: tuple[Any, ...]) -> bool: + return not any(len(bad) <= len(word) and word[len(word) - len(bad) :] == bad for bad in forbidden) + + blocks: list[tuple[Any, ...]] = [()] + for _ in range(length): + blocks = [ + block + (symbol,) for block in blocks for symbol in shift.symbol_alphabet if allowed(block + (symbol,)) + ] + vertices = {block[:-1] for block in blocks} | {block[1:] for block in blocks} + index = {vertex: i for i, vertex in enumerate(vertices)} + matrix = np.zeros((len(vertices), len(vertices)), dtype=object) + for block in blocks: + matrix[index[block[:-1]], index[block[1:]]] += 1 + return matrix + + +def signed_subset_matrices(model: SymbolicModel) -> list[tuple[np.ndarray, int]]: + """Return the signed subset matrices ``(A_j, (-1)^(j+1))`` of a right-resolving presentation. + + Following :cite:`LindMarcus1995` (§6.4), the ``j``-th auxiliary graph has + the ``j``-element state subsets as vertices; a symbol ``a`` defined on + every state of ``I = {I_1 < ... < I_j}`` with distinct images gives an + edge to ``J = a(I)`` signed by the parity of the permutation + ``(a(I_1), ..., a(I_j))`` of ``J``. ``model`` is first right-resolved + (:func:`~sofic.shifts.wheeler.right_resolve`) and trimmed. Each matrix is + restricted to the strongly connected parts of its graph that carry a + cycle, which leaves every ``tr A_j^n`` and ``det(I - tA_j)`` unchanged; + empty matrices are dropped. + """ + from sofic.shifts.sofic import SoficShift + from sofic.shifts.wheeler import right_resolve + + presentation = ( + model + if isinstance(model, SoficShift) + else SoficShift(graph=model.graph.copy(), symbol_alphabet=model.symbol_alphabet) + ) + resolved = trim_transient(right_resolve(presentation)) + states = tuple(resolved.states()) + index = {state: i for i, state in enumerate(states)} + step: list[dict[Any, int]] = [{} for _ in states] + for transition in resolved.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is not None: + step[index[transition.source]][symbol] = index[transition.target] + symbols = {symbol for moves in step for symbol in moves} + + matrices: list[tuple[np.ndarray, int]] = [] + for size in range(1, len(states) + 1): + entries: dict[tuple[tuple[int, ...], tuple[int, ...]], int] = defaultdict(int) + for subset in combinations(range(len(states)), size): + for symbol in symbols: + images = [step[i].get(symbol) for i in subset] + if None in images or len(set(images)) < size: + continue + entries[subset, tuple(sorted(images))] += _permutation_sign(images) + matrix = _cyclic_part(entries) + if matrix.shape[0]: + matrices.append((matrix, (-1) ** (size + 1))) + return matrices + + +def _permutation_sign(sequence: list[Any]) -> int: + inversions = sum(1 for i, j in combinations(range(len(sequence)), 2) if sequence[i] > sequence[j]) + return -1 if inversions % 2 else 1 + + +def _cyclic_part(entries: dict[tuple[Hashable, Hashable], int]) -> np.ndarray: + graph = nx.DiGraph((source, target) for (source, target), value in entries.items() if value) + nodes = [ + node + for component in nx.strongly_connected_components(graph) + if len(component) > 1 or graph.has_edge(next(iter(component)), next(iter(component))) + for node in component + ] + index = {node: i for i, node in enumerate(nodes)} + matrix = np.zeros((len(nodes), len(nodes)), dtype=object) + for (source, target), value in entries.items(): + if source in index and target in index: + matrix[index[source], index[target]] += value + return matrix diff --git a/sofic/shifts/invariants.py b/sofic/shifts/invariants.py new file mode 100644 index 0000000..b914f2f --- /dev/null +++ b/sofic/shifts/invariants.py @@ -0,0 +1,147 @@ +"""Algebraic conjugacy invariants of edge shifts (Lind & Marcus §7.4).""" + +from __future__ import annotations + +from typing import Any, NamedTuple + +import numpy as np + +from sofic.shifts.algorithms import integer_adjacency_matrix, reciprocal_characteristic_polynomial +from sofic.shifts.base import SymbolicModel + + +class BowenFranksGroup(NamedTuple): + """``BF(A) = Z^r / Z^r (I - A)`` together with the sign of ``det(I - A)``. + + ``invariant_factors`` are the Smith-form diagonal entries of ``I - A`` + other than ``1``, in divisibility order, with ``0`` standing for a free + ``Z`` summand; ``()`` is the trivial group. + """ + + invariant_factors: tuple[int, ...] + det_sign: int + + +def bowen_franks_group(matrix: np.ndarray | SymbolicModel) -> BowenFranksGroup: + """Return the Bowen-Franks group of ``A`` and the sign of ``det(I - A)``. + + ``BF(A) ≅ Z_{d_1} ⊕ ... ⊕ Z_{d_r}`` where ``d_j`` are the diagonal entries + of the Smith form of ``I - A`` over ``Z`` :cite:`LindMarcus1995` + (Definition 7.4.15); it is an invariant of shift equivalence over ``Z`` + and hence of conjugacy (Theorem 7.4.17). ``det(I - A) = 1/zeta(1)`` is a + conjugacy invariant too, and its sign together with ``BF(A)`` is the + Parry-Sullivan / Bowen-Franks flow-equivalence invariant + :cite:`LindMarcus1995` (§13.6). Arithmetic is exact. + + Parameters + ---------- + matrix : array_like or TopologicalMarkovChain + Square integer matrix, or a shift whose adjacency matrix (with edge + multiplicities) is used. + + Returns + ------- + BowenFranksGroup + ``(invariant_factors, det_sign)``. + """ + a = _integer_matrix(matrix) + size = a.shape[0] + difference = [[int(i == j) - int(a[i, j]) for j in range(size)] for i in range(size)] + diagonal = _smith_diagonal(difference) + det = sum(reciprocal_characteristic_polynomial(a)) + return BowenFranksGroup(tuple(d for d in diagonal if d != 1), (det > 0) - (det < 0)) + + +def jordan_form_away_from_zero(matrix: np.ndarray | SymbolicModel) -> tuple[tuple[Any, int], ...]: + """Return ``J^x(A)``, the Jordan blocks of ``A`` with nonzero eigenvalue. + + ``J^x(A)`` is the Jordan form of the invertible part of ``A``, an invariant + of shift equivalence over ``Z`` :cite:`LindMarcus1995` (Definition 7.4.9, + Theorem 7.4.10). Block sizes come from exact ranks: for each irreducible + factor ``q != t`` of ``chi_A`` over ``Q``, ``d_k = nullity(q(A)^k) / deg q`` + and the number of blocks of size ``k`` is ``2 d_k - d_(k-1) - d_(k+1)`` + for each root of ``q`` (Exercise 7.4.5). Requires sympy + (``sofic[symbolic]``). + + Parameters + ---------- + matrix : array_like or TopologicalMarkovChain + Square integer matrix, or a shift whose adjacency matrix is used. + + Returns + ------- + tuple of (sympy.Expr, int) + ``(eigenvalue, block size)`` pairs, one per block, sorted canonically. + """ + import sympy as sp + + a = _integer_matrix(matrix) + size = a.shape[0] + m = sp.Matrix(size, size, lambda i, j: int(a[i, j])) + x = sp.Symbol("x") + _, factors = sp.Poly(m.charpoly(x).as_expr(), x).factor_list() + blocks: list[tuple[Any, int]] = [] + for q, multiplicity in factors: + if q.eval(0) == 0: + continue + q_at_a = sp.zeros(size, size) + for coefficient in q.all_coeffs(): + q_at_a = q_at_a * m + coefficient * sp.eye(size) + nullities = [0] + power = sp.eye(size) + while nullities[-1] < multiplicity: + power = power * q_at_a + nullities.append((size - power.rank()) // q.degree()) + nullities.append(nullities[-1]) + counts = {k: 2 * nullities[k] - nullities[k - 1] - nullities[k + 1] for k in range(1, len(nullities) - 1)} + for root in q.all_roots(): + blocks.extend((root, k) for k, count in counts.items() for _ in range(count)) + return tuple(sorted(blocks, key=lambda block: (sp.default_sort_key(block[0]), block[1]))) + + +def _integer_matrix(matrix: np.ndarray | SymbolicModel) -> np.ndarray: + if isinstance(matrix, SymbolicModel): + return integer_adjacency_matrix(matrix)[0] + a = np.asarray(matrix) + if a.size == 0: + return np.zeros((0, 0), dtype=object) + if a.ndim != 2 or a.shape[0] != a.shape[1]: + raise ValueError("matrix must be square") + if any(int(v) != v for v in a.flat): + raise ValueError("matrix must have integer entries") + return np.vectorize(int, otypes=[object])(a) + + +def _smith_diagonal(matrix: list[list[int]]) -> list[int]: + """Smith-form diagonal ``d_1 | d_2 | ... | d_r`` (``d_j >= 0``) of a square integer matrix.""" + m = [row[:] for row in matrix] + size = len(m) + diagonal: list[int] = [] + for t in range(size): + entries = [(abs(m[i][j]), i, j) for i in range(t, size) for j in range(t, size) if m[i][j]] + if not entries: + break + while True: + _, i, j = min(entries) + m[t], m[i] = m[i], m[t] + for row in m: + row[t], row[j] = row[j], row[t] + pivot = m[t][t] + for i in range(t + 1, size): + q = m[i][t] // pivot + m[i] = [a - q * b for a, b in zip(m[i], m[t], strict=True)] + for j in range(t + 1, size): + q = m[t][j] // pivot + for row in m: + row[j] -= q * row[t] + residue = [(i, j) for i in range(t + 1, size) for j in range(t + 1, size) if m[i][j] % pivot] + if any(m[i][t] for i in range(t + 1, size)) or any(m[t][j] for j in range(t + 1, size)): + pass + elif residue: + i = residue[0][0] + m[t] = [a + b for a, b in zip(m[t], m[i], strict=True)] + else: + diagonal.append(abs(pivot)) + break + entries = [(abs(m[i][j]), i, j) for i in range(t, size) for j in range(t, size) if m[i][j]] + return diagonal + [0] * (size - len(diagonal)) diff --git a/sofic/shifts/sft.py b/sofic/shifts/sft.py index e73c225..09bd45d 100644 --- a/sofic/shifts/sft.py +++ b/sofic/shifts/sft.py @@ -53,6 +53,26 @@ def topological_entropy(self) -> float: return sofic_topological_entropy(self) + def periodic_points(self, n: int) -> int: + """Number of points of period ``n``. + + Counted on the higher-block edge graph of the forbidden-word + specification, or by Manning's formula on a presentation-only SFT. See + :func:`~sofic.shifts.algorithms.periodic_points` :cite:`LindMarcus1995`. + """ + from sofic.shifts.algorithms import periodic_points + + return periodic_points(self, n) + + def zeta_function(self, t: Any = None) -> Any: + """Zeta function, the reciprocal of a polynomial, as a sympy expression. + + See :func:`~sofic.shifts.algorithms.zeta_function` :cite:`LindMarcus1995`. + """ + from sofic.shifts.algorithms import zeta_function + + return zeta_function(self, t) + def forbidden_words( self, *, diff --git a/sofic/shifts/sofic.py b/sofic/shifts/sofic.py index 94c0baa..937d316 100644 --- a/sofic/shifts/sofic.py +++ b/sofic/shifts/sofic.py @@ -22,6 +22,24 @@ def topological_entropy(self) -> float: return sofic_topological_entropy(self) + def periodic_points(self, n: int) -> int: + """Number of points of period ``n`` by Manning's signed-subset formula. + + See :func:`~sofic.shifts.algorithms.periodic_points` :cite:`LindMarcus1995`. + """ + from sofic.shifts.algorithms import periodic_points + + return periodic_points(self, n) + + def zeta_function(self, t: Any = None) -> Any: + """Zeta function ``prod_j det(I - tA_j)^((-1)^j)`` as a sympy rational function. + + See :func:`~sofic.shifts.algorithms.zeta_function` :cite:`LindMarcus1995`. + """ + from sofic.shifts.algorithms import zeta_function + + return zeta_function(self, t) + def parry_measure(self) -> MealyHMM: """Measure of maximal entropy (Parry measure) as a labeled Mealy HMM. diff --git a/sofic/shifts/state_splitting.py b/sofic/shifts/state_splitting.py new file mode 100644 index 0000000..911b37b --- /dev/null +++ b/sofic/shifts/state_splitting.py @@ -0,0 +1,287 @@ +"""State splitting and amalgamation of edge shifts (Lind & Marcus §2.4). + +Edges of a :class:`~sofic.shifts.tmc.TopologicalMarkovChain` are counted with +multiplicity, so an individual edge is named ``(source, target, key, copy)``: +the transition ``key`` between ``source`` and ``target`` and the copy index +``0 <= copy < multiplicity``. :func:`out_edges` and :func:`in_edges` list them. + +Splitting a state ``I`` into ``m`` parts names the new states ``(I, 0), ..., +(I, m - 1)``; states left out of the partition keep their names. +""" + +from __future__ import annotations + +from collections import Counter +from collections.abc import Collection, Hashable, Mapping, Sequence +from typing import Any + +import numpy as np + +from sofic.graph import ATTR_MULTIPLICITY, ATTR_SYMBOL +from sofic.shifts.algorithms import edge_multiplicity, integer_adjacency_matrix +from sofic.shifts.tmc import TopologicalMarkovChain + +Edge = tuple[Hashable, Hashable, int, int] +Partition = Mapping[Hashable, Sequence[Collection[Edge]]] + + +def out_edges(tmc: TopologicalMarkovChain, state: Hashable) -> list[Edge]: + """Return the edges leaving ``state``, one per unit of multiplicity.""" + return [ + (t.source, t.target, t.key, copy) + for t in tmc.transitions() + if t.source == state + for copy in range(edge_multiplicity(t.data)) + ] + + +def in_edges(tmc: TopologicalMarkovChain, state: Hashable) -> list[Edge]: + """Return the edges entering ``state``, one per unit of multiplicity.""" + return [ + (t.source, t.target, t.key, copy) + for t in tmc.transitions() + if t.target == state + for copy in range(edge_multiplicity(t.data)) + ] + + +def out_split(tmc: TopologicalMarkovChain, partitions: Partition) -> TopologicalMarkovChain: + """Return the out-split graph ``G[P]`` :cite:`LindMarcus1995` (Definition 2.4.3). + + Each state ``I`` in ``partitions`` has its out-edges partitioned into + nonempty parts ``E_I^0, ..., E_I^(m-1)``. An edge ``e`` from ``I`` to + ``J`` in part ``i`` becomes one edge from ``(I, i)`` to every copy of + ``J``: out-edges are split, in-edges are copied. Every copy keeps the + symbol of ``e``, so erasing the split is a 1-block conjugacy of edge shifts + :cite:`LindMarcus1995` (Theorem 2.4.10). + + Parameters + ---------- + tmc : TopologicalMarkovChain + The graph ``G``. + partitions : Mapping + ``{state: [part, ...]}`` where each part is a collection of + :func:`out_edges` of ``state``; the parts must partition them. + + Returns + ------- + TopologicalMarkovChain + The split graph. + """ + names, part_of = _split_names(tmc, partitions, out_edges) + return _build_split(tmc, names, part_of, outgoing=True) + + +def in_split(tmc: TopologicalMarkovChain, partitions: Partition) -> TopologicalMarkovChain: + """Return the in-split graph ``G[P]`` :cite:`LindMarcus1995` (Definition 2.4.7). + + Each state ``J`` in ``partitions`` has its in-edges partitioned into + nonempty parts. An edge ``e`` from ``I`` to ``J`` in part ``j`` becomes + one edge from every copy of ``I`` to ``(J, j)``: in-edges are split, + out-edges are copied. The edge shift is conjugate to that of ``tmc``. + + Parameters + ---------- + tmc : TopologicalMarkovChain + The graph ``G``. + partitions : Mapping + ``{state: [part, ...]}`` where each part is a collection of + :func:`in_edges` of ``state``; the parts must partition them. + + Returns + ------- + TopologicalMarkovChain + The split graph. + """ + names, part_of = _split_names(tmc, partitions, in_edges) + return _build_split(tmc, names, part_of, outgoing=False) + + +def out_split_matrices(tmc: TopologicalMarkovChain, partitions: Partition) -> tuple[np.ndarray, np.ndarray]: + """Return the division and edge matrices ``(D, E)`` of an out-splitting. + + ``D(I, J^j) = [I = J]`` and ``E(I^i, J) = |E_I^i ∩ E^J|``, so that + ``A_G = DE`` and ``A_{G[P]} = ED`` :cite:`LindMarcus1995` (Definition + 2.4.11, Theorem 2.4.12). Rows and columns follow ``tmc.states()`` and + ``out_split(tmc, partitions).states()``. + """ + names, part_of = _split_names(tmc, partitions, out_edges) + return _split_matrices(tmc, names, part_of, outgoing=True) + + +def in_split_matrices(tmc: TopologicalMarkovChain, partitions: Partition) -> tuple[np.ndarray, np.ndarray]: + """Return ``(D, E)`` with ``A_G = DE`` and ``A_{G[P]} = ED`` for an in-splitting. + + In-splitting ``G`` is out-splitting its transpose + :cite:`LindMarcus1995` (Exercise 2.4.9), so here ``E`` is the transposed + division matrix, ``E(J^j, J) = 1``, and ``D(I, J^j) = |E^j_J ∩ E_I|`` counts + the edges from ``I`` in part ``j`` of ``J``. Rows and columns follow + ``tmc.states()`` and ``in_split(tmc, partitions).states()``. + """ + names, part_of = _split_names(tmc, partitions, in_edges) + return _split_matrices(tmc, names, part_of, outgoing=False) + + +def out_amalgamate( + tmc: TopologicalMarkovChain, classes: Mapping[Hashable, Sequence[Hashable]] +) -> TopologicalMarkovChain: + """Merge each class of states into one state, inverting an out-splitting. + + A class can be out-amalgamated when its states have identical columns in + the adjacency matrix (every state sends the same number of edges to each + of them); then ``A' = ED`` and ``A = DE`` with ``D`` the class division + matrix :cite:`LindMarcus1995` (Definition 2.4.9, Theorem 2.4.14). The + merged state keeps the out-edges of all members and the in-edges of the + first listed member. + + Parameters + ---------- + tmc : TopologicalMarkovChain + The split graph. + classes : Mapping + ``{new_state: [state, ...]}``; unlisted states are kept as is. + + Returns + ------- + TopologicalMarkovChain + The amalgamated graph. + """ + return _amalgamate(tmc, classes, outgoing=True) + + +def in_amalgamate( + tmc: TopologicalMarkovChain, classes: Mapping[Hashable, Sequence[Hashable]] +) -> TopologicalMarkovChain: + """Merge each class of states into one state, inverting an in-splitting. + + Dual of :func:`out_amalgamate`: members must have identical rows, and the + merged state keeps the in-edges of all members and the out-edges of the + first listed member :cite:`LindMarcus1995` (Definition 2.4.9, Exercise 2.4.9). + """ + return _amalgamate(tmc, classes, outgoing=False) + + +def _split_names( + tmc: TopologicalMarkovChain, partitions: Partition, edges_of: Any +) -> tuple[dict[Hashable, list[Hashable]], dict[Edge, int]]: + states = list(tmc.states()) + names: dict[Hashable, list[Hashable]] = {state: [state] for state in states} + part_of: dict[Edge, int] = {} + for state, parts in partitions.items(): + if state not in names: + raise ValueError(f"unknown state {state!r}") + expected = set(edges_of(tmc, state)) + seen: set[Edge] = set() + for index, part in enumerate(parts): + part = set(part) + if not part: + raise ValueError(f"partition of {state!r} has an empty part") + if part & seen or not part <= expected: + raise ValueError(f"parts for {state!r} must be disjoint edges of that state") + seen |= part + part_of.update(dict.fromkeys(part, index)) + if seen != expected: + raise ValueError(f"parts for {state!r} do not cover all of its edges") + names[state] = [(state, index) for index in range(len(parts))] + new_names = [name for state in states for name in names[state]] + if len(set(new_names)) != len(new_names): + raise ValueError("split state names collide with existing states") + return names, part_of + + +def _part_counts(tmc: TopologicalMarkovChain, part_of: dict[Edge, int], transition: Any) -> Counter[int]: + return Counter( + part_of.get((transition.source, transition.target, transition.key, copy), 0) + for copy in range(edge_multiplicity(transition.data)) + ) + + +def _build_split( + tmc: TopologicalMarkovChain, + names: dict[Hashable, list[Hashable]], + part_of: dict[Edge, int], + *, + outgoing: bool, +) -> TopologicalMarkovChain: + result = type(tmc)(symbol_alphabet=tmc.symbol_alphabet) + for state in tmc.states(): + for name in names[state]: + result.graph.add_state(name) + for transition in tmc.transitions(): + attrs = {k: v for k, v in transition.data.items() if k not in (ATTR_SYMBOL, ATTR_MULTIPLICITY)} + symbol = transition.data.get(ATTR_SYMBOL) + for part, count in _part_counts(tmc, part_of, transition).items(): + if outgoing: + pairs = [(names[transition.source][part], target) for target in names[transition.target]] + else: + pairs = [(source, names[transition.target][part]) for source in names[transition.source]] + for source, target in pairs: + result.add_transition(source, target, symbol, **{ATTR_MULTIPLICITY: count, **attrs}) + return result + + +def _split_matrices( + tmc: TopologicalMarkovChain, + names: dict[Hashable, list[Hashable]], + part_of: dict[Edge, int], + *, + outgoing: bool, +) -> tuple[np.ndarray, np.ndarray]: + states = list(tmc.states()) + row = {state: i for i, state in enumerate(states)} + column = {name: j for j, name in enumerate(name for state in states for name in names[state])} + division = np.zeros((len(states), len(column)), dtype=int) + for state in states: + for name in names[state]: + division[row[state], column[name]] = 1 + edge = np.zeros((len(column), len(states)), dtype=int) + for transition in tmc.transitions(): + owner = transition.source if outgoing else transition.target + other = transition.target if outgoing else transition.source + for part, count in _part_counts(tmc, part_of, transition).items(): + edge[column[names[owner][part]], row[other]] += count + if outgoing: + return division, edge + return edge.T, division.T + + +def _amalgamate( + tmc: TopologicalMarkovChain, classes: Mapping[Hashable, Sequence[Hashable]], *, outgoing: bool +) -> TopologicalMarkovChain: + matrix, states = integer_adjacency_matrix(tmc) + index = {state: i for i, state in enumerate(states)} + merged = {state: state for state in states} + representative: dict[Hashable, Hashable] = {} + assigned: set[Hashable] = set() + for new_state, members in classes.items(): + members = list(members) + if not members: + raise ValueError(f"class {new_state!r} is empty") + for member in members: + if member not in index: + raise ValueError(f"unknown state {member!r}") + if member in assigned: + raise ValueError(f"state {member!r} is in more than one class") + assigned.add(member) + merged[member] = new_state + representative[new_state] = members[0] + lines = [matrix[:, index[m]] if outgoing else matrix[index[m], :] for m in members] + if any(list(line) != list(lines[0]) for line in lines): + kind = "columns" if outgoing else "rows" + raise ValueError(f"class {new_state!r} cannot be amalgamated: its {kind} differ") + if (set(states) - assigned) & set(classes): + raise ValueError("amalgamated state names collide with existing states") + + result = type(tmc)(symbol_alphabet=tmc.symbol_alphabet) + for state in states: + if not result.graph.has_state(merged[state]): + result.graph.add_state(merged[state]) + for transition in tmc.transitions(): + kept = transition.target if outgoing else transition.source + if merged[kept] in representative and representative[merged[kept]] != kept: + continue + attrs = {k: v for k, v in transition.data.items() if k != ATTR_SYMBOL} + result.add_transition( + merged[transition.source], merged[transition.target], transition.data.get(ATTR_SYMBOL), **attrs + ) + return result diff --git a/sofic/shifts/tmc.py b/sofic/shifts/tmc.py index b25c556..53a5009 100644 --- a/sofic/shifts/tmc.py +++ b/sofic/shifts/tmc.py @@ -41,6 +41,24 @@ def topological_entropy(self) -> float: return topological_entropy(self) + def periodic_points(self, n: int) -> int: + """Number of points of period ``n``: ``tr A^n`` counting edge multiplicities. + + See :func:`~sofic.shifts.algorithms.periodic_points` :cite:`LindMarcus1995`. + """ + from sofic.shifts.algorithms import periodic_points + + return periodic_points(self, n) + + def zeta_function(self, t: Any = None) -> Any: + """Zeta function ``1 / det(I - tA)`` as a sympy rational function. + + See :func:`~sofic.shifts.algorithms.zeta_function` :cite:`LindMarcus1995`. + """ + from sofic.shifts.algorithms import zeta_function + + return zeta_function(self, t) + def validate(self) -> None: super().validate() for transition in self.transitions(): diff --git a/tests/test_automata_decisions.py b/tests/test_automata_decisions.py new file mode 100644 index 0000000..4a8e001 --- /dev/null +++ b/tests/test_automata_decisions.py @@ -0,0 +1,146 @@ +"""Emptiness, shortest-witness, universality, and inclusion on finite automata.""" + +from __future__ import annotations + +from itertools import product + +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.automata import DFA, NFA +from sofic.graph import EPSILON +from sofic.testing.strategies import dfas, nfas + +ALPHABET = ("a", "b") +MAX_LENGTH = 6 + +small_nfas = nfas(alphabet=ALPHABET, max_states=4, max_transitions=8) +small_automata = st.one_of(small_nfas, dfas(alphabet=ALPHABET, max_states=3)) + + +def _words(max_length: int = MAX_LENGTH): + for length in range(max_length + 1): + yield from product(ALPHABET, repeat=length) + + +def _language(aut, max_length: int = MAX_LENGTH) -> set[tuple[str, ...]]: + return {word for word in _words(max_length) if aut.recognizes(word)} + + +def _nfa(initial, accepting, edges, alphabet=ALPHABET) -> NFA: + nfa = NFA( + input_alphabet=frozenset(alphabet), + initial_states=frozenset(initial), + accepting_states=frozenset(accepting), + ) + for source, target, _symbol in edges: + nfa.graph.add_state(source) + nfa.graph.add_state(target) + for state in set(initial) | set(accepting): + nfa.graph.add_state(state) + for source, target, symbol in edges: + nfa.add_transition(source, target, symbol) + return nfa + + +def test_empty_language_has_no_witness(): + nfa = _nfa({0}, {1}, [(0, 0, "a")]) + assert nfa.is_empty() + assert nfa.accepted_word() is None + + +def test_epsilon_reaches_acceptance_with_empty_word(): + nfa = _nfa({0}, {2}, [(0, 1, EPSILON), (1, 2, EPSILON)]) + assert nfa.accepted_word() == () + + +def test_shortest_witness_ignores_epsilon_cost(): + nfa = _nfa({0}, {3}, [(0, 1, EPSILON), (1, 3, "b"), (0, 2, "a"), (2, 3, "a")]) + assert nfa.accepted_word() == ("b",) + + +def test_universal_nfa_without_determinization(): + nfa = _nfa({0}, {0}, [(0, 0, "a"), (0, 0, "b")]) + assert nfa.is_universal() + assert not nfa.is_universal(alphabet=frozenset({"a", "b", "c"})) + + +def test_includes_detects_missing_word(): + sigma_star = _nfa({0}, {0}, [(0, 0, "a"), (0, 0, "b")]) + only_a = _nfa({0}, {0}, [(0, 0, "a")]) + assert sigma_star.includes(only_a) + assert not only_a.includes(sigma_star) + + +def test_dfa_methods(): + dfa = DFA(input_alphabet=frozenset(ALPHABET), initial_states=frozenset({0}), accepting_states=frozenset({1})) + dfa.graph.add_state(0) + dfa.graph.add_state(1) + dfa.add_transition(0, 1, "a") + dfa.add_transition(1, 1, "a") + assert dfa.accepted_word() == ("a",) + assert not dfa.is_universal() + assert dfa.includes(_nfa({0}, {1}, [(0, 1, "a"), (1, 1, "a")])) + + +@pytest.mark.parametrize("states", [3, 4]) +def test_universality_exponential_witness(states): + """``(a|b)* a (a|b)^{k}`` complements need words of length ``k + 1`` to reject.""" + edges = [(0, 0, "a"), (0, 0, "b"), (0, 1, "a")] + edges += [(i, i + 1, symbol) for i in range(1, states) for symbol in ALPHABET] + nfa = _nfa({0}, {states}, edges) + assert not nfa.is_universal() + assert nfa.union(nfa.complement()).is_universal() + + +@settings(max_examples=150, deadline=None) +@given(small_automata) +def test_is_empty_matches_brute_force(aut): + # With at most 4 states a non-empty language has a word of length < 4. + assert aut.is_empty() == (not _language(aut, max_length=4)) + + +@settings(max_examples=150, deadline=None) +@given(small_automata) +def test_accepted_word_is_accepted_and_shortest(aut): + word = aut.accepted_word() + language = _language(aut) + if word is None: + assert not language + return + assert aut.recognizes(word) + assert len(word) == min(len(w) for w in language) + + +@settings(max_examples=150, deadline=None) +@given(small_automata) +def test_is_universal_matches_complement_and_brute_force(aut): + universal = aut.is_universal() + assert universal == aut.complement().is_empty() + assert universal == (len(_language(aut)) == len(list(_words()))) + assert universal == aut.minimize(alphabet=frozenset(ALPHABET)).is_universal() + + +@settings(max_examples=150, deadline=None) +@given(small_automata, small_automata) +def test_includes_matches_brute_force_and_difference(left, right): + included = left.includes(right) + assert included == right.difference(left, alphabet=frozenset(ALPHABET)).is_empty() + assert included == (_language(right) <= _language(left)) + + +@settings(max_examples=100, deadline=None) +@given(small_automata, small_automata, small_automata) +def test_includes_is_a_preorder(a, b, c): + assert a.includes(a) + if a.includes(b) and b.includes(c): + assert a.includes(c) + + +@settings(max_examples=100, deadline=None) +@given(small_automata, small_automata) +def test_union_includes_both_operands(a, b): + union = a.union(b) + assert union.includes(a) + assert union.includes(b) diff --git a/tests/test_buchi.py b/tests/test_buchi.py index 1f8ea61..f7a5172 100644 --- a/tests/test_buchi.py +++ b/tests/test_buchi.py @@ -1,9 +1,11 @@ """Tests for Büchi automata.""" import pytest +from hypothesis import given from sofic.automata.buchi import BuchiAutomaton from sofic.graph import ATTR_SYMBOL +from sofic.testing.strategies import buchi_automata def _accepting_loop_ba() -> BuchiAutomaton: @@ -149,3 +151,40 @@ def test_lasso_acceptance_matches_brute_force(): for p in itertools.product("ab", repeat=pl): for loop in itertools.product("ab", repeat=ll): assert ba.accepts_lasso(p, loop) == _brute_lasso(ba, p, loop) + + +def test_buchi_emptiness_uses_omega_semantics(): + accepting = _accepting_loop_ba() + assert not accepting.is_empty() + prefix, loop = accepting.accepted_lasso() + assert accepting.accepts_lasso(prefix, loop) + + rejecting = _rejecting_loop_ba() + assert rejecting.accepted_lasso() is None + assert rejecting.is_empty() + + +def test_buchi_finite_word_decisions_are_refused(): + ba = _accepting_loop_ba() + for call in (ba.accepted_word, ba.is_universal, lambda: ba.includes(ba)): + with pytest.raises(NotImplementedError): + call() + + +def _short_lassos(alphabet, max_prefix=2, max_loop=3): + from itertools import product + + for p in range(max_prefix + 1): + for prefix in product(sorted(alphabet), repeat=p): + for q in range(1, max_loop + 1): + for loop in product(sorted(alphabet), repeat=q): + yield prefix, loop + + +@given(buchi_automata(allow_epsilon=True)) +def test_accepted_lasso_matches_lasso_oracle(ba): + witness = ba.accepted_lasso() + if witness is not None: + assert ba.accepts_lasso(*witness) + else: + assert not any(ba.accepts_lasso(p, l) for p, l in _short_lassos(ba.input_alphabet)) diff --git a/tests/test_complexity_dimension.py b/tests/test_complexity_dimension.py new file mode 100644 index 0000000..4d42bc2 --- /dev/null +++ b/tests/test_complexity_dimension.py @@ -0,0 +1,109 @@ +"""Statistical complexity dimension (Jurgens & Crutchfield, Chaos 2021).""" + +import math + +import numpy as np +import pytest + +from sofic.examples import even_process, golden_mean, nonunifilar_golden_mean, sns +from sofic.generators.complexity_dimension import ifs_lyapunov_dimension, statistical_complexity_dimension +from sofic.generators.mealy import MealyHMM + + +def _from_matrices(matrices: dict[str, np.ndarray]) -> MealyHMM: + hmm = MealyHMM(observation_alphabet=frozenset(matrices)) + for symbol, matrix in matrices.items(): + for i in range(matrix.shape[0]): + hmm.graph.add_state(i) + for j in range(matrix.shape[1]): + if matrix[i, j] > 0: + hmm.add_transition(i, j, symbol, float(matrix[i, j])) + return hmm + + +def cantor_machine(s: float) -> MealyHMM: + """Fig. 1/2: the mixed states form the middle-1/s Cantor set; both maps contract by 1/s.""" + return _from_matrices( + { + "0": np.array([[0.5, 0.0], [(s - 1) / (2 * s), 1 / (2 * s)]]), + "1": np.array([[1 / (2 * s), (s - 1) / (2 * s)], [0.0, 0.5]]), + } + ) + + +def sierpinski_machine(s: float = 2.0, a: float = 1 / 6) -> MealyHMM: + """Appendix D, Eq. (S1): the mixed-state attractor is the Sierpinski triangle.""" + b = 1 - a * s + return _from_matrices( + { + "0": np.array([[a, 0, a * (s - 1)], [0, a, a * (s - 1)], [0, 0, a * s]]), + "1": np.array( + [[b / 2, 0, 0], [b * (s - 1) / (2 * s), b / (2 * s), 0], [b * (s - 1) / (2 * s), 0, b / (2 * s)]] + ), + "2": np.array( + [[b / (2 * s), b * (s - 1) / (2 * s), 0], [0, b / 2, 0], [0, b * (s - 1) / (2 * s), b / (2 * s)]] + ), + } + ) + + +@pytest.mark.parametrize("s", [3.0, 4.0]) +def test_cantor_machine_dimension(s): + result = statistical_complexity_dimension(cantor_machine(s), n_samples=5_000, seed=0) + assert not result.finite + assert result.entropy_rate == pytest.approx(1.0) + np.testing.assert_allclose(result.lyapunov_exponents, [-math.log2(s)], rtol=1e-9) + assert result.dimension == pytest.approx(math.log(2) / math.log(s), rel=1e-9) + + +def test_sierpinski_machine_dimension(): + result = statistical_complexity_dimension(sierpinski_machine(), n_samples=5_000, seed=0) + assert result.entropy_rate == pytest.approx(math.log2(3)) + np.testing.assert_allclose(result.lyapunov_exponents, [-1.0, -1.0], rtol=1e-9) + assert result.dimension == pytest.approx(math.log2(3), rel=1e-9) + + +@pytest.mark.parametrize("factory", [golden_mean, even_process]) +def test_unifilar_processes_have_zero_dimension(factory): + result = statistical_complexity_dimension(factory(), n_samples=2_000, seed=0) + assert result.finite + assert result.dimension == 0.0 + assert result.entropy_rate == pytest.approx(factory().entropy_rate()) + + +def test_nonunifilar_presentation_with_finite_epsilon_machine_has_zero_dimension(): + result = statistical_complexity_dimension(nonunifilar_golden_mean(), n_samples=2_000, seed=0) + assert result.finite + assert result.dimension == 0.0 + assert result.entropy_rate == pytest.approx(2 / 3) + + +def test_countable_mixed_states_have_zero_dimension(): + result = statistical_complexity_dimension(sns(), n_samples=5_000, seed=0, max_states=200) + assert not result.finite + assert result.lyapunov_exponents[0] == -math.inf + assert result.dimension == 0.0 + + +def test_seed_determinism(): + hmm = cantor_machine(3.0) + first = statistical_complexity_dimension(hmm, n_samples=1_000, seed=5, max_states=None) + second = statistical_complexity_dimension(hmm, n_samples=1_000, seed=np.random.default_rng(5), max_states=None) + assert first.dimension == second.dimension + np.testing.assert_array_equal(first.lyapunov_exponents, second.lyapunov_exponents) + + +@pytest.mark.parametrize( + ("rate", "exponents", "expected"), + [ + (1.0, [-2.0], 0.5), + (math.log2(3), [-1.0, -1.0], math.log2(3)), + (3.0, [-1.0, -1.0], 2.0), + (0.0, [-1.0], 0.0), + (1.0, [-math.inf], 0.0), + (1.5, [-1.0, -math.inf], 1.0), + (1.0, [], 0.0), + ], +) +def test_ifs_lyapunov_dimension(rate, exponents, expected): + assert ifs_lyapunov_dimension(rate, np.array(exponents)) == pytest.approx(expected) diff --git a/tests/test_correlations.py b/tests/test_correlations.py new file mode 100644 index 0000000..7c0216b --- /dev/null +++ b/tests/test_correlations.py @@ -0,0 +1,192 @@ +"""Tests for closed-form correlations, power spectra and mutual information functions.""" + +from __future__ import annotations + +import math +from collections import defaultdict + +import numpy as np +import pytest +from hypothesis import assume, given, settings + +from sofic.examples import even_process, golden_mean, iid, nonunifilar_golden_mean, periodic +from sofic.generators import EpsilonMachine, MealyHMM +from sofic.generators.correlations import autocorrelation, mutual_information_function, power_spectrum +from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.testing.strategies import mealy_hmms +from tests.oracles import word_distribution + +OMEGAS = np.linspace(0.0, np.pi, 9) + + +def _relabeled(hmm: MealyHMM, symbol_map: dict | None = None) -> MealyHMM: + """Copy ``hmm`` with states renamed (and symbols remapped by ``symbol_map``).""" + symbol_map = symbol_map or {} + rename = {state: ("relabeled", i) for i, state in enumerate(reversed(list(hmm.states())))} + out = MealyHMM( + observation_alphabet=frozenset(symbol_map.get(s, s) for s in hmm.observation_alphabet), + initial_distribution={rename[s]: mass for s, mass in hmm.initial_distribution.items()}, + ) + for state in rename.values(): + out.graph.add_state(state) + for transition in hmm.transitions(): + symbol = transition.data[ATTR_EMISSION] + out.add_transition( + rename[transition.source], + rename[transition.target], + symbol_map.get(symbol, symbol), + transition.data[ATTR_PROB], + ) + return out + + +def _brute_force_lagged(hmm: MealyHMM, lag: int) -> dict[tuple, float]: + """``P(X_0, X_lag)`` by marginalizing the oracle's length ``lag + 1`` word law.""" + joint: dict[tuple, float] = defaultdict(float) + for word, prob in word_distribution(hmm, lag + 1).items(): + joint[(word[0], word[-1])] += prob + return joint + + +def _mi_bits(joint: dict[tuple, float]) -> float: + left: dict = defaultdict(float) + right: dict = defaultdict(float) + for (x, y), p in joint.items(): + left[x] += p + right[y] += p + return sum(p * math.log2(p / (left[x] * right[y])) for (x, y), p in joint.items() if p > 0) + + +# --------------------------------------------------------------------------- i.i.d. + + +def test_iid_correlations_are_trivial() -> None: + process = iid(2) + gamma = autocorrelation(process, 5) + assert gamma[0] == pytest.approx(0.5) + np.testing.assert_allclose(gamma[1:], 0.25) + np.testing.assert_allclose(autocorrelation(process, 5, centered=True)[1:], 0.0, atol=1e-15) + np.testing.assert_allclose(mutual_information_function(process, 5), [1.0, 0, 0, 0, 0, 0], atol=1e-12) + np.testing.assert_allclose(power_spectrum(process, OMEGAS), 0.25) + + +def test_iid_with_custom_observable_has_flat_spectrum_at_variance() -> None: + process = iid(("a", "b", "c")) + values = {"a": -1.0, "b": 0.0, "c": 2.0} + mean = sum(values.values()) / 3 + var = sum(v**2 for v in values.values()) / 3 - mean**2 + np.testing.assert_allclose(autocorrelation(process, 3, observable=values)[1:], mean**2) + np.testing.assert_allclose(power_spectrum(process, OMEGAS, observable=values), var) + np.testing.assert_allclose(power_spectrum(process, OMEGAS, observable=values.__getitem__), var) + + +def test_non_numeric_symbols_need_observable() -> None: + with pytest.raises(ValueError, match="not numeric"): + autocorrelation(iid(("a", "b")), 2) + with pytest.raises(ValueError, match="no value"): + autocorrelation(iid(("a", "b")), 2, observable={"a": 1.0}) + with pytest.raises(ValueError, match="nonnegative"): + mutual_information_function(iid(2), -1) + + +# --------------------------------------------------------------------------- periodic + + +def test_period_two() -> None: + process = periodic("01") + np.testing.assert_allclose(autocorrelation(process, 4), [0.5, 0.0, 0.5, 0.0, 0.5]) + np.testing.assert_allclose(mutual_information_function(process, 4), 1.0) + np.testing.assert_allclose(power_spectrum(process, OMEGAS), 0.0, atol=1e-12) + + +@pytest.mark.parametrize("word", ["001", "0011", "00101"]) +def test_periodic_processes_have_no_continuous_spectrum(word: str) -> None: + np.testing.assert_allclose(power_spectrum(periodic(word), OMEGAS), 0.0, atol=1e-12) + + +# --------------------------------------------------------------------------- golden mean + + +@pytest.mark.parametrize("p", [0.3, 0.5, 0.8]) +def test_golden_mean_closed_form(p: float) -> None: + process = golden_mean(p) + lam = p - 1.0 + mean = (1.0 - p) / (2.0 - p) + var = mean * (1.0 - mean) + lags = np.arange(10) + np.testing.assert_allclose(autocorrelation(process, 9, centered=True), var * lam**lags, atol=1e-14) + np.testing.assert_allclose(autocorrelation(process, 9), mean**2 + var * lam**lags, atol=1e-14) + lorentzian = var * (1.0 - lam**2) / (1.0 - 2.0 * lam * np.cos(OMEGAS) + lam**2) + np.testing.assert_allclose(power_spectrum(process, OMEGAS), lorentzian) + + +def test_golden_mean_mi_decays_with_squared_subdominant_eigenvalue() -> None: + p = 0.6 + info = mutual_information_function(golden_mean(p), 10) + ratios = info[6:] / info[5:-1] + np.testing.assert_allclose(ratios, (p - 1.0) ** 2, rtol=1e-2) + + +# --------------------------------------------------------------------------- oracles + + +@pytest.mark.parametrize("make", [lambda: golden_mean(0.4), even_process, lambda: nonunifilar_golden_mean()]) +def test_against_brute_force_word_enumeration(make) -> None: + process = make() + process.initial_distribution = dict(zip(process.reindex().states, process.stationary_distribution(), strict=True)) + gamma = autocorrelation(process, 4) + info = mutual_information_function(process, 4) + for lag in range(1, 5): + joint = _brute_force_lagged(process, lag) + assert gamma[lag] == pytest.approx(sum(float(x) * float(y) * p for (x, y), p in joint.items())) + assert info[lag] == pytest.approx(_mi_bits(joint), abs=1e-12) + + +@pytest.mark.parametrize("make", [lambda: golden_mean(0.3), even_process, nonunifilar_golden_mean]) +def test_wiener_khinchin(make) -> None: + process = make() + lags = 512 + cov = autocorrelation(process, lags, centered=True) + symmetric = np.concatenate([cov, cov[-2:0:-1]]) + estimate = np.real(np.fft.fft(symmetric)) + omegas = 2.0 * np.pi * np.arange(len(symmetric)) / len(symmetric) + np.testing.assert_allclose(power_spectrum(process, omegas), estimate, atol=1e-9) + + +# --------------------------------------------------------------------------- invariants + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(max_states=3)) +def test_mi_bounds(hmm: MealyHMM) -> None: + info = mutual_information_function(hmm, 6) + assert np.all(info >= 0.0) + assert np.all(info[1:] <= info[0] + 1e-9) + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(max_states=3)) +def test_invariant_under_state_and_symbol_relabeling(hmm: MealyHMM) -> None: + assume(hmm.is_irreducible()) + relabeled = _relabeled(hmm) + np.testing.assert_allclose(autocorrelation(relabeled, 6), autocorrelation(hmm, 6), atol=1e-10) + np.testing.assert_allclose(power_spectrum(relabeled, OMEGAS), power_spectrum(hmm, OMEGAS), atol=1e-8) + swapped = _relabeled(hmm, {"0": "1", "1": "0"}) + np.testing.assert_allclose(mutual_information_function(swapped, 6), mutual_information_function(hmm, 6), atol=1e-10) + + +@settings(max_examples=30, deadline=None) +@given(mealy_hmms(max_states=3)) +def test_invariant_under_epsilon_machine_presentation(hmm: MealyHMM) -> None: + assume(hmm.is_irreducible() and hmm.is_unifilar()) + machine = EpsilonMachine.from_hmm(hmm) + np.testing.assert_allclose(autocorrelation(machine, 6), autocorrelation(hmm, 6), atol=1e-10) + np.testing.assert_allclose(mutual_information_function(machine, 6), mutual_information_function(hmm, 6), atol=1e-10) + np.testing.assert_allclose(power_spectrum(machine, OMEGAS), power_spectrum(hmm, OMEGAS), atol=1e-8) + + +def test_nonunifilar_and_epsilon_machine_agree() -> None: + hmm = nonunifilar_golden_mean() + machine = EpsilonMachine.from_hmm(hmm) + np.testing.assert_allclose(autocorrelation(machine, 8), autocorrelation(hmm, 8), atol=1e-10) + np.testing.assert_allclose(power_spectrum(machine, OMEGAS), power_spectrum(hmm, OMEGAS), atol=1e-8) diff --git a/tests/test_epsilon_machine.py b/tests/test_epsilon_machine.py index c9f4a31..ce23850 100644 --- a/tests/test_epsilon_machine.py +++ b/tests/test_epsilon_machine.py @@ -260,3 +260,13 @@ def test_numerically_truncated_reverse_machine_warns(): assert not forward.reverse_epsilon_machine_is_finite() with pytest.warns(RuntimeWarning, match="finite approximation"): EpsilonMachine.from_time_reversed(forward) + + +def test_from_hmm_honors_max_states(): + import pytest + + from sofic.examples import sns + from sofic.exceptions import MixedStateExplosionError + + with pytest.raises(MixedStateExplosionError): + EpsilonMachine.from_hmm(sns(), max_states=3) diff --git a/tests/test_k_testable.py b/tests/test_k_testable.py new file mode 100644 index 0000000..3e6f354 --- /dev/null +++ b/tests/test_k_testable.py @@ -0,0 +1,127 @@ +import itertools + +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.automata import learn_dfa_k_testable +from sofic.automata.learning import learn_dfa_k_testable as learning_export + + +def _words(alphabet, max_length): + for length in range(max_length + 1): + yield from itertools.product(alphabet, repeat=length) + + +def _k_tss_member(word, k, prefixes, suffixes, factors, short): + window = k - 1 + if len(word) < window: + return word in short + return ( + word[:window] in prefixes + and word[len(word) - window :] in suffixes + and all(word[index : index + k] in factors for index in range(len(word) - k + 1)) + ) + + +def _definition(samples, k): + words = [tuple(word) for word in samples] + window = k - 1 + prefixes = {word[:window] for word in words if len(word) >= window} + suffixes = {word[len(word) - window :] for word in words if len(word) >= window} + factors = {word[index : index + k] for word in words for index in range(len(word) - k + 1)} + short = {word for word in words if len(word) < window} + return prefixes, suffixes, factors, short + + +_samples = st.lists(st.lists(st.sampled_from("ab"), max_size=6).map(tuple), max_size=6) + + +def test_exported_from_learning_package(): + assert learn_dfa_k_testable is learning_export + + +@settings(max_examples=100, deadline=None) +@given(_samples, st.integers(1, 4)) +def test_agrees_with_definition(samples, k): + dfa = learn_dfa_k_testable(samples, k, alphabet="ab") + sets = _definition(samples, k) + for word in _words("ab", 7): + assert dfa.recognizes(word) == _k_tss_member(word, k, *sets), word + + +@settings(max_examples=100, deadline=None) +@given(_samples, st.integers(1, 4)) +def test_contains_samples(samples, k): + dfa = learn_dfa_k_testable(samples, k) + assert all(dfa.recognizes(word) for word in samples) + + +@settings(max_examples=100, deadline=None) +@given(_samples, st.integers(1, 3)) +def test_monotone_in_k(samples, k): + coarse = learn_dfa_k_testable(samples, k, alphabet="ab") + fine = learn_dfa_k_testable(samples, k + 1, alphabet="ab") + for word in _words("ab", 7): + assert not fine.recognizes(word) or coarse.recognizes(word), word + + +@settings(max_examples=60, deadline=None) +@given( + st.integers(2, 3).flatmap( + lambda k: st.tuples( + st.just(k), + st.frozensets(st.tuples(*[st.sampled_from("ab")] * (k - 1))), + st.frozensets(st.tuples(*[st.sampled_from("ab")] * (k - 1))), + st.frozensets(st.tuples(*[st.sampled_from("ab")] * k)), + st.frozensets(st.lists(st.sampled_from("ab"), max_size=k - 2).map(tuple)), + ) + ) +) +def test_identifies_target_in_the_limit(target): + k, prefixes, suffixes, factors, short = target + language = [word for word in _words("ab", 11) if _k_tss_member(word, k, prefixes, suffixes, factors, short)] + dfa = learn_dfa_k_testable(language, k, alphabet="ab") + for word in _words("ab", 8): + assert dfa.recognizes(word) == _k_tss_member(word, k, prefixes, suffixes, factors, short), word + + +def test_characteristic_sample_for_ab_star(): + dfa = learn_dfa_k_testable(["", "ab", "abab"], 2) + for word in _words("ab", 8): + assert dfa.recognizes(word) == ("".join(word) in {"ab" * n for n in range(5)}) + + +def test_k_equal_one_is_star_of_occurring_symbols(): + dfa = learn_dfa_k_testable(["a", "aa"], 1, alphabet="ab") + assert dfa.recognizes(()) + assert dfa.recognizes(("a",) * 5) + assert not dfa.recognizes(("b",)) + assert len(list(dfa.states())) == 1 + + +def test_k_equal_one_empty_word_only(): + dfa = learn_dfa_k_testable([""], 1, alphabet="ab") + assert dfa.recognizes(()) + assert not dfa.recognizes(("a",)) + + +@pytest.mark.parametrize("k", [1, 2, 3]) +def test_empty_sample_gives_empty_language(k): + dfa = learn_dfa_k_testable([], k, alphabet="ab") + assert not any(dfa.recognizes(word) for word in _words("ab", 5)) + + +def test_short_words_are_kept_exactly(): + dfa = learn_dfa_k_testable(["a", "abba"], 3) + assert dfa.recognizes(("a",)) + assert not dfa.recognizes(()) + assert not dfa.recognizes(("b",)) + assert dfa.recognizes(tuple("abba")) + + +def test_invalid_arguments(): + with pytest.raises(ValueError, match="at least 1"): + learn_dfa_k_testable(["a"], 0) + with pytest.raises(ValueError, match="not in the alphabet"): + learn_dfa_k_testable(["ac"], 2, alphabet="ab") diff --git a/tests/test_mealy_hmm.py b/tests/test_mealy_hmm.py index d95ea62..3d61150 100644 --- a/tests/test_mealy_hmm.py +++ b/tests/test_mealy_hmm.py @@ -70,7 +70,7 @@ def test_unifilarity_check(): bad.validate() -def test_entropy_rate_requires_unifilar_presentation(): +def test_entropy_rate_hmm_requires_unifilar_but_method_dispatches(): hmm = MealyHMM( initial_distribution={"q0": 1.0}, observation_alphabet=frozenset({"0", "1"}), @@ -81,8 +81,11 @@ def test_entropy_rate_requires_unifilar_presentation(): hmm.add_transition("q0", "q1", "0", 0.5) hmm.add_transition("q1", "q1", "1", 1.0) + from sofic.generators.measures import entropy_rate_hmm + with pytest.raises(NotImplementedError, match="unifilar"): - hmm.entropy_rate() + entropy_rate_hmm(hmm) + assert hmm.entropy_rate() == pytest.approx(0.0) def test_mealy_without_alphabet_or_initial_distribution(): diff --git a/tests/test_nonunifilar_entropy_rate.py b/tests/test_nonunifilar_entropy_rate.py new file mode 100644 index 0000000..379dfbb --- /dev/null +++ b/tests/test_nonunifilar_entropy_rate.py @@ -0,0 +1,132 @@ +"""Entropy rates of non-unifilar HMM presentations: bounds, Blackwell estimates, and dispatch.""" + +import math +import warnings + +import numpy as np +import pytest +from hypothesis import given, settings + +from sofic.examples import golden_mean, nonunifilar_golden_mean, sns +from sofic.exceptions import MixedStateExplosionError +from sofic.generators.measures import entropy_rate, entropy_rate_blackwell, entropy_rate_bounds +from sofic.testing.strategies import markov_chains, mealy_hmms + +pytestmark = pytest.mark.filterwarnings("ignore:stationary distribution is not unique:RuntimeWarning") + + +def _sns_entropy_rate() -> float: + """SNS is a renewal process: ``0`` resynchronizes, and ``P(L = n) = (n - 1) / 2^n`` for ``n >= 2``.""" + n = np.arange(2, 400) + p = (n - 1) / 2.0**n + return float(-(p * np.log2(p)).sum() / (n * p).sum()) + + +SNS_RATE = _sns_entropy_rate() + + +def test_sns_reference_value(): + assert pytest.approx(0.677867181055153, abs=1e-12) == SNS_RATE + + +@pytest.mark.parametrize("n", [0, 1, 3, 8, 14]) +def test_bounds_bracket_sns_rate(n): + lower, upper = entropy_rate_bounds(sns(), n) + assert lower - 1e-12 <= SNS_RATE <= upper + 1e-12 + + +def test_bounds_are_monotone_and_converge_for_sns(): + pairs = [entropy_rate_bounds(sns(), n) for n in range(16)] + lowers, uppers = np.array(pairs).T + assert np.all(np.diff(lowers) >= -1e-12) + assert np.all(np.diff(uppers) <= 1e-12) + assert uppers[-1] - lowers[-1] < 1e-7 + + +def test_bounds_for_finite_epsilon_machine_presentation(): + lower, upper = entropy_rate_bounds(nonunifilar_golden_mean(), 10) + assert lower == pytest.approx(2 / 3, abs=1e-9) + assert upper == pytest.approx(2 / 3, abs=1e-9) + + +def test_bounds_reject_negative_length(): + with pytest.raises(ValueError, match="nonnegative"): + entropy_rate_bounds(sns(), -1) + + +def test_exact_method_uses_finite_mixed_states(): + hmm = nonunifilar_golden_mean() + assert not hmm.is_unifilar() + assert hmm.entropy_rate(method="exact") == pytest.approx(golden_mean(0.5).entropy_rate()) + assert hmm.entropy_rate() == pytest.approx(2 / 3) + + +def test_exact_method_raises_when_mixed_states_explode(): + with pytest.raises(MixedStateExplosionError): + sns().entropy_rate(method="exact", max_states=200) + + +def test_auto_falls_back_to_converged_bounds(): + with warnings.catch_warnings(): + warnings.simplefilter("error") + value = sns().entropy_rate(max_states=200) + assert value == pytest.approx(SNS_RATE, abs=1e-6) + + +def test_bounds_method_warns_when_not_converged(): + with pytest.warns(RuntimeWarning, match="did not converge"): + value = entropy_rate(sns(), "bounds", max_length=3) + lower, upper = entropy_rate_bounds(sns(), 3) + assert value == pytest.approx(0.5 * (lower + upper)) + + +def test_unknown_method_raises(): + with pytest.raises(ValueError, match="unknown entropy-rate method"): + sns().entropy_rate(method="nope") + + +def test_unifilar_dispatch_is_unchanged(): + eps = golden_mean(0.5) + assert eps.entropy_rate() == eps.entropy_rate(method="exact") == pytest.approx(2 / 3) + + +def test_blackwell_estimate_within_standard_errors(): + estimate, stderr = entropy_rate_blackwell(sns(), n_samples=40_000, burn_in=500, seed=1) + assert stderr > 0.0 + assert abs(estimate - SNS_RATE) < 4 * stderr + + +def test_blackwell_is_seed_deterministic(): + first = entropy_rate_blackwell(sns(), n_samples=2_000, seed=7) + second = entropy_rate_blackwell(sns(), n_samples=2_000, seed=np.random.default_rng(7)) + assert first == second + + +def test_blackwell_method_returns_estimate(): + value = sns().entropy_rate(method="blackwell", n_samples=2_000, seed=3) + assert value == entropy_rate_blackwell(sns(), n_samples=2_000, seed=3).estimate + + +def test_blackwell_on_unifilar_process_is_exact_per_step(): + estimate, _stderr = entropy_rate_blackwell(golden_mean(0.5), n_samples=20_000, burn_in=100, seed=0) + assert estimate == pytest.approx(2 / 3, abs=0.02) + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(max_states=3)) +def test_bounds_sandwich_property(hmm): + pairs = [entropy_rate_bounds(hmm, n) for n in range(6)] + lowers, uppers = np.array(pairs).T + assert np.all(lowers <= uppers + 1e-9) + assert np.all(np.diff(lowers) >= -1e-9) + assert np.all(np.diff(uppers) <= 1e-9) + + +@settings(max_examples=30, deadline=None) +@given(markov_chains(max_states=3)) +def test_markov_chain_bounds_match_exact_rate(chain): + exact = chain.entropy_rate() + lower, upper = entropy_rate_bounds(chain, 2) + assert lower == pytest.approx(exact, abs=1e-9) + assert upper == pytest.approx(exact, abs=1e-9) + assert not math.isnan(entropy_rate(chain)) diff --git a/tests/test_predictive_rd.py b/tests/test_predictive_rd.py new file mode 100644 index 0000000..f8747e7 --- /dev/null +++ b/tests/test_predictive_rd.py @@ -0,0 +1,110 @@ +"""Tests for predictive rate-distortion (causal information bottleneck) curves.""" + +from __future__ import annotations + +import numpy as np +import pytest + +from sofic.examples import even_process, golden_mean, iid, nemo_process, nonunifilar_golden_mean, periodic +from sofic.generators import EpsilonMachine +from sofic.generators.predictive_rd import PredictiveRateDistortionCurve, predictive_rate_distortion + +PROCESSES = { + "golden_mean": lambda: golden_mean(0.5), + "even": even_process, + "period3": lambda: periodic("001"), + "nemo": nemo_process, +} + + +@pytest.fixture(scope="module", params=sorted(PROCESSES)) +def curve(request) -> tuple[EpsilonMachine, PredictiveRateDistortionCurve]: + machine = EpsilonMachine.from_hmm(PROCESSES[request.param]()) + return machine, predictive_rate_distortion(machine, 30, beta_range=(0.25, 2000.0)) + + +def test_endpoints(curve) -> None: + machine, result = curve + assert result.statistical_complexity == pytest.approx(machine.statistical_complexity()) + assert result.predictive_information <= machine.excess_entropy() + 1e-9 + assert result.predictive_information == pytest.approx(machine.excess_entropy(), abs=1e-6) + small = result.beta <= 1.0 + np.testing.assert_allclose(result.rate[small], 0.0, atol=1e-9) + np.testing.assert_allclose(result.relevant_information[small], 0.0, atol=1e-9) + assert result.rate[-1] == pytest.approx(result.statistical_complexity, abs=1e-4) + assert result.relevant_information[-1] == pytest.approx(result.predictive_information, abs=1e-6) + np.testing.assert_allclose(result.distortion, result.predictive_information - result.relevant_information) + + +def test_monotone_and_bounded(curve) -> None: + _machine, result = curve + assert np.all(np.diff(result.beta) > 0) + assert np.all(np.diff(result.rate) >= -1e-7) + assert np.all(np.diff(result.relevant_information) >= -1e-7) + assert np.all(result.relevant_information <= result.rate + 1e-9) + assert np.all(result.rate <= result.statistical_complexity + 1e-9) + assert np.all(result.distortion >= -1e-12) + + +def test_each_point_optimizes_its_lagrangian(curve) -> None: + """Every computed point is feasible, so the one chosen at ``β`` must minimize ``I[S;R] − β I[R;Y]``.""" + _machine, result = curve + for beta, rate, info in zip(result.beta, result.rate, result.relevant_information, strict=True): + lagrangians = result.rate - beta * result.relevant_information + assert rate - beta * info <= lagrangians.min() + 1e-6 + + +def test_curve_is_concave(curve) -> None: + _machine, result = curve + rate, order = np.unique(np.round(result.rate, 6), return_index=True) + info = result.relevant_information[order] + if len(rate) < 3: + return + slopes = np.diff(info) / np.diff(rate) + assert np.all(np.diff(slopes) <= 1e-3) + assert np.all(slopes <= 1.0 + 1e-6) + + +def test_iid_curve_is_zero() -> None: + result = predictive_rate_distortion(iid(3), 10) + assert result.statistical_complexity == 0.0 + assert result.future_length == 0 + for values in (result.rate, result.relevant_information, result.distortion): + np.testing.assert_allclose(values, 0.0, atol=1e-12) + + +def test_hmm_is_converted_to_epsilon_machine() -> None: + from_hmm = predictive_rate_distortion(nonunifilar_golden_mean(), 20) + direct = predictive_rate_distortion(golden_mean(0.5), 20) + np.testing.assert_allclose(from_hmm.rate, direct.rate, atol=1e-8) + np.testing.assert_allclose(from_hmm.relevant_information, direct.relevant_information, atol=1e-8) + + +def test_seed_determinism() -> None: + machine = nemo_process() + first = predictive_rate_distortion(machine, 12, seed=7) + second = predictive_rate_distortion(machine, 12, seed=7) + np.testing.assert_array_equal(first.rate, second.rate) + np.testing.assert_array_equal(first.relevant_information, second.relevant_information) + + +def test_explicit_betas_and_future_length() -> None: + result = predictive_rate_distortion(golden_mean(0.5), [100.0, 0.5, 5.0], future_length=3) + np.testing.assert_array_equal(result.beta, [0.5, 5.0, 100.0]) + assert result.future_length == 3 + with pytest.raises(ValueError, match="nonnegative"): + predictive_rate_distortion(golden_mean(0.5), [-1.0]) + with pytest.raises(ValueError, match="future_length"): + predictive_rate_distortion(golden_mean(0.5), 3, future_length=-1) + + +def test_short_future_captures_less_than_excess_entropy() -> None: + machine = even_process() + short = predictive_rate_distortion(machine, [1000.0], future_length=2) + assert short.predictive_information < machine.excess_entropy() - 1e-3 + + +def test_warns_when_future_length_cap_is_hit() -> None: + with pytest.warns(RuntimeWarning, match="did not reach E"): + result = predictive_rate_distortion(even_process(), 3, max_future_length=2) + assert result.future_length == 2 diff --git a/tests/test_regex_parser.py b/tests/test_regex_parser.py new file mode 100644 index 0000000..0328c40 --- /dev/null +++ b/tests/test_regex_parser.py @@ -0,0 +1,161 @@ +import itertools +import re + +import pytest +from hypothesis import given, settings +from hypothesis import strategies as st + +from sofic.automata import DFA, NFA, Regex, RegexSyntaxError, automaton_to_regex, equivalent, parse_regex, regex_to_nfa +from sofic.testing.strategies import dfas, nfas + + +def _words(alphabet, max_length): + for length in range(max_length + 1): + yield from itertools.product(alphabet, repeat=length) + + +def _assert_same_words(left, right, alphabet, max_length): + for word in _words(sorted(alphabet, key=repr), max_length): + assert left.recognizes(word) == right.recognizes(word), word + + +@pytest.mark.parametrize( + "text", + ["(a", "a)", "(a|b", "a|b)", "*a", "a|*", "+", "?a", "a|", "|a", "|", "a||b", "a**", "a+?", "", "()", "(|a)"], +) +def test_parse_errors_raise_regex_syntax_error(text): + with pytest.raises(RegexSyntaxError) as info: + parse_regex(text) + assert isinstance(info.value, ValueError) + assert 0 <= info.value.position <= len(text) + + +@pytest.mark.parametrize("text", ["'ab", "a\\", ".", "[ab]", "a{2}", "^a$", "(?=a)"]) +def test_unsupported_or_unterminated_syntax_raises(text): + with pytest.raises(RegexSyntaxError): + parse_regex(text) + + +def test_parse_tree_is_unsimplified(): + assert parse_regex("ab|c*") == Regex( + "union", + ( + Regex("concat", (Regex("literal", ("a",)), Regex("literal", ("b",)))), + Regex("star", (Regex("literal", ("c",)),)), + ), + ) + assert parse_regex("'ab'") == Regex("literal", ("ab",)) + assert parse_regex("ε") == parse_regex("(?:)") == Regex("epsilon") + assert parse_regex("∅") == parse_regex("(?!)") == Regex("empty") + + +@pytest.mark.parametrize( + ("text", "python"), + [ + ("a", "a"), + ("ab|ba", "ab|ba"), + ("(ab)*", "(ab)*"), + ("(?:a|b)*abb", "(?:a|b)*abb"), + ("a+b?", "a+b?"), + ("(a|ε)b", "(a|)b"), + ("ε", ""), + ("∅", "(?!)"), + ("a∅|b", "a(?!)|b"), + ("(a*b*)*c", "(a*b*)*c"), + ("((a|b)(a|b))*", "((a|b)(a|b))*"), + ("a(?:)b", "ab"), + ], +) +def test_known_languages_against_python_re(text, python): + nfa = regex_to_nfa(text, alphabet={"a", "b", "c"}) + compiled = re.compile(python) + for word in _words("abc", 6): + assert nfa.recognizes(word) == bool(compiled.fullmatch("".join(word))), (text, word) + + +def test_multi_character_and_escaped_symbols(): + nfa = regex_to_nfa("'ab'(a|b)*\\*'it\\'s'") + assert nfa.input_alphabet == frozenset({"ab", "a", "b", "*", "it's"}) + assert nfa.recognizes(("ab", "*", "it's")) + assert nfa.recognizes(("ab", "a", "b", "*", "it's")) + assert not nfa.recognizes(("a", "b", "*", "it's")) + + +def test_thompson_shape(): + nfa = regex_to_nfa("(a|b)*c") + assert len(nfa.initial_states) == 1 + assert len(nfa.accepting_states) == 1 + assert len(list(nfa.states())) <= 2 * 6 + + +def test_alphabet_resolves_non_string_symbols(): + nfa = regex_to_nfa("0(1|'10')*", alphabet={0, 1, 10}) + assert nfa.input_alphabet == frozenset({0, 1, 10}) + assert nfa.recognizes((0, 10, 1)) + with pytest.raises(ValueError, match="not in"): + regex_to_nfa("2", alphabet={0, 1}) + + +def test_nfa_from_regex_delegates(): + assert equivalent(NFA.from_regex("(ab)*"), regex_to_nfa("(ab)*")) + + +def test_automaton_to_regex_quotes_multi_character_symbols(): + def single(symbol): + dfa = DFA(input_alphabet=frozenset({"a", "b", "ab"}), initial_states=frozenset({0})) + dfa.graph.add_state(0) + dfa.graph.add_state(1) + dfa.accepting_states = frozenset({1}) + if symbol == "a b": + dfa.graph.add_state(2) + dfa.add_transition(0, 2, "a") + dfa.add_transition(2, 1, "b") + else: + dfa.add_transition(0, 1, symbol) + return dfa + + assert automaton_to_regex(single("ab")) == "'ab'" + assert automaton_to_regex(single("a b")) == "ab" + + +def test_automaton_to_regex_single_character_output_is_python_re(): + dfa = DFA(input_alphabet=frozenset({"a", "b"}), initial_states=frozenset({0}), accepting_states=frozenset({0})) + dfa.graph.add_state(0) + dfa.graph.add_state(1) + dfa.add_transition(0, 1, "a") + dfa.add_transition(1, 0, "b") + pattern = automaton_to_regex(dfa) + assert pattern == "(?:(?:)|a(?:ba)*b)" + for word in _words("ab", 6): + assert bool(re.fullmatch(pattern, "".join(word))) == dfa.recognizes(word) + + +_SYMBOL_ALPHABETS = [ + ("0", "1"), + ("a", "ab", "b"), + ("ε", "∅", "'"), + ("(", "|", "x y", "\\", "a'b"), +] + + +@settings(max_examples=60, deadline=None) +@given(st.sampled_from(_SYMBOL_ALPHABETS).flatmap(lambda alphabet: nfas(alphabet=alphabet, max_states=3))) +def test_round_trip_nfa(nfa): + rebuilt = regex_to_nfa(automaton_to_regex(nfa), alphabet=nfa.input_alphabet) + _assert_same_words(nfa, rebuilt, nfa.input_alphabet, 4) + assert equivalent(nfa, rebuilt, nfa.input_alphabet) + + +@settings(max_examples=60, deadline=None) +@given(st.sampled_from(_SYMBOL_ALPHABETS[:3]).flatmap(lambda alphabet: dfas(alphabet=alphabet, max_states=3))) +def test_round_trip_dfa(dfa): + rebuilt = regex_to_nfa(automaton_to_regex(dfa), alphabet=dfa.input_alphabet) + _assert_same_words(dfa, rebuilt, dfa.input_alphabet, 4) + assert equivalent(dfa, rebuilt, dfa.input_alphabet) + + +@settings(max_examples=40, deadline=None) +@given(dfas(alphabet=(0, 1, 10), max_states=3)) +def test_round_trip_non_string_symbols(dfa): + rebuilt = regex_to_nfa(automaton_to_regex(dfa), alphabet=dfa.input_alphabet) + assert equivalent(dfa, rebuilt, dfa.input_alphabet) diff --git a/tests/test_relative_entropy_rate.py b/tests/test_relative_entropy_rate.py new file mode 100644 index 0000000..9a2cd6c --- /dev/null +++ b/tests/test_relative_entropy_rate.py @@ -0,0 +1,163 @@ +"""Relative entropy rate between stationary hidden Markov processes.""" + +from __future__ import annotations + +import math +from itertools import product + +import numpy as np +import pytest +from hypothesis import assume, 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 ( + _support_includes, + _support_nfa, + relative_entropy_rate, + relative_entropy_rate_bounds, +) +from sofic.testing.strategies import epsilon_machines, mealy_hmms + +TOL = 1e-9 + + +def _kl(p, q) -> float: + return sum(a * math.log2(a / b) for a, b in zip(p, q, strict=True) if a > 0) + + +def _noisy_hmm(eps: float = 0.2) -> MealyHMM: + """Non-unifilar two-state HMM in which every state can emit every symbol.""" + 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.6 * (1 - eps)) + hmm.add_transition("A", "B", "0", 0.4 * (1 - eps)) + hmm.add_transition("A", "A", "1", 0.6 * eps) + hmm.add_transition("A", "B", "1", 0.4 * eps) + hmm.add_transition("B", "A", "1", 0.3 * (1 - eps)) + hmm.add_transition("B", "B", "1", 0.7 * (1 - eps)) + hmm.add_transition("B", "A", "0", 0.3 * eps) + hmm.add_transition("B", "B", "0", 0.7 * eps) + return hmm + + +@pytest.mark.parametrize(("a", "b"), [(0.3, 0.6), (0.5, 0.1), (0.9, 0.5)]) +def test_iid_equals_per_symbol_kl(a, b): + expected = _kl([1 - a, a], [1 - b, b]) + assert relative_entropy_rate(bernoulli(a), bernoulli(b)) == pytest.approx(expected, abs=TOL) + + +@pytest.mark.parametrize(("p", "q"), [((0.4, 0.3), (0.2, 0.6)), ((0.1, 0.8), (0.5, 0.45))]) +def test_markov_chain_is_stationary_weighted_row_kl(p, q): + """Gray (1990), Lemma 3.10 with k = 1: ``sum_i pi_i D(P(.|i) || Q(.|i))``.""" + (p01, p10), (q01, q10) = p, q + pi = np.array([p10, p01]) / (p01 + p10) + rows_p = [[1 - p01, p01], [p10, 1 - p10]] + rows_q = [[1 - q01, q01], [q10, 1 - q10]] + expected = sum(pi[i] * _kl(rows_p[i], rows_q[i]) for i in range(2)) + actual = relative_entropy_rate(binary_markov_chain(*p), binary_markov_chain(*q)) + assert actual == pytest.approx(expected, abs=TOL) + + +@pytest.mark.parametrize("model", [bernoulli(0.3), golden_mean(0.4), even_process(0.6), binary_markov_chain(0.2, 0.7)]) +def test_self_divergence_is_zero(model): + assert relative_entropy_rate(model, model) == pytest.approx(0.0, abs=TOL) + + +def test_golden_mean_against_fair_coin(): + # Fair coin cross entropy is 1 bit; golden mean at p = 1/2 has h = 2/3. + assert relative_entropy_rate(golden_mean(0.5), bernoulli(0.5)) == pytest.approx(1 / 3, abs=TOL) + + +def test_infinite_when_support_not_included(): + assert math.isfinite(relative_entropy_rate(golden_mean(0.5), bernoulli())) + assert relative_entropy_rate(bernoulli(), golden_mean(0.5)) == math.inf + assert relative_entropy_rate(even_process(), golden_mean(0.5)) == math.inf + assert relative_entropy_rate(golden_mean(0.5), even_process()) == math.inf + bounds = relative_entropy_rate_bounds(bernoulli(), golden_mean(0.5), 3) + assert bounds.lower == bounds.upper == math.inf + + +def test_support_nfa_matches_positive_probability_words(): + for model in (golden_mean(0.4), even_process(0.3), _noisy_hmm()): + nfa = _support_nfa(model) + mealy = model.to_mealy() + pi = mealy.stationary_distribution() + for length in range(6): + for word in product(("0", "1"), repeat=length): + positive = model.word_probability(word, start=pi) > 1e-12 + assert nfa.recognizes(word) == positive + + +def test_support_includes_orders_golden_mean_below_full_shift(): + assert _support_includes(golden_mean(0.5), bernoulli()) + assert not _support_includes(bernoulli(), golden_mean(0.5)) + assert _support_includes(even_process(), bernoulli()) + + +def test_non_unifilar_reference_requires_bounds(): + with pytest.raises(NotImplementedError): + relative_entropy_rate(golden_mean(0.5), _noisy_hmm()) + + +@pytest.mark.parametrize( + ("p", "q"), [(golden_mean(0.4), binary_markov_chain(0.3, 0.6)), (even_process(), bernoulli(0.3))] +) +def test_bounds_bracket_exact_value_for_unifilar_reference(p, q): + exact = relative_entropy_rate(p, q) + for n in (1, 3, 6): + bounds = relative_entropy_rate_bounds(p, q, n) + assert bounds.lower - TOL <= exact <= bounds.upper + TOL + + +def test_bounds_converge_for_non_unifilar_reference(): + p, q = even_process(0.4), _noisy_hmm() + gaps = [] + for n in (2, 5, 9): + bounds = relative_entropy_rate_bounds(p, q, n) + assert bounds.lower <= bounds.upper + TOL + assert bounds.lower - TOL <= bounds.conditional_estimate + gaps.append(bounds.upper - bounds.lower) + assert gaps[0] > gaps[1] > gaps[2] + + +def test_self_divergence_bounds_contain_zero(): + q = _noisy_hmm() + bounds = relative_entropy_rate_bounds(q, q, 6) + assert bounds.lower == 0.0 + assert bounds.upper >= -TOL + + +def test_block_estimate_converges(): + p, q = golden_mean(0.4), binary_markov_chain(0.3, 0.6) + exact = relative_entropy_rate(p, q) + errors = [abs(relative_entropy_rate_bounds(p, q, n).block_estimate - exact) for n in (2, 4, 8)] + assert errors[0] > errors[1] > errors[2] + # A 1-step Markov reference makes the conditional estimate exact from n = 2 on. + assert relative_entropy_rate_bounds(p, q, 4).conditional_estimate == pytest.approx(exact, abs=1e-9) + + +def test_block_length_must_be_positive(): + with pytest.raises(ValueError, match="block_length"): + relative_entropy_rate_bounds(bernoulli(), bernoulli(), 0) + + +@settings(max_examples=60, deadline=None) +@given(epsilon_machines(max_states=3), epsilon_machines(max_states=3)) +def test_nonnegative_and_bracketed(p, q): + try: + exact = relative_entropy_rate(p, q) + except NotImplementedError: + assume(False) + assert exact >= 0.0 + if math.isfinite(exact): + bounds = relative_entropy_rate_bounds(p, q, 4) + assert bounds.lower - 1e-9 <= exact <= bounds.upper + 1e-9 + + +@settings(max_examples=40, deadline=None) +@given(mealy_hmms(max_states=2), mealy_hmms(max_states=2)) +def test_general_bounds_are_ordered(p, q): + bounds = relative_entropy_rate_bounds(p, q, 4) + assert 0.0 <= bounds.lower <= bounds.upper + 1e-9 diff --git a/tests/test_state_splitting.py b/tests/test_state_splitting.py new file mode 100644 index 0000000..8ab4902 --- /dev/null +++ b/tests/test_state_splitting.py @@ -0,0 +1,231 @@ +"""State splitting, amalgamation, and conjugacy invariants (Lind & Marcus §2.4, §7.4).""" + +from __future__ import annotations + +import math + +import numpy as np +import pytest +from hypothesis import HealthCheck, given, settings +from hypothesis import strategies as st + +from sofic.shifts import ( + TopologicalMarkovChain, + bowen_franks_group, + in_amalgamate, + in_edges, + in_split, + in_split_matrices, + jordan_form_away_from_zero, + out_amalgamate, + out_edges, + out_split, + out_split_matrices, +) +from sofic.shifts.algorithms import integer_adjacency_matrix + +SETTINGS = settings(max_examples=40, deadline=None, suppress_health_check=[HealthCheck.too_slow]) + + +def adjacency(tmc): + return np.array(integer_adjacency_matrix(tmc)[0], dtype=int) + + +def chain(rows): + return TopologicalMarkovChain.from_adjacency(np.array(rows, dtype=int)) + + +@st.composite +def split_cases(draw, *, outgoing=True): + """A small TMC plus a random partition of the out- (or in-) edges of some states.""" + n = draw(st.integers(1, 3)) + rows = draw(st.lists(st.lists(st.integers(0, 2), min_size=n, max_size=n), min_size=n, max_size=n)) + tmc = chain(rows) + edges_of = out_edges if outgoing else in_edges + partitions = {} + for state in tmc.states(): + edges = edges_of(tmc, state) + if not edges or not draw(st.booleans()): + continue + labels = draw(st.lists(st.integers(0, len(edges) - 1), min_size=len(edges), max_size=len(edges))) + parts = {} + for edge, label in zip(edges, labels, strict=True): + parts.setdefault(label, []).append(edge) + partitions[state] = list(parts.values()) + return tmc, partitions + + +def invariants(tmc, *, symbolic): + values = [tmc.periodic_points(n) for n in range(1, 7)] + values += [round(tmc.topological_entropy(), 9), bowen_franks_group(tmc)] + if symbolic: + values += [tmc.zeta_function(), jordan_form_away_from_zero(tmc)] + return values + + +# --- Known values ---------------------------------------------------------- + + +def example_2_4_6(): + """Lind & Marcus Figure 2.4.5(a): a, b: I -> J; c: I -> K; d: J -> I; e: K -> J; f: K -> I.""" + tmc = chain([[0, 2, 1], [1, 0, 0], [1, 1, 0]]) + i, j, k = tmc.states() + a, b = [e for e in out_edges(tmc, i) if e[1] == j] + (c,) = [e for e in out_edges(tmc, i) if e[1] == k] + e, f = sorted(out_edges(tmc, k), key=lambda edge: edge[1] != j) + return tmc, {i: [[a], [b, c]], k: [[e], [f]]} + + +def test_division_and_edge_matrices_example_2_4_6(): + tmc, partitions = example_2_4_6() + d, e = out_split_matrices(tmc, partitions) + assert d.tolist() == [[1, 1, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 1, 1]] + assert e.tolist() == [[0, 1, 0], [0, 1, 1], [1, 0, 0], [0, 1, 0], [1, 0, 0]] + split = out_split(tmc, partitions) + assert (d @ e == adjacency(tmc)).all() + assert (e @ d == adjacency(split)).all() + assert adjacency(split).tolist() == [ + [0, 0, 1, 0, 0], + [0, 0, 1, 1, 1], + [1, 1, 0, 0, 0], + [0, 0, 1, 0, 0], + [1, 1, 0, 0, 0], + ] + + +def test_complete_out_splitting_of_full_two_shift_example_2_4_5(): + tmc = chain([[2]]) + (state,) = tmc.states() + split = out_split(tmc, {state: [[edge] for edge in out_edges(tmc, state)]}) + assert adjacency(split).tolist() == [[1, 1], [1, 1]] + + +def test_split_copies_symbols(): + tmc = TopologicalMarkovChain(symbol_alphabet=frozenset("ab")) + tmc.add_transition("I", "I", "a") + tmc.add_transition("I", "I", "b") + split = out_split(tmc, {"I": [[e] for e in out_edges(tmc, "I")]}) + symbols = sorted((t.source, t.target, t.data["symbol"]) for t in split.transitions()) + assert symbols == sorted((("I", s), ("I", d), sym) for s, sym in enumerate("ab") for d in range(2)) + + +@pytest.mark.parametrize( + ("rows", "factors"), + [ + ([[2]], ()), + ([[1, 1], [1, 0]], ()), + ([[4, 1], [1, 0]], (4,)), + ([[3, 2], [2, 1]], (2, 2)), + ([[5, 3], [3, 5]], (7,)), + ([[6, 8], [1, 4]], (7,)), + ([[1]], (0,)), + ([[1, 0], [0, 1]], (0, 0)), + ], +) +def test_bowen_franks_known_values(rows, factors): + # Lind & Marcus Example 7.4.16 and Exercise 7.4.12. + assert bowen_franks_group(np.array(rows)).invariant_factors == factors + + +def test_bowen_franks_sign_and_inputs(): + assert bowen_franks_group(np.array([[2]])) == ((), -1) + assert bowen_franks_group(np.array([[1]])).det_sign == 0 + assert bowen_franks_group(np.array([[0, 1], [1, 0]])) == ((0,), 0) + assert bowen_franks_group(chain([[4, 1], [1, 0]])) == bowen_franks_group(np.array([[4, 1], [1, 0]])) + with pytest.raises(ValueError): + bowen_franks_group(np.array([[0.5]])) + + +def test_jordan_form_away_from_zero_known_values(): + sp = pytest.importorskip("sympy") + # Lind & Marcus Examples 7.4.7 and 7.4.8. + assert jordan_form_away_from_zero(np.array([[1, 1], [1, 1]])) == ((2, 1),) + assert jordan_form_away_from_zero(np.array([[3, 1, 1], [2, 2, 1], [1, 2, 2]])) == ((1, 2), (5, 1)) + assert jordan_form_away_from_zero(np.array([[3, 1, 1], [2, 2, 1], [2, 1, 2]])) == ((1, 1), (1, 1), (5, 1)) + assert jordan_form_away_from_zero(np.array([[0, 1], [0, 0]])) == () + golden = jordan_form_away_from_zero(np.array([[1, 1], [1, 0]])) + assert {root for root, _ in golden} == {(1 + sp.sqrt(5)) / 2, (1 - sp.sqrt(5)) / 2} + + +def test_invalid_partitions_and_amalgamations(): + tmc = chain([[1, 1], [1, 0]]) + i, j = tmc.states() + edges = out_edges(tmc, i) + with pytest.raises(ValueError, match="cover"): + out_split(tmc, {i: [edges[:1]]}) + with pytest.raises(ValueError, match="empty"): + out_split(tmc, {i: [edges, []]}) + with pytest.raises(ValueError, match="disjoint"): + out_split(tmc, {i: [edges, edges[:1]]}) + with pytest.raises(ValueError, match="columns"): + out_amalgamate(tmc, {"merged": [i, j]}) + with pytest.raises(ValueError, match="rows"): + in_amalgamate(tmc, {"merged": [i, j]}) + + +# --- Properties ------------------------------------------------------------- + + +@SETTINGS +@given(split_cases(outgoing=True)) +def test_out_split_matrices_and_round_trip(case): + tmc, partitions = case + split = out_split(tmc, partitions) + d, e = out_split_matrices(tmc, partitions) + assert ((d == 0) | (d == 1)).all() and (d.sum(axis=0) == 1).all() + assert (d @ e == adjacency(tmc)).all() + assert (e @ d == adjacency(split)).all() + classes = {state: [(state, i) for i in range(len(parts))] for state, parts in partitions.items()} + assert (adjacency(out_amalgamate(split, classes)) == adjacency(tmc)).all() + + +@SETTINGS +@given(split_cases(outgoing=False)) +def test_in_split_matrices_and_round_trip(case): + tmc, partitions = case + split = in_split(tmc, partitions) + d, e = in_split_matrices(tmc, partitions) + assert (d @ e == adjacency(tmc)).all() + assert (e @ d == adjacency(split)).all() + classes = {state: [(state, i) for i in range(len(parts))] for state, parts in partitions.items()} + assert (adjacency(in_amalgamate(split, classes)) == adjacency(tmc)).all() + + +@SETTINGS +@given(split_cases(outgoing=True)) +def test_out_splitting_preserves_conjugacy_invariants(case): + pytest.importorskip("sympy") + tmc, partitions = case + assert invariants(out_split(tmc, partitions), symbolic=True) == invariants(tmc, symbolic=True) + + +@SETTINGS +@given(split_cases(outgoing=False)) +def test_in_splitting_preserves_conjugacy_invariants(case): + pytest.importorskip("sympy") + tmc, partitions = case + assert invariants(in_split(tmc, partitions), symbolic=True) == invariants(tmc, symbolic=True) + + +@SETTINGS +@given( + st.integers(1, 4).flatmap( + lambda n: st.lists(st.lists(st.integers(-3, 3), min_size=n, max_size=n), min_size=n, max_size=n) + ) +) +def test_bowen_franks_matches_sympy_smith_form(rows): + sp = pytest.importorskip("sympy") + from sympy.matrices.normalforms import smith_normal_form + + difference = sp.eye(len(rows)) - sp.Matrix(rows) + diagonal = [ + abs(int(difference_entry)) for difference_entry in smith_normal_form(difference, domain=sp.ZZ).diagonal() + ] + group = bowen_franks_group(np.array(rows)) + assert sorted(group.invariant_factors) == sorted(d for d in diagonal if d != 1) + nonzero = [d for d in group.invariant_factors if d] + assert all(b % a == 0 for a, b in zip(nonzero, nonzero[1:], strict=False)) + det = int(difference.det()) + assert group.det_sign == (det > 0) - (det < 0) + if det: + assert math.prod(group.invariant_factors) == abs(det) diff --git a/tests/test_zeta.py b/tests/test_zeta.py new file mode 100644 index 0000000..57fa87a --- /dev/null +++ b/tests/test_zeta.py @@ -0,0 +1,230 @@ +"""Periodic points and zeta functions (Lind & Marcus §6.4).""" + +from __future__ import annotations + +import math +from itertools import product + +import numpy as np +import pytest +from hypothesis import HealthCheck, assume, given, settings +from hypothesis import strategies as st + +from sofic.graph import ATTR_SYMBOL, TransitionGraph +from sofic.shifts import ShiftOfFiniteType, SoficShift, TopologicalMarkovChain, higher_block_presentation +from sofic.shifts.algorithms import reciprocal_characteristic_polynomial, signed_subset_matrices +from sofic.testing.strategies import sfts, sofic_shifts + +sp = pytest.importorskip("sympy") +t = sp.Symbol("t") + +SETTINGS = settings(max_examples=40, deadline=None, suppress_health_check=[HealthCheck.too_slow]) + + +def brute_force_periodic_points(shift, n): + """Count words ``w`` of length ``n`` with ``w^infinity`` in the shift. + + ``w^infinity`` is a bi-infinite path label iff some state returns to itself + reading ``w^k`` for some ``1 <= k <= |V|`` (pigeonhole on the states at + multiples of ``n``), which is weaker than a cycle labeled ``w`` itself. + """ + outgoing = {state: [] for state in shift.states()} + for transition in shift.transitions(): + outgoing[transition.source].append((transition.data[ATTR_SYMBOL], transition.target)) + states = list(outgoing) + + def read(current, word): + for symbol in word: + current = {target for state in current for label, target in outgoing[state] if label == symbol} + return current + + count = 0 + for word in product(sorted(shift.symbol_alphabet), repeat=n): + if any(_returns_within(read, state, word, len(states)) for state in states): + count += 1 + return count + + +def _returns_within(read, state, word, max_repeats): + """Whether reading ``word^k`` from ``state`` can return to it for some ``1 <= k <= max_repeats``.""" + current = {state} + for _ in range(max_repeats): + current = read(current, word) + if state in current: + return True + return False + + +def series_from_zeta(zeta, order): + """``p_n = n [t^n] log zeta`` for ``n = 1..order``.""" + series = sp.series(sp.log(zeta), t, 0, order + 1).removeO() + return [int(n * series.coeff(t, n)) for n in range(1, order + 1)] + + +def block_recoding(shift, order): + """The higher block shift ``X^[order+1]``, conjugate to ``shift``.""" + refined = higher_block_presentation(shift, order) + graph = TransitionGraph() + for state in refined.states(): + graph.add_state(state) + labels = set() + for transition in refined.transitions(): + label = transition.source[0] + (transition.data[ATTR_SYMBOL],) + labels.add(label) + graph.add_transition(transition.source, transition.target, **{ATTR_SYMBOL: label}) + return SoficShift(graph=graph, symbol_alphabet=frozenset(labels)) + + +def as_sofic(shift): + return SoficShift(graph=shift.graph.copy(), symbol_alphabet=shift.symbol_alphabet) + + +def even_shift(): + shift = SoficShift(symbol_alphabet=frozenset("ab")) + shift.add_transition(0, 0, "a") + shift.add_transition(0, 1, "b") + shift.add_transition(1, 0, "b") + return shift + + +def same_rational(f, g): + return sp.cancel(f - g) == 0 + + +tmc_matrices = st.integers(1, 3).flatmap( + lambda n: st.lists(st.lists(st.integers(0, 2), min_size=n, max_size=n), min_size=n, max_size=n) +) + + +# --- Known values ---------------------------------------------------------- + + +def test_full_two_shift(): + shift = TopologicalMarkovChain.from_adjacency(np.array([[2]])) + assert [shift.periodic_points(n) for n in range(1, 8)] == [2**n for n in range(1, 8)] + assert same_rational(shift.zeta_function(t), 1 / (1 - 2 * t)) + + +def test_golden_mean_lucas_numbers(): + shift = TopologicalMarkovChain.from_adjacency(np.array([[1, 1], [1, 0]])) + assert [shift.periodic_points(n) for n in range(1, 8)] == [1, 3, 4, 7, 11, 18, 29] + assert same_rational(shift.zeta_function(t), 1 / (1 - t - t**2)) + + +def test_even_shift_example_6_4_10(): + shift = even_shift() + lucas = [1, 3, 4, 7, 11, 18, 29] + assert [shift.periodic_points(n) for n in range(1, 8)] == [lucas[n - 1] - (-1) ** n for n in range(1, 8)] + assert same_rational(shift.zeta_function(t), (1 + t) / (1 - t - t**2)) + matrices = signed_subset_matrices(shift) + assert [(m.tolist(), sign) for m, sign in matrices] == [([[1, 1], [1, 0]], 1), ([[-1]], -1)] + + +def test_golden_mean_sft_from_forbidden_words(): + shift = ShiftOfFiniteType.from_forbidden_words({(1, 1)}, frozenset({0, 1})) + assert [shift.periodic_points(n) for n in range(1, 8)] == [1, 3, 4, 7, 11, 18, 29] + assert same_rational(shift.zeta_function(), 1 / (1 - t - t**2)) + + +@pytest.mark.parametrize( + "rows", + [ + [["", "a", "b"], ["a", "", "b"], ["a", "b", ""]], + [["a", "b", ""], ["", "", "b"], ["b", "", "a"]], + ], +) +def test_exercise_6_4_2_matches_brute_force(rows): + shift = SoficShift(symbol_alphabet=frozenset("ab")) + for i, row in enumerate(rows): + for j, label in enumerate(row): + if label: + shift.add_transition(i, j, label) + zeta = shift.zeta_function(t) + expected = [brute_force_periodic_points(shift, n) for n in range(1, 8)] + assert [shift.periodic_points(n) for n in range(1, 8)] == expected + assert series_from_zeta(zeta, 7) == expected + + +def test_nondeterministic_cycle_counts_once(): + # 0^infinity is presented only by the 2-cycle 0 -> 1 -> 0 and the self-loop at 2. + shift = SoficShift(symbol_alphabet=frozenset("0")) + shift.add_transition(0, 1, "0") + shift.add_transition(1, 0, "0") + shift.add_transition(2, 2, "0") + assert [shift.periodic_points(n) for n in range(1, 5)] == [1, 1, 1, 1] + assert same_rational(shift.zeta_function(t), 1 / (1 - t)) + + +def test_empty_shift_and_bad_period(): + shift = SoficShift(symbol_alphabet=frozenset("a")) + shift.add_transition(0, 1, "a") + assert shift.periodic_points(3) == 0 + assert shift.zeta_function(t) == 1 + with pytest.raises(ValueError): + shift.periodic_points(0) + empty = ShiftOfFiniteType.from_forbidden_words({("a",)}, frozenset("a")) + assert empty.periodic_points(2) == 0 + + +def test_exact_for_large_periods(): + shift = TopologicalMarkovChain.from_adjacency(np.array([[3]])) + assert shift.periodic_points(90) == 3**90 + + +def test_reciprocal_characteristic_polynomial(): + a = np.array([[3, 1, 1], [2, 2, 1], [1, 2, 2]], dtype=object) + # chi_A = (t - 5)(t - 1)^2 (Lind & Marcus Example 7.4.8) + assert reciprocal_characteristic_polynomial(a) == [1, -7, 11, -5] + + +# --- Properties ------------------------------------------------------------- + + +@SETTINGS +@given(sofic_shifts(max_states=3)) +def test_sofic_matches_brute_force(shift): + for n in range(1, 9): + assert shift.periodic_points(n) == brute_force_periodic_points(shift, n) + + +@SETTINGS +@given(sofic_shifts(max_states=3)) +def test_zeta_series_is_periodic_points(shift): + zeta = shift.zeta_function(t) + assert series_from_zeta(zeta, 6) == [shift.periodic_points(n) for n in range(1, 7)] + + +@SETTINGS +@given(sofic_shifts(max_states=3), st.integers(1, 2)) +def test_higher_block_recoding_is_invariant(shift, order): + recoded = block_recoding(shift, order) + assert [recoded.periodic_points(n) for n in range(1, 7)] == [shift.periodic_points(n) for n in range(1, 7)] + assert same_rational(recoded.zeta_function(t), shift.zeta_function(t)) + + +@SETTINGS +@given(sfts(max_states=12)) +def test_sft_block_matrix_agrees_with_manning(sft): + presentation = as_sofic(sft) + for n in range(1, 7): + assert sft.periodic_points(n) == presentation.periodic_points(n) + assert same_rational(sft.zeta_function(t), presentation.zeta_function(t)) + + +@SETTINGS +@given(tmc_matrices) +def test_tmc_trace_and_zeta(rows): + matrix = np.array(rows, dtype=int) + shift = TopologicalMarkovChain.from_adjacency(matrix) + for n in range(1, 6): + assert shift.periodic_points(n) == int(np.trace(np.linalg.matrix_power(matrix, n))) + assert same_rational(shift.zeta_function(t), 1 / sp.Matrix(sp.eye(len(rows)) - t * sp.Matrix(rows)).det()) + + +@SETTINGS +@given(sofic_shifts(max_states=3)) +def test_smallest_pole_gives_entropy(shift): + _, denominator = sp.fraction(shift.zeta_function(t)) + assume(sp.degree(denominator, t) > 0) + radius = min(abs(complex(root)) for root in sp.Poly(denominator, t).nroots()) + assert math.isclose(math.log2(1 / radius), shift.topological_entropy(), abs_tol=1e-6)