From f131da703ba4022f83c21fae229bf3bac1537cba Mon Sep 17 00:00:00 2001
From: goulart-paul <paul.goulart@eng.ox.ac.uk>
Date: Tue, 7 Jan 2025 13:00:01 +0000
Subject: [PATCH] pytest for basic qp/sdp

---
 python/clarabel/tests/test_basic_qp.py  | 86 +++++++++++++++++++++++++
 python/clarabel/tests/test_basic_sdp.py | 81 +++++++++++++++++++++++
 tests/basic_sdp.rs                      | 43 ++++++-------
 3 files changed, 185 insertions(+), 25 deletions(-)
 create mode 100644 python/clarabel/tests/test_basic_qp.py
 create mode 100644 python/clarabel/tests/test_basic_sdp.py

diff --git a/python/clarabel/tests/test_basic_qp.py b/python/clarabel/tests/test_basic_qp.py
new file mode 100644
index 00000000..daccc3b5
--- /dev/null
+++ b/python/clarabel/tests/test_basic_qp.py
@@ -0,0 +1,86 @@
+import clarabel
+import pytest
+import numpy as np
+from scipy import sparse
+
+
+@pytest.fixture
+def basic_qp_data():
+
+    P = sparse.csc_matrix([[4., 1.], [1., 2.]])
+    P = sparse.triu(P).tocsc()
+
+    A = sparse.csc_matrix(
+        [[-1., -1.],
+         [-1.,  0.],
+         [0.,  -1.],
+         [1.,   1.],
+         [1.,   0.],
+         [0.,   1.]])
+
+    q = np.array([1., 1.])
+    b = np.array([-1., 0., 0., 1., 0.7, 0.7])
+
+    cones = [clarabel.NonnegativeConeT(3), clarabel.NonnegativeConeT(3)]
+    settings = clarabel.DefaultSettings()
+    return P, q, A, b, cones, settings
+
+
+@pytest.fixture
+def basic_qp_data_dual_inf():
+
+    P = sparse.csc_matrix([[1., 1.], [1., 1.]])
+    P = sparse.triu(P).tocsc()
+
+    A = sparse.csc_matrix(
+        [[1.,   1.],
+         [1.,   0.]])
+
+    q = np.array([1., -1.])
+    b = np.array([1., 1.])
+
+    cones = [clarabel.NonnegativeConeT(2)]
+    settings = clarabel.DefaultSettings()
+    return P, q, A, b, cones, settings
+
+
+def test_qp_feasible(basic_qp_data):
+
+    P, q, A, b, cones, settings = basic_qp_data
+
+    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
+    solution = solver.solve()
+
+    refsol = np.array([0.3, 0.7])
+    refobj = 1.8800000298331538
+
+    assert solution.status == clarabel.SolverStatus.Solved
+    assert np.allclose(solution.x, refsol)
+    assert np.allclose(solution.obj_val, refobj)
+    assert np.allclose(solution.obj_val_dual, refobj)
+
+
+def test_qp_primal_infeasible(basic_qp_data):
+
+    P, q, A, b, cones, settings = basic_qp_data
+    b[0] = -1.
+    b[3] = -1.
+
+    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
+    solution = solver.solve()
+
+    assert solution.status == clarabel.SolverStatus.PrimalInfeasible
+    assert np.isnan(solution.obj_val)
+    assert np.isnan(solution.obj_val_dual)
+
+
+def test_qp_dual_infeasible(basic_qp_data_dual_inf):
+
+    P, q, A, b, cones, settings = basic_qp_data_dual_inf
+
+    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
+    solution = solver.solve()
+
+    assert solution.status == clarabel.SolverStatus.DualInfeasible
+    assert np.isnan(solution.obj_val)
+    assert np.isnan(solution.obj_val_dual)
diff --git a/python/clarabel/tests/test_basic_sdp.py b/python/clarabel/tests/test_basic_sdp.py
new file mode 100644
index 00000000..74308c69
--- /dev/null
+++ b/python/clarabel/tests/test_basic_sdp.py
@@ -0,0 +1,81 @@
+import clarabel
+import pytest
+import numpy as np
+from scipy import sparse
+from scipy.sparse import vstack
+
+
+@pytest.fixture
+def basic_sdp_data():
+
+    P = sparse.eye(6).tocsc()
+    A = sparse.eye(6).tocsc()
+
+    q = np.zeros(6)
+    b = np.array([-3., 1., 4., 1., 2., 5.])
+
+    cones = [clarabel.PSDTriangleConeT(3)]
+    settings = clarabel.DefaultSettings()
+    return P, q, A, b, cones, settings
+
+
+@pytest.fixture
+def basic_sdp_solution():
+
+    refsol = np.array([
+        -3.0729833267361095,
+        0.3696004167288786,
+        -0.022226685581313674,
+        0.31441213129613066,
+        -0.026739700851545107,
+        -0.016084530571308823,
+    ])
+    refobj = 4.840076866013861
+
+    return refsol, refobj
+
+
+def test_sdp_feasible(basic_sdp_data, basic_sdp_solution):
+
+    P, q, A, b, cones, settings = basic_sdp_data
+    refsol, refobj = basic_sdp_solution
+
+    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
+    solution = solver.solve()
+
+    assert solution.status == clarabel.SolverStatus.Solved
+    assert np.allclose(solution.x, refsol)
+    assert np.allclose(solution.obj_val, refobj)
+    assert np.allclose(solution.obj_val_dual, refobj)
+
+
+def test_sdp_empty_cone(basic_sdp_data, basic_sdp_solution):
+
+    P, q, A, b, cones, settings = basic_sdp_data
+    refsol, refobj = basic_sdp_solution
+
+    cones = np.append(cones, clarabel.PSDTriangleConeT(0))
+
+    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
+    solution = solver.solve()
+
+    assert solution.status == clarabel.SolverStatus.Solved
+    assert np.allclose(solution.x, refsol)
+    assert np.allclose(solution.obj_val, refobj)
+    assert np.allclose(solution.obj_val_dual, refobj)
+
+
+def test_sdp_primal_infeasible(basic_sdp_data):
+
+    P, q, A, b, cones, settings = basic_sdp_data
+
+    A = vstack((A, -A))
+    b = np.pad(b, (0, len(b)))
+    cones = np.concatenate((cones, cones))
+
+    solver = clarabel.DefaultSolver(P, q, A, b, cones, settings)
+    solution = solver.solve()
+
+    assert solution.status == clarabel.SolverStatus.PrimalInfeasible
+    assert np.isnan(solution.obj_val)
+    assert np.isnan(solution.obj_val_dual)
diff --git a/tests/basic_sdp.rs b/tests/basic_sdp.rs
index 815ce3ab..683ce77a 100644
--- a/tests/basic_sdp.rs
+++ b/tests/basic_sdp.rs
@@ -26,9 +26,24 @@ fn basic_sdp_data() -> (
     (P, c, A, b, cones)
 }
 
