Literature gaps: zeta functions, conjugacy invariants, decisions, regex, entropy rates, correlations, predictive RD (3/4)

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
 

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 <sofic.automata.nfa.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

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

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

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)
 ========================================
 

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

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

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