feat(sinks): gurobi, and the two families it revealed

CONTRIBUTING.md

@@ -12,8 +12,11 @@ uv run pre-commit install  # once per clone
 ```
 
 `uv sync` installs the `[linopy]` extra too, because the differential test
-suite needs a second lane to compare against. The engine itself never imports
-linopy, xarray or pandas — see *the bare install* below.
+suite needs a second lane to compare against, and `[gurobi]`, because the
+second solver sink needs a second solver to compare against — gurobipy's wheel
+carries a size-limited licence, so those tests run on a plain checkout with no
+licence of your own. They skip where it is absent. The engine itself never
+imports linopy, xarray, pandas or gurobipy — see *the bare install* below.
 
 ## The loop
 

README.md

@@ -167,6 +167,7 @@ To see it rather than read it, `python examples/walkthrough.py` runs one small m
 pip install lpspec  # the relational engine (polars, highspy)
 pip install "lpspec[linopy]"  # adds linopy + xarray + pandas: the shim, the
                               # oracle, and to_pandas / to_dataarray
+pip install "lpspec[gurobi]"  # adds the gurobi sink: solver_name='gurobi'
 ```
 
 Not a solver wrapper, not a domain package, not a data-loading layer — bring

bench/README.md

@@ -44,7 +44,7 @@ is the hard part.
 
 | | `lp` | `highs` |
 |---|---|---|
-| `lpspec` | `lps.build(...)` then `ex.write_lp(...)` | `lps.build(...)` then `build_highs(...)` |
+| `lpspec` | `lps.build(...)` then `ex.write(...)` | `lps.build(...)` then `build_highs(...)` |
 | `linopy` | `lpspec.linopy.build(...)` then `Model.to_file(io_api='lp-polars')` | `lpspec.linopy.build(...)` then `Model.to_highspy()` |
 
 **The `highs` sink stops at the handoff — `run()` is never called.** That is the
@@ -98,7 +98,7 @@ the other never does. The boundaries are therefore explicit:
 | **before the clock** | splitting parquet paths into parameters vs dimensions (harness bookkeeping — it re-parses the YAML only because the *runner* decides which file is which) | — |
 | `import` | `import lpspec` | `import lpspec.linopy` → linopy, xarray |
 | `build` | `lps.build(...)` — the engine scans the parquet itself | `read_parquet` + reshape + `lpspec.linopy.build(...)` |
-| `emit` | `ex.write_lp(path)` / `build_highs(ex._tables())` | `Model.to_file(path, io_api='lp-polars')` / `Model.to_highspy()` |
+| `emit` | `ex.write(path)` / `build_highs(ex._tables())` | `Model.to_file(path, io_api='lp-polars')` / `Model.to_highspy()` |
 | `teardown` | `ex.close()` — releases the built model | — (nothing to release) |
 | **after the clock** | row, column and nonzero counts off the built frames | `nvars` / `ncons` |
 
@@ -194,7 +194,7 @@ Data is generated deterministically (a blake2b digest of the shape seeds the
 RNG — `hash()` is salted per process and would give the two arms different
 numbers), cached under `bench/.cache/`, and feasible by construction.
 
-Storage (`roll`, the bounded-halo self-join) and a MILP through `solver_direct`
+Storage (`roll`, the bounded-halo self-join) and a MILP through the `highs` solver
 are the next rungs — see docs/benchmarks.md.
 
 ## The other harness: regression benchmarks

bench/_run_case.py

@@ -60,7 +60,7 @@ def _run_lpspec(
     case: Case, paths: dict[str, str], lp: Path, phases: Phases, opts: argparse.Namespace
 ) -> dict[str, Any]:
     import lpspec as lps
-    from lpspec.relational.sinks.highs import build_highs
+    from lpspec.relational.sinks.solvers.highs import build_highs
 
     # The parameter/dimension split is harness bookkeeping — it re-parses the
     # YAML only because the runner, not lpspec, decides which parquet file is
