Correctness fixes: bits everywhere, exact Fischer/Krieger covers, and other bugs (1/5)

docs/inference/model_selection.rst

@@ -8,9 +8,11 @@ Model Selection
 Classical model-selection criteria for choosing the order (state count) of a
 fitted generator when a fully-Bayesian evidence is unavailable or undesirable.
 The module scores any fitted :class:`~sofic.generators.base.HiddenMarkovModel`
-using the natural-log likelihood from
+using the log-likelihood (in bits) from
 :func:`sofic.generators.hmm_inference.log_likelihood` and a free-parameter count
-read off the transition graph:
+read off the transition graph. Log-likelihoods, log scores, and the MDL code
+length are in bits; AIC, AICc, BIC, and WAIC are computed from the natural
+log-likelihood so they keep their standard deviance scale:
 
 * **AIC** :cite:`Akaike1974` and the small-sample-corrected **AICc**
   :cite:`HurvichTsai1989`,

docs/quickstart.rst.txt

@@ -77,7 +77,7 @@ topological entropy (the log of the golden ratio) :cite:`LindMarcus1995`:
 
    @doctest float
    In [4]: tmc.topological_entropy()
-   Out[4]: 0.4812118250596035
+   Out[4]: 0.6942419136306174
 
 The Parry measure turns the shift into the maximum-entropy ``MealyHMM`` on it
 (``tmc.parry_measure()``), the measure that attains this topological entropy.

docs/shifts/covers.rst

@@ -5,18 +5,36 @@
 Covers
 ******
 
-Fischer and Krieger covers convert a Sofic shift into unifilar presentations
-:cite:`Fischer1975,Krieger1984,LindMarcus1995`.
-
-* :class:`LeftFischerCover`, :class:`RightFischerCover` — implemented.
-* :class:`LeftKriegerCover`, :class:`RightKriegerCover` — construction raises
-  :exc:`NotImplementedError`.
+Fischer and Krieger covers are canonical resolving presentations of a sofic
+shift :cite:`Fischer1975,Krieger1984,LindMarcus1995`. Following Lind and Marcus,
+*right* means right-resolving (deterministic reading forward):
+
+* :class:`RightFischerCover` -- the unique minimal right-resolving presentation
+  of an irreducible sofic shift. Reducible shifts raise
+  :exc:`~sofic.exceptions.SoficValidationError`, since their minimal
+  right-resolving presentation need not be unique.
+* :class:`RightKriegerCover` -- the future cover, with one vertex per follower
+  set of a left-infinite ray. Defined for every sofic shift; for an irreducible
+  shift the Fischer cover is its unique terminal component.
+* :class:`LeftFischerCover`, :class:`LeftKriegerCover` -- the left-resolving
+  mirror images, built from the reversed shift.
+
+Both constructions are exact: the Fischer cover is the terminal component of the
+follower-merged subset construction, and the Krieger cover's vertices are the
+images of path relations lying on cycles of the finite relation monoid.
 
 .. ipython::
 
-   In [1]: from sofic.examples import golden_mean_shift_parry; parry = golden_mean_shift_parry()
+   In [1]: from sofic.graph import ATTR_SYMBOL
+
+   In [2]: from sofic.shifts import RightKriegerCover, SoficShift
+
+   In [3]: even = SoficShift(symbol_alphabet=frozenset("01"))
 
-   In [2]: parry.validate()
+   In [4]: for s, t, a in [("A", "A", "0"), ("A", "B", "1"), ("B", "A", "1")]:
+      ...:     even.graph.add_transition(s, t, **{ATTR_SYMBOL: a})
+
+   In [5]: len(list(RightKriegerCover.from_sofic(even).states()))
 
 API
 ===
@@ -26,7 +44,11 @@ API
 .. autoclass:: RightFischerCover
    :members: from_sofic
 .. autoclass:: LeftKriegerCover
+   :members: from_sofic
 .. autoclass:: RightKriegerCover
+   :members: from_sofic
 
 .. autofunction:: sofic.shifts.cover_construction.left_fischer_from_sofic
 .. autofunction:: sofic.shifts.cover_construction.right_fischer_from_sofic
+.. autofunction:: sofic.shifts.cover_construction.left_krieger_from_sofic
+.. autofunction:: sofic.shifts.cover_construction.right_krieger_from_sofic

docs/shifts/topological_anatomy.rst

@@ -67,8 +67,8 @@ entropy across both parts, while the sofic even shift is purely bound
    Out[10]: 0.5527864045001022
 
 The parts add up to :math:`h_\mathrm{top}`, which equals
