Adding missing tests for the C++ wrapper (#40)

* Adding missing tests for the C++ wrapper

The tests are adapted directly from the rust tests.

* Remove unnecessary spacing

tests/CMakeLists.txt

@@ -20,9 +20,16 @@ add_executable(clarabel_cpp_tests
     api_dimension_checks.cpp
     basic_lp.cpp
     basic_qp.cpp
+    basic_eq_constrained.cpp
+    basic_powcone.cpp
+    basic_genpowcone.cpp
+    basic_expcone.cpp
+    basic_sdp.cpp
+    basic_socp.cpp
     basic_unconstrained.cpp
     mixed_conic.cpp
     data_updating.cpp
+    sdp_chordal.cpp
 )
 target_link_libraries(clarabel_cpp_tests 
     libclarabel_c_shared

tests/basic_eq_constrained.cpp

@@ -0,0 +1,108 @@
+#include <Clarabel>
+#include <Eigen/Eigen>
+#include <cmath>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+using namespace std;
+using namespace clarabel;
+using namespace Eigen;
+
+SparseMatrix<double> eq_constrained_A1() {
+    // A =
+    //[ 0. 1.  1.;
+    //  0. 1. -1.]
+    int colptr[] = { 0, 0, 2, 4 }; // start index per column + 1 for last col
+    int rowval[] = { 0, 1, 0, 1 }; // nonzero row indices
+    double nzval[] = { 1., 1., 1., -1. };
+    return SparseMatrix<double>::Map(
+        2 /*rows*/, 3 /*cols*/, 4 /*nonzeros*/, colptr, rowval, nzval);
+}
+
+SparseMatrix<double> eq_constrained_A2() {
+    // A = [
+    // 0    1.0   1.0;
+    // 0    1.0  -1.0;
+    //1.0   2.0  -1.0l
+    //2.0  -1.0   3.0l
+    //]
+    int colptr[] = { 0, 2, 6, 10 }; // start index per column + 1 for last col
+    int rowval[] = { 2, 3, 0, 1, 2, 3, 0, 1, 2, 3}; // nonzero row indices
+    double nzval[] = {1., 2., 1., 1., 2., -1., 1., -1., -1., 3. };
+    return SparseMatrix<double>::Map(
+        4 /*rows*/, 3 /*cols*/, 10 /*nonzeros*/, colptr, rowval, nzval);
+}
+
+TEST(BasicEqConstrainedTest, Feasible)
+{
+    SparseMatrix<double> P = MatrixXd::Identity(3, 3).sparseView();
+    P.makeCompressed();
+
+    Vector<double, 3> c = { 0., 0., 0. };
+
+    SparseMatrix<double> A = eq_constrained_A1();
+
+    Vector<double, 2> b = {2., 0.};
+
+    vector<SupportedConeT<double>> cones = {ZeroConeT<double>(2)};
+
+    DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::Solved);
+
+    // Compare the solution to the reference solution
+    Vector3d ref_solution{ 0., 1., 1. };
+    ASSERT_EQ(solution.x.size(), 3);
+    ASSERT_TRUE(solution.x.isApprox(ref_solution, 1e-6));
+}
+
+TEST(BasicEqConstrainedTest, PrimalInfeasible)
+{
+    SparseMatrix<double> P = MatrixXd::Identity(3, 3).sparseView();
+    P.makeCompressed();
+
+    Vector<double, 3> c = { 0., 0., 0. };
+
+    SparseMatrix<double> A = eq_constrained_A2();
+
+    Vector<double, 4> b = {1., 1., 1., 1.};
+
+    vector<SupportedConeT<double>> cones = {ZeroConeT<double>(4)};
+
+    DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::PrimalInfeasible);
+}
+
+TEST(BasicEqConstrainedTest, DualInfeasible)
+{
+    SparseMatrix<double> P = MatrixXd::Identity(3, 3).sparseView();
+    P.makeCompressed();
+    P.coeffRef(0,0) = 0.0;
+
+    Vector<double, 3> c = { 1., 1., 1. };
+
+    SparseMatrix<double> A = eq_constrained_A1();
+
+    Vector<double, 2> b = {2., 0.};
+
+    vector<SupportedConeT<double>> cones = {ZeroConeT<double>(2)};
+
+    DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::DualInfeasible);
+}

