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