added 305

Former-commit-id: 477339d410

python/305_QuadraticProgramming.py

@@ -0,0 +1,86 @@
+import igl
+
+b = igl.eigen.MatrixXi()
+B = igl.eigen.MatrixXd()
+bc = igl.eigen.MatrixXd()
+lx = igl.eigen.MatrixXd()
+ux = igl.eigen.MatrixXd()
+Beq = igl.eigen.MatrixXd()
+Bieq = igl.eigen.MatrixXd()
+Z = igl.eigen.MatrixXd()
+
+Q = igl.eigen.SparseMatrixd()
+Aeq = igl.eigen.SparseMatrixd()
+Aieq = igl.eigen.SparseMatrixd()
+
+def solve(viewer):
+    global Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,Z
+    params = igl.active_set_params()
+    params.max_iter = 8
+
+    igl.active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,params,Z)
+
+    C = igl.eigen.MatrixXd()
+    igl.jet(Z,0,1,C)
+    viewer.data.set_colors(C)
+
+def key_down(viewer, key, mod):
+    global Beq,solve
+    if key == ord('.'):
+         Beq[0,0] = Beq[0,0] * 2.0
+         solve(viewer)
+         return True
+    elif key == ord(','):
+         Beq[0,0] = Beq[0,0] / 2.0
+         solve(viewer)
+         return True
+    elif key == ord(' '):
+         solve(viewer)
+         return True
+    return False;
+
+
+V = igl.eigen.MatrixXd()
+F = igl.eigen.MatrixXi()
+
+igl.readOFF("../tutorial/shared/cheburashka.off",V,F)
+
+# Plot the mesh
+viewer = igl.viewer.Viewer()
+viewer.data.set_mesh(V, F)
+viewer.core.show_lines = False
+viewer.callback_key_down = key_down
+
+# One fixed point on belly
+b  = igl.eigen.MatrixXi([[2556]])
+bc = igl.eigen.MatrixXd([[1]])
+
+# Construct Laplacian and mass matrix
+L = igl.eigen.SparseMatrixd()
+M = igl.eigen.SparseMatrixd()
+Minv = igl.eigen.SparseMatrixd()
+
+igl.cotmatrix(V,F,L)
+igl.massmatrix(V,F,igl.MASSMATRIX_TYPE_VORONOI,M);
+igl.invert_diag(M,Minv)
+
+# Bi-Laplacian
+Q = L.transpose() * (Minv * L)
+
+# Zero linear term
+B = igl.eigen.MatrixXd.Zero(V.rows(),1)
+
+# Lower and upper bound
+lx = igl.eigen.MatrixXd.Zero(V.rows(),1)
+ux = igl.eigen.MatrixXd.Ones(V.rows(),1)
+
+# Equality constraint constrain solution to sum to 1
+Beq = igl.eigen.MatrixXd([[0.08]])
+Aeq = M.diagonal().transpose().sparseView()
+
+# (Empty inequality constraints)
+solve(viewer)
+print("Press '.' to increase scale and resolve.")
+print("Press ',' to decrease scale and resolve.")
+
+viewer.launch()

python/py_igl.cpp

@@ -29,6 +29,8 @@
 #include <igl/unique.h>
 #include <igl/setdiff.h>
 #include <igl/min_quad_with_fixed.h>
+#include <igl/SolverStatus.h>
+#include <igl/active_set.h>
 
 void python_export_igl(py::module &m)
 {
@@ -59,5 +61,7 @@ void python_export_igl(py::module &m)
 #include "py_igl/py_unique.cpp"
 #include "py_igl/py_setdiff.cpp"
 #include "py_igl/py_min_quad_with_fixed.cpp"
+#include "py_igl/py_SolverStatus.cpp"
+#include "py_igl/py_active_set.cpp"
 
 }

python/py_igl/py_SolverStatus.cpp

@@ -0,0 +1,6 @@
+py::enum_<igl::SolverStatus>(m, "SolverStatus")
+    .value("SOLVER_STATUS_CONVERGED", igl::SOLVER_STATUS_CONVERGED)
+    .value("SOLVER_STATUS_MAX_ITER", igl::SOLVER_STATUS_MAX_ITER)
+    .value("SOLVER_STATUS_ERROR", igl::SOLVER_STATUS_ERROR)
+    .value("NUM_SOLVER_STATUSES", igl::NUM_SOLVER_STATUSES)
+    .export_values();

python/py_igl/py_active_set.cpp