+fn basic_sdp_solution() -> (Vec<f64>, f64) {
+    let refsol = vec![
+        -3.0729833267361095,
+        0.3696004167288786,
+        -0.022226685581313674,
+        0.31441213129613066,
+        -0.026739700851545107,
+        -0.016084530571308823,
+    ];
+    let refobj = 4.840076866013861;
+
+    (refsol, refobj)
+}
+
 #[test]
 fn test_sdp_feasible() {
     let (P, c, A, b, cones) = basic_sdp_data();
+    let (refsol, refobj) = basic_sdp_solution();
 
     let settings = DefaultSettings::default();
 
@@ -37,25 +52,15 @@ fn test_sdp_feasible() {
     solver.solve();
 
     assert_eq!(solver.solution.status, SolverStatus::Solved);
-
-    let refsol = vec![
-        -3.0729833267361095,
-        0.3696004167288786,
-        -0.022226685581313674,
-        0.31441213129613066,
-        -0.026739700851545107,
-        -0.016084530571308823,
-    ];
-
     assert!(solver.solution.x.dist(&refsol) <= 1e-6);
-
-    let refobj = 4.840076866013861;
     assert!(f64::abs(solver.info.cost_primal - refobj) <= 1e-6);
 }
 
 #[test]
-fn empty_sdp_cone() {
+fn test_sdp_empty_cone() {
     let (P, c, A, b, mut cones) = basic_sdp_data();
+    let (refsol, refobj) = basic_sdp_solution();
+
     cones.append(&mut vec![PSDTriangleConeT(0)]);
 
     let settings = DefaultSettings::default();
@@ -65,19 +70,7 @@ fn empty_sdp_cone() {
     solver.solve();
 
     assert_eq!(solver.solution.status, SolverStatus::Solved);
-
-    let refsol = vec![
-        -3.0729833267361095,
-        0.3696004167288786,
-        -0.022226685581313674,
-        0.31441213129613066,
-        -0.026739700851545107,
-        -0.016084530571308823,
-    ];
-
     assert!(solver.solution.x.dist(&refsol) <= 1e-6);
-
-    let refobj = 4.840076866013861;
     assert!(f64::abs(solver.info.cost_primal - refobj) <= 1e-6);
 }