One INIT: globmp/init.py serves the grid start and the synthetic study

- globmp/init.py: INIT (Algorithm 3) as one function, the shortest path over the joints for given
per-joint state sets, with exact piece tests, budget pruning (cost-to-go rule by default, plain
budget rule as an argument) and back-pointers.
- globmp/gridstart: instance.py supplies the grid values per joint, the piece cost and the piece
test (GridStates, grid_start); plain_grid.py holds the plain grid; init.py removed. reference/
reproduces (same z_h and test counts in 36 of 36 cases).
- globmp/synthetic: joints.py supplies the per-pair states (positions on a grid, zero derivatives),
piece cost, piece test and route label; init.py removed. Same class and joints as before on all
80 instances; INIT takes 0.5 to 1.3 s per instance (median per robot count).
- The name INIT is used everywhere; the ranking file is N*_i*_init.json.
- gridstart run.py --first-order calls the current Algorithm 2 again.

README.md

@@ -28,8 +28,9 @@ globmp/                      the method (pip install -e .)
   certificate.py             contacts, multipliers, global-optimality test
   check.py                   exact continuous-time feasibility check
   seed.py                    initial curves from waypoints
+  init.py                    Algorithm 3 (INIT): shortest path over the joints for given per-joint state sets
   synthetic/                 Section 4.1: robots exchanging places among spheres (2 to 200 robots)
-  gridstart/                 Section 3, Algorithm 3 (INIT): the grid start of Theorem 1 by a shortest path over the joints
+  gridstart/                 Section 3: INIT on the grid of Theorem 1 (the grid start), against the plain grid
   tests/                     fast unit tests
 applications/
   multi_robot/               Section 4.2 and appendix: multi-robot planning
@@ -73,7 +74,7 @@ continuity `C^ell_j` at junctions.
 
 | paper | code | check against the paper's runs |
 |---|---|---|
-| Sec. 4.1, Figs. 2-3 (pairs of robots among spheres, 2 to 200 robots; RRT-Connect, RRT*, GCS, SDP-MP) | `globmp/synthetic/` | default: INIT-P picks one route class, then one first-order run of Algorithms 1-2 (median 47 / 66 / 78 / 78 s per instance for 2 / 4 / 6 / 8 robots); 77 of 80 instances reach the cost of the all-class run within 1e-5, the other 3 are 0.2-3.0% above. All-class run (`--routes all`, Fig. 2): 2472 route classes, all feasible, median 58 / 76 / 79 / 95 s per class. Baselines bit-identical. Fig. 3 (scaling) is from the paper's runs, not rerun; 200 robots on CPU: 65 restoration iterations in 70 s (paper 71 s) |
+| Sec. 4.1, Figs. 2-3 (pairs of robots among spheres, 2 to 200 robots; RRT-Connect, RRT*, GCS, SDP-MP) | `globmp/synthetic/` | default: INIT picks one route class, then one first-order run of Algorithms 1-2 (median 47 / 66 / 78 / 78 s per instance for 2 / 4 / 6 / 8 robots); 77 of 80 instances reach the cost of the all-class run within 1e-5, the other 3 are 0.2-3.0% above. All-class run (`--routes all`, Fig. 2): 2472 route classes, all feasible, median 58 / 76 / 79 / 95 s per class. Baselines bit-identical. Fig. 3 (scaling) is from the paper's runs, not rerun; 200 robots on CPU: 65 restoration iterations in 70 s (paper 71 s) |
 | Sec. 3, Algorithm 3 (INIT, the grid start of Theorem 1) and its appendix table | `globmp/gridstart/` | start of the same cost as the plain grid in 27 of 27 cases (same curve except one tie); 16 free coordinates, 5.9e18 grid points: 208,929 exact piece tests |
 | Sec. 4.2 and appendix (Robotarium arena; simulated suite) | `applications/multi_robot/robotarium/`, `applications/multi_robot/multi_robot_extension/` | suite cells: same success, path length to 3e-5; Robotarium: saved plans rescored to the paper's numbers on all 20 cells |
 | appendix (Maze2D) | `applications/maze2d/` | score, success and violations identical; plans match the figure data to 1e-11 |

globmp/gridstart/README.md

@@ -5,6 +5,8 @@ points. INIT (Algorithm 3, Section 3) computes the same `z_h` without listing th
 per joint, a piece is determined by the values at its two ends, the cost is a sum over the pieces, and a curve is
 feasible exactly when every piece is. INIT is the shortest path over the joints; it tests single pieces exactly and
 skips a piece when the cost so far plus the piece cost exceeds a budget `U` (the cost of a known feasible curve).
+The default rule of `globmp.init` also adds a lower bound on the cost of the remaining pieces (cost-to-go): same
+curve, fewer piece tests.
 
 Instance: one robot in the plane passes two discs, acceleration cost (k = 2), endpoint positions fixed, `N` cubic
 pieces joined C^1, so 4N free coordinates. Budget `U = 4 E*_N`, where `E*_N` is the best certified cost over the
@@ -14,8 +16,9 @@ four route classes.
 
 | file | content |
 |---|---|
-| `init.py` | INIT (`init`; comments give the line numbers of Algorithm 3) and the plain grid (`plain_grid`) |
-| `instance.py` | the instance, the chart of free control points, the cost, the exact test of one piece |
+| `../init.py` | INIT (Algorithm 3), the package's one implementation: shortest path over the joints for given per-joint state sets |
+| `instance.py` | the instance, the chart of free control points, the cost, the exact test of one piece; `GridStates` hands INIT the grid values per joint, the piece cost and the piece test (`grid_start`) |
+| `plain_grid.py` | the plain grid (`plain_grid`), for comparison |
 | `exact_repair.py` | the exact repair step (5) as a semidefinite program (Clarabel), and GlobMP with that repair |
 | `run.py` | sanity checks, `E*_N`, one job per `(N, h)`, summary table |
 | `reference/` | `optimum.json` and `summary.json` of the run below |
@@ -60,8 +63,9 @@ Counts are deterministic; times depend on the machine.
 ## Limits
 
 - With one piece INIT is the plain grid. The work is at most `N S^2` piece tests, `S` the number of grid values at
-  one joint: linear in `N` and exponential in the number of coordinates per joint. It is not run on the team
-  instances of `globmp/synthetic/`, where a joint has 16 to 64 coordinates.
+  one joint: linear in `N` and exponential in the number of coordinates per joint. On the team instances of
+  `globmp/synthetic/` a joint has 16 to 64 coordinates, and INIT runs there on smaller per-joint state sets
+  (`globmp/synthetic/joints.py`).
 - `E*_N` is the best certified cost found, not a proved optimum.
-- `--first-order` also runs the released Algorithms 1-2 from `z_h`. At N = 2 they need more than the default
-  6000 iterations to come within 1e-5 of `E*_2` (`--T 100000`).
+- `--first-order` also runs the released Algorithms 1-2 (`globmp.synthetic.run.alg2`) from `z_h`. Checked at
+  `h = 0.5`: they end within 1e-5 of `E*_N` for N = 2 and 3.

globmp/gridstart/instance.py

@@ -10,9 +10,11 @@ Piece i has control points (Q_i, o_i, a_{i+1}, Q_{i+1}) with Q_0 = START, Q_N =
 Time scale: the package gives every piece unit duration, so cost(total duration 1) = N^3 * cost(package).
 
 Contents: the program (`build`), the chart, the cost as a quadratic (`Quad`), the exact test of one piece
