added bbw to mosek package and fixed bug in boundary conditions

Former-commit-id: a6677894489c5689d29373ea6c86ee972f9155ab

include/igl/boundary_conditions.cpp

@@ -6,6 +6,7 @@
 
 #include <vector>
 #include <map>
+#include <iostream>
 
 IGL_INLINE bool igl::boundary_conditions(
   const Eigen::MatrixXd & V  ,
@@ -41,9 +42,20 @@ IGL_INLINE bool igl::boundary_conditions(
     {
       // Find samples just on pos
       //Vec3 vi(V(i,0),V(i,1),V(i,2));
+      // EIGEN GOTCHA:
+      // double sqrd = (V.row(i)-pos).array().pow(2).sum();
+      // Must first store in temporary
+      VectorXd vi = V.row(i)
       double sqrd = (V.row(i)-pos).array().pow(2).sum();
       if(sqrd <= FLOAT_EPS)
       {
+        cout<<"sum((["<<
+          V(i,0)<<" "<<
+          V(i,1)<<" "<<
+          V(i,2)<<"] - ["<<
+          pos(0)<<" "<<
+          pos(1)<<" "<<
+          pos(2)<<"]).^2) = "<<sqrd<<endl;
         bci.push_back(i);
         bcj.push_back(p);
         bcv.push_back(1.0);

include/igl/list_to_matrix.cpp

@@ -17,8 +17,10 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Mat
   int m = V.size();
   if(m == 0)
   {
-    fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
-    return false;
+    //fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
+    //return false;
+    M.resize(0,0);
+    return true;
   }
   // number of columns
   int n = igl::min_size(V);
@@ -52,8 +54,10 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Mat & M)
   int m = V.size();
   if(m == 0)
   {
-    fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
-    return false;
+    //fprintf(stderr,"Error: list_to_matrix() list is empty()\n");
+    //return false;
+    M.resize(0,0);
+    return true;
   }
   // Resize output
   M.resize(m,1);

include/igl/mosek/bbw.cpp

@@ -0,0 +1,114 @@
+#define VERBOSE
+#include "bbw.h"
+
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/invert_diag.h>
+#include <igl/speye.h>
+#include <igl/slice_into.h>
+#include <igl/normalize_row_sums.h>
+#include <igl/verbose.h>
+
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Sparse>
+
+#include <iostream>
+#include <cstdio>
+
+
+igl::BBWData::BBWData():
+  partition_unity(false)
+{}
+
+void igl::BBWData::print()
+{
+  using namespace std;
+  cout<<"partition_unity: "<<partition_unity<<endl;
+  cout<<"W0=["<<endl<<W0<<endl<<"];"<<endl;
+}
+
+
+template <
+  typename DerivedV, 
+  typename DerivedEle, 
+  typename Derivedb,
+  typename Derivedbc, 
+  typename DerivedW>
+IGL_INLINE bool igl::bbw(
+  const Eigen::MatrixBase<DerivedV> & V, 
+  const Eigen::MatrixBase<DerivedEle> & Ele, 
+  const Eigen::MatrixBase<Derivedb> & b, 
+  const Eigen::MatrixBase<Derivedbc> & bc, 
+  igl::BBWData & data,
+  Eigen::MatrixBase<DerivedW> & W
+  )
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+
+  // number of domain vertices
+  int n = V.rows();
+  // number of handles
+  int m = bc.cols();
+
+  SparseMatrix<typename DerivedW::Scalar> L;
+  cotmatrix(V,Ele,L);
+  MassMatrixType mmtype = MASSMATRIX_VORONOI;
+  if(Ele.cols() == 4)
+  {
+    mmtype = MASSMATRIX_BARYCENTRIC;
+  }
+  SparseMatrix<typename DerivedW::Scalar> M;
+  SparseMatrix<typename DerivedW::Scalar> Mi;
+  massmatrix(V,Ele,mmtype,M);
+
+  invert_diag(M,Mi);
+
+  // Biharmonic operator
+  SparseMatrix<typename DerivedW::Scalar> Q = L.transpose() * Mi * L;
+
+  W.derived().resize(n,m);
+  if(data.partition_unity)
+  {
+    // Not yet implemented
+    assert(false);
+  }else
+  {
+    // No linear terms
+    VectorXd c = VectorXd::Zero(n);
+    // No linear constraints
+    SparseMatrix<typename DerivedW::Scalar> A(0,n);
+    VectorXd uc(0,1);
+    VectorXd lc(0,1);
+    // 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++)
+    {
+      verbose("\n^%s: Computing weight for handle %d out of %d.\n\n",
+        __FUNCTION__,i,m);
+      // impose boundary conditions
+      VectorXd bci = bc.col(i);
+      slice_into(bci,b,ux);
+      slice_into(bci,b,lx);
+      VectorXd Wi;
+      bool r = igl::mosek_quadprog(Q,c,0,A,lc,uc,lx,ux,data.mosek_data,Wi);
+      if(!r)
+      {
+        return false;
+      }
+      W.col(i) = Wi;
+    }
+    // Need to normalize
+    igl::normalize_row_sums(W,W); 
+  }
+
+  return true;
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template bool igl::bbw<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 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::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::BBWData&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/mosek/bbw.h

@@ -0,0 +1,66 @@
+#ifndef BBW_H
+#define BBW_H
+#include "../igl_inline.h"
+
+#include <Eigen/Dense>
+#include "mosek_quadprog.h"
+
+namespace igl
+{
+  // Container for BBW computation related data and flags
+  class BBWData
+  {
+    public:
+      // Enforce partition of unity during optimization (optimize all weight
+      // simultaneously)
+      bool partition_unity;
+      // TODO: Is it safe if this is a reference?
+      // Initial guess
+      Eigen::MatrixXd W0;
+      // TODO: Mosek options
+      igl::MosekData mosek_data;
+      // TODO: Active set options
+    public:
+      BBWData();
+      // Print current state of object
+      void print();
+  };
+
+  // Compute Bounded Biharmonic Weights on a given domain (V,Ele) with a given
+  // set of boundary conditions
+  //
+  // Templates
+  //   DerivedV  derived type of eigen matrix for V (e.g. MatirxXd)
+  //   DerivedF  derived type of eigen matrix for F (e.g. MatirxXi)
+  //   Derivedb  derived type of eigen matrix for b (e.g. VectorXi)
+  //   Derivedbc  derived type of eigen matrix for bc (e.g. MatirxXd)
+  //   DerivedW  derived type of eigen matrix for W (e.g. MatirxXd)
+  // Inputs:
+  //   V  #V by dim vertex positions
+  //   Ele  #Elements by simplex-size list of element indices
+  //   b  #b boundary indices into V
+  //   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
+  // Returns true on success, false on failure
+  template <
+    typename DerivedV, 
+    typename DerivedEle, 
+    typename Derivedb,
+    typename Derivedbc, 
+    typename DerivedW>
+  IGL_INLINE bool bbw(
+    const Eigen::MatrixBase<DerivedV> & V, 
+    const Eigen::MatrixBase<DerivedEle> & Ele, 
+    const Eigen::MatrixBase<Derivedb> & b, 
+    const Eigen::MatrixBase<Derivedbc> & bc, 
+    BBWData & data,
+    Eigen::MatrixBase<DerivedW> & W);
+}
+  
+#ifdef IGL_HEADER_ONLY
+#  include "bbw.cpp"
+#endif
+
+#endif