better separation of bbw solvers, openmp for active set and unproject eigen wrapper

Former-commit-id: 00152cbe796abcfed1703eb859351b942614e6a7

examples/bbw/main.cpp

@@ -11,6 +11,7 @@
 #include <igl/mosek/bbw.h>
 #include <igl/writeDMAT.h>
 #include <igl/writeMESH.h>
+#include <igl/normalize_row_sums.h>
 
 #include <Eigen/Dense>
 
@@ -255,6 +256,8 @@ int main(int argc, char * argv[])
   {
     return 1;
   }
+  // Normalize weights to sum to one
+  normalize_row_sums(W,W);
   // Save output
   save_output(argv[1],argv[2],V,F,VV,TT,FF,W);
   return 0;

include/igl/min_quad_with_fixed.cpp

@@ -513,5 +513,6 @@ IGL_INLINE bool igl::min_quad_with_fixed_solve(
 // Explicit template specialization
 template bool igl::min_quad_with_fixed_solve<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::Matrix<double, -1, 1, 0, -1, 1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template bool igl::min_quad_with_fixed_precompute<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int> const&, bool, igl::min_quad_with_fixed_data<double>&);
+template bool igl::min_quad_with_fixed_solve<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::Matrix<double, -1, -1, 0, -1, -1> >(igl::min_quad_with_fixed_data<double> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif
 

include/igl/mosek/Makefile

@@ -5,8 +5,9 @@ all: libiglmosek
 debug: libiglmosek
 
 include ../../../Makefile.conf
-all: CFLAGS += -O3 -DNDEBUG
-debug: CFLAGS += -g -Wall -Werror
+all: OPTFLAGS += -O3 -DNDEBUG $(OPENMP)
+debug: OPTFLAGS += -g -Wall -Werror
+CFLAGS += $(OPTFLAGS)
 
 .PHONY: libiglmosek
 libiglmosek: obj ../../../lib/libiglmosek.a

include/igl/mosek/bbw.cpp

@@ -8,6 +8,7 @@
 #include <igl/slice_into.h>
 #include <igl/normalize_row_sums.h>
 #include <igl/verbose.h>
+#include <igl/min_quad_with_fixed.h>
 
 #define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
 #include <Eigen/Sparse>
@@ -85,19 +86,36 @@ IGL_INLINE bool igl::bbw(
     // Upper and lower box constraints (Constant bounds)
     VectorXd ux = VectorXd::Ones(n);
     VectorXd lx = VectorXd::Zero(n);
-    // Loop over handles
-    for(int i = 0;i<m;i++)
+    active_set_params eff_params = data.active_set_params;
+    switch(data.qp_solver)
     {
-      verbose("\n^%s: Computing weight for handle %d out of %d.\n\n",
-        __FUNCTION__,i+1,m);
-      VectorXd bci = bc.col(i);
-      VectorXd Wi;
-      switch(data.qp_solver)
+      case QP_SOLVER_IGL_ACTIVE_SET:
       {
-        case QP_SOLVER_IGL_ACTIVE_SET:
+        verbose("\n^%s: Computing initial weights for %d handle%s.\n\n",
+          __FUNCTION__,m,(m!=1?"s":""));
+        min_quad_with_fixed_data<typename DerivedW::Scalar > mqwf;
+        min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
+        min_quad_with_fixed_solve(mqwf,c,bc,Beq,W);
+        // decrement
+        eff_params.max_iter--;
+        bool error = false;
+        // Loop over handles
+#pragma omp parallel for
+        for(int i = 0;i<m;i++)
         {
+          // Quicker exit for openmp
+          if(error)
+          {
+            continue;
+          }
+          verbose("\n^%s: Computing weight for handle %d out of %d.\n\n",
+              __FUNCTION__,i+1,m);
+          VectorXd bci = bc.col(i);
+          VectorXd Wi;
+          // use initial guess
+          Wi = W.col(i);
           SolverStatus ret = active_set(
-            Q,c,b,bci,Aeq,Beq,Aieq,Bieq,lx,ux,data.active_set_params,Wi);
+              Q,c,b,bci,Aeq,Beq,Aieq,Bieq,lx,ux,eff_params,Wi);
           switch(ret)
           {
             case SOLVER_STATUS_CONVERGED:
@@ -108,12 +126,25 @@ IGL_INLINE bool igl::bbw(
             case SOLVER_STATUS_ERROR:
             default:
               cout<<"active_set error."<<endl;
-              return false;
+              error = true;
           }
-          break;
+          W.col(i) = Wi;
         }
-        case QP_SOLVER_MOSEK:
+        if(error)
+        {
+          return false;
+        }
+        break;
+      }
+      case QP_SOLVER_MOSEK:
+      {
+        // Loop over handles
+        for(int i = 0;i<m;i++)
         {
+          verbose("\n^%s: Computing weight for handle %d out of %d.\n\n",
+              __FUNCTION__,i+1,m);
+          VectorXd bci = bc.col(i);
+          VectorXd Wi;
           // impose boundary conditions via bounds
           slice_into(bci,b,ux);
           slice_into(bci,b,lx);
@@ -122,24 +153,24 @@ IGL_INLINE bool igl::bbw(
           {
             return false;
           }
-          break;
-        }
-        default:
-        {
-          assert(false && "Unknown qp_solver");
-          return false;
+          W.col(i) = Wi;
         }
+        break;
+      }
+      default:
+      {
+        assert(false && "Unknown qp_solver");
+        return false;
       }
-      W.col(i) = Wi;
     }
+#ifndef NDEBUG
     const double min_rowsum = W.rowwise().sum().array().abs().minCoeff();
     if(min_rowsum < 0.1)
     {
-      cerr<<"bbw.cpp: Warning, row sum is very low. Consider more iterations "
-        "or enforcing partition of unity."<<endl;
+      cerr<<"bbw.cpp: Warning, minimum row sum is very low. Consider more "
+        "active set iterations or enforcing partition of unity."<<endl;
     }
-    // Need to normalize
-    igl::normalize_row_sums(W,W); 
+#endif
   }
 
   return true;

include/igl/mosek/bbw.h

@@ -57,7 +57,8 @@ namespace igl
   //   bc #b by #W list of boundary values
   //   data  object containing options, intial guess --> solution and results
   // Outputs:
-  //   W  #V by #W list of weights
+  //   W  #V by #W list of *unnormalized* weights to normalize use 
+  //    igl::normalize_row_sums(W,W); 
   // Returns true on success, false on failure
   template <
     typename DerivedV, 

include/igl/unproject.cpp

@@ -20,4 +20,21 @@ IGL_INLINE int igl::unproject(
   glGetIntegerv(GL_VIEWPORT, VP);
   return gluUnProject(winX,winY,winZ,MV,P,VP,objX,objY,objZ);
 }
+
+template <typename Derivedwin, typename Derivedobj>
+IGL_INLINE int igl::unproject(
+  const Eigen::PlainObjectBase<Derivedwin> & win,
+  Eigen::PlainObjectBase<Derivedobj> & obj)
+{
+  return unproject(win(0),win(1),win(2),
+      &obj.data()[0],
+      &obj.data()[1],
+      &obj.data()[2]);
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template instanciation
+template int igl::unproject<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
+#endif
+
 #endif

include/igl/unproject.h

@@ -1,6 +1,7 @@
 #ifndef IGL_UNPROJECT_H
 #define IGL_UNPROJECT_H
 #include "igl_inline.h"
+#include <Eigen/Core>
 namespace igl
 {
   // Wrapper for gluUnproject that uses the current GL_MODELVIEW_MATRIX,
@@ -17,6 +18,10 @@ namespace igl
     double* objX,
     double* objY,
     double* objZ);
+  template <typename Derivedwin, typename Derivedobj>
+  IGL_INLINE int unproject(
+    const Eigen::PlainObjectBase<Derivedwin> & win,
+    Eigen::PlainObjectBase<Derivedobj> & obj);
 }
 
 #ifdef IGL_HEADER_ONLY