GIC.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: Mathieu Carriere
4  *
5  * Copyright (C) 2017 Inria
6  *
7  * Modification(s):
8  * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
9  * - YYYY/MM Author: Description of the modification
10  */
11 
12 #ifndef GIC_H_
13 #define GIC_H_
14 
15 #ifdef GUDHI_USE_TBB
16 #include <tbb/parallel_for.h>
17 #include <mutex>
18 #endif
19 
20 #if __has_include(<CGAL/version.h>)
21 # define GUDHI_GIC_USE_CGAL 1
22 # include <gudhi/Bottleneck.h>
23 #elif __has_include(<hera/bottleneck.h>)
24 # define GUDHI_GIC_USE_HERA 1
25 # ifdef _MSC_VER
26 // https://github.com/grey-narn/hera/issues/3
27 // ssize_t is a non-standard type (well, posix)
28 # include <type_traits>
29 # include <cstdlib>
30 using ssize_t = std::make_signed_t<std::size_t>;
31 # endif
32 # include <hera/bottleneck.h>
33 #endif
34 
35 #include <gudhi/Debug_utils.h>
37 #include <gudhi/reader_utils.h>
38 #include <gudhi/Simplex_tree.h>
39 #include <gudhi/Rips_complex.h>
40 #include <gudhi/Points_off_io.h>
42 #include <gudhi/Persistent_cohomology.h>
43 
44 #include <boost/config.hpp>
45 #include <boost/graph/graph_traits.hpp>
46 #include <boost/graph/adjacency_list.hpp>
47 #include <boost/graph/connected_components.hpp>
48 #include <boost/graph/dijkstra_shortest_paths.hpp>
49 #include <boost/graph/subgraph.hpp>
50 #include <boost/graph/graph_utility.hpp>
51 
52 #include <iostream>
53 #include <vector>
54 #include <map>
55 #include <string>
56 #include <limits> // for numeric_limits
57 #include <utility> // for std::pair<>
58 #include <algorithm> // for (std::max)
59 #include <random>
60 #include <cassert>
61 #include <cmath>
62 
63 namespace Gudhi {
64 
65 namespace cover_complex {
66 
70 using Persistence_diagram = std::vector<std::pair<double, double> >;
71 using Graph = boost::subgraph<
72  boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
73  boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
74 using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
75 using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
76 using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
77 
98 template <typename Point>
100  private:
101  bool verbose = false; // whether to display information.
102  std::string type; // Nerve or GIC
103 
104  std::vector<Point> point_cloud; // input point cloud.
105  std::vector<std::vector<double> > distances; // all pairwise distances.
106  int maximal_dim; // maximal dimension of output simplicial complex.
107  int data_dimension; // dimension of input data.
108  int num_points; // number of points.
109 
110  std::vector<double> func; // function used to compute the output simplicial complex.
111  std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex.
112  bool functional_cover = false; // whether we use a cover with preimages of a function or not.
113 
114  Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists).
115  Graph one_skeleton; // one-skeleton used to compute the connected components.
116  std::vector<Vertex_t> vertices; // vertices of one_skeleton.
117 
118  std::vector<std::vector<int> > simplices; // simplices of output simplicial complex.
119  std::vector<int> voronoi_subsamples; // Voronoi germs (in case of Voronoi cover).
120 
121  Persistence_diagram PD;
122  std::vector<double> distribution;
123 
124  std::vector<std::vector<int> >
125  cover; // function associating to each data point the vector of cover elements to which it belongs.
126  std::map<int, std::vector<int> >
127  cover_back; // inverse of cover, in order to get the data points associated to a specific cover element.
128  std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence
129  // diagram of the output simplicial complex.
130  std::map<int, int>
131  cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not.
132  std::map<int, std::pair<int, double> >
133  cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex.
134 
135  int resolution_int = -1;
136  double resolution_double = -1;
137  double gain = -1;
138  double rate_constant = 10; // Constant in the subsampling.
139  double rate_power = 0.001; // Power in the subsampling.
140  int mask = 0; // Ignore nodes containing less than mask points.
141 
142  std::map<int, int> name2id, name2idinv;
143 
144  std::string cover_name;
145  std::string point_cloud_name;
146  std::string color_name;
147 
148  // Remove all edges of a graph.
149  void remove_edges(Graph& G) {
150  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
151  for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
152  }
153 
154  // Find random number in [0,1].
155  double GetUniform() {
156  thread_local std::default_random_engine re;
157  std::uniform_real_distribution<double> Dist(0, 1);
158  return Dist(re);
159  }
160 
161  // Subsample points.
162  void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector<int>& samples) {
163  int t = 0;
164  int m = 0;
165  double u;
166  while (m < sampleSize) {
167  u = GetUniform();
168  if ((populationSize - t) * u >= sampleSize - m) {
169  t++;
170  } else {
171  samples[m] = t;
172  t++;
173  m++;
174  }
175  }
176  }
177 
178  // *******************************************************************************************************************
179  // Utils.
180  // *******************************************************************************************************************
181 
182  public:
188  void set_type(const std::string& t) { type = t; }
189 
190  public:
196  void set_verbose(bool verb = false) { verbose = verb; }
197 
198  public:
206  void set_subsampling(double constant, double power) {
207  rate_constant = constant;
208  rate_power = power;
209  }
210 
211  public:
219  void set_mask(int nodemask) { mask = nodemask; }
220 
221  public:
222 
223 
229  void set_point_cloud_from_range(const std::vector<std::vector<double> > & point_cloud) {
230  this->num_points = point_cloud.size(); data_dimension = point_cloud[0].size();
231  point_cloud_name = "cloud"; cover.resize(this->num_points);
232  for(int i = 0; i < this->num_points; i++){
233  boost::add_vertex(one_skeleton_OFF);
234  vertices.push_back(boost::add_vertex(one_skeleton));
235  }
236  this->point_cloud = point_cloud;
237  }
238 
244  bool read_point_cloud(const std::string& off_file_name) {
245  point_cloud_name = off_file_name;
246  std::ifstream input(off_file_name);
247  std::string line;
248 
249  char comment = '#';
250  while (comment == '#') {
251  std::getline(input, line);
252  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
253  comment = line[line.find_first_not_of(' ')];
254  }
255  if (strcmp((char*)line.c_str(), "nOFF") == 0) {
256  comment = '#';
257  while (comment == '#') {
258  std::getline(input, line);
259  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
260  comment = line[line.find_first_not_of(' ')];
261  }
262  std::stringstream stream(line);
263  stream >> data_dimension;
264  } else {
265  data_dimension = 3;
266  }
267 
268  comment = '#';
269  int numedges, numfaces, i, dim;
270  while (comment == '#') {
271  std::getline(input, line);
272  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
273  comment = line[line.find_first_not_of(' ')];
274  }
275  std::stringstream stream(line);
276  stream >> this->num_points;
277  stream >> numfaces;
278  stream >> numedges;
279 
280  i = 0;
281  while (i < this->num_points) {
282  std::getline(input, line);
283  if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
284  !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
285  std::stringstream iss(line);
286  std::vector<double> point;
287  point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
288  point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
289  boost::add_vertex(one_skeleton_OFF);
290  vertices.push_back(boost::add_vertex(one_skeleton));
291  cover.emplace_back();
292  i++;
293  }
294  }
295 
296  i = 0;
297  while (i < numfaces) {
298  std::getline(input, line);
299  if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
300  !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
301  std::vector<int> simplex;
302  std::stringstream iss(line);
303  simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>());
304  dim = simplex[0];
305  for (int j = 1; j <= dim; j++)
306  for (int k = j + 1; k <= dim; k++)
307  boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF);
308  i++;
309  }
310  }
311 
312  return input.is_open();
313  }
314 
315  // *******************************************************************************************************************
316  // Graphs.
317  // *******************************************************************************************************************
318 
319  public: // Set graph from file.
327  void set_graph_from_file(const std::string& graph_file_name) {
328  remove_edges(one_skeleton);
329  int neighb;
330  std::ifstream input(graph_file_name);
331  std::string line;
332  int source;
333  while (std::getline(input, line)) {
334  std::stringstream stream(line);
335  stream >> source;
336  while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
337  }
338  }
339 
340  public: // Set graph from OFF file.
345  remove_edges(one_skeleton);
346  if (num_edges(one_skeleton_OFF))
347  one_skeleton = one_skeleton_OFF;
348  else
349  std::cerr << "No triangulation read in OFF file!" << std::endl;
350  }
351 
352  public: // Set graph from Rips complex.
359  template <typename Distance>
360  void set_graph_from_rips(double threshold, Distance distance) {
361  remove_edges(one_skeleton);
362  if (distances.size() == 0) compute_pairwise_distances(distance);
363  for (int i = 0; i < this->num_points; i++) {
364  for (int j = i + 1; j < this->num_points; j++) {
365  if (distances[i][j] <= threshold) {
366  boost::add_edge(vertices[i], vertices[j], one_skeleton);
367  boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first,
368  distances[i][j]);
369  }
370  }
371  }
372  }
373 
374  public:
375  void set_graph_weights() {
376  Index_map index = boost::get(boost::vertex_index, one_skeleton);
377  Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
378  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
379  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
380  boost::put(weight, *ei,
381  distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]);
382  }
383 
384  public:
390  void set_distances_from_range(const std::vector<std::vector<double> > & distance_matrix) {
391  this->num_points = distance_matrix.size(); data_dimension = 0; point_cloud_name = "matrix";
392  cover.resize(this->num_points); point_cloud.resize(this->num_points);
393  for(int i = 0; i < this->num_points; i++){
394  boost::add_vertex(one_skeleton_OFF);
395  vertices.push_back(boost::add_vertex(one_skeleton));
396  }
397  distances = distance_matrix;
398  }
399 
400  public: // Pairwise distances.
403  template <typename Distance>
404  void compute_pairwise_distances(Distance ref_distance) {
405  std::vector<double> zeros(this->num_points);
406  for (int i = 0; i < this->num_points; i++) distances.push_back(zeros);
407  if (verbose) std::clog << "Computing distances..." << std::endl;
408  for (int i = 0; i < this->num_points; i++) {
409  int state = 100 * (i + 1) / this->num_points;
410  if (verbose && state % 10 == 0) std::clog << "\r" << state << "%" << std::flush;
411  for (int j = i; j < this->num_points; j++) {
412  double dis = ref_distance(point_cloud[i], point_cloud[j]);
413  distances[i][j] = dis;
414  distances[j][i] = dis;
415  }
416  }
417  if (verbose) std::clog << std::endl;
418  }
419 
420  public: // Automatic tuning of Rips complex.
430  template <typename Distance>
431  double set_graph_from_automatic_rips(Distance distance, int N = 100) {
432  int m = floor(this->num_points / std::exp((1 + rate_power) * std::log(std::log(this->num_points) / std::log(rate_constant))));
433  m = (std::min)(m, this->num_points - 1);
434  double delta = 0;
435 
436  if (verbose) std::clog << this->num_points << " points in R^" << data_dimension << std::endl;
437  if (verbose) std::clog << "Subsampling " << m << " points" << std::endl;
438 
439  if (distances.size() == 0) compute_pairwise_distances(distance);
440 
441  #ifdef GUDHI_USE_TBB
442  std::mutex deltamutex;
443  tbb::parallel_for(0, N, [&](int i){
444  std::vector<int> samples(m);
445  SampleWithoutReplacement(this->num_points, m, samples);
446  double hausdorff_dist = 0;
447  for (int j = 0; j < this->num_points; j++) {
448  double mj = distances[j][samples[0]];
449  for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
450  hausdorff_dist = (std::max)(hausdorff_dist, mj);
451  }
452  deltamutex.lock();
453  delta += hausdorff_dist / N;
454  deltamutex.unlock();
455  });
456  #else
457  for (int i = 0; i < N; i++) {
458  std::vector<int> samples(m);
459  SampleWithoutReplacement(this->num_points, m, samples);
460  double hausdorff_dist = 0;
461  for (int j = 0; j < this->num_points; j++) {
462  double mj = distances[j][samples[0]];
463  for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
464  hausdorff_dist = (std::max)(hausdorff_dist, mj);
465  }
466  delta += hausdorff_dist / N;
467  }
468  #endif
469 
470  if (verbose) std::clog << "delta = " << delta << std::endl;
471  set_graph_from_rips(delta, distance);
472  return delta;
473  }
474 
475  // *******************************************************************************************************************
476  // Functions.
477  // *******************************************************************************************************************
478 
479  public: // Set function from file.
485  void set_function_from_file(const std::string& func_file_name) {
486  std::ifstream input(func_file_name);
487  std::string line;
488  double f;
489  while (std::getline(input, line)) {
490  std::stringstream stream(line);
491  stream >> f;
492  func.push_back(f);
493  }
494  functional_cover = true;
495  cover_name = func_file_name;
496  }
497 
498  public: // Set function from kth coordinate
505  if(point_cloud[0].size() > 0){
506  for (int i = 0; i < this->num_points; i++) func.push_back(point_cloud[i][k]);
507  functional_cover = true;
508  cover_name = "coordinate " + std::to_string(k);
509  }
510  else{
511  std::cerr << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl;
512  for (int i = 0; i < this->num_points; i++) func.push_back(0.0);
513  functional_cover = true;
514  cover_name = "null";
515  }
516  }
517 
518  public: // Set function from vector.
524  template <class InputRange>
525  void set_function_from_range(InputRange const& function) {
526  for (int i = 0; i < this->num_points; i++) func.push_back(function[i]);
527  functional_cover = true;
528  }
529 
530  // *******************************************************************************************************************
531  // Covers.
532  // *******************************************************************************************************************
533 
534  public: // Automatic tuning of resolution.
543  if (!functional_cover) {
544  std::cerr << "Cover needs to come from the preimages of a function." << std::endl;
545  return 0;
546  }
547  if (type != "Nerve" && type != "GIC") {
548  std::cerr << "Type of complex needs to be specified." << std::endl;
549  return 0;
550  }
551 
552  double reso = 0;
553  Index_map index = boost::get(boost::vertex_index, one_skeleton);
554 
555  if (type == "GIC") {
556  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
557  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
558  reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
559  func[index[boost::target(*ei, one_skeleton)]]));
560  if (verbose) std::clog << "resolution = " << reso << std::endl;
561  resolution_double = reso;
562  }
563 
564  if (type == "Nerve") {
565  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
566  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
567  reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
568  func[index[boost::target(*ei, one_skeleton)]]) /
569  gain);
570  if (verbose) std::clog << "resolution = " << reso << std::endl;
571  resolution_double = reso;
572  }
573 
574  return reso;
575  }
576 
577  public:
583  void set_resolution_with_interval_length(double reso) { resolution_double = reso; }
589  void set_resolution_with_interval_number(int reso) { resolution_int = reso; }
595  void set_gain(double g = 0.3) { gain = g; }
596 
597  public: // Set cover with preimages of function.
602  if (resolution_double == -1 && resolution_int == -1) {
603  std::cerr << "Number and/or length of intervals not specified" << std::endl;
604  return;
605  }
606  if (gain == -1) {
607  std::cerr << "Gain not specified" << std::endl;
608  return;
609  }
610 
611  // Read function values and compute min and max
612  double minf = (std::numeric_limits<float>::max)();
613  double maxf = std::numeric_limits<float>::lowest();
614  for (int i = 0; i < this->num_points; i++) {
615  minf = (std::min)(minf, func[i]);
616  maxf = (std::max)(maxf, func[i]);
617  }
618  if (verbose) std::clog << "Min function value = " << minf << " and Max function value = " << maxf << std::endl;
619 
620  // Compute cover of im(f)
621  std::vector<std::pair<double, double> > intervals;
622  int res;
623 
624  if (resolution_double == -1) { // Case we use an integer for the number of intervals.
625  double incr = (maxf - minf) / resolution_int;
626  double x = minf;
627  double alpha = (incr * gain) / (2 - 2 * gain);
628  double y = minf + incr + alpha;
629  std::pair<double, double> interm(x, y);
630  intervals.push_back(interm);
631  for (int i = 1; i < resolution_int - 1; i++) {
632  x = minf + i * incr - alpha;
633  y = minf + (i + 1) * incr + alpha;
634  std::pair<double, double> inter(x, y);
635  intervals.push_back(inter);
636  }
637  x = minf + (resolution_int - 1) * incr - alpha;
638  y = maxf;
639  std::pair<double, double> interM(x, y);
640  intervals.push_back(interM);
641  res = intervals.size();
642  if (verbose) {
643  for (int i = 0; i < res; i++)
644  std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
645  << std::endl;
646  }
647  } else {
648  if (resolution_int == -1) { // Case we use a double for the length of the intervals.
649  double x = minf;
650  double y = x + resolution_double;
651  while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
652  std::pair<double, double> inter(x, y);
653  intervals.push_back(inter);
654  x = y - gain * resolution_double;
655  y = x + resolution_double;
656  }
657  std::pair<double, double> interM(x, maxf);
658  intervals.push_back(interM);
659  res = intervals.size();
660  if (verbose) {
661  for (int i = 0; i < res; i++)
662  std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
663  << std::endl;
664  }
665  } else { // Case we use an integer and a double for the length of the intervals.
666  double x = minf;
667  double y = x + resolution_double;
668  int count = 0;
669  while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
670  std::pair<double, double> inter(x, y);
671  intervals.push_back(inter);
672  count++;
673  x = y - gain * resolution_double;
674  y = x + resolution_double;
675  }
676  res = intervals.size();
677  if (verbose) {
678  for (int i = 0; i < res; i++)
679  std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
680  << std::endl;
681  }
682  }
683  }
684 
685  // Sort points according to function values
686  std::vector<int> points(this->num_points);
687  for (int i = 0; i < this->num_points; i++) points[i] = i;
688  std::sort(points.begin(), points.end(), [this](int p1, int p2){return (this->func[p1] < this->func[p2]);});
689 
690  int id = 0;
691  int pos = 0;
692  Index_map index = boost::get(boost::vertex_index, one_skeleton); // int maxc = -1;
693  std::map<int, std::vector<int> > preimages;
694  std::map<int, double> funcstd;
695 
696  if (verbose) std::clog << "Computing preimages..." << std::endl;
697  for (int i = 0; i < res; i++) {
698  // Find points in the preimage
699  std::pair<double, double> inter1 = intervals[i];
700  int tmp = pos;
701  double u, v;
702 
703  if (i != res - 1) {
704  if (i != 0) {
705  std::pair<double, double> inter3 = intervals[i - 1];
706  while (func[points[tmp]] < inter3.second && tmp != this->num_points) {
707  preimages[i].push_back(points[tmp]);
708  tmp++;
709  }
710  u = inter3.second;
711  } else {
712  u = inter1.first;
713  }
714 
715  std::pair<double, double> inter2 = intervals[i + 1];
716  while (func[points[tmp]] < inter2.first && tmp != this->num_points) {
717  preimages[i].push_back(points[tmp]);
718  tmp++;
719  }
720  v = inter2.first;
721  pos = tmp;
722  while (func[points[tmp]] < inter1.second && tmp != this->num_points) {
723  preimages[i].push_back(points[tmp]);
724  tmp++;
725  }
726 
727  } else {
728  std::pair<double, double> inter3 = intervals[i - 1];
729  while (func[points[tmp]] < inter3.second && tmp != this->num_points) {
730  preimages[i].push_back(points[tmp]);
731  tmp++;
732  }
733  while (tmp != this->num_points) {
734  preimages[i].push_back(points[tmp]);
735  tmp++;
736  }
737  u = inter3.second;
738  v = inter1.second;
739  }
740 
741  funcstd[i] = 0.5 * (u + v);
742  }
743 
744  #ifdef GUDHI_USE_TBB
745  if (verbose) std::clog << "Computing connected components (parallelized)..." << std::endl;
746  std::mutex covermutex, idmutex;
747  tbb::parallel_for(0, res, [&](int i){
748  // Compute connected components
749  Graph G = one_skeleton.create_subgraph();
750  int num = preimages[i].size();
751  std::vector<int> component(num);
752  for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
753  boost::connected_components(G, &component[0]);
754  int max = 0;
755 
756  // For each point in preimage
757  for (int j = 0; j < num; j++) {
758  // Update number of components in preimage
759  if (component[j] > max) max = component[j];
760 
761  // Identify component with Cantor polynomial N^2 -> N
762  int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
763 
764  // Update covers
765  covermutex.lock();
766  cover[preimages[i][j]].push_back(identifier);
767  cover_back[identifier].push_back(preimages[i][j]);
768  cover_fct[identifier] = i;
769  cover_std[identifier] = funcstd[i];
770  cover_color[identifier].second += func_color[preimages[i][j]];
771  cover_color[identifier].first += 1;
772  covermutex.unlock();
773  }
774 
775  // Maximal dimension is total number of connected components
776  idmutex.lock();
777  id += max + 1;
778  idmutex.unlock();
779  });
780  #else
781  if (verbose) std::clog << "Computing connected components..." << std::endl;
782  for (int i = 0; i < res; i++) {
783  // Compute connected components
784  Graph G = one_skeleton.create_subgraph();
785  int num = preimages[i].size();
786  std::vector<int> component(num);
787  for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
788  boost::connected_components(G, &component[0]);
789  int max = 0;
790 
791  // For each point in preimage
792  for (int j = 0; j < num; j++) {
793  // Update number of components in preimage
794  if (component[j] > max) max = component[j];
795 
796  // Identify component with Cantor polynomial N^2 -> N
797  int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
798 
799  // Update covers
800  cover[preimages[i][j]].push_back(identifier);
801  cover_back[identifier].push_back(preimages[i][j]);
802  cover_fct[identifier] = i;
803  cover_std[identifier] = funcstd[i];
804  cover_color[identifier].second += func_color[preimages[i][j]];
805  cover_color[identifier].first += 1;
806  }
807 
808  // Maximal dimension is total number of connected components
809  id += max + 1;
810  }
811  #endif
812 
813  maximal_dim = id - 1;
814  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
815  iit->second.second /= iit->second.first;
816  }
817 
818  public: // Set cover from file.
825  void set_cover_from_file(const std::string& cover_file_name) {
826  int i = 0;
827  int cov;
828  std::vector<int> cov_elts, cov_number;
829  std::ifstream input(cover_file_name);
830  std::string line;
831  while (std::getline(input, line)) {
832  cov_elts.clear();
833  std::stringstream stream(line);
834  while (stream >> cov) {
835  cov_elts.push_back(cov);
836  cov_number.push_back(cov);
837  cover_fct[cov] = cov;
838  cover_color[cov].second += func_color[i];
839  cover_color[cov].first++;
840  cover_back[cov].push_back(i);
841  }
842  cover[i] = cov_elts;
843  i++;
844  }
845 
846  std::sort(cov_number.begin(), cov_number.end());
847  std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
848  cov_number.resize(std::distance(cov_number.begin(), it));
849 
850  maximal_dim = cov_number.size() - 1;
851  for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
852  cover_name = cover_file_name;
853  }
854 
855  public: // Set cover from range.
862  template <class AssignmentRange>
863  void set_cover_from_range(AssignmentRange const& assignments) {
864  std::vector<int> cov_elts, cov_number;
865  for(int i=0; i < static_cast<int>(assignments.size()); i++){
866  cov_elts.clear();
867  for (int cov : assignments[i]){
868  cov_elts.push_back(cov);
869  cov_number.push_back(cov);
870  cover_fct[cov] = cov;
871  auto& cc = cover_color[cov];
872  cc.second += func_color[i];
873  cc.first++;
874  cover_back[cov].push_back(i);
875  }
876  cover[i] = cov_elts;
877  }
878 
879  std::sort(cov_number.begin(), cov_number.end());
880  std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
881  cov_number.resize(std::distance(cov_number.begin(), it));
882 
883  maximal_dim = cov_number.size() - 1;
884  for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
885  }
886 
887  public: // Set cover from Voronoi
894  template <typename Distance>
895  void set_cover_from_Voronoi(Distance distance, int m = 100) {
896  voronoi_subsamples.resize(m);
897  SampleWithoutReplacement(this->num_points, m, voronoi_subsamples);
898  if (distances.size() == 0) compute_pairwise_distances(distance);
899  set_graph_weights();
900  Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
901  Index_map index = boost::get(boost::vertex_index, one_skeleton);
902  std::vector<double> mindist(this->num_points);
903  for (int j = 0; j < this->num_points; j++) mindist[j] = (std::numeric_limits<double>::max)();
904 
905  // Compute the geodesic distances to subsamples with Dijkstra
906  #ifdef GUDHI_USE_TBB
907  if (verbose) std::clog << "Computing geodesic distances (parallelized)..." << std::endl;
908  std::mutex coverMutex; std::mutex mindistMutex;
909  tbb::parallel_for(0, m, [&](int i){
910  int seed = voronoi_subsamples[i];
911  std::vector<double> dmap(this->num_points);
912  boost::dijkstra_shortest_paths(
913  one_skeleton, vertices[seed],
914  boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
915 
916  coverMutex.lock(); mindistMutex.lock();
917  for (int j = 0; j < this->num_points; j++)
918  if (mindist[j] > dmap[j]) {
919  mindist[j] = dmap[j];
920  if (cover[j].size() == 0)
921  cover[j].push_back(i);
922  else
923  cover[j][0] = i;
924  }
925  coverMutex.unlock(); mindistMutex.unlock();
926  });
927  #else
928  for (int i = 0; i < m; i++) {
929  if (verbose) std::clog << "Computing geodesic distances to seed " << i << "..." << std::endl;
930  int seed = voronoi_subsamples[i];
931  std::vector<double> dmap(this->num_points);
932  boost::dijkstra_shortest_paths(
933  one_skeleton, vertices[seed],
934  boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
935 
936  for (int j = 0; j < this->num_points; j++)
937  if (mindist[j] > dmap[j]) {
938  mindist[j] = dmap[j];
939  if (cover[j].size() == 0)
940  cover[j].push_back(i);
941  else
942  cover[j][0] = i;
943  }
944  }
945  #endif
946 
947  for (int i = 0; i < this->num_points; i++) {
948  cover_back[cover[i][0]].push_back(i);
949  cover_color[cover[i][0]].second += func_color[i];
950  cover_color[cover[i][0]].first++;
951  }
952  for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first;
953  maximal_dim = m - 1;
954  cover_name = "Voronoi";
955  }
956 
957  public: // return subset of data corresponding to a node
964  const std::vector<int>& subpopulation(int c) { return cover_back[c]; }
965 
966 
967  public:
974  double subcolor(int c) { return cover_color[c].second; }
975 
976  // *******************************************************************************************************************
977  // Visualization.
978  // *******************************************************************************************************************
979 
980  public: // Set color from file.
987  void set_color_from_file(const std::string& color_file_name) {
988  std::ifstream input(color_file_name);
989  std::string line;
990  double f;
991  while (std::getline(input, line)) {
992  std::stringstream stream(line);
993  stream >> f;
994  func_color.push_back(f);
995  }
996  color_name = color_file_name;
997  }
998 
999  public: // Set color from kth coordinate
1005  void set_color_from_coordinate(int k = 0) {
1006  if(point_cloud[0].size() > 0){
1007  for (int i = 0; i < this->num_points; i++) func_color.push_back(point_cloud[i][k]);
1008  color_name = "coordinate ";
1009  color_name.append(std::to_string(k));
1010  }
1011  else{
1012  std::cerr << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl;
1013  for (int i = 0; i < this->num_points; i++) func.push_back(0.0);
1014  functional_cover = true;
1015  cover_name = "null";
1016  }
1017  }
1018 
1019  public: // Set color from vector.
1025  void set_color_from_range(std::vector<double> color) {
1026  for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
1027  }
1028 
1029  public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
1034  void plot_DOT() {
1035  std::string mapp = point_cloud_name + "_sc.dot";
1036  std::ofstream graphic(mapp);
1037 
1038  double maxv = std::numeric_limits<double>::lowest();
1039  double minv = (std::numeric_limits<double>::max)();
1040  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1041  maxv = (std::max)(maxv, iit->second.second);
1042  minv = (std::min)(minv, iit->second.second);
1043  }
1044 
1045  std::vector<int> nodes;
1046  nodes.clear();
1047 
1048  graphic << "graph GIC {" << std::endl;
1049  int id = 0;
1050  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1051  if (iit->second.first > mask) {
1052  nodes.push_back(iit->first);
1053  name2id[iit->first] = id;
1054  name2idinv[id] = iit->first;
1055  id++;
1056  graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
1057  << ":" << iit->second.first << "\" style=filled fillcolor=\""
1058  << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl;
1059  }
1060  }
1061  int num_simplices = simplices.size();
1062  for (int i = 0; i < num_simplices; i++)
1063  if (simplices[i].size() == 2) {
1064  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
1065  graphic << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];"
1066  << std::endl;
1067  }
1068  }
1069  graphic << "}";
1070  graphic.close();
1071  std::clog << mapp << " file generated. It can be visualized with e.g. neato." << std::endl;
1072  }
1073 
1074  public: // Create a .txt file that can be compiled with KeplerMapper.
1078  void write_info() {
1079  int num_simplices = simplices.size();
1080  int num_edges = 0;
1081  std::string mapp = point_cloud_name + "_sc.txt";
1082  std::ofstream graphic(mapp);
1083 
1084  for (int i = 0; i < num_simplices; i++)
1085  if (simplices[i].size() == 2)
1086  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++;
1087 
1088  graphic << point_cloud_name << std::endl;
1089  graphic << cover_name << std::endl;
1090  graphic << color_name << std::endl;
1091  graphic << resolution_double << " " << gain << std::endl;
1092  graphic << cover_color.size() << " " << num_edges << std::endl;
1093 
1094  int id = 0;
1095  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1096  graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl;
1097  name2id[iit->first] = id;
1098  name2idinv[id] = iit->first;
1099  id++;
1100  }
1101 
1102  for (int i = 0; i < num_simplices; i++)
1103  if (simplices[i].size() == 2)
1104  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
1105  graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl;
1106  graphic.close();
1107  std::clog << mapp
1108  << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
1109  << std::endl;
1110  }
1111 
1112  public: // Create a .off file that can be visualized (e.g. with Geomview).
1117  void plot_OFF() {
1118  assert(cover_name == "Voronoi");
1119 
1120  int m = voronoi_subsamples.size();
1121  int numedges = 0;
1122  int numfaces = 0;
1123  std::vector<std::vector<int> > edges, faces;
1124  int numsimplices = simplices.size();
1125 
1126  std::string mapp = point_cloud_name + "_sc.off";
1127  std::ofstream graphic(mapp);
1128 
1129  graphic << "OFF" << std::endl;
1130  for (int i = 0; i < numsimplices; i++) {
1131  if (simplices[i].size() == 2) {
1132  numedges++;
1133  edges.push_back(simplices[i]);
1134  }
1135  if (simplices[i].size() == 3) {
1136  numfaces++;
1137  faces.push_back(simplices[i]);
1138  }
1139  }
1140  graphic << m << " " << numedges + numfaces << std::endl;
1141  for (int i = 0; i < m; i++) {
1142  if (data_dimension <= 3) {
1143  for (int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1144  for (int j = data_dimension; j < 3; j++) graphic << 0 << " ";
1145  graphic << std::endl;
1146  } else {
1147  for (int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1148  }
1149  }
1150  for (int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl;
1151  for (int i = 0; i < numfaces; i++)
1152  graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
1153  graphic.close();
1154  std::clog << mapp << " generated. It can be visualized with e.g. geomview." << std::endl;
1155  }
1156 
1157  // *******************************************************************************************************************
1158  // Extended Persistence Diagrams.
1159  // *******************************************************************************************************************
1160 
1161  public:
1165  Persistence_diagram compute_PD() {
1166  Simplex_tree st;
1167 
1168  // Compute max and min
1169  double maxf = std::numeric_limits<double>::lowest();
1170  double minf = (std::numeric_limits<double>::max)();
1171  for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1172  maxf = (std::max)(maxf, it->second);
1173  minf = (std::min)(minf, it->second);
1174  }
1175 
1176  // Build filtration
1177  for (auto const& simplex : simplices) {
1178  std::vector<int> splx = simplex;
1179  splx.push_back(-2);
1180  st.insert_simplex_and_subfaces(splx, -3);
1181  }
1182 
1183  for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1184  int vertex = it->first; float val = it->second;
1185  int vert[] = {vertex}; int edge[] = {vertex, -2};
1186  if(st.find(vert) != st.null_simplex()){
1187  st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf));
1188  st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf));
1189  }
1190  }
1192 
1193  // Compute PD
1196 
1197  // Output PD
1198  int max_dim = st.dimension();
1199  for (int i = 0; i < max_dim; i++) {
1200  std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i);
1201  int num_bars = bars.size(); if(i == 0) num_bars -= 1;
1202  if(verbose) std::clog << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
1203  for (int j = 0; j < num_bars; j++) {
1204  double birth = bars[j].first;
1205  double death = bars[j].second;
1206  if (i == 0 && std::isinf(death)) continue;
1207  if (birth < 0)
1208  birth = minf + (birth + 2) * (maxf - minf);
1209  else
1210  birth = minf + (2 - birth) * (maxf - minf);
1211  if (death < 0)
1212  death = minf + (death + 2) * (maxf - minf);
1213  else
1214  death = minf + (2 - death) * (maxf - minf);
1215  PD.push_back(std::pair<double, double>(birth, death));
1216  if (verbose) std::clog << " [" << birth << ", " << death << "]" << std::endl;
1217  }
1218  }
1219  return PD;
1220  }
1221 
1222  public:
1228  void compute_distribution(unsigned int N = 100) {
1229  unsigned int sz = distribution.size();
1230  if (sz < N) {
1231  for (unsigned int i = 0; i < N - sz; i++) {
1232  if (verbose) std::clog << "Computing " << i << "th bootstrap, bottleneck distance = ";
1233 
1234  Cover_complex Cboot; Cboot.num_points = this->num_points; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true;
1235 
1236  std::vector<int> boot(this->num_points);
1237  for (int j = 0; j < this->num_points; j++) {
1238  double u = GetUniform();
1239  int id = std::floor(u * (this->num_points)); boot[j] = id;
1240  Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]);
1241  boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
1242  }
1243  Cboot.set_color_from_range(Cboot.func);
1244 
1245  for (int j = 0; j < this->num_points; j++) {
1246  std::vector<double> dist(this->num_points);
1247  for (int k = 0; k < this->num_points; k++) dist[k] = distances[boot[j]][boot[k]];
1248  Cboot.distances.push_back(dist);
1249  }
1250 
1252  Cboot.set_gain();
1253  Cboot.set_automatic_resolution();
1254  Cboot.set_cover_from_function();
1255  Cboot.find_simplices();
1256  Cboot.compute_PD();
1257 #ifdef GUDHI_GIC_USE_CGAL
1258  double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD);
1259 #elif defined GUDHI_GIC_USE_HERA
1260  double db = hera::bottleneckDistExact(this->PD, Cboot.PD);
1261 #else
1262  double db;
1263  throw std::logic_error("This function requires CGAL or Hera for the bottleneck distance.");
1264 #endif
1265  if (verbose) std::clog << db << std::endl;
1266  distribution.push_back(db);
1267  }
1268 
1269  std::sort(distribution.begin(), distribution.end());
1270  }
1271  }
1272 
1273  public:
1280  unsigned int N = distribution.size();
1281  double d = distribution[std::floor(alpha * N)];
1282  if (verbose) std::clog << "Distance corresponding to confidence " << alpha << " is " << d << std::endl;
1283  return d;
1284  }
1285 
1286  public:
1293  unsigned int N = distribution.size();
1294  double level = 1;
1295  for (unsigned int i = 0; i < N; i++)
1296  if (distribution[i] >= d){ level = i * 1.0 / N; break; }
1297  if (verbose) std::clog << "Confidence level of distance " << d << " is " << level << std::endl;
1298  return level;
1299  }
1300 
1301  public:
1306  double compute_p_value() {
1307  double distancemin = (std::numeric_limits<double>::max)(); int N = PD.size();
1308  for (int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first));
1309  double p_value = 1 - compute_confidence_level_from_distance(distancemin);
1310  if (verbose) std::clog << "p value = " << p_value << std::endl;
1311  return p_value;
1312  }
1313 
1314  // *******************************************************************************************************************
1315  // Computation of simplices.
1316  // *******************************************************************************************************************
1317 
1318  public:
1324  template <typename SimplicialComplex>
1325  void create_complex(SimplicialComplex& complex) {
1326  unsigned int dimension = 0;
1327  for (auto const& simplex : simplices) {
1328  int numvert = simplex.size();
1329  double filt = std::numeric_limits<double>::lowest();
1330  for (int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt);
1331  complex.insert_simplex_and_subfaces(simplex, filt);
1332  if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
1333  }
1334  }
1335 
1336  public:
1340  if (type != "Nerve" && type != "GIC") {
1341  std::cerr << "Type of complex needs to be specified." << std::endl;
1342  return;
1343  }
1344 
1345  if (type == "Nerve") {
1346  simplices = cover;
1347  std::sort(simplices.begin(), simplices.end());
1348  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1349  simplices.resize(std::distance(simplices.begin(), it));
1350  }
1351 
1352  if (type == "GIC") {
1353  Index_map index = boost::get(boost::vertex_index, one_skeleton);
1354 
1355  if (functional_cover) {
1356  // Computes the simplices in the GIC by looking at all the edges of the graph and adding the
1357  // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals.
1358 
1359  if (gain >= 0.5)
1360  throw std::invalid_argument(
1361  "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
1362 
1363  // Loop on all edges.
1364  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1365  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
1366  int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
1367  for (int i = 0; i < nums; i++) {
1368  int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
1369  int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
1370  for (int j = 0; j < numt; j++) {
1371  int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
1372  if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
1373  std::vector<int> edge(2);
1374  edge[0] = (std::min)(vs, vt);
1375  edge[1] = (std::max)(vs, vt);
1376  simplices.push_back(edge);
1377  goto afterLoop;
1378  }
1379  }
1380  }
1381  afterLoop:;
1382  }
1383  std::sort(simplices.begin(), simplices.end());
1384  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1385  simplices.resize(std::distance(simplices.begin(), it));
1386 
1387  } else {
1388  // Find edges to keep
1389  Simplex_tree st;
1390  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1391  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
1392  if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
1393  cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
1394  std::vector<int> edge(2);
1395  edge[0] = index[boost::source(*ei, one_skeleton)];
1396  edge[1] = index[boost::target(*ei, one_skeleton)];
1397  st.insert_simplex_and_subfaces(edge);
1398  }
1399 
1400  // st.insert_graph(one_skeleton);
1401 
1402  // Build the Simplex Tree corresponding to the graph
1403  st.expansion(maximal_dim);
1404 
1405  // Find simplices of GIC
1406  simplices.clear();
1407  for (auto simplex : st.complex_simplex_range()) {
1408  if (!st.has_children(simplex)) {
1409  std::vector<int> simplx;
1410  for (auto vertex : st.simplex_vertex_range(simplex)) {
1411  unsigned int sz = cover[vertex].size();
1412  for (unsigned int i = 0; i < sz; i++) {
1413  simplx.push_back(cover[vertex][i]);
1414  }
1415  }
1416  std::sort(simplx.begin(), simplx.end());
1417  std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
1418  simplx.resize(std::distance(simplx.begin(), it));
1419  simplices.push_back(simplx);
1420  }
1421  }
1422  std::sort(simplices.begin(), simplices.end());
1423  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1424  simplices.resize(std::distance(simplices.begin(), it));
1425  }
1426  }
1427  }
1428 };
1429 
1430 } // namespace cover_complex
1431 
1432 } // namespace Gudhi
1433 
1434 #endif // GIC_H_
Compute the Euclidean distance between two Points given by a range of coordinates....
Definition: distance_functions.h:32
Options::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Simplex_tree.h:102
void assign_filtration(Simplex_handle sh, Filtration_value fv)
Sets the filtration value of a simplex.
Definition: Simplex_tree.h:625
std::pair< Simplex_handle, bool > insert_simplex_and_subfaces(const InputVertexRange &Nsimplex, Filtration_value filtration=0)
Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of Vertex_handles,...
Definition: Simplex_tree.h:942
bool make_filtration_non_decreasing()
This function ensures that each simplex has a higher filtration value than its faces by increasing th...
Definition: Simplex_tree.h:1895
Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const
Returns a range over the vertices of a simplex.
Definition: Simplex_tree.h:349
bool has_children(SimplexHandle sh) const
Returns true if the node in the simplex tree pointed by the given simplex handle has children.
Definition: Simplex_tree.h:742
Simplex_handle find(const InputVertexRange &s)
Given a range of Vertex_handles, returns the Simplex_handle of the simplex in the simplicial complex ...
Definition: Simplex_tree.h:767
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1360
static Simplex_handle null_simplex()
Returns a Simplex_handle different from all Simplex_handles associated to the simplices in the simpli...
Definition: Simplex_tree.h:635
Complex_simplex_range complex_simplex_range()
Returns a range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:304
Cover complex data structure.
Definition: GIC.h:99
void set_function_from_file(const std::string &func_file_name)
Creates the function f from a file containing the function values.
Definition: GIC.h:485
double set_automatic_resolution()
Computes the optimal length of intervals (i.e. the smallest interval length avoiding discretization a...
Definition: GIC.h:542
void set_cover_from_Voronoi(Distance distance, int m=100)
Creates the cover C from the Voronoï cells of a subsampling of the point cloud.
Definition: GIC.h:895
void set_resolution_with_interval_number(int reso)
Sets a number of intervals from a value stored in memory.
Definition: GIC.h:589
void set_mask(int nodemask)
Sets the mask, which is a threshold integer such that nodes in the complex that contain a number of d...
Definition: GIC.h:219
Persistence_diagram compute_PD()
Computes the extended persistence diagram of the complex.
Definition: GIC.h:1165
double compute_distance_from_confidence_level(double alpha)
Computes the bottleneck distance threshold corresponding to a specific confidence level.
Definition: GIC.h:1279
void set_graph_from_rips(double threshold, Distance distance)
Creates a graph G from a Rips complex.
Definition: GIC.h:360
void set_cover_from_file(const std::string &cover_file_name)
Creates the cover C from a file containing the cover elements of each point (the order has to be the ...
Definition: GIC.h:825
void set_graph_from_file(const std::string &graph_file_name)
Creates a graph G from a file containing the edges.
Definition: GIC.h:327
void set_cover_from_range(AssignmentRange const &assignments)
Creates the cover C from a vector of assignments stored in memory. The assignments,...
Definition: GIC.h:863
void create_complex(SimplicialComplex &complex)
Creates the simplicial complex.
Definition: GIC.h:1325
void set_type(const std::string &t)
Specifies whether the type of the output simplicial complex.
Definition: GIC.h:188
void set_distances_from_range(const std::vector< std::vector< double > > &distance_matrix)
Reads and stores the distance matrices from vector stored in memory.
Definition: GIC.h:390
void find_simplices()
Computes the simplices of the simplicial complex.
Definition: GIC.h:1339
void set_function_from_range(InputRange const &function)
Creates the function f from a vector stored in memory.
Definition: GIC.h:525
void set_cover_from_function()
Creates a cover C from the preimages of the function f.
Definition: GIC.h:601
double subcolor(int c)
Returns the mean color corresponding to a specific node of the created complex.
Definition: GIC.h:974
void set_color_from_file(const std::string &color_file_name)
Computes the function used to color the nodes of the simplicial complex from a file containing the fu...
Definition: GIC.h:987
void write_info()
Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e....
Definition: GIC.h:1078
void plot_OFF()
Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC....
Definition: GIC.h:1117
void plot_DOT()
Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial comple...
Definition: GIC.h:1034
void set_color_from_range(std::vector< double > color)
Computes the function used to color the nodes of the simplicial complex from a vector stored in memor...
Definition: GIC.h:1025
void set_subsampling(double constant, double power)
Sets the constants used to subsample the data set. These constants are explained in .
Definition: GIC.h:206
void set_point_cloud_from_range(const std::vector< std::vector< double > > &point_cloud)
Reads and stores the input point cloud from vector stored in memory.
Definition: GIC.h:229
double set_graph_from_automatic_rips(Distance distance, int N=100)
Creates a graph G from a Rips complex whose threshold value is automatically tuned with subsampling—s...
Definition: GIC.h:431
void set_resolution_with_interval_length(double reso)
Sets a length of intervals from a value stored in memory.
Definition: GIC.h:583
void set_gain(double g=0.3)
Sets a gain from a value stored in memory (default value 0.3).
Definition: GIC.h:595
void set_function_from_coordinate(int k)
Creates the function f from the k-th coordinate of the point cloud P.
Definition: GIC.h:504
const std::vector< int > & subpopulation(int c)
Returns the data subset corresponding to a specific node of the created complex.
Definition: GIC.h:964
void set_color_from_coordinate(int k=0)
Computes the function used to color the nodes of the simplicial complex from the k-th coordinate.
Definition: GIC.h:1005
void compute_distribution(unsigned int N=100)
Computes bootstrapped distances distribution.
Definition: GIC.h:1228
double compute_confidence_level_from_distance(double d)
Computes the confidence level of a specific bottleneck distance threshold.
Definition: GIC.h:1292
void set_graph_from_OFF()
Creates a graph G from the triangulation given by the input .OFF file.
Definition: GIC.h:344
double compute_p_value()
Computes the p-value, i.e. the opposite of the confidence level of the largest bottleneck distance pr...
Definition: GIC.h:1306
bool read_point_cloud(const std::string &off_file_name)
Reads and stores the input point cloud from .(n)OFF file.
Definition: GIC.h:244
void set_verbose(bool verb=false)
Specifies whether the program should display information or not.
Definition: GIC.h:196
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:54
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:677
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:168
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:152
Rips complex data structure.
Definition: Rips_complex.h:45
Global distance functions.
Graph simplicial complex methods.
double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=(std::numeric_limits< double >::min)())
Function to compute the Bottleneck distance between two persistence diagrams.
Definition: Bottleneck.h:116
Gudhi namespace.
Definition: SimplicialComplexForAlpha.h:14
This file includes common file reader for GUDHI.
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20