Use more efficient colmatrix-based reduction

- Use colmatrix as boundary matrix.
- Use faster set_symmetric_difference() for colmatrix::add_column().
- Move diagram sanity checks to unit tests.
- Make lowestones a col-index vector indexed by row-indices.
- Test lowestones in unit tests, which are independt of the reduction
  algorithm.

diff --git a/include/libstick-0.1/booleanmatrix.h b/include/libstick-0.1/booleanmatrix.h
index 8626f53d58eca6c755737b29f16e37bebd316399..2b64a4ecb79bcd4200e3e26ab3c16bfa2fd10025 100644 (file)
--- a/include/libstick-0.1/booleanmatrix.h
+++ b/include/libstick-0.1/booleanmatrix.h
@@ -15,6 +15,10 @@
 namespace libstick {
 
 
 namespace libstick {
 
 

+template<class T>
+std::ostream& operator<<(std::ostream &, const std::vector<T> &);
+
+

 /** The base class of boolean_colmatrix and boolean_colrowmatrix which implements
  * the common logic of both. */
 template<class IT, class D>
 /** The base class of boolean_colmatrix and boolean_colrowmatrix which implements
  * the common logic of both. */
 template<class IT, class D>


@@ -32,11 +36,26 @@ class boolean_colmatrix_base {
         }
 
     public:
         }
 
     public:

-        /** Get height resp. width of the matrix. */
+        /** A casting constructor for any colmatrix with same entries type. */
+        template<class D2>
+        boolean_colmatrix_base(const boolean_colmatrix_base<IT, D2> &mat) :
+            cols(mat.get_columns()) {
+        }
+
+        /** Get width of the matrix. */

         size_t width() const {
             return cols.size();
         }
 
         size_t width() const {
             return cols.size();
         }
 

+        /** Get height of the matrix, i.e. maximum row-index + 1 among all columns. */
+        size_t height() const {
+            IT h = 0;
+            for (unsigned c=0; c < width(); ++c)
+                if (cols[c].size() > 0)
+                    h = std::max(h, cols[c].back());
+            return h+1;
+        }
+

         /** Get the matrix entry at row 'r' and column 'c'. */
         bool get(index_type r, index_type c) const {
             assert(c < width());
         /** Get the matrix entry at row 'r' and column 'c'. */
         bool get(index_type r, index_type c) const {
             assert(c < width());


@@ -61,6 +80,11 @@ class boolean_colmatrix_base {
             return cols[c];
         }
 
             return cols[c];
         }
 

+        /** Get all columns */
+        const std::vector<column_type>& get_columns() const {
+            return cols;
+        }
+

         /** Add the column-vector 'col' to the c-th column. Note that 'col'
          * actually contains the list of row-indices that are 1. */
         void add_column(index_type c, const column_type &col) {
         /** Add the column-vector 'col' to the c-th column. Note that 'col'
          * actually contains the list of row-indices that are 1. */
         void add_column(index_type c, const column_type &col) {


@@ -72,12 +96,14 @@ class boolean_colmatrix_base {
         }
 
         /** Two matrices are equal iff they have the same entries */
         }
 
         /** Two matrices are equal iff they have the same entries */

-        bool operator==(const boolean_colmatrix_base<IT, D> &m) const {
-            return cols == m.cols;
+        template<class D2>
+        bool operator==(const boolean_colmatrix_base<IT, D2> &m) const {
+            return cols == m.get_columns();

         }
 
         /** Two matrices are equal iff they have the same entries */
         }
 
         /** Two matrices are equal iff they have the same entries */

-        bool operator!=(const boolean_colmatrix_base<IT, D> &m) const {
+        template<class D2>
+        bool operator!=(const boolean_colmatrix_base<IT, D2> &m) const {

             return !(*this == m);
         }
 
             return !(*this == m);
         }
 


@@ -131,7 +157,7 @@ class boolean_colmatrix_base {
 #endif
         }
 
 #endif
         }
 

-    private:
+    protected:

         /** The matrix is the set of columns. */
         std::vector<column_type> cols;
 
         /** The matrix is the set of columns. */
         std::vector<column_type> cols;
 


@@ -151,6 +177,7 @@ class boolean_colmatrix : public boolean_colmatrix_base<IT, boolean_colmatrix<IT
     public:
         typedef IT index_type;
         typedef boolean_colmatrix_base<IT, boolean_colmatrix<IT> > base;
     public:
         typedef IT index_type;
         typedef boolean_colmatrix_base<IT, boolean_colmatrix<IT> > base;

+        typedef typename base::column_type column_type;

 
         /** Create a matrix with 'width' columns, initalized with zero entries. */
         boolean_colmatrix(size_t columns) :
 
         /** Create a matrix with 'width' columns, initalized with zero entries. */
         boolean_colmatrix(size_t columns) :


@@ -161,6 +188,37 @@ class boolean_colmatrix : public boolean_colmatrix_base<IT, boolean_colmatrix<IT
         void _set(index_type r, index_type c, bool value) {
             base::_set(r, c, value);
         }
         void _set(index_type r, index_type c, bool value) {
             base::_set(r, c, value);
         }

+
+        /** A faster implementation of boolean_colmatrix_base::add_column().
+         * Assumes that 'col' is sorted. */
+        void add_column(index_type c, const column_type &col) {
+            assert(c < base::width());
+
+#ifndef NDEBUG
+            for (unsigned i=1; i < col.size(); ++i)
+                assert(col[i-1] < col[i]);
+#endif
+
+            // The original column
+            column_type &orig_col = base::cols[c];
+
+            // Make target column large enough
+            const size_t maxsize = orig_col.size() + col.size();
+            if (tcol.size() < maxsize)
+                tcol.resize(maxsize);
+
+            // Compute symmetric difference
+            typename column_type::iterator it = std::set_symmetric_difference(
+                    orig_col.begin(), orig_col.end(), col.begin(), col.end(), tcol.begin());
+
+            // Copy back to the original column
+            orig_col.resize(it - tcol.begin());
+            std::copy(tcol.begin(), it, orig_col.begin());
+        }
+
+    private:
+        /** A temporary container to speed up add_column() */
+        column_type tcol;

 };
 
 
 };
 
 


