diff --git a/cmake/LibiglDownloadExternal.cmake b/cmake/LibiglDownloadExternal.cmake index 600c35d8c..7e37dcc8a 100644 --- a/cmake/LibiglDownloadExternal.cmake +++ b/cmake/LibiglDownloadExternal.cmake @@ -156,7 +156,7 @@ endfunction() function(igl_download_catch2) igl_download_project(catch2 GIT_REPOSITORY https://github.com/catchorg/Catch2.git - GIT_TAG 03d122a35c3f5c398c43095a87bc82ed44642516 + GIT_TAG v2.11.0 ) endfunction() diff --git a/include/igl/slice.cpp b/include/igl/slice.cpp index ced1965be..8b0fb3843 100644 --- a/include/igl/slice.cpp +++ b/include/igl/slice.cpp @@ -9,7 +9,6 @@ #include "colon.h" #include -#include template < typename TX, @@ -18,136 +17,8 @@ template < typename DerivedC> IGL_INLINE void igl::slice( const Eigen::SparseMatrix &X, - const Eigen::MatrixBase &R, - const Eigen::MatrixBase &C, - Eigen::SparseMatrix &Y) -{ -#if 1 - int xm = X.rows(); - int xn = X.cols(); - int ym = R.size(); - int yn = C.size(); - - // special case when R or C is empty - if (ym == 0 || yn == 0) - { - Y.resize(ym, yn); - return; - } - - assert(R.minCoeff() >= 0); - assert(R.maxCoeff() < xm); - assert(C.minCoeff() >= 0); - assert(C.maxCoeff() < xn); - - // Build reindexing maps for columns and rows, -1 means not in map - std::vector> RI; - RI.resize(xm); - for (int i = 0; i < ym; i++) - { - RI[R(i)].push_back(i); - } - std::vector> CI; - CI.resize(xn); - // initialize to -1 - for (int i = 0; i < yn; i++) - { - CI[C(i)].push_back(i); - } - // Resize output - Eigen::DynamicSparseMatrix dyn_Y(ym, yn); - // Take a guess at the number of nonzeros (this assumes uniform distribution - // not banded or heavily diagonal) - dyn_Y.reserve((X.nonZeros() / (X.rows() * X.cols())) * (ym * yn)); - // Iterate over outside - for (int k = 0; k < X.outerSize(); ++k) - { - // Iterate over inside - for (typename Eigen::SparseMatrix::InnerIterator it(X, k); it; ++it) - { - typename std::vector::iterator rit; - typename std::vector::iterator cit; - for (rit = RI[it.row()].begin(); rit != RI[it.row()].end(); rit++) - { - for (cit = CI[it.col()].begin(); cit != CI[it.col()].end(); cit++) - { - dyn_Y.coeffRef(*rit, *cit) = it.value(); - } - } - } - } - Y = Eigen::SparseMatrix(dyn_Y); -#else - - // Alec: This is _not_ valid for arbitrary R,C since they don't necessary - // representation a strict permutation of the rows and columns: rows or - // columns could be removed or replicated. The removal of rows seems to be - // handled here (although it's not clear if there is a performance gain when - // the #removals >> #remains). If this is sufficiently faster than the - // correct code above, one could test whether all entries in R and C are - // unique and apply the permutation version if appropriate. - // - - int xm = X.rows(); - int xn = X.cols(); - int ym = R.size(); - int yn = C.size(); - - // special case when R or C is empty - if (ym == 0 || yn == 0) - { - Y.resize(ym, yn); - return; - } - - assert(R.minCoeff() >= 0); - assert(R.maxCoeff() < xm); - assert(C.minCoeff() >= 0); - assert(C.maxCoeff() < xn); - - // initialize row and col permutation vectors - Eigen::VectorXi rowIndexVec = igl::LinSpaced(xm, 0, xm - 1); - Eigen::VectorXi rowPermVec = igl::LinSpaced(xm, 0, xm - 1); - for (int i = 0; i < ym; i++) - { - int pos = rowIndexVec.coeffRef(R(i)); - if (pos != i) - { - int &val = rowPermVec.coeffRef(i); - std::swap(rowIndexVec.coeffRef(val), rowIndexVec.coeffRef(R(i))); - std::swap(rowPermVec.coeffRef(i), rowPermVec.coeffRef(pos)); - } - } - Eigen::PermutationMatrix rowPerm(rowIndexVec); - - Eigen::VectorXi colIndexVec = igl::LinSpaced(xn, 0, xn - 1); - Eigen::VectorXi colPermVec = igl::LinSpaced(xn, 0, xn - 1); - for (int i = 0; i < yn; i++) - { - int pos = colIndexVec.coeffRef(C(i)); - if (pos != i) - { - int &val = colPermVec.coeffRef(i); - std::swap(colIndexVec.coeffRef(val), colIndexVec.coeffRef(C(i))); - std::swap(colPermVec.coeffRef(i), colPermVec.coeffRef(pos)); - } - } - Eigen::PermutationMatrix colPerm(colPermVec); - - Eigen::SparseMatrix M = (rowPerm * X); - Y = (M * colPerm).block(0, 0, ym, yn); -#endif -} - -template < - typename TX, - typename TY, - typename DerivedR, - typename DerivedC> -IGL_INLINE void igl::slice( - const Eigen::SparseMatrix &X, - const Eigen::ArrayBase &R, - const Eigen::MatrixBase &C, + const Eigen::DenseBase &R, + const Eigen::DenseBase &C, Eigen::SparseMatrix &Y) { int xm = X.rows(); @@ -167,7 +38,7 @@ IGL_INLINE void igl::slice( assert(C.minCoeff() >= 0); assert(C.maxCoeff() < xn); - // Build reindexing maps for columns and rows, -1 means not in map + // Build reindexing maps for columns and rows std::vector> RI; RI.resize(xm); for (int i = 0; i < ym; i++) @@ -176,141 +47,39 @@ IGL_INLINE void igl::slice( } std::vector> CI; CI.resize(xn); - // initialize to -1 for (int i = 0; i < yn; i++) { CI[C(i)].push_back(i); } - // Resize output - Eigen::DynamicSparseMatrix dyn_Y(ym, yn); + // Take a guess at the number of nonzeros (this assumes uniform distribution // not banded or heavily diagonal) - dyn_Y.reserve((X.nonZeros() / (X.rows() * X.cols())) * (ym * yn)); + std::vector> entries; + entries.reserve((X.nonZeros()/(X.rows()*X.cols())) * (ym*yn)); + // Iterate over outside for (int k = 0; k < X.outerSize(); ++k) { // Iterate over inside for (typename Eigen::SparseMatrix::InnerIterator it(X, k); it; ++it) { - typename std::vector::iterator rit; - typename std::vector::iterator cit; - for (rit = RI[it.row()].begin(); rit != RI[it.row()].end(); rit++) + for (auto rit = RI[it.row()].begin(); rit != RI[it.row()].end(); rit++) { - for (cit = CI[it.col()].begin(); cit != CI[it.col()].end(); cit++) + for (auto cit = CI[it.col()].begin(); cit != CI[it.col()].end(); cit++) { - dyn_Y.coeffRef(*rit, *cit) = it.value(); + entries.emplace_back(*rit, *cit, it.value()); } } } } - Y = Eigen::SparseMatrix(dyn_Y); -} - -template < - typename TX, - typename TY, - typename DerivedR, - typename DerivedC> -IGL_INLINE void igl::slice( - const Eigen::SparseMatrix &X, - const Eigen::MatrixBase &R, - const Eigen::ArrayBase &C, - Eigen::SparseMatrix &Y) -{ - int xm = X.rows(); - int xn = X.cols(); - int ym = R.size(); - int yn = C.size(); - - // special case when R or C is empty - if (ym == 0 || yn == 0) - { - Y.resize(ym, yn); - return; - } - - assert(R.minCoeff() >= 0); - assert(R.maxCoeff() < xm); - assert(C.minCoeff() >= 0); - assert(C.maxCoeff() < xn); - - // Build reindexing maps for columns and rows, -1 means not in map - std::vector> RI; - RI.resize(xm); - for (int i = 0; i < ym; i++) - { - RI[R(i)].push_back(i); - } - std::vector> CI; - CI.resize(xn); - // initialize to -1 - for (int i = 0; i < yn; i++) - { - CI[C(i)].push_back(i); - } - // Resize output - Eigen::DynamicSparseMatrix dyn_Y(ym, yn); - // Take a guess at the number of nonzeros (this assumes uniform distribution - // not banded or heavily diagonal) - dyn_Y.reserve((X.nonZeros() / (X.rows() * X.cols())) * (ym * yn)); - // Iterate over outside - for (int k = 0; k < X.outerSize(); ++k) - { - // Iterate over inside - for (typename Eigen::SparseMatrix::InnerIterator it(X, k); it; ++it) - { - typename std::vector::iterator rit; - typename std::vector::iterator cit; - for (rit = RI[it.row()].begin(); rit != RI[it.row()].end(); rit++) - { - for (cit = CI[it.col()].begin(); cit != CI[it.col()].end(); cit++) - { - dyn_Y.coeffRef(*rit, *cit) = it.value(); - } - } - } - } - Y = Eigen::SparseMatrix(dyn_Y); + Y.resize(ym, yn); + Y.setFromTriplets(entries.begin(), entries.end()); } template IGL_INLINE void igl::slice( const MatX &X, - const Eigen::MatrixBase &R, - const int dim, - MatY &Y) -{ - Eigen::Matrix C; - switch (dim) - { - case 1: - // boring base case - if (X.cols() == 0) - { - Y.resize(R.size(), 0); - return; - } - igl::colon(0, X.cols() - 1, C); - return slice(X, R, C, Y); - case 2: - // boring base case - if (X.rows() == 0) - { - Y.resize(0, R.size()); - return; - } - igl::colon(0, X.rows() - 1, C); - return slice(X, C, R, Y); - default: - assert(false && "Unsupported dimension"); - return; - } -} - -template -IGL_INLINE void igl::slice( - const MatX &X, - const Eigen::ArrayBase &R, + const Eigen::DenseBase &R, const int dim, MatY &Y) { @@ -347,9 +116,9 @@ template < typename DerivedC, typename DerivedY> IGL_INLINE void igl::slice( - const Eigen::ArrayBase &X, - const Eigen::MatrixBase &R, - const Eigen::MatrixBase &C, + const Eigen::DenseBase &X, + const Eigen::DenseBase &R, + const Eigen::DenseBase &C, Eigen::PlainObjectBase &Y) { #ifndef NDEBUG @@ -383,52 +152,10 @@ IGL_INLINE void igl::slice( } } -template < - typename DerivedX, - typename DerivedR, - typename DerivedC, - typename DerivedY> -IGL_INLINE void igl::slice( - const Eigen::MatrixBase &X, - const Eigen::MatrixBase &R, - const Eigen::MatrixBase &C, - Eigen::PlainObjectBase &Y) -{ -#ifndef NDEBUG - int xm = X.rows(); - int xn = X.cols(); -#endif - int ym = R.size(); - int yn = C.size(); - - // special case when R or C is empty - if (ym == 0 || yn == 0) - { - Y.resize(ym, yn); - return; - } - - assert(R.minCoeff() >= 0); - assert(R.maxCoeff() < xm); - assert(C.minCoeff() >= 0); - assert(C.maxCoeff() < xn); - - // Resize output - Y.resize(ym, yn); - // loop over output rows, then columns - for (int i = 0; i < ym; i++) - { - for (int j = 0; j < yn; j++) - { - Y(i, j) = X(R(i,0), C(j,0)); - } - } -} - template IGL_INLINE void igl::slice( - const Eigen::MatrixBase &X, - const Eigen::MatrixBase &R, + const Eigen::DenseBase &X, + const Eigen::DenseBase &R, Eigen::PlainObjectBase &Y) { // phony column indices @@ -440,8 +167,8 @@ IGL_INLINE void igl::slice( template IGL_INLINE DerivedX igl::slice( - const Eigen::MatrixBase &X, - const Eigen::MatrixBase &R) + const Eigen::DenseBase &X, + const Eigen::DenseBase &R) { DerivedX Y; igl::slice(X, R, Y); @@ -450,8 +177,8 @@ IGL_INLINE DerivedX igl::slice( template IGL_INLINE DerivedX igl::slice( - const Eigen::MatrixBase &X, - const Eigen::MatrixBase &R, + const Eigen::DenseBase &X, + const Eigen::DenseBase &R, const int dim) { DerivedX Y; @@ -460,64 +187,61 @@ IGL_INLINE DerivedX igl::slice( } #ifdef IGL_STATIC_LIBRARY -// Explicit template instantiation -// generated by autoexplicit.sh -template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, int, Eigen::PlainObjectBase >&); -// generated by autoexplicit.sh -template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, int, Eigen::PlainObjectBase >&); -template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::Matrix &); -template void igl::slice>, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::Matrix &); -template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::MatrixBase> const &, int, Eigen::SparseMatrix &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::Matrix const &, Eigen::MatrixBase> const &, int, Eigen::Matrix &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::MatrixBase> const &, int, Eigen::SparseMatrix &); -template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::MatrixBase> const &, int, Eigen::SparseMatrix &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Block const, -1, 1, true>>(Eigen::MatrixBase> const &, Eigen::MatrixBase const, -1, 1, true>> const &, Eigen::PlainObjectBase> &); -template Eigen::Matrix igl::slice, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::MatrixBase> const &, int, Eigen::SparseMatrix &); -template Eigen::Matrix igl::slice, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::Matrix const &, Eigen::MatrixBase> const &, int, Eigen::Matrix &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); + +template Eigen::Matrix igl::slice, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, int); +template Eigen::Matrix igl::slice, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::Matrix const &, Eigen::DenseBase> const &, int, Eigen::Matrix &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Block const, -1, 1, true>>(Eigen::DenseBase> const &, Eigen::DenseBase const, -1, 1, true>> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::DenseBase> const &, Eigen::DenseBase> const &, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::Matrix>(Eigen::Matrix const &, Eigen::DenseBase> const &, int, Eigen::Matrix &); +template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::Matrix const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice >, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::Matrix&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::Matrix&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::Matrix&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::DenseBase > const&, int, Eigen::PlainObjectBase >&); +template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::DenseBase> const &, int, Eigen::SparseMatrix &); +template void igl::slice, Eigen::Array, Eigen::SparseMatrix >(Eigen::SparseMatrix const&, Eigen::DenseBase > const&, int, Eigen::SparseMatrix&); +template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::DenseBase> const &, int, Eigen::SparseMatrix &); +template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::DenseBase> const &, int, Eigen::SparseMatrix &); +template void igl::slice, Eigen::Matrix, Eigen::SparseMatrix>(Eigen::SparseMatrix const &, Eigen::DenseBase> const &, int, Eigen::SparseMatrix &); + #ifdef WIN32 -template void igl::slice, Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase>>(Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice >,class Eigen::Matrix<__int64,-1,1,0,-1,1>,class Eigen::PlainObjectBase > >(class Eigen::MatrixBase > const &,class Eigen::MatrixBase > const &,int,class Eigen::PlainObjectBase > &); +template void igl::slice >,class Eigen::Matrix<__int64,-1,1,0,-1,1>,class Eigen::PlainObjectBase > >(class Eigen::DenseBase > const &,class Eigen::DenseBase > const &,int,class Eigen::PlainObjectBase > &); +template void igl::slice >,class Eigen::Matrix<__int64,-1,1,0,-1,1>,class Eigen::PlainObjectBase > >(class Eigen::MatrixBase > const &,class Eigen::DenseBase > const &,int,class Eigen::PlainObjectBase > &); +template void igl::slice, Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase>>(Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); +template void igl::slice>, Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::DenseBase> const &, int, Eigen::PlainObjectBase> &); #endif -#if EIGEN_VERSION_AT_LEAST(3, 3, 7) -#else -template void igl::slice, Eigen::Matrix const>>, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase, Eigen::Matrix const>> const &, Eigen::MatrixBase> const &, int, Eigen::Matrix &); -#endif -template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, int, Eigen::PlainObjectBase >&); -template void igl::slice, Eigen::Matrix, Eigen::Matrix, Eigen::Matrix>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::PlainObjectBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, int, Eigen::PlainObjectBase >&); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice >, Eigen::Matrix, Eigen::PlainObjectBase > >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, int, Eigen::PlainObjectBase >&); -template void igl::slice>, Eigen::Matrix, Eigen::PlainObjectBase>>(Eigen::MatrixBase> const &, Eigen::MatrixBase> const &, int, Eigen::PlainObjectBase> &); -template void igl::slice, Eigen::Array, Eigen::SparseMatrix >(Eigen::SparseMatrix const&, Eigen::ArrayBase > const&, int, Eigen::SparseMatrix&); -template void igl::slice >, Eigen::Matrix, Eigen::Matrix >(Eigen::MatrixBase > const&, Eigen::MatrixBase > const&, int, Eigen::Matrix&); + #endif diff --git a/include/igl/slice.h b/include/igl/slice.h index d3cfd9573..fcb6bbbeb 100644 --- a/include/igl/slice.h +++ b/include/igl/slice.h @@ -30,29 +30,10 @@ namespace igl typename DerivedC> IGL_INLINE void slice( const Eigen::SparseMatrix& X, - const Eigen::MatrixBase & R, - const Eigen::MatrixBase & C, - Eigen::SparseMatrix& Y); - template < - typename TX, - typename TY, - typename DerivedR, - typename DerivedC> - IGL_INLINE void slice( - const Eigen::SparseMatrix& X, - const Eigen::ArrayBase & R, - const Eigen::MatrixBase & C, - Eigen::SparseMatrix& Y); - template < - typename TX, - typename TY, - typename DerivedR, - typename DerivedC> - IGL_INLINE void slice( - const Eigen::SparseMatrix& X, - const Eigen::MatrixBase & R, - const Eigen::ArrayBase & C, + const Eigen::DenseBase & R, + const Eigen::DenseBase & C, Eigen::SparseMatrix& Y); + // Wrapper to only slice in one direction // // Inputs: @@ -65,45 +46,27 @@ namespace igl typename MatY> IGL_INLINE void slice( const MatX& X, - const Eigen::MatrixBase & R, + const Eigen::DenseBase & R, const int dim, MatY& Y); - template < - typename MatX, - typename DerivedR, - typename MatY> - IGL_INLINE void slice( - const MatX& X, - const Eigen::ArrayBase & R, - const int dim, - MatY& Y); - template < - typename DerivedX, - typename DerivedR, - typename DerivedC, - typename DerivedY> - IGL_INLINE void slice( - const Eigen::MatrixBase & X, - const Eigen::MatrixBase & R, - const Eigen::MatrixBase & C, - Eigen::PlainObjectBase & Y); template < typename DerivedX, typename DerivedR, typename DerivedC, typename DerivedY> -IGL_INLINE void slice( - const Eigen::ArrayBase &X, - const Eigen::MatrixBase &R, - const Eigen::MatrixBase &C, - Eigen::PlainObjectBase &Y); + IGL_INLINE void slice( + const Eigen::DenseBase & X, + const Eigen::DenseBase & R, + const Eigen::DenseBase & C, + Eigen::PlainObjectBase & Y); template IGL_INLINE void slice( - const Eigen::MatrixBase & X, - const Eigen::MatrixBase & R, + const Eigen::DenseBase & X, + const Eigen::DenseBase & R, Eigen::PlainObjectBase & Y); + // VectorXi Y = slice(X,R); // // This templating is bad because the return type might not have the same @@ -112,12 +75,12 @@ IGL_INLINE void slice( // the number of rows in `DerivedX`. template IGL_INLINE DerivedX slice( - const Eigen::MatrixBase & X, - const Eigen::MatrixBase & R); + const Eigen::DenseBase & X, + const Eigen::DenseBase & R); template IGL_INLINE DerivedX slice( - const Eigen::MatrixBase& X, - const Eigen::MatrixBase & R, + const Eigen::DenseBase& X, + const Eigen::DenseBase & R, const int dim); } diff --git a/include/igl/slice_mask.cpp b/include/igl/slice_mask.cpp index 669f4e231..9fd5b45db 100644 --- a/include/igl/slice_mask.cpp +++ b/include/igl/slice_mask.cpp @@ -1,12 +1,13 @@ // This file is part of libigl, a simple c++ geometry processing library. -// +// // Copyright (C) 2015 Alec Jacobson -// -// 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 +// +// 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_mask.h" #include "slice.h" +#include "slice_sorted.h" #include "find.h" #include @@ -143,7 +144,7 @@ IGL_INLINE void igl::slice_mask( find(R,Ri); Eigen::VectorXi Ci; find(C,Ci); - return slice(X,Ri,Ci,Y); + return slice_sorted(X,Ri,Ci,Y); } #ifdef IGL_STATIC_LIBRARY diff --git a/include/igl/slice_sorted.cpp b/include/igl/slice_sorted.cpp new file mode 100644 index 000000000..000a20aa9 --- /dev/null +++ b/include/igl/slice_sorted.cpp @@ -0,0 +1,88 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2013 Alec Jacobson +// +// 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_sorted.h" + +#include + +// TODO: Write a version that works for row-major sparse matrices as well. +template +IGL_INLINE void igl::slice_sorted(const Eigen::SparseMatrix &X, + const Eigen::DenseBase &R, + const Eigen::DenseBase &C, + Eigen::SparseMatrix &Y) +{ + int xm = X.rows(); + int xn = X.cols(); + int ym = R.size(); + int yn = C.size(); + + // Special case when R or C is empty + if (ym == 0 || yn == 0) + { + Y.resize(ym, yn); + return; + } + + assert(R.minCoeff() >= 0); + assert(R.maxCoeff() < xm); + assert(C.minCoeff() >= 0); + assert(C.maxCoeff() < xn); + + // Multiplicity count for each row/col + using RowIndexType = typename DerivedR::Scalar; + using ColIndexType = typename DerivedC::Scalar; + std::vector slicedRowStart(xm); + std::vector rowRepeat(xm, 0); + for (int i = 0; i < ym; ++i) + { + if (rowRepeat[R(i)] == 0) + { + slicedRowStart[R(i)] = i; + } + rowRepeat[R(i)]++; + } + std::vector columnRepeat(xn, 0); + for (int i = 0; i < yn; i++) + { + columnRepeat[C(i)]++; + } + // Count number of nnz per outer row/col + Eigen::VectorXi nnz(yn); + for (int k = 0, c = 0; k < X.outerSize(); ++k) + { + int cnt = 0; + for (typename Eigen::SparseMatrix::InnerIterator it(X, k); it; ++it) + { + cnt += rowRepeat[it.row()]; + } + for (int i = 0; i < columnRepeat[k]; ++i, ++c) + { + nnz(c) = cnt; + } + } + Y.resize(ym, yn); + Y.reserve(nnz); + // Insert values + for (int k = 0, c = 0; k < X.outerSize(); ++k) + { + for (int i = 0; i < columnRepeat[k]; ++i, ++c) + { + for (typename Eigen::SparseMatrix::InnerIterator it(X, k); it; ++it) + { + for (int j = 0, r = slicedRowStart[it.row()]; j < rowRepeat[it.row()]; ++j, ++r) + { + Y.insert(r, c) = it.value(); + } + } + } + } +} + +#ifdef IGL_STATIC_LIBRARY +template void igl::slice_sorted, Eigen::Matrix >(Eigen::SparseMatrix const&, Eigen::DenseBase > const&, Eigen::DenseBase > const&, Eigen::SparseMatrix&); +#endif diff --git a/include/igl/slice_sorted.h b/include/igl/slice_sorted.h new file mode 100644 index 000000000..47826b927 --- /dev/null +++ b/include/igl/slice_sorted.h @@ -0,0 +1,42 @@ +// This file is part of libigl, a simple c++ geometry processing library. +// +// Copyright (C) 2019 Jérémie Dumas +// +// 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_SORTED_H +#define IGL_SLICE_SORTED_H + +#include "igl_inline.h" +#include +#include + +namespace igl +{ + // Act like the matlab X(row_indices,col_indices) operator, where row_indices, + // col_indices are non-negative integer indices. This version is about 2x faster + // than igl::slice, but it assumes that the indices to slice with are already sorted. + // + // Inputs: + // X m by n matrix + // R list of row indices + // C list of column indices + // + // Output: + // Y #R by #C matrix + // + template + IGL_INLINE void slice_sorted(const Eigen::SparseMatrix &X, + const Eigen::DenseBase &R, + const Eigen::DenseBase &C, + Eigen::SparseMatrix &Y); + +} // namespace igl + +#ifndef IGL_STATIC_LIBRARY +#include "slice_sorted.cpp" +#endif + +#endif diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 86c487ce0..273324d7a 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -27,6 +27,7 @@ add_subdirectory(${LIBIGL_EXTERNAL}/catch2 catch2) add_executable(libigl_tests main.cpp test_common.h) target_link_libraries(libigl_tests PUBLIC igl::core Catch2::Catch2) target_include_directories(libigl_tests PUBLIC ${CMAKE_CURRENT_LIST_DIR}) +target_compile_definitions(libigl_tests PUBLIC CATCH_CONFIG_ENABLE_BENCHMARKING) # Set DATA_DIR definition set(DATA_DIR "${CMAKE_CURRENT_SOURCE_DIR}/data/") diff --git a/tests/include/igl/slice_sorted.cpp b/tests/include/igl/slice_sorted.cpp new file mode 100644 index 000000000..026c070ed --- /dev/null +++ b/tests/include/igl/slice_sorted.cpp @@ -0,0 +1,79 @@ +#include +#include +#include +#include + +namespace +{ + Eigen::SparseMatrix generate_random_sparse_matrix(int rows, int cols) + { + std::mt19937 gen; + std::uniform_real_distribution dist(0.0, 1.0); + + using T = Eigen::Triplet; + std::vector tripletList; + for (int i = 0; i < rows; ++i) + { + for (int j = 0; j < cols; ++j) + { + auto v_ij = dist(gen); // generate random number + if (v_ij < 0.1) + { + tripletList.push_back(T(i, j, v_ij)); // if larger than treshold, insert it + } + } + } + Eigen::SparseMatrix mat(rows, cols); + mat.setFromTriplets(tripletList.begin(), tripletList.end()); // create the matrix + return mat; + } + +} // namespace + +TEST_CASE("slice_sorted: correctness", "[igl]") +{ + constexpr int rows = 1e3; + constexpr int cols = 1e3; + + Eigen::SparseMatrix M = generate_random_sparse_matrix(rows, cols); + + Eigen::Matrix R(rows / 2); + Eigen::Matrix C(cols / 2); + for (int i = 0; i < rows; i += 2) R[i / 2] = i; + for (int i = 0; i < cols; i += 2) C[i / 2] = i; + + SECTION("correctness") + { + // Check for correctness + Eigen::SparseMatrix A, B; + igl::slice(M, R, C, A); + igl::slice_sorted(M, R, C, B); + REQUIRE((A - B).norm() == 0); + } +} + +TEST_CASE("slice_sorted: benchmark", "[igl]" IGL_DEBUG_OFF) +{ + constexpr int rows = 1e3; + constexpr int cols = 1e3; + + Eigen::SparseMatrix M = generate_random_sparse_matrix(rows, cols); + + Eigen::Matrix R(rows / 2); + Eigen::Matrix C(cols / 2); + for (int i = 0; i < rows; i += 2) R[i / 2] = i; + for (int i = 0; i < cols; i += 2) C[i / 2] = i; + + BENCHMARK("igl::slice") { + Eigen::SparseMatrix A; + igl::slice(M, R, C, A); + return A.norm(); + }; + + BENCHMARK("igl::slice_sorted") { + Eigen::SparseMatrix A; + igl::slice_sorted(M, R, C, A); + return A.norm(); + }; +} + diff --git a/tests/test_common.h b/tests/test_common.h index e15203324..d7bd8c9ff 100644 --- a/tests/test_common.h +++ b/tests/test_common.h @@ -1,8 +1,7 @@ #pragma once - // These are not directly used but would otherwise be included in most files. -// Leaving them included here. +// Leaving them included here. #include #include @@ -17,6 +16,13 @@ #include #include +// Disable lengthy tests in debug mode +#ifdef NDEBUG +#define IGL_DEBUG_OFF "" +#else +#define IGL_DEBUG_OFF "[!hide]" +#endif + namespace test_common { template