CSSR and transCSSR learn from independent Trials

CHANGELOG.md

@@ -115,6 +115,11 @@ brute-force reference or a closed-form value and has a regression test.
 
 ### New
 
+- `sofic.inference.Trials`: `learn_epsilon_machine_cssr`,
+  `learn_epsilon_transducer_cssr`, and `suggest_max_history` learn from
+  independent trials, pooling counts without any history spanning two trials
+  (`SuffixCounts.from_sequences`; `JointSuffixCounts.from_sequences` accepts
+  paired `Trials`). A `dit.inference.Trials` is accepted too.
 - Hypothesis strategies in `sofic.testing` for NFAs, Büchi automata and lassos,
   Wheeler NFAs, Markov chains, Mealy HMMs, sofic shifts, SFTs, VPAs, NWAs, and
   Mealy transducers; `ci` and `nightly` Hypothesis profiles (`HYPOTHESIS_PROFILE`).

docs/inference/cssr.rst

@@ -94,10 +94,24 @@ apply directly.
        observations, method="cssr", max_history="auto", test="exact", correction="bonferroni"
    )
 
+Independent recordings of the same process go in a :class:`Trials`. Counts are
+pooled across trials, but no suffix spans two trials, so the boundaries add no
+spurious histories. :func:`learn_epsilon_transducer_cssr` and
+:func:`suggest_max_history` accept :class:`Trials` too, as does a
+``dit.inference.Trials``.
+
+.. code-block:: python
+
+   from sofic.inference import Trials, learn_epsilon_machine_cssr
+
+   machine = learn_epsilon_machine_cssr(Trials([run_0, run_1, run_2]), max_history=3)
+
 .. autofunction:: learn_epsilon_machine_cssr
 
 .. autofunction:: suggest_max_history
 
+.. autoclass:: Trials
+
 After reconstruction, check the result with
 :func:`~sofic.inference.diagnostics.goodness_of_fit` and
 :func:`~sofic.inference.diagnostics.structure_stability`

sofic/inference/__init__.py

