diff --git a/python/305_QuadraticProgramming.py b/python/305_QuadraticProgramming.py new file mode 100644 index 000000000..df941644c --- /dev/null +++ b/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() diff --git a/python/py_igl.cpp b/python/py_igl.cpp index fafecc74f..596281053 100644 --- a/python/py_igl.cpp +++ b/python/py_igl.cpp @@ -29,6 +29,8 @@ #include #include #include +#include +#include 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" } diff --git a/python/py_igl/py_SolverStatus.cpp b/python/py_igl/py_SolverStatus.cpp new file mode 100644 index 000000000..1e5ccc5d1 --- /dev/null +++ b/python/py_igl/py_SolverStatus.cpp @@ -0,0 +1,6 @@ +py::enum_(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(); diff --git a/python/py_igl/py_active_set.cpp b/python/py_igl/py_active_set.cpp new file mode 100644 index 000000000..e9cba466e --- /dev/null +++ b/python/py_igl/py_active_set.cpp @@ -0,0 +1,48 @@ +// Wrap the params struct +py::class_ 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& A, + const Eigen::MatrixXd& B, + const Eigen::MatrixXi& known, + const Eigen::MatrixXd& Y, + const Eigen::SparseMatrix& Aeq, + const Eigen::MatrixXd& Beq, + const Eigen::SparseMatrix& 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")); diff --git a/python/py_igl/py_min_quad_with_fixed.cpp b/python/py_igl/py_min_quad_with_fixed.cpp index 01bdff0f2..3780c7492 100644 --- a/python/py_igl/py_min_quad_with_fixed.cpp +++ b/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 & 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")); diff --git a/python/py_vector.cpp b/python/py_vector.cpp index 3a7931fbc..ecaaffa63 100644 --- a/python/py_vector.cpp +++ b/python/py_vector.cpp @@ -146,6 +146,8 @@ py::class_ 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(m.asDiagonal());} ) + .def("sparseView", [](Type &m) { return Eigen::SparseMatrix(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_ bind_eigen_sparse_2(py::module &m, const char *name, oss << v; return oss.str(); }) - // .def("__getitem__", [](const Type &m, std::pair 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 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::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_ bind_eigen_sparse_2(py::module &m, const char *name, return m.makeCompressed(); }) + .def("diagonal", [](const Type &m) {return Eigen::Matrix(m.diagonal());} ) + ; return matrix; } @@ -532,10 +511,12 @@ py::class_ 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; } diff --git a/python/python.h b/python/python.h index e5a0136f8..84dd86c46 100644 --- a/python/python.h +++ b/python/python.h @@ -14,6 +14,9 @@ template void assert_is_VectorX(const std::string name, const Eigen::PlainObjectBase& 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 void assert_is_RowVectorX(const std::string name, const Eigen::PlainObjectBase& 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 void assert_is_Vector3(const std::string name, const Eigen::PlainObjectBase& 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 void assert_is_RowVector3(const std::string name, const Eigen::PlainObjectBase& 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 void assert_is_Vector4(const std::string name, const Eigen::PlainObjectBase& 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 void assert_is_RowVector4(const std::string name, const Eigen::PlainObjectBase& 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 void assert_is_Matrix4(const std::string name, const Eigen::PlainObjectBase& v) { + if (v.size() == 0) + return; + if ((v.cols() != 4) || (v.rows() != 4)) throw std::runtime_error(name + " must be a 4x4 matrix."); }