better tetgen, matlab workspace is header only, bbw with active set option, eigen project(),

Former-commit-id: e093344745305205d2e60cd42509b93a80c05c3b

include/igl/active_set.cpp

@@ -36,6 +36,7 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   Eigen::PlainObjectBase<DerivedZ> & Z
   )
 {
+#define ACTIVE_SET_CPP_DEBUG
   using namespace igl;
   using namespace Eigen;
   using namespace std;
@@ -105,8 +106,10 @@ IGL_INLINE igl::SolverStatus igl::active_set(
   int iter = 0;
   while(true)
   {
-    //cout<<iter<<":"<<endl;
-    //cout<<"  pre"<<endl;
+#ifdef ACTIVE_SET_CPP_DEBUG
+    cout<<iter<<":"<<endl;
+    cout<<"  pre"<<endl;
+#endif
     // FIND BREACHES OF CONSTRAINTS
     int new_as_lx = 0;
     int new_as_ux = 0;
@@ -138,10 +141,14 @@ IGL_INLINE igl::SolverStatus igl::active_set(
           as_ieq(a) = TRUE;
         }
       }
-      //cout<<"new_as_lx: "<<new_as_lx<<endl;
-      //cout<<"new_as_ux: "<<new_as_ux<<endl;
+#ifdef ACTIVE_SET_CPP_DEBUG
+      cout<<"  new_as_lx: "<<new_as_lx<<endl;
+      cout<<"  new_as_ux: "<<new_as_ux<<endl;
+#endif
       const double diff = (Z-old_Z).squaredNorm();
-      //cout<<"diff: "<<diff<<endl;
+#ifdef ACTIVE_SET_CPP_DEBUG
+      cout<<"diff: "<<diff<<endl;
+#endif
       if(diff < params.solution_diff_threshold)
       {
         ret = SOLVER_STATUS_CONVERGED;
@@ -241,7 +248,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
     }
 #endif
     
-    //cout<<"  min_quad_with_fixed_precompute"<<endl;
+#ifdef ACTIVE_SET_CPP_DEBUG
+    cout<<"  min_quad_with_fixed_precompute"<<endl;
+#endif
     if(!min_quad_with_fixed_precompute(A,known_i,Aeq_i,params.Auu_pd,data))
     {
       cerr<<"Error: min_quad_with_fixed precomputation failed."<<endl;
@@ -253,7 +262,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       ret = SOLVER_STATUS_ERROR;
       break;
     }
-    //cout<<"  min_quad_with_fixed_solve"<<endl;
+#ifdef ACTIVE_SET_CPP_DEBUG
+    cout<<"  min_quad_with_fixed_solve"<<endl;
+#endif
     Eigen::PlainObjectBase<DerivedZ> sol;
     if(!min_quad_with_fixed_solve(data,B,Y_i,Beq_i,Z,sol))
     {
@@ -261,7 +272,9 @@ IGL_INLINE igl::SolverStatus igl::active_set(
       ret = SOLVER_STATUS_ERROR;
       break;
     }
-    //cout<<"  post"<<endl;
+#ifdef ACTIVE_SET_CPP_DEBUG
+    cout<<"  post"<<endl;
+#endif
 
     // Compute Lagrange multiplier values for known_i
     // This needs to be adjusted slightly if A is not symmetric

include/igl/matlab/MatlabWorkspace.cpp

@@ -1,342 +0,0 @@
-#include "igl/matlab/MatlabWorkspace.h"
-
-// IGL
-#include "igl/list_to_matrix.h"
-
-// MATLAB
-#include "mat.h"
-
-// STL
-#include <iostream>
-#include <algorithm>
-#include <vector>
-
-IGL_INLINE igl::MatlabWorkspace::MatlabWorkspace()
-{
-}
-
-IGL_INLINE igl::MatlabWorkspace::~MatlabWorkspace()
-{
-  // clean up data
-  clear();
-}
-
-IGL_INLINE void igl::MatlabWorkspace::clear()
-{
-  for_each(data.begin(),data.end(),&mxDestroyArray);
-  data.clear();
-  names.clear();
-}
-
-IGL_INLINE bool igl::MatlabWorkspace::write(const std::string & path) const
-{
-  using namespace std;
-  MATFile * mat_file = matOpen(path.c_str(), "w");
-  assert(names.size() == data.size());
-  // loop over names and data
-  for(int i = 0;i < (int)names.size(); i++)
-  {
-    // Put variable as LOCAL variable
-    int status = matPutVariable(mat_file,names[i].c_str(), data[i]);
-    if(status != 0) 
-    {
-      cerr<<"^MatlabWorkspace::save Error: matPutVariable ("<<names[i]<<
-        ") failed"<<endl;
-      return false;
-    } 
-  }
-  if(matClose(mat_file) != 0)
-  {
-    fprintf(stderr,"Error closing file %s\n",path.c_str());
-    return false;
-  }
-  return true;
-}
-
-IGL_INLINE bool igl::MatlabWorkspace::read(const std::string & path)
-{
-  using namespace std;
-
-  MATFile * mat_file;
-
-  mat_file = matOpen(path.c_str(), "r");
-  if (mat_file == NULL) 
-  {
-    cerr<<"Error: failed to open "<<path<<endl;
-    return false;
-  }
-
-  int ndir;
-  const char ** dir = (const char **)matGetDir(mat_file, &ndir);
-  if (dir == NULL) {
-    cerr<<"Error reading directory of file "<< path<<endl;
-    return false;
-  }
-  mxFree(dir);
-
-  // Must close and reopen
-  if(matClose(mat_file) != 0)
-  {
-    cerr<<"Error: failed to close file "<<path<<endl;
-    return false;
-  }
-  mat_file = matOpen(path.c_str(), "r");
-  if (mat_file == NULL) 
-  {
-    cerr<<"Error: failed to open "<<path<<endl;
-    return false;
-  }
-  
-
-  /* Read in each array. */
-  for (int i=0; i<ndir; i++) 
-  {
-    const char * name;
-    mxArray * mx_data = matGetNextVariable(mat_file, &name);
-    if (mx_data == NULL) 
-    {
-      cerr<<"Error: matGetNextVariable failed in "<<path<<endl;
-      return false;
-    } 
-    const int dims = mxGetNumberOfDimensions(mx_data);
-    assert(dims == 2);
-    if(dims != 2)
-    {
-      fprintf(stderr,"Variable '%s' has %d ≠ 2 dimensions. Skipping\n",
-          name,dims);
-      mxDestroyArray(mx_data);
-      continue;
-    }
-    // don't destroy
-    names.push_back(name);
-    data.push_back(mx_data);
-  }
-
-  if(matClose(mat_file) != 0)
-  {
-    cerr<<"Error: failed to close file "<<path<<endl;
-    return false;
-  }
-
-  return true;
-}
-
-// Treat everything as a double
-template <typename DerivedM>
-IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
-  const Eigen::PlainObjectBase<DerivedM>& M,
-  const std::string & name)
-{
-  using namespace std;
-  const int m = M.rows();
-  const int n = M.cols();
-  mxArray * mx_data = mxCreateDoubleMatrix(m,n,mxREAL);
-  data.push_back(mx_data);
-  names.push_back(name);
-  // Copy data immediately
-  // Q: Won't this be incorrect for integers?
-  copy(M.data(),M.data()+M.size(),mxGetPr(mx_data));
-  return *this;
-}
-
-// Treat everything as a double
-template <typename MT>
-IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
-  const Eigen::SparseMatrix<MT>& M,
-  const std::string & name)
-{
-  using namespace std;
-  const int m = M.rows();
-  const int n = M.cols();
-  // THIS WILL NOT WORK FOR ROW-MAJOR
-  assert(n==M.outerSize());
-  const int nzmax = M.nonZeros();
-  mxArray * mx_data = mxCreateSparse(m, n, nzmax, mxREAL);
-  data.push_back(mx_data);
-  names.push_back(name);
-  // Copy data immediately
-  double * pr = mxGetPr(mx_data);
-  mwIndex * ir = mxGetIr(mx_data);
-  mwIndex * jc = mxGetJc(mx_data);
-
-  // Iterate over outside
-  int k = 0;
-  for(int j=0; j<M.outerSize();j++)
-  {
-    jc[j] = k;
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<MT>::InnerIterator it (M,j); it; ++it)
-    {
-      pr[k] = it.value();
-      ir[k] = it.row();
-      k++;
-    }
-  }
-  jc[M.outerSize()] = k;
-
-  return *this;
-}
-
-template <typename ScalarM>
-IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
-  const std::vector<std::vector<ScalarM> > & vM,
-  const std::string & name)
-{
-  Eigen::MatrixXd M;
-  list_to_matrix(vM,M);
-  return this->save(M,name);
-}
-
-template <typename ScalarV>
-IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
-  const std::vector<ScalarV> & vV,
-  const std::string & name)
-{
-  Eigen::MatrixXd V;
-  list_to_matrix(vV,V);
-  return this->save(V,name);
-}
-
-template <typename DerivedM>
-IGL_INLINE igl::MatlabWorkspace& 
-  igl::MatlabWorkspace::save_index(
-    const Eigen::PlainObjectBase<DerivedM>& M,
-    const std::string & name)
-{
-  DerivedM Mp1 = M;
-  Mp1.array() += 1;
-  return this->save(Mp1,name);
-}
-
-template <typename ScalarM>
-IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save_index(
-  const std::vector<std::vector<ScalarM> > & vM,
-  const std::string & name)
-{
-  Eigen::MatrixXd M;
-  list_to_matrix(vM,M);
-  return this->save_index(M,name);
-}
-
-template <typename ScalarV>
-IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save_index(
-  const std::vector<ScalarV> & vV,
-  const std::string & name)
-{
-  Eigen::MatrixXd V;
-  list_to_matrix(vV,V);
-  return this->save_index(V,name);
-}
-
-template <typename DerivedM>
-IGL_INLINE bool igl::MatlabWorkspace::find( 
-  const std::string & name,
-  Eigen::PlainObjectBase<DerivedM>& M)
-{
-  using namespace std;
-  const int i = std::find(names.begin(), names.end(), name)-names.begin();
-  if(i>=(int)names.size())
-  {
-    return false;
-  }
-  assert(i<=(int)data.size());
-  mxArray * mx_data = data[i];
-  assert(!mxIsSparse(mx_data));
-  assert(mxGetNumberOfDimensions(mx_data) == 2);
-  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
-  const int m = mxGetM(mx_data);
-  const int n = mxGetN(mx_data);
-
-  // Handle vectors
-  if(DerivedM::IsVectorAtCompileTime)
-  {
-    assert(m==1 || n==1 || (m==0 && n==0));
-    M.resize(m*n);
-  }else
-  {
-    M.resize(m,n);
-  }
-  copy(
-    mxGetPr(mx_data),
-    mxGetPr(mx_data)+mxGetNumberOfElements(mx_data),
-    M.data());
-  return true;
-}
-
-template <typename MT>
-IGL_INLINE bool igl::MatlabWorkspace::find( 
-  const std::string & name,
-  Eigen::SparseMatrix<MT>& M)
-{
-  using namespace std;
-  using namespace Eigen;
-  const int i = std::find(names.begin(), names.end(), name)-names.begin();
-  if(i>=(int)names.size())
-  {
-    return false;
-  }
-  assert(i<=(int)data.size());
-  mxArray * mx_data = data[i];
-  // Handle boring case where matrix is actually an empty dense matrix
-  if(mxGetNumberOfElements(mx_data) == 0)
-  {
-    M.resize(0,0);
-    return true;
-  }
-  assert(mxIsSparse(mx_data));
-  assert(mxGetNumberOfDimensions(mx_data) == 2);
-  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
-  const int m = mxGetM(mx_data);
-  const int n = mxGetN(mx_data);
-
-  // Copy data immediately
-  double * pr = mxGetPr(mx_data);
-  mwIndex * ir = mxGetIr(mx_data);
-  mwIndex * jc = mxGetJc(mx_data);
-
-  vector<Triplet<MT> > MIJV;
-  MIJV.reserve(mxGetNumberOfElements(mx_data));
-
-  // Iterate over outside
-  int k = 0;
-  for(int j=0; j<n;j++)
-  {
-    // Iterate over inside
-    while(k<(int)jc[j+1])
-    {
-      //cout<<ir[k]<<" "<<j<<" "<<pr[k]<<endl;
-      assert((int)ir[k]<m);
-      assert((int)j<n);
-      MIJV.push_back(Triplet<MT >(ir[k],j,pr[k]));
-      k++;
-    }
-  }
-  M.resize(m,n);
-  M.setFromTriplets(MIJV.begin(),MIJV.end());
-  return true;
-
-}
-
-template <typename DerivedM>
-IGL_INLINE bool igl::MatlabWorkspace::find_index( 
-  const std::string & name,
-  Eigen::PlainObjectBase<DerivedM>& M)
-{
-  if(!find(name,M))
-  {
-    return false;
-  }
-  M.array() -= 1;
-  return true;
-}
-
-
-//template <typename Data>
-//bool igl::MatlabWorkspace::save(const Data & M, const std::string & name)
-//{
-//  using namespace std;
-//  // If I don't know the type then I can't save it
-//  cerr<<"^MatlabWorkspace::save Error: Unknown data type. "<<
-//    name<<" not saved."<<endl;
-//  return false;
-//}

