diff --git a/CMakeLists.txt b/CMakeLists.txt index 0701bb4a..b1876839 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/src/osqp/interface.py b/src/osqp/interface.py index a800b2ac..1d558745 100644 --- a/src/osqp/interface.py +++ b/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): """ diff --git a/src/osqp/tests/update_matrices_test.py b/src/osqp/tests/update_matrices_test.py index 76aa890b..46843a91 100644 --- a/src/osqp/tests/update_matrices_test.py +++ b/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 diff --git a/src/osqp/tests/utils.py b/src/osqp/tests/utils.py index c787f5a8..747862dc 100644 --- a/src/osqp/tests/utils.py +++ b/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'] \ No newline at end of file + 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())