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..5274d3f --- /dev/null +++ b/docs/automata/wheeler.rst @@ -0,0 +1,150 @@ +.. 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. 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. + +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 227ccb6..c4fa183 100644 --- a/docs/references.bib +++ b/docs/references.bib @@ -902,6 +902,107 @@ @inproceedings{Volkov2008 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}, +} + +@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}, +} + @misc{TetrisWikiNES, author = {{Tetris Wiki}}, title = {Tetris ({NES})}, 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 0000000..055b6da Binary files /dev/null and b/sofic/.DS_Store differ 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..45b28ee --- /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 :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. + """ + 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 2900ba2..7f57f01 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 @@ -98,6 +99,7 @@ "tetris_nes", "tetris_tgm", "tetris_tgm2", + "wheeler_infinite_order_process", ] __all__ += _process_all __all__ = [name for name in __all__ if name not in {"processes", "shifts", "tetris"}] diff --git a/sofic/examples/epsilon_machines.py b/sofic/examples/epsilon_machines.py index c999185..b30e777 100644 --- a/sofic/examples/epsilon_machines.py +++ b/sofic/examples/epsilon_machines.py @@ -1062,6 +1062,51 @@ 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, 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: + 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()) 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" },