bug fix in active set, qp tutorial example

Former-commit-id: 987fd623cc5aac9cf951ad2fb10d80942555a618

include/igl/active_set.cpp

@@ -224,6 +224,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     // Gather active constraints and resp. rhss
     PlainObjectBase<DerivedBeq> Beq_i;
     Beq_i.resize(Beq.rows()+as_ieq_count,1);
+    Beq_i.head(Beq.rows()) = Beq;
     {
       int k =0;
       for(int a=0;a<as_ieq.size();a++)
@@ -282,6 +283,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       ret = SOLVER_STATUS_ERROR;
       break;
     }
+    //cout<<matlab_format((Aeq*Z-Beq).eval(),"cr")<<endl;
     //cout<<matlab_format(Z,"Z")<<endl;
 #ifdef ACTIVE_SET_CPP_DEBUG
     cout<<"  post"<<endl;

include/igl/find.cpp

@@ -66,4 +66,5 @@ IGL_INLINE void igl::find(
 template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::find<double, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
+template void igl::find<double, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/is_symmetric.cpp

@@ -50,6 +50,7 @@ IGL_INLINE bool igl::is_symmetric(
 {
   using namespace Eigen;
   using namespace igl;
+  using namespace std;
   if(A.rows() != A.cols())
   {
     return false;
@@ -60,6 +61,10 @@ IGL_INLINE bool igl::is_symmetric(
   VectorXi AmATI,AmATJ;
   Matrix<AType,Dynamic,1> AmATV;
   find(AmAT,AmATI,AmATJ,AmATV);
+  if(AmATI.size() == 0)
+  {
+    return true;
+  }
   
   return AmATV.maxCoeff() < epsilon && AmATV.minCoeff() > -epsilon;
 }

include/igl/min_quad_with_fixed.cpp

@@ -12,7 +12,7 @@
 #include "find.h"
 #include "sparse.h"
 #include "repmat.h"
-#include "lu_lagrange.h"
+//#include "lu_lagrange.h"
 #include "full.h"
 #include "matlab_format.h"
 #include "EPS.h"
@@ -109,7 +109,10 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   }else
   {
     // determine if A(unknown,unknown) is symmetric and/or positive definite
-    data.Auu_sym = is_symmetric(Auu,EPS<double>());
+    VectorXi AuuI,AuuJ;
+    MatrixXd AuuV;
+    find(Auu,AuuI,AuuJ,AuuV);
+    data.Auu_sym = is_symmetric(Auu,EPS<double>()*AuuV.maxCoeff());
   }
 
   // Determine number of linearly independent constraints
@@ -144,6 +147,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     nc = data.AeqTQR.rank();
     assert(nc<=neq);
     data.Aeq_li = nc == neq;
+    //cout<<"data.Aeq_li: "<<data.Aeq_li<<endl;
   }else
   {
     data.Aeq_li = true;
@@ -378,6 +382,8 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
   using namespace std;
   using namespace Eigen;
   using namespace igl;
+  typedef Matrix<T,Dynamic,1> VectorXT;
+  typedef Matrix<T,Dynamic,Dynamic> MatrixXT;
   // number of known rows
   int kr = data.known.size();
   if(kr!=0)
@@ -402,56 +408,55 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 
   if(data.Aeq_li)
   {
+    // number of lagrange multipliers aka linear equality constraints
+    int neq = data.lagrange.size();
+    // append lagrange multiplier rhs's
+    VectorXT BBeq(B.size() + Beq.size());
+    BBeq << B, (Beq*-2.0);
+    // Build right hand side
+    VectorXT BBequl;
+    igl::slice(BBeq,data.unknown_lagrange,BBequl);
+    MatrixXT BBequlcols;
+    igl::repmat(BBequl,1,cols,BBequlcols);
+    MatrixXT NB;
+    if(kr == 0)
+    {
+      NB = BBequlcols;
+    }else
+    {
+      NB = data.preY * Y + BBequlcols;
+    }
 
-  // number of lagrange multipliers aka linear equality constraints
-  int neq = data.lagrange.size();
-  // append lagrange multiplier rhs's
-  Eigen::Matrix<T,Eigen::Dynamic,1> BBeq(B.size() + Beq.size());
-  BBeq << B, (Beq*-2.0);
-  // Build right hand side
-  Eigen::Matrix<T,Eigen::Dynamic,1> BBequl;
-  igl::slice(BBeq,data.unknown_lagrange,BBequl);
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> BBequlcols;
-  igl::repmat(BBequl,1,cols,BBequlcols);
-  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> NB;
-  if(kr == 0)
-  {
-    NB = BBequlcols;
-  }else
-  {
-    NB = data.preY * Y + BBequlcols;
-  }
-
-  //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
-  //cout<<matlab_format(NB,"NB")<<endl;
-  switch(data.solver_type)
-  {
-    case igl::min_quad_with_fixed_data<T>::LLT:
-      sol = data.llt.solve(NB);
-      break;
-    case igl::min_quad_with_fixed_data<T>::LDLT:
-      sol = data.ldlt.solve(NB);
-      break;
-    case igl::min_quad_with_fixed_data<T>::LU:
-      // Not a bottleneck
-      sol = data.lu.solve(NB);
-      break;
-    default:
-      cerr<<"Error: invalid solver type"<<endl;
-      return false;
-  }
-  //std::cout<<"sol=["<<std::endl<<sol<<std::endl<<"];"<<std::endl;
-  // Now sol contains sol/-0.5
-  sol *= -0.5;
-  // Now sol contains solution
-  // Place solution in Z
-  for(int i = 0;i<(sol.rows()-neq);i++)
-  {
-    for(int j = 0;j<sol.cols();j++)
+    //std::cout<<"NB=["<<std::endl<<NB<<std::endl<<"];"<<std::endl;
+    //cout<<matlab_format(NB,"NB")<<endl;
+    switch(data.solver_type)
     {
-      Z(data.unknown_lagrange(i),j) = sol(i,j);
+      case igl::min_quad_with_fixed_data<T>::LLT:
+        sol = data.llt.solve(NB);
+        break;
+      case igl::min_quad_with_fixed_data<T>::LDLT:
+        sol = data.ldlt.solve(NB);
+        break;
+      case igl::min_quad_with_fixed_data<T>::LU:
+        // Not a bottleneck
+        sol = data.lu.solve(NB);
+        break;
+      default:
+        cerr<<"Error: invalid solver type"<<endl;
+        return false;
+    }
+    //std::cout<<"sol=["<<std::endl<<sol<<std::endl<<"];"<<std::endl;
+    // Now sol contains sol/-0.5
+    sol *= -0.5;
+    // Now sol contains solution
+    // Place solution in Z
+    for(int i = 0;i<(sol.rows()-neq);i++)
+    {
+      for(int j = 0;j<sol.cols();j++)
+      {
+        Z(data.unknown_lagrange(i),j) = sol(i,j);
+      }
     }
-  }
   }else
   {
     assert(data.solver_type == min_quad_with_fixed_data<T>::QR_LLT);

tutorial/305_QuadraticProgramming/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(305_QuadraticProgramming)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/305_QuadraticProgramming/main.cpp

@@ -0,0 +1,94 @@
+#include <igl/active_set.h>
+#include <igl/boundary_faces.h>
+#include <igl/cotmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/jet.h>
+#include <igl/massmatrix.h>
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+
+  
+Eigen::VectorXi b;
+Eigen::VectorXd B,bc,lx,ux,Beq,Bieq,Z;
+Eigen::SparseMatrix<double> Q,Aeq,Aieq;
+
+void solve(igl::Viewer &viewer)
+{
+  using namespace std;
+  igl::active_set_params as;
+  as.max_iter = 8;
+  igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
+  // Pseudo-color based on solution
+  Eigen::MatrixXd C;
+  igl::jet(Z,0,1,C);
+  viewer.set_colors(C);
+}
+
+bool key_down(igl::Viewer &viewer, unsigned char key, int mod)
+{
+  switch(key)
+  {
+    case '.':
+      Beq(0) *= 2.0;
+      solve(viewer);
+      return true;
+    case ',':
+      Beq(0) /= 2.0;
+      solve(viewer);
+      return true;
+    case ' ':
+      solve(viewer);
+      return true;
+    default:
+      return false;
+  }
+}
+
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/cheburashka.off",V,F);
+
+  // Plot the mesh
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.options.show_lines = false;
+  viewer.callback_key_down = &key_down;
+
+  // One fixed point
+  b.resize(1,1);
+  // point on belly.
+  b<<2556;
+  bc.resize(1,1);
+  bc<<1;
+
+  // Construct Laplacian and mass matrix
+  SparseMatrix<double> L,M,Minv;
+  igl::cotmatrix(V,F,L);
+  igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
+  //M = (M/M.diagonal().maxCoeff()).eval();
+  igl::invert_diag(M,Minv);
+  // Bi-Laplacian
+  Q = L.transpose() * (Minv * L);
+  // Zero linear term
+  B = VectorXd::Zero(V.rows(),1);
+
+  // Lower and upper bound
+  lx = VectorXd::Zero(V.rows(),1);
+  ux = VectorXd::Ones(V.rows(),1);
+
+  // Equality constraint constrain solution to sum to 1
+  Beq.resize(1,1);
+  Beq(0) = 0.08;
+  Aeq = M.diagonal().transpose().sparseView();
+  // (Empty inequality constraints)
+  solve(viewer);
+
+  viewer.launch();
+}

tutorial/images/cheburashka-multiscale-biharmonic-kernels.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+e81844bdd47fc074d33e369495ea84569e66856d

tutorial/readme.md

@@ -37,7 +37,7 @@ of these lecture notes links to a cross-platform example application.
     * [204 Laplacian](#laplacian)
         * [Mass matrix](#massmatrix)
         * [Alternative construction of
-          Laplacian](#alternativeconstuctionoflaplacian)
+          Laplacian](#alternativeconstructionoflaplacian)
 * [Chapter 3: Matrices and Linear Algebra](#chapter3:matricesandlinearalgebra)
     * [301 Slice](#slice)
     * [302 Sort](#sort)
@@ -45,6 +45,15 @@ of these lecture notes links to a cross-platform example application.
     * [303 Laplace Equation](#laplaceequation)
         * [Quadratic energy minimization](#quadraticenergyminimization)
     * [304 Linear Equality Constraints](#linearequalityconstraints)
+    * [305 Quadratic Programming](#quadraticprogramming)
+* [Chapter 4: Shape Deformation](#chapter4:shapedeformation)
+    * [401 Biharmonic Deformation](#biharmonicdeformation)
+    * [402 Bounded Biharmonic Weights](#boundedbiharmonicweights)
+    * [403 Dual Quaternion Skinning](#dualquaternionskinning)
+    * [404 As-rigid-as-possible](#as-rigid-as-possible)
+    * [405 Fast automatic skinning
+      transformations](#fastautomaticskinningtransformations)
+
 
 # Compilation Instructions
 
@@ -370,7 +379,7 @@ An alternative construction of the discrete cotangent Laplacian is by
 "squaring" the discrete gradient operator. This may be derived by applying
 Green's identity (ignoring boundary conditions for the moment):
 
-  $\int_S \nabla f \nabla f dA = \int_S f \Delta f dA$
+  $\int_S \|\nabla f\|^2 dA = \int_S f \Delta f dA$
 
 Or in matrix form which is immediately translatable to code:
 
@@ -647,8 +656,9 @@ Notice that we can rewrite the last constraint in the familiar form from above:
 
  $z_{c} - z_{d} = 0.$
 
-Now we can assembly `Aeq` as a $1 \times n$ sparse matrix with a coefficient 1
-in the column corresponding to vertex $c$ and a -1 at $d$. The right-hand side
+Now we can assembly `Aeq` as a $1 \times n$ sparse matrix with a coefficient
+$1$
+in the column corresponding to vertex $c$ and a $-1$ at $d$. The right-hand side
 `Beq` is simply zero.
 
 Internally, `min_quad_with_fixed_*` solves using the Lagrange Multiplier
@@ -710,6 +720,48 @@ hand and foot constrained to be equal).](images/cheburashka-biharmonic-leq.jpg)
 
 ## Quadratic Programming
 
+We can generalize the quadratic optimization in the previous section even more
+by allowing inequality constraints. Specifically box constraints (lower and
+upper bounds):
+
+ $\mathbf{l} \le \mathbf{z} \le \mathbf{u},$
+
+where $\mathbf{l},\mathbf{u}$ are $n \times 1$ vectors of lower and upper
+bounds
+and general linear inequality constraints:
+
+ $\mathbf{A}_{ieq} \mathbf{z} \le \mathbf{B}_{ieq},$
+
+where $\mathbf{A}_{ieq}$ is a $k \times n$ matrix of linear coefficients and
+$\mathbf{B}_{ieq}$ is a $k \times 1$ matrix of constraint right-hand sides.
+
+Again, we are overly general as the box constraints could be written as
+rows of the linear inequality constraints, but bounds appear frequently enough
+to merit a dedicated api.
+
+Libigl implements its own active set routine for solving _quadratric programs_
+(QPs). This algorithm works by iteratively "activating" violated inequality
+constraints by enforcing them as equalities and "deactivating" constraints
+which are no longer needed.
+
+After deciding which constraints are active each iteration simple reduces to a
+quadratic minimization subject to linear _equality_ constraints, and the method
+from the previous section is invoked. This is repeated until convergence.
+
+Currently the implementation is efficient for box constraints and sparse
+non-overlapping linear inequality constraints.
+
+Unlike alternative interior-point methods, the active set method benefits from
+a warm-start (initial guess for the solution vector $\mathbf{z}$).
+
+```cpp
+igl::active_set_params as;
+// Z is optional initial guess and output
+igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
+```
+
+![The example `QuadraticProgramming` uses an active set solver to optimize
+discrete biharmonic kernels at multiple scales [#rustamov_2011][].](images/cheburashka-multiscale-biharmonic-kernels.jpg)
 
 [#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
 "Discrete Differential-Geometry Operators for Triangulated
@@ -721,3 +773,4 @@ _Algorithms and Interfaces for Real-Time Deformation of 2D and 3D Shapes_,
 2013.
 [#kazhdan_2012]: Michael Kazhdan, Jake Solomon, Mirela Ben-Chen,
 "Can Mean-Curvature Flow Be Made Non-Singular," 2012.
+[#rustamov_2011]: Raid M. Rustamov, "Multiscale Biharmonic Kernels", 2011.