tests/basic_expcone.cpp

@@ -0,0 +1,107 @@
+#include <Clarabel>
+#include <Eigen/Eigen>
+#include <cmath>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+using namespace std;
+using namespace clarabel;
+using namespace Eigen;
+
+class BasicExpConeTest : public testing::Test
+{
+    // produces data for the following exponential cone problem
+    // max  x
+    // s.t. y * exp(x / y) <= z
+    //      y == 1, z == exp(5)
+  protected:
+    SparseMatrix<double> P, A;
+    Vector<double, 3> c = {-1., 0., 0.};
+    Vector<double, 5> b = {0., 0., 0., 1.0, exp(5.)};
+    vector<SupportedConeT<double>> cones = {
+        ExponentialConeT<double>(),
+        ZeroConeT<double>(2),
+    };
+    DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
+
+    BasicExpConeTest()
+    {
+        P = MatrixXd::Zero(3, 3).sparseView();
+        P.makeCompressed();
+
+        MatrixXd A1 = -MatrixXd::Identity(3,3);
+        MatrixXd A2(2,3);
+        A2 << 0, 1, 0, // y = 1
+              0, 0, 1; // z = exp(5)
+
+        MatrixXd A_dense = MatrixXd::Zero(5,3);
+        A_dense << A1, A2;
+        A = A_dense.sparseView();
+        A.makeCompressed();
+    }
+};
+
+TEST_F(BasicExpConeTest, Feasible)
+{
+    // solve the following exponential cone problem
+    // max  x
+    // s.t. y * exp(x / y) <= z
+    //      y == 1, z == exp(5)
+    //
+    // This is just the default problem data above
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::Solved);
+
+    // Check the solution
+    Vector3d ref_solution{ 5.0, 1.0, exp(5.0) };
+    ASSERT_EQ(solution.x.size(), 3);
+    ASSERT_TRUE(solution.x.isApprox(ref_solution, 1e-6));
+
+    double ref_obj = -5.0;
+    ASSERT_NEAR(solution.obj_val, ref_obj, 1e-6);
+    ASSERT_NEAR(solution.obj_val_dual, ref_obj, 1e-6);
+}
+
+TEST_F(BasicExpConeTest, PrimalInfeasible)
+{
+    // solve the following exponential cone problem
+    // max  x
+    // s.t. y * exp(x / y) <= z
+    //      y == 1, z == -1
+    //
+    // Same as default, but last element of b is different
+    b[4] = -1.;
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::PrimalInfeasible);
+}
+
+
+TEST_F(BasicExpConeTest, DualInfeasible)
+{
+    // solve the following exponential cone problem
+    // max  x
+    // s.t. y * exp(x / y) <= z
+    //
+    // Same as default, but no equality constraint
+
+    Vector<double, 3> b = {0., 0., 0.};
+    vector<SupportedConeT<double>> cones = {ExponentialConeT<double>()};
+    MatrixXd A1 = -MatrixXd::Identity(3,3);
+    SparseMatrix<double> A = A1.sparseView();
+    A.makeCompressed();
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::DualInfeasible);
+}

tests/basic_genpowcone.cpp

