diff --git a/doc/release_notes.rst b/doc/release_notes.rst index ea9416ad..4cc8b881 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -59,6 +59,7 @@ Upcoming Version **Performance** +* Direct solver APIs no longer set variable and constraint names by default: ``Model.set_names_in_solver_io`` now defaults to ``False``, and ``to_highspy``, ``to_gurobipy``, ``to_mosek`` and ``to_xpress`` follow it when ``set_names`` is not given. The default names ``x{label}``/``c{label}`` carry no information, and the solution is mapped back by position. ``Model.solve`` still sets them when ``warmstart_fn`` or ``basis_fn`` is given, as basis files refer to names. Setting them made ``to_highspy`` about 2x slower and cost ~100 ms in ``to_gurobipy`` on a model with 200k variables. **Behaviour change:** pass ``set_names=True`` or ``Model(set_names_in_solver_io=True)`` to keep names on the native solver model, e.g. to write it to a file. (`#978 `__) * ``@``/``dot`` against a constant matrix that holds zeros no longer densifies the result to one term per contracted member. The zero-coefficient terms are dropped, so the term dimension shrinks to the widest non-zero cell. On PyPSA's Kirchhoff Voltage Law constraint (a cycle matrix with ~3 branches per cycle) this cuts the expression from 852 to 3 terms — 284x fewer cells — which in turn shrinks the downstream ``merge``. A constant without zeros is unaffected. (`#748 `__) * ``densify_terms`` (used by ``sum(drop_zeros=True)`` and the sparse ``@`` path) is now fully vectorised. It previously counted the non-zero positions with a Python loop that scaled quadratically in the number of non-zero terms — 127 s for a (2000 x 60) expression, now 3 ms — and allocated the compacted output at the full original term width. It now allocates only the compacted width and returns the expression unchanged when it holds no zeros. * Under v1, ``@``/``dot`` against a constant now runs as sparse linear algebra instead of building the dense broadcast intermediate (``self * other`` then ``.sum()``): peak memory scales with ``nnz(C) x nterm`` rather than with the full broadcast shape. A ``CSRLinearExpression``-backed operand (from ``groupby(...).sum(sparse=True)``) stays CSR-backed through ``@``, and the result is the compact canonical form (duplicate variables summed, terms label-ordered, explicit zeros pruned — cell activeness is carried by ``const`` alone). (`#748 `__, `#756 `__, `#925 `__) @@ -74,6 +75,7 @@ Upcoming Version **Bug fixes** +* ``Variable.get_solver_attribute`` for Gurobi now maps values by position for the direct API. It parsed labels from the solver's variable names, which returned the values of the wrong variables when names were off and a variable was masked. (`#978 `__) * ``densify_terms`` no longer raises on expressions without coordinate dimensions (``expr.sum(drop_zeros=True)`` over all dimensions) and now works on ``QuadraticExpression``, where it previously indexed the ``_factor`` axis as the term axis. * ``sum()`` over a dimension no longer raises when another dimension of the expression has size 0; it returns an expression without terms over the kept coordinates, as summing over the empty dimension itself already did. (https://github.com/PyPSA/linopy/issues/906) * A multi-key ``groupby`` now returns its groups sorted by key tuple, like the single-key path. The key combinations were numbered by iterating a ``set``, so the group order was arbitrary and changed between processes with ``PYTHONHASHSEED``. diff --git a/linopy/io.py b/linopy/io.py index d7cd0063..3f2bb8f4 100644 --- a/linopy/io.py +++ b/linopy/io.py @@ -798,7 +798,7 @@ def to_file( # Use very fast highspy implementation # Might be replaced by custom writer, however needs C/Rust bindings for performance h = solvers.Highs._build_solver_model( - m, explicit_coordinate_names=explicit_coordinate_names + m, explicit_coordinate_names=explicit_coordinate_names, set_names=True ) h.writeModel(str(fn)) else: @@ -813,26 +813,26 @@ def to_mosek( m: Model, task: Any | None = None, explicit_coordinate_names: bool = False, - set_names: bool = True, + set_names: bool | None = None, ) -> Any: """Build the MOSEK task for `m`.""" import mosek - if task is None: - task = mosek.Task() - return solvers.Mosek._build_solver_model( + solver = solvers.Mosek.from_model( m, - task, + io_api="direct", explicit_coordinate_names=explicit_coordinate_names, set_names=set_names, + task=mosek.Task() if task is None else task, ) + return solver._detach_solver_model() def to_gurobipy( m: Model, env: Any | None = None, explicit_coordinate_names: bool = False, - set_names: bool = True, + set_names: bool | None = None, ) -> Any: """Build the gurobipy.Model for `m`.""" solver = solvers.Gurobi.from_model( @@ -848,7 +848,7 @@ def to_gurobipy( def to_highspy( m: Model, explicit_coordinate_names: bool = False, - set_names: bool = True, + set_names: bool | None = None, ) -> Highs: """Build the highspy.Highs instance for `m`.""" solver = solvers.Highs.from_model( @@ -863,14 +863,16 @@ def to_highspy( def to_xpress( m: Model, explicit_coordinate_names: bool = False, - set_names: bool = True, + set_names: bool | None = None, ) -> Any: """Build the xpress.problem instance for `m`.""" - return solvers.Xpress._build_solver_model( + solver = solvers.Xpress.from_model( m, + io_api="direct", explicit_coordinate_names=explicit_coordinate_names, set_names=set_names, ) + return solver._detach_solver_model() def to_cupdlpx(m: Model) -> cupdlpxModel: diff --git a/linopy/model.py b/linopy/model.py index 1e4f005b..80aef147 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -237,7 +237,7 @@ def __init__( force_dim_names: bool = False, auto_mask: bool = False, freeze_constraints: bool | None = None, - set_names_in_solver_io: bool = True, + set_names_in_solver_io: bool = False, dtypes: Mapping[DtypeKey, type[np.signedinteger]] | None = None, sparse: bool = False, ) -> None: @@ -270,7 +270,8 @@ def __init__( ``sparse=True``. The default is False. set_names_in_solver_io : bool Whether direct solver exports should include variable and - constraint names by default. The default is True. + constraint names by default. Names cost build time and are not + needed to map the solution back. The default is False. dtypes : mapping, optional Integer dtypes for the model's data, exposed read-only as ``Model.dtypes``. Only ``"labels"`` is supported, e.g. @@ -2115,8 +2116,10 @@ def solve( set_names : bool, optional Whether to set variable and constraint names when using the direct solver API (io_api='direct'). Setting to False can significantly - speed up model export. If None, uses the model default - ``Model.set_names_in_solver_io`` setting (default True). + speed up model export. If None, names are set when + ``warmstart_fn`` or ``basis_fn`` is given, as basis files refer + to names, and otherwise the model default + ``Model.set_names_in_solver_io`` (default False) applies. problem_fn : path_like, optional Path of the lp file or output file/directory which is written out during the process. The default None results in a temporary file. @@ -2291,8 +2294,8 @@ def solve( try: self.solver = None # closes any previous solver if io_api == "direct": - if set_names is None: - set_names = self.set_names_in_solver_io + if set_names is None and (warmstart_fn or basis_fn): + set_names = True build_kwargs: dict[str, Any] = { "explicit_coordinate_names": explicit_coordinate_names, "set_names": set_names, @@ -2664,22 +2667,16 @@ def compute_infeasibilities(self) -> list[int]: def _compute_infeasibilities_gurobi(self, solver_model: Any) -> list[int]: """Compute infeasibilities for Gurobi solver.""" + solver = self.solver + assert solver is not None solver_model.computeIIS() - f = NamedTemporaryFile(suffix=".ilp", prefix="linopy-iis-", delete=False) - solver_model.write(f.name) - labels = [] - pattern = re.compile(r"^ [^:]+#([0-9]+):") - for line in f.readlines(): - line_decoded = line.decode() - try: - if line_decoded.startswith(" c"): - labels.append(int(line_decoded.split(":")[0][2:])) - except ValueError as _: - match = pattern.match(line_decoded) - if match: - labels.append(int(match.group(1))) - f.close() - return labels + constrs = solver_model.getConstrs() + in_iis = np.asarray(solver_model.getAttr("IISConstr", constrs), dtype=bool) + if solver.io_api == "direct": + clabels = self.constraints.label_index.clabels + else: + clabels = solvers._names_to_labels([c.ConstrName for c in constrs]) + return sorted(int(label) for label in clabels[in_iis] if label >= 0) def _compute_infeasibilities_xpress(self, solver_model: Any) -> list[int]: """ diff --git a/linopy/solvers.py b/linopy/solvers.py index c1641649..6e1ec74d 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -714,6 +714,8 @@ def _build(self, **build_kwargs: Any) -> None: raise RuntimeError("Solver has no model attached; cannot build.") self._validate_model() if self.io_api == "direct": + if build_kwargs.get("set_names") is None: + build_kwargs["set_names"] = self.model.set_names_in_solver_io self._build_direct(**build_kwargs) if self.track_updates: self.snapshot = ModelSnapshot.capture(self.model) @@ -756,7 +758,7 @@ def _validate_model(self) -> None: "Use a solver that supports them." ) - def _build_direct(self, **build_kwargs: Any) -> None: + def _build_direct(self, *, set_names: bool, **build_kwargs: Any) -> None: """Build the native solver model from ``self.model``. Override per-solver.""" raise NotImplementedError( f"Solver {self.solver_name.value} does not support direct API model export." @@ -1116,6 +1118,17 @@ def _cache_model_sizes(self, model: Model) -> None: def update_solver_model(self, model: Model, **kwargs: Any) -> None: raise NotImplementedError + def _detach_solver_model(self) -> Any: + """ + Hand ownership of the native solver model to the caller: unregister it + from the solver's teardown so ``close()`` does not dispose it. + """ + m = self.solver_model + if self._env_stack is not None: + self._env_stack.pop_all() + self.close() + return m + def close(self) -> None: """ Dispose the native solver model and env, releasing any held license. @@ -1636,7 +1649,8 @@ def _apply_obj_sense(self, ctx: Any, sense: str) -> None: def _build_direct( self, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, log_fn: Path | None = None, **kwargs: Any, ) -> None: @@ -1670,7 +1684,8 @@ def _build_direct( def _build_solver_model( model: Model, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, ) -> highspy.Highs: """Build a highspy.Highs instance that mirrors the linopy `model`.""" if model.variables.sos: @@ -1973,23 +1988,12 @@ def _register_solver_model(self, m: gurobipy.Model) -> gurobipy.Model: assert self._env_stack is not None return self._env_stack.enter_context(m) - def _detach_solver_model(self) -> gurobipy.Model: - """ - Hand ownership of the gurobipy model to the caller: unregister it from - the solver's teardown so ``close()`` does not dispose it. gurobipy - frees the underlying env once the caller drops the model. - """ - m = self.solver_model - if self._env_stack is not None: - self._env_stack.pop_all() - self.close() - return m - def _build_direct( self, explicit_coordinate_names: bool = False, env: gurobipy.Env | dict[str, Any] | None = None, - set_names: bool = True, + *, + set_names: bool, **kwargs: Any, ) -> None: model = self.model @@ -2011,7 +2015,8 @@ def _build_solver_model( model: Model, env: gurobipy.Env | None = None, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, ) -> gurobipy.Model: """Build a gurobipy.Model that mirrors the linopy `model`.""" model.constraints.sanitize_missings() @@ -2741,7 +2746,8 @@ def _apply_obj_sense(self, ctx: Any, sense: str) -> None: def _build_direct( self, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, **kwargs: Any, ) -> None: model = self.model @@ -2763,7 +2769,8 @@ def _build_direct( def _build_solver_model( model: Model, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, ) -> xpress.problem: """ Build an ``xpress.problem`` that mirrors the linopy ``model`` via ``loadproblem``. @@ -3465,15 +3472,18 @@ def _run_direct( def _build_direct( self, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, + task: mosek.Task | None = None, **kwargs: Any, ) -> None: model = self.model assert model is not None self.close() self._env_stack = contextlib.ExitStack() - env = self._env_stack.enter_context(mosek.Env()) - task = self._env_stack.enter_context(env.Task(0, 0)) + if task is None: + env = self._env_stack.enter_context(mosek.Env()) + task = self._env_stack.enter_context(env.Task(0, 0)) m = self._build_solver_model( model, task, @@ -3490,7 +3500,8 @@ def _build_solver_model( model: Model, task: mosek.Task, explicit_coordinate_names: bool = False, - set_names: bool = True, + *, + set_names: bool, ) -> mosek.Task: """Populate an empty MOSEK task with the contents of `model`.""" if model.variables.sos: diff --git a/linopy/variables.py b/linopy/variables.py index 4a97057f..0589756e 100644 --- a/linopy/variables.py +++ b/linopy/variables.py @@ -51,7 +51,6 @@ has_optimized_model, iterate_slices, save_join, - set_int_index, to_dataframe, to_polars, ) @@ -1195,6 +1194,7 @@ def get_solver_attribute(self, attr: str) -> DataArray: xr.DataArray """ from linopy.solver_capabilities import SolverFeature, solver_supports + from linopy.solvers import _solution_from_labels, _solution_from_names solver_model = self.model.solver_model if not solver_supports( @@ -1204,17 +1204,18 @@ def get_solver_attribute(self, attr: str) -> DataArray: "Solver attribute getter only supports the Gurobi solver for now." ) - vals = pd.Series( - {v.VarName: getattr(v, attr) for v in solver_model.getVars()}, dtype=float - ) - vals = set_int_index(vals) - - idx = np.ravel(self.labels) - try: - values = vals[idx].to_numpy().reshape(self.labels.shape) - except KeyError: - values = vals.reindex(idx).to_numpy().reshape(self.labels.shape) + solver = self.model.solver + assert solver is not None + gurobi_vars = solver_model.getVars() + vals = solver_model.getAttr(attr, gurobi_vars) + if solver.io_api == "direct": + lookup = _solution_from_labels(vals, solver._vlabels, solver._n_vars) + else: + names = [v.VarName for v in gurobi_vars] + lookup = _solution_from_names(vals, names, solver._n_vars) + labels = self.labels.values + values = np.where(labels == -1, np.nan, lookup[labels]) return DataArray(values, self.coords) @property diff --git a/test/test_infeasibility.py b/test/test_infeasibility.py index 2ba20d33..bee4fb33 100644 --- a/test/test_infeasibility.py +++ b/test/test_infeasibility.py @@ -247,8 +247,13 @@ def test_deprecated_method( assert len(subset) > 0 @pytest.mark.parametrize("solver", ["gurobi", "xpress", "highs"]) + @pytest.mark.parametrize("io_api,set_names", [("lp", None), ("direct", False)]) def test_masked_constraint_infeasibility( - self, solver: str, capsys: pytest.CaptureFixture[str] + self, + solver: str, + io_api: str, + set_names: bool | None, + capsys: pytest.CaptureFixture[str], ) -> None: """ Test infeasibility detection with masked constraints. @@ -275,7 +280,9 @@ def test_masked_constraint_infeasibility( m.add_constraints(x <= 4, name="x_upper", mask=mask) m.add_objective(x.sum() + y.sum()) - status, condition = m.solve(solver_name=solver) + status, condition = m.solve( + solver_name=solver, io_api=io_api, set_names=set_names + ) assert status == "warning" assert "infeasible" in condition diff --git a/test/test_io.py b/test/test_io.py index 825ca16a..67919367 100644 --- a/test/test_io.py +++ b/test/test_io.py @@ -10,6 +10,7 @@ import pickle from collections.abc import Callable from pathlib import Path +from typing import Any import numpy as np import pandas as pd @@ -608,29 +609,12 @@ def test_to_gurobipy(model: Model) -> None: assert gm.NumVars > 0 -@pytest.mark.skipif("gurobi" not in available_solvers, reason="Gurobipy not installed") -def test_to_gurobipy_no_names(model: Model) -> None: - m_with = model.to_gurobipy(set_names=True) - m_without = model.to_gurobipy(set_names=False) - names_with = [v.VarName for v in m_with.getVars()] - names_without = [v.VarName for v in m_without.getVars()] - assert names_with != names_without - - @pytest.mark.skipif("highs" not in available_solvers, reason="Highspy not installed") def test_to_highspy(model: Model) -> None: h = model.to_highspy() assert h.getLp().num_col_ > 0 -@pytest.mark.skipif("highs" not in available_solvers, reason="Highspy not installed") -def test_to_highspy_no_names(model: Model) -> None: - h = model.to_highspy(set_names=False) - lp = h.getLp() - assert len(lp.col_names_) == 0 - assert len(lp.row_names_) == 0 - - @pytest.mark.skipif("mosek" not in available_solvers, reason="Mosek not installed") def test_to_mosek(model: Model) -> None: task = model.to_mosek() @@ -644,13 +628,60 @@ def test_to_xpress(model: Model) -> None: assert p.attributes.rows > 0 -@pytest.mark.skipif("xpress" not in available_solvers, reason="Xpress not installed") -def test_to_xpress_no_names(model: Model) -> None: - p_with = model.to_xpress(set_names=True) - p_without = model.to_xpress(set_names=False) - names_with = [v.name for v in p_with.getVariable()] - names_without = [v.name for v in p_without.getVariable()] - assert names_with != names_without +def _gurobi_names(gm: Any) -> tuple[list[str], list[str]]: + gm.update() + return [v.VarName for v in gm.getVars()], [c.ConstrName for c in gm.getConstrs()] + + +SOLVER_IO: dict[ + str, tuple[Callable[..., Any], Callable[[Any], tuple[list[str], list[str]]]] +] = { + "highs": ( + Model.to_highspy, + lambda h: (list(h.getLp().col_names_), list(h.getLp().row_names_)), + ), + "gurobi": (Model.to_gurobipy, _gurobi_names), + "xpress": ( + Model.to_xpress, + lambda p: ( + [v.name for v in p.getVariable()], + [c.name for c in p.getConstraint()], + ), + ), + "mosek": ( + Model.to_mosek, + lambda task: ( + [task.getvarname(i) for i in range(task.getnumvar())], + [task.getconname(i) for i in range(task.getnumcon())], + ), + ), +} + + +@pytest.mark.parametrize("solver", SOLVER_IO) +@pytest.mark.parametrize( + "model_default,set_names,expected", + [ + (False, None, False), + (True, None, True), + (False, True, True), + (True, False, False), + ], +) +def test_to_solver_set_names( + model: Model, + solver: str, + model_default: bool, + set_names: bool | None, + expected: bool, +) -> None: + if solver not in available_solvers: + pytest.skip(f"{solver} not installed") + to_solver, names = SOLVER_IO[solver] + named = to_solver(model, set_names=True) + model.set_names_in_solver_io = model_default + built = to_solver(model, set_names=set_names) + assert (names(built) == names(named)) == expected @pytest.mark.skipif("cupdlpx" not in available_solvers, reason="cuPDLPx not installed") @@ -667,7 +698,7 @@ def test_to_cuopt(model: Model) -> None: def test_model_set_names_in_solver_io_default() -> None: - assert Model().set_names_in_solver_io is True + assert Model().set_names_in_solver_io is False @pytest.mark.skipif("highs" not in available_solvers, reason="Highspy not installed") @@ -675,9 +706,12 @@ def test_model_set_names_in_solver_io(model: Model) -> None: model.solve(solver_name="highs", io_api="direct") expected_obj = model.objective.value - model.set_names_in_solver_io = False + assert len(model.solver_model.getLp().col_names_) == 0 + + model.set_names_in_solver_io = True status, _ = model.solve(solver_name="highs", io_api="direct") assert status == "ok" + assert len(model.solver_model.getLp().col_names_) > 0 assert model.objective.value == pytest.approx(expected_obj) diff --git a/test/test_model.py b/test/test_model.py index 1b93c5c6..d3f32928 100644 --- a/test/test_model.py +++ b/test/test_model.py @@ -59,10 +59,10 @@ def test_model_config_defaults(sparse: bool) -> None: @SPARSE_CONFIG def test_model_copy_preserves_config(sparse: bool) -> None: - copied = Model(sparse=sparse, set_names_in_solver_io=False).copy() + copied = Model(sparse=sparse, set_names_in_solver_io=True).copy() assert copied.sparse is sparse assert copied.freeze_constraints is sparse - assert copied.set_names_in_solver_io is False + assert copied.set_names_in_solver_io is True def test_model_is_weakrefable() -> None: diff --git a/test/test_optimization.py b/test/test_optimization.py index 79a80db8..33184320 100644 --- a/test/test_optimization.py +++ b/test/test_optimization.py @@ -1140,6 +1140,16 @@ def test_basis_and_warmstart( ) +@pytest.mark.parametrize("solver", set_names_direct_solvers) +def test_warmstart_direct_from_lp_basis( + tmp_path: Any, model: Model, solver: str +) -> None: + basis_fn = tmp_path / "basis.bas" + model.solve(solver, basis_fn=basis_fn, io_api="lp") + status, _ = model.solve(solver, warmstart_fn=basis_fn, io_api="direct") + assert status == "ok" + + @pytest.mark.parametrize("solver,io_api,explicit_coordinate_names", params) def test_solution_fn_parent_dir_doesnt_exist( model: Model, @@ -1179,6 +1189,29 @@ def test_solver_attribute_getter( assert set(rc) == set(model.variables) +@pytest.mark.skipif("gurobi" not in direct_solvers, reason="Gurobi not available") +@pytest.mark.parametrize( + "io_api,set_names", [("lp", None), ("direct", True), ("direct", False)] +) +def test_solver_attribute_getter_maps_labels( + io_api: str, set_names: bool | None +) -> None: + m = Model() + x = m.add_variables( + 0, 10, coords=[range(4)], name="x", mask=pd.Series([True, False, True, True]) + ) + y = m.add_variables(0, 10, name="y") + m.add_constraints(x.sum() + y >= 1) + m.add_objective(x.sum() + 2 * y) + m.solve("gurobi", io_api=io_api, set_names=set_names) + obj = m.variables.get_solver_attribute("Obj") + np.testing.assert_array_equal(obj.x, [1.0, np.nan, 1.0, 1.0]) + assert obj.y.item() == 2.0 + np.testing.assert_array_equal( + x.isel(dim_0=[3, 1]).get_solver_attribute("Obj"), [1.0, np.nan] + ) + + def assert_semantically_equal_direct_solves( solved_with_names: Model, solved_without_names: Model, solver: str ) -> None: