]> git.sthu.org Git - libstick.git/blob - include/libstick-0.1/persistence.h
White space editing in booleanmatrix.h
[libstick.git] / include / libstick-0.1 / persistence.h
1 #ifndef persistence_h_aPeinoXeiFeethuz
2 #define persistence_h_aPeinoXeiFeethuz
3
4 #include <iostream>
5 #include <ostream>
6 #include <vector>
7 #include <cassert>
8
9 #include "simplicialcomplex.h"
10 #include "booleanmatrix.h"
11
12
13 namespace libstick {
14
15
16 /** Persistence computers various persistent homology information on a
17 * simplex_order of a simplicial_complex. */
18 template<int MAXDIM, class IT, class VT>
19 class persistence {
20
21 public:
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;
27
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. */
33 IT birth;
34 /** 'death'-th simplex in 'order' killed this cycle. If 'death' is
35 * zero, the cycle never dies. Otherwise, 'death' is always larger
36 * than 'birth'. */
37 IT death;
38 };
39
40 /** Save the births and deaths of simplices of a fixed dimension, say p. */
41 struct diagram {
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;
46
47
48 /** Gives the p-th Betti number, i.e., the number of immortal
49 * p-dimensional homology classes. */
50 size_t betti() const {
51 size_t b = 0;
52
53 for (unsigned i=0; i < births.size(); ++i)
54 if (births[i].death == 0)
55 b++;
56
57 return b;
58 }
59
60 /** Gives the number of p-dimensional homology classes that are
61 * born by simplex 'from' or earlier and die after simplex 'to' is
62 * inserted. */
63 size_t persistent_betti(IT from, IT to) const {
64 assert(from <= to);
65 size_t b = 0;
66
67 for (unsigned i=0; i < births.size(); ++i) {
68 // All following simplices are born too late.
69 if (births[i].birth > from)
70 break;
71 if (births[i].death > to || births[i].death == 0)
72 b++;
73 }
74
75 return b;
76 }
77 };
78
79 public:
80 /** Create a new peristence object on the given simplex_order */
81 persistence(const simplex_order &order) :
82 order(order),
83 bm(0),
84 rbm(0),
85 lowestones(0),
86 tm(0)
87 {
88 reset();
89 }
90
91 /** Reset all results gained from 'order'. */
92 void reset() {
93 done_matrices = false;
94 done_diagrams = false;
95 }
96
97 /** Get simplicial order of this persistence object. */
98 const simplex_order& get_order() const {
99 return order;
100 }
101
102 /** Get boundary matrix 'bm' of 'order', compute reduces boundary
103 * matrix 'rbm', and the transformation matrix 'tm' such that rbm = bm *
104 * tm. */
105 void compute_matrices() {
106 if (done_matrices)
107 return;
108 done_matrices = true;
109
110 bm = order.get_boundary_matrix();
111 rbm = bm;
112 tm = create_unit_matrix<transformation_matrix>(bm.width());
113 lowestones = lowestones_matrix(bm.width());
114
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) {
117 //if (c % 100 == 0)
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];
125 assert(c2 < c);
126
127 // There is no valid c2, hence we have completely reduced
128 if (c2 == 0)
129 break;
130 // Reduce col by column c2
131 else {
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));
136 }
137 }
138
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);
143 assert(r < c);
144 lowestones[r] = c;
145 }
146 }
147 }
148
149 /** Get the lowest-one matrix of 'order'. */
150 const lowestones_matrix& get_lowestones_matrix() const {
151 assert(done_matrices);
152 return lowestones;
153 }
154
155 /** Get the boundary matrix of 'order'. */
156 const boundary_matrix& get_boundary_matrix() const {
157 assert(done_matrices);
158 return bm;
159 }
160
161 /** Get the reduced boundary matrix of 'order'. */
162 const boundary_matrix& get_reduced_boundary_matrix() const {
163 assert(done_matrices);
164 return rbm;
165 }
166
167 /** Get the transformation matrix of 'order', which transforms the
168 * boundary matrix to the reduced boundary matrix, when multiplied from
169 * the right. */
170 const transformation_matrix& get_transformation_matrix() const {
171 assert(done_matrices);
172 return tm;
173 }
174
175 /** Print the sequence of births and deaths of cycles (homology classes). */
176 void interpret_reduction(std::ostream &os) const {
177 assert(done_matrices);
178
179 for (unsigned c=0; c < bm.width(); ++c) {
180 os << c << ". inserting ";
181
182 switch (bm.get_column(c).size()) {
183 case 0:
184 os << "dummy vertex of dim -1";
185 break;
186 case 1:
187 os << "vertex";
188 break;
189 case 2:
190 os << "edge";
191 break;
192 case 3:
193 os << "triangle";
194 break;
195 case 4:
196 os << "tetrahedron";
197 break;
198 default:
199 os << "simplex";
200 break;
201 }
202
203 os << "(value: " << order.get_simplex(c).value << ")" << std::endl;
204 os << " boundary: " << bm.get_column(c) << std::endl;
205
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;
209 } else {
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;
212 }
213 os << std::endl;
214 }
215 }
216
217 /** Get the 'dim'-dimensional peristence diagram. */
218 const diagram& get_persistence_diagram(int dim) const{
219 assert(-1 <= dim);
220 assert(dim <= MAXDIM);
221 assert(done_diagrams);
222
223 return diagrams[dim+1];
224 }
225
226 /** Compute persistence diagrams */
227 void compute_diagrams() {
228 if (done_diagrams)
229 return;
230 done_diagrams = true;
231
232 compute_matrices();
233
234 #ifndef NDEBUG
235 std::set<IT> deaths, births;
236 #endif
237 //Clear all diagrams
238 for (int d=0; d < MAXDIM + 2; ++d)
239 diagrams[d] = diagram();
240
241 for (unsigned c=0; c < lowestones.size(); ++c) {
242 // A cycle was born
243 if (rbm.get_column(c).size() == 0) {
244 const int dim = order.get_simplex(c).dim;
245
246 // Create a diagram point
247 diagram_point p;
248 p.birth = c;
249 p.death = 0;
250 // If the class actually dies, get the index of the killing simplex
251 if (lowestones[c] > 0)
252 p.death = lowestones[c];
253
254 _get_persistence_diagram(dim).births.push_back(p);
255 #ifndef NDEBUG
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);
260 #endif
261 }
262 }
263 #ifndef NDEBUG
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);
270 else
271 assert(births.count(c) > 0);
272 }
273 #endif
274 }
275
276 /** Print birth and death information in the persistence diagrams. */
277 void interpret_persistence(std::ostream &os) const {
278 assert(done_diagrams);
279
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) {
284 os << " ";
285 if (dia.births[i].death > 0)
286 os << dia.births[i].birth << "\033[1;31m -> " << dia.births[i].death;
287 else
288 os << "\033[1;32m" << dia.births[i].birth;
289 os << "\033[0;m" << std::endl;
290 }
291 }
292 }
293
294 private:
295
296 diagram& _get_persistence_diagram(int dim) {
297 assert(-1 <= dim);
298 assert(dim <= MAXDIM);
299 assert(done_diagrams);
300
301 return diagrams[dim+1];
302 }
303
304 /** The underlying simplex order of this diagram. */
305 const simplex_order &order;
306
307 bool done_matrices;
308 bool done_diagrams;
309
310 /** The boundary matrix of 'order' */
311 boundary_matrix bm;
312 /** The reduced boundary matrix of 'order' */
313 boundary_matrix rbm;
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;
320
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];
324 };
325
326
327 } // namespace libstick
328
329 #endif