Clp-core: true Forrest-Tomlin (sparse Hessenberg) + solution-writer & unbounded fixes

[andig]

Clp-faithful engine core: true Forrest-Tomlin + solution-writer / unbounded fixes

Branch delivers the Clp-core rewrite (docs/rewrite-clp-core.md): Forrest-Tomlin
factorization, sparse factorize, D1 model reduction, 9/9 Cgl generators, DSE
dual — all behind env gates, default engine byte-identical. This PR's headline
work is the true Forrest-Tomlin column update plus two correctness fixes
surfaced by the ../optimizer battery models. Consolidated results in
Benchmarks at the end.

020 vs CBC — parity status (honest)

Node/tree parity: achieved — 774 nodes vs CBC's 833. Wall parity (3.6s):
not achieved, and it is a floor, not a missing fix. The campaign took 020
from 60s → 12.6s → 10.9s (gap to CBC ~16× → ~3×); every remaining lever is
either the Go-vs-C++ per-pivot kernel or a tree-perturbing change that re-rolls
020's roundoff-fragile search.

The 10.9-vs-3.6s wall factors exactly into two locked terms:

factor	cbcgo	real CBC	ratio
pivot volume	182,619	99,465	1.84×
per-pivot	59.7 µs	36.2 µs	1.65×
wall	10.9s	3.6s	3.03×

• Per-pivot 1.65× — Go sparse-triangular-solve (s -= cv[k]·y[r] random
  gather; ~45% of the solve is the dual BTRAN) vs Clp's C++ kernel. A codegen
  floor: bounds-check elimination on the hot loops measured a no-op (the cost
  is the memory gather, not the checks).
• Volume 1.84× — strong-branch probes (81748 pivots vs CBC's 55266
  strong-branch iterations, 1.5×) + the root incumbent dive. faceWalk fires once at the
  root and spends 37711 pivots to find the first feasible point where CBC's
  DiveCoefficient needs 706 (53×): the model forces incremental per-variable
  diving (RENS batch-rounding fails on 020), so faceWalk re-solves once per
  fixed integer. CBC's cheaper (crunched/hot-started) re-solves and coefficient
  selection are exactly what re-rolls 020's fragile tree. Every volume-cutting
  lever confirms it: faceWalk pivot budget → timeout (loses the root incumbent),
  CBC_CRUNCH → 846 nodes, CBC_FREECOL → 3400 nodes, probe stall/cap/DSE →
  all worse.

Reaching 3.6s needs a scope decision, not a search/correctness fix: either
accept node-parity + this floor, or authorize a cgo/C-kernel (or SIMD/assembly)
sparse BTRAN/FTRAN to close the 1.65× per-pivot. The tree-safe wins that stay
on by default (gap-scheduled diving, slices.SortFunc ratio test, pooled
propagate) took 020 from 12.6→10.9s without moving a single node.

Simplex: faster dual ratio-test sort (e9d9386, 9aedff3)

The Clp long-step / DSE ratio test sorts entering candidates by ascending
ratio every dual pivot. sort.Slice did this by reflection — 740ms / 6.8% of
the 020 solve. Microbenchmark (200-elem tie-heavy sort): sort.Slice
6743ns/3allocs → concrete sort.Sort 3269ns/1alloc → slices.SortFunc
3080ns/0allocs. Shipped slices.SortFunc (2.2× over sort.Slice,
zero-alloc). Order is byte-identical to sort.Slice — verified across 200
tie-heavy seeds and by 020 staying at 774 nodes / 182619 pivots. Numerics-
neutral, no tree re-roll.

MIP: three CBC-parity levers — diving, presolve, alloc (1b5829f, 554b68a, 74f91ba)

Measured the three remaining 020 levers against real CBC 2.10.3 and shipped
what wins on the default; gated what re-rolls 020's proof-fragile tree.

