From 1bd07f2827c73fdb61ae498a75219f3d20f04f73 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Fri, 31 Jul 2026 19:16:51 +0200 Subject: [PATCH 1/2] perf(sinks): a row is seated the way a column is, for both solvers MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `dense_columns` exists because `col` is dense `0..n-1` and so *is* the position a value belongs at. `row` is dense in exactly the same way, and neither solver was treating it that way: both sorted the row frame and then filtered it once per chunk, so `fleet/l` read the same 6M rows nine times over to hand each of them over once. So the row half joins the column half in the contract, where tables.py already argues it belongs: two sinks computing the same projection separately could disagree about the model they loaded. ModelTables.dense_rows(infinity) -> (sense, rhs) over the solver's row index ModelTables.row_blocks(budget) -> each chunk of rows with the entries it owns `dense_rows` stops at the sense rather than at bounds because that is where the two solvers part — HiGHS wants a row's lower/upper, Gurobi its comparison and right-hand side, and both are that pair spelled differently. `SENSE_CODES` is the spelling neither owns; it was already inside gurobi.py as a per-chunk `replace_strict`. `row_blocks` orders the matrix once and locates each chunk by binary search on the label column. Searched in polars rather than through numpy: pulling the column out to search it there is marginally faster and holds a second copy of it for the whole loop, 0.11 GB at `transport/l`, which is the wrong trade in a pass that exists to stay bounded. Paired in-process, alternating, min of five, same frames: highs fleet/l 6.0M rows / 14.0M nz 0.722s -> 0.644s -10.8% transport/l 1.4M rows / 12.6M nz 0.354s -> 0.308s -13.0% nodal/l 1.0M rows / 3.0M nz 0.140s -> 0.131s -6.7% sector/l 1.0M rows / 1.0M nz 0.058s -> 0.054s -6.9% dispatch/l 0.1M rows / 10.0M nz 0.403s -> 0.379s -6.0% fleet/m 0.063s -> 0.059s -6.9% **Gurobi does not move**, and that is worth stating rather than hiding: its own `addMVar`/`addMConstr` take 1.7 s where HiGHS takes 0.06 s for the same model, so the polars work this removes is invisible inside it. The sharing is justified by the duplication, not by a second speedup. Peak is unchanged on both — the dense row vectors replace a sorted copy of the row frame, and what is left is inside the run-to-run spread. Co-Authored-By: Claude Opus 5 (1M context) --- src/lpspec/relational/sinks/solvers/gurobi.py | 19 ++--- src/lpspec/relational/sinks/solvers/highs.py | 36 ++++----- src/lpspec/relational/sinks/tables.py | 78 +++++++++++++++++-- tests/test_relational.py | 31 ++++++++ 4 files changed, 125 insertions(+), 39 deletions(-) diff --git a/src/lpspec/relational/sinks/solvers/gurobi.py b/src/lpspec/relational/sinks/solvers/gurobi.py index d191d75a6..6ba766f5d 100644 --- a/src/lpspec/relational/sinks/solvers/gurobi.py +++ b/src/lpspec/relational/sinks/solvers/gurobi.py @@ -161,23 +161,16 @@ def _load( 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()) + for lo, hi, a, starts in model.row_blocks(batch): block = scipy.sparse.csr_matrix( (a['coeff'].to_numpy(), a['col'].to_numpy(), np.append(starts, a.height)), - shape=(rows.height, model.column_count), + shape=(hi - lo, model.column_count), ) - blocks.append(m.addMConstr(block, x, senses[rows['op'].to_numpy()], rows['rhs'].to_numpy())) + blocks.append(m.addMConstr(block, x, spelling[sense[lo:hi]], rhs[lo:hi])) if model.objective_sense == 'max': m.ModelSense = gurobipy.GRB.MAXIMIZE diff --git a/src/lpspec/relational/sinks/solvers/highs.py b/src/lpspec/relational/sinks/solvers/highs.py index cd8e86e6e..9b7ec1835 100644 --- a/src/lpspec/relational/sinks/solvers/highs.py +++ b/src/lpspec/relational/sinks/solvers/highs.py @@ -6,8 +6,9 @@ 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 ``'<='`` is made -in polars and only the bound it selects is handed over — the same rule +boxing every value as a Python object, so a row's ``'<='`` becomes a +:data:`~lpspec.relational.sinks.tables.SENSE_CODES` byte before it is read here +— the same rule :meth:`~lpspec.relational.sinks.tables.ModelTables.dense_columns` states for the column vectors, where it was measured. @@ -22,6 +23,7 @@ import polars as pl from lpspec.errors import LpspecError +from lpspec.relational.sinks.tables import SENSE_CODES from lpspec.relational.status import SolveStatus if TYPE_CHECKING: @@ -39,9 +41,9 @@ #: #: Deliberately small, and the opposite of what the same budget settled at on #: the duckdb engine. There a wider chunk was worth a third of the hand-off, -#: because every chunk re-ran an ordered query; here the columns are numpy -#: slices and the rows a filter over an already-sorted frame, so more chunks -#: cost almost nothing and only residency scales with the budget. Measured at +#: because every chunk re-ran an ordered query; here both columns and rows are +#: numpy slices, so more chunks cost almost nothing and only residency scales +#: with the budget. Measured at #: `l`: 2e6 against 1e5 is 0.50s against 0.59s on `dispatch` and 0.71 against #: 0.75 on `transport`, for 0.6 GB and 0.2 GB more peak. A tenth of a second #: on a hand-off that precedes a minute of simplex is not worth half a @@ -112,26 +114,18 @@ def build_highs( integrality = np.full(len(noncontinuous), int(highspy.HighsVarType.kInteger), dtype=np.uint8) h.changeColsIntegrality(len(noncontinuous), noncontinuous, integrality) - ordered_rows = model.rows.sort('row') - ordered_matrix = model.matrix.sort('row') - for lo, hi in model.row_chunks_by_nonzeros(batch): - rows = ordered_rows.filter(pl.col('row').is_between(lo, hi, closed='left')).select( - 'row', - pl.when(pl.col('sense') == '<=').then(-inf).otherwise(pl.col('rhs')).alias('rlb'), - pl.when(pl.col('sense') == '>=').then(inf).otherwise(pl.col('rhs')).alias('rub'), - ) - a = ordered_matrix.filter(pl.col('row').is_between(lo, hi, closed='left')) - rlb = rows['rlb'].to_numpy() - rub = rows['rub'].to_numpy() - starts = np.searchsorted(a['row'].to_numpy(), rows['row'].to_numpy()).astype(np.int32) + sense, rhs = model.dense_rows(inf) + rlb = np.where(sense == SENSE_CODES['<='], -inf, rhs) + rub = np.where(sense == SENSE_CODES['>='], inf, rhs) + for lo, hi, a, starts in model.row_blocks(batch): _loaded( h, h.addRows( - len(rlb), - rlb, - rub, + hi - lo, + rlb[lo:hi], + rub[lo:hi], a.height, - starts, + starts.astype(np.int32), a['col'].to_numpy().astype(np.int32), a['coeff'].to_numpy(), ), diff --git a/src/lpspec/relational/sinks/tables.py b/src/lpspec/relational/sinks/tables.py index 8bf0d80dd..a411db21d 100644 --- a/src/lpspec/relational/sinks/tables.py +++ b/src/lpspec/relational/sinks/tables.py @@ -3,10 +3,10 @@ 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. +Also the *projections* of them more than one sink needs — the dense column and +row vectors, and the matrix a block at a time. They belong to the contract +rather than to either solver, because two sinks computing them separately could +disagree about the model they loaded — the one thing neither may do. """ from __future__ import annotations @@ -30,6 +30,19 @@ npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.float64], npt.NDArray[np.bool_] ] + #: A sense code and a right-hand side, each as long as the model has rows. + DenseRows = tuple[npt.NDArray[np.uint8], npt.NDArray[np.float64]] + + #: One chunk of rows: the half-open label range, the matrix entries those + #: rows own, and the offset of each row's entries within them. + RowBlock = tuple[int, int, pl.DataFrame, npt.NDArray[np.int64]] + + +#: ``sense`` as a number, so a row's comparison crosses into numpy as one byte +#: rather than as a boxed Python string. The order is arbitrary and shared: a +#: solver indexes its own spelling with these, so the two cannot drift. +SENSE_CODES = {'<=': 0, '>=': 1, '==': 2} + @dataclass(frozen=True) class ModelTables: @@ -107,9 +120,64 @@ def dense_columns(self, infinity: float) -> DenseColumns: np.nan_to_num(ub, copy=False, neginf=-infinity, posinf=infinity) return lb, ub, cost, integral + def dense_rows(self, infinity: float) -> DenseRows: + """``(sense, rhs)`` as numpy vectors over the solver's row index. + + The row half of :meth:`dense_columns`, and for the same reason: ``row`` + is dense ``0..n-1``, so it *is* the position a value belongs at, and a + chunk of rows is then a slice rather than a search. Sorting the row + frame and filtering it once per chunk read the same 6M rows nine times + over on `fleet/l`, to hand each of them over once. + + It stops at the sense rather than at bounds because that is where the + two solvers part: HiGHS wants a row's ``lower``/``upper`` and Gurobi its + comparison and right-hand side, and both are this pair spelled + differently. :data:`SENSE_CODES` is the spelling neither owns. + + A row with no entry is left with a comparison nothing can fail — + ``>=`` against ``-infinity`` — rather than the ``== 0`` that would be a + real equality the model never stated. + """ + sided = self.rows.select( + 'row', + pl.col('sense').replace_strict(SENSE_CODES, return_dtype=pl.UInt8).alias('op'), + 'rhs', + ) + at = sided['row'].to_numpy() + sense = _scattered(self.row_count, at, sided['op'].to_numpy(), SENSE_CODES['>=']) + rhs = _scattered(self.row_count, at, sided['rhs'].to_numpy(), -infinity) + return sense, rhs + + def row_blocks(self, budget: int) -> Iterator[RowBlock]: + """Each chunk of rows with the matrix entries it owns. + + The matrix is ordered once, and a chunk is then a ``slice`` of it + located by binary search on the label column — the range is contiguous + because ``row`` is dense and the frame is sorted, so scanning for it + would re-read the whole model once per chunk. + + **Searched in polars rather than through numpy.** Pulling the label + column out to search it there is marginally faster and holds a second + copy of one column of the model for the whole loop, which is 0.11 GB at + `transport/l` — the wrong trade in a pass that exists to stay bounded. + + ``starts`` is each row's offset within the block, which is what both + solvers' matrix APIs ask for. A row with no entries takes the next + row's offset, and so occupies no span. + """ + import numpy as np + + ordered = self.matrix.sort('row') + label = ordered['row'] + for lo, hi in self.row_chunks_by_nonzeros(budget): + first = int(label.search_sorted(lo, 'left')) + last = int(label.search_sorted(hi, 'left')) + entries = ordered.slice(first, last - first) + yield lo, hi, entries, np.searchsorted(entries['row'].to_numpy(), np.arange(lo, hi)) + def _scattered(count: int, at: Any, values: Any, absent: Any) -> Any: - """*values* written at the column each one belongs to, *absent* elsewhere.""" + """*values* written at the label each one belongs to, *absent* elsewhere.""" import numpy as np dense = np.full(count, absent, dtype=values.dtype) diff --git a/tests/test_relational.py b/tests/test_relational.py index 483e9a7f6..227466e56 100644 --- a/tests/test_relational.py +++ b/tests/test_relational.py @@ -37,6 +37,7 @@ Variable, VariableDeclaration, ) +from lpspec.relational.sinks import SOLVERS from tests.conftest import by_coord, solve_lp_file from tests.differential import RTOL, differential from tests.oracle import linopy, pd, transport_eager_objective, xr @@ -715,6 +716,36 @@ def test_infinite_bounds_survive_the_handoff(dispatch_data): assert ex.solve().is_ok +@pytest.mark.parametrize('solver_name', sorted(SOLVERS)) +@pytest.mark.parametrize('batch_rows', [1, 2, 7, 100_000], ids=['one', 'two', 'odd', 'whole']) +def test_a_row_with_no_terms_keeps_its_seat_at_any_chunking(solver_name, batch_rows): + """Rows reach a solver by position, so an empty one still occupies one. + + `where: "t > 0"` leaves `balance` at `t = 0` with nothing to sum, and this + lane keeps the row: `0 == 5` is infeasible and says so. Both solvers now + take the row bounds from a dense vector and slice it per chunk, which is + the same hand-off only if every label in the range has a seat — a row that + fell out of the frame would take a comparison nothing can fail, leave the + constraint unenforced, and the model would come back solved against a + model that cannot be. + + Ragged batches because the range loop is where a seat would be lost, and a + round number is the one split that hides an off-by-one. Both solvers, + because the seating is now theirs jointly rather than either one's. + """ + model = { + 'dimensions': {'t': {'dtype': 'int', 'values': [0, 1, 2]}, 'g': {'values': ['a', 'b']}}, + 'parameters': {'load': {'dims': ['t']}}, + 'variables': {'p': {'foreach': ['t', 'g'], 'where': 't > 0', 'bounds': {'lower': 0, 'upper': 100}}}, + 'constraints': {'balance': {'foreach': ['t'], 'expression': 'sum(p, over=g) == load'}}, + 'objectives': {'o': {'sense': 'minimize', 'expression': 'sum(sum(p, over=g), over=t)'}}, + } + with lps.build(model, {'load': pl.DataFrame({'t': [0, 1, 2], 'value': [5.0, 4.0, 6.0]})}) as ex: + assert sorted(set(ex._tables().matrix['row'].to_list())) == [1, 2], 'row 0 is the orphan under test' + solution = ex.solve(batch_rows=batch_rows, solver_name=solver_name) + assert solution.termination_condition == 'infeasible' + + def test_row_chunks_are_bounded_by_nonzeros_not_by_rows(): """A chunk of rows is a chunk of *entries*, and only entries are residency. From 0a3068405a2d1b562cf4e149291aceb8a74b2626 Mon Sep 17 00:00:00 2001 From: FBumann <117816358+FBumann@users.noreply.github.com> Date: Fri, 31 Jul 2026 19:32:06 +0200 Subject: [PATCH 2/2] fix(gurobi): the operator table is derived from SENSE_CODES, not retyped The array `addMConstr` indexes was written out in `SENSE_CODES`' order, so the two had to be kept in step by hand. Nothing checks that: a wrong order is a model whose comparisons are silently permuted, and every solver would answer it confidently. Built from the mapping instead, keyed by name. A sense added to `SENSE_CODES` and not to `_GUROBI_SENSE` raises where it is built rather than mislabelling a row. Raised by CodeRabbit on #421. Co-Authored-By: Claude Opus 5 (1M context) --- src/lpspec/relational/sinks/solvers/gurobi.py | 33 ++++++++++++++++--- 1 file changed, 28 insertions(+), 5 deletions(-) diff --git a/src/lpspec/relational/sinks/solvers/gurobi.py b/src/lpspec/relational/sinks/solvers/gurobi.py index 6ba766f5d..402ec79d8 100644 --- a/src/lpspec/relational/sinks/solvers/gurobi.py +++ b/src/lpspec/relational/sinks/solvers/gurobi.py @@ -1,8 +1,8 @@ """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 same hand-off as :mod:`~lpspec.relational.sinks.solvers.highs`, reading the +same ``dense_columns``, ``dense_rows`` and ``row_blocks``, 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 @@ -24,6 +24,7 @@ import polars as pl from lpspec.errors import LpspecError +from lpspec.relational.sinks.tables import SENSE_CODES from lpspec.relational.status import SolveStatus if TYPE_CHECKING: @@ -162,8 +163,7 @@ def _load( x = m.addMVar(model.column_count, lb=lb, ub=ub, obj=cost, vtype=np.where(integral, 'I', 'C')) sense, rhs = model.dense_rows(gurobipy.GRB.INFINITY) - # indexed by SENSE_CODES, which is the order that pairing depends on - spelling = np.array([gurobipy.GRB.LESS_EQUAL, gurobipy.GRB.GREATER_EQUAL, gurobipy.GRB.EQUAL], dtype='=': 'GREATER_EQUAL', '==': 'EQUAL'} + + +def _spelled(gurobipy: Any) -> Any: + """:data:`SENSE_CODES` as the characters ``addMConstr`` wants, by code. + + Built from the mapping rather than written out in its order, so the two + cannot come to disagree: a wrong order here is a model whose comparisons + are silently permuted, which every solver would answer confidently. A + sense added to :data:`SENSE_CODES` and not to :data:`_GUROBI_SENSE` raises + instead. + """ + import numpy as np + + spelling = np.empty(len(SENSE_CODES), dtype=' Any: """The optional dependency, or a message naming the extra. Both halves are named, because the missing one is as often scipy."""