From 31990c8eb8b8ae58710424d0eb33fe2fb74039fb Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Wed, 1 Jul 2026 11:27:00 -0700 Subject: [PATCH 1/2] fix: widen frozen constraint CSR to current active-variable count Frozen (CSR-backed) constraints cache their matrix at freeze time, sized against the active-variable set as it was then. Adding variables afterwards (e.g. an auxiliary variable introduced by a constraint added after this one was frozen, common under Model(freeze_constraints=True) when building incrementally) grows the active-variable count, but CSRConstraint.to_matrix and to_matrix_with_rhs returned the cached CSR verbatim, ignoring the VariableLabelIndex passed in. The earlier-frozen blocks stay narrower than later ones, so stacking them fails with "ValueError: incompatible dimensions for axis 1". This breaks the invariant documented in #630 -- to_matrix on both the mutable and frozen constraint "always returns a csr matrix of shape (n_active_cons, n_active_vars)" -- which the VariableLabelIndex exists precisely to uphold. Widen the cached CSR to the current n_active_vars in a shared _csr_for_index helper used by both to_matrix (matrix-accessor path) and to_matrix_with_rhs (direct-solve path). Variable positions are assigned in encounter order and are stable under additions, so widening only appends empty trailing columns. Adds differential tests (freeze vs mutable must agree) covering the baseline, variables added after freeze, and constraints linking later-added variables. --- doc/release_notes.rst | 1 + linopy/constraints.py | 40 ++++++- test/test_frozen_constraint_active_vars.py | 121 +++++++++++++++++++++ 3 files changed, 159 insertions(+), 3 deletions(-) create mode 100644 test/test_frozen_constraint_active_vars.py diff --git a/doc/release_notes.rst b/doc/release_notes.rst index 80f9076e..c8bb46c7 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -31,6 +31,7 @@ Upcoming Version * LP file export now honors bounds tightened below ``[0, 1]`` on a binary variable via the ``.lower``/``.upper`` setters after creation (e.g. ``upper = 0``). Previously such bounds were written only by ``io_api="direct"`` and dropped by ``io_api="lp"``. (https://github.com/PyPSA/linopy/issues/776) * Freezing an empty constraint group (e.g. an empty ``isel`` slice) no longer raises ``ValueError: cannot reshape array of size 0``. ``Model(freeze_constraints=True)`` and ``Constraint.freeze()`` now round-trip zero-row constraints losslessly. +* Adding variables after freezing a constraint (e.g. under ``Model(freeze_constraints=True)`` while building incrementally) no longer fails with ``ValueError: incompatible dimensions for axis 1`` when assembling the constraint matrix. Frozen (CSR-backed) constraints now widen their cached CSR to the current active-variable count, matching the ``shape (n_active_cons, n_active_vars)`` invariant of the mutable path. * ``Variable.where`` no longer raises ``ValueError: exact match required for all data variable names`` once a solution is attached (after ``Model.solve``) or the variable is fixed. The fill value now covers auxiliary data variables (``solution``, stashed bounds) instead of only ``labels``/``lower``/``upper``. Version 0.8.0 diff --git a/linopy/constraints.py b/linopy/constraints.py index 0b9dbb0a..72258dc6 100644 --- a/linopy/constraints.py +++ b/linopy/constraints.py @@ -852,11 +852,40 @@ def row_expr(row: int) -> str: return "\n".join(lines) + def _csr_for_index( + self, label_index: VariableLabelIndex | None + ) -> scipy.sparse.csr_array: + """ + Return the cached CSR widened to the current active-variable count. + + The CSR is cached at freeze time with as many columns as there were + active variables then. If variables were added afterwards (e.g. an + auxiliary variable introduced by a constraint added after this one was + frozen), the active-variable count has grown and the cached CSR is too + narrow to be stacked with constraints frozen later, so + ``Constraints.to_matrix`` fails with an axis-1 dimension mismatch. + + Variable positions are assigned in encounter order and are stable under + additions, so widening the column dimension only appends empty trailing + columns and keeps every per-constraint block consistent with + ``label_index.n_active_vars``. + """ + csr = self._csr + if label_index is None: + return csr + n_active_vars = label_index.n_active_vars + if csr.shape[1] < n_active_vars: + csr = scipy.sparse.csr_array( + (csr.data, csr.indices, csr.indptr), + shape=(csr.shape[0], n_active_vars), + ) + return csr + def to_matrix( self, label_index: VariableLabelIndex | None = None ) -> tuple[scipy.sparse.csr_array, np.ndarray]: """Return the stored CSR matrix and con_labels.""" - return self._csr, self._con_labels + return self._csr_for_index(label_index), self._con_labels def to_netcdf_ds(self) -> Dataset: """Return a Dataset with raw CSR components for netcdf serialization.""" @@ -945,12 +974,17 @@ def has_variable(self, variable: variables.Variable) -> bool: def to_matrix_with_rhs( self, label_index: VariableLabelIndex ) -> tuple[scipy.sparse.csr_array, np.ndarray, np.ndarray, np.ndarray]: - """Return (csr, con_labels, b, sense) — all pre-stored, no recomputation.""" + """ + Return (csr, con_labels, b, sense) — pre-stored, widened if needed. + + See ``_csr_for_index`` for why the cached CSR may need widening to the + current active-variable count before it can be stacked. + """ if isinstance(self._sign, str): sense = np.full(len(self._rhs), self._sign[0]) else: sense = np.array([s[0] for s in self._sign]) - return self._csr, self._con_labels, self._rhs, sense + return self._csr_for_index(label_index), self._con_labels, self._rhs, sense def active_labels(self) -> np.ndarray: return self._con_labels diff --git a/test/test_frozen_constraint_active_vars.py b/test/test_frozen_constraint_active_vars.py new file mode 100644 index 00000000..1c6413d3 --- /dev/null +++ b/test/test_frozen_constraint_active_vars.py @@ -0,0 +1,121 @@ +""" +Differential tests for frozen (CSR-backed) constraints: adding variables. + +The invariant under test (stated in PR #630): ``to_matrix`` on both the mutable +``Constraint`` and the frozen ``CSRConstraint`` "always returns a csr matrix of +shape (n_active_constraints, n_active_variables)". A frozen constraint caches +its CSR at freeze time, so adding variables *after* freezing must still yield a +matrix consistent with the mutable path. + +Strategy: build the same model twice under an identical construction order -- +once with ``freeze_constraints=False`` (mutable, treated as ground truth) and +once with ``freeze_constraints=True`` -- and assert the assembled constraint +matrix, shape, and solved objective agree. Any divergence is a frozen-path bug. +""" + +import numpy as np +import pandas as pd +import pytest + +import linopy +from linopy import Model + + +def _assemble(freeze: bool) -> Model: + """ + Return a model built in a fixed order that freezes before adding a variable. + + Order: add x, add a constraint on x, THEN add y, add a constraint on y. + With freeze_constraints=True the first constraint is frozen while only x + exists. + """ + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(3, name="i")], name="x") + m.add_constraints(x >= 1, name="cx") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="j")], name="y") + m.add_constraints(y >= 2, name="cy") + m.add_objective(x.sum() + y.sum()) + return m + + +def _dense_A(m: Model) -> np.ndarray: + A = m.matrices.A + assert A is not None + return np.asarray(A.todense()) + + +# --------------------------------------------------------------------------- +# Scenario A: baseline -- all variables added before any frozen constraint. +# This is what upstream fixtures do; should pass with or without the fix. +# --------------------------------------------------------------------------- +def test_A_all_vars_first() -> None: + m_mut = Model() + x = m_mut.add_variables(lower=0, coords=[pd.RangeIndex(3, name="i")], name="x") + y = m_mut.add_variables(lower=0, coords=[pd.RangeIndex(2, name="j")], name="y") + m_mut.add_constraints(x >= 1, name="cx", freeze=False) + m_mut.add_constraints(y >= 2, name="cy", freeze=False) + + m_frz = Model() + x = m_frz.add_variables(lower=0, coords=[pd.RangeIndex(3, name="i")], name="x") + y = m_frz.add_variables(lower=0, coords=[pd.RangeIndex(2, name="j")], name="y") + m_frz.add_constraints(x >= 1, name="cx", freeze=True) + m_frz.add_constraints(y >= 2, name="cy", freeze=True) + + np.testing.assert_array_equal(_dense_A(m_mut), _dense_A(m_frz)) + + +# --------------------------------------------------------------------------- +# Scenario B: interleaved -- variable added AFTER a frozen constraint. +# This is the widening case. Expected to FAIL on unpatched master. +# --------------------------------------------------------------------------- +def test_B_var_added_after_freeze_matrix_shape() -> None: + m = _assemble(freeze=True) + n_vars = m.variables.label_index.n_active_vars + n_cons = m.constraints.label_index.n_active_cons + A = m.matrices.A + assert A is not None + assert A.shape == (n_cons, n_vars) + + +def test_B_interleaved_matrix_matches_mutable() -> None: + A_mut = _dense_A(_assemble(freeze=False)) + A_frz = _dense_A(_assemble(freeze=True)) + assert A_mut.shape == A_frz.shape + np.testing.assert_array_equal(A_mut, A_frz) + + +def test_B_interleaved_solves_equal() -> None: + if "highs" not in linopy.available_solvers: + pytest.skip("highs unavailable") + m_mut = _assemble(freeze=False) + m_frz = _assemble(freeze=True) + m_mut.solve(solver_name="highs") + m_frz.solve(solver_name="highs") + obj_mut = m_mut.objective.value + obj_frz = m_frz.objective.value + assert obj_mut is not None and obj_frz is not None + assert obj_frz == pytest.approx(obj_mut) + + +# --------------------------------------------------------------------------- +# Scenario C: a frozen constraint that links a variable added later, plus an +# extra late variable that widens the index further. +# --------------------------------------------------------------------------- +def _assemble_linking(freeze: bool) -> Model: + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="x") + m.add_constraints(x >= 1, name="cx") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="y") + # constraint linking x and y, frozen while both exist + m.add_constraints(x + y >= 5, name="cxy") + z = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="z") + m.add_constraints(z >= 3, name="cz") + m.add_objective(x.sum() + y.sum() + z.sum()) + return m + + +def test_C_linking_matrix_matches_mutable() -> None: + A_mut = _dense_A(_assemble_linking(freeze=False)) + A_frz = _dense_A(_assemble_linking(freeze=True)) + assert A_mut.shape == A_frz.shape + np.testing.assert_array_equal(A_mut, A_frz) From b2b622ccb327b9a65fae6ed4433781653ca34e90 Mon Sep 17 00:00:00 2001 From: Koen van Greevenbroek Date: Wed, 1 Jul 2026 11:27:37 -0700 Subject: [PATCH 2/2] fix: remap frozen constraint columns when a variable is removed Frozen (CSR-backed) constraints store absolute dense variable positions from freeze time. Removing a variable renumbers the positions of every later variable, but Model.remove_variables left surviving frozen constraints with those stale positions -- both the wrong width and the wrong columns -- producing a corrupt constraint matrix (e.g. shape (2, 5) instead of (2, 2)) or an IndexError at solve time. The mutable Constraint path handles removal correctly, so this was a frozen-vs-mutable parity gap. remove_variables now captures the pre-removal variable ordering and, once the variable (and any constraints referencing it) are gone, remaps each surviving frozen constraint's CSR columns through variable labels via the new CSRConstraint._remap_columns. The old ordering is still available at the mutation site, so no per-constraint label storage is needed; a frozen constraint that referenced the removed variable has already been dropped, so every remaining column maps to a valid new position. Extends the differential tests with removal scenarios: first/middle-variable removal, straddling references, sequential removals, remove-then-add, and the drop-referencing-constraint parity case. --- doc/release_notes.rst | 1 + linopy/constraints.py | 33 +++++ linopy/model.py | 15 ++ test/test_frozen_constraint_active_vars.py | 152 ++++++++++++++++++++- 4 files changed, 198 insertions(+), 3 deletions(-) diff --git a/doc/release_notes.rst b/doc/release_notes.rst index c8bb46c7..b5ed587a 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -32,6 +32,7 @@ Upcoming Version * LP file export now honors bounds tightened below ``[0, 1]`` on a binary variable via the ``.lower``/``.upper`` setters after creation (e.g. ``upper = 0``). Previously such bounds were written only by ``io_api="direct"`` and dropped by ``io_api="lp"``. (https://github.com/PyPSA/linopy/issues/776) * Freezing an empty constraint group (e.g. an empty ``isel`` slice) no longer raises ``ValueError: cannot reshape array of size 0``. ``Model(freeze_constraints=True)`` and ``Constraint.freeze()`` now round-trip zero-row constraints losslessly. * Adding variables after freezing a constraint (e.g. under ``Model(freeze_constraints=True)`` while building incrementally) no longer fails with ``ValueError: incompatible dimensions for axis 1`` when assembling the constraint matrix. Frozen (CSR-backed) constraints now widen their cached CSR to the current active-variable count, matching the ``shape (n_active_cons, n_active_vars)`` invariant of the mutable path. +* Removing a variable after freezing constraints now remaps the cached CSR columns of the surviving frozen constraints. Previously ``Model.remove_variables`` left their absolute variable positions stale after the removal renumbered them, producing a wrong-width constraint matrix or an ``IndexError`` at solve time. * ``Variable.where`` no longer raises ``ValueError: exact match required for all data variable names`` once a solution is attached (after ``Model.solve``) or the variable is fixed. The fill value now covers auxiliary data variables (``solution``, stashed bounds) instead of only ``labels``/``lower``/``upper``. Version 0.8.0 diff --git a/linopy/constraints.py b/linopy/constraints.py index 72258dc6..2b1cca47 100644 --- a/linopy/constraints.py +++ b/linopy/constraints.py @@ -881,6 +881,39 @@ def _csr_for_index( ) return csr + def _remap_columns( + self, old_vlabels: np.ndarray, new_label_index: VariableLabelIndex + ) -> None: + """ + Rewrite cached CSR column positions after variables were removed. + + The cached CSR stores absolute dense variable positions from freeze + time. Removing a variable renumbers the positions of every later + variable, so those cached positions go stale (wrong columns and wrong + width). ``old_vlabels`` gives the variable label at each column position + currently used by the CSR; ``new_label_index`` provides the post-removal + label -> position mapping. A frozen constraint that referenced a removed + variable has already been dropped, so every remaining column maps to a + valid new position. + """ + csr = self._csr + n_active_vars = new_label_index.n_active_vars + if csr.nnz == 0: + self._csr = scipy.sparse.csr_array( + (csr.shape[0], n_active_vars), dtype=csr.dtype + ) + return + new_positions = new_label_index.label_to_pos[old_vlabels[csr.indices]] + if (new_positions < 0).any(): + raise ValueError( + "Frozen constraint references a removed variable while remapping " + "columns; this should not happen." + ) + self._csr = scipy.sparse.csr_array( + (csr.data, new_positions, csr.indptr), + shape=(csr.shape[0], n_active_vars), + ) + def to_matrix( self, label_index: VariableLabelIndex | None = None ) -> tuple[scipy.sparse.csr_array, np.ndarray]: diff --git a/linopy/model.py b/linopy/model.py index de5c089f..544c1310 100644 --- a/linopy/model.py +++ b/linopy/model.py @@ -1281,8 +1281,23 @@ def remove_variables(self, name: str) -> None: for k in to_remove: self.constraints.remove(k) + # Frozen (CSR-backed) constraints cache absolute variable positions. + # Removing a variable renumbers the positions of every later variable, + # so capture the pre-removal ordering to remap the surviving frozen + # constraints once the variable is gone. + frozen_cons = [ + c for _, c in self.constraints.items() if isinstance(c, CSRConstraint) + ] + old_vlabels = self.variables.label_index.vlabels.copy() if frozen_cons else None + self.variables.remove(name) + if frozen_cons: + assert old_vlabels is not None # captured whenever frozen_cons is non-empty + new_label_index = self.variables.label_index + for con in frozen_cons: + con._remap_columns(old_vlabels, new_label_index) + self.objective = self.objective.sel( {TERM_DIM: ~self.objective.vars.isin(variable.labels)} ) diff --git a/test/test_frozen_constraint_active_vars.py b/test/test_frozen_constraint_active_vars.py index 1c6413d3..697868ab 100644 --- a/test/test_frozen_constraint_active_vars.py +++ b/test/test_frozen_constraint_active_vars.py @@ -1,11 +1,12 @@ """ -Differential tests for frozen (CSR-backed) constraints: adding variables. +Differential tests for frozen (CSR-backed) constraints: changing the variables. The invariant under test (stated in PR #630): ``to_matrix`` on both the mutable ``Constraint`` and the frozen ``CSRConstraint`` "always returns a csr matrix of shape (n_active_constraints, n_active_variables)". A frozen constraint caches -its CSR at freeze time, so adding variables *after* freezing must still yield a -matrix consistent with the mutable path. +its CSR at freeze time, so any change to the active-variable set *after* freezing +(adding or removing variables) must still yield a matrix consistent with the +mutable path. Strategy: build the same model twice under an identical construction order -- once with ``freeze_constraints=False`` (mutable, treated as ground truth) and @@ -119,3 +120,148 @@ def test_C_linking_matrix_matches_mutable() -> None: A_frz = _dense_A(_assemble_linking(freeze=True)) assert A_mut.shape == A_frz.shape np.testing.assert_array_equal(A_mut, A_frz) + + +# --------------------------------------------------------------------------- +# Scenario D: variable removal AFTER freezing. Removing a variable renumbers +# dense positions for all later variables. A frozen constraint caches absolute +# positions, so this is the silent-corruption hypothesis. +# +# We remove the FIRST variable (lowest positions) so that a later frozen +# constraint referencing a higher-positioned variable has its cached column +# indices invalidated by the renumbering. +# --------------------------------------------------------------------------- +def _assemble_for_removal(freeze: bool) -> Model: + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(3, name="i")], name="x") + m.add_constraints(x >= 1, name="cx") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="j")], name="y") + m.add_constraints(y >= 7, name="cy") + m.add_objective(x.sum() + y.sum()) + return m + + +def test_D_remove_first_var_matrix_matches_mutable() -> None: + m_mut = _assemble_for_removal(freeze=False) + m_frz = _assemble_for_removal(freeze=True) + # Removing x also removes cx (references x); cy (on y) must survive intact. + m_mut.remove_variables("x") + m_frz.remove_variables("x") + A_mut = _dense_A(m_mut) + A_frz = _dense_A(m_frz) + assert A_mut.shape == A_frz.shape + np.testing.assert_array_equal(A_mut, A_frz) + + +def test_D_remove_first_var_solves_equal() -> None: + if "highs" not in linopy.available_solvers: + pytest.skip("highs unavailable") + m_mut = _assemble_for_removal(freeze=False) + m_frz = _assemble_for_removal(freeze=True) + m_mut.remove_variables("x") + m_frz.remove_variables("x") + m_mut.solve(solver_name="highs") + m_frz.solve(solver_name="highs") + obj_mut = m_mut.objective.value + obj_frz = m_frz.objective.value + assert obj_mut is not None and obj_frz is not None + assert obj_frz == pytest.approx(obj_mut) + + +# --------------------------------------------------------------------------- +# Scenario E: a surviving frozen constraint whose terms STRADDLE the removed +# variable (references a variable before it and one after it). Exercises the +# column remap within a single CSR row and the sorted-indices assumption. +# --------------------------------------------------------------------------- +def _assemble_straddle(freeze: bool) -> Model: + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="x") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="y") + z = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="z") + # constraint links x and z (straddles y), frozen while all three exist + m.add_constraints(x + z >= 4, name="cxz") + m.add_objective(x.sum() + y.sum() + z.sum()) + return m + + +def test_E_straddle_remove_middle_matches_mutable() -> None: + m_mut = _assemble_straddle(freeze=False) + m_frz = _assemble_straddle(freeze=True) + m_mut.remove_variables("y") # y does not appear in cxz -> cxz survives + m_frz.remove_variables("y") + A_mut = _dense_A(m_mut) + A_frz = _dense_A(m_frz) + assert A_mut.shape == A_frz.shape + np.testing.assert_array_equal(A_mut, A_frz) + # sanity: CSR indices remain sorted within rows after remap + A_frz_sp = m_frz.matrices.A + assert A_frz_sp is not None + dense_before = np.asarray(A_frz_sp.todense()) + A_frz_sp.sort_indices() + np.testing.assert_array_equal(np.asarray(A_frz_sp.todense()), dense_before) + + +# --------------------------------------------------------------------------- +# Scenario F: multiple sequential removals, each renumbering again. +# --------------------------------------------------------------------------- +def _assemble_three(freeze: bool) -> Model: + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="x") + m.add_constraints(x >= 1, name="cx") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="y") + m.add_constraints(y >= 2, name="cy") + z = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="z") + m.add_constraints(z >= 3, name="cz") + m.add_objective(x.sum() + y.sum() + z.sum()) + return m + + +def test_F_multiple_sequential_removals_match_mutable() -> None: + m_mut = _assemble_three(freeze=False) + m_frz = _assemble_three(freeze=True) + for name in ("x", "y"): # remove two, leaving only z + cz + m_mut.remove_variables(name) + m_frz.remove_variables(name) + np.testing.assert_array_equal(_dense_A(m_mut), _dense_A(m_frz)) + + +# --------------------------------------------------------------------------- +# Scenario G: interleave removal and addition (remove, then add a new var and +# a new frozen constraint that widens again). +# --------------------------------------------------------------------------- +def test_G_remove_then_add_matches_mutable() -> None: + def build(freeze: bool) -> Model: + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="x") + m.add_constraints(x >= 1, name="cx") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="y") + m.add_constraints(y >= 2, name="cy") + m.remove_variables("x") # renumber: cy remapped + w = m.add_variables(lower=0, coords=[pd.RangeIndex(3, name="i")], name="w") + m.add_constraints(w >= 5, name="cw") # widen again + m.add_objective(y.sum() + w.sum()) + return m + + np.testing.assert_array_equal(_dense_A(build(False)), _dense_A(build(True))) + + +# --------------------------------------------------------------------------- +# Scenario H: removing a variable also drops a frozen constraint referencing it +# even when that constraint references surviving variables too (parity check). +# --------------------------------------------------------------------------- +def test_H_removal_drops_referencing_frozen_constraint() -> None: + def build(freeze: bool) -> Model: + m = Model(freeze_constraints=freeze) + x = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="x") + y = m.add_variables(lower=0, coords=[pd.RangeIndex(2, name="i")], name="y") + m.add_constraints(x + y >= 3, name="cxy") + m.add_constraints(y >= 1, name="cy") + m.add_objective(x.sum() + y.sum()) + return m + + m_mut = build(False) + m_frz = build(True) + m_mut.remove_variables("x") # drops cxy in both; cy survives + m_frz.remove_variables("x") + assert set(m_mut.constraints) == set(m_frz.constraints) + np.testing.assert_array_equal(_dense_A(m_mut), _dense_A(m_frz))