1 #ifndef persistence_h_aPeinoXeiFeethuz
2 #define persistence_h_aPeinoXeiFeethuz
9 #include "simplicialcomplex.h"
10 #include "booleanmatrix.h"
16 /** Persistence computers various persistent homology information on a
17 * simplex_order of a simplicial_complex. */
18 template<int MAXDIM
, class IT
, class VT
>
22 typedef simplicial_complex
<MAXDIM
, IT
, VT
> scomplex
;
23 typedef typename
scomplex::simplex_order simplex_order
;
24 typedef boolean_colmatrix
<IT
> boundary_matrix
;
25 typedef boolean_colmatrix
<IT
> transformation_matrix
;
26 typedef std::vector
<IT
> lowestones_matrix
;
28 /** A point in a persistence diagram comprises birth- and death-index.
29 * These are the indices of the simplices of 'order' that give birth
30 * and death to the class. */
31 struct diagram_point
{
32 /** 'birth'-th simplex in 'order' gave birth to this cycle. */
34 /** 'death'-th simplex in 'order' killed this cycle. If 'death' is
35 * zero, the cycle never dies. Otherwise, 'death' is always larger
40 /** Save the births and deaths of simplices of a fixed dimension, say p. */
42 /** Saves the sequence of birth/death-points of p-simplices, sorted
43 * by birth-index. If death-index is zero, the corresponding
44 * homology class never dies. */
45 std::vector
<diagram_point
> births
;
48 /** Gives the p-th Betti number, i.e., the number of immortal
49 * p-dimensional homology classes. */
50 size_t betti() const {
53 for (unsigned i
=0; i
< births
.size(); ++i
)
54 if (births
[i
].death
== 0)
60 /** Gives the number of p-dimensional homology classes that are
61 * born by simplex 'from' or earlier and die after simplex 'to' is
63 size_t persistent_betti(IT from
, IT to
) const {
67 for (unsigned i
=0; i
< births
.size(); ++i
) {
68 // All following simplices are born too late.
69 if (births
[i
].birth
> from
)
71 if (births
[i
].death
> to
|| births
[i
].death
== 0)
80 /** Create a new peristence object on the given simplex_order */
81 persistence(const simplex_order
&order
) :
91 /** Reset all results gained from 'order'. */
93 done_matrices
= false;
94 done_diagrams
= false;
97 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
98 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
100 void compute_matrices() {
103 done_matrices
= true;
105 bm
= order
.get_boundary_matrix();
107 tm
= create_unit_matrix
<transformation_matrix
>(bm
.width());
108 lowestones
= lowestones_matrix(bm
.width());
110 // Make every column reduced, i.e., it is a zero-column or contains only one 1.
111 for (unsigned c
=0; c
< rbm
.width(); ++c
) {
113 //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
114 // Reduce as long as we need to reduce
115 while (rbm
.get_column(c
).size() > 0) {
116 // (r, c) is the lowest one of col
117 const IT r
= rbm
.get_column(c
).back();
118 // This column contains its lowest one on the same row
119 const IT c2
= lowestones
[r
];
122 // There is no valid c2, hence we have completely reduced
125 // Reduce col by column c2
127 assert(rbm
.get(r
, c
));
128 rbm
.add_column(c
, rbm
.get_column(c2
));
129 tm
.add_column(c
, tm
.get_column(c2
));
130 assert(!rbm
.get(r
, c
));
134 // A lowest one remained, recall it
135 if (rbm
.get_column(c
).size() > 0) {
136 const IT r
= rbm
.get_column(c
).back();
137 assert(lowestones
[r
] == 0);
144 /** Get the lowest-one matrix of 'order'. */
145 const lowestones_matrix
& get_lowestones_matrix() const {
146 assert(done_matrices
);
150 /** Get the boundary matrix of 'order'. */
151 const boundary_matrix
& get_boundary_matrix() const {
152 assert(done_matrices
);
156 /** Get the reduced boundary matrix of 'order'. */
157 const boundary_matrix
& get_reduced_boundary_matrix() const {
158 assert(done_matrices
);
162 /** Get the transformation matrix of 'order', which transforms the
163 * boundary matrix to the reduced boundary matrix, when multiplied from
165 const transformation_matrix
& get_transformation_matrix() const {
166 assert(done_matrices
);
170 /** Print the sequence of births and deaths of cycles (homology classes). */
171 void interpret_reduction(std::ostream
&os
) const {
172 assert(done_matrices
);
174 for (unsigned c
=0; c
< bm
.width(); ++c
) {
175 os
<< c
<< ". inserting ";
177 switch (bm
.get_column(c
).size()) {
179 os
<< "dummy vertex of dim -1";
199 os
<< " boundary: " << bm
.get_column(c
) << std::endl
;
201 typename boolean_colmatrix
<IT
>::column_type col
= rbm
.get_column(c
);
202 if (col
.size() == 0) {
203 os
<< " \033[1;32mbirth\033[0;m of a cycle: " << tm
.get_column(c
) << std::endl
;
205 os
<< " \033[1;31mdeath\033[0;m of a cycle: " << rbm
.get_column(c
) << std::endl
;
206 os
<< " boundary of: " << tm
.get_column(c
) << std::endl
;
212 /** Get the 'dim'-dimensional peristence diagram. */
213 const diagram
& get_persistence_diagram(int dim
) const{
215 assert(dim
<= MAXDIM
);
216 assert(done_diagrams
);
218 return diagrams
[dim
+1];
221 /** Compute persistence diagrams */
222 void compute_diagrams() {
225 done_diagrams
= true;
230 std::set
<IT
> deaths
, births
;
233 for (int d
=0; d
< MAXDIM
+ 2; ++d
)
234 diagrams
[d
] = diagram();
236 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
238 if (rbm
.get_column(c
).size() == 0) {
239 const int dim
= order
.get_simplex(c
).dim
;
241 // Create a diagram point
245 // If the class actually dies, get the index of the killing simplex
246 if (lowestones
[c
] > 0)
247 p
.death
= lowestones
[c
];
249 _get_persistence_diagram(dim
).births
.push_back(p
);
251 assert(births
.count(p
.birth
) == 0);
252 assert(p
.death
== 0 || deaths
.count(p
.death
) == 0);
253 births
.insert(p
.birth
);
254 deaths
.insert(p
.death
);
259 // Shakespeare principle: A simplex either gives birth or kills, be
260 // or not be, tertium non datur.
261 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
262 // A class died with c
263 if (rbm
.get_column(c
).size() != 0)
264 assert(c
== 0 || deaths
.count(c
) > 0);
266 assert(births
.count(c
) > 0);
271 /** Print birth and death information in the persistence diagrams. */
272 void interpret_persistence(std::ostream
&os
) const {
273 assert(done_diagrams
);
275 for (int d
=-1; d
<= MAXDIM
; ++d
) {
276 os
<< "Dimension " << d
<< std::endl
;
277 const diagram
&dia
= get_persistence_diagram(d
);
278 for (unsigned i
=0; i
< dia
.births
.size(); ++i
) {
280 if (dia
.births
[i
].death
> 0)
281 os
<< dia
.births
[i
].birth
<< "\033[1;31m -> " << dia
.births
[i
].death
;
283 os
<< "\033[1;32m" << dia
.births
[i
].birth
;
284 os
<< "\033[0;m" << std::endl
;
291 diagram
& _get_persistence_diagram(int dim
) {
293 assert(dim
<= MAXDIM
);
294 assert(done_diagrams
);
296 return diagrams
[dim
+1];
299 /** The underlying simplex order of this diagram. */
300 const simplex_order
&order
;
305 /** The boundary matrix of 'order' */
307 /** The reduced boundary matrix of 'order' */
309 /** Only the lowest ones per column of the reduced boundary matrix.
310 * lowestones[i] contains the column-index at which the lowest one at
311 * row i appears, and zero if no such column exists. */
312 lowestones_matrix lowestones
;
313 /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
314 transformation_matrix tm
;
316 /** The container for all p-dimensional diagrams. The p-th diagram is
317 * save at index p+1, with -1 <= p <= MAXDIM. */
318 diagram diagrams
[MAXDIM
+2];
322 } // namespace libstick