From ab0c000eef5a3ca5b8abcad7a3263441ee44116e Mon Sep 17 00:00:00 2001 From: Ryan James Date: Thu, 17 Sep 2026 18:02:31 -0600 Subject: [PATCH 1/5] Add co-lexicographic (Wheeler) ordering Wheeler order is the total order on states in which the words reaching each state form an interval of the co-lexicographically sorted prefixes of the language. Co-lex compares words from the last symbol backwards, so a Wheeler presentation is one whose states partition history space into bands of recency. Topological core in automata/wheeler.py: axiom checking, order search by partition refinement then linearization of the maximum co-lex relation, co-lex width, minimum WDFA, and WNFA determinization. Every graph-backed model inherits is_wheeler/wheeler_order/colex_width from StateMachine. Path coherence gives synchronization.py an interval power automaton: n(n+1)/2 states instead of 2^n, so the Markov order, cryptic order, and reset threshold stay polynomial. automata/wheeler_index.py stores the GMS arrays for O(|w| log |A|) membership and co-lex rank/unrank. shifts/wheeler.py searches higher block presentations for a Wheeler cover; generators/wheeler_epsilon.py proposes C_W, the entropy of the states of a Wheeler presentation, bounded below by both C_mu and H[X_0]. That layer is original and every docstring says no source exists. Wheeler is not finite memory: wheeler_infinite_order_process is Wheeler with infinite Markov order, while the even process is non-Wheeler in every presentation because Wheeler languages are star-free. --- docs/.DS_Store | Bin 0 -> 8196 bytes docs/automata/automata.rst | 1 + docs/automata/wheeler.rst | 148 ++++++ docs/generators/generators.rst | 1 + docs/generators/wheeler_epsilon.rst | 103 ++++ docs/references.bib | 89 ++++ docs/shifts/shifts.rst | 1 + docs/shifts/wheeler.rst | 65 +++ sofic/.DS_Store | Bin 0 -> 8196 bytes sofic/__init__.py | 20 + sofic/automata/__init__.py | 29 ++ sofic/automata/wheeler.py | 697 ++++++++++++++++++++++++++++ sofic/automata/wheeler_index.py | 291 ++++++++++++ sofic/base.py | 36 ++ sofic/examples/__init__.py | 2 + sofic/examples/epsilon_machines.py | 44 ++ sofic/generators/__init__.py | 16 + sofic/generators/epsilon_machine.py | 26 ++ sofic/generators/synchronization.py | 76 ++- sofic/generators/wheeler_epsilon.py | 247 ++++++++++ sofic/serialization.py | 9 +- sofic/shifts/__init__.py | 16 + sofic/shifts/covers.py | 16 + sofic/shifts/wheeler.py | 242 ++++++++++ tests/test_wheeler.py | 542 +++++++++++++++++++++ 25 files changed, 2715 insertions(+), 2 deletions(-) create mode 100644 docs/.DS_Store create mode 100644 docs/automata/wheeler.rst create mode 100644 docs/generators/wheeler_epsilon.rst create mode 100644 docs/shifts/wheeler.rst create mode 100644 sofic/.DS_Store create mode 100644 sofic/automata/wheeler.py create mode 100644 sofic/automata/wheeler_index.py create mode 100644 sofic/generators/wheeler_epsilon.py create mode 100644 sofic/shifts/wheeler.py create mode 100644 tests/test_wheeler.py diff --git a/docs/.DS_Store b/docs/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..3f31ffd658bdf9d36fe56b49bf56c4e0b66df78d GIT binary patch literal 8196 zcmeI1y=xRf7{=ex_-L6my!>gR(=FJz^*` z9Ql~#h>4AviDoz{Gn|xNS=kec(yL=!tmUK<6Qvag#DSs%e0HzVW9ra^w!`zg`u=bw z%n#F6t20PD=&#?~e|zO!V|u%&>CM%Qw(|D`>gOS*pmX5Zr9En+Zd2EPld$l*sBl>h zp4A>q|9oYZV4h)`(^RWm9(xbF0G>VS(~w4RhN=3r>r<&mSg%y)vfXOL_bpq6HYc2? zE#a(#Z34O!q)gfN8K3mPSkSn9@4x)IcyL+oV~H+jJ5J|unY(aAT`-NQ@AE8v1J5ky za(q~AHm^sT!^!9Cm*u^6X&qBzKQ^&bJE&vr2Qviw_WNU%&gFUeJ5)$jjL(pQ%6A`XZHfdeXPt+n`vzCBw9X>-bt(YtOltjkal^$;}r mIF1;x1abC{e5t(WE$JMb5(+DRt> literal 0 HcmV?d00001 diff --git a/docs/automata/automata.rst b/docs/automata/automata.rst index dd380af..7c38409 100644 --- a/docs/automata/automata.rst +++ b/docs/automata/automata.rst @@ -46,6 +46,7 @@ Canonical structures atomaton observation_table learning + wheeler Advanced ======== diff --git a/docs/automata/wheeler.rst b/docs/automata/wheeler.rst new file mode 100644 index 0000000..e957b90 --- /dev/null +++ b/docs/automata/wheeler.rst @@ -0,0 +1,148 @@ +.. wheeler.rst +.. py:module:: sofic.automata.wheeler + +******************************** +Co-lexicographic (Wheeler) order +******************************** + +A labeled graph is *Wheeler* when its states admit a total order in which the +states with no incoming edges come first and, for edges ``(u, v)`` labeled +``a`` and ``(u', v')`` labeled ``a'``, + +* ``a < a'`` implies ``v < v'``, and +* ``a == a'`` and ``u < u'`` imply ``v <= v'`` + +:cite:`Gagie2017`. Equivalently, the words reaching each state form an +interval of the co-lexicographically sorted prefixes of the language +:cite:`Alanko2020`. Co-lex order compares words from the last symbol backwards, +so a Wheeler presentation is one whose states partition history space into +*bands of recency* rather than arbitrary sets. + +Wheeler orders are the width-one case of the co-lexicographic partial orders of +:cite:`CotumaccioPrezza2021`. :func:`colex_width` measures how far a machine is +from being Wheeler; the width bounds the cost of indexing it, storing it, and +determinizing it. + +Presentations versus languages +============================== + +:func:`is_wheeler` asks whether *this presentation* is Wheeler, which is a +sorting question answerable in polynomial time. Whether the *language* is +Wheeler — whether any equivalent automaton is — is much harder: ``O(mn)`` for a +DFA :cite:`Becker2023`, PSPACE-complete for an NFA :cite:`Alanko2021`, and +NP-hard to recognize as a graph property :cite:`GibneyThankachan2019`. +:func:`~sofic.shifts.wheeler.wheeler_cover` and +:func:`~sofic.generators.wheeler_epsilon.wheeler_presentation` search over +presentations for the process-level question. + +Wheeler is not a restatement of finite memory. The golden mean process is +Wheeler; so is +:func:`~sofic.examples.epsilon_machines.wheeler_infinite_order_process`, whose +Markov order is infinite. Conversely the even process is not Wheeler in any +presentation, because Wheeler languages are star-free +:cite:`ShyrThierrin1974` :cite:`Alanko2021` and the even process counts +``1``\ s modulo two. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean, even_process, wheeler_infinite_order_process + + In [2]: golden_mean().wheeler_order().states + + In [3]: even_process().is_wheeler(), even_process().colex_width() + + In [4]: machine = wheeler_infinite_order_process(); machine.is_wheeler(), machine.markov_order() + +Every graph-backed model inherits :meth:`~sofic.base.StateMachine.is_wheeler`, +:meth:`~sofic.base.StateMachine.wheeler_order`, and +:meth:`~sofic.base.StateMachine.colex_width`, so automata, shifts, and +ε-machines all answer the same questions. + +What the order buys +=================== + +* A canonical state numbering, hence canonical transition matrices and + ``O(m)`` equality by comparing Burrows-Wheeler strings + (:func:`wheeler_canonical_form`, :func:`wheeler_isomorphic`). +* Path coherence: the states reachable from an interval on a given symbol are + again an interval. :func:`~sofic.generators.synchronization.power_automaton` + uses this to replace ``2^n`` subsets with ``n(n+1)/2`` intervals, which makes + the Markov order, cryptic order, and reset threshold polynomial. +* Determinization to at most ``2n - 1 - |Sigma|`` states + (:func:`wnfa_to_wdfa`) and a unique minimal WDFA (:func:`minimum_wdfa`), + both of which fail for general automata. +* A succinct index — see below. + +Order +===== + +.. autoclass:: WheelerOrder + :members: validate +.. autoclass:: LabeledGraph +.. autoexception:: WheelerError + +.. autofunction:: labeled_graph +.. autofunction:: wheeler_order +.. autofunction:: wheeler_order_of_graph +.. autofunction:: is_wheeler +.. autofunction:: is_input_consistent +.. autofunction:: check_wheeler_axioms + +Width +===== + +.. autofunction:: colex_width +.. autofunction:: maximum_colex_relation +.. autofunction:: maximum_colex_relation_of_graph + +Canonical forms and minimization +================================ + +.. autofunction:: minimum_wdfa +.. autofunction:: wnfa_to_wdfa +.. autofunction:: wheeler_canonical_form +.. autofunction:: wheeler_isomorphic +.. autofunction:: wheeler_state_index + +Burrows-Wheeler index +===================== + +.. py:currentmodule:: sofic.automata.wheeler_index + +:class:`WheelerIndex` stores a Wheeler machine as the out-degree, in-degree, +and label arrays of :cite:`Gagie2017`. Because the Wheeler axioms make the +edges sorted by ``(label, source)`` coincide with the edges sorted by target, +following a symbol maps one node interval onto another — FM-index backward +search, generalized to labeled graphs. Rank is served by binary search over +per-symbol position arrays, which costs a logarithmic factor but adds no +dependency beyond numpy. + +A shift presents its factor language with every state both initial and +accepting, the case :cite:`Gagie2017` Theorem 6 covers explicitly, so the same +index answers membership and enumeration queries for shifts. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean + + In [2]: from sofic.automata.wheeler_index import WheelerIndex + + In [3]: index = WheelerIndex.from_model(golden_mean()); index.bits() + + In [4]: index.contains((0, 1, 0)), index.contains((1, 1)) + + In [5]: list(index.words_of_length(3)) + + In [6]: index.unrank_word(index.rank_word((0, 1, 0)), 3) + +Words are listed in co-lexicographic order, so :meth:`WheelerIndex.rank_word` +and :meth:`WheelerIndex.unrank_word` invert one another and +:meth:`WheelerIndex.sample_word` draws uniformly from the words of a length +without enumerating them. + +.. autoclass:: WheelerIndex + :members: from_model, bits, step, forward_search, contains, count_states, + states_reached, count_words, words_of_length, rank_word, + unrank_word, sample_word + +.. autofunction:: wheeler_index diff --git a/docs/generators/generators.rst b/docs/generators/generators.rst index 468463e..7f66ba2 100644 --- a/docs/generators/generators.rst +++ b/docs/generators/generators.rst @@ -38,6 +38,7 @@ Computational mechanics generative_models directional_flow alternative_complexity + wheeler_epsilon Constructions and conversions ============================= diff --git a/docs/generators/wheeler_epsilon.rst b/docs/generators/wheeler_epsilon.rst new file mode 100644 index 0000000..f9f91db --- /dev/null +++ b/docs/generators/wheeler_epsilon.rst @@ -0,0 +1,103 @@ +.. wheeler_epsilon.rst +.. py:module:: sofic.generators.wheeler_epsilon + +********************************** +Co-lexicographic epsilon-machines +********************************** + +.. warning:: + + Everything on this page is original and uncited. Wheeler automata are a + purely topological theory :cite:`Gagie2017` :cite:`Alanko2020`; a literature + search turned up no treatment of weighted, probabilistic, or + information-theoretic Wheeler automata, so no canonical source exists for + the quantities defined here. They are proposals, and every docstring says + so. + +Co-lexicographic order compares words from the last symbol backwards, which is +the *recency* order on pasts. A presentation is Wheeler exactly when its +partition of history space is an interval partition under recency: each state +owns a contiguous band of pasts rather than an arbitrary set. Two consequences +follow, one measure-theoretic and one information-theoretic. + +The stationary distribution becomes a genuine cumulative distribution over +pasts. Under an arbitrary state numbering the partial sums of the stationary +vector mean nothing; in Wheeler order they sweep history space from the most +remote pasts to the most recent, which is what arithmetic coding over pasts +needs. :func:`colex_cdf` returns that sweep, :func:`cylinder_measure` the mass +of a rank interval, and :func:`word_cylinder_measure` the mass of the interval +a word selects — pairing directly with +:meth:`~sofic.automata.wheeler_index.WheelerIndex.forward_search`. + +The constraint has a price in bits. A process whose causal states are not +recency intervals must split them to get one, and the split states cost +entropy. :func:`wheeler_statistical_complexity` measures the result and +:func:`wheeler_complexity_gap` the excess over :math:`C_\mu`. + +.. math:: + + C_W = \operatorname{H}[\text{Wheeler states}] \ge \max(C_\mu, \operatorname{H}[X_0]), + +the first bound because every Wheeler presentation refines the causal-state +partition, with equality exactly when the ε-machine is already Wheeler; the +second because a Wheeler presentation is input consistent, so its state +determines the symbol that entered it. The fair coin is the extreme case of the +second bound: nothing at all to predict, :math:`C_\mu = 0`, yet sortability +costs a full bit of memory for a symbol the process will never reuse. + +.. ipython:: + + In [1]: from sofic.examples import golden_mean, fair_coin + + In [2]: from sofic.generators.wheeler_epsilon import colex_cdf, wheeler_complexity_gap + + In [3]: machine = golden_mean(); machine.is_wheeler() + + In [4]: colex_cdf(machine) + + In [5]: machine.wheeler_statistical_complexity(), machine.statistical_complexity() + + In [6]: coin = fair_coin(); coin.is_wheeler(), coin.statistical_complexity() + + In [7]: coin.wheeler_statistical_complexity(), wheeler_complexity_gap(coin) + +Finding a presentation +====================== + +:func:`wheeler_presentation` returns the machine itself when it is already +Wheeler. Otherwise, when the Markov order ``R`` is finite, it returns the +order-``R`` de Bruijn presentation, whose states *are* the length-``R`` words +and so sort co-lexicographically by construction — every finite-order process +is therefore a Wheeler *language* even when its ε-machine is not a Wheeler +*presentation*. Failing that it tries edge-machine refinements, whose order-``k`` +states are length-``k`` transition paths and are therefore entered on a single +symbol. + +The search can fail for good reason: Wheeler languages are star-free +:cite:`Alanko2021`, so the even process has no Wheeler presentation at any +order and :func:`wheeler_presentation` raises +:class:`~sofic.automata.wheeler.WheelerError`. + +The converse separation also holds. +:func:`~sofic.examples.epsilon_machines.wheeler_infinite_order_process` is a +five-state Wheeler ε-machine of infinite Markov order, so Wheelerness is not a +restatement of finite memory in either direction. + +API +=== + +.. autofunction:: wheeler_presentation +.. autofunction:: debruijn_presentation +.. autofunction:: wheeler_statistical_complexity +.. autofunction:: wheeler_complexity_gap +.. autofunction:: colex_cdf +.. autofunction:: cylinder_measure +.. autofunction:: word_cylinder_measure + +:class:`~sofic.generators.epsilon_machine.EpsilonMachine` also exposes +:meth:`~sofic.generators.epsilon_machine.EpsilonMachine.wheeler_presentation` +and +:meth:`~sofic.generators.epsilon_machine.EpsilonMachine.wheeler_statistical_complexity`, +alongside the :meth:`~sofic.base.StateMachine.is_wheeler`, +:meth:`~sofic.base.StateMachine.wheeler_order`, and +:meth:`~sofic.base.StateMachine.colex_width` methods every model inherits. diff --git a/docs/references.bib b/docs/references.bib index a794925..44fb110 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -901,3 +901,92 @@ @inproceedings{Volkov2008 year = {2008}, doi = {10.1007/978-3-540-88282-4_4}, } + +@article{Gagie2017, + author = {Gagie, Travis and Manzini, Giovanni and Sir{\'e}n, Jouni}, + title = {Wheeler Graphs: A Framework for {BWT}-Based Data Structures}, + journal = {Theoretical Computer Science}, + volume = {698}, + pages = {67--78}, + year = {2017}, + doi = {10.1016/j.tcs.2017.06.016}, +} + +@inproceedings{Alanko2020, + author = {Alanko, Jarno and D'Agostino, Giovanna and Policriti, Alberto and Prezza, Nicola}, + title = {Regular Languages Meet Prefix Sorting}, + booktitle = {Proceedings of the Thirty-First Annual {ACM-SIAM} Symposium on Discrete Algorithms ({SODA})}, + pages = {911--930}, + publisher = {SIAM}, + year = {2020}, + doi = {10.1137/1.9781611975994.55}, +} + +@article{Alanko2021, + author = {Alanko, Jarno and D'Agostino, Giovanna and Policriti, Alberto and Prezza, Nicola}, + title = {Wheeler Languages}, + journal = {Information and Computation}, + volume = {281}, + pages = {104820}, + year = {2021}, + doi = {10.1016/j.ic.2021.104820}, +} + +@inproceedings{CotumaccioPrezza2021, + author = {Cotumaccio, Nicola and Prezza, Nicola}, + title = {On Indexing and Compressing Finite Automata}, + booktitle = {Proceedings of the Thirty-Second Annual {ACM-SIAM} Symposium on Discrete Algorithms ({SODA})}, + pages = {2585--2599}, + publisher = {SIAM}, + year = {2021}, + doi = {10.1137/1.9781611976465.153}, +} + +@article{Cotumaccio2023, + author = {Cotumaccio, Nicola and D'Agostino, Giovanna and Policriti, Alberto and Prezza, Nicola}, + title = {Co-lexicographically Ordering Automata and Regular Languages --- Part {I}}, + journal = {Journal of the {ACM}}, + volume = {70}, + number = {4}, + pages = {27:1--27:73}, + year = {2023}, + doi = {10.1145/3607471}, + eprint = {2208.04931}, + archivePrefix = {arXiv}, +} + +@inproceedings{GibneyThankachan2019, + author = {Gibney, Daniel and Thankachan, Sharma V.}, + title = {On the Hardness and Inapproximability of Recognizing Wheeler Graphs}, + booktitle = {27th Annual European Symposium on Algorithms ({ESA})}, + series = {Leibniz International Proceedings in Informatics ({LIPIcs})}, + volume = {144}, + pages = {51:1--51:16}, + publisher = {Schloss Dagstuhl--Leibniz-Zentrum f{\"u}r Informatik}, + year = {2019}, + doi = {10.4230/LIPIcs.ESA.2019.51}, +} + +@inproceedings{Becker2023, + author = {Becker, Ruben and Cenzato, Davide and Kim, Sung-Hwan and Kodric, Bojana and Policriti, Alberto and Prezza, Nicola}, + title = {Optimal {W}heeler Language Recognition}, + booktitle = {String Processing and Information Retrieval ({SPIRE})}, + series = {Lecture Notes in Computer Science}, + volume = {14240}, + pages = {62--74}, + publisher = {Springer}, + year = {2023}, + doi = {10.1007/978-3-031-43980-3_6}, + eprint = {2306.04737}, + archivePrefix = {arXiv}, +} + +@article{ShyrThierrin1974, + author = {Shyr, H. J. and Thierrin, Gabriel}, + title = {Ordered Automata and Associated Languages}, + journal = {Tamkang Journal of Mathematics}, + volume = {5}, + number = {1}, + pages = {9--20}, + year = {1974}, +} diff --git a/docs/shifts/shifts.rst b/docs/shifts/shifts.rst index f83b897..2f07371 100644 --- a/docs/shifts/shifts.rst +++ b/docs/shifts/shifts.rst @@ -22,3 +22,4 @@ type, Sofic shifts, Dyck shifts, topological Markov chains, and covers textile dyck_enumeration covers + wheeler diff --git a/docs/shifts/wheeler.rst b/docs/shifts/wheeler.rst new file mode 100644 index 0000000..37f3905 --- /dev/null +++ b/docs/shifts/wheeler.rst @@ -0,0 +1,65 @@ +.. wheeler.rst +.. py:module:: sofic.shifts.wheeler + +********************** +Wheeler presentations +********************** + +A shift presents its factor language with every state both initial and +accepting, which is the "all states initial" case of :cite:`Gagie2017` +Theorem 6. The minimal right-resolving presentation is rarely Wheeler as it +stands — a state entered on two different symbols already breaks the axioms — +but remembering the last few symbols restores input consistency, and for many +shifts that is enough. + +:func:`wheeler_cover` right-resolves the presentation, walks up the higher +block presentations until one is Wheeler, then merges the Wheeler-consecutive +states that share an incoming label and a follower language. The result is the +smallest Wheeler presentation *reachable that way*, not provably the smallest +Wheeler presentation of the language. + +Not every sofic shift has one at any order. Wheeler languages are star-free +:cite:`ShyrThierrin1974` :cite:`Alanko2021`, so a shift whose syntactic monoid +contains a nontrivial group — the even shift being the standard example — is +excluded outright, and :func:`wheeler_cover` raises +:class:`~sofic.automata.wheeler.WheelerError`. + +.. ipython:: + + In [1]: from sofic.shifts.sofic import SoficShift + + In [2]: from sofic.shifts.wheeler import wheeler_cover, wheeler_order_of_shift + + In [3]: from sofic.graph import ATTR_SYMBOL + + In [4]: shift = SoficShift(symbol_alphabet=frozenset({0, 1})) + + In [5]: for source, symbol, target in ((0, 0, 0), (0, 1, 1), (1, 0, 0)): + ...: shift.graph.add_transition(source, target, **{ATTR_SYMBOL: symbol}) + + In [6]: wheeler_order_of_shift(shift) + + In [7]: cover = wheeler_cover(shift); cover.wheeler_order().states + +The cover is a :class:`~sofic.shifts.covers.WheelerCover`, a +:class:`~sofic.shifts.sofic.SoficShift` subclass that sits beside the Fischer +and Krieger covers. It does not fill the Krieger stubs in +:mod:`sofic.shifts.cover_construction`: a Krieger cover's states are *all* +sets of pasts closed under the follower relation, whereas a Wheeler cover +carries only those that happen to be recency intervals. + +Once a shift has a Wheeler cover, :func:`wheeler_index_of_shift` gives +``O(|w| log |A|)`` factor-language membership in place of scanning +:meth:`~sofic.shifts.base.SymbolicModel.factor_language`. + +API +=== + +.. autofunction:: wheeler_cover +.. autofunction:: is_wheeler_shift +.. autofunction:: wheeler_order_of_shift +.. autofunction:: wheeler_index_of_shift +.. autofunction:: higher_block_presentation +.. autofunction:: right_resolve + +.. autoclass:: sofic.shifts.covers.WheelerCover diff --git a/sofic/.DS_Store b/sofic/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..055b6da0147b3cbe585cf94b5ee211b2de1295e3 GIT binary patch literal 8196 zcmeI1L2DC16vy9~T^p^S9!#qstcRcn5p(jomh=`~-Ry z@A?hAh!XV zKpgll9pL9fgw8RxF*8v=Ixy%G05QO6-EfRLK<7j@#x`aq>Qt<0st09Ll`S!pNymJ| za>Uri%tVt;%A}LBJ1bkED7`!KA}uEsn<%3=APy8Ah;V8i)q0CY^w_K4Tf=^1rJE)r zom6>nr(Bt&&F0!LU4x?f>3jXz%i!JJS?cw_P>%)=5(`k9@mx3aP{aBh!S`?#5a&lk1MS2Xs!v3%XRVR+rj=J zS3#F^&C+>X;U06Hr9Noxqd(MsrHDS;wf7&&bS~e)<(=*G3);7-Yr5Y3N$YXee=+A7 z+5qJ}&B?r7>UjleVkUa4?g5Br&jJyAf+S9{($mY4%q&VyC7k<@$j9_JX|FiSb#DU{GP%){CT>meZ zfB!%J29hZbhyy3y0Tnb?noGd6wRHiQYwZZVjn0jAnTa|D155b4F2@m%{$YsY$j-$! UW+q|=O+Ex@8)Og%{;31M0dxr*OaK4? literal 0 HcmV?d00001 diff --git a/sofic/__init__.py b/sofic/__init__.py index 6a176f2..9d50338 100644 --- a/sofic/__init__.py +++ b/sofic/__init__.py @@ -33,17 +33,23 @@ UnifilarAutomaton, VisiblyPushdownAutomaton, WeightedFiniteStateTransducer, + WheelerError, + WheelerIndex, + WheelerOrder, automaton_to_regex, cartesian_product_gg, cartesian_product_tt, + colex_width, complete, compose_tg, compose_tt, determinize, equivalent, + is_wheeler, minimize, transduce_generator, trim, + wheeler_order, ) from sofic.core import EPSILON, StateIndex, StateMachine, Transition, TransitionGraph from sofic.generators import ( @@ -75,6 +81,8 @@ is_lumpable, lump, minimal_generative_model, + wheeler_presentation, + wheeler_statistical_complexity, wyner_generative_model, ) from sofic.operations import reverse @@ -93,6 +101,8 @@ SymbolicModel, TextileSystem, TopologicalMarkovChain, + WheelerCover, + wheeler_cover, ) __all__ = [ @@ -158,15 +168,21 @@ "UnifilarAutomaton", "VisiblyPushdownAutomaton", "WeightedFiniteStateTransducer", + "WheelerCover", + "WheelerError", + "WheelerIndex", + "WheelerOrder", "WynerGenerativeModel", "automaton_to_regex", "cartesian_product_gg", "cartesian_product_tt", + "colex_width", "complete", "compose_tg", "compose_tt", "determinize", "equivalent", + "is_wheeler", "minimize", "from_yaml", "functional_generative_model", @@ -180,6 +196,10 @@ "reverse", "transduce_generator", "trim", + "wheeler_cover", + "wheeler_order", + "wheeler_presentation", + "wheeler_statistical_complexity", "wyner_generative_model", "__version__", ] diff --git a/sofic/automata/__init__.py b/sofic/automata/__init__.py index 6f1f822..720bf4e 100644 --- a/sofic/automata/__init__.py +++ b/sofic/automata/__init__.py @@ -106,6 +106,21 @@ kleene_star_vpa, union_vpa, ) +from sofic.automata.wheeler import ( + WheelerError, + WheelerOrder, + colex_width, + is_input_consistent, + is_wheeler, + maximum_colex_relation, + minimum_wdfa, + wheeler_canonical_form, + wheeler_isomorphic, + wheeler_order, + wheeler_state_index, + wnfa_to_wdfa, +) +from sofic.automata.wheeler_index import WheelerIndex, wheeler_index __all__ = [ "Atomaton", @@ -150,8 +165,12 @@ "Transducer", "UnifilarAutomaton", "WeightedFiniteStateTransducer", + "WheelerError", + "WheelerIndex", + "WheelerOrder", "VisiblyPushdownAutomaton", "automaton_to_regex", + "colex_width", "cartesian_product_gg", "cartesian_product_tt", "complement_vpa", @@ -172,7 +191,9 @@ "flags_from_string", "icdfa_string_to_dfa", "intersection_vpa", + "is_input_consistent", "is_well_matched", + "is_wheeler", "iter_icdfa", "iter_icdfa_empty_strings", "iter_idfa_strings", @@ -189,7 +210,9 @@ "learn_mealy_lstar", "learn_pfa_alergia", "learn_sofic_dyck_shift_papni", + "maximum_colex_relation", "minimize", + "minimum_wdfa", "next_flags", "next_icdfa_empty_string", "papni_encode", @@ -204,4 +227,10 @@ "union_vpa", "validate_idfa_string", "validate_icdfa_empty_string", + "wheeler_canonical_form", + "wheeler_index", + "wheeler_isomorphic", + "wheeler_order", + "wheeler_state_index", + "wnfa_to_wdfa", ] diff --git a/sofic/automata/wheeler.py b/sofic/automata/wheeler.py new file mode 100644 index 0000000..7427cc2 --- /dev/null +++ b/sofic/automata/wheeler.py @@ -0,0 +1,697 @@ +"""Co-lexicographic (Wheeler) ordering of labeled state machines. + +A labeled graph is *Wheeler* when its states admit a total order in which +states with no incoming edges come first and, for edges ``(u, v)`` labeled +``a`` and ``(u', v')`` labeled ``a'``, ``a < a'`` implies ``v < v'`` while +``a == a'`` and ``u < u'`` imply ``v <= v'`` :cite:`Gagie2017`. Equivalently, +the set of words reaching each state is an interval of the co-lexicographically +sorted prefixes of the recognized language :cite:`Alanko2020`. + +Wheeler orders are the width-one case of the co-lexicographic partial orders of +:cite:`CotumaccioPrezza2021`; :func:`colex_width` measures how far an arbitrary +machine is from being Wheeler. +""" + +from __future__ import annotations + +import itertools +import math +from collections import defaultdict, deque +from collections.abc import Callable, Hashable, Iterable, Sequence +from dataclasses import dataclass +from functools import partial +from typing import Any + +import networkx as nx + +from sofic.exceptions import SoficValidationError +from sofic.graph import ATTR_EMISSION, ATTR_SYMBOL, EPSILON + +#: Cap on the number of tie-breaking permutations tried by :func:`wheeler_order`. +MAX_TIE_PERMUTATIONS = 20_000 + +LabeledEdge = tuple[Hashable, Any, Hashable] + + +class WheelerError(SoficValidationError): + """Raised when a Wheeler invariant is violated or cannot be established.""" + + +@dataclass(frozen=True, slots=True) +class LabeledGraph: + """Symbol-labeled view of a model, shared by every Wheeler routine. + + ``alphabet`` is ordered; ``edges`` are ``(source, symbol, target)`` triples. + """ + + states: tuple[Hashable, ...] + alphabet: tuple[Any, ...] + edges: tuple[LabeledEdge, ...] + + def symbol_rank(self) -> dict[Any, int]: + return {symbol: index for index, symbol in enumerate(self.alphabet)} + + def in_labels(self) -> dict[Hashable, frozenset[Any]]: + """Map each state to the set of labels on its incoming edges.""" + labels: dict[Hashable, set[Any]] = {state: set() for state in self.states} + for _source, symbol, target in self.edges: + labels[target].add(symbol) + return {state: frozenset(symbols) for state, symbols in labels.items()} + + def predecessors(self) -> dict[Hashable, tuple[Hashable, ...]]: + """Map each state to the sources of its incoming edges.""" + preds: dict[Hashable, list[Hashable]] = {state: [] for state in self.states} + for source, _symbol, target in self.edges: + preds[target].append(source) + return {state: tuple(sources) for state, sources in preds.items()} + + +def labeled_graph(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> LabeledGraph: + """Extract the symbol-labeled graph of ``model``. + + Reads :data:`~sofic.graph.ATTR_SYMBOL` for automata and shift presentations + and :data:`~sofic.graph.ATTR_EMISSION` for Mealy HMMs and epsilon-machines. + Epsilon transitions have no Wheeler semantics and raise. + """ + key = repr if symbol_key is None else symbol_key + edges: list[LabeledEdge] = [] + symbols: set[Any] = set() + for transition in model.transitions(): + symbol = transition.data.get(ATTR_SYMBOL, transition.data.get(ATTR_EMISSION)) + if symbol is EPSILON: + raise WheelerError("epsilon transitions have no Wheeler order; remove them first") + if symbol is None: + continue + edges.append((transition.source, symbol, transition.target)) + symbols.add(symbol) + + for attribute in ("input_alphabet", "symbol_alphabet", "observation_alphabet"): + declared = getattr(model, attribute, None) + if declared: + symbols.update(declared) + + return LabeledGraph( + states=tuple(model.states()), + alphabet=tuple(sorted(symbols, key=key)), + edges=tuple(edges), + ) + + +@dataclass(frozen=True, slots=True) +class WheelerOrder: + """A total order on states witnessing the Wheeler property.""" + + states: tuple[Hashable, ...] + rank: dict[Hashable, int] + + @classmethod + def from_sequence(cls, states: Sequence[Hashable]) -> WheelerOrder: + ordered = tuple(states) + return cls(states=ordered, rank={state: index for index, state in enumerate(ordered)}) + + def __len__(self) -> int: + return len(self.states) + + def __iter__(self) -> Any: + return iter(self.states) + + def validate(self, graph: LabeledGraph) -> None: + """Raise :class:`WheelerError` unless the axioms hold for ``graph``.""" + violation = _first_axiom_violation(graph, self.rank) + if violation is not None: + raise WheelerError(violation) + + +def is_input_consistent(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> bool: + """Return whether every state's incoming edges share a single label. + + A necessary condition for the Wheeler property: axiom one forces states + with distinct incoming labels apart, so a state entered on two different + symbols can never be placed in a total co-lex order. This rejects the even + process, Nemo, and butterfly epsilon-machines outright. + """ + graph = labeled_graph(model, symbol_key=symbol_key) + return all(len(labels) <= 1 for labels in graph.in_labels().values()) + + +def _first_axiom_violation(graph: LabeledGraph, rank: dict[Hashable, int]) -> str | None: + """Return a description of the first Wheeler axiom violation, else ``None``.""" + symbol_rank = graph.symbol_rank() + in_labels = graph.in_labels() + + sourceless = [rank[state] for state, labels in in_labels.items() if not labels] + targeted = [rank[state] for state, labels in in_labels.items() if labels] + if sourceless and targeted and max(sourceless) > min(targeted): + return "a state with no incoming edges is ordered after a state with incoming edges" + + # Axiom 1: distinct incoming labels force the order, so every state entered + # on symbol ``a`` must precede every state entered on a later symbol. + by_symbol: dict[Any, list[int]] = defaultdict(list) + for state, labels in in_labels.items(): + if len(labels) > 1: + return f"state {state!r} is entered on more than one symbol: {sorted(labels, key=repr)}" + for symbol in labels: + by_symbol[symbol].append(rank[state]) + ordered_symbols = sorted(by_symbol, key=lambda symbol: symbol_rank[symbol]) + for earlier, later in zip(ordered_symbols, ordered_symbols[1:], strict=False): + if max(by_symbol[earlier]) > min(by_symbol[later]): + return f"states entered on {earlier!r} are not all before states entered on {later!r}" + + # Axiom 2: among equally labeled edges the target order follows the source + # order, so sorting one symbol's edges by source rank must leave the target + # ranks non-decreasing. Edges sharing a source are mutually unconstrained, + # so they are compared as a group against the groups before them. + state_at = {index: state for state, index in rank.items()} + ranked: dict[Any, list[tuple[int, int]]] = defaultdict(list) + for source, symbol, target in graph.edges: + ranked[symbol].append((rank[source], rank[target])) + for symbol, pairs in ranked.items(): + pairs.sort() + ceiling, ceiling_source = -1, -1 + position = 0 + while position < len(pairs): + source_rank = pairs[position][0] + targets = [] + while position < len(pairs) and pairs[position][0] == source_rank: + targets.append(pairs[position][1]) + position += 1 + if min(targets) < ceiling: + return ( + f"edges {state_at[ceiling_source]!r} -{symbol!r}-> {state_at[ceiling]!r} and " + f"{state_at[source_rank]!r} -{symbol!r}-> {state_at[min(targets)]!r} invert the order" + ) + if max(targets) > ceiling: + ceiling, ceiling_source = max(targets), source_rank + return None + + +def check_wheeler_axioms(graph: LabeledGraph, order: Sequence[Hashable]) -> bool: + """Return whether ``order`` satisfies the Wheeler axioms on ``graph``.""" + rank = {state: index for index, state in enumerate(order)} + return _first_axiom_violation(graph, rank) is None + + +def _predecessor_range( + state: Hashable, + *, + predecessors: dict[Hashable, tuple[Hashable, ...]], + block_of: dict[Hashable, int], +) -> tuple[int, int]: + """Lowest and highest block index among ``state``'s predecessors.""" + sources = [block_of[source] for source in predecessors[state]] + return (min(sources), max(sources)) if sources else (-1, -1) + + +def _refine_colex_blocks(graph: LabeledGraph) -> list[list[Hashable]]: + """Refine states into co-lex ordered blocks by predecessor block range. + + Seeds an ordered partition from the incoming label (axiom one) and then + repeatedly splits each block by the lowest and highest block index among + its members' predecessors. In a deterministic automaton two states sharing + an incoming label must have *separated* predecessor sets -- ``u < v`` + forces every predecessor of ``u`` below every predecessor of ``v`` -- so + both ends of that range move monotonically with the Wheeler order. + """ + symbol_rank = graph.symbol_rank() + in_labels = graph.in_labels() + predecessors = graph.predecessors() + + def seed_key(state: Hashable) -> tuple[int, int]: + labels = in_labels[state] + if not labels: + return (0, -1) + return (1, min(symbol_rank[symbol] for symbol in labels)) + + blocks: list[list[Hashable]] = [] + for _key, group in itertools.groupby(sorted(graph.states, key=seed_key), key=seed_key): + blocks.append(list(group)) + + for _round in range(len(graph.states) + 1): + block_of = {state: index for index, block in enumerate(blocks) for state in block} + split_key = partial(_predecessor_range, predecessors=predecessors, block_of=block_of) + refined: list[list[Hashable]] = [] + for block in blocks: + if len(block) == 1: + refined.append(block) + continue + for _key, group in itertools.groupby(sorted(block, key=split_key), key=split_key): + refined.append(list(group)) + if len(refined) == len(blocks): + return refined + blocks = refined + return blocks + + +def wheeler_order( + model: Any, + *, + symbol_key: Callable[[Any], Any] | None = None, + max_tie_permutations: int = MAX_TIE_PERMUTATIONS, +) -> WheelerOrder | None: + """Return a Wheeler order for ``model``, or ``None`` if it is not Wheeler.""" + return wheeler_order_of_graph( + labeled_graph(model, symbol_key=symbol_key), + max_tie_permutations=max_tie_permutations, + ) + + +def wheeler_order_of_graph( + graph: LabeledGraph, + *, + max_tie_permutations: int = MAX_TIE_PERMUTATIONS, +) -> WheelerOrder | None: + """Return a Wheeler order for ``graph``, or ``None`` if it is not Wheeler. + + Tries three increasingly expensive steps. Co-lex partition refinement + usually pins the order outright. Otherwise the maximum co-lex relation is + linearized, which also *disproves* Wheelerness whenever it leaves a pair + incomparable, since it contains every co-lex order. Only states the + relation ranks as mutually comparable are permuted, bounded by + ``max_tie_permutations`` -- deciding Wheelerness is NP-complete for NFAs + :cite:`GibneyThankachan2019`, so some input must stay expensive. + """ + if not graph.states: + return WheelerOrder.from_sequence(()) + if any(len(labels) > 1 for labels in graph.in_labels().values()): + return None + + blocks = _refine_colex_blocks(graph) + candidate = [state for block in blocks for state in block] + if check_wheeler_axioms(graph, candidate): + return WheelerOrder.from_sequence(candidate) + + if _is_deterministic(graph) and all(len(block) == 1 for block in blocks): + # Separated predecessor sets make the refinement order the only + # candidate for a deterministic graph, so a failure here is decisive + # and the maximum co-lex relation has nothing left to contribute. + return None + + chain = _linearize_colex_relation(graph) + if chain is None: + return None + groups = [list(group) for group in chain] + ordered = [state for group in groups for state in group] + if check_wheeler_axioms(graph, ordered): + return WheelerOrder.from_sequence(ordered) + + tied = [index for index, group in enumerate(groups) if len(group) > 1] + total = math.prod(math.factorial(len(groups[index])) for index in tied) + if total > max_tie_permutations: + raise WheelerError( + f"{total} tie-breaking permutations exceed max_tie_permutations=" + f"{max_tie_permutations}; Wheelerness is undecided for this model" + ) + for arrangement in itertools.product(*(itertools.permutations(groups[index]) for index in tied)): + replacement = dict(zip(tied, arrangement, strict=True)) + attempt: list[Hashable] = [] + for index, group in enumerate(groups): + attempt.extend(replacement.get(index, tuple(group))) + if check_wheeler_axioms(graph, attempt): + return WheelerOrder.from_sequence(attempt) + return None + + +def _linearize_colex_relation(graph: LabeledGraph) -> list[tuple[Hashable, ...]] | None: + """Order the states into a chain of mutually comparable groups, or ``None``. + + ``None`` means some pair is incomparable in the maximum co-lex relation, so + no co-lex order compares them and no *total* one exists: the graph is not + Wheeler. This keeps the hard case polynomial instead of leaving unrelated + states, sourceless ones above all, to a factorial search. + """ + relation = maximum_colex_relation_of_graph(graph) + strict = nx.DiGraph() + strict.add_nodes_from(relation) + for left, rights in relation.items(): + for right in rights: + if left != right and left not in relation[right]: + strict.add_edge(left, right) + + # A DAG's reachability order is total exactly when its topological order is + # unique, i.e. Kahn's algorithm never has two ready nodes at once. That is + # linear, where materializing the transitive closure is not. + condensation = nx.condensation(strict) + remaining = dict(condensation.in_degree()) + ready = [node for node, degree in remaining.items() if degree == 0] + chain: list[int] = [] + while len(ready) == 1: + node = ready.pop() + chain.append(node) + for successor in condensation.successors(node): + remaining[successor] -= 1 + if remaining[successor] == 0: + ready.append(successor) + if len(chain) != condensation.number_of_nodes(): + return None + return [tuple(condensation.nodes[node]["members"]) for node in chain] + + +def _is_deterministic(graph: LabeledGraph) -> bool: + """Whether each state has at most one out-edge per symbol.""" + seen: dict[tuple[Hashable, Any], Hashable] = {} + for source, symbol, target in graph.edges: + key = (source, symbol) + if seen.setdefault(key, target) != target: + return False + return True + + +def is_wheeler(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> bool: + """Return whether ``model`` admits a Wheeler order. + + A property of the presentation. Whether the *language* is Wheeler -- whether + some equivalent automaton is -- is a strictly weaker and much more expensive + question, decidable in ``O(mn)`` for a DFA and PSPACE-complete for an NFA + :cite:`Becker2023` :cite:`Alanko2021`. See + :func:`~sofic.shifts.wheeler.wheeler_cover` and + :func:`~sofic.generators.wheeler_epsilon.wheeler_presentation` for the + search over presentations. + """ + return wheeler_order(model, symbol_key=symbol_key) is not None + + +def maximum_colex_relation( + model: Any, + *, + symbol_key: Callable[[Any], Any] | None = None, +) -> dict[Hashable, frozenset[Hashable]]: + """Return the maximum co-lexicographic relation of ``model``. + + Maps each state ``u`` to the states ``v`` with ``u <= v``. Computed as the + greatest fixpoint of the co-lex axioms of :cite:`CotumaccioPrezza2021`: + start from every pair permitted by axiom one, then discard ``u <= v`` + whenever some pair of incoming edges makes axiom two fail. Unlike a minimum + width co-lex *order*, whose computation is NP-hard for NFAs, this relation + always exists and is computable in polynomial time :cite:`Cotumaccio2023`. + """ + return maximum_colex_relation_of_graph(labeled_graph(model, symbol_key=symbol_key)) + + +def maximum_colex_relation_of_graph(graph: LabeledGraph) -> dict[Hashable, frozenset[Hashable]]: + """Greatest fixpoint of the co-lex axioms on ``graph``. + + Every co-lex order of ``graph`` is contained in this relation, since the + fixpoint only ever discards pairs that violate an axiom. A pair left + incomparable is therefore a pair *no* co-lex order can compare. + """ + symbol_rank = graph.symbol_rank() + in_labels = graph.in_labels() + successors: dict[Hashable, set[Hashable]] = {state: set() for state in graph.states} + for source, _symbol, target in graph.edges: + successors[source].add(target) + + def label_precedes(left: Hashable, right: Hashable) -> bool: + """Whether every incoming label of ``left`` is below every one of ``right``.""" + left_labels, right_labels = in_labels[left], in_labels[right] + if not left_labels: + return bool(right_labels) + if not right_labels: + return False + return max(symbol_rank[s] for s in left_labels) < min(symbol_rank[s] for s in right_labels) + + relation: set[tuple[Hashable, Hashable]] = set() + dropped: deque[tuple[Hashable, Hashable]] = deque() + for left in graph.states: + for right in graph.states: + if left == right or not label_precedes(right, left): + relation.add((left, right)) + else: + dropped.append((left, right)) + + # Axiom two fails for a pair exactly when one of its predecessor pairs has + # already failed, so push removals forward along edges rather than + # rescanning every pair on every round. + while dropped: + left_source, right_source = dropped.popleft() + for left in successors[left_source]: + for right in successors[right_source]: + if left == right or (left, right) not in relation: + continue + if in_labels[left] != in_labels[right]: + continue + relation.discard((left, right)) + dropped.append((left, right)) + + return {left: frozenset(right for right in graph.states if (left, right) in relation) for left in graph.states} + + +def colex_width(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> int: + """Return the co-lexicographic width of ``model``. + + The width is the size of the largest antichain of the partial order carried + by the maximum co-lex relation, computed by Dilworth duality as the state + count minus a maximum bipartite matching on the strict order. Width one + means the order is total, i.e. the machine is Wheeler; larger widths bound + the cost of indexing, encoding, and determinizing it + :cite:`CotumaccioPrezza2021`. + + A co-lex order must be antisymmetric, so a pair the relation orients in + both directions is a pair no co-lex order can compare and is counted here + as incomparable. Minimum width over co-lex orders is NP-hard to compute for + NFAs and the relation-based estimate can differ from it + :cite:`Cotumaccio2023`. + + A Wheeler machine always has width one, but the converse needs + :func:`is_input_consistent`: axiom one compares the incoming labels of two + *distinct* states, so it cannot see a single state entered on two different + symbols. The one-state presentation of the full shift is the smallest + example -- width one, yet not Wheeler. Test Wheelerness with + :func:`is_wheeler` rather than ``colex_width(...) == 1``. + """ + relation = maximum_colex_relation(model, symbol_key=symbol_key) + if not relation: + return 0 + + strict = nx.DiGraph() + strict.add_nodes_from(relation) + for left, rights in relation.items(): + for right in rights: + if left != right and left not in relation[right]: + strict.add_edge(left, right) + + # The strict part of a transitive relation is acyclic; condense defensively + # so a non-transitive fixpoint still yields a DAG to close and match on. + condensation = nx.condensation(strict) + if condensation.number_of_nodes() <= 1: + return condensation.number_of_nodes() + + closure = nx.transitive_closure_dag(condensation) + bipartite = nx.Graph() + bipartite.add_nodes_from((("out", node) for node in closure), bipartite=0) + bipartite.add_nodes_from((("in", node) for node in closure), bipartite=1) + bipartite.add_edges_from((("out", source), ("in", target)) for source, target in closure.edges) + matching = nx.bipartite.hopcroft_karp_matching(bipartite, top_nodes=[("out", node) for node in closure]) + matched = sum(1 for node in matching if node[0] == "out") + return closure.number_of_nodes() - matched + + +def _order_of(model: Any, symbol_key: Callable[[Any], Any] | None) -> tuple[LabeledGraph, WheelerOrder]: + graph = labeled_graph(model, symbol_key=symbol_key) + order = wheeler_order(model, symbol_key=symbol_key) + if order is None: + raise WheelerError(f"{type(model).__qualname__} does not admit a Wheeler order") + return graph, order + + +def minimum_wdfa( + dfa: Any, + *, + symbol_key: Callable[[Any], Any] | None = None, +) -> Any: + """Return the minimum Wheeler DFA equivalent to the Wheeler DFA ``dfa``. + + Merges maximal runs of states that are consecutive in the Wheeler order, + share an incoming label, and are Myhill-Nerode equivalent. Merging a range + whose incoming edges carry one label preserves Wheelerness + :cite:`Gagie2017`, and the resulting automaton is the unique smallest + Wheeler DFA for the language :cite:`Alanko2020`. + """ + from sofic.automata.algorithms import minimize + from sofic.automata.dfa import DFA + + graph, order = _order_of(dfa, symbol_key) + in_labels = graph.in_labels() + classes = _nerode_classes(dfa, minimize(dfa)) + + runs: list[list[Hashable]] = [] + for state in order.states: + signature = (in_labels[state], classes[state]) + if runs and (in_labels[runs[-1][-1]], classes[runs[-1][-1]]) == signature: + runs[-1].append(state) + else: + runs.append([state]) + + merged: dict[Hashable, Hashable] = {} + for run in runs: + for state in run: + merged[state] = run[0] + + result = DFA( + input_alphabet=frozenset(dfa.input_alphabet), + initial_states=frozenset(merged[state] for state in dfa.initial_states), + accepting_states=frozenset(merged[state] for state in dfa.accepting_states), + ) + for run in runs: + result.graph.add_state(run[0]) + seen: set[tuple[Hashable, Any, Hashable]] = set() + for source, symbol, target in graph.edges: + edge = (merged[source], symbol, merged[target]) + if edge in seen: + continue + seen.add(edge) + result.graph.add_transition(edge[0], edge[2], **{ATTR_SYMBOL: symbol}) + return result + + +def _nerode_classes(dfa: Any, minimal: Any) -> dict[Hashable, Hashable]: + """Map each state of ``dfa`` to the ``minimal`` state it is equivalent to.""" + from sofic.graph import ATTR_SYMBOL as SYMBOL + + def step(automaton: Any, state: Hashable, symbol: Any) -> Hashable | None: + for transition in automaton.graph.out_transitions(state): + if transition.data.get(SYMBOL) == symbol: + return transition.target + return None + + start = next(iter(dfa.initial_states), None) + minimal_start = next(iter(minimal.initial_states), None) + if start is None or minimal_start is None: + return dict.fromkeys(dfa.states(), 0) + + classes: dict[Hashable, Hashable] = {start: minimal_start} + queue = [(start, minimal_start)] + alphabet = sorted(dfa.input_alphabet, key=repr) + while queue: + state, image = queue.pop() + for symbol in alphabet: + target = step(dfa, state, symbol) + if target is None or target in classes: + continue + classes[target] = step(minimal, image, symbol) + queue.append((target, classes[target])) + + # States unreachable from the start form their own singleton classes. + for state in dfa.states(): + classes.setdefault(state, ("unreachable", repr(state))) + return classes + + +def wnfa_to_wdfa( + nfa: Any, + *, + symbol_key: Callable[[Any], Any] | None = None, +) -> Any: + """Determinize a Wheeler NFA into an equivalent Wheeler DFA. + + Path coherence makes every reachable subset an interval of the Wheeler + order, so the subset construction ranges over intervals rather than subsets + and yields at most ``2n - 1 - |Sigma|`` states :cite:`Alanko2020`. + """ + from sofic.automata.dfa import DFA + + graph, order = _order_of(nfa, symbol_key) + rank = order.rank + by_symbol: dict[Any, list[tuple[int, int]]] = defaultdict(list) + for source, symbol, target in graph.edges: + by_symbol[symbol].append((rank[source], rank[target])) + + def successor(interval: tuple[int, int], symbol: Any) -> tuple[int, int] | None: + targets = [target for source, target in by_symbol[symbol] if interval[0] <= source <= interval[1]] + if not targets: + return None + return (min(targets), max(targets)) + + initial_ranks = [rank[state] for state in nfa.initial_states] + if not initial_ranks: + return DFA(input_alphabet=frozenset(nfa.input_alphabet)) + start = (min(initial_ranks), max(initial_ranks)) + + intervals = {start} + queue = [start] + edges: list[tuple[tuple[int, int], Any, tuple[int, int]]] = [] + while queue: + interval = queue.pop() + for symbol in graph.alphabet: + target = successor(interval, symbol) + if target is None: + continue + edges.append((interval, symbol, target)) + if target not in intervals: + intervals.add(target) + queue.append(target) + + accepting_ranks = {rank[state] for state in nfa.accepting_states} + + def label(interval: tuple[int, int]) -> Hashable: + return tuple(order.states[index] for index in range(interval[0], interval[1] + 1)) + + result = DFA( + input_alphabet=frozenset(nfa.input_alphabet), + initial_states=frozenset({label(start)}), + accepting_states=frozenset( + label(interval) + for interval in intervals + if any(index in accepting_ranks for index in range(interval[0], interval[1] + 1)) + ), + ) + for interval in intervals: + result.graph.add_state(label(interval)) + for source, symbol, target in edges: + result.graph.add_transition(label(source), label(target), **{ATTR_SYMBOL: symbol}) + return result + + +def wheeler_canonical_form(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> tuple[Any, ...]: + """Return a canonical fingerprint of a Wheeler model. + + The Wheeler order is intrinsic, so serializing the machine in that order + gives a representation independent of how its states happen to be named or + inserted. Two Wheeler models share a fingerprint exactly when they are + isomorphic as labeled graphs. + """ + graph, order = _order_of(model, symbol_key) + rank = order.rank + symbol_rank = graph.symbol_rank() + transitions = sorted((rank[source], symbol_rank[symbol], rank[target]) for source, symbol, target in graph.edges) + marked = _marked_states(model, rank) + return (len(order), len(graph.alphabet), tuple(transitions), marked) + + +def wheeler_isomorphic( + left: Any, + right: Any, + *, + symbol_key: Callable[[Any], Any] | None = None, +) -> bool: + """Return whether two Wheeler models are isomorphic as labeled graphs. + + Because the Wheeler order is intrinsic, this compares canonical forms in + time linear in the transition count, where + :func:`~sofic.automata.algorithms.equivalent` has to minimize both inputs. + It is strictly finer than language equivalence: two Wheeler DFAs for the + same language differ here unless both are minimal. Raises + :class:`WheelerError` if either model lacks a Wheeler order. + """ + return wheeler_canonical_form(left, symbol_key=symbol_key) == wheeler_canonical_form(right, symbol_key=symbol_key) + + +def wheeler_state_index(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> Any: + """Return a :class:`~sofic.indexing.StateIndex` in Wheeler order. + + Indexing states this way makes derived matrices canonical: the same machine + yields the same transition matrix no matter how its states were named or + inserted, unlike the insertion order of :meth:`~sofic.base.StateMachine.reindex`. + """ + from sofic.indexing import StateIndex + + _graph, order = _order_of(model, symbol_key) + return StateIndex(order.states) + + +def _marked_states(model: Any, rank: dict[Hashable, int]) -> tuple[tuple[int, ...], ...]: + """Ranks of initial and accepting states, when the model distinguishes them.""" + marks: list[tuple[int, ...]] = [] + for attribute in ("initial_states", "accepting_states"): + states: Iterable[Hashable] = getattr(model, attribute, ()) or () + marks.append(tuple(sorted(rank[state] for state in states if state in rank))) + return tuple(marks) diff --git a/sofic/automata/wheeler_index.py b/sofic/automata/wheeler_index.py new file mode 100644 index 0000000..b49dff0 --- /dev/null +++ b/sofic/automata/wheeler_index.py @@ -0,0 +1,291 @@ +"""Burrows-Wheeler index over a Wheeler-ordered machine. + +Stores a Wheeler graph as the arrays of :cite:`Gagie2017`: out-degrees and +in-degrees in node order, plus the edge labels listed in node order. Because +the Wheeler axioms make the edges sorted by ``(label, source)`` coincide with +the edges sorted by target, following a label maps one node interval onto +another -- the generalization of FM-index backward search to labeled graphs. + +Rank is served by binary search over per-symbol position arrays rather than a +succinct bitvector, which costs a logarithmic factor but keeps the dependency +surface at numpy. +""" + +from __future__ import annotations + +import random +from collections.abc import Callable, Hashable, Iterator, Sequence +from dataclasses import dataclass, field +from typing import Any + +import numpy as np + +from sofic.automata.wheeler import LabeledGraph, WheelerError, WheelerOrder, labeled_graph, wheeler_order_of_graph + +Interval = tuple[int, int] + + +@dataclass(frozen=True, slots=True, eq=False) +class WheelerIndex: + """Searchable Burrows-Wheeler representation of a Wheeler machine.""" + + order: WheelerOrder + alphabet: tuple[Any, ...] + #: Prefix sums of out-degrees in Wheeler order; ``out_start[-1]`` is the edge count. + out_start: np.ndarray + #: Prefix sums of in-degrees in Wheeler order. + in_start: np.ndarray + #: Symbol index of each out-edge, listed in node order then symbol order. + labels: np.ndarray + #: Number of edges whose symbol precedes each symbol index. + symbol_start: np.ndarray + #: Positions in :attr:`labels` carrying each symbol index. + label_positions: tuple[np.ndarray, ...] + initial: Interval | None + accepting: frozenset[int] + #: Reverse adjacency, used for co-lex word ranking. + _backward: dict[tuple[Hashable, Any], frozenset[Hashable]] + _symbol_index: dict[Any, int] + _word_counts: dict[tuple[frozenset[Hashable], int], int] = field(default_factory=dict) + + @classmethod + def from_model( + cls, + model: Any, + *, + symbol_key: Callable[[Any], Any] | None = None, + initial_states: frozenset[Hashable] | None = None, + accepting_states: frozenset[Hashable] | None = None, + ) -> WheelerIndex: + """Build an index for ``model``. + + Initial and accepting states default to the model's own, or to every + state for presentations that mark neither -- the "all states initial" + case of :cite:`Gagie2017` Theorem 6, which is how a shift presents its + factor language. + """ + graph = labeled_graph(model, symbol_key=symbol_key) + order = wheeler_order_of_graph(graph) + if order is None: + raise WheelerError(f"{type(model).__qualname__} does not admit a Wheeler order") + + every = frozenset(graph.states) + if initial_states is None: + initial_states = getattr(model, "initial_states", None) or every + if accepting_states is None: + accepting_states = getattr(model, "accepting_states", None) or every + return cls._build(graph, order, initial_states, accepting_states) + + @classmethod + def _build( + cls, + graph: LabeledGraph, + order: WheelerOrder, + initial_states: frozenset[Hashable], + accepting_states: frozenset[Hashable], + ) -> WheelerIndex: + rank = order.rank + symbol_rank = graph.symbol_rank() + size = len(order.states) + arity = len(graph.alphabet) + + # Edges in node order carry the labels; edges in (label, source) order + # are the same multiset listed by target, which is what makes the + # label-to-target hop a pair of array lookups. + by_source = sorted((rank[source], symbol_rank[symbol], rank[target]) for source, symbol, target in graph.edges) + out_degree = np.zeros(size, dtype=np.int64) + in_degree = np.zeros(size, dtype=np.int64) + for source, _symbol, target in by_source: + out_degree[source] += 1 + in_degree[target] += 1 + + labels = np.fromiter((symbol for _source, symbol, _target in by_source), dtype=np.int64, count=len(by_source)) + symbol_count = np.zeros(arity + 1, dtype=np.int64) + for symbol in labels: + symbol_count[symbol + 1] += 1 + + backward: dict[tuple[Hashable, Any], set[Hashable]] = {} + for source, symbol, target in graph.edges: + backward.setdefault((target, symbol), set()).add(source) + + return cls( + order=order, + alphabet=graph.alphabet, + out_start=np.concatenate(([0], np.cumsum(out_degree))), + in_start=np.concatenate(([0], np.cumsum(in_degree))), + labels=labels, + symbol_start=np.cumsum(symbol_count)[:-1], + label_positions=tuple(np.flatnonzero(labels == symbol) for symbol in range(arity)), + initial=_interval_of(initial_states, rank), + accepting=frozenset(rank[state] for state in accepting_states), + _backward={key: frozenset(value) for key, value in backward.items()}, + _symbol_index=dict(symbol_rank), + ) + + def __len__(self) -> int: + return len(self.order.states) + + def bits(self) -> int: + """Size of the Burrows-Wheeler encoding in bits. + + The ``2(e + n) + e log |A|`` bound of :cite:`Gagie2017` Theorem 6, + excluding the lower-order rank structures. + """ + edges, nodes = int(self.labels.size), len(self) + arity = max(len(self.alphabet), 1) + return 2 * (edges + nodes) + int(np.ceil(edges * np.log2(arity))) + + def step(self, interval: Interval, symbol: Any) -> Interval | None: + """Follow every ``symbol`` edge out of ``interval``, or ``None`` if there are none.""" + index = self._symbol_index.get(symbol) + if index is None: + return None + low = int(self.out_start[interval[0]]) + high = int(self.out_start[interval[1] + 1]) + positions = self.label_positions[index] + first = int(np.searchsorted(positions, low)) + last = int(np.searchsorted(positions, high)) + if first == last: + return None + base = int(self.symbol_start[index]) + return (self._node_of(base + first), self._node_of(base + last - 1)) + + def _node_of(self, in_edge: int) -> int: + """Node owning the given position in the by-target edge ordering.""" + return int(np.searchsorted(self.in_start, in_edge, side="right")) - 1 + + def forward_search(self, word: Sequence[Any], interval: Interval | None = None) -> Interval | None: + """Return the node interval reached by reading ``word``, or ``None``. + + Runs in ``O(|word| log |A|)`` regardless of how many states the word + actually reaches, because path coherence keeps that set an interval + :cite:`Gagie2017`. + """ + current = self.initial if interval is None else interval + if current is None: + return None + for symbol in word: + current = self.step(current, symbol) + if current is None: + return None + return current + + def contains(self, word: Sequence[Any]) -> bool: + """Return whether ``word`` is accepted. + + The indexed replacement for scanning + :meth:`~sofic.shifts.base.SymbolicModel.factor_language`. + """ + interval = self.forward_search(word) + if interval is None: + return False + return any(node in self.accepting for node in range(interval[0], interval[1] + 1)) + + def count_states(self, word: Sequence[Any]) -> int: + """Number of states reachable by reading ``word``.""" + interval = self.forward_search(word) + return 0 if interval is None else interval[1] - interval[0] + 1 + + def states_reached(self, word: Sequence[Any]) -> tuple[Hashable, ...]: + """The states reachable by reading ``word``, in Wheeler order.""" + interval = self.forward_search(word) + if interval is None: + return () + return tuple(self.order.states[node] for node in range(interval[0], interval[1] + 1)) + + # -- Co-lexicographic word ranking ------------------------------------- + # + # Co-lex compares words from the last symbol backwards, which is exactly + # the order the Wheeler order imposes on states via the words reaching + # them. Counting therefore walks the transition graph backwards, subset + # determinizing as it goes so that each branch is a distinct word. + + def _predecessors(self, states: frozenset[Hashable], symbol: Any) -> frozenset[Hashable]: + sources: set[Hashable] = set() + for state in states: + sources.update(self._backward.get((state, symbol), ())) + return frozenset(sources) + + def _count_from(self, states: frozenset[Hashable], length: int) -> int: + """Distinct words of ``length`` symbols that end at some state in ``states``.""" + if length == 0: + return 1 if states else 0 + cached = self._word_counts.get((states, length)) + if cached is not None: + return cached + total = sum( + self._count_from(sources, length - 1) + for symbol in self.alphabet + if (sources := self._predecessors(states, symbol)) + ) + self._word_counts[(states, length)] = total + return total + + def _accepting_states(self) -> frozenset[Hashable]: + return frozenset(self.order.states[node] for node in sorted(self.accepting)) + + def count_words(self, length: int) -> int: + """Number of distinct accepted words of ``length`` symbols.""" + return self._count_from(self._accepting_states(), length) + + def words_of_length(self, length: int) -> Iterator[tuple[Any, ...]]: + """Yield the accepted words of ``length`` symbols in co-lexicographic order.""" + for index in range(self.count_words(length)): + yield self.unrank_word(index, length) + + def rank_word(self, word: Sequence[Any]) -> int: + """Position of ``word`` among equally long accepted words, co-lex ordered.""" + states = self._accepting_states() + position = 0 + for depth, symbol in enumerate(reversed(tuple(word))): + remaining = len(word) - depth - 1 + for candidate in self.alphabet: + sources = self._predecessors(states, candidate) + if candidate == symbol: + if not sources: + raise ValueError(f"{tuple(word)!r} is not an accepted word") + states = sources + break + position += self._count_from(sources, remaining) + else: + raise ValueError(f"symbol {symbol!r} is not in the alphabet") + return position + + def unrank_word(self, index: int, length: int) -> tuple[Any, ...]: + """Inverse of :meth:`rank_word`: the ``index``-th co-lex accepted word.""" + total = self.count_words(length) + if not 0 <= index < total: + raise IndexError(f"rank {index} out of range for {total} words of length {length}") + states = self._accepting_states() + suffix: list[Any] = [] + remaining = index + for depth in range(length): + for candidate in self.alphabet: + sources = self._predecessors(states, candidate) + block = self._count_from(sources, length - depth - 1) + if remaining < block: + suffix.append(candidate) + states = sources + break + remaining -= block + return tuple(reversed(suffix)) + + def sample_word(self, length: int, rng: random.Random | None = None) -> tuple[Any, ...]: + """Sample uniformly from the accepted words of ``length`` symbols.""" + total = self.count_words(length) + if total == 0: + raise ValueError(f"no accepted words of length {length}") + chooser = random.Random() if rng is None else rng + return self.unrank_word(chooser.randrange(total), length) + + +def _interval_of(states: frozenset[Hashable], rank: dict[Hashable, int]) -> Interval | None: + ranks = [rank[state] for state in states if state in rank] + if not ranks: + return None + return (min(ranks), max(ranks)) + + +def wheeler_index(model: Any, **kwargs: Any) -> WheelerIndex: + """Build a :class:`WheelerIndex` for ``model``.""" + return WheelerIndex.from_model(model, **kwargs) diff --git a/sofic/base.py b/sofic/base.py index bfa14da..3509c3e 100644 --- a/sofic/base.py +++ b/sofic/base.py @@ -83,6 +83,42 @@ def reverse(self) -> Self: result.graph = self.graph.reverse() return result + def is_wheeler(self) -> bool: + """Return whether this presentation admits a Wheeler order. + + Wheelerness is a property of a *presentation*, not of the language it + generates: a process whose minimal presentation is not Wheeler may + still have a larger one that is. Nor does it coincide with finite + Markov order -- every definite presentation has a Wheeler order-``R`` + de Bruijn form, yet Wheeler presentations of infinite Markov order also + exist :cite:`Gagie2017` :cite:`Alanko2020`. + """ + from sofic.automata.wheeler import is_wheeler + + return is_wheeler(self) + + def wheeler_order(self) -> Any: + """Return a :class:`~sofic.automata.wheeler.WheelerOrder`, or ``None``. + + Sorts states by the co-lexicographic rank of the words reaching them, + so each state owns an interval of the sorted prefixes + :cite:`Gagie2017`. + """ + from sofic.automata.wheeler import wheeler_order + + return wheeler_order(self) + + def colex_width(self) -> int: + """Co-lexicographic width of this presentation; Wheeler is width one. + + Bounds the cost of indexing, encoding, and determinizing the machine + :cite:`CotumaccioPrezza2021`. Width one implies Wheelerness only for + input-consistent presentations -- prefer :meth:`is_wheeler`. + """ + from sofic.automata.wheeler import colex_width + + return colex_width(self) + def __repr__(self) -> str: states = list(self.states()) transitions = list(self.transitions()) diff --git a/sofic/examples/__init__.py b/sofic/examples/__init__.py index b7b2156..e274ddd 100644 --- a/sofic/examples/__init__.py +++ b/sofic/examples/__init__.py @@ -31,6 +31,7 @@ tent_map_misiurewicz_partition_information_expected, tent_map_misiurewicz_partition_symbol_matrices, tent_map_misiurewicz_reverse, + wheeler_infinite_order_process, ) from sofic.examples.processes import * from sofic.examples.processes import __all__ as _process_all @@ -79,6 +80,7 @@ "tent_map_misiurewicz_partition_information_expected", "tent_map_misiurewicz_partition_symbol_matrices", "tent_map_misiurewicz_reverse", + "wheeler_infinite_order_process", ] __all__ += _process_all __all__ = [name for name in __all__ if name not in {"processes", "shifts"}] diff --git a/sofic/examples/epsilon_machines.py b/sofic/examples/epsilon_machines.py index c999185..d41c905 100644 --- a/sofic/examples/epsilon_machines.py +++ b/sofic/examples/epsilon_machines.py @@ -1062,6 +1062,50 @@ def tent_map_misiurewicz_information_expected(a: Any | None = None) -> dict[str, } +def wheeler_infinite_order_process(p: float = 0.5, q: float = 0.5, r: float = 0.5) -> EpsilonMachine: + """Five-state Wheeler ε-machine of infinite Markov order. + + The discriminating example separating the Wheeler property from finite + memory. Its causal states admit the Wheeler order ``A < E < B < C < D``, so + every state owns an interval of the recency-ordered pasts, yet no bounded + window of symbols fixes the state: :meth:`~EpsilonMachine.markov_order` is + infinite. Wheelerness is therefore *not* a restatement of definiteness. + + Found by exhaustive search over binary topological ε-machines with + :func:`~sofic.generators.topological_epsilon_enumeration.iter_topological_epsilon_machines`; + twenty of the 35186 five-state machines share both properties. No prior + source states this example. + """ + for name, value in (("p", p), ("q", q), ("r", r)): + if not 0.0 < value < 1.0: + raise ValueError(f"{name} must be in (0, 1)") + states = ("A", "B", "C", "D", "E") + return from_symbol_matrices( + states, + (0, 1), + { + 0: np.array( + [ + [p, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0], + [q, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0], + [r, 0.0, 0.0, 0.0, 0.0], + ] + ), + 1: np.array( + [ + [0.0, 1.0 - p, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0 - q, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0 - r, 0.0, 0.0], + ] + ), + }, + ) + + def ellison_fig9_reverse() -> EpsilonMachine: """Reverse ε-machine from Ellison et al., arXiv:1107.2168, Fig.~9. diff --git a/sofic/generators/__init__.py b/sofic/generators/__init__.py index 07d2d97..35a7388 100644 --- a/sofic/generators/__init__.py +++ b/sofic/generators/__init__.py @@ -55,6 +55,15 @@ iter_topological_epsilon_machines, iter_topological_epsilon_strings, ) +from sofic.generators.wheeler_epsilon import ( + colex_cdf, + cylinder_measure, + debruijn_presentation, + wheeler_complexity_gap, + wheeler_presentation, + wheeler_statistical_complexity, + word_cylinder_measure, +) __all__ = [ "BidirectionalEpsilonMachine", @@ -99,7 +108,10 @@ "QuasiStochasticModel", "StochasticModel", "WynerGenerativeModel", + "colex_cdf", "count_topological_epsilon_machines", + "cylinder_measure", + "debruijn_presentation", "epsilon_machine_to_idfa_string", "functional_generative_model", "gacs_korner_generative_model", @@ -111,5 +123,9 @@ "iter_topological_epsilon_machines", "iter_topological_epsilon_strings", "minimal_generative_model", + "wheeler_complexity_gap", + "wheeler_presentation", + "wheeler_statistical_complexity", + "word_cylinder_measure", "wyner_generative_model", ] diff --git a/sofic/generators/epsilon_machine.py b/sofic/generators/epsilon_machine.py index dcd0709..3832e6f 100644 --- a/sofic/generators/epsilon_machine.py +++ b/sofic/generators/epsilon_machine.py @@ -441,6 +441,32 @@ def is_markov(self) -> bool: order = self.markov_order() return not isinstance(order, float) or math.isfinite(order) + def wheeler_presentation(self, **kwargs: Any) -> Any: + """Return a Wheeler presentation of this process. + + Returns ``self`` when :meth:`~sofic.base.StateMachine.is_wheeler`, + otherwise the smallest edge-machine refinement whose states sort + co-lexicographically. See + :func:`sofic.generators.wheeler_epsilon.wheeler_presentation` -- + original, uncited work. + """ + from sofic.generators.wheeler_epsilon import wheeler_presentation + + return wheeler_presentation(self, **kwargs) + + def wheeler_statistical_complexity(self, **kwargs: Any) -> float: + """``C_W``: state entropy of a Wheeler presentation, in bits. + + Always at least :meth:`statistical_complexity`, with equality exactly + when the epsilon-machine is already Wheeler. The gap is the memory + spent making the causal states intervals of the recency order on pasts. + See :func:`sofic.generators.wheeler_epsilon.wheeler_statistical_complexity` + -- original, uncited work. + """ + from sofic.generators.wheeler_epsilon import wheeler_statistical_complexity + + return wheeler_statistical_complexity(self, **kwargs) + def cryptic_order(self) -> int | float: """Cryptic order ``k_chi``: retrodiction depth after synchronization. diff --git a/sofic/generators/synchronization.py b/sofic/generators/synchronization.py index d3e59df..f0b669b 100644 --- a/sofic/generators/synchronization.py +++ b/sofic/generators/synchronization.py @@ -7,6 +7,7 @@ from __future__ import annotations +import bisect import math from collections.abc import Hashable, Mapping from dataclasses import dataclass, field @@ -14,6 +15,7 @@ import networkx as nx +from sofic.automata.wheeler import LabeledGraph, WheelerOrder, wheeler_order_of_graph from sofic.exceptions import UnifilarityError from sofic.graph import ATTR_EMISSION, ATTR_SYMBOL @@ -61,8 +63,31 @@ def build_topological_graph_from_transitions( return TopologicalUnifilarGraph(states=states, alphabet=alphabet, transitions=dict(transitions)) +def labeled_graph_of(graph: TopologicalUnifilarGraph) -> LabeledGraph: + """View a topological graph as a :class:`~sofic.automata.wheeler.LabeledGraph`.""" + return LabeledGraph( + states=tuple(sorted(graph.states, key=repr)), + alphabet=tuple(sorted(graph.alphabet, key=repr)), + edges=tuple((source, symbol, target) for (source, symbol), target in graph.transitions.items()), + ) + + def power_automaton(graph: TopologicalUnifilarGraph) -> PowerAutomaton: - """Build the power automaton via subset construction from the full state set.""" + """Build the power automaton by subset construction from the full state set. + + When ``graph`` is Wheeler the construction runs over co-lex *intervals* + instead of arbitrary subsets. Path coherence guarantees the two agree + :cite:`Gagie2017`, but the interval form visits at most ``n(n+1)/2`` states + and steps in ``O(log n)`` rather than ``O(n)``, so the synchronization + orders built on top of it stay polynomial. + """ + order = wheeler_order_of_graph(labeled_graph_of(graph)) + if order is not None: + return _interval_power_automaton(graph, order) + return _subset_power_automaton(graph) + + +def _subset_power_automaton(graph: TopologicalUnifilarGraph) -> PowerAutomaton: start = frozenset(graph.states) pa = PowerAutomaton(graph=graph, start=start) queue = [start] @@ -81,6 +106,55 @@ def power_automaton(graph: TopologicalUnifilarGraph) -> PowerAutomaton: return pa +def _interval_power_automaton(graph: TopologicalUnifilarGraph, order: WheelerOrder) -> PowerAutomaton: + """Subset construction restricted to intervals of a Wheeler order. + + Axiom two makes the target rank non-decreasing in the source rank for each + symbol, so the image of a rank interval is bracketed by the first and last + edges whose sources fall inside it -- two binary searches per step. + """ + edges_by_symbol: dict[Any, tuple[list[int], list[int]]] = {} + for symbol in graph.alphabet: + pairs = sorted( + (order.rank[source], order.rank[target]) + for (source, edge_symbol), target in graph.transitions.items() + if edge_symbol == symbol + ) + if pairs: + edges_by_symbol[symbol] = ([source for source, _ in pairs], [target for _, target in pairs]) + + def image(interval: tuple[int, int], symbol: Any) -> tuple[int, int] | None: + found = edges_by_symbol.get(symbol) + if found is None: + return None + sources, targets = found + low = bisect.bisect_left(sources, interval[0]) + high = bisect.bisect_right(sources, interval[1]) - 1 + if low > high: + return None + return (targets[low], targets[high]) + + def as_subset(interval: tuple[int, int]) -> frozenset[Hashable]: + return frozenset(order.states[rank] for rank in range(interval[0], interval[1] + 1)) + + start = (0, len(order.states) - 1) + pa = PowerAutomaton(graph=graph, start=as_subset(start)) + queue = [start] + seen = {start} + while queue: + current = queue.pop(0) + out_map = pa.transitions.setdefault(as_subset(current), {}) + for symbol in graph.alphabet: + successor = image(current, symbol) + if successor is None: + continue + out_map[symbol] = as_subset(successor) + if successor not in seen: + seen.add(successor) + queue.append(successor) + return pa + + def _power_automaton_digraph(pa: PowerAutomaton) -> nx.MultiDiGraph: """Materialize the power automaton as a networkx graph. diff --git a/sofic/generators/wheeler_epsilon.py b/sofic/generators/wheeler_epsilon.py new file mode 100644 index 0000000..707222b --- /dev/null +++ b/sofic/generators/wheeler_epsilon.py @@ -0,0 +1,247 @@ +"""Co-lexicographic structure of epsilon-machines. + +.. warning:: + + Everything in this module is original and uncited. Wheeler automata are a + purely topological theory :cite:`Gagie2017` :cite:`Alanko2020`; a literature + search turned up no treatment of weighted, probabilistic, or + information-theoretic Wheeler automata, so no canonical source exists for + the quantities defined here. They are offered as proposals, and each + docstring says so. + +The bridge from Wheeler theory to computational mechanics is that +co-lexicographic order compares words from the last symbol backwards, which is +the *recency* order on pasts. A presentation is Wheeler exactly when its state +partition of history space is an interval partition under recency, so each +state owns a contiguous band of pasts rather than an arbitrary set. That makes +the stationary distribution a genuine cumulative distribution over pasts, and +it prices the constraint: a process whose causal states are not recency +intervals must split them to get one, and the extra states cost entropy. +""" + +from __future__ import annotations + +from collections.abc import Sequence +from typing import Any + +import numpy as np + +from sofic.automata.wheeler import WheelerError, WheelerOrder, wheeler_order +from sofic.generators.mealy import MealyHMM + +#: How many edge-machine refinements :func:`wheeler_presentation` will try. +DEFAULT_MAX_REFINEMENT = 4 + +#: Refinements past this many states are abandoned rather than searched. +DEFAULT_MAX_STATES = 2048 + + +def wheeler_presentation( + machine: Any, + *, + max_refinement: int = DEFAULT_MAX_REFINEMENT, + max_states: int = DEFAULT_MAX_STATES, +) -> MealyHMM: + """Return a Wheeler presentation of the process ``machine`` generates. + + Returns ``machine`` itself when it is already Wheeler. Otherwise, when the + Markov order ``R`` is finite, returns the order-``R`` de Bruijn + presentation of :func:`debruijn_presentation`, whose states *are* the + length-``R`` words and so sort co-lexicographically by construction. + Failing that it tries :func:`~sofic.generators.edge_machine.hmm_to_edge_machine` + refinements, whose order-``k`` states are length-``k`` transition paths and + are therefore entered on a single symbol. + + Raises :class:`~sofic.automata.wheeler.WheelerError` if no refinement in + range works. Wheeler languages are star-free :cite:`Alanko2021`, so a + process that counts modulo anything -- the even process is the standard + example -- has no Wheeler presentation at any order. + + .. note:: Original, uncited: no literature treats probabilistic Wheeler + presentations. + """ + from sofic.generators.edge_machine import hmm_to_edge_machine + + if wheeler_order(machine) is not None: + return machine + + order = machine.markov_order() if hasattr(machine, "markov_order") else float("inf") + budget = max_refinement if order == float("inf") else min(int(order), max_refinement) + if order != float("inf"): + debruijn = debruijn_presentation(machine, max(int(order), 1)) + if debruijn is not None and wheeler_order(debruijn) is not None: + return debruijn + for iterations in range(1, budget + 1): + refined = hmm_to_edge_machine(machine, iterations=iterations) + size = len(list(refined.states())) + if size > max_states: + raise WheelerError( + f"refinement {iterations} needs {size} states, over max_states={max_states}; " + "raise the cap if the process really is Wheeler at this depth" + ) + 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" + ) + + +def debruijn_presentation(machine: Any, order: int) -> MealyHMM | None: + """Return the order-``order`` de Bruijn presentation, or ``None``. + + States are the length-``order`` words of the process, which is well defined + only once ``order`` reaches the Markov order ``R``, so that each word pins + down one causal state; below that this returns ``None``. Because the states + *are* words, sorting them co-lexicographically satisfies the Wheeler axioms + outright, which is why every finite-Markov-order process has a Wheeler + presentation even when its epsilon-machine has none. + + Distinct from :func:`~sofic.generators.edge_machine.hmm_to_edge_machine`, + whose states are transition paths: a path remembers where it started, so + the edge machine can stay unsortable where the de Bruijn form is not. + + .. note:: Original, uncited as a Wheeler construction, though the de Bruijn + presentation itself is standard. + """ + from sofic.generators.synchronization import graph_from_epsilon_machine + from sofic.graph import ATTR_EMISSION, ATTR_PROB + + if order < 1: + raise ValueError("order must be at least one") + graph = graph_from_epsilon_machine(machine) + alphabet = sorted(graph.alphabet, key=repr) + + emissions: dict[tuple[Any, Any], tuple[Any, Any]] = {} + for transition in machine.transitions(): + symbol = transition.data.get(ATTR_EMISSION) + if symbol is not None: + emissions[(transition.source, symbol)] = (transition.data[ATTR_PROB], transition.target) + + reached: dict[tuple[Any, ...], frozenset[Any]] = {(): frozenset(graph.states)} + for _step in range(order): + extended: dict[tuple[Any, ...], frozenset[Any]] = {} + for word, states in reached.items(): + for symbol in alphabet: + image = graph.delta_set(states, symbol) + if image: + extended[(*word, symbol)] = image + reached = extended + if not reached or any(len(states) != 1 for states in reached.values()): + return None + + result = MealyHMM(observation_alphabet=frozenset(alphabet)) + for word in reached: + result.graph.add_state(word) + for word, states in reached.items(): + state = next(iter(states)) + for symbol in alphabet: + found = emissions.get((state, symbol)) + if found is None: + continue + probability, _target = found + successor = (*word, symbol)[-order:] + if successor in reached: + result.add_transition(word, successor, symbol, probability) + + index = result.reindex() + stationary = result.stationary_distribution() + result.initial_distribution = {index.state(i): float(stationary[i]) for i in range(len(index))} + return result + + +def wheeler_statistical_complexity(machine: Any, **kwargs: Any) -> float: + """``C_W``: the entropy of the states of a Wheeler presentation, in bits. + + Every Wheeler presentation refines the causal-state partition, so + ``C_W >= C_mu`` always, with equality exactly when the epsilon-machine is + itself Wheeler. The gap is the entropic price of co-lex sortability: the + extra memory a process must carry for its states to be intervals of the + recency order on pasts, over and above the memory needed to predict it. + + A second bound comes free. A Wheeler presentation is input consistent, so + every state determines the symbol that entered it, making the last symbol a + function of the state: ``C_W >= H[X_0]`` for a stationary process. A fair + coin therefore has ``C_mu = 0`` but ``C_W = 1``, since sortability forces it + to remember a bit it does not need. Together, ``C_W >= max(C_mu, H[X_0])``. + + .. note:: Original, uncited: proposed here, with no canonical source. + """ + from sofic.generators.measures import state_entropy + + return float(state_entropy(wheeler_presentation(machine, **kwargs))) + + +def wheeler_complexity_gap(machine: Any, **kwargs: Any) -> float: + """``C_W - C_mu``: bits of memory spent purely on co-lex sortability. + + Zero exactly when the epsilon-machine is already Wheeler. + + .. note:: Original, uncited: proposed here, with no canonical source. + """ + return wheeler_statistical_complexity(machine, **kwargs) - float(machine.statistical_complexity()) + + +def colex_cdf(machine: Any) -> tuple[tuple[Any, ...], np.ndarray]: + """Cumulative stationary distribution over states in Wheeler order. + + Returns the ordered states and the running total of their stationary + probabilities, so entry ``i`` is the chance the current past falls at or + below state ``i`` in the recency order. The Wheeler order is what makes + this cumulative sum mean anything: under any other state numbering the + partial sums are arbitrary, whereas here they sweep history space from the + most remote pasts to the most recent, which is what arithmetic coding over + pasts needs. + + Raises :class:`~sofic.automata.wheeler.WheelerError` if ``machine`` is not + Wheeler; refine it with :func:`wheeler_presentation` first. + + .. note:: Original, uncited: proposed here, with no canonical source. + """ + order = _require_order(machine) + index = machine.reindex() + stationary = np.asarray(machine.stationary_distribution(), dtype=float) + masses = np.array([stationary[index.index(state)] for state in order.states]) + return order.states, np.cumsum(masses) + + +def cylinder_measure(machine: Any, low: int, high: int) -> float: + """Stationary probability that the current state lies in Wheeler ranks ``[low, high]``. + + Because the Wheeler order sorts states by the co-lex rank of the words + reaching them, a rank interval is a *cylinder of pasts* under the recency + order, and this is its measure. Pair it with + :meth:`~sofic.automata.wheeler_index.WheelerIndex.forward_search`, which + returns exactly such an interval for a word. + + .. note:: Original, uncited: proposed here, with no canonical source. + """ + _states, cumulative = colex_cdf(machine) + if not 0 <= low <= high < len(cumulative): + raise IndexError(f"interval ({low}, {high}) outside 0..{len(cumulative) - 1}") + below = cumulative[low - 1] if low > 0 else 0.0 + return float(cumulative[high] - below) + + +def word_cylinder_measure(machine: Any, word: Sequence[Any]) -> float: + """Stationary probability that the current state is one reachable by ``word``. + + Not the probability of seeing ``word``: it is the mass of the co-lex + interval that ``word`` selects, i.e. how much of history space is + consistent with having just read it. + + .. note:: Original, uncited: proposed here, with no canonical source. + """ + from sofic.automata.wheeler_index import WheelerIndex + + interval = WheelerIndex.from_model(machine).forward_search(word) + if interval is None: + return 0.0 + return cylinder_measure(machine, *interval) + + +def _require_order(machine: Any) -> WheelerOrder: + order = wheeler_order(machine) + if order is None: + raise WheelerError(f"{type(machine).__qualname__} is not Wheeler; call wheeler_presentation() first") + return order diff --git a/sofic/serialization.py b/sofic/serialization.py index 3ac8d71..dac0024 100644 --- a/sofic/serialization.py +++ b/sofic/serialization.py @@ -331,7 +331,13 @@ def _specs() -> tuple[_ModelSpec, ...]: from sofic.generators.quasi_realization import QuasiRealization from sofic.generators.stack_hmm import HiddenMarkovStackModel from sofic.shifts.base import SymbolicModel - from sofic.shifts.covers import LeftFischerCover, LeftKriegerCover, RightFischerCover, RightKriegerCover + from sofic.shifts.covers import ( + LeftFischerCover, + LeftKriegerCover, + RightFischerCover, + RightKriegerCover, + WheelerCover, + ) from sofic.shifts.markov_dyck import MarkovDyckShift from sofic.shifts.sft import ShiftOfFiniteType from sofic.shifts.sofic import SoficShift @@ -447,4 +453,5 @@ def _specs() -> tuple[_ModelSpec, ...]: _spec(RightFischerCover, symbolic), _spec(LeftKriegerCover, symbolic), _spec(RightKriegerCover, symbolic), + _spec(WheelerCover, symbolic), ) diff --git a/sofic/shifts/__init__.py b/sofic/shifts/__init__.py index dfdc2db..157ad34 100644 --- a/sofic/shifts/__init__.py +++ b/sofic/shifts/__init__.py @@ -6,6 +6,7 @@ LeftKriegerCover, RightFischerCover, RightKriegerCover, + WheelerCover, ) from sofic.shifts.dyck_enumeration import ( DyckGraphString, @@ -23,6 +24,14 @@ from sofic.shifts.sofic_relation import SoficRelation from sofic.shifts.textile import TextileSystem from sofic.shifts.tmc import TopologicalMarkovChain +from sofic.shifts.wheeler import ( + higher_block_presentation, + is_wheeler_shift, + right_resolve, + wheeler_cover, + wheeler_index_of_shift, + wheeler_order_of_shift, +) __all__ = [ "DyckGraphString", @@ -44,5 +53,12 @@ "SymbolicModel", "TextileSystem", "TopologicalMarkovChain", + "WheelerCover", "full_shift", + "higher_block_presentation", + "is_wheeler_shift", + "right_resolve", + "wheeler_cover", + "wheeler_index_of_shift", + "wheeler_order_of_shift", ] diff --git a/sofic/shifts/covers.py b/sofic/shifts/covers.py index 617467a..839e883 100644 --- a/sofic/shifts/covers.py +++ b/sofic/shifts/covers.py @@ -45,3 +45,19 @@ def from_sofic(cls, shift: SoficShift, **kwargs: Any) -> RightKriegerCover: from sofic.shifts.cover_construction import right_krieger_from_sofic return right_krieger_from_sofic(shift) + + +class WheelerCover(SoficShift): + """Wheeler presentation of a sofic shift, merged down. + + Unlike the Fischer and Krieger covers this one need not exist: Wheeler + languages are star-free :cite:`Alanko2021`, and even for a shift that has + one the presentation may need extra symbols of memory. See + :func:`sofic.shifts.wheeler.wheeler_cover`. + """ + + @classmethod + def from_sofic(cls, shift: SoficShift, **kwargs: Any) -> WheelerCover: + from sofic.shifts.wheeler import wheeler_cover + + return wheeler_cover(shift, **kwargs) diff --git a/sofic/shifts/wheeler.py b/sofic/shifts/wheeler.py new file mode 100644 index 0000000..2f351bf --- /dev/null +++ b/sofic/shifts/wheeler.py @@ -0,0 +1,242 @@ +"""Wheeler presentations of sofic shifts. + +A shift presents its factor language with every state both initial and +accepting, which is the "all states initial" case of :cite:`Gagie2017` +Theorem 6. The minimal right-resolving presentation of a shift is rarely +Wheeler on its own -- a state reachable on two different symbols already +violates the axioms -- but remembering the last few symbols restores input +consistency, and for many shifts that is enough. + +Not every sofic shift has a Wheeler presentation at any order: Wheeler +languages are star-free :cite:`ShyrThierrin1974` :cite:`Alanko2021`, so a shift +whose syntactic monoid contains a nontrivial group, such as the even shift, +is excluded outright. +""" + +from __future__ import annotations + +from collections import defaultdict, deque +from collections.abc import Hashable +from typing import Any + +from sofic.automata.wheeler import WheelerError, WheelerOrder, labeled_graph, wheeler_order +from sofic.graph import ATTR_SYMBOL, TransitionGraph +from sofic.shifts.covers import WheelerCover +from sofic.shifts.sofic import SoficShift + +#: How many symbols of memory :func:`wheeler_cover` will add before giving up. +DEFAULT_MAX_ORDER = 6 + +#: Word length used to compare follower languages when merging states. +DEFAULT_FOLLOWER_DEPTH = 8 + +BlockState = tuple[tuple[Any, ...], Hashable] + + +def higher_block_presentation(shift: SoficShift, order: int) -> SoficShift: + """Return the order-``order`` refinement of ``shift``. + + States become ``(last order symbols, original state)``. Every edge into + such a state carries the last symbol of its word, so any refinement of + order at least one is input-consistent -- the cheapest necessary condition + for the Wheeler axioms. + """ + if order < 0: + raise ValueError("order must be nonnegative") + if order == 0: + return shift + + outgoing: dict[Hashable, list[tuple[Any, Hashable]]] = defaultdict(list) + for transition in shift.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is not None: + outgoing[transition.source].append((symbol, transition.target)) + + current: set[BlockState] = {((), state) for state in shift.states()} + for _step in range(order): + current = { + ((word + (symbol,))[-order:], target) for word, state in current for symbol, target in outgoing[state] + } + + # Close under transitions so the refinement is a genuine presentation. + reachable = set(current) + queue = deque(current) + edges: set[tuple[BlockState, Any, BlockState]] = set() + while queue: + word, state = queue.popleft() + for symbol, target in outgoing[state]: + successor = ((word + (symbol,))[-order:], target) + edges.add(((word, state), symbol, successor)) + if successor not in reachable: + reachable.add(successor) + queue.append(successor) + + graph = TransitionGraph() + for block in reachable: + graph.add_state(block) + for source, symbol, target in edges: + graph.add_transition(source, target, **{ATTR_SYMBOL: symbol}) + return SoficShift(graph=graph, symbol_alphabet=shift.symbol_alphabet).trim_transient() + + +def right_resolve(shift: SoficShift) -> SoficShift: + """Return a right-resolving presentation of the same factor language. + + Subset construction seeded with every state, since a shift presents its + factor language with all states initial. Returns ``shift`` unchanged when + it is already right-resolving. + """ + if shift.is_unifilar(): + return shift + + outgoing: dict[Hashable, list[tuple[Any, Hashable]]] = defaultdict(list) + for transition in shift.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is not None: + outgoing[transition.source].append((symbol, transition.target)) + + start = frozenset(shift.states()) + subsets = {start} + queue = deque([start]) + edges: set[tuple[frozenset[Hashable], Any, frozenset[Hashable]]] = set() + while queue: + current = queue.popleft() + successors: dict[Any, set[Hashable]] = defaultdict(set) + for state in current: + for symbol, target in outgoing[state]: + successors[symbol].add(target) + for symbol, targets in successors.items(): + successor = frozenset(targets) + edges.add((current, symbol, successor)) + if successor not in subsets: + subsets.add(successor) + queue.append(successor) + + graph = TransitionGraph() + for subset in subsets: + graph.add_state(subset) + for source, symbol, target in edges: + graph.add_transition(source, target, **{ATTR_SYMBOL: symbol}) + return SoficShift(graph=graph, symbol_alphabet=shift.symbol_alphabet).trim_transient() + + +def wheeler_cover( + shift: SoficShift, + *, + max_order: int = DEFAULT_MAX_ORDER, + follower_depth: int = DEFAULT_FOLLOWER_DEPTH, +) -> WheelerCover: + """Return the smallest Wheeler presentation found for ``shift``. + + Right-resolves the presentation, then tries it as given followed by its + order-1, order-2, ... refinements up to ``max_order``, merging the first + Wheeler one down. Raises :class:`~sofic.automata.wheeler.WheelerError` when + no refinement in range is Wheeler; that is evidence, not proof, that the + shift is non-Wheeler. + """ + trimmed = right_resolve(shift.trim_transient()) + for order in range(max_order + 1): + candidate = higher_block_presentation(trimmed, order) + 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" + ) + + +def is_wheeler_shift(shift: SoficShift, *, max_order: int = DEFAULT_MAX_ORDER) -> bool: + """Return whether some refinement of ``shift`` up to ``max_order`` is Wheeler.""" + try: + wheeler_cover(shift, max_order=max_order) + except WheelerError: + return False + return True + + +def wheeler_order_of_shift(shift: SoficShift) -> int | None: + """Smallest refinement order at which ``shift`` becomes Wheeler, or ``None``. + + Zero means the right-resolving presentation is already Wheeler. This is a + property of the presentation as much as of the shift, since refining a + non-minimal presentation can need more memory than refining a minimal one. + """ + trimmed = right_resolve(shift.trim_transient()) + for order in range(DEFAULT_MAX_ORDER + 1): + if wheeler_order(higher_block_presentation(trimmed, order)) is not None: + return order + return None + + +def _follower_signature(shift: SoficShift, depth: int) -> dict[Hashable, frozenset[tuple[Any, ...]]]: + """Bounded follower language of every state, the shift analogue of Nerode classes.""" + outgoing: dict[Hashable, list[tuple[Any, Hashable]]] = defaultdict(list) + for transition in shift.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is not None: + outgoing[transition.source].append((symbol, transition.target)) + + signatures: dict[Hashable, frozenset[tuple[Any, ...]]] = {} + for state in shift.states(): + words: set[tuple[Any, ...]] = set() + queue: deque[tuple[Hashable, tuple[Any, ...]]] = deque([(state, ())]) + while queue: + current, prefix = queue.popleft() + if len(prefix) >= depth: + continue + for symbol, target in outgoing[current]: + word = prefix + (symbol,) + words.add(word) + queue.append((target, word)) + signatures[state] = frozenset(words) + return signatures + + +def _merge_wheeler_runs(shift: SoficShift, order: WheelerOrder, *, follower_depth: int) -> WheelerCover: + """Collapse Wheeler-consecutive states that share an incoming label and followers. + + Merging a range of states whose incoming edges all carry one label keeps the + graph Wheeler :cite:`Gagie2017`, and equal follower languages make the merge + language-preserving; together these are the merges of the minimum-WDFA + construction :cite:`Alanko2020`. Only *adjacent* runs collapse, so a + presentation carrying states the Wheeler order separates -- a source state + with no incoming edges, for one -- keeps them. + """ + in_labels = labeled_graph(shift).in_labels() + followers = _follower_signature(shift, follower_depth) + + runs: list[list[Hashable]] = [] + for state in order.states: + signature = (in_labels[state], followers[state]) + if runs and (in_labels[runs[-1][-1]], followers[runs[-1][-1]]) == signature: + runs[-1].append(state) + else: + runs.append([state]) + + representative = {state: run[0] for run in runs for state in run} + graph = TransitionGraph() + for run in runs: + graph.add_state(run[0]) + seen: set[tuple[Hashable, Any, Hashable]] = set() + for transition in shift.transitions(): + symbol = transition.data.get(ATTR_SYMBOL) + if symbol is None: + continue + edge = (representative[transition.source], symbol, representative[transition.target]) + if edge in seen: + continue + seen.add(edge) + graph.add_transition(edge[0], edge[2], **{ATTR_SYMBOL: symbol}) + return WheelerCover(graph=graph, symbol_alphabet=shift.symbol_alphabet) + + +def wheeler_index_of_shift(shift: SoficShift, **kwargs: Any) -> Any: + """Build a :class:`~sofic.automata.wheeler_index.WheelerIndex` over a Wheeler cover. + + Gives ``O(|w| log |A|)`` factor-language membership in place of scanning + :meth:`~sofic.shifts.base.SymbolicModel.factor_language`. + """ + from sofic.automata.wheeler_index import WheelerIndex + + return WheelerIndex.from_model(wheeler_cover(shift, **kwargs)) diff --git a/tests/test_wheeler.py b/tests/test_wheeler.py new file mode 100644 index 0000000..564a4b0 --- /dev/null +++ b/tests/test_wheeler.py @@ -0,0 +1,542 @@ +"""Tests for co-lexicographic (Wheeler) ordering.""" + +import itertools +import random +from collections import defaultdict + +import numpy as np +import pytest +from hypothesis import given, settings + +from sofic.automata.algorithms import equivalent +from sofic.automata.dfa import DFA +from sofic.automata.nfa import NFA +from sofic.automata.wheeler import ( + WheelerError, + check_wheeler_axioms, + colex_width, + is_input_consistent, + is_wheeler, + labeled_graph, + minimum_wdfa, + wheeler_canonical_form, + wheeler_isomorphic, + wheeler_order, + wheeler_state_index, + wnfa_to_wdfa, +) +from sofic.automata.wheeler_index import WheelerIndex +from sofic.examples.epsilon_machines import ( + butterfly_process, + even_process, + fair_coin, + golden_mean, + golden_mean_markov, + nemo_process, + wheeler_infinite_order_process, +) +from sofic.generators.synchronization import ( + _interval_power_automaton, + _subset_power_automaton, + cryptic_order_from_graph, + graph_from_epsilon_machine, + labeled_graph_of, + markov_order_from_graph, + reset_threshold_from_graph, +) +from sofic.generators.topological_epsilon_enumeration import iter_topological_epsilon_machines +from sofic.generators.wheeler_epsilon import ( + colex_cdf, + cylinder_measure, + debruijn_presentation, + wheeler_presentation, + wheeler_statistical_complexity, + word_cylinder_measure, +) +from sofic.shifts.sft import ShiftOfFiniteType +from sofic.shifts.sofic import SoficShift +from sofic.shifts.wheeler import ( + higher_block_presentation, + is_wheeler_shift, + wheeler_cover, + wheeler_order_of_shift, +) +from sofic.testing.strategies import dfas, epsilon_machines + +BINARY = frozenset({0, 1}) + + +def shift_from_edges(edges, alphabet=(0, 1)): + shift = SoficShift(symbol_alphabet=frozenset(alphabet)) + for state in {end for edge in edges for end in (edge[0], edge[2])}: + shift.graph.add_state(state) + for source, symbol, target in edges: + shift.add_transition(source, target, symbol) + return shift + + +GOLDEN_MEAN_SHIFT = ((0, 0, 0), (0, 1, 1), (1, 0, 0)) +EVEN_SHIFT = ((0, 0, 0), (0, 1, 1), (1, 1, 0)) + + +def sigma_star_dfa(order): + """Order-``order`` de Bruijn presentation of ``{a, b}*`` as a Wheeler DFA.""" + words = ["".join(w) for length in range(order + 1) for w in itertools.product("ab", repeat=length)] + dfa = DFA( + input_alphabet=frozenset("ab"), + initial_states=frozenset({""}), + accepting_states=frozenset(words), + ) + for word in words: + dfa.graph.add_state(word) + for word in words: + for symbol in "ab": + dfa.add_transition(word, (word + symbol)[-order:], symbol) + return dfa + + +def substring_wnfa(text): + """Substring automaton of ``text``: every state initial and accepting.""" + nfa = NFA( + input_alphabet=frozenset(text), + initial_states=frozenset(range(len(text) + 1)), + accepting_states=frozenset(range(len(text) + 1)), + ) + for position in range(len(text) + 1): + nfa.graph.add_state(position) + for position, symbol in enumerate(text): + nfa.add_transition(position, position + 1, symbol) + return nfa + + +# -- Known answers --------------------------------------------------------- + + +def test_golden_mean_is_wheeler_with_a_before_b(): + order = wheeler_order(golden_mean()) + assert order is not None + assert order.states == ("A", "B") + assert order.rank == {"A": 0, "B": 1} + + +@pytest.mark.parametrize("machine", [even_process(), nemo_process(), butterfly_process()]) +def test_canonical_non_wheeler_machines(machine): + assert not is_input_consistent(machine) + assert not is_wheeler(machine) + assert colex_width(machine) > 1 + + +@pytest.mark.parametrize("machine", [golden_mean(), golden_mean_markov()]) +def test_canonical_wheeler_machines(machine): + assert is_input_consistent(machine) + assert is_wheeler(machine) + assert colex_width(machine) == 1 + + +def test_wheeler_does_not_imply_finite_markov_order(): + machine = wheeler_infinite_order_process() + assert machine.wheeler_order().states == ("A", "E", "B", "C", "D") + assert machine.colex_width() == 1 + assert machine.markov_order() == float("inf") + + +def test_input_consistency_is_necessary_but_not_sufficient(): + # One state entered on two symbols cannot sit in any co-lex order, yet + # axiom one only compares *distinct* states, so width misses it. + full_shift = shift_from_edges(((0, 0, 0), (0, 1, 0))) + assert not is_input_consistent(full_shift) + assert not is_wheeler(full_shift) + assert colex_width(full_shift) == 1 + + +def test_wheeler_order_recovers_colex_order_of_words(): + shift = ShiftOfFiniteType.from_forbidden_words({(1, 1)}, BINARY).trim_transient() + order = wheeler_order(shift) + assert order is not None + assert list(order.states) == sorted(order.states, key=lambda word: tuple(reversed(word))) + + +# -- Enumeration regression ------------------------------------------------ + + +@pytest.mark.parametrize(("states", "expected"), [(1, 2), (2, 3), (3, 12), (4, 49)]) +def test_binary_wheeler_counts(states, expected): + found = sum(1 for m in iter_topological_epsilon_machines(2, states) if is_wheeler(m)) + assert found == expected + + +@pytest.mark.slow +def test_binary_wheeler_count_five_states(): + found = sum(1 for m in iter_topological_epsilon_machines(2, 5) if is_wheeler(m)) + assert found == 256 + + +@pytest.mark.slow +def test_wheeler_order_agrees_with_brute_force_search(): + for alphabet, states in ((2, 3), (2, 4), (3, 2), (3, 3)): + for machine in iter_topological_epsilon_machines(alphabet, states): + graph = labeled_graph(machine) + brute = any(check_wheeler_axioms(graph, p) for p in itertools.permutations(graph.states)) + assert brute == (wheeler_order(machine) is not None) + + +# -- Properties ------------------------------------------------------------ + + +@settings(deadline=None, max_examples=50) +@given(epsilon_machines(max_states=4)) +def test_returned_order_satisfies_the_axioms(machine): + order = wheeler_order(machine) + if order is not None: + assert check_wheeler_axioms(labeled_graph(machine), order.states) + assert sorted(order.states, key=repr) == sorted(machine.states(), key=repr) + + +@settings(deadline=None, max_examples=50) +@given(epsilon_machines(max_states=4)) +def test_width_one_iff_wheeler_for_input_consistent_machines(machine): + if is_input_consistent(machine): + assert (colex_width(machine) == 1) == is_wheeler(machine) + else: + assert not is_wheeler(machine) + + +@settings(deadline=None, max_examples=50) +@given(dfas(alphabet=(0, 1), max_states=4)) +def test_wheeler_implies_width_one(dfa): + if is_wheeler(dfa): + assert colex_width(dfa) == 1 + + +@settings(deadline=None, max_examples=30) +@given(epsilon_machines(max_states=4)) +def test_canonical_form_is_invariant_under_copying(machine): + if is_wheeler(machine): + assert wheeler_isomorphic(machine, machine.copy()) + assert wheeler_canonical_form(machine) == wheeler_canonical_form(machine.copy()) + + +# -- Interval power automaton ---------------------------------------------- + + +def wheeler_order_for(graph): + from sofic.automata.wheeler import wheeler_order_of_graph + + return wheeler_order_of_graph(labeled_graph_of(graph)) + + +def synchronization_orders(graph): + return ( + markov_order_from_graph(graph), + cryptic_order_from_graph(graph), + reset_threshold_from_graph(graph), + ) + + +def force_subset_power_automaton(monkeypatch): + """Make :func:`power_automaton` take the unrestricted subset branch.""" + monkeypatch.setattr( + "sofic.generators.synchronization.wheeler_order_of_graph", + lambda _graph, **_kwargs: None, + ) + + +@pytest.mark.parametrize("machine", [golden_mean(), golden_mean_markov(), wheeler_infinite_order_process()]) +def test_interval_and_subset_power_automata_agree(machine): + graph = graph_from_epsilon_machine(machine) + order = wheeler_order_for(graph) + assert order is not None + interval = _interval_power_automaton(graph, order) + subset = _subset_power_automaton(graph) + assert interval.start == subset.start + assert interval.transitions == subset.transitions + + +@pytest.mark.parametrize("machine", [golden_mean(), golden_mean_markov(), wheeler_infinite_order_process()]) +def test_synchronization_orders_are_unchanged_by_the_interval_fast_path(machine, monkeypatch): + graph = graph_from_epsilon_machine(machine) + assert wheeler_order_for(graph) is not None + with_intervals = synchronization_orders(graph) + force_subset_power_automaton(monkeypatch) + assert synchronization_orders(graph) == with_intervals + + +@settings(deadline=None, max_examples=40) +@given(epsilon_machines(max_states=4)) +def test_power_automaton_paths_agree_on_every_wheeler_machine(machine): + graph = graph_from_epsilon_machine(machine) + order = wheeler_order_for(graph) + if order is None: + return + interval = _interval_power_automaton(graph, order) + subset = _subset_power_automaton(graph) + assert interval.start == subset.start + assert interval.transitions == subset.transitions + + +def test_interval_power_automaton_is_polynomially_bounded(): + machine = wheeler_infinite_order_process() + graph = graph_from_epsilon_machine(machine) + size = len(graph.states) + pa = _subset_power_automaton(graph) + assert len(pa.transitions) <= size * (size + 1) // 2 + + +# -- Minimization and determinization -------------------------------------- + + +def test_minimum_wdfa_collapses_a_de_bruijn_presentation(): + dfa = sigma_star_dfa(2) + assert is_wheeler(dfa) + minimal = minimum_wdfa(dfa) + assert len(list(minimal.states())) == 3 + assert is_wheeler(minimal) + assert equivalent(dfa, minimal, frozenset("ab")) + # Already minimal, so minimizing again is a no-op. + assert len(list(minimum_wdfa(minimal).states())) == 3 + + +def test_minimum_wdfa_rejects_non_wheeler_input(): + dfa = DFA(input_alphabet=BINARY, initial_states=frozenset({"A"}), accepting_states=frozenset({"A"})) + for state in ("A", "B"): + dfa.graph.add_state(state) + dfa.add_transition("A", "A", 0) + dfa.add_transition("A", "B", 1) + dfa.add_transition("B", "A", 1) + with pytest.raises(WheelerError): + minimum_wdfa(dfa) + + +def test_wnfa_to_wdfa_stays_within_the_interval_bound(): + nfa = substring_wnfa("abra") + assert is_wheeler(nfa) + dfa = wnfa_to_wdfa(nfa) + states, arity = len(list(nfa.states())), len(nfa.input_alphabet) + assert len(list(dfa.states())) <= 2 * states - 1 - arity + assert is_wheeler(dfa) + assert equivalent(nfa, dfa, nfa.input_alphabet) + + +def test_wheeler_state_index_is_in_wheeler_order(): + index = wheeler_state_index(golden_mean()) + assert tuple(index.states) == ("A", "B") + + +# -- Succinct index -------------------------------------------------------- + + +def test_index_membership_agrees_with_the_factor_language(): + shift = ShiftOfFiniteType.from_forbidden_words({(1, 1)}, BINARY) + index = WheelerIndex.from_model(shift) + for length in range(1, 8): + expected = set(shift.factor_language(length)) + for word in itertools.product((0, 1), repeat=length): + assert index.contains(word) == (word in expected) + + +def test_index_membership_agrees_with_dfa_recognition(): + dfa = sigma_star_dfa(1) + dfa.accepting_states = frozenset({"a"}) + index = WheelerIndex.from_model(dfa) + for length in range(5): + for word in itertools.product("ab", repeat=length): + assert index.contains(word) == dfa.recognizes(word) + + +def test_index_finds_substrings(): + text = "abracadabra" + index = WheelerIndex.from_model(substring_wnfa(text)) + substrings = {text[i:j] for i in range(len(text) + 1) for j in range(i, len(text) + 1)} + for length in range(4): + for word in itertools.product("abcdrx", repeat=length): + assert index.contains(word) == ("".join(word) in substrings) + + +def test_index_ranks_and_unranks_in_colex_order(): + shift = ShiftOfFiniteType.from_forbidden_words({(1, 1)}, BINARY) + index = WheelerIndex.from_model(shift) + for length in range(1, 8): + words = sorted(shift.factor_language(length)) + assert index.count_words(length) == len(words) + assert list(index.words_of_length(length)) == sorted(words, key=lambda w: tuple(reversed(w))) + for word in words: + assert index.unrank_word(index.rank_word(word), length) == word + + +def test_index_word_counts_follow_the_golden_mean_fibonacci(): + index = WheelerIndex.from_model(shift_from_edges(GOLDEN_MEAN_SHIFT)) + assert [index.count_words(length) for length in range(1, 9)] == [2, 3, 5, 8, 13, 21, 34, 55] + + +def test_index_rejects_out_of_range_ranks_and_unknown_words(): + index = WheelerIndex.from_model(shift_from_edges(GOLDEN_MEAN_SHIFT)) + with pytest.raises(IndexError): + index.unrank_word(index.count_words(3), 3) + with pytest.raises(ValueError): + index.rank_word((1, 1)) + + +def test_index_samples_only_language_words(): + index = WheelerIndex.from_model(shift_from_edges(GOLDEN_MEAN_SHIFT)) + rng = random.Random(0) + for _draw in range(25): + assert index.contains(index.sample_word(6, rng)) + + +def test_index_forward_search_returns_reachable_states(): + index = WheelerIndex.from_model(shift_from_edges(GOLDEN_MEAN_SHIFT)) + assert index.states_reached((1,)) == (1,) + assert index.forward_search((1, 1)) is None + assert index.count_states(()) == 2 + + +# -- Shifts ---------------------------------------------------------------- + + +def test_wheeler_cover_of_the_golden_mean_shift(): + shift = shift_from_edges(GOLDEN_MEAN_SHIFT) + assert wheeler_order_of_shift(shift) == 0 + cover = wheeler_cover(shift) + assert len(list(cover.states())) == 2 + assert cover.is_wheeler() + + +def test_wheeler_cover_adds_memory_for_the_full_shift(): + full_shift = shift_from_edges(((0, 0, 0), (0, 1, 0))) + assert wheeler_order_of_shift(full_shift) == 1 + cover = wheeler_cover(full_shift) + assert len(list(cover.states())) == 2 + assert cover.is_wheeler() + + +def test_the_even_shift_has_no_wheeler_presentation(): + # Wheeler languages are star-free; the even shift counts 1s modulo two. + even = shift_from_edges(EVEN_SHIFT) + assert wheeler_order_of_shift(even) is None + assert not is_wheeler_shift(even) + with pytest.raises(WheelerError): + wheeler_cover(even) + + +@pytest.mark.parametrize("edges", [GOLDEN_MEAN_SHIFT, ((0, 0, 0), (0, 1, 0))]) +def test_wheeler_cover_preserves_the_factor_language(edges): + shift = shift_from_edges(edges) + cover = wheeler_cover(shift) + for length in range(1, 8): + assert sorted(cover.factor_language(length)) == sorted(shift.factor_language(length)) + + +def test_wheeler_cover_right_resolves_a_left_resolving_presentation(): + reversed_golden_mean = shift_from_edges(tuple((t, a, s) for s, a, t in GOLDEN_MEAN_SHIFT)) + cover = wheeler_cover(reversed_golden_mean) + assert cover.is_wheeler() + + +def test_higher_block_presentation_is_input_consistent(): + even = shift_from_edges(EVEN_SHIFT) + for order in (1, 2, 3): + assert is_input_consistent(higher_block_presentation(even, order)) + + +# -- Stochastic layer ------------------------------------------------------ + + +@pytest.mark.parametrize("machine", [golden_mean(), golden_mean_markov(), wheeler_infinite_order_process()]) +def test_wheeler_presentation_is_the_machine_itself_when_wheeler(machine): + assert wheeler_presentation(machine) is machine + assert wheeler_statistical_complexity(machine) == pytest.approx(machine.statistical_complexity()) + + +@settings(deadline=None, max_examples=30) +@given(epsilon_machines(max_states=4)) +def test_wheeler_complexity_is_never_below_statistical_complexity(machine): + try: + found = wheeler_statistical_complexity(machine) + except WheelerError: + return + assert found >= machine.statistical_complexity() - 1e-9 + + +def stationary_symbol_entropy(presentation): + """``H[X_0]`` read off an input-consistent presentation's incoming labels.""" + index = presentation.reindex() + stationary = np.asarray(presentation.stationary_distribution(), dtype=float) + in_labels = labeled_graph(presentation).in_labels() + mass = defaultdict(float) + for state in presentation.states(): + for symbol in in_labels[state]: + mass[symbol] += stationary[index.index(state)] + probabilities = np.array(list(mass.values())) + probabilities = probabilities[probabilities > 0] + return float(-(probabilities * np.log2(probabilities)).sum()) + + +@settings(deadline=None, max_examples=30) +@given(epsilon_machines(max_states=4)) +def test_wheeler_complexity_is_never_below_the_single_symbol_entropy(machine): + # A Wheeler presentation is input consistent, so its state determines the + # symbol that entered it and C_W >= H[X_0] on top of C_W >= C_mu. + if not machine.is_irreducible(): + return + try: + presentation = wheeler_presentation(machine) + except WheelerError: + return + found = wheeler_statistical_complexity(machine) + assert found >= stationary_symbol_entropy(presentation) - 1e-9 + + +def test_a_fair_coin_pays_a_whole_bit_for_sortability(): + coin = fair_coin() + assert coin.statistical_complexity() == pytest.approx(0.0) + assert wheeler_statistical_complexity(coin) == pytest.approx(1.0) + + +def test_every_finite_order_process_has_a_wheeler_presentation(): + # Finite Markov order R means the length-R words pin down the causal state, + # so the de Bruijn presentation exists and its states sort by construction. + # Machines whose own causal states are not sortable pay for it in C_W. + checked = 0 + for states in (2, 3): + for machine in iter_topological_epsilon_machines(2, states): + if is_wheeler(machine) or machine.markov_order() == float("inf"): + continue + presentation = wheeler_presentation(machine) + assert is_wheeler(presentation) + assert wheeler_statistical_complexity(machine) > machine.statistical_complexity() + checked += 1 + assert checked > 0 + + +def test_debruijn_presentation_needs_the_markov_order(): + machine = next(m for m in iter_topological_epsilon_machines(2, 3) if not is_wheeler(m) and m.markov_order() == 3) + assert debruijn_presentation(machine, 1) is None + assert is_wheeler(debruijn_presentation(machine, 3)) + + +@pytest.mark.parametrize("machine", [even_process(), nemo_process(), butterfly_process()]) +def test_non_star_free_processes_have_no_wheeler_presentation(machine): + with pytest.raises(WheelerError): + wheeler_presentation(machine) + + +def test_colex_cdf_is_the_cumulative_stationary_distribution(): + machine = golden_mean() + states, cumulative = colex_cdf(machine) + assert states == ("A", "B") + assert cumulative[0] == pytest.approx(2 / 3) + assert cumulative[-1] == pytest.approx(1.0) + assert cylinder_measure(machine, 0, 1) == pytest.approx(1.0) + assert cylinder_measure(machine, 1, 1) == pytest.approx(1 / 3) + + +def test_word_cylinder_measure_weighs_the_reachable_interval(): + machine = golden_mean() + assert word_cylinder_measure(machine, (0,)) == pytest.approx(2 / 3) + assert word_cylinder_measure(machine, (1,)) == pytest.approx(1 / 3) + assert word_cylinder_measure(machine, (1, 1)) == 0.0 + + +def test_colex_cdf_rejects_non_wheeler_machines(): + with pytest.raises(WheelerError): + colex_cdf(even_process()) From 4d0a927895f2e2350bb5e8eb444808750b07c81d Mon Sep 17 00:00:00 2001 From: Ryan James Date: Thu, 17 Sep 2026 18:02:55 -0600 Subject: [PATCH 2/5] Untrack .DS_Store files committed by mistake --- docs/.DS_Store | Bin 8196 -> 0 bytes sofic/.DS_Store | Bin 8196 -> 0 bytes 2 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 docs/.DS_Store delete mode 100644 sofic/.DS_Store diff --git a/docs/.DS_Store b/docs/.DS_Store deleted file mode 100644 index 3f31ffd658bdf9d36fe56b49bf56c4e0b66df78d..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8196 zcmeI1y=xRf7{=ex_-L6my!>gR(=FJz^*` z9Ql~#h>4AviDoz{Gn|xNS=kec(yL=!tmUK<6Qvag#DSs%e0HzVW9ra^w!`zg`u=bw z%n#F6t20PD=&#?~e|zO!V|u%&>CM%Qw(|D`>gOS*pmX5Zr9En+Zd2EPld$l*sBl>h zp4A>q|9oYZV4h)`(^RWm9(xbF0G>VS(~w4RhN=3r>r<&mSg%y)vfXOL_bpq6HYc2? zE#a(#Z34O!q)gfN8K3mPSkSn9@4x)IcyL+oV~H+jJ5J|unY(aAT`-NQ@AE8v1J5ky za(q~AHm^sT!^!9Cm*u^6X&qBzKQ^&bJE&vr2Qviw_WNU%&gFUeJ5)$jjL(pQ%6A`XZHfdeXPt+n`vzCBw9X>-bt(YtOltjkal^$;}r mIF1;x1abC{e5t(WE$JMb5(+DRt> diff --git a/sofic/.DS_Store b/sofic/.DS_Store deleted file mode 100644 index 055b6da0147b3cbe585cf94b5ee211b2de1295e3..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8196 zcmeI1L2DC16vy9~T^p^S9!#qstcRcn5p(jomh=`~-Ry z@A?hAh!XV zKpgll9pL9fgw8RxF*8v=Ixy%G05QO6-EfRLK<7j@#x`aq>Qt<0st09Ll`S!pNymJ| za>Uri%tVt;%A}LBJ1bkED7`!KA}uEsn<%3=APy8Ah;V8i)q0CY^w_K4Tf=^1rJE)r zom6>nr(Bt&&F0!LU4x?f>3jXz%i!JJS?cw_P>%)=5(`k9@mx3aP{aBh!S`?#5a&lk1MS2Xs!v3%XRVR+rj=J zS3#F^&C+>X;U06Hr9Noxqd(MsrHDS;wf7&&bS~e)<(=*G3);7-Yr5Y3N$YXee=+A7 z+5qJ}&B?r7>UjleVkUa4?g5Br&jJyAf+S9{($mY4%q&VyC7k<@$j9_JX|FiSb#DU{GP%){CT>meZ zfB!%J29hZbhyy3y0Tnb?noGd6wRHiQYwZZVjn0jAnTa|D155b4F2@m%{$YsY$j-$! UW+q|=O+Ex@8)Og%{;31M0dxr*OaK4? From fab49a05e1853bfe9ca19169ddb86e5b3cb27685 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Thu, 17 Sep 2026 18:10:08 -0600 Subject: [PATCH 3/5] Note that five states is the smallest Wheeler/infinite-order size --- sofic/.DS_Store | Bin 0 -> 8196 bytes sofic/examples/epsilon_machines.py | 5 +++-- 2 files changed, 3 insertions(+), 2 deletions(-) create mode 100644 sofic/.DS_Store diff --git a/sofic/.DS_Store b/sofic/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..055b6da0147b3cbe585cf94b5ee211b2de1295e3 GIT binary patch literal 8196 zcmeI1L2DC16vy9~T^p^S9!#qstcRcn5p(jomh=`~-Ry z@A?hAh!XV zKpgll9pL9fgw8RxF*8v=Ixy%G05QO6-EfRLK<7j@#x`aq>Qt<0st09Ll`S!pNymJ| za>Uri%tVt;%A}LBJ1bkED7`!KA}uEsn<%3=APy8Ah;V8i)q0CY^w_K4Tf=^1rJE)r zom6>nr(Bt&&F0!LU4x?f>3jXz%i!JJS?cw_P>%)=5(`k9@mx3aP{aBh!S`?#5a&lk1MS2Xs!v3%XRVR+rj=J zS3#F^&C+>X;U06Hr9Noxqd(MsrHDS;wf7&&bS~e)<(=*G3);7-Yr5Y3N$YXee=+A7 z+5qJ}&B?r7>UjleVkUa4?g5Br&jJyAf+S9{($mY4%q&VyC7k<@$j9_JX|FiSb#DU{GP%){CT>meZ zfB!%J29hZbhyy3y0Tnb?noGd6wRHiQYwZZVjn0jAnTa|D155b4F2@m%{$YsY$j-$! UW+q|=O+Ex@8)Og%{;31M0dxr*OaK4? literal 0 HcmV?d00001 diff --git a/sofic/examples/epsilon_machines.py b/sofic/examples/epsilon_machines.py index d41c905..b30e777 100644 --- a/sofic/examples/epsilon_machines.py +++ b/sofic/examples/epsilon_machines.py @@ -1073,8 +1073,9 @@ def wheeler_infinite_order_process(p: float = 0.5, q: float = 0.5, r: float = 0. Found by exhaustive search over binary topological ε-machines with :func:`~sofic.generators.topological_epsilon_enumeration.iter_topological_epsilon_machines`; - twenty of the 35186 five-state machines share both properties. No prior - source states this example. + twenty of the 35186 five-state machines share both properties, and five + states is the smallest size at which any does. No prior source states this + example. """ for name, value in (("p", p), ("q", q), ("r", r)): if not 0.0 < value < 1.0: From b64e9e6879e98a7400d27d69ae24eba5121d5377 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Thu, 17 Sep 2026 18:34:07 -0600 Subject: [PATCH 4/5] Attribute Wheeler language recognition to its actual sources Becker et al. is 2023, not 2024, and the PSPACE-completeness of deciding whether an NFA's language is Wheeler is D'Agostino, Martincigh, and Policriti, not Alanko et al. Add the missing DAgostino2023 entry and correct the claim in both the is_wheeler docstring and the automata docs. Co-authored-by: Cursor --- docs/automata/wheeler.rst | 12 +++++++----- docs/references.bib | 12 ++++++++++++ sofic/automata/wheeler.py | 4 ++-- 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/docs/automata/wheeler.rst b/docs/automata/wheeler.rst index e957b90..5274d3f 100644 --- a/docs/automata/wheeler.rst +++ b/docs/automata/wheeler.rst @@ -26,11 +26,13 @@ determinizing it. Presentations versus languages ============================== -:func:`is_wheeler` asks whether *this presentation* is Wheeler, which is a -sorting question answerable in polynomial time. Whether the *language* is -Wheeler — whether any equivalent automaton is — is much harder: ``O(mn)`` for a -DFA :cite:`Becker2023`, PSPACE-complete for an NFA :cite:`Alanko2021`, and -NP-hard to recognize as a graph property :cite:`GibneyThankachan2019`. +:func:`is_wheeler` asks whether *this presentation* is Wheeler. For a +deterministic presentation that is a sorting question, answerable in polynomial +time; for a general labeled graph, recognizing Wheelerness is NP-complete +:cite:`GibneyThankachan2019`. Whether the *language* is Wheeler — whether any +equivalent automaton is — is harder still: ``O(mn)`` for a DFA +:cite:`Becker2023`, improving the first polynomial algorithm +:cite:`Alanko2021`, and PSPACE-complete for an NFA :cite:`DAgostino2023`. :func:`~sofic.shifts.wheeler.wheeler_cover` and :func:`~sofic.generators.wheeler_epsilon.wheeler_presentation` search over presentations for the process-level question. diff --git a/docs/references.bib b/docs/references.bib index 44fb110..8cab925 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -990,3 +990,15 @@ @article{ShyrThierrin1974 pages = {9--20}, year = {1974}, } + +@article{DAgostino2023, + author = {D'Agostino, Giovanna and Martincigh, Davide and Policriti, Alberto}, + title = {Ordering Regular Languages and Automata: Complexity}, + journal = {Theoretical Computer Science}, + volume = {949}, + pages = {113709}, + year = {2023}, + doi = {10.1016/j.tcs.2023.113709}, + eprint = {2203.12534}, + archivePrefix = {arXiv}, +} diff --git a/sofic/automata/wheeler.py b/sofic/automata/wheeler.py index 7427cc2..45b28ee 100644 --- a/sofic/automata/wheeler.py +++ b/sofic/automata/wheeler.py @@ -361,8 +361,8 @@ def is_wheeler(model: Any, *, symbol_key: Callable[[Any], Any] | None = None) -> A property of the presentation. Whether the *language* is Wheeler -- whether some equivalent automaton is -- is a strictly weaker and much more expensive - question, decidable in ``O(mn)`` for a DFA and PSPACE-complete for an NFA - :cite:`Becker2023` :cite:`Alanko2021`. See + question, decidable in ``O(mn)`` for a DFA :cite:`Becker2023` and + PSPACE-complete for an NFA :cite:`DAgostino2023`. See :func:`~sofic.shifts.wheeler.wheeler_cover` and :func:`~sofic.generators.wheeler_epsilon.wheeler_presentation` for the search over presentations. From 926a3d85836c32296420c77d039d95bc271d8452 Mon Sep 17 00:00:00 2001 From: Ryan James Date: Thu, 17 Sep 2026 18:55:45 -0600 Subject: [PATCH 5/5] Relock after matplotlib moved into the viz extra pyproject declares matplotlib>=3.5 under viz, but uv.lock predated that, so every uv invocation re-resolved the project before running. No package versions change; the lock only records the viz dependency and its specifier. Co-authored-by: Cursor --- uv.lock | 3 +++ 1 file changed, 3 insertions(+) diff --git a/uv.lock b/uv.lock index 2d87508..3f01177 100644 --- a/uv.lock +++ b/uv.lock @@ -2170,6 +2170,7 @@ test = [ ] viz = [ { name = "graphviz" }, + { name = "matplotlib" }, ] [package.metadata] @@ -2189,7 +2190,9 @@ requires-dist = [ { name = "ipython", marker = "extra == 'dev'" }, { name = "ipython", marker = "extra == 'docs'" }, { name = "matplotlib", marker = "extra == 'dev'" }, + { name = "matplotlib", marker = "extra == 'dev'", specifier = ">=3.5" }, { name = "matplotlib", marker = "extra == 'docs'" }, + { name = "matplotlib", marker = "extra == 'viz'", specifier = ">=3.5" }, { name = "networkx", specifier = ">=2.6" }, { name = "numpy", specifier = ">=1.22" }, { name = "pymc", marker = "extra == 'bayes'", specifier = ">=5" },