@@ -72,7 +72,7 @@ def _run_lpspec(
     ex = lps.build(case.model, sources, coords=coords)
     phases.mark('build')
     if opts.sink == 'lp':
-        ex.write_lp(lp)
+        ex.write(lp)
     else:
         # the hand-off, and nothing past it. `run()` is never called: the
         # simplex is the same work whoever filled the model, so timing it would

bench/profile_build.py

@@ -106,7 +106,7 @@ def main() -> None:
             build = time.perf_counter() - started
             phase['now'] = 'emit'
             started = time.perf_counter()
-            executor.write_lp(Path(tmp) / 'model.lp')
+            executor.write(Path(tmp) / 'model.lp')
             emit = time.perf_counter() - started
 
     print(f'\n{args.case}/{args.size}: build {build:.2f}s, emit {emit:.2f}s')

bench/regressions/test_build.py

@@ -73,7 +73,7 @@ def build_and_write(case_name: str, size: str, sources: dict[str, str], coords:
         tempfile.TemporaryDirectory(prefix='lpspec-bench-') as tmp,
         lps.build(CASES[case_name].model, sources, coords=coords) as ex,
     ):
-        ex.write_lp(Path(tmp) / 'model.lp')
+        ex.write(Path(tmp) / 'model.lp')
         return ex._tables().column_count
 
 
@@ -86,7 +86,7 @@ def build_and_hand_over(case_name: str, size: str, sources: dict[str, str], coor
     be watching a number nothing in this repository can move.
     """
     import lpspec as lps
-    from lpspec.relational.sinks.highs import build_highs
+    from lpspec.relational.sinks.solvers.highs import build_highs
 
     with lps.build(CASES[case_name].model, sources, coords=coords) as ex:
         tables = ex._tables()

docs/ARCHITECTURE.md

@@ -80,8 +80,8 @@ flowchart TB
             BIND["binding.py<br/>→ BoundSources, frozen"] --> EXEC["executor.py + labels.py<br/>assemble the model frames"]
         end
         ENG --> TABLES["sinks/tables.py<br/>cols · obj · rows · A"]
-        TABLES --> LPS["sinks/lp_file.py<br/>(mps planned)"]
-        TABLES --> DIRECT["sinks/highs.py<br/>COO batches → HiGHS"]
+        TABLES --> LPS["sinks/writers/<br/>a file, chosen by suffix<br/>lp_file (mps planned)"]
+        TABLES --> DIRECT["sinks/solvers/<br/>COO batches → the solver, chosen by name<br/>highs (ships) · gurobi (extra)"]
         DIRECT --> SOL["result.py<br/>label join, never dense"]
     end
 
@@ -162,17 +162,17 @@ costs seconds, and needs nothing but the file.
 members** — including the ones nobody has built, which is the point: none of
 them is a rewrite. Each reads the same AST the engine reads, so a renderer is a
 tree walk, a check is a pass with no data bound, and a new output format is one
-function in `relational/sinks/`.
+module in `relational/sinks/writers/`.
 
 `typeset/` is that claim cashed — a **spike** that typesets any model the lanes
 can build, in one walk of the resolved AST, holding no opinion the lanes do not
 already hold: a `piecewise:` block prints as the λ-formulation it expands to,
 not as the sugar it was written as. How names *print* is the one thing it does
 not read off the model, since a symbol table is presentation — hence a sidecar
 file (`examples/symbols/`) rather than keys on `MathSchema`, and a model with no
-table still renders. It splits the way `relational/sinks/` does, one module per
-output format, so a format is a spelling table rather than a second walk that
-could disagree. `python -m lpspec <format>` is its shell front, one verb per
+table still renders. It splits the way `relational/sinks/writers/` does, one
+module per output format, so a format is a spelling table rather than a second
+walk that could disagree. `python -m lpspec <format>` is its shell front, one verb per
 entry in `typeset.FORMATS`: a consumer that needs no data needs no runner.
 
 That claim is enforced twice, because a renderer that imports only `language/`
@@ -258,7 +258,7 @@ choice load-bearing in the language's rulebook.
    consumer that binds no data must reach it without reaching a runner; `api.py`
    re-exports it so callers keep saying `lps.load_schema`.
 2. **The engine knows nothing about linopy, xarray or YAML.** `relational/` goes
-   plan → engine → highspy → solver, with linopy's semantics as a spec to match
+   plan → engine → a solver sink → solver, with linopy's semantics as a spec to match
    rather than code to share; it never sees the schema, the AST, or the eager
    builder. **The engine is a directory, not a convention:** `engines/polars/`
    is one implementation, and everything above it — `plan.py`, `sinks/`,
@@ -304,7 +304,8 @@ That split is what makes the ceiling's admissibility test something you can
 *perform* rather than reason about: build a `PolarsCompiler`, hand it a node,
 read `.explain()` — `tests/test_compiler.py` does exactly that over empty
 frames, since a schema is all it takes to compile a query. It is also why a new
-sink is a function in one file rather than another method on the executor.
+sink is a module in one of two families rather than another method on the
+executor.
 
 **What binding produces is a value.** `BoundSources` is frozen — parameters,
 dimensions, their cardinalities, and which parameters are boolean — because a
@@ -369,6 +370,24 @@ per sink (see "Capability is not the ceiling"); that unevenness is what
 [Track 3](ROADMAP.md#track-3--capabilities-and-the-degree-line) exists to make declared rather
 than discovered at solve time.
 
+**A sink is one of two things, and the directory says which.** A **solver**
+runs the tables and returns an answer, chosen by **name** at the call
+(`solver_name='gurobi'`); a **writer** renders them to a file, chosen by the
+output's **suffix** — because a file's format is a property of the file, while
+which solver runs is a property of nothing but the call. Both sets are closed
+dict literals (`SOLVERS`, `WRITERS`): no YAML key names a solver, and nothing
+installed may change what either resolves to.
+
+The split is a directory rather than a convention for the reason `engines/` is:
+**how many solvers there are will change, and what a solver has to answer will
+not.** A new one is a module named for it and a line in `SOLVERS` — no method
+on the executor, no branch in `api.py`, no name on the Python surface. Members
+share the projection of `cols` and `obj` onto the solver's column index, which
+lives on `ModelTables` so two solvers cannot drift into loading different
+models; they never share hand-off code, because the currencies differ (HiGHS
+takes the three CSR arrays, gurobipy a matrix object) and because `gurobipy`
+must stay off the import path of a caller who does not use it.
+
 ## Module map
 
 | Module | Role |
@@ -400,8 +419,8 @@ than discovered at solve time.
 | `relational/engines/polars/executor.py` | assemble the model frames from the bound data |
 | `relational/result.py` | what a solve returned: status, objective, and the label joins that read values back |
 | `relational/engines/polars/data_validation.py` | is the bound data usable — one row per coordinate, labels that exist, single-valued coords |
-| `relational/sinks/tables.py` | what every sink reads and no more — the four frames plus the batching scalars; what an engine produces |
-| `relational/sinks/` | how a built model leaves: `lp_file`, `solver_direct` (one module each, [README](https://github.com/FBumann/lpspec/blob/main/src/lpspec/relational/sinks/README.md)) |
+| `relational/sinks/tables.py` | what every sink reads and no more — the four frames plus the batching scalars, and their projection onto the solver's column index; what an engine produces |
+| `relational/sinks/` | how a built model leaves, in two families: `solvers/` (one module per solver, chosen by name) and `writers/` (one per format, chosen by suffix) — [README](https://github.com/FBumann/lpspec/blob/main/src/lpspec/relational/sinks/README.md) |
 | `linopy/__init__.py` | opt-in shim: `build` / `extend` on a `linopy.Model` |
 | `linopy/loader.py` | data coercion to `xr.Dataset`, master coords |
 | `linopy/builder.py` | eager backend: core AST → `linopy.Model` |
@@ -497,6 +516,14 @@ read back by joining labels to coordinates.
 
 **Add a macro or named expression:** edit YAML. Nothing else.
 
+**Add a sink:** a module in `relational/sinks/solvers/` named for the solver
+(`solve_<name>`, `build_<name>`, one line in `SOLVERS`, its dependency behind an
+extra and imported inside the function), or one in `writers/` keyed by suffix in
+`WRITERS`. Nothing above it changes — no method on the executor, no branch in
+`api.py`, no name on the Python surface. The
+[README](https://github.com/FBumann/lpspec/blob/main/src/lpspec/relational/sinks/README.md)
+is the full list, and `tests/test_architecture.py` checks the shape off the path.
+
 **Add a consumer of the AST** (a renderer, a checker, a report): a directory
 beside `typeset/`, a fence test naming what it may import, and a walk. It reads
 `language.load_schema` and stops there — if it needs the plan it is a lane, not
@@ -506,8 +533,8 @@ a consumer, and the ceiling doc is the conversation to have first.
 signature in `helpers.BUILTINS` (arity and which arguments name dimensions —
 resolution, validation and lowering all read it from there, so the shape is
 declared once) → eager helper → plan node + locality class → executor →
-lowering case → differential test on both sinks → SPEC §5/§7, and this file if
-structural.
+lowering case → differential test through a solver *and* the LP writer → SPEC
+§5/§7, and this file if structural.
 
 Three things are deliberately *not* per-primitive work, because they are one
 implementation each: a primitive's dim rule lives only in `language/dimensions.py` —

docs/ROADMAP.md

@@ -69,7 +69,8 @@ already materialise.
 - **Reading results.** Duals have shipped; reduced costs and slacks ride the
   same join and have not. Derived results (LCOE, curtailment, emissions by
   group) are SQL over the solution tables.
-- **Infeasibility.** HiGHS has no IIS, so the answer is **elastic relaxation**
+- **Infeasibility.** HiGHS has no IIS — Gurobi does, but an answer only the
+  opt-in sink can give is not the answer — so it is **elastic relaxation**
   ([#80](https://github.com/FBumann/lpspec/issues/80)) — slacks with penalty
   costs, then a query of nonzero slacks grouped by block. Needs no solver
   feature and works on every sink. Taxed like a primitive, since new variables
@@ -119,7 +120,7 @@ That unblocks three things, in order of effort:
 |---|---|---|
 | **Semi-continuous** | nothing | HiGHS has `kSemiContinuous` natively and linopy has the oracle — [#383](https://github.com/FBumann/lpspec/issues/383) |
 | **SOS / indicator** | the capability model | `lp_file` carries SOS as a text section, Gurobi natively, HiGHS not at all — [#23](https://github.com/FBumann/lpspec/issues/23) |
-| **Quadratic** | the capability model, **then** a second solver | below |
+| **Quadratic** | the capability model — the second solver ships | below |
 
 **Quadratic is planned, not refused.** The cost side is settled and small, and
 the oracle is free — linopy's `QuadraticExpression` builds the comparison, which
@@ -129,9 +130,11 @@ The blocker is *where it can land*: HiGHS returns `kError` for
 `Hessian + integrality`, and `binary:`, `integer:` and nonconvex `piecewise:` all
 ship today, so on the default path quadratic conflicts with features already in
 the language. That is a capability finding, not a reason to refuse the math — so
-it needs the table above **and** a second solver on `solver_direct` without the
-exclusion (Gurobi, [#106](https://github.com/FBumann/lpspec/issues/106)).
-Landing the primitive before either ships math the default solver refuses.
+it needs the table above **and** a solver without the exclusion. The second
+solver is no longer the blocking half: the `gurobi` solver ships
+([#106](https://github.com/FBumann/lpspec/issues/106)), and Gurobi takes a
+Hessian alongside integrality. What is left is the capability table, without
+which landing the primitive would ship math the *default* solver refuses.
 
 Until then `piecewise: {convex: true}` and the epigraph pattern are the answer
 for convex 1-D curves — and they keep the LP duals, warm starts and MILP

docs/SPEC.md

@@ -502,7 +502,7 @@ language: nothing there changes what a file means.
 | Not here | Instead |
 |---|---|
 | time-series processing (resample, cluster, interpolate, align), file IO, units | data prep; pass a parameter |
-| solver breadth | HiGHS via `solver_direct`, Gurobi planned on the same path, LP files for everything else ([#106](https://github.com/FBumann/lpspec/issues/106)) |
+| solver breadth | two solver sinks — HiGHS, which ships, and Gurobi via the `[gurobi]` extra — chosen with `solver_name` at the call, never in the file; LP files for everything else ([#106](https://github.com/FBumann/lpspec/issues/106)) |
 | SOS and indicator constraints | planned, as a *sink capability* rather than a language question — the default solver has no such concept, `lp_file` and Gurobi do ([#23](https://github.com/FBumann/lpspec/issues/23), [ROADMAP Track 3](ROADMAP.md#track-3--capabilities-and-the-degree-line)). `piecewise:` (§4) covers SOS2's usual purpose today |
 | multi-objective | one objective — declaring a second is a load error (§2); weight them into one expression |
 | schema migrations | — |

docs/api.md

@@ -18,6 +18,7 @@ lps.check('model.yaml')  # parse → validate → lower, no data bound
 schema = lps.load_schema('model.yaml')  # MathSchema
 
 result = lps.solve('model.yaml', sources, solver_options={'time_limit': 60})
+# ...or solver_name='gurobi', the other solver sink — same model either way
 result.status, result.termination_condition, result.objective
 result.is_ok  # rolled-up verdict: not an error, abort or refusal
 result.has_primal  # narrower: are there values to read
@@ -39,7 +40,7 @@ executor when one build should feed more than one sink:
 
 ```python
 ex = lps.build('model.yaml', sources)
-ex.write_lp('model.lp')
+ex.write('model.lp')
 result = ex.solve()
 ```
 
@@ -50,16 +51,25 @@ pandas / xarray, which ship with the `[linopy]` extra. The only build knob is
 `coords`; **`solver_options` is not a build knob** and is forwarded verbatim to
 the solver.
 
+**Which solver is a caller's choice, not the file's.** `solver_name` is
+`highs` (ships with the package) or `gurobi` (needs the `[gurobi]` extra), and
+nothing in the YAML names one — the same file means the same model whichever
+takes it. Options travel in the chosen solver's own vocabulary,
+`{'time_limit': 60}` for HiGHS against `{'TimeLimit': 60}` for Gurobi, because
+forwarding verbatim is the contract and translating names would mean holding an
+opinion about every option either one has. A name outside the two is an error
+listing them, never a quiet fallback to the default.
+
 Reading a result:
 
 | Rule | |
 |---|---|
 | **`is_ok` is not `has_primal`** | `is_ok` rolls up the termination condition; `has_primal` adds the solver's verdict on whether an incumbent exists, and is what every reader gates on. A MIP that hits `time_limit` before finding a feasible point is `ok` with nothing to read |
 | reading anyway | `NoSolutionError`; `objective` is `nan` |
 | `dual` **raises rather than zero-filling** | no values at all is `NoSolutionError`; values but no duals — any integer or binary variable makes them undefined — is `LpspecError`, because only this quantity is missing |
-| duals exist only on `solver_direct` | a model written to LP and solved elsewhere never passes back through here. Reduced costs and slacks ride the same join and are not exposed yet |
+| duals exist only where a solver ran | either solver sink hands them back through the same join; a model written to LP and solved elsewhere never passes back through here. Reduced costs and slacks ride that join too and are not exposed yet |
 | `to_dataset` costs what it says | each variable arrives dense over its own dims — name a subset, or use `to_parquet` |
-| `write` | `.lp` only today; `.mps` raises `NotImplementedError` |
+| `write` | the **suffix** picks the writer — `.lp` today, `.mps` a `NotImplementedError` naming it as planned, anything else a `ValueError` listing both sets. Checked before the build, so a format nothing can write costs no model |
 
 **The linopy shim** (`lpspec.linopy.build` / `.extend`, `[linopy]` extra) puts
 the same YAML math on a `linopy.Model` that already exists in memory. It is

docs/index.md

@@ -274,6 +274,7 @@ Only the notation is a choice, and **How** shows the one that was made here.
 pip install lpspec  # the relational engine (polars, highspy)
 pip install "lpspec[linopy]"  # adds linopy + xarray + pandas: the shim, the
                               # oracle, and to_pandas / to_dataarray
+pip install "lpspec[gurobi]"  # adds the gurobi sink: solver_name='gurobi'
 ```
 
 !!! warning "Alpha, pre-1.0"