@@ -182,7 +240,7 @@ class boolean_rowmatrix_base {
         }
 
     public:
         }
 
     public:

-        /** Get height resp. width of the matrix. */
+        /** Get height of the matrix. */

         size_t height() const {
             return rows.size();
         }
         size_t height() const {
             return rows.size();
         }


@@ -263,7 +321,7 @@ class boolean_rowmatrix_base {
 #endif
         }
 
 #endif
         }
 

-    private:
+    protected:

         derived* get_derived() {
             return static_cast<derived*>(this);
         }
         derived* get_derived() {
             return static_cast<derived*>(this);
         }


@@ -314,6 +372,23 @@ class boolean_colrowmatrix : public boolean_colmatrix_base<IT, boolean_colrowmat
             rowbase(size) {
         }
 
             rowbase(size) {
         }
 

+        /** Casting a colmatrix into a colrow matrix */
+        template<class D>
+        boolean_colrowmatrix(const boolean_colmatrix_base<IT, D>& mat) :
+            colbase(std::max(mat.width(), mat.height())),
+            rowbase(std::max(mat.width(), mat.height())) {
+            for (unsigned c=0; c < mat.width(); ++c) {
+                const typename colbase::column_type &col = mat.get_column(c);
+                for (unsigned i=0; i < col.size(); ++i)
+                    set(col[i], c, true);
+            }
+            for (unsigned r=0; r < size(); ++r) {
+                const typename rowbase::row_type &row = rowbase::get_row(r);
+                for (unsigned i=0; i < row.size(); ++i)
+                    assert(get(r, row[i]) == true);
+            }
+        }
+

         /** Override implementation. */
         void _set(index_type r, index_type c, bool value) {
             colbase::_set(r, c, value);
         /** Override implementation. */
         void _set(index_type r, index_type c, bool value) {
             colbase::_set(r, c, value);


@@ -352,6 +427,12 @@ class boolean_colrowmatrix : public boolean_colmatrix_base<IT, boolean_colrowmat
             return colbase::operator==(m);
         }
 
             return colbase::operator==(m);
         }
 

