Merge commit 'f239034d5eb0790ee7639b2523367a71e298fd21'

include/igl/AtA_cached.cpp

@@ -0,0 +1,130 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "AtA_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <utility>
+
+template <typename Scalar>
+IGL_INLINE void igl::AtA_cached_precompute(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    igl::AtA_cached_data& data)
+{
+  // 1 Compute At (this could be avoided, but performance-wise it will not make a difference)
+  std::vector<std::vector<int> > Col_RowPtr;
+  std::vector<std::vector<int> > Col_IndexPtr;
+
+  Col_RowPtr.resize(A.cols());
+  Col_IndexPtr.resize(A.cols());
+
+  for (unsigned k=0; k<A.outerSize(); ++k)
+  {
+    unsigned outer_index = *(A.outerIndexPtr()+k);
+    unsigned next_outer_index = (k+1 == A.outerSize()) ? A.nonZeros() : *(A.outerIndexPtr()+k+1); 
+
+    for (unsigned l=outer_index; l<next_outer_index; ++l)
+    {
+      int col = k;
+      int row = *(A.innerIndexPtr()+l);
+      int value_index = l;
+      assert(col < A.cols());
+      assert(col >= 0);
+      assert(row < A.rows());
+      assert(row >= 0);
+      assert(value_index >= 0);
+      assert(value_index < A.nonZeros());
+
+      Col_RowPtr[col].push_back(row);
+      Col_IndexPtr[col].push_back(value_index);
+    }
+  }
+
+  Eigen::SparseMatrix<Scalar> At = A.transpose();
+  At.makeCompressed();
+  AtA = At * A;
+  AtA.makeCompressed();
+
+  assert(AtA.isCompressed());
+
+  // If weights are not provided, use 1
+  if (data.W.size() == 0)
+    data.W = Eigen::VectorXd::Ones(A.rows());
+  assert(data.W.size() == A.rows());
+
+  data.I_outer.reserve(AtA.outerSize());
+  data.I_row.reserve(2*AtA.nonZeros());
+  data.I_col.reserve(2*AtA.nonZeros());
+  data.I_w.reserve(2*AtA.nonZeros());
+
+  // 2 Construct the rules
+  for (unsigned k=0; k<AtA.outerSize(); ++k)
+  {
+    unsigned outer_index = *(AtA.outerIndexPtr()+k);
+    unsigned next_outer_index = (k+1 == AtA.outerSize()) ? AtA.nonZeros() : *(AtA.outerIndexPtr()+k+1); 
+
+    for (unsigned l=outer_index; l<next_outer_index; ++l)
+    {
+      int col = k;
+      int row = *(AtA.innerIndexPtr()+l);
+      int value_index = l;
+      assert(col < AtA.cols());
+      assert(col >= 0);
+      assert(row < AtA.rows());
+      assert(row >= 0);
+      assert(value_index >= 0);
+      assert(value_index < AtA.nonZeros());
+
+      data.I_outer.push_back(data.I_row.size());
+
+      // Find correspondences
+      unsigned i=0;
+      unsigned j=0;
+      while (i<Col_RowPtr[row].size() && j<Col_RowPtr[col].size())
+      {
+          if (Col_RowPtr[row][i] == Col_RowPtr[col][j])
+          {
+            data.I_row.push_back(Col_IndexPtr[row][i]);
+            data.I_col.push_back(Col_IndexPtr[col][j]);
+            data.I_w.push_back(Col_RowPtr[col][j]);
+            ++i;
+            ++j;
+          } else 
+          if (Col_RowPtr[row][i] > Col_RowPtr[col][j])
+            ++j;
+          else
+            ++i;
+
+      }
+    }
+  }
+  data.I_outer.push_back(data.I_row.size()); // makes it more efficient to iterate later on
+
+  igl::AtA_cached(A,AtA,data);
+}
+
+template <typename Scalar>
+IGL_INLINE void igl::AtA_cached(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    const igl::AtA_cached_data& data)
+{
+  for (unsigned i=0; i<data.I_outer.size()-1; ++i)
+  {
+    *(AtA.valuePtr() + i) = 0;
+    for (unsigned j=data.I_outer[i]; j<data.I_outer[i+1]; ++j)
+      *(AtA.valuePtr() + i) += *(A.valuePtr() + data.I_row[j]) * data.W[data.I_w[j]] * *(A.valuePtr() + data.I_col[j]);
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::AtA_cached<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, igl::AtA_cached_data const&);
+template void igl::AtA_cached_precompute<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, igl::AtA_cached_data&);
+#endif

include/igl/AtA_cached.h

