Persistent_cohomology_interface.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author(s): Vincent Rouvreau
4 *
5 * Copyright (C) 2016 Inria
6 *
7 * Modification(s):
8 * - YYYY/MM Author: Description of the modification
9 */
10
11#ifndef INCLUDE_PERSISTENT_COHOMOLOGY_INTERFACE_H_
12#define INCLUDE_PERSISTENT_COHOMOLOGY_INTERFACE_H_
13
14#include <gudhi/Persistent_cohomology.h>
15
16#include <cstdlib>
17#include <vector>
18#include <utility> // for std::pair
19#include <algorithm> // for sort
20#include <unordered_map>
21
22namespace Gudhi {
23
24template<class FilteredComplex>
25class Persistent_cohomology_interface : public
26persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp> {
27 private:
28 typedef persistent_cohomology::Persistent_cohomology<FilteredComplex, persistent_cohomology::Field_Zp> Base;
29 /*
30 * Compare two intervals by dimension, then by length.
31 */
32 struct cmp_intervals_by_dim_then_length {
33 template<typename Persistent_interval>
34 bool operator()(const Persistent_interval & p1, const Persistent_interval & p2) {
35 if (std::get<0>(p1) == std::get<0>(p2)) {
36 auto& i1 = std::get<1>(p1);
37 auto& i2 = std::get<1>(p2);
38 return std::get<1>(i1) - std::get<0>(i1) > std::get<1>(i2) - std::get<0>(i2);
39 }
40 else
41 return (std::get<0>(p1) > std::get<0>(p2));
42 // Why does this sort by decreasing dimension?
43 }
44 };
45
46 public:
47 Persistent_cohomology_interface(FilteredComplex* stptr, bool persistence_dim_max=false)
48 : Base(*stptr, persistence_dim_max),
49 stptr_(stptr) { }
50
51 // TODO: move to the constructors?
52 void compute_persistence(int homology_coeff_field, double min_persistence) {
53 Base::init_coefficients(homology_coeff_field);
55 }
56
57 std::vector<std::pair<int, std::pair<double, double>>> get_persistence() {
58 std::vector<std::pair<int, std::pair<double, double>>> persistence;
59 auto const& persistent_pairs = Base::get_persistent_pairs();
60 persistence.reserve(persistent_pairs.size());
61 for (auto pair : persistent_pairs) {
62 persistence.emplace_back(stptr_->dimension(get<0>(pair)),
63 std::make_pair(stptr_->filtration(get<0>(pair)),
64 stptr_->filtration(get<1>(pair))));
65 }
66 // Custom sort and output persistence
67 cmp_intervals_by_dim_then_length cmp;
68 std::sort(std::begin(persistence), std::end(persistence), cmp);
69 return persistence;
70 }
71
72 // This function computes the top-dimensional cofaces associated to the positive and negative
73 // simplices of a cubical complex. The output format is a vector of vectors of three integers,
74 // which are [homological dimension, index of top-dimensional coface of positive simplex,
75 // index of top-dimensional coface of negative simplex]. If the topological feature is essential,
76 // then the index of top-dimensional coface of negative simplex is arbitrarily set to -1.
77 std::vector<std::vector<int>> cofaces_of_cubical_persistence_pairs() {
78
79 // Warning: this function is meant to be used with CubicalComplex only!!
80
81 auto&& pairs = persistent_cohomology::Persistent_cohomology<FilteredComplex,
82 persistent_cohomology::Field_Zp>::get_persistent_pairs();
83
84 // Gather all top-dimensional cells and store their simplex handles
85 std::vector<std::size_t> max_splx;
86 for (auto splx : stptr_->top_dimensional_cells_range())
87 max_splx.push_back(splx);
88 // Sort these simplex handles and compute the ordering function
89 // This function allows to go directly from the simplex handle to the position of the corresponding top-dimensional cell in the input data
90 std::unordered_map<std::size_t, int> order;
91 //std::sort(max_splx.begin(), max_splx.end());
92 for (unsigned int i = 0; i < max_splx.size(); i++) order.emplace(max_splx[i], i);
93
94 std::vector<std::vector<int>> persistence_pairs;
95 for (auto pair : pairs) {
96 int h = stptr_->dimension(get<0>(pair));
97 // Recursively get the top-dimensional cell / coface associated to the persistence generator
98 std::size_t face0 = stptr_->get_top_dimensional_coface_of_a_cell(get<0>(pair));
99 // Retrieve the index of the corresponding top-dimensional cell in the input data
100 int splx0 = order[face0];
101
102 int splx1 = -1;
103 if (get<1>(pair) != stptr_->null_simplex()){
104 // Recursively get the top-dimensional cell / coface associated to the persistence generator
105 std::size_t face1 = stptr_->get_top_dimensional_coface_of_a_cell(get<1>(pair));
106 // Retrieve the index of the corresponding top-dimensional cell in the input data
107 splx1 = order[face1];
108 }
109 persistence_pairs.push_back({ h, splx0, splx1 });
110 }
111 return persistence_pairs;
112 }
113
114 std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs() {
115 std::vector<std::pair<std::vector<int>, std::vector<int>>> persistence_pairs;
116 auto const& pairs = Base::get_persistent_pairs();
117 persistence_pairs.reserve(pairs.size());
118 std::vector<int> birth;
119 std::vector<int> death;
120 for (auto pair : pairs) {
121 birth.clear();
122 if (get<0>(pair) != stptr_->null_simplex()) {
123 for (auto vertex : stptr_->simplex_vertex_range(get<0>(pair))) {
124 birth.push_back(vertex);
125 }
126 }
127
128 death.clear();
129 if (get<1>(pair) != stptr_->null_simplex()) {
130 death.reserve(birth.size()+1);
131 for (auto vertex : stptr_->simplex_vertex_range(get<1>(pair))) {
132 death.push_back(vertex);
133 }
134 }
135
136 persistence_pairs.emplace_back(birth, death);
137 }
138 return persistence_pairs;
139 }
140
141 // TODO: (possibly at the python level)
142 // - an option to return only some of those vectors?
143 typedef std::pair<std::vector<std::vector<int>>, std::vector<std::vector<int>>> Generators;
144
145 Generators lower_star_generators() {
146 Generators out;
147 // diags[i] should be interpreted as vector<array<int,2>>
148 auto& diags = out.first;
149 // diagsinf[i] should be interpreted as vector<int>
150 auto& diagsinf = out.second;
151 for (auto pair : Base::get_persistent_pairs()) {
152 auto s = std::get<0>(pair);
153 auto t = std::get<1>(pair);
154 int dim = stptr_->dimension(s);
155 auto v = stptr_->vertex_with_same_filtration(s);
156 if(t == stptr_->null_simplex()) {
157 while(diagsinf.size() < dim+1) diagsinf.emplace_back();
158 diagsinf[dim].push_back(v);
159 } else {
160 while(diags.size() < dim+1) diags.emplace_back();
161 auto w = stptr_->vertex_with_same_filtration(t);
162 auto& d = diags[dim];
163 d.insert(d.end(), { v, w });
164 }
165 }
166 return out;
167 }
168
169 // An alternative, to avoid those different sizes, would be to "pad" vertex generator v as (v, v) or (v, -1). When using it as index, this corresponds to adding the vertex filtration values either on the diagonal of the distance matrix, or as an extra row or column.
170 // We could also merge the vectors for different dimensions into a single one, with an extra column for the dimension (converted to type double).
171 Generators flag_generators() {
172 Generators out;
173 // diags[0] should be interpreted as vector<array<int,3>> and other diags[i] as vector<array<int,4>>
174 auto& diags = out.first;
175 // diagsinf[0] should be interpreted as vector<int> and other diagsinf[i] as vector<array<int,2>>
176 auto& diagsinf = out.second;
177 for (auto pair : Base::get_persistent_pairs()) {
178 auto s = std::get<0>(pair);
179 auto t = std::get<1>(pair);
180 int dim = stptr_->dimension(s);
181 bool infinite = t == stptr_->null_simplex();
182 if(infinite) {
183 if(dim == 0) {
184 auto v = *std::begin(stptr_->simplex_vertex_range(s));
185 if(diagsinf.size()==0)diagsinf.emplace_back();
186 diagsinf[0].push_back(v);
187 } else {
188 auto e = stptr_->edge_with_same_filtration(s);
189 auto&& e_vertices = stptr_->simplex_vertex_range(e);
190 auto i = std::begin(e_vertices);
191 auto v1 = *i;
192 auto v2 = *++i;
193 GUDHI_CHECK(++i==std::end(e_vertices), "must be an edge");
194 while(diagsinf.size() < dim+1) diagsinf.emplace_back();
195 auto& d = diagsinf[dim];
196 d.insert(d.end(), { v1, v2 });
197 }
198 } else {
199 auto et = stptr_->edge_with_same_filtration(t);
200 auto&& et_vertices = stptr_->simplex_vertex_range(et);
201 auto it = std::begin(et_vertices);
202 auto w1 = *it;
203 auto w2 = *++it;
204 GUDHI_CHECK(++it==std::end(et_vertices), "must be an edge");
205 if(dim == 0) {
206 auto v = *std::begin(stptr_->simplex_vertex_range(s));
207 if(diags.size()==0)diags.emplace_back();
208 auto& d = diags[0];
209 d.insert(d.end(), { v, w1, w2 });
210 } else {
211 auto es = stptr_->edge_with_same_filtration(s);
212 auto&& es_vertices = stptr_->simplex_vertex_range(es);
213 auto is = std::begin(es_vertices);
214 auto v1 = *is;
215 auto v2 = *++is;
216 GUDHI_CHECK(++is==std::end(es_vertices), "must be an edge");
217 while(diags.size() < dim+1) diags.emplace_back();
218 auto& d = diags[dim];
219 d.insert(d.end(), { v1, v2, w1, w2 });
220 }
221 }
222 }
223 return out;
224 }
225
226 private:
227 // A copy
228 FilteredComplex* stptr_;
229};
230
231} // namespace Gudhi
232
233#endif // INCLUDE_PERSISTENT_COHOMOLOGY_INTERFACE_H_
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:172
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:156
const std::vector< Persistent_interval > & get_persistent_pairs() const
Returns a list of persistence birth and death FilteredComplex::Simplex_handle pairs.
Definition: Persistent_cohomology.h:672
The concept FilteredComplex describes the requirements for a type to implement a filtered cell comple...
Definition: FilteredComplex.h:17
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Tue Aug 16 2022 14:01:50 for GUDHIdev by Doxygen 1.9.4