@@ -0,0 +1,48 @@
+// Wrap the params struct
+py::class_<igl::active_set_params > active_set_params(m, "active_set_params");
+
+active_set_params
+.def("__init__", [](igl::active_set_params &m)
+{
+    new (&m) igl::active_set_params();
+    m.Auu_pd = false;
+    m.max_iter = 100;
+    m.inactive_threshold = igl::DOUBLE_EPS;
+    m.constraint_threshold = igl::DOUBLE_EPS;
+    m.solution_diff_threshold = igl::DOUBLE_EPS;
+})
+.def_readwrite("Auu_pd", &igl::active_set_params::Auu_pd)
+.def_readwrite("max_iter", &igl::active_set_params::max_iter)
+.def_readwrite("inactive_threshold", &igl::active_set_params::inactive_threshold)
+.def_readwrite("constraint_threshold", &igl::active_set_params::constraint_threshold)
+.def_readwrite("solution_diff_threshold", &igl::active_set_params::solution_diff_threshold)
+.def_readwrite("Auu_pd", &igl::active_set_params::Auu_pd)
+;
+
+m.def("active_set", []
+(
+  const Eigen::SparseMatrix<double>& A,
+  const Eigen::MatrixXd& B,
+  const Eigen::MatrixXi& known,
+  const Eigen::MatrixXd& Y,
+  const Eigen::SparseMatrix<double>& Aeq,
+  const Eigen::MatrixXd& Beq,
+  const Eigen::SparseMatrix<double>& Aieq,
+  const Eigen::MatrixXd& Bieq,
+  const Eigen::MatrixXd& lx,
+  const Eigen::MatrixXd& ux,
+  const igl::active_set_params& params,
+  Eigen::MatrixXd& Z
+)
+{
+  assert_is_VectorX("B",B);
+  assert_is_VectorX("known",known);
+  assert_is_VectorX("Y",Y);
+  assert_is_VectorX("Beq",Beq);
+  assert_is_VectorX("Bieq",Bieq);
+  assert_is_VectorX("Z",Z);
+
+  return igl::active_set(A,B,known,Y,Aeq,Eigen::VectorXd(Beq),Aieq,Bieq,lx,ux,params,Z);
+}, __doc_igl_active_set,
+py::arg("A"), py::arg("B"), py::arg("known"), py::arg("Y"), py::arg("Aeq"), py::arg("Beq")
+, py::arg("Aieq"), py::arg("Bieq"), py::arg("lx"), py::arg("ux"), py::arg("params"), py::arg("Z"));

python/py_igl/py_min_quad_with_fixed.cpp

@@ -14,6 +14,7 @@ m.def("min_quad_with_fixed_precompute", []
   igl::min_quad_with_fixed_data<double> & data
 )
 {
+  assert_is_VectorX("known",known);
   return igl::min_quad_with_fixed_precompute(A,known,Aeq,pd,data);
 }, __doc_igl_min_quad_with_fixed,
 py::arg("A"), py::arg("known"), py::arg("Aeq"), py::arg("pd"), py::arg("data"));
@@ -28,6 +29,9 @@ m.def("min_quad_with_fixed_solve", []
   Eigen::MatrixXd& sol
 )
 {
+  assert_is_VectorX("B",B);
+  assert_is_VectorX("Y",Y);
+  assert_is_VectorX("Beq",Beq);
   return igl::min_quad_with_fixed_solve(data,B,Y,Beq,Z,sol);
 }, __doc_igl_min_quad_with_fixed,
 py::arg("data"), py::arg("B"), py::arg("Y"), py::arg("Beq"), py::arg("Z"), py::arg("sol"));
@@ -41,6 +45,9 @@ m.def("min_quad_with_fixed_solve", []
   Eigen::MatrixXd& Z
 )
 {
+  assert_is_VectorX("B",B);
+  assert_is_VectorX("Y",Y);
+  assert_is_VectorX("Beq",Beq);
   return igl::min_quad_with_fixed_solve(data,B,Y,Beq,Z);
 }, __doc_igl_min_quad_with_fixed,
 py::arg("data"), py::arg("B"), py::arg("Y"), py::arg("Beq"), py::arg("Z"));
@@ -57,6 +64,10 @@ m.def("min_quad_with_fixed", []
   Eigen::MatrixXd& Z
 )
 {
+  assert_is_VectorX("B",B);
+  assert_is_VectorX("known",known);
+  assert_is_VectorX("Y",Y);
+  assert_is_VectorX("Beq",Beq);
   return igl::min_quad_with_fixed(A,B,known,Y,Aeq,Beq,pd,Z);
 }, __doc_igl_min_quad_with_fixed,
 py::arg("A"), py::arg("B"), py::arg("known"), py::arg("Y"), py::arg("Aeq"), py::arg("Beq"), py::arg("pd"), py::arg("Z"));

python/py_vector.cpp

@@ -146,6 +146,8 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name,
         .def("replicate", [](const Type &m, const int& r, const int& c) {return Type(m.replicate(r,c));} )
         .def("asDiagonal", [](const Type &m) {return Eigen::DiagonalMatrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(m.asDiagonal());} )
 
