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