Vb/pybind11 partial updates

CMakeLists.txt

@@ -7,7 +7,7 @@ list(APPEND CMAKE_MESSAGE_INDENT "  ")
 FetchContent_Declare(
   osqp
   GIT_REPOSITORY https://github.com/osqp/osqp.git
-  GIT_TAG vb/develop1-fix)
+  GIT_TAG develop-1.0)
 
 list(POP_BACK CMAKE_MESSAGE_INDENT)
 FetchContent_MakeAvailable(osqp)

src/osqp/interface.py

@@ -169,15 +169,15 @@ def update(self, q=None, l=None, u=None,
             self._model.update_bounds(l, u)
 
         # update matrix P
-        if Px is not None and Ax is None:
+        if (Px is not None and Px.size > 0) and (Ax is None or Ax.size == 0):
             self._model.update_P(Px, Px_idx, len(Px))
 
         # update matrix A
-        if Ax is not None and Px is None:
+        if (Ax is not None and Ax.size > 0) and (Px is None or Px.size == 0):
             self._model.update_A(Ax, Ax_idx, len(Ax))
 
         # update matrices P and A
-        if Px is not None and Ax is not None:
+        if Px is not None and Px.size > 0 and Ax is not None and Ax.size > 0:
             self._model.update_P_A(Px, Px_idx, len(Px), Ax, Ax_idx, len(Ax))
 
 
@@ -209,10 +209,12 @@ def update(self, q=None, l=None, u=None,
 
         # delete results from self._derivative_cache to prohibit
         # taking the derivative of unsolved problems
-        if "results" in self._derivative_cache.keys():
-            del self._derivative_cache["results"]
-            del self._derivative_cache["M"]
-            del self._derivative_cache["solver"]
+        if 'results' in self._derivative_cache.keys():
+            del self._derivative_cache['results']
+            if 'M' in self._derivative_cache:
+                del self._derivative_cache['M']
+            if 'solver' in self._derivative_cache:
+                del self._derivative_cache['solver']
 
     def update_settings(self, **kwargs):
         """

src/osqp/tests/update_matrices_test.py

@@ -1,6 +1,6 @@
 from types import SimpleNamespace
 import osqp
-from osqp.tests.utils import load_high_accuracy, rel_tol, abs_tol, decimal_tol, SOLVER_TYPES
+from osqp.tests.utils import load_high_accuracy, rel_tol, abs_tol, decimal_tol, Random, SOLVER_TYPES
 import numpy as np
 import scipy as sp
 from scipy import sparse
@@ -25,7 +25,7 @@ def self(request):
 
     self.P = (Pt.T.dot(Pt) + sparse.eye(self.n)).tocsc()
     self.P_new = (Pt_new.T.dot(Pt_new) + sparse.eye(self.n)).tocsc()
-    self.P_triu = sparse.triu(self.P)
+    self.P_triu = sparse.triu(self.P).tocsc()
     self.P_triu_new = sparse.triu(self.P_new)
     self.q = np.random.randn(self.n)
     self.A = sparse.random(self.m, self.n, density=p, format='csc')
@@ -71,6 +71,135 @@ def test_update_P(self):
         res.info.obj_val, obj_sol, decimal=decimal_tol)
 
 
+def test_update_P_partial(self):
+
+    # Run a bunch of tests in a loop - for 1000 iterations, we expect ~200 to yield a converged solution, for which we
+    # test that the partial update solution matches with the case where we supply the perturbed P matrix to begin with.
+    for i in range(1000):
+        model1 = osqp.OSQP()
+        # Note: the supplied P to the model will change between iterations through model.update(..), so it's important
+        # to initialize with a copy of self.P so we're starting afresh
+        model1.setup(P=self.P.copy(), q=self.q, A=self.A, l=self.l, u=self.u, **self.opts)
+
+        with Random(i):
+            n_changed = np.random.randint(self.P_triu.nnz)
+            changed_data = np.random.random(n_changed)
+            changed_indices = np.sort(np.random.choice(np.arange(self.P_triu.nnz), n_changed, replace=False))
+        changed_P_triu_data = self.P_triu.data.copy()
+        changed_P_triu_data[changed_indices] = changed_data
+        changed_P_triu_csc = sparse.csc_matrix((changed_P_triu_data, self.P_triu.indices, self.P_triu.indptr))
+        changed_P_triu_dense = changed_P_triu_csc.todense()
+        changed_P_dense = changed_P_triu_dense + np.tril(changed_P_triu_dense.T, -1)
+
+        if not np.all(np.linalg.eigvals(changed_P_dense) >= 0):  # not positive semi-definite?
+            continue
+
+        model1.update(Px=changed_data, Px_idx=changed_indices)
+        # The results we obtain should be the same as if we were solving a new problem with the new P
+        model2 = osqp.OSQP()
+
+        try:
+            res1 = model1.solve()
+            model2.setup(P=changed_P_triu_csc, q=self.q, A=self.A, l=self.l, u=self.u, **self.opts)
+            res2 = model2.solve()
+        except ValueError:  # ignore any setup errors because of non-convexity etc.
+            continue
+
+        # We only care about solved instances, and for which the solution/objective is not too close to 0.
+        # A zero objective value hints that the solution is not unique, in which case the duals may not be the same.
+        if res1.info.status != 'solved' or res2.info.status != 'solved' or np.max(np.abs(res1.x) < 1e-6):
+            continue
+
+        assert np.allclose(res1.x, res2.x, atol=abs_tol, rtol=rel_tol)
+        assert np.allclose(res1.y, res2.y, atol=abs_tol, rtol=rel_tol)
+        assert np.allclose(res1.info.obj_val, res2.info.obj_val, atol=abs_tol, rtol=rel_tol)
+
+
+def test_update_A_partial(self):
+
+    for i in range(1000):
+
+        model1 = osqp.OSQP()
+        # Note: the supplied A to the model will change between iterations through model.update(..), so it's important
+        # to initialize with a copy of self.A so we're starting afresh
+        model1.setup(P=self.P.copy(), q=self.q, A=self.A.copy(), l=self.l, u=self.u, **self.opts)
+
+        with Random(i):
+            n_changed = np.random.randint(self.A.nnz)
+            changed_data = np.random.random(n_changed)
+            changed_indices = np.sort(np.random.choice(np.arange(self.A.nnz), n_changed, replace=False))
+        changed_A_data = self.A.data.copy()
+        changed_A_data[changed_indices] = changed_data
+        changed_A = sparse.csc_matrix((changed_A_data, self.A.indices, self.A.indptr))
+
+        model1.update(Ax=changed_data, Ax_idx=changed_indices)
+        res1 = model1.solve()
+
+        # The results we obtain should be the same as if we were solving a new problem with the new A
+        model2 = osqp.OSQP()
+        model2.setup(P=self.P.copy(), q=self.q, A=changed_A, l=self.l, u=self.u, **self.opts)
+        res2 = model2.solve()
+
+        # We only care about solved instances, and for which the solution/objective is not too close to 0.
+        # A zero objective value hints that the solution is not unique, in which case the duals may not be the same.
+        if res1.info.status != 'solved' or res2.info.status != 'solved' or np.max(np.abs(res1.x) < 1e-6):
+            continue
+
+        assert np.allclose(res1.x, res2.x, atol=abs_tol, rtol=rel_tol)
+        assert np.allclose(res1.y, res2.y, atol=abs_tol, rtol=rel_tol)
+        assert np.allclose(res1.info.obj_val, res2.info.obj_val, atol=abs_tol, rtol=rel_tol)
+
+
+def test_update_P_A_partial(self):
+
+    for i in range(41, 1000):
+        model1 = osqp.OSQP()
+        # Note: the supplied P/A to the model will change between iterations through model.update(..), so it's important
+        # to initialize with a copy of self.P/self.A so we're starting afresh
+        model1.setup(P=self.P.copy(), q=self.q, A=self.A.copy(), l=self.l, u=self.u, **self.opts)
+
+        with Random(i):
+            n_P_changed = np.random.randint(self.P_triu.nnz)
+            _changed_P_data = np.random.random(n_P_changed)
+            changed_P_indices = np.sort(np.random.choice(np.arange(self.P_triu.nnz), n_P_changed, replace=False))
+            n_A_changed = np.random.randint(self.A.nnz)
+            _changed_A_data = np.random.random(n_A_changed)
+            changed_A_indices = np.sort(np.random.choice(np.arange(self.A.nnz), n_A_changed, replace=False))
+
+        changed_P_triu_data = self.P_triu.data.copy()
+        changed_P_triu_data[changed_P_indices] = _changed_P_data
+        changed_P_triu_csc = sparse.csc_matrix((changed_P_triu_data, self.P_triu.indices, self.P_triu.indptr))
+        changed_P_triu_dense = changed_P_triu_csc.todense()
+        changed_P_dense = changed_P_triu_dense + np.tril(changed_P_triu_dense.T, -1)
+
+        changed_A_data = self.A.data.copy()
+        changed_A_data[changed_A_indices] = _changed_A_data
+        changed_A = sparse.csc_matrix((changed_A_data, self.A.indices, self.A.indptr))
+
+        if not np.all(np.linalg.eigvals(changed_P_dense) >= 0):  # not positive semi-definite?
+            continue
+
+        model1.update(Px=_changed_P_data, Px_idx=changed_P_indices, Ax=_changed_A_data, Ax_idx=changed_A_indices)
+        # The results we obtain should be the same as if we were solving a new problem with the new P/A
+        model2 = osqp.OSQP()
+
+        try:
+            res1 = model1.solve()
+            model2.setup(P=changed_P_triu_csc, q=self.q, A=changed_A, l=self.l, u=self.u, **self.opts)
+            res2 = model2.solve()
+        except ValueError:  # ignore any setup errors because of non-convexity etc.
+            continue
+
+        # We only care about solved instances, and for which the solution/objective is not too close to 0.
+        # A zero objective value hints that the solution is not unique, in which case the duals may not be the same.
+        if res1.info.status != 'solved' or res2.info.status != 'solved' or np.max(np.abs(res1.x) < 1e-6):
+            continue
+
+        assert np.allclose(res1.x, res2.x, atol=abs_tol, rtol=rel_tol)
+        assert np.allclose(res1.y, res2.y, atol=abs_tol, rtol=rel_tol)
+        assert np.allclose(res1.info.obj_val, res2.info.obj_val, atol=abs_tol, rtol=rel_tol)
+
+
 def test_update_P_allind(self):
     # Update matrix P
     Px = self.P_triu_new.data

src/osqp/tests/utils.py

@@ -21,4 +21,29 @@
 def load_high_accuracy(test_name):
     npz = os.path.join(os.path.dirname(__file__), 'solutions', f'{test_name}.npz')
     npzfile = np.load(npz)
-    return npzfile['x_val'], npzfile['y_val'], npzfile['obj']
+    return npzfile['x_val'], npzfile['y_val'], npzfile['obj']
+
+
+# A list of random states, used as a stack
+random_states = []
+
+
+class Random:
+    """
+    A context manager that pushes a random seed to the stack for reproducible results,
+    and pops it on exit.
+    """
+
+    def __init__(self, seed=None):
+        self.seed = seed
+
+    def __enter__(self):
+        if self.seed is not None:
+            # Push current state on stack
+            random_states.append(np.random.get_state())
+            new_state = np.random.RandomState(self.seed)
+            np.random.set_state(new_state.get_state())
+
+    def __exit__(self, *args):
+        if self.seed is not None:
+            np.random.set_state(random_states.pop())