summaryrefslogtreecommitdiff
path: root/src/python/include/Simplex_tree_interface.h
blob: 629f60832899c3e8fe8d21fe6c432e2b2ca9b012 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
 *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
 *    Author(s):       Vincent Rouvreau
 *
 *    Copyright (C) 2016 Inria
 *
 *    Modification(s):
 *      - YYYY/MM Author: Description of the modification
 */

#ifndef INCLUDE_SIMPLEX_TREE_INTERFACE_H_
#define INCLUDE_SIMPLEX_TREE_INTERFACE_H_

#include <gudhi/graph_simplicial_complex.h>
#include <gudhi/distance_functions.h>
#include <gudhi/Simplex_tree.h>
#include <gudhi/Points_off_io.h>
#ifdef GUDHI_USE_EIGEN3
#include <gudhi/Flag_complex_edge_collapser.h>
#endif

#include <iostream>
#include <vector>
#include <utility>  // std::pair
#include <tuple>
#include <iterator>  // for std::distance

namespace Gudhi {

template<typename SimplexTreeOptions = Simplex_tree_options_full_featured>
class Simplex_tree_interface : public Simplex_tree<SimplexTreeOptions> {
 public:
  using Base = Simplex_tree<SimplexTreeOptions>;
  using Filtration_value = typename Base::Filtration_value;
  using Vertex_handle = typename Base::Vertex_handle;
  using Simplex_handle = typename Base::Simplex_handle;
  using Insertion_result = typename std::pair<Simplex_handle, bool>;
  using Simplex = std::vector<Vertex_handle>;
  using Simplex_and_filtration = std::pair<Simplex, Filtration_value>;
  using Filtered_simplices = std::vector<Simplex_and_filtration>;
  using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator;
  using Complex_simplex_iterator = typename Base::Complex_simplex_iterator;
  using Extended_filtration_data = typename Base::Extended_filtration_data;
  using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator;

 public:

  Extended_filtration_data efd;
  
  bool find_simplex(const Simplex& simplex) {
    return (Base::find(simplex) != Base::null_simplex());
  }

  void assign_simplex_filtration(const Simplex& simplex, Filtration_value filtration) {
    Base::assign_filtration(Base::find(simplex), filtration);
    Base::clear_filtration();
  }

  bool insert(const Simplex& simplex, Filtration_value filtration = 0) {
    Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration);
    if (result.first != Base::null_simplex())
      Base::clear_filtration();
    return (result.second);
  }

  // Do not interface this function, only used in alpha complex interface for complex creation
  bool insert_simplex(const Simplex& simplex, Filtration_value filtration = 0) {
    Insertion_result result = Base::insert_simplex(simplex, filtration);
    return (result.second);
  }