-:meth:`~sofic.shifts.sofic.SoficShift.topological_entropy` divided by
-:math:`\ln 2` on a right-resolving presentation.
+:meth:`~sofic.shifts.sofic.SoficShift.topological_entropy` (in bits) on a
+right-resolving presentation.
 
 API
 ===

docs/shifts/topological_markov_chain.rst

@@ -15,7 +15,7 @@ adjacency matrix. Its Parry measure is the maximum-entropy stochastic generator
 
    @doctest float
    In [2]: tmc.topological_entropy()
-   Out[2]: 0.48121182505960347
+   Out[2]: 0.6942419136306174
 
 API
 ===

docs/shifts/wheeler.rst

@@ -43,10 +43,9 @@ excluded outright, and :func:`wheeler_cover` raises
 
 The cover is a :class:`~sofic.shifts.covers.WheelerCover`, a
 :class:`~sofic.shifts.sofic.SoficShift` subclass that sits beside the Fischer
-and Krieger covers. It does not fill the Krieger stubs in
-:mod:`sofic.shifts.cover_construction`: a Krieger cover's states are *all*
-sets of pasts closed under the follower relation, whereas a Wheeler cover
-carries only those that happen to be recency intervals.
+and Krieger covers. It is not a Krieger cover: a Krieger cover's states are
+*all* follower sets of left-infinite pasts, whereas a Wheeler cover carries
+only those that happen to be recency intervals.
 
 Once a shift has a Wheeler cover, :func:`wheeler_index_of_shift` gives
 ``O(|w| log |A|)`` factor-language membership in place of scanning

sofic/automata/algorithms.py

@@ -244,9 +244,15 @@ def minimize_hopcroft(dfa: DFA, *, alphabet: frozenset[Any] | None = None) -> DF
 def equivalent(
     aut1: LabeledAutomaton,
     aut2: LabeledAutomaton,
-    alphabet: frozenset[Any],
+    alphabet: frozenset[Any] | None = None,
 ) -> bool:
-    """Return whether two automata recognize the same language over ``alphabet``."""
+    """Return whether two automata recognize the same language.
+
+    The comparison runs over ``alphabet`` together with every symbol either
+    automaton declares or uses, so a too-small ``alphabet`` cannot hide a
+    difference on the omitted symbols.
+    """
+    alphabet = frozenset(alphabet or ()) | _transition_alphabet(aut1) | _transition_alphabet(aut2)
     d1 = minimize(_to_nfa(aut1), alphabet=alphabet, algorithm="hopcroft")
     d2 = minimize(_to_nfa(aut2), alphabet=alphabet, algorithm="hopcroft")
     return _isomorphic_minimal_dfa(d1, d2, alphabet)
@@ -281,6 +287,15 @@ def _effective_alphabet(aut: LabeledAutomaton) -> frozenset[Any]:
     return frozenset(symbols)
 
 
+def _transition_alphabet(aut: LabeledAutomaton) -> frozenset[Any]:
+    symbols = {symbol for symbol in aut.input_alphabet if symbol is not EPSILON}
+    for transition in aut.transitions():
+        symbol = transition.data.get(ATTR_SYMBOL)
+        if symbol is not None and symbol is not EPSILON:
+            symbols.add(symbol)
+    return frozenset(symbols)
+
+
 def _forward_reachable(aut: LabeledAutomaton) -> set[Hashable]:
     if not aut.initial_states:
         return set()

sofic/automata/buchi_simulation.py

@@ -9,7 +9,12 @@
 
 
 def accepts_lasso_buchi(ba: BuchiAutomaton, prefix: Sequence[Any], loop: Sequence[Any]) -> bool:
-    """Accept if repeating ``loop`` after ``prefix`` visits an accepting state infinitely often."""
+    """Accept if repeating ``loop`` after ``prefix`` visits an accepting state infinitely often.
+
+    ``loop`` must be non-empty: ``prefix loop^omega`` is an infinite word only then.
+    """
+    if not loop:
+        raise ValueError("an ultimately periodic omega-word needs a non-empty loop")
     post = ba._run_nfa(prefix)
     if not post:
         return False

sofic/examples/epsilon_machines.py

@@ -27,6 +27,7 @@
 
 import numpy as np
 
+from sofic.exceptions import SoficError
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.mealy import MealyHMM
 from sofic.graph import ATTR_EMISSION, ATTR_FUTURE_SYMBOL, ATTR_PROB, TransitionGraph
@@ -1004,7 +1005,7 @@ def tent_map_misiurewicz_bidirectional_fig8(a: Any | None = None):
     else:
         try:
             reverse = EpsilonMachine.from_hmm(reverse_raw)
-        except Exception:
+        except SoficError:
             from sofic.generators.epsilon_machine import _row_normalized_presentation
 
             reverse = _row_normalized_presentation(reverse_raw)

