diff --git a/linopy/persistent/diff.py b/linopy/persistent/diff.py index 52058d433..ca38c4509 100644 --- a/linopy/persistent/diff.py +++ b/linopy/persistent/diff.py @@ -2,11 +2,12 @@ import enum from collections.abc import Iterable -from dataclasses import dataclass +from dataclasses import dataclass, replace from functools import cached_property from typing import TYPE_CHECKING import numpy as np +from scipy.sparse import csr_array from linopy.constants import short_GREATER_EQUAL, short_LESS_EQUAL from linopy.constraints import Constraint @@ -460,17 +461,38 @@ def diff_con( self.con_coords[name] = new_coords if not _coords_equal(base_coords, new_coords, self.ignored): return RebuildReason.COORD_REINDEX - if not _same(new_buf.indptr, base_buf.indptr): - return RebuildReason.SPARSITY - if not _same(new_buf.indices, base_buf.indices): - return RebuildReason.SPARSITY - n_rows = new_buf.active_labels.size if n_rows == 0: return None + coef_buf = new_buf changed_rows = None - if not (skip_coef_compare or new_buf.data is base_buf.data): + same_pattern = _same(new_buf.indptr, base_buf.indptr) and _same( + new_buf.indices, base_buf.indices + ) + if not same_pattern: + # Native coefficient setters insert nonzeros and remove zeros. + # Compare by matrix position so culled terms are explicitly cleared. + shape = (n_rows, len(self.var_label_index.vlabels)) + before = csr_array( + (base_buf.data, base_buf.indices, base_buf.indptr), shape=shape + ) + after = csr_array( + (new_buf.data, new_buf.indices, new_buf.indptr), shape=shape + ) + changes = (after - before).tocsr() + changes.eliminate_zeros() + if changes.nnz: + rows = np.repeat(np.arange(n_rows), np.diff(changes.indptr)) + values = np.asarray(after[rows, changes.indices]).ravel() + coef_buf = replace( + new_buf, + indptr=changes.indptr, + indices=changes.indices, + data=values, + ) + changed_rows = np.flatnonzero(np.diff(changes.indptr)) + elif not (skip_coef_compare or new_buf.data is base_buf.data): data_diff = new_buf.data != base_buf.data if data_diff.any(): nnz_per_row = np.diff(new_buf.indptr) @@ -496,10 +518,10 @@ def diff_con( row_positions = self.con_l2p[new_buf.active_labels[changed_rows]].astype( np.int32, copy=False ) - indptr = new_buf.indptr + indptr = coef_buf.indptr nnz = int((indptr[changed_rows + 1] - indptr[changed_rows]).sum()) self.coef_deltas.append( - _CoefDelta(new_buf, changed_rows, row_positions, nnz) + _CoefDelta(coef_buf, changed_rows, row_positions, nnz) ) self._cc_cur += nnz if rhs_idx is not None: diff --git a/linopy/solvers.py b/linopy/solvers.py index 6e1ec74d0..df4df6e95 100644 --- a/linopy/solvers.py +++ b/linopy/solvers.py @@ -818,8 +818,8 @@ def solve( between solves. Pass ``disallow_rebuild=True`` to guarantee that an existing solver - model is updated in place — any condition that would force a rebuild - (structural change, sparsity change, backend rejection, …) raises + model is updated in place. Any condition that would force a rebuild + (structural change or backend rejection) raises :class:`RebuildRequiredError` instead. The initial build on the first ``solve(model, ...)`` is still allowed. diff --git a/test/test_persistent_snapshot_buffers.py b/test/test_persistent_snapshot_buffers.py index fcacb1b8e..95c65cc6c 100644 --- a/test/test_persistent_snapshot_buffers.py +++ b/test/test_persistent_snapshot_buffers.py @@ -57,16 +57,16 @@ def baseline_model() -> Model: return m -def test_shape_mismatch_triggers_sparsity_rebuild(baseline_model: Model) -> None: +def test_added_terms_produce_coefficient_updates(baseline_model: Model) -> None: snap = ModelSnapshot.capture(baseline_model) x = baseline_model.variables["x"] y = baseline_model.variables["y"] baseline_model.constraints["c1"].lhs = 2 * x + 1 * y.sum() diff = ModelDiff.from_snapshot(snap, baseline_model) - assert diff in { - RebuildReason.SPARSITY, - RebuildReason.STRUCTURAL_LABELS, - } + assert isinstance(diff, ModelDiff) + np.testing.assert_array_equal(diff.con_coef_rows, np.repeat([0, 1, 2], 2)) + np.testing.assert_array_equal(diff.con_coef_cols, np.tile([3, 4], 3)) + np.testing.assert_array_equal(diff.con_coef_vals, np.ones(6)) def test_zero_coefficient_term_needs_no_rebuild(baseline_model: Model) -> None: diff --git a/test/test_persistent_snapshot_diff.py b/test/test_persistent_snapshot_diff.py index f9135cb76..8ea2231f5 100644 --- a/test/test_persistent_snapshot_diff.py +++ b/test/test_persistent_snapshot_diff.py @@ -148,7 +148,8 @@ def test_coef_sparsity_change(baseline: Model) -> None: x = baseline.variables["x"] baseline.constraints["c2"].lhs = 2 * x.sum() diff = ModelDiff.from_snapshot(snap, baseline) - assert diff is RebuildReason.SPARSITY + assert isinstance(diff, ModelDiff) + np.testing.assert_array_equal(diff.con_coef_vals, [2, 2, 2, 0, 0]) def test_deep_copy_invariant(baseline: Model) -> None: @@ -431,7 +432,7 @@ def build(aggregate: bool) -> Model: assert diff is RebuildReason.STRUCTURAL_LABELS -def test_indices_change_triggers_sparsity() -> None: +def test_indices_change_clears_old_coefficient() -> None: def build(on: int) -> Model: m = Model() x = m.add_variables(0, 10, coords=[range(2)], name="x") @@ -439,7 +440,11 @@ def build(on: int) -> Model: m.add_objective(x.sum()) return m - assert ModelDiff.from_models(build(0), build(1)) is RebuildReason.SPARSITY + diff = ModelDiff.from_models(build(0), build(1)) + assert isinstance(diff, ModelDiff) + np.testing.assert_array_equal(diff.con_coef_rows, [0, 0]) + np.testing.assert_array_equal(diff.con_coef_cols, [0, 1]) + np.testing.assert_array_equal(diff.con_coef_vals, [0, 1]) def test_sign_only_mutation(baseline: Model) -> None: diff --git a/test/test_persistent_solver_extras.py b/test/test_persistent_solver_extras.py index 3fc28d18d..69ff24949 100644 --- a/test/test_persistent_solver_extras.py +++ b/test/test_persistent_solver_extras.py @@ -128,7 +128,7 @@ def test_cross_model_scenario_sweep(solver_name: str) -> None: @pytest.mark.parametrize("solver_name", SOLVER_PARAMS) -def test_cross_model_sparsity_change_rebuilds(solver_name: str) -> None: +def test_cross_model_sparsity_change_updates_in_place(solver_name: str) -> None: def build(include_y_in_c1: bool) -> Model: m = Model() x = m.add_variables(0, 10, coords=[range(3)], name="x") @@ -148,12 +148,16 @@ def build(include_y_in_c1: bool) -> Model: m2 = build(include_y_in_c1=False) s.solve(m2, assign=True) - assert s._rebuilds == 1 - assert s._last_rebuild_reason in { - RebuildReason.SPARSITY, - RebuildReason.STRUCTURAL_LABELS, - RebuildReason.STRUCTURAL_CONTAINERS, - } + assert s._rebuilds == 0 + assert s._in_place_updates == 1 + fresh = build(include_y_in_c1=False) + fresh_solver = _built(solver_name, fresh) + try: + fresh_solver.solve(assign=True) + np.testing.assert_allclose(_obj(m2), _obj(fresh)) + finally: + fresh_solver.close() + s.close() @pytest.mark.parametrize("solver_name", SOLVER_PARAMS) diff --git a/test/test_persistent_zero_crossing.py b/test/test_persistent_zero_crossing.py new file mode 100644 index 000000000..7d1eee881 --- /dev/null +++ b/test/test_persistent_zero_crossing.py @@ -0,0 +1,92 @@ +""" +Persistent coefficient updates must handle zero crossings. + +Failure modes: crossing zero causes an unnecessary rebuild, a removed term +remains active in the solver, a restored term is not inserted, same-model +updates miss a dirty coefficient, or a fully zero row retains old coefficients. +""" + +from __future__ import annotations + +import numpy as np +import pytest + +from linopy import Model +from linopy.solvers import Solver + + +def _activation_model(ratio: float, frozen: bool) -> Model: + model = Model() + activation = model.add_variables(0, 1, name="activation") + reservation = model.add_variables(0, 1, name="reservation") + model.add_constraints( + activation <= ratio * reservation, name="limit", freeze=frozen + ) + model.add_objective(-activation) + return model + + +@pytest.mark.parametrize("ratios", [(1.0, 0.0, 1.0), (0.0, 1.0, 0.0)]) +@pytest.mark.parametrize( + "frozen,same_model", [(False, False), (True, False), (False, True)] +) +def test_zero_crossing_updates_match_fresh_solve( + ratios: tuple[float, ...], + frozen: bool, + same_model: bool, +) -> None: + solver = Solver.from_name( + "highs", + model=_activation_model(ratios[0], frozen), + io_api="direct", + track_updates=True, + set_names=False, + ) + solver.solve(assign=True) + for ratio in ratios[1:]: + if same_model: + updated = solver.model + assert updated is not None + activation = updated.variables["activation"] + reservation = updated.variables["reservation"] + updated.constraints["limit"].update(lhs=activation - ratio * reservation) + else: + updated = _activation_model(ratio, frozen) + solver.update(updated) + solver.solve(assign=True) + fresh = _activation_model(ratio, frozen) + fresh.solve(solver_name="highs", io_api="direct") + updated_value = updated.objective.value + fresh_value = fresh.objective.value + assert updated_value is not None + assert fresh_value is not None + np.testing.assert_allclose(updated_value, fresh_value) + np.testing.assert_allclose(updated_value, -ratio) + assert solver._rebuilds == 0 + solver.close() + + +@pytest.mark.parametrize("frozen", [False, True]) +def test_entire_row_zero_crossing_updates_in_place(frozen: bool) -> None: + def build(ratio: float) -> Model: + model = Model() + x = model.add_variables(0, 1, name="x") + model.add_constraints(ratio * x <= 0, name="limit", freeze=frozen) + model.add_objective(-x) + return model + + solver = Solver.from_name( + "highs", + model=build(1.0), + io_api="direct", + track_updates=True, + ) + solver.solve(assign=True) + for ratio in (0.0, 1.0): + model = build(ratio) + solver.solve(model, assign=True) + value = model.objective.value + assert value is not None + np.testing.assert_allclose(value, ratio - 1) + assert solver._rebuilds == 0 + solver.close()