Generator caps: forward joint-moment stationarity rows

CHANGELOG.md

@@ -12,6 +12,15 @@ brute-force reference or a closed-form value and has a regression test.
   truncated reverse machines failed with "joint masses from ... sum to 0.86". Restriction to the
   stationary support and to states with complete rows now alternate until neither removes anything.
 
+### Generator caps
+
+- `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.
+  Same cost; two-state caps drop (e.g. 0.918 -> 0.907 at 600 nodes) and the coverage bound on the
+  survey processes `((1,1),(2,0),(1,0))` and `((1,·),(2,0),(1,0))` rises 1.261 -> 1.277 and
+  1.175 -> 1.184.
+
 ### State labels
 
 - `sofic.states.sequential_labels` no longer stops at 26 states: after ``A``–``Z`` it continues

sofic/generators/generator_caps.py

@@ -6,7 +6,11 @@
 :mod:`~sofic.generators.coupled_generators` (edge law, forward stationarity, reverse consistency)
 constrain them further. :func:`generator_weight_cap` keeps ``k`` states explicit, aggregates every
 other state into one group, and writes those equations for per-group and per-group-pair moments
-(masses, belief moments, flows). Products of an explicit state's beliefs with flows are relaxed by
+(masses, belief moments, flows). Stationarity and reverse consistency are imposed both on belief
+moments and on the joint moments :math:`\sum_w \pi_w Q_w R_w^\top`: forward stationarity of each
+state, multiplied by that state's own :math:`R_{w'}`, is an identity between the edge moments and the
+target group's joint moment, the mirror of the reverse one. Products of an explicit state's beliefs
+with flows are relaxed by
 McCormick envelopes :cite:`mccormick1976computability` over a box of beliefs, and the boxes are
 branched. Every generator gives a feasible point (its ``k`` heaviest states explicit, the rest
 aggregated), so the largest open bound caps the total weight of the ``k`` heaviest states at any
@@ -106,6 +110,10 @@ def total(indices, minus) -> list[tuple[int, float]]:
         for h, a2 in itertools.product(groups, range(na)):
             sources = [(a, x) for a in range(na) for x in symbols if problem.F[a, x, a2]]
             eq([(int(Fl[x, g, h, a]), 1.0) for a, x in sources for g in groups] + [(int(y[h, a2]), -1.0)], 0.0)
+            # forward stationarity of each state times its own R_{w'}(b), the mirror of the reverse rows below
+            for b in range(nb):
+                terms = [(int(Gt[x, g, h, a, b]), 1.0) for a, x in sources for g in groups]
+                eq(terms + [(int(T[h, a2, b]), -1.0)], 0.0)
         for g, b in itertools.product(groups, range(nb)):
             targets = [(x, b2) for x in symbols for b2 in range(nb) if problem.G[x, b2, b]]
             eq([(int(Hr[x, g, h, b2]), 1.0) for x, b2 in targets for h in groups] + [(int(z[g, b]), -1.0)], 0.0)

tests/test_generator_caps.py

@@ -54,6 +54,15 @@ def test_generator_coverage_meets_three_halves_on_rll() -> None:
     assert cov.value == pytest.approx(1.5, abs=1e-6)
 
 
+def test_forward_joint_stationarity_tightens_the_two_state_cap() -> None:
+    problem = cg._GeneratorProblem(_survey_machine())
+    cover = functional_cover_search(problem)
+    weights = np.sort(np.bincount(cover.labels, weights=problem.p_pair))[::-1]
+    cap2 = generator_weight_cap(problem, 2, max_nodes=600)
+    # 0.918 without the forward joint-moment rows at this budget (0.907 with them on macOS)
+    assert weights[:2].sum() <= cap2 + 1e-9 < 0.915
+
+
 def test_rll_caps_are_exactly_one_half_and_three_quarters() -> None:
     rll = cg._GeneratorProblem(_rll())
     for k, nodes, cap in ((1, 200, Fraction(1, 2)), (2, 400, Fraction(3, 4))):