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 
22 namespace Gudhi {
23 
24 template<class FilteredComplex>
25 class Persistent_cohomology_interface : public
26 persistent_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);
54  Base::compute_persistent_cohomology(min_persistence);
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.1