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 <tbb/mutex.h>
18 #endif
19 
20 #include <gudhi/Debug_utils.h>
21 #include <gudhi/graph_simplicial_complex.h>
22 #include <gudhi/reader_utils.h>
23 #include <gudhi/Simplex_tree.h>
24 #include <gudhi/Rips_complex.h>
25 #include <gudhi/Points_off_io.h>
27 #include <gudhi/Persistent_cohomology.h>
28 #include <gudhi/Bottleneck.h>
29 
30 #include <boost/config.hpp>
31 #include <boost/graph/graph_traits.hpp>
32 #include <boost/graph/adjacency_list.hpp>
33 #include <boost/graph/connected_components.hpp>
34 #include <boost/graph/dijkstra_shortest_paths.hpp>
35 #include <boost/graph/subgraph.hpp>
36 #include <boost/graph/graph_utility.hpp>
37 
38 #include <CGAL/version.h> // for CGAL_VERSION_NR
39 
40 #include <iostream>
41 #include <vector>
42 #include <map>
43 #include <string>
44 #include <limits> // for numeric_limits
45 #include <utility> // for std::pair<>
46 #include <algorithm> // for (std::max)
47 #include <random>
48 #include <cassert>
49 #include <cmath>
50 
51 // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10
52 #if CGAL_VERSION_NR < 1041101000
53 # error Alpha_complex_3d is only available for CGAL >= 4.11
54 #endif
55 
56 namespace Gudhi {
57 
58 namespace cover_complex {
59 
63 using Persistence_diagram = std::vector<std::pair<double, double> >;
64 using Graph = boost::subgraph<
65  boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
66  boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
67 using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
68 using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
69 using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
70 
91 template <typename Point>
93  private:
94  bool verbose = false; // whether to display information.
95  std::string type; // Nerve or GIC
96 
97  std::vector<Point> point_cloud; // input point cloud.
98  std::vector<std::vector<double> > distances; // all pairwise distances.
99  int maximal_dim; // maximal dimension of output simplicial complex.
100  int data_dimension; // dimension of input data.
101  int n; // number of points.
102 
103  std::vector<double> func; // function used to compute the output simplicial complex.
104  std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex.
105  bool functional_cover = false; // whether we use a cover with preimages of a function or not.
106 
107  Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists).
108  Graph one_skeleton; // one-skeleton used to compute the connected components.
109  std::vector<Vertex_t> vertices; // vertices of one_skeleton.
110 
111  std::vector<std::vector<int> > simplices; // simplices of output simplicial complex.
112  std::vector<int> voronoi_subsamples; // Voronoi germs (in case of Voronoi cover).
113 
114  Persistence_diagram PD;
115  std::vector<double> distribution;
116 
117  std::vector<std::vector<int> >
118  cover; // function associating to each data point the vector of cover elements to which it belongs.
119  std::map<int, std::vector<int> >
120  cover_back; // inverse of cover, in order to get the data points associated to a specific cover element.
121  std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence
122  // diagram of the output simplicial complex.
123  std::map<int, int>
124  cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not.
125  std::map<int, std::pair<int, double> >
126  cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex.
127 
128  int resolution_int = -1;
129  double resolution_double = -1;
130  double gain = -1;
131  double rate_constant = 10; // Constant in the subsampling.
132  double rate_power = 0.001; // Power in the subsampling.
133  int mask = 0; // Ignore nodes containing less than mask points.
134 
135  std::map<int, int> name2id, name2idinv;
136 
137  std::string cover_name;
138  std::string point_cloud_name;
139  std::string color_name;
140 
141  // Remove all edges of a graph.
142  void remove_edges(Graph& G) {
143  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
144  for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
145  }
146 
147  // Thread local is not available on XCode version < V.8
148  // If not available, random engine is a class member.
149 #ifndef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
150  std::default_random_engine re;
151 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
152 
153  // Find random number in [0,1].
154  double GetUniform() {
155  // Thread local is not available on XCode version < V.8
156  // If available, random engine is defined for each thread.
157 #ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
158  thread_local std::default_random_engine re;
159 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
160  std::uniform_real_distribution<double> Dist(0, 1);
161  return Dist(re);
162  }
163 
164  // Subsample points.
165  void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector<int>& samples) {
166  int t = 0;
167  int m = 0;
168  double u;
169  while (m < sampleSize) {
170  u = GetUniform();
171  if ((populationSize - t) * u >= sampleSize - m) {
172  t++;
173  } else {
174  samples[m] = t;
175  t++;
176  m++;
177  }
178  }
179  }
180 
181  // *******************************************************************************************************************
182  // Utils.
183  // *******************************************************************************************************************
184 
185  public:
191  void set_type(const std::string& t) { type = t; }
192 
193  public:
199  void set_verbose(bool verb = false) { verbose = verb; }
200 
201  public:
209  void set_subsampling(double constant, double power) {
210  rate_constant = constant;
211  rate_power = power;
212  }
213 
214  public:
222  void set_mask(int nodemask) { mask = nodemask; }
223 
224  public:
225 
226 
232  void set_point_cloud_from_range(const std::vector<std::vector<double> > & point_cloud) {
233  n = point_cloud.size(); data_dimension = point_cloud[0].size();
234  point_cloud_name = "matrix"; cover.resize(n);
235  for(int i = 0; i < n; i++){
236  boost::add_vertex(one_skeleton_OFF);
237  vertices.push_back(boost::add_vertex(one_skeleton));
238  }
239  this->point_cloud = point_cloud;
240  }
241 
247  bool read_point_cloud(const std::string& off_file_name) {
248  point_cloud_name = off_file_name;
249  std::ifstream input(off_file_name);
250  std::string line;
251 
252  char comment = '#';
253  while (comment == '#') {
254  std::getline(input, line);
255  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
256  comment = line[line.find_first_not_of(' ')];
257  }
258  if (strcmp((char*)line.c_str(), "nOFF") == 0) {
259  comment = '#';
260  while (comment == '#') {
261  std::getline(input, line);
262  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
263  comment = line[line.find_first_not_of(' ')];
264  }
265  std::stringstream stream(line);
266  stream >> data_dimension;
267  } else {
268  data_dimension = 3;
269  }
270 
271  comment = '#';
272  int numedges, numfaces, i, dim;
273  while (comment == '#') {
274  std::getline(input, line);
275  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
276  comment = line[line.find_first_not_of(' ')];
277  }
278  std::stringstream stream(line);
279  stream >> n;
280  stream >> numfaces;
281  stream >> numedges;
282 
283  i = 0;
284  while (i < n) {
285  std::getline(input, line);
286  if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
287  !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
288  std::stringstream iss(line);
289  std::vector<double> point;
290  point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
291  point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
292  boost::add_vertex(one_skeleton_OFF);
293  vertices.push_back(boost::add_vertex(one_skeleton));
294  cover.emplace_back();
295  i++;
296  }
297  }
298 
299  i = 0;
300  while (i < numfaces) {
301  std::getline(input, line);
302  if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
303  !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
304  std::vector<int> simplex;
305  std::stringstream iss(line);
306  simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>());
307  dim = simplex[0];
308  for (int j = 1; j <= dim; j++)
309  for (int k = j + 1; k <= dim; k++)
310  boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF);
311  i++;
312  }
313  }
314 
315  return input.is_open();
316  }
317 
318  // *******************************************************************************************************************
319  // Graphs.
320  // *******************************************************************************************************************
321 
322  public: // Set graph from file.
330  void set_graph_from_file(const std::string& graph_file_name) {
331  remove_edges(one_skeleton);
332  int neighb;
333  std::ifstream input(graph_file_name);
334  std::string line;
335  int source;
336  while (std::getline(input, line)) {
337  std::stringstream stream(line);
338  stream >> source;
339  while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
340  }
341  }
342 
343  public: // Set graph from OFF file.
348  remove_edges(one_skeleton);
349  if (num_edges(one_skeleton_OFF))
350  one_skeleton = one_skeleton_OFF;
351  else
352  std::cout << "No triangulation read in OFF file!" << std::endl;
353  }
354 
355  public: // Set graph from Rips complex.
362  template <typename Distance>
363  void set_graph_from_rips(double threshold, Distance distance) {
364  remove_edges(one_skeleton);
365  if (distances.size() == 0) compute_pairwise_distances(distance);
366  for (int i = 0; i < n; i++) {
367  for (int j = i + 1; j < n; j++) {
368  if (distances[i][j] <= threshold) {
369  boost::add_edge(vertices[i], vertices[j], one_skeleton);
370  boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first,
371  distances[i][j]);
372  }
373  }
374  }
375  }
376 
377  public:
378  void set_graph_weights() {
379  Index_map index = boost::get(boost::vertex_index, one_skeleton);
380  Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
381  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
382  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
383  boost::put(weight, *ei,
384  distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]);
385  }
386 
387  public:
393  void set_distances_from_range(const std::vector<std::vector<double> > & distance_matrix) {
394  n = distance_matrix.size(); data_dimension = 0; point_cloud_name = "matrix";
395  cover.resize(n); point_cloud.resize(n);
396  for(int i = 0; i < n; i++){
397  boost::add_vertex(one_skeleton_OFF);
398  vertices.push_back(boost::add_vertex(one_skeleton));
399  }
400  distances = distance_matrix;
401  }
402 
403  public: // Pairwise distances.
406  template <typename Distance>
407  void compute_pairwise_distances(Distance ref_distance) {
408  double d;
409  std::vector<double> zeros(n);
410  for (int i = 0; i < n; i++) distances.push_back(zeros);
411  std::string distance = point_cloud_name + "_dist";
412  std::ifstream input(distance, std::ios::out | std::ios::binary);
413 
414  if (input.good()) {
415  if (verbose) std::cout << "Reading distances..." << std::endl;
416  for (int i = 0; i < n; i++) {
417  for (int j = i; j < n; j++) {
418  input.read((char*)&d, 8);
419  distances[i][j] = d;
420  distances[j][i] = d;
421  }
422  }
423  input.close();
424  } else {
425  if (verbose) std::cout << "Computing distances..." << std::endl;
426  input.close();
427  std::ofstream output(distance, std::ios::out | std::ios::binary);
428  for (int i = 0; i < n; i++) {
429  int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10;
430  if (state == 0 && verbose) std::cout << "\r" << state << "%" << std::flush;
431  for (int j = i; j < n; j++) {
432  double dis = ref_distance(point_cloud[i], point_cloud[j]);
433  distances[i][j] = dis;
434  distances[j][i] = dis;
435  output.write((char*)&dis, 8);
436  }
437  }
438  output.close();
439  if (verbose) std::cout << std::endl;
440  }
441  }
442 
443  public: // Automatic tuning of Rips complex.
453  template <typename Distance>
454  double set_graph_from_automatic_rips(Distance distance, int N = 100) {
455  int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant))));
456  m = (std::min)(m, n - 1);
457  double delta = 0;
458 
459  if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl;
460  if (verbose) std::cout << "Subsampling " << m << " points" << std::endl;
461 
462  if (distances.size() == 0) compute_pairwise_distances(distance);
463 
464  // This cannot be parallelized if thread_local is not defined
465  // thread_local is not defined for XCode < v.8
466  #if defined(GUDHI_USE_TBB) && defined(GUDHI_CAN_USE_CXX11_THREAD_LOCAL)
467  tbb::mutex deltamutex;
468  tbb::parallel_for(0, N, [&](int i){
469  std::vector<int> samples(m);
470  SampleWithoutReplacement(n, m, samples);
471  double hausdorff_dist = 0;
472  for (int j = 0; j < n; j++) {
473  double mj = distances[j][samples[0]];
474  for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
475  hausdorff_dist = (std::max)(hausdorff_dist, mj);
476  }
477  deltamutex.lock();
478  delta += hausdorff_dist / N;
479  deltamutex.unlock();
480  });
481  #else
482  for (int i = 0; i < N; i++) {
483  std::vector<int> samples(m);
484  SampleWithoutReplacement(n, m, samples);
485  double hausdorff_dist = 0;
486  for (int j = 0; j < n; j++) {
487  double mj = distances[j][samples[0]];
488  for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
489  hausdorff_dist = (std::max)(hausdorff_dist, mj);
490  }
491  delta += hausdorff_dist / N;
492  }
493  #endif
494 
495  if (verbose) std::cout << "delta = " << delta << std::endl;
496  set_graph_from_rips(delta, distance);
497  return delta;
498  }
499 
500  // *******************************************************************************************************************
501  // Functions.
502  // *******************************************************************************************************************
503 
504  public: // Set function from file.
510  void set_function_from_file(const std::string& func_file_name) {
511  int i = 0;
512  std::ifstream input(func_file_name);
513  std::string line;
514  double f;
515  while (std::getline(input, line)) {
516  std::stringstream stream(line);
517  stream >> f;
518  func.push_back(f);
519  i++;
520  }
521  functional_cover = true;
522  cover_name = func_file_name;
523  }
524 
525  public: // Set function from kth coordinate
532  if(point_cloud[0].size() > 0){
533  for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]);
534  functional_cover = true;
535  cover_name = "coordinate " + std::to_string(k);
536  }
537  else{
538  std::cout << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl;
539  for (int i = 0; i < n; i++) func.push_back(0.0);
540  functional_cover = true;
541  cover_name = "null";
542  }
543  }
544 
545  public: // Set function from vector.
551  template <class InputRange>
552  void set_function_from_range(InputRange const& function) {
553  for (int i = 0; i < n; i++) func.push_back(function[i]);
554  functional_cover = true;
555  }
556 
557  // *******************************************************************************************************************
558  // Covers.
559  // *******************************************************************************************************************
560 
561  public: // Automatic tuning of resolution.
570  if (!functional_cover) {
571  std::cout << "Cover needs to come from the preimages of a function." << std::endl;
572  return 0;
573  }
574  if (type != "Nerve" && type != "GIC") {
575  std::cout << "Type of complex needs to be specified." << std::endl;
576  return 0;
577  }
578 
579  double reso = 0;
580  Index_map index = boost::get(boost::vertex_index, one_skeleton);
581 
582  if (type == "GIC") {
583  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
584  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
585  reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
586  func[index[boost::target(*ei, one_skeleton)]]));
587  if (verbose) std::cout << "resolution = " << reso << std::endl;
588  resolution_double = reso;
589  }
590 
591  if (type == "Nerve") {
592  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
593  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
594  reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
595  func[index[boost::target(*ei, one_skeleton)]]) /
596  gain);
597  if (verbose) std::cout << "resolution = " << reso << std::endl;
598  resolution_double = reso;
599  }
600 
601  return reso;
602  }
603 
604  public:
610  void set_resolution_with_interval_length(double reso) { resolution_double = reso; }
616  void set_resolution_with_interval_number(int reso) { resolution_int = reso; }
622  void set_gain(double g = 0.3) { gain = g; }
623 
624  public: // Set cover with preimages of function.
629  if (resolution_double == -1 && resolution_int == -1) {
630  std::cout << "Number and/or length of intervals not specified" << std::endl;
631  return;
632  }
633  if (gain == -1) {
634  std::cout << "Gain not specified" << std::endl;
635  return;
636  }
637 
638  // Read function values and compute min and max
639  double minf = (std::numeric_limits<float>::max)();
640  double maxf = std::numeric_limits<float>::lowest();
641  for (int i = 0; i < n; i++) {
642  minf = (std::min)(minf, func[i]);
643  maxf = (std::max)(maxf, func[i]);
644  }
645  if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl;
646 
647  // Compute cover of im(f)
648  std::vector<std::pair<double, double> > intervals;
649  int res;
650 
651  if (resolution_double == -1) { // Case we use an integer for the number of intervals.
652  double incr = (maxf - minf) / resolution_int;
653  double x = minf;
654  double alpha = (incr * gain) / (2 - 2 * gain);
655  double y = minf + incr + alpha;
656  std::pair<double, double> interm(x, y);
657  intervals.push_back(interm);
658  for (int i = 1; i < resolution_int - 1; i++) {
659  x = minf + i * incr - alpha;
660  y = minf + (i + 1) * incr + alpha;
661  std::pair<double, double> inter(x, y);
662  intervals.push_back(inter);
663  }
664  x = minf + (resolution_int - 1) * incr - alpha;
665  y = maxf;
666  std::pair<double, double> interM(x, y);
667  intervals.push_back(interM);
668  res = intervals.size();
669  if (verbose) {
670  for (int i = 0; i < res; i++)
671  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
672  << std::endl;
673  }
674  } else {
675  if (resolution_int == -1) { // Case we use a double for the length of the intervals.
676  double x = minf;
677  double y = x + resolution_double;
678  while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
679  std::pair<double, double> inter(x, y);
680  intervals.push_back(inter);
681  x = y - gain * resolution_double;
682  y = x + resolution_double;
683  }
684  std::pair<double, double> interM(x, maxf);
685  intervals.push_back(interM);
686  res = intervals.size();
687  if (verbose) {
688  for (int i = 0; i < res; i++)
689  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
690  << std::endl;
691  }
692  } else { // Case we use an integer and a double for the length of the intervals.
693  double x = minf;
694  double y = x + resolution_double;
695  int count = 0;
696  while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
697  std::pair<double, double> inter(x, y);
698  intervals.push_back(inter);
699  count++;
700  x = y - gain * resolution_double;
701  y = x + resolution_double;
702  }
703  res = intervals.size();
704  if (verbose) {
705  for (int i = 0; i < res; i++)
706  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
707  << std::endl;
708  }
709  }
710  }
711 
712  // Sort points according to function values
713  std::vector<int> points(n);
714  for (int i = 0; i < n; i++) points[i] = i;
715  std::sort(points.begin(), points.end(), [=](const int & p1, const int & p2){return (this->func[p1] < this->func[p2]);});
716 
717  int id = 0;
718  int pos = 0;
719  Index_map index = boost::get(boost::vertex_index, one_skeleton); // int maxc = -1;
720  std::map<int, std::vector<int> > preimages;
721  std::map<int, double> funcstd;
722 
723  if (verbose) std::cout << "Computing preimages..." << std::endl;
724  for (int i = 0; i < res; i++) {
725  // Find points in the preimage
726  std::pair<double, double> inter1 = intervals[i];
727  int tmp = pos;
728  double u, v;
729 
730  if (i != res - 1) {
731  if (i != 0) {
732  std::pair<double, double> inter3 = intervals[i - 1];
733  while (func[points[tmp]] < inter3.second && tmp != n) {
734  preimages[i].push_back(points[tmp]);
735  tmp++;
736  }
737  u = inter3.second;
738  } else {
739  u = inter1.first;
740  }
741 
742  std::pair<double, double> inter2 = intervals[i + 1];
743  while (func[points[tmp]] < inter2.first && tmp != n) {
744  preimages[i].push_back(points[tmp]);
745  tmp++;
746  }
747  v = inter2.first;
748  pos = tmp;
749  while (func[points[tmp]] < inter1.second && tmp != n) {
750  preimages[i].push_back(points[tmp]);
751  tmp++;
752  }
753 
754  } else {
755  std::pair<double, double> inter3 = intervals[i - 1];
756  while (func[points[tmp]] < inter3.second && tmp != n) {
757  preimages[i].push_back(points[tmp]);
758  tmp++;
759  }
760  while (tmp != n) {
761  preimages[i].push_back(points[tmp]);
762  tmp++;
763  }
764  u = inter3.second;
765  v = inter1.second;
766  }
767 
768  funcstd[i] = 0.5 * (u + v);
769  }
770 
771  #ifdef GUDHI_USE_TBB
772  if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl;
773  tbb::mutex covermutex, idmutex;
774  tbb::parallel_for(0, res, [&](int i){
775  // Compute connected components
776  Graph G = one_skeleton.create_subgraph();
777  int num = preimages[i].size();
778  std::vector<int> component(num);
779  for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
780  boost::connected_components(G, &component[0]);
781  int max = 0;
782 
783  // For each point in preimage
784  for (int j = 0; j < num; j++) {
785  // Update number of components in preimage
786  if (component[j] > max) max = component[j];
787 
788  // Identify component with Cantor polynomial N^2 -> N
789  int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
790 
791  // Update covers
792  covermutex.lock();
793  cover[preimages[i][j]].push_back(identifier);
794  cover_back[identifier].push_back(preimages[i][j]);
795  cover_fct[identifier] = i;
796  cover_std[identifier] = funcstd[i];
797  cover_color[identifier].second += func_color[preimages[i][j]];
798  cover_color[identifier].first += 1;
799  covermutex.unlock();
800  }
801 
802  // Maximal dimension is total number of connected components
803  idmutex.lock();
804  id += max + 1;
805  idmutex.unlock();
806  });
807  #else
808  if (verbose) std::cout << "Computing connected components..." << std::endl;
809  for (int i = 0; i < res; i++) {
810  // Compute connected components
811  Graph G = one_skeleton.create_subgraph();
812  int num = preimages[i].size();
813  std::vector<int> component(num);
814  for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
815  boost::connected_components(G, &component[0]);
816  int max = 0;
817 
818  // For each point in preimage
819  for (int j = 0; j < num; j++) {
820  // Update number of components in preimage
821  if (component[j] > max) max = component[j];
822 
823  // Identify component with Cantor polynomial N^2 -> N
824  int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
825 
826  // Update covers
827  cover[preimages[i][j]].push_back(identifier);
828  cover_back[identifier].push_back(preimages[i][j]);
829  cover_fct[identifier] = i;
830  cover_std[identifier] = funcstd[i];
831  cover_color[identifier].second += func_color[preimages[i][j]];
832  cover_color[identifier].first += 1;
833  }
834 
835  // Maximal dimension is total number of connected components
836  id += max + 1;
837  }
838  #endif
839 
840  maximal_dim = id - 1;
841  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
842  iit->second.second /= iit->second.first;
843  }
844 
845  public: // Set cover from file.
852  void set_cover_from_file(const std::string& cover_file_name) {
853  int i = 0;
854  int cov;
855  std::vector<int> cov_elts, cov_number;
856  std::ifstream input(cover_file_name);
857  std::string line;
858  while (std::getline(input, line)) {
859  cov_elts.clear();
860  std::stringstream stream(line);
861  while (stream >> cov) {
862  cov_elts.push_back(cov);
863  cov_number.push_back(cov);
864  cover_fct[cov] = cov;
865  cover_color[cov].second += func_color[i];
866  cover_color[cov].first++;
867  cover_back[cov].push_back(i);
868  }
869  cover[i] = cov_elts;
870  i++;
871  }
872 
873  std::sort(cov_number.begin(), cov_number.end());
874  std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
875  cov_number.resize(std::distance(cov_number.begin(), it));
876 
877  maximal_dim = cov_number.size() - 1;
878  for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
879  cover_name = cover_file_name;
880  }
881 
882  public: // Set cover from Voronoi
889  template <typename Distance>
890  void set_cover_from_Voronoi(Distance distance, int m = 100) {
891  voronoi_subsamples.resize(m);
892  SampleWithoutReplacement(n, m, voronoi_subsamples);
893  if (distances.size() == 0) compute_pairwise_distances(distance);
894  set_graph_weights();
895  Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
896  Index_map index = boost::get(boost::vertex_index, one_skeleton);
897  std::vector<double> mindist(n);
898  for (int j = 0; j < n; j++) mindist[j] = (std::numeric_limits<double>::max)();
899 
900  // Compute the geodesic distances to subsamples with Dijkstra
901  #ifdef GUDHI_USE_TBB
902  if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl;
903  tbb::mutex coverMutex; tbb::mutex mindistMutex;
904  tbb::parallel_for(0, m, [&](int i){
905  int seed = voronoi_subsamples[i];
906  std::vector<double> dmap(n);
907  boost::dijkstra_shortest_paths(
908  one_skeleton, vertices[seed],
909  boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
910 
911  coverMutex.lock(); mindistMutex.lock();
912  for (int j = 0; j < n; j++)
913  if (mindist[j] > dmap[j]) {
914  mindist[j] = dmap[j];
915  if (cover[j].size() == 0)
916  cover[j].push_back(i);
917  else
918  cover[j][0] = i;
919  }
920  coverMutex.unlock(); mindistMutex.unlock();
921  });
922  #else
923  for (int i = 0; i < m; i++) {
924  if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
925  int seed = voronoi_subsamples[i];
926  std::vector<double> dmap(n);
927  boost::dijkstra_shortest_paths(
928  one_skeleton, vertices[seed],
929  boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
930 
931  for (int j = 0; j < n; j++)
932  if (mindist[j] > dmap[j]) {
933  mindist[j] = dmap[j];
934  if (cover[j].size() == 0)
935  cover[j].push_back(i);
936  else
937  cover[j][0] = i;
938  }
939  }
940  #endif
941 
942  for (int i = 0; i < n; i++) {
943  cover_back[cover[i][0]].push_back(i);
944  cover_color[cover[i][0]].second += func_color[i];
945  cover_color[cover[i][0]].first++;
946  }
947  for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first;
948  maximal_dim = m - 1;
949  cover_name = "Voronoi";
950  }
951 
952  public: // return subset of data corresponding to a node
959  const std::vector<int>& subpopulation(int c) { return cover_back[name2idinv[c]]; }
960 
961  // *******************************************************************************************************************
962  // Visualization.
963  // *******************************************************************************************************************
964 
965  public: // Set color from file.
972  void set_color_from_file(const std::string& color_file_name) {
973  int i = 0;
974  std::ifstream input(color_file_name);
975  std::string line;
976  double f;
977  while (std::getline(input, line)) {
978  std::stringstream stream(line);
979  stream >> f;
980  func_color.push_back(f);
981  i++;
982  }
983  color_name = color_file_name;
984  }
985 
986  public: // Set color from kth coordinate
992  void set_color_from_coordinate(int k = 0) {
993  if(point_cloud[0].size() > 0){
994  for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]);
995  color_name = "coordinate ";
996  color_name.append(std::to_string(k));
997  }
998  else{
999  std::cout << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl;
1000  for (int i = 0; i < n; i++) func.push_back(0.0);
1001  functional_cover = true;
1002  cover_name = "null";
1003  }
1004  }
1005 
1006  public: // Set color from vector.
1012  void set_color_from_range(std::vector<double> color) {
1013  for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
1014  }
1015 
1016  public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
1021  void plot_DOT() {
1022  std::string mapp = point_cloud_name + "_sc.dot";
1023  std::ofstream graphic(mapp);
1024 
1025  double maxv = std::numeric_limits<double>::lowest();
1026  double minv = (std::numeric_limits<double>::max)();
1027  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1028  maxv = (std::max)(maxv, iit->second.second);
1029  minv = (std::min)(minv, iit->second.second);
1030  }
1031 
1032  int k = 0;
1033  std::vector<int> nodes;
1034  nodes.clear();
1035 
1036  graphic << "graph GIC {" << std::endl;
1037  int id = 0;
1038  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1039  if (iit->second.first > mask) {
1040  nodes.push_back(iit->first);
1041  name2id[iit->first] = id;
1042  name2idinv[id] = iit->first;
1043  id++;
1044  graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
1045  << ":" << iit->second.first << "\" style=filled fillcolor=\""
1046  << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl;
1047  k++;
1048  }
1049  }
1050  int ke = 0;
1051  int num_simplices = simplices.size();
1052  for (int i = 0; i < num_simplices; i++)
1053  if (simplices[i].size() == 2) {
1054  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
1055  graphic << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];"
1056  << std::endl;
1057  ke++;
1058  }
1059  }
1060  graphic << "}";
1061  graphic.close();
1062  std::cout << mapp << " file generated. It can be visualized with e.g. neato." << std::endl;
1063  }
1064 
1065  public: // Create a .txt file that can be compiled with KeplerMapper.
1069  void write_info() {
1070  int num_simplices = simplices.size();
1071  int num_edges = 0;
1072  std::string mapp = point_cloud_name + "_sc.txt";
1073  std::ofstream graphic(mapp);
1074 
1075  for (int i = 0; i < num_simplices; i++)
1076  if (simplices[i].size() == 2)
1077  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++;
1078 
1079  graphic << point_cloud_name << std::endl;
1080  graphic << cover_name << std::endl;
1081  graphic << color_name << std::endl;
1082  graphic << resolution_double << " " << gain << std::endl;
1083  graphic << cover_color.size() << " " << num_edges << std::endl;
1084 
1085  int id = 0;
1086  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1087  graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl;
1088  name2id[iit->first] = id;
1089  name2idinv[id] = iit->first;
1090  id++;
1091  }
1092 
1093  for (int i = 0; i < num_simplices; i++)
1094  if (simplices[i].size() == 2)
1095  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
1096  graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl;
1097  graphic.close();
1098  std::cout << mapp
1099  << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
1100  << std::endl;
1101  }
1102 
1103  public: // Create a .off file that can be visualized (e.g. with Geomview).
1108  void plot_OFF() {
1109  assert(cover_name == "Voronoi");
1110 
1111  int m = voronoi_subsamples.size();
1112  int numedges = 0;
1113  int numfaces = 0;
1114  std::vector<std::vector<int> > edges, faces;
1115  int numsimplices = simplices.size();
1116 
1117  std::string mapp = point_cloud_name + "_sc.off";
1118  std::ofstream graphic(mapp);
1119 
1120  graphic << "OFF" << std::endl;
1121  for (int i = 0; i < numsimplices; i++) {
1122  if (simplices[i].size() == 2) {
1123  numedges++;
1124  edges.push_back(simplices[i]);
1125  }
1126  if (simplices[i].size() == 3) {
1127  numfaces++;
1128  faces.push_back(simplices[i]);
1129  }
1130  }
1131  graphic << m << " " << numedges + numfaces << std::endl;
1132  for (int i = 0; i < m; i++) {
1133  if (data_dimension <= 3) {
1134  for (int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1135  for (int j = data_dimension; j < 3; j++) graphic << 0 << " ";
1136  graphic << std::endl;
1137  } else {
1138  for (int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1139  }
1140  }
1141  for (int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl;
1142  for (int i = 0; i < numfaces; i++)
1143  graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
1144  graphic.close();
1145  std::cout << mapp << " generated. It can be visualized with e.g. geomview." << std::endl;
1146  }
1147 
1148  // *******************************************************************************************************************
1149  // Extended Persistence Diagrams.
1150  // *******************************************************************************************************************
1151 
1152  public:
1156  Persistence_diagram compute_PD() {
1157  Simplex_tree st;
1158 
1159  // Compute max and min
1160  double maxf = std::numeric_limits<double>::lowest();
1161  double minf = (std::numeric_limits<double>::max)();
1162  for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1163  maxf = (std::max)(maxf, it->second);
1164  minf = (std::min)(minf, it->second);
1165  }
1166 
1167  // Build filtration
1168  for (auto const& simplex : simplices) {
1169  std::vector<int> splx = simplex;
1170  splx.push_back(-2);
1171  st.insert_simplex_and_subfaces(splx, -3);
1172  }
1173 
1174  for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1175  int vertex = it->first; float val = it->second;
1176  int vert[] = {vertex}; int edge[] = {vertex, -2};
1177  if(st.find(vert) != st.null_simplex()){
1178  st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf));
1179  st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf));
1180  }
1181  }
1183 
1184  // Compute PD
1187 
1188  // Output PD
1189  int max_dim = st.dimension();
1190  for (int i = 0; i < max_dim; i++) {
1191  std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i);
1192  int num_bars = bars.size(); if(i == 0) num_bars -= 1;
1193  if(verbose) std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
1194  for (int j = 0; j < num_bars; j++) {
1195  double birth = bars[j].first;
1196  double death = bars[j].second;
1197  if (i == 0 && std::isinf(death)) continue;
1198  if (birth < 0)
1199  birth = minf + (birth + 2) * (maxf - minf);
1200  else
1201  birth = minf + (2 - birth) * (maxf - minf);
1202  if (death < 0)
1203  death = minf + (death + 2) * (maxf - minf);
1204  else
1205  death = minf + (2 - death) * (maxf - minf);
1206  PD.push_back(std::pair<double, double>(birth, death));
1207  if (verbose) std::cout << " [" << birth << ", " << death << "]" << std::endl;
1208  }
1209  }
1210  return PD;
1211  }
1212 
1213  public:
1219  void compute_distribution(unsigned int N = 100) {
1220  unsigned int sz = distribution.size();
1221  if (sz >= N) {
1222  std::cout << "Already done!" << std::endl;
1223  } else {
1224  for (unsigned int i = 0; i < N - sz; i++) {
1225  if (verbose) std::cout << "Computing " << i << "th bootstrap, bottleneck distance = ";
1226 
1227  Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true;
1228 
1229  std::vector<int> boot(this->n);
1230  for (int j = 0; j < this->n; j++) {
1231  double u = GetUniform();
1232  int id = std::floor(u * (this->n)); boot[j] = id;
1233  Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]);
1234  boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
1235  }
1236  Cboot.set_color_from_range(Cboot.func);
1237 
1238  for (int j = 0; j < n; j++) {
1239  std::vector<double> dist(n);
1240  for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]];
1241  Cboot.distances.push_back(dist);
1242  }
1243 
1245  Cboot.set_gain();
1246  Cboot.set_automatic_resolution();
1247  Cboot.set_cover_from_function();
1248  Cboot.find_simplices();
1249  Cboot.compute_PD();
1250  double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD);
1251  if (verbose) std::cout << db << std::endl;
1252  distribution.push_back(db);
1253  }
1254 
1255  std::sort(distribution.begin(), distribution.end());
1256  }
1257  }
1258 
1259  public:
1266  unsigned int N = distribution.size();
1267  double d = distribution[std::floor(alpha * N)];
1268  if (verbose) std::cout << "Distance corresponding to confidence " << alpha << " is " << d << std::endl;
1269  return d;
1270  }
1271 
1272  public:
1279  unsigned int N = distribution.size();
1280  double level = 1;
1281  for (unsigned int i = 0; i < N; i++)
1282  if (distribution[i] >= d){ level = i * 1.0 / N; break; }
1283  if (verbose) std::cout << "Confidence level of distance " << d << " is " << level << std::endl;
1284  return level;
1285  }
1286 
1287  public:
1292  double compute_p_value() {
1293  double distancemin = (std::numeric_limits<double>::max)(); int N = PD.size();
1294  for (int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first));
1295  double p_value = 1 - compute_confidence_level_from_distance(distancemin);
1296  if (verbose) std::cout << "p value = " << p_value << std::endl;
1297  return p_value;
1298  }
1299 
1300  // *******************************************************************************************************************
1301  // Computation of simplices.
1302  // *******************************************************************************************************************
1303 
1304  public:
1310  template <typename SimplicialComplex>
1311  void create_complex(SimplicialComplex& complex) {
1312  unsigned int dimension = 0;
1313  for (auto const& simplex : simplices) {
1314  int numvert = simplex.size();
1315  double filt = std::numeric_limits<double>::lowest();
1316  for (int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt);
1317  complex.insert_simplex_and_subfaces(simplex, filt);
1318  if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
1319  }
1320  }
1321 
1322  public:
1326  if (type != "Nerve" && type != "GIC") {
1327  std::cout << "Type of complex needs to be specified." << std::endl;
1328  return;
1329  }
1330 
1331  if (type == "Nerve") {
1332  for(int i = 0; i < n; i++) simplices.push_back(cover[i]);
1333  std::sort(simplices.begin(), simplices.end());
1334  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1335  simplices.resize(std::distance(simplices.begin(), it));
1336  }
1337 
1338  if (type == "GIC") {
1339  Index_map index = boost::get(boost::vertex_index, one_skeleton);
1340 
1341  if (functional_cover) {
1342  // Computes the simplices in the GIC by looking at all the edges of the graph and adding the
1343  // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals.
1344 
1345  if (gain >= 0.5)
1346  throw std::invalid_argument(
1347  "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
1348 
1349  // Loop on all edges.
1350  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1351  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
1352  int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
1353  for (int i = 0; i < nums; i++) {
1354  int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
1355  int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
1356  for (int j = 0; j < numt; j++) {
1357  int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
1358  if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
1359  std::vector<int> edge(2);
1360  edge[0] = (std::min)(vs, vt);
1361  edge[1] = (std::max)(vs, vt);
1362  simplices.push_back(edge);
1363  goto afterLoop;
1364  }
1365  }
1366  }
1367  afterLoop:;
1368  }
1369  std::sort(simplices.begin(), simplices.end());
1370  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1371  simplices.resize(std::distance(simplices.begin(), it));
1372 
1373  } else {
1374  // Find edges to keep
1375  Simplex_tree st;
1376  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1377  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
1378  if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
1379  cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
1380  std::vector<int> edge(2);
1381  edge[0] = index[boost::source(*ei, one_skeleton)];
1382  edge[1] = index[boost::target(*ei, one_skeleton)];
1383  st.insert_simplex_and_subfaces(edge);
1384  }
1385 
1386  // st.insert_graph(one_skeleton);
1387 
1388  // Build the Simplex Tree corresponding to the graph
1389  st.expansion(maximal_dim);
1390 
1391  // Find simplices of GIC
1392  simplices.clear();
1393  for (auto simplex : st.complex_simplex_range()) {
1394  if (!st.has_children(simplex)) {
1395  std::vector<int> simplx;
1396  for (auto vertex : st.simplex_vertex_range(simplex)) {
1397  unsigned int sz = cover[vertex].size();
1398  for (unsigned int i = 0; i < sz; i++) {
1399  simplx.push_back(cover[vertex][i]);
1400  }
1401  }
1402  std::sort(simplx.begin(), simplx.end());
1403  std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
1404  simplx.resize(std::distance(simplx.begin(), it));
1405  simplices.push_back(simplx);
1406  }
1407  }
1408  std::sort(simplices.begin(), simplices.end());
1409  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1410  simplices.resize(std::distance(simplices.begin(), it));
1411  }
1412  }
1413  }
1414 };
1415 
1416 } // namespace cover_complex
1417 
1418 } // namespace Gudhi
1419 
1420 #endif // GIC_H_
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:156
void plot_DOT()
Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial comple...
Definition: GIC.h:1021
bool read_point_cloud(const std::string &off_file_name)
Reads and stores the input point cloud from .(n)OFF file.
Definition: GIC.h:247
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1102
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:222
void set_function_from_range(InputRange const &function)
Creates the function f from a vector stored in memory.
Definition: GIC.h:552
void set_graph_from_file(const std::string &graph_file_name)
Creates a graph G from a file containing the edges.
Definition: GIC.h:330
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:52
const std::vector< int > & subpopulation(int c)
Returns the data subset corresponding to a specific node of the created complex.
Definition: GIC.h:959
Simplex Tree data structure for representing simplicial complexes.
Definition: Simplex_tree.h:60
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:750
double compute_distance_from_confidence_level(double alpha)
Computes the bottleneck distance threshold corresponding to a specific confidence level...
Definition: GIC.h:1265
void set_type(const std::string &t)
Specifies whether the type of the output simplicial complex.
Definition: GIC.h:191
void find_simplices()
Computes the simplices of the simplicial complex.
Definition: GIC.h:1325
void write_info()
Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e...
Definition: GIC.h:1069
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:498
Compute the Euclidean distance between two Points given by a range of coordinates. The points are assumed to have the same dimension.
Definition: distance_functions.h:34
void set_function_from_coordinate(int k)
Creates the function f from the k-th coordinate of the point cloud P.
Definition: GIC.h:531
Persistence_diagram compute_PD()
Computes the extended persistence diagram of the complex.
Definition: GIC.h:1156
Definition: SimplicialComplexForAlpha.h:14
Rips complex data structure.
Definition: Rips_complex.h:45
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:992
void create_complex(SimplicialComplex &complex)
Creates the simplicial complex.
Definition: GIC.h:1311
void set_graph_from_OFF()
Creates a graph G from the triangulation given by the input .OFF file.
Definition: GIC.h:347
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:510
void set_subsampling(double constant, double power)
Sets the constants used to subsample the data set. These constants are explained in ...
Definition: GIC.h:209
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:890
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:1012
void set_verbose(bool verb=false)
Specifies whether the program should display information or not.
Definition: GIC.h:199
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:1321
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:693
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:1292
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:584
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20
Complex_simplex_range complex_simplex_range()
Returns a range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:202
void assign_filtration(Simplex_handle sh, Filtration_value fv)
Sets the filtration value of a simplex.
Definition: Simplex_tree.h:488
Simplex_vertex_range simplex_vertex_range(Simplex_handle sh)
Returns a range over the vertices of a simplex.
Definition: Simplex_tree.h:249
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:393
Global distance functions.
void set_graph_from_rips(double threshold, Distance distance)
Creates a graph G from a Rips complex.
Definition: GIC.h:363
bool has_children(SimplexHandle sh) const
Returns true if the node in the simplex tree pointed by sh has children.
Definition: Simplex_tree.h:571
void compute_distribution(unsigned int N=100)
Computes bootstrapped distances distribution.
Definition: GIC.h:1219
void plot_OFF()
Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC...
Definition: GIC.h:1108
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:112
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Definition: Simplex_tree.h:543
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—...
Definition: GIC.h:454
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:972
double compute_confidence_level_from_distance(double d)
Computes the confidence level of a specific bottleneck distance threshold.
Definition: GIC.h:1278
void set_resolution_with_interval_number(int reso)
Sets a number of intervals from a value stored in memory.
Definition: GIC.h:616
void set_resolution_with_interval_length(double reso)
Sets a length of intervals from a value stored in memory.
Definition: GIC.h:610
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:852
void set_cover_from_function()
Creates a cover C from the preimages of the function f.
Definition: GIC.h:628
Options::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Simplex_tree.h:67
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:172
Cover complex data structure.
Definition: GIC.h:92
double set_automatic_resolution()
Computes the optimal length of intervals (i.e. the smallest interval length avoiding discretization a...
Definition: GIC.h:569
This file includes common file reader for GUDHI.
void set_gain(double g=0.3)
Sets a gain from a value stored in memory (default value 0.3).
Definition: GIC.h:622
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:232
GUDHI  Version 3.0.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : MIT Generated on Tue Dec 10 2019 13:33:51 for GUDHI by Doxygen 1.8.13