Handle zero-crossing coefficients in persistent updates

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:

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.
 

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:

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,15 +432,19 @@ 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")
         m.add_constraints(x.loc[on] >= 1, name="c1")
         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:

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)

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()