include/igl/matlab/MatlabWorkspace.h

@@ -120,10 +120,351 @@ namespace igl
   };
 }
 
+// Implementation
+
 // Be sure that this is not compiled into libigl.a
-#ifdef IGL_HEADER_ONLY
-#  include "MatlabWorkspace.cpp"
-#endif
+// http://stackoverflow.com/a/3318993/148668
+
+// IGL
+#include "igl/list_to_matrix.h"
+
+// MATLAB
+#include "mat.h"
+
+// STL
+#include <iostream>
+#include <algorithm>
+#include <vector>
+
+IGL_INLINE igl::MatlabWorkspace::MatlabWorkspace()
+{
+}
+
+IGL_INLINE igl::MatlabWorkspace::~MatlabWorkspace()
+{
+  // clean up data
+  clear();
+}
+
+IGL_INLINE void igl::MatlabWorkspace::clear()
+{
+  for_each(data.begin(),data.end(),&mxDestroyArray);
+  data.clear();
+  names.clear();
+}
+
+IGL_INLINE bool igl::MatlabWorkspace::write(const std::string & path) const
+{
+  using namespace std;
+  MATFile * mat_file = matOpen(path.c_str(), "w");
+  assert(names.size() == data.size());
+  // loop over names and data
+  for(int i = 0;i < (int)names.size(); i++)
+  {
+    // Put variable as LOCAL variable
+    int status = matPutVariable(mat_file,names[i].c_str(), data[i]);
+    if(status != 0) 
+    {
+      cerr<<"^MatlabWorkspace::save Error: matPutVariable ("<<names[i]<<
+        ") failed"<<endl;
+      return false;
+    } 
+  }
+  if(matClose(mat_file) != 0)
+  {
+    fprintf(stderr,"Error closing file %s\n",path.c_str());
+    return false;
+  }
+  return true;
+}
+
+IGL_INLINE bool igl::MatlabWorkspace::read(const std::string & path)
+{
+  using namespace std;
+
+  MATFile * mat_file;
+
+  mat_file = matOpen(path.c_str(), "r");
+  if (mat_file == NULL) 
+  {
+    cerr<<"Error: failed to open "<<path<<endl;
+    return false;
+  }
+
+  int ndir;
+  const char ** dir = (const char **)matGetDir(mat_file, &ndir);
+  if (dir == NULL) {
+    cerr<<"Error reading directory of file "<< path<<endl;
+    return false;
+  }
+  mxFree(dir);
+
+  // Must close and reopen
+  if(matClose(mat_file) != 0)
+  {
+    cerr<<"Error: failed to close file "<<path<<endl;
+    return false;
+  }
+  mat_file = matOpen(path.c_str(), "r");
+  if (mat_file == NULL) 
+  {
+    cerr<<"Error: failed to open "<<path<<endl;
+    return false;
+  }
+  
+
+  /* Read in each array. */
+  for (int i=0; i<ndir; i++) 
+  {
+    const char * name;
+    mxArray * mx_data = matGetNextVariable(mat_file, &name);
+    if (mx_data == NULL) 
+    {
+      cerr<<"Error: matGetNextVariable failed in "<<path<<endl;
+      return false;
+    } 
+    const int dims = mxGetNumberOfDimensions(mx_data);
+    assert(dims == 2);
+    if(dims != 2)
+    {
+      fprintf(stderr,"Variable '%s' has %d ≠ 2 dimensions. Skipping\n",
+          name,dims);
+      mxDestroyArray(mx_data);
+      continue;
+    }
+    // don't destroy
+    names.push_back(name);
+    data.push_back(mx_data);
+  }
+
+  if(matClose(mat_file) != 0)
+  {
+    cerr<<"Error: failed to close file "<<path<<endl;
+    return false;
+  }
+
+  return true;
+}
+
+// Treat everything as a double
+template <typename DerivedM>
+IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
+  const Eigen::PlainObjectBase<DerivedM>& M,
+  const std::string & name)
+{
+  using namespace std;
+  const int m = M.rows();
+  const int n = M.cols();
+  mxArray * mx_data = mxCreateDoubleMatrix(m,n,mxREAL);
+  data.push_back(mx_data);
+  names.push_back(name);
+  // Copy data immediately
+  // Q: Won't this be incorrect for integers?
+  copy(M.data(),M.data()+M.size(),mxGetPr(mx_data));
+  return *this;
+}
+
+// Treat everything as a double
+template <typename MT>
+IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
+  const Eigen::SparseMatrix<MT>& M,
+  const std::string & name)
+{
+  using namespace std;
+  const int m = M.rows();
+  const int n = M.cols();
+  // THIS WILL NOT WORK FOR ROW-MAJOR
+  assert(n==M.outerSize());
+  const int nzmax = M.nonZeros();
+  mxArray * mx_data = mxCreateSparse(m, n, nzmax, mxREAL);
+  data.push_back(mx_data);
+  names.push_back(name);
+  // Copy data immediately
+  double * pr = mxGetPr(mx_data);
+  mwIndex * ir = mxGetIr(mx_data);
+  mwIndex * jc = mxGetJc(mx_data);
+
+  // Iterate over outside
+  int k = 0;
+  for(int j=0; j<M.outerSize();j++)
+  {
+    jc[j] = k;
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<MT>::InnerIterator it (M,j); it; ++it)
+    {
+      pr[k] = it.value();
+      ir[k] = it.row();
+      k++;
+    }
+  }
+  jc[M.outerSize()] = k;
+
+  return *this;
+}
+
+template <typename ScalarM>
+IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
+  const std::vector<std::vector<ScalarM> > & vM,
+  const std::string & name)
+{
+  Eigen::MatrixXd M;
+  list_to_matrix(vM,M);
+  return this->save(M,name);
+}
+
+template <typename ScalarV>
+IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save(
+  const std::vector<ScalarV> & vV,
+  const std::string & name)
+{
+  Eigen::MatrixXd V;
+  list_to_matrix(vV,V);
+  return this->save(V,name);
+}
+
+template <typename DerivedM>
+IGL_INLINE igl::MatlabWorkspace& 
+  igl::MatlabWorkspace::save_index(
+    const Eigen::PlainObjectBase<DerivedM>& M,
+    const std::string & name)
+{
+  DerivedM Mp1 = M;
+  Mp1.array() += 1;
+  return this->save(Mp1,name);
+}
+
+template <typename ScalarM>
+IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save_index(
+  const std::vector<std::vector<ScalarM> > & vM,
+  const std::string & name)
+{
+  Eigen::MatrixXd M;
+  list_to_matrix(vM,M);
+  return this->save_index(M,name);
+}
+
+template <typename ScalarV>
+IGL_INLINE igl::MatlabWorkspace& igl::MatlabWorkspace::save_index(
+  const std::vector<ScalarV> & vV,
+  const std::string & name)
+{
+  Eigen::MatrixXd V;
+  list_to_matrix(vV,V);
+  return this->save_index(V,name);
+}
+
+template <typename DerivedM>
+IGL_INLINE bool igl::MatlabWorkspace::find( 
+  const std::string & name,
+  Eigen::PlainObjectBase<DerivedM>& M)
+{
+  using namespace std;
+  const int i = std::find(names.begin(), names.end(), name)-names.begin();
+  if(i>=(int)names.size())
+  {
+    return false;
+  }
+  assert(i<=(int)data.size());
+  mxArray * mx_data = data[i];
+  assert(!mxIsSparse(mx_data));
+  assert(mxGetNumberOfDimensions(mx_data) == 2);
+  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
+  const int m = mxGetM(mx_data);
+  const int n = mxGetN(mx_data);
+
+  // Handle vectors
+  if(DerivedM::IsVectorAtCompileTime)
+  {
+    assert(m==1 || n==1 || (m==0 && n==0));
+    M.resize(m*n);
+  }else
+  {
+    M.resize(m,n);
+  }
+  copy(
+    mxGetPr(mx_data),
+    mxGetPr(mx_data)+mxGetNumberOfElements(mx_data),
+    M.data());
+  return true;
+}
+
+template <typename MT>
+IGL_INLINE bool igl::MatlabWorkspace::find( 
+  const std::string & name,
+  Eigen::SparseMatrix<MT>& M)
+{
+  using namespace std;
+  using namespace Eigen;
+  const int i = std::find(names.begin(), names.end(), name)-names.begin();
+  if(i>=(int)names.size())
+  {
+    return false;
+  }
+  assert(i<=(int)data.size());
+  mxArray * mx_data = data[i];
+  // Handle boring case where matrix is actually an empty dense matrix
+  if(mxGetNumberOfElements(mx_data) == 0)
+  {
+    M.resize(0,0);
+    return true;
+  }
+  assert(mxIsSparse(mx_data));
+  assert(mxGetNumberOfDimensions(mx_data) == 2);
+  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
+  const int m = mxGetM(mx_data);
+  const int n = mxGetN(mx_data);
+
+  // Copy data immediately
+  double * pr = mxGetPr(mx_data);
+  mwIndex * ir = mxGetIr(mx_data);
+  mwIndex * jc = mxGetJc(mx_data);
+
+  vector<Triplet<MT> > MIJV;
+  MIJV.reserve(mxGetNumberOfElements(mx_data));
+
+  // Iterate over outside
+  int k = 0;
+  for(int j=0; j<n;j++)
+  {
+    // Iterate over inside
+    while(k<(int)jc[j+1])
+    {
+      //cout<<ir[k]<<" "<<j<<" "<<pr[k]<<endl;
+      assert((int)ir[k]<m);
+      assert((int)j<n);
+      MIJV.push_back(Triplet<MT >(ir[k],j,pr[k]));
+      k++;
+    }
+  }
+  M.resize(m,n);
+  M.setFromTriplets(MIJV.begin(),MIJV.end());
+  return true;
+
+}
+
+template <typename DerivedM>
+IGL_INLINE bool igl::MatlabWorkspace::find_index( 
+  const std::string & name,
+  Eigen::PlainObjectBase<DerivedM>& M)
+{
+  if(!find(name,M))
+  {
+    return false;
+  }
+  M.array() -= 1;
+  return true;
+}
+
+
+//template <typename Data>
+//bool igl::MatlabWorkspace::save(const Data & M, const std::string & name)
+//{
+//  using namespace std;
+//  // If I don't know the type then I can't save it
+//  cerr<<"^MatlabWorkspace::save Error: Unknown data type. "<<
+//    name<<" not saved."<<endl;
+//  return false;
+//}
 
 #endif
 

