From 842f936c429ddacf161ffe081e050d195308b901 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Sun, 25 Oct 2026 15:29:21 -0600 Subject: [PATCH 1/2] Add sofic.algebra: Krohn-Rhodes theory for deterministic presentations Transformation semigroups, Green's relations, syntactic monoids and star-freeness, the holonomy decomposition as a verifiable cascade product, Krohn-Rhodes complexity bounds and prime factors, and hooks on automata, unifilar HMMs / epsilon-machines, and sofic shifts. Following DeDeo (2011): effective theories from truncated cascades, random-map representations of probabilistic automata, walks over semigroup elements, and an irreversibility / dissipation information split. Also language varieties, minimal ideals and synchronization, Schutzenberger groups, and belief-map semigroups of HMMs. Co-authored-by: Cursor --- CHANGELOG.md | 26 ++ docs/algebra/algebra.rst | 39 +++ docs/algebra/complexity.rst | 123 +++++++ docs/algebra/effective_theories.rst | 59 ++++ docs/algebra/holonomy.rst | 77 +++++ docs/algebra/noise.rst | 118 +++++++ docs/algebra/semigroups.rst | 169 +++++++++ docs/examples.rst | 29 ++ docs/index.rst | 8 +- docs/references.bib | 192 +++++++++++ sofic/algebra/__init__.py | 104 ++++++ sofic/algebra/_groups.py | 109 ++++++ sofic/algebra/adapters.py | 108 ++++++ sofic/algebra/belief.py | 99 ++++++ sofic/algebra/cascade.py | 209 ++++++++++++ sofic/algebra/complexity.py | 190 +++++++++++ sofic/algebra/green.py | 207 +++++++++++ sofic/algebra/holonomy.py | 512 ++++++++++++++++++++++++++++ sofic/algebra/information.py | 75 ++++ sofic/algebra/random_maps.py | 390 +++++++++++++++++++++ sofic/algebra/syntactic.py | 148 ++++++++ sofic/algebra/transformation.py | 511 +++++++++++++++++++++++++++ sofic/algebra/varieties.py | 141 ++++++++ sofic/automata/base.py | 22 ++ sofic/automata/wheeler.py | 27 ++ sofic/exceptions.py | 9 + sofic/generators/mealy.py | 22 ++ sofic/generators/wheeler_epsilon.py | 5 +- sofic/shifts/__init__.py | 8 +- sofic/shifts/invariants.py | 29 ++ sofic/shifts/sofic.py | 21 ++ sofic/shifts/wheeler.py | 13 +- tests/test_algebra_belief.py | 38 +++ tests/test_algebra_dedeo.py | 189 ++++++++++ tests/test_algebra_holonomy.py | 280 +++++++++++++++ tests/test_algebra_random_maps.py | 186 ++++++++++ tests/test_algebra_semigroups.py | 246 +++++++++++++ tests/test_algebra_varieties.py | 121 +++++++ 38 files changed, 4848 insertions(+), 11 deletions(-) create mode 100644 docs/algebra/algebra.rst create mode 100644 docs/algebra/complexity.rst create mode 100644 docs/algebra/effective_theories.rst create mode 100644 docs/algebra/holonomy.rst create mode 100644 docs/algebra/noise.rst create mode 100644 docs/algebra/semigroups.rst create mode 100644 sofic/algebra/__init__.py create mode 100644 sofic/algebra/_groups.py create mode 100644 sofic/algebra/adapters.py create mode 100644 sofic/algebra/belief.py create mode 100644 sofic/algebra/cascade.py create mode 100644 sofic/algebra/complexity.py create mode 100644 sofic/algebra/green.py create mode 100644 sofic/algebra/holonomy.py create mode 100644 sofic/algebra/information.py create mode 100644 sofic/algebra/random_maps.py create mode 100644 sofic/algebra/syntactic.py create mode 100644 sofic/algebra/transformation.py create mode 100644 sofic/algebra/varieties.py create mode 100644 tests/test_algebra_belief.py create mode 100644 tests/test_algebra_dedeo.py create mode 100644 tests/test_algebra_holonomy.py create mode 100644 tests/test_algebra_random_maps.py create mode 100644 tests/test_algebra_semigroups.py create mode 100644 tests/test_algebra_varieties.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 48e6012..652a406 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -135,6 +135,32 @@ brute-force reference or a closed-form value and has a regression test. ### New +- `sofic.algebra`: algebraic automata theory for deterministic presentations. + `TransformationSemigroup` (lazy closure, shortest words, sink completion of + partial maps), Green's relations and maximal subgroups, `syntactic_monoid` + and `is_star_free`, and the Krohn-Rhodes holonomy decomposition + (`holonomy_decomposition`) with an explicit `CascadeProduct` whose + `emulates()` checks the covering. `krohn_rhodes_complexity` returns + `ComplexityBounds(lower, upper, exact)`, and `prime_decomposition` lists the + simple-group factors. `holonomy_information` splits `C_mu` across the + holonomy levels. DFAs, unifilar HMMs / ε-machines and sofic shifts gain + `transition_semigroup()` and `krohn_rhodes_complexity()`, and the shift + invariants gain `factor_syntactic_monoid` and + `factor_krohn_rhodes_complexity`. Wheeler-presentation failures now report + whether the language is star-free. +- `sofic.algebra`, following DeDeo's "Effective Theories for Circuits and + Automata" (2011): `HolonomyDecomposition.effective_theory(k)` and + `CascadeProduct.truncate(k)` give coarse-grained machines on overlapping + superstates; `independent_random_maps` / `coupled_random_maps` / + `random_maps` rewrite probabilistic automata as random choices of + deterministic maps; `element_walk` and `classify_elements` follow the walk + over semigroup elements; and `walk_information` splits lost information into + irreversibility and dissipation. Also new: language varieties + (`is_piecewise_testable`, `is_locally_testable`, `is_definite`, ...), + `minimal_ideal` / `minimal_rank` / `is_synchronizing` / + `schutzenberger_group`, belief-map semigroups (`projective_semigroup`, + `mixed_state_semigroup`), and `is_full_transformation_monoid`, which makes + the complexity of `T_n` exact without enumerating it. - `sofic.generators.capacity`: capacity bounds for finite ε-transducers (unifilar finite-state channels). `EpsilonTransducer.feedback_capacity()` brackets the feedback capacity with Q-graph bounds (convex upper bound, diff --git a/docs/algebra/algebra.rst b/docs/algebra/algebra.rst new file mode 100644 index 0000000..c68f2eb --- /dev/null +++ b/docs/algebra/algebra.rst @@ -0,0 +1,39 @@ +.. algebra.rst +.. py:module:: sofic.algebra + +******* +Algebra +******* + +The :mod:`sofic.algebra` package brings algebraic automata theory to every +deterministic presentation in ``sofic``. A DFA, a unifilar automaton, a +unifilar HMM or ε-machine, or a right-resolving shift presentation has one +partial map of its states per symbol. Those maps generate a finite +*transformation semigroup*, and the Krohn-Rhodes theorem breaks that +semigroup into a cascade of simple groups and flip-flops +:cite:`KrohnRhodes1965` :cite:`RhodesNehanivHirsch2010`. + +* :doc:`semigroups` covers transformation semigroups, Green's relations, + syntactic monoids, and star-freeness. +* :doc:`holonomy` covers the holonomy decomposition as an explicit cascade + product that can be checked against the original machine. +* :doc:`complexity` covers Krohn-Rhodes complexity bounds, prime factors, + how statistical complexity splits across the holonomy levels, and the + semigroup of an HMM's belief updates. +* :doc:`effective_theories` covers truncating the cascade into coarse-grained + machines on overlapping superstates :cite:`DeDeo2011`. +* :doc:`noise` covers probabilistic automata: rewriting them as random choices + of deterministic maps, the walk over semigroup elements, and how much + information is erased versus dissipated. + +The group-theoretic outputs need the optional ``sympy`` dependency, which the +``sofic[symbolic]`` extra installs. + +.. toctree:: + :maxdepth: 1 + + semigroups + holonomy + complexity + effective_theories + noise diff --git a/docs/algebra/complexity.rst b/docs/algebra/complexity.rst new file mode 100644 index 0000000..eda0549 --- /dev/null +++ b/docs/algebra/complexity.rst @@ -0,0 +1,123 @@ +.. complexity.rst +.. py:module:: sofic.algebra.complexity + +*********************** +Krohn-Rhodes Complexity +*********************** + +The *group complexity* of a finite semigroup is the least number of group +levels, separated by aperiodic levels, in a cascade that emulates it +:cite:`KrohnRhodes1965` :cite:`RhodesSteinberg2009`. It is zero exactly for +aperiodic semigroups, which for syntactic monoids means star-free languages. +It equals ``n - 1`` for the full transformation monoid ``T_n``. Whether it is +computable in general is open; only complexity one is known to be decidable +:cite:`HenckellRhodesSteinberg2010`. + +:func:`krohn_rhodes_complexity` therefore returns +:class:`ComplexityBounds`: + +* The upper bound counts the holonomy levels with a nontrivial group. +* The lower bound is ``0`` for an aperiodic semigroup and at least ``1`` + otherwise. It is raised to ``k - 1`` when the elements that map some + ``k``-element image set into itself act on it as all of ``T_k``, because + complexity cannot increase under division. + +``exact`` is filled in when the two bounds meet. That happens for every +aperiodic semigroup, for every semigroup with a single nontrivial holonomy +level, and for ``T_n``. The sharper Rhodes-Tilson lower bounds are not +implemented. + +.. ipython:: + + In [1]: from sofic.algebra import TransformationSemigroup, krohn_rhodes_complexity, prime_decomposition + + In [2]: t4 = TransformationSemigroup({"t": (1, 0, 2, 3), "c": (1, 2, 3, 0), "e": (0, 0, 2, 3)}) + + @doctest + In [3]: krohn_rhodes_complexity(t4) + Out[3]: ComplexityBounds(lower=3, upper=3, exact=3) + + @doctest + In [4]: prime_decomposition(t4).simple_group_orders() + Out[4]: (2, 3) + +:func:`prime_decomposition` lists the Jordan-Hölder factors of every +holonomy group. These are the simple groups the Krohn-Rhodes theorem +requires, and every simple group divisor of ``S`` divides one of them. It +also reports whether flip-flops (resets) are needed. A nonabelian simple +factor is reported by its order only, because order alone does not identify +it. + +.. autofunction:: krohn_rhodes_complexity +.. autoclass:: ComplexityBounds +.. autofunction:: prime_decomposition +.. autoclass:: PrimeDecomposition + :members: +.. autoclass:: PrimeFactor + +Holonomy levels and statistical complexity +========================================== + +.. py:module:: sofic.algebra.information + +Lifting each causal state to its holonomy coordinates turns the stationary +state distribution into a joint distribution over the level variables +``L_1, ..., L_h``. Distinct states get distinct coordinates, so +``H[L_1, ..., L_h] = C_mu``. The chain rule then splits the statistical +complexity into ``H[L_i | L_ eta T^(x) / (eta T^(x) 1)``. +:func:`projective_semigroup` enumerates the word matrices ``T^(w)`` up to a +positive scalar, up to a cap, and returns the result as a transformation +monoid, so the rest of :mod:`sofic.algebra` applies. If the enumeration +closes, every initial belief has a finite orbit. Finiteness of matrix +semigroups is decidable :cite:`MandelSimon1977` :cite:`Jacob1977`, but this +is only a capped enumeration. It is also only a sufficient condition: the +mixed-state presentation from one particular belief can be finite even when +this semigroup is infinite. :func:`mixed_state_semigroup` takes the +transition semigroup of the mixed-state presentation directly. + +.. ipython:: + + In [1]: from sofic.algebra import projective_semigroup + + @doctest + In [2]: len(projective_semigroup(even_process()).states) + Out[2]: 7 + +For a unifilar machine like the even process this matches the transition +monoid. The simple nonunifilar source, whose mixed-state presentation is +infinite, never closes and raises +:class:`~sofic.exceptions.SemigroupTooLargeError`. + +.. autofunction:: projective_semigroup +.. autofunction:: mixed_state_semigroup diff --git a/docs/algebra/effective_theories.rst b/docs/algebra/effective_theories.rst new file mode 100644 index 0000000..c130fa8 --- /dev/null +++ b/docs/algebra/effective_theories.rst @@ -0,0 +1,59 @@ +.. effective_theories.rst + +****************** +Effective Theories +****************** + +Each level of a holonomy cascade reads only the levels above it. So keeping +the top ``k`` levels gives a well-defined machine of its own, a +coarse-grained "effective theory" that tracks the large-scale logic of the +computation while ignoring the details below :cite:`DeDeo2011`. Its states +are *superstates*: sets of original states that agree on the top ``k`` +coordinates of some lift. Superstates can overlap. + +:meth:`~sofic.algebra.holonomy.HolonomyDecomposition.effective_theory` +returns the truncated machine as a DFA, together with its superstates and +the projection from full coordinates. The projection commutes with every +input. :meth:`~sofic.algebra.cascade.CascadeProduct.truncate` does the same +for a cascade. + +The example is Fig. 4 of :cite:`DeDeo2011`: an eight-state machine with two +inputs, built from XOR gates on a three-bit shift register, here transcribed +from the paper's TeX source. + +.. ipython:: + + In [1]: from sofic.algebra import TransformationSemigroup, holonomy_decomposition + + In [2]: a = {"A": "B", "B": "E", "C": "F", "D": "A", "E": "H", "F": "C", "G": "D", "H": "G"} + + In [3]: b = {"A": "B", "B": "C", "C": "F", "D": "B", "E": "G", "F": "C", "G": "F", "H": "G"} + + In [4]: hd = holonomy_decomposition(TransformationSemigroup.from_maps({"a": a, "b": b}, states="ABCDEFGH")) + + @doctest + In [5]: [level.order for level in hd.levels] + Out[5]: [3, 2, 2] + + In [6]: top = hd.effective_theory(1) + + @doctest + In [7]: sorted("".join(sorted(s)) for s in top.superstates.values()) + Out[7]: ['ACFH', 'BCFG', 'CDEF'] + +The top level is a counter modulo three on ``a`` and a reset on ``b``. Its +three superstates are AH, BG, and DE, and each one also contains C and F, +exactly as in the paper's Fig. 4b. The two lower levels are cyclic groups of +order two. + +Given a stochastic machine, ``effective_theory`` also returns the stationary +mass of each superstate, grouping states by their canonical lift. It reports +whether that grouping is strongly lumpable +(:func:`~sofic.generators.lumping.is_lumpable`), that is, whether the +coarse-grained process is itself Markov on the superstates. + +.. autoclass:: sofic.algebra.holonomy.EffectiveTheory +.. automethod:: sofic.algebra.holonomy.HolonomyDecomposition.effective_theory + :noindex: +.. automethod:: sofic.algebra.cascade.CascadeProduct.truncate + :noindex: diff --git a/docs/algebra/holonomy.rst b/docs/algebra/holonomy.rst new file mode 100644 index 0000000..22604f7 --- /dev/null +++ b/docs/algebra/holonomy.rst @@ -0,0 +1,77 @@ +.. holonomy.rst +.. py:module:: sofic.algebra.holonomy + +********************** +Holonomy Decomposition +********************** + +The holonomy method of Zeiger :cite:`Zeiger1967` and Eilenberg +:cite:`Eilenberg1976` builds a cascade product that emulates a +transformation semigroup. Each level of the cascade is a permutation group +(a *holonomy group*) together with constant maps (resets). The implementation +follows the constructive proof of Egri-Nagy & Nehaniv +:cite:`EgriNagyNehaniv2015`: + +1. The **extended image set** contains the images ``Q · s`` of the full state + set, plus all singletons. A set ``P`` approximates any state ``x ∈ P``. +2. Sets that are images of each other are **equivalent**, and the classes are + the strongly connected components of ``S`` acting on image sets. + **Subduction** (``P ⊆ R · s`` for some ``s``) orders the classes. The + **height** of a set is the length of the longest strict subduction chain of + non-singletons ending at it, and the cascade has one level per height, + numbered by depth from the top. +3. The **tiles** of a set are its maximal proper subsets in the image set. A + class representative's **holonomy group** is the permutation group its + setwise stabilizer induces on its tiles. It is generated by round trips + through the class, so it is computed from the generators alone, even when + ``S`` is too large to enumerate. + +A state is lifted to a maximal chain of tiles from ``Q`` down to the state's +singleton. Each tile in the chain, recorded relative to its parent's class +representative, becomes one cascade coordinate. An input acts on a chain +level by level, from the top. Each level either applies a holonomy +permutation or resets, and the choice depends only on the coordinates above +it. + +.. ipython:: + + In [1]: from sofic.algebra import TransformationSemigroup, holonomy_decomposition + + In [2]: sg = TransformationSemigroup({"s1": (1, 0, 2), "s2": (0, 1, 1)}) + + In [3]: hd = holonomy_decomposition(sg) + + @doctest + In [4]: [level.order for level in hd.levels] + Out[4]: [1, 2] + + In [5]: cascade = hd.to_cascade() + + @doctest + In [6]: cascade.emulates(sg) + Out[6]: True + +This is Example 7.1 of :cite:`EgriNagyNehaniv2015`. ``{0, 1}`` is swapped by +``s1``, which gives a cyclic group of order two at depth 2. The top level +only records whether the state lies in ``{0, 1}`` or ``{2}``. +:meth:`HolonomyDecomposition.describe` summarizes the levels. +:meth:`~sofic.algebra.cascade.CascadeProduct.emulates` checks that decoding +commutes with every generator, so each decomposition verifies itself. +:meth:`~sofic.algebra.cascade.CascadeProduct.as_dfa` returns the cascade as +an ordinary :class:`~sofic.automata.dfa.DFA` on coordinate tuples. + +.. autofunction:: holonomy_decomposition +.. autoclass:: HolonomyDecomposition + :members: +.. autoclass:: HolonomyLevel + :members: + +Cascade products +================ + +.. py:module:: sofic.algebra.cascade + +.. autoclass:: CascadeProduct + :members: +.. autoclass:: LevelAction + :members: diff --git a/docs/algebra/noise.rst b/docs/algebra/noise.rst new file mode 100644 index 0000000..c05c1f6 --- /dev/null +++ b/docs/algebra/noise.rst @@ -0,0 +1,118 @@ +.. noise.rst +.. py:module:: sofic.algebra.random_maps + +******************** +Noisy Automata +******************** + +A probabilistic automaton gives each input letter a row-stochastic matrix. +It can be rewritten as a random choice of a deterministic map of the states, +followed by applying that map :cite:`Maler1995`. The deterministic part is an +ordinary transformation semigroup, so it has a holonomy decomposition. How +the random choice is made is a *failure model* for the underlying mechanism +:cite:`DeDeo2011`: + +* :func:`independent_random_maps` lets every state fail independently: + ``pi(m) = prod_q P[q, m(q)]`` (Eq. 2 of :cite:`DeDeo2011`). +* :func:`coupled_random_maps` uses one shared random draw for every state, + the random mapping representation of :cite:`LevinPeresWilmer2009`. Failures + are then perfectly correlated. +* :func:`random_maps` takes any decomposition and checks it against the + matrices. + +The example is the two-bit XOR register of :cite:`DeDeo2011` with a 1% +bit-flip rate. Each state appends a 1 with the probability given in the list +``[eps, 1 - eps, 1 - eps, eps]``: + +.. ipython:: + + In [1]: import numpy as np + + In [2]: from sofic.algebra import independent_random_maps, coupled_random_maps + + In [3]: eps = 0.01 + + In [4]: P = np.array([[1 - eps, eps, 0, 0], [0, 0, eps, 1 - eps], [eps, 1 - eps, 0, 0], [0, 0, 1 - eps, eps]]) + + In [5]: independent, coupled = independent_random_maps(P), coupled_random_maps(P) + + @doctest + In [6]: len(independent.maps("t")), len(coupled.maps("t")) + Out[6]: (16, 2) + + @doctest + In [7]: len(independent.semigroup()), len(coupled.semigroup()) + Out[7]: (256, 12) + +Independent failures give all sixteen maps compatible with the register, +and together they generate the full transformation monoid on four states. +Correlated failures give just two maps, XOR and NOT-XOR, which generate the +alternating group ``A_4``. + +Walks on the semigroup +====================== + +:func:`element_walk` returns the distribution of the composed random map +``M_t = m_1 ... m_t`` over semigroup elements at each time ``t``, a Markov +chain on the right Cayley graph. :func:`classify_elements` labels the +identity, the permutations, the constant maps (resets), and the minimal +ideal. Together they reproduce Figs. 6 and 8 of :cite:`DeDeo2011`. Under +independent noise the mass drains into the resets, which form the minimal +ideal, after a brief window in which other permutations outweigh them. Under +correlated noise no reset ever appears, and the walk spreads over ``A_4``. +Here it cycles with period three through the cosets of the Klein group, so +it is uniform on ``A_4`` only on average. + +.. ipython:: + + In [1]: from sofic.algebra import classify_elements, element_walk + + In [2]: sg = independent.semigroup(); resets = sorted(classify_elements(sg)["resets"]) + + @doctest + In [3]: round(float(element_walk(independent, 1000)[-1, resets].sum()), 3) + Out[3]: 0.998 + +Erased versus dissipated information +==================================== + +DeDeo separates *irreversibility*, where a map deletes information, from +*dissipation*, where different histories become indistinguishable because +nobody knows which random operation ran. :func:`walk_information` measures +both. Take ``X_0`` uniform and independent of the noise, and let +``X_t = X_0 · M_t``. Then: + +* ``irreversibility = H[X_0] - I[X_0 : X_t | M_t]`` is the information gone + even if the random choices are known. +* ``dissipation = I[X_0 : X_t | M_t] - I[X_0 : X_t]`` is the information lost + only because the random choices are unknown. + +.. note:: Original, uncited: this information-theoretic split is not in + :cite:`DeDeo2011`. + +.. ipython:: + + In [1]: from sofic.algebra import walk_information + + In [2]: ind, cpl = walk_information(independent, 1000), walk_information(coupled, 5000) + + @doctest + In [3]: round(float(ind.irreversibility[-1]), 2), round(float(cpl.irreversibility[-1]), 2) + Out[3]: (2.0, 0.0) + + @doctest + In [4]: round(float(cpl.dissipation[-1]), 2) + Out[4]: 2.0 + +Independent noise erases the two bits of the initial state. Correlated noise +erases nothing, but loses both bits to ignorance of the noise. + +.. autoclass:: RandomMapRepresentation + :members: +.. autofunction:: independent_random_maps +.. autofunction:: coupled_random_maps +.. autofunction:: random_maps +.. autofunction:: element_walk +.. autofunction:: classify_elements +.. autofunction:: walk_information +.. autoclass:: WalkInformation diff --git a/docs/algebra/semigroups.rst b/docs/algebra/semigroups.rst new file mode 100644 index 0000000..08213a2 --- /dev/null +++ b/docs/algebra/semigroups.rst @@ -0,0 +1,169 @@ +.. semigroups.rst +.. py:module:: sofic.algebra.transformation + +************************* +Transformation Semigroups +************************* + +A :class:`TransformationSemigroup` is generated by labeled total maps of a +finite set, acting on the right so that a word is read left to right +:cite:`Eilenberg1976`. Its elements are enumerated lazily by breadth-first +closure, and each one keeps a shortest generating word. A ``max_size`` cap +raises :class:`~sofic.exceptions.SemigroupTooLargeError`, because ``n`` +states can generate up to ``n**n`` maps. Partial maps, such as the undefined +moves of a unifilar presentation, are completed with an absorbing ``SINK`` +state. + +.. ipython:: + + In [1]: from sofic.algebra import TransformationSemigroup + + In [2]: t3 = TransformationSemigroup({"t": (1, 0, 2), "c": (1, 2, 0), "e": (0, 0, 2)}) + + @doctest + In [3]: len(t3), t3.is_aperiodic() + Out[3]: (27, False) + + @doctest + In [4]: t3.word(t3.evaluate("ec")) + Out[4]: ('e', 'c') + +These three generators give the full transformation monoid ``T_3``. Any sofic +model with symbol-labeled edges converts with +:func:`~sofic.algebra.adapters.transformation_semigroup`, or with the +``transition_semigroup()`` method on automata, ε-machines, and sofic shifts: + +.. ipython:: + + In [1]: from sofic.examples import even_process + + In [2]: sg = even_process().transition_semigroup() + + @doctest + In [3]: sg.states, len(sg) + Out[3]: (('A', 'B', SINK), 7) + +.. autoclass:: Transformation + :members: + +.. autoclass:: TransformationSemigroup + :members: + +.. autodata:: SINK + +.. autofunction:: sofic.algebra.adapters.transformation_semigroup +.. autofunction:: sofic.algebra.adapters.semiautomaton + +Green's relations +================= + +:func:`~sofic.algebra.green.green_structure` computes the R-, L-, H-, and +D-classes (D = J for finite semigroups) as strongly connected components of +the Cayley graphs :cite:`Green1951` :cite:`Pin1986`. Each D-class is laid out +as an egg box of H-classes, and the H-class of an idempotent is a maximal +subgroup, returned as a sympy permutation group on the idempotent's image. + +.. ipython:: + + In [1]: green = t3.green() + + @doctest + In [2]: [d.size for d in green.d_classes] + Out[2]: [6, 18, 3] + + @doctest + In [3]: [int(group.order()) for _, group in green.maximal_subgroups()] + Out[3]: [6, 2, 1] + +Every finite semigroup has a unique minimal ideal: the J-class below all +others, which for transformations consists of the elements of least rank +:cite:`Howie1995`. ``S`` contains a constant map, so the automaton is +synchronizing, exactly when that rank is one. +:meth:`~sofic.algebra.transformation.TransformationSemigroup.is_synchronizing` +decides this from the generators by merging pairs of states +:cite:`Volkov2008`. :meth:`~sofic.algebra.green.GreenStructure.schutzenberger_group` +gives the Schützenberger group of any H-class, regular or not. + +.. ipython:: + + @doctest + In [1]: green.minimal_rank(), t3.is_synchronizing() + Out[1]: (1, True) + +.. autofunction:: sofic.algebra.green.green_structure +.. autoclass:: sofic.algebra.green.GreenStructure + :members: +.. autoclass:: sofic.algebra.green.DClass + :members: + +Syntactic monoids and star-free languages +========================================= + +The syntactic monoid of a regular language is the transition monoid of its +minimal complete DFA :cite:`Pin1986`. Because it is the smallest monoid that +recognizes the language, its properties belong to the language rather than to +any particular automaton. A language is star-free iff its syntactic monoid is +aperiodic :cite:`Schutzenberger1965`, which is also iff it is first-order +definable :cite:`McNaughtonPapert1971`. + +.. ipython:: + + In [1]: from sofic.algebra import is_star_free, syntactic_monoid + + @doctest + In [2]: is_star_free("(ab)*", alphabet="ab"), is_star_free("(aa)*", alphabet="a") + Out[2]: (True, False) + +:func:`~sofic.algebra.syntactic.factor_monoid` makes every state initial and +accepting, so it describes the words labeling paths in a machine. For a shift +presentation that is the factor language ``B(X)``; +:func:`~sofic.shifts.invariants.factor_syntactic_monoid` and +:func:`~sofic.shifts.invariants.factor_krohn_rhodes_complexity` package it as +presentation-independent shift invariants. The Wheeler constructions use the +same test to say whether a failed search was doomed from the start, since +Wheeler languages are star-free :cite:`Alanko2021`. + +Language varieties +================== + +Eilenberg's variety theorem matches classes of regular languages with +classes of finite monoids, so a language's class can be read off identities +in its syntactic monoid or semigroup :cite:`Eilenberg1976` :cite:`Pin1986`. +:mod:`sofic.algebra.varieties` decides the classic ones: + +* Piecewise testable (scattered subwords up to some length) iff the syntactic + monoid is J-trivial :cite:`Simon1975`. +* Locally testable (prefixes, suffixes, and factors up to some length) iff + the syntactic semigroup is locally idempotent and commutative + :cite:`BrzozowskiSimon1973` :cite:`McNaughton1974`. +* Definite or reverse definite (a bounded suffix or prefix decides + membership) iff ``se = e`` or ``es = e`` for every idempotent ``e`` + :cite:`Pin1986`. +* Group languages and commutative languages, plus R-, L-, and J-triviality. + +.. ipython:: + + In [1]: from sofic.algebra import is_locally_testable, is_piecewise_testable, is_definite + + @doctest + In [2]: is_locally_testable("(ab)*", alphabet="ab"), is_piecewise_testable("(ab)*", alphabet="ab") + Out[2]: (True, False) + + @doctest + In [3]: is_definite("(a|b)*a", alphabet="ab") + Out[3]: True + +Shifts of finite type have locally testable factor languages, which +``is_locally_testable(factor_monoid(shift, monoid=False))`` confirms. The even +shift's language is not locally testable. + +.. automodule:: sofic.algebra.varieties + :members: + :noindex: + +.. autofunction:: sofic.algebra.syntactic.syntactic_monoid +.. autofunction:: sofic.algebra.syntactic.transition_monoid +.. autofunction:: sofic.algebra.syntactic.factor_monoid +.. autofunction:: sofic.algebra.syntactic.is_star_free +.. autofunction:: sofic.shifts.invariants.factor_syntactic_monoid +.. autofunction:: sofic.shifts.invariants.factor_krohn_rhodes_complexity diff --git a/docs/examples.rst b/docs/examples.rst index ba4c923..fa0f72e 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -258,6 +258,35 @@ Example In [3]: golden_mean(0.5).entropy_rate() Out[3]: 0.6666666666666665 +Algebraic structure +~~~~~~~~~~~~~~~~~~~ + +The transition semigroup of an ε-machine depends only on its support. The +golden mean process is a shift of finite type, so its semigroup is aperiodic +and its Krohn-Rhodes complexity is zero. The even process has to count ``1``\ s +modulo two, so it needs one group level (see :doc:`algebra/algebra`): + +.. ipython:: + + In [1]: from sofic.algebra import TransformationSemigroup, krohn_rhodes_complexity + + @doctest + In [2]: golden_mean().krohn_rhodes_complexity() + Out[2]: ComplexityBounds(lower=0, upper=0, exact=0) + + @doctest + In [3]: even_process().krohn_rhodes_complexity() + Out[3]: ComplexityBounds(lower=1, upper=1, exact=1) + + In [4]: mod3 = TransformationSemigroup({"0": (0, 1, 2), "1": (1, 2, 0)}) + + @doctest + In [5]: [level.order for level in mod3.holonomy().levels] + Out[5]: [3] + +A counter modulo three is a single cyclic group of order three: one level, +complexity one. + API --- diff --git a/docs/index.rst b/docs/index.rst index 1fbd238..adc4327 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -20,8 +20,8 @@ a consistent, composable Python API built on NumPy, SciPy, and NetworkX. Every model is a graph-backed state machine, so the same objects support construction, validation, simulation, visualization, (de)serialization, and a large library of structural and information-theoretic measures. The model -families group into three presentations of a process and a set of tools for -inferring them from data: +families group into three presentations of a process, plus tools for their +algebraic structure and for inferring them from data: * **Stochastic generators** (:mod:`sofic.generators`) — Markov chains, hidden Markov models (Mealy and Moore presentations), ε-machines and their @@ -35,6 +35,9 @@ inferring them from data: * **Symbolic shifts** (:mod:`sofic.shifts`) — shifts of finite type, sofic shifts, topological Markov chains, Dyck and sofic-Dyck shifts, and their covers. These describe the *support* (set of allowed sequences) of a process. +* **Algebra** (:mod:`sofic.algebra`) — transformation semigroups, syntactic + monoids, Green's relations, and the Krohn-Rhodes holonomy decomposition of + any deterministic presentation, with complexity bounds and prime factors. * **Inference** (:mod:`sofic.inference`) — Bayesian structural inference for Markov chains, ε-machines, and stack HMMs via exact conjugate Dirichlet evidences, with optional PyMC backends. @@ -56,6 +59,7 @@ Contents: generators/generators automata/automata shifts/shifts + algebra/algebra inference/inference viz examples diff --git a/docs/references.bib b/docs/references.bib index 0f31ad3..049780c 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -1421,3 +1421,195 @@ @article{Grigorescu2024 year = {2024}, doi = {10.1109/TIT.2024.3411919}, } + +@article{KrohnRhodes1965, + author = {Krohn, Kenneth and Rhodes, John}, + title = {Algebraic Theory of Machines. {I}. Prime Decomposition Theorem for Finite Semigroups and Machines}, + journal = {Transactions of the American Mathematical Society}, + volume = {116}, + pages = {450--464}, + year = {1965}, + doi = {10.1090/S0002-9947-1965-0188316-1}, +} + +@article{Zeiger1967, + author = {Zeiger, H. Paul}, + title = {Cascade Synthesis of Finite-State Machines}, + journal = {Information and Control}, + volume = {10}, + number = {4}, + pages = {419--433}, + year = {1967}, +} + +@book{Eilenberg1976, + author = {Eilenberg, Samuel}, + title = {Automata, Languages, and Machines}, + volume = {B}, + publisher = {Academic Press}, + year = {1976}, +} + +@article{Schutzenberger1965, + author = {Sch{\"u}tzenberger, Marcel-Paul}, + title = {On Finite Monoids Having Only Trivial Subgroups}, + journal = {Information and Control}, + volume = {8}, + number = {2}, + pages = {190--194}, + year = {1965}, +} + +@book{McNaughtonPapert1971, + author = {McNaughton, Robert and Papert, Seymour}, + title = {Counter-Free Automata}, + publisher = {MIT Press}, + year = {1971}, +} + +@book{RhodesSteinberg2009, + author = {Rhodes, John and Steinberg, Benjamin}, + title = {The q-theory of Finite Semigroups}, + series = {Springer Monographs in Mathematics}, + publisher = {Springer}, + year = {2009}, +} + +@misc{EgriNagyNehaniv2015, + author = {Egri-Nagy, Attila and Nehaniv, Chrystopher L.}, + title = {Computational Holonomy Decomposition of Transformation Semigroups}, + year = {2015}, + eprint = {1508.06345}, + archivePrefix = {arXiv}, + primaryClass = {math.GR}, +} + +@book{RhodesNehanivHirsch2010, + author = {Rhodes, John}, + editor = {Nehaniv, Chrystopher L.}, + title = {Applications of Automata Theory and Algebra: Via the Mathematical Theory of Complexity to Biology, Physics, Psychology, Philosophy, and Games}, + note = {With a foreword by Morris W. Hirsch}, + publisher = {World Scientific}, + year = {2010}, +} + +@article{Green1951, + author = {Green, J. A.}, + title = {On the Structure of Semigroups}, + journal = {Annals of Mathematics}, + volume = {54}, + number = {1}, + pages = {163--172}, + year = {1951}, +} + +@book{Pin1986, + author = {Pin, Jean-{\'E}ric}, + title = {Varieties of Formal Languages}, + publisher = {North Oxford Academic}, + year = {1986}, +} + +@article{HenckellRhodesSteinberg2010, + author = {Henckell, Karsten and Rhodes, John and Steinberg, Benjamin}, + title = {A Profinite Approach to Stable Pairs}, + journal = {International Journal of Algebra and Computation}, + volume = {20}, + number = {2}, + pages = {269--285}, + year = {2010}, +} + +@article{DeDeo2011, + author = {DeDeo, Simon}, + title = {Effective Theories for Circuits and Automata}, + journal = {Chaos}, + volume = {21}, + number = {3}, + pages = {037106}, + year = {2011}, + doi = {10.1063/1.3640747}, +} + +@article{Maler1995, + author = {Maler, Oded}, + title = {A Decomposition Theorem for Probabilistic Transition Systems}, + journal = {Theoretical Computer Science}, + volume = {145}, + pages = {391}, + year = {1995}, +} + +@book{LevinPeresWilmer2009, + author = {Levin, David A. and Peres, Yuval and Wilmer, Elizabeth L.}, + title = {Markov Chains and Mixing Times}, + publisher = {American Mathematical Society}, + year = {2009}, +} + +@inproceedings{Simon1975, + author = {Simon, Imre}, + title = {Piecewise Testable Events}, + booktitle = {Automata Theory and Formal Languages}, + series = {Lecture Notes in Computer Science}, + volume = {33}, + pages = {214--222}, + publisher = {Springer}, + year = {1975}, +} + +@article{BrzozowskiSimon1973, + author = {Brzozowski, Janusz A. and Simon, Imre}, + title = {Characterizations of Locally Testable Events}, + journal = {Discrete Mathematics}, + volume = {4}, + number = {3}, + pages = {243--271}, + year = {1973}, +} + +@article{McNaughton1974, + author = {McNaughton, Robert}, + title = {Algebraic Decision Procedures for Local Testability}, + journal = {Mathematical Systems Theory}, + volume = {8}, + number = {1}, + pages = {60--76}, + year = {1974}, +} + +@article{Linton1998, + author = {Linton, Steve A. and Pfeiffer, G{\"o}tz and Robertson, Edmund F. and Ru{\v{s}}kuc, Nik}, + title = {Groups and Actions in Transformation Semigroups}, + journal = {Mathematische Zeitschrift}, + volume = {228}, + pages = {435--450}, + year = {1998}, +} + +@article{MandelSimon1977, + author = {Mandel, Arnaldo and Simon, Imre}, + title = {On Finite Semigroups of Matrices}, + journal = {Theoretical Computer Science}, + volume = {5}, + number = {2}, + pages = {101--111}, + year = {1977}, +} + +@article{Jacob1977, + author = {Jacob, G{\'e}rard}, + title = {Un algorithme calculant le cardinal, fini ou infini, des demi-groupes de matrices}, + journal = {Theoretical Computer Science}, + volume = {5}, + number = {2}, + pages = {183--204}, + year = {1977}, +} + +@book{Howie1995, + author = {Howie, John M.}, + title = {Fundamentals of Semigroup Theory}, + publisher = {Oxford University Press}, + year = {1995}, +} \ No newline at end of file diff --git a/sofic/algebra/__init__.py b/sofic/algebra/__init__.py new file mode 100644 index 0000000..47fd22e --- /dev/null +++ b/sofic/algebra/__init__.py @@ -0,0 +1,104 @@ +"""Algebraic automata theory: transformation semigroups and Krohn-Rhodes decompositions. + +* :mod:`~sofic.algebra.transformation` -- transformations and the semigroups + they generate. +* :mod:`~sofic.algebra.green` -- Green's relations, egg boxes, maximal + subgroups. +* :mod:`~sofic.algebra.syntactic` -- syntactic monoids and star-freeness. +* :mod:`~sofic.algebra.holonomy` and :mod:`~sofic.algebra.cascade` -- the + holonomy decomposition as a verifiable cascade product. +* :mod:`~sofic.algebra.complexity` -- Krohn-Rhodes complexity bounds and + prime factors. +* :mod:`~sofic.algebra.information` -- the stationary state entropy split by + holonomy level. + +Group-theoretic outputs (holonomy and maximal subgroups as permutation +groups, composition factors) need the optional ``sympy`` dependency, which +the ``sofic[symbolic]`` extra installs. +""" + +from sofic.algebra.adapters import semiautomaton, transformation_semigroup +from sofic.algebra.belief import mixed_state_semigroup, projective_semigroup +from sofic.algebra.cascade import CascadeProduct, LevelAction +from sofic.algebra.complexity import ( + ComplexityBounds, + PrimeDecomposition, + PrimeFactor, + is_full_transformation_monoid, + krohn_rhodes_complexity, + prime_decomposition, +) +from sofic.algebra.green import DClass, GreenStructure, green_structure +from sofic.algebra.holonomy import EffectiveTheory, HolonomyDecomposition, HolonomyLevel, holonomy_decomposition +from sofic.algebra.information import HolonomyInformation, holonomy_information +from sofic.algebra.random_maps import ( + RandomMapRepresentation, + WalkInformation, + classify_elements, + coupled_random_maps, + element_walk, + independent_random_maps, + random_maps, + walk_information, +) +from sofic.algebra.syntactic import factor_monoid, is_star_free, syntactic_monoid, transition_monoid +from sofic.algebra.transformation import SINK, Transformation, TransformationSemigroup +from sofic.algebra.varieties import ( + is_commutative_language, + is_definite, + is_group_language, + is_j_trivial, + is_l_trivial, + is_locally_testable, + is_piecewise_testable, + is_r_trivial, + is_reverse_definite, +) + +__all__ = [ + "SINK", + "CascadeProduct", + "ComplexityBounds", + "DClass", + "EffectiveTheory", + "GreenStructure", + "HolonomyDecomposition", + "HolonomyInformation", + "HolonomyLevel", + "LevelAction", + "PrimeDecomposition", + "PrimeFactor", + "RandomMapRepresentation", + "Transformation", + "TransformationSemigroup", + "WalkInformation", + "classify_elements", + "coupled_random_maps", + "element_walk", + "factor_monoid", + "green_structure", + "holonomy_decomposition", + "holonomy_information", + "independent_random_maps", + "is_commutative_language", + "is_definite", + "is_full_transformation_monoid", + "is_group_language", + "is_j_trivial", + "is_l_trivial", + "is_locally_testable", + "is_piecewise_testable", + "is_r_trivial", + "is_reverse_definite", + "is_star_free", + "krohn_rhodes_complexity", + "mixed_state_semigroup", + "prime_decomposition", + "projective_semigroup", + "random_maps", + "semiautomaton", + "syntactic_monoid", + "transformation_semigroup", + "transition_monoid", + "walk_information", +] diff --git a/sofic/algebra/_groups.py b/sofic/algebra/_groups.py new file mode 100644 index 0000000..cb37c3b --- /dev/null +++ b/sofic/algebra/_groups.py @@ -0,0 +1,109 @@ +"""Permutation-group helpers backed by :mod:`sympy.combinatorics` (optional).""" + +from __future__ import annotations + +from collections.abc import Iterable +from typing import Any + + +def require_sympy() -> Any: + try: + import sympy.combinatorics as combinatorics + except ImportError as exc: # pragma: no cover - exercised only without sympy + raise ImportError( + "group-theoretic parts of sofic.algebra need sympy; install with `pip install sofic[symbolic]`" + ) from exc + return combinatorics + + +def permutation_group(perms: Iterable[tuple[int, ...]], degree: int) -> Any: + """Return the sympy ``PermutationGroup`` generated by ``perms`` on ``degree`` points. + + Generators already in the group so far are skipped, so a large redundant + generating set (for example every element of a subgroup) collapses to at + most ``log2 |G|`` generators. + """ + combinatorics = require_sympy() + permutation = combinatorics.Permutation + identity = permutation(list(range(degree)), size=degree) + gens: list[Any] = [] + group = combinatorics.PermutationGroup([identity]) + for perm in dict.fromkeys(perms): + p = permutation(list(perm), size=degree) + if p.is_Identity or group.contains(p): + continue + gens.append(p) + group = combinatorics.PermutationGroup(gens) + return group + + +def composition_factor_orders(group: Any) -> list[int]: + """Orders of the Jordan-Hölder composition factors, ascending, with multiplicity. + + Solvable groups have cyclic factors of prime order, read off the order. + Otherwise abelian layers of the derived series are peeled off; a perfect + group is simple iff every nontrivial conjugacy class has normal closure + the whole group, and else splits through a proper normal subgroup and + the quotient (built as the action on cosets). + """ + order = int(group.order()) + if order == 1: + return [] + if group.is_solvable: + return _prime_factors(order) + derived = group.derived_subgroup() + if derived.order() < order: + return sorted(_prime_factors(order // int(derived.order())) + composition_factor_orders(derived)) + normal = _proper_normal_subgroup(group) + if normal is None: + return [order] + return sorted(composition_factor_orders(normal) + composition_factor_orders(_quotient(group, normal))) + + +def _prime_factors(n: int) -> list[int]: + factors = [] + p = 2 + while p * p <= n: + while n % p == 0: + factors.append(p) + n //= p + p += 1 + if n > 1: + factors.append(n) + return factors + + +def _proper_normal_subgroup(group: Any) -> Any: + order = group.order() + for cls in group.conjugacy_classes(): + x = next(iter(cls)) + if x.is_Identity: + continue + closure = group.normal_closure(x) + if closure.order() < order: + return closure + return None + + +def _quotient(group: Any, normal: Any) -> Any: + """``G / N`` as the permutation group of ``G`` acting on the right cosets of ``N``.""" + combinatorics = require_sympy() + reps = [group.identity] + + def coset(y: Any) -> int: + for i, r in enumerate(reps): + if normal.contains(y * ~r): + return i + reps.append(y) + return len(reps) - 1 + + images: dict[int, list[int]] = {} + i = 0 + while i < len(reps): + for k, g in enumerate(group.generators): + images.setdefault(k, []).append(coset(reps[i] * g)) + i += 1 + degree = len(reps) + return combinatorics.PermutationGroup( + [combinatorics.Permutation(images[k], size=degree) for k in range(len(group.generators))] + ) diff --git a/sofic/algebra/adapters.py b/sofic/algebra/adapters.py new file mode 100644 index 0000000..acf4d7f --- /dev/null +++ b/sofic/algebra/adapters.py @@ -0,0 +1,108 @@ +"""Turn sofic models into transformation semigroups. + +Every deterministic (unifilar) symbol-labeled presentation -- a DFA, a +:class:`~sofic.automata.unifilar.UnifilarAutomaton`, a unifilar HMM or +ε-machine, or a right-resolving shift presentation -- has an underlying +semiautomaton: one partial map of the states per symbol. Completing the maps +with an absorbing sink gives the transition semigroup. +""" + +from __future__ import annotations + +from collections.abc import Hashable, Iterable +from typing import Any + +from sofic.algebra.transformation import DEFAULT_MAX_SIZE, TransformationSemigroup +from sofic.graph import ATTR_EMISSION_EDGE, ATTR_PROB, ATTR_SYMBOL, EPSILON + +__all__ = ["semiautomaton", "transformation_semigroup"] + + +def semiautomaton( + machine: Any, + *, + alphabet: Iterable[Hashable] | None = None, +) -> tuple[list[Hashable], dict[Hashable, dict[Hashable, Hashable]]]: + """Return ``(states, maps)`` for a deterministic symbol-labeled machine. + + ``maps[symbol][state]`` is the unique successor. Edges with probability + zero are ignored, so a stochastic generator contributes its support. + Raises ``ValueError`` if some state has two successors on one symbol or + an ε-edge. + """ + states = list(machine.states()) + symbols: dict[Hashable, None] = dict.fromkeys(alphabet) if alphabet is not None else {} + if alphabet is None: + for attr in ("input_alphabet", "symbol_alphabet", "observation_alphabet"): + declared = getattr(machine, attr, None) + if declared: + symbols.update(dict.fromkeys(_sorted(declared))) + break + maps: dict[Hashable, dict[Hashable, Hashable]] = {s: {} for s in symbols} + for transition in machine.transitions(): + symbol = transition.data.get(ATTR_SYMBOL, transition.data.get(ATTR_EMISSION_EDGE)) + if symbol is None: + continue + if symbol is EPSILON: + raise ValueError("ε-transitions have no transition semigroup; remove them first") + prob = transition.data.get(ATTR_PROB) + if prob is not None and _is_zero(prob): + continue + table = maps.setdefault(symbol, {}) + if table.get(transition.source, transition.target) != transition.target: + raise ValueError(f"not unifilar: state {transition.source!r} has two successors on {symbol!r}") + table[transition.source] = transition.target + return states, maps + + +def transformation_semigroup( + machine: Any, + *, + alphabet: Iterable[Hashable] | None = None, + monoid: bool = True, + max_size: int = DEFAULT_MAX_SIZE, +) -> TransformationSemigroup: + """Return the transition semigroup of a deterministic presentation. + + Accepts a :class:`TransformationSemigroup` (returned unchanged), a DFA or + other deterministic labeled automaton, a unifilar HMM or ε-machine + (probabilities dropped, support kept), or a shift presentation. A + non-right-resolving :class:`~sofic.shifts.sofic.SoficShift` is replaced by + its right Fischer cover, the canonical right-resolving presentation + :cite:`LindMarcus1995`; other presentations are used as given, so the + result describes *that* presentation. For a presentation-independent + invariant of a language use :func:`~sofic.algebra.syntactic.syntactic_monoid`. + + Undefined moves go to an added absorbing ``SINK`` state. + """ + if isinstance(machine, TransformationSemigroup): + return machine + from sofic.shifts.sofic import SoficShift + + if isinstance(machine, SoficShift) and not machine.is_unifilar(): + from sofic.shifts.cover_construction import right_fischer_cover + + machine = right_fischer_cover(machine) + elif hasattr(machine, "is_unifilar") and not machine.is_unifilar(): + raise ValueError(f"{type(machine).__name__} is not unifilar; determinize it first") + states, maps = semiautomaton(machine, alphabet=alphabet) + if not states: + raise ValueError("machine has no states") + if not maps: + raise ValueError(f"{type(machine).__name__} has no symbol-labeled edges") + return TransformationSemigroup.from_maps(maps, states=states, monoid=monoid, max_size=max_size) + + +def _is_zero(prob: Any) -> bool: + try: + return bool(prob == 0) + except TypeError: + return False + + +def _sorted(items: Iterable[Hashable]) -> list[Hashable]: + values = list(items) + try: + return sorted(values) # type: ignore[type-var] + except TypeError: + return sorted(values, key=repr) diff --git a/sofic/algebra/belief.py b/sofic/algebra/belief.py new file mode 100644 index 0000000..176e226 --- /dev/null +++ b/sofic/algebra/belief.py @@ -0,0 +1,99 @@ +"""Semigroups of belief (mixed-state) maps of hidden Markov models. + +An observer of an HMM with symbol matrices ``T^(x)`` updates a belief +``eta`` over hidden states by ``eta -> eta T^(x) / (eta T^(x) 1)``. The +matrices ``T^(w)`` of words, taken up to a positive scalar, form a semigroup +that contains every such update. If that semigroup is finite, every initial +belief has a finite orbit, so the mixed-state presentation is finite from any +start. The converse fails: distinct matrices can induce the same belief map +(two rank-one matrices with the same row space both give one constant map), +and the orbit of one particular belief can be finite when the semigroup is +not. Finiteness of matrix semigroups is decidable :cite:`MandelSimon1977` +:cite:`Jacob1977`; this module only enumerates with a cap, so it proves +finiteness when the enumeration closes and is otherwise inconclusive. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Hashable +from typing import Any + +import numpy as np + +from sofic.algebra.transformation import DEFAULT_MAX_SIZE, TransformationSemigroup +from sofic.exceptions import SemigroupTooLargeError + +__all__ = ["mixed_state_semigroup", "projective_semigroup"] + + +def projective_semigroup( + hmm: Any, + *, + max_size: int = 10_000, + decimals: int = 10, +) -> TransformationSemigroup: + """Return the monoid of word matrices ``T^(w)`` up to scaling, as a transformation monoid. + + Each element is a word matrix normalized to total mass one (or the zero + matrix, for words no state can emit), rounded to ``decimals`` places to + compare. The monoid acts on itself by right multiplication by each + symbol, which is faithful, so the result is a + :class:`~sofic.algebra.transformation.TransformationSemigroup` whose + states are shortest representative words (``()`` is the identity) and + whose generators are the symbols. All of :mod:`sofic.algebra` then + applies. + + Raises :class:`~sofic.exceptions.SemigroupTooLargeError` when more than + ``max_size`` distinct elements appear: the semigroup may be infinite. + """ + from sofic.generators.matrices import symbol_matrices + + matrices = {symbol: np.asarray(m, dtype=float) for symbol, m in symbol_matrices(hmm).items()} + symbols = list(matrices) + n = next(iter(matrices.values())).shape[0] + + def normalize(matrix: np.ndarray) -> np.ndarray: + total = matrix.sum() + return matrix / total if total > 0 else np.zeros_like(matrix) + + def key(matrix: np.ndarray) -> tuple[float, ...]: + return tuple(np.round(matrix, decimals).ravel() + 0.0) + + identity = np.eye(n) / n + elements = [identity] + words: list[tuple[Hashable, ...]] = [()] + index = {key(identity): 0} + right: list[list[int]] = [] + queue = deque([0]) + while queue: + i = queue.popleft() + row = [] + for symbol in symbols: + product = normalize(elements[i] @ matrices[symbol]) + k = key(product) + if k not in index: + if len(elements) >= max_size: + raise SemigroupTooLargeError( + f"more than {max_size} belief-map matrices; the semigroup may be infinite" + ) + index[k] = len(elements) + elements.append(product) + words.append((*words[i], symbol)) + queue.append(index[k]) + row.append(index[k]) + right.append(row) + generators = {symbol: tuple(right[i][s] for i in range(len(elements))) for s, symbol in enumerate(symbols)} + return TransformationSemigroup(generators, states=words, max_size=max(DEFAULT_MAX_SIZE, len(elements))) + + +def mixed_state_semigroup(hmm: Any, **kwargs: Any) -> TransformationSemigroup: + """Return the transition semigroup of ``hmm``'s mixed-state presentation. + + The mixed-state presentation is unifilar, so its symbols act as partial + maps of the reachable beliefs; keyword arguments go to + :meth:`~sofic.generators.mealy.MealyHMM.mixed_state_presentation`. + """ + from sofic.algebra.adapters import transformation_semigroup + + return transformation_semigroup(hmm.mixed_state_presentation(**kwargs)) diff --git a/sofic/algebra/cascade.py b/sofic/algebra/cascade.py new file mode 100644 index 0000000..4f87f6d --- /dev/null +++ b/sofic/algebra/cascade.py @@ -0,0 +1,209 @@ +"""Cascade products with explicit dependency functions. + +A cascade over components ``(X_1, S_1), ..., (X_n, S_n)`` acts on tuples +``(x_1, ..., x_n)`` coordinatewise: an input applies ``d_i(x_1, ..., x_{i-1}) +∈ S_i`` to ``x_i``, so information only flows down the hierarchy +:cite:`EgriNagyNehaniv2015`. It is the finite, explicitly tabulated part of +the iterated wreath product that a particular decomposition actually uses. +""" + +from __future__ import annotations + +from collections.abc import Hashable, Iterable, Mapping, Sequence +from dataclasses import dataclass +from typing import TYPE_CHECKING, Literal, cast + +from sofic.exceptions import SemigroupTooLargeError + +if TYPE_CHECKING: + from sofic.algebra.holonomy import HolonomyDecomposition + from sofic.algebra.transformation import TransformationSemigroup + from sofic.automata.dfa import DFA + +__all__ = ["CascadeProduct", "LevelAction", "build_holonomy_cascade"] + +Coordinates = tuple[int | None, ...] + + +@dataclass(frozen=True) +class LevelAction: + """A component transformation: a reset to ``value`` or a permutation. + + ``mapping`` lists the moved states of a permutation; every other state is + fixed. A reset to ``None`` clears a coordinate the new chain does not use. + """ + + kind: Literal["constant", "permutation"] + value: int | None = None + mapping: tuple[tuple[int, int], ...] = () + + @classmethod + def constant(cls, value: int | None) -> LevelAction: + return cls("constant", value=value) + + @classmethod + def permutation(cls, mapping: Mapping[int, int]) -> LevelAction: + return cls("permutation", mapping=tuple(sorted((a, b) for a, b in mapping.items() if a != b))) + + def apply(self, x: int | None) -> int | None: + if self.kind == "constant": + return self.value + return dict(self.mapping).get(x, x) if x is not None else None + + @property + def is_identity(self) -> bool: + return self.kind == "permutation" and not self.mapping + + +class CascadeProduct: + """A finite cascade product given by dependency tables. + + Parameters + ---------- + level_sizes : sequence of int + Number of (non-``None``) states of each component, top level first. + dependencies : mapping + ``dependencies[label][i]`` maps a prefix ``coords[:i]`` to the + :class:`LevelAction` applied at level ``i``. + decoding : mapping + Valid coordinate tuple to the original state it emulates. + encoding : mapping + Original state to a chosen coordinate tuple (a lift). + """ + + def __init__( + self, + level_sizes: Sequence[int], + dependencies: Mapping[Hashable, Sequence[Mapping[Coordinates, LevelAction]]], + decoding: Mapping[Coordinates, Hashable], + encoding: Mapping[Hashable, Coordinates], + ) -> None: + self.level_sizes = tuple(level_sizes) + self.dependencies = {label: tuple(dict(t) for t in tables) for label, tables in dependencies.items()} + self.decoding = dict(decoding) + self.encoding = dict(encoding) + + @property + def n_levels(self) -> int: + return len(self.level_sizes) + + @property + def states(self) -> tuple[Coordinates, ...]: + return tuple(self.decoding) + + @property + def labels(self) -> tuple[Hashable, ...]: + return tuple(self.dependencies) + + def __repr__(self) -> str: + return f"CascadeProduct(levels={list(self.level_sizes)}, states={len(self.decoding)}, inputs={len(self.dependencies)})" + + def act(self, coords: Sequence[int | None], label: Hashable) -> Coordinates: + """Apply input ``label``: level ``i`` sees only the old ``coords[:i]``.""" + coords = tuple(coords) + tables = self.dependencies[label] + return tuple(tables[i][coords[:i]].apply(coords[i]) for i in range(self.n_levels)) + + def run(self, word: Iterable[Hashable], start: Hashable) -> Coordinates: + """Read ``word`` from ``start`` (an original state or a coordinate tuple).""" + coords = cast("Coordinates", start) if start in self.decoding else self.encoding[start] + for label in word: + coords = self.act(coords, label) + return coords + + def truncate(self, k: int) -> CascadeProduct: + """Keep the top ``k`` levels: the coarse-grained "effective theory". + + Level ``i`` reads only ``coords[:i]``, so the first ``k`` levels are a + cascade in their own right and the projection onto them commutes with + every input :cite:`DeDeo2011`. In the result, ``decoding`` maps each + prefix to the frozenset of original states it covers (a *superstate*; + superstates may overlap) and ``encoding`` maps each state to the + prefix of its lift. + """ + if not 0 <= k <= self.n_levels: + raise ValueError(f"k must be between 0 and {self.n_levels}") + covers: dict[Coordinates, set[Hashable]] = {} + for coords, state in self.decoding.items(): + covers.setdefault(coords[:k], set()).add(state) + return CascadeProduct( + self.level_sizes[:k], + {label: tables[:k] for label, tables in self.dependencies.items()}, + {prefix: frozenset(states) for prefix, states in covers.items()}, + {state: coords[:k] for state, coords in self.encoding.items()}, + ) + + def as_dfa(self, initial: Hashable | None = None) -> DFA: + """Return the cascade as a complete DFA on its coordinate tuples.""" + from sofic.automata.dfa import DFA + + dfa = DFA( + input_alphabet=frozenset(self.labels), + initial_states=frozenset() if initial is None else frozenset({self.encoding[initial]}), + ) + for coords in self.decoding: + dfa.graph.add_state(coords) + for coords in self.decoding: + for label in self.labels: + dfa.add_transition(coords, self.act(coords, label), label) + return dfa + + def emulates(self, semigroup: TransformationSemigroup) -> bool: + """Return whether the cascade covers ``semigroup``'s generators. + + Checks that the valid coordinate tuples are closed under every input, + that decoding is onto the states, and that decoding commutes with the + action: ``decode(act(c, a)) == decode(c) · a``. + """ + if set(self.decoding.values()) != set(semigroup.states): + return False + gens = semigroup.generators + if set(gens) != set(self.labels): + return False + position = {s: i for i, s in enumerate(semigroup.states)} + for coords, state in self.decoding.items(): + for label, g in gens.items(): + try: + image = self.act(coords, label) + except KeyError: + return False + if image not in self.decoding: + return False + if position[self.decoding[image]] != g(position[state]): + return False + return True + + +def build_holonomy_cascade(decomposition: HolonomyDecomposition, *, max_states: int) -> CascadeProduct: + """Tabulate the holonomy cascade over all maximal tile chains. + + Raises ``RuntimeError`` if a level's action is not determined by the + coordinates above it, or does not reproduce the lifted action -- either + would mean the lift is not consistent with chain structure. + """ + sg = decomposition.semigroup + gens = sg.generators + n_levels = decomposition.n_levels + dependencies: dict[Hashable, list[dict[Coordinates, LevelAction]]] = { + label: [{} for _ in range(n_levels)] for label in gens + } + decoding: dict[Coordinates, Hashable] = {} + for chain in decomposition.chains(): + if len(decoding) >= max_states: + raise SemigroupTooLargeError(f"holonomy cascade exceeds max_states={max_states} coordinate tuples") + coords = decomposition.encode(chain) + (x,) = chain[-1] + decoding[coords] = sg.states[x] + for label, g in gens.items(): + actions = decomposition.level_actions(chain, g) + expected = decomposition.encode(decomposition.step(chain, g)) + for i, action in enumerate(actions): + table = dependencies[label][i] + previous = table.setdefault(coords[:i], action) + if previous != action: + raise RuntimeError(f"level {i} action on {label!r} depends on more than the levels above it") + if action.apply(coords[i]) != expected[i]: + raise RuntimeError(f"level {i} action on {label!r} does not reproduce the lifted chain") + encoding = {state: decomposition.coordinates(state) for state in sg.states} + sizes = [len(level.states) for level in decomposition.levels] + return CascadeProduct(sizes, dependencies, decoding, encoding) diff --git a/sofic/algebra/complexity.py b/sofic/algebra/complexity.py new file mode 100644 index 0000000..853541c --- /dev/null +++ b/sofic/algebra/complexity.py @@ -0,0 +1,190 @@ +"""Krohn-Rhodes complexity bounds and prime decompositions. + +By the Krohn-Rhodes theorem every finite transformation semigroup divides a +cascade of simple groups and flip-flops (the three-element monoid of two +resets and the identity), and the simple groups needed are exactly the simple +group divisors of ``S`` :cite:`KrohnRhodes1965`. The *group complexity* +``c(S)`` is the least number of group levels, separated by aperiodic levels, +in such a cascade. It is monotone under division, ``c(S) = 0`` exactly when +``S`` is aperiodic, and ``c(T_n) = n - 1`` for the full transformation monoid +:cite:`RhodesSteinberg2009`. Whether ``c(S)`` is computable in general is +open; only complexity one is known to be decidable +:cite:`HenckellRhodesSteinberg2010`, and that procedure is not implemented +here. This module returns bounds instead, exact when they meet. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any, NamedTuple + +from sofic.algebra._groups import composition_factor_orders +from sofic.exceptions import SemigroupTooLargeError + +if TYPE_CHECKING: + from sofic.algebra.holonomy import HolonomyDecomposition, Subset + +__all__ = [ + "ComplexityBounds", + "PrimeDecomposition", + "PrimeFactor", + "is_full_transformation_monoid", + "krohn_rhodes_complexity", + "prime_decomposition", +] + + +class ComplexityBounds(NamedTuple): + """Lower and upper bounds on Krohn-Rhodes group complexity. + + ``exact`` is the complexity when ``lower == upper`` and ``None`` otherwise. + """ + + lower: int + upper: int + exact: int | None + + +def krohn_rhodes_complexity(machine: Any, *, decomposition: HolonomyDecomposition | None = None) -> ComplexityBounds: + """Return bounds on the Krohn-Rhodes complexity of a semigroup or machine. + + * **Upper bound:** the number of holonomy levels with a nontrivial group. + The holonomy cascade emulates ``S`` and each level is a permutation + group with resets :cite:`EgriNagyNehaniv2015`. + * **Lower bound:** ``0`` when ``S`` is aperiodic, which holds iff every + holonomy group is trivial, and otherwise at least ``1``. If the + elements mapping an image set ``A`` into itself restrict to all + ``|A|**|A|`` maps of ``A``, then ``T_|A|`` divides ``S`` and so + ``c(S) >= |A| - 1`` :cite:`RhodesSteinberg2009`. This second test + enumerates ``S`` and is skipped if ``S`` is too large to enumerate, + except that ``S = T_n`` itself is recognized from its generators + (:func:`is_full_transformation_monoid`). + + Sharper lower bounds from the literature (Rhodes-Tilson chains of + essential J-classes) are not implemented. + + Examples + -------- + >>> from sofic.algebra import TransformationSemigroup + >>> krohn_rhodes_complexity(TransformationSemigroup({"a": (1, 2, 0)})) + ComplexityBounds(lower=1, upper=1, exact=1) + """ + from sofic.algebra.holonomy import holonomy_decomposition + + hd = decomposition if decomposition is not None else holonomy_decomposition(machine) + upper = sum(not level.is_trivial for level in hd.levels) + lower = 0 + if upper: + lower = max(1, _full_transformation_bound(hd)) + return ComplexityBounds(lower, upper, lower if lower == upper else None) + + +def is_full_transformation_monoid(semigroup: Any) -> bool: + """Return whether ``S`` is the full transformation monoid ``T_n``, from its generators alone. + + A product involving a non-permutation is not a permutation, and rank + never increases under composition. So ``S = T_n`` exactly when the + permutation generators generate ``S_n`` and some generator has rank + ``n - 1``; conversely ``S_n`` and any rank ``n - 1`` map generate ``T_n`` + :cite:`Howie1995`. + """ + from math import factorial + + from sofic.algebra._groups import permutation_group + + n = semigroup.degree + gens = list(semigroup.generators.values()) + if n == 1: + return True + if not any(g.rank() == n - 1 for g in gens): + return False + perms = [g.images for g in gens if g.is_permutation()] + return bool(perms) and int(permutation_group(perms, n).order()) == factorial(n) + + +def _full_transformation_bound(hd: HolonomyDecomposition) -> int: + """Largest ``k - 1`` such that ``T_k`` appears as a local action on an image set.""" + sg = hd.semigroup + if is_full_transformation_monoid(sg): + return sg.degree - 1 + try: + elements = [g.images for g in sg] + except SemigroupTooLargeError: + return 0 + best = 0 + candidates: list[Subset] = sorted((s for s in hd.image_set if len(s) > 2), key=len, reverse=True) + for subset in candidates: + k = len(subset) + if k - 1 <= best or k**k > len(elements): + continue + ordered = sorted(subset) + position = {x: i for i, x in enumerate(ordered)} + local = {tuple(position[e[x]] for x in ordered) for e in elements if all(e[x] in position for x in ordered)} + if len(local) == k**k: + best = k - 1 + return best + + +@dataclass(frozen=True) +class PrimeFactor: + """A simple group occurring as a composition factor of a holonomy group. + + ``kind`` is ``"cyclic"`` (order prime) or ``"nonabelian simple"``. A + nonabelian simple group is reported by order only, since order does not + always determine it (``|A_8| = |PSL(3, 4)|``). + """ + + kind: str + order: int + depth: int + representative: frozenset[int] + + +@dataclass(frozen=True) +class PrimeDecomposition: + """The Krohn-Rhodes primes read off a holonomy decomposition. + + ``groups`` lists the composition factors of every holonomy group, with + multiplicity and position. ``needs_flip_flop`` is true when ``S`` + contains a non-permutation, so resets (flip-flops) are required. + """ + + groups: tuple[PrimeFactor, ...] + needs_flip_flop: bool + + def simple_group_orders(self) -> tuple[int, ...]: + """Distinct orders of the simple group factors, ascending.""" + return tuple(sorted({f.order for f in self.groups})) + + +def prime_decomposition(machine: Any, *, decomposition: HolonomyDecomposition | None = None) -> PrimeDecomposition: + """Return the simple-group and flip-flop components needed to build ``S``. + + The Jordan-Hölder factors of the holonomy groups are simple group + divisors of ``S`` and, by the prime decomposition theorem, every simple + group divisor of ``S`` divides one of them :cite:`KrohnRhodes1965` + :cite:`EgriNagyNehaniv2015`. + + Examples + -------- + >>> from sofic.algebra import TransformationSemigroup + >>> pd = prime_decomposition(TransformationSemigroup({"t": (1, 0, 2), "c": (1, 2, 0), "e": (0, 0, 2)})) + >>> pd.simple_group_orders(), pd.needs_flip_flop + ((2, 3), True) + """ + from sofic.algebra.holonomy import holonomy_decomposition + + hd = decomposition if decomposition is not None else holonomy_decomposition(machine) + factors = [] + for level in hd.levels: + if level.is_trivial: + continue + for rep, group in zip(level.representatives, level.groups, strict=True): + for order in composition_factor_orders(group): + kind = "cyclic" if _is_prime(order) else "nonabelian simple" + factors.append(PrimeFactor(kind, order, level.depth, rep)) + return PrimeDecomposition(tuple(factors), not hd.semigroup.is_group()) + + +def _is_prime(n: int) -> bool: + return n > 1 and all(n % p for p in range(2, int(n**0.5) + 1)) diff --git a/sofic/algebra/green.py b/sofic/algebra/green.py new file mode 100644 index 0000000..1d10eb6 --- /dev/null +++ b/sofic/algebra/green.py @@ -0,0 +1,207 @@ +"""Green's relations of a finite transformation semigroup. + +Two elements are R-related when they generate the same right ideal +(``sS¹ = tS¹``), L-related for left ideals, J-related for two-sided ideals, +and H = R ∩ L. In a finite semigroup D = J, every D-class is a rectangular +"egg box" of H-classes, and the H-class of an idempotent is a maximal subgroup +:cite:`Green1951` :cite:`Pin1986`. Here R- and L-classes are the strongly +connected components of the right and left Cayley graphs and J-classes those +of their union, which is exactly the "reachable by multiplying on that side" +characterization. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import networkx as nx + +from sofic.algebra._groups import permutation_group + +if TYPE_CHECKING: + from sofic.algebra.transformation import TransformationSemigroup + +__all__ = ["DClass", "GreenStructure", "green_structure"] + + +@dataclass(frozen=True) +class DClass: + """One D-class laid out as an egg box. + + ``cells[r][c]`` is the H-class at the intersection of the ``r``-th R-class + and the ``c``-th L-class, as a tuple of element indices. + """ + + r_classes: tuple[tuple[int, ...], ...] + l_classes: tuple[tuple[int, ...], ...] + cells: tuple[tuple[tuple[int, ...], ...], ...] + idempotents: frozenset[int] + + @property + def is_regular(self) -> bool: + return bool(self.idempotents) + + @property + def size(self) -> int: + return sum(len(r) for r in self.r_classes) + + @property + def h_class_size(self) -> int: + return len(self.cells[0][0]) + + def elements(self) -> tuple[int, ...]: + return tuple(sorted(i for r in self.r_classes for i in r)) + + +class GreenStructure: + """Green's R, L, H and D (= J) classes of a :class:`TransformationSemigroup`. + + ``j_order`` is the condensation DAG of the two-sided Cayley graph on + D-class ids, with an edge from a class to each class strictly below it in + one multiplication step. + """ + + def __init__(self, semigroup: TransformationSemigroup) -> None: + self.semigroup = semigroup + right = semigroup.right_cayley() + left = semigroup.left_cayley() + n = len(right) + self.r_of, self.r_classes = _scc_labels(n, right) + self.l_of, self.l_classes = _scc_labels(n, left) + both = [r + l for r, l in zip(right, left, strict=True)] + self.d_of, d_members = _scc_labels(n, both) + + h_map: dict[tuple[int, int], list[int]] = {} + for i in range(n): + h_map.setdefault((self.r_of[i], self.l_of[i]), []).append(i) + self.h_classes = [tuple(v) for v in h_map.values()] + self.h_of = [0] * n + for k, members in enumerate(self.h_classes): + for i in members: + self.h_of[i] = k + + idempotents = set(semigroup.idempotents()) + self.d_classes: list[DClass] = [] + for members in d_members: + r_ids = list(dict.fromkeys(self.r_of[i] for i in members)) + l_ids = list(dict.fromkeys(self.l_of[i] for i in members)) + cells = tuple(tuple(tuple(h_map[(r, c)]) for c in l_ids) for r in r_ids) + self.d_classes.append( + DClass( + r_classes=tuple(self.r_classes[r] for r in r_ids), + l_classes=tuple(self.l_classes[c] for c in l_ids), + cells=cells, + idempotents=frozenset(i for i in members if i in idempotents), + ) + ) + + graph = nx.DiGraph() + graph.add_nodes_from(range(len(d_members))) + for i, row in enumerate(both): + for j in row: + if self.d_of[i] != self.d_of[j]: + graph.add_edge(self.d_of[i], self.d_of[j]) + self.j_order = graph + + def is_aperiodic(self) -> bool: + """Return whether every H-class is trivial.""" + return all(len(h) == 1 for h in self.h_classes) + + def maximal_subgroup(self, idempotent: int) -> Any: + """Return the maximal subgroup ``H_e`` as a sympy group on ``im(e)``. + + Each element of ``H_e`` restricts to a permutation of the image of + ``e``, and the restriction is faithful, so the group is returned as a + permutation group on ``len(im(e))`` points (the image in sorted order). + """ + e = self.semigroup[idempotent] + if not e.is_idempotent(): + raise ValueError(f"element {idempotent} is not idempotent") + image = sorted(e.image()) + members = self.h_classes[self.h_of[idempotent]] + return permutation_group((self.semigroup[i].restrict(image) for i in members), len(image)) + + def maximal_subgroups(self) -> list[tuple[int, Any]]: + """One ``(idempotent, group)`` pair per regular D-class. + + All maximal subgroups in a D-class are isomorphic, so one + representative idempotent per regular D-class suffices. + """ + return [(min(d.idempotents), self.maximal_subgroup(min(d.idempotents))) for d in self.d_classes if d.is_regular] + + def minimal_ideal(self) -> DClass: + """Return the minimal ideal ``K(S)``, the unique J-class below all others. + + Every finite semigroup has one; it is a completely simple semigroup + and, for a transformation semigroup, consists of the elements of + minimal rank :cite:`Howie1995`. + """ + (bottom,) = [d for d in self.j_order.nodes if self.j_order.out_degree(d) == 0] + return self.d_classes[bottom] + + def minimal_rank(self) -> int: + """Rank of the elements of the minimal ideal (``1`` iff ``S`` contains a constant map).""" + return self.semigroup[self.minimal_ideal().elements()[0]].rank() + + def schutzenberger_group(self, element: int) -> Any: + """Return the Schützenberger group of an element's H-class, regular or not. + + The right Schützenberger group of ``H = H_x`` is induced by + ``{u ∈ S¹ : x u ∈ H}``; such ``u`` permute ``im(x)``, and ``xu = xu'`` + iff they agree there, so the group is faithful on ``im(x)`` with order + ``|H|`` :cite:`Howie1995`. It is returned as a sympy group on the + sorted image. For an idempotent it is the maximal subgroup + (:meth:`maximal_subgroup`). Compare the *permutator group* of + ``im(x)``, induced by its whole setwise stabilizer, which can be larger + :cite:`Linton1998`. + """ + sg = self.semigroup + x = sg[element] + image = sorted(x.image()) + target = self.h_of[element] + perms = [tuple(range(len(image)))] + for u in sg: + product = x * u + if product in sg and self.h_of[sg.index(product)] == target: + perms.append(u.restrict(image)) + return permutation_group(perms, len(image)) + + def egg_box(self) -> str: + """Render every D-class as a grid of H-classes labeled by shortest words. + + Idempotents are starred; the empty word is written ``1``. + """ + sg = self.semigroup + idempotents = {i for d in self.d_classes for i in d.idempotents} + + def label(i: int) -> str: + word = sg.word(i) + text = "".join(map(str, word)) if word else "1" + return f"*{text}" if i in idempotents else text + + blocks = [] + for k, d in enumerate(self.d_classes): + rows = [" | ".join(",".join(label(i) for i in cell) for cell in row) for row in d.cells] + header = f"D{k} (size {d.size}, {'regular' if d.is_regular else 'non-regular'})" + blocks.append("\n".join([header, *rows])) + return "\n\n".join(blocks) + + +def green_structure(semigroup: TransformationSemigroup) -> GreenStructure: + """Return Green's relations of ``semigroup``; enumerates its elements.""" + return GreenStructure(semigroup) + + +def _scc_labels(n: int, adjacency: list[list[int]]) -> tuple[list[int], list[tuple[int, ...]]]: + graph = nx.DiGraph() + graph.add_nodes_from(range(n)) + for i, row in enumerate(adjacency): + graph.add_edges_from((i, j) for j in row) + labels = [0] * n + components = sorted((sorted(int(i) for i in c) for c in nx.strongly_connected_components(graph)), key=min) + classes = [tuple(component) for component in components] + for k, component in enumerate(components): + for i in component: + labels[i] = k + return labels, classes diff --git a/sofic/algebra/holonomy.py b/sofic/algebra/holonomy.py new file mode 100644 index 0000000..be1315b --- /dev/null +++ b/sofic/algebra/holonomy.py @@ -0,0 +1,512 @@ +"""Holonomy decomposition of a finite transformation semigroup. + +The holonomy method of Zeiger :cite:`Zeiger1967`, streamlined by Eilenberg +:cite:`Eilenberg1976`, emulates ``(Q, S)`` by a cascade (iterated wreath +product) whose components are permutation groups -- the *holonomy groups* -- +each augmented with constant maps (resets). This module follows the +constructive presentation of Egri-Nagy & Nehaniv :cite:`EgriNagyNehaniv2015`: + +1. The *extended image set* ``I`` is ``{Q · s : s ∈ S¹}`` together with all + singletons. +2. Two sets are equivalent when each is the image of the other; the classes + are the strongly connected components of ``S`` acting on ``I``. + *Subduction* ``P ≤ R`` (``P ⊆ R · s`` for some ``s ∈ S¹``) orders the + classes, and the *height* of a set is the length of the longest strict + subduction chain of non-singletons ending at it. The cascade has one level + per height ``1, ..., h(Q)``; levels are numbered by *depth* + ``d = h(Q) - h + 1`` so the top level is depth 1. +3. The *tiles* of a set are its maximal proper subsets in ``I``. The holonomy + group of a class representative ``R`` is the permutation group its + setwise stabilizer induces on its tiles. It is generated by round trips + ``R -> P -s-> P' -> R`` through the class (Schreier generators), so it is + computed from the generators of ``S`` alone, without enumerating ``S``. + +States are lifted to maximal tile chains ``Q = P_0 ⊃ P_1 ⊃ ... ⊃ {x}``; the +tile below ``P_i`` sits at level ``d(P_i)`` and is recorded relative to the +representative of ``P_i``. Those relative tiles are the cascade coordinates. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Hashable, Iterator, Sequence +from dataclasses import dataclass +from math import prod +from typing import TYPE_CHECKING, Any, NamedTuple, cast + +import networkx as nx + +from sofic.algebra.transformation import ( + DEFAULT_MAX_SIZE, + Transformation, + TransformationSemigroup, + _compose, + _permutation_order, +) + +if TYPE_CHECKING: + from sofic.algebra.cascade import CascadeProduct, LevelAction + +__all__ = ["EffectiveTheory", "HolonomyDecomposition", "HolonomyLevel", "holonomy_decomposition"] + +Subset = frozenset[int] +Chain = tuple[Subset, ...] +Coordinates = tuple[int | None, ...] + + +def _set_key(subset: Subset) -> tuple[int, tuple[int, ...]]: + return (-len(subset), tuple(sorted(subset))) + + +class EffectiveTheory(NamedTuple): + """The coarse-grained machine on the top levels of a holonomy cascade. + + ``automaton`` is a DFA on coordinate prefixes, ``superstates`` maps each + prefix to the original states it covers, and ``projection`` maps every + full coordinate tuple to its prefix. ``distribution`` and ``lumpable`` + are filled in only when a stochastic machine is supplied. + """ + + automaton: Any + superstates: dict[Coordinates, frozenset[Hashable]] + projection: dict[Coordinates, Coordinates] + distribution: dict[Coordinates, float] | None + lumpable: bool | None + + +@dataclass(frozen=True) +class HolonomyLevel: + """One level of the holonomy cascade. + + ``representatives[k]`` is a class representative of this depth, + ``tiles[k]`` its tiles, and ``group_generators[k]`` permutations of + ``range(len(tiles[k]))`` generating its holonomy group. The component's + state set is ``states``: pairs ``(k, t)`` meaning tile ``t`` of + representative ``k`` (the disjoint union of the tile sets). + """ + + depth: int + height: int + representatives: tuple[Subset, ...] + tiles: tuple[tuple[Subset, ...], ...] + group_generators: tuple[tuple[tuple[int, ...], ...], ...] + states: tuple[tuple[int, int], ...] + + @property + def groups(self) -> tuple[Any, ...]: + """Holonomy groups as sympy ``PermutationGroup`` objects on tile positions.""" + from sofic.algebra._groups import permutation_group + + return tuple( + permutation_group(gens, len(tiles)) for gens, tiles in zip(self.group_generators, self.tiles, strict=True) + ) + + @property + def order(self) -> int: + """Order of the level group, the direct product of its holonomy groups.""" + if self.is_trivial: + return 1 + return prod(int(g.order()) for g in self.groups) + + @property + def is_trivial(self) -> bool: + return all(all(p == tuple(range(len(p))) for p in gens) for gens in self.group_generators) + + +class HolonomyDecomposition: + """The holonomy decomposition of a :class:`TransformationSemigroup`. + + Attributes + ---------- + semigroup : TransformationSemigroup + image_set : tuple of frozenset + The extended image set, largest sets first. + height : dict + Height of every set in ``image_set`` (singletons have height 0). + levels : tuple of HolonomyLevel + Top (depth 1) first; ``len(levels) == height[Q]``. + """ + + def __init__(self, semigroup: TransformationSemigroup) -> None: + self.semigroup = semigroup + n = semigroup.degree + self._full: Subset = frozenset(range(n)) + self._gens = [(label, g) for label, g in semigroup.generators.items()] + + nodes, edges = self._orbit() + self.image_set: tuple[Subset, ...] = tuple(sorted(nodes, key=_set_key)) + graph = nx.DiGraph() + graph.add_nodes_from(nodes) + for source, moves in edges.items(): + graph.add_edges_from((source, target) for _, _, target in moves) + + classes = [frozenset(c) for c in nx.strongly_connected_components(graph) if len(next(iter(c))) > 1] + classes.sort(key=lambda c: _set_key(min(c, key=_set_key))) + self.classes: tuple[tuple[Subset, ...], ...] = tuple(tuple(sorted(c, key=_set_key)) for c in classes) + self.representative: dict[Subset, Subset] = {} + for members in self.classes: + rep = min(members, key=lambda s: tuple(sorted(s))) + for member in members: + self.representative[member] = rep + + self.height = self._heights(graph) + self.n_levels = self.height[self._full] + self.tiles = self._tiles(nodes) + self._tile_sets = {p: frozenset(ts) for p, ts in self.tiles.items()} + self.to_rep, self.from_rep = self._translations(edges) + self.levels = self._levels(edges) + self._level_index = [{state: i for i, state in enumerate(self._level_keys(level))} for level in self.levels] + + # -- construction --------------------------------------------------------- + + def _orbit(self) -> tuple[set[Subset], dict[Subset, list[tuple[Hashable, Transformation, Subset]]]]: + singletons = {frozenset({x}) for x in range(self.semigroup.degree)} + nodes = {self._full} | singletons + edges: dict[Subset, list[tuple[Hashable, Transformation, Subset]]] = {} + queue = deque(nodes) + while queue: + subset = queue.popleft() + moves = [] + for label, g in self._gens: + image = g.act(subset) + moves.append((label, g, image)) + if image not in nodes: + nodes.add(image) + queue.append(image) + edges[subset] = moves + return nodes, edges + + def _heights(self, graph: nx.DiGraph) -> dict[Subset, int]: + reps = sorted(set(self.representative.values()), key=_set_key) + reach = {r: nx.descendants(graph, r) | {r} for r in reps} + below = {a: [b for b in reps if b != a and any(b <= c for c in reach[a]) and len(b) > 1] for a in reps} + heights: dict[Subset, int] = {} + + def height(rep: Subset) -> int: + if rep not in heights: + heights[rep] = 1 + max((height(b) for b in below[rep]), default=0) + return heights[rep] + + result = {s: 0 for s in graph.nodes if len(s) == 1} + for member, rep in self.representative.items(): + result[member] = height(rep) + return result + + def _tiles(self, nodes: set[Subset]) -> dict[Subset, tuple[Subset, ...]]: + tiles: dict[Subset, tuple[Subset, ...]] = {} + for parent in nodes: + if len(parent) == 1: + continue + proper = [s for s in nodes if s < parent] + maximal = [s for s in proper if not any(s < t for t in proper)] + tiles[parent] = tuple(sorted(maximal, key=_set_key)) + return tiles + + def _translations( + self, edges: dict[Subset, list[tuple[Hashable, Transformation, Subset]]] + ) -> tuple[dict[Subset, Transformation], dict[Subset, Transformation]]: + """Maps ``P -> rep(P)`` and back, inverse to each other on ``P`` and ``rep(P)``.""" + n = self.semigroup.degree + identity = tuple(range(n)) + to_rep: dict[Subset, tuple[int, ...]] = {} + from_rep: dict[Subset, tuple[int, ...]] = {} + incoming: dict[Subset, list[tuple[tuple[int, ...], Subset]]] = {} + for source, moves in edges.items(): + for _, g, target in moves: + incoming.setdefault(target, []).append((g.images, source)) + for members in self.classes: + rep = self.representative[members[0]] + inside = set(members) + from_rep[rep] = identity + queue = deque([rep]) + while queue: + p = queue.popleft() + for _, g, q in edges[p]: + if q in inside and q not in from_rep: + from_rep[q] = _compose(from_rep[p], g.images) + queue.append(q) + to_rep[rep] = identity + queue = deque([rep]) + while queue: + q = queue.popleft() + for g_images, p in incoming.get(q, []): + if p in inside and p not in to_rep: + to_rep[p] = _compose(g_images, to_rep[q]) + queue.append(p) + for p in members: + roundtrip = _compose(to_rep[p], from_rep[p]) + order = _permutation_order(roundtrip, p) + from_rep[p] = _compose(from_rep[p], Transformation(roundtrip).power(order - 1).images) + return ( + {p: Transformation(m) for p, m in to_rep.items()}, + {p: Transformation(m) for p, m in from_rep.items()}, + ) + + def _levels(self, edges: dict[Subset, list[tuple[Hashable, Transformation, Subset]]]) -> tuple[HolonomyLevel, ...]: + by_depth: dict[int, list[Subset]] = {} + for rep in sorted(set(self.representative.values()), key=lambda s: tuple(sorted(s))): + by_depth.setdefault(self.depth(rep), []).append(rep) + levels = [] + for depth in range(1, self.n_levels + 1): + reps = tuple(by_depth.get(depth, [])) + tiles = tuple(self.tiles[r] for r in reps) + gens = tuple(self._holonomy_generators(r, edges) for r in reps) + states = tuple((k, t) for k, ts in enumerate(tiles) for t in range(len(ts))) + levels.append( + HolonomyLevel( + depth=depth, + height=self.n_levels - depth + 1, + representatives=reps, + tiles=tiles, + group_generators=gens, + states=states, + ) + ) + return tuple(levels) + + def _holonomy_generators( + self, rep: Subset, edges: dict[Subset, list[tuple[Hashable, Transformation, Subset]]] + ) -> tuple[tuple[int, ...], ...]: + tiles = self.tiles[rep] + position = {t: i for i, t in enumerate(tiles)} + perms: dict[tuple[int, ...], None] = {} + for p, r in self.representative.items(): + if r != rep: + continue + for _, g, q in edges[p]: + if self.representative.get(q) != rep: + continue + loop = self.from_rep[p] * g * self.to_rep[q] + perms[tuple(position[loop.act(t)] for t in tiles)] = None + identity = tuple(range(len(tiles))) + return tuple(p for p in perms if p != identity) or (identity,) + + @staticmethod + def _level_keys(level: HolonomyLevel) -> Iterator[tuple[Subset, Subset]]: + for k, t in level.states: + yield (level.representatives[k], level.tiles[k][t]) + + # -- geometry of sets ----------------------------------------------------- + + def depth(self, subset: Subset) -> int: + """Depth ``h(Q) - h(P) + 1`` of a set; singletons sit below the last level.""" + return self.n_levels - self.height[subset] + 1 + + def equivalence_class(self, subset: Subset) -> tuple[Subset, ...]: + rep = self.representative[subset] + return next(c for c in self.classes if rep in c) + + # -- chains and coordinates ----------------------------------------------- + + def lift(self, state: Hashable) -> Chain: + """Return the canonical maximal tile chain from ``Q`` down to ``{state}``.""" + x = self.semigroup.state_index(state) + chain = [self._full] + while len(chain[-1]) > 1: + chain.append(next(t for t in self.tiles[chain[-1]] if x in t)) + return tuple(chain) + + def chains(self) -> Iterator[Chain]: + """Yield every maximal tile chain from ``Q`` to a singleton.""" + stack: list[Chain] = [(self._full,)] + while stack: + chain = stack.pop() + if len(chain[-1]) == 1: + yield chain + continue + for tile in reversed(self.tiles[chain[-1]]): + stack.append((*chain, tile)) + + def encode(self, chain: Sequence[Subset]) -> Coordinates: + """Return the cascade coordinates of a maximal tile chain. + + Coordinate ``i`` is an index into ``levels[i].states`` or ``None`` + when no set of the chain sits at depth ``i + 1``. + """ + coords: list[int | None] = [None] * self.n_levels + for parent, tile in zip(chain, chain[1:], strict=False): + rep = self.representative[parent] + level = self.depth(parent) - 1 + coords[level] = self._level_index[level][(rep, self.to_rep[parent].act(tile))] + return tuple(coords) + + def decode(self, coords: Sequence[int | None]) -> Chain: + """Inverse of :meth:`encode`; raises ``ValueError`` on invalid coordinates.""" + if len(coords) != self.n_levels: + raise ValueError(f"expected {self.n_levels} coordinates") + chain = [self._full] + for i, c in enumerate(coords): + parent = chain[-1] + on_level = len(parent) > 1 and self.depth(parent) == i + 1 + if c is None: + if on_level: + raise ValueError(f"missing coordinate at level {i}") + continue + if not on_level: + raise ValueError(f"unexpected coordinate at level {i}") + k, t = self.levels[i].states[c] + level = self.levels[i] + if level.representatives[k] != self.representative[parent]: + raise ValueError(f"coordinate at level {i} belongs to another representative") + chain.append(self.from_rep[parent].act(level.tiles[k][t])) + if len(chain[-1]) != 1: + raise ValueError("coordinates do not reach a singleton") + return tuple(chain) + + def coordinates(self, state: Hashable) -> Coordinates: + """Return the coordinates of the canonical lift of ``state``.""" + return self.encode(self.lift(state)) + + def state(self, coords: Sequence[int | None]) -> Hashable: + """Return the original state a coordinate tuple stands for.""" + (x,) = self.decode(coords)[-1] + return self.semigroup.states[x] + + def step(self, chain: Sequence[Subset], s: Transformation) -> Chain: + """Lift the action of ``s`` to maximal tile chains. + + The images ``P_i · s`` of the chain form a chain of image sets; this + returns a maximal tile chain containing all of them (a *dominating* + chain) built top down. Below a set ``C`` it takes the first image that + is a proper subset of ``C``: that image itself if it is a tile of + ``C``, otherwise the first tile of ``C`` containing it. The choice down + to ``P · s`` depends only on the chain down to ``P``, which is what + makes the cascade's dependency functions well defined + :cite:`EgriNagyNehaniv2015`. + """ + images = [s.act(p) for p in chain] + current = self._full + result = [current] + while len(current) > 1: + target = next(d for d in images if d < current) + if target in self._tile_sets[current]: + current = target + else: + current = next(t for t in self.tiles[current] if target <= t) + result.append(current) + return tuple(result) + + def act(self, coords: Sequence[int | None], label: Hashable) -> Coordinates: + """Apply generator ``label`` to coordinates.""" + g = self.semigroup.generators[label] + return self.encode(self.step(self.decode(coords), g)) + + def level_actions(self, chain: Sequence[Subset], s: Transformation) -> tuple[LevelAction, ...]: + """Return the component action of ``s`` at each level for ``chain``. + + At level ``i`` let ``P`` and ``P'`` be the sets the old and new chains + are refining there. If ``P'`` does not sit at depth ``i + 1`` the + coordinate is cleared; if ``P · s = P'`` the action is the holonomy + permutation ``m(R -> P) s m(P' -> R)`` of the representative's tiles; + otherwise it is a reset to the new coordinate. + """ + from sofic.algebra.cascade import LevelAction + + new = self.step(chain, s) + new_coords = self.encode(new) + old_frames = self._frames(chain) + new_frames = self._frames(new) + actions = [] + for i in range(self.n_levels): + p, q = old_frames[i], new_frames[i] + if len(q) == 1 or self.depth(q) != i + 1: + actions.append(LevelAction.constant(None)) + elif s.act(p) == q: + rep = self.representative[q] + loop = self.from_rep[p] * s * self.to_rep[q] + level = self.levels[i] + k = level.representatives.index(rep) + index = self._level_index[i] + mapping = {index[(rep, tile)]: index[(rep, loop.act(tile))] for tile in level.tiles[k]} + actions.append(LevelAction.permutation(mapping)) + else: + actions.append(LevelAction.constant(new_coords[i])) + return tuple(actions) + + def _frames(self, chain: Sequence[Subset]) -> list[Subset]: + """``frames[i]`` is the deepest set of ``chain`` placed above level ``i``.""" + frames = [] + current = self._full + members = iter(chain[1:]) + nxt = next(members, None) + for i in range(self.n_levels): + frames.append(current) + if len(current) > 1 and self.depth(current) == i + 1 and nxt is not None: + current = nxt + nxt = next(members, None) + return frames + + # -- outputs -------------------------------------------------------------- + + def to_cascade(self, *, max_states: int = DEFAULT_MAX_SIZE) -> CascadeProduct: + """Return the holonomy cascade product with explicit dependency functions.""" + from sofic.algebra.cascade import build_holonomy_cascade + + return build_holonomy_cascade(self, max_states=max_states) + + def effective_theory(self, k: int, *, machine: Any = None, max_states: int = DEFAULT_MAX_SIZE) -> EffectiveTheory: + """Return the automaton on the top ``k`` levels of the holonomy cascade. + + Truncating the cascade gives a coarse-grained machine whose states + are *superstates*: sets of original states that agree on the top + ``k`` coordinates of some lift. Superstates may overlap. This is the + "effective theory" of :cite:`DeDeo2011`. + + With a stochastic ``machine`` (whose transition semigroup is this + decomposition's), states are grouped by the prefix of their canonical + lift. The result reports the stationary mass of each group and + whether that partition is strongly lumpable (see + :func:`~sofic.generators.lumping.is_lumpable`), that is, whether the + coarse-grained process is itself Markov on the superstates. + """ + full = self.to_cascade(max_states=max_states) + cascade = full.truncate(k) + automaton = cascade.as_dfa() + projection = {coords: coords[:k] for coords in full.decoding} + distribution = lumpable = None + if machine is not None: + from sofic.generators.lumping import is_lumpable + + index = machine.reindex() + pi = machine.stationary_distribution() + partition = {index.state(i): self.coordinates(index.state(i))[:k] for i in range(len(index))} + distribution = {} + for i in range(len(index)): + key = partition[index.state(i)] + distribution[key] = distribution.get(key, 0.0) + float(pi[i]) + lumpable = is_lumpable(machine, partition) + superstates = {prefix: cast("frozenset[Hashable]", states) for prefix, states in cascade.decoding.items()} + return EffectiveTheory(automaton, superstates, projection, distribution, lumpable) + + def describe(self) -> str: + """Human-readable summary: one line per level with representatives and groups.""" + labels = self.semigroup.states + + def show(subset: Subset) -> str: + return "{" + ", ".join(str(labels[x]) for x in sorted(subset)) + "}" + + lines = [f"holonomy decomposition: {len(self.image_set)} image sets, {self.n_levels} levels"] + for level in self.levels: + parts = [] + for rep, tiles, group in zip(level.representatives, level.tiles, level.groups, strict=True): + parts.append(f"{show(rep)} [{len(tiles)} tiles, |H|={group.order()}]") + lines.append(f" depth {level.depth}: " + "; ".join(parts)) + return "\n".join(lines) + + +def holonomy_decomposition(machine: Any, **kwargs: Any) -> HolonomyDecomposition: + """Return the holonomy decomposition of a semigroup or deterministic machine. + + ``machine`` is anything :func:`~sofic.algebra.adapters.transformation_semigroup` + accepts; keyword arguments are passed through to it. + + Examples + -------- + >>> from sofic.algebra import TransformationSemigroup + >>> hd = holonomy_decomposition(TransformationSemigroup({"t": (1, 0, 2), "c": (0, 1, 1)})) + >>> [level.order for level in hd.levels] + [1, 2] + """ + from sofic.algebra.adapters import transformation_semigroup + + return HolonomyDecomposition(transformation_semigroup(machine, **kwargs)) diff --git a/sofic/algebra/information.py b/sofic/algebra/information.py new file mode 100644 index 0000000..1a22446 --- /dev/null +++ b/sofic/algebra/information.py @@ -0,0 +1,75 @@ +"""Information carried by each level of a holonomy decomposition. + +.. note:: Original, uncited: the literature does not pair holonomy + coordinates with a stationary measure. +""" + +from __future__ import annotations + +from typing import Any, NamedTuple + +__all__ = ["HolonomyInformation", "holonomy_information"] + + +class HolonomyInformation(NamedTuple): + """Per-level split of the stationary state entropy. + + ``distribution`` is the joint law of the level coordinates + ``(L_1, ..., L_h)``, top level first, with ``"*"`` for an unused level. + ``level_entropies[i]`` is ``H[L_{i+1} | L_1, ..., L_i]``; they sum to + ``total``, the entropy of the state distribution. + """ + + distribution: Any + level_entropies: tuple[float, ...] + total: float + + +def holonomy_information(machine: Any) -> HolonomyInformation: + """Split the stationary state entropy across holonomy levels. + + Each state of a unifilar stochastic generator is lifted to its canonical + maximal tile chain (:meth:`~sofic.algebra.holonomy.HolonomyDecomposition.lift`) + and so to cascade coordinates ``(L_1, ..., L_h)``. Distinct states get + distinct coordinates, so ``H[L_1, ..., L_h] = H[S]`` and the chain rule + gives ``H[S] = sum_i H[L_i | L_>> from sofic.examples import even_process + >>> info = holonomy_information(even_process()) + >>> round(info.total, 6) == round(even_process().statistical_complexity(), 6) + True + """ + import dit + from dit.shannon import conditional_entropy, entropy + + from sofic.algebra.adapters import transformation_semigroup + from sofic.algebra.holonomy import HolonomyDecomposition + + hd = HolonomyDecomposition(transformation_semigroup(machine)) + index = machine.reindex() + pi = machine.stationary_distribution() + outcomes = [] + pmf = [] + for i in range(len(index)): + p = float(pi[i]) + if p <= 0: + continue + coords = hd.coordinates(index.state(i)) + outcomes.append(tuple("*" if c is None else str(c) for c in coords)) + pmf.append(p) + if hd.n_levels == 0: + return HolonomyInformation(None, (), 0.0) + dist = dit.Distribution(outcomes, pmf) + entropies = tuple(float(conditional_entropy(dist, [i], list(range(i)))) for i in range(hd.n_levels)) + return HolonomyInformation(dist, entropies, float(entropy(dist))) diff --git a/sofic/algebra/random_maps.py b/sofic/algebra/random_maps.py new file mode 100644 index 0000000..27586e5 --- /dev/null +++ b/sofic/algebra/random_maps.py @@ -0,0 +1,390 @@ +"""Random-map representations of probabilistic automata. + +A probabilistic automaton assigns each input letter ``s`` a row-stochastic +matrix ``P_s``. It can be rewritten as a cascade: a Bernoulli process picks a +deterministic map ``m`` of the states with probability ``pi(m | s)``, and a +deterministic machine applies it, so that ``P_s[q, q'] = sum_{m(q) = q'} +pi(m | s)`` :cite:`Maler1995`. The deterministic part is an ordinary +transformation semigroup and has a Krohn-Rhodes decomposition. Different +choices of ``pi`` are different *failure models* for the mechanism +:cite:`DeDeo2011`: + +* :func:`independent_random_maps` -- each state fails independently, + ``pi(m | s) = prod_q P_s[q, m(q)]`` (the ε-independent model, Eq. 2 of + :cite:`DeDeo2011`). +* :func:`coupled_random_maps` -- one shared uniform draw decides every + state's move (the random mapping representation of + :cite:`LevinPeresWilmer2009`), a model of correlated failures. +* :func:`random_maps` -- any decomposition given explicitly. + +:func:`element_walk` follows the distribution of the composed random map +over time, and :func:`walk_information` splits the information the state +loses about its past into an irreversible and a dissipative part. +""" + +from __future__ import annotations + +import itertools +from collections.abc import Hashable, Iterable, Mapping, Sequence +from math import prod +from typing import Any, NamedTuple + +import numpy as np + +from sofic.algebra.transformation import DEFAULT_MAX_SIZE, Transformation, TransformationSemigroup + +__all__ = [ + "RandomMapRepresentation", + "WalkInformation", + "classify_elements", + "coupled_random_maps", + "element_walk", + "independent_random_maps", + "random_maps", + "walk_information", +] + +#: Letter used when a source has a single input (the passage of time). +TIME = "t" + + +class RandomMapRepresentation: + """A distribution over deterministic maps of the states for each input letter. + + Parameters + ---------- + weights : mapping + ``weights[letter][m]`` is the probability of applying map ``m`` on + ``letter``; ``m`` is a :class:`~sofic.algebra.transformation.Transformation` + or a tuple of images. Zero weights are dropped. + states : sequence, optional + Labels for the points ``0, ..., n - 1``. + """ + + def __init__( + self, + weights: Mapping[Hashable, Mapping[Any, float]], + *, + states: Sequence[Hashable] | None = None, + atol: float = 1e-9, + ) -> None: + self._weights: dict[Hashable, dict[Transformation, float]] = {} + degree = None + for letter, table in weights.items(): + maps: dict[Transformation, float] = {} + for m, w in table.items(): + t = m if isinstance(m, Transformation) else Transformation(tuple(m)) + if w < -atol: + raise ValueError(f"negative weight on {letter!r}") + if w > atol: + maps[t] = maps.get(t, 0.0) + float(w) + degree = t.degree if degree is None else degree + if t.degree != degree: + raise ValueError("all maps must act on the same number of states") + total = sum(maps.values()) + if abs(total - 1.0) > 1e-6: + raise ValueError(f"weights on {letter!r} sum to {total}, not 1") + self._weights[letter] = maps + if degree is None: + raise ValueError("need at least one letter with one map") + self.degree = degree + self.states = tuple(states) if states is not None else tuple(range(degree)) + self._semigroups: dict[int, TransformationSemigroup] = {} + + @property + def letters(self) -> tuple[Hashable, ...]: + return tuple(self._weights) + + def maps(self, letter: Hashable) -> dict[Transformation, float]: + return dict(self._weights[letter]) + + def matrices(self) -> dict[Hashable, np.ndarray]: + """Rebuild ``P_s[q, q'] = sum_{m(q) = q'} pi(m | s)`` for every letter.""" + result = {} + for letter, table in self._weights.items(): + matrix = np.zeros((self.degree, self.degree)) + for m, w in table.items(): + matrix[np.arange(self.degree), m.images] += w + result[letter] = matrix + return result + + def generator_weights(self) -> dict[Hashable, dict[Hashable, float]]: + """``{letter: {generator label: probability}}`` for :meth:`semigroup`'s generators.""" + return { + letter: {(letter, k): w for k, w in enumerate(table.values())} for letter, table in self._weights.items() + } + + def semigroup(self, *, max_size: int = DEFAULT_MAX_SIZE) -> TransformationSemigroup: + """The monoid generated by every map in the support, labeled ``(letter, k)``. + + This is the deterministic machine beneath the Bernoulli process; its + holonomy decomposition is the Krohn-Rhodes decomposition of the noisy + automaton :cite:`DeDeo2011`. + """ + if max_size not in self._semigroups: + gens: dict[Hashable, Transformation] = { + (letter, k): m for letter, table in self._weights.items() for k, m in enumerate(table) + } + self._semigroups[max_size] = TransformationSemigroup(gens, states=self.states, max_size=max_size) + return self._semigroups[max_size] + + def __repr__(self) -> str: + sizes = {letter: len(table) for letter, table in self._weights.items()} + return f"RandomMapRepresentation(degree={self.degree}, maps={sizes})" + + +def _as_matrices(source: Any) -> tuple[dict[Hashable, np.ndarray], tuple[Hashable, ...] | None]: + if isinstance(source, Mapping): + return {letter: np.asarray(m, dtype=float) for letter, m in source.items()}, None + if isinstance(source, (np.ndarray, list, tuple)): + return {TIME: np.asarray(source, dtype=float)}, None + if hasattr(source, "reindex") and hasattr(source, "transitions"): + from sofic.graph import ATTR_PROB + from sofic.properties import transition_matrix + + index = source.reindex() + matrix, states = transition_matrix(source, attr=ATTR_PROB, states=index.states) + return {TIME: np.asarray(matrix, dtype=float)}, tuple(states) + raise TypeError(f"cannot read stochastic matrices from {type(source).__name__}") + + +def _check_stochastic(matrices: Mapping[Hashable, np.ndarray]) -> int: + degrees = {m.shape for m in matrices.values()} + if len(degrees) != 1: + raise ValueError("all matrices must have the same shape") + (shape,) = degrees + if len(shape) != 2 or shape[0] != shape[1]: + raise ValueError("matrices must be square") + for letter, m in matrices.items(): + if (m < -1e-12).any() or not np.allclose(m.sum(axis=1), 1.0): + raise ValueError(f"matrix for {letter!r} is not row-stochastic") + return shape[0] + + +def independent_random_maps( + source: Any, + *, + states: Sequence[Hashable] | None = None, + max_maps: int = 100_000, +) -> RandomMapRepresentation: + """Each state picks its successor independently: ``pi(m | s) = prod_q P_s[q, m(q)]``. + + This is Eq. 2 of :cite:`DeDeo2011` (the ε-independent failure model) + and the decomposition of :cite:`Maler1995`. The support is every map + compatible with the nonzero entries, so its size is the product of the + row support sizes; ``max_maps`` caps it. + + ``source`` is a mapping from letter to a row-stochastic matrix, a single + matrix (one letter, ``"t"``), or a Markov chain / HMM, whose state + transition matrix is used as the single letter. + """ + matrices, model_states = _as_matrices(source) + n = _check_stochastic(matrices) + weights: dict[Hashable, dict[Transformation, float]] = {} + for letter, matrix in matrices.items(): + supports = [np.flatnonzero(matrix[q] > 0) for q in range(n)] + count = prod(len(s) for s in supports) + if count > max_maps: + raise ValueError(f"{count} maps on {letter!r} exceed max_maps={max_maps}") + weights[letter] = { + Transformation(tuple(int(x) for x in choice)): float(prod(matrix[q, x] for q, x in enumerate(choice))) + for choice in itertools.product(*supports) + } + return RandomMapRepresentation(weights, states=states or model_states) + + +def coupled_random_maps( + source: Any, + *, + states: Sequence[Hashable] | None = None, + decimals: int = 12, +) -> RandomMapRepresentation: + """All states share one uniform draw ``U``: ``m_U(q)`` is the ``U``-quantile of row ``q``. + + Each row's outcomes are ordered by decreasing probability (ties by + index) and ``m_U(q)`` is the outcome whose cumulative interval contains + ``U``. The breakpoints of all rows cut ``[0, 1)`` into at most + ``1 + sum_q (|supp_q| - 1)`` intervals, each giving one map weighted by + its length. This is the random mapping representation of + :cite:`LevinPeresWilmer2009` and models perfectly correlated failures: + for the noisy XOR register it gives exactly two maps, XOR and NOT-XOR + :cite:`DeDeo2011`. ``source`` is read as in :func:`independent_random_maps`. + """ + matrices, model_states = _as_matrices(source) + n = _check_stochastic(matrices) + weights: dict[Hashable, dict[Transformation, float]] = {} + for letter, matrix in matrices.items(): + orders = [] + cuts = {0.0, 1.0} + for q in range(n): + support = sorted(np.flatnonzero(matrix[q] > 0), key=lambda x, q=q: (-matrix[q, x], x)) + cumulative = np.round(np.cumsum([matrix[q, x] for x in support]), decimals) + orders.append((support, cumulative)) + cuts.update(float(c) for c in cumulative[:-1]) + bounds = sorted(cuts) + table: dict[Transformation, float] = {} + for lo, hi in itertools.pairwise(bounds): + if hi - lo <= 0: + continue + u = (lo + hi) / 2 + images = tuple( + int(support[min(int(np.searchsorted(cumulative, u, side="right")), len(support) - 1)]) + for support, cumulative in orders + ) + m = Transformation(images) + table[m] = table.get(m, 0.0) + (hi - lo) + weights[letter] = table + return RandomMapRepresentation(weights, states=states or model_states) + + +def random_maps( + weights: Mapping[Hashable, Mapping[Any, float]], + *, + matrices: Any = None, + states: Sequence[Hashable] | None = None, +) -> RandomMapRepresentation: + """An explicit decomposition, optionally checked against ``matrices``. + + Raises ``ValueError`` if the maps do not reproduce the given stochastic + matrices (read as in :func:`independent_random_maps`). + """ + rep = RandomMapRepresentation(weights, states=states) + if matrices is not None: + expected, _ = _as_matrices(matrices) + rebuilt = rep.matrices() + if set(expected) != set(rebuilt) or any(not np.allclose(expected[s], rebuilt[s]) for s in expected): + raise ValueError("the maps do not reproduce the given matrices") + return rep + + +def element_walk( + rep: RandomMapRepresentation, + steps: int | None = None, + *, + inputs: Sequence[Hashable] | None = None, + max_size: int = DEFAULT_MAX_SIZE, +) -> np.ndarray: + """Distribution of the composed random map ``M_t = m_1 m_2 ... m_t`` over semigroup elements. + + Row ``t`` of the result is indexed like :meth:`RandomMapRepresentation.semigroup` + (index ``0`` is the identity, ``M_0``). Each step multiplies on the right + by a random generator, a Markov chain on the right Cayley graph + :cite:`DeDeo2011`. With several letters pass ``inputs``, the letter read + at each step; ``steps`` then defaults to ``len(inputs)``. + """ + sg = rep.semigroup(max_size=max_size) + if inputs is None: + if len(rep.letters) != 1: + raise ValueError("pass inputs when there is more than one letter") + if steps is None: + raise ValueError("pass steps or inputs") + inputs = [rep.letters[0]] * steps + inputs = list(inputs) + steps = len(inputs) if steps is None else steps + if steps > len(inputs): + raise ValueError("fewer inputs than steps") + right = np.asarray(sg.right_cayley()) + labels = sg.generator_labels + weights = rep.generator_weights() + column = {label: k for k, label in enumerate(labels)} + dist = np.zeros((steps + 1, len(sg))) + dist[0, 0] = 1.0 + for t in range(steps): + current = dist[t] + nxt = dist[t + 1] + for label, w in weights[inputs[t]].items(): + np.add.at(nxt, right[:, column[label]], w * current) + return dist + + +def classify_elements(semigroup: TransformationSemigroup) -> dict[str, frozenset[int]]: + """Index sets of the identity, units (permutations), resets (constants), and minimal ideal.""" + units, resets, identity = set(), set(), set() + for i, element in enumerate(semigroup): + if element.is_permutation(): + units.add(i) + if element == Transformation.identity(semigroup.degree): + identity.add(i) + if element.is_constant(): + resets.add(i) + minimal = frozenset(semigroup.green().minimal_ideal().elements()) + return { + "identity": frozenset(identity), + "units": frozenset(units), + "resets": frozenset(resets), + "minimal_ideal": minimal, + } + + +class WalkInformation(NamedTuple): + """Information about the initial state retained by the noisy machine, per step. + + All fields are arrays of length ``steps + 1`` in bits. ``element_entropy`` + is ``H[M_t]``; ``mutual`` is ``I[X_0 : X_t]``; ``conditional`` is + ``I[X_0 : X_t | M_t]``. ``irreversibility = H[X_0] - conditional`` is + information erased by maps that merge states, even if the random choices + are known; ``dissipation = conditional - mutual`` is information lost only + because the choices are unknown. Both are nonnegative. + """ + + element_entropy: np.ndarray + mutual: np.ndarray + conditional: np.ndarray + irreversibility: np.ndarray + dissipation: np.ndarray + + +def _entropy(p: np.ndarray) -> float: + p = p[p > 0] + return float(-(p * np.log2(p)).sum()) + + +def walk_information( + rep: RandomMapRepresentation, + steps: int | None = None, + *, + initial: Iterable[float] | None = None, + inputs: Sequence[Hashable] | None = None, + max_size: int = DEFAULT_MAX_SIZE, +) -> WalkInformation: + """Split the information ``X_t`` keeps about ``X_0`` into irreversible and dissipative losses. + + ``X_0`` has law ``initial`` (uniform by default) and is independent of + the random maps, and ``X_t = X_0 · M_t``. Given ``M_t = m`` the state + is a deterministic function of ``X_0``, so + ``I[X_0 : X_t | M_t] = sum_m P(m) H[m(X_0)]``, and independence gives + ``I[X_0 : X_t | M_t] >= I[X_0 : X_t]``. DeDeo distinguishes + *irreversibility* (a map that deletes information) from *dissipation* + (histories merged by ignorance of which random operation ran) + :cite:`DeDeo2011`; this function measures each. + + .. note:: Original, uncited: the information-theoretic split is not in + :cite:`DeDeo2011`. + """ + sg = rep.semigroup(max_size=max_size) + walk = element_walk(rep, steps, inputs=inputs, max_size=max_size) + n = sg.degree + p0 = np.full(n, 1.0 / n) if initial is None else np.asarray(list(initial), dtype=float) + images = np.asarray([e.images for e in sg]) + pushed = np.zeros((len(sg), n)) + for i in range(len(sg)): + np.add.at(pushed[i], images[i], p0) + branch_entropy = np.asarray([_entropy(row) for row in pushed]) + h0 = _entropy(p0) + element_entropy, mutual, conditional = [], [], [] + for q in walk: + element_entropy.append(_entropy(q)) + conditional.append(float(q @ branch_entropy)) + joint = np.zeros((n, n)) + for i in np.flatnonzero(q > 0): + joint[np.arange(n), images[i]] += q[i] * p0 + mutual.append(h0 + _entropy(joint.sum(axis=0)) - _entropy(joint.ravel())) + conditional_arr = np.asarray(conditional) + mutual_arr = np.asarray(mutual) + return WalkInformation( + np.asarray(element_entropy), + mutual_arr, + conditional_arr, + h0 - conditional_arr, + conditional_arr - mutual_arr, + ) diff --git a/sofic/algebra/syntactic.py b/sofic/algebra/syntactic.py new file mode 100644 index 0000000..c10b3cc --- /dev/null +++ b/sofic/algebra/syntactic.py @@ -0,0 +1,148 @@ +"""Syntactic monoids and star-freeness of regular languages. + +The syntactic monoid of a regular language ``L`` is the quotient of ``A*`` +by the syntactic congruence ``u ~ v iff xuy ∈ L ⇔ xvy ∈ L`` for all ``x, y``; +it is isomorphic to the transition monoid of the minimal complete DFA of +``L`` :cite:`Pin1986`. It is the smallest monoid recognizing ``L``, so its +properties are properties of the language rather than of a presentation. +""" + +from __future__ import annotations + +from collections.abc import Hashable, Iterable +from typing import Any + +from sofic.algebra.adapters import semiautomaton +from sofic.algebra.transformation import DEFAULT_MAX_SIZE, SINK, TransformationSemigroup + +__all__ = ["factor_monoid", "is_star_free", "syntactic_monoid", "transition_monoid"] + + +def transition_monoid( + automaton: Any, + *, + alphabet: Iterable[Hashable] | None = None, + max_size: int = DEFAULT_MAX_SIZE, + monoid: bool = True, +) -> TransformationSemigroup: + """Return the transition monoid of a deterministic automaton. + + Missing transitions go to an absorbing ``SINK``; symbols of ``alphabet`` + that never appear map every state to the sink. With ``monoid=False`` the + result is the transition *semigroup* of nonempty words. + """ + states, maps = semiautomaton(automaton, alphabet=alphabet) + if not states: + states = [SINK] + sink_maps: dict[Hashable, dict[Hashable, Hashable]] = {symbol: {SINK: SINK} for symbol in maps} + maps = sink_maps + if not maps: + return TransformationSemigroup({}, states=states, monoid=True, max_size=max_size) + return TransformationSemigroup.from_maps(maps, states=states, monoid=monoid, max_size=max_size) + + +def syntactic_monoid( + language: Any, + *, + alphabet: Iterable[Hashable] | None = None, + max_size: int = DEFAULT_MAX_SIZE, + monoid: bool = True, +) -> TransformationSemigroup: + """Return the syntactic monoid of a regular language. + + ``language`` may be a labeled automaton (DFA, NFA, ...), an + :class:`~sofic.automata.languages.base.AutomatonLanguage`, or a regular + expression string accepted by :func:`~sofic.automata.regex.regex_to_nfa`. + The alphabet defaults to the automaton's input alphabet; it matters, + because a symbol that never occurs acts as a zero. With ``monoid=False`` + the result is the syntactic *semigroup* (images of nonempty words), which + several variety tests need :cite:`Pin1986`. + + Examples + -------- + >>> m = syntactic_monoid("(aa)*", alphabet="a") + >>> len(m), m.is_aperiodic() + (2, False) + """ + from sofic.automata.algorithms import _to_nfa, _transition_alphabet, minimize + from sofic.automata.base import LabeledAutomaton + from sofic.automata.dfa import DFA + from sofic.automata.languages.base import AutomatonLanguage + from sofic.automata.nfa import NFA + + if isinstance(language, str): + from sofic.automata.regex import regex_to_nfa + + language = regex_to_nfa(language, alphabet=alphabet) + if isinstance(language, AutomatonLanguage): + language = language.automaton + if not isinstance(language, LabeledAutomaton): + raise TypeError(f"cannot build a syntactic monoid from {type(language).__name__}") + symbols = frozenset(alphabet) if alphabet is not None else _transition_alphabet(language) + automaton = language if isinstance(language, (DFA, NFA)) else _to_nfa(language) + minimal = minimize(automaton, alphabet=symbols) + return transition_monoid(minimal, alphabet=_ordered(symbols), max_size=max_size, monoid=monoid) + + +def factor_monoid( + machine: Any, + *, + alphabet: Iterable[Hashable] | None = None, + max_size: int = DEFAULT_MAX_SIZE, + monoid: bool = True, +) -> TransformationSemigroup: + """Return the syntactic monoid of the words labeling paths in ``machine``. + + Every state is initial and accepting, so for a shift presentation this + is the factor language ``B(X)`` and for a stochastic generator the + support of its word distribution from any start state. Edges may carry + ``symbol`` or ``emission`` labels; zero-probability edges are ignored, and + the machine need not be unifilar. + """ + from sofic.automata.nfa import NFA + from sofic.graph import ATTR_EMISSION_EDGE, ATTR_PROB, ATTR_SYMBOL + + states = frozenset(machine.states()) + nfa = NFA(initial_states=states, accepting_states=states) + for state in states: + nfa.graph.add_state(state) + symbols: set[Hashable] = set() + for transition in machine.transitions(): + label = transition.data.get(ATTR_SYMBOL, transition.data.get(ATTR_EMISSION_EDGE)) + prob = transition.data.get(ATTR_PROB) + if label is None or (prob is not None and prob == 0): + continue + nfa.add_transition(transition.source, transition.target, label) + symbols.add(label) + alphabet = frozenset(alphabet) if alphabet is not None else frozenset(symbols) + nfa.input_alphabet = alphabet | symbols + return syntactic_monoid(nfa, alphabet=alphabet, max_size=max_size, monoid=monoid) + + +def is_star_free( + language: Any, + *, + alphabet: Iterable[Hashable] | None = None, + max_size: int = DEFAULT_MAX_SIZE, +) -> bool: + """Return whether a regular language is star-free. + + A regular language is star-free (expressible with union, concatenation + and complement but no Kleene star) iff its syntactic monoid is aperiodic + :cite:`Schutzenberger1965`, iff it is definable in first-order logic + :cite:`McNaughtonPapert1971`, iff its Krohn-Rhodes complexity is zero. + + Examples + -------- + >>> is_star_free("(ab)*", alphabet="ab"), is_star_free("(aa)*", alphabet="a") + (True, False) + """ + return syntactic_monoid(language, alphabet=alphabet, max_size=max_size).is_aperiodic() + + +def _ordered(symbols: Iterable[Hashable]) -> list[Hashable]: + values = list(symbols) + try: + return sorted(values) # type: ignore[type-var] + except TypeError: + return sorted(values, key=repr) diff --git a/sofic/algebra/transformation.py b/sofic/algebra/transformation.py new file mode 100644 index 0000000..cc39534 --- /dev/null +++ b/sofic/algebra/transformation.py @@ -0,0 +1,511 @@ +"""Finite transformation semigroups. + +A transformation semigroup ``(Q, S)`` is a finite state set ``Q`` together +with a semigroup ``S`` of total maps ``Q -> Q`` closed under composition. Maps +act on the right, ``x · (st) = (x · s) · t``, so a word read left to right is +the product of its letters in the same order :cite:`Eilenberg1976`. A +deterministic semiautomaton (a DFA without initial or accepting states) is the +same data with a chosen generating set: one transformation per input symbol. +""" + +from __future__ import annotations + +from collections import deque +from collections.abc import Hashable, Iterable, Iterator, Mapping, Sequence +from dataclasses import dataclass +from math import lcm +from typing import TYPE_CHECKING, Any, Literal + +from sofic.exceptions import SemigroupTooLargeError + +if TYPE_CHECKING: + import networkx as nx + + from sofic.algebra.complexity import ComplexityBounds + from sofic.algebra.green import GreenStructure + from sofic.algebra.holonomy import HolonomyDecomposition + +DEFAULT_MAX_SIZE = 100_000 + +__all__ = [ + "DEFAULT_MAX_SIZE", + "SINK", + "Transformation", + "TransformationSemigroup", +] + + +class _Sink: + """Absorbing state added when completing partial transformations.""" + + _instance: _Sink | None = None + + def __new__(cls) -> _Sink: + if cls._instance is None: + cls._instance = super().__new__(cls) + return cls._instance + + def __repr__(self) -> str: + return "SINK" + + def __reduce__(self) -> tuple[Any, ...]: + return (_Sink, ()) + + +SINK = _Sink() + + +@dataclass(frozen=True, slots=True) +class Transformation: + """A total map on ``{0, ..., n - 1}`` acting on the right. + + ``x · s = images[x]``, and ``s * t`` is "``s`` then ``t``", so + ``x · (s * t) = (x · s) · t``. + + Examples + -------- + >>> s = Transformation((1, 0, 2)) + >>> t = Transformation((2, 2, 2)) + >>> (s * t).images + (2, 2, 2) + >>> s.index_period() + (1, 2) + """ + + images: tuple[int, ...] + + @classmethod + def identity(cls, degree: int) -> Transformation: + return cls(tuple(range(degree))) + + @classmethod + def constant(cls, degree: int, value: int) -> Transformation: + return cls((value,) * degree) + + @property + def degree(self) -> int: + return len(self.images) + + def __call__(self, x: int) -> int: + return self.images[x] + + def __mul__(self, other: object) -> Transformation: + if not isinstance(other, Transformation): + return NotImplemented + return Transformation(_compose(self.images, other.images)) + + def act(self, subset: Iterable[int]) -> frozenset[int]: + """Return the setwise image ``P · s``.""" + return frozenset(self.images[x] for x in subset) + + def image(self) -> frozenset[int]: + return frozenset(self.images) + + def rank(self) -> int: + return len(set(self.images)) + + def kernel(self) -> tuple[frozenset[int], ...]: + """Return the kernel partition: blocks of states with a common image.""" + blocks: dict[int, set[int]] = {} + for x, y in enumerate(self.images): + blocks.setdefault(y, set()).add(x) + return tuple(sorted((frozenset(b) for b in blocks.values()), key=sorted)) + + def is_permutation(self) -> bool: + return self.rank() == self.degree + + def is_idempotent(self) -> bool: + return all(self.images[y] == y for y in self.images) + + def is_constant(self) -> bool: + return self.rank() <= 1 + + def power(self, exponent: int) -> Transformation: + """Return ``s**exponent`` for ``exponent >= 0`` (``s**0`` is the identity).""" + if exponent < 0: + raise ValueError("exponent must be nonnegative") + result = tuple(range(self.degree)) + base = self.images + while exponent: + if exponent & 1: + result = _compose(result, base) + base = _compose(base, base) + exponent >>= 1 + return Transformation(result) + + def __pow__(self, exponent: int) -> Transformation: + return self.power(exponent) + + def index_period(self) -> tuple[int, int]: + """Return ``(m, r)``, the least pair with ``s**(m + r) == s**m``. + + The element generates an aperiodic cyclic semigroup exactly when + ``r == 1``. + """ + seen: dict[tuple[int, ...], int] = {} + current = self.images + k = 1 + while current not in seen: + seen[current] = k + current = _compose(current, self.images) + k += 1 + index = seen[current] + return index, k - index + + def restrict(self, subset: Sequence[int]) -> tuple[int, ...]: + """Return the action on the ordered ``subset`` as positions in it. + + Requires ``subset · s`` to be contained in ``subset``. + """ + position = {x: i for i, x in enumerate(subset)} + return tuple(position[self.images[x]] for x in subset) + + +def _compose(left: tuple[int, ...], right: tuple[int, ...]) -> tuple[int, ...]: + """Return ``left`` then ``right`` as raw image tuples.""" + return tuple(right[y] for y in left) + + +def _permutation_order(images: tuple[int, ...], support: Iterable[int]) -> int: + """Order of ``images`` restricted to ``support`` (which it must permute).""" + order = 1 + seen: set[int] = set() + for start in support: + if start in seen: + continue + length = 0 + x = start + while x not in seen: + seen.add(x) + x = images[x] + length += 1 + order = lcm(order, length) + return order + + +class TransformationSemigroup: + """The transformation semigroup generated by labeled maps of a finite set. + + Elements are enumerated lazily by breadth-first closure the first time + they are needed, so every element carries a shortest generating word. + Constructions that only need the generators -- notably the holonomy + decomposition -- never trigger the enumeration. + + Parameters + ---------- + generators : mapping + Generator label (typically an input symbol) to a :class:`Transformation` + or a sequence of images on ``range(degree)``. + states : sequence, optional + Labels for the points ``0, ..., degree - 1``. Defaults to the integers. + monoid : bool + Adjoin the identity (the empty word). This is the convention for the + transition monoid of an automaton. + max_size : int + Raise :class:`~sofic.exceptions.SemigroupTooLargeError` if the closure + exceeds this many elements. + + Examples + -------- + >>> flip = TransformationSemigroup({"a": (1, 0)}) + >>> len(flip), flip.is_aperiodic() + (2, False) + """ + + def __init__( + self, + generators: Mapping[Hashable, Transformation | Sequence[int]], + *, + states: Sequence[Hashable] | None = None, + monoid: bool = True, + max_size: int = DEFAULT_MAX_SIZE, + ) -> None: + gens = { + label: g if isinstance(g, Transformation) else Transformation(tuple(int(x) for x in g)) + for label, g in generators.items() + } + if states is not None: + degree = len(states) + elif gens: + degree = next(iter(gens.values())).degree + else: + raise ValueError("need at least one generator or an explicit state list") + if degree == 0: + raise ValueError("a transformation semigroup needs at least one state") + for label, g in gens.items(): + if g.degree != degree or any(not 0 <= y < degree for y in g.images): + raise ValueError(f"generator {label!r} is not a map on {degree} states") + if not gens and not monoid: + raise ValueError("a semigroup with no generators must be a monoid") + self._generators = gens + self._labels = tuple(gens) + self._states = tuple(states) if states is not None else tuple(range(degree)) + if len(set(self._states)) != degree: + raise ValueError("state labels must be distinct") + self._degree = degree + self._monoid = monoid + self._max_size = max_size + self._elements: list[tuple[int, ...]] | None = None + self._words: list[tuple[Hashable, ...]] = [] + self._index: dict[tuple[int, ...], int] = {} + self._right: list[list[int]] = [] + self._left: list[list[int]] | None = None + + # -- construction helpers ------------------------------------------------- + + @classmethod + def from_maps( + cls, + maps: Mapping[Hashable, Mapping[Hashable, Hashable]], + *, + states: Iterable[Hashable] | None = None, + complete: Literal["sink", "error"] = "sink", + monoid: bool = True, + max_size: int = DEFAULT_MAX_SIZE, + ) -> TransformationSemigroup: + """Build from per-label state maps, completing partial maps. + + ``maps[label][state]`` is the successor of ``state`` on ``label``; a + missing entry is undefined. With ``complete="sink"`` undefined moves go + to an absorbing :data:`SINK` state, added only when needed; with + ``complete="error"`` they raise ``ValueError``. + """ + if states is None: + seen: dict[Hashable, None] = {} + for table in maps.values(): + for source, target in table.items(): + seen.setdefault(source) + seen.setdefault(target) + ordered = _sorted_states(seen) + else: + ordered = list(dict.fromkeys(states)) + position = {s: i for i, s in enumerate(ordered)} + partial = any(s not in table for table in maps.values() for s in ordered) + if partial: + if complete == "error": + raise ValueError("partial transformation; pass complete='sink' to add an absorbing state") + if SINK not in position: + position[SINK] = len(ordered) + ordered.append(SINK) + generators: dict[Hashable, tuple[int, ...]] = {} + for label, table in maps.items(): + images: list[int] = [] + for state in ordered: + target = table.get(state, SINK) if state is not SINK else SINK + if target not in position: + raise ValueError(f"unknown target state {target!r} on {label!r}") + images.append(position[target]) + generators[label] = tuple(images) + return cls(generators, states=ordered, monoid=monoid, max_size=max_size) + + # -- basic data ----------------------------------------------------------- + + @property + def degree(self) -> int: + return self._degree + + @property + def states(self) -> tuple[Hashable, ...]: + return self._states + + @property + def generators(self) -> dict[Hashable, Transformation]: + return dict(self._generators) + + @property + def generator_labels(self) -> tuple[Hashable, ...]: + return self._labels + + @property + def is_monoid(self) -> bool: + return self._monoid + + def state_index(self, state: Hashable) -> int: + return self._states.index(state) + + # -- enumeration ---------------------------------------------------------- + + def _closure(self) -> list[tuple[int, ...]]: + if self._elements is not None: + return self._elements + try: + return self._enumerate() + except SemigroupTooLargeError: + self._index, self._words, self._right = {}, [], [] + raise + + def _enumerate(self) -> list[tuple[int, ...]]: + elements: list[tuple[int, ...]] = [] + gens = [self._generators[label].images for label in self._labels] + + def add(element: tuple[int, ...], word: tuple[Hashable, ...]) -> int: + index = self._index.get(element) + if index is not None: + return index + if len(elements) >= self._max_size: + raise SemigroupTooLargeError( + f"semigroup exceeds max_size={self._max_size} elements; raise the cap if intended" + ) + self._index[element] = len(elements) + elements.append(element) + self._words.append(word) + self._right.append([]) + return len(elements) - 1 + + if self._monoid: + add(tuple(range(self._degree)), ()) + for label, g in zip(self._labels, gens, strict=True): + add(g, (label,)) + queue = deque(range(len(elements))) + while queue: + i = queue.popleft() + row = [] + for label, g in zip(self._labels, gens, strict=True): + before = len(elements) + j = add(_compose(elements[i], g), (*self._words[i], label)) + if j >= before: + queue.append(j) + row.append(j) + self._right[i] = row + self._elements = elements + return elements + + def __len__(self) -> int: + return len(self._closure()) + + def __iter__(self) -> Iterator[Transformation]: + return (Transformation(e) for e in self._closure()) + + def __getitem__(self, index: int) -> Transformation: + return Transformation(self._closure()[index]) + + def __contains__(self, item: object) -> bool: + if not isinstance(item, Transformation): + return False + self._closure() + return item.images in self._index + + def __repr__(self) -> str: + size = len(self._elements) if self._elements is not None else "?" + return ( + f"{type(self).__name__}(degree={self._degree}, generators={len(self._labels)}, " + f"size={size}, monoid={self._monoid})" + ) + + def index(self, element: Transformation) -> int: + self._closure() + return self._index[element.images] + + def word(self, element: Transformation | int) -> tuple[Hashable, ...]: + """Return a shortest generator word evaluating to ``element``.""" + index = element if isinstance(element, int) else self.index(element) + self._closure() + return self._words[index] + + def evaluate(self, word: Iterable[Hashable]) -> Transformation: + """Return the transformation of a generator word (identity for ``()``).""" + images = tuple(range(self._degree)) + for label in word: + images = _compose(images, self._generators[label].images) + return Transformation(images) + + def multiply(self, i: int, j: int) -> int: + """Return the index of ``self[i] * self[j]``.""" + elements = self._closure() + return self._index[_compose(elements[i], elements[j])] + + def right_cayley(self) -> list[list[int]]: + """``right[i][g]`` is the index of ``self[i] * generator g``.""" + self._closure() + return self._right + + def left_cayley(self) -> list[list[int]]: + """``left[i][g]`` is the index of ``generator g * self[i]``.""" + elements = self._closure() + if self._left is None: + gens = [self._generators[label].images for label in self._labels] + self._left = [[self._index[_compose(g, e)] for g in gens] for e in elements] + return self._left + + def cayley_graph(self, side: Literal["right", "left"] = "right") -> nx.MultiDiGraph: + """Return the Cayley graph with edges labeled by generator.""" + import networkx as nx + + table = self.right_cayley() if side == "right" else self.left_cayley() + graph = nx.MultiDiGraph() + graph.add_nodes_from(range(len(table))) + for i, row in enumerate(table): + for label, j in zip(self._labels, row, strict=True): + graph.add_edge(i, j, symbol=label) + return graph + + # -- element-level properties --------------------------------------------- + + def idempotents(self) -> list[int]: + return [i for i, e in enumerate(self._closure()) if all(e[y] == y for y in e)] + + def is_aperiodic(self) -> bool: + """Return whether every cyclic subsemigroup has period one. + + Equivalently every subgroup is trivial, every H-class is a singleton + :cite:`Pin1986`, and -- for a syntactic monoid -- the language is + star-free :cite:`Schutzenberger1965`. + """ + return all(Transformation(e).index_period()[1] == 1 for e in self._closure()) + + def is_synchronizing(self) -> bool: + """Return whether some word sends every state to one state (``S`` contains a constant). + + Equivalently the minimal ideal has rank one. A word reaches a constant + iff every pair of states can be merged, so this is decided on pairs + without enumerating ``S`` :cite:`Volkov2008`. + """ + if self._degree == 1: + return True + gens = [g.images for g in self._generators.values()] + pairs = {(x, y) for x in range(self._degree) for y in range(x + 1, self._degree)} + mergeable: set[tuple[int, int]] = set() + changed = True + while changed: + changed = False + for pair in pairs - mergeable: + x, y = pair + for g in gens: + a, b = sorted((g[x], g[y])) + if a == b or (a, b) in mergeable: + mergeable.add(pair) + changed = True + break + return mergeable == pairs + + def is_group(self) -> bool: + """Return whether ``S`` is a permutation group (all generators are permutations).""" + return all(g.is_permutation() for g in self._generators.values()) + + # -- delegated constructions ---------------------------------------------- + + def green(self) -> GreenStructure: + """Return Green's relations; see :func:`~sofic.algebra.green.green_structure`.""" + from sofic.algebra.green import green_structure + + return green_structure(self) + + def holonomy(self, **kwargs: Any) -> HolonomyDecomposition: + """Return the holonomy decomposition; see :func:`~sofic.algebra.holonomy.holonomy_decomposition`.""" + from sofic.algebra.holonomy import holonomy_decomposition + + return holonomy_decomposition(self, **kwargs) + + def krohn_rhodes_complexity(self) -> ComplexityBounds: + """Return complexity bounds; see :func:`~sofic.algebra.complexity.krohn_rhodes_complexity`.""" + from sofic.algebra.complexity import krohn_rhodes_complexity + + return krohn_rhodes_complexity(self) + + +def _sorted_states(states: Iterable[Hashable]) -> list[Hashable]: + items = list(states) + try: + return sorted(items) # type: ignore[type-var] + except TypeError: + return sorted(items, key=repr) diff --git a/sofic/algebra/varieties.py b/sofic/algebra/varieties.py new file mode 100644 index 0000000..a30a3f9 --- /dev/null +++ b/sofic/algebra/varieties.py @@ -0,0 +1,141 @@ +"""Varieties of regular languages, decided on the syntactic monoid. + +Eilenberg's variety theorem matches classes of regular languages closed +under Boolean operations and quotients with pseudovarieties of finite monoids +(or semigroups), so membership is decided by checking an identity on the +syntactic monoid or semigroup :cite:`Eilenberg1976` :cite:`Pin1986`. +Products here follow word order: ``s * t`` is the image of ``uv`` when ``s`` +and ``t`` are the images of ``u`` and ``v``. + +Each test accepts a language (anything +:func:`~sofic.algebra.syntactic.syntactic_monoid` accepts, plus an +``alphabet`` keyword) or an already-built +:class:`~sofic.algebra.transformation.TransformationSemigroup`, which is then +used as is. +""" + +from __future__ import annotations + +import itertools +from typing import Any + +from sofic.algebra.transformation import TransformationSemigroup + +__all__ = [ + "is_commutative_language", + "is_definite", + "is_group_language", + "is_j_trivial", + "is_l_trivial", + "is_locally_testable", + "is_piecewise_testable", + "is_r_trivial", + "is_reverse_definite", +] + + +def _monoid(language: Any, **kwargs: Any) -> TransformationSemigroup: + if isinstance(language, TransformationSemigroup): + return language + from sofic.algebra.syntactic import syntactic_monoid + + return syntactic_monoid(language, **kwargs) + + +def _semigroup(language: Any, **kwargs: Any) -> TransformationSemigroup: + if isinstance(language, TransformationSemigroup): + return language + from sofic.algebra.syntactic import syntactic_monoid + + return syntactic_monoid(language, monoid=False, **kwargs) + + +def is_j_trivial(language: Any, **kwargs: Any) -> bool: + """Return whether every J-class of the syntactic monoid is a singleton.""" + return all(d.size == 1 for d in _monoid(language, **kwargs).green().d_classes) + + +def is_r_trivial(language: Any, **kwargs: Any) -> bool: + """Return whether every R-class of the syntactic monoid is a singleton.""" + return all(len(r) == 1 for r in _monoid(language, **kwargs).green().r_classes) + + +def is_l_trivial(language: Any, **kwargs: Any) -> bool: + """Return whether every L-class of the syntactic monoid is a singleton.""" + return all(len(c) == 1 for c in _monoid(language, **kwargs).green().l_classes) + + +def is_piecewise_testable(language: Any, **kwargs: Any) -> bool: + """Return whether membership depends only on the scattered subwords up to some length. + + Simon's theorem: a language is piecewise testable iff its syntactic monoid + is J-trivial :cite:`Simon1975`. + + Examples + -------- + >>> is_piecewise_testable("a*b*", alphabet="ab"), is_piecewise_testable("(a|b)*a", alphabet="ab") + (True, False) + """ + return is_j_trivial(language, **kwargs) + + +def is_locally_testable(language: Any, **kwargs: Any) -> bool: + """Return whether membership depends only on prefixes, suffixes, and factors of bounded length. + + Brzozowski-Simon and McNaughton: a language is locally testable iff its + syntactic *semigroup* is locally idempotent and commutative, that is, + ``eSe`` is an idempotent commutative monoid for every idempotent ``e`` + :cite:`BrzozowskiSimon1973` :cite:`McNaughton1974`. + + Examples + -------- + >>> is_locally_testable("(ab)*", alphabet="ab"), is_locally_testable("a*ba*", alphabet="ab") + (True, False) + """ + sg = _semigroup(language, **kwargs) + elements = list(range(len(sg))) + for e in sg.idempotents(): + local = {sg.multiply(sg.multiply(e, s), e) for s in elements} + for x in local: + if sg.multiply(x, x) != x: + return False + for x, y in itertools.combinations(local, 2): + if sg.multiply(x, y) != sg.multiply(y, x): + return False + return True + + +def is_definite(language: Any, **kwargs: Any) -> bool: + """Return whether membership depends only on a bounded-length suffix. + + A language is definite iff its syntactic semigroup satisfies ``s e = e`` + for every element ``s`` and idempotent ``e`` :cite:`Pin1986`. + + Examples + -------- + >>> is_definite("(a|b)*a", alphabet="ab"), is_definite("a(a|b)*", alphabet="ab") + (True, False) + """ + sg = _semigroup(language, **kwargs) + return all(sg.multiply(s, e) == e for e in sg.idempotents() for s in range(len(sg))) + + +def is_reverse_definite(language: Any, **kwargs: Any) -> bool: + """Return whether membership depends only on a bounded-length prefix. + + A language is reverse definite iff its syntactic semigroup satisfies + ``e s = e`` for every element ``s`` and idempotent ``e`` :cite:`Pin1986`. + """ + sg = _semigroup(language, **kwargs) + return all(sg.multiply(e, s) == e for e in sg.idempotents() for s in range(len(sg))) + + +def is_group_language(language: Any, **kwargs: Any) -> bool: + """Return whether the syntactic monoid is a group (the minimal DFA is a permutation automaton).""" + return _monoid(language, **kwargs).is_group() + + +def is_commutative_language(language: Any, **kwargs: Any) -> bool: + """Return whether membership is invariant under permuting letters (commutative syntactic monoid).""" + gens = list(_monoid(language, **kwargs).generators.values()) + return all(a * b == b * a for a, b in itertools.combinations(gens, 2)) diff --git a/sofic/automata/base.py b/sofic/automata/base.py index ae76d07..ec3a83b 100644 --- a/sofic/automata/base.py +++ b/sofic/automata/base.py @@ -10,6 +10,8 @@ from sofic.graph import ATTR_SYMBOL, EPSILON if TYPE_CHECKING: + from sofic.algebra.complexity import ComplexityBounds + from sofic.algebra.transformation import TransformationSemigroup from sofic.automata.dfa import DFA from sofic.automata.nfa import NFA @@ -148,6 +150,26 @@ def is_deterministic(self) -> bool: return is_deterministic_automaton(self) + def transition_semigroup(self, **kwargs: Any) -> TransformationSemigroup: + """Return the transition monoid of this deterministic automaton. + + See :func:`~sofic.algebra.adapters.transformation_semigroup`. + """ + from sofic.algebra.adapters import transformation_semigroup + + return transformation_semigroup(self, **kwargs) + + def krohn_rhodes_complexity(self) -> ComplexityBounds: + """Krohn-Rhodes complexity bounds of the *language*, via its syntactic monoid. + + Zero exactly when the language is star-free :cite:`Schutzenberger1965`. + See :func:`~sofic.algebra.complexity.krohn_rhodes_complexity`. + """ + from sofic.algebra.complexity import krohn_rhodes_complexity + from sofic.algebra.syntactic import syntactic_monoid + + return krohn_rhodes_complexity(syntactic_monoid(self)) + def delta(self, state: Hashable, symbol: Any) -> set[Hashable]: targets: set[Hashable] = set() for transition in self.graph.out_transitions(state): diff --git a/sofic/automata/wheeler.py b/sofic/automata/wheeler.py index f03689f..4f0a619 100644 --- a/sofic/automata/wheeler.py +++ b/sofic/automata/wheeler.py @@ -37,6 +37,33 @@ class WheelerError(SoficValidationError): """Raised when a Wheeler invariant is violated or cannot be established.""" +def star_free_obstruction(machine: Any, *, max_size: int = 10_000) -> str: + """Explain whether star-freeness rules out every Wheeler presentation of ``machine``. + + Wheeler languages are star-free :cite:`Alanko2021`, and a language is + star-free iff its syntactic monoid is aperiodic :cite:`Schutzenberger1965`. + The path language of ``machine`` is tested with + :func:`~sofic.algebra.syntactic.factor_monoid`; the test is skipped when + the monoid exceeds ``max_size`` elements. + """ + from sofic.algebra.syntactic import factor_monoid + from sofic.exceptions import SemigroupTooLargeError + + try: + aperiodic = factor_monoid(machine, max_size=max_size).is_aperiodic() + except SemigroupTooLargeError: + return ( + "Wheeler languages are star-free, so a language with a nontrivial syntactic group " + "has no Wheeler presentation at any order" + ) + if aperiodic: + return "the language is star-free, so a deeper refinement may still be Wheeler" + return ( + "the language is not star-free (its syntactic monoid contains a nontrivial group) and " + "Wheeler languages are star-free, so there is no Wheeler presentation at any order" + ) + + @dataclass(frozen=True, slots=True) class LabeledGraph: """Symbol-labeled view of a model, shared by every Wheeler routine. diff --git a/sofic/exceptions.py b/sofic/exceptions.py index 49259c1..0d246c9 100644 --- a/sofic/exceptions.py +++ b/sofic/exceptions.py @@ -41,6 +41,15 @@ class LumpabilityError(SoficValidationError): """Raised when a partition is not strongly lumpable for a model.""" +class SemigroupTooLargeError(SoficError): + """Raised when a transformation semigroup or cascade exceeds its size cap. + + A semigroup on ``n`` states can have up to ``n**n`` elements, so the caps + are a guard against accidental exponential blow-up, not a sign that the + input is malformed. + """ + + class InfiniteTransductionError(SoficError): """Raised when a finite input has infinitely many transducer outputs.""" diff --git a/sofic/generators/mealy.py b/sofic/generators/mealy.py index 2777b56..667f33b 100644 --- a/sofic/generators/mealy.py +++ b/sofic/generators/mealy.py @@ -11,6 +11,8 @@ from sofic.graph import ATTR_EMISSION, ATTR_PROB if TYPE_CHECKING: + from sofic.algebra.complexity import ComplexityBounds + from sofic.algebra.transformation import TransformationSemigroup from sofic.generators.lumping import LabelsLike, PartitionLike from sofic.generators.mixed_state import MixedState, MixedStatePresentation @@ -48,6 +50,26 @@ def add_transition(self, source: Hashable, target: Hashable, symbol: Any, prob: **{ATTR_EMISSION: symbol, ATTR_PROB: as_prob(prob), **attrs}, ) + def transition_semigroup(self, **kwargs: Any) -> TransformationSemigroup: + """Transition semigroup of the support of this unifilar presentation. + + Probabilities are dropped. For an ε-machine the presentation is the + canonical causal-state one, so the semigroup is a property of the + process. See :func:`~sofic.algebra.adapters.transformation_semigroup`. + """ + from sofic.algebra.adapters import transformation_semigroup + + return transformation_semigroup(self, **kwargs) + + def krohn_rhodes_complexity(self) -> ComplexityBounds: + """Krohn-Rhodes complexity bounds of :meth:`transition_semigroup`. + + See :func:`~sofic.algebra.complexity.krohn_rhodes_complexity`. + """ + from sofic.algebra.complexity import krohn_rhodes_complexity + + return krohn_rhodes_complexity(self.transition_semigroup()) + def validate_stochastic(self) -> None: super().validate_stochastic() validate_stochastic_edge_emissions( diff --git a/sofic/generators/wheeler_epsilon.py b/sofic/generators/wheeler_epsilon.py index 707222b..2b061de 100644 --- a/sofic/generators/wheeler_epsilon.py +++ b/sofic/generators/wheeler_epsilon.py @@ -26,7 +26,7 @@ import numpy as np -from sofic.automata.wheeler import WheelerError, WheelerOrder, wheeler_order +from sofic.automata.wheeler import WheelerError, WheelerOrder, star_free_obstruction, wheeler_order from sofic.generators.mealy import MealyHMM #: How many edge-machine refinements :func:`wheeler_presentation` will try. @@ -82,8 +82,7 @@ def wheeler_presentation( if wheeler_order(refined) is not None: return refined raise WheelerError( - f"no Wheeler presentation within {budget} edge-machine refinements; Wheeler languages " - "are star-free, so a process with a nontrivial syntactic group has none at any order" + f"no Wheeler presentation within {budget} edge-machine refinements; {star_free_obstruction(machine)}" ) diff --git a/sofic/shifts/__init__.py b/sofic/shifts/__init__.py index 4b9eb1c..a7bec3f 100644 --- a/sofic/shifts/__init__.py +++ b/sofic/shifts/__init__.py @@ -17,7 +17,13 @@ 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.invariants import ( + BowenFranksGroup, + bowen_franks_group, + factor_krohn_rhodes_complexity, + factor_syntactic_monoid, + 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 diff --git a/sofic/shifts/invariants.py b/sofic/shifts/invariants.py index b914f2f..a080226 100644 --- a/sofic/shifts/invariants.py +++ b/sofic/shifts/invariants.py @@ -99,6 +99,35 @@ def jordan_form_away_from_zero(matrix: np.ndarray | SymbolicModel) -> tuple[tupl return tuple(sorted(blocks, key=lambda block: (sp.default_sort_key(block[0]), block[1]))) +def factor_syntactic_monoid(shift: SymbolicModel, *, max_size: int | None = None) -> Any: + """Return the syntactic monoid of the shift's factor language ``B(X)``. + + ``B(X)`` is the set of words labeling finite paths in any presentation, + so this monoid depends only on the shift, not on the presentation + :cite:`LindMarcus1995`. It is a language invariant, not a conjugacy + invariant: recoding the alphabet by a sliding block code changes the + language. See :func:`~sofic.algebra.syntactic.syntactic_monoid`. + """ + from sofic.algebra.syntactic import factor_monoid + from sofic.algebra.transformation import DEFAULT_MAX_SIZE + + cap = DEFAULT_MAX_SIZE if max_size is None else max_size + return factor_monoid(shift, alphabet=shift.symbol_alphabet or None, max_size=cap) + + +def factor_krohn_rhodes_complexity(shift: SymbolicModel) -> Any: + """Krohn-Rhodes complexity bounds of :func:`factor_syntactic_monoid`. + + Zero exactly when ``B(X)`` is star-free :cite:`Schutzenberger1965`; every + shift of finite type has complexity zero, while the even shift, whose + language counts the symbol ``1`` modulo two, has complexity one. See + :func:`~sofic.algebra.complexity.krohn_rhodes_complexity`. + """ + from sofic.algebra.complexity import krohn_rhodes_complexity + + return krohn_rhodes_complexity(factor_syntactic_monoid(shift)) + + def _integer_matrix(matrix: np.ndarray | SymbolicModel) -> np.ndarray: if isinstance(matrix, SymbolicModel): return integer_adjacency_matrix(matrix)[0] diff --git a/sofic/shifts/sofic.py b/sofic/shifts/sofic.py index 937d316..34af6fa 100644 --- a/sofic/shifts/sofic.py +++ b/sofic/shifts/sofic.py @@ -7,6 +7,8 @@ from sofic.shifts.base import SymbolicModel if TYPE_CHECKING: + from sofic.algebra.complexity import ComplexityBounds + from sofic.algebra.transformation import TransformationSemigroup from sofic.generators.mealy import MealyHMM @@ -65,6 +67,25 @@ def topological_anatomy(self) -> dict[str, float]: return topological_anatomy(self) + def transition_semigroup(self, **kwargs: Any) -> TransformationSemigroup: + """Transition semigroup of this presentation, or of its right Fischer cover if not right-resolving. + + See :func:`~sofic.algebra.adapters.transformation_semigroup`. + """ + from sofic.algebra.adapters import transformation_semigroup + + return transformation_semigroup(self, **kwargs) + + def krohn_rhodes_complexity(self) -> ComplexityBounds: + """Krohn-Rhodes complexity bounds of the factor language's syntactic monoid. + + Presentation-independent; see + :func:`~sofic.shifts.invariants.factor_krohn_rhodes_complexity`. + """ + from sofic.shifts.invariants import factor_krohn_rhodes_complexity + + return factor_krohn_rhodes_complexity(self) + def markov_order(self) -> int | float: """Markov order ``R`` when the presentation is right-resolving (unifilar).""" from sofic.generators.synchronization import graph_from_sofic_shift, markov_order_from_graph diff --git a/sofic/shifts/wheeler.py b/sofic/shifts/wheeler.py index 2f351bf..ce7cd9a 100644 --- a/sofic/shifts/wheeler.py +++ b/sofic/shifts/wheeler.py @@ -19,7 +19,13 @@ from collections.abc import Hashable from typing import Any -from sofic.automata.wheeler import WheelerError, WheelerOrder, labeled_graph, wheeler_order +from sofic.automata.wheeler import ( + WheelerError, + WheelerOrder, + labeled_graph, + star_free_obstruction, + wheeler_order, +) from sofic.graph import ATTR_SYMBOL, TransitionGraph from sofic.shifts.covers import WheelerCover from sofic.shifts.sofic import SoficShift @@ -140,10 +146,7 @@ def wheeler_cover( found = wheeler_order(candidate) if found is not None: return _merge_wheeler_runs(candidate, found, follower_depth=follower_depth) - raise WheelerError( - f"no Wheeler presentation of order <= {max_order}; Wheeler languages are star-free, " - "so a shift that counts modulo anything (the even shift, for one) has none at any order" - ) + raise WheelerError(f"no Wheeler presentation of order <= {max_order}; {star_free_obstruction(trimmed)}") def is_wheeler_shift(shift: SoficShift, *, max_order: int = DEFAULT_MAX_ORDER) -> bool: diff --git a/tests/test_algebra_belief.py b/tests/test_algebra_belief.py new file mode 100644 index 0000000..ade35e4 --- /dev/null +++ b/tests/test_algebra_belief.py @@ -0,0 +1,38 @@ +"""Tests for belief-map semigroups of hidden Markov models.""" + +import pytest + +from sofic.algebra import krohn_rhodes_complexity, transformation_semigroup +from sofic.algebra.belief import mixed_state_semigroup, projective_semigroup +from sofic.examples.epsilon_machines import even_process, golden_mean +from sofic.examples.processes import nonunifilar_golden_mean, sns +from sofic.exceptions import SemigroupTooLargeError + + +@pytest.mark.parametrize("factory", [golden_mean, even_process]) +def test_unifilar_machines_have_finite_projective_semigroups(factory): + machine = factory() + projective = projective_semigroup(machine) + assert len(projective.states) == len(transformation_semigroup(machine)) + assert projective.states[0] == () + assert krohn_rhodes_complexity(projective) == krohn_rhodes_complexity(machine) + + +def test_projective_semigroup_acts_faithfully_on_itself(): + projective = projective_semigroup(even_process()) + assert len(projective) == len(projective.states) + + +def test_nonunifilar_presentations_can_still_be_finite(): + assert len(projective_semigroup(nonunifilar_golden_mean()).states) == 6 + + +def test_simple_nonunifilar_source_does_not_close(): + with pytest.raises(SemigroupTooLargeError): + projective_semigroup(sns(), max_size=500) + + +def test_mixed_state_semigroup_uses_the_mixed_state_presentation(): + sg = mixed_state_semigroup(even_process()) + assert len(sg.states) == 5 + assert krohn_rhodes_complexity(sg).exact == 1 diff --git a/tests/test_algebra_dedeo.py b/tests/test_algebra_dedeo.py new file mode 100644 index 0000000..4b5c292 --- /dev/null +++ b/tests/test_algebra_dedeo.py @@ -0,0 +1,189 @@ +"""Reproductions of DeDeo, "Effective Theories for Circuits and Automata" (Chaos 21, 037106, 2011).""" + +import itertools + +import numpy as np +import pytest + +from sofic.algebra import ( + TransformationSemigroup, + classify_elements, + coupled_random_maps, + element_walk, + holonomy_decomposition, + independent_random_maps, + is_full_transformation_monoid, + krohn_rhodes_complexity, + prime_decomposition, + walk_information, +) + + +def register(n, rule): + """Left-shift register on ``n`` bits: state ``(t_{i-n}, ..., t_{i-1})`` appends ``rule(state)``.""" + states = list(itertools.product((0, 1), repeat=n)) + index = {s: i for i, s in enumerate(states)} + return tuple(index[s[1:] + (rule(s),)] for s in states) + + +def truth_table_rule(table): + return lambda s: table[int("".join(map(str, s)), 2)] + + +def noisy_register(n, p_one): + """Stochastic matrix of an ``n``-bit register appending 1 with probability ``p_one[state]``.""" + states = list(itertools.product((0, 1), repeat=n)) + index = {s: i for i, s in enumerate(states)} + matrix = np.zeros((len(states), len(states))) + for s in states: + matrix[index[s], index[s[1:] + (1,)]] += p_one[index[s]] + matrix[index[s], index[s[1:] + (0,)]] += 1 - p_one[index[s]] + return matrix + + +#: Fig. 4a, transcribed from the arXiv TeX source (the 101 -> 010 edge, drawn +#: unlabeled, is that state's only edge and so carries both letters). +FIG_4A = [ + ("000", "001", "ab"), + ("001", "010", "b"), + ("100", "001", "b"), + ("100", "000", "a"), + ("010", "101", "ab"), + ("101", "010", "ab"), + ("001", "011", "a"), + ("110", "100", "a"), + ("011", "110", "b"), + ("110", "101", "b"), + ("011", "111", "a"), + ("111", "110", "ab"), +] +FIG_4A_NAMES = {"000": "A", "001": "B", "010": "C", "100": "D", "011": "E", "101": "F", "110": "G", "111": "H"} + + +def left_shift_powerset(n, **kwargs): + tables = itertools.product((0, 1), repeat=2**n) + return TransformationSemigroup({t: register(n, truth_table_rule(t)) for t in tables}, **kwargs) + + +# -- deterministic results (Secs. III-VI, appendix) ---------------------------------- + + +def test_two_bit_xor_is_a_mod_three_counter(): + xor = TransformationSemigroup({"a": register(2, lambda s: s[0] ^ s[1])}) + assert len(xor) == 3 + assert [level.order for level in holonomy_decomposition(xor).levels] == [3] + assert krohn_rhodes_complexity(xor).exact == 1 + + +def test_left_shift_powerset_is_the_full_transformation_monoid(): + lsp2 = left_shift_powerset(2) + assert len(lsp2) == 4**4 + assert is_full_transformation_monoid(lsp2) + + +def test_three_bit_left_shift_powerset_complexity_without_enumeration(): + lsp3 = left_shift_powerset(3, max_size=1000) + assert is_full_transformation_monoid(lsp3) + assert krohn_rhodes_complexity(lsp3).exact == 7 + assert lsp3._elements is None + + +def test_table_one_three_bit_rows(): + one_letter = TransformationSemigroup({"a": register(3, lambda s: s[0] ^ s[1])}) + assert prime_decomposition(one_letter).simple_group_orders() == (7,) + two_letters = TransformationSemigroup( + {"a": register(3, lambda s: s[0] ^ s[1]), "b": register(3, lambda s: s[0] ^ s[2])} + ) + (factor,) = prime_decomposition(two_letters).groups + assert (factor.kind, factor.order) == ("nonabelian simple", 168) + + +def test_figure_four_effective_theory(): + maps = {"a": {}, "b": {}} + for source, target, letters in FIG_4A: + for letter in letters: + maps[letter][FIG_4A_NAMES[source]] = FIG_4A_NAMES[target] + sg = TransformationSemigroup.from_maps(maps, states="ABCDEFGH") + hd = holonomy_decomposition(sg) + assert [level.order for level in hd.levels] == [3, 2, 2] + top = hd.effective_theory(1) + names = {prefix: "".join(sorted(states)) for prefix, states in top.superstates.items()} + # Fig. 4b: superstates AH, BG, DE, each also containing C and F. + assert sorted(names.values()) == ["ACFH", "BCFG", "CDEF"] + by_name = {name: prefix for prefix, name in names.items()} + + def step(name, letter): + (target,) = top.automaton.delta(by_name[name], letter) + return names[target] + + assert [step("ACFH", "a"), step("BCFG", "a"), step("CDEF", "a")] == ["BCFG", "CDEF", "ACFH"] + assert {step(name, "b") for name in names.values()} == {"BCFG"} + + +# -- probabilistic automata (Secs. VII-IX) --------------------------------------------- + +XOR_MAP = (0, 3, 1, 2) + + +def noisy_xor(eps): + return noisy_register(2, [eps, 1 - eps, 1 - eps, eps]) + + +def test_figure_five_decomposition(): + rep = independent_random_maps(noisy_register(2, [0.5, 1, 1, 0.75])) + assert sorted(rep.maps("t").values()) == pytest.approx([1 / 8, 1 / 8, 3 / 8, 3 / 8]) + levels = holonomy_decomposition(rep.semigroup()).levels + assert len(levels) == 4 + assert {level.order for level in levels} == {2, 3} + + +def _z3(sg): + xor = next(g for g in sg.generators.values() if g.images == XOR_MAP) + return {sg.index(xor**k) for k in range(3)} + + +def test_figure_six_independent_noise_decays_to_resets(): + rep = independent_random_maps(noisy_xor(0.01)) + assert len(rep.maps("t")) == 16 + sg = rep.semigroup() + assert is_full_transformation_monoid(sg) + classes = classify_elements(sg) + z3 = _z3(sg) + other = sorted(classes["units"] - z3) + walk = element_walk(rep, 1000) + on_z3 = walk[:, sorted(z3)].sum(axis=1) + on_resets = walk[:, sorted(classes["resets"])].sum(axis=1) + on_other = walk[:, other].sum(axis=1) + assert on_z3[0] == 1.0 and on_z3[-1] < 1e-6 + assert on_resets[-1] > 0.99 + # The brief window in which other permutations outweigh the resets. + assert (on_other[1:10] > on_resets[1:10]).all() + assert on_other[-1] < 1e-6 + + +def test_figure_eight_correlated_noise_stays_reversible(): + rep = coupled_random_maps(noisy_xor(0.01)) + weights = {m.images: w for m, w in rep.maps("t").items()} + assert weights == pytest.approx({XOR_MAP: 0.99, (1, 2, 0, 3): 0.01}) + sg = rep.semigroup() + assert sg.is_group() and len(sg) == 12 + walk = element_walk(rep, 5000) + # Both maps are 3-cycles in one coset of the Klein group, so the walk has + # period three: uniform on a coset at each step, uniform on A_4 on average. + assert np.count_nonzero(walk[-1] > 1e-9) == 4 + assert walk[-3:].mean(axis=0) == pytest.approx(np.full(12, 1 / 12), abs=1e-3) + assert walk[-1, sorted(_z3(sg))].sum() == pytest.approx(0.25, abs=1e-3) + + +def test_irreversibility_versus_dissipation(): + independent = walk_information(independent_random_maps(noisy_xor(0.01)), 1000) + correlated = walk_information(coupled_random_maps(noisy_xor(0.01)), 5000) + assert independent.irreversibility[-1] == pytest.approx(2.0, abs=0.01) + assert np.allclose(correlated.irreversibility, 0.0) + assert correlated.dissipation[-1] == pytest.approx(2.0, abs=0.01) + + +def test_full_transformation_recognition_needs_a_rank_deficient_generator(): + assert not is_full_transformation_monoid(TransformationSemigroup({"t": (1, 0, 2), "c": (1, 2, 0)})) + assert not is_full_transformation_monoid(TransformationSemigroup({"c": (1, 2, 0), "e": (0, 0, 2)})) + assert is_full_transformation_monoid(TransformationSemigroup({"t": (1, 0, 2), "c": (1, 2, 0), "e": (0, 0, 2)})) diff --git a/tests/test_algebra_holonomy.py b/tests/test_algebra_holonomy.py new file mode 100644 index 0000000..782e203 --- /dev/null +++ b/tests/test_algebra_holonomy.py @@ -0,0 +1,280 @@ +"""Tests for the holonomy decomposition, cascades, and Krohn-Rhodes complexity.""" + +import math + +import pytest +from hypothesis import given +from hypothesis import strategies as st +from sympy.combinatorics.named_groups import AlternatingGroup, DirectProduct, SymmetricGroup + +from sofic.algebra import ( + TransformationSemigroup, + holonomy_decomposition, + holonomy_information, + krohn_rhodes_complexity, + prime_decomposition, +) +from sofic.algebra._groups import composition_factor_orders +from sofic.automata.wheeler import WheelerError +from sofic.examples.epsilon_machines import ( + butterfly_process, + even_process, + golden_mean, + noisy_random_phase_slip, +) +from sofic.shifts import factor_krohn_rhodes_complexity, wheeler_cover +from sofic.shifts.sofic import SoficShift + +EXAMPLE_7_1 = {"s1": (1, 0, 2), "s2": (0, 1, 1)} +EXAMPLE_7_2 = { + "s1": (0, 1, 2, 0, 0, 0), + "s2": (3, 3, 3, 4, 3, 5), + "s3": (3, 3, 3, 4, 5, 3), + "s4": (3, 3, 3, 3, 4, 4), + "s5": (3, 3, 3, 0, 1, 2), + "s6": (1, 2, 0, 3, 3, 3), +} + + +def full_transformation_monoid(n, **kwargs): + rest = tuple(range(2, n)) + return TransformationSemigroup({"t": (1, 0, *rest), "c": (*range(1, n), 0), "e": (0, 0, *rest)}, **kwargs) + + +def shift(edges): + s = SoficShift(symbol_alphabet=frozenset(symbol for _, _, symbol in edges)) + for source, target, symbol in edges: + s.add_transition(source, target, symbol) + return s + + +EVEN_SHIFT = [("A", "A", "0"), ("A", "B", "1"), ("B", "A", "1")] +GOLDEN_MEAN_SHIFT = [("A", "A", "0"), ("A", "B", "1"), ("B", "A", "0")] + + +@st.composite +def semigroups(draw, max_degree=5, max_generators=3): + degree = draw(st.integers(1, max_degree)) + count = draw(st.integers(1, max_generators)) + gens = { + f"g{i}": tuple(draw(st.lists(st.integers(0, degree - 1), min_size=degree, max_size=degree))) + for i in range(count) + } + return TransformationSemigroup(gens) + + +# -- holonomy decomposition ------------------------------------------------------ + + +def test_egri_nagy_nehaniv_example_7_1(): + hd = holonomy_decomposition(TransformationSemigroup(EXAMPLE_7_1)) + assert len(hd.image_set) == 5 + assert [level.order for level in hd.levels] == [1, 2] + + +def test_egri_nagy_nehaniv_example_7_2(): + hd = holonomy_decomposition(TransformationSemigroup(EXAMPLE_7_2)) + assert len(hd.image_set) == 19 + assert hd.n_levels == 5 + groups = [sorted(int(g.order()) for g in level.groups) for level in hd.levels] + assert groups == [[1], [3], [2, 6], [2, 3], [1]] + + +@pytest.mark.parametrize("n", [2, 3, 4, 5]) +def test_full_transformation_monoid_has_symmetric_holonomy(n): + hd = holonomy_decomposition(full_transformation_monoid(n)) + assert [level.order for level in hd.levels] == [math.factorial(k) for k in range(n, 1, -1)] + + +def test_holonomy_needs_only_generators(): + hd = holonomy_decomposition(full_transformation_monoid(5, max_size=10)) + assert hd.n_levels == 4 + assert hd.semigroup._elements is None + + +@given(semigroups()) +def test_encode_and_decode_are_inverse(sg): + hd = holonomy_decomposition(sg) + for chain in hd.chains(): + assert hd.decode(hd.encode(chain)) == chain + for state in sg.states: + assert hd.state(hd.coordinates(state)) == state + + +@given(semigroups()) +def test_heights_respect_strict_subduction(sg): + hd = holonomy_decomposition(sg) + for parent, tiles in hd.tiles.items(): + assert all(hd.height[tile] < hd.height[parent] for tile in tiles) + assert frozenset().union(*tiles) == parent + + +def test_invalid_coordinates_are_rejected(): + hd = holonomy_decomposition(TransformationSemigroup(EXAMPLE_7_1)) + with pytest.raises(ValueError): + hd.decode((None, None)) + with pytest.raises(ValueError): + hd.decode((0,)) + + +# -- cascades ------------------------------------------------------------------------- + + +@given(semigroups()) +def test_the_holonomy_cascade_emulates_the_semigroup(sg): + hd = holonomy_decomposition(sg) + cascade = hd.to_cascade() + assert cascade.emulates(sg) + for coords in cascade.states: + for label in sg.generator_labels: + assert cascade.act(coords, label) == hd.act(coords, label) + + +@given(semigroups(), st.lists(st.integers(0, 2), max_size=6)) +def test_cascade_runs_track_the_original_machine(sg, picks): + labels = sg.generator_labels + word = [labels[i % len(labels)] for i in picks] + cascade = holonomy_decomposition(sg).to_cascade() + g = sg.evaluate(word) + for x, state in enumerate(sg.states): + assert cascade.decoding[cascade.run(word, state)] == sg.states[g(x)] + + +@given(semigroups(), st.integers(0, 6)) +def test_truncation_commutes_with_every_input(sg, k): + cascade = holonomy_decomposition(sg).to_cascade() + k = min(k, cascade.n_levels) + coarse = cascade.truncate(k) + for coords, state in cascade.decoding.items(): + superstate = coarse.decoding[coords[:k]] + assert isinstance(superstate, frozenset) and state in superstate + for label in sg.generator_labels: + assert cascade.act(coords, label)[:k] == coarse.act(coords[:k], label) + + +def test_effective_theory_superstates_overlap(): + hd = holonomy_decomposition(TransformationSemigroup(EXAMPLE_7_2)) + top = hd.effective_theory(1) + assert sorted(map(sorted, top.superstates.values())) == [[0, 1, 2, 3], [3, 4, 5]] + assert len(list(top.automaton.states())) == 2 + assert hd.effective_theory(0).superstates == {(): frozenset(range(6))} + assert len(hd.effective_theory(hd.n_levels).superstates) == len(hd.to_cascade().states) + + +def test_effective_theory_of_a_process_reports_lumpability(): + machine = noisy_random_phase_slip() + top = holonomy_decomposition(machine).effective_theory(1, machine=machine) + assert sum(top.distribution.values()) == pytest.approx(1.0) + assert top.lumpable is False + finest = holonomy_decomposition(machine) + assert finest.effective_theory(finest.n_levels, machine=machine).lumpable is True + + +def test_cascade_as_dfa(): + sg = TransformationSemigroup(EXAMPLE_7_2) + cascade = holonomy_decomposition(sg).to_cascade() + dfa = cascade.as_dfa(initial=0) + assert len(list(dfa.states())) == len(cascade.states) + assert dfa.is_deterministic() + + +def test_wrong_semigroup_is_not_emulated(): + cascade = holonomy_decomposition(TransformationSemigroup(EXAMPLE_7_1)).to_cascade() + assert not cascade.emulates(TransformationSemigroup({"s1": (0, 1, 2), "s2": (0, 1, 1)})) + + +# -- complexity ------------------------------------------------------------------------- + + +@pytest.mark.parametrize("n", [2, 3, 4, 5]) +def test_full_transformation_monoid_complexity(n): + assert krohn_rhodes_complexity(full_transformation_monoid(n)) == (n - 1, n - 1, n - 1) + + +@pytest.mark.parametrize("n", [2, 3, 5, 6, 7]) +def test_cyclic_groups_have_complexity_one(n): + sg = TransformationSemigroup({"a": (*range(1, n), 0)}) + assert krohn_rhodes_complexity(sg).exact == 1 + expected = [p for p in range(2, n + 1) if n % p == 0 and all(p % q for q in range(2, p))] + assert set(prime_decomposition(sg).simple_group_orders()) == set(expected) + + +def test_flip_flop_is_aperiodic(): + flip_flop = TransformationSemigroup({"r0": (0, 0), "r1": (1, 1)}) + assert krohn_rhodes_complexity(flip_flop).exact == 0 + decomposition = prime_decomposition(flip_flop) + assert decomposition.groups == () and decomposition.needs_flip_flop + + +def test_lower_bound_degrades_gracefully_when_enumeration_is_capped(): + assert krohn_rhodes_complexity(full_transformation_monoid(5, max_size=10)).exact == 4 + rank_three = TransformationSemigroup( + {"t": (1, 0, 2, 3, 4), "c": (1, 2, 3, 4, 0), "e": (0, 0, 0, 3, 4)}, max_size=10 + ) + bounds = krohn_rhodes_complexity(rank_three) + assert bounds.lower == 1 and bounds.upper > 1 and bounds.exact is None + + +@given(semigroups()) +def test_complexity_bounds_are_consistent(sg): + bounds = krohn_rhodes_complexity(sg) + assert 0 <= bounds.lower <= bounds.upper + assert (bounds.upper == 0) == sg.is_aperiodic() + assert bounds.exact == (bounds.lower if bounds.lower == bounds.upper else None) + + +def test_a5_is_a_nonabelian_simple_prime(): + sg = TransformationSemigroup({"a": (1, 2, 3, 4, 0), "b": (1, 2, 0, 3, 4)}) + (factor,) = prime_decomposition(sg).groups + assert (factor.kind, factor.order) == ("nonabelian simple", 60) + assert not prime_decomposition(sg).needs_flip_flop + + +@pytest.mark.parametrize( + ("group", "orders"), + [ + (AlternatingGroup(5), [60]), + (SymmetricGroup(5), [2, 60]), + (SymmetricGroup(6), [2, 360]), + (SymmetricGroup(4), [2, 2, 2, 3]), + (DirectProduct(AlternatingGroup(5), AlternatingGroup(5)), [60, 60]), + (DirectProduct(AlternatingGroup(5), SymmetricGroup(3)), [2, 3, 60]), + ], +) +def test_composition_factor_orders(group, orders): + assert composition_factor_orders(group) == orders + + +# -- processes and shifts ---------------------------------------------------------------- + + +def test_even_and_golden_mean_processes(): + assert even_process().krohn_rhodes_complexity().exact == 1 + assert golden_mean().krohn_rhodes_complexity().exact == 0 + + +def test_shift_complexity_ignores_the_presentation(): + redundant = shift([*EVEN_SHIFT, ("A", "C", "1"), ("C", "A", "1")]) + assert not redundant.is_unifilar() + assert factor_krohn_rhodes_complexity(shift(EVEN_SHIFT)).exact == 1 + assert redundant.krohn_rhodes_complexity().exact == 1 + assert shift(GOLDEN_MEAN_SHIFT).krohn_rhodes_complexity().exact == 0 + + +def test_non_right_resolving_shifts_use_the_fischer_cover(): + redundant = shift([*EVEN_SHIFT, ("A", "C", "1"), ("C", "A", "1")]) + assert len(redundant.transition_semigroup()) == len(shift(EVEN_SHIFT).transition_semigroup()) + + +def test_wheeler_failure_names_the_obstruction(): + with pytest.raises(WheelerError, match="not star-free"): + wheeler_cover(shift(EVEN_SHIFT)) + + +@pytest.mark.parametrize("factory", [even_process, golden_mean, noisy_random_phase_slip, butterfly_process]) +def test_holonomy_information_splits_statistical_complexity(factory): + machine = factory() + info = holonomy_information(machine) + assert info.total == pytest.approx(machine.statistical_complexity()) + assert sum(info.level_entropies) == pytest.approx(info.total) + assert all(h >= -1e-12 for h in info.level_entropies) diff --git a/tests/test_algebra_random_maps.py b/tests/test_algebra_random_maps.py new file mode 100644 index 0000000..b6d6920 --- /dev/null +++ b/tests/test_algebra_random_maps.py @@ -0,0 +1,186 @@ +"""Tests for random-map representations, semigroup walks, and minimal ideals.""" + +import dit +import numpy as np +import pytest +from dit.shannon import mutual_information +from hypothesis import given +from hypothesis import strategies as st + +from sofic.algebra import TransformationSemigroup +from sofic.algebra.random_maps import ( + classify_elements, + coupled_random_maps, + element_walk, + independent_random_maps, + random_maps, + walk_information, +) +from sofic.automata.unifilar import UnifilarAutomaton +from sofic.examples.epsilon_machines import even_process + + +@st.composite +def stochastic_matrices(draw, max_states=4): + n = draw(st.integers(1, max_states)) + rows = [] + for _ in range(n): + support = draw(st.lists(st.integers(0, n - 1), min_size=1, max_size=min(n, 3), unique=True)) + raw = [draw(st.integers(1, 5)) for _ in support] + row = np.zeros(n) + row[support] = np.asarray(raw) / sum(raw) + rows.append(row) + return np.asarray(rows) + + +@st.composite +def semigroups(draw, max_degree=5, max_generators=3): + degree = draw(st.integers(1, max_degree)) + count = draw(st.integers(1, max_generators)) + gens = { + f"g{i}": tuple(draw(st.lists(st.integers(0, degree - 1), min_size=degree, max_size=degree))) + for i in range(count) + } + return TransformationSemigroup(gens) + + +# -- representations ------------------------------------------------------------------ + + +@given(stochastic_matrices()) +def test_both_constructions_reproduce_the_matrix(matrix): + for rep in (independent_random_maps(matrix), coupled_random_maps(matrix)): + assert np.allclose(rep.matrices()["t"], matrix) + + +@given(stochastic_matrices()) +def test_coupling_uses_few_maps(matrix): + supports = (matrix > 0).sum(axis=1) + assert len(coupled_random_maps(matrix).maps("t")) <= 1 + int((supports - 1).sum()) + + +def test_independent_support_is_the_product_of_row_supports(): + matrix = np.array([[0.5, 0.5, 0.0], [0.0, 0.2, 0.8], [1.0, 0.0, 0.0]]) + assert len(independent_random_maps(matrix).maps("t")) == 4 + with pytest.raises(ValueError, match="max_maps"): + independent_random_maps(matrix, max_maps=3) + + +def test_explicit_maps_are_checked(): + matrix = np.array([[0.5, 0.5], [0.5, 0.5]]) + random_maps({"t": {(0, 1): 0.5, (1, 0): 0.5}}, matrices=matrix) + with pytest.raises(ValueError, match="reproduce"): + random_maps({"t": {(0, 0): 0.5, (0, 1): 0.5}}, matrices=matrix) + with pytest.raises(ValueError, match="sum"): + random_maps({"t": {(0, 1): 0.4}}) + + +def test_markov_models_are_read_through_their_state_matrix(): + rep = independent_random_maps(even_process()) + assert rep.states == ("A", "B") + assert len(rep.maps("t")) == 2 + + +def test_rejects_non_stochastic_input(): + with pytest.raises(ValueError, match="row-stochastic"): + independent_random_maps(np.array([[0.5, 0.2], [0.0, 1.0]])) + + +# -- walks ------------------------------------------------------------------------------ + + +@given(stochastic_matrices(), st.integers(0, 6)) +def test_walk_marginals_are_matrix_powers(matrix, steps): + for rep in (independent_random_maps(matrix), coupled_random_maps(matrix)): + sg = rep.semigroup() + walk = element_walk(rep, steps) + assert np.allclose(walk.sum(axis=1), 1.0) + n = len(matrix) + reached = np.zeros((n, n)) + for i, element in enumerate(sg): + reached[np.arange(n), element.images] += walk[steps, i] + assert np.allclose(reached, np.linalg.matrix_power(matrix, steps)) + + +def test_multi_letter_walks_need_inputs(): + rep = random_maps({"a": {(1, 0): 1.0}, "b": {(0, 0): 1.0}}) + with pytest.raises(ValueError, match="inputs"): + element_walk(rep, 3) + walk = element_walk(rep, inputs=["a", "a", "b"]) + sg = rep.semigroup() + assert sg[int(np.argmax(walk[2]))] == sg.evaluate([]) + assert sg[int(np.argmax(walk[3]))].images == (0, 0) + + +@given(stochastic_matrices(), st.integers(0, 5)) +def test_information_losses_are_nonnegative_and_add_up(matrix, steps): + info = walk_information(coupled_random_maps(matrix), steps) + h0 = np.log2(len(matrix)) + assert (info.irreversibility >= -1e-9).all() + assert (info.dissipation >= -1e-9).all() + assert np.allclose(info.irreversibility + info.dissipation + info.mutual, h0) + + +def test_walk_information_matches_dit(): + matrix = np.array([[0.7, 0.3, 0.0], [0.0, 0.4, 0.6], [0.5, 0.0, 0.5]]) + rep = independent_random_maps(matrix) + info = walk_information(rep, 3) + p_t = np.linalg.matrix_power(matrix, 3) / 3 + joint = dit.Distribution.from_ndarray(p_t) + assert info.mutual[3] == pytest.approx(mutual_information(joint, [0], [1])) + + +def test_permutation_noise_never_irreversible(): + rep = coupled_random_maps(np.array([[0.9, 0.1], [0.1, 0.9]])) + info = walk_information(rep, 50) + assert np.allclose(info.irreversibility, 0.0) + assert info.dissipation[-1] == pytest.approx(1.0, abs=1e-6) + + +# -- minimal ideal, synchronization, Schutzenberger groups ------------------------------- + + +@given(semigroups()) +def test_synchronizing_iff_minimal_rank_one(sg): + green = sg.green() + assert sg.is_synchronizing() == (green.minimal_rank() == 1) + ideal = set(green.minimal_ideal().elements()) + lowest = min(e.rank() for e in sg) + assert ideal == {i for i, e in enumerate(sg) if e.rank() == lowest} + + +@given(semigroups()) +def test_schutzenberger_group_order_is_h_class_size(sg): + green = sg.green() + for members in green.h_classes: + assert int(green.schutzenberger_group(members[0]).order()) == len(members) + for e, group in green.maximal_subgroups(): + assert int(green.schutzenberger_group(e).order()) == int(group.order()) + + +def unifilar(alphabet, edges): + aut = UnifilarAutomaton(input_alphabet=frozenset(alphabet)) + for source, target, _ in edges: + aut.graph.add_state(source) + aut.graph.add_state(target) + for source, target, symbol in edges: + aut.graph.add_transition(source, target, symbol=symbol) + return aut + + +def test_synchronization_agrees_with_reset_threshold(): + cerny = unifilar("ab", [(0, 1, "a"), (1, 2, "a"), (2, 0, "a"), (0, 0, "b"), (1, 1, "b"), (2, 0, "b")]) + assert cerny.transition_semigroup().is_synchronizing() + assert cerny.reset_threshold() == 4 + flip = unifilar("a", [(0, 1, "a"), (1, 0, "a")]) + assert not flip.transition_semigroup().is_synchronizing() + assert flip.reset_threshold() == float("inf") + + +def test_classify_elements_of_t3(): + t3 = TransformationSemigroup({"t": (1, 0, 2), "c": (1, 2, 0), "e": (0, 0, 2)}) + classes = classify_elements(t3) + assert len(classes["units"]) == 6 + assert classes["resets"] == classes["minimal_ideal"] + assert len(classes["resets"]) == 3 + assert len(classes["identity"]) == 1 diff --git a/tests/test_algebra_semigroups.py b/tests/test_algebra_semigroups.py new file mode 100644 index 0000000..6748c5e --- /dev/null +++ b/tests/test_algebra_semigroups.py @@ -0,0 +1,246 @@ +"""Tests for transformation semigroups, Green's relations, and syntactic monoids.""" + +import pytest +from hypothesis import given +from hypothesis import strategies as st + +from sofic.algebra import ( + SINK, + Transformation, + TransformationSemigroup, + is_star_free, + syntactic_monoid, + transformation_semigroup, +) +from sofic.automata.dfa import DFA +from sofic.automata.nfa import NFA +from sofic.examples.epsilon_machines import even_process, golden_mean +from sofic.exceptions import SemigroupTooLargeError + + +def full_transformation_monoid(n): + """T_n from a transposition, an n-cycle, and an elementary collapse.""" + rest = tuple(range(2, n)) + return TransformationSemigroup( + {"t": (1, 0, *rest), "c": (*range(1, n), 0), "e": (0, 0, *rest)}, + ) + + +@st.composite +def semigroups(draw, max_degree=5, max_generators=3): + degree = draw(st.integers(1, max_degree)) + count = draw(st.integers(1, max_generators)) + gens = { + f"g{i}": tuple(draw(st.lists(st.integers(0, degree - 1), min_size=degree, max_size=degree))) + for i in range(count) + } + return TransformationSemigroup(gens, monoid=draw(st.booleans())) + + +# -- transformations ---------------------------------------------------------- + + +def test_composition_acts_on_the_right(): + s = Transformation((1, 2, 0)) + t = Transformation((0, 0, 2)) + assert (s * t).images == tuple(t(s(x)) for x in range(3)) + + +def test_transformation_basics(): + t = Transformation((2, 2, 0, 1)) + assert t.rank() == 3 + assert t.image() == {0, 1, 2} + assert t.kernel() == (frozenset({0, 1}), frozenset({2}), frozenset({3})) + assert not t.is_permutation() + assert Transformation.identity(3).is_idempotent() + assert Transformation.constant(3, 1).is_constant() + assert t.power(0) == Transformation.identity(4) + assert t.restrict([0, 2]) == (1, 0) + + +@given(st.lists(st.integers(0, 5), min_size=6, max_size=6)) +def test_index_period_is_the_first_repeat(images): + s = Transformation(tuple(images)) + m, r = s.index_period() + assert s ** (m + r) == s**m + assert all(s ** (k + r) != s**k for k in range(1, m)) + + +# -- semigroups ----------------------------------------------------------------- + + +@pytest.mark.parametrize("n", [2, 3, 4, 5]) +def test_full_transformation_monoid_has_n_to_the_n_elements(n): + assert len(full_transformation_monoid(n)) == n**n + + +def test_monoid_flag_adjoins_the_identity(): + # Example 7.2 of Egri-Nagy & Nehaniv: |S| = 138 without the identity. + gens = { + "s1": (0, 1, 2, 0, 0, 0), + "s2": (3, 3, 3, 4, 3, 5), + "s3": (3, 3, 3, 4, 5, 3), + "s4": (3, 3, 3, 3, 4, 4), + "s5": (3, 3, 3, 0, 1, 2), + "s6": (1, 2, 0, 3, 3, 3), + } + assert len(TransformationSemigroup(gens, monoid=False)) == 138 + assert len(TransformationSemigroup(gens)) == 139 + + +def test_max_size_guard_raises_and_resets(): + t5 = full_transformation_monoid(5) + t5._max_size = 100 + with pytest.raises(SemigroupTooLargeError): + len(t5) + t5._max_size = 10_000 + assert len(t5) == 3125 + + +def test_bad_generators_are_rejected(): + with pytest.raises(ValueError): + TransformationSemigroup({"a": (0, 1), "b": (0, 1, 2)}) + with pytest.raises(ValueError): + TransformationSemigroup({"a": (0, 3)}) + with pytest.raises(ValueError): + TransformationSemigroup({}, monoid=False) + + +def test_from_maps_adds_a_sink_only_when_needed(): + total = TransformationSemigroup.from_maps({"a": {0: 1, 1: 0}}) + assert SINK not in total.states + partial = TransformationSemigroup.from_maps({"a": {0: 1}}) + assert partial.states == (0, 1, SINK) + assert partial.generators["a"].images == (1, 2, 2) + with pytest.raises(ValueError): + TransformationSemigroup.from_maps({"a": {0: 1}}, complete="error") + + +@given(semigroups()) +def test_words_evaluate_to_their_elements(sg): + for i, element in enumerate(sg): + assert sg.evaluate(sg.word(i)) == element + assert sg.index(element) == i + + +@given(semigroups()) +def test_cayley_tables_are_products(sg): + right, left = sg.right_cayley(), sg.left_cayley() + gens = [sg.generators[label] for label in sg.generator_labels] + for i, element in enumerate(sg): + for k, g in enumerate(gens): + assert sg[right[i][k]] == element * g + assert sg[left[i][k]] == g * element + + +def test_cayley_graph_edges_are_labeled(): + graph = TransformationSemigroup({"a": (1, 0)}).cayley_graph() + assert {d["symbol"] for _, _, d in graph.edges(data=True)} == {"a"} + + +# -- Green's relations ------------------------------------------------------------ + + +def test_green_structure_of_t3(): + green = full_transformation_monoid(3).green() + assert sorted(d.size for d in green.d_classes) == [3, 6, 18] + assert sorted(int(g.order()) for _, g in green.maximal_subgroups()) == [1, 2, 6] + assert not green.is_aperiodic() + assert "D0" in green.egg_box() + + +@given(semigroups()) +def test_green_classes_partition_and_form_egg_boxes(sg): + green = sg.green() + n = len(sg) + for classes in (green.r_classes, green.l_classes, green.h_classes): + assert sorted(i for c in classes for i in c) == list(range(n)) + assert sorted(i for d in green.d_classes for i in d.elements()) == list(range(n)) + for d in green.d_classes: + sizes = {len(cell) for row in d.cells for cell in row} + assert len(sizes) == 1 and 0 not in sizes + + +@given(semigroups()) +def test_aperiodicity_agrees_with_trivial_h_classes(sg): + assert sg.is_aperiodic() == sg.green().is_aperiodic() + + +@given(semigroups()) +def test_maximal_subgroup_orders_match_their_h_classes(sg): + green = sg.green() + for e, group in green.maximal_subgroups(): + assert int(group.order()) == len(green.h_classes[green.h_of[e]]) + + +# -- syntactic monoids -------------------------------------------------------------- + + +def test_syntactic_monoid_known_values(): + assert len(syntactic_monoid("(aa)*", alphabet="a")) == 2 + assert len(syntactic_monoid("(a|b)*", alphabet="ab")) == 1 + assert not syntactic_monoid("(aaa)*", alphabet="a").is_aperiodic() + + +def test_syntactic_monoid_ignores_the_presentation(): + nfa = NFA.from_regex("(ab)*", alphabet="ab") + dfa = nfa.determinize() + assert len(syntactic_monoid(nfa)) == len(syntactic_monoid(dfa)) == len(syntactic_monoid("(ab)*", alphabet="ab")) + + +def test_unused_symbols_act_as_zero(): + assert len(syntactic_monoid("(aa)*", alphabet="a")) == 2 + assert len(syntactic_monoid("(aa)*", alphabet="ab")) == 3 + + +@st.composite +def star_free_regexes(draw, depth=2): + """Expressions from words and ``A* w A*`` closed under union and concatenation.""" + word = draw(st.text(alphabet="ab", min_size=1, max_size=3)) + atom = draw(st.sampled_from([word, f"(a|b)*{word}(a|b)*"])) + if depth == 0 or draw(st.booleans()): + return atom + other = draw(star_free_regexes(depth=depth - 1)) + return draw(st.sampled_from([f"({atom})({other})", f"({atom})|({other})"])) + + +@given(star_free_regexes()) +def test_star_free_constructions_are_star_free(regex): + assert is_star_free(regex, alphabet="ab") + + +@pytest.mark.parametrize("period", [2, 3, 4]) +def test_counting_languages_are_not_star_free(period): + assert not is_star_free("(" + "(a|b)" * period + ")*", alphabet="ab") + + +# -- adapters ----------------------------------------------------------------------- + + +def test_epsilon_machines_convert_through_their_support(): + sg = even_process().transition_semigroup() + assert sg.states == ("A", "B", SINK) + assert len(sg) == 7 + assert not sg.is_aperiodic() + assert golden_mean().transition_semigroup().is_aperiodic() + + +def test_dfa_transition_semigroup(): + dfa = DFA(input_alphabet=frozenset("ab"), 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, "a") + dfa.add_transition(0, 0, "b") + dfa.add_transition(1, 1, "b") + sg = dfa.transition_semigroup() + assert len(sg) == 2 + assert transformation_semigroup(sg) is sg + + +def test_nondeterministic_machines_are_rejected(): + nfa = NFA(input_alphabet=frozenset("a"), initial_states=frozenset({0}), accepting_states=frozenset({0})) + nfa.add_transition(0, 0, "a") + nfa.add_transition(0, 1, "a") + with pytest.raises(ValueError, match="not unifilar"): + transformation_semigroup(nfa) diff --git a/tests/test_algebra_varieties.py b/tests/test_algebra_varieties.py new file mode 100644 index 0000000..8d8a077 --- /dev/null +++ b/tests/test_algebra_varieties.py @@ -0,0 +1,121 @@ +"""Tests for language varieties decided on syntactic monoids.""" + +import itertools + +import pytest +from hypothesis import given +from hypothesis import strategies as st + +from sofic.algebra import factor_monoid, syntactic_monoid +from sofic.algebra.varieties import ( + is_commutative_language, + is_definite, + is_group_language, + is_l_trivial, + is_locally_testable, + is_piecewise_testable, + is_r_trivial, + is_reverse_definite, +) +from sofic.automata.dfa import DFA +from sofic.shifts.sofic import SoficShift + +#: (regex, piecewise testable, locally testable, definite, reverse definite, group, commutative) +KNOWN = [ + ("(a|b)*a", False, True, True, False, False, False), + ("a(a|b)*", False, True, False, True, False, False), + ("(a|b)*ab(a|b)*", True, True, False, False, False, False), + ("(ab)*", False, True, False, False, False, False), + ("a*b*", True, True, False, False, False, False), + ("a*ba*", True, False, False, False, False, True), + ("((a|b)(a|b))*", False, False, False, False, True, True), + ("(a|b)*", True, True, True, True, True, True), +] + + +@pytest.mark.parametrize(("regex", "pt", "lt", "definite", "reverse", "group", "commutative"), KNOWN) +def test_known_languages(regex, pt, lt, definite, reverse, group, commutative): + assert is_piecewise_testable(regex, alphabet="ab") is pt + assert is_locally_testable(regex, alphabet="ab") is lt + assert is_definite(regex, alphabet="ab") is definite + assert is_reverse_definite(regex, alphabet="ab") is reverse + assert is_group_language(regex, alphabet="ab") is group + assert is_commutative_language(regex, alphabet="ab") is commutative + + +def test_green_trivialities_bracket_piecewise_testability(): + assert is_l_trivial("(a|b)*a", alphabet="ab") and not is_r_trivial("(a|b)*a", alphabet="ab") + assert is_r_trivial("a(a|b)*", alphabet="ab") and not is_l_trivial("a(a|b)*", alphabet="ab") + + +@st.composite +def complete_dfas(draw, max_states=4): + n = draw(st.integers(1, max_states)) + dfa = DFA( + input_alphabet=frozenset("ab"), + initial_states=frozenset({0}), + accepting_states=frozenset(draw(st.sets(st.integers(0, n - 1)))), + ) + for q in range(n): + dfa.graph.add_state(q) + for q in range(n): + for symbol in "ab": + dfa.add_transition(q, draw(st.integers(0, n - 1)), symbol) + return dfa + + +def _words(max_length): + for length in range(max_length + 1): + yield from itertools.product("ab", repeat=length) + + +@given(complete_dfas()) +def test_definite_languages_depend_only_on_a_suffix(dfa): + if not is_definite(dfa): + return + k = len(syntactic_monoid(dfa)) + for word in _words(k + 3): + if len(word) >= k: + assert dfa.recognizes(word) == dfa.recognizes(word[-k:]) + + +@given(complete_dfas()) +def test_reverse_definite_languages_depend_only_on_a_prefix(dfa): + if not is_reverse_definite(dfa): + return + k = len(syntactic_monoid(dfa)) + for word in _words(k + 3): + if len(word) >= k: + assert dfa.recognizes(word) == dfa.recognizes(word[:k]) + + +@given(complete_dfas()) +def test_inclusions_between_varieties(dfa): + if is_definite(dfa) or is_reverse_definite(dfa): + assert is_locally_testable(dfa) + if is_piecewise_testable(dfa): + assert is_r_trivial(dfa) and is_l_trivial(dfa) + + +def shift(edges): + s = SoficShift(symbol_alphabet=frozenset(symbol for _, _, symbol in edges)) + for source, target, symbol in edges: + s.add_transition(source, target, symbol) + return s + + +@pytest.mark.parametrize( + "edges", + [ + [("A", "A", "0"), ("A", "B", "1"), ("B", "A", "0")], + [("A", "A", "0"), ("A", "A", "1")], + [("A", "A", "0"), ("A", "B", "1"), ("B", "C", "1"), ("C", "A", "0")], + ], +) +def test_shifts_of_finite_type_have_locally_testable_languages(edges): + assert is_locally_testable(factor_monoid(shift(edges), monoid=False)) + + +def test_the_even_shift_is_not_locally_testable(): + even = shift([("A", "A", "0"), ("A", "B", "1"), ("B", "A", "1")]) + assert not is_locally_testable(factor_monoid(even, monoid=False)) From c559c7ae8d97149395fbcc0008fee696d3d689f9 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Sun, 25 Oct 2026 17:16:43 -0600 Subject: [PATCH 2/2] Keep sofic.algebra usable without sympy Compute holonomy group orders and recognize T_n by closing small permutation groups directly, falling back to sympy only for very large groups. Skip the tests that need sympy group objects when it is not installed. Co-authored-by: Cursor --- sofic/algebra/_groups.py | 24 ++++++++++++++++++++++++ sofic/algebra/complexity.py | 4 ++-- sofic/algebra/holonomy.py | 17 ++++++++++++----- tests/test_algebra_dedeo.py | 1 + tests/test_algebra_holonomy.py | 25 +++++++++++++++---------- tests/test_algebra_random_maps.py | 1 + tests/test_algebra_semigroups.py | 2 ++ 7 files changed, 57 insertions(+), 17 deletions(-) diff --git a/sofic/algebra/_groups.py b/sofic/algebra/_groups.py index cb37c3b..e21cc97 100644 --- a/sofic/algebra/_groups.py +++ b/sofic/algebra/_groups.py @@ -37,6 +37,30 @@ def permutation_group(perms: Iterable[tuple[int, ...]], degree: int) -> Any: return group +def group_order(perms: Iterable[tuple[int, ...]], degree: int, *, cap: int = 1_000_000) -> int: + """Order of the permutation group generated by ``perms``, without sympy when small. + + Enumerates the group by breadth-first closure; only a group larger than + ``cap`` falls back to sympy's Schreier-Sims. + """ + identity = tuple(range(degree)) + gens = [p for p in dict.fromkeys(perms) if p != identity] + seen = {identity} + frontier = [identity] + while frontier: + nxt = [] + for p in frontier: + for g in gens: + q = tuple(g[x] for x in p) + if q not in seen: + if len(seen) >= cap: + return int(permutation_group(gens, degree).order()) + seen.add(q) + nxt.append(q) + frontier = nxt + return len(seen) + + def composition_factor_orders(group: Any) -> list[int]: """Orders of the Jordan-Hölder composition factors, ascending, with multiplicity. diff --git a/sofic/algebra/complexity.py b/sofic/algebra/complexity.py index 853541c..40bb8be 100644 --- a/sofic/algebra/complexity.py +++ b/sofic/algebra/complexity.py @@ -90,7 +90,7 @@ def is_full_transformation_monoid(semigroup: Any) -> bool: """ from math import factorial - from sofic.algebra._groups import permutation_group + from sofic.algebra._groups import group_order n = semigroup.degree gens = list(semigroup.generators.values()) @@ -99,7 +99,7 @@ def is_full_transformation_monoid(semigroup: Any) -> bool: if not any(g.rank() == n - 1 for g in gens): return False perms = [g.images for g in gens if g.is_permutation()] - return bool(perms) and int(permutation_group(perms, n).order()) == factorial(n) + return bool(perms) and group_order(perms, n) == factorial(n) def _full_transformation_bound(hd: HolonomyDecomposition) -> int: diff --git a/sofic/algebra/holonomy.py b/sofic/algebra/holonomy.py index be1315b..f05a637 100644 --- a/sofic/algebra/holonomy.py +++ b/sofic/algebra/holonomy.py @@ -101,12 +101,19 @@ def groups(self) -> tuple[Any, ...]: permutation_group(gens, len(tiles)) for gens, tiles in zip(self.group_generators, self.tiles, strict=True) ) + @property + def orders(self) -> tuple[int, ...]: + """Order of each representative's holonomy group (no sympy needed for small groups).""" + from sofic.algebra._groups import group_order + + return tuple( + group_order(gens, len(tiles)) for gens, tiles in zip(self.group_generators, self.tiles, strict=True) + ) + @property def order(self) -> int: """Order of the level group, the direct product of its holonomy groups.""" - if self.is_trivial: - return 1 - return prod(int(g.order()) for g in self.groups) + return prod(self.orders) @property def is_trivial(self) -> bool: @@ -488,8 +495,8 @@ def show(subset: Subset) -> str: lines = [f"holonomy decomposition: {len(self.image_set)} image sets, {self.n_levels} levels"] for level in self.levels: parts = [] - for rep, tiles, group in zip(level.representatives, level.tiles, level.groups, strict=True): - parts.append(f"{show(rep)} [{len(tiles)} tiles, |H|={group.order()}]") + for rep, tiles, order in zip(level.representatives, level.tiles, level.orders, strict=True): + parts.append(f"{show(rep)} [{len(tiles)} tiles, |H|={order}]") lines.append(f" depth {level.depth}: " + "; ".join(parts)) return "\n".join(lines) diff --git a/tests/test_algebra_dedeo.py b/tests/test_algebra_dedeo.py index 4b5c292..f7642a0 100644 --- a/tests/test_algebra_dedeo.py +++ b/tests/test_algebra_dedeo.py @@ -89,6 +89,7 @@ def test_three_bit_left_shift_powerset_complexity_without_enumeration(): def test_table_one_three_bit_rows(): + pytest.importorskip("sympy") one_letter = TransformationSemigroup({"a": register(3, lambda s: s[0] ^ s[1])}) assert prime_decomposition(one_letter).simple_group_orders() == (7,) two_letters = TransformationSemigroup( diff --git a/tests/test_algebra_holonomy.py b/tests/test_algebra_holonomy.py index 782e203..030f8fe 100644 --- a/tests/test_algebra_holonomy.py +++ b/tests/test_algebra_holonomy.py @@ -5,7 +5,6 @@ import pytest from hypothesis import given from hypothesis import strategies as st -from sympy.combinatorics.named_groups import AlternatingGroup, DirectProduct, SymmetricGroup from sofic.algebra import ( TransformationSemigroup, @@ -76,7 +75,7 @@ def test_egri_nagy_nehaniv_example_7_2(): hd = holonomy_decomposition(TransformationSemigroup(EXAMPLE_7_2)) assert len(hd.image_set) == 19 assert hd.n_levels == 5 - groups = [sorted(int(g.order()) for g in level.groups) for level in hd.levels] + groups = [sorted(level.orders) for level in hd.levels] assert groups == [[1], [3], [2, 6], [2, 3], [1]] @@ -195,6 +194,7 @@ def test_full_transformation_monoid_complexity(n): def test_cyclic_groups_have_complexity_one(n): sg = TransformationSemigroup({"a": (*range(1, n), 0)}) assert krohn_rhodes_complexity(sg).exact == 1 + pytest.importorskip("sympy") expected = [p for p in range(2, n + 1) if n % p == 0 and all(p % q for q in range(2, p))] assert set(prime_decomposition(sg).simple_group_orders()) == set(expected) @@ -224,6 +224,7 @@ def test_complexity_bounds_are_consistent(sg): def test_a5_is_a_nonabelian_simple_prime(): + pytest.importorskip("sympy") sg = TransformationSemigroup({"a": (1, 2, 3, 4, 0), "b": (1, 2, 0, 3, 4)}) (factor,) = prime_decomposition(sg).groups assert (factor.kind, factor.order) == ("nonabelian simple", 60) @@ -231,17 +232,21 @@ def test_a5_is_a_nonabelian_simple_prime(): @pytest.mark.parametrize( - ("group", "orders"), + ("factors", "orders"), [ - (AlternatingGroup(5), [60]), - (SymmetricGroup(5), [2, 60]), - (SymmetricGroup(6), [2, 360]), - (SymmetricGroup(4), [2, 2, 2, 3]), - (DirectProduct(AlternatingGroup(5), AlternatingGroup(5)), [60, 60]), - (DirectProduct(AlternatingGroup(5), SymmetricGroup(3)), [2, 3, 60]), + ((("A", 5),), [60]), + ((("S", 5),), [2, 60]), + ((("S", 6),), [2, 360]), + ((("S", 4),), [2, 2, 2, 3]), + ((("A", 5), ("A", 5)), [60, 60]), + ((("A", 5), ("S", 3)), [2, 3, 60]), ], ) -def test_composition_factor_orders(group, orders): +def test_composition_factor_orders(factors, orders): + named = pytest.importorskip("sympy.combinatorics.named_groups") + build = {"A": named.AlternatingGroup, "S": named.SymmetricGroup} + groups = [build[kind](n) for kind, n in factors] + group = groups[0] if len(groups) == 1 else named.DirectProduct(*groups) assert composition_factor_orders(group) == orders diff --git a/tests/test_algebra_random_maps.py b/tests/test_algebra_random_maps.py index b6d6920..7288a4c 100644 --- a/tests/test_algebra_random_maps.py +++ b/tests/test_algebra_random_maps.py @@ -151,6 +151,7 @@ def test_synchronizing_iff_minimal_rank_one(sg): @given(semigroups()) def test_schutzenberger_group_order_is_h_class_size(sg): + pytest.importorskip("sympy") green = sg.green() for members in green.h_classes: assert int(green.schutzenberger_group(members[0]).order()) == len(members) diff --git a/tests/test_algebra_semigroups.py b/tests/test_algebra_semigroups.py index 6748c5e..1fc3d11 100644 --- a/tests/test_algebra_semigroups.py +++ b/tests/test_algebra_semigroups.py @@ -142,6 +142,7 @@ def test_cayley_graph_edges_are_labeled(): def test_green_structure_of_t3(): + pytest.importorskip("sympy") green = full_transformation_monoid(3).green() assert sorted(d.size for d in green.d_classes) == [3, 6, 18] assert sorted(int(g.order()) for _, g in green.maximal_subgroups()) == [1, 2, 6] @@ -168,6 +169,7 @@ def test_aperiodicity_agrees_with_trivial_h_classes(sg): @given(semigroups()) def test_maximal_subgroup_orders_match_their_h_classes(sg): + pytest.importorskip("sympy") green = sg.green() for e, group in green.maximal_subgroups(): assert int(group.order()) == len(green.h_classes[green.h_of[e]])