reductions for sparse matrices

Former-commit-id: ff21ae6a99cd54d4a8461e6b60bc93e2e3e7abde

include/igl/all.cpp

@@ -6,6 +6,7 @@
 // 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"
+#include "redux.h"
 
 
 template <typename AType, typename DerivedB>

include/igl/all.h

@@ -8,6 +8,8 @@
 #ifndef IGL_ALL_H
 #define IGL_ALL_H
 #include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
 namespace igl
 {
   // For Dense matrices use: A.rowwise().all() or A.colwise().all()

include/igl/any.cpp

@@ -6,6 +6,7 @@
 // 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"
+#include "redux.h"
 
 
 template <typename AType, typename DerivedB>

include/igl/any.h

@@ -8,6 +8,8 @@
 #ifndef IGL_ANY_H
 #define IGL_ANY_H
 #include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
 namespace igl
 {
   // For Dense matrices use: A.rowwise().any() or A.colwise().any()

include/igl/find_zero.cpp

@@ -0,0 +1,42 @@
+#include "find_zero.h"
+#include "for_each.h"
+#include "any.h"
+
+template <typename AType, typename DerivedI>
+IGL_INLINE void igl::find_zero(
+  const Eigen::SparseMatrix<AType> & A,
+  const int dim,
+  Eigen::PlainObjectBase<DerivedI> & I)
+{
+  assert((dim == 1 || dim == 2) && "dim must be 2 or 1");
+  // Get size of input
+  int m = A.rows();
+  int n = A.cols();
+  // I starts by containing guess where 0 might be
+  I = DerivedI::Zero(dim==1?n:m);
+  Eigen::Array<bool,Eigen::Dynamic,1> found = 
+    Eigen::Array<bool,Eigen::Dynamic,1>::Zero(dim==1?n:m);
+  const auto func = [&I,&found,&dim](int i, int j, const int v)
+  {
+    if(dim == 2)
+    {
+      std::swap(i,j);
+    }
+    // Coded as if dim == 1, assuming swap for dim == 2
+    // Have we already found a zero?
+    if(!found(j))
+    {
+      if(I(j) != i || v == 0)
+      {
+        // either there was an implicit zero between the last element and this
+        // one, or this one is zero
+        found(j) = true;
+      }else
+      {
+        // If not found, then guess that next element will be zero
+        I(j) = I(j)+1;
+      }
+    }
+  };
+  for_each(A,func);
+}

include/igl/find_zero.h

@@ -0,0 +1,28 @@
+#ifndef IGL_FIND_ZERO_H
+#define IGL_FIND_ZERO_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Find the first zero (whether implicit or explicitly stored) in the
+  // rows/columns of a matrix.
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   dim  dimension along which to check for any (1 or 2)
+  // Output:
+  //   I  n-long vector (if dim == 1)  {m means no zeros found}
+  //   or
+  //   I  m-long vector (if dim == 2)  {n means no zeros found}
+  //
+  template <typename AType, typename DerivedI>
+  IGL_INLINE void find_zero(
+    const Eigen::SparseMatrix<AType> & A,
+    const int dim,
+    Eigen::PlainObjectBase<DerivedI> & I);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "find_zero.cpp"
+#endif
+#endif
+

include/igl/max.cpp

@@ -0,0 +1,40 @@
+#include "max.h"
+#include "for_each.h"
+#include "find_zero.h"
+
+template <typename AType, typename DerivedB, typename DerivedI>
+IGL_INLINE void igl::max(
+  const Eigen::SparseMatrix<AType> & A,
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedI> & I)
+{
+  const int n = A.cols();
+  const int m = A.rows();
+  B.resize(dim==1?n:m);
+  B.setConstant(std::numeric_limits<typename DerivedB::Scalar>::lowest());
+  I.resize(dim==1?n:m);
+  for_each(A,[&B,&I,&dim](int i, int j,const typename DerivedB::Scalar v)
+    {
+      if(dim == 2)
+      {
+        std::swap(i,j);
+      }
+      // Coded as if dim == 1, assuming swap for dim == 2
+      if(v > B(j))
+      {
+        B(j) = v;
+        I(j) = i;
+      }
+    });
+  Eigen::VectorXi Z;
+  find_zero(A,dim,Z);
+  for(int j = 0;j<I.size();j++)
+  {
+    if(Z(j) != (dim==1?m:n) && 0 > B(j))
+    {
+      B(j) = 0;
+      I(j) = Z(j);
+    }
+  }
+}

include/igl/max.h

@@ -0,0 +1,27 @@
+#ifndef IGL_MAX_H
+#define IGL_MAX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Inputs:
+  //   X  m by n matrix
+  //   dim  dimension along which to take max
+  // Outputs:
+  //   Y  n-long vector (if dim == 1) 
+  //   or
+  //   Y  m-long vector (if dim == 2)
+  //   I  vector the same size as Y containing the indices along dim of maximum
+  //     entries
+  template <typename AType, typename DerivedB, typename DerivedI>
+  IGL_INLINE void max(
+    const Eigen::SparseMatrix<AType> & A,
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedI> & I);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "max.cpp"
+#endif
+#endif

include/igl/min.cpp

@@ -0,0 +1,41 @@
+#include "min.h"
+#include "for_each.h"
+#include "find_zero.h"
+
+template <typename AType, typename DerivedB, typename DerivedI>
+IGL_INLINE void igl::min(
+  const Eigen::SparseMatrix<AType> & A,
+  const int dim,
+  Eigen::PlainObjectBase<DerivedB> & B,
+  Eigen::PlainObjectBase<DerivedI> & I)
+{
+  const int n = A.cols();
+  const int m = A.rows();
+  B.resize(dim==1?n:m);
+  B.setConstant(std::numeric_limits<typename DerivedB::Scalar>::max());
+  I.resize(dim==1?n:m);
+  for_each(A,[&B,&I,&dim](int i, int j,const typename DerivedB::Scalar v)
+    {
+      if(dim == 2)
+      {
+        std::swap(i,j);
+      }
+      // Coded as if dim == 1, assuming swap for dim == 2
+      if(v < B(j))
+      {
+        B(j) = v;
+        I(j) = i;
+      }
+    });
+  Eigen::VectorXi Z;
+  find_zero(A,dim,Z);
+  for(int j = 0;j<I.size();j++)
+  {
+    if(Z(j) != (dim==1?m:n) && 0 < B(j))
+    {
+      B(j) = 0;
+      I(j) = Z(j);
+    }
+  }
+}
+

include/igl/min.h

@@ -0,0 +1,28 @@
+#ifndef IGL_MIN_H
+#define IGL_MIN_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Inputs:
+  //   X  m by n matrix
+  //   dim  dimension along which to take min
+  // Outputs:
+  //   Y  n-long vector (if dim == 1) 
+  //   or
+  //   Y  m-long vector (if dim == 2)
+  //   I  vector the same size as Y containing the indices along dim of minimum
+  //     entries
+  template <typename AType, typename DerivedB, typename DerivedI>
+  IGL_INLINE void min(
+    const Eigen::SparseMatrix<AType> & A,
+    const int dim,
+    Eigen::PlainObjectBase<DerivedB> & B,
+    Eigen::PlainObjectBase<DerivedI> & I);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "min.cpp"
+#endif
+#endif
+