Model.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): David Salinas
4  *
5  * Copyright (C) 2014 Inria
6  *
7  * Modification(s):
8  * - YYYY/MM Author: Description of the modification
9  * - 2019/07 Vincent Rouvreau: homsimpl (from CHOMP) is no more delivered.
10  * Error message when not previously installed.
11  */
12 
13 #ifndef MODEL_MODEL_H_
14 #define MODEL_MODEL_H_
15 
16 #include <gudhi/Clock.h>
17 #include <gudhi/Skeleton_blocker/Skeleton_blocker_simple_geometric_traits.h>
18 #include <gudhi/Skeleton_blocker_geometric_complex.h>
19 #include <gudhi/Off_reader.h>
20 
21 #include <CGAL/Euclidean_distance.h>
22 
23 #include <fstream>
24 #include <limits>
25 #include <string>
26 #include <vector>
27 
28 #include "utils/UI_utils.h"
29 #include "utils/Lloyd_builder.h"
30 #include "utils/Rips_builder.h"
31 #include "utils/K_nearest_builder.h"
32 #include "utils/Vertex_collapsor.h"
33 #include "utils/Edge_collapsor.h"
34 #include "utils/Edge_contractor.h"
35 #include "utils/Persistence_compute.h"
36 #include "utils/Critical_points.h"
37 #include "utils/Is_manifold.h"
38 
39 #include "Complex_typedefs.h"
40 
41 template<typename Complex>
42 class CGAL_geometric_flag_complex_wrapper {
43  Complex& complex_;
44  typedef typename Complex::Vertex_handle Vertex_handle;
45  typedef typename Complex::Point Point;
46 
47  const bool load_only_points_;
48 
49  public:
50  CGAL_geometric_flag_complex_wrapper(Complex& complex, bool load_only_points = false) :
51  complex_(complex),
52  load_only_points_(load_only_points) { }
53 
54  void init(int dim, int num_vertices, int num_max_faces, int num_edges) const { }
55 
56  void point(const std::vector<double>& coords) {
57  Point p(coords.size(), coords.begin(), coords.end());
58  complex_.add_vertex(p);
59  }
60 
61  void maximal_face(std::vector<int> vertices) {
62  if (!load_only_points_) {
63  // std::clog << "size:" << vertices.size() << std::endl;
64  for (std::size_t i = 0; i < vertices.size(); ++i)
65  for (std::size_t j = i + 1; j < vertices.size(); ++j)
66  complex_.add_edge_without_blockers(Vertex_handle(vertices[i]), Vertex_handle(vertices[j]));
67  }
68  }
69 
70  void done() const { }
71 };
72 
73 class Model {
74  public:
75  Complex complex_;
77 
78  Model() : complex_() { }
79 
80  public:
81  void off_file_open(const std::string& name_file) {
82  UIDBGMSG("load off file", name_file);
83  complex_.clear();
84  CGAL_geometric_flag_complex_wrapper<Complex> read_wraper(complex_, false);
85  Gudhi::read_off(name_file, read_wraper);
86  }
87 
88  void off_points_open(const std::string& name_file) {
89  UIDBGMSG("load off points", name_file);
90  complex_.clear();
91  CGAL_geometric_flag_complex_wrapper<Complex> read_wraper(complex_, true);
92  Gudhi::read_off(name_file, read_wraper);
93  }
94 
95  void off_file_save(const std::string& name_file) {
96  UIDBG("save off file");
97  UIDBG("save off off_points_save");
98  std::ofstream file(name_file);
99  if (file.is_open()) {
100  file << "OFF\n";
101  file << complex_.num_vertices() << " " << complex_.num_edges() << " 0\n";
102  for (auto v : complex_.vertex_range()) {
103  const auto& pt(complex_.point(v));
104  for (auto it = pt.cartesian_begin(); it != pt.cartesian_end(); ++it)
105  file << *it << " ";
106  file << std::endl;
107  }
108  for (auto e : complex_.edge_range())
109  file << "2 " << complex_.first_vertex(e) << " " << complex_.second_vertex(e) << "\n";
110  file.close();
111  } else {
112  std::cerr << "Could not open file " << name_file << std::endl;
113  }
114  }
115 
116  void off_points_save(const std::string& name_file) {
117  UIDBG("save off off_points_save");
118  std::ofstream file(name_file);
119  if (file.is_open()) {
120  file << "OFF\n";
121  file << complex_.num_vertices() << " 0 0\n";
122  for (auto v : complex_.vertex_range()) {
123  const auto& pt(complex_.point(v));
124  for (auto it = pt.cartesian_begin(); it != pt.cartesian_end(); ++it)
125  file << *it << " ";
126  file << std::endl;
127  }
128  file.close();
129  } else {
130  std::cerr << "Could not open file " << name_file << std::endl;
131  }
132  }
133 
134  // point sets operations
135  void uniform_noise(double amplitude) {
136  UIDBG("unif noise");
137  for (auto v : complex_.vertex_range())
138  complex_.point(v) = add_uniform_noise(complex_.point(v), amplitude);
139  }
140 
141  private:
142  Point add_uniform_noise(const Point& point, double amplitude) {
143  std::vector<double> new_point(point.dimension());
144  for (int i = 0; i < point.dimension(); ++i) {
145  new_point[i] = point[i] + (rand() % 2 - .5) * amplitude;
146  }
147  return Point(point.dimension(), new_point.begin(), new_point.end());
148  }
149 
150  public:
151  void lloyd(int num_iterations, int num_closest_neighbors) {
152  UIDBG("lloyd");
153  Lloyd_builder<Complex> lloyd_builder(complex_, 1);
154  }
155 
156  double squared_eucl_distance(const Point& p1, const Point& p2) const {
157  return Geometry_trait::Squared_distance_d()(p1, p2);
158  }
159 
160  // complex operations from points
161 
162  void build_rips(double alpha) {
163  UIDBG("build_rips");
164  Rips_builder<Complex> rips_builder(complex_, alpha);
165  }
166 
167  void build_k_nearest_neighbors(unsigned k) {
168  UIDBG("build_k_nearest");
169  complex_.keep_only_vertices();
170  K_nearest_builder<Complex> k_nearest_builder(complex_, k);
171  }
172 
173  void build_delaunay() {
174  UIDBG("build_delaunay");
175  complex_.keep_only_vertices();
176  }
177 
178  void contract_edges(unsigned num_contractions) {
179  Gudhi::Clock c;
180  Edge_contractor<Complex> contractor(complex_, num_contractions);
181  std::clog << "Time to simplify: " << c.num_seconds() << "s" << std::endl;
182  }
183 
184  void collapse_vertices(unsigned num_collapses) {
185  auto old_num_vertices = complex_.num_vertices();
186  Vertex_collapsor<Complex> collapsor(complex_, complex_.num_vertices());
187  UIDBGMSG("num vertices collapsed:", old_num_vertices - complex_.num_vertices());
188  }
189 
190  void collapse_edges(unsigned num_collapses) {
191  Edge_collapsor<Complex> collapsor(complex_, num_collapses);
192  }
193 
194  void show_graph_stats() {
195  std::clog << "++++++ Graph stats +++++++" << std::endl;
196  std::clog << "Num vertices : " << complex_.num_vertices() << std::endl;
197  std::clog << "Num edges : " << complex_.num_edges() << std::endl;
198  std::clog << "Num connected components : " << complex_.num_connected_components() << std::endl;
199  std::clog << "Min/avg/max degree : " << min_degree() << "/" << avg_degree() << "/" << max_degree() << std::endl;
200  std::clog << "Num connected components : " << complex_.num_connected_components() << std::endl;
201  std::clog << "Num connected components : " << complex_.num_connected_components() << std::endl;
202  std::clog << "+++++++++++++++++++++++++" << std::endl;
203  }
204 
205  private:
206  int min_degree() const {
207  int res = (std::numeric_limits<int>::max)();
208  for (auto v : complex_.vertex_range())
209  res = (std::min)(res, complex_.degree(v));
210  return res;
211  }
212 
213  int max_degree() const {
214  int res = 0;
215  for (auto v : complex_.vertex_range())
216  res = (std::max)(res, complex_.degree(v));
217  return res;
218  }
219 
220  int avg_degree() const {
221  int res = 0;
222  for (auto v : complex_.vertex_range())
223  res += complex_.degree(v);
224  return res / complex_.num_vertices();
225  }
226 
227  public:
228  void show_complex_stats() {
229  std::clog << "++++++ Mesh stats +++++++" << std::endl;
230  std::clog << "Num vertices : " << complex_.num_vertices() << std::endl;
231  std::clog << "Num edges : " << complex_.num_edges() << std::endl;
232  std::clog << "Num connected components : " << complex_.num_connected_components() << std::endl;
233  std::clog << "+++++++++++++++++++++++++" << std::endl;
234  }
235 
236  void show_complex_dimension() {
237  unsigned num_simplices = 0;
238  int euler = 0;
239  int dimension = 0;
240  Gudhi::Clock clock;
241  for (const auto &s : complex_.complex_simplex_range()) {
242  num_simplices++;
243  dimension = (std::max)(s.dimension(), dimension);
244  if (s.dimension() % 2 == 0)
245  euler += 1;
246  else
247  euler -= 1;
248  }
249  clock.end();
250  std::clog << "++++++ Mesh dimension +++++++" << std::endl;
251  std::clog << "Dimension : " << dimension << std::endl;
252  std::clog << "Euler characteristic : " << euler << std::endl;
253  std::clog << "Num simplices : " << num_simplices << std::endl;
254  std::clog << "Total time: " << clock << std::endl;
255  std::clog << "Time per simplex: " << clock.num_seconds() / num_simplices << " s" << std::endl;
256  std::clog << "+++++++++++++++++++++++++" << std::endl;
257  }
258 
259  void show_homology_group() {
260 #ifdef _WIN32
261  std::clog << "Works only on linux x64 for the moment\n";
262 #else
263  Gudhi::Clock clock;
264  run_chomp();
265  clock.end();
266 #endif
267  }
268 
269  void show_euler_characteristic() {
270  unsigned num_simplices = 0;
271  int euler = 0;
272  int dimension = 0;
273  for (const auto &s : complex_.complex_simplex_range()) {
274  num_simplices++;
275  dimension = (std::max)(s.dimension(), dimension);
276  if (s.dimension() % 2 == 0)
277  euler += 1;
278  else
279  euler -= 1;
280  }
281  std::clog << "Saw " << num_simplices << " simplices with maximum dimension " << dimension << std::endl;
282  std::clog << "The euler characteristic is : " << euler << std::endl;
283  }
284 
285  void show_persistence(int p, double threshold, int max_dim, double min_pers) {
286  Persistence_compute<Complex> persistence(complex_, std::clog, Persistence_params(p, threshold, max_dim, min_pers));
287  }
288 
289  void show_critical_points(double max_distance) {
290  Critical_points<Complex> critical_points(complex_, std::clog, max_distance);
291  }
292 
293  void show_is_manifold() {
294  unsigned dim;
295  bool is_manifold;
296  Is_manifold<Complex> test_manifold(complex_, dim, is_manifold);
297 
298  if (is_manifold) {
299  std::clog << "The complex is a " << dim << "-manifold\n";
300  } else {
301  if (dim < 4) {
302  std::clog << "The complex has dimension greater than " << dim << " and is not a manifold\n";
303  } else {
304  std::clog << "The complex has dimension>=4 and may or may not be a manifold\n";
305  }
306  }
307  }
308 
309  private:
310  void run_chomp() {
311  save_complex_in_file_for_chomp();
312  std::clog << "Call CHOMP library\n";
313  int returnValue = system("homsimpl chomp.sim");
314  if (returnValue != 0) {
315  std::cerr << "homsimpl (from CHOMP) failed. Please check it is installed or available in the PATH."
316  << std::endl;
317  }
318  }
319 
320  void save_complex_in_file_for_chomp() {
321  std::ofstream file;
322  file.open("chomp.sim");
323  for (const auto &s : complex_.complex_simplex_range()) {
324  bool first = true;
325  file << "(";
326  for (auto x : s) {
327  if (first)
328  first = false;
329  else
330  file << ",";
331  file << x;
332  }
333  file << ")\n";
334  }
335  }
336 
337  public:
338  unsigned num_vertices() const {
339  return complex_.num_vertices();
340  }
341 
342  unsigned num_edges() const {
343  return complex_.num_edges();
344  }
345 };
346 
347 #endif // MODEL_MODEL_H_
Definition: Critical_points.h:26
Definition: Edge_collapsor.h:21
Definition: Edge_contractor.h:24
Edge_handle add_edge_without_blockers(Vertex_handle a, Vertex_handle b)
Adds an edge between vertices a and b without blockers.
Definition: Skeleton_blocker_complex.h:566
void keep_only_vertices()
The complex is reduced to its set of vertices. All the edges and blockers are removed.
Definition: Skeleton_blocker_complex.h:623
Complex_edge_range edge_range() const
Returns a Complex_edge_range over all edges of the simplicial complex.
Definition: Skeleton_blocker_complex.h:1318
Vertex_handle first_vertex(Edge_handle edge_handle) const
returns the first vertex of an edge
Definition: Skeleton_blocker_complex.h:512
int degree(Vertex_handle local) const
return the graph degree of a vertex.
Definition: Skeleton_blocker_complex.h:468
Vertex_handle second_vertex(Edge_handle edge_handle) const
returns the first vertex of an edge
Definition: Skeleton_blocker_complex.h:520
int num_connected_components() const
returns the number of connected components in the graph of the 1-skeleton.
Definition: Skeleton_blocker_complex.h:1029
virtual void clear()
Definition: Skeleton_blocker_complex.h:311
Complex_vertex_range vertex_range() const
Returns a Complex_vertex_range over all vertices of the complex.
Definition: Skeleton_blocker_complex.h:1284
Complex_simplex_range complex_simplex_range() const
Returns a Complex_simplex_range over all the simplices of the complex.
Definition: Skeleton_blocker_complex.h:1431
Vertex_handle add_vertex()
Add a vertex to the complex with a default constructed associated point.
Definition: Skeleton_blocker_geometric_complex.h:85
const Point & point(Vertex_handle v) const
Returns the Point associated to the vertex v.
Definition: Skeleton_blocker_geometric_complex.h:101
Definition: Is_manifold.h:24
Definition: Lloyd_builder.h:20
Definition: Persistence_compute.h:36
Definition: Vertex_collapsor.h:23
Handle type for the vertices of a cell complex.
Definition: VertexHandle.h:15
GUDHIdev  Version 3.5.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Wed Apr 6 2022 19:26:28 for GUDHIdev by Doxygen 1.9.1