-(`pieces_free`), route classes, and the exact rational check of a whole curve.
+(`pieces_free`), what INIT needs (`GridStates`: the grid values per joint, the piece cost, the piece test; `grid_start`
+calls globmp.init), route classes, and the exact rational check of a whole curve.
 """
 import math
+import time
 
 import numpy as np
 import torch
@@ -20,6 +22,7 @@ import torch
 from .. import Problem, DT
 from .. import spline as bern
 from ..check import block_minima
+from ..init import init
 
 torch.set_num_threads(1)
 
@@ -187,6 +190,69 @@ def pieces_free(B):
     return ok, nfb
 
 
+# ------------------------------------------------------------------ INIT on the grid h Z^{4N}
+def lattice(lo, hi, h):
+    ax = [np.arange(l, u + 1) * h for l, u in zip(lo, hi)]
+    return np.stack(np.meshgrid(*ax, indexing='ij'), -1).reshape(-1, len(ax))
+
+
+class GridStates:
+    """The grid h Z^{4N} inside the bounding box of {f <= U}, as INIT's input (globmp.init): the state at joint j
+    is the grid value of (a_j, Q_j) (c0 at joint 0, a_N at joint N); piece i has control points
+    (Q_i, o_i, a_{i+1}, Q_{i+1}); its cost is the acceleration cost, its test is `pieces_free`."""
+
+    def __init__(self, N, h, U):
+        self.q = q = Quad(N); self.N = self.n_pieces = N; self.lo, self.hi = lo, hi = q.box(U, h)
+        self.T = [lattice(lo[idx_joint(j, N)], hi[idx_joint(j, N)], h) for j in range(N + 1)]
+        self.root_fallback = 0; self._piece = (None,)
+
+    def n_states(self, j): return len(self.T[j])
+
+    def _source(self, i, s):
+        """(Q_i, o_i) of the states s at joint i."""
+        S = self.T[i][s]
+        return (np.broadcast_to(START, S.shape), S) if i == 0 else (S[:, 2:], 2 * S[:, 2:] - S[:, :2])
+
+    def _target(self, i):
+        """(a_{i+1}, Q_{i+1}) of every state at joint i+1, the target part of the piece cost and the target factor
+        of its cross term."""
+        if self._piece[0] != i:
+            T = self.T[i + 1]; Ta, Tq = (T[:, :2], T[:, 2:]) if i + 1 < self.N else (T, np.broadcast_to(GOAL, T.shape))
+            v = np.stack([Ta, Tq], 1)
+            self._piece = (i, Ta, Tq, np.einsum('pq,kpc,kqc->k', K4[2:, 2:], v, v), v.reshape(len(v), 4))
+        return self._piece[1:]
+
+    def cost(self, i, s, V):
+        Ta, Tq, wt, vf = self._target(i); u = np.stack(self._source(i, s), 1)
+        ws = np.einsum('pq,kpc,kqc->k', K4[:2, :2], u, u); uK = np.einsum('kpc,pq->kqc', u, K4[:2, 2:]).reshape(len(u), 4)
+        tot = 2 * uK @ vf.T; tot += (V + ws)[:, None] + wt[None, :]
+        return tot
+
+    def free_pairs(self, i, s, t):
+        Sq, So = self._source(i, s); Ta, Tq = self._target(i)[:2]
+        ok, nfb = pieces_free(np.stack([Sq, So, Ta[t], Tq[t]], 1)); self.root_fallback += nfb
+        return ok
+
+    def togo(self, j): return self.q.togo(j, self.T[j])
+
+
+def grid_start(N, h, U, rule='togo', improve_only=True):
+    """z_h by INIT (globmp.init) on the grid of spacing h with budget U.  Returns the cheapest feasible grid point x
+    (or None = FAIL: no feasible grid point with cost <= U), its cost, and the counts."""
+    t0 = time.time(); g = GridStates(N, h, U)
+    r = init(g, U * (1 + 1e-12), rule=rule, improve_only=improve_only)
+    st = dict(layers=r['layers'], fail=not r['curves'], box=[(int(a), int(b)) for a, b in zip(g.lo, g.hi)],
+              root_fallback=g.root_fallback)
+    if not r['curves']: return None, None, dict(st, time=time.time() - t0)
+    cost, path = r['curves'][0]; x = np.zeros(4 * N)
+    for j in range(N + 1): x[idx_joint(j, N)] = g.T[j][path[j]]
+    npts = 1
+    for a, b in zip(g.lo, g.hi): npts *= int(b - a + 1)
+    st.update(time=time.time() - t0, grid_points_box=npts, tests=sum(l['tested'] for l in r['layers']),
+              pairs=sum(l['pairs'] for l in r['layers']), budget_pass=sum(l['budget_pass'] for l in r['layers']))
+    return x, cost, st
+
+
 # ------------------------------------------------------------------ route classes and checks
 def side_of(p, Z):
     """route class: side of each disc at the curve's closest approach ('p' above the centre)."""

globmp/gridstart/{init.py → plain_grid.py}

@@ -1,88 +1,9 @@
-"""INIT (Algorithm 3 of the paper): the grid start z_h of Theorem 1 by a shortest path over the joints, and the
-plain grid it replaces.
-
-The grid h Z^{4N} in the free control points is a product of finite sets, one per joint.  A piece is determined by
-the values at its two ends, the cost is a sum over the pieces, and a curve is feasible exactly when every piece
-is.  `init` keeps, for every grid value at a joint, the cheapest feasible way to reach it.
-"""
+"""The plain grid that INIT replaces: every grid point with cost <= U gets a whole-curve test."""
 import time
 
 import numpy as np
 
