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-0.1/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, {1, 6},       6.12},
  90                 /* 23 */ {2, {12, 13, 22}, 6.1201},
  91                 /* 24 */ {2, {8, 16, 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_persistence(std::cout);
 211             //std::cout << std::endl;
 212             //std::cout << std::endl;
 213         }
 214
 215         void test_betti_numbers() {
 216             pers ringp(oring);
 217             ringp.compute_diagrams();
 218
 219             TEST_ASSERT(ringp.get_persistence_diagram(-1).betti() == 0);
 220             TEST_ASSERT(ringp.get_persistence_diagram(0).betti() == 0);
 221             TEST_ASSERT(ringp.get_persistence_diagram(1).betti() == 1);
 222             TEST_ASSERT(ringp.get_persistence_diagram(2).betti() == 0);
 223             TEST_ASSERT(ringp.get_persistence_diagram(3).betti() == 0);
 224             for (unsigned d=0; d < 3; ++d)
 225                 TEST_ASSERT(ringp.get_persistence_diagram(d).persistent_betti(oring.size()-1, oring.size()-1) ==
 226                         ringp.get_persistence_diagram(d).betti());
 227
 228
 229             pers torusp(otorus);
 230             torusp.compute_diagrams();
 231
 232             TEST_ASSERT(torusp.get_persistence_diagram(-1).betti() == 0);
 233             TEST_ASSERT(torusp.get_persistence_diagram(0).betti() == 0);
 234             TEST_ASSERT(torusp.get_persistence_diagram(1).betti() == 2);
 235             TEST_ASSERT(torusp.get_persistence_diagram(2).betti() == 1);
 236             TEST_ASSERT(torusp.get_persistence_diagram(3).betti() == 0);
 237             for (unsigned d=0; d < 3; ++d)
 238                 TEST_ASSERT(torusp.get_persistence_diagram(d).persistent_betti(otorus.size()-1, otorus.size()-1) ==
 239                         torusp.get_persistence_diagram(d).betti());
 240
 241             // Torus was made from ring. Hence, persistent betti numbers must correspond to betti numbers of ring.
 242             for (unsigned d=0; d < 3; ++d)
 243                 TEST_ASSERT(torusp.get_persistence_diagram(d).persistent_betti(oring.size()-1, oring.size()-1) ==
 244                         ringp.get_persistence_diagram(d).betti());
 245
 246
 247             pers ballp(oball);
 248             ballp.compute_diagrams();
 249
 250             //std::cout << std::endl;
 251             //ballp.interpret_reduction(std::cout);
 252             //std::cout << std::endl;
 253
 254             TEST_ASSERT(ballp.get_persistence_diagram(-1).betti() == 0);
 255             TEST_ASSERT(ballp.get_persistence_diagram(0).betti() == 0);
 256             TEST_ASSERT(ballp.get_persistence_diagram(1).betti() == 0);
 257             TEST_ASSERT(ballp.get_persistence_diagram(2).betti() == 0);
 258             TEST_ASSERT(ballp.get_persistence_diagram(3).betti() == 0);
 259             for (unsigned d=0; d < 3; ++d)
 260                 TEST_ASSERT(ballp.get_persistence_diagram(d).persistent_betti(oball.size()-1, oball.size()-1) ==
 261                         ballp.get_persistence_diagram(d).betti());
 262
 263             // Before we inserted the two tetrahedra and the interior triangle, we had a sphere.
 264             TEST_ASSERT(ballp.get_persistence_diagram(-1).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 265             TEST_ASSERT(ballp.get_persistence_diagram(0).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 266             TEST_ASSERT(ballp.get_persistence_diagram(1).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 267             TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-4, oball.size()-4) == 1);
 268             TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-4, oball.size()-4) == 0);
 269
 270             // Before we inserted the two tetrahedra, we had a double-sphere.
 271             TEST_ASSERT(ballp.get_persistence_diagram(-1).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 272             TEST_ASSERT(ballp.get_persistence_diagram(0).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 273             TEST_ASSERT(ballp.get_persistence_diagram(1).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 274             TEST_ASSERT(ballp.get_persistence_diagram(2).persistent_betti(oball.size()-3, oball.size()-3) == 2);
 275             TEST_ASSERT(ballp.get_persistence_diagram(3).persistent_betti(oball.size()-3, oball.size()-3) == 0);
 276         }
 277
 278         void test_lowestones() {
 279             test_lowestones_impl(oring);
 280             test_lowestones_impl(otorus);
 281             test_lowestones_impl(oball);
 282         }
 283
 284         void test_lowestones_impl(const scomplex::simplex_order& so) {
 285             pers p(so);
 286             p.compute_matrices();
 287             const lom &lowestones = p.get_lowestones_matrix();
 288
 289             // Check that lowestones contains every column-index only once,
 290             // except the zeros.
 291             lom cpy = lowestones;
 292             sort(cpy.begin(), cpy.end());
 293             lom::iterator it = cpy.begin();
 294             // Skip the zeros
 295             while (it != cpy.end() && *it == 0)
 296                 ++it;
 297             // From now on, no column index appears twice
 298             for (; it+1 != cpy.end(); ++it)
 299                 TEST_ASSERT(*it < *(it+1));
 300
 301             for (unsigned r=0; r < lowestones.size(); ++r) {
 302                 // If (r,c) is a one then
 303                 //   (i) a cycle dies with c --> row c is zero
 304                 //  (ii) a cycle is born with r --> column r is zero
 305                 //Hence
 306                 //   (i) the column r is a zero-column
 307                 //   (i) the row c is a zero-column
 308
 309                 const uint32_t c = lowestones[r];
 310                 if (c > 0) {
 311                     TEST_ASSERT(p.get_reduced_boundary_matrix().get(r,c) == true);
 312                     TEST_ASSERT(p.get_reduced_boundary_matrix().get_column(r).size() == 0);
 313                     TEST_ASSERT(lowestones[c] == 0);
 314                 }
 315             }
 316         }
 317
 318         void test_diagram() {
 319             test_diagram_impl(oring);
 320             test_diagram_impl(otorus);
 321             test_diagram_impl(oball);
 322         }
 323
 324         void test_diagram_impl(const scomplex::simplex_order& so) {
 325             pers p(so);
 326             p.compute_diagrams();
 327
 328             // Even Jesus died after his birth. Homology classes are nothing better.
 329             for (int d=-1; d <= 3; ++d) {
 330                 const pers::diagram &dia = p.get_persistence_diagram(d);
 331                 for (unsigned i=0; i < dia.births.size(); ++i)
 332                     if (dia.births[i].death > 0)
 333                         assert(dia.births[i].birth < dia.births[i].death);
 334             }
 335
 336         }
 337 };
 338
 339 #endif