X-Git-Url: https://git.sthu.org/?a=blobdiff_plain;f=tests%2Fpersistence.h;h=aac47234bb7bd07f18b3042b7007319fcfb2dc22;hb=39f867a091cbcdbf67499403d08fbe28103f06b7;hp=1a6600c50139cbcf7641602990bfb9452eb2bee5;hpb=2dbba971939606ec3cc2bbd07978ff61e9566d01;p=libstick.git diff --git a/tests/persistence.h b/tests/persistence.h index 1a6600c..aac4723 100644 --- a/tests/persistence.h +++ b/tests/persistence.h @@ -4,7 +4,7 @@ #include #include -#include +#include 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: @@ -61,68 +66,68 @@ 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} + // 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); 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}, + // 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); torus = ring; @@ -130,30 +135,30 @@ class persistence_TestSuite: public Test::Suite { 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}, + // 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); ball.add_simplices(ssball, cntssball); @@ -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(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(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(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