From 660fc5318392485f86f699539f51ef3e9d6c6ff7 Mon Sep 17 00:00:00 2001
From: Ryan James <ryan.james@reddit.com>
Date: Fri, 13 Nov 2026 14:22:40 -0700
Subject: [PATCH 1/6] Anytime pipeline on large joints: LP time limits, dit
 gate, memory-lift deadline

Three stages ran far past their deadlines on 4-state survey joints with
95-144 causal-state pairs:
- generator caps: one degenerate node LP ran 8+ minutes; child LPs now get
  the time left as a HiGHS time_limit (a cut-off LP keeps its parent's bound);
- dit exact/Wyner fallbacks take no deadline; skipped above 30 pairs (every
  3-state survey joint has at most 21);
- memory_lift_search ran its whole first start regardless of the deadline;
  a start in progress now stops at it.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md                           | 16 ++++++++++++++++
 sofic/generators/coupled_generators.py | 14 ++++++++++++--
 sofic/generators/generative_search.py  | 10 ++++++++--
 sofic/generators/generator_caps.py     | 16 +++++++++++++---
 tests/test_coupled_generators.py       | 21 +++++++++++++++++++++
 tests/test_generative_model.py         | 24 ++++++++++++++++++++++++
 tests/test_generator_caps.py           | 26 ++++++++++++++++++++++++++
 7 files changed, 120 insertions(+), 7 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a85ade8..9b8aaf3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,6 +20,11 @@ brute-force reference or a closed-form value and has a regression test.
 
 ### Generator caps
 
+- With a deadline, each branch-and-bound child LP of `generator_weight_cap` gets the time left as a
+  HiGHS `time_limit`; a cut-off LP keeps its parent's bound, so the cap stays valid. On a 4-state
+  survey process (4 × 36 causal states) one degenerate node LP had run more than 8 minutes past a
+  90 s stage deadline; the stage now ends on time (92.7 s).
+
 - `generator_weight_cap` adds the forward analogue of the reverse joint-moment rows: per-state forward
   stationarity times the state's own reverse belief gives
   `T_h(a', b) = sum_{x, g, a : f(a, x) = a'} Gt[x, g, h, a, b]`, which every generator satisfies.
@@ -27,6 +32,17 @@ brute-force reference or a closed-form value and has a regression test.
   survey processes `((1,1),(2,0),(1,0))` and `((1,·),(2,0),(1,0))` rises 1.261 -> 1.277 and
   1.175 -> 1.184.
 
+### Anytime pipeline on large joints
+
+- `generative_model` skips the dit exact- and Wyner-common-information fallbacks when the joint has
+  more than 30 causal-state pairs (logged as a skipped `dit fallbacks` stage). dit's optimizers take no
+  deadline: on finished 4-state survey cases with 30-70 pairs they used 190-760 s of a 900 s budget,
+  and on a 95-pair case one call ran past the survey's hard timeout. Every 3-state survey joint has at
+  most 21 pairs, so the benchmark is unaffected.
+- `memory_lift_search` stops a start in progress at its `deadline` (keeping the last finished penalty
+  stage) instead of only refusing to begin new starts. Its first start used to run all four
+  finite-difference L-BFGS-B stages whatever the deadline, more than 14 minutes on a 95-pair joint.
+
 ### Support-type bound speed
 
 - Exact vertex enumeration of the support-type polytope (`support_type_bound`, the certified flow
diff --git a/sofic/generators/coupled_generators.py b/sofic/generators/coupled_generators.py
index 6978be4..aacf7e0 100644
--- a/sofic/generators/coupled_generators.py
+++ b/sofic/generators/coupled_generators.py
@@ -407,6 +407,10 @@ def memory_lift(bidir: BidirectionalEpsilonMachine | _GeneratorProblem, nu: np.n
     return pi, Q, R, K
 
 
+class _DeadlineReached(Exception):
+    """Raised inside an objective to stop an optimizer at a deadline."""
+
+
 def memory_lift_search(
     bidir: BidirectionalEpsilonMachine,
     m: int,
@@ -420,7 +424,8 @@ def memory_lift_search(
 
     A penalty continuation over the channel's logits (numerical gradients; meant for small problems).
     Returns ``(entropy, residual, nu)`` for every start, best first. No new start begins after
-    ``deadline`` (a :func:`time.monotonic` value).
+    ``deadline`` (a :func:`time.monotonic` value), and a start in progress stops at it, keeping the logits
+    of its last finished penalty stage.
     """
     rng = rng or np.random.default_rng()
     problem = _GeneratorProblem(bidir)
@@ -432,6 +437,8 @@ def channel(z: np.ndarray) -> np.ndarray:
         return z / z.sum(axis=1, keepdims=True)
 
     def value(z: np.ndarray, mu: float) -> float:
+        if deadline is not None and time.monotonic() > deadline:
+            raise _DeadlineReached
         pi, Q, R, K = memory_lift(problem, channel(z))
         p = pi[pi > 1e-15]
         penalty = sum(float((r**2).sum()) for r in problem.residuals(pi, Q, R, K))
@@ -443,7 +450,10 @@ def value(z: np.ndarray, mu: float) -> float:
             break
         z = rng.normal(scale=2.0, size=n * m)
         for mu in (1e2, 1e4, 1e6, 1e8):
-            z = minimize(value, z, args=(mu,), method="L-BFGS-B", options={"maxiter": maxiter}).x
+            try:
+                z = minimize(value, z, args=(mu,), method="L-BFGS-B", options={"maxiter": maxiter}).x
+            except _DeadlineReached:
+                break
         nu = channel(z)
         pi, Q, R, K = memory_lift(problem, nu)
         p = pi[pi > 1e-15]
diff --git a/sofic/generators/generative_search.py b/sofic/generators/generative_search.py
index c8cf38d..c7d4b05 100644
--- a/sofic/generators/generative_search.py
+++ b/sofic/generators/generative_search.py
@@ -75,6 +75,9 @@
 _TARGETED_GAP = 1e-2
 # ``I[S+ : S- | W]`` (bits) above which a channel is not treated as an exact decomposition.
 _SEPARATION_ATOL = 1e-9
+# Causal-state pairs above which the dit fallbacks are skipped: dit's optimizers take no deadline, and
+# on 4-state survey joints with 30-70 pairs they spent 190-760 s of a 900 s budget (95 pairs: over 10 min).
+_DIT_MAX_PAIRS = 30
 
 
 def _entropy(p: np.ndarray) -> float:
@@ -500,7 +503,10 @@ def maybe_return_early() -> GenerativeModelResult | None:
         best, optima = separating_optima(problem, starts=exact_starts, deadline=share(0.25), rng=rng)
         if np.isfinite(best):
             g_estimate = best
-        if dit_fallback and not out_of_time("dit exact common information"):
+        dit_small = len(lifted.pairs) <= _DIT_MAX_PAIRS
+        if dit_fallback and not dit_small:
+            log.append({"stage": "dit fallbacks", "skipped": f"more than {_DIT_MAX_PAIRS} causal-state pairs"})
+        if dit_fallback and dit_small and not out_of_time("dit exact common information"):
             try:
                 joint_n = _normalized_joint_distribution(bidir, cutoff=1e-10)
                 c = _unconstrained_channel(
@@ -526,7 +532,7 @@ def maybe_return_early() -> GenerativeModelResult | None:
                             optima.append(c)
             except (StochasticValidationError, ImportError):
                 pass
-        if wyner and dit_fallback and not out_of_time("dit Wyner common information"):
+        if wyner and dit_fallback and dit_small and not out_of_time("dit Wyner common information"):
             try:
                 c = _unconstrained_channel(
                     lifted,
diff --git a/sofic/generators/generator_caps.py b/sofic/generators/generator_caps.py
index 1b32824..287a665 100644
--- a/sofic/generators/generator_caps.py
+++ b/sofic/generators/generator_caps.py
@@ -174,8 +174,11 @@ def envelope(prod: int, factor: int, e: int, is_q: bool, i: int) -> None:
         self.c = np.zeros(self.n)
         self.c[m[:k]] = -1.0
 
-    def bound(self, boxes: tuple[tuple[np.ndarray, ...], ...]) -> float | None:
-        """Max total weight of the explicit states with beliefs in ``boxes``; ``None`` if the box is infeasible."""
+    def bound(self, boxes: tuple[tuple[np.ndarray, ...], ...], time_limit: float | None = None) -> float | None:
+        """Max total weight of the explicit states with beliefs in ``boxes``; ``None`` if the box is infeasible.
+
+        An LP stopped by ``time_limit`` (seconds) returns the trivial bound 1.
+        """
         res = linprog(
             self.c,
             A_ub=self.inequalities(boxes),
@@ -184,6 +187,7 @@ def bound(self, boxes: tuple[tuple[np.ndarray, ...], ...]) -> float | None:
             b_eq=self.b_eq,
             bounds=(0, 1),
             method="highs",
+            options=None if time_limit is None else {"time_limit": max(time_limit, 1e-3)},
         )
         if res.status == _HIGHS_INFEASIBLE:
             return None
@@ -271,9 +275,15 @@ def _branch_and_bound(
     Returns the open leaves ``(bound, boxes, path)`` and the ``(boxes, path)`` dropped because their LP is
     infeasible; a path lists the ``(state, part, coordinate, side)`` splits from the root. Boxes dropped
     by :func:`_tighten` miss the simplex, so the open and infeasible leaves cover every generator.
+    With a ``deadline`` each child LP is also limited to the time left, since a single degenerate LP can
+    otherwise run far past it; a cut-off LP keeps its parent's bound. The root LP always runs to the end.
     """
     k = relaxation.k
     root = _root(relaxation)
+
+    def left() -> float | None:
+        return None if deadline is None else deadline - time.monotonic()
+
     first = relaxation.bound(root)
     if first is None:
         return [], [(root, ())]
@@ -294,7 +304,7 @@ def _branch_and_bound(
             if child is None:
                 continue
             child_path = (*path, (e, part, i, side))
-            value = relaxation.bound(child)
+            value = relaxation.bound(child, left())
             if value is None:
                 infeasible.append((child, child_path))
             else:
diff --git a/tests/test_coupled_generators.py b/tests/test_coupled_generators.py
index 5509286..c6aa454 100644
--- a/tests/test_coupled_generators.py
+++ b/tests/test_coupled_generators.py
@@ -389,6 +389,27 @@ def test_separating_optima_and_bilevel_respect_an_expired_deadline() -> None:
     assert len(out) == 1
 
 
+def test_memory_lift_search_stops_inside_a_start_at_the_deadline(monkeypatch: pytest.MonkeyPatch) -> None:
+    """The first start used to run all four penalty stages to the end whatever the deadline."""
+    clock = {"now": 0.0}
+    calls = {"n": 0}
+    real = cg.memory_lift
+
+    def tick() -> float:
+        clock["now"] += 1.0
+        return clock["now"]
+
+    def counted(*args, **kwargs):
+        calls["n"] += 1
+        return real(*args, **kwargs)
+
+    monkeypatch.setattr(cg.time, "monotonic", tick)
+    monkeypatch.setattr(cg, "memory_lift", counted)
+    lifts = cg.memory_lift_search(_rll(), 2, starts=5, deadline=5.0, rng=np.random.default_rng(0))
+    assert len(lifts) == 1
+    assert calls["n"] <= 6
+
+
 def test_coupled_memory_and_targeted_searches_stop_at_an_expired_deadline() -> None:
     import time
 
diff --git a/tests/test_generative_model.py b/tests/test_generative_model.py
index cc4c218..f8b1029 100644
--- a/tests/test_generative_model.py
+++ b/tests/test_generative_model.py
@@ -68,3 +68,27 @@ def test_zero_time_budget_returns_a_valid_bracket_from_the_machines() -> None:
 def test_unlimited_budget_is_not_reported_as_exhausted() -> None:
     result = generative_model(golden_mean(0.5), rng=np.random.default_rng(0))
     assert not result.budget_exhausted
+
+
+def test_dit_fallbacks_are_skipped_on_large_joints(monkeypatch) -> None:
+    import sofic.generators.generative_search as gs
+    import sofic.generators.lifted_common_information as lci
+    from tests.test_generator_caps import _survey_machine
+
+    def fail(*args, **kwargs):
+        raise AssertionError("dit fallback ran")
+
+    monkeypatch.setattr(gs, "_DIT_MAX_PAIRS", 0)
+    monkeypatch.setattr(lci, "_unconstrained_channel", fail)
+    result = generative_model_from_bidirectional(
+        _survey_machine(),
+        stages=("machines", "exact"),
+        exact_starts=2,
+        coverage="off",
+        generator_caps=0,
+        early_certified=False,
+        time_budget=5.0,
+        rng=np.random.default_rng(0),
+    )
+    skipped = [s for s in result.stages if s["stage"] == "dit fallbacks"]
+    assert skipped and skipped[0]["skipped"] == "more than 0 causal-state pairs"
diff --git a/tests/test_generator_caps.py b/tests/test_generator_caps.py
index fe3d386..031b98a 100644
--- a/tests/test_generator_caps.py
+++ b/tests/test_generator_caps.py
@@ -119,3 +119,29 @@ def test_expired_deadline_still_gives_a_valid_cap() -> None:
     assert 0.5 - 1e-9 <= cap <= 1.0
     cov = generator_coverage_bound(rll, explicit=2, deadline=time.monotonic())
     assert cov.value <= 1.5 + 1e-6
+
+
+def test_child_lps_get_the_time_left_and_keep_the_parent_bound(monkeypatch: pytest.MonkeyPatch) -> None:
+    """With a deadline, a child LP that HiGHS stops at its time limit keeps the parent's bound."""
+    import time
+    from types import SimpleNamespace
+
+    import sofic.generators.generator_caps as gc
+
+    real = gc.linprog
+    seen: list[dict | None] = []
+
+    def fake(*args, options=None, **kwargs):
+        seen.append(options)
+        if options is not None:
+            return SimpleNamespace(status=1, success=False, fun=None)
+        return real(*args, **kwargs)
+
+    monkeypatch.setattr(gc, "linprog", fake)
+    relaxation = gc._HeaviestStatesRelaxation(cg._GeneratorProblem(golden_mean_bidirectional(0.5)), 1)
+    leaves, infeasible = gc._branch_and_bound(relaxation, 4, deadline=time.monotonic() + 60)
+    assert seen[0] is None
+    assert len(seen) > 1 and all(o is not None and 0 < o["time_limit"] <= 60 for o in seen[1:])
+    root = relaxation.bound(gc._root(relaxation))
+    assert not infeasible
+    assert max(b for b, _, _ in leaves) == pytest.approx(root)

From 867545123bb1e90a0fc3498f3897071d1e1e68f5 Mon Sep 17 00:00:00 2001
From: Ryan James <ryan.james@reddit.com>
Date: Fri, 13 Nov 2026 14:45:16 -0700
Subject: [PATCH 2/6] Line-split ALM polish stops at the pipeline deadline

_optimize and _optimize_coupled take a deadline; every L-BFGS-B objective
raises once it passes and the run returns its last finished step. The
pipeline's line-split stage passes its deadline and checks it between
floors (one polish ran past the budget for over 10 minutes on a 144-pair
4-state joint). Rename the private exception to _DeadlineReachedError (N818).

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 sofic/generators/coupled_generators.py | 106 ++++++++++++-------------
 sofic/generators/generative_search.py  |   4 +-
 tests/test_coupled_generators.py       |  16 ++++
 3 files changed, 69 insertions(+), 57 deletions(-)

diff --git a/sofic/generators/coupled_generators.py b/sofic/generators/coupled_generators.py
index aacf7e0..155318c 100644
--- a/sofic/generators/coupled_generators.py
+++ b/sofic/generators/coupled_generators.py
@@ -407,7 +407,7 @@ def memory_lift(bidir: BidirectionalEpsilonMachine | _GeneratorProblem, nu: np.n
     return pi, Q, R, K
 
 
-class _DeadlineReached(Exception):
+class _DeadlineReachedError(Exception):
     """Raised inside an objective to stop an optimizer at a deadline."""
 
 
@@ -438,7 +438,7 @@ def channel(z: np.ndarray) -> np.ndarray:
 
     def value(z: np.ndarray, mu: float) -> float:
         if deadline is not None and time.monotonic() > deadline:
-            raise _DeadlineReached
+            raise _DeadlineReachedError
         pi, Q, R, K = memory_lift(problem, channel(z))
         p = pi[pi > 1e-15]
         penalty = sum(float((r**2).sum()) for r in problem.residuals(pi, Q, R, K))
@@ -452,7 +452,7 @@ def value(z: np.ndarray, mu: float) -> float:
         for mu in (1e2, 1e4, 1e6, 1e8):
             try:
                 z = minimize(value, z, args=(mu,), method="L-BFGS-B", options={"maxiter": maxiter}).x
-            except _DeadlineReached:
+            except _DeadlineReachedError:
                 break
         nu = channel(z)
         pi, Q, R, K = memory_lift(problem, nu)
@@ -1007,8 +1007,12 @@ def _optimize(
     maxiter: int,
     *,
     alternating: bool = False,
+    deadline: float | None = None,
 ) -> np.ndarray:
-    """Feasibility pass, then block ALM; joint beliefs+kernel by default."""
+    """Feasibility pass, then block ALM; joint beliefs+kernel by default.
+
+    At ``deadline`` (a :func:`time.monotonic` value) the run stops and returns its last finished step.
+    """
     block_mus = (_ALM_MU, _ALM_MU, _ALM_MU)
     pi, Q, R, K = problem.unpack(theta, k)
     if problem.residual(pi, Q, R, K) <= _ALM_TOL:
@@ -1017,31 +1021,45 @@ def _optimize(
     def _block_max(res: tuple[np.ndarray, ...]) -> tuple[float, float, float]:
         return float(np.abs(res[0]).max()), float(np.abs(res[1]).max()), float(np.abs(res[2]).max())
 
+    def lbfgs(fun, x: np.ndarray, args: tuple) -> np.ndarray:
+        def timed(*a):
+            if deadline is not None and time.monotonic() > deadline:
+                raise _DeadlineReachedError
+            return fun(*a)
+
+        return minimize(timed, x, args=args, jac=True, method="L-BFGS-B", options={"maxiter": maxiter}).x
+
     if alternating:
         theta_b = problem.logits_beliefs(pi, Q, R)
         lam = [np.zeros_like(r) for r in problem.residuals(pi, Q, R, K)]
-        theta_b = minimize(
-            problem.loss_and_grad_beliefs,
-            theta_b,
-            args=(k, K, 1.0, None, 0.0, block_mus),
-            jac=True,
-            method="L-BFGS-B",
-            options={"maxiter": maxiter},
-        ).x
-        pi, Q, R = problem.unpack_beliefs(theta_b, k)
+        try:
+            theta_b = lbfgs(problem.loss_and_grad_beliefs, theta_b, (k, K, 1.0, None, 0.0, block_mus))
+            pi, Q, R = problem.unpack_beliefs(theta_b, k)
+            for _ in range(_ALM_ROUNDS):
+                _viol, K = _kernel_lp_beliefs(problem, pi, Q, R)
+                shifts = tuple(l / m for l, m in zip(lam, block_mus, strict=False))
+                theta_b = lbfgs(problem.loss_and_grad_beliefs, theta_b, (k, K, 1.0, shifts, 0.0, block_mus))
+                pi, Q, R = problem.unpack_beliefs(theta_b, k)
+                res = problem.residuals(pi, Q, R, K)
+                lam = [l + m * r for l, r, m in zip(lam, res, block_mus, strict=False)]
+                mE, mF, mB = _block_max(res)
+                if max(mE, mF, mB) < _ALM_TOL:
+                    break
+                block_mus = tuple(
+                    min(m * (2.0 if v > _ALM_TOL else 1.0), _ALM_MU_MAX)
+                    for m, v in zip(block_mus, (mE, mF, mB), strict=False)
+                )
+        except _DeadlineReachedError:
+            pass
+        return problem.logits(pi, Q, R, K)
+
+    try:
+        theta = lbfgs(problem.loss_and_grad, theta, (k, 1.0, None, 0.0, block_mus))
+        lam = [np.zeros_like(r) for r in problem.residuals(*problem.unpack(theta, k))]
         for _ in range(_ALM_ROUNDS):
-            _viol, K = _kernel_lp_beliefs(problem, pi, Q, R)
             shifts = tuple(l / m for l, m in zip(lam, block_mus, strict=False))
-            theta_b = minimize(
-                problem.loss_and_grad_beliefs,
-                theta_b,
-                args=(k, K, 1.0, shifts, 0.0, block_mus),
-                jac=True,
-                method="L-BFGS-B",
-                options={"maxiter": maxiter},
-            ).x
-            pi, Q, R = problem.unpack_beliefs(theta_b, k)
-            res = problem.residuals(pi, Q, R, K)
+            theta = lbfgs(problem.loss_and_grad, theta, (k, 1.0, shifts, 0.0, block_mus))
+            res = problem.residuals(*problem.unpack(theta, k))
             lam = [l + m * r for l, r, m in zip(lam, res, block_mus, strict=False)]
             mE, mF, mB = _block_max(res)
             if max(mE, mF, mB) < _ALM_TOL:
@@ -1050,35 +1068,8 @@ def _block_max(res: tuple[np.ndarray, ...]) -> tuple[float, float, float]:
                 min(m * (2.0 if v > _ALM_TOL else 1.0), _ALM_MU_MAX)
                 for m, v in zip(block_mus, (mE, mF, mB), strict=False)
             )
-        return problem.logits(pi, Q, R, K)
-
-    theta = minimize(
-        problem.loss_and_grad,
-        theta,
-        args=(k, 1.0, None, 0.0, block_mus),
-        jac=True,
-        method="L-BFGS-B",
-        options={"maxiter": maxiter},
-    ).x
-    lam = [np.zeros_like(r) for r in problem.residuals(*problem.unpack(theta, k))]
-    for _ in range(_ALM_ROUNDS):
-        shifts = tuple(l / m for l, m in zip(lam, block_mus, strict=False))
-        theta = minimize(
-            problem.loss_and_grad,
-            theta,
-            args=(k, 1.0, shifts, 0.0, block_mus),
-            jac=True,
-            method="L-BFGS-B",
-            options={"maxiter": maxiter},
-        ).x
-        res = problem.residuals(*problem.unpack(theta, k))
-        lam = [l + m * r for l, r, m in zip(lam, res, block_mus, strict=False)]
-        mE, mF, mB = _block_max(res)
-        if max(mE, mF, mB) < _ALM_TOL:
-            break
-        block_mus = tuple(
-            min(m * (2.0 if v > _ALM_TOL else 1.0), _ALM_MU_MAX) for m, v in zip(block_mus, (mE, mF, mB), strict=False)
-        )
+    except _DeadlineReachedError:
+        pass
     return theta
 
 
@@ -1091,18 +1082,21 @@ def _optimize_coupled(
     alternating: bool = False,
     dual_strategy: bool = True,
     residual_atol: float = 1e-6,
+    deadline: float | None = None,
 ) -> np.ndarray:
     """Run block ALM; optionally retry with alternating ``K`` LP if joint ALM stays infeasible."""
     pi, Q, R, K = problem.unpack(theta, k)
     if problem.residual(pi, Q, R, K) <= residual_atol:
         return theta
     if alternating:
-        return _optimize(problem, theta, k, maxiter, alternating=True)
-    theta_joint = _optimize(problem, theta, k, maxiter, alternating=False)
+        return _optimize(problem, theta, k, maxiter, alternating=True, deadline=deadline)
+    theta_joint = _optimize(problem, theta, k, maxiter, alternating=False, deadline=deadline)
     r_joint = problem.residual(*problem.unpack(theta_joint, k))
     if r_joint <= residual_atol or not dual_strategy:
         return theta_joint
-    theta_alt = _optimize(problem, theta, k, maxiter, alternating=True)
+    if deadline is not None and time.monotonic() > deadline:
+        return theta_joint
+    theta_alt = _optimize(problem, theta, k, maxiter, alternating=True, deadline=deadline)
     r_alt = problem.residual(*problem.unpack(theta_alt, k))
     if r_alt + 1e-15 < r_joint:
         return theta_alt
diff --git a/sofic/generators/generative_search.py b/sofic/generators/generative_search.py
index c7d4b05..f5a554a 100644
--- a/sofic/generators/generative_search.py
+++ b/sofic/generators/generative_search.py
@@ -667,8 +667,10 @@ def maybe_return_early() -> GenerativeModelResult | None:
                 break
             k = len(seed[0])
             for floor in _LINE_SPLIT_FLOORS:
+                if out_of_time("line-split seeds"):
+                    break
                 theta = problem.logits(*seed, floor=floor)
-                theta = _optimize_coupled(problem, theta, k, 3000, residual_atol=residual_atol)
+                theta = _optimize_coupled(problem, theta, k, 3000, residual_atol=residual_atol, deadline=deadline)
                 offer(f"{name}, floor {floor:g}", problem.unpack(theta, k))
 
     if "memory" in stages and not out_of_time("memory lifts"):
diff --git a/tests/test_coupled_generators.py b/tests/test_coupled_generators.py
index c6aa454..d8a2aaf 100644
--- a/tests/test_coupled_generators.py
+++ b/tests/test_coupled_generators.py
@@ -410,6 +410,22 @@ def counted(*args, **kwargs):
     assert calls["n"] <= 6
 
 
+def test_alm_polish_returns_at_an_expired_deadline() -> None:
+    import time
+
+    problem = cg._GeneratorProblem(_rll())
+    pi = np.full(3, 1 / 3)
+    Q = np.full((3, problem.n_plus), 1 / problem.n_plus)
+    R = np.full((3, problem.n_minus), 1 / problem.n_minus)
+    K = np.full((3, problem.n_sym, 3), 1 / (3 * problem.n_sym))
+    theta = problem.logits(pi, Q, R, K)
+    assert problem.residual(*problem.unpack(theta, 3)) > 1e-6
+    start = time.monotonic()
+    out = cg._optimize_coupled(problem, theta, 3, 3000, deadline=start - 1.0)
+    assert time.monotonic() - start < 5
+    assert np.allclose(problem.unpack(out, 3)[0], problem.unpack(theta, 3)[0])
+
+
 def test_coupled_memory_and_targeted_searches_stop_at_an_expired_deadline() -> None:
     import time
 

From 8de705c5bda7bbe7493733cd6cf24da76973d7e1 Mon Sep 17 00:00:00 2001
From: Ryan James <ryan.james@reddit.com>
Date: Fri, 13 Nov 2026 15:26:24 -0700
Subject: [PATCH 3/6] Changelog: line-split ALM deadline

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9b8aaf3..49c29df 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -42,6 +42,9 @@ brute-force reference or a closed-form value and has a regression test.
 - `memory_lift_search` stops a start in progress at its `deadline` (keeping the last finished penalty
   stage) instead of only refusing to begin new starts. Its first start used to run all four
   finite-difference L-BFGS-B stages whatever the deadline, more than 14 minutes on a 95-pair joint.
+- The line-split seeds' ALM polish (`_optimize`, `_optimize_coupled`) stops at the pipeline deadline
+  and returns its last finished step; one polish had run more than 10 minutes past the budget on a
+  144-pair joint.
 
 ### Support-type bound speed
 

From da69554f51813e1e01fe19b685fb2d323e233e52 Mon Sep 17 00:00:00 2001
From: Ryan James <ryan.james@reddit.com>
Date: Fri, 13 Nov 2026 16:51:45 -0700
Subject: [PATCH 4/6] Sparse stationary solve for large chains; ladder reuses
 the preflight build

Above 2000 states, stationary_distribution_hmm builds a sparse transition
matrix from the edges and solves pi (P - I) = 0 with one equation replaced
by the normalization. The dense path's matrix_rank and eig are O(n^3) and
never finished on a 50k-state reverse machine in the ladder. Symbolic
weights, several closed classes or a failed residual check fall back to the
dense path.

The ladder and try_build_bidirectional reuse the machines built by the
preflight (_preflight) instead of inferring and building them twice.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md                            |  8 ++++
 sofic/generators/bidirectional_build.py | 29 ++++++++------
 sofic/generators/stationary.py          | 53 +++++++++++++++++++++++++
 tests/test_bidirectional_build.py       |  6 +--
 tests/test_stationary.py                | 46 +++++++++++++++++++++
 5 files changed, 126 insertions(+), 16 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 49c29df..ac01735 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,14 @@ brute-force reference or a closed-form value and has a regression test.
 
 ### Bidirectional builds
 
+- `stationary_distribution_hmm` solves chains with more than 2,000 states sparsely (sparse transition
+  matrix from the edges, one balance equation replaced by the normalization, `spsolve`). The dense
+  path built the full matrix and ran `matrix_rank` and `eig` on it, $O(n^3)$, which never finished on
+  a 50k-state reverse machine. Symbolic weights, several closed classes, or a failed residual check
+  fall back to the dense path and its warnings.
+- The retry ladder and `try_build_bidirectional` reuse the reverse and bidirectional machines built by
+  the preflight instead of inferring and building them a second time (about half the build time).
+
 - `build_bidirectional_epsilon_machine` prunes joint states to a fixed point: dropping states
   outside the stationary support could leave others with incomplete outgoing rows, so builds from
   truncated reverse machines failed with "joint masses from ... sum to 0.86". Restriction to the
diff --git a/sofic/generators/bidirectional_build.py b/sofic/generators/bidirectional_build.py
index 5870c98..b54bc23 100644
--- a/sofic/generators/bidirectional_build.py
+++ b/sofic/generators/bidirectional_build.py
@@ -63,11 +63,18 @@ def _pe_max_drift(a: np.ndarray, b: np.ndarray) -> float:
 
 def preflight_bidir(forward: EpsilonMachine, *, reverse_max_states: int = 10_000) -> BidirPreflightResult:
     """Classify reverse presentation and whether a bidirectional machine builds at ``reverse_max_states``."""
+    return _preflight(forward, reverse_max_states)[0]
+
+
+def _preflight(
+    forward: EpsilonMachine, reverse_max_states: int
+) -> tuple[BidirPreflightResult, BidirectionalEpsilonMachine | None]:
+    """:func:`preflight_bidir` plus the bidirectional machine it built (``None`` unless it built)."""
     twins = forward.reverse_is_finite()
     recurrence = forward.reverse_epsilon_machine_is_finite()
     try:
         reverse = infer_reverse_epsilon_machine(forward, max_states=reverse_max_states)
-        build_bidirectional_epsilon_machine(forward, reverse)
+        bidir = build_bidirectional_epsilon_machine(forward, reverse)
     except MixedStateExplosionError as exc:
         return BidirPreflightResult(
             "explosion",
@@ -75,7 +82,7 @@ def preflight_bidir(forward: EpsilonMachine, *, reverse_max_states: int = 10_000
             reverse_recurrence_finite=recurrence,
             reverse_twins_finite=twins,
             error=str(exc)[:300],
-        )
+        ), None
     except (StochasticValidationError, ValueError) as exc:
         return BidirPreflightResult(
             "failed_build",
@@ -83,20 +90,20 @@ def preflight_bidir(forward: EpsilonMachine, *, reverse_max_states: int = 10_000
             reverse_recurrence_finite=recurrence,
             reverse_twins_finite=twins,
             error=str(exc)[:300],
-        )
+        ), None
     if recurrence:
         return BidirPreflightResult(
             "finite",
             reverse_max_states=reverse_max_states,
             reverse_recurrence_finite=True,
             reverse_twins_finite=twins,
-        )
+        ), bidir
     return BidirPreflightResult(
         "truncated",
         reverse_max_states=reverse_max_states,
         reverse_recurrence_finite=False,
         reverse_twins_finite=twins,
-    )
+    ), bidir
 
 
 def build_bidirectional_with_ladder(
@@ -128,13 +135,11 @@ def build_bidirectional_with_ladder(
                 exhausted = True
                 break
         attempt = time.monotonic()
-        pre = preflight_bidir(forward, reverse_max_states=cap)
+        pre, bidir = _preflight(forward, cap)
         last_preflight = pre
-        if pre.bidir_class in ("explosion", "failed_build"):
+        if bidir is None:
             previous = (cap, time.monotonic() - attempt)
             continue
-        reverse = infer_reverse_epsilon_machine(forward, max_states=cap)
-        bidir = build_bidirectional_epsilon_machine(forward, reverse)
         pe = _edge_tensor(bidir)
         previous = (cap, time.monotonic() - attempt)
         if last_pe is not None:
@@ -179,14 +184,12 @@ def try_build_bidirectional(
     :func:`build_bidirectional_with_ladder`.
     """
     start = time.monotonic()
-    pre = preflight_bidir(forward, reverse_max_states=reverse_max_states)
+    pre, bidir = _preflight(forward, reverse_max_states)
     if pre.bidir_class == "explosion" and use_ladder:
         remaining = None if time_budget is None else max(0.0, time_budget - (time.monotonic() - start))
         return build_bidirectional_with_ladder(forward, caps=caps, edge_tol=edge_tol, time_budget=remaining)
-    if pre.bidir_class in ("explosion", "failed_build"):
+    if bidir is None:
         return BidirBuildResult(None, pre)
-    reverse = infer_reverse_epsilon_machine(forward, max_states=reverse_max_states)
-    bidir = build_bidirectional_epsilon_machine(forward, reverse)
     return BidirBuildResult(bidir, pre, edge_stabilized=True, pe_drift=0.0)
 
 
diff --git a/sofic/generators/stationary.py b/sofic/generators/stationary.py
index 539dc6e..1ac576d 100644
--- a/sofic/generators/stationary.py
+++ b/sofic/generators/stationary.py
@@ -15,6 +15,10 @@
 )
 from sofic.graph import ATTR_PROB
 
+# Above this many states a numeric chain is solved with sparse matrices: the dense path builds the full
+# matrix and computes its rank and eigendecomposition (O(n^2) memory, O(n^3) time).
+_SPARSE_STATES = 2000
+
 
 def stationary_distribution_hmm(hmm: HiddenMarkovModel) -> np.ndarray:
     from sofic.properties import transition_matrix
@@ -22,10 +26,59 @@ def stationary_distribution_hmm(hmm: HiddenMarkovModel) -> np.ndarray:
     idx = hmm.reindex()
     if len(idx) == 0:
         return np.array([], dtype=float)
+    if len(idx) > _SPARSE_STATES:
+        pi = _stationary_distribution_sparse(hmm, idx.states)
+        if pi is not None:
+            return pi
     transition, _states = transition_matrix(hmm, attr=ATTR_PROB, states=idx.states)
     return stationary_distribution_from_transition(transition)
 
 
+def _stationary_distribution_sparse(hmm: HiddenMarkovModel, states: list) -> np.ndarray | None:
+    """Sparse solve of ``π P = π``, ``sum(π) = 1`` for a large numeric chain.
+
+    One balance equation is replaced by the normalization, which leaves a nonsingular system exactly when
+    the stationary distribution is unique. Returns ``None`` (so the caller falls back to the dense path,
+    with its warnings) for symbolic weights, a singular system, or a solution that is not a stationary
+    distribution.
+    """
+    import scipy.sparse as sp
+    from scipy.sparse.linalg import MatrixRankWarning, spsolve
+
+    index = {state: i for i, state in enumerate(states)}
+    rows, cols, values = [], [], []
+    for state in states:
+        for transition in hmm.graph.out_transitions(state):
+            j = index.get(transition.target)
+            if j is None:
+                continue
+            value = transition.data.get(ATTR_PROB, 0.0)
+            if has_symbolic([value]):
+                return None
+            rows.append(index[state])
+            cols.append(j)
+            values.append(float(value))
+    n = len(states)
+    matrix = sp.csr_matrix((values, (rows, cols)), shape=(n, n))
+    system = (matrix.T - sp.identity(n, format="csr")).tolil()
+    system[n - 1, :] = np.ones(n)
+    target = np.zeros(n)
+    target[-1] = 1.0
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore", MatrixRankWarning)
+        pi = spsolve(system.tocsc(), target)
+    if not np.all(np.isfinite(pi)) or np.any(pi < -1e-10):
+        return None
+    pi = np.maximum(pi, 0.0)
+    total = float(pi.sum())
+    if total <= 0.0:
+        return None
+    pi = _clean_stationary_distribution(pi / total)
+    if not np.allclose(matrix.T @ pi, pi, rtol=1e-8, atol=1e-10):
+        return None
+    return pi
+
+
 def stationary_distribution_from_transition(transition: np.ndarray) -> np.ndarray:
     """Return a normalized left eigenvector of ``transition`` for eigenvalue one.
 
diff --git a/tests/test_bidirectional_build.py b/tests/test_bidirectional_build.py
index 5e356fc..29516d0 100644
--- a/tests/test_bidirectional_build.py
+++ b/tests/test_bidirectional_build.py
@@ -67,12 +67,12 @@ def test_ladder_budget_skips_caps_forecast_to_overrun(monkeypatch: pytest.Monkey
     clock = {"now": 0.0}
     attempted: list[int] = []
 
-    def fake_preflight(forward, *, reverse_max_states):
+    def fake_preflight(forward, reverse_max_states):
         attempted.append(reverse_max_states)
         clock["now"] += 10.0 * (reverse_max_states / 10_000) ** 2
-        return bb.BidirPreflightResult("explosion", reverse_max_states=reverse_max_states)
+        return bb.BidirPreflightResult("explosion", reverse_max_states=reverse_max_states), None
 
-    monkeypatch.setattr(bb, "preflight_bidir", fake_preflight)
+    monkeypatch.setattr(bb, "_preflight", fake_preflight)
     monkeypatch.setattr(bb.time, "monotonic", lambda: clock["now"])
     built = bb.build_bidirectional_with_ladder(None, caps=(10_000, 20_000, 40_000), time_budget=100.0)
     assert attempted == [10_000, 20_000]
diff --git a/tests/test_stationary.py b/tests/test_stationary.py
index 0ef5579..798ba09 100644
--- a/tests/test_stationary.py
+++ b/tests/test_stationary.py
@@ -6,7 +6,11 @@
 
 import numpy as np
 import pytest
+from hypothesis import given, settings
+from hypothesis import strategies as st
 
+import sofic.generators.stationary as stationary
+from sofic.generators.mealy import MealyHMM
 from sofic.generators.stationary import stationary_distribution_from_transition
 
 
@@ -31,3 +35,45 @@ def test_unique_stationary_distribution_does_not_warn(transition):
         warnings.simplefilter("error", RuntimeWarning)
         pi = stationary_distribution_from_transition(transition)
     assert pi @ transition == pytest.approx(pi)
+
+
+def _chain(n: int, edges: list[tuple[int, int, float]]) -> MealyHMM:
+    hmm = MealyHMM(observation_alphabet=frozenset({0}))
+    for i in range(n):
+        hmm.graph.add_state(i)
+    for i, j, p in edges:
+        hmm.add_transition(i, j, 0, p)
+    return hmm
+
+
+def _edges(T: np.ndarray) -> list[tuple[int, int, float]]:
+    return [(i, j, float(T[i, j])) for i, j in zip(*np.nonzero(T), strict=True)]
+
+
+@settings(max_examples=40, deadline=None)
+@given(st.integers(2, 12), st.integers(0, 2**32 - 1), st.floats(0.0, 0.7))
+def test_sparse_stationary_matches_dense(n: int, seed: int, sparsity: float) -> None:
+    rng = np.random.default_rng(seed)
+    T = rng.random((n, n)) * (rng.random((n, n)) >= sparsity)
+    T[np.arange(n), (np.arange(n) + 1) % n] += 0.1  # a cycle keeps the chain irreducible
+    T /= T.sum(axis=1, keepdims=True)
+    pi = stationary._stationary_distribution_sparse(_chain(n, _edges(T)), list(range(n)))
+    assert pi is not None
+    assert np.allclose(pi, stationary_distribution_from_transition(T), atol=1e-10)
+
+
+def test_sparse_stationary_declines_two_closed_classes() -> None:
+    T = np.array([[0.5, 0.5, 0, 0], [0.5, 0.5, 0, 0], [0, 0, 0.5, 0.5], [0, 0, 0.5, 0.5]])
+    assert stationary._stationary_distribution_sparse(_chain(4, _edges(T)), [0, 1, 2, 3]) is None
+
+
+def test_large_chain_uses_the_sparse_solve() -> None:
+    import time
+
+    n = 20_000
+    edges = [(i, (i + 1) % n, 0.75) for i in range(n)] + [(i, 0, 0.25) for i in range(n)]
+    start = time.monotonic()
+    pi = stationary.stationary_distribution_hmm(_chain(n, edges))
+    assert time.monotonic() - start < 30
+    expected = 0.25 * 0.75 ** np.arange(n) / (1 - 0.75**n)
+    assert np.allclose(pi, expected, atol=1e-12)

From 52a44baa3741f0645c6fce87bb2d98665460a90d Mon Sep 17 00:00:00 2001
From: Ryan James <ryan.james@reddit.com>
Date: Fri, 13 Nov 2026 17:47:57 -0700
Subject: [PATCH 5/6] Ladder: max_reverse_states classifies oversized truncated
 builds as explosions

A ladder build whose reverse machine exceeds max_reverse_states (off by
default) is reported as an explosion so surveys use window bounds. The last
4-state survey timeout built an 8,949-state truncated reverse machine at the
50k cap, and the pipeline then spent its whole budget setting up the problem.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md                            |  4 ++++
 sofic/generators/bidirectional_build.py | 19 +++++++++++++++----
 tests/test_bidirectional_build.py       | 17 +++++++++++++++++
 3 files changed, 36 insertions(+), 4 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index ac01735..5799443 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -14,6 +14,10 @@ brute-force reference or a closed-form value and has a regression test.
   fall back to the dense path and its warnings.
 - The retry ladder and `try_build_bidirectional` reuse the reverse and bidirectional machines built by
   the preflight instead of inferring and building them a second time (about half the build time).
+- `build_bidirectional_with_ladder` and `try_build_bidirectional` take `max_reverse_states`: a ladder
+  build whose truncated reverse machine is larger is classified as an explosion, so surveys fall back
+  to window bounds. Off by default. A 4-state survey process built at the 50k cap with tens of
+  thousands of reverse states, and the pipeline then spent its whole budget setting up the problem.
 
 - `build_bidirectional_epsilon_machine` prunes joint states to a fixed point: dropping states
   outside the stationary support could leave others with incomplete outgoing rows, so builds from
diff --git a/sofic/generators/bidirectional_build.py b/sofic/generators/bidirectional_build.py
index b54bc23..605f531 100644
--- a/sofic/generators/bidirectional_build.py
+++ b/sofic/generators/bidirectional_build.py
@@ -8,7 +8,7 @@
 from __future__ import annotations
 
 import time
-from dataclasses import dataclass
+from dataclasses import dataclass, replace
 from typing import Literal
 
 import numpy as np
@@ -112,6 +112,7 @@ def build_bidirectional_with_ladder(
     caps: tuple[int, ...] = _DEFAULT_CAPS,
     edge_tol: float = 1e-8,
     time_budget: float | None = None,
+    max_reverse_states: int | None = None,
 ) -> BidirBuildResult:
     """Build a bidirectional machine, increasing reverse ``max_states`` until ``Pe`` stabilizes.
 
@@ -119,6 +120,8 @@ def build_bidirectional_with_ladder(
     the budget, and ``budget_exhausted`` is set on the result. A single mixed-state construction
     cannot be interrupted, so the forecast scales the previous cap's time by the square of the cap
     ratio (construction time grows faster than linearly in the cap); the first cap always runs.
+    A build whose reverse machine has more than ``max_reverse_states`` states is classified as an
+    explosion: such truncated machines are far too large for any generator search.
     """
     start = time.monotonic()
     last_pe: np.ndarray | None = None
@@ -136,6 +139,11 @@ def build_bidirectional_with_ladder(
                 break
         attempt = time.monotonic()
         pre, bidir = _preflight(forward, cap)
+        if bidir is not None and max_reverse_states is not None:
+            size = len(list(bidir.reverse_machine.states()))
+            if size > max_reverse_states:
+                error = f"reverse machine has {size} states, more than max_reverse_states={max_reverse_states}"
+                pre, bidir = replace(pre, bidir_class="explosion", error=error), None
         last_preflight = pre
         if bidir is None:
             previous = (cap, time.monotonic() - attempt)
@@ -177,17 +185,20 @@ def try_build_bidirectional(
     caps: tuple[int, ...] = _DEFAULT_CAPS,
     edge_tol: float = 1e-8,
     time_budget: float | None = None,
+    max_reverse_states: int | None = None,
 ) -> BidirBuildResult:
     """Single-cap build, or ladder when ``use_ladder`` and the default cap explodes.
 
-    ``time_budget`` (seconds, including the first preflight) is passed on to
-    :func:`build_bidirectional_with_ladder`.
+    ``time_budget`` (seconds, including the first preflight) and ``max_reverse_states`` are passed on
+    to :func:`build_bidirectional_with_ladder`.
     """
     start = time.monotonic()
     pre, bidir = _preflight(forward, reverse_max_states)
     if pre.bidir_class == "explosion" and use_ladder:
         remaining = None if time_budget is None else max(0.0, time_budget - (time.monotonic() - start))
-        return build_bidirectional_with_ladder(forward, caps=caps, edge_tol=edge_tol, time_budget=remaining)
+        return build_bidirectional_with_ladder(
+            forward, caps=caps, edge_tol=edge_tol, time_budget=remaining, max_reverse_states=max_reverse_states
+        )
     if bidir is None:
         return BidirBuildResult(None, pre)
     return BidirBuildResult(bidir, pre, edge_stabilized=True, pe_drift=0.0)
diff --git a/tests/test_bidirectional_build.py b/tests/test_bidirectional_build.py
index 29516d0..92b1cdc 100644
--- a/tests/test_bidirectional_build.py
+++ b/tests/test_bidirectional_build.py
@@ -87,6 +87,23 @@ def fake_preflight(forward, reverse_max_states):
     assert not built.budget_exhausted
 
 
+def test_ladder_treats_an_oversized_reverse_machine_as_an_explosion(monkeypatch: pytest.MonkeyPatch) -> None:
+    from types import SimpleNamespace
+
+    import sofic.generators.bidirectional_build as bb
+
+    huge = SimpleNamespace(reverse_machine=SimpleNamespace(states=lambda: range(5000)))
+
+    def fake_preflight(forward, reverse_max_states):
+        return bb.BidirPreflightResult("truncated", reverse_max_states=reverse_max_states), huge
+
+    monkeypatch.setattr(bb, "_preflight", fake_preflight)
+    built = bb.build_bidirectional_with_ladder(None, caps=(50_000,), max_reverse_states=1000)
+    assert built.bidir is None
+    assert built.preflight.bidir_class == "explosion"
+    assert "5000 states" in (built.preflight.error or "")
+
+
 def test_truncated_reverse_machine_prunes_to_complete_rows() -> None:
     """Regression: dropping sub-tolerance joint states left (2, 'CE') with 86% of its outgoing mass."""
     mats = {x: np.zeros((4, 4)) for x in (0, 1)}

From d471c9e1c379826f2c61277f025fba20f3dcc853 Mon Sep 17 00:00:00 2001
From: Ryan James <ryan.james@reddit.com>
Date: Sat, 14 Nov 2026 08:17:03 -0700
Subject: [PATCH 6/6] Stationary test: check the sparse path is taken instead
 of timing it (Ubuntu CI took 31 s)

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 tests/test_stationary.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/tests/test_stationary.py b/tests/test_stationary.py
index 798ba09..1f2b715 100644
--- a/tests/test_stationary.py
+++ b/tests/test_stationary.py
@@ -67,13 +67,13 @@ def test_sparse_stationary_declines_two_closed_classes() -> None:
     assert stationary._stationary_distribution_sparse(_chain(4, _edges(T)), [0, 1, 2, 3]) is None
 
 
-def test_large_chain_uses_the_sparse_solve() -> None:
-    import time
+def test_large_chain_uses_the_sparse_solve(monkeypatch: pytest.MonkeyPatch) -> None:
+    def dense(*args, **kwargs):
+        raise AssertionError("dense stationary solve on a large chain")
 
+    monkeypatch.setattr(stationary, "_stationary_distribution_numeric", dense)
     n = 20_000
     edges = [(i, (i + 1) % n, 0.75) for i in range(n)] + [(i, 0, 0.25) for i in range(n)]
-    start = time.monotonic()
     pi = stationary.stationary_distribution_hmm(_chain(n, edges))
-    assert time.monotonic() - start < 30
     expected = 0.25 * 0.75 ** np.arange(n) / (1 - 0.75**n)
     assert np.allclose(pi, expected, atol=1e-12)