diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index d215e333a..05bd91273 100644 --- a/CONTRIBUTING.md +++ b/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 diff --git a/README.md b/README.md index 397cc9e23..b8eff2754 100644 --- a/README.md +++ b/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 diff --git a/bench/README.md b/bench/README.md index 913d2bb1d..245b6ab5f 100644 --- a/bench/README.md +++ b/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 diff --git a/bench/_run_case.py b/bench/_run_case.py index 03171ed73..11fec376f 100644 --- a/bench/_run_case.py +++ b/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 diff --git a/bench/profile_build.py b/bench/profile_build.py index 59fb880a0..08a5d84af 100644 --- a/bench/profile_build.py +++ b/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') diff --git a/bench/regressions/test_build.py b/bench/regressions/test_build.py index 51b401403..f08157cd2 100644 --- a/bench/regressions/test_build.py +++ b/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() diff --git a/docs/ARCHITECTURE.md b/docs/ARCHITECTURE.md index 65aa5add3..7aa485d0d 100644 --- a/docs/ARCHITECTURE.md +++ b/docs/ARCHITECTURE.md @@ -80,8 +80,8 @@ flowchart TB BIND["binding.py
→ BoundSources, frozen"] --> EXEC["executor.py + labels.py
assemble the model frames"] end ENG --> TABLES["sinks/tables.py
cols · obj · rows · A"] - TABLES --> LPS["sinks/lp_file.py
(mps planned)"] - TABLES --> DIRECT["sinks/highs.py
COO batches → HiGHS"] + TABLES --> LPS["sinks/writers/
a file, chosen by suffix
lp_file (mps planned)"] + TABLES --> DIRECT["sinks/solvers/
COO batches → the solver, chosen by name
highs (ships) · gurobi (extra)"] DIRECT --> SOL["result.py
label join, never dense"] end @@ -162,7 +162,7 @@ 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 @@ -170,9 +170,9 @@ 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 ` 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 ` 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_`, `build_`, 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` — diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index fccf9b753..de3063932 100644 --- a/docs/ROADMAP.md +++ b/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 diff --git a/docs/SPEC.md b/docs/SPEC.md index 02dd846c3..46b33f8da 100644 --- a/docs/SPEC.md +++ b/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 | — | diff --git a/docs/api.md b/docs/api.md index d7441cd70..55a9bdd8a 100644 --- a/docs/api.md +++ b/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,6 +51,15 @@ 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 | | @@ -57,9 +67,9 @@ Reading a result: | **`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 diff --git a/docs/index.md b/docs/index.md index a046d6dc2..1d31db9d1 100644 --- a/docs/index.md +++ b/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" diff --git a/examples/walkthrough.py b/examples/walkthrough.py index 8dc5d2a02..f794ca44d 100644 --- a/examples/walkthrough.py +++ b/examples/walkthrough.py @@ -143,11 +143,11 @@ def main() -> None: # -------------------------------------------------------------- banner(5, 'sink: stream the frames to an LP file', 'relational/sinks/lp_file.py') - # Same frames, second sink. The other one (solver_direct, stage 6) + # Same frames, second sink. The other one (the `highs` solver, stage 6) # hands COO batches to highspy without ever forming a full CSR here. with tempfile.TemporaryDirectory() as tmp: lp = Path(tmp) / 'model.lp' - ex.write_lp(lp) + ex.write(lp) text = lp.read_text().splitlines() print(_indent('\n'.join(text[:12]))) print(f' ... ({len(text)} lines total)') diff --git a/pyproject.toml b/pyproject.toml index 33e5f70ac..c50b7a9db 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -73,6 +73,21 @@ linopy = [ # solves and reads results back without it. "pandas>=2.2", ] +# The `gurobi` sink (relational/sinks/gurobi.py). Opt-in, and the only sink +# with a licence question: gurobipy's own wheel carries a size-limited licence +# that solves small models out of the box, which is what the tests run on. +# +# scipy is not a convenience here. gurobipy's matrix API takes a *matrix +# object* where HiGHS takes the three CSR arrays, so `scipy.sparse.csr_matrix` +# is the only way to hand over a block of rows without one Python call per +# row. Only the long-stable `csr_matrix` constructor is used, so its floor is +# the oldest release we would test against rather than a feature; gurobipy's +# is 11, the release whose status codes the sink's condition map is written +# against — 17 is `MEM_LIMIT` from that version on. +gurobi = [ + "gurobipy>=11", + "scipy>=1.11", +] [project.urls] repository = "https://github.com/FBumann/lpspec" @@ -119,6 +134,11 @@ dev = [ # the shim and the oracle, by reference: the pins live in the extra above, # so they cannot drift from what a user installing `lpspec[linopy]` gets "lpspec[linopy]", + # the second solver sink, by reference for the same reason. In `dev` so CI + # exercises it: gurobipy's wheel carries a size-limited licence, and the + # sink's tests stay under it. Its own tests skip where it is absent, so a + # contributor who does not want it can drop it without failing the suite. + "lpspec[gurobi]", # how pandas reads parquet, which the benchmark harness's linopy arm does # as part of the cost it is there to measure. Not needed by the engine — # the bare-install job runs the suite without it. diff --git a/src/lpspec/api.py b/src/lpspec/api.py index 4594206f9..dcda84868 100644 --- a/src/lpspec/api.py +++ b/src/lpspec/api.py @@ -34,6 +34,7 @@ from lpspec.language.validation import load_schema from lpspec.lowering import lower_program from lpspec.relational.engines.polars.executor import PolarsExecutor +from lpspec.relational.sinks import solver, writer from lpspec.sources import tidy_sources if TYPE_CHECKING: @@ -71,7 +72,7 @@ def build( ``sources`` maps parameter names to parquet paths or in-memory tables (and optionally dimension names to index tables). One build can feed more than - one sink: call ``ex.solve()`` and ``ex.write_lp(path)`` on the same object. + one sink: call ``ex.solve()`` and ``ex.write(path)`` on the same object. Raises ------ @@ -94,22 +95,30 @@ def solve( model: str | Path | dict[str, Any] | MathSchema, sources: Mapping[str, Any], solver_options: Mapping[str, Any] | None = None, + solver_name: str = 'highs', **build_kwargs: Any, ) -> Result: """Build and solve in one call. - ``solver_options`` is forwarded verbatim to the solver — the same shape - linopy takes, e.g. ``{'time_limit': 60, 'mip_rel_gap': 0.01}``. Build - options stay separate, because they govern *construction* and never reach - the solver. + ``solver_name`` is which solver sink to hand the built model to — + ``highs``, which ships with the package, or ``gurobi``, which needs the + ``[gurobi]`` extra. linopy's spelling, and a decision the *caller* makes: + the same file solves the same model either way, so nothing in the YAML + names a solver. + + ``solver_options`` is forwarded verbatim to it — the same shape linopy + takes, e.g. ``{'time_limit': 60, 'mip_rel_gap': 0.01}``, in whichever + solver's vocabulary was chosen. Build options stay separate, because they + govern *construction* and never reach the solver. The executor stays attached to the returned :class:`Result`, whose label frames back ``result.primal(...)``. Nothing has to be released, though ``result.close()`` drops a large model early if you want the memory back. """ + solver(solver_name) # before the build, for the reason `write` checks the suffix first ex = build(model, sources, **build_kwargs) try: - return ex.solve(solver_options=solver_options) + return ex.solve(solver_options=solver_options, solver_name=solver_name) except BaseException: ex.close() raise @@ -123,17 +132,15 @@ def write( ) -> Path: """Build and stream the model to a file; format from the suffix. - ``.lp`` is supported today; ``.mps`` is planned (same streaming - mechanics, see docs/ARCHITECTURE.md sinks). + ``.lp`` is supported today and ``.mps`` is planned, both answered by the + writer family rather than by a branch here — this verb owns *when* to + build, not what can be written. + + The suffix is checked **before** the build, because a caller who named a + format nothing can write should not pay for a model first. """ out = Path(out) - suffix = out.suffix.lower() - if suffix == '.lp': - with build(model, sources, **build_kwargs) as ex: - ex.write_lp(out) - return out - if suffix == '.mps': - msg = 'the mps sink is planned but not implemented yet (docs/ARCHITECTURE.md, sinks)' - raise NotImplementedError(msg) - msg = f"unsupported output format '{suffix}' — supported: .lp (planned: .mps)" - raise ValueError(msg) + writer(out.suffix.lower()) + with build(model, sources, **build_kwargs) as ex: + ex.write(out) + return out diff --git a/src/lpspec/relational/__init__.py b/src/lpspec/relational/__init__.py index 339e30ed6..67037532a 100644 --- a/src/lpspec/relational/__init__.py +++ b/src/lpspec/relational/__init__.py @@ -21,16 +21,12 @@ ``lpspec.relational.plan`` and are imported from there — so adding a node no longer silently widens something that reads like public API, and the import site says which layer the caller is reaching into. - -``RelationalBuildError`` is kept as a deprecated alias only. The engine now -raises ``lpspec.errors.LanguageError`` and ``DataError``; catch those. """ -from lpspec.relational.engines.polars import PolarsExecutor, RelationalBuildError +from lpspec.relational.engines.polars import PolarsExecutor from lpspec.relational.result import Result __all__ = [ 'PolarsExecutor', - 'RelationalBuildError', 'Result', ] diff --git a/src/lpspec/relational/chunking.py b/src/lpspec/relational/chunking.py index ed2a42403..b0f1a616b 100644 --- a/src/lpspec/relational/chunking.py +++ b/src/lpspec/relational/chunking.py @@ -8,7 +8,7 @@ costs one. The width is the part that gets forgotten, and forgetting it does not look -like a bug. ``solver_direct`` chunked the constraint matrix by rows with no +like a bug. The HiGHS solver chunked the constraint matrix by rows with no width at all, which reads as bounded and is not: a row is nine entries in one model and a hundred in another, so what the pass held tracked the model's shape rather than the budget — the one thing a *batched* pass exists to stop diff --git a/src/lpspec/relational/engines/polars/__init__.py b/src/lpspec/relational/engines/polars/__init__.py index 8309c786e..731a2633c 100644 --- a/src/lpspec/relational/engines/polars/__init__.py +++ b/src/lpspec/relational/engines/polars/__init__.py @@ -10,6 +10,6 @@ which of the eleven modules under `relational/` were which. """ -from lpspec.relational.engines.polars.executor import PolarsExecutor, RelationalBuildError +from lpspec.relational.engines.polars.executor import PolarsExecutor -__all__ = ['PolarsExecutor', 'RelationalBuildError'] +__all__ = ['PolarsExecutor'] diff --git a/src/lpspec/relational/engines/polars/executor.py b/src/lpspec/relational/engines/polars/executor.py index 5bb942bf4..30867225d 100644 --- a/src/lpspec/relational/engines/polars/executor.py +++ b/src/lpspec/relational/engines/polars/executor.py @@ -16,6 +16,7 @@ from __future__ import annotations +from pathlib import Path from typing import TYPE_CHECKING, Any, Literal, get_args import polars as pl @@ -23,7 +24,6 @@ from lpspec.errors import ( DataError, LanguageError, - LpspecError, null_bounds_message, sparse_divisor_message, ) @@ -35,7 +35,6 @@ if TYPE_CHECKING: from collections.abc import Mapping, Sequence - from pathlib import Path #: The four frames a sink reads, as schemas. Stated here because the executor @@ -61,14 +60,6 @@ } # fmt: skip -#: Deprecated. The engine's failures are now split between -#: :class:`~lpspec.errors.LanguageError` (the program says something the -#: engine cannot build) and :class:`~lpspec.errors.DataError` (a source is -#: missing or the wrong shape). This alias is their common base, so an existing -#: ``except RelationalBuildError`` keeps catching everything it used to. -RelationalBuildError = LpspecError - - class PolarsExecutor: """Build a :class:`Program` into polars frames, then sink it.""" @@ -309,23 +300,39 @@ def _tables(self) -> sinks.ModelTables: objective_constant=self._obj_const, ) - def write_lp(self, path: str | Path) -> None: - """Sink the built model to an LP file.""" - sinks.write_lp_file(self._tables(), path) + def write(self, path: str | Path) -> None: + """Sink the built model to a file; the **suffix** picks the writer. + + ``.lp`` today, ``.mps`` planned — an unknown suffix is an error naming + both sets. The caller names an output rather than a writer, which is + the one place this differs from :meth:`solve`: a file's format is a + property of the file, where which solver runs is not a property of + anything but the call. + """ + path = Path(path) + sinks.writer(path.suffix.lower())(self._tables(), path) def solve( self, batch_rows: int | None = None, solver_options: Mapping[str, Any] | None = None, + solver_name: str = 'highs', ) -> Result: - """Sink the built model straight into HiGHS and solve it. - - ``solver_options`` is forwarded verbatim to the solver, the way - linopy's is — ``{'time_limit': 60, 'mip_rel_gap': 0.01}``. - ``batch_rows`` is the hand-off budget in elements, and defaults to the - sink's own — see :data:`~lpspec.relational.sinks.highs.HANDOFF_BUDGET`. + """Sink the built model straight into a solver and solve it. + + ``solver_name`` picks the sink — ``highs``, which ships with the + package, or ``gurobi``, which needs the ``[gurobi]`` extra. Spelled + the way linopy spells it, and a *caller's* choice at the call: no YAML + file can express it, because a model means the same thing whoever + solves it. + + ``solver_options`` is forwarded verbatim to that solver, the way + linopy's is — ``{'time_limit': 60, 'mip_rel_gap': 0.01}``, and so + named in the solver's own vocabulary. ``batch_rows`` is the hand-off + budget in elements, and defaults to the sink's own — see + :data:`~lpspec.relational.sinks.highs.HANDOFF_BUDGET`. """ - status, objective, primal, dual = sinks.solve_direct(self._tables(), batch_rows, solver_options) + status, objective, primal, dual = sinks.solver(solver_name)(self._tables(), batch_rows, solver_options) return Result( _status=status, _objective=objective, diff --git a/src/lpspec/relational/sinks/README.md b/src/lpspec/relational/sinks/README.md index b3bab38a6..afd91648b 100644 --- a/src/lpspec/relational/sinks/README.md +++ b/src/lpspec/relational/sinks/README.md @@ -1,13 +1,23 @@ # Sinks -How a built model leaves the engine. These are the two boxes downstream of the -executor in [docs/ARCHITECTURE.md](../../../../docs/ARCHITECTURE.md)'s pipeline. +How a built model leaves the engine — the boxes downstream of the executor in +[docs/ARCHITECTURE.md](../../../../docs/ARCHITECTURE.md)'s pipeline, which +carries the argument for the split. This page is the membership list. -| Module | Sink | Needs | +**Two families.** A **solver** takes the tables and runs them; a **writer** +takes the tables and renders them to a file. Everything else follows. + +| | solvers/ | writers/ | |---|---|---| -| `tables.py` | — | the contract every sink reads | -| `lp_file.py` | `lp_file` — LP text | nothing beyond polars | -| `highs.py` | `solver_direct` — COO batches → HiGHS | `highspy` | +| answers | `(tables, batch_rows, options) -> (status, objective, primal, dual)` | `(tables, path) -> None` | +| chosen by | **name**, at the call — `solver_name='gurobi'` | **suffix**, from the output — `model.lp` | +| registry | `SOLVERS`, closed | `WRITERS` + `PLANNED_WRITERS`, closed | +| members | `highs.py` (`highspy`, ships), `gurobi.py` (`[gurobi]`: `gurobipy`, `scipy`) | `lp_file.py` (nothing beyond polars) | + +`tables.py` is what both read. Neither family imports the other and no member +imports a sibling — `tests/test_architecture.py` reads all of that off the +path, which is what keeps `gurobipy` off the import path of a caller who +solves with HiGHS. ## The contract @@ -18,46 +28,37 @@ constant — those last two live outside the tables because a constant has no column to attach to. A sink never learns how the tables were filled, and the executor never learns -how they are drained. That is the point: adding `mps` is a new module here, not -another method on `PolarsExecutor`. - -## Adding a sink - -1. New module, named for the sink. Take `ModelTables`, return whatever the sink - naturally returns (`None` for a writer, `(status, objective)` for a solver). -2. Re-export it from `__init__.py`. -3. Thin delegation on `PolarsExecutor` — three lines, no logic. -4. If it needs an optional dependency, import it **inside the function**. The - module boundary is the fence; the lazy import is what keeps importing this - package free for callers who will never use that sink. -5. Stream. Nothing here may materialise the model a second time. Sink a lazy - frame, or hand the solver slices of what the labels already laid out. - -## Why one module per sink - -Because that is where the fences are. `highspy` is an optional dependency of -`solver_direct` alone, and a caller that only writes LP files should not import -it. Splitting by *kind* — all writers together, all solvers together — would -put two optional imports in one module and a function that branches on which -solver you meant. - -When `mps` lands it may well belong beside `lp_file` (both sink rendered lines -into part files and concatenate them bytewise; they would share `_cat` and the -float formatting). That is a decision to take with -the code in hand, not now — a `text.py` holding one function today would be a -guess about a sink that does not exist. - -## When Track 4 lands - -[Track 4](../../../../docs/ROADMAP.md#track-4--sink-capabilities) gives each sink a -declared capability table so `check(model, sink=...)` can answer "will this -sink take it". Two notes for whoever writes it: - -- The capability table belongs next to the sink it describes, in that sink's - module. `__init__.py` collects them; it does not own them. -- Make the set of sinks **closed**, like `helpers.BUILTINS` — not an open - `register_sink()`. An installed plugin that can change the answer to - `lps.check(model, sink=...)` is hard rule 5's failure mode one level down. +how they are drained. That is the point: adding `mps` is a new module in +`writers/`, not another method on `PolarsExecutor`. + +The one thing sinks may share is a *projection* of those frames, never a step +of the work — `ModelTables.dense_columns`, which both solvers read. + +## Adding one + +**A solver:** `solvers/.py` named for the solver, defining `solve_` +and the `build_` seam `bench/` measures, plus one line in `SOLVERS`. +Import the solver **inside the function** and declare an extra for it — the +module boundary is the fence, the lazy import is what keeps this package free +to import for callers who will never use it. Copy linopy's status map for it +and pin the copy in `tests/test_solve_status.py`, including anywhere you +deliberately diverge. + +**A writer:** `writers/.py`, one line in `WRITERS` keyed by suffix, +moved out of `PLANNED_WRITERS` if it was there. + +Either way: stream — nothing here may materialise the model a second time — +and nothing above changes. No method on the executor, no branch in `api.py`, +no name on the Python surface. + +## When Track 3 lands + +[Track 3](../../../../docs/ROADMAP.md#track-3--capabilities-and-the-degree-line) +gives each sink a declared capability table so `check(model, sink=...)` can +answer "will this sink take it". The table belongs in the sink's own module, +collected by the family `__init__` rather than owned by it — and it stops +being uniform at exactly the seam this directory draws: SOS is native in +`gurobi`, a text section in `lp_file`, absent in `highs`. ## Stable output @@ -66,5 +67,6 @@ Two runs of one model produce the same bytes is easy to lose: a parallel join hands back a group in whatever order it finished it, so a sink that gathers a row's terms and *then* orders the rows has already lost the order within one. `lp_file` emits one frame of lines -carrying its own sort key instead, and sorts once. `solver_direct` is ordered -for a different reason — `searchsorted` requires it. +carrying its own sort key instead, and sorts once. The solvers are ordered for +a different reason — `searchsorted` requires it, and the CSR `indptr` it +produces is only a row's extent if the rows it indexes are sorted. diff --git a/src/lpspec/relational/sinks/__init__.py b/src/lpspec/relational/sinks/__init__.py index 4b446615d..f34aa8751 100644 --- a/src/lpspec/relational/sinks/__init__.py +++ b/src/lpspec/relational/sinks/__init__.py @@ -1,23 +1,24 @@ -"""Sinks: how a built model leaves the engine. +"""Sinks: how a built model leaves the engine. See README.md. -docs/ARCHITECTURE.md's pipeline draws two boxes downstream of the executor. This -package is those boxes — one module each, because that is where the fences -are: ``highspy`` is an optional dependency of ``solver_direct`` alone, and a -caller that only writes LP files should never import it. +**Two families, and that is the whole mental model.** A *solver* takes the +tables and runs them (``solvers/``, chosen by name); a *writer* takes the +tables and renders them to a file (``writers/``, chosen by suffix). The +families are directories rather than a convention, so +``tests/test_architecture.py`` reads membership off the path. -A sink reads :class:`ModelTables` and nothing else. Neither sink knows how the -tables were filled, and the executor does not know how they are drained, which -is what makes the planned ``mps`` sink a new module here rather than another -method on the executor. +``tables.py`` is what both read, and neither family imports the other. This +module is the seam a caller uses: the contract, and the two lookups. """ -from lpspec.relational.sinks.highs import build_highs, solve_direct -from lpspec.relational.sinks.lp_file import write_lp_file +from lpspec.relational.sinks.solvers import SOLVERS, solver from lpspec.relational.sinks.tables import ModelTables +from lpspec.relational.sinks.writers import PLANNED_WRITERS, WRITERS, writer __all__ = [ + 'PLANNED_WRITERS', + 'SOLVERS', + 'WRITERS', 'ModelTables', - 'build_highs', - 'solve_direct', - 'write_lp_file', + 'solver', + 'writer', ] diff --git a/src/lpspec/relational/sinks/solvers/__init__.py b/src/lpspec/relational/sinks/solvers/__init__.py new file mode 100644 index 000000000..2213c6e4e --- /dev/null +++ b/src/lpspec/relational/sinks/solvers/__init__.py @@ -0,0 +1,54 @@ +"""The solver family: the tables in, an answer out. See ../README.md. + +One module per solver, **named for the solver** — nothing here is named for +the mechanism, because every member uses the same one. Each answers:: + + (tables, batch_rows, solver_options) -> (status, objective, primal, dual) + +plus a ``build_`` that loads the model and stops, which is the seam +`bench/` measures. ``tests/test_architecture.py`` checks that off the path. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from lpspec.errors import LpspecError +from lpspec.relational.sinks.solvers.gurobi import solve_gurobi +from lpspec.relational.sinks.solvers.highs import solve_highs + +if TYPE_CHECKING: + from collections.abc import Callable, Mapping + from typing import Any + + import polars as pl + + from lpspec.relational.sinks.tables import ModelTables + from lpspec.relational.status import SolveStatus + + Solve = Callable[ + [ModelTables, int | None, Mapping[str, Any] | None], + tuple[SolveStatus, float, pl.DataFrame | None, pl.DataFrame | None], + ] + +__all__ = ['SOLVERS', 'solver'] + +#: Every solver a caller may name, and **closed** — a dict literal, not a +#: registry something installed can add to. Which solver runs is the caller's +#: choice at the call and never the file's, so an installed package that could +#: change what ``solver_name='x'`` resolves to is hard rule 5 one level down. +SOLVERS: Mapping[str, Solve] = { + 'highs': solve_highs, + 'gurobi': solve_gurobi, +} + + +def solver(name: str) -> Solve: + """The solver called *name*, or an error listing every alternative.""" + try: + return SOLVERS[name] + except KeyError: + raise LpspecError( + f'unknown solver {name!r} — this build solves with {", ".join(sorted(SOLVERS))}. ' + 'HiGHS ships with the package and is the default; gurobi needs the [gurobi] extra.' + ) from None diff --git a/src/lpspec/relational/sinks/solvers/gurobi.py b/src/lpspec/relational/sinks/solvers/gurobi.py new file mode 100644 index 000000000..d191d75a6 --- /dev/null +++ b/src/lpspec/relational/sinks/solvers/gurobi.py @@ -0,0 +1,253 @@ +"""The ``gurobi`` solver: COO blocks straight into gurobipy. + +The same hand-off as :mod:`~lpspec.relational.sinks.solvers.highs`, reading +the same ``dense_columns``, so the two cannot disagree about the model they +load. Two things differ: + +- **The matrix's currency.** HiGHS takes the three CSR arrays; gurobipy's + matrix API takes a matrix *object*, so a block is wrapped in a + ``scipy.sparse.csr_matrix`` over them — a view, not a copy. That wrapper is + why the ``[gurobi]`` extra carries scipy: the alternative is a Python call + per row. +- **Columns arrive in one call.** ``addMConstr`` writes into one ``MVar`` + spanning the model, so only the matrix is chunked. + +``gurobipy`` and ``scipy`` are imported inside the functions, so importing +this module stays free for a caller who never solves with it. +""" + +from __future__ import annotations + +import weakref +from typing import TYPE_CHECKING, Any + +import polars as pl + +from lpspec.errors import LpspecError +from lpspec.relational.status import SolveStatus + +if TYPE_CHECKING: + from collections.abc import Mapping + + from lpspec.relational.sinks.tables import ModelTables + + +#: Nonzeros per row block, spent through :mod:`~lpspec.relational.chunking`. +#: Not independently tuned — the measurement is +#: :data:`~lpspec.relational.sinks.solvers.highs.HANDOFF_BUDGET`'s, and holds +#: because a block is the same filter over the same sorted frame; wrapping +#: three arrays in a CSR view is constant work on top. +HANDOFF_BUDGET = 100_000 + +#: Gurobi status -> termination condition. Copied from linopy's own +#: ``Gurobi.CONDITION_MAP`` bar three entries (:data:`_LINOPY_DIVERGENCES`); +#: ``tests/test_solve_status.py`` asserts both halves, so linopy moving — or +#: fixing — fails here rather than silently. +_CONDITION_OF_GUROBI_STATUS = { + 1: 'unknown', + 2: 'optimal', + 3: 'infeasible', + 4: 'infeasible_or_unbounded', + 5: 'unbounded', + 6: 'other', + 7: 'iteration_limit', + 8: 'terminated_by_limit', + 9: 'time_limit', + 10: 'terminated_by_limit', + 11: 'user_interrupt', + 12: 'other', + 13: 'suboptimal', + 14: 'unknown', + 15: 'terminated_by_limit', + 16: 'terminated_by_limit', + 17: 'resource_interrupt', +} + +#: Where the table above does not copy linopy's, and why: each contradicts a +#: status Gurobi documents, so copying it would import a wrong answer rather +#: than a shared vocabulary — the trade +#: :attr:`~lpspec.relational.status.SolveStatus.is_readable` already refused +#: once. The words stay linopy's; only the verdicts differ. +_LINOPY_DIVERGENCES = { + 10: 'SOLUTION_LIMIT stopped early after n incumbents; linopy calls it optimal', + 16: 'WORK_LIMIT is a limit, not a solver failure; linopy calls it internal_solver_error', + 17: 'MEM_LIMIT is the resource_interrupt linopy itself maps kMemoryLimit to on HiGHS', +} + + +def build_gurobi( + model: ModelTables, + batch_rows: int | None = None, + solver_options: Mapping[str, Any] | None = None, +) -> Any: + """Load the model into a :class:`gurobipy.Model` and stop there. + + :func:`~lpspec.relational.sinks.solvers.highs.build_highs`'s seam, drawn + for its reason: the search is the same work whoever filled the model. + ``batch_rows`` is the budget in *nonzeros*, and stays a parameter so tests + can force ragged blocks. + + **The caller owns the model, so the environment follows it.** gurobipy has + no ``Model.getEnv()``, so a caller handed only the model could never + release the licence it holds; a finalizer disposes the environment when + the model is collected, which under refcounting is when the caller drops + it. That is one thing to own rather than two, and it is why this returns a + model rather than a pair. (linopy's equivalent — a solver model detached + from its `Solver` — leaves the environment to gurobipy to free whenever.) + """ + m, _x, _blocks, environment = _load(model, batch_rows, solver_options) + weakref.finalize(m, environment.dispose) + return m + + +def solve_gurobi( + model: ModelTables, + batch_rows: int | None = None, + solver_options: Mapping[str, Any] | None = None, +) -> tuple[SolveStatus, float, pl.DataFrame | None, pl.DataFrame | None]: + """Feed the model to Gurobi and solve it. + + The family's shape, and the two ``None`` cases mean what they mean in + :func:`~lpspec.relational.sinks.solvers.highs.solve_highs`: no primal, no + dual. Gurobi refuses the attribute in both cases rather than handing back + zeros, which is the one place it makes this easier than HiGHS. + + Model and environment are disposed before returning, so the licence is + released when the solve ends rather than whenever the collector gets to + it. Everything read back is a numpy array taken before the dispose. + :func:`build_gurobi` cannot do this — its caller owns the model. + """ + m, x, blocks, environment = _load(model, batch_rows, solver_options) + try: + m.optimize() + status = _status_of(m) + if not status.is_readable: + return status, float('nan'), None, None + return status, m.ObjVal, _labelled('col', x.X), _duals(model.row_count, blocks) + finally: + m.dispose() + environment.dispose() + + +def _load( + model: ModelTables, + batch_rows: int | None, + solver_options: Mapping[str, Any] | None, +) -> tuple[Any, Any, list[Any], Any]: + """The model, the handles to read it back, and the environment to release. + + ``x.X`` and ``block.Pi`` are numpy arrays; ``getVars()`` / ``getConstrs()`` + would build one Python object per column and row to reach the same numbers. + The environment comes back because gurobipy has no ``Model.getEnv()``, and + whoever disposes the model has to dispose it too. + + **Options go on the environment, not the model.** A licence parameter — + ``WLSAccessID``, ``ComputeServer``, ``TokenServer`` — can only be set + before an environment starts, and ``setParam`` on the model refuses it + with *unable to modify parameter after environment started*. So a + Compute-Server or WLS user could not reach this sink at all. Everything + else is unaffected: an environment's parameters are the defaults of every + model built on it, and an unknown name still raises at the same point. + ``OutputFlag`` leads so a caller can put the log back by passing their own. + """ + gurobipy = _gurobipy() + import numpy as np + import scipy.sparse + + batch = HANDOFF_BUDGET if batch_rows is None else batch_rows + environment = gurobipy.Env(params={'OutputFlag': 0, **dict(solver_options or {})}) + m = gurobipy.Model(env=environment) + + lb, ub, cost, integral = model.dense_columns(gurobipy.GRB.INFINITY) + x = m.addMVar(model.column_count, lb=lb, ub=ub, obj=cost, vtype=np.where(integral, 'I', 'C')) + + ordered_rows = model.rows.sort('row') + ordered_matrix = model.matrix.sort('row') + senses = np.array([gurobipy.GRB.LESS_EQUAL, gurobipy.GRB.GREATER_EQUAL, gurobipy.GRB.EQUAL], dtype='=': 1, '==': 2}, return_dtype=pl.UInt8).alias('op'), + 'rhs', + ) + a = ordered_matrix.filter(pl.col('row').is_between(lo, hi, closed='left')) + starts = np.searchsorted(a['row'].to_numpy(), rows['row'].to_numpy()) + block = scipy.sparse.csr_matrix( + (a['coeff'].to_numpy(), a['col'].to_numpy(), np.append(starts, a.height)), + shape=(rows.height, model.column_count), + ) + blocks.append(m.addMConstr(block, x, senses[rows['op'].to_numpy()], rows['rhs'].to_numpy())) + + if model.objective_sense == 'max': + m.ModelSense = gurobipy.GRB.MAXIMIZE + m.ObjCon = model.objective_constant + m.update() + return m, x, blocks, environment + + +def _gurobipy() -> Any: + """The optional dependency, or a message naming the extra. Both halves are + named, because the missing one is as often scipy.""" + try: + import gurobipy + import scipy.sparse # noqa: F401 — guarded here so the message covers it + except ModuleNotFoundError as exc: + msg = 'The gurobi sink requires the [gurobi] extra (gurobipy, scipy): pip install "lpspec[gurobi]"' + raise ModuleNotFoundError(msg) from exc + return gurobipy + + +def _status_of(m: Any) -> SolveStatus: + """What the solve concluded, on both axes. ``SolCount`` answers "is there + anything here", which the termination condition does not: a run stopped at + a limit may or may not hold an incumbent.""" + code = int(m.Status) + return SolveStatus( + termination_condition=_CONDITION_OF_GUROBI_STATUS.get(code, 'unknown'), + solver_wording=_wording(code), + has_primal=m.SolCount > 0, + ) + + +def _wording(code: int) -> str: + """Gurobi's own name for a status code, read off ``GRB.Status`` rather than + tabulated — so one this package has never heard of still arrives + searchable.""" + gurobipy = _gurobipy() + names = {getattr(gurobipy.GRB.Status, name): name for name in dir(gurobipy.GRB.Status) if not name.startswith('_')} + return names.get(code, str(code)) + + +def _duals(row_count: int, blocks: list[Any]) -> pl.DataFrame | None: + """Shadow prices in row order, or ``None`` where the model has none. + + Blocks were added in ascending row ranges, so concatenating their slices + reproduces the row index without a sort. Gurobi refuses ``Pi`` on a + mixed-integer model, and that refusal *is* the answer — no zero vector to + test. + """ + import numpy as np + + gurobipy = _gurobipy() + try: + slices = [block.Pi for block in blocks] + except (AttributeError, gurobipy.GurobiError): + return None + values = np.concatenate(slices) if slices else np.empty(0, dtype=np.float64) + if len(values) != row_count: + raise LpspecError( + f'the solver returned {len(values)} duals for {row_count} rows. The join that reads ' + f'them back is positional, so a short vector would drop rows silently. This is an ' + f'engine bug rather than a problem with the model — please report it.' + ) + return _labelled('row', values) + + +def _labelled(label: str, values: Any) -> pl.DataFrame: + """``(label, value)`` over the solver's own dense index, which *is* our + label — so the join column is an ``arange``, not something read back.""" + import numpy as np + + values = np.asarray(values, dtype=np.float64) + return pl.DataFrame({label: np.arange(len(values), dtype=np.int64), 'value': values}) diff --git a/src/lpspec/relational/sinks/highs.py b/src/lpspec/relational/sinks/solvers/highs.py similarity index 81% rename from src/lpspec/relational/sinks/highs.py rename to src/lpspec/relational/sinks/solvers/highs.py index 4165c687e..cd8e86e6e 100644 --- a/src/lpspec/relational/sinks/highs.py +++ b/src/lpspec/relational/sinks/solvers/highs.py @@ -1,14 +1,15 @@ -"""The ``solver_direct`` sink: COO batches straight into HiGHS. +"""The ``highs`` solver: COO batches straight into HiGHS. -No float→text→parse round trip — that is why this exists beside -:mod:`~lpspec.relational.sinks.lp_file`. Columns and rows arrive as numpy -slices, in batches. +The default, and the only one whose dependency ships with the package. No +float→text→parse round trip — that is why this exists beside +:mod:`~lpspec.relational.sinks.writers.lp_file`. Columns and rows arrive as +numpy slices, in batches. **Nothing textual crosses into numpy.** A polars ``String`` column converts by -boxing every value as a Python object, so a comparison against ``'continuous'`` -or ``'<='`` is made in polars and only its answer — a bool, or the bound it -selects — is handed over. At 10M columns the same test costs 0.95 s across the -boundary and 0.04 s on this side of it. +boxing every value as a Python object, so a comparison against ``'<='`` is made +in polars and only the bound it selects is handed over — the same rule +:meth:`~lpspec.relational.sinks.tables.ModelTables.dense_columns` states for +the column vectors, where it was measured. ``highspy`` is imported inside the function: it is an optional dependency, and importing this module must stay free for callers that only write LP files. @@ -81,7 +82,7 @@ def build_highs( ) -> Any: """Load the model into a :class:`highspy.Highs` and stop there. - The hand-off without the simplex. :func:`solve_direct` is this plus + The hand-off without the simplex. :func:`solve_highs` is this plus ``run()``, and the seam exists because the two are different questions: the simplex is the same work whoever filled the model, so a measurement that includes it says nothing about the lane that filled it. `bench/` ends here, @@ -103,14 +104,7 @@ def build_highs( empty_i = np.empty(0, dtype=np.int32) empty_f = np.empty(0, dtype=np.float64) - count = model.column_count - at = model.cols['col'].to_numpy() - lb = _by_column(count, at, model.cols['lb'].to_numpy(), -inf) - ub = _by_column(count, at, model.cols['ub'].to_numpy(), inf) - integral = _by_column(count, at, model.cols.select(pl.col('vtype') != 'continuous').to_series().to_numpy(), False) - cost = _by_column(count, model.obj['col'].to_numpy(), model.obj['coeff'].to_numpy(), 0.0) - np.nan_to_num(lb, copy=False, neginf=-inf, posinf=inf) - np.nan_to_num(ub, copy=False, neginf=-inf, posinf=inf) + lb, ub, cost, integral = model.dense_columns(inf) for lo, hi in model.col_chunks(batch): _loaded(h, h.addCols(hi - lo, cost[lo:hi], lb[lo:hi], ub[lo:hi], 0, empty_i, empty_i, empty_f), 'columns') noncontinuous = np.flatnonzero(integral[lo:hi]).astype(np.int32) + np.int32(lo) @@ -149,7 +143,7 @@ def build_highs( return h -def solve_direct( +def solve_highs( model: ModelTables, batch_rows: int | None = None, solver_options: Mapping[str, Any] | None = None, @@ -181,26 +175,6 @@ def solve_direct( return status, objective, primal, dual -def _by_column(count: int, at: Any, values: Any, absent: Any) -> Any: - """*values* written at the column each one belongs to, *absent* elsewhere. - - ``col`` is dense ``0..n-1`` (:class:`ModelTables`), so it *is* the position - a value has to end up at: lining a frame up with the solver's index is a - scatter, and neither the join that fills the objective's gaps nor the sort - that puts the bounds in order has anything to do that this does not. The - frame those two produced had to be collected whole before the first batch - could be handed over, which cost more than the model does. - - A column the table somehow has no row for is left free rather than left - holding whatever the allocator returned. - """ - import numpy as np - - dense = np.full(count, absent, dtype=values.dtype) - dense[at] = values - return dense - - def _loaded(h: Any, status: Any, what: str) -> None: """Raise unless the solver accepted the batch. diff --git a/src/lpspec/relational/sinks/tables.py b/src/lpspec/relational/sinks/tables.py index f43281d9c..8bf0d80dd 100644 --- a/src/lpspec/relational/sinks/tables.py +++ b/src/lpspec/relational/sinks/tables.py @@ -2,19 +2,33 @@ Four frames plus the scalars a writer needs to size its batching. A sink that needs a fifth thing states it here, where both sides can see it. + +Also the one *projection* of them more than one sink needs +(:meth:`ModelTables.dense_columns`). It belongs to the contract rather than to +either solver, because two sinks computing it separately could disagree about +the model they loaded — the one thing neither may do. """ from __future__ import annotations from dataclasses import dataclass -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Any + +import polars as pl from lpspec.relational import chunking if TYPE_CHECKING: from collections.abc import Iterator - import polars as pl + import numpy as np + import numpy.typing as npt + + #: What a solver sink is handed: three float vectors and an integrality + #: mask, each as long as the model has columns. + DenseColumns = tuple[ + npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.bool_] + ] @dataclass(frozen=True) @@ -58,3 +72,46 @@ def col_chunks(self, budget: int) -> Iterator[tuple[int, int]]: :mod:`~lpspec.relational.chunking` asks for. """ return chunking.ranges(self.column_count, budget, 1.0) + + def dense_columns(self, infinity: float) -> DenseColumns: + """``(lb, ub, cost, integral)`` as numpy vectors over the solver's index. + + *infinity* is the solver's own spelling of an absent bound — the one + thing the two disagree on — so it is asked for rather than assumed and + the vectors come back ready to hand over unedited. + + ``col`` is dense ``0..n-1``, so it *is* the position a value has to end + up at: lining a frame up with the solver's index is a scatter, and + neither the join that fills the objective's gaps nor the sort that + orders the bounds does anything this does not. That pair had to be + collected whole before the first batch could go over, which cost more + than the model does. A column with no row is left free rather than + holding whatever the allocator returned. + + **Nothing textual crosses into numpy.** A polars ``String`` converts by + boxing every value as a Python object, so the test against + ``'continuous'`` is made in polars and only its answer crosses: 0.04 s + against 0.95 s at 10M columns. + """ + import numpy as np + + count = self.column_count + at = self.cols['col'].to_numpy() + lb = _scattered(count, at, self.cols['lb'].to_numpy(), -infinity) + ub = _scattered(count, at, self.cols['ub'].to_numpy(), infinity) + integral = _scattered( + count, at, self.cols.select(pl.col('vtype') != 'continuous').to_series().to_numpy(), False + ) + cost = _scattered(count, self.obj['col'].to_numpy(), self.obj['coeff'].to_numpy(), 0.0) + np.nan_to_num(lb, copy=False, neginf=-infinity, posinf=infinity) + np.nan_to_num(ub, copy=False, neginf=-infinity, posinf=infinity) + return lb, ub, cost, integral + + +def _scattered(count: int, at: Any, values: Any, absent: Any) -> Any: + """*values* written at the column each one belongs to, *absent* elsewhere.""" + import numpy as np + + dense = np.full(count, absent, dtype=values.dtype) + dense[at] = values + return dense diff --git a/src/lpspec/relational/sinks/writers/__init__.py b/src/lpspec/relational/sinks/writers/__init__.py new file mode 100644 index 000000000..eec82ad50 --- /dev/null +++ b/src/lpspec/relational/sinks/writers/__init__.py @@ -0,0 +1,44 @@ +"""The writer family: the tables in, a file out. See ../README.md. + +One module per format, chosen by the output's **suffix** — the caller names an +output, not a writer, because a file's format is a property of the file. Each +answers ``(tables, path) -> None``, and streams. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +from lpspec.relational.sinks.writers.lp_file import write_lp_file + +if TYPE_CHECKING: + from collections.abc import Callable, Mapping + from pathlib import Path + + from lpspec.relational.sinks.tables import ModelTables + + Write = Callable[[ModelTables, Path], None] + +__all__ = ['PLANNED_WRITERS', 'WRITERS', 'write_lp_file', 'writer'] + +#: What can be written today, by suffix. Closed, for +#: :data:`~lpspec.relational.sinks.solvers.SOLVERS`' reason. +WRITERS: Mapping[str, Write] = {'.lp': write_lp_file} + +#: Formats with a module coming. Separate from :data:`WRITERS` because "not +#: yet" and "no" are different answers and a caller acts differently on each. +PLANNED_WRITERS: Mapping[str, str] = { + '.mps': 'the mps writer is planned but not implemented yet (docs/ARCHITECTURE.md, sinks)', +} + + +def writer(suffix: str) -> Write: + """The writer for *suffix* — or a ``NotImplementedError`` naming a format + that is coming, or a ``ValueError`` listing what can be written.""" + if suffix in WRITERS: + return WRITERS[suffix] + if suffix in PLANNED_WRITERS: + raise NotImplementedError(PLANNED_WRITERS[suffix]) + supported = ', '.join(sorted(WRITERS)) + planned = ', '.join(sorted(PLANNED_WRITERS)) + raise ValueError(f'unsupported output format {suffix!r} — supported: {supported} (planned: {planned})') diff --git a/src/lpspec/relational/sinks/lp_file.py b/src/lpspec/relational/sinks/writers/lp_file.py similarity index 100% rename from src/lpspec/relational/sinks/lp_file.py rename to src/lpspec/relational/sinks/writers/lp_file.py diff --git a/tests/conftest.py b/tests/conftest.py index 0fc6d41c2..ada55a035 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -103,7 +103,7 @@ def raw_of(source: str | Path | dict[str, Any]) -> dict[str, Any]: def solve_lp_file(path: Path | str) -> float: """Objective HiGHS reaches reading the written LP file back from disk. - The third opinion in a differential: ``solver_direct`` builds the model + The third opinion in a differential: the ``highs`` solver builds the model through the HiGHS API, this one round-trips it through text, and a sink that writes a wrong file is otherwise invisible. Lives here rather than in ``tests.differential`` because highspy is a core dependency — a bare diff --git a/tests/differential.py b/tests/differential.py index 1ce8dfb15..82db2c50f 100644 --- a/tests/differential.py +++ b/tests/differential.py @@ -109,7 +109,7 @@ def differential( lp_path = None if lp: lp_path = work / 'model.lp' - ex.write_lp(lp_path) + ex.write(lp_path) assert solve_lp_file(lp_path) == pytest.approx(oracle, rel=rel) yield Agreement(oracle=oracle, model=m, result=result, schema=schema, executor=ex, lp=lp_path) diff --git a/tests/test_api.py b/tests/test_api.py index 21477597b..c76ca5c7d 100644 --- a/tests/test_api.py +++ b/tests/test_api.py @@ -35,7 +35,7 @@ def test_solve(dispatch_yaml, dispatch_frame_inputs): result.close() -def test_build_context_manager_and_write_lp(dispatch_yaml, dispatch_frame_inputs, tmp_path): +def test_build_context_manager_and_write(dispatch_yaml, dispatch_frame_inputs, tmp_path): sources, coords = dispatch_frame_inputs with lps.build(dispatch_yaml, sources, coords=coords) as ex: result = ex.solve() @@ -149,8 +149,6 @@ def test_check_reports_language_errors_before_any_data_is_bound( def test_error_hierarchy_is_one_catchable_tree(): """One ``except`` covers the package, and the model/run split is real.""" - from lpspec.relational import RelationalBuildError - for cls in (lps.LanguageError, lps.DataError): assert issubclass(cls, lps.LpspecError) for cls in (lps.SchemaError, lps.DimensionError, lps.PiecewiseExpansionError): @@ -158,8 +156,19 @@ def test_error_hierarchy_is_one_catchable_tree(): assert not issubclass(lps.DataError, lps.LanguageError) assert issubclass(lps.LpspecError, ValueError) - # the retired name still catches everything it used to - assert RelationalBuildError is lps.LpspecError + +def test_an_unknown_solver_is_refused_with_the_alternatives(dispatch_yaml, dispatch_frame_inputs): + """The set of solvers is closed, and a name outside it never falls back to + the default — solving with a solver other than the one asked for is the one + answer that cannot be right. Here rather than in ``test_gurobi_sink.py``, + which skips without the extra: the closed set is a property of the package, + not of gurobi. Refused before the build, as an unwritable suffix is.""" + from lpspec.relational.sinks import SOLVERS + + sources, coords = dispatch_frame_inputs + with pytest.raises(lps.LpspecError, match='unknown solver'): + lps.solve(dispatch_yaml, sources, solver_name='cplex', coords=coords) + assert set(SOLVERS) == {'highs', 'gurobi'} def test_multi_file_composition_reserved(dispatch_yaml): diff --git a/tests/test_architecture.py b/tests/test_architecture.py index 23a28d3b3..bf57bb903 100644 --- a/tests/test_architecture.py +++ b/tests/test_architecture.py @@ -379,6 +379,66 @@ def test_expansion_has_no_mutable_module_state(): ) +#: The two sink families. The directory *is* the family, so a member cannot +#: land in the wrong one by being spelled differently. +SINKS = PKG / 'relational' / 'sinks' + + +def _family(name: str) -> set[str]: + return {p.stem for p in (SINKS / name).glob('*.py') if p.stem != '__init__'} + + +def test_each_sink_family_is_its_directory_and_its_registry(): + """One shape per family, checked off the path. + + A solver is four things that must agree: a module under ``solvers/`` named + for the solver, a ``solve_``, a ``build_`` seam for `bench/`, + and the ``SOLVERS`` key. Writers are keyed by suffix instead, but the rule + is the same. Agreement is what makes adding one mechanical — nothing above + the module to teach, nothing to remember but the name. + """ + import importlib + + from lpspec.relational.sinks import PLANNED_WRITERS, SOLVERS, WRITERS + + solvers = _family('solvers') + assert set(SOLVERS) == solvers, f'solver modules and SOLVERS keys disagree: {solvers ^ set(SOLVERS)}' + for name in sorted(solvers): + module = importlib.import_module(f'lpspec.relational.sinks.solvers.{name}') + assert SOLVERS[name] is getattr(module, f'solve_{name}'), f'SOLVERS[{name!r}] is not {name}.solve_{name}' + assert hasattr(module, f'build_{name}'), f'{name} has no build_{name}: the load-only seam `bench/` measures' + + assert {w.__module__.rsplit('.', 1)[-1] for w in WRITERS.values()} == _family('writers') + assert all(s.startswith('.') for s in (*WRITERS, *PLANNED_WRITERS)), 'writers are keyed by file suffix' + assert not set(WRITERS) & set(PLANNED_WRITERS), 'a format is either written or planned, never both' + + +def test_no_sink_reaches_a_sibling(): + """The fence that keeps an optional dependency optional. + + ``gurobipy`` stays the ``gurobi`` module's alone only because no other + sink imports it, directly or by importing the module that does. Each leaf + reads ``tables.py`` and nothing else in the family. + """ + offenders = {} + for family in ('solvers', 'writers'): + for path in sorted((SINKS / family).glob('*.py')): + reached = { + node.module + for node in ast.walk(ast.parse(path.read_text())) + if isinstance(node, ast.ImportFrom) + and node.module + and node.module.startswith('lpspec.relational.sinks.') + and not node.module.endswith('.tables') + } + if reached and path.stem != '__init__': + offenders[f'{family}/{path.name}'] = sorted(reached) + assert not offenders, ( + f'sink modules reaching a sibling: {offenders} — a sink reads tables.py and its own ' + f'dependency; anything shared belongs on ModelTables, where both families can see it' + ) + + def test_every_plan_node_is_handled_by_the_compiler(): """Two-tier economy: a primitive is not done until the engine consumes it. @@ -431,6 +491,11 @@ def test_every_module_is_documented_somewhere(): that map exists *not* to be. A ``README.md`` beside the code counts instead: it is what you read when you open the directory, and it stays next to the thing it describes. + + "Beside" reaches *up* as well as across, because a family may be a + directory of its own — ``sinks/solvers/highs.py`` is documented by + ``sinks/README.md``, which is the page describing both families. One + README per tree, not one per level. """ architecture = (REPO / 'docs/ARCHITECTURE.md').read_text() missing = [] @@ -440,8 +505,8 @@ def test_every_module_is_documented_somewhere(): continue # private plumbing (_notes) needs no doc entry if name == '__init__.py': continue - local_readme = path.parent / 'README.md' - documented = name in architecture or (local_readme.exists() and name in local_readme.read_text()) + readmes = [d / 'README.md' for d in path.parents if PKG in d.parents or d == PKG] + documented = name in architecture or any(r.exists() and name in r.read_text() for r in readmes) if not documented: missing.append(str(path.relative_to(PKG))) assert not missing, ( diff --git a/tests/test_group_sum.py b/tests/test_group_sum.py index e9a2ba75d..dde1f01db 100644 --- a/tests/test_group_sum.py +++ b/tests/test_group_sum.py @@ -2,7 +2,7 @@ Three-way differential on examples/transport.yaml: 1. eager lpspec_linopy.build + solve (group_sum via linopy groupby) - 2. lowered Program -> PolarsExecutor solver_direct, plus the LP file + 2. lowered Program -> PolarsExecutor -> the `highs` solver, plus the LP file 3. hand-built indicator-matrix linopy model (an independent oracle that involves no group_sum at all) """ diff --git a/tests/test_gurobi_sink.py b/tests/test_gurobi_sink.py new file mode 100644 index 000000000..3e4f372a2 --- /dev/null +++ b/tests/test_gurobi_sink.py @@ -0,0 +1,237 @@ +"""The ``gurobi`` sink, against the sink that was already here. + +Two sinks loading one :class:`ModelTables` must produce the same model, so +HiGHS is the oracle for Gurobi the way linopy is the oracle for the math: the +interesting assertions here are agreements, not values. Where a value *is* +asserted it comes from ``examples/ports/references.json`` — somebody else's +published optimum, which neither sink can talk the other into. + +Every test skips without ``gurobipy``. It ships a size-limited licence in its +own wheel, which is what makes this runnable in CI at all, so the models here +stay small enough for it — a few hundred columns, where the limit is 2000. +""" + +from __future__ import annotations + +import builtins +import gc +import weakref +from typing import Any + +import polars as pl +import pytest + +import lpspec as lps +from lpspec.errors import LpspecError, NoSolutionError +from lpspec.relational.sinks.solvers.gurobi import build_gurobi +from tests.test_ports import PORTS, REFERENCES, sources + +gurobipy = pytest.importorskip('gurobipy', reason='the gurobi sink needs the [gurobi] extra') + + +LP = { + 'dimensions': {'t': {'dtype': 'int', 'values': [0, 1, 2]}}, + 'parameters': {'load': {'dims': ['t']}, 'price': {'dims': ['t']}}, + 'variables': {'p': {'foreach': ['t'], 'bounds': {'lower': 0, 'upper': 100}}}, + 'constraints': {'meet': {'foreach': ['t'], 'expression': 'p >= load'}}, + 'objectives': {'cost': {'sense': 'minimize', 'expression': 'sum(p * price, over=t)'}}, +} + +#: Maximisation *and* an objective constant, which are the two things the +#: sink states outside the frames: ``ModelSense`` and ``ObjCon``. +MAX = { + 'dimensions': {'t': {'dtype': 'int', 'values': [0, 1]}}, + 'parameters': {'cap': {'dims': ['t']}}, + 'variables': {'p': {'foreach': ['t'], 'bounds': {'lower': 0, 'upper': 10}}}, + 'constraints': {'lim': {'foreach': ['t'], 'expression': 'p <= cap'}}, + 'objectives': {'profit': {'sense': 'maximize', 'expression': 'sum(p, over=t) + 5'}}, +} + +MIP = { + 'dimensions': {'i': {'dtype': 'int', 'values': [0, 1, 2]}, 'one': {'dtype': 'int', 'values': [0]}}, + 'parameters': {'w': {'dims': ['i']}, 'cap': {'dims': ['one']}}, + 'variables': {'x': {'foreach': ['i'], 'binary': True}}, + 'constraints': {'budget': {'foreach': ['one'], 'expression': 'sum(x * w, over=i) <= cap'}}, + 'objectives': {'o': {'sense': 'maximize', 'expression': 'sum(x * w, over=i)'}}, +} + +INFEASIBLE = { + 'dimensions': {'t': {'dtype': 'int', 'values': [0]}}, + 'parameters': {'load': {'dims': ['t']}}, + 'variables': {'p': {'foreach': ['t'], 'bounds': {'lower': 0, 'upper': 1}}}, + 'constraints': {'meet': {'foreach': ['t'], 'expression': 'p == load'}}, + 'objectives': {'c': {'sense': 'minimize', 'expression': 'p'}}, +} + +DATA: dict[str, dict[str, Any]] = { + 'LP': { + 'load': pl.DataFrame({'t': [0, 1, 2], 'value': [1.0, 2.0, 3.0]}), + 'price': pl.DataFrame({'t': [0, 1, 2], 'value': [10.0, 20.0, 30.0]}), + }, + 'MAX': {'cap': pl.DataFrame({'t': [0, 1], 'value': [3.0, 4.0]})}, + 'MIP': { + 'w': pl.DataFrame({'i': [0, 1, 2], 'value': [2.0, 3.0, 4.0]}), + 'cap': pl.DataFrame({'one': [0], 'value': [5.0]}), + }, + 'INFEASIBLE': {'load': pl.DataFrame({'t': [0], 'value': [99.0]})}, +} + +MODELS = {'LP': LP, 'MAX': MAX, 'MIP': MIP} + + +# --------------------------------------------------------------------------- +# the two sinks answer the same question the same way +# --------------------------------------------------------------------------- + + +@pytest.mark.parametrize( + ('name', 'variable', 'constraint'), [('LP', 'p', 'meet'), ('MAX', 'p', 'lim'), ('MIP', 'x', None)] +) +def test_gurobi_and_highs_agree(name: str, variable: str, constraint: str | None) -> None: + """The claim the second solver has to earn, on all three quantities. + + Coordinates as well as values, since a sink that loaded the columns in a + different order would still reach the same objective on these models — and + duals under ``maximize``, where a sign convention could differ and nothing + else in the suite would notice. + """ + with lps.solve(MODELS[name], DATA[name]) as highs, lps.solve(MODELS[name], DATA[name], solver_name='gurobi') as gb: + assert gb.termination_condition == highs.termination_condition + assert gb.objective == pytest.approx(highs.objective) + + expected, got = highs.primal(variable), gb.primal(variable) + assert got.columns == expected.columns + assert got.drop('value').equals(expected.drop('value')) + assert got['value'].to_list() == pytest.approx(expected['value'].to_list()) + + if constraint: + assert gb.dual(constraint)['value'].to_list() == pytest.approx(highs.dual(constraint)['value'].to_list()) + + +@pytest.mark.parametrize('port', sorted(REFERENCES), ids=str) +def test_every_port_reaches_its_reference_optimum_on_gurobi(port: str) -> None: + """``test_ports.py``'s corpus, solved by the other solver. + + The one assertion here no part of this package produced. A sink that + mis-loads the matrix — a block boundary off by a row, a sense inverted — + still reaches *a* number; this is what that number is checked against. + """ + reference = REFERENCES[port] + with lps.solve(PORTS / f'{port}.yaml', sources(port), solver_name='gurobi') as solution: + assert solution.is_ok, f'{port} did not solve: {solution.status}' + assert solution.objective == pytest.approx(reference['objective'], rel=reference['rtol']) + + +def test_block_boundaries_do_not_move_the_answer() -> None: + """``batch_rows=1`` forces one block per row, so every CSR view is built at + a boundary — where an off-by-one in ``indptr`` shifts coefficients into the + neighbouring row rather than dropping them.""" + with lps.build(LP, DATA['LP']) as ex: + whole = ex.solve(solver_name='gurobi') + ragged = ex.solve(batch_rows=1, solver_name='gurobi') + assert ragged.objective == pytest.approx(whole.objective) + assert ragged.primal('p')['value'].to_list() == pytest.approx(whole.primal('p')['value'].to_list()) + + +# --------------------------------------------------------------------------- +# what the sink says when there is nothing to read +# --------------------------------------------------------------------------- + + +def test_an_infeasible_solve_reports_both_axes_in_gurobis_wording() -> None: + with lps.solve(INFEASIBLE, DATA['INFEASIBLE'], solver_name='gurobi') as solution: + assert solution.status == 'warning' + assert solution.termination_condition == 'infeasible' + assert not solution.has_primal + assert solution.objective != solution.objective # nan, not 0.0 + with pytest.raises(NoSolutionError, match='INFEASIBLE'): + solution.primal('p') + + +def test_a_mixed_integer_model_has_no_duals() -> None: + """Gurobi refuses ``Pi`` rather than returning zeros; the sink passes the + refusal on as the ``None`` that makes ``dual`` explain itself.""" + with lps.solve(MIP, DATA['MIP'], solver_name='gurobi') as solution: + assert solution.has_primal + with pytest.raises(LpspecError, match='mixed-integer'): + solution.dual('budget') + + +def test_solver_options_reach_gurobi() -> None: + """Verbatim, in Gurobi's own vocabulary — ``TimeLimit``, not HiGHS' + ``time_limit``. Forwarding is the contract; translating names is not, and + an option the solver does not know reaches the caller as the solver's own + complaint rather than as a guess at what was meant.""" + with lps.solve(MIP, DATA['MIP'], solver_options={'TimeLimit': 0.0}, solver_name='gurobi') as solution: + assert solution.termination_condition == 'time_limit' + with pytest.raises(gurobipy.GurobiError, match='no_such_parameter'): + lps.solve(MIP, DATA['MIP'], solver_options={'no_such_parameter': 1}, solver_name='gurobi') + + +# --------------------------------------------------------------------------- +# the seams: the build without the search, and choosing a sink at all +# --------------------------------------------------------------------------- + + +def test_solver_options_land_on_the_environment() -> None: + """Where a licence parameter has to go. + + ``WLSAccessID`` / ``ComputeServer`` / ``TokenServer`` can only be set + before an environment starts, so applying options to the *model* — as this + sink first did — locks out every Compute-Server and WLS user. Asserted + through an ordinary parameter, since a licence one would need a licence: + the model sees it as its default, which is what environment-level means. + """ + with lps.build(MIP, DATA['MIP']) as ex: + assert build_gurobi(ex._tables(), solver_options={'TimeLimit': 5.0}).Params.TimeLimit == 5.0 + + +def test_build_gurobi_loads_the_model_and_stops() -> None: + """`bench/`'s seam: the hand-off with no search behind it, so what it + reports is what was loaded rather than what was solved.""" + with lps.build(MIP, DATA['MIP']) as ex: + tables = ex._tables() + m = build_gurobi(tables) + assert (m.NumVars, m.NumConstrs) == (tables.column_count, tables.row_count) + assert m.NumIntVars == tables.cols.filter(pl.col('vtype') != 'continuous').height + assert m.ModelSense == gurobipy.GRB.MAXIMIZE + assert m.SolCount == 0 + + +def test_nothing_keeps_a_built_model_alive() -> None: + """The precondition for releasing the licence a built model holds. + + :func:`build_gurobi` hands ownership over and disposes the environment + through a finalizer on the model, so anything in this package still + referencing that model would hold a Gurobi licence open for the life of + the process. The one-shot :func:`solve_gurobi` path does not depend on + this — it disposes both in a ``finally``. + """ + with lps.build(MIP, DATA['MIP']) as ex: + reference = weakref.ref(build_gurobi(ex._tables())) + gc.collect() + assert reference() is None, 'a built gurobi model outlived its caller — its environment cannot be released' + + +def test_the_objective_constant_rides_on_the_model_not_the_answer() -> None: + """Gurobi has ``ObjCon``, so the constant is part of the model it holds — + which makes the build seam a complete hand-off rather than a model plus a + number to remember.""" + with lps.build(MAX, DATA['MAX']) as ex: + assert build_gurobi(ex._tables()).ObjCon == pytest.approx(5.0) + + +def test_the_missing_extra_is_named() -> None: + """What a caller without gurobipy meets — both halves named, since the + absent one is as often scipy.""" + real_import = builtins.__import__ + + def refuse(name: str, *args: Any, **kwargs: Any) -> Any: + if name in {'gurobipy', 'scipy.sparse'}: + raise ModuleNotFoundError(f'No module named {name!r}') + return real_import(name, *args, **kwargs) + + with lps.build(LP, DATA['LP']) as ex, pytest.MonkeyPatch.context() as patch: + patch.setattr(builtins, '__import__', refuse) + with pytest.raises(ModuleNotFoundError, match=r'\[gurobi\] extra \(gurobipy, scipy\)'): + ex.solve(solver_name='gurobi') diff --git a/tests/test_lp_text.py b/tests/test_lp_text.py index e10656dc2..8348ffd60 100644 --- a/tests/test_lp_text.py +++ b/tests/test_lp_text.py @@ -22,7 +22,7 @@ import pytest import lpspec as lps -from lpspec.relational.sinks.lp_file import _number, _signed +from lpspec.relational.sinks.writers.lp_file import _number, _signed from tests.conftest import DISPATCH_MODEL, override #: Doubles that break naive formatters: repeating binary fractions, the @@ -115,7 +115,7 @@ def test_written_bounds_are_bit_exact() -> None: with tempfile.TemporaryDirectory() as tmp: lp = Path(tmp) / 'model.lp' with lps.build(DISPATCH_MODEL, data) as ex: - ex.write_lp(lp) + ex.write(lp) text = lp.read_text() section = text.split('bounds\n')[1].split('\nend')[0] @@ -150,7 +150,7 @@ def test_one_model_writes_the_same_bytes_every_time(tmp_path: Path) -> None: with lps.build(schema, data) as ex: for attempt in range(3): lp = tmp_path / f'{attempt}.lp' - ex.write_lp(lp) + ex.write(lp) written.append(hashlib.sha256(lp.read_bytes()).hexdigest()) assert len(set(written)) == 1, 'the same model wrote different bytes' @@ -177,7 +177,7 @@ def test_section_keywords_survive_sections_far_larger_than_a_buffer(tmp_path: Pa lp = tmp_path / 'model.lp' with lps.build(schema, data) as ex: - ex.write_lp(lp) + ex.write(lp) lines = lp.read_text().splitlines() keywords = ['min', 'obj:', 's.t.', 'bounds', 'end'] diff --git a/tests/test_milp.py b/tests/test_milp.py index 13e1a7510..4ecce3885 100644 --- a/tests/test_milp.py +++ b/tests/test_milp.py @@ -2,7 +2,7 @@ Binary commitment variables u with p <= p_max * u and a fixed commitment cost. Verifies the relational backend's vtype path end to end: cols vtype -column, HiGHS changeColsIntegrality in solver_direct, and the LP binary +column, HiGHS changeColsIntegrality in the `highs` solver, and the LP binary section. """ @@ -90,7 +90,7 @@ def test_commitment_milp_agrees_and_stays_integral(commitment_inputs): @pytest.mark.parametrize('batch_rows', [7, 13, 100_000], ids=['tiny-chunks', 'odd-chunks', 'one-chunk']) -def test_solver_direct_ingests_columns_in_order_whatever_the_chunking(commitment_inputs, batch_rows): +def test_the_highs_solver_ingests_columns_in_order_whatever_the_chunking(commitment_inputs, batch_rows): """Columns reach HiGHS in label order however the range loop splits them. ``addCols`` appends, so column *k* must be the *k*-th row handed over. The diff --git a/tests/test_relational.py b/tests/test_relational.py index b10ec91da..483e9a7f6 100644 --- a/tests/test_relational.py +++ b/tests/test_relational.py @@ -1,7 +1,7 @@ """Phase-2 gate: two real models round-trip through solve on the relational backend. Each model is built three ways and must agree on the objective: - 1. relational executor -> solver_direct (HiGHS via batched addCols/addRows) + 1. relational executor -> the `highs` solver (batched addCols/addRows) 2. relational executor -> lp_file sink -> HiGHS reads and solves the file 3. eager linopy build (the correctness oracle) """ @@ -133,7 +133,7 @@ def test_dispatch_roundtrip(dispatch_data, tmp_path): assert result.objective == pytest.approx(oracle, rel=RTOL) lp = tmp_path / 'dispatch.lp' - ex.write_lp(lp) + ex.write(lp) assert solve_lp_file(lp) == pytest.approx(oracle, rel=RTOL) # masked variable rows are absent, and primal joins back to coords @@ -223,7 +223,7 @@ def test_transport_roundtrip(transport_data, tmp_path): assert result.objective == pytest.approx(oracle, rel=RTOL) lp = tmp_path / 'transport.lp' - ex.write_lp(lp) + ex.write(lp) assert solve_lp_file(lp) == pytest.approx(oracle, rel=RTOL) # flows respect line capacity bounds diff --git a/tests/test_solve_status.py b/tests/test_solve_status.py index 3b7fac124..a31381e4a 100644 --- a/tests/test_solve_status.py +++ b/tests/test_solve_status.py @@ -1,10 +1,14 @@ """The solve outcome, and linopy as its oracle. `relational/status.py` copies linopy's status vocabulary spelling for -spelling, and `sinks/highs.py` copies linopy's own HiGHS mapping. Copies rot. -These tests import linopy and compare, so a divergence — ours drifting, or a -linopy release moving — fails here instead of being discovered by a user who -knows one vocabulary and is handed another. +spelling, and each solver sink copies linopy's own mapping for its solver. +Copies rot. These tests import linopy and compare, so a divergence — ours +drifting, or a linopy release moving — fails here instead of being discovered +by a user who knows one vocabulary and is handed another. + +The gurobi map diverges from linopy's in three declared places, and that is +checked in both directions: a copy is only honest if the exceptions are as +pinned as the agreements. The engine itself never imports linopy (docs/ARCHITECTURE.md, hard rule 2). Tests may, and this is the same oracle arrangement the differential tests use for @@ -13,12 +17,16 @@ from __future__ import annotations +import ast +from typing import Any + import polars as pl import pytest import lpspec as lps from lpspec.errors import NoSolutionError -from lpspec.relational.sinks.highs import _CONDITION_OF_HIGHS_STATUS +from lpspec.relational.sinks.solvers.gurobi import _CONDITION_OF_GUROBI_STATUS, _LINOPY_DIVERGENCES +from lpspec.relational.sinks.solvers.highs import _CONDITION_OF_HIGHS_STATUS from lpspec.relational.status import STATUS_TO_TERMINATION_CONDITIONS, SolveStatus INFEASIBLE = { @@ -49,35 +57,69 @@ def test_the_status_rollup_matches_linopy(): def test_the_highs_mapping_matches_linopy(): - """linopy builds its map inside a method, so it is read from the source. + assert _linopy_condition_map('Highs', ast.Attribute, 'attr') == _CONDITION_OF_HIGHS_STATUS + + +def test_the_gurobi_mapping_matches_linopy_where_it_claims_to(): + """The same copy, and the same brittleness — plus three declared exceptions. - Brittle to a linopy refactor, deliberately: this table is a copy, and a - copy nobody checks is a copy that rots. If linopy moves it and this stops - finding it, the assertion below says so rather than passing vacuously. + linopy's Gurobi map contradicts Gurobi's own documented status codes in + three places, so copying it whole would import a wrong answer rather than + a shared vocabulary. Each is listed in ``_LINOPY_DIVERGENCES`` with its + reason, and this asserts **both** directions: everything else still + matches, and every declared divergence is still a divergence — so if + linopy fixes one, this fails and the exception goes away. + """ + theirs = _linopy_condition_map('Gurobi', ast.Constant, 'value') + assert set(theirs) == set(_CONDITION_OF_GUROBI_STATUS), ( + 'linopy and this package no longer cover the same Gurobi status codes' + ) + for code, condition in theirs.items(): + if code in _LINOPY_DIVERGENCES: + assert _CONDITION_OF_GUROBI_STATUS[code] != condition, ( + f'linopy now agrees with us on status {code} — drop the entry from _LINOPY_DIVERGENCES' + ) + else: + assert _CONDITION_OF_GUROBI_STATUS[code] == condition + + +def test_every_gurobi_divergence_stays_inside_linopys_vocabulary(): + """Diverging on a verdict is not licence to invent a word for it. Every + condition this package reports is one linopy also defines, which is what + keeps `status`, `is_ok` and the rollup meaningful across both.""" + assert set(_CONDITION_OF_GUROBI_STATUS.values()) <= set().union(*STATUS_TO_TERMINATION_CONDITIONS.values()) + + +def _linopy_condition_map(solver: str, node: type[ast.expr], attribute: str) -> dict[Any, Any]: + """linopy's ``CONDITION_MAP`` for *solver*, read out of its source. + + Each solver spells the map differently — HiGHS keys it by + ``HighsModelStatus`` attributes, Gurobi by integer literals — so the node + type and the attribute holding the value are arguments. + + Brittle to a linopy refactor, deliberately: the map is a local inside a + method, so there is nothing to import, and a copy nobody checks is a copy + that rots. If linopy moves it, the assertions say so rather than passing + vacuously. """ - import ast import inspect solvers = pytest.importorskip('linopy.solvers') tree = ast.parse(inspect.getsource(solvers)) - # several solver classes define a CONDITION_MAP; only the HiGHS one is ours - highs = next((n for n in ast.walk(tree) if isinstance(n, ast.ClassDef) and n.name == 'Highs'), None) - assert highs is not None, 'linopy no longer has a Highs solver class — re-verify the copy by hand' + cls = next((n for n in ast.walk(tree) if isinstance(n, ast.ClassDef) and n.name == solver), None) + assert cls is not None, f'linopy no longer has a {solver} solver class — re-verify the copy by hand' literals = [ - node.value - for node in ast.walk(highs) - if isinstance(node, (ast.AnnAssign, ast.Assign)) - and 'CONDITION_MAP' in ast.dump(node) - and isinstance(node.value, ast.Dict) + n.value + for n in ast.walk(cls) + if isinstance(n, (ast.AnnAssign, ast.Assign)) and 'CONDITION_MAP' in ast.dump(n) + if isinstance(n.value, ast.Dict) ] - assert literals, 'linopy no longer defines CONDITION_MAP as a dict literal — re-verify the copy by hand' - - theirs = { - key.attr: value.attr + assert literals, f'linopy no longer defines {solver}.CONDITION_MAP as a dict literal — re-verify by hand' + return { + getattr(key, attribute): getattr(value, attribute) for key, value in zip(literals[0].keys, literals[0].values, strict=True) - if isinstance(key, ast.Attribute) and isinstance(value, ast.Attribute) + if isinstance(key, node) and isinstance(value, node) } - assert theirs == _CONDITION_OF_HIGHS_STATUS # ---------------------------------------------------------------------------