include/igl/min_quad_with_fixed.cpp

@@ -28,10 +28,13 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
   min_quad_with_fixed_data<T> & data
   )
 {
+#define MIN_QUAD_WITH_FIXED_CPP_DEBUG
   using namespace Eigen;
   using namespace std;
   using namespace igl;
-  //cout<<"    pre"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+  cout<<"    pre"<<endl;
+#endif
   // number of rows
   int n = A.rows();
   // cache problem size
@@ -100,7 +103,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     data.Auu_sym = is_symmetric(Auu,EPS<double>());
   }
 
-  //cout<<"    qr"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+  cout<<"    qr"<<endl;
+#endif
   // Determine number of linearly independent constraints
   int nc = 0;
   if(neq>0)
@@ -134,7 +139,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
 
   if(data.Aeq_li)
   {
-    //cout<<"    Aeq_li=true"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    Aeq_li=true"<<endl;
+#endif
     // Append lagrange multiplier quadratic terms
     SparseMatrix<T> new_A;
     SparseMatrix<T> AeqT = Aeq.transpose();
@@ -168,9 +175,14 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
 
     // Positive definite and no equality constraints (Postive definiteness
     // implies symmetric)
-    //cout<<"    factorize"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    factorize"<<endl;
+#endif
     if(data.Auu_pd && neq == 0)
     {
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    llt"<<endl;
+#endif
       data.llt.compute(Auu);
       switch(data.llt.info())
       {
@@ -186,6 +198,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
       data.solver_type = min_quad_with_fixed_data<T>::LLT;
     }else
     {
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    ldlt"<<endl;
+#endif
       // Either not PD or there are equality constraints
       SparseMatrix<T> NA;
       slice(new_A,data.unknown_lagrange,data.unknown_lagrange,NA);
@@ -210,6 +225,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
         data.solver_type = min_quad_with_fixed_data<T>::LDLT;
       }else
       {
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    lu"<<endl;
+#endif
         // Resort to LU
         // Bottleneck >1/2
         data.lu.compute(NA);
@@ -233,7 +251,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     }
   }else
   {
-    //cout<<"    Aeq_li=false"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    Aeq_li=false"<<endl;
+#endif
     data.neq = neq;
     const int nu = data.unknown.size();
     //cout<<"nu: "<<nu<<endl;
@@ -244,7 +264,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     AeqTR = data.AeqTQR.matrixR();
     // This shouldn't be necessary
     AeqTR.prune(0.0);
-    //cout<<"    matrixQ"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    matrixQ"<<endl;
+#endif
     // THIS IS ESSENTIALLY DENSE AND THIS IS BY FAR THE BOTTLENECK
     // http://forum.kde.org/viewtopic.php?f=74&t=117500
     AeqTQ = data.AeqTQR.matrixQ();
@@ -271,11 +293,15 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
     // Null space
     data.AeqTQ2 = AeqTQ.bottomRightCorner(nu,nu-nc);
     data.AeqTQ2T = data.AeqTQ2.transpose().eval();
-    //cout<<"    proj"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    proj"<<endl;
+#endif
     // Projected hessian
     SparseMatrix<T> QRAuu = data.AeqTQ2T * Auu * data.AeqTQ2;
     {
-      //cout<<"    factorize"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+      cout<<"    factorize"<<endl;
+#endif
       // QRAuu should always be PD
       data.llt.compute(QRAuu);
       switch(data.llt.info())
@@ -291,7 +317,9 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
       }
       data.solver_type = min_quad_with_fixed_data<T>::QR_LLT;
     }
-    //cout<<"    smash"<<endl;
+#ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
+    cout<<"    smash"<<endl;
+#endif
     // Known value multiplier
     SparseMatrix<T> Auk;
     slice(A,data.unknown,data.known,Auk);

include/igl/mosek/bbw.cpp

@@ -17,7 +17,8 @@
 
 
 igl::BBWData::BBWData():
-  partition_unity(false)
+  partition_unity(false),
+  qp_solver(QP_SOLVER_IGL_ACTIVE_SET)
 {}
 
 void igl::BBWData::print()
@@ -25,6 +26,7 @@ void igl::BBWData::print()
   using namespace std;
   cout<<"partition_unity: "<<partition_unity<<endl;
   cout<<"W0=["<<endl<<W0<<endl<<"];"<<endl;
+  cout<<"qp_solver: "<<QPSolverNames[qp_solver]<<endl;
 }
 
 
