GIC.h
1 /* This file is part of the Gudhi Library. The Gudhi library
2  * (Geometric Understanding in Higher Dimensions) is a generic C++
3  * library for computational topology.
4  *
5  * Author: Mathieu Carriere
6  *
7  * Copyright (C) 2017 Inria
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #ifndef GIC_H_
24 #define GIC_H_
25 
26 #ifdef GUDHI_USE_TBB
27 #include <tbb/parallel_for.h>
28 #include <tbb/mutex.h>
29 #endif
30 
31 #include <gudhi/Debug_utils.h>
32 #include <gudhi/graph_simplicial_complex.h>
33 #include <gudhi/reader_utils.h>
34 #include <gudhi/Simplex_tree.h>
35 #include <gudhi/Rips_complex.h>
36 #include <gudhi/Points_off_io.h>
38 #include <gudhi/Persistent_cohomology.h>
39 #include <gudhi/Bottleneck.h>
40 
41 #include <boost/config.hpp>
42 #include <boost/graph/graph_traits.hpp>
43 #include <boost/graph/adjacency_list.hpp>
44 #include <boost/graph/connected_components.hpp>
45 #include <boost/graph/dijkstra_shortest_paths.hpp>
46 #include <boost/graph/subgraph.hpp>
47 #include <boost/graph/graph_utility.hpp>
48 
49 #include <iostream>
50 #include <vector>
51 #include <map>
52 #include <string>
53 #include <limits> // for numeric_limits
54 #include <utility> // for std::pair<>
55 #include <algorithm> // for std::max
56 #include <random>
57 #include <cassert>
58 #include <cmath>
59 
60 namespace Gudhi {
61 
62 namespace cover_complex {
63 
67 using Persistence_diagram = std::vector<std::pair<double, double> >;
68 using Graph = boost::subgraph<
69  boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
70  boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
71 using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
72 using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
73 using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
74 
95 template <typename Point>
97  private:
98  bool verbose = false; // whether to display information.
99  std::string type; // Nerve or GIC
100 
101  std::vector<Point> point_cloud; // input point cloud.
102  std::vector<std::vector<double> > distances; // all pairwise distances.
103  int maximal_dim; // maximal dimension of output simplicial complex.
104  int data_dimension; // dimension of input data.
105  int n; // number of points.
106 
107  std::vector<double> func; // function used to compute the output simplicial complex.
108  std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex.
109  bool functional_cover = false; // whether we use a cover with preimages of a function or not.
110 
111  Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists).
112  Graph one_skeleton; // one-skeleton used to compute the connected components.
113  std::vector<Vertex_t> vertices; // vertices of one_skeleton.
114 
115  std::vector<std::vector<int> > simplices; // simplices of output simplicial complex.
116  std::vector<int> voronoi_subsamples; // Voronoi germs (in case of Voronoi cover).
117 
118  Persistence_diagram PD;
119  std::vector<double> distribution;
120 
121  std::vector<std::vector<int> >
122  cover; // function associating to each data point the vector of cover elements to which it belongs.
123  std::map<int, std::vector<int> >
124  cover_back; // inverse of cover, in order to get the data points associated to a specific cover element.
125  std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence
126  // diagram of the output simplicial complex.
127  std::map<int, int>
128  cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not.
129  std::map<int, std::pair<int, double> >
130  cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex.
131 
132  int resolution_int = -1;
133  double resolution_double = -1;
134  double gain = -1;
135  double rate_constant = 10; // Constant in the subsampling.
136  double rate_power = 0.001; // Power in the subsampling.
137  int mask = 0; // Ignore nodes containing less than mask points.
138 
139  std::map<int, int> name2id, name2idinv;
140 
141  std::string cover_name;
142  std::string point_cloud_name;
143  std::string color_name;
144 
145  // Remove all edges of a graph.
146  void remove_edges(Graph& G) {
147  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
148  for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
149  }
150 
151  // Thread local is not available on XCode version < V.8
152  // If not available, random engine is a class member.
153 #ifndef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
154  std::default_random_engine re;
155 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
156 
157  // Find random number in [0,1].
158  double GetUniform() {
159  // Thread local is not available on XCode version < V.8
160  // If available, random engine is defined for each thread.
161 #ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL
162  thread_local std::default_random_engine re;
163 #endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
164  std::uniform_real_distribution<double> Dist(0, 1);
165  return Dist(re);
166  }
167 
168  // Subsample points.
169  void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector<int>& samples) {
170  int t = 0;
171  int m = 0;
172  double u;
173  while (m < sampleSize) {
174  u = GetUniform();
175  if ((populationSize - t) * u >= sampleSize - m) {
176  t++;
177  } else {
178  samples[m] = t;
179  t++;
180  m++;
181  }
182  }
183  }
184 
185  // *******************************************************************************************************************
186  // Utils.
187  // *******************************************************************************************************************
188 
189  public:
195  void set_type(const std::string& t) { type = t; }
196 
197  public:
203  void set_verbose(bool verb = false) { verbose = verb; }
204 
205  public:
213  void set_subsampling(double constant, double power) {
214  rate_constant = constant;
215  rate_power = power;
216  }
217 
218  public:
226  void set_mask(int nodemask) { mask = nodemask; }
227 
228  public:
234  bool read_point_cloud(const std::string& off_file_name) {
235  point_cloud_name = off_file_name;
236  std::ifstream input(off_file_name);
237  std::string line;
238 
239  char comment = '#';
240  while (comment == '#') {
241  std::getline(input, line);
242  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
243  comment = line[line.find_first_not_of(' ')];
244  }
245  if (strcmp((char*)line.c_str(), "nOFF") == 0) {
246  comment = '#';
247  while (comment == '#') {
248  std::getline(input, line);
249  if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
250  comment = line[line.find_first_not_of(' ')];
251  }
252  std::stringstream stream(line);
253  stream >> data_dimension;
254  } else {
255  data_dimension = 3;
256  }
257 
258  comment = '#';
259  int numedges, numfaces, i, dim;
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 >> n;
267  stream >> numfaces;
268  stream >> numedges;
269 
270  i = 0;
271  while (i < n) {
272  std::getline(input, line);
273  if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
274  !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
275  std::stringstream iss(line);
276  std::vector<double> point;
277  point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
278  point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
279  boost::add_vertex(one_skeleton_OFF);
280  vertices.push_back(boost::add_vertex(one_skeleton));
281  cover.emplace_back();
282  i++;
283  }
284  }
285 
286  i = 0;
287  while (i < numfaces) {
288  std::getline(input, line);
289  if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
290  !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
291  std::vector<int> simplex;
292  std::stringstream iss(line);
293  simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>());
294  dim = simplex[0];
295  for (int j = 1; j <= dim; j++)
296  for (int k = j + 1; k <= dim; k++)
297  boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF);
298  i++;
299  }
300  }
301 
302  return input.is_open();
303  }
304 
305  // *******************************************************************************************************************
306  // Graphs.
307  // *******************************************************************************************************************
308 
309  public: // Set graph from file.
317  void set_graph_from_file(const std::string& graph_file_name) {
318  remove_edges(one_skeleton);
319  int neighb;
320  std::ifstream input(graph_file_name);
321  std::string line;
322  int source;
323  while (std::getline(input, line)) {
324  std::stringstream stream(line);
325  stream >> source;
326  while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
327  }
328  }
329 
330  public: // Set graph from OFF file.
335  remove_edges(one_skeleton);
336  if (num_edges(one_skeleton_OFF))
337  one_skeleton = one_skeleton_OFF;
338  else
339  std::cout << "No triangulation read in OFF file!" << std::endl;
340  }
341 
342  public: // Set graph from Rips complex.
349  template <typename Distance>
350  void set_graph_from_rips(double threshold, Distance distance) {
351  remove_edges(one_skeleton);
352  if (distances.size() == 0) compute_pairwise_distances(distance);
353  for (int i = 0; i < n; i++) {
354  for (int j = i + 1; j < n; j++) {
355  if (distances[i][j] <= threshold) {
356  boost::add_edge(vertices[i], vertices[j], one_skeleton);
357  boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first,
358  distances[i][j]);
359  }
360  }
361  }
362  }
363 
364  public:
365  void set_graph_weights() {
366  Index_map index = boost::get(boost::vertex_index, one_skeleton);
367  Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
368  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
369  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
370  boost::put(weight, *ei,
371  distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]);
372  }
373 
374  public: // Pairwise distances.
377  template <typename Distance>
378  void compute_pairwise_distances(Distance ref_distance) {
379  double d;
380  std::vector<double> zeros(n);
381  for (int i = 0; i < n; i++) distances.push_back(zeros);
382  std::string distance = point_cloud_name + "_dist";
383  std::ifstream input(distance, std::ios::out | std::ios::binary);
384 
385  if (input.good()) {
386  if (verbose) std::cout << "Reading distances..." << std::endl;
387  for (int i = 0; i < n; i++) {
388  for (int j = i; j < n; j++) {
389  input.read((char*)&d, 8);
390  distances[i][j] = d;
391  distances[j][i] = d;
392  }
393  }
394  input.close();
395  } else {
396  if (verbose) std::cout << "Computing distances..." << std::endl;
397  input.close();
398  std::ofstream output(distance, std::ios::out | std::ios::binary);
399  for (int i = 0; i < n; i++) {
400  int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10;
401  if (state == 0 && verbose) std::cout << "\r" << state << "%" << std::flush;
402  for (int j = i; j < n; j++) {
403  double dis = ref_distance(point_cloud[i], point_cloud[j]);
404  distances[i][j] = dis;
405  distances[j][i] = dis;
406  output.write((char*)&dis, 8);
407  }
408  }
409  output.close();
410  if (verbose) std::cout << std::endl;
411  }
412  }
413 
414  public: // Automatic tuning of Rips complex.
424  template <typename Distance>
425  double set_graph_from_automatic_rips(Distance distance, int N = 100) {
426  int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant))));
427  m = std::min(m, n - 1);
428  double delta = 0;
429 
430  if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl;
431  if (verbose) std::cout << "Subsampling " << m << " points" << std::endl;
432 
433  if (distances.size() == 0) compute_pairwise_distances(distance);
434 
435  // This cannot be parallelized if thread_local is not defined
436  // thread_local is not defined for XCode < v.8
437  #if defined(GUDHI_USE_TBB) && defined(GUDHI_CAN_USE_CXX11_THREAD_LOCAL)
438  tbb::mutex deltamutex;
439  tbb::parallel_for(0, N, [&](int i){
440  std::vector<int> samples(m);
441  SampleWithoutReplacement(n, m, samples);
442  double hausdorff_dist = 0;
443  for (int j = 0; j < n; j++) {
444  double mj = distances[j][samples[0]];
445  for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
446  hausdorff_dist = std::max(hausdorff_dist, mj);
447  }
448  deltamutex.lock();
449  delta += hausdorff_dist / N;
450  deltamutex.unlock();
451  });
452  #else
453  for (int i = 0; i < N; i++) {
454  std::vector<int> samples(m);
455  SampleWithoutReplacement(n, m, samples);
456  double hausdorff_dist = 0;
457  for (int j = 0; j < n; j++) {
458  double mj = distances[j][samples[0]];
459  for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
460  hausdorff_dist = std::max(hausdorff_dist, mj);
461  }
462  delta += hausdorff_dist / N;
463  }
464  #endif
465 
466  if (verbose) std::cout << "delta = " << delta << std::endl;
467  set_graph_from_rips(delta, distance);
468  return delta;
469  }
470 
471  // *******************************************************************************************************************
472  // Functions.
473  // *******************************************************************************************************************
474 
475  public: // Set function from file.
481  void set_function_from_file(const std::string& func_file_name) {
482  int i = 0;
483  std::ifstream input(func_file_name);
484  std::string line;
485  double f;
486  while (std::getline(input, line)) {
487  std::stringstream stream(line);
488  stream >> f;
489  func.push_back(f);
490  i++;
491  }
492  functional_cover = true;
493  cover_name = func_file_name;
494  }
495 
496  public: // Set function from kth coordinate
503  for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]);
504  functional_cover = true;
505  cover_name = "coordinate " + std::to_string(k);
506  }
507 
508  public: // Set function from vector.
514  template <class InputRange>
515  void set_function_from_range(InputRange const& function) {
516  for (int i = 0; i < n; i++) func.push_back(function[i]);
517  functional_cover = true;
518  }
519 
520  // *******************************************************************************************************************
521  // Covers.
522  // *******************************************************************************************************************
523 
524  public: // Automatic tuning of resolution.
533  if (!functional_cover) {
534  std::cout << "Cover needs to come from the preimages of a function." << std::endl;
535  return 0;
536  }
537  if (type != "Nerve" && type != "GIC") {
538  std::cout << "Type of complex needs to be specified." << std::endl;
539  return 0;
540  }
541 
542  double reso = 0;
543  Index_map index = boost::get(boost::vertex_index, one_skeleton);
544 
545  if (type == "GIC") {
546  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
547  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
548  reso = std::max(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
549  func[index[boost::target(*ei, one_skeleton)]]));
550  if (verbose) std::cout << "resolution = " << reso << std::endl;
551  resolution_double = reso;
552  }
553 
554  if (type == "Nerve") {
555  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
556  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
557  reso = std::max(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
558  func[index[boost::target(*ei, one_skeleton)]]) /
559  gain);
560  if (verbose) std::cout << "resolution = " << reso << std::endl;
561  resolution_double = reso;
562  }
563 
564  return reso;
565  }
566 
567  public:
573  void set_resolution_with_interval_length(double reso) { resolution_double = reso; }
579  void set_resolution_with_interval_number(int reso) { resolution_int = reso; }
585  void set_gain(double g = 0.3) { gain = g; }
586 
587  public: // Set cover with preimages of function.
592  if (resolution_double == -1 && resolution_int == -1) {
593  std::cout << "Number and/or length of intervals not specified" << std::endl;
594  return;
595  }
596  if (gain == -1) {
597  std::cout << "Gain not specified" << std::endl;
598  return;
599  }
600 
601  // Read function values and compute min and max
602  double minf = std::numeric_limits<float>::max();
603  double maxf = std::numeric_limits<float>::lowest();
604  for (int i = 0; i < n; i++) {
605  minf = std::min(minf, func[i]);
606  maxf = std::max(maxf, func[i]);
607  }
608  if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl;
609 
610  // Compute cover of im(f)
611  std::vector<std::pair<double, double> > intervals;
612  int res;
613 
614  if (resolution_double == -1) { // Case we use an integer for the number of intervals.
615  double incr = (maxf - minf) / resolution_int;
616  double x = minf;
617  double alpha = (incr * gain) / (2 - 2 * gain);
618  double y = minf + incr + alpha;
619  std::pair<double, double> interm(x, y);
620  intervals.push_back(interm);
621  for (int i = 1; i < resolution_int - 1; i++) {
622  x = minf + i * incr - alpha;
623  y = minf + (i + 1) * incr + alpha;
624  std::pair<double, double> inter(x, y);
625  intervals.push_back(inter);
626  }
627  x = minf + (resolution_int - 1) * incr - alpha;
628  y = maxf;
629  std::pair<double, double> interM(x, y);
630  intervals.push_back(interM);
631  res = intervals.size();
632  if (verbose) {
633  for (int i = 0; i < res; i++)
634  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
635  << std::endl;
636  }
637  } else {
638  if (resolution_int == -1) { // Case we use a double for the length of the intervals.
639  double x = minf;
640  double y = x + resolution_double;
641  while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
642  std::pair<double, double> inter(x, y);
643  intervals.push_back(inter);
644  x = y - gain * resolution_double;
645  y = x + resolution_double;
646  }
647  std::pair<double, double> interM(x, maxf);
648  intervals.push_back(interM);
649  res = intervals.size();
650  if (verbose) {
651  for (int i = 0; i < res; i++)
652  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
653  << std::endl;
654  }
655  } else { // Case we use an integer and a double for the length of the intervals.
656  double x = minf;
657  double y = x + resolution_double;
658  int count = 0;
659  while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
660  std::pair<double, double> inter(x, y);
661  intervals.push_back(inter);
662  count++;
663  x = y - gain * resolution_double;
664  y = x + resolution_double;
665  }
666  res = intervals.size();
667  if (verbose) {
668  for (int i = 0; i < res; i++)
669  std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
670  << std::endl;
671  }
672  }
673  }
674 
675  // Sort points according to function values
676  std::vector<int> points(n);
677  for (int i = 0; i < n; i++) points[i] = i;
678  std::sort(points.begin(), points.end(), [=](const int & p1, const int & p2){return (this->func[p1] < this->func[p2]);});
679 
680  int id = 0;
681  int pos = 0;
682  Index_map index = boost::get(boost::vertex_index, one_skeleton); // int maxc = -1;
683  std::map<int, std::vector<int> > preimages;
684  std::map<int, double> funcstd;
685 
686  if (verbose) std::cout << "Computing preimages..." << std::endl;
687  for (int i = 0; i < res; i++) {
688  // Find points in the preimage
689  std::pair<double, double> inter1 = intervals[i];
690  int tmp = pos;
691  double u, v;
692 
693  if (i != res - 1) {
694  if (i != 0) {
695  std::pair<double, double> inter3 = intervals[i - 1];
696  while (func[points[tmp]] < inter3.second && tmp != n) {
697  preimages[i].push_back(points[tmp]);
698  tmp++;
699  }
700  u = inter3.second;
701  } else {
702  u = inter1.first;
703  }
704 
705  std::pair<double, double> inter2 = intervals[i + 1];
706  while (func[points[tmp]] < inter2.first && tmp != n) {
707  preimages[i].push_back(points[tmp]);
708  tmp++;
709  }
710  v = inter2.first;
711  pos = tmp;
712  while (func[points[tmp]] < inter1.second && tmp != n) {
713  preimages[i].push_back(points[tmp]);
714  tmp++;
715  }
716 
717  } else {
718  std::pair<double, double> inter3 = intervals[i - 1];
719  while (func[points[tmp]] < inter3.second && tmp != n) {
720  preimages[i].push_back(points[tmp]);
721  tmp++;
722  }
723  while (tmp != n) {
724  preimages[i].push_back(points[tmp]);
725  tmp++;
726  }
727  u = inter3.second;
728  v = inter1.second;
729  }
730 
731  funcstd[i] = 0.5 * (u + v);
732  }
733 
734  #ifdef GUDHI_USE_TBB
735  if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl;
736  tbb::mutex covermutex, idmutex;
737  tbb::parallel_for(0, res, [&](int i){
738  // Compute connected components
739  Graph G = one_skeleton.create_subgraph();
740  int num = preimages[i].size();
741  std::vector<int> component(num);
742  for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
743  boost::connected_components(G, &component[0]);
744  int max = 0;
745 
746  // For each point in preimage
747  for (int j = 0; j < num; j++) {
748  // Update number of components in preimage
749  if (component[j] > max) max = component[j];
750 
751  // Identify component with Cantor polynomial N^2 -> N
752  int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
753 
754  // Update covers
755  covermutex.lock();
756  cover[preimages[i][j]].push_back(identifier);
757  cover_back[identifier].push_back(preimages[i][j]);
758  cover_fct[identifier] = i;
759  cover_std[identifier] = funcstd[i];
760  cover_color[identifier].second += func_color[preimages[i][j]];
761  cover_color[identifier].first += 1;
762  covermutex.unlock();
763  }
764 
765  // Maximal dimension is total number of connected components
766  idmutex.lock();
767  id += max + 1;
768  idmutex.unlock();
769  });
770  #else
771  if (verbose) std::cout << "Computing connected components..." << std::endl;
772  for (int i = 0; i < res; i++) {
773  // Compute connected components
774  Graph G = one_skeleton.create_subgraph();
775  int num = preimages[i].size();
776  std::vector<int> component(num);
777  for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
778  boost::connected_components(G, &component[0]);
779  int max = 0;
780 
781  // For each point in preimage
782  for (int j = 0; j < num; j++) {
783  // Update number of components in preimage
784  if (component[j] > max) max = component[j];
785 
786  // Identify component with Cantor polynomial N^2 -> N
787  int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
788 
789  // Update covers
790  cover[preimages[i][j]].push_back(identifier);
791  cover_back[identifier].push_back(preimages[i][j]);
792  cover_fct[identifier] = i;
793  cover_std[identifier] = funcstd[i];
794  cover_color[identifier].second += func_color[preimages[i][j]];
795  cover_color[identifier].first += 1;
796  }
797 
798  // Maximal dimension is total number of connected components
799  id += max + 1;
800  }
801  #endif
802 
803  maximal_dim = id - 1;
804  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
805  iit->second.second /= iit->second.first;
806  }
807 
808  public: // Set cover from file.
815  void set_cover_from_file(const std::string& cover_file_name) {
816  int i = 0;
817  int cov;
818  std::vector<int> cov_elts, cov_number;
819  std::ifstream input(cover_file_name);
820  std::string line;
821  while (std::getline(input, line)) {
822  cov_elts.clear();
823  std::stringstream stream(line);
824  while (stream >> cov) {
825  cov_elts.push_back(cov);
826  cov_number.push_back(cov);
827  cover_fct[cov] = cov;
828  cover_color[cov].second += func_color[i];
829  cover_color[cov].first++;
830  cover_back[cov].push_back(i);
831  }
832  cover[i] = cov_elts;
833  i++;
834  }
835 
836  std::sort(cov_number.begin(), cov_number.end());
837  std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
838  cov_number.resize(std::distance(cov_number.begin(), it));
839 
840  maximal_dim = cov_number.size() - 1;
841  for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
842  cover_name = cover_file_name;
843  }
844 
845  public: // Set cover from Voronoi
852  template <typename Distance>
853  void set_cover_from_Voronoi(Distance distance, int m = 100) {
854  voronoi_subsamples.resize(m);
855  SampleWithoutReplacement(n, m, voronoi_subsamples);
856  if (distances.size() == 0) compute_pairwise_distances(distance);
857  set_graph_weights();
858  Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
859  Index_map index = boost::get(boost::vertex_index, one_skeleton);
860  std::vector<double> mindist(n);
861  for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max();
862 
863  // Compute the geodesic distances to subsamples with Dijkstra
864  #ifdef GUDHI_USE_TBB
865  if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl;
866  tbb::mutex coverMutex; tbb::mutex mindistMutex;
867  tbb::parallel_for(0, m, [&](int i){
868  int seed = voronoi_subsamples[i];
869  std::vector<double> dmap(n);
870  boost::dijkstra_shortest_paths(
871  one_skeleton, vertices[seed],
872  boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
873 
874  coverMutex.lock(); mindistMutex.lock();
875  for (int j = 0; j < n; j++)
876  if (mindist[j] > dmap[j]) {
877  mindist[j] = dmap[j];
878  if (cover[j].size() == 0)
879  cover[j].push_back(i);
880  else
881  cover[j][0] = i;
882  }
883  coverMutex.unlock(); mindistMutex.unlock();
884  });
885  #else
886  for (int i = 0; i < m; i++) {
887  if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
888  int seed = voronoi_subsamples[i];
889  std::vector<double> dmap(n);
890  boost::dijkstra_shortest_paths(
891  one_skeleton, vertices[seed],
892  boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
893 
894  for (int j = 0; j < n; j++)
895  if (mindist[j] > dmap[j]) {
896  mindist[j] = dmap[j];
897  if (cover[j].size() == 0)
898  cover[j].push_back(i);
899  else
900  cover[j][0] = i;
901  }
902  }
903  #endif
904 
905  for (int i = 0; i < n; i++) {
906  cover_back[cover[i][0]].push_back(i);
907  cover_color[cover[i][0]].second += func_color[i];
908  cover_color[cover[i][0]].first++;
909  }
910  for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first;
911  maximal_dim = m - 1;
912  cover_name = "Voronoi";
913  }
914 
915  public: // return subset of data corresponding to a node
922  const std::vector<int>& subpopulation(int c) { return cover_back[name2idinv[c]]; }
923 
924  // *******************************************************************************************************************
925  // Visualization.
926  // *******************************************************************************************************************
927 
928  public: // Set color from file.
935  void set_color_from_file(const std::string& color_file_name) {
936  int i = 0;
937  std::ifstream input(color_file_name);
938  std::string line;
939  double f;
940  while (std::getline(input, line)) {
941  std::stringstream stream(line);
942  stream >> f;
943  func_color.push_back(f);
944  i++;
945  }
946  color_name = color_file_name;
947  }
948 
949  public: // Set color from kth coordinate
955  void set_color_from_coordinate(int k = 0) {
956  for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]);
957  color_name = "coordinate ";
958  color_name.append(std::to_string(k));
959  }
960 
961  public: // Set color from vector.
967  void set_color_from_vector(std::vector<double> color) {
968  for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
969  }
970 
971  public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
976  void plot_DOT() {
977  std::string mapp = point_cloud_name + "_sc.dot";
978  std::ofstream graphic(mapp);
979 
980  double maxv = std::numeric_limits<double>::lowest();
981  double minv = std::numeric_limits<double>::max();
982  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
983  maxv = std::max(maxv, iit->second.second);
984  minv = std::min(minv, iit->second.second);
985  }
986 
987  int k = 0;
988  std::vector<int> nodes;
989  nodes.clear();
990 
991  graphic << "graph GIC {" << std::endl;
992  int id = 0;
993  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
994  if (iit->second.first > mask) {
995  nodes.push_back(iit->first);
996  name2id[iit->first] = id;
997  name2idinv[id] = iit->first;
998  id++;
999  graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
1000  << ":" << iit->second.first << "\" style=filled fillcolor=\""
1001  << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl;
1002  k++;
1003  }
1004  }
1005  int ke = 0;
1006  int num_simplices = simplices.size();
1007  for (int i = 0; i < num_simplices; i++)
1008  if (simplices[i].size() == 2) {
1009  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
1010  graphic << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];"
1011  << std::endl;
1012  ke++;
1013  }
1014  }
1015  graphic << "}";
1016  graphic.close();
1017  std::cout << mapp << " file generated. It can be visualized with e.g. neato." << std::endl;
1018  }
1019 
1020  public: // Create a .txt file that can be compiled with KeplerMapper.
1024  void write_info() {
1025  int num_simplices = simplices.size();
1026  int num_edges = 0;
1027  std::string mapp = point_cloud_name + "_sc.txt";
1028  std::ofstream graphic(mapp);
1029 
1030  for (int i = 0; i < num_simplices; i++)
1031  if (simplices[i].size() == 2)
1032  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++;
1033 
1034  graphic << point_cloud_name << std::endl;
1035  graphic << cover_name << std::endl;
1036  graphic << color_name << std::endl;
1037  graphic << resolution_double << " " << gain << std::endl;
1038  graphic << cover_color.size() << " " << num_edges << std::endl;
1039 
1040  int id = 0;
1041  for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1042  graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl;
1043  name2id[iit->first] = id;
1044  name2idinv[id] = iit->first;
1045  id++;
1046  }
1047 
1048  for (int i = 0; i < num_simplices; i++)
1049  if (simplices[i].size() == 2)
1050  if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
1051  graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl;
1052  graphic.close();
1053  std::cout << mapp
1054  << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
1055  << std::endl;
1056  }
1057 
1058  public: // Create a .off file that can be visualized (e.g. with Geomview).
1063  void plot_OFF() {
1064  assert(cover_name == "Voronoi");
1065 
1066  int m = voronoi_subsamples.size();
1067  int numedges = 0;
1068  int numfaces = 0;
1069  std::vector<std::vector<int> > edges, faces;
1070  int numsimplices = simplices.size();
1071 
1072  std::string mapp = point_cloud_name + "_sc.off";
1073  std::ofstream graphic(mapp);
1074 
1075  graphic << "OFF" << std::endl;
1076  for (int i = 0; i < numsimplices; i++) {
1077  if (simplices[i].size() == 2) {
1078  numedges++;
1079  edges.push_back(simplices[i]);
1080  }
1081  if (simplices[i].size() == 3) {
1082  numfaces++;
1083  faces.push_back(simplices[i]);
1084  }
1085  }
1086  graphic << m << " " << numedges + numfaces << std::endl;
1087  for (int i = 0; i < m; i++) {
1088  if (data_dimension <= 3) {
1089  for (int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1090  for (int j = data_dimension; j < 3; j++) graphic << 0 << " ";
1091  graphic << std::endl;
1092  } else {
1093  for (int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1094  }
1095  }
1096  for (int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl;
1097  for (int i = 0; i < numfaces; i++)
1098  graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
1099  graphic.close();
1100  std::cout << mapp << " generated. It can be visualized with e.g. geomview." << std::endl;
1101  }
1102 
1103  // *******************************************************************************************************************
1104  // Extended Persistence Diagrams.
1105  // *******************************************************************************************************************
1106 
1107  public:
1111  void compute_PD() {
1112  Simplex_tree st;
1113 
1114  // Compute max and min
1115  double maxf = std::numeric_limits<double>::lowest();
1116  double minf = std::numeric_limits<double>::max();
1117  for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1118  maxf = std::max(maxf, it->second);
1119  minf = std::min(minf, it->second);
1120  }
1121 
1122  // Build filtration
1123  for (auto const& simplex : simplices) {
1124  std::vector<int> splx = simplex; splx.push_back(-2);
1125  st.insert_simplex_and_subfaces(splx, -3);
1126  }
1127 
1128  for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1129  int vertex = it->first; float val = it->second;
1130  int vert[] = {vertex}; int edge[] = {vertex, -2};
1131  st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf));
1132  st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf));
1133  }
1135 
1136  // Compute PD
1139 
1140  // Output PD
1141  int max_dim = st.dimension();
1142  for (int i = 0; i < max_dim; i++) {
1143  std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i);
1144  int num_bars = bars.size(); if(i == 0) num_bars -= 1;
1145  if(verbose) std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
1146  for (int j = 0; j < num_bars; j++) {
1147  double birth = bars[j].first;
1148  double death = bars[j].second;
1149  if (i == 0 && std::isinf(death)) continue;
1150  if (birth < 0)
1151  birth = minf + (birth + 2) * (maxf - minf);
1152  else
1153  birth = minf + (2 - birth) * (maxf - minf);
1154  if (death < 0)
1155  death = minf + (death + 2) * (maxf - minf);
1156  else
1157  death = minf + (2 - death) * (maxf - minf);
1158  PD.push_back(std::pair<double, double>(birth, death));
1159  if (verbose) std::cout << " [" << birth << ", " << death << "]" << std::endl;
1160  }
1161  }
1162  }
1163 
1164  public:
1170  void compute_distribution(unsigned int N = 100) {
1171  unsigned int sz = distribution.size();
1172  if (sz >= N) {
1173  std::cout << "Already done!" << std::endl;
1174  } else {
1175  for (unsigned int i = 0; i < N - sz; i++) {
1176  if (verbose) std::cout << "Computing " << i << "th bootstrap, bottleneck distance = ";
1177 
1178  Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true;
1179 
1180  std::vector<int> boot(this->n);
1181  for (int j = 0; j < this->n; j++) {
1182  double u = GetUniform();
1183  int id = std::floor(u * (this->n)); boot[j] = id;
1184  Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]);
1185  boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
1186  }
1187  Cboot.set_color_from_vector(Cboot.func);
1188 
1189  for (int j = 0; j < n; j++) {
1190  std::vector<double> dist(n);
1191  for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]];
1192  Cboot.distances.push_back(dist);
1193  }
1194 
1196  Cboot.set_automatic_resolution();
1197  Cboot.set_gain();
1198  Cboot.set_cover_from_function();
1199  Cboot.find_simplices();
1200  Cboot.compute_PD();
1201  double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD);
1202  if (verbose) std::cout << db << std::endl;
1203  distribution.push_back(db);
1204  }
1205 
1206  std::sort(distribution.begin(), distribution.end());
1207  }
1208  }
1209 
1210  public:
1217  unsigned int N = distribution.size();
1218  return distribution[std::floor(alpha * N)];
1219  }
1220 
1221  public:
1228  unsigned int N = distribution.size();
1229  double level = 1;
1230  for (unsigned int i = 0; i < N; i++)
1231  if (distribution[i] > d){ level = i * 1.0 / N; break; }
1232  return level;
1233  }
1234 
1235  public:
1240  double compute_p_value() {
1241  double distancemin = -std::numeric_limits<double>::lowest();
1242  int N = PD.size();
1243  for (int i = 0; i < N; i++) distancemin = std::min(distancemin, 0.5 * (PD[i].second - PD[i].first));
1244  double p_value = 1 - compute_confidence_level_from_distance(distancemin);
1245  if (verbose) std::cout << "p value = " << p_value << std::endl;
1246  return p_value;
1247  }
1248 
1249  // *******************************************************************************************************************
1250  // Computation of simplices.
1251  // *******************************************************************************************************************
1252 
1253  public:
1259  template <typename SimplicialComplex>
1260  void create_complex(SimplicialComplex& complex) {
1261  unsigned int dimension = 0;
1262  for (auto const& simplex : simplices) {
1263  int numvert = simplex.size();
1264  double filt = std::numeric_limits<double>::lowest();
1265  for (int i = 0; i < numvert; i++) filt = std::max(cover_color[simplex[i]].second, filt);
1266  complex.insert_simplex_and_subfaces(simplex, filt);
1267  if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
1268  }
1269  }
1270 
1271  public:
1275  if (type != "Nerve" && type != "GIC") {
1276  std::cout << "Type of complex needs to be specified." << std::endl;
1277  return;
1278  }
1279 
1280  if (type == "Nerve") {
1281  for(int i = 0; i < n; i++) simplices.push_back(cover[i]);
1282  std::sort(simplices.begin(), simplices.end());
1283  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1284  simplices.resize(std::distance(simplices.begin(), it));
1285  }
1286 
1287  if (type == "GIC") {
1288  Index_map index = boost::get(boost::vertex_index, one_skeleton);
1289 
1290  if (functional_cover) {
1291  // Computes the simplices in the GIC by looking at all the edges of the graph and adding the
1292  // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals.
1293 
1294  if (gain >= 0.5)
1295  throw std::invalid_argument(
1296  "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
1297 
1298  // Loop on all edges.
1299  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1300  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
1301  int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
1302  for (int i = 0; i < nums; i++) {
1303  int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
1304  int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
1305  for (int j = 0; j < numt; j++) {
1306  int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
1307  if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
1308  std::vector<int> edge(2);
1309  edge[0] = std::min(vs, vt);
1310  edge[1] = std::max(vs, vt);
1311  simplices.push_back(edge);
1312  goto afterLoop;
1313  }
1314  }
1315  }
1316  afterLoop:;
1317  }
1318  std::sort(simplices.begin(), simplices.end());
1319  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1320  simplices.resize(std::distance(simplices.begin(), it));
1321 
1322  } else {
1323  // Find edges to keep
1324  Simplex_tree st;
1325  boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1326  for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
1327  if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
1328  cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
1329  std::vector<int> edge(2);
1330  edge[0] = index[boost::source(*ei, one_skeleton)];
1331  edge[1] = index[boost::target(*ei, one_skeleton)];
1332  st.insert_simplex_and_subfaces(edge);
1333  }
1334 
1335  // st.insert_graph(one_skeleton);
1336 
1337  // Build the Simplex Tree corresponding to the graph
1338  st.expansion(maximal_dim);
1339 
1340  // Find simplices of GIC
1341  simplices.clear();
1342  for (auto simplex : st.complex_simplex_range()) {
1343  if (!st.has_children(simplex)) {
1344  std::vector<int> simplx;
1345  for (auto vertex : st.simplex_vertex_range(simplex)) {
1346  unsigned int sz = cover[vertex].size();
1347  for (unsigned int i = 0; i < sz; i++) {
1348  simplx.push_back(cover[vertex][i]);
1349  }
1350  }
1351  std::sort(simplx.begin(), simplx.end());
1352  std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
1353  simplx.resize(std::distance(simplx.begin(), it));
1354  simplices.push_back(simplx);
1355  }
1356  }
1357  std::sort(simplices.begin(), simplices.end());
1358  std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1359  simplices.resize(std::distance(simplices.begin(), it));
1360  }
1361  }
1362  }
1363 };
1364 
1365 } // namespace cover_complex
1366 
1367 } // namespace Gudhi
1368 
1369 #endif // GIC_H_
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:168
void set_color_from_vector(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:967
void plot_DOT()
Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial comple...
Definition: GIC.h:976
bool read_point_cloud(const std::string &off_file_name)
Reads and stores the input point cloud.
Definition: GIC.h:234
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1037
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:226
void set_function_from_range(InputRange const &function)
Creates the function f from a vector stored in memory.
Definition: GIC.h:515
void set_graph_from_file(const std::string &graph_file_name)
Creates a graph G from a file containing the edges.
Definition: GIC.h:317
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:64
const std::vector< int > & subpopulation(int c)
Returns the data subset corresponding to a specific node of the created complex.
Definition: GIC.h:922
Simplex Tree data structure for representing simplicial complexes.
Definition: Simplex_tree.h:72
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:683
double compute_distance_from_confidence_level(double alpha)
Computes the bottleneck distance threshold corresponding to a specific confidence level...
Definition: GIC.h:1216
void set_type(const std::string &t)
Specifies whether the type of the output simplicial complex.
Definition: GIC.h:195
void find_simplices()
Computes the simplices of the simplicial complex.
Definition: GIC.h:1274
void write_info()
Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e...
Definition: GIC.h:1024
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:46
void set_function_from_coordinate(int k)
Creates the function f from the k-th coordinate of the point cloud P.
Definition: GIC.h:502
Definition: SimplicialComplexForAlpha.h:26
Rips complex data structure.
Definition: Rips_complex.h:57
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:955
void create_complex(SimplicialComplex &complex)
Creates the simplicial complex.
Definition: GIC.h:1260
void set_graph_from_OFF()
Creates a graph G from the triangulation given by the input .OFF file.
Definition: GIC.h:334
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:481
void set_subsampling(double constant, double power)
Sets the constants used to subsample the data set. These constants are explained in ...
Definition: GIC.h:213
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:853
void set_verbose(bool verb=false)
Specifies whether the program should display information or not.
Definition: GIC.h:203
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:1252
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:703
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:1240
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:517
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:32
Complex_simplex_range complex_simplex_range()
Returns a range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:214
void assign_filtration(Simplex_handle sh, Filtration_value fv)
Sets the filtration value of a simplex.
Definition: Simplex_tree.h:421
Simplex_vertex_range simplex_vertex_range(Simplex_handle sh)
Returns a range over the vertices of a simplex.
Definition: Simplex_tree.h:261
Global distance functions.
void set_graph_from_rips(double threshold, Distance distance)
Creates a graph G from a Rips complex.
Definition: GIC.h:350
bool has_children(SimplexHandle sh) const
Returns true if the node in the simplex tree pointed by sh has children.
Definition: Simplex_tree.h:504
void compute_distribution(unsigned int N=100)
Computes bootstrapped distances distribution.
Definition: GIC.h:1170
void plot_OFF()
Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC...
Definition: GIC.h:1063
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:111
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Definition: Simplex_tree.h:476
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:425
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:935
double compute_confidence_level_from_distance(double d)
Computes the confidence level of a specific bottleneck distance threshold.
Definition: GIC.h:1227
void set_resolution_with_interval_number(int reso)
Sets a number of intervals from a value stored in memory.
Definition: GIC.h:579
void set_resolution_with_interval_length(double reso)
Sets a length of intervals from a value stored in memory.
Definition: GIC.h:573
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:815
void set_cover_from_function()
Creates a cover C from the preimages of the function f.
Definition: GIC.h:591
Options::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Simplex_tree.h:79
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:184
Cover complex data structure.
Definition: GIC.h:96
void compute_PD()
Computes the extended persistence diagram of the complex.
Definition: GIC.h:1111
double set_automatic_resolution()
Computes the optimal length of intervals (i.e. the smallest interval length avoiding discretization a...
Definition: GIC.h:532
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:585
GUDHI  Version 2.2.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Thu Jun 14 2018 15:00:54 for GUDHI by Doxygen 1.8.13