rm DynamicSparseMatrix; rm commented code; fix bug in slice_into; +tests (#2016)

include/igl/cat.cpp

@@ -39,81 +39,7 @@ IGL_INLINE void igl::cat(
     return;
   }
 
-#if false
-  // This **must** be DynamicSparseMatrix, otherwise this implementation is
-  // insanely slow
-  DynamicSparseMatrix<Scalar, RowMajor> dyn_C;
-  if(dim == 1)
-  {
-    assert(A.cols() == B.cols());
-    dyn_C.resize(A.rows()+B.rows(),A.cols());
-  }else if(dim == 2)
-  {
-    assert(A.rows() == B.rows());
-    dyn_C.resize(A.rows(),A.cols()+B.cols());
-  }else
-  {
-    fprintf(stderr,"cat.h: Error: Unsupported dimension %d\n",dim);
-  }
-
-  dyn_C.reserve(A.nonZeros()+B.nonZeros());
-
-  // Iterate over outside of A
-  for(int k=0; k<A.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
-    {
-      dyn_C.coeffRef(it.row(),it.col()) += it.value();
-    }
-  }
-
-  // Iterate over outside of B
-  for(int k=0; k<B.outerSize(); ++k)
-  {
-    // Iterate over inside
-    for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
-    {
-      int r = (dim == 1 ? A.rows()+it.row() : it.row());
-      int c = (dim == 2 ? A.cols()+it.col() : it.col());
-      dyn_C.coeffRef(r,c) += it.value();
-    }
-  }
-
-  C = SparseMatrix<Scalar>(dyn_C);
-#elif false
-  std::vector<Triplet<Scalar> > CIJV;
-  CIJV.reserve(A.nonZeros() + B.nonZeros());
-  {
-    // Iterate over outside of A
-    for(int k=0; k<A.outerSize(); ++k)
-    {
-      // Iterate over inside
-      for(typename SparseMatrix<Scalar>::InnerIterator it (A,k); it; ++it)
-      {
-        CIJV.emplace_back(it.row(),it.col(),it.value());
-      }
-    }
-    // Iterate over outside of B
-    for(int k=0; k<B.outerSize(); ++k)
-    {
-      // Iterate over inside
-      for(typename SparseMatrix<Scalar>::InnerIterator it (B,k); it; ++it)
-      {
-        int r = (dim == 1 ? A.rows()+it.row() : it.row());
-        int c = (dim == 2 ? A.cols()+it.col() : it.col());
-        CIJV.emplace_back(r,c,it.value());
-      }
-    }
-
-  }
-
-  C = SparseMatrix<Scalar>(
-      dim == 1 ? A.rows()+B.rows() : A.rows(),
-      dim == 1 ? A.cols()          : A.cols()+B.cols());
-  C.reserve(A.nonZeros() + B.nonZeros());
-  C.setFromTriplets(CIJV.begin(),CIJV.end());
-#else
+  // This is faster than using DynamicSparseMatrix or setFromTriplets
   C = SparseMatrix<Scalar>(
       dim == 1 ? A.rows()+B.rows() : A.rows(),
       dim == 1 ? A.cols()          : A.cols()+B.cols());
@@ -181,9 +107,6 @@ IGL_INLINE void igl::cat(
     }
   }
   C.makeCompressed();
-
-#endif
-
 }
 
 template <typename Derived, class MatC>

include/igl/diag.cpp

