Use scomplex throughout rather than simplcompltype
[libstick.git] / tests / persistence.h
index f48c3867c32efdf55603b960cb5c31e07ab94235..aac47234bb7bd07f18b3042b7007319fcfb2dc22 100644 (file)
@@ -4,7 +4,7 @@
 #include <cpptest.h>
 #include <cpptest-suite.h>
 
-#include <libstick-0.1/persistence.h>
+#include <libstick/persistence.h>
 
 using namespace libstick;
 
@@ -12,9 +12,11 @@ using namespace libstick;
 class persistence_TestSuite: public Test::Suite {
 
     private:
-        typedef simplicial_complex<3, uint32_t, double> scomplex;
-        typedef persistence<3, uint32_t, double> pers;
+        typedef simplicial_function<3, uint32_t, uint32_t> sfunction;
+        typedef sfunction::scomplex scomplex;
+        typedef persistence<3, uint32_t> pers;
         typedef pers::boundary_matrix bm;
+        typedef pers::lowestones_matrix lom;
         typedef pers::transformation_matrix tm;
 
         bool setupcalled;
@@ -30,6 +32,9 @@ 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_diagram);
+            TEST_ADD(persistence_TestSuite::test_diagrampoint_value);
         }
 
     protected:
@@ -60,102 +65,102 @@ class persistence_TestSuite: public Test::Suite {
             //                      \ /
             //                       9
 
-            scomplex::Simplex ssring[] = {
-                // dimension, faces, value...
-                /*  1 */ {0, {},           1},
-                /*  2 */ {0, {},           2},
-                /*  3 */ {0, {},           3},
-                /*  4 */ {0, {},           4},
-                /*  5 */ {0, {},           5},
-                /*  6 */ {0, {},           6},
-                /*  7 */ {1, {2, 3},       6.01},
-                /*  8 */ {1, {3, 1},       6.02},
-                /*  9 */ {1, {1, 2},       6.03},
-                /* 10 */ {1, {4, 5},       6.04},
-                /* 11 */ {1, {5, 6},       6.05},
-                /* 12 */ {1, {6, 4},       6.06},
-                /* 13 */ {1, {1, 4},       6.07},
-                /* 14 */ {1, {1, 5},       6.08},
-                /* 15 */ {2, {13, 14, 10}, 6.0801},
-                /* 16 */ {1, {3, 6},       6.09},
-                /* 17 */ {1, {2, 5},       6.10},
-                /* 18 */ {1, {2, 6},       6.11},
-                /* 19 */ {2, {9, 14, 17},  6.1101},
-                /* 20 */ {2, {7, 16, 18},  6.1102},
-                /* 21 */ {2, {11, 17, 18}, 6.1103},
-                /* 22 */ {1, {1, 6},       6.12},
-                /* 23 */ {2, {12, 13, 22}, 6.1201},
-                /* 24 */ {2, {8, 16, 22},  6.1202}
+            scomplex::simplex ssring[] = {
+                // dimension, faces
+                /*  1 */ {0, {}},
+                /*  2 */ {0, {}},
+                /*  3 */ {0, {}},
+                /*  4 */ {0, {}},
+                /*  5 */ {0, {}},
+                /*  6 */ {0, {}},
+                /*  7 */ {1, {2, 3}},
+                /*  8 */ {1, {3, 1}},
+                /*  9 */ {1, {1, 2}},
+                /* 10 */ {1, {4, 5}},
+                /* 11 */ {1, {5, 6}},
+                /* 12 */ {1, {6, 4}},
+                /* 13 */ {1, {1, 4}},
+                /* 14 */ {1, {1, 5}},
+                /* 15 */ {2, {13, 14, 10}},
+                /* 16 */ {1, {3, 6}},
+                /* 17 */ {1, {2, 5}},
+                /* 18 */ {1, {2, 6}},
+                /* 19 */ {2, {9, 14, 17}},
+                /* 20 */ {2, {7, 16, 18}},
+                /* 21 */ {2, {11, 17, 18}},
+                /* 22 */ {1, {3, 4}},
+                /* 23 */ {2, {12, 16, 22}},
+                /* 24 */ {2, {8, 13, 22}},
             };
-            const size_t cntssring = sizeof(ssring)/sizeof(scomplex::Simplex);
+            const size_t cntssring = sizeof(ssring)/sizeof(scomplex::simplex);
             ring.add_simplices(ssring, cntssring);
             oring.reset();
 
-            scomplex::Simplex sstorus[] = {
-                // dimension, faces, value...
-                /* 25 */ {0, {},           7},
-                /* 26 */ {0, {},           8},
-                /* 27 */ {0, {},           9},
-                /* 28 */ {1, {25, 26},     9.01},
-                /* 29 */ {1, {26, 27},     9.02},
-                /* 30 */ {1, {25, 27},     9.03},
-                /* 31 */ {1, {25, 1},      9.04},
-                /* 32 */ {1, {26, 1},      9.05},
-                /* 33 */ {2, {28, 31, 32}, 9.0501},
-                /* 34 */ {1, {27, 1},      9.06},
-                /* 35 */ {2, {30, 31, 34}, 9.0601},
-                /* 36 */ {1, {27, 3},      9.07},
-                /* 37 */ {2, {36, 34, 8},  9.0701},
-                /* 38 */ {1, {27, 2},      9.08},
-                /* 39 */ {1, {26, 2},      9.09},
-                /* 40 */ {2, {9, 32, 39},  9.0901},
-                /* 41 */ {2, {29, 38, 39}, 9.0902},
-                /* 42 */ {2, {7, 38, 36},  9.0903},
-                /* 43 */ {1, {4, 25},      9.10},
-                /* 44 */ {1, {5, 26},      9.11},
-                /* 45 */ {1, {6, 27},      9.12},
-                /* 46 */ {1, {4, 26},      9.13},
-                /* 47 */ {1, {4, 27},      9.14},
-                /* 48 */ {1, {5, 27},      9.15},
-                /* 49 */ {2, {43, 47, 30}, 9.1501},
-                /* 50 */ {2, {12, 45, 47}, 9.1502},
-                /* 51 */ {2, {29, 44, 48}, 9.1503},
-                /* 52 */ {2, {48, 11, 45}, 9.1504},
-                /* 53 */ {2, {10, 46, 44}, 9.1505},
-                /* 54 */ {2, {43, 46, 28}, 9.1506},
+            scomplex::simplex sstorus[] = {
+                // dimension, faces
+                /* 25 */ {0, {}},
+                /* 26 */ {0, {}},
+                /* 27 */ {0, {}},
+                /* 28 */ {1, {25, 26}},
+                /* 29 */ {1, {26, 27}},
+                /* 30 */ {1, {25, 27}},
+                /* 31 */ {1, {25, 1}},
+                /* 32 */ {1, {26, 1}},
+                /* 33 */ {2, {28, 31, 32}},
+                /* 34 */ {1, {27, 1}},
+                /* 35 */ {2, {30, 31, 34}},
+                /* 36 */ {1, {27, 3}},
+                /* 37 */ {2, {36, 34, 8}},
+                /* 38 */ {1, {27, 2}},
+                /* 39 */ {1, {26, 2}},
+                /* 40 */ {2, {9, 32, 39}},
+                /* 41 */ {2, {29, 38, 39}},
+                /* 42 */ {2, {7, 38, 36}},
+                /* 43 */ {1, {4, 25}},
+                /* 44 */ {1, {5, 26}},
+                /* 45 */ {1, {6, 27}},
+                /* 46 */ {1, {4, 26}},
+                /* 47 */ {1, {4, 27}},
+                /* 48 */ {1, {5, 27}},
+                /* 49 */ {2, {43, 47, 30}},
+                /* 50 */ {2, {12, 45, 47}},
+                /* 51 */ {2, {29, 44, 48}},
+                /* 52 */ {2, {48, 11, 45}},
+                /* 53 */ {2, {10, 46, 44}},
+                /* 54 */ {2, {43, 46, 28}},
             };
-            const size_t cntsstorus = sizeof(sstorus)/sizeof(scomplex::Simplex);
+            const size_t cntsstorus = sizeof(sstorus)/sizeof(scomplex::simplex);
             torus = ring;
             torus.add_simplices(sstorus, cntsstorus);
             otorus.reset();
 
-            scomplex::Simplex ssball[] = {
-                // dimension, faces, value...
-                {0, {},                1},
-                {0, {},                2},
-                {0, {},                3},
-                {0, {},                4},
-                {0, {},                5},
-                {1, {1, 2},            6},
-                {1, {2, 3},            7},
-                {1, {3, 4},            8},
-                {1, {4, 1},            9},
-                {1, {1, 5},           10},
-                {1, {2, 5},           11},
-                {1, {3, 5},           12},
-                {1, {4, 5},           13},
-                {2, {6, 10, 11},      14},
-                {2, {7, 11, 12},      15},
-                {2, {8, 12, 13},      16},
-                {2, {9, 13, 10},      17},
-                {1, {1, 3},           18},
-                {2, {18, 6, 7},       19},
-                {2, {18, 8, 9},       20},
-                {2, {18, 10, 12},     21},
-                {3, {21, 14, 15, 19}, 22},
-                {3, {21, 16, 17, 20}, 23},
+            scomplex::simplex ssball[] = {
+                // dimension, faces
+                {0, {}},
+                {0, {}},
+                {0, {}},
+                {0, {}},
+                {0, {}},
+                {1, {1, 2}},
+                {1, {2, 3}},
+                {1, {3, 4}},
+                {1, {4, 1}},
+                {1, {1, 5}},
+                {1, {2, 5}},
+                {1, {3, 5}},
+                {1, {4, 5}},
+                {2, {6, 10, 11}},
+                {2, {7, 11, 12}},
+                {2, {8, 12, 13}},
+                {2, {9, 13, 10}},
+                {1, {1, 3}},
+                {2, {18, 6, 7}},
+                {2, {18, 8, 9}},
+                {2, {18, 10, 12}},
+                {3, {21, 14, 15, 19}},
+                {3, {21, 16, 17, 20}},
             };
-            const size_t cntssball = sizeof(ssball)/sizeof(scomplex::Simplex);
+            const size_t cntssball = sizeof(ssball)/sizeof(scomplex::simplex);
             ball.add_simplices(ssball, cntssball);
             oball.reset();
         }
@@ -169,30 +174,23 @@ 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();
-            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();
-            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);
@@ -200,25 +198,50 @@ class persistence_TestSuite: public Test::Suite {
             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);
+
+#if 0
+            sfunction ftorus(torus);
+            ftorus.make_complete();
+            for (unsigned i=1; i < torus.size(); ++i)
+                ftorus.set_value(i, i);
+
+            std::cout << std::endl;
+            ballp.interpret_reduction(std::cout);
+            torusp.compute_diagrams();
+            torusp.interpret_reduction(std::cout, &ftorus);
+            torusp.interpret_persistence(std::cout, &ftorus);
+            std::cout << std::endl;
+            std::cout << std::endl;
+#endif
+        }
 
-            //std::cout << std::endl;
-            //ballp.interpret_reduction(std::cout);
-            //torusp.compute_diagrams();
-            //torusp.interpret_persistence(std::cout);
-            //std::cout << std::endl;
-            //std::cout << std::endl;
+        void test_diagrampoint_value() {
+            pers torusp(otorus);
+            torusp.compute_diagrams();
+
+            // We make a simplicial function that assigns to each simplex twice
+            // the index as value.
+            sfunction ftorus(torus);
+            ftorus.make_complete();
+            for (unsigned i=1; i < torus.size(); ++i)
+                ftorus.set_value(i, 2*i);
+
+            for (unsigned d=0; d <= 3; ++d) {
+                const pers::diagram &dia = torusp.get_persistence_diagram(d);
+                for (unsigned i=0; i < dia.births.size(); ++i) {
+                    TEST_ASSERT(2 * dia.births[i].birth == dia.births[i].get_birth_value(otorus, ftorus));
+                    TEST_ASSERT(2 * dia.births[i].death == dia.births[i].get_death_value(otorus, ftorus));
+                    TEST_ASSERT(2 * (dia.births[i].death - dia.births[i].birth) ==
+                            dia.births[i].get_persistence(otorus, ftorus));
+                }
+            }
         }
 
         void test_betti_numbers() {
@@ -283,6 +306,66 @@ 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 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
+
+                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