@@ -9,6 +9,7 @@
 
 from sofic.inference import bayesian, cssr, diagnostics, hmm, model_selection, spectral
 from sofic.inference.cssr import (
+    Trials,
     learn_epsilon_machine_cssr,
     learn_epsilon_machine_subtree,
     learn_epsilon_transducer_cssr,
@@ -69,6 +70,7 @@
     "hmm",
     "model_selection",
     "spectral",
+    "Trials",
     "learn_epsilon_machine_cssr",
     "GoodnessOfFit",
     "StructureStability",

sofic/inference/cssr/__init__.py

@@ -13,6 +13,7 @@
     JointSuffixCounts,
     StackSuffixCounts,
     SuffixCounts,
+    Trials,
 )
 from sofic.inference.cssr.process import learn_epsilon_machine_cssr, suggest_max_history
 from sofic.inference.cssr.significance import (
@@ -40,6 +41,7 @@
     "StackSuffixCounts",
     "SuffixCounts",
     "TableTest",
+    "Trials",
     "aggregates_differ",
     "learn_epsilon_machine_cssr",
     "learn_stack_hmm_mle",

sofic/inference/cssr/counts.py

@@ -22,7 +22,35 @@
 ConfigurationHistory = tuple[tuple[Any, ...], tuple[Any, ...]]
 
 
-def infer_alphabet(tokens: Sequence[Any], alphabet: Sequence[Any] | None) -> tuple[Any, ...]:
+class Trials(list):
+    """Independent realizations (trials) of the same process.
+
+    Wrap a list of sequences in :class:`Trials` to tell the CSSR learners that they
+    are separate recordings: counts are pooled across trials, but no history ever
+    spans the boundary between two trials. Trials may have different lengths. A
+    ``dit.inference.Trials`` is accepted wherever this class is.
+    """
+
+
+def is_trials(data: Any) -> bool:
+    """Whether ``data`` is a :class:`Trials` (or a ``dit.inference.Trials``)."""
+    if isinstance(data, Trials):
+        return True
+    try:
+        import dit.inference
+    except ImportError:  # pragma: no cover - dit is a dependency
+        return False
+    # Older dit releases have no Trials.
+    dit_trials = getattr(dit.inference, "Trials", None)
+    return dit_trials is not None and isinstance(data, dit_trials)
+
+
+def split_trials(data: Any) -> list[tuple[Any, ...]]:
+    """The trials of ``data`` as tuples: each trial of a :class:`Trials`, else ``data`` as one trial."""
+    return [tuple(trial) for trial in data] if is_trials(data) else [tuple(data)]
+
+
+def infer_alphabet(tokens: Iterable[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)
 
@@ -56,17 +84,32 @@ def from_sequence(
         alphabet: Sequence[Any] | None = None,
         max_length: int | None = None,
     ) -> SuffixCounts:
-        seq = tuple(sequence)
-        if not seq:
+        if not is_trials(sequence) and not tuple(sequence):
             raise ValueError("sequence must be non-empty")
-        alphabet = infer_alphabet(seq, alphabet)
-        unknown = set(seq) - set(alphabet)
+        return cls.from_sequences(split_trials(sequence), alphabet=alphabet, max_length=max_length)
+
+    @classmethod
+    def from_sequences(
+        cls,
+        sequences: Iterable[Sequence[Any]],
+        *,
+        alphabet: Sequence[Any] | None = None,
+        max_length: int | None = None,
+    ) -> SuffixCounts:
+        """Counts pooled over independent ``sequences``; no history spans two of them."""
+        trials = [tuple(seq) for seq in sequences]
+        symbols = {symbol for seq in trials for symbol in seq}
+        if not symbols:
+            raise ValueError("sequences must be non-empty")
+        alphabet = infer_alphabet(symbols, alphabet)
+        unknown = symbols - set(alphabet)
         if unknown:
             raise ValueError(f"symbols {unknown!r} not in alphabet")
         counts = cls(alphabet=alphabet)
-        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
+        for seq in trials:
+            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]:
@@ -138,18 +181,18 @@ def from_sequences(
         output_alphabet: Sequence[Any] | None = None,
         max_length: int,
     ) -> JointSuffixCounts:
-        xs = tuple(inputs)
-        ys = tuple(outputs)
-        if len(xs) != len(ys):
-            raise ValueError("inputs and outputs must have equal length")
-        if not xs:
+        """Counts of ``(inputs, outputs)``; both may be :class:`Trials` of matching lengths."""
+        x_trials, y_trials = split_paired_trials(inputs, outputs)
+        if not any(x_trials):
             raise ValueError("sequences must be non-empty")
         counts = cls(
-            input_alphabet=infer_alphabet(xs, input_alphabet), output_alphabet=infer_alphabet(ys, output_alphabet)
+            input_alphabet=infer_alphabet((x for xs in x_trials for x in xs), input_alphabet),
+            output_alphabet=infer_alphabet((y for ys in y_trials for y in 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
+        for xs, ys in zip(x_trials, y_trials, strict=True):
+            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]:
@@ -172,6 +215,16 @@ 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 split_paired_trials(inputs: Any, outputs: Any) -> tuple[list[tuple[Any, ...]], list[tuple[Any, ...]]]:
+    """The trials of paired ``inputs`` / ``outputs``: both :class:`Trials`, or both single sequences."""
+    if is_trials(inputs) != is_trials(outputs):
+        raise ValueError("inputs and outputs must both be Trials, or both be single sequences")
+    x_trials, y_trials = split_trials(inputs), split_trials(outputs)
+    if len(x_trials) != len(y_trials) or any(len(xs) != len(ys) for xs, ys in zip(x_trials, y_trials, strict=True)):
+        raise ValueError("inputs and outputs must have equal length")
+    return x_trials, y_trials
+
+
 #: Aggregated output counts of a state: ``agg[input_symbol]`` is a Counter over outputs.
 StateAggregate = dict[Any, Counter[Any]]
 

sofic/inference/cssr/process.py

@@ -14,7 +14,7 @@
 from sofic.exceptions import StochasticValidationError
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.graph import ATTR_EMISSION, ATTR_PROB, TransitionGraph
-from sofic.inference.cssr.counts import History, SuffixCounts
+from sofic.inference.cssr.counts import History, SuffixCounts, is_trials, split_trials
 from sofic.inference.cssr.significance import (
     MorphTest,
     _bonferroni_alpha,
@@ -197,7 +197,7 @@ def _suffix_determinize(
 def _suffix_machine(
     states: list[set[History]],
     counts: SuffixCounts,
-    sequence: Sequence[Any],
+    trials: Sequence[tuple[Any, ...]],
     alive: set[int],
     *,
     max_history: int,
@@ -210,13 +210,13 @@ def _suffix_machine(
     history_to_state = {h: index for index in alive for h in states[index]}
 
     visits: Counter[int] = Counter()
-    seq = tuple(sequence)
-    for t in range(len(seq) + 1):
-        for length in range(min(t, max_history), -1, -1):
-            state = history_to_state.get(seq[t - length : t])
-            if state is not None:
-                visits[state] += 1
-                break
+    for seq in trials:
+        for t in range(len(seq) + 1):
+            for length in range(min(t, max_history), -1, -1):
+                state = history_to_state.get(seq[t - length : t])
+                if state is not None:
+                    visits[state] += 1
+                    break
 
     classes = _recurrent_states(edges)
     if not classes:
@@ -265,13 +265,16 @@ def _suffix_machine(
 def _suffix_reconstruct(
     states: list[set[History]],
     counts: SuffixCounts,
-    sequence: Sequence[Any],
+    trials: Sequence[tuple[Any, ...]],
     *,
     max_history: int,
     alpha: float,
     test: MorphTest,
 ) -> EpsilonMachine:
-    """Prune transient states, determinize, and build the machine from homogeneous ``states``."""
+    """Prune transient states, determinize, and build the machine from homogeneous ``states``.
+
+    ``trials`` are the observed sequences, one tuple per independent trial.
+    """
 
     def reconstruct(resolve: bool) -> EpsilonMachine:
         everything = set(range(len(states)))
@@ -283,13 +286,13 @@ def reconstruct(resolve: bool) -> EpsilonMachine:
             states, counts, alive, max_history=max_history, alpha=alpha, test=test, resolve=resolve
         )
         return _suffix_machine(
-            split, counts, sequence, alive, max_history=max_history, alpha=alpha, test=test, resolve=resolve
+            split, counts, trials, alive, max_history=max_history, alpha=alpha, test=test, resolve=resolve
         )
 
     machine = reconstruct(resolve=True)
     # Resolving truncated successors needs max_history at least the synchronization length. When it
     # is shorter, resolution can close off a state that never emits some observed symbol.
-    if {t.data[ATTR_EMISSION] for t in machine.transitions()} < set(sequence):
+    if {t.data[ATTR_EMISSION] for t in machine.transitions()} < {symbol for seq in trials for symbol in seq}:
         machine = reconstruct(resolve=False)
     return machine
 
@@ -314,7 +317,7 @@ def suggest_max_history(
     Parameters
     ----------
     sequence
-        Observed symbols.
+        Observed symbols, or independent :class:`~sofic.inference.cssr.Trials` of them.
     alpha
         Significance level of each order test.
     max_order
@@ -341,11 +344,19 @@ def suggest_max_history(
     if select_markov_order is None:  # pragma: no cover - depends on the installed dit
         raise ImportError("suggest_max_history requires a dit release with dit.inference.select_markov_order")
     codes: dict[Any, str] = {}
-    seq = [codes.setdefault(symbol, str(len(codes))) for symbol in sequence]
+    trials = [[codes.setdefault(symbol, str(len(codes))) for symbol in seq] for seq in split_trials(sequence)]
+    n = sum(len(seq) for seq in trials)
     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))
-    order = select_markov_order(seq, max_order, method=method, alpha=alpha, n_surrogates=n_surrogates, prng=seed)
+        k = max(2, len(codes))
+        max_order = max(1, min(10, int(np.log(max(n, 1) / 5) / np.log(k)) - 1))
+    if is_trials(sequence):
+        DitTrials = getattr(dit.inference, "Trials", None)
+        if DitTrials is None:  # pragma: no cover - depends on the installed dit
+            raise ImportError("suggest_max_history on Trials requires a dit release with dit.inference.Trials")
+        data: Any = DitTrials(trials)
+    else:
+        data = trials[0]
+    order = select_markov_order(data, max_order, method=method, alpha=alpha, n_surrogates=n_surrogates, prng=seed)
     return max(1, int(order))
 
 
@@ -369,11 +380,12 @@ def learn_epsilon_machine_cssr(
     Parameters
     ----------
     sequence
-        Observed symbols.
+        Observed symbols, or independent :class:`~sofic.inference.cssr.Trials` of them:
+        counts are pooled across trials, but no suffix spans two trials.
     alphabet
         Symbol alphabet; defaults to the symbols in ``sequence``.
     max_history
-        Longest suffix considered; by default a third of ``log_k len(sequence)``,
+        Longest suffix considered; by default a third of ``log_k n`` for ``n`` symbols in total,
         between 1 and 10. It should be at least the synchronization length of the
         source (for a Markov source, its order). Larger values run many more
         significance tests, and some of them split states by chance. ``"auto"``
@@ -403,21 +415,22 @@ def learn_epsilon_machine_cssr(
     reconstruction. CSSR then returns more states than the ε-machine, with an
     entropy rate that approaches the true one from above as ``max_history`` grows.
     """
-    seq = tuple(sequence)
-    if len(seq) < 2:
+    trials = split_trials(sequence)
+    n = sum(len(seq) for seq in trials)
+    if n < 2:
         raise ValueError("sequence must contain at least two symbols")
-    alphabet_size = len(set(seq)) if alphabet is None else len(tuple(alphabet))
+    alphabet_size = len({symbol for seq in trials for symbol in seq}) if alphabet is None else len(tuple(alphabet))
     if max_history == "auto":
-        max_length = suggest_max_history(seq, alpha=alpha)
+        max_length = suggest_max_history(sequence if is_trials(sequence) else trials[0], alpha=alpha)
     else:
-        max_length = max_history if max_history is not None else _default_max_history(len(seq), alphabet_size)
+        max_length = max_history if max_history is not None else _default_max_history(n, alphabet_size)
     if max_length < 0:
         raise ValueError("max_history must be non-negative")
-    counts = SuffixCounts.from_sequence(seq, alphabet=alphabet, max_length=max_length + 1)
+    counts = SuffixCounts.from_sequences(trials, alphabet=alphabet, max_length=max_length + 1)
     if correction == "bonferroni":
         alpha = _bonferroni_alpha(counts, alpha, max_length=max_length, min_count=min_count, resolution_tests=True)
     elif correction is not None:
         raise ValueError(f"unknown correction {correction!r}")
 
     homogeneous = _suffix_homogenize(counts, max_history=max_length, alpha=alpha, test=test, min_count=min_count)
-    return _suffix_reconstruct(homogeneous, counts, seq, max_history=max_length, alpha=alpha, test=test)
+    return _suffix_reconstruct(homogeneous, counts, trials, max_history=max_length, alpha=alpha, test=test)

sofic/inference/cssr/subtree.py

@@ -122,4 +122,4 @@ def learn_epsilon_machine_subtree(
         raise ValueError(f"unknown correction {correction!r}")
     states = list(_cluster_histories_by_morph(counts, histories, delta=delta, alpha=alpha, test=test).values())
 
-    return _suffix_reconstruct(states, counts, seq, max_history=max_history, alpha=alpha, test=test)
+    return _suffix_reconstruct(states, counts, [seq], max_history=max_history, alpha=alpha, test=test)