@@ -0,0 +1,68 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_ATA_CACHED_H
+#define IGL_ATA_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{  
+  struct AtA_cached_data
+  {
+    // Weights
+    Eigen::VectorXd W;
+
+    // Flatten composition rules
+    std::vector<int> I_row;
+    std::vector<int> I_col;
+    std::vector<int> I_w;
+
+    // For each entry of AtA, points to the beginning
+    // of the composition rules
+    std::vector<int> I_outer;
+  };
+
+  // Computes At * W * A, where A is sparse and W is diagonal. Divides the 
+  // construction in two phases, one
+  // for fixing the sparsity pattern, and one to populate it with values. Compared to
+  // evaluating it directly, this version is slower for the first time (since it requires a
+  // precomputation), but faster to the subsequent evaluations.
+  //
+  // Input:
+  //   A m x n sparse matrix
+  //   data stores the precomputed sparsity pattern, data.W contains the optional diagonal weights (stored as a dense vector). If W is not provided, it is replaced by the identity.
+  // Outputs:
+  //   AtA  m by m matrix computed as AtA * W * A
+  //
+  // Example:
+  // AtA_data = igl::AtA_cached_data();
+  // AtA_data.W = W;
+  // if (s.AtA.rows() == 0)
+  //   igl::AtA_cached_precompute(s.A,s.AtA,s.AtA_data);
+  // else
+  //   igl::AtA_cached(s.A,s.AtA,s.AtA_data);
+  template <typename Scalar>
+  IGL_INLINE void AtA_cached_precompute(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    AtA_cached_data& data);
+
+  template <typename Scalar>
+  IGL_INLINE void AtA_cached(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    const AtA_cached_data& data);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "AtA_cached.cpp"
+#endif
+
+#endif

include/igl/slice.cpp

@@ -350,6 +350,7 @@ template void igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<in
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 3, 0, -1, 3>&);
 template void igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::slice<unsigned int, unsigned int>(Eigen::SparseMatrix<unsigned int, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<unsigned int, 0, int>&);
 #ifdef WIN32
 template void igl::slice<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>>>(class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> &);
 template void igl::slice<class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>>(class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> &);

include/igl/slice_cached.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "slice_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <igl/slice.h>
+
+template <typename TX, typename TY>
+IGL_INLINE void igl::slice_cached_precompute(
+  const Eigen::SparseMatrix<TX>& X,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+  Eigen::SparseMatrix<TY>& Y,
+  Eigen::VectorXi& data)
+{
+  // Create a sparse matrix whose entries are the ids
+  Eigen::SparseMatrix<unsigned> TS = X.template cast<unsigned>();
+
+  TS.makeCompressed();
+  for (unsigned i=0;i<TS.nonZeros();++i)
+    *(TS.valuePtr() + i) = i;
+
+  Eigen::SparseMatrix<unsigned> TS_sliced;
+  igl::slice(TS,R,C,TS_sliced);
+  Y = TS_sliced.cast<TY>();
+
+  data.resize(TS_sliced.nonZeros());
+  for (unsigned i=0;i<data.size();++i)
+  {
+    data[i] = *(TS_sliced.valuePtr() + i);
+    *(Y.valuePtr() + i) = *(X.valuePtr() + data[i]);
+  }
+}
+
+template <typename TX, typename TY>
+IGL_INLINE void igl::slice_cached(
+  const Eigen::SparseMatrix<TX>& X,
+  Eigen::SparseMatrix<TY>& Y,
+  const Eigen::VectorXi& data)
+{
+  for (unsigned i=0; i<data.size(); ++i)
+    *(Y.valuePtr() + i) = *(X.valuePtr() + data[i]);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::slice_cached_precompute<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::slice_cached<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
+#endif

include/igl/slice_cached.h

@@ -0,0 +1,66 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SLICE_CACHED_H
+#define IGL_SLICE_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{  
+
+  // Act like the matlab X(row_indices,col_indices) operator, where
+  // row_indices, col_indices are non-negative integer indices. This is a fast version
+  // of igl::slice that can analyze and store the sparsity structure. It is slower at the // first evaluation (slice_cached_precompute), but faster on the subsequent ones.
+  // 
+  // Inputs:
+  //   X  m by n matrix
+  //   R  list of row indices
+  //   C  list of column indices
+  //   
+  // Output:
+  //   Y  #R by #C matrix
+  //   data Temporary data used by slice_cached to repeat this operation
+  //
+  // Usage:
+  //
+  // // Construct and slice up Laplacian
+  // SparseMatrix<double> L,L_sliced;
+  // igl::cotmatrix(V,F,L);
+
+  // // Normal igl::slice call
+  // igl::slice(L,in,in,L_in_in);
+
+  // // Fast version
+  // static VectorXi data; // static or saved in a global state
+  // if (data.size() == 0)
+  //   igl::slice_cached_precompute(L,in,in,L_sliced,data);
+  // else
+  //   igl::slice_cached(L,L_sliced,temp);
+
+  template <typename TX, typename TY>
+  IGL_INLINE void slice_cached_precompute(
+    const Eigen::SparseMatrix<TX>& X,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+    Eigen::SparseMatrix<TY>& Y,
+    Eigen::VectorXi& data);
+
+  template <typename TX, typename TY>
+  IGL_INLINE void slice_cached(
+    const Eigen::SparseMatrix<TX>& X,
+    Eigen::SparseMatrix<TY>& Y,
+    const Eigen::VectorXi& data);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slice_cached.cpp"
+#endif
+
+#endif