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::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
;
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
);
40 virtual void setup() {
45 // 7 ------------------------------------- 8
47 // \ 4 ------------------------- 5 /
49 // \ \ 1 ------- 2 / /
67 scomplex::simplex ssring
[] = {
83 /* 15 */ {2, {13, 14, 10}},
87 /* 19 */ {2, {9, 14, 17}},
88 /* 20 */ {2, {7, 16, 18}},
89 /* 21 */ {2, {11, 17, 18}},
91 /* 23 */ {2, {12, 16, 22}},
92 /* 24 */ {2, {8, 13, 22}},
94 const size_t cntssring
= sizeof(ssring
)/sizeof(scomplex::simplex
);
95 ring
.add_simplices(ssring
, cntssring
);
98 scomplex::simplex sstorus
[] = {
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}},
131 const size_t cntsstorus
= sizeof(sstorus
)/sizeof(scomplex::simplex
);
133 torus
.add_simplices(sstorus
, cntsstorus
);
136 scomplex::simplex ssball
[] = {
159 {3, {21, 14, 15, 19}},
160 {3, {21, 16, 17, 20}},
162 const size_t cntssball
= sizeof(ssball
)/sizeof(scomplex::simplex
);
163 ball
.add_simplices(ssball
, cntssball
);
167 virtual void tear_down() {
170 void test_matrix_reduction() {
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
);
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
);
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
);
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
);
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
);
209 sfunction
ftorus(torus
);
210 ftorus
.make_complete();
211 for (unsigned i
=1; i
< torus
.size(); ++i
)
212 ftorus
.set_value(i
, i
);
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
;
224 void test_betti_numbers() {
226 ringp
.compute_diagrams();
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());
239 torusp
.compute_diagrams();
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());
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());
257 ballp
.compute_diagrams();
259 //std::cout << std::endl;
260 //ballp.interpret_reduction(std::cout);
261 //std::cout << std::endl;
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());
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);
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);
287 void test_lowestones() {
288 test_lowestones_impl(oring
);
289 test_lowestones_impl(otorus
);
290 test_lowestones_impl(oball
);
293 void test_lowestones_impl(const scomplex::simplex_order
& so
) {
295 p
.compute_matrices();
296 const lom
&lowestones
= p
.get_lowestones_matrix();
298 // Check that lowestones contains every column-index only once,
300 lom cpy
= lowestones
;
301 sort(cpy
.begin(), cpy
.end());
302 lom::iterator it
= cpy
.begin();
304 while (it
!= cpy
.end() && *it
== 0)
306 // From now on, no column index appears twice
307 for (; it
+1 != cpy
.end(); ++it
)
308 TEST_ASSERT(*it
< *(it
+1));
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
315 // (i) the column r is a zero-column
316 // (i) the row c is a zero-column
318 const uint32_t c
= lowestones
[r
];
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);
327 void test_diagram() {
328 test_diagram_impl(oring
);
329 test_diagram_impl(otorus
);
330 test_diagram_impl(oball
);
333 void test_diagram_impl(const scomplex::simplex_order
& so
) {
335 p
.compute_diagrams();
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
);