From 52cf07d2660099c5c05a5d183991e920eac63895 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Wed, 21 Jan 2009 18:46:04 +0000
Subject: [PATCH] sparse module: * add row(i), col(i) functions * add prune()
 function to remove small coefficients

---
 Eigen/src/Core/util/StaticAssert.h     |  1 +
 Eigen/src/Sparse/CompressedStorage.h   | 17 +++++++++++
 Eigen/src/Sparse/DynamicSparseMatrix.h | 10 +++++--
 Eigen/src/Sparse/SparseBlock.h         | 34 +++++++++++++++++++++
 Eigen/src/Sparse/SparseMatrix.h        | 24 ++++++++++++++-
 Eigen/src/Sparse/SparseMatrixBase.h    |  4 +++
 Eigen/src/Sparse/SparseVector.h        | 24 ++++++++-------
 test/sparse_basic.cpp                  | 41 ++++++++++++++++++++++++--
 8 files changed, 138 insertions(+), 17 deletions(-)

diff --git a/Eigen/src/Core/util/StaticAssert.h b/Eigen/src/Core/util/StaticAssert.h
index 350630c0a..60055683e 100644
--- a/Eigen/src/Core/util/StaticAssert.h
+++ b/Eigen/src/Core/util/StaticAssert.h
@@ -72,6 +72,7 @@
         INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION,
         YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY,
         THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES,
+        THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES,
         INVALID_MATRIX_TEMPLATE_PARAMETERS
       };
     };
diff --git a/Eigen/src/Sparse/CompressedStorage.h b/Eigen/src/Sparse/CompressedStorage.h
index d66cb5c94..4e8074b98 100644
--- a/Eigen/src/Sparse/CompressedStorage.h
+++ b/Eigen/src/Sparse/CompressedStorage.h
@@ -31,6 +31,7 @@
 template<typename Scalar>
 class CompressedStorage
 {
+    typedef typename NumTraits<Scalar>::Real RealScalar;
   public:
     CompressedStorage()
       : m_values(0), m_indices(0), m_size(0), m_allocatedSize(0)
@@ -183,6 +184,22 @@ class CompressedStorage
       }
       return m_values[id];
     }
+    
+    void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
+    {
+      int k = 0;
+      int n = size();
+      for (int i=0; i<n; ++i)
+      {
+        if (!ei_isMuchSmallerThan(value(i), reference, epsilon))
+        {
+          value(k) = value(i);
+          index(k) = index(i);
+          ++k;
+        }
+      }
+      resize(k,0);
+    }
 
   protected:
 
diff --git a/Eigen/src/Sparse/DynamicSparseMatrix.h b/Eigen/src/Sparse/DynamicSparseMatrix.h
index 30e54a981..5e7ce9327 100644
--- a/Eigen/src/Sparse/DynamicSparseMatrix.h
+++ b/Eigen/src/Sparse/DynamicSparseMatrix.h
@@ -74,7 +74,7 @@ class DynamicSparseMatrix
     std::vector<CompressedStorage<Scalar> > m_data;
 
   public:
-
+  
     inline int rows() const { return IsRowMajor ? outerSize() : m_innerSize; }
     inline int cols() const { return IsRowMajor ? m_innerSize : outerSize(); }
     inline int innerSize() const { return m_innerSize; }
@@ -102,8 +102,6 @@ class DynamicSparseMatrix
       return m_data[outer].atWithInsertion(inner);
     }
 
-  public:
-
     class InnerIterator;
 
     inline void setZero()
@@ -176,6 +174,12 @@ class DynamicSparseMatrix
 
     /** Does nothing. Provided for compatibility with SparseMatrix. */
     inline void endFill() {}
+    
+    void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
+    {
+      for (int j=0; j<outerSize(); ++j)
+        m_data[j].prune(reference,epsilon);
+    }
 
     /** Resize the matrix without preserving the data (the matrix is set to zero)
       */
diff --git a/Eigen/src/Sparse/SparseBlock.h b/Eigen/src/Sparse/SparseBlock.h
index 48be3754f..cbec36e61 100644
--- a/Eigen/src/Sparse/SparseBlock.h
+++ b/Eigen/src/Sparse/SparseBlock.h
@@ -78,6 +78,40 @@ public:
   }
 };
 
