From 26a783cf57acb79e6fe509f7b42a3a8dbdaf1a85 Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Mon, 14 Mar 2022 15:13:56 -0400 Subject: [PATCH 1/2] partial update tests --- src/osqp/tests/update_matrices_test.py | 86 +++++++++++++++++++++++++- src/osqp/tests/utils.py | 27 +++++++- 2 files changed, 110 insertions(+), 3 deletions(-) diff --git a/src/osqp/tests/update_matrices_test.py b/src/osqp/tests/update_matrices_test.py index dae0790..f0a26fd 100644 --- a/src/osqp/tests/update_matrices_test.py +++ b/src/osqp/tests/update_matrices_test.py @@ -1,6 +1,8 @@ # Test osqp python module +import pytest + import osqp -from osqp.tests.utils import load_high_accuracy, rel_tol, abs_tol, decimal_tol +from osqp.tests.utils import load_high_accuracy, rel_tol, abs_tol, decimal_tol, Random import numpy as np import scipy as sp from scipy import sparse @@ -20,7 +22,7 @@ def setUp(self): self.m = 8 p = 0.7 - Pt = sparse.random(self.n, self.n, density=p) + self.Pt = Pt = sparse.random(self.n, self.n, density=p) Pt_new = Pt.copy() Pt_new.data += 0.1 * np.random.randn(Pt.nnz) @@ -67,6 +69,86 @@ def test_update_P(self): nptest.assert_almost_equal( res.info.obj_val, obj_sol, decimal=decimal_tol) + def test_update_P_partial(self): + with Random(4234): + n_changed = np.random.randint(self.P_triu.nnz) + changed_data = np.random.random(n_changed) + changed_indices = np.random.choice(np.arange(self.P_triu.nnz), n_changed) + changed_P_triu_data = self.P_triu.data.copy() + changed_P_triu_data[changed_indices] = changed_data + changed_P_triu = np.array(sparse.coo_matrix((changed_P_triu_data, (self.P_triu.row, self.P_triu.col))).todense()) + changed_P = np.triu(changed_P_triu, 1) + np.tril(changed_P_triu.T) + if not np.all(np.linalg.eigvals(changed_P) > 0): + pytest.skip("Perturbed P not positive semi-definite") + + self.model.update(Px=changed_data, Px_idx=changed_indices) + res1 = self.model.solve() + + # The results we obtain should be the same as if we were solving a new problem with the new P + model = osqp.OSQP() + model.setup(P=changed_P, q=self.q, A=self.A, l=self.l, u=self.u, + **self.opts) + res2 = model.solve() + + 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): + with Random(60023): + n_changed = np.random.randint(self.A.nnz) + changed_data = np.random.random(n_changed) + changed_indices = np.random.choice(np.arange(self.A.nnz), n_changed) + 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)) + + self.model.update(Ax=changed_data, Ax_idx=changed_indices) + res1 = self.model.solve() + + # The results we obtain should be the same as if we were solving a new problem with the new A + model = osqp.OSQP() + model.setup(P=self.P, q=self.q, A=changed_A, l=self.l, u=self.u, + **self.opts) + res2 = model.solve() + + 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): + with Random(54355): + n_P_changed = np.random.randint(self.P_triu.nnz) + _changed_P_data = np.random.random(n_P_changed) + changed_P_indices = np.random.choice(np.arange(self.P_triu.nnz), n_P_changed) + n_A_changed = np.random.randint(self.A.nnz) + _changed_A_data = np.random.random(n_A_changed) + changed_A_indices = np.random.choice(np.arange(self.A.nnz), n_A_changed) + + changed_P_triu_data = self.P_triu.data.copy() + changed_P_triu_data[changed_P_indices] = _changed_P_data + changed_P_triu = np.array(sparse.coo_matrix((changed_P_triu_data, (self.P_triu.row, self.P_triu.col))).todense()) + changed_P = np.triu(changed_P_triu, 1) + np.tril(changed_P_triu.T) + if not np.all(np.linalg.eigvals(changed_P) > 0): + pytest.skip("Perturbed P not positive semi-definite") + + 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)) + + self.model.update(Px=_changed_P_data, Px_idx=changed_P_indices, Ax=_changed_A_data, Ax_idx=changed_A_indices) + res1 = self.model.solve() + + # The results we obtain should be the same as if we were solving a new problem with the new P/A + model = osqp.OSQP() + model.setup(P=changed_P, q=self.q, A=changed_A, l=self.l, u=self.u, + **self.opts) + res2 = model.solve() + + 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 e5e2f9a..a6b06fc 100644 --- a/src/osqp/tests/utils.py +++ b/src/osqp/tests/utils.py @@ -10,4 +10,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()) From 889fae750b19f1afd8675e31a35ecc130830002f Mon Sep 17 00:00:00 2001 From: Vineet Bansal Date: Wed, 16 Mar 2022 15:54:47 -0400 Subject: [PATCH 2/2] more robust partial update tests --- CMakeLists.txt | 2 +- src/osqp/interface.py | 16 ++- src/osqp/tests/update_matrices_test.py | 190 +++++++++++++++---------- 3 files changed, 128 insertions(+), 80 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 0701bb4..b187683 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 a800b2a..1d55874 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 d7a0f10..46843a9 100644 --- a/src/osqp/tests/update_matrices_test.py +++ b/src/osqp/tests/update_matrices_test.py @@ -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') @@ -72,86 +72,132 @@ def test_update_P(self): def test_update_P_partial(self): - with Random(4234): - n_changed = np.random.randint(self.P_triu.nnz) - changed_data = np.random.random(n_changed) - changed_indices = np.random.choice(np.arange(self.P_triu.nnz), n_changed) - changed_P_triu_data = self.P_triu.data.copy() - changed_P_triu_data[changed_indices] = changed_data - changed_P_triu = np.array(sparse.coo_matrix((changed_P_triu_data, (self.P_triu.row, self.P_triu.col))).todense()) - changed_P = sparse.csc_matrix(np.triu(changed_P_triu, 1) + np.tril(changed_P_triu.T)) - if not np.all(np.linalg.eigvals(changed_P.todense()) > 0): - pytest.skip("Perturbed P not positive semi-definite") - - self.model.update(Px=changed_data, Px_idx=changed_indices) - res1 = self.model.solve() - - # The results we obtain should be the same as if we were solving a new problem with the new P - model = osqp.OSQP() - model.setup(P=changed_P, q=self.q, A=self.A, l=self.l, u=self.u, - **self.opts) - res2 = model.solve() - 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) + # 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): - with Random(60023): - n_changed = np.random.randint(self.A.nnz) - changed_data = np.random.random(n_changed) - changed_indices = np.random.choice(np.arange(self.A.nnz), n_changed) - 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)) - - self.model.update(Ax=changed_data, Ax_idx=changed_indices) - res1 = self.model.solve() - - # The results we obtain should be the same as if we were solving a new problem with the new A - model = osqp.OSQP() - model.setup(P=self.P, q=self.q, A=changed_A, l=self.l, u=self.u, - **self.opts) - res2 = model.solve() - 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) + 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): - with Random(54355): - n_P_changed = np.random.randint(self.P_triu.nnz) - _changed_P_data = np.random.random(n_P_changed) - changed_P_indices = np.random.choice(np.arange(self.P_triu.nnz), n_P_changed) - n_A_changed = np.random.randint(self.A.nnz) - _changed_A_data = np.random.random(n_A_changed) - changed_A_indices = np.random.choice(np.arange(self.A.nnz), n_A_changed) - - changed_P_triu_data = self.P_triu.data.copy() - changed_P_triu_data[changed_P_indices] = _changed_P_data - changed_P_triu = np.array(sparse.coo_matrix((changed_P_triu_data, (self.P_triu.row, self.P_triu.col))).todense()) - changed_P = sparse.csc_matrix(np.triu(changed_P_triu, 1) + np.tril(changed_P_triu.T)) - if not np.all(np.linalg.eigvals(changed_P.todense()) > 0): - pytest.skip("Perturbed P not positive semi-definite") - - 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)) - - self.model.update(Px=_changed_P_data, Px_idx=changed_P_indices, Ax=_changed_A_data, Ax_idx=changed_A_indices) - res1 = self.model.solve() - - # The results we obtain should be the same as if we were solving a new problem with the new P/A - model = osqp.OSQP() - model.setup(P=changed_P, q=self.q, A=changed_A, l=self.l, u=self.u, - **self.opts) - res2 = model.solve() - 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) + 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):