From f21f6033a54cf5157d49964b29deaa9b660a9e2d Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 12:38:37 +0200 Subject: [PATCH 1/8] perf(matrices): cache Model.matrices for frozen models and trim the export path (#1008) --- benchmark/benchmark_sparse_export.py | 147 +++++++++++++++++++++++++++ doc/release_notes.rst | 1 + linopy/matrices.py | 44 ++++---- linopy/model.py | 44 +++++++- test/test_matrices.py | 66 ++++++++++++ 5 files changed, 279 insertions(+), 23 deletions(-) create mode 100644 benchmark/benchmark_sparse_export.py diff --git a/benchmark/benchmark_sparse_export.py b/benchmark/benchmark_sparse_export.py new file mode 100644 index 000000000..03071a2d7 --- /dev/null +++ b/benchmark/benchmark_sparse_export.py @@ -0,0 +1,147 @@ +#!/usr/bin/env python3 +""" +Benchmark the build and export of a PyPSA-like dispatch model. + +Run as ``python benchmark/benchmark_sparse_export.py {sparse,dense,frozen}``. +``sparse`` uses ``Model(sparse=True)``, ``dense`` keeps mutable constraints and +``frozen`` freezes each constraint of a dense model. Every phase reports its +wall time, its peak traced memory and the operations that densified. +""" + +from __future__ import annotations + +import argparse +import gc +import tempfile +import time +import tracemalloc +import warnings +from collections.abc import Iterator +from contextlib import contextmanager +from pathlib import Path + +import numpy as np +import pandas as pd +import xarray as xr + +import linopy +from linopy.constants import PerformanceWarning +from linopy.io import to_highspy + +MATRIX_ATTRS = ("A", "b", "c", "lb", "ub", "sense", "vlabels", "clabels") + + +class Phases: + def __init__(self) -> None: + self.rows: list[tuple[str, float, float, list[str]]] = [] + + @contextmanager + def __call__(self, name: str) -> Iterator[None]: + gc.collect() + tracemalloc.start() + start = time.perf_counter() + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + yield + elapsed = time.perf_counter() - start + peak = tracemalloc.get_traced_memory()[1] + tracemalloc.stop() + densified = [ + str(w.message)[:90] + for w in caught + if issubclass(w.category, PerformanceWarning) + ] + self.rows.append((name, elapsed, peak / 1e6, densified)) + + def report(self) -> None: + for name, elapsed, peak, densified in self.rows: + note = "DENSIFY: " + "; ".join(densified) if densified else "" + print(f"{name:28s} {elapsed * 1e3:9.0f} ms peak {peak:8.1f} MB {note}") + total = sum(row[1] for row in self.rows) + print(f"{'TOTAL':28s} {total * 1e3:9.0f} ms") + + +def run(mode: str, n_bus: int, n_t: int, n_gen: int, n_line: int, n_sto: int) -> None: + linopy.options["semantics"] = "v1" + linopy.options["warn_on_densify"] = True + rng = np.random.default_rng(0) + + bus = pd.RangeIndex(n_bus, name="bus") + t = pd.RangeIndex(n_t, name="t") + gen = pd.RangeIndex(n_gen, name="gen") + line = pd.RangeIndex(n_line, name="line") + sto = pd.RangeIndex(n_sto, name="sto") + + gen_bus = rng.integers(0, n_bus, n_gen) + gen_bus[rng.random(n_gen) < 0.4] = 0 + gen_bus = pd.Series(gen_bus, index=gen, name="bus") + sto_bus = pd.Series(np.arange(n_sto) % n_bus, index=sto, name="bus") + bus0 = np.arange(n_line) % n_bus + bus1 = (bus0 + 1 + rng.integers(0, 3, n_line)) % n_bus + incidence = np.zeros((n_line, n_bus)) + incidence[np.arange(n_line), bus0] = -1 + incidence[np.arange(n_line), bus1] = 1 + incidence = xr.DataArray(incidence, coords=[line, bus]) + demand = xr.DataArray(rng.uniform(10, 100, (n_bus, n_t)), coords=[bus, t]) + pmax = xr.DataArray(rng.uniform(0, 1, (n_gen, n_t)), coords=[gen, t]) + + freeze = {"sparse": None, "dense": False, "frozen": True}[mode] + m = linopy.Model(sparse=mode == "sparse") + phase = Phases() + + with phase("vars"): + p = m.add_variables(0, coords=[gen, t], name="p") + flow = m.add_variables(-100, 100, coords=[line, t], name="flow") + soc = m.add_variables(0, 100, coords=[sto, t], name="soc") + ch = m.add_variables(0, 50, coords=[sto, t], name="ch") + dis = m.add_variables(0, 50, coords=[sto, t], name="dis") + with phase("expr: groupby nodal supply"): + supply = p.groupby(gen_bus).sum() + with phase("expr: sto groupby"): + sto_net = (dis - ch).groupby(sto_bus).sum() + with phase("expr: flow @ incidence"): + net_flow = flow @ incidence + with phase("expr: balance sum"): + balance = supply + sto_net + net_flow + with phase("cons: balance"): + m.add_constraints(balance == demand, name="balance", freeze=freeze) + with phase("expr: soc shift"): + soc_expr = soc - 0.99 * soc.shift(t=1) - 0.95 * ch + dis / 0.95 + with phase("cons: soc"): + m.add_constraints(soc_expr == 0, name="soc", freeze=freeze) + with phase("expr+cons: p <= pmax"): + m.add_constraints(p <= pmax * 100, name="p_max", freeze=freeze) + with phase("objective"): + m.add_objective((p * 2).sum() + (ch + dis).sum()) + with phase("matrices"): + matrices = m.matrices + for attr in MATRIX_ATTRS: + getattr(matrices, attr) + with phase("to_highspy"): + to_highspy(m) + with phase("to_file(lp)"): + m.to_file(Path(tempfile.mkdtemp()) / "model.lp", progress=False) + + print( + f"=== {mode} buses={n_bus} snapshots={n_t} generators={n_gen} " + f"lines={n_line} storage={n_sto} nvars={m._xCounter} ncons={m._cCounter}" + ) + phase.report() + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("mode", choices=["sparse", "dense", "frozen"]) + parser.add_argument("--buses", type=int, default=200) + parser.add_argument("--snapshots", type=int, default=720) + parser.add_argument("--generators", type=int, default=2000) + parser.add_argument("--lines", type=int, default=300) + parser.add_argument("--storage", type=int, default=200) + args = parser.parse_args() + run( + args.mode, args.buses, args.snapshots, args.generators, args.lines, args.storage + ) + + +if __name__ == "__main__": + main() diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 42e8bee49..d7c166789 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -41,6 +41,7 @@ Upcoming Version * ``@``/``dot`` against a sparse constant no longer scales with the total number of variables in the model: each chunk of the sparse product now runs on only the variables it uses, which removes several seconds of allocation overhead on models with millions of variables. (`#990 `__) * Deprecated in favour of ``Model(sparse=True)``, each with a ``FutureWarning`` and to be removed with the legacy semantics: ``Model(freeze_constraints=...)`` and the ``Model.freeze_constraints`` setter, ``groupby(...).sum(sparse=...)`` and ``linopy.options["sparse_groupby"]``. They keep their current behaviour until then, except that ``@`` ignores ``sparse_groupby``, and netcdf files that store ``freeze_constraints`` still load. (`#976 `__) * Adding a frozen constraint from a sparse expression no longer copies the lhs matrix when every row stays active, and picks the mask and the row scaling at the active rows without expanding them over the full coordinate grid. This roughly halves the peak memory of ``add_constraints`` on a sparse model. (`#977 `__) +* ``Model.matrices`` is cached while every constraint is frozen, which is always the case in a sparse model. Any change to the variables, constraints, objective or solution rebuilds it on the next access. Building the matrices no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen blocks. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per export. (`#1008 `__) *Other* diff --git a/linopy/matrices.py b/linopy/matrices.py index bca66d70f..8203626ef 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -15,14 +15,13 @@ from numpy import ndarray from linopy import expressions -from linopy.constraints import CSRConstraint -from linopy.scaling import constraint_scaling_lookup, variable_scaling_lookup +from linopy.constraints import ConstraintBase, CSRConstraint if TYPE_CHECKING: from linopy.model import Model -def _stack(csrs: list) -> scipy.sparse.csr_array | None: +def _stack(csrs: list, pruned: bool) -> scipy.sparse.csr_array | None: """ Vertically stack CSR blocks, or None when there are none. @@ -35,7 +34,8 @@ def _stack(csrs: list) -> scipy.sparse.csr_array | None: if not csrs: return None stacked = cast(scipy.sparse.csr_array, scipy.sparse.vstack(csrs, format="csr")) - stacked.eliminate_zeros() + if not pruned: + stacked.eliminate_zeros() return stacked @@ -44,6 +44,13 @@ def _concat(arrays: list, dtype: type | None = None) -> ndarray: return np.concatenate(arrays) if arrays else np.array([], dtype=dtype) +def _row_scaling(c: ConstraintBase) -> ndarray: + """Scaling of the active rows of a constraint, in row order.""" + if isinstance(c, CSRConstraint): + return c._scaling + return c.scaling.values.ravel()[c.active_row_mask()] + + def _binval_per_row(binval: int | np.ndarray, n: int) -> ndarray: """Broadcast an indicator triggering value to one entry per active row.""" if np.ndim(binval) == 0: @@ -72,18 +79,12 @@ def __init__(self, model: Model) -> None: def _build_vars(self) -> None: m = self._parent - label_index = m.variables.label_index - self.vlabels: ndarray = label_index.vlabels - var_scaling_by_label = variable_scaling_lookup(m) - self.var_scaling: ndarray = ( - var_scaling_by_label[self.vlabels] - if len(self.vlabels) - else np.array([], dtype=float) - ) + self.vlabels: ndarray = m.variables.label_index.vlabels lb_list = [] ub_list = [] vtypes_list = [] + scaling_list = [] for name, var in m.variables.items(): labels = var.labels.values.ravel() @@ -101,7 +102,9 @@ def _build_vars(self) -> None: lb_list.append(var.lower.values.ravel()[mask]) ub_list.append(var.upper.values.ravel()[mask]) vtypes_list.append(np.full(mask.sum(), vtype)) + scaling_list.append(var.solver_scaling.values.ravel()[mask]) + self.var_scaling: ndarray = _concat(scaling_list, dtype=float) if lb_list: self.lb: ndarray = np.concatenate(lb_list) * self.var_scaling self.ub: ndarray = np.concatenate(ub_list) * self.var_scaling @@ -115,15 +118,13 @@ def _build_cons(self) -> None: m = self._parent label_index = m.variables.label_index label_to_pos = label_index.label_to_pos - con_scaling_by_label = constraint_scaling_lookup(m) unit_cols = bool((self.var_scaling == 1).all()) def scale_rows_and_cols( - csr: scipy.sparse.csr_array, con_labels: np.ndarray, b: np.ndarray + csr: scipy.sparse.csr_array, row_scaling: np.ndarray, b: np.ndarray ) -> tuple[scipy.sparse.csr_array, np.ndarray]: if csr.shape[0] == 0: return csr, b - row_scaling = con_scaling_by_label[con_labels] unit_rows = bool((row_scaling == 1).all()) if unit_rows and unit_cols: return csr, b @@ -144,8 +145,8 @@ def scale_rows_and_cols( for c in m.constraints.data.values(): if c.is_indicator: cc = c if isinstance(c, CSRConstraint) else c.freeze() - csr, con_labels, b, sense = cc.to_matrix_with_rhs(label_index) - csr, b = scale_rows_and_cols(csr, con_labels, b) + csr, _, b, sense = cc.to_matrix_with_rhs(label_index) + csr, b = scale_rows_and_cols(csr, cc._scaling, b) ind_csrs.append(csr) ind_b.append(b) ind_sense.append(sense) @@ -153,17 +154,18 @@ def scale_rows_and_cols( binval = cast("int | np.ndarray", cc._binval) ind_binval.append(_binval_per_row(binval, len(b))) else: - csr, con_labels, b, sense = c.to_matrix_with_rhs(label_index) - csr, b = scale_rows_and_cols(csr, con_labels, b) + csr, _, b, sense = c.to_matrix_with_rhs(label_index) + csr, b = scale_rows_and_cols(csr, _row_scaling(c), b) reg_csrs.append(csr) reg_b.append(b) reg_sense.append(sense) self.clabels: ndarray = m.constraints.label_index.clabels - self.A: scipy.sparse.csr_array | None = _stack(reg_csrs) + frozen = all(isinstance(c, CSRConstraint) for c in m.constraints.data.values()) + self.A: scipy.sparse.csr_array | None = _stack(reg_csrs, frozen) self.b: ndarray = _concat(reg_b) self.sense: ndarray = _concat(reg_sense, dtype=object) - self.indicator_A: scipy.sparse.csr_array | None = _stack(ind_csrs) + self.indicator_A: scipy.sparse.csr_array | None = _stack(ind_csrs, True) self.indicator_b: ndarray = _concat(ind_b) self.indicator_sense: ndarray = _concat(ind_sense, dtype=object) self.indicator_binvar: ndarray = _concat(ind_binvar, dtype=np.intp) diff --git a/linopy/model.py b/linopy/model.py index c873f6846..94a3fde7d 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -135,6 +135,10 @@ def _check_infinities(sign: Any, rhs: Any, name: str) -> None: raise ValueError(f"Constraint {name} contains incorrect infinite values.") +def _same(a: tuple, b: tuple) -> bool: + return len(a) == len(b) and all(x is y for x, y in zip(a, b)) + + class Model: """ Linear optimization model. @@ -208,6 +212,7 @@ class Model: "_piecewise_formulations", "_solver", "_sos_reformulation_state", + "_matrices", "__weakref__", ) @@ -326,6 +331,7 @@ def __init__( ) self._solver: solvers.Solver | None = None self._sos_reformulation_state: SOSReformulationResult | None = None + self._matrices: tuple[tuple, tuple, MatrixAccessor] | None = None @property def solver(self) -> solvers.Solver | None: @@ -357,10 +363,44 @@ def solver_name(self, value: str | None) -> None: raise AttributeError("solver state is managed via model.solver") self.solver = None + def _matrices_key(self) -> tuple[tuple, tuple] | None: + data = self.constraints.data + constraints = [c for c in data.values() if isinstance(c, CSRConstraint)] + if len(constraints) < len(data): + return None + objective = self.objective + identities = ( + objective, + objective._expression, + *(v._data for v in self.variables.data.values()), + *(a for c in constraints for a in (c, c._csr, c._rhs, c._dual)), + ) + values = ( + self._status, + objective._sense, + objective._scaling, + tuple(self._relaxed_registry.items()), + ) + return identities, values + @property def matrices(self) -> MatrixAccessor: - """Matrix representation of the model, computed fresh on each access.""" - return MatrixAccessor(self) + """ + Matrix representation of the model. + + The accessor is cached while every constraint is frozen and the + model is unchanged; otherwise it is rebuilt on each access. + """ + key = self._matrices_key() + if key is None: + self._matrices = None + return MatrixAccessor(self) + identities, values = key + cached = self._matrices + if cached is None or cached[1] != values or not _same(cached[0], identities): + cached = (identities, values, MatrixAccessor(self)) + self._matrices = cached + return cached[2] @property def variables(self) -> Variables: diff --git a/test/test_matrices.py b/test/test_matrices.py index 6da3eaf1d..1f22f94ee 100644 --- a/test/test_matrices.py +++ b/test/test_matrices.py @@ -5,11 +5,16 @@ @author: fabian """ +from collections.abc import Callable +from typing import Any + import numpy as np import pandas as pd +import pytest import xarray as xr from linopy import EQUAL, GREATER_EQUAL, Model +from linopy.matrices import MatrixAccessor def test_basic_matrices() -> None: @@ -98,3 +103,64 @@ def test_matrices_float_c() -> None: c = m.matrices.c assert np.all(c == np.array([1.5, 1.5])) + + +def _sparse_model() -> Model: + m = Model(sparse=True) + i = pd.RangeIndex(4, name="i") + x = m.add_variables(0, 1, coords=[i], name="x") + n = m.add_variables(0, 5, coords=[i], name="n", integer=True) + m.add_constraints(x + n >= 1, name="c") + m.add_constraints(x - n <= 1, name="d") + m.add_objective((x + 2 * n).sum()) + return m + + +MUTATIONS = { + "add_variables": lambda m: m.add_variables(coords=[m.variables["x"].indexes["i"]]), + "add_constraints": lambda m: m.add_constraints(m.variables["x"] <= 0.5), + "remove_constraints": lambda m: m.remove_constraints("c"), + "update_bounds": lambda m: m.variables["x"].update(upper=2), + "variable_scaling": lambda m: setattr(m.variables["x"], "scaling", 2), + "relax": lambda m: m.variables["n"].relax(), + "objective": lambda m: m.add_objective(-m.variables["x"].sum(), overwrite=True), + "objective_scaling": lambda m: setattr(m.objective, "scaling", 2), +} + + +@pytest.mark.v1 +def test_matrices_cached_for_frozen_model() -> None: + m = _sparse_model() + assert m.matrices is m.matrices + + +@pytest.mark.v1 +def test_matrices_not_cached_with_mutable_constraint() -> None: + m = _sparse_model() + m.add_constraints(m.variables["x"] >= 0, name="mutable", freeze=False) + assert m.matrices is not m.matrices + + +@pytest.mark.v1 +@pytest.mark.parametrize("mutate", MUTATIONS.values(), ids=MUTATIONS.keys()) +def test_matrices_cache_invalidated(mutate: Callable[[Model], Any]) -> None: + m = _sparse_model() + cached = m.matrices + mutate(m) + fresh = MatrixAccessor(m) + assert m.matrices is not cached + assert m.matrices is m.matrices + for attr in ("vlabels", "clabels", "lb", "ub", "vtypes", "b", "sense", "c"): + np.testing.assert_array_equal(getattr(m.matrices, attr), getattr(fresh, attr)) + assert m.matrices.A is not None and fresh.A is not None + np.testing.assert_array_equal(m.matrices.A.toarray(), fresh.A.toarray()) + + +@pytest.mark.v1 +def test_matrices_cache_refreshes_solution() -> None: + m = _sparse_model() + m._mock_solve() + first = m.matrices.sol + for var in m.variables.data.values(): + var.solution = var.solution + 1 + np.testing.assert_array_equal(m.matrices.sol, first + 1) From f9058a6a873458b2fa7aeea51fe9d4ce5a2e2ca1 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 12:51:01 +0200 Subject: [PATCH 2/8] fix(xpress): stop editing cached matrix bounds in place; review follow-ups --- benchmark/benchmark_sparse_export.py | 15 ++++++------- doc/release_notes.rst | 2 +- linopy/matrices.py | 24 ++++++--------------- linopy/model.py | 24 +++++++++++++++------ linopy/solvers.py | 6 ++---- test/test_matrices.py | 16 ++++++++++---- test/test_scaling.py | 32 ++++++++++++++++++++++++++++ 7 files changed, 77 insertions(+), 42 deletions(-) diff --git a/benchmark/benchmark_sparse_export.py b/benchmark/benchmark_sparse_export.py index 03071a2d7..a3e051e36 100644 --- a/benchmark/benchmark_sparse_export.py +++ b/benchmark/benchmark_sparse_export.py @@ -5,7 +5,8 @@ Run as ``python benchmark/benchmark_sparse_export.py {sparse,dense,frozen}``. ``sparse`` uses ``Model(sparse=True)``, ``dense`` keeps mutable constraints and ``frozen`` freezes each constraint of a dense model. Every phase reports its -wall time, its peak traced memory and the operations that densified. +wall time, its peak traced memory and the operations that densified. In the +``sparse`` and ``frozen`` modes, ``to_highspy`` reuses the cached matrices. """ from __future__ import annotations @@ -28,8 +29,6 @@ from linopy.constants import PerformanceWarning from linopy.io import to_highspy -MATRIX_ATTRS = ("A", "b", "c", "lb", "ub", "sense", "vlabels", "clabels") - class Phases: def __init__(self) -> None: @@ -114,17 +113,15 @@ def run(mode: str, n_bus: int, n_t: int, n_gen: int, n_line: int, n_sto: int) -> with phase("objective"): m.add_objective((p * 2).sum() + (ch + dis).sum()) with phase("matrices"): - matrices = m.matrices - for attr in MATRIX_ATTRS: - getattr(matrices, attr) + m.matrices.c with phase("to_highspy"): to_highspy(m) - with phase("to_file(lp)"): - m.to_file(Path(tempfile.mkdtemp()) / "model.lp", progress=False) + with phase("to_file(lp)"), tempfile.TemporaryDirectory() as tmp: + m.to_file(Path(tmp) / "model.lp", progress=False) print( f"=== {mode} buses={n_bus} snapshots={n_t} generators={n_gen} " - f"lines={n_line} storage={n_sto} nvars={m._xCounter} ncons={m._cCounter}" + f"lines={n_line} storage={n_sto} nvars={m.nvars} ncons={m.ncons}" ) phase.report() diff --git a/doc/release_notes.rst b/doc/release_notes.rst index d7c166789..5dc7c6e5f 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -41,7 +41,7 @@ Upcoming Version * ``@``/``dot`` against a sparse constant no longer scales with the total number of variables in the model: each chunk of the sparse product now runs on only the variables it uses, which removes several seconds of allocation overhead on models with millions of variables. (`#990 `__) * Deprecated in favour of ``Model(sparse=True)``, each with a ``FutureWarning`` and to be removed with the legacy semantics: ``Model(freeze_constraints=...)`` and the ``Model.freeze_constraints`` setter, ``groupby(...).sum(sparse=...)`` and ``linopy.options["sparse_groupby"]``. They keep their current behaviour until then, except that ``@`` ignores ``sparse_groupby``, and netcdf files that store ``freeze_constraints`` still load. (`#976 `__) * Adding a frozen constraint from a sparse expression no longer copies the lhs matrix when every row stays active, and picks the mask and the row scaling at the active rows without expanding them over the full coordinate grid. This roughly halves the peak memory of ``add_constraints`` on a sparse model. (`#977 `__) -* ``Model.matrices`` is cached while every constraint is frozen, which is always the case in a sparse model. Any change to the variables, constraints, objective or solution rebuilds it on the next access. Building the matrices no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen blocks. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per export. (`#1008 `__) +* ``Model.matrices`` is cached while every constraint is frozen, which is always the case in a sparse model. Any change to the variables, constraints, objective or solution rebuilds it on the next access. Building the matrices no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen blocks. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per export. The cached matrices stay in memory until the model changes. (`#1008 `__) *Other* diff --git a/linopy/matrices.py b/linopy/matrices.py index 8203626ef..08e328e41 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -21,22 +21,11 @@ from linopy.model import Model -def _stack(csrs: list, pruned: bool) -> scipy.sparse.csr_array | None: - """ - Vertically stack CSR blocks, or None when there are none. - - Explicit zeros are dropped: expressions that broadcast against a dense - coordinate store one coefficient per pair, most of them zero, and a zero - coefficient never changes a constraint. Keeping them only inflates the - stored nnz handed to the solvers/writers (e.g. ``highspy.addRows`` scales - with stored nnz), so we prune them once, centrally, for every backend. - """ +def _stack(csrs: list[scipy.sparse.csr_array]) -> scipy.sparse.csr_array | None: + """Vertically stack CSR blocks, or None when there are none.""" if not csrs: return None - stacked = cast(scipy.sparse.csr_array, scipy.sparse.vstack(csrs, format="csr")) - if not pruned: - stacked.eliminate_zeros() - return stacked + return cast(scipy.sparse.csr_array, scipy.sparse.vstack(csrs, format="csr")) def _concat(arrays: list, dtype: type | None = None) -> ndarray: @@ -155,17 +144,18 @@ def scale_rows_and_cols( ind_binval.append(_binval_per_row(binval, len(b))) else: csr, _, b, sense = c.to_matrix_with_rhs(label_index) + if not isinstance(c, CSRConstraint): + csr.eliminate_zeros() csr, b = scale_rows_and_cols(csr, _row_scaling(c), b) reg_csrs.append(csr) reg_b.append(b) reg_sense.append(sense) self.clabels: ndarray = m.constraints.label_index.clabels - frozen = all(isinstance(c, CSRConstraint) for c in m.constraints.data.values()) - self.A: scipy.sparse.csr_array | None = _stack(reg_csrs, frozen) + self.A: scipy.sparse.csr_array | None = _stack(reg_csrs) self.b: ndarray = _concat(reg_b) self.sense: ndarray = _concat(reg_sense, dtype=object) - self.indicator_A: scipy.sparse.csr_array | None = _stack(ind_csrs, True) + self.indicator_A: scipy.sparse.csr_array | None = _stack(ind_csrs) self.indicator_b: ndarray = _concat(ind_b) self.indicator_sense: ndarray = _concat(ind_sense, dtype=object) self.indicator_binvar: ndarray = _concat(ind_binvar, dtype=np.intp) diff --git a/linopy/model.py b/linopy/model.py index 94a3fde7d..f72731c2d 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -14,7 +14,7 @@ from pathlib import Path from tempfile import NamedTemporaryFile, gettempdir from types import MappingProxyType -from typing import TYPE_CHECKING, Any, Literal, get_args, overload +from typing import TYPE_CHECKING, Any, Literal, NamedTuple, get_args, overload from warnings import warn import numpy as np @@ -135,7 +135,13 @@ def _check_infinities(sign: Any, rhs: Any, name: str) -> None: raise ValueError(f"Constraint {name} contains incorrect infinite values.") -def _same(a: tuple, b: tuple) -> bool: +class _MatricesCache(NamedTuple): + identities: tuple[object, ...] + values: tuple[object, ...] + accessor: MatrixAccessor + + +def _same(a: tuple[object, ...], b: tuple[object, ...]) -> bool: return len(a) == len(b) and all(x is y for x, y in zip(a, b)) @@ -331,7 +337,7 @@ def __init__( ) self._solver: solvers.Solver | None = None self._sos_reformulation_state: SOSReformulationResult | None = None - self._matrices: tuple[tuple, tuple, MatrixAccessor] | None = None + self._matrices: _MatricesCache | None = None @property def solver(self) -> solvers.Solver | None: @@ -363,7 +369,7 @@ def solver_name(self, value: str | None) -> None: raise AttributeError("solver state is managed via model.solver") self.solver = None - def _matrices_key(self) -> tuple[tuple, tuple] | None: + def _matrices_key(self) -> tuple[tuple[object, ...], tuple[object, ...]] | None: data = self.constraints.data constraints = [c for c in data.values() if isinstance(c, CSRConstraint)] if len(constraints) < len(data): @@ -397,10 +403,14 @@ def matrices(self) -> MatrixAccessor: return MatrixAccessor(self) identities, values = key cached = self._matrices - if cached is None or cached[1] != values or not _same(cached[0], identities): - cached = (identities, values, MatrixAccessor(self)) + if ( + cached is None + or cached.values != values + or not _same(cached.identities, identities) + ): + cached = _MatricesCache(identities, values, MatrixAccessor(self)) self._matrices = cached - return cached[2] + return cached.accessor @property def variables(self) -> Variables: diff --git a/linopy/solvers.py b/linopy/solvers.py index 6e1ec74d0..600e3d6d6 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -2799,10 +2799,8 @@ def _build_solver_model( rowind = np.empty(0, dtype=np.int64) rowcoef = np.empty(0, dtype=float) - lb = np.asarray(M.lb, dtype=float) - ub = np.asarray(M.ub, dtype=float) - np.place(lb, np.isneginf(lb), -xpress.infinity) - np.place(ub, np.isposinf(ub), xpress.infinity) + lb = np.where(np.isneginf(M.lb), -xpress.infinity, M.lb) + ub = np.where(np.isposinf(M.ub), xpress.infinity, M.ub) rowtype: np.ndarray rhs: np.ndarray diff --git a/test/test_matrices.py b/test/test_matrices.py index 1f22f94ee..918af5845 100644 --- a/test/test_matrices.py +++ b/test/test_matrices.py @@ -73,7 +73,8 @@ def test_matrices_duplicated_variables() -> None: assert np.isin(np.unique(np.array(A)), [0.0, 2.0]).all() -def test_matrices_drops_explicit_zeros() -> None: +@pytest.mark.parametrize("freeze", [False, True]) +def test_matrices_drops_explicit_zeros(freeze: bool) -> None: # https://github.com/PyPSA/linopy/issues/814 # Expressions that broadcast against a dense coordinate store one coefficient # per pair, most of them structurally zero. Those must not reach A, whose @@ -84,7 +85,9 @@ def test_matrices_drops_explicit_zeros() -> None: coeff = xr.DataArray( np.eye(4), dims=["j", "i"], coords={"j": range(4), "i": range(4)} ) - m.add_constraints((coeff * x.rename(dim_0="i")).sum("i") <= 1, name="c") + m.add_constraints( + (coeff * x.rename(dim_0="i")).sum("i") <= 1, name="c", freeze=freeze + ) A = m.matrices.A assert A is not None @@ -159,8 +162,13 @@ def test_matrices_cache_invalidated(mutate: Callable[[Model], Any]) -> None: @pytest.mark.v1 def test_matrices_cache_refreshes_solution() -> None: m = _sparse_model() + cached = m.matrices m._mock_solve() - first = m.matrices.sol + assert m.matrices is not cached + sol, dual = m.matrices.sol, m.matrices.dual for var in m.variables.data.values(): var.solution = var.solution + 1 - np.testing.assert_array_equal(m.matrices.sol, first + 1) + for con in m.constraints.data.values(): + con.dual = con.dual + 1 + np.testing.assert_array_equal(m.matrices.sol, sol + 1) + np.testing.assert_array_equal(m.matrices.dual, dual + 1) diff --git a/test/test_scaling.py b/test/test_scaling.py index 16c318421..709172c99 100644 --- a/test/test_scaling.py +++ b/test/test_scaling.py @@ -19,6 +19,38 @@ def _dense(matrix: Any) -> np.ndarray: return matrix.toarray() +@pytest.mark.parametrize("freeze", [False, True]) +def test_masked_scaling_in_matrices(freeze: bool) -> None: + m = Model() + i = pd.RangeIndex(3, name="i") + mask = xr.DataArray([True, False, True], coords=[i]) + x = m.add_variables( + 0, + 1, + coords=[i], + name="x", + mask=mask, + scaling=xr.DataArray([2.0, 3.0, 4.0], coords=[i]), + ) + y = m.add_variables(0, 1, coords=[i], name="y") + m.add_constraints( + x.where(mask) + y >= 1, + name="c", + mask=mask, + scaling=xr.DataArray([5.0, 6.0, 7.0], coords=[i]), + freeze=freeze, + ) + + matrices = m.matrices + + np.testing.assert_allclose(matrices.var_scaling, [2.0, 4.0, 1.0, 1.0, 1.0]) + np.testing.assert_allclose(matrices.b, [5.0, 7.0]) + np.testing.assert_allclose( + _dense(matrices.A), + [[5 / 2, 0.0, 5.0, 0.0, 0.0], [0.0, 7 / 4, 0.0, 0.0, 7.0]], + ) + + def test_variable_constraint_and_objective_scaling_in_matrices() -> None: m = Model() i = pd.Index(["a", "b"], name="i") From 2c9dee42db8f3abdd758c22eb2c1427b54368443 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 14:14:07 +0200 Subject: [PATCH 3/8] docs: note the shared Model.matrices contract in the release notes --- doc/release_notes.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 5dc7c6e5f..2303c33af 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -41,7 +41,7 @@ Upcoming Version * ``@``/``dot`` against a sparse constant no longer scales with the total number of variables in the model: each chunk of the sparse product now runs on only the variables it uses, which removes several seconds of allocation overhead on models with millions of variables. (`#990 `__) * Deprecated in favour of ``Model(sparse=True)``, each with a ``FutureWarning`` and to be removed with the legacy semantics: ``Model(freeze_constraints=...)`` and the ``Model.freeze_constraints`` setter, ``groupby(...).sum(sparse=...)`` and ``linopy.options["sparse_groupby"]``. They keep their current behaviour until then, except that ``@`` ignores ``sparse_groupby``, and netcdf files that store ``freeze_constraints`` still load. (`#976 `__) * Adding a frozen constraint from a sparse expression no longer copies the lhs matrix when every row stays active, and picks the mask and the row scaling at the active rows without expanding them over the full coordinate grid. This roughly halves the peak memory of ``add_constraints`` on a sparse model. (`#977 `__) -* ``Model.matrices`` is cached while every constraint is frozen, which is always the case in a sparse model. Any change to the variables, constraints, objective or solution rebuilds it on the next access. Building the matrices no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen blocks. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per export. The cached matrices stay in memory until the model changes. (`#1008 `__) +* ``Model.matrices`` is cached while every constraint is frozen, which is always the case in a sparse model. Any change to the variables, constraints, objective or solution rebuilds it on the next access. Building the matrices no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen blocks. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per export. The cached matrices stay in memory until the model changes. In a fully frozen model, ``m.matrices`` returns the same object on repeated access, so its arrays must not be edited in place. (`#1008 `__) *Other* From fbf7c7497859a89b76b2eb8e9bdf37e259fd54b7 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 21:02:17 +0200 Subject: [PATCH 4/8] fix(matrices): key the cache on the objective expression data --- linopy/model.py | 6 +++++- test/test_matrices.py | 6 ++++++ 2 files changed, 11 insertions(+), 1 deletion(-) diff --git a/linopy/model.py b/linopy/model.py index f72731c2d..e00cddaca 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -375,9 +375,13 @@ def _matrices_key(self) -> tuple[tuple[object, ...], tuple[object, ...]] | None: if len(constraints) < len(data): return None objective = self.objective + expression = objective._expression + csr = expression._csr if isinstance(expression, LinearExpression) else None identities = ( objective, - objective._expression, + expression, + expression._data, + csr, *(v._data for v in self.variables.data.values()), *(a for c in constraints for a in (c, c._csr, c._rhs, c._dual)), ) diff --git a/test/test_matrices.py b/test/test_matrices.py index 918af5845..533b27ac2 100644 --- a/test/test_matrices.py +++ b/test/test_matrices.py @@ -128,6 +128,12 @@ def _sparse_model() -> Model: "relax": lambda m: m.variables["n"].relax(), "objective": lambda m: m.add_objective(-m.variables["x"].sum(), overwrite=True), "objective_scaling": lambda m: setattr(m.objective, "scaling", 2), + "objective_coeffs": lambda m: setattr( + m.objective.expression, "coeffs", m.objective.expression.coeffs * 2 + ), + "objective_vars": lambda m: setattr( + m.objective.expression, "vars", m.objective.expression.vars.roll(_term=1) + ), } From 087292b7924c5d76721958d29c45219a6929bc79 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 21:55:33 +0200 Subject: [PATCH 5/8] refactor(matrices): drop the Model.matrices cache, build the matrices once per caller Model.matrices returns a fresh accessor again. Dualization binds it once and expression solutions map labels directly, without building the matrices. --- benchmark/benchmark_sparse_export.py | 4 +- doc/release_notes.rst | 2 +- linopy/dualization.py | 9 ++-- linopy/expressions.py | 8 +-- linopy/model.py | 60 ++-------------------- test/test_expressions.py | 13 +++++ test/test_matrices.py | 76 ---------------------------- 7 files changed, 28 insertions(+), 144 deletions(-) diff --git a/benchmark/benchmark_sparse_export.py b/benchmark/benchmark_sparse_export.py index a3e051e36..9b0590f39 100644 --- a/benchmark/benchmark_sparse_export.py +++ b/benchmark/benchmark_sparse_export.py @@ -5,8 +5,8 @@ Run as ``python benchmark/benchmark_sparse_export.py {sparse,dense,frozen}``. ``sparse`` uses ``Model(sparse=True)``, ``dense`` keeps mutable constraints and ``frozen`` freezes each constraint of a dense model. Every phase reports its -wall time, its peak traced memory and the operations that densified. In the -``sparse`` and ``frozen`` modes, ``to_highspy`` reuses the cached matrices. +wall time, its peak traced memory and the operations that densified. The +``matrices``, ``to_highspy`` and ``to_file(lp)`` phases each build the matrices. """ from __future__ import annotations diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 2303c33af..97ae738fe 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -41,7 +41,7 @@ Upcoming Version * ``@``/``dot`` against a sparse constant no longer scales with the total number of variables in the model: each chunk of the sparse product now runs on only the variables it uses, which removes several seconds of allocation overhead on models with millions of variables. (`#990 `__) * Deprecated in favour of ``Model(sparse=True)``, each with a ``FutureWarning`` and to be removed with the legacy semantics: ``Model(freeze_constraints=...)`` and the ``Model.freeze_constraints`` setter, ``groupby(...).sum(sparse=...)`` and ``linopy.options["sparse_groupby"]``. They keep their current behaviour until then, except that ``@`` ignores ``sparse_groupby``, and netcdf files that store ``freeze_constraints`` still load. (`#976 `__) * Adding a frozen constraint from a sparse expression no longer copies the lhs matrix when every row stays active, and picks the mask and the row scaling at the active rows without expanding them over the full coordinate grid. This roughly halves the peak memory of ``add_constraints`` on a sparse model. (`#977 `__) -* ``Model.matrices`` is cached while every constraint is frozen, which is always the case in a sparse model. Any change to the variables, constraints, objective or solution rebuilds it on the next access. Building the matrices no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen blocks. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per export. The cached matrices stay in memory until the model changes. In a fully frozen model, ``m.matrices`` returns the same object on repeated access, so its arrays must not be edited in place. (`#1008 `__) +* Building ``Model.matrices`` no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen constraints. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per build. (`#1008 `__) *Other* diff --git a/linopy/dualization.py b/linopy/dualization.py index edcfffdf0..27bf68e72 100644 --- a/linopy/dualization.py +++ b/linopy/dualization.py @@ -479,12 +479,13 @@ def _add_dual_feasibility_constraints( dual_vars : dict ``{constraint_name: dual_variable}`` as returned by ``_add_dual_variables()``. """ - A = m.matrices.A + M = m.matrices + A = M.A if A is None: raise ValueError("Constraint matrix is None, model has no constraints.") - vlabels = np.asarray(m.matrices.vlabels, dtype=np.int64) - clabels = np.asarray(m.matrices.clabels, dtype=np.int64) + vlabels = np.asarray(M.vlabels, dtype=np.int64) + clabels = np.asarray(M.clabels, dtype=np.int64) flat_con_to_dual = _build_flat_con_to_dual_label_lookup(m, dual_vars) if not len(flat_con_to_dual): @@ -496,7 +497,7 @@ def _add_dual_feasibility_constraints( flat_v, flat_d, nnz_data = _extract_dual_feas_entries( A, vlabels, clabels, flat_con_to_dual ) - c_lookup = _build_obj_coeff_lookup(vlabels, m.matrices.c) + c_lookup = _build_obj_coeff_lookup(vlabels, M.c) logger.debug("Building dual feasibility constraints for each primal variable.") for var_name, var in m.variables.items(): diff --git a/linopy/expressions.py b/linopy/expressions.py index 083e673f2..881b457e1 100644 --- a/linopy/expressions.py +++ b/linopy/expressions.py @@ -1743,11 +1743,11 @@ def _map_solution(self) -> DataArray: Replace variable labels by solution values. """ m = self.model - M = m.matrices - sol = pd.Series(M.sol, M.vlabels) + sol = np.full(m._xCounter + 1, np.nan) + for var in m.variables.data.values(): + sol[var.labels.values] = var.solution.values sol[-1] = np.nan - idx = np.ravel(self.vars) - values = np.asarray(sol[idx]).reshape(self.vars.shape) + values = sol[self.vars.values] return xr.DataArray(values, dims=self.vars.dims, coords=self.vars.coords) @property diff --git a/linopy/model.py b/linopy/model.py index e00cddaca..c873f6846 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -14,7 +14,7 @@ from pathlib import Path from tempfile import NamedTemporaryFile, gettempdir from types import MappingProxyType -from typing import TYPE_CHECKING, Any, Literal, NamedTuple, get_args, overload +from typing import TYPE_CHECKING, Any, Literal, get_args, overload from warnings import warn import numpy as np @@ -135,16 +135,6 @@ def _check_infinities(sign: Any, rhs: Any, name: str) -> None: raise ValueError(f"Constraint {name} contains incorrect infinite values.") -class _MatricesCache(NamedTuple): - identities: tuple[object, ...] - values: tuple[object, ...] - accessor: MatrixAccessor - - -def _same(a: tuple[object, ...], b: tuple[object, ...]) -> bool: - return len(a) == len(b) and all(x is y for x, y in zip(a, b)) - - class Model: """ Linear optimization model. @@ -218,7 +208,6 @@ class Model: "_piecewise_formulations", "_solver", "_sos_reformulation_state", - "_matrices", "__weakref__", ) @@ -337,7 +326,6 @@ def __init__( ) self._solver: solvers.Solver | None = None self._sos_reformulation_state: SOSReformulationResult | None = None - self._matrices: _MatricesCache | None = None @property def solver(self) -> solvers.Solver | None: @@ -369,52 +357,10 @@ def solver_name(self, value: str | None) -> None: raise AttributeError("solver state is managed via model.solver") self.solver = None - def _matrices_key(self) -> tuple[tuple[object, ...], tuple[object, ...]] | None: - data = self.constraints.data - constraints = [c for c in data.values() if isinstance(c, CSRConstraint)] - if len(constraints) < len(data): - return None - objective = self.objective - expression = objective._expression - csr = expression._csr if isinstance(expression, LinearExpression) else None - identities = ( - objective, - expression, - expression._data, - csr, - *(v._data for v in self.variables.data.values()), - *(a for c in constraints for a in (c, c._csr, c._rhs, c._dual)), - ) - values = ( - self._status, - objective._sense, - objective._scaling, - tuple(self._relaxed_registry.items()), - ) - return identities, values - @property def matrices(self) -> MatrixAccessor: - """ - Matrix representation of the model. - - The accessor is cached while every constraint is frozen and the - model is unchanged; otherwise it is rebuilt on each access. - """ - key = self._matrices_key() - if key is None: - self._matrices = None - return MatrixAccessor(self) - identities, values = key - cached = self._matrices - if ( - cached is None - or cached.values != values - or not _same(cached.identities, identities) - ): - cached = _MatricesCache(identities, values, MatrixAccessor(self)) - self._matrices = cached - return cached.accessor + """Matrix representation of the model, computed fresh on each access.""" + return MatrixAccessor(self) @property def variables(self) -> Variables: diff --git a/test/test_expressions.py b/test/test_expressions.py index 53c5a57ce..64308ae0b 100644 --- a/test/test_expressions.py +++ b/test/test_expressions.py @@ -141,3 +141,16 @@ def test_expressions_solution() -> None: assert isinstance(sol, xr.Dataset) assert "double_x" in sol assert (sol["double_x"] == 4).all() + + +def test_expression_solution_maps_variable_values() -> None: + m = Model() + y = m.add_variables(coords=[pd.RangeIndex(3, name="i")], name="y") + z = m.add_variables(coords=[pd.RangeIndex(2, name="j")], name="z") + m._mock_solve() + y.solution = xr.DataArray([10.0, 20.0, 30.0], coords=y.coords) + z.solution = xr.DataArray([1.0, 2.0], coords=z.coords) + expected = 2 * y.solution + z.solution + xr.testing.assert_equal((2 * y + z).solution, expected.rename("solution")) + expected = y.solution * z.solution + xr.testing.assert_equal((y * z).solution, expected.rename("solution")) diff --git a/test/test_matrices.py b/test/test_matrices.py index 533b27ac2..fb93c5a37 100644 --- a/test/test_matrices.py +++ b/test/test_matrices.py @@ -5,16 +5,12 @@ @author: fabian """ -from collections.abc import Callable -from typing import Any - import numpy as np import pandas as pd import pytest import xarray as xr from linopy import EQUAL, GREATER_EQUAL, Model -from linopy.matrices import MatrixAccessor def test_basic_matrices() -> None: @@ -106,75 +102,3 @@ def test_matrices_float_c() -> None: c = m.matrices.c assert np.all(c == np.array([1.5, 1.5])) - - -def _sparse_model() -> Model: - m = Model(sparse=True) - i = pd.RangeIndex(4, name="i") - x = m.add_variables(0, 1, coords=[i], name="x") - n = m.add_variables(0, 5, coords=[i], name="n", integer=True) - m.add_constraints(x + n >= 1, name="c") - m.add_constraints(x - n <= 1, name="d") - m.add_objective((x + 2 * n).sum()) - return m - - -MUTATIONS = { - "add_variables": lambda m: m.add_variables(coords=[m.variables["x"].indexes["i"]]), - "add_constraints": lambda m: m.add_constraints(m.variables["x"] <= 0.5), - "remove_constraints": lambda m: m.remove_constraints("c"), - "update_bounds": lambda m: m.variables["x"].update(upper=2), - "variable_scaling": lambda m: setattr(m.variables["x"], "scaling", 2), - "relax": lambda m: m.variables["n"].relax(), - "objective": lambda m: m.add_objective(-m.variables["x"].sum(), overwrite=True), - "objective_scaling": lambda m: setattr(m.objective, "scaling", 2), - "objective_coeffs": lambda m: setattr( - m.objective.expression, "coeffs", m.objective.expression.coeffs * 2 - ), - "objective_vars": lambda m: setattr( - m.objective.expression, "vars", m.objective.expression.vars.roll(_term=1) - ), -} - - -@pytest.mark.v1 -def test_matrices_cached_for_frozen_model() -> None: - m = _sparse_model() - assert m.matrices is m.matrices - - -@pytest.mark.v1 -def test_matrices_not_cached_with_mutable_constraint() -> None: - m = _sparse_model() - m.add_constraints(m.variables["x"] >= 0, name="mutable", freeze=False) - assert m.matrices is not m.matrices - - -@pytest.mark.v1 -@pytest.mark.parametrize("mutate", MUTATIONS.values(), ids=MUTATIONS.keys()) -def test_matrices_cache_invalidated(mutate: Callable[[Model], Any]) -> None: - m = _sparse_model() - cached = m.matrices - mutate(m) - fresh = MatrixAccessor(m) - assert m.matrices is not cached - assert m.matrices is m.matrices - for attr in ("vlabels", "clabels", "lb", "ub", "vtypes", "b", "sense", "c"): - np.testing.assert_array_equal(getattr(m.matrices, attr), getattr(fresh, attr)) - assert m.matrices.A is not None and fresh.A is not None - np.testing.assert_array_equal(m.matrices.A.toarray(), fresh.A.toarray()) - - -@pytest.mark.v1 -def test_matrices_cache_refreshes_solution() -> None: - m = _sparse_model() - cached = m.matrices - m._mock_solve() - assert m.matrices is not cached - sol, dual = m.matrices.sol, m.matrices.dual - for var in m.variables.data.values(): - var.solution = var.solution + 1 - for con in m.constraints.data.values(): - con.dual = con.dual + 1 - np.testing.assert_array_equal(m.matrices.sol, sol + 1) - np.testing.assert_array_equal(m.matrices.dual, dual + 1) From d757470c2659464aa7c36c5d7a1e1bca2a159351 Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 22:08:26 +0200 Subject: [PATCH 6/8] fix(expressions): raise when a solution maps labels missing from the model --- linopy/expressions.py | 6 +++++- test/test_expressions.py | 3 +++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/linopy/expressions.py b/linopy/expressions.py index 881b457e1..d2f084635 100644 --- a/linopy/expressions.py +++ b/linopy/expressions.py @@ -1743,11 +1743,15 @@ def _map_solution(self) -> DataArray: Replace variable labels by solution values. """ m = self.model + labels = self.vars.values + known = np.append(m.variables.label_index.label_to_pos != -1, True) + if not known[labels].all(): + raise KeyError("Expression references variables missing from the model.") sol = np.full(m._xCounter + 1, np.nan) for var in m.variables.data.values(): sol[var.labels.values] = var.solution.values sol[-1] = np.nan - values = sol[self.vars.values] + values = sol[labels] return xr.DataArray(values, dims=self.vars.dims, coords=self.vars.coords) @property diff --git a/test/test_expressions.py b/test/test_expressions.py index 64308ae0b..90cc15386 100644 --- a/test/test_expressions.py +++ b/test/test_expressions.py @@ -154,3 +154,6 @@ def test_expression_solution_maps_variable_values() -> None: xr.testing.assert_equal((2 * y + z).solution, expected.rename("solution")) expected = y.solution * z.solution xr.testing.assert_equal((y * z).solution, expected.rename("solution")) + m.remove_variables("z") + with pytest.raises(KeyError): + (2 * y + z).solution From 89560f39e13567577bdfc17e40506506d54f595f Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 22:19:02 +0200 Subject: [PATCH 7/8] docs: restate the matrices build saving without the cache --- doc/release_notes.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 97ae738fe..45661be79 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -41,7 +41,7 @@ Upcoming Version * ``@``/``dot`` against a sparse constant no longer scales with the total number of variables in the model: each chunk of the sparse product now runs on only the variables it uses, which removes several seconds of allocation overhead on models with millions of variables. (`#990 `__) * Deprecated in favour of ``Model(sparse=True)``, each with a ``FutureWarning`` and to be removed with the legacy semantics: ``Model(freeze_constraints=...)`` and the ``Model.freeze_constraints`` setter, ``groupby(...).sum(sparse=...)`` and ``linopy.options["sparse_groupby"]``. They keep their current behaviour until then, except that ``@`` ignores ``sparse_groupby``, and netcdf files that store ``freeze_constraints`` still load. (`#976 `__) * Adding a frozen constraint from a sparse expression no longer copies the lhs matrix when every row stays active, and picks the mask and the row scaling at the active rows without expanding them over the full coordinate grid. This roughly halves the peak memory of ``add_constraints`` on a sparse model. (`#977 `__) -* Building ``Model.matrices`` no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen constraints. On a model with 2M variables this saves about 200 ms and 185 MB of peak memory per build. (`#1008 `__) +* Building ``Model.matrices`` no longer allocates label-sized scaling lookups and skips ``eliminate_zeros`` for frozen constraints. On a model with 2M variables this makes each build of a sparse or frozen model about 20 to 35 ms faster. (`#1008 `__) *Other* From eeb63e780f903a5134095ee5fabf5edd110a448f Mon Sep 17 00:00:00 2001 From: Fabian Date: Tue, 6 Oct 2026 22:52:51 +0200 Subject: [PATCH 8/8] test(matrices): cover the solution read-back of the accessor --- test/test_matrices.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/test_matrices.py b/test/test_matrices.py index fb93c5a37..04317f8aa 100644 --- a/test/test_matrices.py +++ b/test/test_matrices.py @@ -102,3 +102,19 @@ def test_matrices_float_c() -> None: c = m.matrices.c assert np.all(c == np.array([1.5, 1.5])) + + +def test_matrices_sol_aligned_with_vlabels() -> None: + m = Model() + i = pd.RangeIndex(3, name="i") + x = m.add_variables(coords=[i], name="x", mask=pd.Series([True, False, True], i)) + y = m.add_variables(coords=[i], name="y") + m.add_constraints(x + y >= 0, name="c") + with pytest.raises(ValueError, match="not optimized"): + m.matrices.sol + m._mock_solve() + x.solution = xr.DataArray([1.0, np.nan, 3.0], coords=[i]) + y.solution = xr.DataArray([4.0, 5.0, 6.0], coords=[i]) + M = m.matrices + np.testing.assert_array_equal(M.vlabels, [0, 2, 3, 4, 5]) + np.testing.assert_array_equal(M.sol, [1.0, 3.0, 4.0, 5.0, 6.0])