+/** \returns the i-th row of the matrix \c *this. For row-major matrix only. */
+template<typename Derived>
+SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i)
+{
+  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  return innerVector(i);
+}
+
+/** \returns the i-th row of the matrix \c *this. For row-major matrix only. 
+  * (read-only version) */
+template<typename Derived>
+const SparseInnerVector<Derived> SparseMatrixBase<Derived>::row(int i) const
+{
+  EIGEN_STATIC_ASSERT(IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  return innerVector(i);
+}
+
+/** \returns the i-th column of the matrix \c *this. For column-major matrix only. */
+template<typename Derived>
+SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i)
+{
+  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  return innerVector(i);
+}
+
+/** \returns the i-th column of the matrix \c *this. For column-major matrix only. 
+  * (read-only version) */
+template<typename Derived>
+const SparseInnerVector<Derived> SparseMatrixBase<Derived>::col(int i) const
+{
+  EIGEN_STATIC_ASSERT(!IsRowMajor,THIS_METHOD_IS_ONLY_FOR_ROW_MAJOR_MATRICES);
+  return innerVector(i);
+}
+
 /** \returns the \a outer -th column (resp. row) of the matrix \c *this if \c *this
   * is col-major (resp. row-major).
   */
diff --git a/Eigen/src/Sparse/SparseMatrix.h b/Eigen/src/Sparse/SparseMatrix.h
index 2ce466428..e0c56beb0 100644
--- a/Eigen/src/Sparse/SparseMatrix.h
+++ b/Eigen/src/Sparse/SparseMatrix.h
@@ -69,11 +69,11 @@ class SparseMatrix
     int* m_outerIndex;
     CompressedStorage<Scalar> m_data;
 
-
   public:
 
     inline int rows() const { return IsRowMajor ? m_outerSize : m_innerSize; }
     inline int cols() const { return IsRowMajor ? m_innerSize : m_outerSize; }
+    
     inline int innerSize() const { return m_innerSize; }
     inline int outerSize() const { return m_outerSize; }
     inline int innerNonZeros(int j) const { return m_outerIndex[j+1]-m_outerIndex[j]; }
@@ -230,6 +230,28 @@ class SparseMatrix
         ++i;
       }
     }