+        .def("sparseView", [](Type &m) { return Eigen::SparseMatrix<Scalar>(m.sparseView()); })
+
         /* Arithmetic operators (def_cast forcefully casts the result back to a
            Type to avoid type issues with Eigen's crazy expression templates) */
         .def_cast(-py::self)
@@ -391,31 +393,6 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
             oss << v;
             return oss.str();
         })
-        // .def("__getitem__", [](const Type &m, std::pair<size_t, size_t> i) {
-        //     if (i.first >= (size_t) m.rows() || i.second >= (size_t) m.cols())
-        //         throw py::index_error();
-        //     return m(i.first, i.second);
-        //  })
-        // .def("__setitem__", [](Type &m, std::pair<size_t, size_t> i, Scalar v) {
-        //     if (i.first >= (size_t) m.rows() || i.second >= (size_t) m.cols())
-        //         throw py::index_error();
-        //     m(i.first, i.second) = v;
-        //  })
-
-        // /* Buffer access for interacting with NumPy */
-        // .def_buffer([](Type &m) -> py::buffer_info {
-        //     return py::buffer_info(
-        //         m.data(),                /* Pointer to buffer */
-        //         sizeof(Scalar),          /* Size of one scalar */
-        //         /* Python struct-style format descriptor */
-        //         py::format_descriptor<Scalar>::value(),
-        //         2,                       /* Number of dimensions */
-        //         { (size_t) m.rows(),     /* Buffer dimensions */
-        //           (size_t) m.cols() },
-        //         { sizeof(Scalar),        /* Strides (in bytes) for each index */
-        //           sizeof(Scalar) * m.rows() }
-        //     );
-        //  })
 
         /* Static initializers */
         // .def_static("Zero", [](size_t n, size_t m) { return Type(Type::Zero(n, m)); })
@@ -460,6 +437,8 @@ py::class_<Type> bind_eigen_sparse_2(py::module &m, const char *name,
           return m.makeCompressed();
         })
 
+        .def("diagonal", [](const Type &m) {return Eigen::Matrix<Scalar,Eigen::Dynamic,Eigen::Dynamic>(m.diagonal());} )
+
         ;
     return matrix;
 }
@@ -532,10 +511,12 @@ py::class_<Type> bind_eigen_diagonal_2(py::module &m, const char *name,
         /* Python protocol implementations */
         .def("__repr__", [](const Type &/*v*/) {
             std::ostringstream oss;
-            oss << "Eigen is not able to plot Diagonal Matrices";
+            oss << "<< operator undefined for diagonal matrices";
             return oss.str();
         })
 
+        /* Other transformations */
+
         ;
     return matrix;
 }

python/python.h

@@ -14,6 +14,9 @@
 template<typename Scalar>
 void assert_is_VectorX(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if (v.cols() != 1)
     throw std::runtime_error(name + " must be a column vector.");
 }
@@ -21,6 +24,9 @@ void assert_is_VectorX(const std::string name, const Eigen::PlainObjectBase<Scal
 template<typename Scalar>
 void assert_is_RowVectorX(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if (v.rows() != 1)
     throw std::runtime_error(name + " must be a row vector.");
 }
@@ -28,6 +34,9 @@ void assert_is_RowVectorX(const std::string name, const Eigen::PlainObjectBase<S
 template<typename Scalar>
 void assert_is_Vector3(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if ((v.cols() != 1) || (v.rows() != 3))
     throw std::runtime_error(name + " must be a column vector with 3 entries.");
 }
@@ -35,6 +44,9 @@ void assert_is_Vector3(const std::string name, const Eigen::PlainObjectBase<Scal
 template<typename Scalar>
 void assert_is_RowVector3(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if ((v.cols() != 3) || (v.rows() != 1))
     throw std::runtime_error(name + " must be a row vector with 3 entries.");
 }
@@ -42,6 +54,9 @@ void assert_is_RowVector3(const std::string name, const Eigen::PlainObjectBase<S
 template<typename Scalar>
 void assert_is_Vector4(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if ((v.cols() != 1) || (v.rows() != 4))
     throw std::runtime_error(name + " must be a column vector with 4 entries.");
 }
@@ -49,6 +64,9 @@ void assert_is_Vector4(const std::string name, const Eigen::PlainObjectBase<Scal
 template<typename Scalar>
 void assert_is_RowVector4(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if ((v.cols() != 4) || (v.rows() != 1))
     throw std::runtime_error(name + " must be a row vector with 4 entries.");
 }
@@ -56,6 +74,9 @@ void assert_is_RowVector4(const std::string name, const Eigen::PlainObjectBase<S
 template<typename Scalar>
 void assert_is_Matrix4(const std::string name, const Eigen::PlainObjectBase<Scalar>& v)
 {
+  if (v.size() == 0)
+    return;
+
   if ((v.cols() != 4) || (v.rows() != 4))
     throw std::runtime_error(name + " must be a 4x4 matrix.");
 }