Move persistence checking code to unit tests

diff --git a/include/libstick-0.1/persistence.h b/include/libstick-0.1/persistence.h
index 3840bd2f8d65a68385948187860b98c241771ea5..0a69be22cb0acf7c974595e4a91a3037aac69c2b 100644 (file)
--- a/include/libstick-0.1/persistence.h
+++ b/include/libstick-0.1/persistence.h
@@ -23,7 +23,8 @@ template<int MAXDIM, class IT, class VT>
 class persistence {
 
     public:
-        typedef typename simplicial_complex<MAXDIM, IT, VT>::simplex_order simplex_order;
+        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> transformation_matrix;
 
@@ -109,32 +110,17 @@ class persistence {
             tm = create_unit_matrix<transformation_matrix>(bm.size());
 
             reduce_boundary_matrix(rbm, tm);
-            assert(rbm == bm * tm);
 
             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);
+        }
 
-#ifndef NDEBUG
-            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) {
-                // 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 IT r = lowestones.get_column(c).back();
-                    assert(lowestones.get_column(r).size() == 0);
-                    assert(lowestones.get_row(c).size() == 0);
-                }
-            }
-#endif
+        /** Get the lowest-one matrix of 'order'. */
+        const boundary_matrix& get_lowestones_matrix() const {
+            assert(done_matrices);
+            return lowestones;
         }
 
         /** Get the boundary matrix of 'order'. */

diff --git a/tests/persistence.h b/tests/persistence.h
index 1a6600c50139cbcf7641602990bfb9452eb2bee5..27b01f2f18c73aa0dbfc3ef79e85bc8d583f1419 100644 (file)
--- a/tests/persistence.h
+++ b/tests/persistence.h
@@ -30,6 +30,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);
         }
 
     protected:
@@ -194,7 +195,6 @@ class persistence_TestSuite: public Test::Suite {
             torusrbme.set_all(torusrbme_coords, sizeof(torusrbme_coords)/(2*sizeof(uint32_t)), true);
             TEST_ASSERT(torusrbme == torusrbm);
 
-
             pers ballp(oball);
             ballp.compute_matrices();
             const bm &ballbm = ballp.get_boundary_matrix();
@@ -283,6 +283,36 @@ class persistence_TestSuite: public Test::Suite {
             TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-3, oball.size()-3) == 2);
             TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-3, oball.size()-3) == 0);
         }
+
+        void test_lowestones() {
+            test_lowestones_impl(oring);
+            test_lowestones_impl(otorus);
+            test_lowestones_impl(oball);
+        }
+
+        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) {
+                // 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);
+                }
+            }
+        }
 };
 
 #endif