  // Do not interface this function, only used in interface for complex creation
  bool insert_simplex_and_subfaces(const Simplex& simplex, Filtration_value filtration = 0) {
    Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration);
    return (result.second);
  }

  // Do not interface this function, only used in strong witness interface for complex creation
  bool insert_simplex(const std::vector<std::size_t>& simplex, Filtration_value filtration = 0) {
    Insertion_result result = Base::insert_simplex(simplex, filtration);
    return (result.second);
  }

  // Do not interface this function, only used in strong witness interface for complex creation
  bool insert_simplex_and_subfaces(const std::vector<std::size_t>& simplex, Filtration_value filtration = 0) {
    Insertion_result result = Base::insert_simplex_and_subfaces(simplex, filtration);
    return (result.second);
  }

  Filtration_value simplex_filtration(const Simplex& simplex) {
    return Base::filtration(Base::find(simplex));
  }

  void remove_maximal_simplex(const Simplex& simplex) {
    Base::remove_maximal_simplex(Base::find(simplex));
    Base::clear_filtration();
  }

  Simplex_and_filtration get_simplex_and_filtration(Simplex_handle f_simplex) {
    Simplex simplex;
    for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
      simplex.insert(simplex.begin(), vertex);
    }
    return std::make_pair(std::move(simplex), Base::filtration(f_simplex));
  }

  Filtered_simplices get_star(const Simplex& simplex) {
    Filtered_simplices star;
    for (auto f_simplex : Base::star_simplex_range(Base::find(simplex))) {
      Simplex simplex_star;
      for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
        simplex_star.insert(simplex_star.begin(), vertex);
      }
      star.push_back(std::make_pair(simplex_star, Base::filtration(f_simplex)));
    }
    return star;
  }

  Filtered_simplices get_cofaces(const Simplex& simplex, int dimension) {
    Filtered_simplices cofaces;
    for (auto f_simplex : Base::cofaces_simplex_range(Base::find(simplex), dimension)) {
      Simplex simplex_coface;
      for (auto vertex : Base::simplex_vertex_range(f_simplex)) {
        simplex_coface.insert(simplex_coface.begin(), vertex);
      }
      cofaces.push_back(std::make_pair(simplex_coface, Base::filtration(f_simplex)));
    }
    return cofaces;
  }

  void compute_extended_filtration() {
    this->efd = this->extend_filtration();
    return;
  }

  std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> compute_extended_persistence_subdiagrams(const std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>& dgm, Filtration_value min_persistence){
    std::vector<std::vector<std::pair<int, std::pair<Filtration_value, Filtration_value>>>> new_dgm(4);
    for (unsigned int i = 0; i < dgm.size(); i++){
      std::pair<Filtration_value, Extended_simplex_type> px = this->decode_extended_filtration(dgm[i].second.first, this->efd);
      std::pair<Filtration_value, Extended_simplex_type> py = this->decode_extended_filtration(dgm[i].second.second, this->efd);
      std::pair<int, std::pair<Filtration_value, Filtration_value>> pd_point = std::make_pair(dgm[i].first, std::make_pair(px.first, py.first));
      if(std::abs(px.first - py.first) > min_persistence){
        //Ordinary
        if (px.second == Extended_simplex_type::UP && py.second == Extended_simplex_type::UP){
          new_dgm[0].push_back(pd_point);
        }
        // Relative
        else if (px.second == Extended_simplex_type::DOWN && py.second == Extended_simplex_type::DOWN){
          new_dgm[1].push_back(pd_point);
        }
        else{
          // Extended+
          if (px.first < py.first){
            new_dgm[2].push_back(pd_point);
          }
          //Extended-
          else{
            new_dgm[3].push_back(pd_point);
          }
        }
      }
    }
    return new_dgm;
  }

  Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) {
#ifdef GUDHI_USE_EIGEN3
    using Filtered_edge = std::tuple<Vertex_handle, Vertex_handle, Filtration_value>;
    std::vector<Filtered_edge> edges;
    for (Simplex_handle sh : Base::skeleton_simplex_range(1)) {
      if (Base::dimension(sh) == 1) {
        typename Base::Simplex_vertex_range rg = Base::simplex_vertex_range(sh);
        auto vit = rg.begin();
        Vertex_handle v = *vit;
        Vertex_handle w = *++vit;
        edges.emplace_back(v, w, Base::filtration(sh));
      }
    }

    for (int iteration = 0; iteration < nb_collapse_iteration; iteration++) {
      edges = Gudhi::collapse::flag_complex_collapse_edges(edges);
    }
    Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface();
    // Copy the original 0-skeleton
    for (Simplex_handle sh : Base::skeleton_simplex_range(0)) {
      collapsed_stree_ptr->insert({*(Base::simplex_vertex_range(sh).begin())}, Base::filtration(sh));
    }
    // Insert remaining edges
    for (auto remaining_edge : edges) {
      collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge));
    }
    return collapsed_stree_ptr;
#else
    throw std::runtime_error("Unable to collapse edges as it requires Eigen3 >= 3.1.0.");
#endif
  }

  // Iterator over the simplex tree
  Complex_simplex_iterator get_simplices_iterator_begin() {
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    return Base::complex_simplex_range().begin();
  }

  Complex_simplex_iterator get_simplices_iterator_end() {
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    return Base::complex_simplex_range().end();
  }

  typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_begin() {
    // Base::initialize_filtration(); already performed in filtration_simplex_range
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    return Base::filtration_simplex_range().begin();
  }

  typename std::vector<Simplex_handle>::const_iterator get_filtration_iterator_end() {
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    return Base::filtration_simplex_range().end();
  }

  Skeleton_simplex_iterator get_skeleton_iterator_begin(int dimension) {
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    return Base::skeleton_simplex_range(dimension).begin();
  }

  Skeleton_simplex_iterator get_skeleton_iterator_end(int dimension) {
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    return Base::skeleton_simplex_range(dimension).end();
  }

  std::pair<Boundary_simplex_iterator, Boundary_simplex_iterator> get_boundary_iterators(const Simplex& simplex) {
    auto bd_sh = Base::find(simplex);
    if (bd_sh == Base::null_simplex())
      throw std::runtime_error("simplex not found - cannot find boundaries");
    // this specific case works because the range is just a pair of iterators - won't work if range was a vector
    auto boundary_srange = Base::boundary_simplex_range(bd_sh);
    return std::make_pair(boundary_srange.begin(), boundary_srange.end());
  }
};

}  // namespace Gudhi

#endif  // INCLUDE_SIMPLEX_TREE_INTERFACE_H_