tests/persistence.h

   1 #ifndef persistence_h_Fahlaewahgaebaqu
   2 #define persistence_h_Fahlaewahgaebaqu
   3
   4 #include <cpptest.h>
   5 #include <cpptest-suite.h>
   6
   7 #include <libstick/persistence.h>
   8
   9 using namespace libstick;
  10
  11
  12 class persistence_TestSuite: public Test::Suite {
  13
  14     private:
  15         typedef simplicial_complex<3, uint32_t, double> scomplex;
  16         typedef persistence<3, uint32_t, double> pers;
  17         typedef pers::boundary_matrix bm;
  18         typedef pers::lowestones_matrix lom;
  19         typedef pers::transformation_matrix tm;
  20
  21         bool setupcalled;
  22         scomplex ball, ring, torus;
  23         scomplex::simplex_order oball, oring, otorus;
  24
  25     public:
  26         persistence_TestSuite() :
  27             setupcalled(false),
  28             oball(ball),
  29             oring(ring),
  30             otorus(torus)
  31             {
  32             TEST_ADD(persistence_TestSuite::test_matrix_reduction);
  33             TEST_ADD(persistence_TestSuite::test_betti_numbers);
  34             TEST_ADD(persistence_TestSuite::test_lowestones);
  35             TEST_ADD(persistence_TestSuite::test_diagram);
  36         }
  37
  38     protected:
  39         virtual void setup() {
  40             if (setupcalled)
  41                 return;
  42             setupcalled = true;
  43
  44             //   7 ------------------------------------- 8
  45             //    \                                     /
  46             //     \    4 ------------------------- 5  /
  47             //      \    \                         /  /
  48             //       \    \      1 ------- 2      /  /
  49             //        \    \      \       /      /  /
  50             //         \    \      \     /      /  /
  51             //          \    \      \   /      /  /
  52             //           \    \      \ /      /  /
  53             //            \    \      3      /  /
  54             //             \    \           /  /
  55             //              \    \         /  /
  56             //               \    \       /  /
  57             //                \    \     /  /
  58             //                 \    \   /  /
  59             //                  \    \ /  /
  60             //                   \    6  /
  61             //                    \     /
  62             //                     \   /
  63             //                      \ /
  64             //                       9
  65
  66             scomplex::simplex ssring[] = {
  67                 // dimension, faces, value...
  68                 /*  1 */ {0, {},           1},
  69                 /*  2 */ {0, {},           2},
  70                 /*  3 */ {0, {},           3},
  71                 /*  4 */ {0, {},           4},
  72                 /*  5 */ {0, {},           5},
  73                 /*  6 */ {0, {},           6},
  74                 /*  7 */ {1, {2, 3},       6.01},
  75                 /*  8 */ {1, {3, 1},       6.02},
  76                 /*  9 */ {1, {1, 2},       6.03},
  77                 /* 10 */ {1, {4, 5},       6.04},
  78                 /* 11 */ {1, {5, 6},       6.05},
  79                 /* 12 */ {1, {6, 4},       6.06},
  80                 /* 13 */ {1, {1, 4},       6.07},
  81                 /* 14 */ {1, {1, 5},       6.08},
  82                 /* 15 */ {2, {13, 14, 10}, 6.0801},
  83                 /* 16 */ {1, {3, 6},       6.09},
  84                 /* 17 */ {1, {2, 5},       6.10},
  85                 /* 18 */ {1, {2, 6},       6.11},
  86                 /* 19 */ {2, {9, 14, 17},  6.1101},
  87                 /* 20 */ {2, {7, 16, 18},  6.1102},
  88                 /* 21 */ {2, {11, 17, 18}, 6.1103},
  89                 /* 22 */ {1, {3, 4},       6.12},
  90                 /* 23 */ {2, {12, 16, 22}, 6.1201},
  91                 /* 24 */ {2, {8, 13, 22},  6.1202}
  92             };
  93             const size_t cntssring = sizeof(ssring)/sizeof(scomplex::simplex);
  94             ring.add_simplices(ssring, cntssring);
  95             oring.reset();
  96
  97             scomplex::simplex sstorus[] = {
  98                 // dimension, faces, value...
  99                 /* 25 */ {0, {},           7},
 100                 /* 26 */ {0, {},           8},
 101                 /* 27 */ {0, {},           9},
 102                 /* 28 */ {1, {25, 26},     9.01},
 103                 /* 29 */ {1, {26, 27},     9.02},
 104                 /* 30 */ {1, {25, 27},     9.03},
 105                 /* 31 */ {1, {25, 1},      9.04},
 106                 /* 32 */ {1, {26, 1},      9.05},
 107                 /* 33 */ {2, {28, 31, 32}, 9.0501},
 108                 /* 34 */ {1, {27, 1},      9.06},
 109                 /* 35 */ {2, {30, 31, 34}, 9.0601},
 110                 /* 36 */ {1, {27, 3},      9.07},
 111                 /* 37 */ {2, {36, 34, 8},  9.0701},
 112                 /* 38 */ {1, {27, 2},      9.08},
 113                 /* 39 */ {1, {26, 2},      9.09},
 114                 /* 40 */ {2, {9, 32, 39},  9.0901},
 115                 /* 41 */ {2, {29, 38, 39}, 9.0902},
 116                 /* 42 */ {2, {7, 38, 36},  9.0903},
 117                 /* 43 */ {1, {4, 25},      9.10},
 118                 /* 44 */ {1, {5, 26},      9.11},
 119                 /* 45 */ {1, {6, 27},      9.12},
 120                 /* 46 */ {1, {4, 26},      9.13},
 121                 /* 47 */ {1, {4, 27},      9.14},
 122                 /* 48 */ {1, {5, 27},      9.15},
 123                 /* 49 */ {2, {43, 47, 30}, 9.1501},
 124                 /* 50 */ {2, {12, 45, 47}, 9.1502},
 125                 /* 51 */ {2, {29, 44, 48}, 9.1503},
 126                 /* 52 */ {2, {48, 11, 45}, 9.1504},
 127                 /* 53 */ {2, {10, 46, 44}, 9.1505},
 128                 /* 54 */ {2, {43, 46, 28}, 9.1506},
 129             };
 130             const size_t cntsstorus = sizeof(sstorus)/sizeof(scomplex::simplex);
 131             torus = ring;
 132             torus.add_simplices(sstorus, cntsstorus);
 133             otorus.reset();
 134
 135             scomplex::simplex ssball[] = {
 136                 // dimension, faces, value...
 137                 {0, {},                1},
 138                 {0, {},                2},
 139                 {0, {},                3},
 140                 {0, {},                4},
 141                 {0, {},                5},
 142                 {1, {1, 2},            6},
 143                 {1, {2, 3},            7},
 144                 {1, {3, 4},            8},
 145                 {1, {4, 1},            9},
 146                 {1, {1, 5},           10},
 147                 {1, {2, 5},           11},
 148                 {1, {3, 5},           12},
 149                 {1, {4, 5},           13},
 150                 {2, {6, 10, 11},      14},
 151                 {2, {7, 11, 12},      15},
 152                 {2, {8, 12, 13},      16},
 153                 {2, {9, 13, 10},      17},
 154                 {1, {1, 3},           18},
 155                 {2, {18, 6, 7},       19},
 156                 {2, {18, 8, 9},       20},
 157                 {2, {18, 10, 12},     21},
 158                 {3, {21, 14, 15, 19}, 22},
 159                 {3, {21, 16, 17, 20}, 23},
 160             };
 161             const size_t cntssball = sizeof(ssball)/sizeof(scomplex::simplex);
 162             ball.add_simplices(ssball, cntssball);
 163             oball.reset();
 164         }
 165
 166         virtual void tear_down() {
 167         }
 168
 169         void test_matrix_reduction() {
 170             pers ringp(oring);
 171             ringp.compute_matrices();
 172             const bm &ringbm = ringp.get_boundary_matrix();
 173             const bm &ringrbm = ringp.get_reduced_boundary_matrix();
 174             const tm &ringtm = ringp.get_transformation_matrix();
 175             TEST_ASSERT(ringrbm == boolean_colrowmatrix<uint32_t>(ringbm) * ringtm);
 176
 177
 178             pers torusp(otorus);
 179             torusp.compute_matrices();
 180             const bm &torusbm = torusp.get_boundary_matrix();
 181             const bm &torusrbm = torusp.get_reduced_boundary_matrix();
 182             const tm &torustm = torusp.get_transformation_matrix();
 183             TEST_ASSERT(torusrbm == boolean_colrowmatrix<uint32_t>(torusbm) * torustm);
 184
 185             const lom &toruslo = torusp.get_lowestones_matrix();
 186             uint32_t torusloe_coords[] = {1, 0, 8, 7, 13, 10, 11, 0, 0, 0, 0,
 187                 0, 24, 0, 15, 0, 21, 19, 20, 0, 0, 0, 23, 0, 0, 31, 28, 29, 0,
 188                 0, 42, 0, 33, 0, 35, 0, 37, 0, 41, 40, 0, 0, 0, 0, 52, 50, 53,
 189                 49, 51, 0, 0, 0, 0, 0, 0};
 190             lom torusloe(torusloe_coords, torusloe_coords + sizeof(torusloe_coords)/sizeof(uint32_t));
 191             TEST_ASSERT(torusloe == toruslo);
 192
 193
 194             pers ballp(oball);
 195             ballp.compute_matrices();
 196             const bm &ballbm = ballp.get_boundary_matrix();
 197             const bm &ballrbm = ballp.get_reduced_boundary_matrix();
 198             const tm &balltm = ballp.get_transformation_matrix();
 199             TEST_ASSERT(ballrbm == boolean_colrowmatrix<uint32_t>(ballbm) * balltm);
 200
 201             const lom &balllo = ballp.get_lowestones_matrix();
 202             uint32_t ballloe_coords[] = {1, 0, 6, 7, 8, 10, 0, 0, 0, 17, 0, 14,
 203                 15, 16, 0, 0, 0, 0, 19, 0, 23, 22, 0, 0};
 204             lom ballloe(ballloe_coords, ballloe_coords + sizeof(ballloe_coords)/sizeof(uint32_t));
 205             TEST_ASSERT(ballloe == balllo);
 206
 207             //std::cout << std::endl;
 208             //ballp.interpret_reduction(std::cout);
 209             //torusp.compute_diagrams();
 210             //torusp.interpret_reduction(std::cout);
 211             //torusp.interpret_persistence(std::cout);
 212             //std::cout << std::endl;
 213             //std::cout << std::endl;
 214         }
 215
 216         void test_betti_numbers() {
 217             pers ringp(oring);
 218             ringp.compute_diagrams();
 219
 220             TEST_ASSERT(ringp.get_persistence_diagram(-1).betti() == 0);
 221             TEST_ASSERT(ringp.get_persistence_diagram(0).betti() == 0);
 222             TEST_ASSERT(ringp.get_persistence_diagram(1).betti() == 1);
 223             TEST_ASSERT(ringp.get_persistence_diagram(2).betti() == 0);
 224             TEST_ASSERT(ringp.get_persistence_diagram(3).betti() == 0);
 225             for (unsigned d=0; d < 3; ++d)
 226                 TEST_ASSERT(ringp.get_persistence_diagram(d).persistent_betti(oring.size()-1, oring.size()-1) ==
 227                         ringp.get_persistence_diagram(d).betti());
 228
 229
 230             pers torusp(otorus);
 231             torusp.compute_diagrams();
 232
 233             TEST_ASSERT(torusp.get_persistence_diagram(-1).betti() == 0);
 234             TEST_ASSERT(torusp.get_persistence_diagram(0).betti() == 0);
 235             TEST_ASSERT(torusp.get_persistence_diagram(1).betti() == 2);
 236             TEST_ASSERT(torusp.get_persistence_diagram(2).betti() == 1);
 237             TEST_ASSERT(torusp.get_persistence_diagram(3).betti() == 0);
 238             for (unsigned d=0; d < 3; ++d)
 239                 TEST_ASSERT(torusp.get_persistence_diagram(d).persistent_betti(otorus.size()-1, otorus.size()-1) ==
 240                         torusp.get_persistence_diagram(d).betti());
 241
 242             // Torus was made from ring. Hence, persistent betti numbers must correspond to betti numbers of ring.
 243             for (unsigned d=0; d < 3; ++d)
 244                 TEST_ASSERT(torusp.get_persistence_diagram(d).persistent_betti(oring.size()-1, oring.size()-1) ==
 245                         ringp.get_persistence_diagram(d).betti());
 246
 247
 248             pers ballp(oball);
 249             ballp.compute_diagrams();
 250
 251             //std::cout << std::endl;
 252             //ballp.interpret_reduction(std::cout);
 253             //std::cout << std::endl;
 254
 255             TEST_ASSERT(ballp.get_persistence_diagram(-1).betti() == 0);
 256             TEST_ASSERT(ballp.get_persistence_diagram(0).betti() == 0);
 257             TEST_ASSERT(ballp.get_persistence_diagram(1).betti() == 0);
 258             TEST_ASSERT(ballp.get_persistence_diagram(2).betti() == 0);
 259             TEST_ASSERT(ballp.get_persistence_diagram(3).betti() == 0);
 260             for (unsigned d=0; d < 3; ++d)
 261                 TEST_ASSERT(ballp.get_persistence_diagram(d).persistent_betti(oball.size()-1, oball.size()-1) ==
 262                         ballp.get_persistence_diagram(d).betti());
 263
 264             // Before we inserted the two tetrahedra and the interior triangle, we had a sphere.
 265             TEST_ASSERT(ballp.get_persistence_diagram(-1).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 266             TEST_ASSERT(ballp.get_persistence_diagram(0).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 267             TEST_ASSERT(ballp.get_persistence_diagram(1).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 268             TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-4, oball.size()-4) == 1);
 269             TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 270
 271             // Before we inserted the two tetrahedra, we had a double-sphere.
 272             TEST_ASSERT(ballp.get_persistence_diagram(-1).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 273             TEST_ASSERT(ballp.get_persistence_diagram(0).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 274             TEST_ASSERT(ballp.get_persistence_diagram(1).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 275             TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-3, oball.size()-3) == 2);
 276             TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 277         }
 278
 279         void test_lowestones() {
 280             test_lowestones_impl(oring);
 281             test_lowestones_impl(otorus);
 282             test_lowestones_impl(oball);
 283         }
 284
 285         void test_lowestones_impl(const scomplex::simplex_order& so) {
 286             pers p(so);
 287             p.compute_matrices();
 288             const lom &lowestones = p.get_lowestones_matrix();
 289
 290             // Check that lowestones contains every column-index only once,
 291             // except the zeros.
 292             lom cpy = lowestones;
 293             sort(cpy.begin(), cpy.end());
 294             lom::iterator it = cpy.begin();
 295             // Skip the zeros
 296             while (it != cpy.end() && *it == 0)
 297                 ++it;
 298             // From now on, no column index appears twice
 299             for (; it+1 != cpy.end(); ++it)
 300                 TEST_ASSERT(*it < *(it+1));
 301
 302             for (unsigned r=0; r < lowestones.size(); ++r) {
 303                 // If (r,c) is a one then
 304                 //   (i) a cycle dies with c --> row c is zero
 305                 //  (ii) a cycle is born with r --> column r is zero
 306                 //Hence
 307                 //   (i) the column r is a zero-column
 308                 //   (i) the row c is a zero-column
 309
 310                 const uint32_t c = lowestones[r];
 311                 if (c > 0) {
 312                     TEST_ASSERT(p.get_reduced_boundary_matrix().get(r,c) == true);
 313                     TEST_ASSERT(p.get_reduced_boundary_matrix().get_column(r).size() == 0);
 314                     TEST_ASSERT(lowestones[c] == 0);
 315                 }
 316             }
 317         }
 318
 319         void test_diagram() {
 320             test_diagram_impl(oring);
 321             test_diagram_impl(otorus);
 322             test_diagram_impl(oball);
 323         }
 324
 325         void test_diagram_impl(const scomplex::simplex_order& so) {
 326             pers p(so);
 327             p.compute_diagrams();
 328
 329             // Even Jesus died after his birth. Homology classes are nothing better.
 330             for (int d=-1; d <= 3; ++d) {
 331                 const pers::diagram &dia = p.get_persistence_diagram(d);
 332                 for (unsigned i=0; i < dia.births.size(); ++i)
 333                     if (dia.births[i].death > 0)
 334                         assert(dia.births[i].birth < dia.births[i].death);
 335             }
 336
 337         }
 338 };
 339
 340 #endif