Add active_set function

Author: SE-starshippilot

src/active_set.cpp

@@ -0,0 +1,100 @@
+#include "default_types.h"
+#include <igl/active_set.h>
+#include <nanobind/nanobind.h>
+#include <nanobind/ndarray.h>
+#include <nanobind/eigen/dense.h>
+#include <nanobind/eigen/sparse.h>
+
+namespace nb = nanobind;
+using namespace nb::literals;
+
+namespace pyigl
+{
+  // Wrapper for active_set function
+  auto active_set(
+    const Eigen::SparseMatrixN &A,
+    const nb::DRef<const Eigen::VectorXN> &B,
+    const nb::DRef<const Eigen::VectorXI> &known,
+    const nb::DRef<const Eigen::VectorXN> &Y,
+    const Eigen::SparseMatrixN &Aeq,
+    const nb::DRef<const Eigen::VectorXN> &Beq,
+    const Eigen::SparseMatrixN &Aieq,
+    const nb::DRef<const Eigen::VectorXN> &Bieq,
+    const nb::DRef<const Eigen::VectorXN> &lx,
+    const nb::DRef<const Eigen::VectorXN> &ux,
+    const bool Auu_pd,
+    const int max_iter,
+    const double inactive_threshold,
+    const double constraint_threshold,
+    const double solution_diff_threshold)
+  {
+    igl::active_set_params params;
+    params.Auu_pd = Auu_pd;
+    params.max_iter = max_iter;
+    params.inactive_threshold = inactive_threshold;
+    params.constraint_threshold = constraint_threshold;
+    params.solution_diff_threshold = solution_diff_threshold;
+
+    Eigen::VectorXN Z;
+    igl::SolverStatus status = igl::active_set(
+      A, B, known, Y, Aeq, Beq, Aieq, Bieq, lx, ux, params, Z);
+
+    if(status == igl::SOLVER_STATUS_ERROR)
+    {
+      throw std::runtime_error("active_set: solver failed");
+    }
+    return Z;
+  }
+}
+
+// Bind the wrapper to the Python module
+void bind_active_set(nb::module_ &m)
+{
+  igl::active_set_params params;
+  m.def(
+    "active_set",
+    &pyigl::active_set,
+    "A"_a,
+    "B"_a,
+    "known"_a      = Eigen::VectorXI(),
+    "Y"_a          = Eigen::VectorXN(),
+    "Aeq"_a        = Eigen::SparseMatrixN(),
+    "Beq"_a        = Eigen::VectorXN(),
+    "Aieq"_a       = Eigen::SparseMatrixN(),
+    "Bieq"_a       = Eigen::VectorXN(),
+    "lx"_a         = Eigen::VectorXN(),
+    "ux"_a         = Eigen::VectorXN(),
+    "Auu_pd"_a     = params.Auu_pd,
+    "max_iter"_a   = params.max_iter,
+    "inactive_threshold"_a      = params.inactive_threshold,
+    "constraint_threshold"_a    = params.constraint_threshold,
+    "solution_diff_threshold"_a = params.solution_diff_threshold,
+R"(Minimize a convex quadratic energy subject to linear constraints.
+
+Solves:
+    min  0.5 * Z' * A * Z + Z' * B
+     Z
+    subject to
+        Z(known) = Y
+        Aeq * Z = Beq
+        Aieq * Z <= Bieq
+        lx <= Z <= ux
+
+@param[in] A  n by n matrix of quadratic coefficients
+@param[in] B  n by 1 column of linear coefficients
+@param[in] known  list of indices to known rows in Z
+@param[in] Y  list of fixed values corresponding to known rows in Z
+@param[in] Aeq  meq by n list of linear equality constraint coefficients
+@param[in] Beq  meq by 1 list of linear equality constraint constant values
+@param[in] Aieq  mieq by n list of linear inequality constraint coefficients
+@param[in] Bieq  mieq by 1 list of linear inequality constraint constant values
+@param[in] lx  n by 1 list of lower bounds [] implies -Inf
+@param[in] ux  n by 1 list of upper bounds [] implies Inf
+@param[in] Auu_pd  whether Auu is positive definite {false}
+@param[in] max_iter  maximum number of iterations {100}
+@param[in] inactive_threshold  threshold on Lagrange multiplier values {1e-14}
+@param[in] constraint_threshold  threshold on constraint violation {1e-14}
+@param[in] solution_diff_threshold  threshold on solution difference {1e-14}
+@return Z  n by 1 solution vector)"
+  );
+}

tests/test_all.py

@@ -258,7 +258,25 @@ def test_min_quad():
     data = igl.min_quad_with_fixed_data()
     igl.min_quad_with_fixed_precompute(A,known,Aeq,True,data)
     Z = igl.min_quad_with_fixed_solve(data,B,Y,Beq)
-    
+
+def test_active_set():
+    V,F,T = single_tet()
+    n = V.shape[0]
+    A = -igl.cotmatrix(V,F)
+    B = np.zeros(n,dtype=np.float64)
+    known = np.array([0,1],dtype=np.int64)
+    Y = np.array([0.0,1.0],dtype=np.float64)
+    Aeq = scipy.sparse.csc_matrix((0,n),dtype=np.float64)
+    Beq = np.array([],dtype=np.float64)
+    Aieq = scipy.sparse.csc_matrix((0,n),dtype=np.float64)
+    Bieq = np.array([],dtype=np.float64)
+    lx = np.array([],dtype=np.float64)
+    ux = np.array([],dtype=np.float64)
+    Z = igl.active_set(A,B,known,Y,Aeq,Beq,Aieq,Bieq,lx,ux)
+    assert Z.shape[0] == n
+    assert Z[0] == pytest.approx(0.0)
+    assert Z[1] == pytest.approx(1.0)
+
 def test_volume():
     V = np.array([[0,0,0],[1,0,0],[0,1,0],[0,0,1]],dtype=np.float64)
     T = np.array([[0,1,2,3]],dtype=np.int64)