@@ -35,12 +37,12 @@ template <
   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, 
+  const Eigen::PlainObjectBase<DerivedV> & V, 
+  const Eigen::PlainObjectBase<DerivedEle> & Ele, 
+  const Eigen::PlainObjectBase<Derivedb> & b, 
+  const Eigen::PlainObjectBase<Derivedbc> & bc, 
   igl::BBWData & data,
-  Eigen::MatrixBase<DerivedW> & W
+  Eigen::PlainObjectBase<DerivedW> & W
   )
 {
   using namespace igl;
@@ -78,9 +80,8 @@ IGL_INLINE bool igl::bbw(
     // 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);
+    SparseMatrix<typename DerivedW::Scalar> A(0,n),Aeq(0,n),Aieq(0,n);
+    VectorXd uc(0,1),Beq(0,1),Bieq(0,1),lc(0,1);
     // Upper and lower box constraints (Constant bounds)
     VectorXd ux = VectorXd::Ones(n);
     VectorXd lx = VectorXd::Zero(n);
@@ -89,15 +90,45 @@ IGL_INLINE bool igl::bbw(
     {
       verbose("\n^%s: Computing weight for handle %d out of %d.\n\n",
         __FUNCTION__,i+1,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)
+      switch(data.qp_solver)
       {
-        return false;
+        case QP_SOLVER_IGL_ACTIVE_SET:
+        {
+          SolverStatus ret = active_set(
+            Q,c,b,bci,Aeq,Beq,Aieq,Bieq,lx,ux,data.active_set_params,Wi);
+          switch(ret)
+          {
+            case SOLVER_STATUS_CONVERGED:
+              break;
+            case SOLVER_STATUS_MAX_ITER:
+              cout<<"active_set: max iter without convergence."<<endl;
+              break;
+            case SOLVER_STATUS_ERROR:
+            default:
+              cout<<"active_set error."<<endl;
+              return false;
+          }
+          break;
+        }
+        case QP_SOLVER_MOSEK:
+        {
+          // impose boundary conditions via bounds
+          slice_into(bci,b,ux);
+          slice_into(bci,b,lx);
+          bool r = mosek_quadprog(Q,c,0,A,lc,uc,lx,ux,data.mosek_data,Wi);
+          if(!r)
+          {
+            return false;
+          }
+          break;
+        }
+        default:
+        {
+          assert(false && "Unknown qp_solver");
+          return false;
+        }
       }
       W.col(i) = Wi;
     }
@@ -110,5 +141,5 @@ IGL_INLINE bool igl::bbw(
 
 #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> >&);
+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::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, igl::BBWData&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/mosek/bbw.h

@@ -4,9 +4,21 @@
 
 #include <Eigen/Dense>
 #include "mosek_quadprog.h"
+#include <igl/active_set.h>
 
 namespace igl
 {
+  enum QPSolver
+  {
+    QP_SOLVER_IGL_ACTIVE_SET = 0,
+    QP_SOLVER_MOSEK = 1,
+    NUM_QP_SOLVERS = 2
+  };
+  const char * const QPSolverNames[NUM_QP_SOLVERS] =
+  {
+    "QP_SOLVER_IGL_ACTIVE_SET",
+    "QP_SOLVER_MOSEK"
+  };
   // Container for BBW computation related data and flags
   class BBWData
   {
@@ -20,6 +32,9 @@ namespace igl
       // TODO: Mosek options
       igl::MosekData mosek_data;
       // TODO: Active set options
+      igl::active_set_params active_set_params;
+      // Which solver
+      QPSolver qp_solver;
     public:
       BBWData();
       // Print current state of object
@@ -51,12 +66,12 @@ namespace igl
     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, 
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedEle> & Ele, 
+    const Eigen::PlainObjectBase<Derivedb> & b, 
+    const Eigen::PlainObjectBase<Derivedbc> & bc, 
     BBWData & data,
-    Eigen::MatrixBase<DerivedW> & W);
+    Eigen::PlainObjectBase<DerivedW> & W);
 }
   
 #ifdef IGL_HEADER_ONLY

include/igl/project.cpp

@@ -62,4 +62,18 @@ IGL_INLINE int igl::project(
 
   return gluProject(objX,objY,objZ,MV,P,VP,winX,winY,winZ);
 }
+
+template <typename Derivedobj, typename Derivedwin>
+IGL_INLINE int igl::project(
+  const Eigen::PlainObjectBase<Derivedobj> & obj,
+  Eigen::PlainObjectBase<Derivedwin> & win)
+{
+  return igl::project(obj(0),obj(1),obj(2),&win(0),&win(1),&win(2));
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template instanciations
+template int igl::project<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/project.h

@@ -1,6 +1,7 @@
 #ifndef IGL_PROJECT_H
 #define IGL_PROJECT_H
 #include "igl_inline.h"
+#include <Eigen/Dense>
 namespace igl
 {
   // Wrapper for gluProject that uses the current GL_MODELVIEW_MATRIX,
@@ -17,6 +18,11 @@ namespace igl
     double* winX,
     double* winY,
     double* winZ);
+  // Eigen wrapper
+  template <typename Derivedobj, typename Derivedwin>
+  IGL_INLINE int project(
+    const Eigen::PlainObjectBase<Derivedobj> & obj,
+    Eigen::PlainObjectBase<Derivedwin> & win);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/tetgen/mesh_with_skeleton.cpp

@@ -6,6 +6,9 @@
 #include <igl/writeOFF.h>
 
 #include <iostream>
+// Default settings pq2Y tell tetgen to mesh interior of triangle mesh and
+// to produce a graded tet mesh
+const static std::string DEFAULT_TETGEN_FLAGS = "pq2Y";
 
 bool igl::mesh_with_skeleton(
   const Eigen::MatrixXd & V,
@@ -15,6 +18,7 @@ bool igl::mesh_with_skeleton(
   const Eigen::MatrixXi & BE,
   const Eigen::MatrixXi & CE,
   const int samples_per_bone,
+  const std::string & tetgen_flags,
   Eigen::MatrixXd & VV,
   Eigen::MatrixXi & TT,
   Eigen::MatrixXi & FF)
@@ -22,6 +26,8 @@ bool igl::mesh_with_skeleton(
   using namespace Eigen;
   using namespace igl;
   using namespace std;
+  const string eff_tetgen_flags = 
+    (tetgen_flags.length() == 0?DEFAULT_TETGEN_FLAGS:tetgen_flags);
   // Collect all edges that need samples:
   MatrixXi BECE = cat(1,BE,CE);
   MatrixXd S;
@@ -38,11 +44,9 @@ bool igl::mesh_with_skeleton(
   //   * has consistent orientation
   //   * has no self-intersections
   //   * has no 0-volume pieces
-  // Default settings pq100 tell tetgen to mesh interior of triangle mesh and
-  // to produce a graded tet mesh
   //writeOFF("mesh_with_skeleton.off",VS,F);
   cerr<<"tetgen begin()"<<endl;
-  int status = tetrahedralize( VS,F,"pq100Y",VV,TT,FF);
+  int status = tetrahedralize( VS,F,eff_tetgen_flags,VV,TT,FF);
   cerr<<"tetgen end()"<<endl;
   if(FF.rows() != F.rows())
   {
@@ -75,3 +79,18 @@ bool igl::mesh_with_skeleton(
   return true;
 }
 
+bool igl::mesh_with_skeleton(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & C,
+  const Eigen::VectorXi & P,
+  const Eigen::MatrixXi & BE,
+  const Eigen::MatrixXi & CE,
+  const int samples_per_bone,
+  Eigen::MatrixXd & VV,
+  Eigen::MatrixXi & TT,
+  Eigen::MatrixXi & FF)
+{
+  return igl::mesh_with_skeleton(
+    V,F,C,P,BE,CE,samples_per_bone,DEFAULT_TETGEN_FLAGS,VV,TT,FF);
+}

include/igl/tetgen/mesh_with_skeleton.h

@@ -2,6 +2,7 @@
 #define IGL_MESH_WITH_SKELETON_H
 #include "../igl_inline.h"
 #include <Eigen/Dense>
+#include <string>
 
 namespace igl
 {
@@ -17,11 +18,26 @@ namespace igl
 //  BE #BE by 2 list of bone-edge indices
 //  CE #CE by 2 list of cage-edge indices
 //  samples_per_bone  #samples to add per bone
+//  tetgen_flags  flags to pass to tetgen {""-->"pq2Y"} otherwise you're on
+//    your own and it's your funeral if you pass nonsense flags
 // Outputs:
 //  VV  #VV by 3 list of tet-mesh vertex positions
 //  TT  #TT by 4 list of tetrahedra indices
 //  FF  #FF by 3 list of surface triangle indices
 // Returns true only on success
+bool mesh_with_skeleton(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & C,
+  const Eigen::VectorXi & /*P*/,
+  const Eigen::MatrixXi & BE,
+  const Eigen::MatrixXi & CE,
+  const int samples_per_bone,
+  const std::string & tetgen_flags,
+  Eigen::MatrixXd & VV,
+  Eigen::MatrixXi & TT,
+  Eigen::MatrixXi & FF);
+// Wrapper using default tetgen_flags
 bool mesh_with_skeleton(
   const Eigen::MatrixXd & V,
   const Eigen::MatrixXi & F,