1 #ifndef persistence_h_Fahlaewahgaebaqu
2 #define persistence_h_Fahlaewahgaebaqu
5 #include <cpptest-suite.h>
7 #include <libstick-0.1/persistence.h>
9 using namespace libstick
;
12 class persistence_TestSuite
: public Test::Suite
{
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
;
22 scomplex ball
, ring
, torus
;
23 scomplex::simplex_order oball
, oring
, otorus
;
26 persistence_TestSuite() :
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
);
39 virtual void setup() {
44 // 7 ------------------------------------- 8
46 // \ 4 ------------------------- 5 /
48 // \ \ 1 ------- 2 / /
66 scomplex::simplex ssring
[] = {
67 // dimension, faces, value...
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}
93 const size_t cntssring
= sizeof(ssring
)/sizeof(scomplex::simplex
);
94 ring
.add_simplices(ssring
, cntssring
);
97 scomplex::simplex sstorus
[] = {
98 // dimension, faces, value...
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},
130 const size_t cntsstorus
= sizeof(sstorus
)/sizeof(scomplex::simplex
);
132 torus
.add_simplices(sstorus
, cntsstorus
);
135 scomplex::simplex ssball
[] = {
136 // dimension, faces, value...
150 {2, {6, 10, 11}, 14},
151 {2, {7, 11, 12}, 15},
152 {2, {8, 12, 13}, 16},
153 {2, {9, 13, 10}, 17},
157 {2, {18, 10, 12}, 21},
158 {3, {21, 14, 15, 19}, 22},
159 {3, {21, 16, 17, 20}, 23},
161 const size_t cntssball
= sizeof(ssball
)/sizeof(scomplex::simplex
);
162 ball
.add_simplices(ssball
, cntssball
);
166 virtual void tear_down() {
169 void test_matrix_reduction() {
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
);
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
);
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
);
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
);
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
);
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;
215 void test_betti_numbers() {
217 ringp
.compute_diagrams();
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());
230 torusp
.compute_diagrams();
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());
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());
248 ballp
.compute_diagrams();
250 //std::cout << std::endl;
251 //ballp.interpret_reduction(std::cout);
252 //std::cout << std::endl;
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());
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);
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);
278 void test_lowestones() {
279 test_lowestones_impl(oring
);
280 test_lowestones_impl(otorus
);
281 test_lowestones_impl(oball
);
284 void test_lowestones_impl(const scomplex::simplex_order
& so
) {
286 p
.compute_matrices();
287 const lom
&lowestones
= p
.get_lowestones_matrix();
289 // Check that lowestones contains every column-index only once,
291 lom cpy
= lowestones
;
292 sort(cpy
.begin(), cpy
.end());
293 lom::iterator it
= cpy
.begin();
295 while (it
!= cpy
.end() && *it
== 0)
297 // From now on, no column index appears twice
298 for (; it
+1 != cpy
.end(); ++it
)
299 TEST_ASSERT(*it
< *(it
+1));
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
306 // (i) the column r is a zero-column
307 // (i) the row c is a zero-column
309 const uint32_t c
= lowestones
[r
];
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);
318 void test_diagram() {
319 test_diagram_impl(oring
);
320 test_diagram_impl(otorus
);
321 test_diagram_impl(oball
);
324 void test_diagram_impl(const scomplex::simplex_order
& so
) {
326 p
.compute_diagrams();
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
);