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 simplicial_complex
<MAXDIM
, IT
, VT
> scomplex
;
27 typedef typename
scomplex::simplex_order simplex_order
;
28 typedef typename
simplex_order::boundary_matrix boundary_matrix
;
29 typedef boolean_colmatrix
<IT
> transformation_matrix
;
31 /** A point in a persistence diagram comprises birth- and death-index.
32 * These are the indices of the simplices of 'order' that give birth
33 * and death to the class. */
34 struct diagram_point
{
35 /** 'birth'-th simplex in 'order' gave birth to this cycle. */
37 /** 'death'-th simplex in 'order' killed this cycle. If 'death' is
38 * zero, the cycle never dies. Otherwise, 'death' is always larger
43 /** Save the births and deaths of simplices of a fixed dimension, say p. */
45 /** Saves the sequence of birth/death-points of p-simplices, sorted
46 * by birth-index. If death-index is zero, the corresponding
47 * homology class never dies. */
48 std::vector
<diagram_point
> births
;
51 /** Gives the p-th Betti number, i.e., the number of immortal
52 * p-dimensional homology classes. */
53 size_t betti() const {
56 for (unsigned i
=0; i
< births
.size(); ++i
)
57 if (births
[i
].death
== 0)
63 /** Gives the number of p-dimensional homology classes that are
64 * born by simplex 'from' or earlier and die after simplex 'to' is
66 size_t persistent_betti(IT from
, IT to
) const {
70 for (unsigned i
=0; i
< births
.size(); ++i
) {
71 // All following simplices are born too late.
72 if (births
[i
].birth
> from
)
74 if (births
[i
].death
> to
|| births
[i
].death
== 0)
83 /** Create a new peristence object on the given simplex_order */
84 persistence(const simplex_order
&order
) :
94 /** Reset all results gained from 'order'. */
96 done_matrices
= false;
97 done_diagrams
= false;
100 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
101 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
103 void compute_matrices() {
106 done_matrices
= true;
108 bm
= order
.get_boundary_matrix();
110 tm
= create_unit_matrix
<transformation_matrix
>(bm
.size());
112 reduce_boundary_matrix(rbm
, 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 /** Get the lowest-one matrix of 'order'. */
121 const boundary_matrix
& get_lowestones_matrix() const {
122 assert(done_matrices
);
126 /** Get the boundary matrix of 'order'. */
127 const boundary_matrix
& get_boundary_matrix() const {
128 assert(done_matrices
);
132 /** Get the reduced boundary matrix of 'order'. */
133 const boundary_matrix
& get_reduced_boundary_matrix() const {
134 assert(done_matrices
);
138 /** Get the transformation matrix of 'order', which transforms the
139 * boundary matrix to the reduced boundary matrix, when multiplied from
141 const transformation_matrix
& get_transformation_matrix() const {
142 assert(done_matrices
);
146 /** Print the sequence of births and deaths of cycles (homology classes). */
147 void interpret_reduction(std::ostream
&os
) const {
148 assert(done_matrices
);
150 for (unsigned c
=0; c
< bm
.size(); ++c
) {
151 os
<< c
<< ". inserting ";
153 switch (bm
.get_column(c
).size()) {
155 os
<< "dummy vertex of dim -1";
175 os
<< " boundary: " << bm
.get_column(c
) << std::endl
;
177 typename boolean_colrowmatrix
<IT
>::colbase::column_type col
= rbm
.get_column(c
);
178 if (col
.size() == 0) {
179 os
<< " \033[1;32mbirth\033[0;m of a cycle: " << tm
.get_column(c
) << std::endl
;
181 os
<< " \033[1;31mdeath\033[0;m of a cycle: " << rbm
.get_column(c
) << std::endl
;
182 os
<< " boundary of: " << tm
.get_column(c
) << std::endl
;
188 /** Get the 'dim'-dimensional peristence diagram. */
189 const diagram
& get_persistence_diagram(int dim
) const{
191 assert(dim
<= MAXDIM
);
192 assert(done_diagrams
);
194 return diagrams
[dim
+1];
197 /** Compute persistence diagrams */
198 void compute_diagrams() {
201 done_diagrams
= true;
206 std::set
<IT
> deaths
, births
;
209 for (int d
=0; d
< MAXDIM
+ 2; ++d
)
210 diagrams
[d
] = diagram();
212 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
213 if (lowestones
.get_column(c
).size() == 0) {
214 const int dim
= order
.get_simplex(c
).dim
;
216 // Create a diagram point
220 // If the class actually dies, get the index of the killing simplex
221 if (lowestones
.get_row(c
).size() > 0)
222 p
.death
= lowestones
.get_row(c
).back();
224 _get_persistence_diagram(dim
).births
.push_back(p
);
226 assert(births
.count(p
.birth
) == 0);
227 assert(p
.death
== 0 || deaths
.count(p
.death
) == 0);
228 births
.insert(p
.birth
);
229 deaths
.insert(p
.death
);
234 // Shakespeare principle: A simplex either gives birth or kills, be
235 // or not be, tertium non datur.
236 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
237 if (lowestones
.get_column(c
).size() != 0)
238 assert(c
== 0 || deaths
.count(c
) > 0);
240 assert(births
.count(c
) > 0);
242 // Even Jesus died after his birth. Homology classes are nothing better.
243 for (int d
=-1; d
<= MAXDIM
; ++d
) {
244 const diagram
&dia
= get_persistence_diagram(d
);
245 for (unsigned i
=0; i
< dia
.births
.size(); ++i
)
246 if (dia
.births
[i
].death
> 0)
247 assert(dia
.births
[i
].birth
< dia
.births
[i
].death
);
252 /** Print birth and death information in the persistence diagrams. */
253 void interpret_persistence(std::ostream
&os
) const {
254 assert(done_diagrams
);
256 for (int d
=-1; d
<= MAXDIM
; ++d
) {
257 os
<< "Dimension " << d
<< std::endl
;
258 const diagram
&dia
= get_persistence_diagram(d
);
259 for (unsigned i
=0; i
< dia
.births
.size(); ++i
) {
261 if (dia
.births
[i
].death
> 0)
262 os
<< dia
.births
[i
].birth
<< "\033[1;31m -> " << dia
.births
[i
].death
;
264 os
<< "\033[1;32m" << dia
.births
[i
].birth
;
265 os
<< "\033[0;m" << std::endl
;
272 diagram
& _get_persistence_diagram(int dim
) {
274 assert(dim
<= MAXDIM
);
275 assert(done_diagrams
);
277 return diagrams
[dim
+1];
280 /** The underlying simplex order of this diagram. */
281 const simplex_order
&order
;
286 /** The boundary matrix of 'order' */
288 /** The reduced boundary matrix of 'order' */
290 /** Only the lowest ones per column of the reduced boundary matrix.
291 * Each one at (r, c) corresponds to the death of a class when
292 * inserting simplex c, which was born when inserting simplex. */
293 boundary_matrix lowestones
;
294 /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
295 transformation_matrix tm
;
297 /** The container for all p-dimensional diagrams. The p-th diagram is
298 * save at index p+1, with -1 <= p <= MAXDIM. */
299 diagram diagrams
[MAXDIM
+2];
303 /** Takes the boundary matrix b and transforms it inplace into a reduced matrix
304 * b. The matrix v is required to be a unit matrix having the same dimensions
305 * as b. It is maintained in order to leave the product b*inverse(v) as an
306 * invariant. Hence, the resulting b is equal to the product of the boundary
308 template<class IT
, class D
>
309 void reduce_boundary_matrix(boolean_colrowmatrix
<IT
> &b
, boolean_colmatrix_base
<IT
, D
> &v
) {
310 assert(b
.size() == v
.width());
312 // Make every column reduced, i.e., it is a zero-column or contains only one 1.
313 for (unsigned c
=0; c
< b
.size(); ++c
) {
314 const typename boolean_colrowmatrix
<IT
>::colbase::column_type
&col
= b
.get_column(c
);
318 // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
322 // Get a copy of the r-th row
323 typename boolean_colrowmatrix
<IT
>::rowbase::row_type
row(b
.get_row(r
).size());
324 copy(b
.get_row(r
).begin(), b
.get_row(r
).end(), row
.begin());
325 assert(row
.size() >= 1);
327 // Get rid of 1s at that row right of column c.
328 typename boolean_colrowmatrix
<IT
>::row_type::const_iterator it
= lower_bound(row
.begin(), row
.end(), c
);
329 for(; it
!= row
.end(); ++it
) {
333 assert(b
.get(r
, *it
));
334 b
.add_column(*it
, col
);
335 v
.add_column(*it
, v
.get_column(c
));
336 assert(!b
.get(r
, *it
));
342 std::ostream
& operator<<(std::ostream
& os
, const std::vector
<T
> &vec
) {
345 typename
std::vector
<T
>::const_iterator it
= vec
.begin();
346 while ( it
!= vec
.end()) {
348 if (++it
!= vec
.end())
357 } // namespace libstick