use sparse matrix redux/for each

Former-commit-id: f0f6ffda94377887e1916b9dde19083517e6c4ee

include/igl/all.cpp

@@ -0,0 +1,25 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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 "all.h"
+
+
+template <typename AType, typename DerivedB>
+IGL_INLINE void igl::all(
+  const Eigen::SparseMatrix<AType> & A, 
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  typedef typename DerivedB::Scalar Scalar;
+  igl::redux(A,dim,[](Scalar a, Scalar b){ return a && b!=0;},B);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif
+
+

include/igl/all.h

@@ -0,0 +1,34 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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_ALL_H
+#define IGL_ALL_H
+#include "igl_inline.h"
+namespace igl
+{
+  // For Dense matrices use: A.rowwise().all() or A.colwise().all()
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   dim  dimension along which to check for all (1 or 2)
+  // Output:
+  //   B  n-long vector (if dim == 1) 
+  //   or
+  //   B  m-long vector (if dim == 2)
+  //
+  template <typename AType, typename DerivedB>
+  IGL_INLINE void all(
+    const Eigen::SparseMatrix<AType> & A, 
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "all.cpp"
+#endif
+#endif
+
+

include/igl/any.cpp

@@ -0,0 +1,24 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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 "any.h"
+
+
+template <typename AType, typename DerivedB>
+IGL_INLINE void igl::any(
+  const Eigen::SparseMatrix<AType> & A, 
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  typedef typename DerivedB::Scalar Scalar;
+  igl::redux(A,dim,[](Scalar a, Scalar b){ return a || b!=0;},B);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif
+

include/igl/any.h

@@ -0,0 +1,33 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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_ANY_H
+#define IGL_ANY_H
+#include "igl_inline.h"
+namespace igl
+{
+  // For Dense matrices use: A.rowwise().any() or A.colwise().any()
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   dim  dimension along which to check for any (1 or 2)
+  // Output:
+  //   B  n-long vector (if dim == 1) 
+  //   or
+  //   B  m-long vector (if dim == 2)
+  //
+  template <typename AType, typename DerivedB>
+  IGL_INLINE void any(
+    const Eigen::SparseMatrix<AType> & A, 
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB>& B);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "any.cpp"
+#endif
+#endif
+

include/igl/count.cpp

@@ -0,0 +1,61 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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 "count.h"
+#include "redux.h"
+
+template <typename XType, typename SType>
+IGL_INLINE void igl::count(
+  const Eigen::SparseMatrix<XType>& X, 
+  const int dim,
+  Eigen::SparseVector<SType>& S)
+{
+  // dim must be 2 or 1
+  assert(dim == 1 || dim == 2);
+  // Get size of input
+  int m = X.rows();
+  int n = X.cols();
+  // resize output
+  if(dim==1)
+  {
+    S = Eigen::SparseVector<SType>(n);
+  }else
+  {
+    S = Eigen::SparseVector<SType>(m);
+  }
+
+  // Iterate over outside
+  for(int k=0; k<X.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<XType>::InnerIterator it (X,k); it; ++it)
+    {
+      if(dim == 1)
+      {
+        S.coeffRef(it.col()) += (it.value() == 0? 0: 1);
+      }else
+      {
+        S.coeffRef(it.row()) += (it.value() == 0? 0: 1);
+      }
+    }
+  }
+
+}
+
+template <typename AType, typename DerivedB>
+IGL_INLINE void igl::count(
+  const Eigen::SparseMatrix<AType>& A, 
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  typedef typename DerivedB::Scalar Scalar;
+  igl::redux(A,dim,[](Scalar a, Scalar b){ return a+(b==0?0:1);},B);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/count.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@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_COUNT_H
+#define IGL_COUNT_H
+#include "igl_inline.h"
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // Note: If your looking for dense matrix matlab like sum for eigen matrics
+  // just use:
+  //   M.colwise().count() or M.rowwise().count()
+  // 
+
+  // Count the number of non-zeros in the columns or rows of a sparse matrix
+  //
+  // Inputs:
+  //   X  m by n sparse matrix
+  //   dim  dimension along which to sum (1 or 2)
+  // Output:
+  //   S  n-long sparse vector (if dim == 1) 
+  //   or
+  //   S  m-long sparse vector (if dim == 2)
+  template <typename XType, typename SType>
+  IGL_INLINE void count(
+    const Eigen::SparseMatrix<XType>& X, 
+    const int dim,
+    Eigen::SparseVector<SType>& S);
+  template <typename XType, typename DerivedS>
+  IGL_INLINE void count(
+    const Eigen::SparseMatrix<XType>& X, 
+    const int dim,
+    Eigen::PlainObjectBase<DerivedS>& S);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "count.cpp"
+#endif
+
+#endif
+

include/igl/for_each.cpp

@@ -0,0 +1,48 @@
+#include "for_each.h"
+
+template <typename AType, typename Func>
+IGL_INLINE void igl::for_each(
+  const Eigen::SparseMatrix<AType> & A,
+  const Func & func)
+{
+  // Can **not** use parallel for because this must be _in order_
+  // Iterate over outside
+  for(int k=0; k<A.outerSize(); ++k)
+  {
+    // Iterate over inside
+    for(typename Eigen::SparseMatrix<AType>::InnerIterator it (A,k); it; ++it)
+    {
+      func(it.row(),it.col(),it.value());
+    }
+  }
+}
+
+template <typename DerivedA, typename Func>
+IGL_INLINE void igl::for_each(
+  const Eigen::DenseBase<DerivedA> & A,
+  const Func & func)
+{
+  // Can **not** use parallel for because this must be _in order_
+  if(A.IsRowMajor)
+  {
+    for(typename DerivedA::Index i = 0;i<A.rows();i++)
+    {
+      for(typename DerivedA::Index j = 0;j<A.cols();j++)
+      {
+        func(i,j,A(i,j));
+      }
+    }
+  }else
+  {
+    for(typename DerivedA::Index j = 0;j<A.cols();j++)
+    {
+      for(typename DerivedA::Index i = 0;i<A.rows();i++)
+      {
+        func(i,j,A(i,j));
+      }
+    }
+  }
+  
+}
+
+#include <functional>

include/igl/for_each.h

@@ -0,0 +1,30 @@
+#ifndef IGL_FOR_EACH_H
+#define IGL_FOR_EACH_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // FOR_EACH  Call a given function for each non-zero (i.e., explicit value
+  // might actually be ==0) in a Sparse Matrix A _in order (of storage)_. This is
+  // useless unless func has _side-effects_.
+  //
+  // Inputs:
+  //   A  m by n SparseMatrix
+  //   func  function handle with prototype "compatible with" `void (Index i,
+  //     Index j, Scalar & v)`. Return values will be ignored.
+  //
+  // See also: std::for_each
+  template <typename AType, typename Func>
+  IGL_INLINE void for_each(
+    const Eigen::SparseMatrix<AType> & A,
+    const Func & func);
+  template <typename DerivedA, typename Func>
+  IGL_INLINE void for_each(
+    const Eigen::DenseBase<DerivedA> & A,
+    const Func & func);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "for_each.cpp"
+#endif
+#endif

include/igl/redux.cpp

@@ -0,0 +1,32 @@
+#include "redux.h"
+#include "for_each.h"
+
+template <typename AType, typename Func, typename DerivedB>
+IGL_INLINE void igl::redux(
+  const Eigen::SparseMatrix<AType> & A,
+  const int dim,
+  const Func & func,
+  Eigen::PlainObjectBase<DerivedB> & B)
+{
+  assert((dim == 1 || dim == 2) && "dim must be 2 or 1");
+  // Get size of input
+  int m = A.rows();
+  int n = A.cols();
+  // resize output
+  B = DerivedB::Zero(dim==1?n:m);
+  const auto func_wrap = [&func,&B,&dim](const int i, const int j, const int v)
+  {
+    if(dim == 1)
+    {
+      B(j) = i == 0? v : func(B(j),v);
+    }else
+    {
+      B(i) = j == 0? v : func(B(i),v);
+    }
+  };
+  for_each(A,func_wrap);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/redux.h

@@ -0,0 +1,38 @@
+#ifndef IGL_REDUX_H
+#define IGL_REDUX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // REDUX Perform reductions on the rows or columns of a SparseMatrix. This is
+  // _similar_ to DenseBase::redux, but different in two important ways:
+  //  1. (unstored) Zeros are **not** "visited", however if the first element
+  //     in the column/row  does not appear in the first row/column then the
+  //     reduction is assumed to start with zero. In this way, "any", "all",
+  //     "count"(non-zeros) work as expected. This means it is **not** possible
+  //     to use this to count (implicit) zeros.
+  //  2. This redux is more powerful in the sense that A and B may have
+  //     different types. This makes it possible to count the number of
+  //     non-zeros in a SparseMatrix<bool> A into a VectorXi B.
+  //
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   dim  dimension along which to sum (1 or 2)
+  //   func  function handle with the prototype `X(Y a, I i, J j, Z b)` where a
+  //     is the running value, b is A(i,j)
+  // Output:
+  //   S  n-long sparse vector (if dim == 1) 
+  //   or
+  //   S  m-long sparse vector (if dim == 2)
+  template <typename AType, typename Func, typename DerivedB>
+  IGL_INLINE void redux(
+    const Eigen::SparseMatrix<AType> & A,
+    const int dim,
+    const Func & func,
+    Eigen::PlainObjectBase<DerivedB> & B);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "redux.cpp"
+#endif
+#endif

include/igl/sum.cpp

@@ -1,11 +1,12 @@
 // This file is part of libigl, a simple c++ geometry processing library.
 // 
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// Copyright (C) 2016 Alec Jacobson <alecjacobson@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 "sum.h"
+#include "redux.h"
 
 template <typename T>
 IGL_INLINE void igl::sum(
@@ -13,8 +14,7 @@ IGL_INLINE void igl::sum(
   const int dim,
   Eigen::SparseVector<T>& S)
 {
-  // dim must be 2 or 1
-  assert(dim == 1 || dim == 2);
+  assert((dim == 1 || dim == 2) && "dim must be 2 or 1");
   // Get size of input
   int m = X.rows();
   int n = X.cols();
@@ -42,7 +42,16 @@ IGL_INLINE void igl::sum(
       }
     }
   }
+}
 
+template <typename AType, typename DerivedB>
+IGL_INLINE void igl::sum(
+  const Eigen::SparseMatrix<AType> & A, 
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB>& B)
+{
+  typedef typename DerivedB::Scalar Scalar;
+  igl::redux(A,dim,[](Scalar a, Scalar b){ return a+b;},B);
 }
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/sum.h

@@ -32,6 +32,12 @@ namespace igl
     const Eigen::SparseMatrix<T>& X, 
     const int dim,
     Eigen::SparseVector<T>& S);
+  // Sum is "conducted" in the type of DerivedB::Scalar 
+  template <typename AType, typename DerivedB>
+  IGL_INLINE void sum(
+    const Eigen::SparseMatrix<AType> & A, 
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB>& B);
 }
 
 #ifndef IGL_STATIC_LIBRARY