Benchmark coverage: solver_direct, storage/MILP cases, sparsity sweep, speed-of-light floor

[FBumann]

bench/ (#143) ships two cases on one sink. These are the gaps, in the order they would change a decision. Each is a small addition — the runner, the parity gate and the report are already case- and arm-agnostic.

1. A solver_direct arm

The biggest gap, because the LP file is not the shipped path. bench/ currently measures write_lp against to_file(io_api='lp-polars'), and the phase split (see the build-perf issue) shows the emit phase is at parity while build is 6-10x slower. On solver_direct there is no emit phase, so the ratio users actually experience is the build ratio — unmeasured today.

Wants a matching linopy arm (io_api="direct" / to_highspy) and a peak that separates the build side from HiGHS's own model, which hard rule 4 exempts. The earlier harness measured 5.76 GB end-to-end at 35.6M variables, dominated by the solver.

2. storage — the roll self-join

roll/shift is the only locality class in the ceiling with no eager analogue in cost: xarray shifts an array, we join on dim.ord - k. It is also the shape that makes real energy models large (temporal coupling), so its absence from the ladder is the most misleading thing about the current numbers. examples/storage.yaml is the model; it needs a scaling twin under bench/models/ with data-driven dimension members.

3. A MILP through solver_direct

Solve time dwarfs build, so this measures the handoff rather than the build — a different question, and the one where being 2x slower to build is irrelevant. tests/test_milp.py has the shape. Also the only case where vtype streaming gets exercised at scale.

4. A where-density sweep

Masks are row absence here and NaN-dense in xarray. That is a claim ARCHITECTURE.md makes and nothing has ever measured. Sweeping mask density (100% / 30% / 5% active) is the one axis where the gap should be qualitative rather than 2x — the eager lane pays for every masked-out cell and we pay for none.

The dispatch model already carries where: p_max > 0; the generator only needs to zero a configurable fraction of p_max, and the ladder a density parameter.

5. A hand-written highspy/CSR arm — the speed-of-light floor

Without it, "2x slower than linopy" has no denominator: there is no way to tell engine overhead from the irreducible cost of emitting 10^7 coefficients. With it, the sentence becomes "we are at Nx the floor and linopy is at Mx", which is a claim about engineering rather than a ranking.

Not a fair modelling comparison — it hardcodes one model — which is exactly why it belongs as a floor rather than as an arm in the headline table.

6. A 10^8-variable capability run

The old harness's 107M-variable run (0.57 GB peak, a model whose dense build does not fit on the machine) is the single most persuasive number this project has ever produced, and it is currently in docs/benchmarks.md under "not reproducible". One rung above l on the existing ladder gets it back, by hand rather than in the default run.

Explicitly not planned

• CI regression gating. Shared runners make peak RSS noisy enough that the gate either fails constantly or is loosened until it catches nothing.
• Pyomo / PuLP / Calliope arms. Beating a Python-object modelling layer on build time proves nothing anyone disputes, and the high-level frameworks would mostly be measuring their own data prep. If it comes up, a paragraph in the doc answers it more honestly than a bar chart.

[FBumann]

Correction to item 1's rationale. I claimed "the phase split shows the emit phase is at parity while build is 6-10x slower" — that comparison is invalid. linopy's constraints_to_file calls con_slice.to_polars(), so the dense-to-long flattening (the work matching our whole build) happens inside its emit phase. See #146 for the correction and the profile.

The item itself stands, and for a better reason: on solver_direct there is no text-writing phase at all, and the remaining gap after #151 is concentrated in exactly that text path (duckdb printf + COPY vs polars concat_str + write_csv). So the solver_direct arm is not just "the shipped path" — it is the measurement that tells us whether the remaining gap matters to anyone.

[FBumann]

Carried over from #153 (closed, re-cut one change per PR) — a measured gap in the ladder, not a prediction.

The published ladder never measured a mask that removes anything. bench/models/dispatch.yaml guards with where: "p_max > 0" against p_max ~ uniform(50, 150), so the predicate is vacuous; transport has no where at all. Every number in docs/benchmarks.md is therefore the fully dense coordinate product — simultaneously the eager lane's best case and ours worst.

bench/models/sparse.yaml at ~20% density, top rung:

case, top rung	wall	peak
dispatch, 10M	1.00x	0.28x
transport, 9.8M	1.21x	0.48x
sparse, 2.1M live	0.57x	0.20x

docs/benchmarks.md listed a where-density sweep under "not measured yet" and predicted the gap would be "qualitative rather than 2x". It is — in our favour.

The case is also worth keeping for a second reason: it reads the leading dim in its mask, so it exercises the general windowed label path rather than either arithmetic fast path (#152, #173). Without it the ladder covers only the frames those two optimise.

[FBumann]

Three of the four gaps are covered on #189, with numbers.

Sparsity sweep — nodal at 12/6/3/1 technologies per node. The result is not the expected one: linopy's peak falls with density too, and at the sparsest rung ends up below ours (0.35 vs 0.40 GB). The sweep holds the coordinate product at 1.2M, where a dense array over it is ~10 MB and the interpreter dominates, so it cannot show the effect whatever the density.

A mixed-density case — new sector: 50 nodes x 12 technologies at 8% installed, crossed with 5 dense carriers and dense snapshots, so sum(p * produces, over=tech) broadcasts a dense dim onto a sparse frame and then reduces the sparse one. Same 8% sparsity at a 12M product, and there it is unmistakable:

variables	peak: farkas	peak: linopy
100k	0.42 GB	0.53 GB
1M	1.27 GB	2.79 GB

Together the two locate the claim: it needs low density and a product large enough for it to cost anything. Either alone is not evidence.

solver_direct / speed-of-light floor — measured by hand rather than in the harness, with Highs.run stubbed so the simplex workspace is excluded and HiGHS's own model is not. dispatch at 10M: build 2.12 GB, build+handoff 3.13 GB, against linopy's 0.75 / 3.38 GB through io_api='direct'. Also worth recording: linopy's io_api='lp' costs it 6.92 GB and 55 s against 3.38 GB and 2.9 s direct, so an LP-file comparison is not a comparison against its real cost.

Folding that into the harness as --mode handoff is the obvious next change and is not done.

Storage / MILP cases — still missing.

[FBumann]

Status check against main @ 8ecc2af. Three of the six landed, mostly in #189; three are open. The issue stays open for those three — but item 1's rationale inverted along the way, which is worth correcting before anyone reads it as a live argument.

Done

1. A solver_direct arm. ✅ bench/run.py takes --sinks lp highs, both by default, and bench/_run_case.py gives it a matching linopy arm — both lanes end holding a populated highspy.Highs with run() never called, ours through build_highs and linopy's through to_highspy(), so the simplex is excluded on both sides. docs/benchmarks.md carries the solver_direct table per case.

4. A where-density sweep. ✅ Case.density is the fraction of the coordinate product a mask keeps, _density_sweep runs 1.0 / 0.5 / 0.25 / 0.083 on nodal at fixed size so a density effect can be told from a size effect, and bench/results/density.jsonl is committed. The runner reports the density actually achieved rather than assuming the requested one.

6. A 10^8-variable capability run. ✅ Back and reproducible — bench/results/scaling.jsonl and docs/benchmarks-scaling.html, out to 120M variables, where the lane builds and loads in 35 s against linopy's 106.

Still open

2. storage — the roll self-join. There is no roll( in any file under bench/models/, so the one locality class with no eager analogue in cost is still absent from the ladder. examples/storage.yaml is the model; it needs a scaling twin with data-driven dimension members. Unchanged, and still the most misleading gap.

3. A MILP through solver_direct. No bench model declares an integer or binary variable, so vtype streaming is still unexercised at scale.

5. The hand-written highspy/CSR floor. Not started.

Correcting item 1's rationale

Item 1 argued the solver_direct arm mattered because "the phase split shows the emit phase is at parity while build is 6-10x slower", so the ratio users actually experience is the build ratio. That was true of the duckdb engine and is not true now — #146 is closed, and #189 inverted it: the lane is ahead of linopy on wall on five of six cases to a loaded solver.

The arm was still worth building, for the reason underneath the number rather than the number itself: the LP file is not the shipped path, so measuring only write_lp measures a path nobody takes. That argument is engine-independent and survives. The "6-10x" figure should not be quoted from this issue.