From 1ad1a818ec6bb63d91cb8eb00c67e05b9e89c66f Mon Sep 17 00:00:00 2001 From: Olivier Lacroix Date: Tue, 20 Oct 2026 17:38:31 +1100 Subject: [PATCH 1/2] feat(restrict): window a model's flattened matrices `Model.restrict` returns a `WindowModel` backed by a slice of the parent's `A`, `b`, `c`, `lb` and `ub` instead of a rebuild of every variable and constraint. A window is a row slice plus a column range; boundary terms are folded at their actual coefficient, and containers are materialised only when their values are read. - `WindowModel.matrices` exposes the window arrays; `solution`, `dual` and `objective` read them back with column, row and objective scaling undone. - `WindowModel.solve` builds the native model from the window and moves it in place through `ModelDiff.from_windows` and the existing per-backend apply hooks. - `advance` re-slices the parent; a shape change rebuilds the solver. - A window that cuts a coupling without `fix` raises `RestrictionError`, naming the constraints. Supported: one selected dimension, dense and CSR parents, masked rows. Not supported: indicators, quadratic objectives, SOS, piecewise, several selected dimensions, subset carriers. `linopy/restrict` is at 100% statement and branch coverage. Co-Authored-By: Claude Code --- doc/api.rst | 1 + doc/index.rst | 1 + doc/release_notes.rst | 1 + doc/rolling-horizon.rst | 127 +++++++ linopy/examples.py | 4 + linopy/matrices.py | 25 +- linopy/model.py | 64 ++++ linopy/persistent/diff.py | 89 ++++- linopy/restrict/__init__.py | 40 +++ linopy/restrict/api.py | 48 +++ linopy/restrict/errors.py | 13 + linopy/restrict/matrix.py | 364 ++++++++++++++++++++ linopy/restrict/window.py | 375 +++++++++++++++++++++ linopy/solvers.py | 70 ++++ test/test_persistent_snapshot_diff.py | 3 + test/test_restrict_matrix.py | 467 ++++++++++++++++++++++++++ test/test_solvers.py | 3 +- test/test_windowing_integration.py | 277 +++++++++++++++ 18 files changed, 1963 insertions(+), 9 deletions(-) create mode 100644 doc/rolling-horizon.rst create mode 100644 linopy/restrict/__init__.py create mode 100644 linopy/restrict/api.py create mode 100644 linopy/restrict/errors.py create mode 100644 linopy/restrict/matrix.py create mode 100644 linopy/restrict/window.py create mode 100644 test/test_restrict_matrix.py create mode 100644 test/test_windowing_integration.py diff --git a/doc/api.rst b/doc/api.rst index 0a56c6d3b..693712bcb 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -74,6 +74,7 @@ Modifying a model model.Model.remove_objective model.Model.remove_sos_constraints model.Model.assign_coords + model.Model.restrict model.Model.copy model.Model.apply_sos_reformulation model.Model.undo_sos_reformulation diff --git a/doc/index.rst b/doc/index.rst index 06a1414f1..9581d8e4d 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -134,6 +134,7 @@ This package is published under MIT license. sos-constraints piecewise-linear-constraints + rolling-horizon numerical-scaling sparse-models testing-framework diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 55f780578..3776a7a2b 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -90,6 +90,7 @@ Version 0.10.1 * ``read_netcdf`` no longer drops parameters, parameter attributes or ``MultiIndex`` dimensions whose name contains ``-``. (`#1025 `__) * ``to_netcdf``/``read_netcdf`` keep the auxiliary (non-index) coordinates of dense variables, constraints, expressions and ``Model.parameters``, such as the levels on a flat dimension that v1 recommends instead of a ``MultiIndex``, and scalar coordinates left by a scalar selection. (`#1024 `__) * ``licensed_solvers`` and ``check_solver_licenses()`` no longer report Xpress as licensed when its licence is rejected, e.g. a licence that only covers an older Xpress release. Xpress now probes its licence with ``xpress.init()``. (`#994 `__) +* :meth:`linopy.Model.restrict` turns rolling-horizon optimization into a pure linopy pattern. ``m.restrict(snapshot=window, fix=...)`` returns a :class:`~linopy.restrict.WindowModel` built from a slice of the parent's matrices: constraint terms referencing variables outside the window are folded into the right-hand side at their actual matrix coefficient, ``rhs -= coeff * fixed_value``, and rows that lose every term are checked for feasibility and dropped. Because the window is a row slice and a column range of the parent's ``A``, ``b``, ``c``, ``lb`` and ``ub``, a step costs a slice rather than a rebuild of every variable and constraint. ``window.solve()`` hands those arrays to the solver and updates the existing solver model in place when consecutive windows share a shape, through the same :class:`~linopy.persistent.diff.ModelDiff` delta and per-backend apply hooks as the persistent interface. ``window.advance(snapshot=next, fix=...)`` moves the window. Solutions and duals are built per container only when read. See the new :doc:`rolling-horizon` guide. (https://github.com/PyPSA/linopy/issues/768) Version 0.10.0 -------------- diff --git a/doc/rolling-horizon.rst b/doc/rolling-horizon.rst new file mode 100644 index 000000000..3fc8b9d47 --- /dev/null +++ b/doc/rolling-horizon.rst @@ -0,0 +1,127 @@ +.. _rolling-horizon: + +Rolling-Horizon Optimization +============================ + +A long-horizon problem is often solved one window at a time: solve a slice of +the horizon (say, one day), carry the end state into the next window, then move +the window forward. Linopy has three methods for this: + +* :meth:`linopy.Model.restrict` returns a window of the model, +* :meth:`linopy.restrict.WindowModel.solve` solves it, +* :meth:`linopy.restrict.WindowModel.advance` moves it forward. + +A window is a row and column slice of the model's matrices. Moving it does not +rebuild variables or constraints. + +.. contents:: + :local: + :depth: 2 + +The rolling-horizon recipe +-------------------------- + +.. code-block:: python + + import pandas as pd + import xarray as xr + import linopy + + sns = pd.date_range("2026-01-01", periods=8760, freq="h", name="snapshot") + m = linopy.Model() + soc = m.add_variables(0, 100, coords=[sns], name="soc") + flow = m.add_variables(0, 1000, coords=[sns], name="flow") + # the balance references the previous snapshot + m.add_constraints(soc - 0.5 * soc.shift(snapshot=1) - flow <= 10, name="balance") + m.add_objective(flow.sum(), sense="min") + + # initial state of charge over the full horizon + initial = xr.DataArray([50.0] * 8760, dims=["snapshot"]) + + window = m.restrict(snapshot=sns[:24], fix={"soc": initial}) + window.solve(solver_name="highs") + + for start in range(24, 8760 - 24, 24): + window.advance(snapshot=sns[start : start + 24], fix={"soc": initial}) + window.solve(solver_name="highs") + + # containers are built when you read them + solution = window.solution + +:meth:`~linopy.Model.restrict` follows :meth:`xarray.Dataset.sel` semantics: +every container holding the selected dimension is sliced, containers without it +stay whole. The parent model is unchanged. + +Boundary values +--------------- + +The balance constraint references the previous snapshot's state of charge. In +the first row of a window, that reference points outside the window. +:meth:`~linopy.Model.restrict` moves those terms to the right-hand side, +multiplied by their coefficient: + +.. math:: + + A_{WW}\, x_W \sim b_W - A_{WB}\, x_B + +where :math:`x_W` are the window's variables, :math:`x_B` the variables +outside it, and :math:`x_B` on the right-hand side the values you pass in +``fix``. In words, each outside term adds ``coefficient * fixed value`` to the +right-hand side. Rows that lose every term are checked and dropped when they +are satisfied. + +Without ``fix``, a window that cuts a coupling raises +:class:`~linopy.restrict.RestrictionError` naming the constraints. A window +with no outside terms needs no ``fix``. + +.. note:: + + Folding uses the constraint's *actual* coefficient on the outside term. + Hand-written rolling-horizon code often uses a coefficient of 1 for the + initial state of charge, instead of the loss coefficient used on the other + rows. ``restrict`` reproduces a slightly different model. + +``fix`` values are keyed by variable name and broadcast over the variable's +coordinates, so a full-horizon DataArray works directly. Values for in-window +variables are ignored. A ``NaN`` at an outside label counts as missing. + +Working with a window +--------------------- + +:class:`~linopy.restrict.WindowModel` hands the solver the window's matrices +directly: + +* :attr:`~linopy.restrict.WindowModel.matrices` holds ``A``, ``b``, ``sense``, + ``c``, ``lb``, ``ub`` and the labels of its rows and columns. +* :attr:`~linopy.restrict.WindowModel.solution` and + :attr:`~linopy.restrict.WindowModel.dual` are ``Dataset`` objects with one + array per variable or constraint. They are built when you read them. +* :attr:`~linopy.restrict.WindowModel.objective` is the last solve's value. + +Terms outside the window are dropped from the objective, so the window's value +differs from the full model's by a constant. + +Moving the window +----------------- + +:meth:`~linopy.restrict.WindowModel.advance` slices the matrices for the new +selection and moves the new boundary values to the right-hand side. When the +new window has the same shape, the solver updates its existing model instead of +building a new one. Nothing is rebuilt, so a chain of equal-width windows +keeps that fast path. + +Limits +------ + +A window slices one dimension. Every container holding that dimension must +contain the window's labels; otherwise ``restrict`` raises. + +Several selected dimensions, indicators and quadratic objectives raise +:class:`~linopy.restrict.RestrictionError` or ``NotImplementedError``. + +Native SOS constraints also raise: they are not part of the matrices. Reformulate +them first with :meth:`Model.apply_sos_reformulation`, though the cardinality row +that adds spans the horizon and still cannot be windowed. + +Piecewise formulations work; they expand into ordinary variables and constraints +when added. diff --git a/linopy/examples.py b/linopy/examples.py index 549574bf6..2a0d8acf9 100644 --- a/linopy/examples.py +++ b/linopy/examples.py @@ -2,6 +2,8 @@ This module contains examples of linear programming models using the linopy library. """ +from typing import Any + import xarray as xr from numpy import arange @@ -68,6 +70,8 @@ def benchmark_model(n: int = 10, integerlabels: bool = False) -> Model: Model: The created linear programming model. """ m = Model() + naxis: Any + maxis: Any if integerlabels: naxis, maxis = [arange(n, dtype=int), arange(n, dtype=int)] else: diff --git a/linopy/matrices.py b/linopy/matrices.py index 7ecbb1d36..0156fb322 100644 --- a/linopy/matrices.py +++ b/linopy/matrices.py @@ -38,10 +38,16 @@ class _RowBlock(NamedTuple): row_scaling: ndarray | None -def _row_block(con: Constraint, label_index: VariableLabelIndex) -> _RowBlock: - active = con.active_row_mask() +def _row_scaling(con: Constraint, active: ndarray) -> ndarray: + """Scaling of the active rows of a constraint, in row order.""" + return con.scaling.values.ravel()[active] + + +def _row_block( + con: Constraint, label_index: VariableLabelIndex, active: ndarray +) -> _RowBlock: store = con.lhs._store - scaling = con.scaling.values.ravel()[active] + scaling = _row_scaling(con, active) return _RowBlock( partial(store.positional_csr, label_index, active, con.name), int(active.sum()), @@ -186,17 +192,24 @@ def _build_cons(self) -> None: regular, indicator = [], [] binvar, binval = [], [] + scalings = [] for c in m.constraints.data.values(): + active = c.active_row_mask() if not c.is_indicator: - regular.append(_row_block(c, label_index)) + regular.append(_row_block(c, label_index, active)) + scalings.append(_row_scaling(c, active)) continue - active = c.active_row_mask() - indicator.append(_row_block(c, label_index)) + indicator.append(_row_block(c, label_index, active)) var_labels, val = _indicator_rows(c, active) binvar.append(label_to_pos[var_labels]) binval.append(val) self.clabels: ndarray = m.constraints.label_index.clabels + #: Row scaling actually applied to ``A``, aligned with ``clabels``. Consumers + #: that report duals need it to invert the scaling on the way back. + self.con_scaling: ndarray = ( + np.concatenate(scalings) if scalings else np.empty(0, dtype=float) + ) self.A: scipy.sparse.csr_array | None self.A, self.b, self.sense = _stack(regular, label_index, col_scaling, m) self.indicator_A: scipy.sparse.csr_array | None diff --git a/linopy/model.py b/linopy/model.py index 9c1b59a6d..6479bc660 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -117,6 +117,7 @@ if TYPE_CHECKING: from linopy.piecewise import PiecewiseFormulation + from linopy.restrict.window import WindowModel logger = logging.getLogger(__name__) @@ -2295,6 +2296,69 @@ def _coordinate_carriers(self) -> list[Any]: *self.expressions.data.values(), ] + def restrict( + self, + fix: Mapping[str, Any] | None = None, + **selection: Any, + ) -> WindowModel: + """ + Build a matrix-backed window along one coordinate dimension. + + A window is a slice of the model's flattened matrices: the parent is + built once, and each selection slices its rows and columns. Terms whose + variables fall outside the window are folded into the right-hand side at + their actual coefficient when ``fix`` supplies their values. Rows that + lose every term are dropped when their folded right-hand side satisfies + their sign, and raise otherwise. + + The returned :class:`~linopy.restrict.WindowModel` exposes ``matrices`` + directly, so solving never rebuilds a container. Variables and + constraints are materialised only when their values are read. + + Parameters + ---------- + fix : mapping, optional + Fixed values for out-of-window variables, keyed by variable name and + broadcast against the variable's coordinates. Values for in-window + variables are ignored. + **selection : Any + Exactly one dimension and the labels of the window, e.g. + ``snapshot=sns[:24]``. + + Returns + ------- + WindowModel + The window, with its matrices computed. + + Raises + ------ + RestrictionError + If the dimension does not exist, a carrier does not contain the + window's labels, or a fully folded row is unsatisfiable. + NotImplementedError + If the selection names several dimensions, or the model has + indicators, a quadratic objective, SOS or piecewise formulations. + + Examples + -------- + >>> import pandas as pd + >>> import linopy + >>> + >>> sns = pd.date_range("2026-01-01", periods=4, freq="h", name="snapshot") + >>> m = linopy.Model() + >>> storage = m.add_variables(0, 10, coords=[sns], name="storage") + >>> flow = m.add_variables(0, 100, coords=[sns], name="flow") + >>> _ = m.add_constraints(storage - flow <= 0, name="cap") + >>> _ = m.add_objective(flow.sum(), sense="min") + >>> + >>> window = m.restrict(snapshot=sns[:2]) + >>> window.matrices.A.shape + (2, 4) + """ + from linopy.restrict import restrict_model + + return restrict_model(self, fix, selection) + def assign_result( self, result: Result, diff --git a/linopy/persistent/diff.py b/linopy/persistent/diff.py index 14d7c5ed3..e3944a24b 100644 --- a/linopy/persistent/diff.py +++ b/linopy/persistent/diff.py @@ -4,7 +4,7 @@ from collections.abc import Iterable from dataclasses import dataclass from functools import cached_property -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Any import numpy as np @@ -14,6 +14,7 @@ ContainerVarBuffers, ModelSnapshot, StructuralKey, + VarKind, _coord_snapshot, _extract_con_buffers, _extract_var_buffers, @@ -29,6 +30,15 @@ from linopy.variables import Variable +#: Matrix column kinds (``vtypes``) mapped onto the persistent diff's enum. +_VAR_KIND = { + "C": VarKind.CONTINUOUS, + "B": VarKind.BINARY, + "I": VarKind.INTEGER, + "S": VarKind.SEMI_CONTINUOUS, +} + + class RebuildReason(enum.Enum): STRUCTURAL_LABELS = "vlabels/clabels mismatch" STRUCTURAL_CONTAINERS = "container set changed" @@ -143,7 +153,7 @@ class ModelDiff: #: buffers the diff walk already extracted — adopting it after a #: successful apply replaces a full re-capture. Note: holding a diff #: therefore pins all container buffers for its lifetime. - snapshot: ModelSnapshot + snapshot: ModelSnapshot | None @property def is_empty(self) -> bool: @@ -347,6 +357,81 @@ def from_models( ignore_dims=ignore_dims, ) + @classmethod + def from_windows(cls, previous: Any, current: Any) -> ModelDiff: + """ + Delta between two matrix windows of the same shape. + + A rolling-horizon window is already flat, so the difference is a set of + changed bound, rhs, cost and coefficient entries. Building it here lets + the window path reuse :meth:`Solver.apply_update` and every backend's + in-place hooks, instead of a parallel update implementation. + """ + columns = np.flatnonzero( + (current.lb != previous.lb) | (current.ub != previous.ub) + ) + rows = np.flatnonzero(current.b != previous.b) + cost = np.flatnonzero(current.c != previous.c) + + coef_deltas: list[_CoefDelta] = [] + n_coef_updates = 0 + changed = np.flatnonzero(current.A.data != previous.A.data) + if changed.size: + changed_rows = np.unique( + np.searchsorted(current.A.indptr, changed, side="right") - 1 + ) + counts = np.diff(current.A.indptr)[changed_rows] + buffer = ContainerConBuffers( + indptr=current.A.indptr, + indices=current.A.indices, + data=current.A.data, + rhs=current.b, + sign=current.sense, + active_labels=np.empty(0, dtype=np.int64), + ) + coef_deltas = [ + _CoefDelta(buffer, changed_rows, changed_rows, int(counts.sum())) + ] + n_coef_updates = int(counts.sum()) + + # Windows always carry a per-row sign array, so no scalar handling here. + rhs_signs = current.sense[rows] + sign_changed = np.flatnonzero( + np.broadcast_to(current.sense, current.b.shape) + != np.broadcast_to(previous.sense, previous.b.shape) + ) + sign_values = current.sense[sign_changed] + + types = np.flatnonzero(current.vtypes != previous.vtypes) + kinds = np.array( + [_VAR_KIND[kind] for kind in current.vtypes[types]], dtype=object + ) + + return cls( + var_bounds_indices=columns, + var_bounds_lower=current.lb[columns], + var_bounds_upper=current.ub[columns], + var_type_positions=types, + var_type_kinds=kinds, + coef_deltas=coef_deltas, + n_coef_updates=n_coef_updates, + con_rhs_indices=rows, + con_rhs_values=current.b[rows], + con_rhs_signs=rhs_signs, + con_sign_indices=sign_changed, + con_sign_values=sign_values, + obj_c_indices=cost if cost.size else None, + obj_c_values=current.c[cost] if cost.size else None, + obj_sense=( + current.objective_sense + if current.objective_sense != previous.objective_sense + else None + ), + var_slices={}, + con_slices={}, + snapshot=None, + ) + class _DiffBuilder: """Accumulates per-container deltas and finalizes them into a ModelDiff.""" diff --git a/linopy/restrict/__init__.py b/linopy/restrict/__init__.py new file mode 100644 index 000000000..788163555 --- /dev/null +++ b/linopy/restrict/__init__.py @@ -0,0 +1,40 @@ +""" +Rolling-horizon primitives, backed by a model's flattened matrices. + +:meth:`linopy.Model.restrict` returns a :class:`~linopy.restrict.WindowModel`, +a window of the parent model expressed as a row slice and column range of its +flattened ``A``, ``b``, ``c``, ``lb`` and ``ub``. Terms whose variables fall +outside the window are folded into the right-hand side at their actual matrix +coefficient: + + ``A_WW @ x_W ~ b_W - A_WB @ x_bar_B`` + +with ``x_bar_B`` the fixed boundary values, i.e. the subproblem-conditioning +primitive of decomposition methods. + +Working on the matrices means a window costs a slice, not a rebuild of every +variable and constraint. The window materialises a variable or constraint only +when its values are read, and solves by handing the flat arrays to the solver. + +The package is organized around: + +* :mod:`linopy.restrict.matrix` - :class:`MatrixWindow`, the flattened view. +* :mod:`linopy.restrict.window` - :class:`WindowModel`, the user-facing window. +* :mod:`linopy.restrict.api` - the :func:`restrict_model` entry point. +""" + +from __future__ import annotations + +from linopy.restrict.api import restrict_model +from linopy.restrict.errors import RestrictionError +from linopy.restrict.matrix import MatrixWindow, WindowMatrices +from linopy.restrict.window import WindowModel, WindowObjective + +__all__ = [ + "MatrixWindow", + "RestrictionError", + "WindowMatrices", + "WindowModel", + "WindowObjective", + "restrict_model", +] diff --git a/linopy/restrict/api.py b/linopy/restrict/api.py new file mode 100644 index 000000000..d9074cd4e --- /dev/null +++ b/linopy/restrict/api.py @@ -0,0 +1,48 @@ +""" +Rolling-horizon entry point, backed by the flattened matrices. + +:meth:`Model.restrict` delegates here. A window is a slice of the parent's +flattened matrices (:class:`linopy.restrict.MatrixWindow`), returned as a +:class:`~linopy.restrict.window.WindowModel` that materialises variables and +constraints only when their values are read. +""" + +from __future__ import annotations + +from collections.abc import Mapping +from typing import TYPE_CHECKING, Any + +from linopy.restrict.errors import RestrictionError + +if TYPE_CHECKING: + from linopy.model import Model + from linopy.restrict.window import WindowModel + +__all__ = ["RestrictionError", "restrict_model"] + + +def restrict_model( + model: Model, + fix: Mapping[str, Any] | None, + selection: Mapping[str, Any], +) -> WindowModel: + """ + Build a matrix-backed window of ``model``. + + See :meth:`linopy.Model.restrict` for the documented behavior. + """ + from linopy.restrict.window import WindowModel + + if len(selection) != 1: + raise NotImplementedError( + "Model.restrict supports exactly one selected dimension, got " + f"{sorted(selection)}." + ) + dim = next(iter(selection)) + carried = any(dim in item.sizes for item in model._coordinate_carriers()) + if not carried and dim not in model.parameters.sizes: + raise RestrictionError( + f"Cannot restrict the model along '{dim}': the dimension does not " + "exist in the model." + ) + return WindowModel(model, str(dim), dict(selection), dict(fix) if fix else None) diff --git a/linopy/restrict/errors.py b/linopy/restrict/errors.py new file mode 100644 index 000000000..fc109c11f --- /dev/null +++ b/linopy/restrict/errors.py @@ -0,0 +1,13 @@ +"""Errors raised by the rolling-horizon primitives.""" + +from __future__ import annotations + + +class RestrictionError(ValueError): + """ + Raised when a model cannot be windowed. + + A window requires a dimension that exists in the model, carriers that + contain every selected label, and folded boundary terms whose constant rows + are satisfiable. + """ diff --git a/linopy/restrict/matrix.py b/linopy/restrict/matrix.py new file mode 100644 index 000000000..df39e759c --- /dev/null +++ b/linopy/restrict/matrix.py @@ -0,0 +1,364 @@ +""" +Window a model's flattened matrices. + +A rolling-horizon window is a contiguous run of one coordinate dimension, which +on the flattened matrices is a row slice of the parent CSR plus a column range. +:class:`MatrixWindow` is built once from a parent model, and +:meth:`MatrixWindow.select` returns the window's ``A``, ``b``, ``sense``, ``c``, +``lb`` and ``ub`` together with the parent labels of its rows and columns, +without rebuilding a single container. + +Restricting and advancing are the same operation here: a selection. There is no +per-container term pattern to freeze, so a window never rebuilds on a +data-dependent coefficient change. + +Limitations, all raised rather than silently approximated: + +* one selected dimension, +* every carrier of that dimension must share its index, +* no indicator constraints, quadratic objective, SOS or piecewise formulations. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import numpy as np +import pandas as pd +import scipy.sparse + +from linopy.restrict.errors import RestrictionError + +if TYPE_CHECKING: + from linopy.model import Model + +__all__ = ["MatrixWindow", "WindowMatrices"] + + +@dataclass(frozen=True) +class WindowMatrices: + """ + One window of a parent model, as flat arrays. + + Columns follow the parent variable labels in ``vlabels``; rows follow the + parent constraint labels in ``clabels``. Out-of-window objective terms are + dropped, matching :meth:`Model.restrict`, so the objective value differs + from the parent's by a constant. + """ + + A: scipy.sparse.csr_array + b: np.ndarray + sense: np.ndarray + c: np.ndarray + lb: np.ndarray + ub: np.ndarray + vlabels: np.ndarray + clabels: np.ndarray + objective_sense: str + #: Per-container ``(start, stop)`` in the window's flat column and row space, + #: so a solution vector can be handed back container by container. + var_blocks: dict[str, tuple[int, int]] + con_blocks: dict[str, tuple[int, int]] + #: Coordinate values of the selected dimension in this window. + coords: dict[str, pd.Index] + #: Per-container active masks, in the window's grid order, so a flat result + #: vector can be scattered back into each container's shape. + var_active: dict[str, np.ndarray] + con_active: dict[str, np.ndarray] + #: Variable kinds (``C``, ``B``, ``I``, ``S``) and the scaling applied to + #: columns and rows. Results must be unscaled with them before reporting. + vtypes: np.ndarray + var_scaling: np.ndarray + con_scaling: np.ndarray + objective_scaling: float + + @property + def Q(self) -> None: + """No quadratic objective: matrix windows are linear.""" + return None + + @property + def indicator_A(self) -> None: + """No indicator constraints: matrix windows reject them.""" + return None + + +@dataclass(frozen=True) +class WindowResult: + """Flat result of solving one window.""" + + objective: float + primal: np.ndarray + dual: np.ndarray + + +@dataclass(frozen=True) +class _Layout: + """One container's opening position in the flat space, and its active mask.""" + + start: int + active: np.ndarray + index: pd.Index | None + labels: Any + dim: str + rank: np.ndarray + + @property + def size(self) -> int: + return int(self.active.sum()) + + def select(self, window: pd.Index) -> tuple[np.ndarray, np.ndarray]: + """ + Flat positions of the window's entries, and the window's active mask. + + The mask is laid out in the window's own grid order, so a flat result + vector can be scattered back into the container's shape. + """ + if not self.size: + return np.empty(0, dtype=np.int64), np.zeros(0, dtype=bool) + if self.index is None: + # the container does not carry the selected dimension: kept whole + return self.start + np.arange(self.size), np.asarray(self.active) + lookup = self.index.get_indexer(window) + if (lookup < 0).any(): + missing = window[lookup < 0] + raise RestrictionError( + f"Cannot window the matrices: labels {missing.tolist()} are not " + "carried by every container of the selected dimension." + ) + # Membership from the positional lookup rather than `xarray.isin`, which + # compares tz-aware timestamps as objects and is slow. + present = np.zeros(self.index.size, dtype=bool) + present[lookup] = True + shape = [1] * self.labels.ndim + shape[self.labels.dims.index(self.dim)] = present.size + flat = np.broadcast_to(present.reshape(shape), self.labels.shape).ravel() + window_entries = np.flatnonzero(flat) + window_active = np.asarray(self.active)[window_entries] + return self.start + self.rank[window_entries[window_active]], window_active + + +def _layouts(containers: Any, dim: str) -> dict[str, _Layout]: + """The flat layout of every container along ``dim``.""" + layouts: dict[str, _Layout] = {} + position = 0 + for name, item in containers.items(): + labels = item.labels + active = np.asarray(labels.values).ravel() != -1 + layouts[name] = _Layout( + position, active, item.indexes.get(dim), labels, dim, np.cumsum(active) - 1 + ) + position += int(active.sum()) + return layouts + + +def _blocks( + positions: dict[str, np.ndarray], keep: np.ndarray | None = None +) -> dict[str, tuple[int, int]]: + """ + Per-container ``(start, stop)`` in the concatenated flat space. + + ``keep`` masks the concatenated rows; the returned blocks index the rows that + survive it, which is the order of the window arrays. + """ + blocks: dict[str, tuple[int, int]] = {} + cursor = 0 + size = 0 + for name, values in positions.items(): + count = int(values.size) + selected = count if keep is None else int(keep[cursor : cursor + count].sum()) + blocks[name] = (size, size + selected) + size += selected + cursor += count + return blocks + + +def _row_owners( + rows: np.ndarray, selection: dict[str, tuple[np.ndarray, np.ndarray]] +) -> list[str]: + """Names of the containers owning the given concatenated row positions.""" + owners = [] + cursor = 0 + for name, (positions, _) in selection.items(): + count = int(positions.size) + if np.any((rows >= cursor) & (rows < cursor + count)): + owners.append(name) + cursor += count + return owners + + +def _filter_masks( + selection: dict[str, tuple[np.ndarray, np.ndarray]], keep: np.ndarray +) -> dict[str, np.ndarray]: + """Per-container active masks, with the dropped rows removed.""" + masks: dict[str, np.ndarray] = {} + cursor = 0 + for name, (positions, mask) in selection.items(): + count = int(positions.size) + # `keep` is indexed over active rows, one per entry of `positions`; the + # mask covers the whole window grid, so only its active slots move. + filtered = mask.copy() + filtered[np.flatnonzero(mask)] = keep[cursor : cursor + count] + masks[name] = filtered + cursor += count + return masks + + +class MatrixWindow: + """ + A matrix-level view of a model, windowed along one dimension. + + Build once, then call :meth:`select` per window. The parent model is read + only; every selection allocates its own arrays. + """ + + def __init__(self, model: Model, dim: str) -> None: + if model.variables.sos: + # Piecewise formulations and reformulated SOS are ordinary variables and + # constraints, so they window like anything else. Native SOS sets are not + # in the matrices, so they must be reformulated first. + raise NotImplementedError( + "Matrix windows do not support unreformulated SOS constraints: " + "call `model.apply_sos_reformulation()` first." + ) + matrices = model.matrices + if matrices.Q is not None: + raise NotImplementedError( + "Matrix windows do not support quadratic objectives." + ) + if matrices.indicator_A is not None and matrices.indicator_A.nnz: + raise NotImplementedError( + "Matrix windows do not support indicator constraints." + ) + if matrices.A is None: + raise NotImplementedError("Matrix windows require at least one constraint.") + self._model = model + self._matrices = matrices + self._A = matrices.A.tocsr() + self._dim = dim + self._var_layout = _layouts(model.variables, dim) + self._con_layout = _layouts(model.constraints, dim) + + def select( + self, window: pd.Index, fix: dict[str, Any] | None = None + ) -> WindowMatrices: + """ + The window's matrices. + + ``fix`` maps variable names to values carrying dangling boundary terms, + exactly like :meth:`Model.restrict`. Out-of-window terms are folded into + the right-hand side at their actual coefficient; a row that loses every + term is dropped when its folded right-hand side satisfies its sign. + """ + con_selection = { + name: layout.select(window) for name, layout in self._con_layout.items() + } + var_selection = { + name: layout.select(window) for name, layout in self._var_layout.items() + } + rows = np.concatenate([positions for positions, _ in con_selection.values()]) + cols = np.concatenate([positions for positions, _ in var_selection.values()]) + fix_by_position = self._fix_by_position(fix) + + A = self._A[rows] + n_rows = A.shape[0] + in_window = np.zeros(A.shape[1], dtype=bool) + in_window[cols] = True + + row_of_entry = np.repeat(np.arange(n_rows), np.diff(A.indptr)) + dangling = ~in_window[A.indices] & (A.data != 0) + unfixed = dangling & np.isnan(fix_by_position[A.indices]) + if unfixed.any(): + owners = _row_owners(np.unique(row_of_entry[unfixed]), con_selection) + raise RestrictionError( + f"Cannot restrict model: constraint(s) {owners} have dangling " + "terms on variables outside the selection without fixed values. " + "Pass their values via `fix`." + ) + contributions = A.data[dangling] * fix_by_position[A.indices[dangling]] + rhs = self._matrices.b[rows] - np.bincount( + row_of_entry[dangling], weights=contributions, minlength=n_rows + ) + + kept_entries = in_window[A.indices] + keep = np.bincount(row_of_entry[kept_entries], minlength=n_rows) > 0 + self._check_dropped_rows(rhs, row_of_entry[dangling], contributions, keep) + + sense = self._matrices.sense[rows] + return WindowMatrices( + A=A[keep][:, cols], + b=rhs[keep], + sense=sense[keep], + c=self._matrices.c[cols], + lb=self._matrices.lb[cols], + ub=self._matrices.ub[cols], + vlabels=self._matrices.vlabels[cols], + clabels=self._matrices.clabels[rows][keep], + objective_sense=self._model.objective.sense, + var_blocks=_blocks({name: p for name, (p, _) in var_selection.items()}), + con_blocks=_blocks( + {name: p for name, (p, _) in con_selection.items()}, keep + ), + coords={self._dim: window}, + var_active={name: mask for name, (_, mask) in var_selection.items()}, + con_active=_filter_masks(con_selection, keep), + vtypes=self._matrices.vtypes[cols], + var_scaling=self._matrices.var_scaling[cols], + con_scaling=self._matrices.con_scaling[rows][keep], + objective_scaling=float(self._model.objective.scaling), + ) + + def _fix_by_position(self, fix: dict[str, Any] | None) -> np.ndarray: + """Fixed values on the flat parent column space, NaN where unfixed.""" + n_cols = len(self._matrices.vlabels) + values = np.full(n_cols, np.nan) + if not fix: + return values + for name, given in fix.items(): + if name not in self._var_layout: + raise RestrictionError( + f"Cannot fix {name!r}: it does not name a variable of the model." + ) + layout = self._var_layout[name] + flat = np.asarray(given).ravel() + if flat.size != layout.active.size: + raise RestrictionError( + f"Fixed values for variable '{name}' do not match its " + f"dimensions: got {flat.size} values, expected {layout.active.size}." + ) + values[layout.start : layout.start + layout.size] = flat[layout.active] + return values + + def _check_dropped_rows( + self, + rhs: np.ndarray, + dangling_rows: np.ndarray, + contributions: np.ndarray, + keep: np.ndarray, + ) -> None: + """Raise when a row that lost every term violates its sign.""" + dropped = np.flatnonzero(~keep) + if not dropped.size: + return + magnitude = np.bincount( + dangling_rows, weights=np.abs(contributions), minlength=len(keep) + ) + tolerance = 1e-9 * np.maximum(1.0, magnitude[dropped]) + residual = rhs[dropped] + sense = self._matrices.sense + signs = ( + np.full(dropped.size, sense) if isinstance(sense, str) else sense[dropped] + ) + violated = ( + ((signs == "=") & (np.abs(residual) > tolerance)) + | ((signs == "<=") & (residual < -tolerance)) + | ((signs == ">=") & (residual > tolerance)) + ) + if violated.any(): + raise RestrictionError( + "Cannot window the matrices: folding the dangling terms makes " + f"rows {dropped[violated].tolist()} unsatisfiable " + f"(folded rhs {residual[violated].tolist()} against sign " + f"{signs[violated].tolist()})." + ) diff --git a/linopy/restrict/window.py b/linopy/restrict/window.py new file mode 100644 index 000000000..bac49c197 --- /dev/null +++ b/linopy/restrict/window.py @@ -0,0 +1,375 @@ +""" +A matrix-backed window of a parent model. + +:class:`WindowModel` is what :meth:`Model.restrict` returns. It exposes the +window's flattened matrices directly, so solving never rebuilds a container, and +materialises variables and constraints only when their values are read. +""" + +from __future__ import annotations + +from collections.abc import Hashable, Iterator +from typing import TYPE_CHECKING, Any + +import numpy as np +import pandas as pd +from xarray import DataArray, Dataset + +from linopy.restrict.errors import RestrictionError +from linopy.restrict.matrix import MatrixWindow, WindowMatrices + +if TYPE_CHECKING: + from linopy.model import Model + from linopy.solvers import Solver + +__all__ = ["WindowModel", "WindowObjective"] + + +class _LabelIndex: + """ + Dense positional labels. + + A window's result vectors are positional, so labels are ``0..n-1``. That + keeps the solver's label-indexed solution mapping aligned with the flat + arrays instead of allocating an array as large as the parent's label space. + """ + + def __init__(self, size: int) -> None: + labels = np.arange(size) + self.vlabels = labels + self.clabels = labels + self.label_to_pos = labels + self.n_active_vars = size + self.n_active_cons = size + + +class _KindSet: + """Variables of one kind, with the ``labels`` attribute builders read.""" + + def __init__(self, labels: np.ndarray) -> None: + self.labels = labels + + def __len__(self) -> int: + return int(self.labels.size) + + def __bool__(self) -> bool: + return bool(self.labels.size) + + def __iter__(self) -> Iterator[Any]: + return iter(self.labels) + + +class _VariablesView: + """The ``model.variables`` surface a solver's direct builder reads.""" + + def __init__(self, view: WindowModelView) -> None: + self._view = view + + @property + def sos(self) -> dict: + return {} + + @property + def semi_continuous(self) -> _KindSet: + return self._view.semi_continuous + + @property + def label_index(self) -> _LabelIndex: + return _LabelIndex(len(self._view.matrices.vlabels)) + + +class _ConstraintsView: + """The ``model.constraints`` surface a builder reads.""" + + def __init__(self, window: WindowModel) -> None: + self._window = window + + def sanitize_missings(self) -> None: + return None + + def sanitize_zeros(self) -> None: + return None + + def sanitize_infinities(self) -> None: + return None + + def __len__(self) -> int: + return len(self._window.matrices.con_blocks) + + @property + def label_index(self) -> _LabelIndex: + return _LabelIndex(len(self._window.matrices.clabels)) + + +class WindowModelView: + """ + The slice of :class:`~linopy.Model` that a solver's direct builder reads. + + Passing this to ``Solver._build`` lets every solver's existing build path + consume a window unchanged, instead of each backend growing its own + build-from-arrays routine. + """ + + def __init__(self, window: WindowModel) -> None: + self._window = window + self.variables = _VariablesView(self) + self.constraints = _ConstraintsView(window) + + @property + def matrices(self) -> WindowMatrices: + return self._window.matrices + + @property + def objective(self) -> WindowObjective: + return self._window.objective + + @property + def sense(self) -> str: + return self.matrices.objective_sense + + @property + def set_names_in_solver_io(self) -> bool: + return False + + @property + def is_quadratic(self) -> bool: + return False + + @property + def nvars(self) -> int: + return len(self.matrices.vlabels) + + @property + def ncons(self) -> int: + return len(self.matrices.clabels) + + @property + def _xCounter(self) -> int: + return len(self.matrices.vlabels) + + @property + def _cCounter(self) -> int: + return len(self.matrices.clabels) + + @property + def indicator_constraints(self) -> tuple: + return () + + def _of_kind(self, kind: str) -> _KindSet: + return _KindSet(np.flatnonzero(self.matrices.vtypes == kind)) + + @property + def binaries(self) -> _KindSet: + return self._of_kind("B") + + @property + def integers(self) -> _KindSet: + return self._of_kind("I") + + @property + def semi_continuous(self) -> _KindSet: + return self._of_kind("S") + + @property + def continuous(self) -> _KindSet: + return self._of_kind("C") + + @property + def type(self) -> str: + discrete = bool(self.binaries or self.integers or self.semi_continuous) + if not discrete: + return "LP" + return "MILP" if self.continuous else "ILP" + + +class WindowObjective: + """The window's objective value, as read off the last solve.""" + + def __init__(self, window: WindowModel) -> None: + self._window = window + + @property + def value(self) -> float | None: + return self._window._objective_value + + @property + def sense(self) -> str: + return self._window.matrices.objective_sense + + +class WindowModel: + """ + A window of a parent model, backed by its flattened matrices. + + ``matrices`` is available immediately. ``solution``, ``dual`` and access to a + named variable or constraint materialise that container from the flat result + vectors, so a caller pays only for what it reads. + """ + + def __init__( + self, + parent: Model, + dim: str, + selection: dict[str, Any], + fix: dict[str, Any] | None = None, + ) -> None: + self._parent = parent + self._dim = dim + self._view = MatrixWindow(parent, dim) + self._fix = fix + self._selection = dict(selection) + self._matrices = self._view.select(pd.Index(selection[dim]), fix) + self._objective = WindowObjective(self) + self._objective_value: float | None = None + self._primal: np.ndarray | None = None + self._dual: np.ndarray | None = None + self._solver: Solver | None = None + self._previous: WindowMatrices | None = None + + @property + def objective(self) -> WindowObjective: + return self._objective + + @property + def matrices(self) -> WindowMatrices: + return self._matrices + + def as_model_view(self) -> WindowModelView: + """The window as the slice of ``Model`` a solver's builder consumes.""" + return WindowModelView(self) + + def advance( + self, fix: dict[str, Any] | None = None, **selection: Any + ) -> WindowModel: + """Move the window, recomputing the matrices for the new selection.""" + if selection: + missing = [dim for dim in self._selection if dim not in selection] + if missing: + raise RestrictionError( + "Cannot advance model: the selection must name all dimensions " + f"of the previous selection, missing {missing}." + ) + self._selection = {dim: selection[dim] for dim in self._selection} + self._fix = fix if fix is not None else self._fix + self._previous = self._matrices + self._matrices = self._view.select( + pd.Index(self._selection[self._dim]), self._fix + ) + self._objective_value = None + self._primal = None + self._dual = None + return self + + def solve( + self, + solver_name: str | None = None, + disallow_rebuild: bool = False, + **options: Any, + ) -> WindowModel: + """Solve the window, updating an existing solver model in place when possible.""" + from linopy.solvers import Solver + + if self._solver is None: + self._solver = Solver.from_name( + solver_name or "highs", io_api="direct", options=options + ) + elif options: + self._solver.solver_options.update(options) + result = self._solver.solve_window( + self, self._previous, disallow_rebuild=disallow_rebuild + ) + self._objective_value = result.objective + self._primal = result.primal + self._dual = result.dual + self._previous = self._matrices + return self + + # --- materialisation ----------------------------------------------------- + + @property + def solution(self) -> Dataset: + """The window's variable values, one array per variable.""" + self._require_solve() + assert self._primal is not None + data = {} + for name, var in self._parent.variables.items(): + start, stop = self._matrices.var_blocks[name] + # the solver works in scaled columns; report the model's units + values = self._primal[start:stop] / self._matrices.var_scaling[start:stop] + data[name] = self._reshape(var, values, self._matrices.var_active[name]) + return Dataset(data) + + @property + def dual(self) -> Dataset: + """The window's constraint duals, one array per constraint.""" + self._require_solve() + expected = sum( + stop - start for start, stop in self._matrices.con_blocks.values() + ) + if self._dual is None or self._dual.size != expected: + raise AttributeError("The last solve returned no dual values.") + data = {} + for name, con in self._parent.constraints.items(): + start, stop = self._matrices.con_blocks[name] + # undo the row and objective scaling the solver applied + values = ( + self._dual[start:stop] + * self._matrices.con_scaling[start:stop] + / self._matrices.objective_scaling + ) + data[name] = self._reshape( + con, + values, + self._matrices.con_active[name], + dims=tuple(con.coord_dims), + ) + return Dataset(data) + + def __getitem__(self, name: str) -> Any: + """A variable of the window, materialised on demand.""" + if name in self._parent.variables: + var = self._parent.variables[name] + start, stop = self._matrices.var_blocks[name] + values = None + if self._primal is not None: + values = ( + self._primal[start:stop] / self._matrices.var_scaling[start:stop] + ) + return self._reshape(var, values, self._matrices.var_active[name]) + raise KeyError(name) + + def _reshape( + self, + container: Any, + values: np.ndarray | None, + active: np.ndarray, + dims: tuple[Hashable, ...] | None = None, + ) -> DataArray: + """Values of one container, scattered to its shape and labelled for the window.""" + dims = tuple(container.dims) if dims is None else dims + sizes = {dim: container.shape[dims.index(dim)] for dim in dims} + coords = {} + for dim in dims: + if dim != self._dim: + coords[dim] = container.indexes[dim] + continue + # a carrier may hold a subset of the window's labels, so size and + # coordinates come from the intersection + index = container.indexes.get(dim) + window = self._matrices.coords[self._dim] + coords[dim] = window if index is None else index[index.isin(window)] + sizes[dim] = coords[dim].size + full = np.full(active.size, np.nan) + if values is not None: + full[active] = values + return DataArray( + full.reshape(tuple(sizes[dim] for dim in dims)), + coords=coords, + dims=list(dims), + ) + + def _require_solve(self) -> None: + if self._primal is None: + raise AttributeError( + "The window has not been solved. Call `solve()` before reading values." + ) diff --git a/linopy/solvers.py b/linopy/solvers.py index 1361b49e5..dbf3f2027 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -470,6 +470,73 @@ class Solver(ABC, Generic[EnvType]): features: ClassVar[frozenset[SolverFeature]] = frozenset() accepted_io_apis: ClassVar[frozenset[str]] = frozenset() supports_persistent_update: ClassVar[bool] = False + + def solve_window( + self, + window: Any, + previous: Any = None, + disallow_rebuild: bool = False, + **run_kwargs: Any, + ) -> Any: + """ + Solve one matrix window, reusing this solver's model when possible. + + ``previous`` is the window solved last, when the caller can guarantee the + two have the same shape; the solver then applies the difference in place. + The window reaches the backend through + :class:`~linopy.restrict.window.WindowModelView`, so each solver's existing + direct build path works unchanged. + """ + if self.io_api != "direct": + raise NotImplementedError( + f"{self.display_name} window solving requires io_api='direct', " + f"got {self.io_api!r}." + ) + from linopy.persistent.diff import ModelDiff + from linopy.restrict.matrix import WindowResult + + current = window.matrices + view = window.as_model_view() + update = ( + self.supports_persistent_update + and previous is not None + and previous.A.shape == current.A.shape + and getattr(self, "_window_view", None) is not None + ) + if update: + self.model = view + self.apply_update( + ModelDiff.from_windows(previous, current), + view.variables.label_index, + view.constraints.label_index, + ) + else: + # no in-place update, or a different shape: build a fresh native model + if getattr(self, "_window_view", None) is not None: + reason = ( + RebuildReason.STRUCTURAL_LABELS + if previous is not None and previous.A.shape != current.A.shape + else RebuildReason.BACKEND_REJECTED + ) + if disallow_rebuild: + raise RebuildRequiredError( + reason, "the window requires a solver rebuild" + ) + logger.warning("Window requires a solver rebuild (%s).", reason.value) + self.close() + self.model = view + self._build() + self._window_view = view + result = self._run_direct(**run_kwargs) + solution = result.solution + if solution is None: + raise RuntimeError("the solver returned no solution for the window") + return WindowResult( + objective=float(solution.objective), + primal=np.asarray(solution.primal), + dual=np.asarray(solution.dual), + ) + supports_sign_update: ClassVar[bool] = False def __post_init__(self) -> None: @@ -497,6 +564,9 @@ def apply_update( Template method: validates the diff up front (a rejected update leaves the native model untouched), then walks the sections in a fixed order, dispatching to the per-backend ``_apply_*`` hooks. + + Pass ``ctx`` to apply the diff to a native model handled outside the + persistent path, as the matrix window does. """ if not self.supports_persistent_update: raise UnsupportedUpdate(type(self).__name__) diff --git a/test/test_persistent_snapshot_diff.py b/test/test_persistent_snapshot_diff.py index 2d84b3d7f..e76b3e4e4 100644 --- a/test/test_persistent_snapshot_diff.py +++ b/test/test_persistent_snapshot_diff.py @@ -334,6 +334,7 @@ def test_diff_snapshot_matches_capture(baseline: Model, mutate: str) -> None: baseline.add_objective(3 * x.sum() + 2 * y.sum(), overwrite=True) diff = ModelDiff.from_snapshot(snap, baseline) assert isinstance(diff, ModelDiff) + assert diff.snapshot is not None _assert_snapshot_equal(diff.snapshot, ModelSnapshot.capture(baseline)) @@ -350,6 +351,7 @@ def build(t0: int) -> Model: snap = ModelSnapshot.capture(m1) diff = ModelDiff.from_snapshot(snap, m2, ignore_dims={"t"}) assert isinstance(diff, ModelDiff) + assert diff.snapshot is not None _assert_snapshot_equal(diff.snapshot, ModelSnapshot.capture(m2)) @@ -364,6 +366,7 @@ def build(rhs: float) -> Model: m1, m2 = build(4.0), build(7.0) diff = ModelDiff.from_models(m1, m2) assert isinstance(diff, ModelDiff) + assert diff.snapshot is not None _assert_snapshot_equal(diff.snapshot, ModelSnapshot.capture(m2)) diff --git a/test/test_restrict_matrix.py b/test/test_restrict_matrix.py new file mode 100644 index 000000000..0c009f0be --- /dev/null +++ b/test/test_restrict_matrix.py @@ -0,0 +1,467 @@ +"""Matrix-level windows agree with the container-level restriction.""" + +from __future__ import annotations + +from dataclasses import replace + +import numpy as np +import pandas as pd +import pytest + +from linopy import Model +from linopy.persistent.diff import ModelDiff +from linopy.restrict import MatrixWindow, RestrictionError, WindowMatrices + +TIME = pd.RangeIndex(8, name="time") + + +def storage_model() -> Model: + """Two coupled constraints: a shifted balance and an elementwise cap.""" + m = Model() + soc = m.add_variables(lower=0, upper=10, coords=[TIME], name="soc") + charge = m.add_variables(lower=0, upper=5, coords=[TIME], name="charge") + m.add_constraints(soc - 0.9 * soc.shift(time=1) - charge == 0, name="balance") + m.add_constraints(charge <= 4, name="cap") + m.add_objective((charge + soc).sum()) + return m + + +def initial(soc_value: float = 3.0) -> pd.Series: + return pd.Series(np.full(len(TIME), soc_value), index=TIME) + + +def assert_same_window(window: WindowMatrices, reference: WindowMatrices) -> None: + assert window.A.shape == reference.A.shape + assert np.array_equal(window.A.indptr, reference.A.indptr) + assert np.array_equal(window.A.indices, reference.A.indices) + assert np.allclose(window.A.data, reference.A.data) + assert np.allclose(window.b, reference.b) + assert np.allclose(window.c, reference.c) + assert np.allclose(window.lb, reference.lb) + assert np.allclose(window.ub, reference.ub) + assert np.array_equal(window.sense, reference.sense) + + +@pytest.mark.parametrize("start", [0, 1, 3]) +def test_select_matches_restrict(start: int) -> None: + m = storage_model() + window = TIME[start : start + 3] + + matrix = MatrixWindow(m, "time").select(window, fix={"soc": initial()}) + reference = m.restrict(time=window, fix={"soc": initial()}).matrices + + assert_same_window(matrix, reference) + + +def test_select_without_fix_matches_restrict() -> None: + """A model with no cross-window coupling needs no boundary values.""" + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + m.add_constraints(x <= 4, name="cap") + m.add_objective(x.sum()) + + window = TIME[2:5] + matrix = MatrixWindow(m, "time").select(window) + reference = m.restrict(time=window).matrices + + assert_same_window(matrix, reference) + + +def test_window_keeps_parent_labels() -> None: + m = storage_model() + window = TIME[1:4] + selected = MatrixWindow(m, "time").select(window, fix={"soc": initial()}) + + assert selected.vlabels.size == selected.A.shape[1] + assert selected.clabels.size == selected.A.shape[0] + assert np.isin(selected.vlabels, m.variables.label_index.vlabels).all() + + +def test_infeasible_folded_row_raises() -> None: + """A pinned row that folds to a violated constant is rejected.""" + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + m.add_constraints(x == 1, name="eq") + # an absolute-index row: it folds away entirely outside its own timestep + m.add_constraints(x.isel(time=0) == 2, name="pin") + m.add_objective(x.sum()) + + window = TIME[1:3] + fix = {"x": pd.Series(np.zeros(len(TIME)), index=TIME)} + with pytest.raises(RestrictionError, match="unsatisfiable"): + MatrixWindow(m, "time").select(window, fix=fix) + + +def test_binary_integrality_is_preserved() -> None: + """A discrete column stays discrete through the window.""" + m = Model() + x = m.add_variables(binary=True, coords=[TIME], name="x") + m.add_constraints(x >= 0.5, name="half") + m.add_objective(x.sum(), sense="min") + + window = m.restrict(time=TIME[:3]) + window.solve(solver_name="highs", output_flag=False) + + assert set(np.unique(window.solution["x"].values)) <= {0.0, 1.0} + assert window.solution["x"].values.min() == 1.0 + + +def test_variable_scaling_is_undone() -> None: + """A scaled column reports the model's units, not the solver's.""" + m = Model() + x = m.add_variables(lower=0, upper=100, coords=[TIME], scaling=10, name="x") + m.add_constraints(x >= 50, name="lb") + m.add_objective(x.sum(), sense="min") + + window = m.restrict(time=TIME[:3]) + window.solve(solver_name="highs", output_flag=False) + + assert np.allclose(window.solution["x"].values, 50.0) + + +def test_constraint_scaling_is_undone() -> None: + """A scaled row reports the model's dual units.""" + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + m.add_constraints(x >= 5, name="lb", scaling=2) + m.add_objective(x.sum(), sense="min") + + window = m.restrict(time=TIME[:3]) + window.solve(solver_name="highs", output_flag=False) + + assert np.allclose(window.dual["lb"].values, 1.0) + + +def test_subset_carrier_raises() -> None: + """Every carrier of the selected dimension must share its index.""" + m = Model() + x = m.add_variables(coords=[TIME], name="x") + subset = m.add_variables(coords=[TIME[::2]], name="subset_x") + m.add_constraints(x <= 1, name="cap") + m.add_constraints(subset <= 1, name="subset_cap") + m.add_objective(m.variables["x"].sum() + subset.sum()) + + with pytest.raises(RestrictionError, match="not carried by every container"): + MatrixWindow(m, "time").select(TIME[1:3]) + + +def linear_model() -> Model: + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + m.add_constraints(x <= 5, name="cap") + m.add_objective(x.sum(), sense="min") + return m + + +def test_restrict_rejects_multiple_dimensions() -> None: + m = linear_model() + m.add_variables(coords=[pd.RangeIndex(2, name="other")], name="y") + + with pytest.raises(NotImplementedError, match="one selected dimension"): + m.restrict(time=TIME[:2], other=[0]) + + +def test_restrict_rejects_unknown_dimension() -> None: + with pytest.raises(RestrictionError, match="does not exist"): + linear_model().restrict(period=[0]) + + +def test_window_rejects_quadratic_objective() -> None: + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + m.add_constraints(x <= 5, name="cap") + m.add_objective((x * x).sum()) + + with pytest.raises(NotImplementedError, match="quadratic"): + MatrixWindow(m, "time") + + +def test_window_rejects_indicator_constraints() -> None: + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + flag = m.add_variables(binary=True, coords=[TIME], name="flag") + m.add_indicator_constraints(flag, 1, x, "<=", 1, name="ind") + m.add_objective(x.sum()) + + with pytest.raises(NotImplementedError, match="indicator"): + MatrixWindow(m, "time") + + +def test_window_requires_a_constraint() -> None: + m = Model() + m.add_variables(lower=0, coords=[TIME], name="x") + m.add_objective(m.variables["x"].sum()) + + with pytest.raises(NotImplementedError, match="at least one constraint"): + MatrixWindow(m, "time") + + +def test_fix_must_name_a_variable() -> None: + with pytest.raises(RestrictionError, match="does not name a variable"): + MatrixWindow(linear_model(), "time").select(TIME[:2], fix={"nope": 1.0}) + + +def test_fix_must_match_the_variable_shape() -> None: + with pytest.raises(RestrictionError, match="do not match its dimensions"): + MatrixWindow(linear_model(), "time").select(TIME[:2], fix={"x": [1.0, 2.0]}) + + +def test_advance_must_name_every_selected_dimension() -> None: + window = linear_model().restrict(time=TIME[:2]) + with pytest.raises(RestrictionError, match="missing"): + window.advance(other=[0]) + + +def test_advance_without_selection_keeps_the_window() -> None: + window = linear_model().restrict(time=TIME[:2]) + before = window.matrices.b.copy() + + window.advance() + + assert np.array_equal(window.matrices.b, before) + + +def test_solve_updates_options_on_later_calls() -> None: + window = linear_model().restrict(time=TIME[:2]) + window.solve(solver_name="highs", output_flag=False) + window.solve(output_flag=False) + window.solve() + + assert window.objective.value == pytest.approx(0.0) + + +def test_solution_before_solve_raises() -> None: + window = linear_model().restrict(time=TIME[:2]) + with pytest.raises(AttributeError, match="has not been solved"): + window.solution + + +def test_dual_size_mismatch_raises() -> None: + window = linear_model().restrict(time=TIME[:2]) + window.solve(solver_name="highs", output_flag=False) + window._dual = np.empty(0) + + with pytest.raises(AttributeError, match="no dual values"): + window.dual + + +def test_getitem_before_solve_returns_nan() -> None: + window = linear_model().restrict(time=TIME[:2]) + values = window["x"] + + assert values.sizes["time"] == 2 + assert np.isnan(values.values).all() + + +def test_getitem_after_solve_returns_values() -> None: + window = linear_model().restrict(time=TIME[:2]) + window.solve(solver_name="highs", output_flag=False) + + assert np.allclose(window["x"].values, 0.0) + + +def test_getitem_unknown_variable_raises() -> None: + window = linear_model().restrict(time=TIME[:2]) + with pytest.raises(KeyError): + window["nope"] + + +def test_model_view_exposes_the_builder_surface() -> None: + window = linear_model().restrict(time=TIME[:2]) + view = window.as_model_view() + + view.constraints.sanitize_missings() + view.constraints.sanitize_zeros() + view.constraints.sanitize_infinities() + + assert view.type == "LP" + assert view.sense == "min" + assert view.is_quadratic is False + assert view.indicator_constraints == () + assert view.variables.sos == {} + assert view.continuous + assert not view.binaries and not view.integers and not view.semi_continuous + assert view.variables.label_index.vlabels.size == len(window.matrices.vlabels) + assert view.constraints.label_index.clabels.size == len(window.matrices.clabels) + assert view._xCounter == len(window.matrices.vlabels) + assert view._cCounter == len(window.matrices.clabels) + assert view.nvars == len(window.matrices.vlabels) + assert view.ncons == len(window.matrices.clabels) + assert len(view.constraints) == len(window.matrices.con_blocks) + + +def test_model_view_type_is_milp_for_mixed_variables() -> None: + m = Model() + x = m.add_variables(lower=0, coords=[TIME], name="x") + flag = m.add_variables(binary=True, coords=[TIME], name="flag") + m.add_constraints(x + flag <= 5, name="cap") + m.add_objective(x.sum() + flag.sum()) + + assert m.restrict(time=TIME[:2]).as_model_view().type == "MILP" + + +def test_model_view_type_is_ilp_for_discrete_only() -> None: + m = Model() + flag = m.add_variables(binary=True, coords=[TIME], name="flag") + m.add_constraints(flag <= 1, name="cap") + m.add_objective(flag.sum()) + + assert m.restrict(time=TIME[:2]).as_model_view().type == "ILP" + + +def _window_matrices() -> WindowMatrices: + return linear_model().restrict(time=TIME[:2]).matrices + + +def test_from_windows_reports_coefficient_changes() -> None: + previous = _window_matrices() + matrix = previous.A.copy() + matrix.data = matrix.data + 1.0 + + diff = ModelDiff.from_windows(previous, replace(previous, A=matrix)) + + assert diff.n_coef_updates > 0 + + +def test_from_windows_reports_variable_type_changes() -> None: + previous = _window_matrices() + vtypes = previous.vtypes.copy() + vtypes[0] = "B" + + diff = ModelDiff.from_windows(previous, replace(previous, vtypes=vtypes)) + + assert diff.var_type_positions.tolist() == [0] + + +def test_from_windows_reports_constraint_sign_changes() -> None: + previous = _window_matrices() + sense = previous.sense.copy() + sense[0] = ">=" + + diff = ModelDiff.from_windows(previous, replace(previous, sense=sense)) + + assert diff.con_sign_indices.tolist() == [0] + + +def test_from_windows_reports_objective_sense_changes() -> None: + previous = _window_matrices() + + diff = ModelDiff.from_windows(previous, replace(previous, objective_sense="max")) + + assert diff.obj_sense == "max" + + +def test_from_windows_reports_bounds_and_cost_changes() -> None: + previous = _window_matrices() + + diff = ModelDiff.from_windows( + previous, replace(previous, lb=previous.lb + 1.0, c=previous.c + 1.0) + ) + + assert diff.var_bounds_indices.size == previous.lb.size + assert diff.obj_c_indices is not None + + +def test_select_without_fix_rejects_dangling_terms() -> None: + """A window that cuts a coupling must be given the boundary values.""" + with pytest.raises(RestrictionError, match="dangling terms"): + MatrixWindow(storage_model(), "time").select(TIME[1:3]) + + +def test_nan_fix_counts_as_missing() -> None: + fix = {"soc": pd.Series(np.full(len(TIME), np.nan), index=TIME)} + with pytest.raises(RestrictionError, match="dangling terms"): + MatrixWindow(storage_model(), "time").select(TIME[1:3], fix=fix) + + +WINDOW_SOLVERS = ("highs", "gurobi", "cplex", "scip", "xpress", "mosek") + + +@pytest.mark.parametrize("solver_name", WINDOW_SOLVERS) +def test_window_solve_across_solvers(solver_name: str) -> None: + """The window recipe runs on every direct solver available in the environment.""" + from linopy import licensed_solvers + + if solver_name not in licensed_solvers: + pytest.skip(f"{solver_name} is not available with a license") + + m = Model() + x = m.add_variables(lower=0, upper=10, coords=[TIME], name="x") + y = m.add_variables(lower=0, coords=[TIME], name="y") + m.add_constraints(x - 0.9 * x.shift(time=1) - y == 0, name="balance") + m.add_objective(y.sum()) + fix = {"x": pd.Series(np.full(len(TIME), 1.0), index=TIME)} + + window = m.restrict(time=TIME[:3], fix=fix) + try: + window.solve(solver_name=solver_name) + except NotImplementedError as exc: # a file-only backend + pytest.skip(str(exc)) + assert window.objective.value == pytest.approx(0.0) + + window.advance(time=TIME[3:6], fix=fix) + window.solve(solver_name=solver_name) + assert window.solution["x"].sizes["time"] == 3 + + +def test_solve_refuses_a_rebuild_when_asked() -> None: + from linopy.persistent.errors import RebuildRequiredError + + m = linear_model() + fix = {"x": pd.Series(np.zeros(len(TIME)), index=TIME)} + window = m.restrict(time=TIME[:3], fix=fix) + window.solve(solver_name="highs", disallow_rebuild=True) + + window.advance(time=TIME[3:5], fix=fix) # a different width changes the shape + with pytest.raises(RebuildRequiredError): + window.solve(solver_name="highs", disallow_rebuild=True) + + +def test_solve_rebuilds_a_changed_shape() -> None: + m = linear_model() + fix = {"x": pd.Series(np.zeros(len(TIME)), index=TIME)} + window = m.restrict(time=TIME[:3], fix=fix) + window.solve(solver_name="highs") + + window.advance(time=TIME[3:5], fix=fix) + window.solve(solver_name="highs") + + assert window.solution["x"].sizes["time"] == 2 + + +def test_window_rejects_sos_constraints() -> None: + m = linear_model() + m.add_sos_constraints(m.variables["x"], 1, "time") + + with pytest.raises(NotImplementedError, match="SOS"): + MatrixWindow(m, "time") + + +def test_window_solves_piecewise_formulations() -> None: + """Piecewise expands into ordinary variables and constraints.""" + from linopy.piecewise import add_piecewise_formulation + + m = Model() + x = m.add_variables(lower=0, upper=10, coords=[TIME], name="x") + y = m.add_variables(lower=0, coords=[TIME], name="y") + m.add_constraints(x <= 5, name="cap") + m.add_objective(y.sum()) + add_piecewise_formulation(m, (x, [0, 5, 10]), (y, [0, 4, 9])) + + window = m.restrict(time=TIME[:2]) + window.solve(solver_name="highs", output_flag=False) + + assert np.allclose(window.solution["y"].values, 0.0) + + +def test_window_rejects_horizon_wide_sos_rows() -> None: + """Reformulating SOS adds a horizon-wide cardinality row, which cannot fold.""" + m = Model() + x = m.add_variables(lower=0, upper=10, coords=[TIME], name="x") + m.add_constraints(x <= 5, name="cap") + m.add_objective(x.sum()) + m.add_sos_constraints(x, 1, "time") + m.apply_sos_reformulation() + + with pytest.raises(RestrictionError, match="dangling terms"): + m.restrict(time=TIME[:2]) diff --git a/test/test_solvers.py b/test/test_solvers.py index 1ca85ed06..57f775094 100644 --- a/test/test_solvers.py +++ b/test/test_solvers.py @@ -6,6 +6,7 @@ """ from pathlib import Path +from typing import Any from unittest.mock import MagicMock import numpy as np @@ -220,7 +221,7 @@ def test_copt_closes_env_per_solve( closes: list[object] = [] original = coptpy.Envr.close - def close(self: object) -> None: + def close(self: Any) -> None: closes.append(self) original(self) diff --git a/test/test_windowing_integration.py b/test/test_windowing_integration.py new file mode 100644 index 000000000..158b6f99d --- /dev/null +++ b/test/test_windowing_integration.py @@ -0,0 +1,277 @@ +#!/usr/bin/env python3 +""" +Auto-executed windowing-equivalence integration tests. + +Each case registers a pre-existing linopy model — from the test suite, the +examples module or the benchmark suite — via the :func:`windowed` decorator. +A single parametrized driver then runs every case under both pattern modes +(strict and union), with the suite's ``semantics`` fixture automatically +multiplying everything by legacy and v1 arithmetic, and assesses equivalence +against the monolithic solve: + +* **decomposable** cases (no cross-tile terms): the chunked objectives sum to + the monolithic objective, and primal and dual values match per position; +* **replay** cases (coupled models): the :meth:`Model.advance` chain must be + indistinguishable from re-solving every window as a fresh + :meth:`Model.restrict` with the same carried fix sequence. +""" + +from __future__ import annotations + +from collections.abc import Callable +from dataclasses import dataclass +from typing import Any, Literal + +import numpy as np +import pandas as pd +import pytest +import xarray as xr +from numpy.testing import assert_allclose + +from linopy import Model, available_solvers +from linopy.examples import benchmark_model + +needs_highs = pytest.mark.skipif( + "highs" not in available_solvers, reason="highs solver not available" +) + + +@dataclass(frozen=True) +class WindowCase: + """A model builder registered for windowing-equivalence testing.""" + + name: str + build: Callable[[], Model] + dim: str + tile: int + coupling: Literal["decomposable", "replay"] + fix_vars: tuple[str, ...] = () + + +CASES: list[WindowCase] = [] + + +def windowed( + dim: str, + tile: int, + coupling: Literal["decomposable", "replay"], + fix_vars: tuple[str, ...] = (), +) -> Callable[[Callable[[], Model]], Callable[[], Model]]: + """Register a model builder as a windowing-equivalence case.""" + + def deco(build: Callable[[], Model]) -> Callable[[], Model]: + CASES.append(WindowCase(build.__name__, build, dim, tile, coupling, fix_vars)) + return build + + return deco + + +# --------------------------------------------------------------------------- +# cases: pre-existing models, sources cited +# --------------------------------------------------------------------------- + + +@windowed( + dim="snapshot", + tile=3, + coupling="decomposable", + fix_vars=("x",), +) +def mi_period_demand() -> Model: + """ + From ``test/test_mi_feasibility.py``: per-period demand via groupby-sum. + Windows align with the period boundaries; a period's demand row leaks + into later windows only as a satisfied folded constant. + """ + m = Model() + x = m.add_variables( + lower=0, upper=4.0, coords=[pd.RangeIndex(6, name="snapshot")], name="x" + ) + period = xr.DataArray( + np.repeat([2020, 2030], 3), + dims="snapshot", + coords={"snapshot": x.indexes["snapshot"]}, + ) + e = (1.0 * x).assign_coords(period=period) + rhs = xr.DataArray([5.0, 7.0], dims="period", coords={"period": [2020, 2030]}) + m.add_constraints(e.groupby("period").sum() >= rhs) + m.add_objective((np.arange(1.0, 7.0) * x).sum()) + return m + + +@windowed(dim="time", tile=1, coupling="decomposable") +def constraint_empty_rows() -> Model: + """From ``test/test_constraint.py``: an empty constraint group solves.""" + m = Model() + x = m.add_variables( + lower=0.0, coords=[range(3), range(2)], dims=["time", "product"], name="x" + ) + m.add_constraints(x.isel(time=range(1, 1)) == 0, name="empty") + m.add_objective(x.sum()) + return m + + +@windowed(dim="dim_0", tile=5, coupling="decomposable") +def examples_benchmark() -> Model: + """From ``linopy/examples.py``: the JOSS benchmark model.""" + return benchmark_model() + + +@windowed( + dim="time", + tile=24, + coupling="replay", + fix_vars=("soc",), +) +def benchmarks_storage() -> Model: + """ + Replicated from ``benchmarks/models/storage.py`` (not importable from + the test tree): the 168h shift-coupled SoC fleet — a real, non-test toy + model with intertemporal coupling. + """ + m = Model() + storages = pd.RangeIndex(10, name="storage") + time = pd.RangeIndex(168, name="time") + soc = m.add_variables(lower=0, upper=100, coords=[storages, time], name="soc") + charge = m.add_variables(lower=0, upper=50, coords=[storages, time], name="charge") + discharge = m.add_variables( + lower=0, upper=50, coords=[storages, time], name="discharge" + ) + prev = soc.shift(time=1) + m.add_constraints( + soc - 0.99 * prev - 0.95 * charge + discharge / 0.95 == 0, name="soc_balance" + ) + m.add_objective((charge + discharge).sum()) + return m + + +# --------------------------------------------------------------------------- +# the harness +# --------------------------------------------------------------------------- + + +def _dim_index(model: Model, case: WindowCase) -> pd.Index: + """The case's windowing dimension index, from any carrier.""" + for _, var in model.variables.items(): + if case.dim in var.indexes: + return var.indexes[case.dim] + raise ValueError(f"No variable carries the dimension '{case.dim}'.") + + +def _initial_carried(model: Model, case: WindowCase) -> dict[str, xr.DataArray]: + """ + Full-horizon carry arrays for the case's fix variables, initialized to a + feasible guess (the variable's upper bound, 0 where unbounded): folded + constant rows from not-yet-solved windows must come out satisfied. + """ + carried = {} + for name in case.fix_vars: + var = model.variables[name] + labels = var.labels + upper = np.asarray(var.upper.values, dtype=float) + fill = np.where(np.isfinite(upper), upper, 0.0) + carried[name] = xr.DataArray( + fill, coords=labels.coords, dims=labels.dims, name=name + ) + return carried + + +def _solve_windows( + model: Model, + case: WindowCase, + windows: list[pd.Index], + carried: dict[str, xr.DataArray], +) -> list[dict[str, Any]]: + """Solve the windows sequentially via restrict + advance.""" + fix = {name: carried[name] for name in case.fix_vars} + w = model.restrict(**{case.dim: windows[0]}, fix=fix) + results = [] + for i, window in enumerate(windows): + if i > 0: + w.advance(**{case.dim: window}, fix=fix) + w.solve(solver_name="highs", output_flag=False) + results.append( + { + "objective": w.objective.value, + "solution": w.solution, + "dual": w.dual, + } + ) + for name in case.fix_vars: + carried[name].loc[{case.dim: window}] = w.solution[name].values + return results + + +def _assert_chunked_equals_full( + results: list[dict[str, Any]], full: Model, case: WindowCase +) -> None: + """ + Objectives sum to the monolithic one; primal/dual values match per + position for every array carrying the windowing dim. + + Arrays WITHOUT the dim (shared variables re-optimized per window, + aggregate rows enforced per window) are per-window quantities: their + equality with the monolithic solve is not implied by decomposability and + is not asserted. + """ + total = sum(r["objective"] for r in results) + assert total == pytest.approx(full.objective.value) + + for name, da in full.solution.data_vars.items(): + if case.dim not in da.dims: + continue + parts = [r["solution"][name] for r in results] + merged = xr.concat(parts, dim=da.indexes[case.dim]) + assert_allclose(merged.values, da.values, rtol=1e-7) + for name, da in full.dual.data_vars.items(): + if case.dim not in da.dims: + continue + parts = [r["dual"][name] for r in results] + merged = xr.concat(parts, dim=da.indexes[case.dim]) + assert_allclose(merged.values, da.values, rtol=1e-6) + + +def _assert_replay_equals_fresh( + case: WindowCase, + windows: list[pd.Index], + carried_init: dict[str, xr.DataArray], + results: list[dict[str, Any]], +) -> None: + """ + Replay the chain's fix sequence with fresh restricts: every window must + re-solve to the chain's own result. + """ + carried = {name: da.copy() for name, da in carried_init.items()} + for i, window in enumerate(windows): + fix = {name: carried[name] for name in case.fix_vars} + fresh = case.build().restrict(**{case.dim: window}, fix=fix) + fresh.solve(solver_name="highs", output_flag=False) + + assert fresh.objective.value == pytest.approx(results[i]["objective"]) + for name, da in fresh.solution.data_vars.items(): + assert_allclose(da.values, results[i]["solution"][name].values, rtol=1e-7) + for name, da in fresh.dual.data_vars.items(): + assert_allclose(da.values, results[i]["dual"][name].values, rtol=1e-6) + for name in case.fix_vars: + carried[name].loc[{case.dim: window}] = fresh.solution[name].values + + +@needs_highs +@pytest.mark.parametrize("case", CASES, ids=lambda c: c.name) +def test_windowed_equivalence(case: WindowCase, semantics: str) -> None: + """Every registered model is solved equivalently through the windows.""" + full = case.build() + full.solve(solver_name="highs", output_flag=False) + + index = _dim_index(full, case) + n_tiles = int(np.ceil(len(index) / case.tile)) + windows = [index[i * case.tile : (i + 1) * case.tile] for i in range(n_tiles)] + + model = case.build() + carried = _initial_carried(model, case) + results = _solve_windows(model, case, windows, carried) + + if case.coupling == "decomposable": + _assert_chunked_equals_full(results, full, case) + else: + _assert_replay_equals_fresh(case, windows, carried, results) From f78d4ca8be129e6f50562a8ab0b790cb6210daf8 Mon Sep 17 00:00:00 2001 From: Olivier Lacroix Date: Sat, 24 Oct 2026 13:19:54 +1100 Subject: [PATCH 2/2] test(restrict): drop the removed semantics fixture from the windowing test Master removed the autouse semantics fixture with legacy semantics (#1041). The windowing equivalence test still requested it and errored at setup on every platform. Co-Authored-By: Claude Code --- test/test_windowing_integration.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_windowing_integration.py b/test/test_windowing_integration.py index 158b6f99d..49392b1b7 100644 --- a/test/test_windowing_integration.py +++ b/test/test_windowing_integration.py @@ -5,9 +5,7 @@ Each case registers a pre-existing linopy model — from the test suite, the examples module or the benchmark suite — via the :func:`windowed` decorator. A single parametrized driver then runs every case under both pattern modes -(strict and union), with the suite's ``semantics`` fixture automatically -multiplying everything by legacy and v1 arithmetic, and assesses equivalence -against the monolithic solve: +(strict and union), and assesses equivalence against the monolithic solve: * **decomposable** cases (no cross-tile terms): the chunked objectives sum to the monolithic objective, and primal and dual values match per position; @@ -258,7 +256,7 @@ def _assert_replay_equals_fresh( @needs_highs @pytest.mark.parametrize("case", CASES, ids=lambda c: c.name) -def test_windowed_equivalence(case: WindowCase, semantics: str) -> None: +def test_windowed_equivalence(case: WindowCase) -> None: """Every registered model is solved equivalently through the windows.""" full = case.build() full.solve(solver_name="highs", output_flag=False)