1 #ifndef persistence_h_aPeinoXeiFeethuz
2 #define persistence_h_aPeinoXeiFeethuz
8 #include "simplicialcomplex.h"
9 #include "booleanmatrix.h"
16 std::ostream
& operator<<(std::ostream
&, const std::vector
<T
> &);
19 /** Persistence computers various persistent homology information on a
20 * simplex_order of a simplicial_complex. */
21 template<int MAXDIM
, class IT
, class VT
>
25 typedef typename simplicial_complex
<MAXDIM
, IT
, VT
>::simplex_order simplex_order
;
26 typedef typename
simplex_order::boundary_matrix boundary_matrix
;
27 typedef boolean_colmatrix
<IT
> transformation_matrix
;
29 /** A point in a diagram comprises birth- and death-index. These are the
30 * indices of the simplices of 'order' that give birth and death to the
32 struct diagram_point
{
39 /** Create a new peristence object on the given simplex_order */
40 persistence(const simplex_order
&order
) :
49 /** Reset all results gained from 'order'. */
51 done_matrices
= false;
54 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
55 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
57 void compute_matrices() {
62 bm
= order
.get_boundary_matrix();
64 tm
= create_unit_matrix
<transformation_matrix
>(bm
.size());
65 reduce_boundary_matrix(rbm
, tm
);
67 assert(rbm
== bm
* tm
);
70 /** Get the boundary matrix of 'order'. */
71 const boundary_matrix
& get_boundary_matrix() {
76 /** Get the reduced boundary matrix of 'order'. */
77 const boundary_matrix
& get_reduced_boundary_matrix() {
82 /** Get the transformation matrix of 'order', which transforms the
83 * boundary matrix to the reduced boundary matrix, when multiplied from
85 const transformation_matrix
& get_transformation_matrix() {
91 /** Print the sequence of births and deaths of cycles (homology classes). */
92 void interpret_reduction(std::ostream
&os
) {
94 for (unsigned c
=0; c
< bm
.size(); ++c
) {
95 os
<< c
<< ". inserting ";
97 switch (bm
.get_column(c
).size()) {
99 os
<< "dummy vertex of dim -1";
119 os
<< " boundary: " << bm
.get_column(c
) << std::endl
;
121 typename boolean_colrowmatrix
<IT
>::colbase::column_type col
= rbm
.get_column(c
);
122 if (col
.size() == 0) {
123 os
<< " \033[1;32mbirth\033[0;m of a cycle: " << tm
.get_column(c
) << std::endl
;
125 os
<< " \033[1;31mdeath\033[0;m of a cycle: " << rbm
.get_column(c
) << std::endl
;
126 os
<< " boundary of: " << tm
.get_column(c
) << std::endl
;
134 /** The underlying simplex order of this diagram. */
135 const simplex_order
&order
;
141 transformation_matrix tm
;
145 /** Takes the boundary matrix b and transforms it inplace into a reduced matrix
146 * b. The matrix v is required to be a unit matrix having the same dimensions
147 * as b. It is maintained in order to leave the product b*inverse(v) as an
148 * invariant. Hence, the resulting b is equal to the product of the boundary
150 template<class IT
, class D
>
151 void reduce_boundary_matrix(boolean_colrowmatrix
<IT
> &b
, boolean_colmatrix_base
<IT
, D
> &v
) {
152 assert(b
.size() == v
.width());
154 // Make every column reduced, i.e., it is a zero-column or contains only one 1.
155 for (unsigned c
=0; c
< b
.size(); ++c
) {
156 const typename boolean_colrowmatrix
<IT
>::colbase::column_type
&col
= b
.get_column(c
);
160 // The column is not a zero-column, take the lowest 1, and reduce the other columns at the row
164 // Get a copy of the r-th row
165 typename boolean_colrowmatrix
<IT
>::rowbase::row_type
row(b
.get_row(r
).size());
166 copy(b
.get_row(r
).begin(), b
.get_row(r
).end(), row
.begin());
167 assert(row
.size() >= 1);
169 // Get rid of 1s at that row right of column c.
170 typename boolean_colrowmatrix
<IT
>::row_type::const_iterator it
= lower_bound(row
.begin(), row
.end(), c
);
171 for(; it
!= row
.end(); ++it
) {
175 assert(b
.get(r
, *it
));
176 b
.add_column(*it
, col
);
177 v
.add_column(*it
, v
.get_column(c
));
178 assert(!b
.get(r
, *it
));
184 std::ostream
& operator<<(std::ostream
& os
, const std::vector
<T
> &vec
) {
187 typename
std::vector
<T
>::const_iterator it
= vec
.begin();
188 while ( it
!= vec
.end()) {
190 if (++it
!= vec
.end())
199 } // namespace libstick