1 #ifndef persistence_h_aPeinoXeiFeethuz
2 #define persistence_h_aPeinoXeiFeethuz
9 #include "simplicialfunction.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
>
22 typedef simplicial_complex
<MAXDIM
, IT
> 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 /** Get the value of the birth index w.r.t. to a simplicial order o
41 * and a simplicial function f. Requires that o and f are defined
42 * on the same simplicial complex. */
44 VT
get_birth_value(const simplex_order
&o
, const simplicial_function
<MAXDIM
, IT
, VT
> &f
) const {
45 assert(&f
.get_complex() == &o
.get_complex());
46 return f
.get_value(o
.order_to_complex(birth
));
49 /** Get the value of the death index w.r.t. to a simplicial order o
50 * and a simplicial function f. Requires that o and f are defined
51 * on the same simplicial complex. */
53 VT
get_death_value(const simplex_order
&o
, const simplicial_function
<MAXDIM
, IT
, VT
> &f
) const {
54 assert(&f
.get_complex() == &o
.get_complex());
55 return f
.get_value(o
.order_to_complex(death
));
58 /** Get the persistence of this point w.r.t. to a simplicial order
59 * o and a simplicial function f. Requires that o and f are defined
60 * on the same simplicial complex. */
62 VT
get_persistence(const simplex_order
&o
, const simplicial_function
<MAXDIM
, IT
, VT
> &f
) const {
63 assert(&f
.get_complex() == &o
.get_complex());
64 return get_death_value
<VT
>(o
, f
) - get_birth_value
<VT
>(o
, f
);
68 /** Save the births and deaths of simplices of a fixed dimension, say p. */
70 /** Saves the sequence of birth/death-points of p-simplices, sorted
71 * by birth-index. If death-index is zero, the corresponding
72 * homology class never dies. */
73 std::vector
<diagram_point
> births
;
76 /** Gives the p-th Betti number, i.e., the number of immortal
77 * p-dimensional homology classes. */
78 size_t betti() const {
81 for (unsigned i
=0; i
< births
.size(); ++i
)
82 if (births
[i
].death
== 0)
88 /** Gives the number of p-dimensional homology classes that are
89 * born by simplex 'from' or earlier and die after simplex 'to' is
91 size_t persistent_betti(IT from
, IT to
) const {
95 for (unsigned i
=0; i
< births
.size(); ++i
) {
96 // All following simplices are born too late.
97 if (births
[i
].birth
> from
)
99 if (births
[i
].death
> to
|| births
[i
].death
== 0)
108 /** Create a new peristence object on the given simplex_order */
109 persistence(const simplex_order
&order
) :
119 /** Reset all results gained from 'order'. */
121 done_matrices
= false;
122 done_diagrams
= false;
125 /** Get simplicial order of this persistence object. */
126 const simplex_order
& get_order() const {
130 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
131 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
133 void compute_matrices() {
136 done_matrices
= true;
138 bm
= order
.get_boundary_matrix();
140 tm
= create_unit_matrix
<transformation_matrix
>(bm
.width());
141 lowestones
= lowestones_matrix(bm
.width());
143 // Make every column reduced, i.e., it is a zero-column or contains only one 1.
144 for (unsigned c
=0; c
< rbm
.width(); ++c
) {
146 //std::cout << "c = " << c << " (" << (100.0*float(c)/rbm.width()) << " %)" << std::endl;
147 // Reduce as long as we need to reduce
148 while (!rbm
.get_column(c
).isempty()) {
149 // (r, c) is the lowest one of col
150 const IT r
= rbm
.get_column(c
).back();
151 // This column contains its lowest one on the same row
152 const IT c2
= lowestones
[r
];
155 // There is no valid c2, hence we have completely reduced
158 // Reduce col by column c2
160 assert(rbm
.get(r
, c
));
161 rbm
.add_column(c
, rbm
.get_column(c2
));
162 tm
.add_column(c
, tm
.get_column(c2
));
163 assert(!rbm
.get(r
, c
));
167 // A lowest one remained, recall it
168 if (!rbm
.get_column(c
).isempty()) {
169 const IT r
= rbm
.get_column(c
).back();
170 assert(lowestones
[r
] == 0);
177 /** Get the lowest-one matrix of 'order'. */
178 const lowestones_matrix
& get_lowestones_matrix() const {
179 assert(done_matrices
);
183 /** Get the boundary matrix of 'order'. */
184 const boundary_matrix
& get_boundary_matrix() const {
185 assert(done_matrices
);
189 /** Get the reduced boundary matrix of 'order'. */
190 const boundary_matrix
& get_reduced_boundary_matrix() const {
191 assert(done_matrices
);
195 /** Get the transformation matrix of 'order', which transforms the
196 * boundary matrix to the reduced boundary matrix, when multiplied from
198 const transformation_matrix
& get_transformation_matrix() const {
199 assert(done_matrices
);
203 /** Print the sequence of births and deaths of cycles (homology
204 * classes). If a non-NULL simplicial function is given, also the
205 * values of the simplices are printed. */
207 void interpret_reduction(std::ostream
&os
, simplicial_function
<MAXDIM
, IT
, VT
> *f
) const {
208 assert(done_matrices
);
210 for (unsigned c
=0; c
< bm
.width(); ++c
) {
211 os
<< c
<< ". inserting ";
213 switch (bm
.get_column(c
).size()) {
215 os
<< "dummy vertex of dim -1";
236 os
<< " value: " << f
->get_value(order
.order_to_complex(c
)) << std::endl
;
237 os
<< " boundary: " << bm
.get_column(c
) << std::endl
;
239 typename boolean_colmatrix
<IT
>::column_type col
= rbm
.get_column(c
);
241 os
<< " \033[1;32mbirth\033[0;m of a cycle: " << tm
.get_column(c
) << std::endl
;
243 os
<< " \033[1;31mdeath\033[0;m of a cycle: " << rbm
.get_column(c
);
244 os
<< ", boundary of: " << tm
.get_column(c
) << std::endl
;
250 /** Print the sequence of births and deaths of cycles (homology classes). */
251 void interpret_reduction(std::ostream
&os
) const {
252 interpret_reduction
<float>(os
, NULL
);
255 /** Get the 'dim'-dimensional peristence diagram. */
256 const diagram
& get_persistence_diagram(int dim
) const{
258 assert(dim
<= MAXDIM
);
259 assert(done_diagrams
);
261 return diagrams
[dim
+1];
264 /** Compute persistence diagrams */
265 void compute_diagrams() {
268 done_diagrams
= true;
273 std::set
<IT
> deaths
, births
;
276 for (int d
=0; d
< MAXDIM
+ 2; ++d
)
277 diagrams
[d
] = diagram();
279 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
281 if (rbm
.get_column(c
).isempty()) {
282 const int dim
= order
.get_simplex(c
).dim
;
284 // Create a diagram point
288 // If the class actually dies, get the index of the killing simplex
289 if (lowestones
[c
] > 0)
290 p
.death
= lowestones
[c
];
292 _get_persistence_diagram(dim
).births
.push_back(p
);
294 assert(births
.count(p
.birth
) == 0);
295 assert(p
.death
== 0 || deaths
.count(p
.death
) == 0);
296 births
.insert(p
.birth
);
297 deaths
.insert(p
.death
);
302 // Shakespeare principle: A simplex either gives birth or kills, be
303 // or not be, tertium non datur.
304 for (unsigned c
=0; c
< lowestones
.size(); ++c
) {
305 // A class died with c
306 if (!rbm
.get_column(c
).isempty())
307 assert(c
== 0 || deaths
.count(c
) > 0);
309 assert(births
.count(c
) > 0);
314 /** Print birth and death information in the persistence diagrams. If a
315 * non-NULL simplicial function is given, also the persistence w.r.t.
316 * this function is printed. */
318 void interpret_persistence(std::ostream
&os
, simplicial_function
<MAXDIM
, IT
, VT
> *f
) const {
319 assert(done_diagrams
);
321 for (int d
=-1; d
<= MAXDIM
; ++d
) {
322 os
<< "Dimension " << d
<< std::endl
;
323 const diagram
&dia
= get_persistence_diagram(d
);
325 for (unsigned i
=0; i
< dia
.births
.size(); ++i
) {
326 const IT birth
= dia
.births
[i
].birth
;
327 const IT death
= dia
.births
[i
].death
;
331 os
<< birth
<< "\033[1;31m -> " << death
<< "\033[0;m";
334 os
<< " \tpers: " << dia
.births
[i
].get_persistence(get_order(), *f
);
336 os
<< "\033[1;32m" << birth
<< "\033[0;m";
338 os
<< " \t\t[essential]";
345 /** Print birth and death information in the persistence diagrams. */
346 void interpret_persistence(std::ostream
&os
) const {
347 interpret_persistence
<float>(os
, NULL
);
352 diagram
& _get_persistence_diagram(int dim
) {
354 assert(dim
<= MAXDIM
);
355 assert(done_diagrams
);
357 return diagrams
[dim
+1];
360 /** The underlying simplex order of this diagram. */
361 const simplex_order
&order
;
366 /** The boundary matrix of 'order' */
368 /** The reduced boundary matrix of 'order' */
370 /** Only the lowest ones per column of the reduced boundary matrix.
371 * lowestones[i] contains the column-index at which the lowest one at
372 * row i appears, and zero if no such column exists. */
373 lowestones_matrix lowestones
;
374 /** The transformation matrix that transforms bm into rbm, i.e. rbm = bm * tm. */
375 transformation_matrix tm
;
377 /** The container for all p-dimensional diagrams. The p-th diagram is
378 * save at index p+1, with -1 <= p <= MAXDIM. */
379 diagram diagrams
[MAXDIM
+2];
383 } // namespace libstick