Test correctness of GPU paths (closes #40)

test/gen_random_cone_prob.py

@@ -1,4 +1,3 @@
-import scs
 import numpy as np
 from scipy import sparse
 
@@ -159,7 +158,7 @@ def proj_sdp(z, n):
   a = np.dot(v, np.dot(np.diag(w), np.transpose(v)))
   a[didx] = a[didx] / np.sqrt(2.)
   z = a[tidx]
-  return z
+  return np.real(z)
 
 
 def proj_pow(v, a):
@@ -221,7 +220,8 @@ def project_exp_bisection(v):
   s = v[1]
   t = v[2]
   # v in cl(Kexp)
-  if (s * np.exp(r / s) <= t and s > 0) or (r <= 0 and s == 0 and t >= 0):
+  if (s > 0 and t > 0 and r <= s * np.log(t / s)) or (r <= 0 and s == 0 and
+                                                      t >= 0):
     return v
   # -v in Kexp^*
   if (-r < 0 and r * np.exp(s / r) <= -np.exp(1) * t) or (-r == 0 and

test/test_solve_random_cone_prob.py

@@ -8,6 +8,7 @@
 #  Uses scs to solve a random cone problem  #
 #############################################
 
+
 def import_error(msg):
   print()
   print('## IMPORT ERROR:' + msg)
@@ -34,33 +35,50 @@ def import_error(msg):
     'f': 10,
     'l': 15,
     'q': [5, 10, 0, 1],
-    's': [3, 4, 0, 0, 1],
+    's': [3, 4, 0, 0, 1, 10],
     'ep': 10,
     'ed': 10,
     'p': [-0.25, 0.5, 0.75, -0.33]
 }
 m = tools.get_scs_cone_dims(K)
-params = {'normalize': True, 'scale': 5, 'cg_rate': 2, 'verbose': False}
+params = {'normalize': True, 'scale': 5, 'cg_rate': 2, 'verbose': True}
+
 
-@pytest.mark.parametrize("use_indirect,gpu", flags)
+@pytest.mark.parametrize('use_indirect,gpu', flags)
 def test_solve_feasible(use_indirect, gpu):
-  data, p_star = tools.gen_feasible(K, n=m // 3, density=0.01)
+  data, p_star = tools.gen_feasible(K, n=m // 3, density=0.1)
 
   sol = scs.solve(data, K, use_indirect=use_indirect, gpu=gpu, **params)
   x = sol['x']
   y = sol['y']
-  print('p*  = ', p_star)
-  print('pri error = ', (np.dot(data['c'], x) - p_star) / p_star)
-  print('dual error = ', (-np.dot(data['b'], y) - p_star) / p_star)
+  s = sol['s']
+  np.testing.assert_almost_equal(np.dot(data['c'], x), p_star, decimal=3)
+  np.testing.assert_almost_equal(np.dot(data['c'], x), p_star, decimal=3)
+  np.testing.assert_array_less(
+      np.linalg.norm(data['A'] @ x - data['b'] + s), 1e-3)
+  np.testing.assert_array_less(
+      np.linalg.norm(data['A'].T @ y + data['c']), 1e-3)
+  np.testing.assert_almost_equal(s.T @ y, 0.)
+  np.testing.assert_almost_equal(s, tools.proj_cone(s, K), decimal=4)
+  np.testing.assert_almost_equal(y, tools.proj_dual_cone(y, K), decimal=4)
 
 
-@pytest.mark.parametrize("use_indirect,gpu", flags)
+@pytest.mark.parametrize('use_indirect,gpu', flags)
 def test_solve_infeasible(use_indirect, gpu):
   data = tools.gen_infeasible(K, n=m // 3)
   sol = scs.solve(data, K, use_indirect=use_indirect, gpu=gpu, **params)
+  y = sol['y']
+  np.testing.assert_array_less(np.linalg.norm(data['A'].T @ y), 1e-3)
+  np.testing.assert_array_less(data['b'].T @ y, -0.1)
+  np.testing.assert_almost_equal(y, tools.proj_dual_cone(y, K), decimal=4)
 
 
-@pytest.mark.parametrize("use_indirect,gpu", flags)
+@pytest.mark.parametrize('use_indirect,gpu', flags)
 def test_solve_unbounded(use_indirect, gpu):
   data = tools.gen_unbounded(K, n=m // 3)
   sol = scs.solve(data, K, use_indirect=use_indirect, gpu=gpu, **params)
+  x = sol['x']
+  s = sol['s']
+  np.testing.assert_array_less(np.linalg.norm(data['A'] @ x + s), 1e-3)
+  np.testing.assert_array_less(data['c'].T @ x, -0.1)
+  np.testing.assert_almost_equal(s, tools.proj_cone(s, K), decimal=4)