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 simplicial order of this persistence object. */
98 const simplex_order
& get_order() const {
102 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
103 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
105 void compute_matrices() {
108 done_matrices
= true;
110 bm
= order
.get_boundary_matrix();
112 tm
= create_unit_matrix
<transformation_matrix
>(bm
.width());
113 lowestones
= lowestones_matrix(bm
.width());
115 // Make every column reduced, i.e., it is a zero-column or contains only one 1.
116 for (unsigned c
=0; c
< rbm
.width(); ++c
) {
118 //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
119 // Reduce as long as we need to reduce
120 while (rbm
.get_column(c
).size() > 0) {
121 // (r, c) is the lowest one of col
122 const IT r
= rbm
.get_column(c
).back();
123 // This column contains its lowest one on the same row
124 const IT c2
= lowestones
[r
];
127 // There is no valid c2, hence we have completely reduced
130 // Reduce col by column c2
132 assert(rbm
.get(r
, c
));
133 rbm
.add_column(c
, rbm
.get_column(c2
));
134 tm
.add_column(c
, tm
.get_column(c2
));
135 assert(!rbm
.get(r
, c
));
139 // A lowest one remained, recall it
140 if (rbm
.get_column(c
).size() > 0) {
141 const IT r
= rbm
.get_column(c
).back();
142 assert(lowestones
[r
] == 0);
149 /** Get the lowest-one matrix of 'order'. */
150 const lowestones_matrix
& get_lowestones_matrix() const {
151 assert(done_matrices
);
155 /** Get the boundary matrix of 'order'. */
156 const boundary_matrix
& get_boundary_matrix() const {
157 assert(done_matrices
);
161 /** Get the reduced boundary matrix of 'order'. */
162 const boundary_matrix
& get_reduced_boundary_matrix() const {
163 assert(done_matrices
);
167 /** Get the transformation matrix of 'order', which transforms the
168 * boundary matrix to the reduced boundary matrix, when multiplied from
170 const transformation_matrix
& get_transformation_matrix() const {
171 assert(done_matrices
);
175 /** Print the sequence of births and deaths of cycles (homology classes). */
176 void interpret_reduction(std::ostream
&os
) const {
177 assert(done_matrices
);
179 for (unsigned c
=0; c
< bm
.width(); ++c
) {
180 os
<< c
<< ". inserting ";
182 switch (bm
.get_column(c
).size()) {
184 os
<< "dummy vertex of dim -1";
203 os
<< "(value: " << order
.get_simplex(c
).value
<< ")" << std::endl
;
204 os
<< " boundary: " << bm
.get_column(c
) << std::endl
;
206 typename boolean_colmatrix
<IT
>::column_type col
= rbm
.get_column(c
);
207 if (col
.size() == 0) {
208 os
<< " \033[1;32mbirth\033[0;m of a cycle: " << tm
.get_column(c
) << std::endl
;
210 os
<< " \033[1;31mdeath\033[0;m of a cycle: " << rbm
.get_column(c
) << std::endl
;
211 os
<< " boundary of: " << tm
.get_column(c
) << std::endl
;
217 /** Get the 'dim'-dimensional peristence diagram. */
218 const diagram
& get_persistence_diagram(int dim
) const{
220 assert(dim
<= MAXDIM
);
221 assert(done_diagrams
);
223 return diagrams
[dim
+1];
226 /** Compute persistence diagrams */
227 void compute_diagrams() {
230 done_diagrams
= true;
235 std::set
<IT
> deaths
, births
;
238 for (int d
=0; d
< MAXDIM
+ 2; ++d
)
239 diagrams
[d
] = diagram();
241 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
243 if (rbm
.get_column(c
).size() == 0) {
244 const int dim
= order
.get_simplex(c
).dim
;
246 // Create a diagram point
250 // If the class actually dies, get the index of the killing simplex
251 if (lowestones
[c
] > 0)
252 p
.death
= lowestones
[c
];
254 _get_persistence_diagram(dim
).births
.push_back(p
);
256 assert(births
.count(p
.birth
) == 0);
257 assert(p
.death
== 0 || deaths
.count(p
.death
) == 0);
258 births
.insert(p
.birth
);
259 deaths
.insert(p
.death
);
264 // Shakespeare principle: A simplex either gives birth or kills, be
265 // or not be, tertium non datur.
266 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
267 // A class died with c
268 if (rbm
.get_column(c
).size() != 0)
269 assert(c
== 0 || deaths
.count(c
) > 0);
271 assert(births
.count(c
) > 0);
276 /** Print birth and death information in the persistence diagrams. */
277 void interpret_persistence(std::ostream
&os
) const {
278 assert(done_diagrams
);
280 for (int d
=-1; d
<= MAXDIM
; ++d
) {
281 os
<< "Dimension " << d
<< std::endl
;
282 const diagram
&dia
= get_persistence_diagram(d
);
283 for (unsigned i
=0; i
< dia
.births
.size(); ++i
) {
285 if (dia
.births
[i
].death
> 0)
286 os
<< dia
.births
[i
].birth
<< "\033[1;31m -> " << dia
.births
[i
].death
;
288 os
<< "\033[1;32m" << dia
.births
[i
].birth
;
289 os
<< "\033[0;m" << std::endl
;
296 diagram
& _get_persistence_diagram(int dim
) {
298 assert(dim
<= MAXDIM
);
299 assert(done_diagrams
);
301 return diagrams
[dim
+1];
304 /** The underlying simplex order of this diagram. */
305 const simplex_order
&order
;
310 /** The boundary matrix of 'order' */
312 /** The reduced boundary matrix of 'order' */
314 /** Only the lowest ones per column of the reduced boundary matrix.
315 * lowestones[i] contains the column-index at which the lowest one at
316 * row i appears, and zero if no such column exists. */
317 lowestones_matrix lowestones
;
318 /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
319 transformation_matrix tm
;
321 /** The container for all p-dimensional diagrams. The p-th diagram is
322 * save at index p+1, with -1 <= p <= MAXDIM. */
323 diagram diagrams
[MAXDIM
+2];
327 } // namespace libstick