@@ -0,0 +1,64 @@
+#include <Clarabel>
+#include <Eigen/Eigen>
+#include <cmath>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+using namespace std;
+using namespace clarabel;
+using namespace Eigen;
+
+TEST(BasicGenPowerConeTest, Feasible)
+{
+    // solve the following power cone problem
+    // max  x1^0.6 y^0.4 + x2^0.1
+    // s.t. x1, y, x2 >= 0
+    //      x1 + 2y  + 3x2 == 3
+    // which is equivalent to
+    // max z1 + z2
+    // s.t. (x1, y, z1) in K_pow(0.6)
+    //      (x2, 1, z2) in K_pow(0.1)
+    //      x1 + 2y + 3x2 == 3
+
+    // x = (x1, y, z1, x2, y2, z2)
+
+    SparseMatrix<double> P = MatrixXd::Zero(6, 6).sparseView();
+    P.makeCompressed();
+
+    Vector<double, 6> c = { 0., 0., -1., 0., 0., -1. };
+
+    // Assembling A
+    MatrixXd A1_dense = -MatrixXd::Identity(6, 6);
+    MatrixXd A2_dense(2,6);
+    A2_dense << 1., 2., 0., 3., 0., 0., //
+                0., 0., 0., 0., 1., 0.; //
+
+    MatrixXd A_dense = MatrixXd::Zero(8,6);
+    A_dense << A1_dense, A2_dense;
+    SparseMatrix<double> A = A_dense.sparseView();
+    A.makeCompressed();
+
+    // Assembling b
+    Vector<double, 8> b = {0., 0., 0., 0., 0., 0., 3., 1.};
+    // Assembling cones
+
+    Vector<double, 2> cone1 = {0.6, 0.4};
+    Vector<double, 2> cone2 = {0.1, 0.9};
+    vector<SupportedConeT<double>> cones = {GenPowerConeT<double>(cone1, 1), 
+                                            GenPowerConeT<double>(cone2, 1), 
+                                            ZeroConeT<double>(2)};
+
+    DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::Solved);
+
+    // Compare the solution to the reference solution
+    double ref_obj = -1.8458;
+    ASSERT_NEAR(solution.obj_val, ref_obj, 1e-3);
+}

tests/basic_powcone.cpp

@@ -0,0 +1,61 @@
+#include <Clarabel>
+#include <Eigen/Eigen>
+#include <cmath>
+#include <gtest/gtest.h>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+using namespace std;
+using namespace clarabel;
+using namespace Eigen;
+
+TEST(BasicPowerConeTest, Feasible)
+{
+    // solve the following power cone problem
+    // max  x1^0.6 y^0.4 + x2^0.1
+    // s.t. x1, y, x2 >= 0
+    //      x1 + 2y  + 3x2 == 3
+    // which is equivalent to
+    // max z1 + z2
+    // s.t. (x1, y, z1) in K_pow(0.6)
+    //      (x2, 1, z2) in K_pow(0.1)
+    //      x1 + 2y + 3x2 == 3
+
+    // x = (x1, y, z1, x2, y2, z2)
+
+    SparseMatrix<double> P = MatrixXd::Zero(6, 6).sparseView();
+    P.makeCompressed();
+
+    Vector<double, 6> c = { 0., 0., -1., 0., 0., -1. };
+
+    // Assembling A
+    MatrixXd A1_dense = -MatrixXd::Identity(6, 6);
+    MatrixXd A2_dense(2,6);
+    A2_dense << 1., 2., 0., 3., 0., 0., //
+                0., 0., 0., 0., 1., 0.; //
+
+    MatrixXd A_dense = MatrixXd::Zero(8,6);
+    A_dense << A1_dense, A2_dense;
+    SparseMatrix<double> A = A_dense.sparseView();
+    A.makeCompressed();
+
+    // Assembling b
+    Vector<double, 8> b = {0., 0., 0., 0., 0., 0., 3., 1.};
+    // Assembling cones
+    vector<SupportedConeT<double>> cones = {PowerConeT<double>(0.6), 
+                                             PowerConeT<double>(0.1), 
+                                             ZeroConeT<double>(2)};
+
+    DefaultSettings<double> settings = DefaultSettings<double>::default_settings();
+
+    DefaultSolver<double> solver(P, c, A, b, cones, settings);
+    solver.solve();
+
+    DefaultSolution<double> solution = solver.solution();
+    ASSERT_EQ(solution.status, SolverStatus::Solved);
+
+    // Compare the solution to the reference solution
+    double ref_obj = -1.8458;
+    ASSERT_NEAR(solution.obj_val, ref_obj, 1e-3);
+}