1 — Gap-scheduled diving (1b5829f, ships on). CBC runs diving heuristics
only while the gap is open. cbcgo bursted at every 256th node unconditionally;
on 020 the lone mid-tree burst fired with the incumbent already within 0.02 of
the optimum, found nothing, and cost ~26k pivots through the full
rens→faceWalk→fpump→dive fallback chain. Gating the burst on !closeGap
(incumbent >10% from the node bound) skips it — the root burst still always
runs (no incumbent yet). 020 208k→183k pivots, 12.9→10.9s, tree identical
(774 nodes), robustness 5/5; 018/021 byte-identical. Rejected by
measurement: faceWalk internal pivot budget (16×/24× rows no-op, 8× loses the
root incumbent → timeout — the root walk needs its full 37k pivots);
burst-failure backoff (the costly burst is the first mid one, nothing later
to disable).

2 — Presolve: the reduction gap is CoinPresolve column elimination, not
doubletons (74f91ba, gated). A structural census of the golden models
finds 0 duplicate rows, 0 duplicate columns, 0 equality doubletons — every
equality row has 3+ nonzeros, all 2-term rows are inequalities — matching CBC,
which prints 0 substitutions. Yet CBC reduces 020 to 1375 rows / 1399 cols
where cbcgo reaches 2032 / 2421; the ~1000-col gap is 260 free columns + 721
continuous singletons. Free-column removal (CoinPresolve empty-column
removal: fix at the cost-optimal bound, exact elimConst postsolve) ships
opt-in CBC_FREECOL=1 — optimum-preserving and ~2× on 018 (31746→13917
pivots, 314→163ms) but the perturbed model re-rolls 020's tree (774→3400
nodes), same fragility class as crunch. Fixing (105 = CBC's 105) and big-M
coefficient strengthening (301 ≈ CBC's 304, Cgl0010I-style) already match
and ship on by default.

3 — Alloc: pooled propagate queue (554b68a, ships on, numerics-neutral).
propagate allocated its worklist per call and used a reslice-front + append
FIFO that reallocated ~O(64)× per call. Rewrote it as a pooled cap-nr ring
(inQ dedups → at most nr rows queued, so the ring is exact); FIFO order is
byte-preserved → results identical (020 774 nodes / 182619 pivots). 020
total allocs 1878MB→1625MB. Wall is flat on the golden cases (GC ran
concurrently) — a memory/GC-pressure reduction, not a speedup.

Simplex/MIP: strong-branch probe cost vs CBC — Clp crunch probes (e60547e)

Head-to-head against real CBC 2.10.3 on 020 (-log 3): CBC spends 2032
strong probes / 55,266 hot-start iterations — 27 per probe and 56% of all
its simplex work — vs cbcgo's 1022 probes / 82k pivots (80/probe, 50% of
total). Probe economics already match CBC; the absolute wall gap is
per-pivot engine constants. Every CBC hot-start mechanism was then built and
measured on 020:

• Clp crunch (CBC_CRUNCH=1, ships opt-in): probes solve a row-subset LP
  with provably-redundant rows dropped (interval activity inside the row
  bounds, Clp's crunch() test behind solveFromHotStart; only slack-basic
  rows drop so the reduced basis stays square by cofactor expansion). Probe
  bounds are exactly equal to full-LP probes — not a relaxation guess (a
  first slack-basic-only version weakened seeds and blew the tree to 57s).
  020 drops ~24% of rows; 021 probe pivots 128→28. Opt-in because the changed
  probe roundoff re-rolls 020's tree lottery: CBC_MAXETAS=24 loses its
  proof. Defaults stay byte-identical (774 nodes, robustness sweep 5/5).
• Rejected by measurement: DSE/perturbed probe duals (65 piv/probe but
  +0.7s wall), probe stall-exit at 15 degenerate pivots (16/probe but weak
  pseudocost seeds → 774→1567 nodes, 31s), probe cap 100→300 (107/probe —
  pure degenerate grind, probes never converge), pre-probe refactorize à la
  markHotStart (byte-identical trees, neutral wall).
• Probe-side waste removed: probes skip Solution() extraction — 312MB
  of the 020 alloc profile was probe x/rowAct/rc/price nobody read.
• Honest pivot accounting: the attribution counters (pivot split,
  heuristic budgets' bases) now include the DSE dual's pivots
  (Stats.Dual2); node/sb numbers were blind to dual2 work before.

MIP: probing coefficient strengthening + CBC cutoff-increment gap (abaebf8)

Two CBC-parity changes that stabilize the search tree and prove 020 for the
first time. (1) CglProbing coefficient strengthening (CBC's Cgl0003I
"strengthened rows"): while probing binary x_j, one-sided rows tighten their
x_j coefficient and rhs against the probe side's propagated activity — the
implication-aware version of the static big-M tightening; multi-pass as
CglPreProcess. 018's pre-cut root bound jumps 18092 → 18291.44; fix/strengthen
counts match CBC's preprocessing on 020/021. (2) Default absolute gap 1e-5
mirroring CBC's default cutoff increment (CbcCutoffIncrement);
-allow/-ratio override. Soundness locked by TestProbeStrengthenSound
(random big-M MIPs, sampled feasible points must survive probing).

MIP: in-tree global cut rounds (3b26247)

Every 256 nodes the tree separates the globally-valid cut families (knapsack
cover, clique, zero-half, flow cover, lift-project) at the current node's
fractional point; on success the LP rebuilds and open nodes cold-restart on
the tightened relaxation. This is the globally-valid subset of CBC's in-tree
cutting — Gomory/probing node cuts are basis-/implication-local and would need
per-subtree row lifecycle management. On 020 the rounds separate nothing new
after the root (measured no-op); kept for cover/clique-separable models.
Raising the root cut pivot budget 10×→30× rows measured strictly worse on 020
(vertex-lottery incumbent regression) and stays at 10×.

MIP: reduced-sub-problem RENS as the burst primary (6fe2392)

CBC runs its heuristics as mini branch and bound on reduced problems; cbcgo's
bursts now do the same — RENS (integral-fix + rounding + node-capped sub-MIP)
leads, faceWalk is the fallback. Five earlier demotion attempts measured
wall-negative; under the substituted model the burst order no longer re-rolls
the tree (020 node counts identical per interval either way), so the
CBC-shaped ordering ships neutrally. 020 attribution: faceWalk 200k → 41k
pivots; strong-branch probes are now the top consumer (82k of 205k).

MIP: equality-chain substitution + node-local probing cuts (506643c)

substituteChains (CglPreProcess / CoinPresolve): implied-free continuous
columns substituted out of equality rows — x = (b − rest)/a drops the column
and the row, folding the expression into every other row containing x, with
exact primal/dual postsolve (x back from the pivot row, rewritten-row
activities restored per record, the dropped row's dual from stationarity;
rowMap remaps the reduced row space; reductions chain with singleton
elimination). Default-on (CBC_SUBST=0 disables). Main models substitute few
columns (the evcc bounds usually bind) but heuristic sub-MIP clones
substitute hundreds, and the change removes 020's last bad basin: 020 now
proves at 5/5 refactorize intervals, 10.1–26.2s (was 4/5 + one timeout);
defaults 020 12.6s / 774 nodes. probingCutsNode: CglProbing implication
cuts derived under node bounds through the local-cut channel; measured still
net-negative on 020 (activation repair at subtree boundaries), stays behind
CBC_LOCALCUTS.

simplex/mip: per-solve IterCap; elimfix + maxEtas 64 default; cut pruning (6d76217)

Takes 020 from 59.8s to 15.6s proven (real CBC 3.6s) with 018/021
unregressed. (1) Per-solve pivot cap (SetIterCap, honored by the primal
loop and the DSE dual; heuristic windows use 2m+200): kills the single-LP
degenerate grinds that made every model transformation a lottery re-roll —
CBC's hot-start/heuristic iteration caps. (2) Fixed-column elimination
default-on + default CBC_MAXETAS 32→64: with the cap in place elimination
is a clean win (018 0.18s, 021 2.2s, 020 proven at 15–32s on 4/5 refactorize
intervals; interval 32 times out, hence the default move — scaling made
higher intervals 021-safe). (3) Local-cut effectiveness pruning: a batch
must lift the node bound or is retracted (rows truncated, activations
stripped, two strikes disable); kept batches drop non-binding cuts. With
pruning CBC_LOCALCUTS=1 now proves 020 (57s — was a timeout) but still
loses to the 16s default, so it stays opt-in.

MIP: local in-tree cuts, RENS sub-MIP fallback, fixed-column elimination (69bd7a2)

The three remaining CBC feature gaps, built and benchmarked; two land gated
where the vertex lottery still punishes them (defaults keep the golden
numbers). (1) Node-local in-tree cuts (CBC_LOCALCUTS=1): GMI cuts from a
node's basis are stored as free rows — vacuous for propagation, heuristics and
the .sol writer — and activated per subtree via bound overrides on their
logical variables, inherited by children: the full local-cut lifecycle with no
LP row add/remove. Correct, but net-negative on 020 (degenerate node grind)
pending cut-effectiveness pruning. (2) Burst heuristics as reduced
sub-problems: RENS's failure path is now a node-capped sub-MIP on the
integral-fixed neighborhood via the shared solveNeighborhood engine
(replacing the LP dive); faceWalk stays the burst primary — five reorderings
measured wall-negative. (3) Fixed-column elimination (CBC_ELIMFIX=1):
LB==UB columns fold into row bounds with exact primal/dual postsolve (020
−364 cols, 021 −643; objectives verified exact); 021 2.9s → 1.7s, but 018/020
regress via a single-LP feasibility-pump grind on the scaled reduced model —
gated until in-solve preemption exists. The pump also gains a
between-iteration pivot budget (no-op on golden paths).

MIP: RINS as a node-capped sub-MIP + polish pivot budget (56c82f5)

RINS now clones the problem, fixes the integers where the node LP and the
incumbent agree (≥half must fix), and solves the reduced neighborhood as a
real sub-MIP (MaxNodes 50, ≤1s, incumbent as MIP start, no nesting) — CBC's
CbcHeuristicRINS design — instead of walking the full LP; it runs only
against a NEW incumbent (CBC semantics). 1-opt polish gets a pivot budget on
top of its time box. New pivot attribution (SOLVER_DEBUG) showed node
re-solves were already cheap on 020 — 12.6 pivots/node vs CBC's 53
iterations/node — and RINS drops 20k → 4.5k pivots; robustness reaches 4/5
intervals proving. Measured and rejected: strongProbeCap 100→25 (fewer
pivots, +30s wall) and three faceWalk budget variants (all wall-negative —
the walk pins the good tree path).

True Forrest-Tomlin — O(nnz) row-spike update (5cdd266, supersedes b561cb4)

Two generations of fixes. b561cb4 replaced the original dense trailing-block
re-triangularization (O(block²) per pivot, 256-row cap → refactorize almost
every pivot) with a sparse Hessenberg sweep. That was still Bartels-Golub in FT
clothing: rows never renumber, so the trailing block's subdiagonal (the old U
diagonals) is never zero and every update appended Θ(m−p) ops to the R-file
— replayed by every subsequent ftran/btran. At maxEtas=32 on 018 (m=725)
that's ~11k replayed ops per solve, 29× the eta path's per-pivot cost.

5cdd266 implements the real Forrest-Tomlin update: the leaving pivot's
row/column cycle to the last step and the row spike is eliminated against the
existing U rows — O(nnz(row)) R-ops per update. To survive step renumbering
without rewriting the factor, L and the R-file are keyed by immutable orig
rows and U by immutable basis cols; only permutation arrays and slice headers
shift per update. A lazy column-wise U index (ucRows) makes the leaving
column's entry deletion O(nnz(col)).

018 under CBC_FT=1: 23.1s → 6.7s; 021 per-pivot cost reaches parity with
the eta path (50µs vs 52µs). The remaining FT-vs-eta wall gap is node-basin
placement, not kernel cost — FT stays opt-in.

Solution-writer fix — unnamed cut rows (14241b9)

Branch-and-bound appends cut rows with empty names; solfile.Write emitted them
as 3-token lines ( 2666 0 0), and PuLP's readsol does l[3] on every
line → IndexError, so large MIP solves failed to parse entirely (the
../optimizer 020/021 cases; 018 was spared only because it needs no cuts).
PuLP maps by name and ignores unknowns, so unnamed internal rows are dropped —
matching real CBC, whose .sol carries only the original model.

simplex/mip: Clp geometric scaling (CBC_SCALE) — 021-safe, big speedups (b4ee1dc)

The evcc models carry a 1e6 constraint-coefficient range (big-M + D1). Clp
scales every model; cbcgo did not, so 021 carried ~1e-3 feasibility slop and
its proven optimum drifted with any perturbation — at higher maxEtas the
unscaled solve falsely proved a suboptimal 8.676 (a soundness bug).

Added Clp-style geometric row/column scaling, verified line-by-line against
Clp source (ClpSimplex::createRim / ClpPackedMatrix::scale): objective
*columnScale, column bounds *inverseColumnScale, row bounds *rowScale,
solution *columnScale. Scaling is internal to the LP; bounds/solution/duals/
tableau unscale at the boundary so mip/cuts work unchanged (CBC/OSI-style).
Integer columns keep colScale=1 (clean rounding); the feasibility pump uses
a new original-space objective interface. Key fix: each cut round
rebuilds+rescales the LP, so WarmSolveExtended on the pre-rebuild
(stale-scaled) basis produced garbage (021 round-1 bound 8.70→10.57, cutting
off the optimum) — cold-solve cut rounds when scaled. Impact when it landed:
018 8.8s→0.4s, 021 15.4s→2.6s and CBC's exact optimum at every maxEtas,
020's incumbent from garbage to CBC's region. Now on by default (CBC_SCALE=0
disables); well-conditioned models stay byte-identical.

simplex: pool the DSE node re-solve workspace (f7d6921)

dual2Run (the DSE dual that drives deep branch-and-bound node re-solves)
allocated ~10 fresh slices per warm re-solve. Clp keeps these as persistent
work arrays on ClpSimplex; pooled them the same way into a per-LP dual2WS
(each simplex.Build makes a fresh LP, so m/nt are frozen for its
lifetime — no size guard). Correctness-neutral. BenchmarkDual2WarmResolve
(warm bound-change re-solve, m=60): 18596 → 11639 ns/op (−37%), 25217 →
8017 B/op (−68%), 40 → 28 allocs/op. On the large ../optimizer models the
node re-solve wall is dominated by the eta FTRAN/BTRAN so the pooling is in
the noise there (018/020/021 unchanged and correct). Adds an env-gated
(CBC_CPUPROFILE) dev CPU-profile hook to cmd/cbc.

MIP: CBC-quality feasibility pump + pump-first ordering (27b01a8)

The old feasibility pump used pure L1 distance and a dual WarmSolve after each
objective swap (mathematically wrong — dual repair assumes dual feasibility a
cost change breaks), so it converged to a bad basin (evcc 020 incumbent ~42),
and faceWalk ran first on a shared deadline, starving it. Rewritten as the
Achterberg-Berthold objective pump: blend the decaying true objective into the
L1 projection, project with a new PrimalResolve (correct after a cost change),
anti-cycle by flipping the most-fractional roundings, and run the pump first
at the root with its own sub-budget (CBC's ordering). On ../optimizer 020 the
incumbent goes 42 → 0.541 (≈ CBC's optimum 0.5583). 018/021 unchanged; full
suite + PuLP green.

MIP: CglPreProcess applied uniformly — CBC/Clp parity (6b2f460)

Preprocessing (binary big-M coefficient strengthening + forcing/redundant-row
removal) was gated on a cbcgo-specific heuristic — the root-LP integer
fractionality. CBC/Clp apply CglPreProcess / presolve uniformly (gated on
structure, not the relaxation). The gate also skipped D1 on the evcc 021 case,
where cbcgo then returned a suboptimal point (8.6837) labelled Optimal — below
its own LP relaxation: a soundness bug. Removed the gate (and the CBC_D1 env,
d1Adaptive, d1FracThreshold); D1 now runs on every model. 018 and 021 then
matched CBC's raw MIP objective (021: 8.70073 vs CBC 8.70087 — was the wrong
8.6837); 020's cut bound (−0.629) matched CBC's region.

Sparser factors — Clp Markowitz least-fill pivoting (d0a7a8c)

rebuild picked the largest entry in the shortest row. Now, among the row's
max-abs entries (numerically strongest — no stability loss), it takes the
sparsest column by active-row count: Clp/CoinFactorization's Markowitz rule.
On ../optimizer 021 (3312 rows) the average factor L+U nnz drops 8429 → 6436
(~24% sparser), wall-clock neutral. Pivot choice never changes the
elimination math, so correctness is unchanged.

FT solve speedups (43cb3ba)

Profile-driven, correctness-neutral: btran's U^-T push skips zero steps (dual
pricing has a near-unit RHS), and refactorize no longer builds the sparse-LU
st.f every time under FT — it drives all solves via the FT factor, so st.f
is built lazily only as a singular-fallback (was a full luFactorize per
refactorize, ~14% of FT solve time). 021 whole-solve under CBC_FT: 9.1s → 7.4s.

Unbounded-ray verification (98c120e)

The primal loop declared Unbounded the instant the ratio test found no blocking
row, trusting a possibly-stale factorization. The optimal path already guards
against this with a refactorize+reprice recheck; this adds the symmetric guard
for unboundedness (as CBC/Clp do) — a ray must survive a clean factorization to
be reported. Resets per pivot; genuine unbounded LPs still resolve.

Safety

Default engine unchanged on the golden path: go test ./... green, PuLP
COIN_CMDTest suite green (78 passed, only the 2 pre-documented failures).
Scaling (CBC_SCALE), the extra Cgl families (CBC_CGL) and the DSE dual
are on by default; Forrest-Tomlin stays opt-in behind CBC_FT=1.

Benchmarks (evcc golden cases, Apple M4; real CBC 2.10.3)

Defaults on both solvers, wall-clock + nodes + objective:

case	main (pre-rewrite)	this branch	real CBC
018	4.9s, 18291.45	0.34s, 7 nodes, 18291.4519	0.04s, 0 nodes, 18291.4519
021	5.8s, 8.6901 (wrong)	2.2s, 3 nodes, 8.70087	0.09s, 0 nodes, 8.70083
020	60s, −140 (garbage)	10.9s, 774 nodes, 0.55835 proven	3.6s, 833 nodes, 0.55835

Tree robustness — nodes (and wall) as solve roundoff is perturbed via the
refactorize interval CBC_MAXETAS ∈ {24, 32, 48, 64, 100}:

case	before strengthening	after	real CBC (any seed)
018	27–336 nodes, 0.3–6.4s	7–10 nodes, 0.20–0.65s	0 nodes
021	3–291 nodes, 2.6–31s	3 nodes every run, 2.7–5.7s	0 nodes
020	never proven	proven 5/5, 10.3–24.6s	833 nodes, 3.6s

020 root parity: preprocessing fixes 105 = CBC's 105 (~350 rows strengthened
vs CBC's 501); root cut bound −0.884 → −0.664 in 4 rounds vs CBC −0.885 →
−0.658 in 10 passes — within 1% of CBC's closed distance. Remaining gap vs
CBC: root-closing power (CBC ends 018/021 at 0 nodes via CglPreProcess +
pump/mini-B&B) and 020 tree cost. Pivot attribution shows 020's node
re-solves are already cheap (12.6 pivots/node vs CBC's 53 iterations/node);
the bulk is faceWalk (~198k) + strong-branch probes (~112k), both measured
to resist capping (three budget variants wall-negative), plus CBC's ~950
node-local probing cuts that cbcgo's global-only in-tree rounds can't match.