+    
+    void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
+    {
+      int k = 0;
+      for (int j=0; j<m_outerSize; ++j)
+      {
+        int previousStart = m_outerIndex[j];
+        m_outerIndex[j] = k;
+        int end = m_outerIndex[j+1];
+        for (int i=previousStart; i<end; ++i)
+        {
+          if (!ei_isMuchSmallerThan(m_data.value(i), reference, epsilon))
+          {
+            m_data.value(k) = m_data.value(i);
+            m_data.index(k) = m_data.index(i);
+            ++k;
+          }
+        }
+      }
+      m_outerIndex[m_outerSize] = k;
+      m_data.resize(k,0);
+    }
 
     void resize(int rows, int cols)
     {
diff --git a/Eigen/src/Sparse/SparseMatrixBase.h b/Eigen/src/Sparse/SparseMatrixBase.h
index 32011afaa..46c0c546a 100644
--- a/Eigen/src/Sparse/SparseMatrixBase.h
+++ b/Eigen/src/Sparse/SparseMatrixBase.h
@@ -327,6 +327,10 @@ template<typename Derived> class SparseMatrixBase
     // void transposeInPlace();
     const AdjointReturnType adjoint() const { return conjugate()/*.nestByValue()*/; }
 
+    SparseInnerVector<Derived> row(int i);
+    const SparseInnerVector<Derived> row(int i) const;
+    SparseInnerVector<Derived> col(int j);
+    const SparseInnerVector<Derived> col(int j) const;
     SparseInnerVector<Derived> innerVector(int outer);
     const SparseInnerVector<Derived> innerVector(int outer) const;
 
diff --git a/Eigen/src/Sparse/SparseVector.h b/Eigen/src/Sparse/SparseVector.h
index ff9d27aa7..9af565351 100644
--- a/Eigen/src/Sparse/SparseVector.h
+++ b/Eigen/src/Sparse/SparseVector.h
@@ -67,20 +67,19 @@ class SparseVector
     CompressedStorage<Scalar> m_data;
     int m_size;
 
-
   public:
 
-    inline int rows() const { return IsColVector ? m_size : 1; }
-    inline int cols() const { return IsColVector ? 1 : m_size; }
-    inline int innerSize() const { return m_size; }
-    inline int outerSize() const { return 1; }
-    inline int innerNonZeros(int j) const { ei_assert(j==0); return m_size; }
+    EIGEN_STRONG_INLINE int rows() const { return IsColVector ? m_size : 1; }
+    EIGEN_STRONG_INLINE int cols() const { return IsColVector ? 1 : m_size; }
+    EIGEN_STRONG_INLINE int innerSize() const { return m_size; }
+    EIGEN_STRONG_INLINE int outerSize() const { return 1; }
+    EIGEN_STRONG_INLINE int innerNonZeros(int j) const { ei_assert(j==0); return m_size; }
 
-    inline const Scalar* _valuePtr() const { return &m_data.value(0); }
-    inline Scalar* _valuePtr() { return &m_data.value(0); }
+    EIGEN_STRONG_INLINE const Scalar* _valuePtr() const { return &m_data.value(0); }
+    EIGEN_STRONG_INLINE Scalar* _valuePtr() { return &m_data.value(0); }
 
-    inline const int* _innerIndexPtr() const { return &m_data.index(0); }
-    inline int* _innerIndexPtr() { return &m_data.index(0); }
+    EIGEN_STRONG_INLINE const int* _innerIndexPtr() const { return &m_data.index(0); }
+    EIGEN_STRONG_INLINE int* _innerIndexPtr() { return &m_data.index(0); }
 
     inline Scalar coeff(int row, int col) const
     {
@@ -145,6 +144,11 @@ class SparseVector
       m_data.value(id+1) = 0;
       return m_data.value(id+1);
     }
+    
+    void prune(Scalar reference, RealScalar epsilon = precision<RealScalar>())
+    {
+      m_data.prune(reference,epsilon);
+    }
 
     void resize(int newSize)
     {
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp
index addd40f9e..5ae038ef8 100644
--- a/test/sparse_basic.cpp
+++ b/test/sparse_basic.cpp
@@ -323,14 +323,49 @@ template<typename SparseMatrixType> void sparse_basic(const SparseMatrixType& re
     VERIFY_IS_APPROX(x=mLo.template marked<LowerTriangular|SelfAdjoint>()*b, refX=refS*b);
     VERIFY_IS_APPROX(x=mS.template marked<SelfAdjoint>()*b, refX=refS*b);
   }
+  
+  // test prune
+  {
+    SparseMatrixType m2(rows, rows);
+    DenseMatrix refM2(rows, rows);
+    refM2.setZero();
+    int countFalseNonZero = 0;
+    int countTrueNonZero = 0;
+    m2.startFill();
+    for (int j=0; j<m2.outerSize(); ++j)
+      for (int i=0; i<m2.innerSize(); ++i)
+      {
+        float x = ei_random<float>(0,1);
+        if (x<0.1)
+        {
+          // do nothing
+        }
+        else if (x<0.5)
+        {
+          countFalseNonZero++;
+          m2.fill(i,j) = Scalar(0);
+        }
+        else
+        {
+          countTrueNonZero++;
+          m2.fill(i,j) = refM2(i,j) = Scalar(1);
+        }
+      }
+    m2.endFill();
+    VERIFY(countFalseNonZero+countTrueNonZero == m2.nonZeros());
+    VERIFY_IS_APPROX(m2, refM2);
+    m2.prune(1);
+    VERIFY(countTrueNonZero==m2.nonZeros());
+    VERIFY_IS_APPROX(m2, refM2);
+  }
 }
 
 void test_sparse_basic()
 {
   for(int i = 0; i < g_repeat; i++) {
-//     CALL_SUBTEST( sparse_basic(SparseMatrix<double>(8, 8)) );
-//     CALL_SUBTEST( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) );
-//     CALL_SUBTEST( sparse_basic(SparseMatrix<double>(33, 33)) );
+    CALL_SUBTEST( sparse_basic(SparseMatrix<double>(8, 8)) );
+    CALL_SUBTEST( sparse_basic(SparseMatrix<std::complex<double> >(16, 16)) );
+    CALL_SUBTEST( sparse_basic(SparseMatrix<double>(33, 33)) );
     
     CALL_SUBTEST( sparse_basic(DynamicSparseMatrix<double>(8, 8)) );
   }