+        /** Two matrices are equal iff they have the same entries */
+        template<class D>
+        bool operator==(const boolean_colmatrix_base<IT, D> &m) const {
+            return colbase::operator==(m);
+        }
+

         /** Two matrices are equal iff they have the same entries */
         bool operator!=(const boolean_colrowmatrix<IT> &m) const {
             return !(*this == m);
         /** Two matrices are equal iff they have the same entries */
         bool operator!=(const boolean_colrowmatrix<IT> &m) const {
             return !(*this == m);


@@ -524,6 +605,21 @@ std::ostream& operator<<(std::ostream &os, boolean_colmatrix_base<IT, D> &mat) {
     return os << m;
 }
 
     return os << m;
 }
 

+template<class T>
+std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
+    os << "[";
+
+    typename std::vector<T>::const_iterator it = vec.begin();
+    while ( it != vec.end()) {
+        os << *it;
+        if (++it != vec.end())
+            os << " ";
+    }
+
+    os << "]";
+    return os;
+}
+

 }
 
 #endif
 }
 
 #endif

diff --git a/include/libstick-0.1/persistence.h b/include/libstick-0.1/persistence.h
index 0a69be22cb0acf7c974595e4a91a3037aac69c2b..6df70abc7c8345b0f11a86d168fdf3bcd9f4f0f2 100644 (file)
--- a/include/libstick-0.1/persistence.h
+++ b/include/libstick-0.1/persistence.h
@@ -13,10 +13,6 @@
 namespace libstick {
 
 
 namespace libstick {
 
 

-template<class T>
-std::ostream& operator<<(std::ostream &, const std::vector<T> &);
-
-

 /** Persistence computers various persistent homology information on a
  * simplex_order of a simplicial_complex. */
 template<int MAXDIM, class IT, class VT>
 /** Persistence computers various persistent homology information on a
  * simplex_order of a simplicial_complex. */
 template<int MAXDIM, class IT, class VT>


@@ -25,8 +21,9 @@ class persistence {
     public:
         typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
         typedef typename scomplex::simplex_order simplex_order;
     public:
         typedef simplicial_complex<MAXDIM, IT, VT> scomplex;
         typedef typename scomplex::simplex_order simplex_order;

-        typedef typename simplex_order::boundary_matrix boundary_matrix;
+        typedef boolean_colmatrix<IT> boundary_matrix;

         typedef boolean_colmatrix<IT> transformation_matrix;
         typedef boolean_colmatrix<IT> transformation_matrix;

+        typedef std::vector<IT> lowestones_matrix;

 
         /** A point in a persistence diagram comprises birth- and death-index.
          * These are the indices of the simplices of 'order' that give birth
 
         /** A point in a persistence diagram comprises birth- and death-index.
          * These are the indices of the simplices of 'order' that give birth


@@ -107,18 +104,45 @@ class persistence {
 
             bm = order.get_boundary_matrix();
             rbm = bm;
 
             bm = order.get_boundary_matrix();
             rbm = bm;

-            tm = create_unit_matrix<transformation_matrix>(bm.size());
-
-            reduce_boundary_matrix(rbm, tm);
+            tm = create_unit_matrix<transformation_matrix>(bm.width());
+            lowestones = lowestones_matrix(bm.width());
+
+            // Make every column reduced, i.e., it is a zero-column or contains only one 1.
+            for (unsigned c=0; c < rbm.width(); ++c) {
+                //if (c % 100 == 0)
+                //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
+                // Reduce as long as we need to reduce
+                while (rbm.get_column(c).size() > 0) {
+                    // (r, c) is the lowest one of col
+                    const IT r = rbm.get_column(c).back();
+                    // This column contains its lowest one on the same row
+                    const IT c2 = lowestones[r];
+                    assert(c2 < c);
+
+                    // There is no valid c2, hence we have completely reduced
+                    if (c2 == 0)
+                        break;
+                    // Reduce col by column c2
+                    else {
+                        assert(rbm.get(r, c));
+                        rbm.add_column(c, rbm.get_column(c2));
+                        tm.add_column(c, tm.get_column(c2));
+                        assert(!rbm.get(r, c));
+                    }
+                }

 
 

-            lowestones = boundary_matrix(rbm.size());
-            for (unsigned c=0; c < rbm.size(); ++c)
-                if (rbm.get_column(c).size() > 0)
-                    lowestones.set(rbm.get_column(c).back(), c, true);
+                // A lowest one remained, recall it
+                if (rbm.get_column(c).size() > 0) {
+                    const IT r = rbm.get_column(c).back();
+                    assert(lowestones[r] == 0);
+                    assert(r < c);
+                    lowestones[r] = c;
+                }
+            }

         }
 
         /** Get the lowest-one matrix of 'order'. */
         }
 
         /** Get the lowest-one matrix of 'order'. */

-        const boundary_matrix& get_lowestones_matrix() const {
+        const lowestones_matrix& get_lowestones_matrix() const {

             assert(done_matrices);
             return lowestones;
         }
             assert(done_matrices);
             return lowestones;
         }


@@ -147,7 +171,7 @@ class persistence {
         void interpret_reduction(std::ostream &os) const {
             assert(done_matrices);
 
         void interpret_reduction(std::ostream &os) const {
             assert(done_matrices);
 

-            for (unsigned c=0; c < bm.size(); ++c) {
+            for (unsigned c=0; c < bm.width(); ++c) {

                 os << c << ". inserting ";
 
                 switch (bm.get_column(c).size()) {
                 os << c << ". inserting ";
 
                 switch (bm.get_column(c).size()) {


@@ -174,7 +198,7 @@ class persistence {
                 os << std::endl;
                 os << "  boundary: " << bm.get_column(c) << std::endl;
 
                 os << std::endl;
                 os << "  boundary: " << bm.get_column(c) << std::endl;
 

-                typename boolean_colrowmatrix<IT>::colbase::column_type col = rbm.get_column(c);
+                typename boolean_colmatrix<IT>::column_type col = rbm.get_column(c);

                 if (col.size() == 0) {
                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
                 } else {
                 if (col.size() == 0) {
                     os << "  \033[1;32mbirth\033[0;m of a cycle: " << tm.get_column(c) << std::endl;
                 } else {


@@ -210,7 +234,8 @@ class persistence {
                 diagrams[d] = diagram();
 
             for (unsigned c=0; c < lowestones.size(); ++c) {
                 diagrams[d] = diagram();
 
             for (unsigned c=0; c < lowestones.size(); ++c) {

-                if (lowestones.get_column(c).size() == 0) {
+                // A cycle was born
+                if (rbm.get_column(c).size() == 0) {

                     const int dim = order.get_simplex(c).dim;
 
                     // Create a diagram point
                     const int dim = order.get_simplex(c).dim;
 
                     // Create a diagram point


@@ -218,8 +243,8 @@ class persistence {
                     p.birth = c;
                     p.death = 0;
                     // If the class actually dies, get the index of the killing simplex
                     p.birth = c;
                     p.death = 0;
                     // If the class actually dies, get the index of the killing simplex

-                    if (lowestones.get_row(c).size() > 0)
-                        p.death = lowestones.get_row(c).back();
+                    if (lowestones[c] > 0)
+                        p.death = lowestones[c];

 
                     _get_persistence_diagram(dim).births.push_back(p);
 #ifndef NDEBUG
 
                     _get_persistence_diagram(dim).births.push_back(p);
 #ifndef NDEBUG


@@ -234,18 +259,12 @@ class persistence {
             // Shakespeare principle: A simplex either gives birth or kills, be
             // or not be, tertium non datur.
             for (unsigned c=0; c < lowestones.size(); ++c) {
             // Shakespeare principle: A simplex either gives birth or kills, be
             // or not be, tertium non datur.
             for (unsigned c=0; c < lowestones.size(); ++c) {

-                if (lowestones.get_column(c).size() != 0)
+                // A class died with c
+                if (rbm.get_column(c).size() != 0)

                     assert(c == 0 || deaths.count(c) > 0);
                 else
                     assert(births.count(c) > 0);
             }
                     assert(c == 0 || deaths.count(c) > 0);
                 else
                     assert(births.count(c) > 0);
             }

-            // Even Jesus died after his birth. Homology classes are nothing better.
-            for (int d=-1; d <= MAXDIM; ++d) {
-                const diagram &dia = get_persistence_diagram(d);
-                for (unsigned i=0; i < dia.births.size(); ++i)
-                    if (dia.births[i].death > 0)
-                        assert(dia.births[i].birth < dia.births[i].death);
-            }

 #endif
         }
 
 #endif
         }
 


@@ -288,9 +307,9 @@ class persistence {
         /** The reduced boundary matrix of 'order' */
         boundary_matrix rbm;
         /** Only the lowest ones per column of the reduced boundary matrix.
         /** The reduced boundary matrix of 'order' */
         boundary_matrix rbm;
         /** Only the lowest ones per column of the reduced boundary matrix.

-         * Each one at (r, c) corresponds to the death of a class when
-         * inserting simplex c, which was born when inserting simplex. */
-        boundary_matrix lowestones;
+         * lowestones[i] contains the column-index at which the lowest one at
+         * row i appears, and zero if no such column exists. */
+        lowestones_matrix lowestones;

         /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
         transformation_matrix tm;
 
         /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
         transformation_matrix tm;
 


@@ -306,51 +325,45 @@ class persistence {
  * invariant.  Hence, the resulting b is equal to the product of the boundary
  * matrix times v. */
 template<class IT, class D>
  * invariant.  Hence, the resulting b is equal to the product of the boundary
  * matrix times v. */
 template<class IT, class D>

-void reduce_boundary_matrix(boolean_colrowmatrix<IT> &b, boolean_colmatrix_base<IT, D> &v) {
-    assert(b.size() == v.width());
+void reduce_boundary_matrix(boolean_colmatrix_base<IT, D> &b, boolean_colmatrix_base<IT, D> &v) {
+    assert(b.width() == v.width());
+
+    // lowestones[i] gives the column whose lowest 1 is at row i. If it contains
+    // zero, there is no such column.
+    std::vector<IT> lowestones(b.width());

 
     // Make every column reduced, i.e., it is a zero-column or contains only one 1.
 
     // Make every column reduced, i.e., it is a zero-column or contains only one 1.

-    for (unsigned c=0; c < b.size(); ++c) {
-        const typename boolean_colrowmatrix<IT>::colbase::column_type &col = b.get_column(c);
-        if (col.size() == 0)
-            continue;
-
-        // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
-        IT r = col.back();
-        assert(b.get(r, c));
-
-        // Get a copy of the r-th row
-        typename boolean_colrowmatrix<IT>::rowbase::row_type row(b.get_row(r).size());
-        copy(b.get_row(r).begin(), b.get_row(r).end(), row.begin());
-        assert(row.size() >= 1);
-
-        // Get rid of 1s at that row right of column c.
-        typename boolean_colrowmatrix<IT>::row_type::const_iterator it = lower_bound(row.begin(), row.end(), c);
-        for(; it != row.end(); ++it) {
-            if (*it <= c)
-                continue;
-
-            assert(b.get(r, *it));
-            b.add_column(*it, col);
-            v.add_column(*it, v.get_column(c));
-            assert(!b.get(r, *it));
+    for (unsigned c=0; c < b.width(); ++c) {
+        //if (c % 100 == 0)
+            //std::cout << "c = " << c << " (" << (100.0*float(c)/b.width()) << " %)" << std::endl;
+        // Reduce as long as we need to reduce
+        while (b.get_column(c).size() > 0) {
+            // (r, c) is the lowest one of col
+            const IT r = b.get_column(c).back();
+            // This column contains its lowest one on the same row
+            const IT c2 = lowestones[r];
+            assert(c2 < c);
+
+            // There is no valid c2, hence we have completely reduced
+            if (c2 == 0)
+                break;
+            // Reduce col by column c2
+            else {
+                assert(b.get(r, c));
+                b.add_column(c, b.get_column(c2));
+                v.add_column(c, v.get_column(c2));
+                assert(!b.get(r, c));
+            }

         }
         }

-    }
-}
-
-template<class T>
-std::ostream& operator<<(std::ostream& os, const std::vector<T> &vec) {
-    os << "[";

 
 

-    typename std::vector<T>::const_iterator it = vec.begin();
-    while ( it != vec.end()) {
-        os << *it;
-        if (++it != vec.end())
-            os << " ";
+        // A lowest one remained, recall it
+        if (b.get_column(c).size() > 0) {
+            const IT r = b.get_column(c).back();
+            assert(lowestones[r] == 0);
+            assert(r < c);
+            lowestones[r] = c;
+        }

     }
     }

-
-    os << "]";
-    return os;

 }
 
 
 }
 
 

diff --git a/include/libstick-0.1/simplicialcomplex.h b/include/libstick-0.1/simplicialcomplex.h
index 1183a92d699dfa896853e8148abb988f3ca06ed4..973d632e548d581fc3b893df6582644f5ce6e021 100644 (file)
--- a/include/libstick-0.1/simplicialcomplex.h
+++ b/include/libstick-0.1/simplicialcomplex.h
@@ -92,7 +92,7 @@ class simplicial_complex {
         class simplex_order {
 
             public:
         class simplex_order {
 
             public:

-                typedef boolean_colrowmatrix<IT> boundary_matrix;
+                typedef boolean_colmatrix<IT> boundary_matrix;

 
                 /** Create a standard order of the complex c, i.e., the identity permutation. */
                 simplex_order(const simplcompltype &c) :
 
                 /** Create a standard order of the complex c, i.e., the identity permutation. */
                 simplex_order(const simplcompltype &c) :

diff --git a/tests/booleanmatrix.h b/tests/booleanmatrix.h
index 71533ffbaf392ef44841d778196e3aaea8390e95..6cbf53847a62a93724a0ed91821e4a5577f7ea10 100644 (file)
--- a/tests/booleanmatrix.h
+++ b/tests/booleanmatrix.h
@@ -63,10 +63,10 @@ class boolean_matrix_TestSuite: public Test::Suite {
             for (unsigned c=0; c < size; ++c)
                 TEST_ASSERT(mat.get_column(c).size() == 0);
 
             for (unsigned c=0; c < size; ++c)
                 TEST_ASSERT(mat.get_column(c).size() == 0);
 

+            col.push_back(0);

             col.push_back(3);
             col.push_back(3);

-            col.push_back(14);

             col.push_back(8);
             col.push_back(8);

-            col.push_back(0);
+            col.push_back(14);

 
             // Add column and test for values
             mat.add_column(5, col);
 
             // Add column and test for values
             mat.add_column(5, col);

diff --git a/tests/image.h b/tests/image.h
index ff7817af6193cc4f2744b1c542c82e5f55b41f27..8633a56b5859e9da8dc91fd41d4406a57a0933fb 100644 (file)
--- a/tests/image.h
+++ b/tests/image.h
@@ -49,7 +49,7 @@ class image_TestSuite: public Test::Suite {
             sorder so(s);
             assert(so.is_filtration());
 
             sorder so(s);
             assert(so.is_filtration());
 

-            typedef sorder::boundary_matrix bmatrix;
+            typedef boolean_colrowmatrix<uint32_t> bmatrix;

             bmatrix bm = so.get_boundary_matrix();
 
             // Check for the right vertex incidences
             bmatrix bm = so.get_boundary_matrix();
 
             // Check for the right vertex incidences

diff --git a/tests/persistence.h b/tests/persistence.h
index 27b01f2f18c73aa0dbfc3ef79e85bc8d583f1419..832457e505fd3d5a6ab81310e1913a741f0b12e2 100644 (file)
--- a/tests/persistence.h
+++ b/tests/persistence.h
@@ -15,6 +15,7 @@ class persistence_TestSuite: public Test::Suite {
         typedef simplicial_complex<3, uint32_t, double> scomplex;
         typedef persistence<3, uint32_t, double> pers;
         typedef pers::boundary_matrix bm;
         typedef simplicial_complex<3, uint32_t, double> scomplex;
         typedef persistence<3, uint32_t, double> pers;
         typedef pers::boundary_matrix bm;

+        typedef pers::lowestones_matrix lom;

         typedef pers::transformation_matrix tm;
 
         bool setupcalled;
         typedef pers::transformation_matrix tm;
 
         bool setupcalled;


@@ -31,6 +32,7 @@ class persistence_TestSuite: public Test::Suite {
             TEST_ADD(persistence_TestSuite::test_matrix_reduction);
             TEST_ADD(persistence_TestSuite::test_betti_numbers);
             TEST_ADD(persistence_TestSuite::test_lowestones);
             TEST_ADD(persistence_TestSuite::test_matrix_reduction);
             TEST_ADD(persistence_TestSuite::test_betti_numbers);
             TEST_ADD(persistence_TestSuite::test_lowestones);

+            TEST_ADD(persistence_TestSuite::test_diagram);

         }
 
     protected:
         }
 
     protected:


@@ -170,48 +172,37 @@ class persistence_TestSuite: public Test::Suite {
             const bm &ringbm = ringp.get_boundary_matrix();
             const bm &ringrbm = ringp.get_reduced_boundary_matrix();
             const tm &ringtm = ringp.get_transformation_matrix();
             const bm &ringbm = ringp.get_boundary_matrix();
             const bm &ringrbm = ringp.get_reduced_boundary_matrix();
             const tm &ringtm = ringp.get_transformation_matrix();

-            TEST_ASSERT(ringrbm == ringbm * ringtm);
+            TEST_ASSERT(ringrbm == boolean_colrowmatrix<uint32_t>(ringbm) * ringtm);
+

 
             pers torusp(otorus);
             torusp.compute_matrices();
             const bm &torusbm = torusp.get_boundary_matrix();
             const bm &torusrbm = torusp.get_reduced_boundary_matrix();
             const tm &torustm = torusp.get_transformation_matrix();
 
             pers torusp(otorus);
             torusp.compute_matrices();
             const bm &torusbm = torusp.get_boundary_matrix();
             const bm &torusrbm = torusp.get_reduced_boundary_matrix();
             const tm &torustm = torusp.get_transformation_matrix();

-            TEST_ASSERT(torusrbm == torusbm * torustm);
-
-            uint32_t torusrbme_coords[][2] = {
-                    {0, 1}, {2, 7}, {3, 7}, {1, 8}, {2, 8}, {4, 10}, {5, 10}, {4, 11}, {6, 11}, {1, 13},
-                    {4, 13}, {10, 15}, {13, 15}, {14, 15}, {9, 19}, {10, 19}, {13, 19}, {17, 19}, {7, 20},
-                    {16, 20}, {18, 20}, {7, 21}, {9, 21}, {10, 21}, {11, 21}, {13, 21}, {16, 21}, {12, 23},
-                    {13, 23}, {22, 23}, {7, 24}, {8, 24}, {9, 24}, {10, 24}, {11, 24}, {12, 24}, {25, 28},
-                    {26, 28}, {25, 29}, {27, 29}, {1, 31}, {25, 31}, {28, 33}, {31, 33}, {32, 33}, {30, 35},
-                    {31, 35}, {34, 35}, {8, 37}, {30, 37}, {31, 37}, {36, 37}, {9, 40}, {28, 40}, {31, 40},
-                    {39, 40}, {9, 41}, {28, 41}, {29, 41}, {31, 41}, {38, 41}, {7, 42}, {8, 42}, {9, 42},
-                    {28, 42}, {29, 42}, {30, 42}, {7, 49}, {8, 49}, {9, 49}, {28, 49}, {29, 49}, {43, 49},
-                    {47, 49}, {10, 50}, {11, 50}, {28, 50}, {29, 50}, {43, 50}, {45, 50}, {29, 51}, {44, 51},
-                    {48, 51}, {10, 52}, {28, 52}, {43, 52}, {44, 52}, {28, 53}, {43, 53}, {46, 53 }
-            };
-            bm torusrbme(torusbm.size());
-            torusrbme.set_all(torusrbme_coords, sizeof(torusrbme_coords)/(2*sizeof(uint32_t)), true);
-            TEST_ASSERT(torusrbme == torusrbm);
+            TEST_ASSERT(torusrbm == boolean_colrowmatrix<uint32_t>(torusbm) * torustm);
+
+            const lom &toruslo = torusp.get_lowestones_matrix();
+            uint32_t torusloe_coords[] = {1, 0, 8, 7, 13, 10, 11, 0, 0, 0, 0,
+                0, 24, 0, 15, 0, 21, 19, 20, 0, 0, 0, 23, 0, 0, 31, 28, 29, 0,
+                0, 42, 0, 33, 0, 35, 0, 37, 0, 41, 40, 0, 0, 0, 0, 52, 50, 53,
+                49, 51, 0, 0, 0, 0, 0, 0};
+            lom torusloe(torusloe_coords, torusloe_coords + sizeof(torusloe_coords)/sizeof(uint32_t));
+            TEST_ASSERT(torusloe == toruslo);
+

 
             pers ballp(oball);
             ballp.compute_matrices();
             const bm &ballbm = ballp.get_boundary_matrix();
             const bm &ballrbm = ballp.get_reduced_boundary_matrix();
             const tm &balltm = ballp.get_transformation_matrix();
 
             pers ballp(oball);
             ballp.compute_matrices();
             const bm &ballbm = ballp.get_boundary_matrix();
             const bm &ballrbm = ballp.get_reduced_boundary_matrix();
             const tm &balltm = ballp.get_transformation_matrix();

-            TEST_ASSERT(ballrbm == ballbm * balltm);
-
-            uint32_t ballrbme_coords[][2] = {
-                    {0, 1}, {1, 6}, {2, 6}, {1, 7}, {3, 7}, {1, 8}, {4, 8}, {1, 10}, {5, 10}, {6, 14},
-                    {10, 14}, {11, 14}, {6, 15}, {7, 15}, {10, 15}, {12, 15}, {6, 16}, {7, 16}, {8, 16},
-                    {10, 16}, {13, 16}, {6, 17}, {7, 17}, {8, 17}, {9, 17}, {6, 19}, {7, 19}, {18, 19},
-                    {14, 22}, {15, 22}, {19, 22}, {21, 22}, {14, 23}, {15, 23}, {16, 23}, {17, 23}, {19, 23},
-                    {20, 23 }
-            };
-            bm ballrbme(ballrbm.size());
-            ballrbme.set_all(ballrbme_coords, sizeof(ballrbme_coords)/(2*sizeof(uint32_t)), true);
-            TEST_ASSERT(ballrbme == ballrbm);
+            TEST_ASSERT(ballrbm == boolean_colrowmatrix<uint32_t>(ballbm) * balltm);
+
+            const lom &balllo = ballp.get_lowestones_matrix();
+            uint32_t ballloe_coords[] = {1, 0, 6, 7, 8, 10, 0, 0, 0, 17, 0, 14,
+                15, 16, 0, 0, 0, 0, 19, 0, 23, 22, 0, 0};
+            lom ballloe(ballloe_coords, ballloe_coords + sizeof(ballloe_coords)/sizeof(uint32_t));
+            TEST_ASSERT(ballloe == balllo);

 
             //std::cout << std::endl;
             //ballp.interpret_reduction(std::cout);
 
             //std::cout << std::endl;
             //ballp.interpret_reduction(std::cout);


@@ -293,26 +284,56 @@ class persistence_TestSuite: public Test::Suite {
         void test_lowestones_impl(const scomplex::simplex_order& so) {
             pers p(so);
             p.compute_matrices();
         void test_lowestones_impl(const scomplex::simplex_order& so) {
             pers p(so);
             p.compute_matrices();

-            const bm &lowestones = p.get_lowestones_matrix();
-
-            for (unsigned c=0; c < lowestones.size(); ++c)
-                assert(lowestones.get_column(c).size() <= 1);
-            for (unsigned r=0; r < lowestones.size(); ++r)
-                assert(lowestones.get_row(r).size() <= 1);
-            for (unsigned c=0; c < lowestones.size(); ++c) {
+            const lom &lowestones = p.get_lowestones_matrix();
+
+            // Check that lowestones contains every column-index only once,
+            // except the zeros.
+            lom cpy = lowestones;
+            sort(cpy.begin(), cpy.end());
+            lom::iterator it = cpy.begin();
+            // Skip the zeros
+            while (it != cpy.end() && *it == 0)
+                ++it;
+            // From now on, no column index appears twice
+            for (; it+1 != cpy.end(); ++it)
+                TEST_ASSERT(*it < *(it+1));
+
+            for (unsigned r=0; r < lowestones.size(); ++r) {

                 // If (r,c) is a one then
                 //   (i) a cycle dies with c --> row c is zero
                 //  (ii) a cycle is born with r --> column r is zero
                 //Hence
                 //   (i) the column r is a zero-column
                 //   (i) the row c is a zero-column
                 // If (r,c) is a one then
                 //   (i) a cycle dies with c --> row c is zero
                 //  (ii) a cycle is born with r --> column r is zero
                 //Hence
                 //   (i) the column r is a zero-column
                 //   (i) the row c is a zero-column

-                if (lowestones.get_column(c).size() > 0) {
-                    const uint32_t r = lowestones.get_column(c).back();
-                    assert(lowestones.get_column(r).size() == 0);
-                    assert(lowestones.get_row(c).size() == 0);
+
+                const uint32_t c = lowestones[r];
+                if (c > 0) {
+                    TEST_ASSERT(p.get_reduced_boundary_matrix().get(r,c) == true);
+                    TEST_ASSERT(p.get_reduced_boundary_matrix().get_column(r).size() == 0);
+                    TEST_ASSERT(lowestones[c] == 0);

                 }
             }
         }
                 }
             }
         }

+
+        void test_diagram() {
+            test_diagram_impl(oring);
+            test_diagram_impl(otorus);
+            test_diagram_impl(oball);
+        }
+
+        void test_diagram_impl(const scomplex::simplex_order& so) {
+            pers p(so);
+            p.compute_diagrams();
+
+            // Even Jesus died after his birth. Homology classes are nothing better.
+            for (int d=-1; d <= 3; ++d) {
+                const pers::diagram &dia = p.get_persistence_diagram(d);
+                for (unsigned i=0; i < dia.births.size(); ++i)
+                    if (dia.births[i].death > 0)
+                        assert(dia.births[i].birth < dia.births[i].death);
+            }
+
+        }

 };
 
 #endif
 };
 
 #endif