1 #ifndef persistence_h_aPeinoXeiFeethuz
2 #define persistence_h_aPeinoXeiFeethuz
9 #include "simplicialcomplex.h"
10 #include "booleanmatrix.h"
17 std::ostream
& operator<<(std::ostream
&, const std::vector
<T
> &);
20 /** Persistence computers various persistent homology information on a
21 * simplex_order of a simplicial_complex. */
22 template<int MAXDIM
, class IT
, class VT
>
26 typedef typename simplicial_complex
<MAXDIM
, IT
, VT
>::simplex_order simplex_order
;
27 typedef typename
simplex_order::boundary_matrix boundary_matrix
;
28 typedef boolean_colmatrix
<IT
> transformation_matrix
;
30 /** A point in a persistence diagram comprises birth- and death-index.
31 * These are the indices of the simplices of 'order' that give birth
32 * and death to the class. */
33 struct diagram_point
{
34 /** 'birth'-th simplex in 'order' gave birth to this cycle. */
36 /** 'death'-th simplex in 'order' killed this cycle. If 'death' is
37 * zero, the cycle never dies. Otherwise, 'death' is always larger
42 /** Save the births and deaths of simplices of a fixed dimension, say p. */
44 /** Saves the sequence of birth/death-points of p-simplices, sorted
45 * by birth-index. If death-index is zero, the corresponding
46 * homology class never dies. */
47 std::vector
<diagram_point
> births
;
50 /** Gives the p-th Betti number, i.e., the number of immortal
51 * p-dimensional homology classes. */
52 size_t betti() const {
55 for (unsigned i
=0; i
< births
.size(); ++i
)
56 if (births
[i
].death
== 0)
62 /** Gives the number of p-dimensional homology classes that are
63 * born by simplex 'from' or earlier and die after simplex 'to' is
65 size_t persistent_betti(IT from
, IT to
) const {
69 for (unsigned i
=0; i
< births
.size(); ++i
) {
70 // All following simplices are born too late.
71 if (births
[i
].birth
> from
)
73 if (births
[i
].death
> to
|| births
[i
].death
== 0)
82 /** Create a new peristence object on the given simplex_order */
83 persistence(const simplex_order
&order
) :
93 /** Reset all results gained from 'order'. */
95 done_matrices
= false;
96 done_diagrams
= false;
99 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
100 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
102 void compute_matrices() {
105 done_matrices
= true;
107 bm
= order
.get_boundary_matrix();
109 tm
= create_unit_matrix
<transformation_matrix
>(bm
.size());
111 reduce_boundary_matrix(rbm
, tm
);
112 assert(rbm
== bm
* tm
);
114 lowestones
= boundary_matrix(rbm
.size());
115 for (unsigned c
=0; c
< rbm
.size(); ++c
)
116 if (rbm
.get_column(c
).size() > 0)
117 lowestones
.set(rbm
.get_column(c
).back(), c
, true);
120 for (unsigned c
=0; c
< lowestones
.size(); ++c
)
121 assert(lowestones
.get_column(c
).size() <= 1);
122 for (unsigned r
=0; r
< lowestones
.size(); ++r
)
123 assert(lowestones
.get_row(r
).size() <= 1);
124 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
125 // If (r,c) is a one then
126 // (i) a cycle dies with c --> row c is zero
127 // (ii) a cycle is born with r --> column r is zero
129 // (i) the column r is a zero-column
130 // (i) the row c is a zero-column
131 if (lowestones
.get_column(c
).size() > 0) {
132 const IT r
= lowestones
.get_column(c
).back();
133 assert(lowestones
.get_column(r
).size() == 0);
134 assert(lowestones
.get_row(c
).size() == 0);
140 /** Get the boundary matrix of 'order'. */
141 const boundary_matrix
& get_boundary_matrix() const {
142 assert(done_matrices
);
146 /** Get the reduced boundary matrix of 'order'. */
147 const boundary_matrix
& get_reduced_boundary_matrix() const {
148 assert(done_matrices
);
152 /** Get the transformation matrix of 'order', which transforms the
153 * boundary matrix to the reduced boundary matrix, when multiplied from
155 const transformation_matrix
& get_transformation_matrix() const {
156 assert(done_matrices
);
160 /** Print the sequence of births and deaths of cycles (homology classes). */
161 void interpret_reduction(std::ostream
&os
) const {
162 assert(done_matrices
);
164 for (unsigned c
=0; c
< bm
.size(); ++c
) {
165 os
<< c
<< ". inserting ";
167 switch (bm
.get_column(c
).size()) {
169 os
<< "dummy vertex of dim -1";
189 os
<< " boundary: " << bm
.get_column(c
) << std::endl
;
191 typename boolean_colrowmatrix
<IT
>::colbase::column_type col
= rbm
.get_column(c
);
192 if (col
.size() == 0) {
193 os
<< " \033[1;32mbirth\033[0;m of a cycle: " << tm
.get_column(c
) << std::endl
;
195 os
<< " \033[1;31mdeath\033[0;m of a cycle: " << rbm
.get_column(c
) << std::endl
;
196 os
<< " boundary of: " << tm
.get_column(c
) << std::endl
;
202 /** Get the 'dim'-dimensional peristence diagram. */
203 const diagram
& get_persistence_diagram(int dim
) const{
205 assert(dim
<= MAXDIM
);
206 assert(done_diagrams
);
208 return diagrams
[dim
+1];
211 /** Compute persistence diagrams */
212 void compute_diagrams() {
215 done_diagrams
= true;
220 std::set
<IT
> deaths
, births
;
223 for (int d
=0; d
< MAXDIM
+ 2; ++d
)
224 diagrams
[d
] = diagram();
226 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
227 if (lowestones
.get_column(c
).size() == 0) {
228 const int dim
= order
.get_simplex(c
).dim
;
230 // Create a diagram point
234 // If the class actually dies, get the index of the killing simplex
235 if (lowestones
.get_row(c
).size() > 0)
236 p
.death
= lowestones
.get_row(c
).back();
238 _get_persistence_diagram(dim
).births
.push_back(p
);
240 assert(births
.count(p
.birth
) == 0);
241 assert(p
.death
== 0 || deaths
.count(p
.death
) == 0);
242 births
.insert(p
.birth
);
243 deaths
.insert(p
.death
);
248 // Shakespeare principle: A simplex either gaves birth or kills, be
249 // or not be, tertium non datur.
250 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
251 if (lowestones
.get_column(c
).size() != 0)
252 assert(c
== 0 || deaths
.count(c
) > 0);
254 assert(births
.count(c
) > 0);
256 // Even Jesus died after his birth. Homology classes are nothing better.
257 for (int d
=-1; d
<= MAXDIM
; ++d
) {
258 const diagram
&dia
= get_persistence_diagram(d
);
259 for (unsigned i
=0; i
< dia
.births
.size(); ++i
)
260 if (dia
.births
[i
].death
> 0)
261 assert(dia
.births
[i
].birth
< dia
.births
[i
].death
);
266 /** Print birth and death information in the persistence diagrams. */
267 void interpret_persistence(std::ostream
&os
) const {
268 assert(done_diagrams
);
270 for (int d
=-1; d
<= MAXDIM
; ++d
) {
271 os
<< "Dimension " << d
<< std::endl
;
272 const diagram
&dia
= get_persistence_diagram(d
);
273 for (unsigned i
=0; i
< dia
.births
.size(); ++i
) {
275 if (dia
.births
[i
].death
> 0)
276 os
<< dia
.births
[i
].birth
<< "\033[1;31m -> " << dia
.births
[i
].death
;
278 os
<< "\033[1;32m" << dia
.births
[i
].birth
;
279 os
<< "\033[0;m" << std::endl
;
286 diagram
& _get_persistence_diagram(int dim
) {
288 assert(dim
<= MAXDIM
);
289 assert(done_diagrams
);
291 return diagrams
[dim
+1];
294 /** The underlying simplex order of this diagram. */
295 const simplex_order
&order
;
300 /** The boundary matrix of 'order' */
302 /** The reduced boundary matrix of 'order' */
304 /** Only the lowest ones per column of the reduced boundary matrix.
305 * Each one at (r, c) corresponds to the death of a class when
306 * inserting simplex c, which was born when inserting simplex. */
307 boundary_matrix lowestones
;
308 /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
309 transformation_matrix tm
;
311 /** The container for all p-dimensional diagrams. The p-th diagram is
312 * save at index p+1, with -1 <= p <= MAXDIM. */
313 diagram diagrams
[MAXDIM
+2];
317 /** Takes the boundary matrix b and transforms it inplace into a reduced matrix
318 * b. The matrix v is required to be a unit matrix having the same dimensions
319 * as b. It is maintained in order to leave the product b*inverse(v) as an
320 * invariant. Hence, the resulting b is equal to the product of the boundary
322 template<class IT
, class D
>
323 void reduce_boundary_matrix(boolean_colrowmatrix
<IT
> &b
, boolean_colmatrix_base
<IT
, D
> &v
) {
324 assert(b
.size() == v
.width());
326 // Make every column reduced, i.e., it is a zero-column or contains only one 1.
327 for (unsigned c
=0; c
< b
.size(); ++c
) {
328 const typename boolean_colrowmatrix
<IT
>::colbase::column_type
&col
= b
.get_column(c
);
332 // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
336 // Get a copy of the r-th row
337 typename boolean_colrowmatrix
<IT
>::rowbase::row_type
row(b
.get_row(r
).size());
338 copy(b
.get_row(r
).begin(), b
.get_row(r
).end(), row
.begin());
339 assert(row
.size() >= 1);
341 // Get rid of 1s at that row right of column c.
342 typename boolean_colrowmatrix
<IT
>::row_type::const_iterator it
= lower_bound(row
.begin(), row
.end(), c
);
343 for(; it
!= row
.end(); ++it
) {
347 assert(b
.get(r
, *it
));
348 b
.add_column(*it
, col
);
349 v
.add_column(*it
, v
.get_column(c
));
350 assert(!b
.get(r
, *it
));
356 std::ostream
& operator<<(std::ostream
& os
, const std::vector
<T
> &vec
) {
359 typename
std::vector
<T
>::const_iterator it
= vec
.begin();
360 while ( it
!= vec
.end()) {
362 if (++it
!= vec
.end())
371 } // namespace libstick