diff --git a/docs/automata/dfa.rst b/docs/automata/dfa.rst index 23a0a49..70b4462 100644 --- a/docs/automata/dfa.rst +++ b/docs/automata/dfa.rst @@ -39,4 +39,4 @@ API === .. autoclass:: DFA - :members: add_transition, recognizes, union, intersection, intersect, complement, difference, concat, concatenate, kleene_star, star, minimize, from_nfa + :members: add_transition, recognizes, union, intersection, complement, difference, concat, kleene_star, minimize, from_nfa diff --git a/docs/automata/nfa.rst b/docs/automata/nfa.rst index 67f6c30..64983e6 100644 --- a/docs/automata/nfa.rst +++ b/docs/automata/nfa.rst @@ -31,4 +31,4 @@ API === .. autoclass:: NFA - :members: add_transition, recognizes, union, intersection, intersect, complement, difference, concat, concatenate, kleene_star, star, determinize, minimize + :members: add_transition, recognizes, union, intersection, complement, difference, concat, kleene_star, determinize, minimize diff --git a/docs/examples.rst b/docs/examples.rst index 566f186..4c749b1 100644 --- a/docs/examples.rst +++ b/docs/examples.rst @@ -181,6 +181,19 @@ but accepting a ``machine_type`` argument: In [3]: gm.entropy_rate() Out[3]: 0.6666666666666665 +Several factories share one implementation with a curated example and differ only +in string symbols (``"0"``/``"1"``), state names, or parametrization: +``BiasedCoin(b)`` is :func:`bernoulli` ``(b)``, ``Even`` is :func:`even_process`, +``Nemo`` is :func:`nemo_process`, ``GoldenMean(b)`` is +:func:`golden_mean_forward` ``(1 - b)`` (and :func:`golden_mean_markov` ``(b)`` +with ``A``/``B`` swapped), ``ABC(p, q)`` is :func:`alternating_biased_coins` +``(1 - p, 1 - q)`` for ``p != q``, ``RestrictedGM`` is +:func:`restricted_golden_mean`, and ``IrreversibleTwoState()`` is +:func:`ellison_fig9_forward`. Similar names do not always mean the same process: +:func:`golden_mean` is the ``0 <-> 1`` mirror of ``GoldenMean``, and ``Butterfly`` +and ``PSB`` differ from :func:`butterfly_process` and +:func:`~sofic.examples.epsilon_machines.phase_slip_backtrack`. + The module also exposes registries — ``processes.process_list`` and ``processes.transducer_list`` — that enumerate every factory, which is handy for parametrized tests and sweeps: diff --git a/docs/generators/epsilon_inference.rst b/docs/generators/epsilon_inference.rst index a9494fb..f4808b3 100644 --- a/docs/generators/epsilon_inference.rst +++ b/docs/generators/epsilon_inference.rst @@ -76,7 +76,11 @@ it as the longest history the data supports, not as a synchronization length. The morph tests also rely on the chi-squared limit, which fails for the sparse counts of long suffixes. ``test="exact"`` compares the G statistic with tables drawn uniformly given the observed margins whenever an expected count is below 5 -(seeded from the table, so reconstruction stays deterministic). +(seeded from the table, so reconstruction stays deterministic). The G-test +(``test="g"``) applies Yates' continuity correction when the table has one degree +of freedom (two observed symbols), matching +``scipy.stats.chi2_contingency(table, lambda_="log-likelihood")``; transCSSR and +stack CSSR share the same implementation. ``correction="bonferroni"`` divides ``alpha`` by the number of suffixes eligible for testing, bounding the chance of any spurious split. Because CSSR chooses each test in light of earlier outcomes, false-discovery-rate step-up procedures do not diff --git a/docs/generators/epsilon_transducer_inference.rst b/docs/generators/epsilon_transducer_inference.rst index 707be51..da3280e 100644 --- a/docs/generators/epsilon_transducer_inference.rst +++ b/docs/generators/epsilon_transducer_inference.rst @@ -13,7 +13,9 @@ channel :cite:`Barnett2015`. Causal states are equivalence classes of joint ``(input, output)`` pasts that induce the same conditional next-output law ``P(y | history, x)`` for every input symbol ``x``. Rare histories inherit their parent's state (controlled by -``min_count``); the split decision uses a G-test at significance ``alpha``. +``min_count``); the split decision uses a G-test at significance ``alpha``, the +same test as :func:`~sofic.generators.epsilon_inference.morphs_differ` (with +Yates' continuity correction at one degree of freedom). As in :func:`~sofic.generators.epsilon_inference.cssr`, ``test="exact"`` uses a Monte Carlo exact G-test for tables with small expected counts, and ``correction="bonferroni"`` divides ``alpha`` by the number of diff --git a/docs/generators/hmm_inference.rst b/docs/generators/hmm_inference.rst index c430f29..f4ae796 100644 --- a/docs/generators/hmm_inference.rst +++ b/docs/generators/hmm_inference.rst @@ -99,3 +99,27 @@ Score and observed information .. autofunction:: observed_information .. autofunction:: free_parameter_labels .. autofunction:: standard_errors + +Transition matrices and start policies +-------------------------------------- + +.. py:module:: sofic.generators.matrices + :no-index: + +Every routine above works with the symbol-labeled joint transition matrices +:math:`T^{(x)}_{ij} = P(S_{t+1} = j, X_t = x \mid S_t = i)` of the model's +Mealy presentation :cite:`Rabiner1989,Ellison2009`, built in one place by +:mod:`sofic.generators.matrices`. The initial state law follows one of two +policies: + +- ``"model"`` (likelihoods, decoding, sampling, word probabilities): the + model's ``initial_distribution``, or its stationary distribution when none is + given. +- ``"stationary"`` (block and window statistics such as + ``joint_block_distribution`` and directional flow): the stationary process + law, preferring the limit reached from ``initial_distribution`` on reducible + chains. + +.. autofunction:: sofic.generators.matrices.start_vector +.. autofunction:: sofic.generators.matrices.emission_tensors +.. autofunction:: sofic.generators.matrices.symbol_matrices diff --git a/sofic/automata/_words.py b/sofic/automata/_words.py new file mode 100644 index 0000000..1f631a9 --- /dev/null +++ b/sofic/automata/_words.py @@ -0,0 +1,23 @@ +"""Shared word enumeration for automata modules.""" + +from __future__ import annotations + +from collections.abc import Iterable, Iterator +from typing import Any + +Word = tuple[Any, ...] + + +def _words_up_to(max_length: int, alphabet: Iterable[Any]) -> Iterator[Word]: + """Yield every word of length at most ``max_length`` in breadth-first order.""" + symbols = tuple(alphabet) + frontier: list[Word] = [()] + yield () + for _ in range(max_length): + nxt: list[Word] = [] + for word in frontier: + for symbol in symbols: + extended = (*word, symbol) + yield extended + nxt.append(extended) + frontier = nxt diff --git a/sofic/automata/active.py b/sofic/automata/active.py index 01393f3..55fc281 100644 --- a/sofic/automata/active.py +++ b/sofic/automata/active.py @@ -24,6 +24,7 @@ import numpy as np +from sofic.automata._words import _words_up_to from sofic.automata.dfa import DFA from sofic.automata.transducers import MealyMachine @@ -135,32 +136,50 @@ def output(self, word: Sequence[Any]) -> Word: return next(iter(outputs)) -def _words_up_to(max_length: int, alphabet: Sequence[Any]) -> Iterator[Word]: - frontier: list[Word] = [()] - yield () - for _ in range(max_length): - nxt: list[Word] = [] - for word in frontier: - for symbol in alphabet: - extended = (*word, symbol) - yield extended - nxt.append(extended) - frontier = nxt +class _BoundedEquivalenceOracle: + """Approximate equivalence test over a finite stream of candidate words. + + Subclasses supply the candidate words (:meth:`_words`) and the comparison of + target versus hypothesis on one word (:meth:`_disagrees`); the first + disagreeing word is the counterexample. + """ + + def __init__(self, alphabet: Iterable[Any]) -> None: + self._alphabet = tuple(sorted(alphabet, key=repr)) + + def _words(self) -> Iterator[Word]: + raise NotImplementedError + + def _disagrees(self, hypothesis: Any, word: Word) -> bool: + raise NotImplementedError + + def find_counterexample(self, hypothesis: Any) -> Word | None: + for word in self._words(): + if self._disagrees(hypothesis, word): + return word + return None + + +class _MembershipBoundedOracle(_BoundedEquivalenceOracle): + """Bounded oracle comparing a membership oracle against a DFA hypothesis.""" + + def __init__(self, membership: MembershipOracle, alphabet: Iterable[Any]) -> None: + super().__init__(alphabet) + self._membership = membership + + def _disagrees(self, hypothesis: DFA, word: Word) -> bool: + return self._membership.member(word) != hypothesis.recognizes(word) -class ExhaustiveEquivalenceOracle: +class ExhaustiveEquivalenceOracle(_MembershipBoundedOracle): """Bounded exhaustive equivalence test for DFA hypotheses.""" def __init__(self, membership: MembershipOracle, alphabet: Iterable[Any], *, max_length: int = 10) -> None: - self._membership = membership - self._alphabet = tuple(sorted(alphabet, key=repr)) + super().__init__(membership, alphabet) self._max_length = int(max_length) - def find_counterexample(self, hypothesis: DFA) -> Word | None: - for word in _words_up_to(self._max_length, self._alphabet): - if self._membership.member(word) != hypothesis.recognizes(word): - return word - return None + def _words(self) -> Iterator[Word]: + return _words_up_to(self._max_length, self._alphabet) class AutomatonEquivalenceOracle: @@ -206,7 +225,7 @@ def _step(aut: Any, states: frozenset[Any], symbol: Any) -> frozenset[Any]: return _closure(aut, targets) -class RandomWalkEquivalenceOracle: +class RandomWalkEquivalenceOracle(_MembershipBoundedOracle): """Randomized equivalence test drawing random input words for DFA hypotheses.""" def __init__( @@ -218,39 +237,33 @@ def __init__( max_steps: int = 30, rng: np.random.Generator | int | None = None, ) -> None: - self._membership = membership - self._alphabet = tuple(sorted(alphabet, key=repr)) + super().__init__(membership, alphabet) self._num_walks = int(num_walks) self._max_steps = int(max_steps) self._rng = rng if isinstance(rng, np.random.Generator) else np.random.default_rng(rng) - def find_counterexample(self, hypothesis: DFA) -> Word | None: + def _words(self) -> Iterator[Word]: n_symbols = len(self._alphabet) for _ in range(self._num_walks): length = int(self._rng.integers(0, self._max_steps + 1)) - word = tuple(self._alphabet[int(self._rng.integers(0, n_symbols))] for _ in range(length)) - if self._membership.member(word) != hypothesis.recognizes(word): - return word - return None + yield tuple(self._alphabet[int(self._rng.integers(0, n_symbols))] for _ in range(length)) -class MealyExhaustiveEquivalenceOracle: +class MealyExhaustiveEquivalenceOracle(_BoundedEquivalenceOracle): """Bounded exhaustive equivalence test for Mealy hypotheses.""" def __init__(self, oracle: MealyMembershipOracle, alphabet: Iterable[Any], *, max_length: int = 10) -> None: + super().__init__(alphabet) self._oracle = oracle - self._alphabet = tuple(sorted(alphabet, key=repr)) self._max_length = int(max_length) - def find_counterexample(self, hypothesis: MealyMachine) -> Word | None: - for word in _words_up_to(self._max_length, self._alphabet): - if not word: - continue - produced = hypothesis.transduce(word) - hyp_out = next(iter(produced)) if produced else None - if self._oracle.output(word) != hyp_out: - return word - return None + def _words(self) -> Iterator[Word]: + return (word for word in _words_up_to(self._max_length, self._alphabet) if word) + + def _disagrees(self, hypothesis: MealyMachine, word: Word) -> bool: + produced = hypothesis.transduce(word) + hyp_out = next(iter(produced)) if produced else None + return self._oracle.output(word) != hyp_out # ------------------------------------------------------------------------------ L* @@ -452,10 +465,6 @@ def build() -> DFA: raise RuntimeError("TTT did not converge within max_rounds; check the equivalence oracle") -def _hypothesis_access(word: Word, sift: Callable[[Word], _DTNode]) -> Word: - return sift(word).access - - def _split_leaf( counterexample: Word, sift: Callable[[Word], _DTNode], @@ -466,7 +475,7 @@ def _split_leaf( length = len(counterexample) def alpha(index: int) -> Word: - return _hypothesis_access(counterexample[:index], sift) + counterexample[index:] + return sift(counterexample[:index]).access + counterexample[index:] base = member(alpha(0)) breakpoint_index = None @@ -477,7 +486,7 @@ def alpha(index: int) -> Word: if breakpoint_index is None: # pragma: no cover - guaranteed by a valid counterexample raise RuntimeError("counterexample analysis found no breakpoint") - state_access = _hypothesis_access(counterexample[:breakpoint_index], sift) + state_access = sift(counterexample[:breakpoint_index]).access symbol = counterexample[breakpoint_index] discriminator = counterexample[breakpoint_index + 1 :] new_access = state_access + (symbol,) diff --git a/sofic/automata/base.py b/sofic/automata/base.py index f107546..d75e0a9 100644 --- a/sofic/automata/base.py +++ b/sofic/automata/base.py @@ -79,10 +79,6 @@ def intersection(self, other: LabeledAutomaton) -> DFA: return intersection_dfa(self, other) - def intersect(self, other: LabeledAutomaton) -> DFA: - """Alias for :meth:`intersection`.""" - return self.intersection(other) - def complement(self, alphabet: frozenset[Any] | None = None) -> DFA: """Return a complete DFA recognizing the complement over ``alphabet``.""" from sofic.automata.languages.automaton_ops import complement_dfa @@ -101,20 +97,12 @@ def concat(self, other: LabeledAutomaton) -> NFA: return concat_nfa(self, other) - def concatenate(self, other: LabeledAutomaton) -> NFA: - """Alias for :meth:`concat`.""" - return self.concat(other) - def kleene_star(self) -> NFA: """Return an NFA recognizing the Kleene star of this language.""" from sofic.automata.languages.automaton_ops import kleene_star_nfa return kleene_star_nfa(self) - def star(self) -> NFA: - """Alias for :meth:`kleene_star`.""" - return self.kleene_star() - def is_deterministic(self) -> bool: """Return whether this automaton is DFA-deterministic.""" from sofic.properties import is_deterministic_automaton @@ -140,14 +128,7 @@ def epsilon_closure(self, states: set[Hashable]) -> set[Hashable]: return closure def _run_nfa(self, word: Sequence[Any], start: set[Hashable] | None = None) -> set[Hashable]: - seed = set(self.initial_states) if start is None else set(start) - current = self.epsilon_closure(seed) - for symbol in word: - next_states: set[Hashable] = set() - for state in current: - next_states.update(self.delta(state, symbol)) - current = self.epsilon_closure(next_states) - return current + return run_nfa(self, word, start=start) def reverse(self) -> NFA: """Return an NFA recognizing the reversed language.""" @@ -159,3 +140,22 @@ def reverse(self) -> NFA: accepting_states=frozenset(self.epsilon_closure(set(self.initial_states))), graph=self.graph.reverse(), ) + + +def run_nfa( + aut: LabeledAutomaton, + word: Sequence[Any], + start: set[Hashable] | frozenset[Hashable] | None = None, +) -> set[Hashable]: + """Return the epsilon-closed state set reached by reading ``word``. + + Simulation starts from ``start`` (default: ``aut.initial_states``). + """ + seed = set(aut.initial_states) if start is None else set(start) + current = aut.epsilon_closure(seed) + for symbol in word: + next_states: set[Hashable] = set() + for state in current: + next_states.update(aut.delta(state, symbol)) + current = aut.epsilon_closure(next_states) + return current diff --git a/sofic/automata/enumeration.py b/sofic/automata/enumeration.py index 3da18f2..c5353f7 100644 --- a/sofic/automata/enumeration.py +++ b/sofic/automata/enumeration.py @@ -6,6 +6,7 @@ from itertools import product from typing import Any +from sofic.automata.algorithms import _effective_alphabet from sofic.automata.base import LabeledAutomaton @@ -36,9 +37,3 @@ def iter_language( while max_length is None or length <= max_length: yield from words_of_length(automaton, length) length += 1 - - -def _effective_alphabet(automaton: LabeledAutomaton) -> frozenset[Any]: - from sofic.automata.algorithms import _effective_alphabet as _shared - - return _shared(automaton) diff --git a/sofic/automata/languages/_quotient_utils.py b/sofic/automata/languages/_quotient_utils.py index ec136da..e4f524d 100644 --- a/sofic/automata/languages/_quotient_utils.py +++ b/sofic/automata/languages/_quotient_utils.py @@ -9,15 +9,6 @@ from sofic.automata.languages.base import AutomatonLanguage, ExplicitLanguage -def _words_up_to(length: int, alphabet: frozenset[Any]) -> list[tuple[Any, ...]]: - if length < 0: - return [] - if length == 0: - return [()] - shorter = _words_up_to(length - 1, alphabet) - return [word + (symbol,) for word in shorter for symbol in alphabet] - - def _residual_from_state(aut: AutomatonLanguage, state) -> AutomatonLanguage: sub = aut.automaton.copy() sub.initial_states = frozenset({state}) diff --git a/sofic/automata/languages/automaton_ops.py b/sofic/automata/languages/automaton_ops.py index 593e472..c119d57 100644 --- a/sofic/automata/languages/automaton_ops.py +++ b/sofic/automata/languages/automaton_ops.py @@ -6,8 +6,8 @@ from collections.abc import Callable, Hashable, Sequence from typing import Any -from sofic.automata.algorithms import complete, determinize, minimize, trim -from sofic.automata.base import LabeledAutomaton +from sofic.automata.algorithms import _effective_alphabet, complete, determinize, minimize, trim +from sofic.automata.base import LabeledAutomaton, run_nfa from sofic.automata.dfa import DFA from sofic.automata.nfa import NFA from sofic.graph import ATTR_SYMBOL, EPSILON, TransitionGraph @@ -202,12 +202,7 @@ def kleene_star_nfa(aut: LabeledAutomaton) -> NFA: def left_quotient_automaton(u: Sequence[Any], aut: LabeledAutomaton) -> NFA: nfa = _to_nfa(aut) - current = nfa.epsilon_closure(set(nfa.initial_states)) - for symbol in u: - next_states: set[Hashable] = set() - for state in current: - next_states.update(nfa.delta(state, symbol)) - current = nfa.epsilon_closure(next_states) + current = run_nfa(nfa, u) return NFA( input_alphabet=nfa.input_alphabet, initial_states=frozenset(current), @@ -235,9 +230,3 @@ def state_residual_languages(dfa: DFA) -> dict[Hashable, DFA]: sub.accepting_states = dfa.accepting_states residuals[state] = minimize(_to_dfa(sub)) return residuals - - -def _effective_alphabet(aut: LabeledAutomaton) -> frozenset[Any]: - from sofic.automata.algorithms import _effective_alphabet as _shared - - return _shared(aut) diff --git a/sofic/automata/languages/quotients.py b/sofic/automata/languages/quotients.py index 070cb16..ad32484 100644 --- a/sofic/automata/languages/quotients.py +++ b/sofic/automata/languages/quotients.py @@ -5,11 +5,11 @@ from collections.abc import Sequence from typing import Any +from sofic.automata._words import _words_up_to from sofic.automata.languages._quotient_utils import ( _prefixes_if_suffix, _residual_from_state, _suffixes_if_prefix, - _words_up_to, ) from sofic.automata.languages.automaton_ops import ( left_quotient_automaton, @@ -54,9 +54,8 @@ def left_quotients(language: RegularLanguage) -> frozenset[RegularLanguage]: if isinstance(lang, ExplicitLanguage): alphabet = lang.alphabet quotients: set[RegularLanguage] = set() - for length in range(4): - for word in _words_up_to(length, alphabet): - quotients.add(left_quotient(word, lang)) + for word in _words_up_to(3, alphabet): + quotients.add(left_quotient(word, lang)) return frozenset(quotients) raise TypeError(f"unsupported language type {type(lang)!r}") diff --git a/sofic/automata/transducer_simulation.py b/sofic/automata/transducer_simulation.py index cd700d5..b6fd9c8 100644 --- a/sofic/automata/transducer_simulation.py +++ b/sofic/automata/transducer_simulation.py @@ -50,12 +50,7 @@ def step(config: tuple[Hashable, tuple[Any, ...]], symbol: Any) -> Iterator[tupl for transition in transducer.graph.out_transitions(state): if transition.data.get(ATTR_SYMBOL) != symbol: continue - if moore: - extended = output_prefix + _state_output(transducer, transition.target) - else: - out = transition.data.get(ATTR_OUTPUT) - extended = output_prefix + ((out,) if out is not None and out is not EPSILON else ()) - yield (transition.target, extended) + yield (transition.target, output_prefix + _emitted(transducer, transition, moore=moore)) stack = simulate_configs( initial, @@ -85,12 +80,7 @@ def _epsilon_closure( for transition in transducer.graph.out_transitions(state): if transition.data.get(ATTR_SYMBOL, EPSILON) is not EPSILON: continue - if moore: - extended = output_prefix + _state_output(transducer, transition.target) - else: - out = transition.data.get(ATTR_OUTPUT) - extended = output_prefix + ((out,) if out is not None and out is not EPSILON else ()) - config = (transition.target, extended) + config = (transition.target, output_prefix + _emitted(transducer, transition, moore=moore)) if config not in closure: closure.add(config) stack.append(config) @@ -131,20 +121,14 @@ def _has_productive_epsilon_cycle( if target not in component: continue for edge_data in eps_graph.get_edge_data(source, target).values(): - transition = edge_data["transition"] - if _epsilon_edge_is_productive(transducer, transition.target, transition.data, moore=moore): + if _emitted(transducer, edge_data["transition"], moore=moore): return True return False -def _epsilon_edge_is_productive( - transducer: Transducer, - target: Hashable, - data: dict[str, Any], - *, - moore: bool, -) -> bool: +def _emitted(transducer: Transducer, transition: Any, *, moore: bool) -> tuple[Any, ...]: + """Output appended when taking ``transition`` (Moore: the entered state's output).""" if moore: - return bool(_state_output(transducer, target)) - output = data.get(ATTR_OUTPUT) - return output is not None and output is not EPSILON + return _state_output(transducer, transition.target) + out = transition.data.get(ATTR_OUTPUT) + return (out,) if out is not None and out is not EPSILON else () diff --git a/sofic/automata/vpa_constructions.py b/sofic/automata/vpa_constructions.py index a49c733..494a39d 100644 --- a/sofic/automata/vpa_constructions.py +++ b/sofic/automata/vpa_constructions.py @@ -99,12 +99,13 @@ def normalize(vpa: VisiblyPushdownAutomaton) -> NormalVPA: A wildcard return (no stack symbol) fires on every stack symbol, and also on the empty stack when ``vpa`` has a bottom symbol; a return guarded by the - bottom symbol fires only on the empty stack. + bottom symbol fires only on the empty stack. A call without a stack symbol + never fires and is dropped. """ bottom = vpa.bottom_stack_symbol gamma = {g for g in vpa.stack_alphabet if g != bottom} for transition in vpa.transitions(): - if transition.data.get(ATTR_KIND) == KIND_CALL: + if transition.data.get(ATTR_KIND) == KIND_CALL and transition.data.get(ATTR_STACK_SYMBOL) is not None: gamma.add(transition.data.get(ATTR_STACK_SYMBOL)) result = NormalVPA( call_alphabet=frozenset(vpa.call_alphabet), @@ -121,7 +122,8 @@ def normalize(vpa: VisiblyPushdownAutomaton) -> NormalVPA: continue source, target = transition.source, transition.target if kind == KIND_CALL: - result.calls.add((source, symbol, target, data.get(ATTR_STACK_SYMBOL))) + if data.get(ATTR_STACK_SYMBOL) is not None: + result.calls.add((source, symbol, target, data.get(ATTR_STACK_SYMBOL))) elif kind == KIND_INTERNAL: result.internals.add((source, symbol, target)) elif kind == KIND_RETURN: diff --git a/sofic/automata/vpa_simulation.py b/sofic/automata/vpa_simulation.py index 55efe6c..f40e1cb 100644 --- a/sofic/automata/vpa_simulation.py +++ b/sofic/automata/vpa_simulation.py @@ -7,47 +7,34 @@ from sofic.automata._config_simulation import simulate_configs from sofic.automata.vpa import VisiblyPushdownAutomaton -from sofic.graph import ATTR_KIND, ATTR_STACK_SYMBOL, ATTR_SYMBOL, KIND_CALL, KIND_INTERNAL, KIND_RETURN - -_BOTTOM = object() +from sofic.automata.vpa_constructions import BOTTOM, normalize Config = tuple[Hashable, tuple[Any, ...]] def recognizes_vpa(vpa: VisiblyPushdownAutomaton, word: Sequence[Any]) -> bool: - """Return whether ``vpa`` accepts ``word``.""" - if vpa.initial_state is None: - return False + """Return whether ``vpa`` accepts ``word``. - bottom = _BOTTOM if vpa.bottom_stack_symbol is None else vpa.bottom_stack_symbol - initial: set[Config] = {(vpa.initial_state, (bottom,))} + Simulates the normalized form (:func:`~sofic.automata.vpa_constructions.normalize`): + the stack holds only pushed symbols, its emptiness stands for the explicit + bottom, and a return guarded by ``BOTTOM`` fires only on the empty stack. + """ + machine = normalize(vpa) + calls, internals, returns = machine.call_map(), machine.internal_map(), machine.return_map() def step(config: Config, symbol: Any) -> Iterator[Config]: state, stack = config - for transition in vpa.graph.out_transitions(state): - data = transition.data - if data.get(ATTR_SYMBOL) != symbol: - continue - kind = data.get(ATTR_KIND) - if kind == KIND_CALL: - stack_sym = data.get(ATTR_STACK_SYMBOL) - if stack_sym is None: - continue - yield (transition.target, stack + (stack_sym,)) - elif kind == KIND_RETURN: - stack_sym = data.get(ATTR_STACK_SYMBOL) - if len(stack) <= 1: - if vpa.bottom_stack_symbol is None: - continue - if stack_sym is not None and stack_sym != stack[-1]: - continue - yield (transition.target, stack) - continue - if stack_sym is not None and stack_sym != stack[-1]: - continue - yield (transition.target, stack[:-1]) - elif kind == KIND_INTERNAL: - yield (transition.target, stack) - + for target, pushed in calls.get((state, symbol), ()): + yield (target, (*stack, pushed)) + for target in internals.get((state, symbol), ()): + yield (target, stack) + for guard, target in returns.get((state, symbol), ()): + if guard == BOTTOM: + if not stack: + yield (target, stack) + elif stack and stack[-1] == guard: + yield (target, stack[:-1]) + + initial: set[Config] = {(state, ()) for state in machine.initial} current = simulate_configs(initial, word, step) - return any(state in vpa.accepting_states for state, _ in current) + return any(state in machine.accepting for state, _ in current) diff --git a/sofic/examples/_construction.py b/sofic/examples/_construction.py new file mode 100644 index 0000000..09f9cb6 --- /dev/null +++ b/sofic/examples/_construction.py @@ -0,0 +1,131 @@ +"""Private helpers shared by the example constructor modules.""" + +from __future__ import annotations + +from collections.abc import Hashable, Iterable, Mapping, Sequence +from typing import Any + +import numpy as np + +from sofic.generators.epsilon_machine import EpsilonMachine +from sofic.generators.mealy import MealyHMM +from sofic.graph import ATTR_EMISSION, ATTR_PROB + + +def _uniform_initial(states: Sequence[Hashable]) -> dict[Hashable, float]: + if not states: + return {} + mass = 1.0 / len(states) + return dict.fromkeys(states, mass) + + +def _stationary_initial(states: Sequence[Hashable], transition: Any) -> dict[Hashable, Any]: + """Stationary distribution of ``transition`` keyed by ``states``. + + Float matrices are row-normalized first and fall back to the uniform + distribution when a row is empty or the solve fails; exact (sympy) matrices + are solved as given and keep exact entries. + """ + from sofic.generators.prob import as_prob, has_symbolic + from sofic.generators.stationary import stationary_distribution_from_transition + + if not states: + return {} + matrix = np.asarray(transition) + if not (matrix.dtype == object and has_symbolic(matrix.ravel())): + matrix = np.asarray(matrix, dtype=float) + row_sums = matrix.sum(axis=1) + if np.any(row_sums <= 0.0): + return _uniform_initial(states) + matrix = matrix / row_sums[:, None] + try: + pi = stationary_distribution_from_transition(matrix) + except (ValueError, np.linalg.LinAlgError): + return _uniform_initial(states) + if pi.dtype == object or has_symbolic(pi.ravel()): + return {state: as_prob(pi[i]) for i, state in enumerate(states)} + return {state: float(pi[i]) for i, state in enumerate(states)} + + +def _normalize_edges( + edges: Sequence[tuple[Hashable, Hashable, Any, float]], +) -> list[tuple[Hashable, Hashable, Any, float]]: + row_totals: dict[Hashable, float] = {} + for source, _target, _symbol, prob in edges: + row_totals[source] = row_totals.get(source, 0.0) + float(prob) + normalized = [] + for source, target, symbol, prob in edges: + total = row_totals[source] + normalized.append((source, target, symbol, float(prob) / total if total else 0.0)) + return normalized + + +def _edge_machine( + edges: Iterable[tuple[Hashable, Hashable, Any, float]], + *, + machine_type: type[MealyHMM] = EpsilonMachine, + name: str | None = None, + initial_distribution: Mapping[Hashable, float] | None = None, + normalize: bool = True, + validate: bool = True, +) -> MealyHMM: + edge_list = list(edges) + if normalize: + edge_list = _normalize_edges(edge_list) + + states = list(dict.fromkeys([source for source, *_ in edge_list] + [target for _source, target, *_ in edge_list])) + symbols = frozenset(symbol for _source, _target, symbol, _prob in edge_list) + if initial_distribution is not None: + initial = dict(initial_distribution) + else: + index = {state: i for i, state in enumerate(states)} + transition = np.zeros((len(states), len(states)), dtype=float) + for source, target, _symbol, prob in edge_list: + transition[index[source], index[target]] += float(prob) + initial = _stationary_initial(states, transition) + + machine = machine_type(initial_distribution=initial, observation_alphabet=symbols) + if name is not None: + machine.name = name + for state in states: + machine.graph.add_state(state) + for source, target, symbol, prob in edge_list: + if prob > 0.0: + machine.graph.add_transition(source, target, **{ATTR_EMISSION: symbol, ATTR_PROB: float(prob)}) + if validate: + machine.validate() + return machine + + +def _relabel( + machine: MealyHMM, + *, + symbols: Mapping[Any, Any] | None = None, + states: Mapping[Hashable, Hashable] | None = None, + machine_type: type[MealyHMM] | None = None, + name: str | None = None, +) -> MealyHMM: + """Copy ``machine`` with emissions and states renamed; probabilities are kept as-is. + + A ``states`` mapping that covers every state also fixes the output state + order (the order of its values). + """ + symbol_map = dict(symbols or {}) + state_map = dict(states or {}) + old_states = list(machine.graph.states()) + if set(state_map) >= set(old_states): + old_states = [old for old in state_map if old in set(old_states)] + out = (machine_type or type(machine))( + initial_distribution={state_map.get(s, s): p for s, p in machine.initial_distribution.items()}, + observation_alphabet=frozenset(symbol_map.get(x, x) for x in machine.observation_alphabet), + ) + if name is not None: + out.name = name + for state in old_states: + out.graph.add_state(state_map.get(state, state), **dict(machine.graph.state_attrs(state))) + for t in machine.transitions(): + data = dict(t.data) + data[ATTR_EMISSION] = symbol_map.get(data[ATTR_EMISSION], data[ATTR_EMISSION]) + out.graph.add_transition(state_map.get(t.source, t.source), state_map.get(t.target, t.target), **data) + out.validate() + return out diff --git a/sofic/examples/epsilon_machines.py b/sofic/examples/epsilon_machines.py index 6673ea6..239039d 100644 --- a/sofic/examples/epsilon_machines.py +++ b/sofic/examples/epsilon_machines.py @@ -27,6 +27,7 @@ import numpy as np +from sofic.examples._construction import _edge_machine, _relabel, _stationary_initial from sofic.exceptions import SoficError from sofic.generators.epsilon_machine import EpsilonMachine from sofic.generators.mealy import MealyHMM @@ -35,20 +36,6 @@ from sofic.states import sequential_labels -def _stationary_distribution( - states: Sequence[Hashable], - symbol_matrices: Mapping[Any, np.ndarray], -) -> dict[Hashable, Any]: - from sofic.generators.prob import as_prob, has_symbolic - from sofic.generators.stationary import stationary_distribution_from_transition - - transition = sum(symbol_matrices.values()) - pi = stationary_distribution_from_transition(transition) - if pi.dtype == object or has_symbolic(pi.ravel()): - return {states[i]: as_prob(pi[i]) for i in range(len(states))} - return {states[i]: float(pi[i]) for i in range(len(states))} - - def from_symbol_matrices( states: Sequence[Hashable], symbols: Sequence[Any], @@ -76,7 +63,9 @@ def from_symbol_matrices( arrays[symbol] = arr pi = ( - dict(initial_distribution) if initial_distribution is not None else _stationary_distribution(state_list, arrays) + dict(initial_distribution) + if initial_distribution is not None + else _stationary_initial(state_list, sum(arrays.values())) ) eps = EpsilonMachine( initial_distribution=pi, @@ -106,8 +95,6 @@ def bernoulli(p: float = 0.5, *, symbols: tuple[Any, Any] = ("0", "1")) -> Epsil """Memoryless (Bernoulli) source with ``P(symbols[0]) = 1 - p``.""" if not 0.0 < p < 1.0: raise ValueError("p must be in (0, 1)") - from sofic.examples.processes import _edge_machine - zero, one = symbols state = sequential_labels(1)[0] return _edge_machine( @@ -150,8 +137,6 @@ def noisy_random_phase_slip() -> EpsilonMachine: emission noise at state ``D``. Prototype for block-convergence figures in *Anatomy of a Bit* :cite:`James2011`. """ - from sofic.examples.processes import _edge_machine - states = sequential_labels(5) a, b, c, d, e = states return _edge_machine( @@ -179,7 +164,9 @@ def golden_mean(p: float = 0.5) -> EpsilonMachine: see :func:`golden_mean_markov`; for the bidirectional machine in Ellison et al., arXiv:0905.3587, Fig.~4, see :func:`golden_mean_forward` and :func:`golden_mean_reverse`; for the Parry max-entropy measure on the - same shift, see :func:`golden_mean_shift_parry`. + same shift, see :func:`golden_mean_shift_parry`. cmpy's + :func:`~sofic.examples.processes.GoldenMean` is the ``0 <-> 1`` mirror + (forbids ``00``) and equals :func:`golden_mean_forward`. """ if not 0.0 < p < 1.0: raise ValueError("p must be in (0, 1)") @@ -244,16 +231,7 @@ def golden_mean_markov(p: float = 0.5) -> EpsilonMachine: For the paper's ``A``/``B`` labeling and bidirectional machine, prefer :func:`golden_mean_forward`. """ - if not 0.0 < p < 1.0: - raise ValueError("p must be in (0, 1)") - return from_symbol_matrices( - sequential_labels(2), - (0, 1), - { - 0: np.array([[0.0, 0.0], [p, 0.0]]), - 1: np.array([[0.0, 1.0], [0.0, 1.0 - p]]), - }, - ) + return _relabel(golden_mean_forward(1 - p), states={"B": "A", "A": "B"}) def golden_mean_shift_parry() -> EpsilonMachine: @@ -295,16 +273,10 @@ def restricted_golden_mean(k: int = 1) -> EpsilonMachine: """ if k < 1: raise ValueError("k must be >= 1") - states = sequential_labels(k + 1) - n = k + 1 - t0 = np.zeros((n, n), dtype=float) - t1 = np.zeros((n, n), dtype=float) - t0[0, 1] = 0.5 - t1[0, 0] = 0.5 - for i in range(1, k): - t1[i, i + 1] = 1.0 - t1[k, 0] = 1.0 - return from_symbol_matrices(states, (0, 1), {0: t0, 1: t1}) + from sofic.examples.processes import RestrictedGM + + states = {str(i): label for i, label in enumerate(sequential_labels(k + 1))} + return _relabel(RestrictedGM(k), symbols={"0": 0, "1": 1}, states=states) def nemo_process(p: float = 0.5, q: float = 0.5) -> EpsilonMachine: @@ -340,7 +312,8 @@ def nemo_process(p: float = 0.5, q: float = 0.5) -> EpsilonMachine: def phase_slip_backtrack(p: float = 0.5, q: float = 0.5) -> EpsilonMachine: """Phase-Slip Backtrack (PSB) Process (``R=3``, ``k_chi=2``). - James, Mahoney, Ellison & Crutchfield, arXiv:1010.5545, Fig.~2. + James, Mahoney, Ellison & Crutchfield, arXiv:1010.5545, Fig.~2. A different + process from cmpy's :func:`~sofic.examples.processes.PSB`. """ if not 0.0 < p < 1.0 or not 0.0 < q < 1.0: raise ValueError("p and q must be in (0, 1)") @@ -374,7 +347,8 @@ def butterfly_process() -> EpsilonMachine: Mahoney et al., arXiv:0906.5099, Fig.~1. Each causal state emits every symbol with probability ``1/8``; synchronizing symbols ``2``--``7`` always - reach the same causal state regardless of the source. + reach the same causal state regardless of the source. Not cmpy's + :func:`~sofic.examples.processes.Butterfly` (two branches per state, ``h_mu = 1``). """ states = ("A", "B", "C", "D", "E") prob = 1.0 / 8.0 @@ -388,8 +362,6 @@ def butterfly_process() -> EpsilonMachine: 6: "C", 7: "E", } - from sofic.examples.processes import _edge_machine - edges = [] for source in states: for symbol in range(8): @@ -405,18 +377,9 @@ def butterfly_process() -> EpsilonMachine: def ellison_fig9_forward() -> EpsilonMachine: """Forward ε-machine from Ellison et al., arXiv:1107.2168, Fig.~9.""" - from sofic.examples.processes import _edge_machine + from sofic.examples.processes import IrreversibleTwoState - return _edge_machine( - [ - ("A", "A", 0, 0.5), - ("A", "B", 1, 0.5), - ("B", "B", 1, 0.5), - ("B", "A", 2, 0.5), - ], - initial_distribution={"A": 0.5, "B": 0.5}, - normalize=False, - ) + return _relabel(IrreversibleTwoState(0.5, 0.5), symbols={"0": 0, "1": 1, "2": 2}) def tent_map_misiurewicz_a(symbolic: bool = False): @@ -879,7 +842,6 @@ def _stationary_distribution_from_joint_graph( graph: TransitionGraph, ) -> dict[tuple[str, str], Any]: from sofic.generators.prob import as_prob, has_symbolic, zeros - from sofic.generators.stationary import stationary_distribution_from_transition states = list(graph.states()) if not states: @@ -894,13 +856,7 @@ def _stationary_distribution_from_joint_graph( transition[source, target] = as_prob(transition[source, target]) + as_prob( transition_edge.data.get(ATTR_PROB, 0.0) ) - if symbolic: - pi = stationary_distribution_from_transition(transition) - return {states[position]: as_prob(pi[position]) for position in range(len(states))} - stationary = np.ones(len(states), dtype=float) / len(states) - for _ in range(20_000): - stationary = stationary @ np.asarray(transition, dtype=float) - return {states[position]: float(stationary[position]) for position in range(len(states))} + return _stationary_initial(states, transition) def _project_bidirectional_side( diff --git a/sofic/examples/processes.py b/sofic/examples/processes.py index b0b11e3..ee4aa43 100644 --- a/sofic/examples/processes.py +++ b/sofic/examples/processes.py @@ -15,6 +15,15 @@ import numpy as np from sofic.automata.transducers import MealyMachine +from sofic.examples._construction import _edge_machine, _relabel +from sofic.examples.epsilon_machines import ( + alternating_biased_coins, + bernoulli, + even_process, + golden_mean_forward, + nemo_process, + noisy_random_phase_slip, +) from sofic.generators.base import QuasiStochasticModel from sofic.generators.epsilon_machine import EpsilonMachine from sofic.generators.mealy import MealyHMM @@ -22,6 +31,8 @@ RecurrentEpsilonMachine = EpsilonMachine +_STR_BITS = {0: "0", 1: "1"} + def _require_machine_type(machine_type: Any, *allowed: type) -> None: """Raise ``NotImplementedError`` unless ``machine_type`` is ``None`` or allowed. @@ -39,79 +50,6 @@ def _as_alphabet(symbols: int | Sequence[Any]) -> tuple[Any, ...]: return tuple(symbols) -def _uniform_initial(states: Sequence[Hashable]) -> dict[Hashable, float]: - if not states: - return {} - mass = 1.0 / len(states) - return dict.fromkeys(states, mass) - - -def _stationary_initial( - states: Sequence[Hashable], - edges: Sequence[tuple[Hashable, Hashable, Any, float]], -) -> dict[Hashable, float]: - from sofic.generators.stationary import stationary_distribution_from_transition - - if not states: - return {} - index = {state: i for i, state in enumerate(states)} - transition = np.zeros((len(states), len(states)), dtype=float) - for source, target, _symbol, prob in edges: - transition[index[source], index[target]] += float(prob) - row_sums = transition.sum(axis=1) - if np.any(row_sums <= 0.0): - return _uniform_initial(states) - transition = transition / row_sums[:, None] - try: - pi = stationary_distribution_from_transition(transition) - except (ValueError, np.linalg.LinAlgError): - return _uniform_initial(states) - return {state: float(pi[i]) for i, state in enumerate(states)} - - -def _normalize_edges( - edges: Sequence[tuple[Hashable, Hashable, Any, float]], -) -> list[tuple[Hashable, Hashable, Any, float]]: - row_totals: dict[Hashable, float] = {} - for source, _target, _symbol, prob in edges: - row_totals[source] = row_totals.get(source, 0.0) + float(prob) - normalized = [] - for source, target, symbol, prob in edges: - total = row_totals[source] - normalized.append((source, target, symbol, float(prob) / total if total else 0.0)) - return normalized - - -def _edge_machine( - edges: Iterable[tuple[Hashable, Hashable, Any, float]], - *, - machine_type: type[MealyHMM] = EpsilonMachine, - name: str | None = None, - initial_distribution: Mapping[Hashable, float] | None = None, - normalize: bool = True, - validate: bool = True, -) -> MealyHMM: - edge_list = list(edges) - if normalize: - edge_list = _normalize_edges(edge_list) - - states = list(dict.fromkeys([source for source, *_ in edge_list] + [target for _source, target, *_ in edge_list])) - symbols = frozenset(symbol for _source, _target, symbol, _prob in edge_list) - initial = dict(initial_distribution) if initial_distribution is not None else _stationary_initial(states, edge_list) - - machine = machine_type(initial_distribution=initial, observation_alphabet=symbols) - if name is not None: - machine.name = name - for state in states: - machine.graph.add_state(state) - for source, target, symbol, prob in edge_list: - if prob > 0.0: - machine.graph.add_transition(source, target, **{ATTR_EMISSION: symbol, ATTR_PROB: float(prob)}) - if validate: - machine.validate() - return machine - - def _from_string( spec: str, *, @@ -173,7 +111,7 @@ def _compatible_machine_type(machine_type: Any, default: type[MealyHMM] = Epsilo def ABC(p: float = 0.75, q: float = 0.25) -> EpsilonMachine: if math.isclose(p, q): return _from_string(f"A A 0 {p}; A A 1 {1 - p}", name="ABC Process") - return _from_string(f"A B 0 {p}; A B 1 {1 - p}; B A 0 {q}; B A 1 {1 - q}", name="ABC Process") + return _relabel(alternating_biased_coins(1 - p, 1 - q), symbols=_STR_BITS, name="ABC Process") def AFC(n: int) -> EpsilonMachine: @@ -288,13 +226,7 @@ def BeforeAfter(machine_type: Any = MealyHMM, style: str = "simple") -> MealyHMM def BiasedCoin(bias: float, machine_type: Any = EpsilonMachine) -> EpsilonMachine: - cls = _compatible_machine_type(machine_type) - return _edge_machine( - [("A", "A", "1", bias), ("A", "A", "0", 1 - bias)], - machine_type=cls, - name=f"Coin, p = {bias}", - normalize=False, - ) + return _relabel(bernoulli(bias), machine_type=_compatible_machine_type(machine_type), name=f"Coin, p = {bias}") def FairCoin(machine_type: Any = EpsilonMachine) -> EpsilonMachine: @@ -436,6 +368,7 @@ def BMC_lohr(p: float) -> MealyHMM: def Butterfly() -> EpsilonMachine: + """cmpy's two-branch Butterfly (``h_mu = 1``); not :func:`~sofic.examples.butterfly_process` (8 symbols, ``h_mu = 3``).""" return _from_string( """ B C 0 .5; D E 0 .5; C A 1 .5; E A 1 .5; A B 2 .5; @@ -556,8 +489,7 @@ def Ehrenfest(p: float = 0.5, N: int = 5, machine_type: Any = EpsilonMachine) -> def Even(machine_type: Any = EpsilonMachine, bias: float = 0.5) -> EpsilonMachine: if machine_type in (EpsilonMachine, RecurrentEpsilonMachine, None): - edges = [("A", "A", "0", bias), ("A", "B", "1", 1 - bias), ("B", "A", "1", 1)] - return _edge_machine(edges, machine_type=EpsilonMachine, name="Even Process", normalize=False) + return _relabel(even_process(bias), symbols=_STR_BITS, name="Even Process") raise NotImplementedError @@ -662,13 +594,9 @@ def Girvan_fig6d( def GoldenMean(bias: float = 0.5, machine_type: Any = EpsilonMachine) -> EpsilonMachine: + """Golden mean forbidding ``00``; :func:`~sofic.examples.golden_mean` is its ``0 <-> 1`` mirror (forbids ``11``).""" _require_machine_type(machine_type, EpsilonMachine, RecurrentEpsilonMachine) - return _edge_machine( - [("A", "A", "1", 1 - bias), ("A", "B", "0", bias), ("B", "A", "1", 1)], - machine_type=EpsilonMachine, - name="Golden Mean Process", - normalize=False, - ) + return _relabel(golden_mean_forward(1 - bias), symbols=_STR_BITS, name="Golden Mean Process") def RestrictedGM(k: int) -> EpsilonMachine: @@ -922,10 +850,7 @@ def MultipleN(n: int, machine_type: Any = MealyHMM, bias: float = 0.5) -> MealyH def Nemo(machine_type: Any = EpsilonMachine, p: float = 0.5, q: float = 0.5) -> EpsilonMachine: _require_machine_type(machine_type, EpsilonMachine, RecurrentEpsilonMachine) - return _from_string( - f"A A 1 {p}; A B 0 {1 - p}; B C 0 1; C A 0 {1 - q}; C A 1 {q};", - name="Nemo Process", - ) + return _relabel(nemo_process(p, q), symbols=_STR_BITS, name="Nemo Process") def NemoRedundant(machine_type: Any = MealyHMM, p: float = 0.5, q: float = 0.5) -> MealyHMM: @@ -950,8 +875,6 @@ def NoisyPeriod2(noise: float = 0.5) -> EpsilonMachine: def NRPS() -> EpsilonMachine: - from sofic.examples.epsilon_machines import noisy_random_phase_slip - return noisy_random_phase_slip() @@ -1077,6 +1000,7 @@ def PerturbedCoin(p: float = 0.2, q: float | None = None, machine_type: Any = Ep def PSB() -> EpsilonMachine: + """cmpy's Phase-Slip Backtrack; a different process from :func:`~sofic.examples.epsilon_machines.phase_slip_backtrack`.""" return _from_string("A B 1; A D 0; B B 1; B C 0; C D 0; D A 1", name="Phase-Slip Backtrack") diff --git a/sofic/examples/tetris.py b/sofic/examples/tetris.py index f83381e..c80a9d1 100644 --- a/sofic/examples/tetris.py +++ b/sofic/examples/tetris.py @@ -22,7 +22,7 @@ import numpy as np -from sofic.examples.processes import _edge_machine +from sofic.examples._construction import _edge_machine from sofic.generators.epsilon_machine import EpsilonMachine TETROMINOES: tuple[str, ...] = ("I", "J", "L", "O", "S", "T", "Z") diff --git a/sofic/generators/_morph_tests.py b/sofic/generators/_morph_tests.py new file mode 100644 index 0000000..518d0a4 --- /dev/null +++ b/sofic/generators/_morph_tests.py @@ -0,0 +1,154 @@ +"""Significance tests on two-row contingency tables of next-symbol counts. + +Shared by process CSSR (:mod:`sofic.generators.epsilon_inference`), stack CSSR, +and transCSSR (:mod:`sofic.generators.epsilon_transducer_inference`). Each table +holds the counts of the symbols following two sets of histories, one row per set. +""" + +from __future__ import annotations + +import zlib +from collections.abc import Mapping, Sequence +from functools import lru_cache +from typing import Any, Literal + +import numpy as np +from scipy import stats + +#: Contingency-table tests: G-test, Pearson chi-squared, or Monte Carlo exact G-test. +TableTest = Literal["g", "chi2", "exact"] + +#: Monte Carlo tables drawn per ``"exact"`` test. +EXACT_DRAWS = 999 + +#: Smallest expected count at which the ``"exact"`` test trusts the chi-squared limit. +EXACT_MIN_EXPECTED = 5.0 + + +def contingency_table(left: Mapping[Any, int], right: Mapping[Any, int], alphabet: Sequence[Any]) -> np.ndarray | None: + """Two-row table of ``left`` and ``right`` counts over the symbols either one observed. + + Returns ``None`` when there is nothing to test: no observed symbol, identical + rows, or a single symbol with an empty row or equal proportions. + """ + active = [symbol for symbol in alphabet if left.get(symbol, 0) + right.get(symbol, 0) > 0] + if not active: + return None + table = np.array( + [[left.get(symbol, 0) for symbol in active], [right.get(symbol, 0) for symbol in active]], + dtype=float, + ) + if np.allclose(table[0], table[1]): + return None + if table.shape[1] < 2: + left_total = table[0].sum() + right_total = table[1].sum() + if left_total == 0.0 or right_total == 0.0: + return None + if np.isclose(table[0, 0] / left_total, table[1, 0] / right_total): + return None + return table + + +def g_statistic(table: np.ndarray) -> float | None: + """G-test statistic of a contingency table, as ``scipy.stats.chi2_contingency`` computes it. + + Applies Yates' continuity correction when the table has one degree of freedom + (``chi2_contingency`` does so for every ``lambda_``, including + ``"log-likelihood"``), and returns ``None`` if an expected count is zero. + Inlined because the test runs once per pair of histories, and the general + scipy routine dominated inference time. + """ + expected = table.sum(axis=1, keepdims=True) * table.sum(axis=0, keepdims=True) / table.sum() + if np.any(expected == 0): + return None + observed = table + if (table.shape[0] - 1) * (table.shape[1] - 1) == 1: + diff = expected - observed + observed = observed + np.sign(diff) * np.minimum(0.5, np.abs(diff)) + with np.errstate(divide="ignore", invalid="ignore"): + terms = np.where(observed > 0, observed * np.log(observed / expected), 0.0) + return 2.0 * float(terms.sum()) + + +@lru_cache(maxsize=256) +def chi2_critical(alpha: float, dof: int) -> float: + """Upper ``alpha`` quantile of the chi-squared distribution with ``dof`` degrees of freedom.""" + return float(stats.chi2.ppf(1.0 - alpha, dof)) + + +def g_significant(table: np.ndarray, alpha: float) -> bool: + """Asymptotic G-test (see :func:`g_statistic`) at level ``alpha``.""" + statistic = g_statistic(table) + if statistic is None or not np.isfinite(statistic): + return False + return statistic > chi2_critical(alpha, max(1, table.shape[1] - 1)) + + +def chi2_significant(table: np.ndarray, alpha: float) -> bool: + """Pearson chi-squared test (with Yates' correction at one degree of freedom) at level ``alpha``.""" + try: + _statistic, p_value, _dof, expected = stats.chi2_contingency(table) + except ValueError: + return False + if np.any(expected == 0): + return False + return float(p_value) < alpha + + +def exact_g_pvalue(table: np.ndarray) -> float: + """Monte Carlo p-value of the G statistic among tables with the same margins. + + Tables are drawn uniformly given both margins (``scipy.stats.random_table``), + the exact null of equal morphs. The generator is seeded from the table itself, + so reconstruction stays deterministic. + """ + counts = table.astype(np.int64) + rows, cols = counts.sum(axis=1), counts.sum(axis=0) + expected = np.outer(rows, cols) / counts.sum() + + def g(observed: np.ndarray) -> np.ndarray: + with np.errstate(divide="ignore", invalid="ignore"): + terms = np.where(observed > 0, observed * np.log(observed / expected), 0.0) + return 2.0 * terms.sum(axis=(-2, -1)) + + rng = np.random.default_rng(zlib.crc32(counts.tobytes())) + draws = stats.random_table(rows, cols, seed=rng).rvs(size=EXACT_DRAWS) + observed = g(counts.astype(float)) + extreme = np.sum(g(draws.astype(float)) >= observed - 1e-9 * max(1.0, observed)) + return float((1 + extreme) / (1 + EXACT_DRAWS)) + + +def exact_significant(table: np.ndarray, alpha: float) -> bool: + """The ``"exact"`` decision for a contingency table of next-symbol counts. + + Uses the Monte Carlo exact G-test when an expected count is below + :data:`EXACT_MIN_EXPECTED` and the asymptotic G-test otherwise. + """ + if np.any(table.sum(axis=1) == 0): + return False + expected = np.outer(table.sum(axis=1), table.sum(axis=0)) / table.sum() + if expected.min() < EXACT_MIN_EXPECTED: + return exact_g_pvalue(table) < alpha + return g_significant(table, alpha) + + +def table_significant(table: np.ndarray, alpha: float, test: TableTest) -> bool: + """Whether the rows of ``table`` differ significantly at level ``alpha`` under ``test``.""" + if test == "exact": + return exact_significant(table, alpha) + if test == "g": + return g_significant(table, alpha) + return chi2_significant(table, alpha) + + +def table_score(table: np.ndarray, test: TableTest) -> float: + """Test statistic of ``table`` (lower is more similar): G for ``"g"``/``"exact"``, else Pearson.""" + if test == "chi2": + try: + statistic, _p, _dof, _expected = stats.chi2_contingency(table) + except ValueError: + return 0.0 + return float(statistic) + statistic = g_statistic(table) + return statistic if statistic is not None and np.isfinite(statistic) else 0.0 diff --git a/sofic/generators/_suffix_counts.py b/sofic/generators/_suffix_counts.py new file mode 100644 index 0000000..ed3779e --- /dev/null +++ b/sofic/generators/_suffix_counts.py @@ -0,0 +1,21 @@ +"""Sliding-window suffix scan shared by the CSSR history counters.""" + +from __future__ import annotations + +from collections.abc import Iterator, Sequence +from typing import Any + + +def infer_alphabet(tokens: Sequence[Any], alphabet: Sequence[Any] | None) -> tuple[Any, ...]: + """``alphabet`` as a tuple, or the distinct ``tokens`` sorted by ``repr``.""" + return tuple(sorted(set(tokens), key=repr)) if alphabet is None else tuple(alphabet) + + +def iter_suffixes(tokens: tuple[Any, ...], max_length: int) -> Iterator[tuple[int, tuple[Any, ...]]]: + """Yield ``(t, tokens[t - L : t])`` for each position ``t`` and each ``L`` in ``0..min(t, max_length)``. + + These are the pasts, up to ``max_length`` long, that precede the token at ``t``. + """ + for t in range(len(tokens)): + for length in range(min(t, max_length) + 1): + yield t, tokens[t - length : t] diff --git a/sofic/generators/base.py b/sofic/generators/base.py index 276ce6b..47325b9 100644 --- a/sofic/generators/base.py +++ b/sofic/generators/base.py @@ -79,19 +79,18 @@ def reverse(self) -> Self: """ from sofic.generators.reversal import is_markov_like, time_reverse_stochastic - if not is_markov_like(self): - from sofic.generators.epsilon_machine import EpsilonMachine - from sofic.generators.mealy import MealyHMM - from sofic.generators.moore import MooreHMM - - if isinstance(self, (MealyHMM, MooreHMM)): - from sofic.generators.epsilon_machine import EpsilonMachine - - if isinstance(self, EpsilonMachine): - return EpsilonMachine.from_time_reversed(self) - return EpsilonMachine.from_hmm(time_reverse_stochastic(self)) - raise NotImplementedError("time-reversed generators with edge emissions require EpsilonMachine.from_hmm") - return time_reverse_stochastic(self) + if is_markov_like(self): + return time_reverse_stochastic(self) + + from sofic.generators.epsilon_machine import EpsilonMachine + from sofic.generators.mealy import MealyHMM + from sofic.generators.moore import MooreHMM + + if isinstance(self, EpsilonMachine): + return EpsilonMachine.from_time_reversed(self) + if isinstance(self, (MealyHMM, MooreHMM)): + return EpsilonMachine.from_hmm(time_reverse_stochastic(self)) + raise NotImplementedError("time-reversed generators with edge emissions require EpsilonMachine.from_hmm") class HiddenMarkovModel(StochasticModel): diff --git a/sofic/generators/block_entropy.py b/sofic/generators/block_entropy.py index 318fee3..a6288f6 100644 --- a/sofic/generators/block_entropy.py +++ b/sofic/generators/block_entropy.py @@ -428,10 +428,10 @@ def _block_entropy_curves( def _stationary_symbol_matrices(machine: EpsilonMachine) -> tuple[np.ndarray, dict[Any, np.ndarray]]: - from sofic.generators.hmm_inference import _emission_transition_tensors + from sofic.generators.matrices import symbol_matrices pi = machine.stationary_distribution() - _, raw_matrices = _emission_transition_tensors(machine) + raw_matrices = symbol_matrices(machine.to_mealy()) n = len(pi) zero = np.zeros((n, n), dtype=float) matrices = { diff --git a/sofic/generators/channel_measures.py b/sofic/generators/channel_measures.py index 0ff3c58..28ac235 100644 --- a/sofic/generators/channel_measures.py +++ b/sofic/generators/channel_measures.py @@ -35,6 +35,8 @@ def channel_statistical_complexity(transducer: MealyMachine, input_process: Hidd ``S`` is the transducer's causal-state component of the driven joint process, weighted by its stationary occupation distribution. """ + from sofic.generators.stochastic import shannon_entropy + joint = driven_joint_generator(transducer, input_process) idx = joint.reindex() if len(idx) == 0: @@ -44,11 +46,7 @@ def channel_statistical_complexity(transducer: MealyMachine, input_process: Hidd for state, weight in zip(idx.states, pi, strict=True): transducer_state = state[1] if isinstance(state, tuple) and len(state) == 2 else state mass[transducer_state] += float(weight) - probs = np.array([value for value in mass.values() if value > 0.0], dtype=float) - if probs.size == 0: - return 0.0 - probs = probs / probs.sum() - return float(-(probs * np.log2(probs)).sum()) + return shannon_entropy(mass.values(), normalize=True) def driven_entropy_rate(transducer: MealyMachine, input_process: HiddenMarkovModel) -> float: diff --git a/sofic/generators/directional_flow.py b/sofic/generators/directional_flow.py index e5e9f6c..6aaf8ba 100644 --- a/sofic/generators/directional_flow.py +++ b/sofic/generators/directional_flow.py @@ -19,36 +19,22 @@ def _require_dit(): def _pair_block_distribution(generator: HiddenMarkovModel, *, history: int) -> Any: """Joint law over flattened ``(x0, y0, x1, y1, ...)`` windows.""" - from sofic.generators.hmm_inference import _stationary_emission_tensors + from sofic.generators.matrices import emission_tensors + from sofic.generators.words import _enumerate_words, _matrix_step dit = _require_dit() # Directional-flow statistics describe the stationary joint process, so weight # the initial state by the stationary law rather than ``initial_distribution``. - pi, joint = _stationary_emission_tensors(generator) + pi, joint = emission_tensors(generator, policy="stationary") + if any(not isinstance(symbol, tuple) or len(symbol) != 2 for symbol in joint): + raise TypeError("generator must emit length-2 tuple symbols") - block_length = history + 1 ones = np.ones(len(pi), dtype=float) outcomes: list[tuple[Any, ...]] = [] probs: list[float] = [] - - def walk( - mass: np.ndarray, - prefix: tuple[Any, ...], - steps_remaining: int, - ) -> None: - if steps_remaining == 0: - outcomes.append(prefix) - probs.append(float(mass @ ones)) - return - for symbol, matrix in joint.items(): - if not isinstance(symbol, tuple) or len(symbol) != 2: - raise TypeError("generator must emit length-2 tuple symbols") - next_mass = mass @ matrix - if next_mass.sum() <= 0.0: - continue - walk(next_mass, prefix + (symbol[0], symbol[1]), steps_remaining - 1) - - walk(pi.copy(), (), block_length) + for word, mass in _enumerate_words(list(joint), history + 1, pi.copy(), _matrix_step(joint, len(pi))): + outcomes.append(tuple(value for pair in word for value in pair)) + probs.append(float(mass @ ones)) total = sum(probs) if total > 0.0: probs = [p / total for p in probs] diff --git a/sofic/generators/epsilon_inference.py b/sofic/generators/epsilon_inference.py index ad3ab3c..294d868 100644 --- a/sofic/generators/epsilon_inference.py +++ b/sofic/generators/epsilon_inference.py @@ -8,17 +8,16 @@ from __future__ import annotations -import zlib from collections import Counter, defaultdict -from collections.abc import Callable, Iterable, Sequence +from collections.abc import Callable, Iterable, Mapping, Sequence from dataclasses import dataclass, field -from functools import lru_cache from typing import Any, ClassVar, Literal import numpy as np -from scipy import stats from sofic.exceptions import StochasticValidationError +from sofic.generators._morph_tests import contingency_table, table_score, table_significant +from sofic.generators._suffix_counts import infer_alphabet, iter_suffixes from sofic.generators.epsilon_machine import EpsilonMachine from sofic.graph import ATTR_EMISSION, ATTR_PROB, TransitionGraph @@ -50,19 +49,14 @@ def from_sequence( seq = tuple(sequence) if not seq: raise ValueError("sequence must be non-empty") - alphabet = tuple(sorted(set(seq), key=repr)) if alphabet is None else tuple(alphabet) + alphabet = infer_alphabet(seq, alphabet) unknown = set(seq) - set(alphabet) if unknown: raise ValueError(f"symbols {unknown!r} not in alphabet") - max_len = max_length if max_length is not None else len(seq) counts = cls(alphabet=alphabet) - n = len(seq) - for t in range(n): - for length in range(0, min(t, max_len) + 1): - history = seq[t - length : t] - counts.history_counts[history] += 1 - nxt = seq[t] - counts.next_counts[history][nxt] += 1 + for t, history in iter_suffixes(seq, max_length if max_length is not None else len(seq)): + counts.history_counts[history] += 1 + counts.next_counts[history][seq[t]] += 1 return counts def morph(self, history: History, *, smoothing: float = 0.0) -> dict[Any, float]: @@ -129,101 +123,11 @@ def _contingency_rows( left_histories: set[History], right_histories: set[History], ) -> np.ndarray | None: - left_obs = _observed_counts_for_morph(counts, left_histories) - right_obs = _observed_counts_for_morph(counts, right_histories) - active = [symbol for symbol in counts.alphabet if left_obs.get(symbol, 0) + right_obs.get(symbol, 0) > 0] - if not active: - return None - table = np.array( - [ - [left_obs.get(symbol, 0) for symbol in active], - [right_obs.get(symbol, 0) for symbol in active], - ], - dtype=float, + return contingency_table( + _observed_counts_for_morph(counts, left_histories), + _observed_counts_for_morph(counts, right_histories), + counts.alphabet, ) - if np.allclose(table[0], table[1]): - return None - if table.shape[1] < 2: - left_total = table[0].sum() - right_total = table[1].sum() - if left_total == 0.0 or right_total == 0.0: - return None - left_prob = table[0, 0] / left_total - right_prob = table[1, 0] / right_total - if np.isclose(left_prob, right_prob): - return None - return table - - -def _g_statistic(table: np.ndarray) -> float | None: - """G-test statistic of a contingency table, as ``scipy.stats.chi2_contingency`` computes it. - - Includes Yates' correction for one degree of freedom, and returns ``None`` if an - expected count is zero. Inlined because the test runs once per pair of - histories, and the general scipy routine dominated inference time. - """ - expected = table.sum(axis=1, keepdims=True) * table.sum(axis=0, keepdims=True) / table.sum() - if np.any(expected == 0): - return None - observed = table - if (table.shape[0] - 1) * (table.shape[1] - 1) == 1: - diff = expected - observed - observed = observed + np.sign(diff) * np.minimum(0.5, np.abs(diff)) - with np.errstate(divide="ignore", invalid="ignore"): - terms = np.where(observed > 0, observed * np.log(observed / expected), 0.0) - return 2.0 * float(terms.sum()) - - -@lru_cache(maxsize=256) -def _chi2_critical(alpha: float, dof: int) -> float: - return float(stats.chi2.ppf(1.0 - alpha, dof)) - - -#: Monte Carlo tables drawn per ``"exact"`` morph test. -_EXACT_DRAWS = 999 - -#: Smallest expected count at which the ``"exact"`` test trusts the chi-squared limit. -_EXACT_MIN_EXPECTED = 5.0 - - -def _exact_g_pvalue(table: np.ndarray) -> float: - """Monte Carlo p-value of the G statistic among tables with the same margins. - - Tables are drawn uniformly given both margins (``scipy.stats.random_table``), - the exact null of equal morphs. The generator is seeded from the table itself, - so reconstruction stays deterministic. - """ - counts = table.astype(np.int64) - rows, cols = counts.sum(axis=1), counts.sum(axis=0) - expected = np.outer(rows, cols) / counts.sum() - - def g(observed: np.ndarray) -> np.ndarray: - with np.errstate(divide="ignore", invalid="ignore"): - terms = np.where(observed > 0, observed * np.log(observed / expected), 0.0) - return 2.0 * terms.sum(axis=(-2, -1)) - - rng = np.random.default_rng(zlib.crc32(counts.tobytes())) - draws = stats.random_table(rows, cols, seed=rng).rvs(size=_EXACT_DRAWS) - observed = g(counts.astype(float)) - extreme = np.sum(g(draws.astype(float)) >= observed - 1e-9 * max(1.0, observed)) - return float((1 + extreme) / (1 + _EXACT_DRAWS)) - - -def _exact_significant(table: np.ndarray, alpha: float) -> bool: - """The ``"exact"`` decision for a contingency table of next-symbol counts. - - Uses the Monte Carlo exact G-test when an expected count is below - ``_EXACT_MIN_EXPECTED`` and the asymptotic G-test otherwise. - """ - if np.any(table.sum(axis=1) == 0): - return False - expected = np.outer(table.sum(axis=1), table.sum(axis=0)) / table.sum() - if expected.min() < _EXACT_MIN_EXPECTED: - return _exact_g_pvalue(table) < alpha - statistic = _g_statistic(table) - if statistic is None or not np.isfinite(statistic): - return False - return statistic > _chi2_critical(alpha, max(1, table.shape[1] - 1)) def _bonferroni_alpha( @@ -254,10 +158,14 @@ def morphs_differ( ) -> bool: """Return whether two history sets have significantly different morphs. - ``"g"`` and ``"chi2"`` use the chi-squared limit, which is unreliable when - expected counts are small. ``"exact"`` instead compares the G statistic with - tables drawn uniformly given the observed margins whenever an expected count is - below 5, and uses the G-test otherwise. + ``"g"`` is the G-test (log-likelihood ratio) with Yates' continuity correction + when the table has one degree of freedom, i.e. two observed symbols; this is + the statistic ``scipy.stats.chi2_contingency(table, lambda_="log-likelihood")`` + reports. transCSSR uses the same test. ``"g"`` and ``"chi2"`` use the + chi-squared limit, which is unreliable when expected counts are small. + ``"exact"`` instead compares the G statistic with tables drawn uniformly given + the observed margins whenever an expected count is below 5, and uses the + G-test otherwise. """ if test == "tv": left = counts.state_morph(left_histories) @@ -268,20 +176,7 @@ def morphs_differ( table = _contingency_rows(counts, left_histories, right_histories) if table is None: return False - if test == "exact": - return _exact_significant(table, alpha) - if test == "g": - statistic = _g_statistic(table) - if statistic is None or not np.isfinite(statistic): - return False - return statistic > _chi2_critical(alpha, max(1, table.shape[1] - 1)) - try: - statistic, p_value, _dof, expected = stats.chi2_contingency(table) - except ValueError: - return False - if np.any(expected == 0): - return False - return float(p_value) < alpha + return table_significant(table, alpha, test) def morph_test_score( @@ -299,231 +194,7 @@ def morph_test_score( table = _contingency_rows(counts, left_histories, right_histories) if table is None: return 0.0 - try: - if test in ("g", "exact"): - statistic = _g_statistic(table) - return statistic if statistic is not None and np.isfinite(statistic) else 0.0 - statistic, _p, _dof, _expected = stats.chi2_contingency(table) - return float(statistic) - except ValueError: - return 0.0 - - -def _default_lmax(n: int, alphabet_size: int, min_count: int) -> int: - if alphabet_size <= 0: - return 1 - return max(1, min(15, n // max(1, alphabet_size * min_count))) - - -def _grow_history(history: History, symbol: Any) -> History: - """Default successor: append the symbol without truncation (flat ε-machine CSSR).""" - return history + (symbol,) - - -def _cssr_homogenize( - counts: SuffixCounts, - *, - Lmax: int, - alpha: float, - test: MorphTest, - successor_fn: Callable[[History, Any], History] = _grow_history, -) -> tuple[dict[int, set[History]], dict[History, int]]: - """Return state id -> histories and history -> state id.""" - states: dict[int, set[History]] = {0: {counts.empty_history}} - history_to_state: dict[History, int] = {counts.empty_history: 0} - next_state_id = 1 - - for _length in range(Lmax + 1): - for state_id in sorted(states): - histories = set(states[state_id]) - for history in list(histories): - for symbol in counts.alphabet: - child = successor_fn(history, symbol) - if child in history_to_state: - continue - if counts.history_counts.get(child, 0) == 0: - continue - child_histories = {child} - if morphs_differ( - counts, - histories, - child_histories, - alpha=alpha, - test=test, - ): - best_state: int | None = None - best_score = float("inf") - for candidate_id, candidate_histories in states.items(): - if morphs_differ( - counts, - candidate_histories, - child_histories, - alpha=alpha, - test=test, - ): - continue - score = morph_test_score( - counts, - candidate_histories, - child_histories, - test=test, - ) - if score < best_score: - best_score = score - best_state = candidate_id - if best_state is None: - best_state = next_state_id - states[next_state_id] = set() - next_state_id += 1 - states[best_state].add(child) - history_to_state[child] = best_state - else: - states[state_id].add(child) - history_to_state[child] = state_id - return states, history_to_state - - -def _cssr_determinize( - states: dict[int, set[History]], - history_to_state: dict[History, int], - counts: SuffixCounts, - *, - length: int, - successor_fn: Callable[[History, Any], History] = _grow_history, -) -> dict[int, set[History]]: - """Split homogeneous states until transitions are unifilar.""" - current = {state_id: set(histories) for state_id, histories in states.items()} - changed = True - next_state_id = max(current) + 1 if current else 0 - - while changed: - changed = False - for state_id in sorted(current): - histories = current[state_id] - if len(histories) <= 1: - continue - for symbol in counts.alphabet: - buckets: dict[int, set[History]] = defaultdict(set) - for history in histories: - if counts.next_counts.get(history, Counter()).get(symbol, 0) == 0: - continue - child = successor_fn(history, symbol) - target = history_to_state.get(child) - if target is None: - continue - buckets[target].add(history) - if len(buckets) <= 1: - continue - # Keep the largest bucket in the original state; split others. - ordered = sorted(buckets.items(), key=lambda item: (-len(item[1]), min(item[1]))) - keep_target, keep_histories = ordered[0] - current[state_id] = keep_histories - for _target, split_histories in ordered[1:]: - new_id = next_state_id - next_state_id += 1 - current[new_id] = split_histories - for history in split_histories: - history_to_state[history] = new_id - changed = True - break - if changed: - break - return current - - -def _merge_similar_states( - states: dict[int, set[History]], - history_to_state: dict[History, int], - counts: SuffixCounts, - *, - alpha: float, - test: MorphTest, -) -> dict[int, set[History]]: - """Merge inferred states whose pooled morphs are statistically indistinguishable. - - Merging on the morph alone can fuse states with incompatible ``symbol -> - successor`` maps, yielding a non-unifilar partition. Callers must re-run - :func:`_cssr_determinize` afterwards to restore unifilarity. - """ - current = {state_id: set(histories) for state_id, histories in states.items()} - changed = True - while changed: - changed = False - state_ids = sorted(current) - for index, left_id in enumerate(state_ids): - if left_id not in current: - continue - for right_id in state_ids[index + 1 :]: - if right_id not in current: - continue - if morphs_differ( - counts, - current[left_id], - current[right_id], - alpha=alpha, - test=test, - ): - continue - current[left_id].update(current.pop(right_id)) - for history in current[left_id]: - history_to_state[history] = left_id - changed = True - break - if changed: - break - return current - - -def _drop_transient_states( - states: dict[int, set[History]], - history_to_state: dict[History, int], - counts: SuffixCounts, - *, - length: int, - successor_fn: Callable[[History, Any], History] = _grow_history, -) -> dict[int, set[History]]: - """Keep only states in bottom strongly connected components.""" - import networkx as nx - - successors: dict[int, dict[Any, set[int]]] = defaultdict(lambda: defaultdict(set)) - for state_id, histories in states.items(): - for history in histories: - for symbol in counts.alphabet: - if counts.next_counts.get(history, Counter()).get(symbol, 0) == 0: - continue - target = history_to_state.get(successor_fn(history, symbol)) - if target is None: - continue - successors[state_id][symbol].add(target) - - graph = nx.DiGraph() - for state_id in states: - graph.add_node(state_id) - for state_id, by_symbol in successors.items(): - for targets in by_symbol.values(): - for target in targets: - graph.add_edge(state_id, target) - - if graph.number_of_edges() == 0: - return states - - recurrent: set[int] = set() - for component in nx.strongly_connected_components(graph): - if not component: - continue - subgraph = graph.subgraph(component) - has_cycle = subgraph.number_of_edges() > 0 and ( - len(component) > 1 or any(subgraph.has_edge(node, node) for node in component) - ) - if not has_cycle: - continue - outgoing = any(graph.has_edge(v, w) for v in component for w in graph.nodes if w not in component) - if not outgoing: - recurrent.update(component) - - if not recurrent: - return states - return {state_id: histories for state_id, histories in states.items() if state_id in recurrent} + return table_score(table, test) def _cssr_default_lmax(n: int, alphabet_size: int) -> int: @@ -638,8 +309,12 @@ def _suffix_edges( return edges -def _recurrent_states(edges: dict[int, dict[Any, dict[int, set[History]]]]) -> list[set[int]]: - """Closed communicating classes (with at least one edge) of the state graph.""" +def _recurrent_states(edges: Mapping[int, Mapping[Any, Iterable[int]]]) -> list[set[int]]: + """Closed communicating classes (with at least one edge) of the state graph. + + ``edges[state][label]`` iterates over the successor states on ``label`` (a + symbol here, an ``(input, output)`` pair in transCSSR). + """ import networkx as nx graph = nx.DiGraph() diff --git a/sofic/generators/epsilon_transducer_inference.py b/sofic/generators/epsilon_transducer_inference.py index 2466b1b..db779fa 100644 --- a/sofic/generators/epsilon_transducer_inference.py +++ b/sofic/generators/epsilon_transducer_inference.py @@ -17,10 +17,10 @@ from dataclasses import dataclass, field from typing import Any, Literal -import numpy as np -from scipy import stats - from sofic.exceptions import StochasticValidationError +from sofic.generators._morph_tests import TableTest, contingency_table, table_score, table_significant +from sofic.generators._suffix_counts import infer_alphabet, iter_suffixes +from sofic.generators.epsilon_inference import _recurrent_states, suggest_lmax from sofic.generators.epsilon_transducer import EpsilonTransducer from sofic.graph import ATTR_OUTPUT, ATTR_PROB, ATTR_SYMBOL, TransitionGraph @@ -53,17 +53,12 @@ def from_sequences( raise ValueError("inputs and outputs must have equal length") if not xs: raise ValueError("sequences must be non-empty") - in_alpha = tuple(sorted(set(xs), key=repr)) if input_alphabet is None else tuple(input_alphabet) - out_alpha = tuple(sorted(set(ys), key=repr)) if output_alphabet is None else tuple(output_alphabet) - counts = cls(input_alphabet=in_alpha, output_alphabet=out_alpha) - pairs = tuple(zip(xs, ys, strict=True)) - n = len(pairs) - for t in range(n): - for length in range(0, min(t, max_length) + 1): - history = pairs[t - length : t] - counts.history_counts[history] += 1 - by_input = counts.next_counts.setdefault(history, {}) - by_input.setdefault(xs[t], Counter())[ys[t]] += 1 + counts = cls( + input_alphabet=infer_alphabet(xs, input_alphabet), output_alphabet=infer_alphabet(ys, output_alphabet) + ) + for t, history in iter_suffixes(tuple(zip(xs, ys, strict=True)), max_length): + counts.history_counts[history] += 1 + counts.next_counts.setdefault(history, {}).setdefault(xs[t], Counter())[ys[t]] += 1 return counts def output_counts(self, histories: set[JointHistory], input_symbol: Any) -> Counter[Any]: @@ -86,26 +81,6 @@ def state_morph(self, histories: set[JointHistory], input_symbol: Any) -> dict[A return {symbol: observed.get(symbol, 0) / total for symbol in self.output_alphabet} -def _output_contingency(left: Counter[Any], right: Counter[Any], alphabet: tuple[Any, ...]) -> np.ndarray | None: - active = [symbol for symbol in alphabet if left.get(symbol, 0) + right.get(symbol, 0) > 0] - if not active: - return None - table = np.array( - [[left.get(symbol, 0) for symbol in active], [right.get(symbol, 0) for symbol in active]], - dtype=float, - ) - if np.allclose(table[0], table[1]): - return None - if table.shape[1] < 2: - left_total = table[0].sum() - right_total = table[1].sum() - if left_total == 0.0 or right_total == 0.0: - return None - if np.isclose(table[0, 0] / left_total, table[1, 0] / right_total): - return None - return table - - #: Aggregated output counts of a state: ``agg[input_symbol]`` is a Counter over outputs. StateAggregate = dict[Any, Counter[Any]] @@ -126,18 +101,23 @@ def aggregates_differ( input_alphabet: tuple[Any, ...], output_alphabet: tuple[Any, ...], alpha: float, - test: Literal["g", "chi2", "exact"] = "g", + test: TableTest = "g", ) -> bool: - """Return whether two aggregated morphs differ on ``P(output | ., input)`` for some input.""" + """Return whether two aggregated morphs differ on ``P(output | ., input)`` for some input. + + Each input symbol's output counts are compared with the same test as process + CSSR (:func:`~sofic.generators.epsilon_inference.morphs_differ`); in particular + ``"g"`` is the G-test with Yates' continuity correction at one degree of freedom. + """ for input_symbol in input_alphabet: - table = _output_contingency( + table = contingency_table( left.get(input_symbol, Counter()), right.get(input_symbol, Counter()), output_alphabet, ) if table is None: continue - if _table_significant(table, alpha=alpha, test=test): + if table_significant(table, alpha, test): return True return False @@ -151,44 +131,16 @@ def _aggregate_score( ) -> float: total = 0.0 for input_symbol in input_alphabet: - table = _output_contingency( + table = contingency_table( left.get(input_symbol, Counter()), right.get(input_symbol, Counter()), output_alphabet, ) - if table is None: - continue - try: - with np.errstate(invalid="ignore", divide="ignore"): - statistic, _p, _dof, _expected = stats.chi2_contingency(table, lambda_="log-likelihood") - if np.isfinite(statistic): - total += float(statistic) - except ValueError: - continue + if table is not None: + total += table_score(table, "g") return total -def _table_significant(table: np.ndarray, *, alpha: float, test: Literal["g", "chi2", "exact"]) -> bool: - if test == "exact": - from sofic.generators.epsilon_inference import _exact_significant - - return _exact_significant(table, alpha) - try: - if test == "g": - with np.errstate(invalid="ignore", divide="ignore"): - statistic, _p, _dof, expected = stats.chi2_contingency(table, lambda_="log-likelihood") - if not np.isfinite(statistic) or np.any(expected == 0): - return False - dof = max(1, table.shape[1] - 1) - return float(statistic) > float(stats.chi2.ppf(1.0 - alpha, dof)) - statistic, p_value, _dof, expected = stats.chi2_contingency(table) - except ValueError: - return False - if np.any(expected == 0): - return False - return float(p_value) < alpha - - def _state_aggregate(counts: JointSuffixCounts, histories: Iterable[JointHistory]) -> StateAggregate: aggregate: StateAggregate = {} for history in histories: @@ -205,7 +157,7 @@ def _homogenize( *, Lmax: int, alpha: float, - test: Literal["g", "chi2", "exact"], + test: TableTest, min_count: int, ) -> list[set[JointHistory]]: """transCSSR homogenization: grow joint suffixes one ``(input, output)`` pair into the past. @@ -272,7 +224,7 @@ def _edges( *, Lmax: int, alpha: float, - test: Literal["g", "chi2", "exact"], + test: TableTest, ) -> dict[int, dict[tuple[Any, Any], dict[int, set[JointHistory]]]]: """Successor states by ``(input, output)`` pair, with the histories that lead there. @@ -330,22 +282,6 @@ def _edges( return edges -def _closed_classes(edges: dict[int, dict[tuple[Any, Any], dict[int, set[JointHistory]]]]) -> list[set[int]]: - import networkx as nx - - graph = nx.DiGraph() - graph.add_nodes_from(edges) - for source, by_pair in edges.items(): - for targets in by_pair.values(): - graph.add_edges_from((source, target) for target in targets) - condensed = nx.condensation(graph) - return [ - set(condensed.nodes[node]["members"]) - for node in condensed - if condensed.out_degree(node) == 0 and graph.subgraph(condensed.nodes[node]["members"]).number_of_edges() > 0 - ] - - def _determinize( states: list[set[JointHistory]], counts: JointSuffixCounts, @@ -353,7 +289,7 @@ def _determinize( *, Lmax: int, alpha: float, - test: Literal["g", "chi2", "exact"], + test: TableTest, ) -> tuple[list[set[JointHistory]], set[int]]: """Split alive states until each ``(input, output)`` pair has one alive successor.""" states = [set(h) for h in states] @@ -388,7 +324,7 @@ def _build_transducer( *, Lmax: int, alpha: float, - test: Literal["g", "chi2", "exact"], + test: TableTest, ) -> EpsilonTransducer: edges = _edges(states, counts, alive, Lmax=Lmax, alpha=alpha, test=test) history_to_state = {h: index for index in alive for h in states[index]} @@ -402,7 +338,7 @@ def _build_transducer( visits[state] += 1 break - classes = _closed_classes(edges) + classes = _recurrent_states(edges) if not classes: raise StochasticValidationError("no recurrent inferred states; the sample is too short for this Lmax") keep = max(classes, key=lambda members: (sum(visits[s] for s in members), -min(members))) @@ -468,7 +404,7 @@ def transcssr( output_alphabet: Sequence[Any] | None = None, Lmax: int | Literal["auto"] | None = None, alpha: float = 0.001, - test: Literal["g", "chi2", "exact"] = "g", + test: TableTest = "g", min_count: int = 5, correction: Literal["bonferroni"] | None = None, ) -> EpsilonTransducer: @@ -498,8 +434,6 @@ def transcssr( else len(tuple(input_alphabet)) * len(tuple(output_alphabet)) ) if Lmax == "auto": - from sofic.generators.epsilon_inference import suggest_lmax - max_length = suggest_lmax(list(zip(xs, ys, strict=True)), alpha=alpha) else: max_length = Lmax if Lmax is not None else _default_lmax(len(xs), joint_alphabet_size, min_count) @@ -523,6 +457,6 @@ def transcssr( states = _homogenize(counts, Lmax=max_length, alpha=alpha, test=test, min_count=min_count) everything = set(range(len(states))) edges = _edges(states, counts, everything, Lmax=max_length, alpha=alpha, test=test) - alive = set().union(*_closed_classes(edges)) or everything + alive = set().union(*_recurrent_states(edges)) or everything states, alive = _determinize(states, counts, alive, Lmax=max_length, alpha=alpha, test=test) return _build_transducer(states, counts, alive, xs, ys, Lmax=max_length, alpha=alpha, test=test) diff --git a/sofic/generators/hmm_inference.py b/sofic/generators/hmm_inference.py index 48f575b..8674f4e 100644 --- a/sofic/generators/hmm_inference.py +++ b/sofic/generators/hmm_inference.py @@ -16,7 +16,7 @@ import numpy as np from sofic.generators.base import HiddenMarkovModel -from sofic.graph import ATTR_EMISSION, ATTR_PROB +from sofic.generators.matrices import emission_tensors, symbol_matrices def _as_mealy_hmm(hmm: HiddenMarkovModel) -> Any: @@ -24,119 +24,6 @@ def _as_mealy_hmm(hmm: HiddenMarkovModel) -> Any: return hmm.to_mealy() -def _emission_transition_tensors_from_mealy( - hmm: Any, -) -> tuple[np.ndarray, dict[Any, np.ndarray]]: - """Return initial vector ``pi`` and symbol -> joint transition matrices. - - Symbols are keyed in a fixed (``repr``-sorted) order, so seeded sampling is - reproducible across interpreter runs. A model without an initial - distribution starts from its stationary distribution. - """ - from sofic.generators.prob import as_prob, has_symbolic, zeros - - idx = hmm.reindex() - n = len(idx) - edge_probs = [transition.data.get(ATTR_PROB, 0.0) for transition in hmm.transitions()] - init_probs = list(hmm.initial_distribution.values()) - symbolic = has_symbolic(edge_probs) or has_symbolic(init_probs) - - pi = zeros((n,), symbolic=symbolic) - if hmm.initial_distribution: - for state, mass in hmm.initial_distribution.items(): - pi[idx.index(state)] = as_prob(mass) - elif n: - pi = np.asarray(hmm.stationary_distribution(), dtype=object if symbolic else float) - - emissions = {transition.data.get(ATTR_EMISSION) for transition in hmm.transitions()} - {None} - symbols = sorted(set(hmm.observation_alphabet) | emissions, key=repr) - joint: dict[Any, np.ndarray] = {symbol: zeros((n, n), symbolic=symbolic) for symbol in symbols} - - for transition in hmm.transitions(): - emission = transition.data.get(ATTR_EMISSION) - if emission is None: - continue - i = idx.index(transition.source) - j = idx.index(transition.target) - joint[emission][i, j] = as_prob(joint[emission][i, j]) + as_prob(transition.data.get(ATTR_PROB, 0.0)) - return pi, joint - - -def _emission_transition_tensors( - hmm: HiddenMarkovModel, -) -> tuple[np.ndarray, dict[Any, np.ndarray]]: - """Return initial vector ``pi`` and symbol -> joint transition matrices.""" - return _emission_transition_tensors_from_mealy(_as_mealy_hmm(hmm)) - - -def _limit_distribution_from_initial(pi_initial: np.ndarray, transition: np.ndarray) -> np.ndarray | None: - """Return the limiting occupation law of ``pi_initial`` under ``transition``. - - On reducible chains the left-eigenvector stationary law is not unique; the - process measure is the limit reached from the model's initial distribution. - """ - pi = np.asarray(pi_initial, dtype=float).copy() - total = float(pi.sum()) - if total <= 0.0: - return None - pi /= total - matrix = np.asarray(transition, dtype=float) - n = len(pi) - for _ in range(max(100, 20 * n)): - nxt = pi @ matrix - mass = float(nxt.sum()) - if mass <= 0.0: - return None - nxt /= mass - if np.allclose(nxt, pi, rtol=1e-12, atol=1e-14): - pi = nxt - break - pi = nxt - pi[np.isclose(pi, 0.0, atol=1e-15)] = 0.0 - mass = float(pi.sum()) - if mass <= 0.0: - return None - return pi / mass - - -def _stationary_emission_tensors( - hmm: HiddenMarkovModel, -) -> tuple[np.ndarray, dict[Any, np.ndarray]]: - """Return the stationary state law and symbol -> joint transition matrices. - - Block/word statistics of a *stationary* process must weight the initial state - by the stationary distribution, not by the model's (possibly transient) - ``initial_distribution``. The stationary vector is recovered directly from the - summed emission-transition matrices so it stays aligned with ``joint``'s state - indexing. - - When the chain is reducible (multiple absorbing classes), the eigenvector - stationary law is not unique — prefer the limiting occupation reached from - ``initial_distribution``. Fall back to the eigenvector solution, then to the - initial vector, only when the limit cannot be formed. - """ - from sofic.generators.prob import zeros - from sofic.generators.stationary import stationary_distribution_from_transition - - pi_initial, joint = _emission_transition_tensors(hmm) - n = len(pi_initial) - if n == 0: - return pi_initial, joint - symbolic = pi_initial.dtype == object or any(matrix.dtype == object for matrix in joint.values()) - transition = zeros((n, n), symbolic=symbolic) - for matrix in joint.values(): - transition = transition + matrix - if not symbolic: - limited = _limit_distribution_from_initial(pi_initial, transition) - if limited is not None and np.allclose(limited @ transition, limited, rtol=1e-8, atol=1e-10): - return limited, joint - try: - pi = stationary_distribution_from_transition(transition) - except (np.linalg.LinAlgError, ValueError): - pi = pi_initial - return pi, joint - - def _forward_scaled( pi: np.ndarray, joint: dict[Any, np.ndarray], @@ -178,7 +65,7 @@ def forward(hmm: HiddenMarkovModel, observations: Sequence[Any], *, scaled: bool With ``scaled=True`` each row is normalized to sum to one (the numerically stable message used for posteriors); otherwise the raw messages are returned. """ - pi, joint = _emission_transition_tensors(hmm) + pi, joint = emission_tensors(hmm) obs = list(observations) if scaled: alpha_hat, _log_scales = _forward_scaled(pi, joint, obs) @@ -202,7 +89,7 @@ def backward(hmm: HiddenMarkovModel, observations: Sequence[Any], *, scaled: boo posterior is then ``normalize(alpha_hat[t] * beta_hat[t])`` (the per-row scaling constants cancel on renormalization). """ - _, joint = _emission_transition_tensors(hmm) + joint = symbol_matrices(_as_mealy_hmm(hmm)) n = next(iter(joint.values())).shape[0] if joint else len(_as_mealy_hmm(hmm).reindex()) obs = list(observations) beta = np.zeros((len(obs) + 1, n), dtype=float) @@ -245,7 +132,7 @@ def log_likelihood(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> float Uses the per-step-scaled forward recursion so the result stays finite for long sequences instead of underflowing to ``-inf``. """ - pi, joint = _emission_transition_tensors(hmm) + pi, joint = emission_tensors(hmm) _alpha_hat, log_scales = _forward_scaled(pi, joint, list(observations)) if not np.all(np.isfinite(log_scales)): return float("-inf") @@ -261,7 +148,7 @@ def smooth(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> np.ndarray: forward-backward smoother of Cappe, Moulines & Ryden (2005, Section 3.2). Rows for observation sequences of zero probability are returned as zeros. """ - pi, joint = _emission_transition_tensors(hmm) + pi, joint = emission_tensors(hmm) obs = list(observations) n_states = len(pi) alpha_hat, log_scales = _forward_scaled(pi, joint, obs) @@ -283,7 +170,7 @@ def two_slice_marginals(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> 2005, Section 3.2). Marginalizing over ``j`` recovers ``gamma[t]`` for ``t < n``. Returns an all-zero tensor for zero-probability sequences. """ - pi, joint = _emission_transition_tensors(hmm) + pi, joint = emission_tensors(hmm) obs = list(observations) n_states = len(pi) xi = np.zeros((len(obs), n_states, n_states), dtype=float) @@ -411,7 +298,7 @@ def baum_welch( alphabet = frozenset(mealy.observation_alphabet) seqs = _as_sequence_list(sequences) - pi, joint = _emission_transition_tensors_from_mealy(mealy) + pi, joint = emission_tensors(mealy) support = { (i, symbol, j) for symbol, matrix in joint.items() @@ -540,7 +427,7 @@ def score(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> dict[tuple[Has """ mealy = hmm.to_mealy() idx = mealy.reindex() - pi, joint = _emission_transition_tensors_from_mealy(mealy) + pi, joint = emission_tensors(mealy) edge_counts, _source_totals, _gamma0, loglik = _expected_edge_counts(pi, joint, list(observations)) if not np.isfinite(loglik): raise ValueError("observations have zero probability under the model; score is undefined") @@ -603,7 +490,7 @@ def observed_information(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> model has no free parameters. """ mealy = hmm.to_mealy() - pi, joint = _emission_transition_tensors_from_mealy(mealy) + pi, joint = emission_tensors(mealy) obs = list(observations) n_states = len(pi) @@ -684,7 +571,7 @@ def free_parameter_labels(hmm: HiddenMarkovModel) -> list[tuple[Hashable, Any, H """ mealy = hmm.to_mealy() idx = mealy.reindex() - _pi, joint = _emission_transition_tensors_from_mealy(mealy) + joint = symbol_matrices(mealy) free_edges, _reference, _source = _free_parameterization(joint, len(idx)) return [(idx.state(i), symbol, idx.state(j)) for i, symbol, j in free_edges] @@ -725,7 +612,7 @@ def _log_probabilities(values: np.ndarray) -> np.ndarray: def viterbi(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> list[Hashable]: mealy = _as_mealy_hmm(hmm) idx = mealy.reindex() - pi, joint = _emission_transition_tensors_from_mealy(mealy) + pi, joint = emission_tensors(mealy) n = len(idx) obs = list(observations) if n == 0: @@ -775,7 +662,7 @@ def sample( generator = rng if rng is not None else np.random.default_rng() mealy = _as_mealy_hmm(hmm) idx = mealy.reindex() - pi, joint = _emission_transition_tensors_from_mealy(mealy) + pi, joint = emission_tensors(mealy) total = float(np.sum(pi)) if not total > 0.0: raise ValueError("cannot sample: the initial state distribution has no mass") diff --git a/sofic/generators/matrices.py b/sofic/generators/matrices.py new file mode 100644 index 0000000..f1bca36 --- /dev/null +++ b/sofic/generators/matrices.py @@ -0,0 +1,252 @@ +"""Symbol-labeled transition matrices and start vectors for stochastic generators. + +For a Mealy-style hidden Markov model with states indexed by ``model.reindex()``, +the symbol-labeled joint transition matrices are + +.. math:: T^{(x)}_{ij} = P(S_{t+1} = j, X_t = x \\mid S_t = i), + +so the probability of a word ``w = x_0 ... x_{L-1}`` from a start vector +:math:`\\eta` is :math:`\\eta T^{(x_0)} \\cdots T^{(x_{L-1})} \\mathbf{1}` +(:cite:`Rabiner1989,Ellison2009`). Their sum +:math:`T = \\sum_x T^{(x)}` is the internal state-to-state transition matrix. + +Start policies +-------------- +Every helper that needs an initial state law takes ``policy``: + +``"model"`` + The model's ``initial_distribution``, or its ``stationary_distribution()`` + when no initial distribution is given. This is the law of the generator *as + specified* and is used for likelihoods, decoding, sampling, and finite-word + probabilities. + +``"stationary"`` + The stationary law of the process. Block and window statistics of a + stationary process must weight the first state by this law rather than by a + possibly transient ``initial_distribution``. On reducible chains the + eigenvector stationary law is not unique, so the limit of the ``"model"`` + start vector under :math:`T` is preferred; the eigenvector solution is used + when that limit cannot be formed (or the model is symbolic), and the + ``"model"`` start vector as a last resort. + +An explicit ``start`` (a state, a state-to-mass mapping, or a dense vector in +``reindex()`` order) overrides either policy. +""" + +from __future__ import annotations + +from collections.abc import Callable, Hashable, Iterable, Mapping, Sequence +from typing import Any, Literal + +import numpy as np + +from sofic.generators.prob import as_prob, has_symbolic, zeros +from sofic.graph import ATTR_EMISSION, ATTR_PROB + +StartPolicy = Literal["model", "stationary"] +StartSpec = Hashable | Mapping[Hashable, Any] | Sequence[Any] | np.ndarray | None + +SKIP = object() + + +def accumulate_matrices( + model: Any, + *, + attr: str = ATTR_PROB, + states: Iterable[Hashable] | None = None, + label: Callable[[Any], Any] = lambda transition: None, + labels: Iterable[Any] = (), + symbolic: bool | None = None, +) -> tuple[dict[Any, np.ndarray], list[Hashable]]: + """Accumulate edge ``attr`` weights into one state-to-state matrix per edge label. + + Rows and columns follow ``states`` when given (edges leaving the set are + ignored), otherwise all model states in iteration order. ``label(transition)`` + keys the matrix an edge contributes to; returning ``SKIP`` drops the edge. + Matrices for ``labels`` are allocated first, in the given order, so they exist + even without edges; other labels follow in order of first appearance. + ``symbolic`` selects object (sympy) storage and is inferred from the retained + edge weights when omitted. Returns the matrices and the ordered state list. + """ + ordered = list(states) if states is not None else list(model.states()) + index = {state: i for i, state in enumerate(ordered)} + n = len(ordered) + edges: list[tuple[int, int, Any, Any]] = [] + for state in ordered: + for transition in model.graph.out_transitions(state): + j = index.get(transition.target) + key = label(transition) + if j is None or key is SKIP: + continue + edges.append((index[state], j, key, transition.data.get(attr, 0.0))) + if symbolic is None: + symbolic = has_symbolic(value for *_, value in edges) + + matrices: dict[Any, np.ndarray] = {key: zeros((n, n), symbolic=symbolic) for key in labels} + for i, j, key, value in edges: + matrix = matrices.get(key) + if matrix is None: + matrix = matrices[key] = zeros((n, n), symbolic=symbolic) + if symbolic: + matrix[i, j] = as_prob(matrix[i, j]) + as_prob(value) + else: + matrix[i, j] += float(value) + return matrices, ordered + + +def emission_label(transition: Any) -> Any: + """:func:`accumulate_matrices` label keying an edge by its emitted symbol (unlabeled edges skipped).""" + emission = transition.data.get(ATTR_EMISSION) + return SKIP if emission is None else emission + + +def _as_mealy(model: Any) -> Any: + to_mealy = getattr(model, "to_mealy", None) + return to_mealy() if to_mealy is not None else model + + +def _is_symbolic(mealy: Any) -> bool: + edge_probs = [transition.data.get(ATTR_PROB, 0.0) for transition in mealy.transitions()] + return has_symbolic(edge_probs) or has_symbolic(mealy.initial_distribution.values()) + + +def symbol_matrices(mealy: Any) -> dict[Any, np.ndarray]: + """Return symbol -> joint transition matrix ``T^(x)`` of a Mealy-style model. + + Symbols (the observation alphabet plus every emitted symbol) are keyed in + ``repr``-sorted order, so iteration and seeded sampling are reproducible. + Matrices are object-dtype when any edge or initial probability is symbolic. + """ + emissions = {transition.data.get(ATTR_EMISSION) for transition in mealy.transitions()} - {None} + symbols = sorted(set(getattr(mealy, "observation_alphabet", ())) | emissions, key=repr) + matrices, _states = accumulate_matrices( + mealy, + attr=ATTR_PROB, + states=mealy.reindex().states, + label=emission_label, + labels=symbols, + symbolic=_is_symbolic(mealy), + ) + return matrices + + +def limit_distribution(pi_initial: np.ndarray, transition: np.ndarray) -> np.ndarray | None: + """Return the limiting occupation law of ``pi_initial`` under ``transition``. + + Power-iterates the normalized law until it is invariant. Returns ``None`` when + the mass vanishes; on periodic chains the result after the iteration cap need + not be invariant, so callers should check it. + """ + pi = np.asarray(pi_initial, dtype=float).copy() + total = float(pi.sum()) + if total <= 0.0: + return None + pi /= total + matrix = np.asarray(transition, dtype=float) + n = len(pi) + for _ in range(max(100, 20 * n)): + nxt = pi @ matrix + mass = float(nxt.sum()) + if mass <= 0.0: + return None + nxt /= mass + if np.allclose(nxt, pi, rtol=1e-12, atol=1e-14): + pi = nxt + break + pi = nxt + pi[np.isclose(pi, 0.0, atol=1e-15)] = 0.0 + mass = float(pi.sum()) + if mass <= 0.0: + return None + return pi / mass + + +def _model_start(mealy: Any) -> np.ndarray: + idx = mealy.reindex() + n = len(idx) + symbolic = _is_symbolic(mealy) + if mealy.initial_distribution: + pi = zeros((n,), symbolic=symbolic) + for state, mass in mealy.initial_distribution.items(): + pi[idx.index(state)] = as_prob(mass) + return pi + if n: + return np.asarray(mealy.stationary_distribution(), dtype=object if symbolic else float) + return zeros((0,), symbolic=symbolic) + + +def _stationary_start(mealy: Any, joint: Mapping[Any, np.ndarray]) -> np.ndarray: + from sofic.generators.stationary import stationary_distribution_from_transition + + pi_initial = _model_start(mealy) + n = len(pi_initial) + if n == 0: + return pi_initial + symbolic = pi_initial.dtype == object or any(matrix.dtype == object for matrix in joint.values()) + transition = zeros((n, n), symbolic=symbolic) + for matrix in joint.values(): + transition = transition + matrix + if not symbolic: + limited = limit_distribution(pi_initial, transition) + if limited is not None and np.allclose(limited @ transition, limited, rtol=1e-8, atol=1e-10): + return limited + try: + return stationary_distribution_from_transition(transition) + except (np.linalg.LinAlgError, ValueError): + return pi_initial + + +def _explicit_start(mealy: Any, start: Hashable | Mapping[Hashable, Any] | Sequence[Any] | np.ndarray) -> np.ndarray: + idx = mealy.reindex() + n = len(idx) + if isinstance(start, Mapping): + vector = zeros((n,), symbolic=has_symbolic(start.values())) + for state, mass in start.items(): + if not mealy.graph.has_state(state): + raise ValueError(f"unknown start state {state!r}") + vector[idx.index(state)] = as_prob(mass) + return vector + + if isinstance(start, Hashable) and mealy.graph.has_state(start): + vector = zeros((n,)) + vector[idx.index(start)] = 1.0 + return vector + + values = list(np.asarray(start, dtype=object).ravel()) + vector = np.asarray(values, dtype=object if has_symbolic(values) else float) + if vector.shape != (n,): + raise ValueError(f"start vector has length {vector.size}, expected {n}") + return vector + + +def start_vector(model: Any, start: StartSpec = None, *, policy: StartPolicy = "model") -> np.ndarray: + """Return the initial state vector of ``model`` in ``model.to_mealy().reindex()`` order. + + ``start`` may be ``None`` (use ``policy``; see the module docstring), a state + (point mass), a state -> mass mapping, or a dense vector in state order. + Explicit starts are not renormalized. The vector is object-dtype when any + supplied or model probability is symbolic. + """ + mealy = _as_mealy(model) + if start is not None: + return _explicit_start(mealy, start) + if policy == "model": + return _model_start(mealy) + if policy == "stationary": + return _stationary_start(mealy, symbol_matrices(mealy)) + raise ValueError(f"unknown start policy {policy!r}; expected 'model' or 'stationary'") + + +def emission_tensors(model: Any, *, policy: StartPolicy = "model") -> tuple[np.ndarray, dict[Any, np.ndarray]]: + """Return ``(pi, T)``: the ``policy`` start vector and symbol -> ``T^(x)`` matrices. + + ``model`` is converted with ``to_mealy()`` when it provides one, and both + outputs share that presentation's ``reindex()`` state order. + """ + mealy = _as_mealy(model) + joint = symbol_matrices(mealy) + if policy == "model": + return _model_start(mealy), joint + if policy == "stationary": + return _stationary_start(mealy, joint), joint + raise ValueError(f"unknown start policy {policy!r}; expected 'model' or 'stationary'") diff --git a/sofic/generators/measures.py b/sofic/generators/measures.py index d6deee6..7220ee7 100644 --- a/sofic/generators/measures.py +++ b/sofic/generators/measures.py @@ -84,24 +84,23 @@ def joint_block_distribution( ``history_length`` counts symbols before the present symbol, so the emitted block length is ``history_length + 1``. """ - from itertools import product - - from sofic.generators.hmm_inference import _stationary_emission_tensors + from sofic.generators.matrices import emission_tensors + from sofic.generators.words import _enumerate_words, _matrix_step dit = _require_dit() # Blocks of a stationary process are weighted by the stationary state law, not # the model's initial distribution (which may describe only the transient). - pi, joint = _stationary_emission_tensors(generator) + pi, joint = emission_tensors(generator, policy="stationary") symbol_list = sorted(generator.observation_alphabet, key=repr) block_length = max(1, history_length + 1) ones = np.ones(len(pi), dtype=float) - outcomes = list(product(symbol_list, repeat=block_length)) + outcomes = [] probs = [] - for outcome in outcomes: - mass = pi.copy() - for symbol in outcome: - mass = mass @ joint.get(symbol, np.zeros((len(pi), len(pi)), dtype=float)) + for outcome, mass in _enumerate_words( + symbol_list, block_length, pi.copy(), _matrix_step(joint, len(pi), prune=False) + ): + outcomes.append(outcome) probs.append(float(mass @ ones)) total = sum(probs) @@ -168,7 +167,7 @@ def entropy_rate_hmm(hmm: HiddenMarkovModel) -> Any: Returns a sympy :class:`~sympy.Expr` when the stationary law or emission tensors are symbolic; otherwise a Python ``float``. """ - from sofic.generators.hmm_inference import _emission_transition_tensors + from sofic.generators.matrices import symbol_matrices from sofic.generators.prob import ( array_sum, as_prob, @@ -186,7 +185,7 @@ def entropy_rate_hmm(hmm: HiddenMarkovModel) -> Any: dit = _require_dit() idx = hmm.reindex() pi = hmm.stationary_distribution() - _, joint = _emission_transition_tensors(hmm) + joint = symbol_matrices(hmm.to_mealy()) symbolic = pi.dtype == object or has_symbolic(pi.ravel()) if not symbolic: diff --git a/sofic/generators/mixed_state_construction.py b/sofic/generators/mixed_state_construction.py index 4968fc0..840bdd0 100644 --- a/sofic/generators/mixed_state_construction.py +++ b/sofic/generators/mixed_state_construction.py @@ -10,6 +10,7 @@ from sofic.exceptions import MixedStateExplosionError from sofic.generators.base import HiddenMarkovModel +from sofic.generators.matrices import start_vector, symbol_matrices from sofic.generators.mealy import MealyHMM from sofic.generators.mixed_state import ( MixedState, @@ -18,12 +19,10 @@ ) from sofic.generators.prob import ( as_prob, - has_symbolic, is_positive_mass, matvec, simplify_prob, sum_probs, - zeros, ) from sofic.graph import ATTR_EMISSION, ATTR_PROB, TransitionGraph @@ -41,14 +40,8 @@ def _resolve_initial_belief( vector = hmm.stationary_distribution() elif isinstance(initial_mixed_state, MixedState): vector = initial_mixed_state.as_array() - elif isinstance(initial_mixed_state, Mapping): - index = {state: i for i, state in enumerate(basis)} - symbolic = has_symbolic(initial_mixed_state.values()) - vector = zeros((len(basis),), symbolic=symbolic) - for state, mass in initial_mixed_state.items(): - vector[index[state]] = as_prob(mass) else: - vector = np.asarray(initial_mixed_state, dtype=object if has_symbolic(initial_mixed_state) else float) + vector = start_vector(hmm, initial_mixed_state) if len(np.asarray(vector).ravel()) != len(basis): raise ValueError(f"initial belief has length {len(np.asarray(vector).ravel())}, expected {len(basis)}") @@ -84,13 +77,11 @@ def build_mixed_state_presentation( if not isinstance(hmm, MealyHMM): raise TypeError(f"mixed-state presentation requires a MealyHMM, not {type(hmm)!r}") - from sofic.generators.hmm_inference import _emission_transition_tensors - constraints = getattr(hmm, "symbol_constraints", None) idx = hmm.reindex() basis = idx.states - _, joint = _emission_transition_tensors(hmm) + joint = symbol_matrices(hmm) symbols = tuple(sorted(joint, key=str)) eta0 = _resolve_initial_belief(hmm, basis, initial_mixed_state) diff --git a/sofic/generators/process_equivalence.py b/sofic/generators/process_equivalence.py index f1c22d8..bd5881f 100644 --- a/sofic/generators/process_equivalence.py +++ b/sofic/generators/process_equivalence.py @@ -14,8 +14,7 @@ import numpy as np from sofic.generators.base import HiddenMarkovModel -from sofic.generators.hmm_inference import _emission_transition_tensors_from_mealy -from sofic.generators.words import _start_vector +from sofic.generators.matrices import start_vector, symbol_matrices _DEFAULT_RTOL = 1e-9 _DEFAULT_ATOL = 1e-12 @@ -77,11 +76,10 @@ def from_hmm( start: Hashable | Mapping[Hashable, float] | Sequence[float] | np.ndarray | None = None, ) -> _HistoryFutureWordList: mealy = hmm.to_mealy() - pi, matrices = _emission_transition_tensors_from_mealy(mealy) return cls( alphabet=tuple(sorted(mealy.observation_alphabet, key=repr)), - matrices=matrices, - start=_start_vector(mealy, pi, start), + matrices=symbol_matrices(mealy), + start=np.asarray(start_vector(mealy, start), dtype=float), ) @property diff --git a/sofic/generators/quasi_inference.py b/sofic/generators/quasi_inference.py index 2a467d2..5674965 100644 --- a/sofic/generators/quasi_inference.py +++ b/sofic/generators/quasi_inference.py @@ -9,21 +9,14 @@ from sofic.exceptions import QuasiStochasticValidationError from sofic.generators.base import QuasiStochasticModel -from sofic.graph import ATTR_EMISSION, ATTR_QUASIPROB +from sofic.generators.matrices import accumulate_matrices, emission_label +from sofic.graph import ATTR_QUASIPROB def transition_matrices(model: QuasiStochasticModel) -> dict[Any, np.ndarray]: - idx = model.reindex() - n = len(idx) - matrices: dict[Any, np.ndarray] = {} - for transition in model.transitions(): - emission = transition.data.get(ATTR_EMISSION) - if emission is None: - continue - matrix = matrices.setdefault(emission, np.zeros((n, n), dtype=float)) - i = idx.index(transition.source) - j = idx.index(transition.target) - matrix[i, j] += float(transition.data.get(ATTR_QUASIPROB, 0.0)) + matrices, _states = accumulate_matrices( + model, attr=ATTR_QUASIPROB, states=model.reindex().states, label=emission_label, symbolic=False + ) return matrices diff --git a/sofic/generators/stack_inference.py b/sofic/generators/stack_inference.py index 6028214..03bd1f6 100644 --- a/sofic/generators/stack_inference.py +++ b/sofic/generators/stack_inference.py @@ -3,11 +3,12 @@ from __future__ import annotations from collections import Counter, defaultdict -from collections.abc import Callable, Hashable, Sequence +from collections.abc import Callable, Hashable, Mapping, Sequence from typing import Any, ClassVar, Literal from sofic.automata.papni import DyckAlphabet, is_well_matched, learn_sofic_dyck_shift_papni from sofic.exceptions import StochasticValidationError +from sofic.generators._suffix_counts import infer_alphabet, iter_suffixes from sofic.generators.epsilon_inference import ( History, MorphTest, @@ -15,9 +16,7 @@ _bonferroni_alpha, _cluster_histories_by_morph, _cssr_default_lmax, - _cssr_determinize, - _drop_transient_states, - _merge_similar_states, + _recurrent_states, morph_test_score, morphs_differ, suggest_lmax, @@ -71,28 +70,30 @@ def from_sequence( # type: ignore[override] seq = tuple(sequence) if not seq: raise ValueError("sequence must be non-empty") - visible_alphabet = tuple(sorted(alphabet.symbol_alphabet, key=repr)) - max_len = max_length if max_length is not None else len(seq) - counts = cls(alphabet=visible_alphabet) - stack: list[Any] = [] - for t, symbol in enumerate(seq): + stacks: list[tuple[Any, ...]] = [] + stack: tuple[Any, ...] = () + for symbol in seq: if symbol not in alphabet.symbol_alphabet: raise ValueError(f"symbol {symbol!r} not in Dyck alphabet") - for length in range(0, min(t, max_len) + 1): - suffix = seq[t - length : t] - history = (suffix, tuple(stack)) - counts.history_counts[history] += 1 - counts.next_counts[history][symbol] += 1 - if symbol in alphabet.call_alphabet: - if len(stack) >= max_stack_depth: - stack = stack[1:] - stack.append(symbol) - elif symbol in alphabet.return_alphabet: - if stack: - stack.pop() + stacks.append(stack) + stack = _push(stack, symbol, alphabet=alphabet, max_stack_depth=max_stack_depth) + counts = cls(alphabet=infer_alphabet(alphabet.symbol_alphabet, None)) + for t, suffix in iter_suffixes(seq, max_length if max_length is not None else len(seq)): + history = (suffix, stacks[t]) + counts.history_counts[history] += 1 + counts.next_counts[history][seq[t]] += 1 return counts +def _push(stack: tuple[Any, ...], symbol: Any, *, alphabet: DyckAlphabet, max_stack_depth: int) -> tuple[Any, ...]: + """The stack after ``symbol``: calls push (dropping the bottom past ``max_stack_depth``), returns pop.""" + if symbol in alphabet.call_alphabet: + return (*(stack[1:] if len(stack) >= max_stack_depth else stack), symbol) + if symbol in alphabet.return_alphabet: + return stack[:-1] + return stack + + def _successor_history( history: ConfigurationHistory, symbol: Any, @@ -109,15 +110,7 @@ def _successor_history( new_suffix = extended else: new_suffix = extended[-length:] - - stack_list = list(stack) - if symbol in alphabet.call_alphabet: - if len(stack_list) >= max_stack_depth: - stack_list = stack_list[1:] - stack_list.append(symbol) - elif symbol in alphabet.return_alphabet and stack_list: - stack_list.pop() - return new_suffix, tuple(stack_list) + return new_suffix, _push(stack, symbol, alphabet=alphabet, max_stack_depth=max_stack_depth) def _stack_successor_fn( @@ -228,14 +221,48 @@ def _stack_determinize( alphabet: DyckAlphabet, max_stack_depth: int, ) -> dict[int, set[ConfigurationHistory]]: - """Split homogeneous states until stack-lifted transitions are unifilar.""" - return _cssr_determinize( - states, - history_to_state, - counts, - length=length, - successor_fn=_stack_successor_fn(alphabet=alphabet, length=length, max_stack_depth=max_stack_depth), - ) + """Split homogeneous states until stack-lifted transitions are unifilar. + + Updates ``history_to_state`` in place for the histories that move. + """ + successor = _stack_successor_fn(alphabet=alphabet, length=length, max_stack_depth=max_stack_depth) + current = {state_id: set(histories) for state_id, histories in states.items()} + changed = True + next_state_id = max(current) + 1 if current else 0 + + while changed: + changed = False + for state_id in sorted(current): + histories = current[state_id] + if len(histories) <= 1: + continue + for symbol in counts.alphabet: + buckets: dict[int, set[ConfigurationHistory]] = defaultdict(set) + for history in histories: + if counts.next_counts.get(history, Counter()).get(symbol, 0) == 0: + continue + target = history_to_state.get(successor(history, symbol)) + if target is None: + continue + buckets[target].add(history) + if len(buckets) <= 1: + continue + # Keep the largest bucket -- and every history that never emits + # ``symbol`` -- in the original state; split the other buckets off. + ordered = sorted(buckets.items(), key=lambda item: (-len(item[1]), min(item[1]))) + moved = set().union(*(split for _target, split in ordered[1:])) + current[state_id] = histories - moved + for _target, split_histories in ordered[1:]: + new_id = next_state_id + next_state_id += 1 + current[new_id] = split_histories + for history in split_histories: + history_to_state[history] = new_id + changed = True + break + if changed: + break + return current def _stack_merge( @@ -247,8 +274,33 @@ def _stack_merge( test: MorphTest, alphabet: DyckAlphabet, ) -> dict[int, set[ConfigurationHistory]]: - proxy = _control_counts(counts, alphabet).restricted_to(set(history_to_state)) - return _merge_similar_states(states, history_to_state, proxy, alpha=alpha, test=test) + """Merge states whose pooled control morphs (see :func:`_control_counts`) are indistinguishable. + + Merging on the morph alone can fuse states with incompatible ``symbol -> + successor`` maps, so the result need not be unifilar. + """ + control = _control_counts(counts, alphabet) + current = {state_id: set(histories) for state_id, histories in states.items()} + changed = True + while changed: + changed = False + state_ids = sorted(current) + for index, left_id in enumerate(state_ids): + if left_id not in current: + continue + for right_id in state_ids[index + 1 :]: + if right_id not in current: + continue + if morphs_differ(control, current[left_id], current[right_id], alpha=alpha, test=test): + continue + current[left_id].update(current.pop(right_id)) + for history in current[left_id]: + history_to_state[history] = left_id + changed = True + break + if changed: + break + return current def _stack_drop_transient( @@ -260,14 +312,32 @@ def _stack_drop_transient( alphabet: DyckAlphabet, max_stack_depth: int, ) -> dict[int, set[ConfigurationHistory]]: - proxy = counts.restricted_to(set(history_to_state)) - return _drop_transient_states( - states, - history_to_state, - proxy, - length=length, - successor_fn=_stack_successor_fn(alphabet=alphabet, length=length, max_stack_depth=max_stack_depth), - ) + """Keep only states in closed communicating classes; all states if there are none. + + States left with no edge into the kept set (their successor configurations + were too rare to be placed) are pruned repeatedly, since a generator state + needs outgoing mass. + """ + successor = _stack_successor_fn(alphabet=alphabet, length=length, max_stack_depth=max_stack_depth) + + def state_edges(kept: Mapping[int, set[ConfigurationHistory]]) -> dict[int, dict[Any, set[int]]]: + edges: dict[int, dict[Any, set[int]]] = {state_id: defaultdict(set) for state_id in kept} + for state_id, histories in kept.items(): + for history in histories: + for symbol, count in counts.next_counts.get(history, Counter()).items(): + target = history_to_state.get(successor(history, symbol)) if count else None + if target in kept: + edges[state_id][symbol].add(target) + return edges + + recurrent = set().union(*_recurrent_states(state_edges(states))) + kept = {state_id: histories for state_id, histories in states.items() if not recurrent or state_id in recurrent} + while True: + edges = state_edges(kept) + dead = {state_id for state_id, by_symbol in edges.items() if not any(by_symbol.values())} + if not dead or len(dead) == len(kept): + return kept + kept = {state_id: histories for state_id, histories in kept.items() if state_id not in dead} def _counts_to_stack_hmm( @@ -282,21 +352,15 @@ def _counts_to_stack_hmm( ) -> HiddenMarkovStackModel: visits: Counter[int] = Counter() seq = tuple(sequence) - stack: list[Any] = [] + stack: tuple[Any, ...] = () for t in range(len(seq)): # Each step occupies one state: the one keyed by its longest available suffix. for hist_len in range(min(t, length), -1, -1): - state = history_to_state.get((seq[t - hist_len : t], tuple(stack))) + state = history_to_state.get((seq[t - hist_len : t], stack)) if state is not None: visits[state] += 1 break - symbol = seq[t] - if symbol in alphabet.call_alphabet: - if len(stack) >= max_stack_depth: - stack = stack[1:] - stack.append(symbol) - elif symbol in alphabet.return_alphabet and stack: - stack.pop() + stack = _push(stack, seq[t], alphabet=alphabet, max_stack_depth=max_stack_depth) if not visits: raise StochasticValidationError("no empirical configuration visits") diff --git a/sofic/generators/words.py b/sofic/generators/words.py index 6fcffec..2cb2ba1 100644 --- a/sofic/generators/words.py +++ b/sofic/generators/words.py @@ -2,26 +2,72 @@ from __future__ import annotations -from collections.abc import Hashable, Mapping, Sequence -from itertools import product +from collections.abc import Callable, Hashable, Iterator, Mapping, Sequence from typing import Any import numpy as np from sofic.generators.base import HiddenMarkovModel, QuasiStochasticModel -from sofic.generators.hmm_inference import _emission_transition_tensors from sofic.generators.markov import MarkovChain +from sofic.generators.matrices import emission_tensors, start_vector, symbol_matrices from sofic.generators.pfa import ProbabilisticFiniteAutomaton from sofic.graph import ATTR_PROB _TOL = 1e-15 +def _enumerate_words( + alphabet: Sequence[Any], + length: int, + start: Any = None, + step: Callable[[Any, Any], Any] | None = None, +) -> Iterator[tuple[tuple[Any, ...], Any]]: + """Yield ``(word, mass)`` for every word of ``length`` over ``alphabet``. + + Words come in lexicographic (``itertools.product``) order of ``alphabet``. + ``mass`` starts at ``start`` and is advanced by ``step(mass, symbol)`` once per + symbol, so prefixes are propagated once rather than per word; a step returning + ``None`` prunes every extension of that prefix. With ``step=None`` the mass is + carried unchanged and every word is yielded. + """ + + def walk(prefix: tuple[Any, ...], mass: Any, remaining: int) -> Iterator[tuple[tuple[Any, ...], Any]]: + if remaining == 0: + yield prefix, mass + return + for symbol in alphabet: + if step is None: + yield from walk(prefix + (symbol,), mass, remaining - 1) + continue + nxt = step(mass, symbol) + if nxt is not None: + yield from walk(prefix + (symbol,), nxt, remaining - 1) + + yield from walk((), start, length) + + +def _matrix_step(matrices: Mapping[Any, np.ndarray], n: int, *, prune: bool = True) -> Callable[[np.ndarray, Any], Any]: + """Return a :func:`_enumerate_words` step multiplying a row mass by ``matrices[symbol]``. + + Symbols without a matrix act as the zero matrix. With ``prune`` a numeric + prefix whose mass vector is identically zero is dropped. + """ + zero = np.zeros((n, n), dtype=float) + + def step(mass: np.ndarray, symbol: Any) -> np.ndarray | None: + nxt = mass @ matrices.get(symbol, zero) + if prune and nxt.dtype != object and not np.any(nxt): + return None + return nxt + + return step + + def hmm_words_of_length(hmm: HiddenMarkovModel, length: int) -> dict[tuple[Any, ...], float]: """Return observed words of ``length`` and their probabilities.""" if length < 0: raise ValueError("length must be nonnegative") - pi, joint = _emission_transition_tensors(hmm) + pi, joint = emission_tensors(hmm) alphabet = sorted(hmm.observation_alphabet, key=repr) if length == 0: return {(): float(pi.sum())} if pi.sum() > _TOL else {} @@ -29,12 +75,8 @@ def hmm_words_of_length(hmm: HiddenMarkovModel, length: int) -> dict[tuple[Any, return {} terminal = np.ones(len(pi), dtype=float) - zero = np.zeros((len(pi), len(pi)), dtype=float) distribution: dict[tuple[Any, ...], float] = {} - for word in product(alphabet, repeat=length): - mass = pi.copy() - for symbol in word: - mass = mass @ joint.get(symbol, zero) + for word, mass in _enumerate_words(alphabet, length, pi.copy(), _matrix_step(joint, len(pi))): probability = float(mass @ terminal) if abs(probability) > _TOL: distribution[word] = probability @@ -53,8 +95,8 @@ def hmm_word_probability( a state-probability mapping, or a dense vector in the model's state order. """ mealy = hmm.to_mealy() - pi, joint = _emission_transition_tensors(mealy) - mass = _start_vector(mealy, pi, start) + joint = symbol_matrices(mealy) + mass = np.asarray(start_vector(mealy, start), dtype=float) if len(word) == 0: return float(mass.sum()) n = len(mass) @@ -101,7 +143,7 @@ def hmm_word_probabilities( if not sparse or abs(probability) > _TOL: distribution[()] = probability continue - for word in product(alphabet, repeat=length): + for word, _ in _enumerate_words(alphabet, length): probability = hmm_word_probability(mealy, word, start=start) if not sparse or abs(probability) > _TOL: distribution[word] = probability @@ -134,7 +176,7 @@ def pfa_words_of_length(pfa: ProbabilisticFiniteAutomaton, length: int) -> dict[ if not alphabet: return {} distribution: dict[tuple[Any, ...], float] = {} - for word in product(alphabet, repeat=length): + for word, _ in _enumerate_words(alphabet, length): probability = pfa.string_probability(word) if probability > _TOL: distribution[word] = probability @@ -152,7 +194,7 @@ def quasi_words_of_length(model: QuasiStochasticModel, length: int) -> dict[tupl if not alphabet: return {} distribution: dict[tuple[Any, ...], float] = {} - for word in product(sorted(alphabet, key=repr), repeat=length): + for word, _ in _enumerate_words(sorted(alphabet, key=repr), length): probability = float(model.word_probability(word)) if abs(probability) > _TOL: distribution[word] = probability @@ -169,7 +211,7 @@ def markov_words_of_length(chain: MarkovChain, length: int) -> dict[tuple[Hashab total = float(sum(start.values())) return {(): total} if total > _TOL else {} distribution: dict[tuple[Hashable, ...], float] = {} - for word in product(states, repeat=length): + for word, _ in _enumerate_words(states, length): probability = _markov_path_probability(chain, word, start) if probability > _TOL: distribution[word] = probability @@ -184,35 +226,6 @@ def _quasi_alphabet(model: QuasiStochasticModel) -> tuple[Any, ...]: return tuple(model.transition_matrices()) -def _start_vector( - hmm: HiddenMarkovModel, - default: np.ndarray, - start: Hashable | Mapping[Hashable, float] | Sequence[float] | np.ndarray | None, -) -> np.ndarray: - if start is None: - return np.array(default, dtype=float) - - idx = hmm.reindex() - n = len(idx) - if isinstance(start, Mapping): - vector = np.zeros(n, dtype=float) - for state, mass in start.items(): - if not hmm.graph.has_state(state): - raise ValueError(f"unknown start state {state!r}") - vector[idx.index(state)] = float(mass) - return vector - - if hmm.graph.has_state(start): - vector = np.zeros(n, dtype=float) - vector[idx.index(start)] = 1.0 - return vector - - vector = np.asarray(start, dtype=float) - if vector.shape != (n,): - raise ValueError(f"start vector must have shape {(n,)}, got {vector.shape}") - return vector.copy() - - def _markov_start(chain: MarkovChain) -> dict[Hashable, float]: """Initial law of ``chain``, or its stationary law when none is given.""" if chain.initial_distribution: diff --git a/sofic/inference/model_selection.py b/sofic/inference/model_selection.py index 912c177..cc4345c 100644 --- a/sofic/inference/model_selection.py +++ b/sofic/inference/model_selection.py @@ -133,9 +133,9 @@ def _smoothed_log_likelihood( computed by forward filtering. After a symbol the model forbids, the belief is propagated without conditioning on it. """ - from sofic.generators.hmm_inference import _emission_transition_tensors_from_mealy + from sofic.generators.matrices import emission_tensors - pi, joint = _emission_transition_tensors_from_mealy(model.to_mealy()) + pi, joint = emission_tensors(model) total_step = sum(joint.values()) belief = np.asarray(pi, dtype=float) belief = belief / belief.sum() diff --git a/sofic/properties.py b/sofic/properties.py index 4edf4c0..bc211f8 100644 --- a/sofic/properties.py +++ b/sofic/properties.py @@ -294,26 +294,10 @@ def transition_matrix( are ignored), otherwise all model states in iteration order. Returns the matrix together with the ordered state list defining its axes. """ - ordered = list(states) if states is not None else list(model.states()) - index = {state: i for i, state in enumerate(ordered)} - n = len(ordered) - from sofic.generators.prob import as_prob, has_symbolic, zeros - - edge_probs = [] - for state in ordered: - for transition in model.graph.out_transitions(state): - if transition.target in index: - edge_probs.append(transition.data.get(attr, 0.0)) - symbolic = has_symbolic(edge_probs) - matrix = zeros((n, n), symbolic=symbolic) - for state in ordered: - i = index[state] - for transition in model.graph.out_transitions(state): - j = index.get(transition.target) - if j is None: - continue - matrix[i, j] = as_prob(matrix[i, j]) + as_prob(transition.data.get(attr, 0.0)) - return matrix, ordered + from sofic.generators.matrices import accumulate_matrices + + matrices, ordered = accumulate_matrices(model, attr=attr, states=states, labels=(None,)) + return matrices[None], ordered def _initial_vector_and_transition(model: StateMachine) -> tuple[np.ndarray, np.ndarray]: @@ -330,10 +314,9 @@ def _labeled_or_internal_matrices(model: StateMachine) -> list[np.ndarray]: observation_alphabet = getattr(model, "observation_alphabet", None) to_mealy = getattr(model, "to_mealy", None) if observation_alphabet is not None and to_mealy is not None: - from sofic.generators.hmm_inference import _emission_transition_tensors_from_mealy + from sofic.generators.matrices import symbol_matrices - _pi, matrices = _emission_transition_tensors_from_mealy(to_mealy()) - return list(matrices.values()) + return list(symbol_matrices(to_mealy()).values()) _initial, transition = _initial_vector_and_transition(model) return [transition] diff --git a/sofic/viz/_format.py b/sofic/viz/_format.py index 7b8ed91..32c9fc4 100644 --- a/sofic/viz/_format.py +++ b/sofic/viz/_format.py @@ -3,10 +3,10 @@ from __future__ import annotations from collections.abc import Mapping, Sequence +from fractions import Fraction from typing import Any from sofic.viz import _labels -from sofic.viz._rational import two_digit_rational def dot_escape(text: str) -> str: @@ -34,42 +34,26 @@ def format_prob(value: float, *, precision: int = 3) -> str: return dot_escape(f"{value:.{precision}g}") +def _dot_fraction(frac: Fraction) -> str: + return dot_escape(f"{frac.numerator}/{frac.denominator}") + + +def _dot_sympy(expr: Any) -> str: + import sympy as sp + + return dot_escape(sp.sstr(expr)) + + def format_prob_rational(value: float, *, precision: int = 3) -> str: """Format a probability as p/q when exact with two-digit numerator and denominator.""" - if value <= 0.0: - return "0" - if value >= 1.0: - return "1" - frac = two_digit_rational(value) - if frac is not None: - if frac.numerator == frac.denominator: - return "1" - return dot_escape(f"{frac.numerator}/{frac.denominator}") - return format_prob(value, precision=precision) + return format_prob_label(float(value), precision=precision) def format_prob_label(value: Any, *, precision: int = 3) -> str: """Format a probability for Graphviz edge/π labels (float or sympy Expr).""" - try: - from sofic.generators.prob import is_symbolic, simplify_prob - except ImportError: # pragma: no cover - - def is_symbolic(_v: Any) -> bool: - return False - - def simplify_prob(v: Any) -> Any: - return v - - if is_symbolic(value): - simplified = simplify_prob(value) - try: - import sympy as sp - - text = sp.sstr(simplified) - except Exception: - text = str(simplified) - return dot_escape(text) - return format_prob_rational(float(value), precision=precision) + return _labels.format_prob( + value, escape=dot_escape, fraction=_dot_fraction, symbolic=_dot_sympy, precision=precision + ) def format_distribution(dist: Mapping[Any, Any], *, precision: int = 3) -> str: diff --git a/sofic/viz/_labels.py b/sofic/viz/_labels.py index 0cf3b3a..8f3c91f 100644 --- a/sofic/viz/_labels.py +++ b/sofic/viz/_labels.py @@ -1,18 +1,37 @@ -"""Backend-agnostic skeletons for state/symbol label formatting. +"""Backend-agnostic skeletons for state/symbol/probability label formatting. The graphviz and LaTeX backends share the same structural recursion for formatting states (tuples recurse, frozensets render as sorted ``{...}``, -``EPSILON`` becomes a glyph) and symbols; they differ only in the leaf escape -function and the epsilon glyph. These helpers capture that shared shape. +``EPSILON`` becomes a glyph), symbols, and probabilities (symbolic +expressions, clamped ``0`` / ``1``, exact two-digit fractions, else +``precision`` significant digits); they differ only in the leaf escape +function, the epsilon glyph, and how fractions / sympy expressions are +spelled. These helpers capture that shared shape. """ from __future__ import annotations from collections.abc import Callable +from fractions import Fraction from typing import Any from sofic.graph import EPSILON +_TWO_DIGIT_RATIONAL_ATOL = 1e-9 +_MAX_TWO_DIGIT_RATIONAL = 99 + + +def two_digit_rational(value: float, *, atol: float = _TWO_DIGIT_RATIONAL_ATOL) -> Fraction | None: + """Return a reduced rational with 1 <= p, q <= 99 when ``value`` matches exactly.""" + if value <= 0.0 or value >= 1.0: + return None + frac = Fraction(value).limit_denominator(_MAX_TWO_DIGIT_RATIONAL) + if abs(float(frac) - value) >= atol: + return None + if not (1 <= frac.numerator <= _MAX_TWO_DIGIT_RATIONAL and 1 <= frac.denominator <= _MAX_TWO_DIGIT_RATIONAL): + return None + return frac + def format_state(state: Any, *, escape: Callable[[str], str], epsilon: str) -> str: """Format a (possibly nested) state using ``escape`` for leaves.""" @@ -32,3 +51,39 @@ def format_symbol(symbol: Any, *, escape: Callable[[str], str], epsilon: str) -> if symbol is EPSILON: return epsilon return escape(str(symbol)) + + +def format_prob( + value: Any, + *, + escape: Callable[[str], str], + fraction: Callable[[Fraction], str], + symbolic: Callable[[Any], str], + precision: int = 3, +) -> str: + """Format a probability (float or sympy Expr). + + ``symbolic`` spells a simplified sympy expression and ``fraction`` an exact + two-digit rational, each returning backend-ready text. Decimals and the + ``str`` fallback (used when ``symbolic`` raises) pass through ``escape``. + """ + from sofic.generators.prob import is_symbolic, simplify_prob + + if is_symbolic(value): + simplified = simplify_prob(value) + try: + return symbolic(simplified) + except Exception: + return escape(str(simplified)) + + numeric = float(value) + if numeric <= 0.0: + return "0" + if numeric >= 1.0: + return "1" + frac = two_digit_rational(numeric) + if frac is not None: + if frac.numerator == frac.denominator: + return "1" + return fraction(frac) + return escape(f"{numeric:.{precision}g}") diff --git a/sofic/viz/_rational.py b/sofic/viz/_rational.py deleted file mode 100644 index e961660..0000000 --- a/sofic/viz/_rational.py +++ /dev/null @@ -1,20 +0,0 @@ -"""Shared rational-number helper for viz label formatting.""" - -from __future__ import annotations - -from fractions import Fraction - -_TWO_DIGIT_RATIONAL_ATOL = 1e-9 -_MAX_TWO_DIGIT_RATIONAL = 99 - - -def two_digit_rational(value: float, *, atol: float = _TWO_DIGIT_RATIONAL_ATOL) -> Fraction | None: - """Return a reduced rational with 1 <= p, q <= 99 when ``value`` matches exactly.""" - if value <= 0.0 or value >= 1.0: - return None - frac = Fraction(value).limit_denominator(_MAX_TWO_DIGIT_RATIONAL) - if abs(float(frac) - value) >= atol: - return None - if not (1 <= frac.numerator <= _MAX_TWO_DIGIT_RATIONAL and 1 <= frac.denominator <= _MAX_TWO_DIGIT_RATIONAL): - return None - return frac diff --git a/sofic/viz/_tikz_format.py b/sofic/viz/_tikz_format.py index 5a5a8f1..74b591f 100644 --- a/sofic/viz/_tikz_format.py +++ b/sofic/viz/_tikz_format.py @@ -3,10 +3,12 @@ from __future__ import annotations from collections.abc import Sequence +from fractions import Fraction from typing import Any from sofic.viz import _labels -from sofic.viz._rational import _TWO_DIGIT_RATIONAL_ATOL, two_digit_rational + +_HALF = Fraction(1, 2) _LATEX_SPECIAL = { "\\": r"\textbackslash{}", @@ -59,64 +61,25 @@ def format_symbol_latex(symbol: Any) -> str: return _labels.format_symbol(symbol, escape=_latex_arg, epsilon=r"\varepsilon") -def format_prob_latex(value: Any, *, precision: int = 3) -> str: - """Format a probability for Vaucanson edge labels (float or sympy Expr).""" - try: - from sofic.generators.prob import is_symbolic, simplify_prob - except ImportError: # pragma: no cover - is_symbolic = lambda _v: False # noqa: E731 - simplify_prob = lambda v: v # noqa: E731 - - if is_symbolic(value): - simplified = simplify_prob(value) - try: - import sympy as sp - - return latex_escape(sp.latex(simplified)) - except Exception: - return latex_escape(str(simplified)) - - numeric = float(value) - if numeric <= 0.0: - return "0" - if numeric >= 1.0: - return "1" - if abs(numeric - 0.5) < _TWO_DIGIT_RATIONAL_ATOL: +def _latex_fraction(frac: Fraction) -> str: + if frac == _HALF: return r"\half" - frac = two_digit_rational(numeric) - if frac is not None: - if frac.numerator == frac.denominator: - return "1" - return rf"\nicefrac{{{frac.numerator}}}{{{frac.denominator}}}" - return latex_escape(f"{numeric:.{precision}g}") - + return rf"\nicefrac{{{frac.numerator}}}{{{frac.denominator}}}" -def format_symbol_macro(symbol: Any) -> str: - return rf"\Symbol{{{format_symbol_latex(symbol)}}}" +def _latex_sympy(expr: Any) -> str: + import sympy as sp -def format_edge_latex(symbol: Any, prob: Any) -> str: - """Return ``$\\Edge{sym}{prob}$`` math content.""" - return rf"$\Edge{{{format_symbol_latex(symbol)}}}{{{format_prob_latex(prob)}}}$" + return sp.latex(expr) -def format_tedge_latex(forward: Any, reverse: Any, prob: Any) -> str: - """Return ``$\\TEdge{f}{r}{prob}$`` math content.""" - return ( - rf"$\TEdge{{{format_symbol_latex(forward)}}}" - rf"{{{format_symbol_latex(reverse)}}}" - rf"{{{format_prob_latex(prob)}}}$" +def format_prob_latex(value: Any, *, precision: int = 3) -> str: + """Format a probability for Vaucanson edge labels (float or sympy Expr).""" + return _labels.format_prob( + value, escape=latex_escape, fraction=_latex_fraction, symbolic=_latex_sympy, precision=precision ) -def format_symbol_only_latex(symbol: Any) -> str: - return rf"$\Symbol{{{format_symbol_latex(symbol)}}}$" - - -def format_transducer_edge_latex(input_symbol: Any, output_symbol: Any | None = None) -> str: - if output_symbol is None: - return format_symbol_only_latex(input_symbol) - return ( - rf"$\Symbol{{{format_symbol_latex(input_symbol)}}}" - rf"\mid\Symbol{{{format_symbol_latex(output_symbol)}}}$" - ) +def format_edge_latex(symbol: Any, prob: Any) -> str: + """Return ``$\\Edge{sym}{prob}$`` math content.""" + return rf"$\Edge{{{format_symbol_latex(symbol)}}}{{{format_prob_latex(prob)}}}$" diff --git a/sofic/viz/tikz.py b/sofic/viz/tikz.py index cc41b84..353499d 100644 --- a/sofic/viz/tikz.py +++ b/sofic/viz/tikz.py @@ -18,18 +18,17 @@ PART_KIND, PART_MATCH_TAG, PART_MULTIPLICITY, - PART_OUTPUT, PART_PROB, PART_QUASIPROB, PART_STACK, PART_SYMBOL, STYLE_DYCK, STYLE_EDGE, - STYLE_PROB_ONLY, - STYLE_SYMBOL_ONLY, + STYLE_FALLBACK, STYLE_TMC, STYLE_TRANSDUCER, STYLE_VPA, + EdgePart, edge_spec, part_value, ) @@ -41,8 +40,6 @@ format_state_latex, format_state_tikz_node, format_symbol_latex, - format_symbol_only_latex, - format_transducer_edge_latex, latex_escape, ) from sofic.viz._tikz_layout import ( @@ -72,77 +69,39 @@ def _format_dyck_match_tag_latex(tag: str) -> str: return latex_escape(tag) +# Symbols render bare in these styles; elsewhere they are wrapped in ``\Symbol{...}``. +_BARE_SYMBOL_STYLES = frozenset({STYLE_VPA, STYLE_TMC}) +_PART_SEPARATOR = {STYLE_TRANSDUCER: r"\mid", STYLE_TMC: r"\mid"} + + +def _render_latex_part(part: EdgePart, style: str) -> str: + if part.kind == PART_KIND: + text = latex_escape(str(part.value)) + return rf"\mathrm{{{text}}}" if style == STYLE_DYCK else text + if part.kind == PART_STACK: + return rf"\uparrow {format_symbol_latex(part.value)}" + if part.kind == PART_MULTIPLICITY: + return rf"\times {part.value}" + if part.kind == PART_MATCH_TAG: + return _format_dyck_match_tag_latex(part.value) + if part.kind in (PART_PROB, PART_QUASIPROB): + return format_prob_latex(part.value) + # symbol / emission / output + symbol = format_symbol_latex(part.value) + return symbol if style in _BARE_SYMBOL_STYLES else rf"\Symbol{{{symbol}}}" + + def _tikz_edge_label(model: StateMachine, transition: Transition) -> str: spec = edge_spec(model, transition) - style = spec.style - - if style == STYLE_SYMBOL_ONLY: - symbol = part_value(spec, PART_SYMBOL) - return format_symbol_only_latex(symbol) if symbol is not None else "" - - if style == STYLE_TRANSDUCER: - symbol = part_value(spec, PART_SYMBOL) - output = part_value(spec, PART_OUTPUT) - if symbol is None and output is None: - return "" - return format_transducer_edge_latex(symbol, output) - - if style == STYLE_PROB_ONLY: - prob = part_value(spec, PART_PROB) - return rf"${format_prob_latex(prob)}$" if prob is not None else "" - - if style == STYLE_EDGE: - label_symbol = part_value(spec, PART_EMISSION, PART_SYMBOL) - value = part_value(spec, PART_PROB, PART_QUASIPROB) - if label_symbol is not None and value is not None: - return format_edge_latex(label_symbol, value) - if label_symbol is not None: - return format_symbol_only_latex(label_symbol) - if value is not None: - return rf"${format_prob_latex(value)}$" + if spec.style in (STYLE_EDGE, STYLE_FALLBACK): + symbol = part_value(spec, PART_EMISSION, PART_SYMBOL) + prob = part_value(spec, PART_PROB, PART_QUASIPROB) + if symbol is not None and prob is not None: + return format_edge_latex(symbol, prob) + parts = [_render_latex_part(part, spec.style) for part in spec.parts] + if not parts: return "" - - if style == STYLE_VPA: - parts: list[str] = [] - for part in spec.parts: - if part.kind == PART_SYMBOL: - parts.append(format_symbol_latex(part.value)) - elif part.kind == PART_KIND: - parts.append(latex_escape(str(part.value))) - elif part.kind == PART_STACK: - parts.append(rf"\uparrow{format_symbol_latex(part.value)}") - return "$" + r"\mid".join(parts) + "$" if parts else "" - - if style == STYLE_DYCK: - parts = [] - for part in spec.parts: - if part.kind == PART_SYMBOL: - parts.append(rf"\Symbol{{{format_symbol_latex(part.value)}}}") - elif part.kind == PART_KIND: - parts.append(rf"\mathrm{{{latex_escape(str(part.value))}}}") - elif part.kind == PART_MATCH_TAG: - parts.append(_format_dyck_match_tag_latex(part.value)) - return "$" + r"\mid ".join(parts) + "$" if parts else "" - - if style == STYLE_TMC: - parts = [] - for part in spec.parts: - if part.kind == PART_SYMBOL: - parts.append(format_symbol_latex(part.value)) - elif part.kind == PART_MULTIPLICITY: - parts.append(latex_escape(f"\\times {part.value}")) - return "$" + r"\mid".join(parts) + "$" if parts else "" - - # STYLE_FALLBACK - label_symbol = part_value(spec, PART_EMISSION, PART_SYMBOL) - prob = part_value(spec, PART_PROB) - if label_symbol is not None and prob is not None: - return format_edge_latex(label_symbol, prob) - if label_symbol is not None: - return format_symbol_only_latex(label_symbol) - if prob is not None: - return rf"${format_prob_latex(prob)}$" - return "" + return "$" + _PART_SEPARATOR.get(spec.style, r"\mid ").join(parts) + "$" def _tikz_display_kwargs(model: StateMachine) -> dict[str, Any]: diff --git a/tests/test_channel_measures.py b/tests/test_channel_measures.py index fd2b5ca..9c5c801 100644 --- a/tests/test_channel_measures.py +++ b/tests/test_channel_measures.py @@ -64,7 +64,7 @@ def test_directed_information_ignores_unreachable_error_class(): """ from sofic.automata.transducer_operations import compose_tg from sofic.generators.directional_flow import directed_information as di_flow - from sofic.generators.hmm_inference import _stationary_emission_tensors + from sofic.generators.matrices import emission_tensors eps = EpsilonTransducer.from_channel(BinaryChannel(0.1, 0.2)) # Force the historical reducible joint even after complete() stops adding an @@ -75,7 +75,7 @@ def test_directed_information_ignores_unreachable_error_class(): completed.add_transition("?", "?", "1", "?", prob=1.0) joint = compose_tg(completed, _iid_input(), joint=True, complete=False) assert ("S", "?") in list(joint.states()) - pi, _tensors = _stationary_emission_tensors(joint) + pi, _tensors = emission_tensors(joint, policy="stationary") idx = joint.reindex() error_index = idx.index(("S", "?")) assert pi[error_index] == pytest.approx(0.0, abs=1e-12) diff --git a/tests/test_epsilon_transducer_inference.py b/tests/test_epsilon_transducer_inference.py index 6e7137f..be5de34 100644 --- a/tests/test_epsilon_transducer_inference.py +++ b/tests/test_epsilon_transducer_inference.py @@ -2,6 +2,9 @@ import numpy as np import pytest +from hypothesis import given, settings +from hypothesis import strategies as st +from hypothesis.extra import numpy as hnp from sofic import EpsilonTransducer, MealyHMM from sofic.automata.transducer_operations import compose_tg @@ -146,3 +149,17 @@ def test_auto_lmax_delay(): xs, ys = _paired_samples(Delay(1), 6000, seed=4) eps = transcssr(xs, ys, input_alphabet=("0", "1"), output_alphabet=("0", "1"), Lmax="auto") assert len(list(eps.states())) == 2 + + +@settings(max_examples=200, deadline=None) +@given( + hnp.arrays(np.int64, st.tuples(st.just(2), st.integers(2, 4)), elements=st.integers(1, 40)), +) +def test_shared_g_statistic_matches_scipy_log_likelihood(table): + """The G-test shared with process CSSR is scipy's log-likelihood statistic, Yates-corrected at dof 1.""" + from scipy import stats + + from sofic.generators._morph_tests import g_statistic + + expected, _p, _dof, _ = stats.chi2_contingency(table, lambda_="log-likelihood") + assert g_statistic(table.astype(float)) == pytest.approx(expected, rel=1e-9, abs=1e-12) diff --git a/tests/test_hmm_inference.py b/tests/test_hmm_inference.py index 390a3a3..99403a6 100644 --- a/tests/test_hmm_inference.py +++ b/tests/test_hmm_inference.py @@ -7,7 +7,6 @@ from sofic.examples import fair_coin, golden_mean from sofic.generators.hmm_inference import ( - _emission_transition_tensors, _forward_scaled, backward, baum_welch, @@ -22,6 +21,7 @@ two_slice_marginals, viterbi, ) +from sofic.generators.matrices import emission_tensors from sofic.generators.mealy import MealyHMM from sofic.graph import ATTR_EMISSION, ATTR_PROB @@ -196,7 +196,7 @@ def test_baum_welch_accepts_single_sequence(): def test_score_matches_finite_difference_gradient(): gm = golden_mean(0.4) obs = [0, 1, 0, 0, 1, 0, 1, 0] - pi, joint = _emission_transition_tensors(gm) + pi, joint = emission_tensors(gm) idx = gm.to_mealy().reindex() a, b = idx.index("A"), idx.index("B") @@ -217,7 +217,7 @@ def loglik_entry(symbol: int, i: int, j: int, value: float) -> float: def test_observed_information_matches_numeric_hessian_scalar(): gm = golden_mean(0.5) obs = [0, 1, 0, 0, 1, 0, 1, 0] - pi, joint = _emission_transition_tensors(gm) + pi, joint = emission_tensors(gm) idx = gm.to_mealy().reindex() a, b = idx.index("A"), idx.index("B") @@ -246,7 +246,7 @@ def test_observed_information_multi_parameter_symmetric_and_matches_hessian(): labels = free_parameter_labels(hmm) assert labels == [("A", 0, "A"), ("A", 1, "A")] - pi, joint = _emission_transition_tensors(hmm) + pi, joint = emission_tensors(hmm) def loglik_free(theta: np.ndarray) -> float: perturbed = {sym: matrix.copy() for sym, matrix in joint.items()} diff --git a/tests/test_stack_inference.py b/tests/test_stack_inference.py index ce3d78e..d2290e6 100644 --- a/tests/test_stack_inference.py +++ b/tests/test_stack_inference.py @@ -261,3 +261,28 @@ def test_stack_cssr_exact_and_bonferroni_run(): assert oracle.word_probability(tuple(observations[:12])) > 0.0 with pytest.raises(ValueError, match="unknown correction"): stack_cssr(observations, alphabet=alphabet, Lmax=2, correction="holm") + + +@pytest.mark.parametrize("seed", [1, 5]) +def test_stack_cssr_keeps_histories_when_splitting(seed): + """Regression: determinization dropped histories that never emitted the splitting + symbol, leaving states whose successors were all unplaced (zero outgoing mass).""" + from sofic.shifts.sofic_dyck import transition_ref + + shift = motzkin_shift() + rng = np.random.default_rng(seed) + refs = [transition_ref(transition) for transition in shift.transitions()] + weights = rng.dirichlet(np.ones(len(refs))) + oracle = HiddenMarkovStackModel.from_sofic_dyck_shift(shift, dict(zip(refs, map(float, weights), strict=True))) + observations, _ = oracle.sample(4000, rng=rng) + alphabet = DyckAlphabet( + call_alphabet=shift.call_alphabet, + return_alphabet=shift.return_alphabet, + internal_alphabet=shift.internal_alphabet, + ) + inferred = stack_cssr(observations, alphabet=alphabet, Lmax=3, max_stack_depth=4, alpha=0.01) + inferred.validate() + held_out, _ = oracle.sample(400, rng=np.random.default_rng(100 + seed)) + windows = [tuple(held_out[i : i + 6]) for i in range(0, 390, 6)] + supported = sum(inferred.word_probability(window) > 0.0 for window in windows) + assert supported >= 0.9 * len(windows) diff --git a/tests/test_tikz.py b/tests/test_tikz.py index 422a990..34cf7d3 100644 --- a/tests/test_tikz.py +++ b/tests/test_tikz.py @@ -61,6 +61,37 @@ def test_sofic_dyck_tikz_marks_matched_edges(): assert tikz.count("m_{1}") == 2 +def test_vpa_tikz_separates_kind_from_mid(): + from sofic.automata.vpa import VisiblyPushdownAutomaton + + vpa = VisiblyPushdownAutomaton(call_alphabet=frozenset({"a"}), stack_alphabet=frozenset({"A"})) + vpa.graph.add_state("q") + vpa.add_call_transition("q", "q", "a", "A") + + tikz = model_to_tikz(vpa) + assert r"$a\mid call\mid \uparrow A$" in tikz + assert r"\midcall" not in tikz + + +def test_tmc_tikz_multiplicity_is_math_times(): + from sofic.graph import ATTR_MULTIPLICITY, ATTR_SYMBOL + from sofic.shifts.tmc import TopologicalMarkovChain + + tmc = TopologicalMarkovChain(symbol_alphabet=frozenset({"a"})) + tmc.graph.add_state("s") + tmc.graph.add_transition("s", "s", **{ATTR_SYMBOL: "a", ATTR_MULTIPLICITY: 2}) + + assert r"$a\mid\times 2$" in model_to_tikz(tmc) + + +def test_symbolic_prob_latex_is_not_escaped(): + sp = pytest.importorskip("sympy") + from sofic.viz._tikz_format import format_prob_latex + + a = sp.Symbol("a", positive=True) + assert format_prob_latex(a / (a + 1)) == r"\frac{a}{a + 1}" + + def test_bidirectional_tikz_uses_edge_labels(): tikz = model_to_tikz(golden_mean_bidirectional(0.5), style="paper") assert r"\Edge{" in tikz