sofic/examples/processes.py

@@ -64,7 +64,7 @@ def _stationary_initial(
     transition = transition / row_sums[:, None]
     try:
         pi = stationary_distribution_from_transition(transition)
-    except Exception:
+    except (ValueError, np.linalg.LinAlgError):
         return _uniform_initial(states)
     return {state: float(pi[i]) for i, state in enumerate(states)}
 

sofic/generators/alternative_complexity.py

@@ -100,5 +100,5 @@ def _stationary_vector(transition: np.ndarray) -> np.ndarray | None:
 
     try:
         return stationary_distribution_from_transition(transition)
-    except Exception:
+    except (ValueError, np.linalg.LinAlgError):
         return None

sofic/generators/block_convergence.py

@@ -8,6 +8,7 @@
 
 import numpy as np
 
+from sofic.exceptions import SoficError
 from sofic.generators._word_measures import (
     _block_caekl,
     _block_coinformation,
@@ -203,7 +204,7 @@ def _anatomy_curves(
 def _exact_anatomy_scalars(machine: EpsilonMachine) -> dict[str, float] | None:
     try:
         bidir = machine.to_bidirectional()
-    except Exception:
+    except (SoficError, ValueError, NotImplementedError, np.linalg.LinAlgError):
         return None
     h_mu = float(bidir.entropy_rate())
     rho_mu = float(bidir.predicted_information())

sofic/generators/block_entropy.py

@@ -9,6 +9,8 @@
 
 import numpy as np
 
+from sofic.exceptions import SoficError
+
 if TYPE_CHECKING:
     from sofic.generators.epsilon_machine import EpsilonMachine
 
@@ -327,7 +329,10 @@ def block_entropy_estimates(
 
     entropy_asymptote = excess_entropy + h_mu_l
     crypticity_estimate = state_block_entropy - block_state_entropy
-    crypticity = float(crypticity_estimate[-1]) if crypticity_estimate.size else 0.0
+    if use_exact:
+        crypticity = statistical_complexity - excess_entropy
+    else:
+        crypticity = float(crypticity_estimate[-1]) if crypticity_estimate.size else 0.0
 
     cm = _cm_extension_curves(
         lengths,
@@ -453,7 +458,7 @@ def _excess_entropy(
 ) -> float:
     try:
         return float(machine.to_bidirectional().excess_entropy())
-    except Exception:
+    except (SoficError, ValueError, NotImplementedError, np.linalg.LinAlgError):
         return _block_entropy_excess_entropy(machine, entropy_rate=entropy_rate, block_entropy=block_entropy)
 
 
@@ -498,7 +503,7 @@ def _estimated_excess_entropy(machine: EpsilonMachine, estimate: np.ndarray, *,
     if use_exact:
         try:
             return float(machine.excess_entropy())
-        except Exception:
+        except (SoficError, ValueError, NotImplementedError, np.linalg.LinAlgError):
             pass
     return float(estimate[-1]) if estimate.size else 0.0
 

sofic/generators/epsilon_inference.py

@@ -832,7 +832,8 @@ def suggest_lmax(
     select_markov_order = getattr(dit.inference, "select_markov_order", None)
     if select_markov_order is None:  # pragma: no cover - depends on the installed dit
         raise ImportError("suggest_lmax requires a dit release with dit.inference.select_markov_order")
-    seq = [repr(symbol) for symbol in sequence]
+    codes: dict[Any, str] = {}
+    seq = [codes.setdefault(symbol, str(len(codes))) for symbol in sequence]
     if max_order is None:
         k = max(2, len(set(seq)))
         max_order = max(1, min(10, int(np.log(max(len(seq), 1) / 5) / np.log(k)) - 1))

sofic/generators/hmm_inference.py

@@ -8,6 +8,7 @@
 
 from __future__ import annotations
 
+import warnings
 from collections import defaultdict
 from collections.abc import Hashable, Iterable, Sequence
 from typing import Any
@@ -131,7 +132,7 @@ def _stationary_emission_tensors(
             return limited, joint
     try:
         pi = stationary_distribution_from_transition(transition)
-    except Exception:
+    except (np.linalg.LinAlgError, ValueError):
         pi = pi_initial
     return pi, joint
 
@@ -143,7 +144,7 @@ def _forward_scaled(
 ) -> tuple[np.ndarray, np.ndarray]:
     """Return per-step-normalized forward messages and log scaling factors.
 
-    ``alpha_hat[t]`` sums to one; ``log P(obs) = log_scales.sum()``. A ``-inf``
+    ``alpha_hat[t]`` sums to one; ``log2 P(obs) = log_scales.sum()``. A ``-inf``
     entry in ``log_scales`` marks an impossible step. Normalizing each step avoids
     the underflow that makes the raw forward product vanish for long sequences.
     """
@@ -155,7 +156,7 @@ def _forward_scaled(
         log_scales[0] = -np.inf
         return alpha_hat, log_scales
     alpha_hat[0] = pi / total0
-    log_scales[0] = float(np.log(total0))
+    log_scales[0] = float(np.log2(total0))
     for t, symbol in enumerate(obs):
         matrix = joint.get(symbol)
         if matrix is None:
@@ -167,7 +168,7 @@ def _forward_scaled(
             log_scales[t + 1] = -np.inf
             continue
         alpha_hat[t + 1] = row / scale
-        log_scales[t + 1] = float(np.log(scale))
+        log_scales[t + 1] = float(np.log2(scale))
     return alpha_hat, log_scales
 
 
@@ -239,7 +240,7 @@ def _backward_scaled(joint: dict[Any, np.ndarray], obs: list[Any], n_states: int
 
 
 def log_likelihood(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> float:
-    """Natural-log likelihood ``log P(observations)``.
+    """Log-likelihood ``log2 P(observations)`` in bits.
 
     Uses the per-step-scaled forward recursion so the result stays finite for long
     sequences instead of underflowing to ``-inf``.
@@ -316,7 +317,7 @@ def _expected_edge_counts(
     - ``source_totals[i] = \sum_{t=0}^{n-1} P(X_t = i \mid Y)`` is the expected
       number of transitions out of state ``i`` (the Baum-Welch denominator);
     - ``gamma0`` is the smoothed marginal of the initial state ``X_0``;
-    - ``loglik`` is the natural-log likelihood of the sequence.
+    - ``loglik`` is the log-likelihood of the sequence in bits.
 
     Only edges present in ``joint`` (structural support) receive mass, so the
     statistics preserve the model topology.
@@ -391,7 +392,7 @@ def baum_welch(
     so the fit is returned as a plain :class:`~sofic.generators.mealy.MealyHMM`.
 
     Returns ``(fitted_model, loglik_trace)`` where ``loglik_trace`` is the
-    non-decreasing sequence of total natural-log likelihoods observed before each
+    non-decreasing sequence of total log-likelihoods (bits) observed before each
     parameter update.
 
     EM converges to a local maximum of the likelihood. With ``n_restarts > 1`` the
@@ -483,15 +484,25 @@ def _baum_welch_run(
         total_source = np.zeros(n_states, dtype=float)
         gamma0_sum = np.zeros(n_states, dtype=float)
         total_ll = 0.0
+        skipped = 0
         for obs in seqs:
             edge_counts, source_totals, gamma0, loglik = _expected_edge_counts(pi, joint, obs)
             if not np.isfinite(loglik):
+                skipped += 1
                 continue
             for key, value in edge_counts.items():
                 total_edge_counts[key] += value
             total_source += source_totals
             gamma0_sum += gamma0
             total_ll += loglik
+        if seqs and skipped == len(seqs):
+            raise ValueError("every observation sequence has zero probability under the model")
+        if skipped and _iteration == 0:
+            warnings.warn(
+                f"{skipped} of {len(seqs)} sequences have zero probability under the model and are ignored",
+                RuntimeWarning,
+                stacklevel=3,
+            )
         loglik_trace.append(total_ll)
         if prev_ll is not None and abs(total_ll - prev_ll) < tol:
             break
@@ -522,6 +533,10 @@ def score(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> dict[tuple[Has
     ``\nabla \log L(\theta) = E[\nabla \log f(X, Y; \theta) \mid Y]`` (Cappe,
     Moulines & Ryden, 2005, Section 10.2.3), evaluated in the raw (unconstrained)
     joint-edge parameters. Keys are ``(source, symbol, target)`` state labels.
+
+    The score is the gradient of the *natural* log-likelihood
+    :math:`\ln P(Y) = \ln 2 \cdot` :func:`log_likelihood`, the convention under
+    which the observed information and standard errors are standard.
     """
     mealy = hmm.to_mealy()
     idx = mealy.reindex()
@@ -761,7 +776,10 @@ def sample(
     mealy = _as_mealy_hmm(hmm)
     idx = mealy.reindex()
     pi, joint = _emission_transition_tensors_from_mealy(mealy)
-    state = int(generator.choice(len(idx), p=pi / pi.sum()))
+    total = float(np.sum(pi))
+    if not total > 0.0:
+        raise ValueError("cannot sample: the initial state distribution has no mass")
+    state = int(generator.choice(len(idx), p=pi / total))
 
     observations: list[Any] = []
     states: list[Hashable] = []