@@ -20,24 +20,6 @@ IGL_INLINE void igl::diag(
 {
   assert(false && "Just call X.diagonal().sparseView() directly");
   V = X.diagonal().sparseView();
-  //// Get size of input
-  //int m = X.rows();
-  //int n = X.cols();
-  //V = Eigen::SparseVector<T>((m>n?n:m));
-  //V.reserve(V.size());
-
-  //// Iterate over outside
-  //for(int k=0; k<X.outerSize(); ++k)
-  //{
-  //  // Iterate over inside
-  //  for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
-  //  {
-  //    if(it.col() == it.row())
-  //    {
-  //      V.coeffRef(it.col()) += it.value();
-  //    }
-  //  }
-  //}
 }
 
 template <typename T,typename DerivedV>
@@ -47,23 +29,6 @@ IGL_INLINE void igl::diag(
 {
   assert(false && "Just call X.diagonal() directly");
   V = X.diagonal();
-  //// Get size of input
-  //int m = X.rows();
-  //int n = X.cols();
-  //V.derived().resize((m>n?n:m),1);
-
-  //// Iterate over outside
-  //for(int k=0; k<X.outerSize(); ++k)
-  //{
-  //  // Iterate over inside
-  //  for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
-  //  {
-  //    if(it.col() == it.row())
-  //    {
-  //      V(it.col()) = it.value();
-  //    }
-  //  }
-  //}
 }
 
 template <typename T>
@@ -72,14 +37,17 @@ IGL_INLINE void igl::diag(
   Eigen::SparseMatrix<T>& X)
 {
   // clear and resize output
-  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
-  dyn_X.reserve(V.size());
+  std::vector<Eigen::Triplet<double> > Xijv;
+  const int n = V.size();
+  const int nnz = V.nonZeros();
+  Xijv.reserve(nnz);
   // loop over non-zeros
   for(typename Eigen::SparseVector<T>::InnerIterator it(V); it; ++it)
   {
-    dyn_X.coeffRef(it.index(),it.index()) += it.value();
+    Xijv.emplace_back(it.index(),it.index(),it.value());
   }
-  X = Eigen::SparseMatrix<T>(dyn_X);
+  X.resize(n,n);
+  X.setFromTriplets(Xijv.begin(),Xijv.end());
 }
 
 template <typename T, typename DerivedV>
@@ -89,14 +57,13 @@ IGL_INLINE void igl::diag(
 {
   assert(V.rows() == 1 || V.cols() == 1);
   // clear and resize output
-  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(V.size(),V.size());
-  dyn_X.reserve(V.size());
+  std::vector<Eigen::Triplet<double> > Xijv;
+  const int n = V.size();
+  Xijv.reserve(n);
   // loop over non-zeros
-  for(int i = 0;i<V.size();i++)
-  {
-    dyn_X.coeffRef(i,i) += V[i];
-  }
-  X = Eigen::SparseMatrix<T>(dyn_X);
+  for(int i = 0;i<n;i++) { Xijv.emplace_back(i,i,V(i)); }
+  X.resize(n,n);
+  X.setFromTriplets(Xijv.begin(),Xijv.end());
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/slice_into.cpp

@@ -7,6 +7,7 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "slice_into.h"
 #include "colon.h"
+#include "sort.h"
 
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
@@ -20,31 +21,42 @@ IGL_INLINE void igl::slice_into(
   Eigen::SparseMatrix<T>& Y)
 {
 
-#ifndef NDEBUG
-  int xm = X.rows();
-  int xn = X.cols();
+  const int xm = X.rows();
+  const int xn = X.cols();
   assert(R.size() == xm);
   assert(C.size() == xn);
-  int ym = Y.rows();
-  int yn = Y.cols();
+  const int ym = Y.rows();
+  const int yn = Y.cols();
   assert(R.minCoeff() >= 0);
   assert(R.maxCoeff() < ym);
   assert(C.minCoeff() >= 0);
   assert(C.maxCoeff() < yn);
-#endif
 
-  // create temporary dynamic sparse matrix
-  Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_Y(Y);
-  // Iterate over outside
-  for(int k=0; k<X.outerSize(); ++k)
+  std::vector<bool> in_R(Y.rows());
+  for(int r = 0;r<R.size();r++) { in_R[R(r)] = true; }
+
+  // Rebuild each column in C
+  for(int c = 0;c<C.size();c++)
   {
-    // Iterate over inside
-    for(typename Eigen::SparseMatrix<T>::InnerIterator it (X,k); it; ++it)
+    int k = C(c);
+    Eigen::SparseVector<T> Yk = Y.col(k);
+    // implicit zeros
+    for(typename Eigen::SparseMatrix<T>::InnerIterator yit(Y, k);yit;++yit)
     {
-      dyn_Y.coeffRef(R(it.row()),C(it.col())) = it.value();
+      if(in_R[yit.row()]){ Yk.coeffRef(yit.row()) = 0; }
     }
+    // explicit values
+    for(typename Eigen::SparseMatrix<T>::InnerIterator xit(X, c);xit;++xit)
+    {
+      // row in X
+      int r = xit.row();
+      // row in Y
+      int s = R(r);
+      Yk.coeffRef(s) = xit.value();
+    }
+    Y.col(k) = Yk;
   }
-  Y = Eigen::SparseMatrix<T>(dyn_Y);
+
 }
 
 template <typename DerivedX, typename DerivedY, typename DerivedR, typename DerivedC>
@@ -57,16 +69,14 @@ IGL_INLINE void igl::slice_into(
 
   int xm = X.rows();
   int xn = X.cols();
-#ifndef NDEBUG
   assert(R.size() == xm);
   assert(C.size() == xn);
-  int ym = Y.rows();
-  int yn = Y.cols();
+  const int ym = Y.rows();
+  const int yn = Y.cols();
   assert(R.minCoeff() >= 0);
   assert(R.maxCoeff() < ym);
   assert(C.minCoeff() >= 0);
   assert(C.maxCoeff() < yn);
-#endif
 
   // Build reindexing maps for columns and rows, -1 means not in map
   Eigen::Matrix<typename DerivedR::Scalar,Eigen::Dynamic,1> RI;
@@ -128,6 +138,9 @@ IGL_INLINE void igl::slice_into(
 }
 
 #ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+// generated by autoexplicit.sh
+template void igl::slice_into<double, 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> > const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 template void igl::slice_into<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::slice_into<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -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::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::slice_into<double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::SparseMatrix<double, 0, int> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::SparseMatrix<double, 0, int>&);

include/igl/sort.cpp

@@ -328,6 +328,8 @@ if(!ascending)
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::DenseBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 // generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2>, Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1>, Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<unsigned int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 0, -1, -1> >&);

include/igl/sparse.cpp

@@ -49,17 +49,6 @@ IGL_INLINE void igl::sparse(
   assert(J.rows() == V.rows());
   assert(I.cols() == J.cols());
   assert(J.cols() == V.cols());
-  //// number of values
-  //int nv = V.size();
-
-  //Eigen::DynamicSparseMatrix<T, Eigen::RowMajor> dyn_X(m,n);
-  //// over estimate the number of entries
-  //dyn_X.reserve(I.size());
-  //for(int i = 0;i < nv;i++)
-  //{
-  //  dyn_X.coeffRef((int)I(i),(int)J(i)) += (T)V(i);
-  //}
-  //X = Eigen::SparseMatrix<T>(dyn_X);
   vector<Triplet<T> > IJV;
   IJV.reserve(I.size());
   for(int x = 0;x<I.size();x++)

tests/include/igl/diag.cpp

@@ -0,0 +1,22 @@
+#include <test_common.h>
+#include <igl/diag.h>
+
+TEST_CASE("diag: dense-vector-to-sparse", "[igl]")
+{
+  const Eigen::VectorXd v = (Eigen::VectorXd(3)<<1,2,3).finished();
+  Eigen::SparseMatrix<double> X;
+  igl::diag(v,X);
+  const Eigen::MatrixXd X_exact = 
+    (Eigen::MatrixXd(3,3)<<1,0,0,0,2,0,0,0,3).finished();
+  test_common::assert_eq(Eigen::MatrixXd(X),X_exact);
+}
+
+TEST_CASE("diag: sparse-vector-to-sparse", "[igl]")
+{
+  const Eigen::SparseVector<double> v = (Eigen::VectorXd(3)<<1,0,3).finished().sparseView();
+  Eigen::SparseMatrix<double> X;
+  igl::diag(v,X);
+  const Eigen::MatrixXd X_exact = 
+    (Eigen::MatrixXd(3,3)<<1,0,0,0,0,0,0,0,3).finished();
+  test_common::assert_eq(Eigen::MatrixXd(X),X_exact);
+}

tests/include/igl/slice_into.cpp

@@ -70,3 +70,36 @@ TEST_CASE("slice_into: sparse_identity", "[igl]")
     test_common::assert_eq(A,B);
   }
 }
+
+#include <iostream>
+#include <igl/matlab_format.h>
+TEST_CASE("slice_into: every-other", "[igl]")
+{
+  Eigen::MatrixXd Af(2,2);
+  Af<<
+    1,0,
+    5,6;
+  Eigen::SparseMatrix<double> As = Af.sparseView();
+  Eigen::MatrixXd Bf(4,4);
+  Bf<<
+    0,5,0,3,
+    0,6,0,4,
+    3,0,1,5,
+    4,8,0,0;
+  Eigen::SparseMatrix<double> Bs = Bf.sparseView();
+
+  Eigen::VectorXi R(2);
+  R<<1,3;
+  Eigen::VectorXi C(2);
+  C<<1,3;
+  igl::slice_into(Af,R,C,Bf);
+  igl::slice_into(As,R,C,Bs);
+  Eigen::MatrixXd res(4,4);
+  res<<
+    0,5,0,3,
+    0,1,0,0,
+    3,0,1,5,
+    4,5,0,6;
+  test_common::assert_eq(Bf,res);
+  test_common::assert_eq(Eigen::MatrixXd(Bs),res);
+}