-from .instance import GOAL, START, K4, Quad, ctrl, idx_joint, pieces_free
-
-
-def _lattice(lo, hi, h):
-    ax = [np.arange(l, u + 1) * h for l, u in zip(lo, hi)]
-    return np.stack(np.meshgrid(*ax, indexing='ij'), -1).reshape(-1, len(ax))
-
-
-def _pair_cost_parts(Sq, So, Ta, Tq):
-    """cost of the piece with control points (Sq, So, Ta, Tq): source part (ns,), target part (nt,), and the
-    factors of the cross term."""
-    Kss = K4[:2, :2]; Ktt = K4[2:, 2:]; Kst = K4[:2, 2:]
-    u = np.stack([Sq, So], 1)                                      # (ns, 2, 2) [which, coord]
-    v = np.stack([Ta, Tq], 1)                                      # (nt, 2, 2)
-    ws = np.einsum('pq,kpc,kqc->k', Kss, u, u); wt = np.einsum('pq,kpc,kqc->k', Ktt, v, v)
-    uK = np.einsum('kpc,pq->kqc', u, Kst)                          # (ns, 2, 2)
-    return ws, wt, uK.reshape(len(u), 4), v.reshape(len(v), 4)
-
-
-def init(N, h, U, rule='budget', improve_only=False, chunk_entries=2e7, verbose=False):
-    """INIT.  rule 'budget' (the rule of the paper's algorithm box): skip a pair when V(s) + w > U.
-    rule 'togo': skip when V(s) + w + togo(s') > U and drop states with V(s) + togo(s) > U, where togo is the
-    obstacle-free cost of the remaining pieces (a lower bound, so the output is the same).  improve_only: also
-    skip pairs that cannot lower the current V(s') (same output).  Returns the cheapest feasible grid point x
-    (or None = FAIL: no feasible grid point with cost <= U), its cost, and the counts."""
-    q = Quad(N); lo, hi = q.box(U, h); Ue = U * (1 + 1e-12)
-    t0 = time.time(); layers = []
-    s0 = _lattice(lo[:2], hi[:2], h)
-    V = np.zeros(len(s0)); Sq = np.broadcast_to(START, s0.shape).copy(); So = s0
-    states = [s0]; preds = []
-    if rule == 'togo':
-        keep = q.togo(0, s0) <= Ue; V = V[keep]; Sq = Sq[keep]; So = So[keep]; states[0] = s0[keep]
-    for i in range(N):
-        tl = time.time(); ij = idx_joint(i + 1, N); T = _lattice(lo[ij], hi[ij], h)
-        if i + 1 < N: Ta, Tq = T[:, :2], T[:, 2:]
-        else: Ta, Tq = T, np.broadcast_to(GOAL, T.shape).copy()
-        tg = q.togo(i + 1, T) if rule == 'togo' else np.zeros(len(T))
-        ws, wt, uK, vf = _pair_cost_parts(Sq, So, Ta, Tq)
-        Vn = np.full(len(T), np.inf); pr = np.full(len(T), -1, int)
-        st = dict(piece=i, sources=len(V), targets_box=len(T), pairs=len(V) * len(T), budget_pass=0, tested=0, free=0, root_fallback=0)
-        order = np.argsort(V, kind='stable'); cs = max(1, int(chunk_entries // max(1, len(T))))
-        for b in range(0, len(order), cs):
-            si = order[b:b + cs]
-            tot = V[si][:, None] + ws[si][:, None] + wt[None, :] + 2 * uK[si] @ vf.T      # V(s) + w        (line 3)
-            mask = tot + tg[None, :] <= Ue                                              # budget rule     (line 4)
-            st['budget_pass'] += int(mask.sum())
-            if improve_only: mask &= tot < Vn[None, :]
-            a, t = np.nonzero(mask)
-            if len(a) == 0: continue
-            B = np.stack([Sq[si][a], So[si][a], Ta[t], Tq[t]], 1)
-            ok, nfb = pieces_free(B)                                                    # exact test      (line 5)
-            st['tested'] += len(a); st['free'] += int(ok.sum()); st['root_fallback'] += nfb
-            a, t = a[ok], t[ok]
-            if len(a) == 0: continue
-            val = tot[a, t]
-            o2 = np.lexsort((val, t)); t2 = t[o2]; first = np.r_[True, t2[1:] != t2[:-1]]
-            tt = t2[first]; vv = val[o2][first]; aa = si[a[o2][first]]
-            better = vv < Vn[tt]; Vn[tt[better]] = vv[better]; pr[tt[better]] = aa[better]   # update V(s')  (line 6)
-        live = np.isfinite(Vn)
-        st.update(live_targets=int(live.sum()), time=time.time() - tl); layers.append(st)
-        if verbose: print(f'  N{N} h{h} piece {i}: {st}', flush=True)
-        if not live.any():
-            return None, None, dict(layers=layers, time=time.time() - t0, fail=True, box=[(int(a), int(b)) for a, b in zip(lo, hi)])
-        states.append(T[live]); preds.append(pr[live]); V = Vn[live]
-        if i + 1 < N: Sq = T[live][:, 2:]; So = 2 * T[live][:, 2:] - T[live][:, :2]
-    j = int(np.argmin(V)); cost = float(V[j]); x = np.zeros(4 * N)
-    for lev in range(N, -1, -1):
-        x[idx_joint(lev, N)] = states[lev][j]
-        if lev > 0: j = int(preds[lev - 1][j])
-    npts = 1
-    for a, b in zip(lo, hi): npts *= int(b - a + 1)
-    return x, cost, dict(layers=layers, time=time.time() - t0, fail=False, grid_points_box=npts,
-                         tests=sum(l['tested'] for l in layers), pairs=sum(l['pairs'] for l in layers),
-                         budget_pass=sum(l['budget_pass'] for l in layers), box=[(int(a), int(b)) for a, b in zip(lo, hi)])
+from .instance import Quad, ctrl, lattice, pieces_free
 
 
 def plain_grid(N, h, U, chunk_entries=2e7):
@@ -90,7 +11,7 @@ def plain_grid(N, h, U, chunk_entries=2e7):
     the cheapest feasible one (or None), its cost, and the counts."""
     q = Quad(N); lo, hi = q.box(U, h); m = 4 * N; t0 = time.time(); Ue = U * (1 + 1e-12)
     inner = list(range(m - 4, m)); outer = list(range(m - 4))
-    Xi = _lattice(lo[inner], hi[inner], h)
+    Xi = lattice(lo[inner], hi[inner], h)
     fi = 2 * Xi @ q.g[inner] + np.einsum('km,mn,kn->k', Xi, q.H[np.ix_(inner, inner)], Xi)
     axes = [np.arange(l, u + 1) * h for l, u in zip(lo[outer], hi[outer])]
     shape = [len(a) for a in axes]; n_outer = int(np.prod(shape)) if outer else 1

globmp/gridstart/run.py

@@ -2,8 +2,8 @@
 
   check     chart, cost, cost-to-go and piece-test sanity checks; INIT against the plain grid on three small cases
   optimum   E*_N for the chosen N: GlobMP with the exact repair from the four route-class seeds
-  job       one run per (N, h): INIT with the budget rule, INIT with the budget + cost-to-go rule, the plain grid
-            when it is small enough, then GlobMP with the exact repair from z_h
+  job       one run per (N, h): INIT (globmp.init) with the cost-to-go rule and with the plain budget rule, the
+            plain grid when it is small enough, then GlobMP with the exact repair from z_h
             (--first-order: also the released Algorithms 1-2 from z_h)
   summary   table of all jobs (markdown on stdout, summary.json in the output directory)
 
@@ -22,8 +22,8 @@ from multiprocessing import Pool
 import numpy as np
 
 from .exact_repair import exact_descent
-from .init import init, plain_grid
 from .instance import *
+from .plain_grid import plain_grid
 
 REL = 1e-4                       # success: final cost <= E*(1 + REL)
 PLAIN_MAX_PAIRS = 1e11           # run the plain budget rule only below this many state pairs
@@ -50,9 +50,9 @@ def descend_exact(args):
 
 
 def descend_first_order(args):
-    from ..synthetic.run import alg2, TOL_SCALE
+    from ..synthetic.run import alg2
     N, G, T, t_max = args; p = build(N, whiten=True); Z = G_to_Z(p, G)
-    Za, ok, log, hist = alg2(p, Z, p.tolvec() * TOL_SCALE, T=T, t_max=t_max)
+    Za, ok, log, hist = alg2(p, Z, p.tolvec(), t_max=t_max, T=T)
     out = dict(start_ok=bool(ok), iters=log.get('iters'), stop=log.get('stop'), time=log.get('time'), n_restore_fail=log.get('n_restore_fail'))
     if ok: out.update(finish(p, Za))
     return out
@@ -75,7 +75,7 @@ def stage_check(a):
     for N, h, U in ((1, 1.0, 600.0), (2, 2.0, 80.0), (2, 1.25, 60.0)):
         xb, cb, sb = plain_grid(N, h, U)
         for rule, imp in (('budget', False), ('togo', True)):
-            xd, cd, sd = init(N, h, U, rule=rule, improve_only=imp)
+            xd, cd, sd = grid_start(N, h, U, rule=rule, improve_only=imp)
             print(f'N={N} h={h} U={U} rule={rule}: grid {cb} ({sb["grid_points_sublevel"]} pts, {sb["time"]:.1f}s) | INIT {cd} '
                   f'({sd.get("tests")} tests, {sd["time"]:.1f}s) same_cost={cb == cd or (cb is not None and abs(cb - cd) < 1e-9 * cb)} '
                   f'same_point={xb is not None and np.allclose(xb, xd)}', flush=True)
@@ -104,9 +104,9 @@ def job(args):
     S = [float(np.prod(cnt[idx_joint(j, N)])) for j in range(N + 1)]
     pairs = sum(S[j] * S[j + 1] for j in range(N)); npts = float(np.prod(cnt))
     res = dict(N=N, h=h, U=U, E_star_routes=Es, states_per_joint=S, pairs=pairs, grid_points_box=npts, grid_points_sublevel_est=q.sublevel_count(U, h))
-    x, c, st = init(N, h, U, rule='togo', improve_only=True); res['init_togo'] = dict(cost=c, x=None if x is None else x.tolist(), **st)
+    x, c, st = grid_start(N, h, U, rule='togo', improve_only=True); res['init_togo'] = dict(cost=c, x=None if x is None else x.tolist(), **st)
     if pairs <= PLAIN_MAX_PAIRS:
-        x2, c2, st2 = init(N, h, U, rule='budget', improve_only=False); res['init_budget'] = dict(cost=c2, x=None if x2 is None else x2.tolist(), **st2)
+        x2, c2, st2 = grid_start(N, h, U, rule='budget', improve_only=False); res['init_budget'] = dict(cost=c2, x=None if x2 is None else x2.tolist(), **st2)
         res['budget_equals_togo'] = bool((x is None and x2 is None) or (x is not None and x2 is not None and abs(c - c2) <= 1e-9 * abs(c)))
     if npts <= BRUTE_MAX_POINTS:
         xb, cb, sb = plain_grid(N, h, U); res['grid'] = dict(cost=cb, x=None if xb is None else xb.tolist(), **sb)

globmp/init.py

@@ -0,0 +1,98 @@
+"""INIT (Algorithm 3 of the paper): the cheapest feasible curve over finite per-joint state sets, by a shortest
+path over the joints.
+
+A candidate curve is one state per joint.  Piece i is determined by the states at joints i and i+1, the cost is a
+sum over the pieces, and a curve is feasible exactly when every piece is.  INIT keeps, for every state at a joint,
+the cheapest feasible way to reach it (V) and a back-pointer.  A pair (s, s') is skipped without a piece test when
+V(s) + w(s, s') exceeds the budget U (the cost of a known feasible curve).  With the cost-to-go rule the test is
+V(s) + w(s, s') + togo(s') > U, where togo(s') is a lower bound on the cost of the remaining pieces; it skips more
+pairs and returns the same curve.
+
+The instance supplies the state sets, the piece cost and the exact piece test (an object `inst`):
+
+    inst.n_pieces         M; the joints are 0..M
+    inst.n_states(j)      number of states at joint j
+    inst.cost(i, s, V)    V(s) + w(s, s') for the states s (indices at joint i, costs so far V) and every state s'
+                          at joint i+1, w the piece cost: float (len(s), n_states(i+1))
+    inst.free(i, s)       exact test of piece i: bool (len(s), n_states(i+1)), True where the piece from s to s' is
+                          feasible (for tests that are table look-ups); or instead
+    inst.free_pairs(i, s, t)   the same test for the listed pairs (s[k], t[k]) only: bool (len(s),) (for tests that
+                          cost work per pair; INIT then lists the pairs that pass the budget)
+    inst.togo(j)          cost-to-go lower bound of every state at joint j, (n_states(j),)        (rule 'togo')
+    inst.n_labels, inst.label(i, s, lab)                                                         (optional)
+                          a state can carry a label fixed by the pieces before it (e.g. the route class so far):
+                          label returns the label of the end state of every pair, int (len(s), n_states(i+1)),
+                          or (len(s), 1) on a piece that changes no label.
+                          INIT then keeps V per (label, state) and returns the cheapest curve per final label.
+"""
+import time
+
+import numpy as np
+
+
+def init(inst, U=np.inf, rule='togo', improve_only=True, ties='first', chunk=2e7):
+    """INIT.  rule 'togo' (default): budget with cost-to-go; 'budget': the plain budget rule V(s) + w > U (same
+    curve, more piece tests).  improve_only: also skip pairs that cannot lower the current V(s') (same curve).
+    Sources are processed cheapest first, `chunk` matrix entries at a time.  ties: among ways of equal cost into a
+    state, keep the first one processed ('first', fewest piece tests) or the one from the lowest (label, state)
+    index ('state', independent of the processing order).
+
+    Returns dict(curves, layers, time): curves = {label: (cost, [state index at joint 0..M])}, the cheapest
+    feasible curve with cost <= U per final label (empty = FAIL); layers = per-piece counts (sources, targets,
+    pairs, budget_pass, tested, free, live)."""
+    M = inst.n_pieces; L = getattr(inst, 'n_labels', 1); label = getattr(inst, 'label', None); t0 = time.time()
+    pairs = hasattr(inst, 'free_pairs'); first = ties == 'first'
+    idx = np.arange(inst.n_states(0)); lab = np.zeros(len(idx), np.int64); V = np.zeros(len(idx))
+    if rule == 'togo':
+        keep = inst.togo(0) <= U; idx, lab, V = idx[keep], lab[keep], V[keep]
+    back = [(idx, lab, None)]; layers = []
+    for i in range(M):
+        tl = time.time(); nt = inst.n_states(i + 1); cols = np.arange(nt)
+        tg = inst.togo(i + 1)[None, :] if rule == 'togo' else 0.0
+        Vn = np.full((L, nt), np.inf); pred = np.full((L, nt), -1, np.int64)
+        st = dict(piece=i, sources=len(V), targets=nt, pairs=len(V) * nt, budget_pass=0, tested=0, free=0)
+        od = np.argsort(V, kind='stable'); cs = max(1, int(chunk // max(1, nt)))
+        for b in range(0, len(od), cs):
+            si = od[b:b + cs] if first else np.sort(od[b:b + cs]); s = idx[si]
+            tot = inst.cost(i, s, V[si])                                          # V(s) + w(s, s')
+            need = tot + tg <= U if U < np.inf else np.ones(tot.shape, bool)      # budget rule
+            st['budget_pass'] += int(np.count_nonzero(need))
+            lb = label(i, s, lab[si]) if label else None
+            if improve_only:
+                cur = Vn[0][None, :] if lb is None else Vn[lb[:, 0]] if lb.shape[1] == 1 else Vn[lb, cols[None, :]]
+                need &= (tot < cur) if first else (tot <= cur)
+            n_need = int(np.count_nonzero(need)); st['tested'] += n_need
+            if n_need == 0: continue
+            if pairs:                                                             # exact piece test, listed pairs
+                a, t = np.nonzero(need); ok = inst.free_pairs(i, s[a], t); a, t = a[ok], t[ok]; st['free'] += len(a)
+            else:                                                                 # exact piece test, whole rows
+                ok = inst.free(i, s) & need; n_ok = int(np.count_nonzero(ok)); st['free'] += n_ok
+                if n_ok * 20 >= ok.size:                                          # many pairs: reduce over the matrix
+                    val = np.where(ok, tot, np.inf)
+                    if lb is not None and lb.shape[1] == 1: lb = np.broadcast_to(lb, ok.shape)
+                    for l in ([0] if lb is None else np.flatnonzero(np.bincount(lb[ok], minlength=L))):
+                        vl = val if lb is None else np.where(lb == l, val, np.inf)
+                        a = vl.argmin(0); v = vl[a, cols]; better = v < Vn[l]
+                        if not first: better |= (v == Vn[l]) & (si[a] < pred[l])
+                        Vn[l, better] = v[better]; pred[l, better] = si[a[better]]
+                    continue
+                a, t = np.nonzero(ok)
+            if len(a) == 0: continue
+            val = tot[a, t]; k = t if lb is None else lb[a, t if lb.shape[1] > 1 else 0] * nt + t   # key: (label, s')
+            o = np.lexsort((val, k)); k = k[o]; head = np.r_[True, k[1:] != k[:-1]]                 # cheapest per key
+            a, k, val = a[o][head], k[head], val[o][head]
+            better = val < Vn.ravel()[k]
+            if not first: better |= (val == Vn.ravel()[k]) & (si[a] < pred.ravel()[k])
+            Vn.ravel()[k[better]] = val[better]; pred.ravel()[k[better]] = si[a[better]]
+        ll, tt = np.nonzero(np.isfinite(Vn))
+        st.update(live=len(tt), time=time.time() - tl); layers.append(st)
+        if len(tt) == 0: return dict(curves={}, layers=layers, time=time.time() - t0)
+        idx, lab, V = tt, ll, Vn[ll, tt]; back.append((idx, lab, pred[ll, tt]))
+    curves = {}
+    for l in np.unique(lab):
+        j = int(np.flatnonzero(lab == l)[np.argmin(V[lab == l])]); cost = float(V[j]); path = [0] * (M + 1)
+        for lev in range(M, -1, -1):
+            path[lev] = int(back[lev][0][j])
+            if lev > 0: j = int(back[lev][2][j])
+        curves[int(l)] = (cost, path)
+    return dict(curves=curves, layers=layers, time=time.time() - t0)

globmp/synthetic/README.md

@@ -4,9 +4,11 @@ Robots exchange places in pairs among static spheres (Problem B of the paper): m
 degree-9 polynomials, C^3 junctions, derivatives 1..4 zero at both ends, inflation margin 0.5 cm, no workspace
 bound. A route class fixes which side each pair and each crossed sphere is passed on.
 
-GlobMP here: INIT-P (a reduced INIT grid search, `init.py`) picks one route class, then one first-order run of
-Algorithms 1-2 from the INIT-P curve. No Newton step, no SQP, no polish. With `--routes all`, GlobMP solves every
-route class from a waypoint path and keeps the cheapest; that run gives `reference/` and Fig. `fig_synth1`.
+GlobMP here: INIT (Algorithm 3, `globmp/init.py`) picks one route class, then one first-order run of
+Algorithms 1-2 from the INIT curve. No Newton step, no SQP, no polish. The per-joint states given to INIT here are
+robot positions on a grid with zero derivatives, per exchanging pair, so the pick is not the exact grid optimum of
+Theorem 1. With `--routes all`, GlobMP solves every route class from a waypoint path and keeps the cheapest; that
+run gives `reference/` and Fig. `fig_synth1`.
 Baselines: RRT-Connect, RRT*, GCS on IRIS regions (their paths fitted by the cheapest collision-free curve of the
 same class) and SDP-MP (Graesdal et al.).
 
@@ -15,8 +17,8 @@ same class) and SDP-MP (Graesdal et al.).
 | file | content |
 |---|---|
 | `instances.py` | instance generators (`instance`: main comparison, M = 7, two spheres; `instance_scale`: scaling study, M = 4, three spheres), the program (`build`), route classes and their waypoint paths, exact evaluation on the true geometry |
-| `init.py` | INIT-P: route-class ranking by a dynamic program over the joints (grid h = 0.05 m) |
-| `run.py` | GlobMP on a main-comparison instance: INIT-P class (default), or every route class (`--routes all`) |
+| `joints.py` | what INIT takes here: the per-joint states (grid h = 0.05 m), the piece cost and the exact piece test of one exchanging pair; the route-class ranking |
+| `run.py` | GlobMP on a main-comparison instance: INIT class (default), or every route class (`--routes all`) |
 | `scale.py` | GlobMP end to end on a scaling instance, timed from process start; data of Fig. `fig_synth2` (GlobMP, panels a-c) |
 | `baselines.py` | RRT-Connect, RRT* (OMPL), GCS on IRIS (Drake), SDP-MP (cvxpy + Clarabel), and the post-processing of planner paths |
 | `report.py`, `figures.py` | tables and quoted numbers; `fig_synth1.pdf/png` and `fig_synth2.pdf/png` |
@@ -35,7 +37,7 @@ when it is installed, else Nimbus Roman or STIX. Every run is single-threaded. R
 ```bash
 export PYTHONPATH=. OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 MKL_NUM_THREADS=1
 
-# GlobMP: INIT-P class, one solve per instance (main way to run)
+# GlobMP: INIT class, one solve per instance (main way to run)
 python -m globmp.synthetic.run --N 6 --inst $(seq 0 19) --jobs 16 --out results/init
 
 # GlobMP on every route class (reproduces reference/ and fig_synth1)
@@ -60,17 +62,17 @@ python -m globmp.synthetic.figures --res $R --out figures
 python -m globmp.synthetic.smoke   --out results/smoke --jobs 8          # add --full for the 4-robot instance
 ```
 
-`--top K` solves the K best INIT-P classes; `--routes ppmp,pppp` solves the listed classes. Use separate output
+`--top K` solves the K best INIT classes; `--routes ppmp,pppp` solves the listed classes. Use separate output
 directories per mode: the answer is the cheapest class found in the directory. Do not run `scale.py` next to other
 heavy jobs: its wall-clock time is a reported quantity.
 
 ## Settings
 
 * Program: `Problem(2N, d=9, k=4, ell=4, M, ell_j=3, delta=0.005)` in the whitened chart; separation radius 1.2 m.
-* INIT-P (`init.py`): joint y-coordinates on a grid of spacing h = 0.05 m within 2 m of the rows, joint x-coordinates
-  from the obstacle-free optimum, zero jets, one DP per exchanging pair with exact piece tests and budget pruning
-  (1.3 x the coarse-grid cost). It drops INIT's guarantee (the reduced family stops at every joint); its curve is
-  feasible and lies in the class it reports.
+* INIT states (`joints.py`): joint y-coordinates on a grid of spacing h = 0.05 m within 2 m of the rows, joint
+  x-coordinates from the obstacle-free optimum, derivatives 1..4 zero at every joint; one INIT run per exchanging
+  pair, with the route class so far as the state label, exact piece tests (segment against disc) and the budget
+  1.3 x the cost on the grid 2h. The returned curve is feasible and lies in the class it reports.
 * Algorithm 2 (`run.py`): `Y = W - (theta / beta) grad f(W)`, `z = Algorithm1(Y)`; stages theta = 0.5, 0.2, 0.05,
   0.01, 0.002; Nesterov extrapolation with a monotone restart in the last two stages; Algorithm-1 tolerance
   continuation (package tolerance x 1e-3, 1e-3, 1e-4, 1e-6, 1e-6); Algorithm-1 steps Polyak-capped (stages 1, 4, 5)
@@ -96,13 +98,14 @@ Times are per instance or per class on one core (AMD EPYC 7502).
 
 | robots | 2 | 4 | 6 | 8 |
 |---|---|---|---|---|
-| default (INIT-P + one solve): instances with cost equal to the all-class best (within 1e-5) | 20/20 | 19/20 | 18/20 | 20/20 |
+| default (INIT + one solve): instances with cost equal to the all-class best (within 1e-5) | 20/20 | 19/20 | 18/20 | 20/20 |
 | default: median (max) time per instance, s | 47 (100) | 66 (111) | 78 (131) | 78 (131) |
 | all classes: route classes, all feasible | 168 | 320 | 640 | 1344 |
 | all classes: median (max) time per class, s | 58 (633) | 76 (639) | 79 (618) | 95 (559) |
 
 * Default: 77 of 80 instances reach the all-class best. The other 3 end in another class, 0.2% to 3.0% above it
-  (4 robots: instance 15, +3.0%; 6 robots: 8, +1.6%; 15, +0.25%). INIT-P takes 1 to 3 s of the instance time.
+  (4 robots: instance 15, +3.0%; 6 robots: 8, +1.6%; 15, +0.25%). INIT took 1 to 3 s of the instance time in that
+  run; it now takes 0.5 to 1.3 s (median per robot count).
   `--top 2` reaches the all-class best on 80/80 (median 48 / 82 / 87 / 86 s).
 * All classes: the per-instance answer over the paper's GlobMP cost is 1.00000 at the median for 2, 4, 6 and 8
   robots (max 1.00007); 9 instances end 5-13% cheaper (6 robots: 10, 13, 16, 18; 8 robots: 1, 9, 10, 11, 12). No
@@ -111,7 +114,7 @@ Times are per instance or per class on one core (AMD EPYC 7502).
 
 | file | content |
 |---|---|
-| `init.csv`, `init_summary.json` | default run per instance (80): INIT-P class, cost, ratio to the all-class best, time (INIT-P + solve); the same for `--top 2` |
+| `init.csv`, `init_summary.json` | default run per instance (80): INIT class, cost, ratio to the all-class best, time (INIT + solve); the same for `--top 2` |
 | `first_order_routes.csv` | every route class (2472): cost, minimum true clearance (mm), feasible, exactly certified, iterations, stop reason |
 | `first_order_instances.csv` | per instance (80): cheapest feasible class and cost, and the cost of the GlobMP runs behind the paper (`paper_cost`, `ratio`) |
 | `first_order_summary.json` | per robot count: ratios, certified classes, run times |
@@ -123,7 +126,7 @@ Times are per instance or per class on one core (AMD EPYC 7502).
 
 Main comparison, `<S> = N<N>_i<i>`:
 `<S>_ours_route<tag>.json/.npy` (tag in {p,m}^A: side of each pair, then of each crossed sphere;
-p = near side / above), `<S>_initp.json` (INIT-P ranking, default mode), `<S>_rrt.json` (+ `_path.npy`, `.npy`)
+p = near side / above), `<S>_init.json` (INIT ranking, default mode), `<S>_rrt.json` (+ `_path.npy`, `.npy`)
 and `<S>_rrt_post.json` (+ `_Z.npy`), the same for `rrtstar`, `<S>_gcs2.json` (+ `_curve.npy`, `_Z.npy`),
 `<S>_sdp.json` (+ `.npy`).
 Scaling, `<S> = scale_N<N>_s<seed>`: `<S>_ourscpu_b3.json` (+ `_alg12_Z.npy`), `<S>_rrtc_b3.json` and
@@ -148,6 +151,8 @@ Bernstein coefficients). `figures.py` and `report.py` also read the original cam
 * `torch.linalg.lstsq` (the fit that starts every run) is not bitwise reproducible, so GlobMP reruns agree to
   about 1e-8 relative on most classes, not bit for bit; a class with a flat end can stop up to about 2e-5 apart.
   `smoke.py` compares at 1e-4 and checks the chosen class.
+* Reruns of the released code against `reference/`: default mode, the same class on 80/80 instances and costs
+  within 7e-7 (relative); all classes, 2472/2472 feasible, class costs within 1.2e-5, per-instance best within 4e-7.
 * RRT-Connect and its post-processing reproduce exactly. GCS uses Drake's convex solvers and SNOPT rounding, and
   SDP-MP an interior-point solver: expect identical GCS results on the same Drake version and about 6 digits for
   SDP-MP.

globmp/synthetic/__init__.py

@@ -3,9 +3,9 @@ pairs among static spheres (Problem B), GlobMP against RRT-Connect, RRT*, GCS an
 from 2 to 200 robots.
 
     instances  instance generator, program, route classes and their waypoint openers, exact evaluation
-    run        GlobMP on a main-comparison instance: INIT-P picks a route class, then first-order Algorithm 2
+    run        GlobMP on a main-comparison instance: INIT picks a route class, then first-order Algorithm 2
                (--routes all: every route class)
-    init       INIT-P, the route-class ranking
+    joints     what INIT (globmp.init) needs here: per-joint states, piece cost, piece test; the route ranking
     certify    optional exact continuous-time check of returned curves (rational arithmetic); not part of the method
     exact      the exact nonnegativity test behind certify
     scale      GlobMP end to end on a scaling instance (time to feasibility, time per iteration)

globmp/synthetic/{init.py → joints.py}

@@ -1,23 +1,22 @@
-"""INIT-P: route-class ranking for the main-comparison instances (2..8 planar robots, minimum snap, M = 7 pieces).
-
-A reduced form of the paper's INIT grid search:
-  * dynamic program over the M+1 joints, one layer per piece, with back-pointers; exact piece test on the whole
-    piece interval; budget pruning (a transition is dropped when V(s) + w(s,s') + togo(s') > U, togo = the exact
-    obstacle-free lower bound on the remaining pieces, U = 1.3 x the cost of the coarse-grid (2h) solution).
-  * joint state: the robots' y-coordinates on the grid y_row + h Z, |y - y_row| <= 2 m; the x-coordinates at the
-    joints are those of the obstacle-free optimum, and the jets (derivatives 1..4) at every joint are 0.  Every
-    piece is then a + (b - a) H(s), H the degree-9 smoothstep (Bernstein coefficients 0,0,0,0,0,1,1,1,1,1): its
-    configuration path is a straight segment, the piece test is an exact segment/disc distance, and its cost is
-    c_H |b - a|^2, so the DP is exact on this family.
-  * one DP per exchanging pair (its two robots, every sphere test and the pair test; the route bits of the pair and
-    of its robots' sphere crossings are carried in the state).  Team classes are ranked by summed DP cost; the
-    separations between pairs are checked exactly on the assembled curve.
-The reduced family (curves that stop at every joint) is a strict subset of the curve space, so INIT's guarantee
-does not carry over: the best-ranked class need not hold the global optimum.  Every INIT-P curve is feasible and
-lies in the class it reports, so its cost is an upper bound.
-
-run_init writes <out>/N<N>_i<i>_initp.json (ranking, time) and returns the ranked classes with their joint
-configurations Q, from which globmp.synthetic.run starts Algorithm 2 (poly_to_Z).
+"""What INIT (globmp.init) needs on the main-comparison instances, and the start it gives (2..8 planar robots,
+minimum snap, M = 7 pieces).
+
+Per-joint states (`PairStates`, one object per exchanging pair): the y-coordinates of the pair's two robots on the
+grid y_row + h Z, |y - y_row| <= 2 m, with the derivatives 1..4 at every joint equal to 0 and the x-coordinates at
+the joints fixed to those of the obstacle-free optimum.  A piece is then a + (b - a) H(s), H the degree-9 smoothstep
+(Bernstein coefficients 0,0,0,0,0,1,1,1,1,1): its configuration path is a straight segment, the piece test is an
+exact segment/disc distance (each robot against every sphere, the pair against R_S), and the piece cost is
+c_H |b - a|^2.  A state also carries the route class so far (INIT's label: the side of the pair and of the spheres
+its robots have crossed), so INIT returns the cheapest curve of every route class of the pair.  The budget is
+1.3 x the cost of the curve INIT returns on the grid 2h.
+
+`rank_routes` adds the pairs' costs per team route class (the cost is a sum over robots), checks the separations
+between pairs exactly on the assembled curve and ranks the classes; `init_start` writes the ranking, and
+`poly_to_Z` gives the chart point from which globmp.synthetic.run starts Algorithm 2.
+
+These states are positions on a grid with zero derivatives, per exchanging pair: a strict subset of the grid of
+Theorem 1, so the pick is not the exact grid optimum z_h.  Every curve returned is feasible and lies in the class
+it reports.
 """
 import itertools
 import json
@@ -29,10 +28,12 @@ import torch
 
 from .. import DT
 from .. import spline as bern
+from ..init import init
 from . import instances as I
 
 BETA = np.array([0.] * 5 + [1.] * 5)
 CH = float(BETA @ bern.cost_gram(I.D, 4).numpy() @ BETA)          # snap energy of the degree-9 smoothstep
+CHUNK = 2e5                                                       # INIT matrix entries per step
 
 
 def setup(N, i):
@@ -65,74 +66,83 @@ def cross_piece(xs, c):
     return (j - 1, 1.0) if j > 0 else (0, 0.0)
 
 
-def dp_cluster(ins, X, yrow, rob, bits, h, Yr, U, delta=I.DELTA):
-    """INIT-P on one cluster (two robots rob = (ri, rj)).  bits: list of ('pair', None) / ('cross', (robot, c)).
-    Returns {mask: (cost_y, ypath (M+1, 2))} for every route mask reachable within the budget, and stats."""
-    M = X.shape[0] - 1; K = int(round(Yr / h)); n = 2 * K + 1; g = h * np.arange(-K, K + 1)
-    ri, rj = rob; balls = ins['balls']; R = I.RS + delta; w = CH * h * h
-    W = w * (np.arange(n)[:, None] - np.arange(n)[None, :]) ** 2   # (a, b) y-cost
-    # per piece tables
-    tabs = []
-    for m in range(M):
-        T = {}
-        for key, r in (('i', ri), ('j', rj)):
-            ya = yrow[r] + g[:, None]; yb = yrow[r] + g[None, :]; ok = np.ones((n, n), bool)
-            for c, rad in balls:
-                ok &= seg_disc_dist(X[m, r], ya, X[m + 1, r], yb, c[0], c[1]) >= rad + delta
-            T['ok' + key] = ok
+class PairStates:
+    """States, piece cost and piece test of one exchanging pair rob = (ri, rj) for globmp.init.
+    bits: the route bits of the pair, ('pair', None) or ('cross', (robot, sphere centre)).
+    State at joints 1..M-1: grid indices (b0, b1) of the two robots' y, numbered b0 n + b1; at joints 0 and M the
+    single state (K, K), the robots' rows."""
+
+    def __init__(self, ins, X, yrow, rob, bits, h, Yr=2.0, delta=I.DELTA):
+        M = X.shape[0] - 1; K = int(round(Yr / h)); n = 2 * K + 1; g = h * np.arange(-K, K + 1)
+        ri, rj = rob; R = I.RS + delta; self.w = w = CH * h * h
+        self.n_pieces = M; self.n_labels = 1 << len(bits); self.K, self.n, self.g, self.rob, self.yrow = K, n, g, rob, yrow
+        self.W = w * (np.arange(n)[:, None] - np.arange(n)[None, :]) ** 2     # piece cost of one robot, [a, b]
+        self.b0 = np.repeat(np.arange(n), n); self.b1 = np.tile(np.arange(n), n); self.db = self.b0 - self.b1 + 2 * K
         e0 = yrow[ri] - yrow[rj]; dg = h * np.arange(-2 * K, 2 * K + 1)
-        rxa = X[m, ri] - X[m, rj]; rxb = X[m + 1, ri] - X[m + 1, rj]
-        T['okp'] = seg_disc_dist(rxa, e0 + dg[:, None], rxb, e0 + dg[None, :], 0.0, 0.0) >= R
-        T['bi'] = np.zeros((n, n), np.int64); T['bj'] = np.zeros((n, n), np.int64); T['bp'] = np.zeros((4 * K + 1, 4 * K + 1), np.int64)
-        tabs.append(T)
-    for q, (kind, arg) in enumerate(bits):
-        if kind == 'pair':
-            rx = X[:, ri] - X[:, rj]; m, lam = cross_piece(rx, 0.0)
-            e0 = yrow[ri] - yrow[rj]; dg = h * np.arange(-2 * K, 2 * K + 1)
-            ry = e0 + dg[:, None] + (dg[None, :] - dg[:, None]) * lam
-            tabs[m]['bp'] |= (ry < 0).astype(np.int64) << q
-        else:
-            r, (cx, cy) = arg; m, lam = cross_piece(X[:, r], cx)
-            y = yrow[r] + g[:, None] + (g[None, :] - g[:, None]) * lam
-            tabs[m]['bi' if r == ri else 'bj'] |= (y > cy).astype(np.int64) << q
-    Bm = 1 << len(bits)
-    V = np.full((Bm, n, n), np.inf); V[0, K, K] = 0.0
-    preds = []; stats = dict(transitions=0, budget_pass=0, feasible=0)
-    for m in range(M):
-        T = tabs[m]; last = m == M - 1
-        tg = np.zeros((n, n)) if last else w * (np.arange(-K, K + 1)[:, None] ** 2 + np.arange(-K, K + 1)[None, :] ** 2) / (M - m - 1)
-        bS, a0S, a1S = np.nonzero(np.isfinite(V)); VS = V[bS, a0S, a1S]
-        dS = a0S - a1S + 2 * K
-        Vn = np.full((Bm, n, n), np.inf); Pn = np.full((Bm, n, n), -1, np.int64)
-        b0s = [K] if last else range(n); b1s = np.array([K]) if last else np.arange(n)
-        db = b1s.copy()                                                 # placeholder
-        for b0 in b0s:
-            dbi = b0 - b1s + 2 * K
-            c = VS[:, None] + W[a0S, b0][:, None] + W[a1S][:, b1s]
-            stats['transitions'] += c.size
-            bud = c + tg[b0, b1s][None, :] <= U
-            stats['budget_pass'] += int(bud.sum())
-            f = bud & T['oki'][a0S, b0][:, None] & T['okj'][a1S][:, b1s] & T['okp'][dS][:, dbi]
-            stats['feasible'] += int(f.sum())
-            if not f.any(): continue
-            c = np.where(f, c, np.inf)
-            nb = bS[:, None] | T['bi'][a0S, b0][:, None] | T['bj'][a1S][:, b1s] | T['bp'][dS][:, dbi]
-            for nm in np.unique(nb[f]):
-                cm = np.where(nb == nm, c, np.inf); j = np.argmin(cm, 0); v = cm[j, np.arange(len(b1s))]
-                cur = Vn[nm, b0, b1s]; better = v < cur
-                Vn[nm, b0, b1s[better]] = v[better]; Pn[nm, b0, b1s[better]] = j[better]
-        preds.append((bS, a0S, a1S, Pn)); V = Vn
-    out = {}
-    for mask in range(Bm):
-        if not np.isfinite(V[mask, K, K]): continue
-        path = np.zeros((M + 1, 2), int); path[M] = (K, K); st = (mask, K, K)
-        for m in range(M - 1, -1, -1):
-            bS, a0S, a1S, Pn = preds[m]; j = Pn[st]; st = (bS[j], a0S[j], a1S[j]); path[m] = (a0S[j], a1S[j])
-        assert path[0, 0] == K and path[0, 1] == K and st[0] == 0
-        ypath = np.stack([yrow[ri] + g[path[:, 0]], yrow[rj] + g[path[:, 1]]], 1)
-        touch = bool((np.abs(path - K) == K).any())
-        out[mask] = (float(V[mask, K, K]), ypath, touch)
-    return out, stats
+        # exact piece tests as tables: each robot against the spheres [a, b], the pair against R_S [a0 - a1, b0 - b1]
+        self.ok0, self.ok1, self.okp = [], [], []
+        for m in range(M):
+            for r, tab in ((ri, self.ok0), (rj, self.ok1)):
+                o = np.ones((n, n), bool)
+                for c, rad in ins['balls']:
+                    o &= seg_disc_dist(X[m, r], yrow[r] + g[:, None], X[m + 1, r], yrow[r] + g[None, :], c[0], c[1]) >= rad + delta
+                tab.append(o)
+            self.okp.append(seg_disc_dist(X[m, ri] - X[m, rj], e0 + dg[:, None], X[m + 1, ri] - X[m + 1, rj], e0 + dg[None, :], 0.0, 0.0) >= R)
+        # route bits: each is decided on one piece, by the side on which the segment passes
+        bi = np.zeros((M, n, n), np.int64); bj = np.zeros((M, n, n), np.int64); bp = np.zeros((M, 4 * K + 1, 4 * K + 1), np.int64)
+        for q, (kind, arg) in enumerate(bits):
+            if kind == 'pair':
+                m, lam = cross_piece(X[:, ri] - X[:, rj], 0.0)
+                bp[m] |= ((e0 + dg[:, None] + (dg[None, :] - dg[:, None]) * lam) < 0).astype(np.int64) << q
+            else:
+                r, (cx, cy) = arg; m, lam = cross_piece(X[:, r], cx)
+                (bi if r == ri else bj)[m] |= ((yrow[r] + g[:, None] + (g[None, :] - g[:, None]) * lam) > cy).astype(np.int64) << q
+        self.bits = [(bi[m], bj[m], bp[m]) if (bi[m].any() or bj[m].any() or bp[m].any()) else None for m in range(M)]
+
+    def _full(self, j): return 0 < j < self.n_pieces
+
+    def n_states(self, j): return self.n * self.n if self._full(j) else 1
+
+    def _ab(self, j, s):
+        """grid indices (a0, a1) of the states s at joint j."""
+        return (self.b0[s], self.b1[s]) if self._full(j) else (np.full(len(s), self.K), np.full(len(s), self.K))
+
+    def _expand(self, i, A, B, op):
+        """A [source, b0], B [source, b1] -> op(A, B) [source, target state of piece i]."""
+        if self._full(i + 1): return op(A[:, :, None], B[:, None, :]).reshape(len(A), -1)
+        return op(A[:, self.K], B[:, self.K])[:, None]
+
+    def _pair(self, i, P):
+        """P [source, b0 - b1] -> [source, target state of piece i]."""
+        return P[:, self.db] if self._full(i + 1) else P[:, 2 * self.K:2 * self.K + 1]
+
+    def cost(self, i, s, V):
+        a0, a1 = self._ab(i, s)
+        return self._expand(i, V[:, None] + self.W[a0], self.W[a1], np.add)
+
+    def free(self, i, s):
+        a0, a1 = self._ab(i, s)
+        return self._expand(i, self.ok0[i][a0], self.ok1[i][a1], np.logical_and) & self._pair(i, self.okp[i][a0 - a1 + 2 * self.K])
+
+    def label(self, i, s, lab):
+        if self.bits[i] is None: return lab[:, None]
+        a0, a1 = self._ab(i, s); bi, bj, bp = self.bits[i]
+        return self._expand(i, lab[:, None] | bi[a0], bj[a1], np.bitwise_or) | self._pair(i, bp[a0 - a1 + 2 * self.K])
+
+    def togo(self, j):
+        if not self._full(j): return np.zeros(1)
+        return self.w * ((self.b0 - self.K) ** 2 + (self.b1 - self.K) ** 2) / (self.n_pieces - j)
+
+    def curves(self, U=np.inf):
+        """INIT on this pair: {route mask: (cost, y of the two robots at the joints (M+1, 2), path touches the
+        grid border)} for every route class with a feasible curve of cost <= U, and INIT's counts."""
+        r = init(self, U, ties='state', chunk=CHUNK); out = {}
+        for mask, (cost, path) in r['curves'].items():
+            mid = np.array(path[1:-1]); ix = np.full((len(path), 2), self.K); ix[1:-1, 0] = self.b0[mid]; ix[1:-1, 1] = self.b1[mid]
+            y = np.stack([self.yrow[self.rob[0]] + self.g[ix[:, 0]], self.yrow[self.rob[1]] + self.g[ix[:, 1]]], 1)
+            out[mask] = (cost, y, bool((np.abs(ix - self.K) == self.K).any()))
+        L = r['layers']
+        return out, dict(pairs=sum(l['pairs'] for l in L), budget_pass=sum(l['budget_pass'] for l in L), tested=sum(l['tested'] for l in L))
 
 
 def clusters(ins):
@@ -154,9 +164,10 @@ def team_poly(X, ypaths, cls, N):
     return Q
 
 
-def init_p(N, i, h=0.1, Yr=2.0, slack=0.3, topk=3):
-    """INIT-P for instance (N, i): coarse pass (2h, no budget) -> U; fine pass (h, budget U) -> per-class costs per
-    cluster; team classes ranked by summed cost; inter-cluster tests on the assembled curve; top-k feasible."""
+def rank_routes(N, i, h=0.05, Yr=2.0, slack=0.3, topk=1):
+    """Route classes of instance (N, i) ranked by INIT: per pair, INIT on the grid 2h without a budget gives U,
+    then INIT on the grid h with budget U gives the cheapest curve per route class; team classes are ranked by
+    summed cost, with the separations between pairs tested on the assembled curve; the top-k feasible are returned."""
     t0 = time.time()
     ins, p, X, yrow = setup(N, i); cls = clusters(ins); A = len(ins['active']); nb = A + len(I.crossings(ins))
     constX = CH * float((np.diff(X, axis=0) ** 2).sum())
@@ -165,7 +176,7 @@ def init_p(N, i, h=0.1, Yr=2.0, slack=0.3, topk=3):
         R = []
         for c in cls:
             Uc = np.inf if UU is not None else (1 + slack) * min(v[0] for v in res[2 * h][len(R)][0].values())
-            out, st = dp_cluster(ins, X, yrow, c['rob'], c['bits'], hh, Yr, Uc); R.append((out, st))
+            out, st = PairStates(ins, X, yrow, c['rob'], c['bits'], hh, Yr).curves(Uc); R.append((out, st))
             st_all.append(dict(h=hh, rob=c['rob'], U=Uc, n_classes=len(out), **st))
         res[hh] = R
     fine = res[h]
@@ -182,7 +193,7 @@ def init_p(N, i, h=0.1, Yr=2.0, slack=0.3, topk=3):
         sig = [0] * nb
         for c, (mask, _) in zip(cls, pick):
             for q, pos in enumerate(c['pos']): sig[pos] = 1 if (mask >> q) & 1 else -1
-        ranked.append(dict(tag=I.route_tag(sig), dp_cost=cost, Q=Q, sep_margin=sep, ball_margin=ball,
+        ranked.append(dict(tag=I.route_tag(sig), init_cost=cost, Q=Q, sep_margin=sep, ball_margin=ball,
                            touch=any(v[1][2] for v in pick)))
         if len(ranked) >= topk: break
     return dict(ins=ins, p=p, ranked=ranked, n_skip=n_skip, n_combos=len(combos), stats=st_all, time=time.time() - t0)
@@ -200,16 +211,16 @@ def poly_to_Z(p, Q):
     return Z
 
 
-def run_init(N, i, out, h=0.05, top=1):
-    """INIT-P on instance (N, i): writes <out>/N<N>_i<i>_initp.json and returns dict(i, time, ranked), where each
-    ranked entry has the class tag, its DP cost and its joint configurations Q (M+1, N, 2)."""
-    r = init_p(N, i, h=h, topk=top)
-    ranked = [dict(tag=d['tag'], dp_cost=d['dp_cost'], Q=d['Q'].tolist()) for d in r['ranked']]
+def init_start(N, i, out, h=0.05, top=1):
+    """INIT on instance (N, i): writes <out>/N<N>_i<i>_init.json and returns dict(i, time, ranked), where each
+    ranked entry has the class tag, its INIT cost and its joint configurations Q (M+1, N, 2)."""
+    r = rank_routes(N, i, h=h, topk=top)
+    ranked = [dict(tag=d['tag'], init_cost=d['init_cost'], Q=d['Q'].tolist()) for d in r['ranked']]
     rec = dict(N=N, i=i, h=h, top=top, time=r['time'], n_combos=r['n_combos'], n_skip=r['n_skip'],
-               ranked=[dict(tag=d['tag'], dp_cost=d['dp_cost'], sep_margin=d['sep_margin'], ball_margin=d['ball_margin'])
+               ranked=[dict(tag=d['tag'], init_cost=d['init_cost'], sep_margin=d['sep_margin'], ball_margin=d['ball_margin'])
                        for d in r['ranked']])
     os.makedirs(out, exist_ok=True)
-    with open(os.path.join(out, f'N{N}_i{i}_initp.json'), 'w') as fh: json.dump(rec, fh, indent=0, default=float)
-    print(f"[INIT-P N{N} i{i}] {r['time']:.1f}s " + ' '.join(f"{d['tag']}:{d['dp_cost']:.3f}" for d in r['ranked']), flush=True)
+    with open(os.path.join(out, f'N{N}_i{i}_init.json'), 'w') as fh: json.dump(rec, fh, indent=0, default=float)
+    print(f"[INIT N{N} i{i}] {r['time']:.1f}s " + ' '.join(f"{d['tag']}:{d['init_cost']:.3f}" for d in r['ranked']), flush=True)
     return dict(i=i, time=r['time'], ranked=ranked)
 

globmp/synthetic/reference/init_summary.json

@@ -1,5 +1,5 @@
 {
- "_about": "INIT-P (h = 0.05 m) then one first-order Algorithm-2 solve per kept class (globmp.synthetic.run defaults: top 1; top 2 for comparison), against the cheapest class of the all-class run. ratio = cost / all-class best; equal = ratio <= 1 + 1e-5. time = INIT-P + solve(s) per instance, one core per solve (AMD EPYC 7502).",
+ "_about": "INIT (h = 0.05 m) then one first-order Algorithm-2 solve per kept class (globmp.synthetic.run defaults: top 1; top 2 for comparison), against the cheapest class of the all-class run. ratio = cost / all-class best; equal = ratio <= 1 + 1e-5. time = INIT + solve(s) per instance, one core per solve (AMD EPYC 7502).",
  "h": 0.05,
  "per_N": {
   "2": {

globmp/synthetic/reference/smoke_reference.json

@@ -1,5 +1,5 @@
 {
- "_about": "Values from the runs behind the reference results (globmp.synthetic.run: --routes all for the route classes, the default INIT-P mode for N2_i0_init; baselines and scaling as in the paper), for globmp.synthetic.smoke.",
+ "_about": "Values from the runs behind the reference results (globmp.synthetic.run: --routes all for the route classes, the default INIT mode for N2_i0_init; baselines and scaling as in the paper), for globmp.synthetic.smoke.",
  "N2_i0_ours_routes": {
   "mmm": {
    "cost": 22.8338505551219,

globmp/synthetic/report.py

@@ -3,7 +3,7 @@
 Main table, per robot count (20 instances each): raw planner path collision-free (%), returned curve feasible at
 every t (%), solved (feasible and within 1% of the lowest cost any method returns, %), median cost / GlobMP's over
 feasible curves, median minimum clearance of feasible curves (cm), median time (baselines: planner plus
-post-processing; GlobMP: INIT-P plus every solved class, slowest class in parentheses).  A dagger marks a method with
+post-processing; GlobMP: INIT plus every solved class, slowest class in parentheses).  A dagger marks a method with
 runs killed at the one-hour cap (counted as not feasible).
 Scaling table, per size: blocks, joint dimension, seeds feasible per method (baselines after post-processing),
 GlobMP's median time to its first feasible curve from process start on CPU (and on GPU where such runs exist).
@@ -43,7 +43,7 @@ def rec(RES, N, i, m):
         fe = [r for r in R if r['eval']['feasible']]
         b = min(fe, key=lambda r: r['eval']['cost']) if fe else None
         e = b['eval'] if b else None
-        ini = J(f'{RES}/N{N}_i{i}_initp.json'); t0 = ini['time'] if ini else 0.0     # INIT-P ranking (default mode)
+        ini = J(f'{RES}/N{N}_i{i}_init.json'); t0 = ini['time'] if ini else 0.0      # INIT ranking (default mode)
         return dict(clear=None, feas=bool(b), cost=e['cost'] if e else None, mincl=_mincl(e) if e else None,
                     time=t0 + sum(r['time'] for r in R), tmax=t0 + max(r['time'] for r in R))
     if m in ('rrt', 'rrtstar'):

globmp/synthetic/run.py

@@ -1,8 +1,8 @@
 """GlobMP on the main-comparison instances: paper Algorithms 1-2, first order.
 
 Route classes solved (--routes):
-  init (default)  INIT-P (init.py, grid h = 0.05 m) picks one route class from a reduced grid search, and one
-                  Algorithm-2 solve runs from its INIT-P curve (--top K solves the K best-ranked classes instead).
+  init (default)  INIT (globmp.init, with the per-joint states of joints.py, grid h = 0.05 m) picks one route
+                  class, and one Algorithm-2 solve runs from its curve (--top K solves the K best-ranked classes).
   all             every route class, each from its waypoint opener (instances.route_wp); the run behind
                   reference/ and fig_synth1.
   <tag>[,<tag>]   the listed route classes, from their waypoint openers.
@@ -29,11 +29,11 @@ Bernstein coefficients on each of 16 equal subintervals are nonnegative is in th
 block coefficients once per trial point.  Same steps, same values.
 
 Writes <out>/N<N>_i<inst>_ours_route<tag>.json and .npy (the chart point Z in the package's default chart,
-Problem.Gamma(Z) gives the Bernstein coefficients), one pair per solved class, plus <out>/N<N>_i<inst>_initp.json
-(the INIT-P ranking) in init mode; the layout that report.py, figures.py and certify.py read.  Use separate output
+Problem.Gamma(Z) gives the Bernstein coefficients), one pair per solved class, plus <out>/N<N>_i<inst>_init.json
+(the INIT ranking) in init mode; the layout that report.py, figures.py and certify.py read.  Use separate output
 directories for the two modes: best_route takes the cheapest class found in the directory.
 
-    python -m globmp.synthetic.run --N 4 --inst 3 --out results/init                  # INIT-P, one class
+    python -m globmp.synthetic.run --N 4 --inst 3 --out results/init                  # INIT, one class
     python -m globmp.synthetic.run --N 4 --inst 3 --routes all --jobs 16 --out results/synthetic
     python -m globmp.synthetic.run --N 4 --inst 3 --routes ppmp --out results/synthetic
 """
@@ -199,7 +199,7 @@ def alg2(p, Z0, tv0, t_max=T_MAX, T=T_ITERS):
                 hist.append((it, c, bool(okr), time.time() - t0))
                 if okr and full and c < best[0]: best = (c, Z.clone())
     log.update(iters=it, n_restore_fail=n_fail, n_restart=n_restart, stop=reason, time=time.time() - t0,
-               stages=stages, inner=int(p.n_it), t_alg1=float(p.t_it))
+               stages=stages, inner=int(getattr(p, 'n_it', 0)), t_alg1=float(getattr(p, 't_it', 0.0)))
     if best[1] is None: return Z, False, log, hist
     return best[1], True, log, hist
 
@@ -215,7 +215,7 @@ def route_stem(N, i, tag): return f"N{N}_i{i}_ours_route{tag}"
 
 
 def run_route(N, i, tag, out, t_max=T_MAX, start=None):
-    """One route class of main-comparison instance (N, i), started from its waypoint opener, or from the INIT-P
+    """One route class of main-comparison instance (N, i), started from its waypoint opener, or from the INIT
     joint configurations start = (rank, Q); writes its json/npy files and returns the record."""
     T0 = time.time(); ins = I.instance(N, i)
     assert len(tag) == len(ins['active']) + len(I.crossings(ins)) and set(tag) <= set('pm'), tag
@@ -223,8 +223,8 @@ def run_route(N, i, tag, out, t_max=T_MAX, start=None):
     if start is None:
         Z0 = init_waypoints(p, I.route_wp(ins, I.sig_of_tag(tag))); how = 'route_wp'
     else:
-        from .init import poly_to_Z
-        Z0 = poly_to_Z(p, np.asarray(start[1])); how = f'initp_rank{start[0]}'
+        from .joints import poly_to_Z
+        Z0 = poly_to_Z(p, np.asarray(start[1])); how = f'init_rank{start[0]}'
     Z, ok, log, hist = alg2(p, Z0, p.tolvec(), t_max)
     ev = I.evaluate(p, ins, Z)
     stem = os.path.join(out, route_stem(N, i, tag)); os.makedirs(out, exist_ok=True)
@@ -252,8 +252,8 @@ def best_route(N, i, res):
 
 def _init_job(args):
     torch.set_num_threads(1)
-    from .init import run_init
-    return run_init(*args)
+    from .joints import init_start
+    return init_start(*args)
 
 
 def _job(args):
@@ -273,9 +273,9 @@ def main(argv=None):
     ap.add_argument('--N', type=int, required=True, help='number of robots (even)')
     ap.add_argument('--inst', type=int, nargs='+', required=True, help='instance indices (the paper uses 0..19)')
     ap.add_argument('--routes', '--route', dest='routes', default='init',
-                    help="'init' (default: the INIT-P class), 'all', or comma-separated route tags")
+                    help="'init' (default: the INIT class), 'all', or comma-separated route tags")
     ap.add_argument('--top', type=int, default=1, help='init mode: number of ranked classes solved (default 1)')
-    ap.add_argument('--h', type=float, default=0.05, help='init mode: INIT-P grid spacing in m (default 0.05)')
+    ap.add_argument('--h', type=float, default=0.05, help='init mode: INIT grid spacing in m (default 0.05)')
     ap.add_argument('--out', default=os.environ.get('GLOBMP_SYNTH_OUT', 'results/synthetic'))
     ap.add_argument('--jobs', type=int, default=1, help='classes solved in parallel processes (one thread each)')
     ap.add_argument('--skip-existing', action='store_true')
@@ -298,7 +298,7 @@ def main(argv=None):
     _pool(_job, jobs, a.jobs)
     for i in a.inst:
         b = best_route(a.N, i, a.out)
-        t = f" (INIT-P {t_init[i]:.1f} s)" if i in t_init else ''
+        t = f" (INIT {t_init[i]:.1f} s)" if i in t_init else ''
         print(f"N{a.N} i{i}: " + (f"best route {b[0]} cost {b[1]['eval']['cost']:.6f}" if b else 'no feasible route') + t,
               flush=True)
 

globmp/synthetic/smoke.py

@@ -1,7 +1,7 @@
 """Smoke test: rerun a few runs of the synthetic study and compare with the reference values
 (reference/smoke_reference.json).
 
-  GlobMP on the 2-robot instance 0: the default (INIT-P class, one solve) and every route class (--routes all,
+  GlobMP on the 2-robot instance 0: the default (INIT class, one solve) and every route class (--routes all,
   8 classes; with --full also the 4-robot instance 12, 16 classes); RRT-Connect plus post-processing on the
   2-robot instance 0; the scaling run at 2 robots, seed 0.
 
@@ -48,8 +48,8 @@ def main(argv=None):
     cases = [(2, 0)] + ([(4, 12)] if a.full else [])
     q = ref['N2_i0_init']; d = os.path.join(a.out, 'init')
     RUN.main(['--N', '2', '--inst', '0', '--out', d, '--t-max', str(a.t_max)])
-    b = RUN.best_route(2, 0, d); print(f"GlobMP default N=2 instance 0: INIT-P class reference {q['cls']}, rerun {b[0] if b else None}")
-    good &= bool(b and b[0] == q['cls']) and _row('INIT-P class', q, b[1]['eval'], a.rtol)
+    b = RUN.best_route(2, 0, d); print(f"GlobMP default N=2 instance 0: INIT class reference {q['cls']}, rerun {b[0] if b else None}")
+    good &= bool(b and b[0] == q['cls']) and _row('INIT class', q, b[1]['eval'], a.rtol)
     for N, i in cases:
         RUN.main(['--N', str(N), '--inst', str(i), '--routes', 'all', '--jobs', str(a.jobs), '--out', a.out, '--t-max', str(a.t_max)])
         R = ref[f'N{N}_i{i}_ours_routes']; print(f'GlobMP N={N} instance {i}:')

globmp/tests/test_core.py

@@ -70,6 +70,49 @@ def test_certificate_on_one_disc():
     assert float(block_minima(p, Zf).min()) > -1e-9
 
 
+def test_init_against_enumeration():
+    """INIT (globmp.init) on random state sets against listing every curve, per label, for both rules."""
+    import itertools
+    from globmp.init import init
+
+    class Toy:
+        n_pieces = 4; n_labels = 2
+
+        def __init__(self, seed):
+            r = np.random.default_rng(seed); self.ns = [1, 6, 5, 6, 1]
+            self.w = [r.integers(1, 9, size=(a, b)).astype(float) for a, b in zip(self.ns[:-1], self.ns[1:])]
+            self.ok = [r.random((a, b)) < 0.6 for a, b in zip(self.ns[:-1], self.ns[1:])]
+            self.flip = r.random((self.ns[2], self.ns[3])) < 0.5          # piece 2 sets the label
+            self.tg = [np.array([min(sum(self.w[i][c[k], c[k + 1]] for k, i in enumerate(range(j, 4)))
+                                     for c in itertools.product(*[[s] if q == j else range(n) for q, n in list(enumerate(self.ns))[j:]]))
+                                 for s in range(self.ns[j])]) for j in range(5)]
+
+        def n_states(self, j): return self.ns[j]
+        def cost(self, i, s, V): return V[:, None] + self.w[i][s]
+        def free(self, i, s): return self.ok[i][s]
+        def togo(self, j): return self.tg[j]
+        def label(self, i, s, lab): return lab[:, None] | self.flip[s] if i == 2 else lab[:, None]
+
+    class ToyPairs(Toy):                                                  # the piece test given pair by pair
+        def free_pairs(self, i, s, t): return self.ok[i][s, t]
+
+    for seed in range(20):
+        t = Toy(seed); tp = ToyPairs(seed); best = {}
+        for c in itertools.product(*[range(n) for n in t.ns]):
+            if all(t.ok[i][c[i], c[i + 1]] for i in range(4)):
+                f = sum(t.w[i][c[i], c[i + 1]] for i in range(4)); l = int(t.flip[c[2], c[3]])
+                if l not in best or f < best[l]: best[l] = f
+        U = min(best.values()) * 1.5 if best else 10.0
+        want = {l: f for l, f in best.items() if f <= U}
+        for rule in ('togo', 'budget'):
+            for ties, inst in (('first', t), ('state', t), ('first', tp), ('state', tp)):
+                r = init(inst, U, rule=rule, ties=ties, chunk=7)
+                assert {l: v[0] for l, v in r['curves'].items()} == want
+                for l, (f, path) in r['curves'].items():
+                    assert all(t.ok[i][path[i], path[i + 1]] for i in range(4)) and int(t.flip[path[2], path[3]]) == l
+                    assert sum(t.w[i][path[i], path[i + 1]] for i in range(4)) == f
+
+
 if __name__ == '__main__':
     for name, f in list(globals().items()):
         if name.startswith('test_'):