1 #ifndef persistence_h_Fahlaewahgaebaqu
2 #define persistence_h_Fahlaewahgaebaqu
5 #include <cpptest-suite.h>
7 #include <libstick/persistence.h>
9 using namespace libstick
;
12 class persistence_TestSuite
: public Test::Suite
{
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
;
23 scomplex ball
, ring
, torus
;
24 scomplex::simplex_order oball
, oring
, otorus
;
27 persistence_TestSuite() :
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
);
41 virtual void setup() {
46 // 7 ------------------------------------- 8
48 // \ 4 ------------------------- 5 /
50 // \ \ 1 ------- 2 / /
68 scomplex::simplex ssring
[] = {
84 /* 15 */ {2, {13, 14, 10}},
88 /* 19 */ {2, {9, 14, 17}},
89 /* 20 */ {2, {7, 16, 18}},
90 /* 21 */ {2, {11, 17, 18}},
92 /* 23 */ {2, {12, 16, 22}},
93 /* 24 */ {2, {8, 13, 22}},
95 const size_t cntssring
= sizeof(ssring
)/sizeof(scomplex::simplex
);
96 ring
.add_simplices(ssring
, cntssring
);
99 scomplex::simplex sstorus
[] = {
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}},
132 const size_t cntsstorus
= sizeof(sstorus
)/sizeof(scomplex::simplex
);
134 torus
.add_simplices(sstorus
, cntsstorus
);
137 scomplex::simplex ssball
[] = {
160 {3, {21, 14, 15, 19}},
161 {3, {21, 16, 17, 20}},
163 const size_t cntssball
= sizeof(ssball
)/sizeof(scomplex::simplex
);
164 ball
.add_simplices(ssball
, cntssball
);
168 virtual void tear_down() {
171 void test_matrix_reduction() {
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
);
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
);
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
);
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
);
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
);
210 sfunction
ftorus(torus
);
211 ftorus
.make_complete();
212 for (unsigned i
=1; i
< torus
.size(); ++i
)
213 ftorus
.set_value(i
, i
);
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
;
225 void test_diagrampoint_value() {
227 torusp
.compute_diagrams();
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
);
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
));
247 void test_betti_numbers() {
249 ringp
.compute_diagrams();
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());
262 torusp
.compute_diagrams();
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());
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());
280 ballp
.compute_diagrams();
282 //std::cout << std::endl;
283 //ballp.interpret_reduction(std::cout);
284 //std::cout << std::endl;
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());
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);
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);
310 void test_lowestones() {
311 test_lowestones_impl(oring
);
312 test_lowestones_impl(otorus
);
313 test_lowestones_impl(oball
);
316 void test_lowestones_impl(const scomplex::simplex_order
& so
) {
318 p
.compute_matrices();
319 const lom
&lowestones
= p
.get_lowestones_matrix();
321 // Check that lowestones contains every column-index only once,
323 lom cpy
= lowestones
;
324 sort(cpy
.begin(), cpy
.end());
325 lom::iterator it
= cpy
.begin();
327 while (it
!= cpy
.end() && *it
== 0)
329 // From now on, no column index appears twice
330 for (; it
+1 != cpy
.end(); ++it
)
331 TEST_ASSERT(*it
< *(it
+1));
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
338 // (i) the column r is a zero-column
339 // (i) the row c is a zero-column
341 const uint32_t c
= lowestones
[r
];
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);
350 void test_diagram() {
351 test_diagram_impl(oring
);
352 test_diagram_impl(otorus
);
353 test_diagram_impl(oball
);
356 void test_diagram_impl(const scomplex::simplex_order
& so
) {
358 p
.compute_diagrams();
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
);