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::transformation_matrix tm;
  19
  20         bool setupcalled;
  21         scomplex ball, ring, torus;
  22         scomplex::simplex_order oball, oring, otorus;
  23
  24     public:
  25         persistence_TestSuite() :
  26             setupcalled(false),
  27             oball(ball),
  28             oring(ring),
  29             otorus(torus)
  30             {
  31             TEST_ADD(persistence_TestSuite::test_matrix_reduction);
  32             TEST_ADD(persistence_TestSuite::test_betti_numbers);
  33             TEST_ADD(persistence_TestSuite::test_lowestones);
  34         }
  35
  36     protected:
  37         virtual void setup() {
  38             if (setupcalled)
  39                 return;
  40             setupcalled = true;
  41
  42             //   7 ------------------------------------- 8
  43             //    \                                     /
  44             //     \    4 ------------------------- 5  /
  45             //      \    \                         /  /
  46             //       \    \      1 ------- 2      /  /
  47             //        \    \      \       /      /  /
  48             //         \    \      \     /      /  /
  49             //          \    \      \   /      /  /
  50             //           \    \      \ /      /  /
  51             //            \    \      3      /  /
  52             //             \    \           /  /
  53             //              \    \         /  /
  54             //               \    \       /  /
  55             //                \    \     /  /
  56             //                 \    \   /  /
  57             //                  \    \ /  /
  58             //                   \    6  /
  59             //                    \     /
  60             //                     \   /
  61             //                      \ /
  62             //                       9
  63
  64             scomplex::simplex ssring[] = {
  65                 // dimension, faces, value...
  66                 /*  1 */ {0, {},           1},
  67                 /*  2 */ {0, {},           2},
  68                 /*  3 */ {0, {},           3},
  69                 /*  4 */ {0, {},           4},
  70                 /*  5 */ {0, {},           5},
  71                 /*  6 */ {0, {},           6},
  72                 /*  7 */ {1, {2, 3},       6.01},
  73                 /*  8 */ {1, {3, 1},       6.02},
  74                 /*  9 */ {1, {1, 2},       6.03},
  75                 /* 10 */ {1, {4, 5},       6.04},
  76                 /* 11 */ {1, {5, 6},       6.05},
  77                 /* 12 */ {1, {6, 4},       6.06},
  78                 /* 13 */ {1, {1, 4},       6.07},
  79                 /* 14 */ {1, {1, 5},       6.08},
  80                 /* 15 */ {2, {13, 14, 10}, 6.0801},
  81                 /* 16 */ {1, {3, 6},       6.09},
  82                 /* 17 */ {1, {2, 5},       6.10},
  83                 /* 18 */ {1, {2, 6},       6.11},
  84                 /* 19 */ {2, {9, 14, 17},  6.1101},
  85                 /* 20 */ {2, {7, 16, 18},  6.1102},
  86                 /* 21 */ {2, {11, 17, 18}, 6.1103},
  87                 /* 22 */ {1, {1, 6},       6.12},
  88                 /* 23 */ {2, {12, 13, 22}, 6.1201},
  89                 /* 24 */ {2, {8, 16, 22},  6.1202}
  90             };
  91             const size_t cntssring = sizeof(ssring)/sizeof(scomplex::simplex);
  92             ring.add_simplices(ssring, cntssring);
  93             oring.reset();
  94
  95             scomplex::simplex sstorus[] = {
  96                 // dimension, faces, value...
  97                 /* 25 */ {0, {},           7},
  98                 /* 26 */ {0, {},           8},
  99                 /* 27 */ {0, {},           9},
 100                 /* 28 */ {1, {25, 26},     9.01},
 101                 /* 29 */ {1, {26, 27},     9.02},
 102                 /* 30 */ {1, {25, 27},     9.03},
 103                 /* 31 */ {1, {25, 1},      9.04},
 104                 /* 32 */ {1, {26, 1},      9.05},
 105                 /* 33 */ {2, {28, 31, 32}, 9.0501},
 106                 /* 34 */ {1, {27, 1},      9.06},
 107                 /* 35 */ {2, {30, 31, 34}, 9.0601},
 108                 /* 36 */ {1, {27, 3},      9.07},
 109                 /* 37 */ {2, {36, 34, 8},  9.0701},
 110                 /* 38 */ {1, {27, 2},      9.08},
 111                 /* 39 */ {1, {26, 2},      9.09},
 112                 /* 40 */ {2, {9, 32, 39},  9.0901},
 113                 /* 41 */ {2, {29, 38, 39}, 9.0902},
 114                 /* 42 */ {2, {7, 38, 36},  9.0903},
 115                 /* 43 */ {1, {4, 25},      9.10},
 116                 /* 44 */ {1, {5, 26},      9.11},
 117                 /* 45 */ {1, {6, 27},      9.12},
 118                 /* 46 */ {1, {4, 26},      9.13},
 119                 /* 47 */ {1, {4, 27},      9.14},
 120                 /* 48 */ {1, {5, 27},      9.15},
 121                 /* 49 */ {2, {43, 47, 30}, 9.1501},
 122                 /* 50 */ {2, {12, 45, 47}, 9.1502},
 123                 /* 51 */ {2, {29, 44, 48}, 9.1503},
 124                 /* 52 */ {2, {48, 11, 45}, 9.1504},
 125                 /* 53 */ {2, {10, 46, 44}, 9.1505},
 126                 /* 54 */ {2, {43, 46, 28}, 9.1506},
 127             };
 128             const size_t cntsstorus = sizeof(sstorus)/sizeof(scomplex::simplex);
 129             torus = ring;
 130             torus.add_simplices(sstorus, cntsstorus);
 131             otorus.reset();
 132
 133             scomplex::simplex ssball[] = {
 134                 // dimension, faces, value...
 135                 {0, {},                1},
 136                 {0, {},                2},
 137                 {0, {},                3},
 138                 {0, {},                4},
 139                 {0, {},                5},
 140                 {1, {1, 2},            6},
 141                 {1, {2, 3},            7},
 142                 {1, {3, 4},            8},
 143                 {1, {4, 1},            9},
 144                 {1, {1, 5},           10},
 145                 {1, {2, 5},           11},
 146                 {1, {3, 5},           12},
 147                 {1, {4, 5},           13},
 148                 {2, {6, 10, 11},      14},
 149                 {2, {7, 11, 12},      15},
 150                 {2, {8, 12, 13},      16},
 151                 {2, {9, 13, 10},      17},
 152                 {1, {1, 3},           18},
 153                 {2, {18, 6, 7},       19},
 154                 {2, {18, 8, 9},       20},
 155                 {2, {18, 10, 12},     21},
 156                 {3, {21, 14, 15, 19}, 22},
 157                 {3, {21, 16, 17, 20}, 23},
 158             };
 159             const size_t cntssball = sizeof(ssball)/sizeof(scomplex::simplex);
 160             ball.add_simplices(ssball, cntssball);
 161             oball.reset();
 162         }
 163
 164         virtual void tear_down() {
 165         }
 166
 167         void test_matrix_reduction() {
 168             pers ringp(oring);
 169             ringp.compute_matrices();
 170             const bm &ringbm = ringp.get_boundary_matrix();
 171             const bm &ringrbm = ringp.get_reduced_boundary_matrix();
 172             const tm &ringtm = ringp.get_transformation_matrix();
 173             TEST_ASSERT(ringrbm == ringbm * ringtm);
 174
 175             pers torusp(otorus);
 176             torusp.compute_matrices();
 177             const bm &torusbm = torusp.get_boundary_matrix();
 178             const bm &torusrbm = torusp.get_reduced_boundary_matrix();
 179             const tm &torustm = torusp.get_transformation_matrix();
 180             TEST_ASSERT(torusrbm == torusbm * torustm);
 181
 182             uint32_t torusrbme_coords[][2] = {
 183                     {0, 1}, {2, 7}, {3, 7}, {1, 8}, {2, 8}, {4, 10}, {5, 10}, {4, 11}, {6, 11}, {1, 13},
 184                     {4, 13}, {10, 15}, {13, 15}, {14, 15}, {9, 19}, {10, 19}, {13, 19}, {17, 19}, {7, 20},
 185                     {16, 20}, {18, 20}, {7, 21}, {9, 21}, {10, 21}, {11, 21}, {13, 21}, {16, 21}, {12, 23},
 186                     {13, 23}, {22, 23}, {7, 24}, {8, 24}, {9, 24}, {10, 24}, {11, 24}, {12, 24}, {25, 28},
 187                     {26, 28}, {25, 29}, {27, 29}, {1, 31}, {25, 31}, {28, 33}, {31, 33}, {32, 33}, {30, 35},
 188                     {31, 35}, {34, 35}, {8, 37}, {30, 37}, {31, 37}, {36, 37}, {9, 40}, {28, 40}, {31, 40},
 189                     {39, 40}, {9, 41}, {28, 41}, {29, 41}, {31, 41}, {38, 41}, {7, 42}, {8, 42}, {9, 42},
 190                     {28, 42}, {29, 42}, {30, 42}, {7, 49}, {8, 49}, {9, 49}, {28, 49}, {29, 49}, {43, 49},
 191                     {47, 49}, {10, 50}, {11, 50}, {28, 50}, {29, 50}, {43, 50}, {45, 50}, {29, 51}, {44, 51},
 192                     {48, 51}, {10, 52}, {28, 52}, {43, 52}, {44, 52}, {28, 53}, {43, 53}, {46, 53 }
 193             };
 194             bm torusrbme(torusbm.size());
 195             torusrbme.set_all(torusrbme_coords, sizeof(torusrbme_coords)/(2*sizeof(uint32_t)), true);
 196             TEST_ASSERT(torusrbme == torusrbm);
 197
 198             pers ballp(oball);
 199             ballp.compute_matrices();
 200             const bm &ballbm = ballp.get_boundary_matrix();
 201             const bm &ballrbm = ballp.get_reduced_boundary_matrix();
 202             const tm &balltm = ballp.get_transformation_matrix();
 203             TEST_ASSERT(ballrbm == ballbm * balltm);
 204
 205             uint32_t ballrbme_coords[][2] = {
 206                     {0, 1}, {1, 6}, {2, 6}, {1, 7}, {3, 7}, {1, 8}, {4, 8}, {1, 10}, {5, 10}, {6, 14},
 207                     {10, 14}, {11, 14}, {6, 15}, {7, 15}, {10, 15}, {12, 15}, {6, 16}, {7, 16}, {8, 16},
 208                     {10, 16}, {13, 16}, {6, 17}, {7, 17}, {8, 17}, {9, 17}, {6, 19}, {7, 19}, {18, 19},
 209                     {14, 22}, {15, 22}, {19, 22}, {21, 22}, {14, 23}, {15, 23}, {16, 23}, {17, 23}, {19, 23},
 210                     {20, 23 }
 211             };
 212             bm ballrbme(ballrbm.size());
 213             ballrbme.set_all(ballrbme_coords, sizeof(ballrbme_coords)/(2*sizeof(uint32_t)), true);
 214             TEST_ASSERT(ballrbme == ballrbm);
 215
 216             //std::cout << std::endl;
 217             //ballp.interpret_reduction(std::cout);
 218             //torusp.compute_diagrams();
 219             //torusp.interpret_persistence(std::cout);
 220             //std::cout << std::endl;
 221             //std::cout << std::endl;
 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 bm &lowestones = p.get_lowestones_matrix();
 297
 298             for (unsigned c=0; c < lowestones.size(); ++c)
 299                 assert(lowestones.get_column(c).size() <= 1);
 300             for (unsigned r=0; r < lowestones.size(); ++r)
 301                 assert(lowestones.get_row(r).size() <= 1);
 302             for (unsigned c=0; c < lowestones.size(); ++c) {
 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                 if (lowestones.get_column(c).size() > 0) {
 310                     const uint32_t r = lowestones.get_column(c).back();
 311                     assert(lowestones.get_column(r).size() == 0);
 312                     assert(lowestones.get_row(c).size() == 0);
 313                 }
 314             }
 315         }
 316 };
 317
 318 #endif