GIC.h
1/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
2 * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
3 * Author: Mathieu Carriere
4 *
5 * Copyright (C) 2017 Inria
6 *
7 * Modification(s):
8 * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL
9 * - YYYY/MM Author: Description of the modification
10 */
11
12#ifndef GIC_H_
13#define GIC_H_
14
15#ifdef GUDHI_USE_TBB
16#include <tbb/parallel_for.h>
17#include <mutex>
18#endif
19
20#include <gudhi/Debug_utils.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
51namespace Gudhi {
52
53namespace cover_complex {
54
58using Persistence_diagram = std::vector<std::pair<double, double> >;
59using Graph = boost::subgraph<
60 boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
61 boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
62using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
63using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
64using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
65
86template <typename Point>
88 private:
89 bool verbose = false; // whether to display information.
90 std::string type; // Nerve or GIC
91
92 std::vector<Point> point_cloud; // input point cloud.
93 std::vector<std::vector<double> > distances; // all pairwise distances.
94 int maximal_dim; // maximal dimension of output simplicial complex.
95 int data_dimension; // dimension of input data.
96 int n; // number of points.
97
98 std::vector<double> func; // function used to compute the output simplicial complex.
99 std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex.
100 bool functional_cover = false; // whether we use a cover with preimages of a function or not.
101
102 Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists).
103 Graph one_skeleton; // one-skeleton used to compute the connected components.
104 std::vector<Vertex_t> vertices; // vertices of one_skeleton.
105
106 std::vector<std::vector<int> > simplices; // simplices of output simplicial complex.
107 std::vector<int> voronoi_subsamples; // Voronoi germs (in case of Voronoi cover).
108
109 Persistence_diagram PD;
110 std::vector<double> distribution;
111
112 std::vector<std::vector<int> >
113 cover; // function associating to each data point the vector of cover elements to which it belongs.
114 std::map<int, std::vector<int> >
115 cover_back; // inverse of cover, in order to get the data points associated to a specific cover element.
116 std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence
117 // diagram of the output simplicial complex.
118 std::map<int, int>
119 cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not.
120 std::map<int, std::pair<int, double> >
121 cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex.
122
123 int resolution_int = -1;
124 double resolution_double = -1;
125 double gain = -1;
126 double rate_constant = 10; // Constant in the subsampling.
127 double rate_power = 0.001; // Power in the subsampling.
128 int mask = 0; // Ignore nodes containing less than mask points.
129
130 std::map<int, int> name2id, name2idinv;
131
132 std::string cover_name;
133 std::string point_cloud_name;
134 std::string color_name;
135
136 // Remove all edges of a graph.
137 void remove_edges(Graph& G) {
138 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
139 for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
140 }
141
142 // Find random number in [0,1].
143 double GetUniform() {
144 thread_local std::default_random_engine re;
145 std::uniform_real_distribution<double> Dist(0, 1);
146 return Dist(re);
147 }
148
149 // Subsample points.
150 void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector<int>& samples) {
151 int t = 0;
152 int m = 0;
153 double u;
154 while (m < sampleSize) {
155 u = GetUniform();
156 if ((populationSize - t) * u >= sampleSize - m) {
157 t++;
158 } else {
159 samples[m] = t;
160 t++;
161 m++;
162 }
163 }
164 }
165
166 // *******************************************************************************************************************
167 // Utils.
168 // *******************************************************************************************************************
169
170 public:
176 void set_type(const std::string& t) { type = t; }
177
178 public:
184 void set_verbose(bool verb = false) { verbose = verb; }
185
186 public:
194 void set_subsampling(double constant, double power) {
195 rate_constant = constant;
196 rate_power = power;
197 }
198
199 public:
207 void set_mask(int nodemask) { mask = nodemask; }
208
209 public:
210
211
217 void set_point_cloud_from_range(const std::vector<std::vector<double> > & point_cloud) {
218 n = point_cloud.size(); data_dimension = point_cloud[0].size();
219 point_cloud_name = "matrix"; cover.resize(n);
220 for(int i = 0; i < n; i++){
221 boost::add_vertex(one_skeleton_OFF);
222 vertices.push_back(boost::add_vertex(one_skeleton));
223 }
224 this->point_cloud = point_cloud;
225 }
226
232 bool read_point_cloud(const std::string& off_file_name) {
233 point_cloud_name = off_file_name;
234 std::ifstream input(off_file_name);
235 std::string line;
236
237 char comment = '#';
238 while (comment == '#') {
239 std::getline(input, line);
240 if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
241 comment = line[line.find_first_not_of(' ')];
242 }
243 if (strcmp((char*)line.c_str(), "nOFF") == 0) {
244 comment = '#';
245 while (comment == '#') {
246 std::getline(input, line);
247 if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
248 comment = line[line.find_first_not_of(' ')];
249 }
250 std::stringstream stream(line);
251 stream >> data_dimension;
252 } else {
253 data_dimension = 3;
254 }
255
256 comment = '#';
257 int numedges, numfaces, i, dim;
258 while (comment == '#') {
259 std::getline(input, line);
260 if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
261 comment = line[line.find_first_not_of(' ')];
262 }
263 std::stringstream stream(line);
264 stream >> n;
265 stream >> numfaces;
266 stream >> numedges;
267
268 i = 0;
269 while (i < n) {
270 std::getline(input, line);
271 if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
272 !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
273 std::stringstream iss(line);
274 std::vector<double> point;
275 point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
276 point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
277 boost::add_vertex(one_skeleton_OFF);
278 vertices.push_back(boost::add_vertex(one_skeleton));
279 cover.emplace_back();
280 i++;
281 }
282 }
283
284 i = 0;
285 while (i < numfaces) {
286 std::getline(input, line);
287 if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
288 !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
289 std::vector<int> simplex;
290 std::stringstream iss(line);
291 simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>());
292 dim = simplex[0];
293 for (int j = 1; j <= dim; j++)
294 for (int k = j + 1; k <= dim; k++)
295 boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF);
296 i++;
297 }
298 }
299
300 return input.is_open();
301 }
302
303 // *******************************************************************************************************************
304 // Graphs.
305 // *******************************************************************************************************************
306
307 public: // Set graph from file.
315 void set_graph_from_file(const std::string& graph_file_name) {
316 remove_edges(one_skeleton);
317 int neighb;
318 std::ifstream input(graph_file_name);
319 std::string line;
320 int source;
321 while (std::getline(input, line)) {
322 std::stringstream stream(line);
323 stream >> source;
324 while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
325 }
326 }
327
328 public: // Set graph from OFF file.
333 remove_edges(one_skeleton);
334 if (num_edges(one_skeleton_OFF))
335 one_skeleton = one_skeleton_OFF;
336 else
337 std::cerr << "No triangulation read in OFF file!" << std::endl;
338 }
339
340 public: // Set graph from Rips complex.
347 template <typename Distance>
348 void set_graph_from_rips(double threshold, Distance distance) {
349 remove_edges(one_skeleton);
350 if (distances.size() == 0) compute_pairwise_distances(distance);
351 for (int i = 0; i < n; i++) {
352 for (int j = i + 1; j < n; j++) {
353 if (distances[i][j] <= threshold) {
354 boost::add_edge(vertices[i], vertices[j], one_skeleton);
355 boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first,
356 distances[i][j]);
357 }
358 }
359 }
360 }
361
362 public:
363 void set_graph_weights() {
364 Index_map index = boost::get(boost::vertex_index, one_skeleton);
365 Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
366 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
367 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
368 boost::put(weight, *ei,
369 distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]);
370 }
371
372 public:
378 void set_distances_from_range(const std::vector<std::vector<double> > & distance_matrix) {
379 n = distance_matrix.size(); data_dimension = 0; point_cloud_name = "matrix";
380 cover.resize(n); point_cloud.resize(n);
381 for(int i = 0; i < n; i++){
382 boost::add_vertex(one_skeleton_OFF);
383 vertices.push_back(boost::add_vertex(one_skeleton));
384 }
385 distances = distance_matrix;
386 }
387
388 public: // Pairwise distances.
391 template <typename Distance>
392 void compute_pairwise_distances(Distance ref_distance) {
393 double d;
394 std::vector<double> zeros(n);
395 for (int i = 0; i < n; i++) distances.push_back(zeros);
396 std::string distance = point_cloud_name + "_dist";
397 std::ifstream input(distance, std::ios::out | std::ios::binary);
398
399 if (input.good()) {
400 if (verbose) std::clog << "Reading distances..." << std::endl;
401 for (int i = 0; i < n; i++) {
402 for (int j = i; j < n; j++) {
403 input.read((char*)&d, 8);
404 distances[i][j] = d;
405 distances[j][i] = d;
406 }
407 }
408 input.close();
409 } else {
410 if (verbose) std::clog << "Computing distances..." << std::endl;
411 input.close();
412 std::ofstream output(distance, std::ios::out | std::ios::binary);
413 for (int i = 0; i < n; i++) {
414 int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10;
415 if (state == 0 && verbose) std::clog << "\r" << state << "%" << std::flush;
416 for (int j = i; j < n; j++) {
417 double dis = ref_distance(point_cloud[i], point_cloud[j]);
418 distances[i][j] = dis;
419 distances[j][i] = dis;
420 output.write((char*)&dis, 8);
421 }
422 }
423 output.close();
424 if (verbose) std::clog << std::endl;
425 }
426 }
427
428 public: // Automatic tuning of Rips complex.
438 template <typename Distance>
439 double set_graph_from_automatic_rips(Distance distance, int N = 100) {
440 int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant))));
441 m = (std::min)(m, n - 1);
442 double delta = 0;
443
444 if (verbose) std::clog << n << " points in R^" << data_dimension << std::endl;
445 if (verbose) std::clog << "Subsampling " << m << " points" << std::endl;
446
447 if (distances.size() == 0) compute_pairwise_distances(distance);
448
449 #ifdef GUDHI_USE_TBB
450 std::mutex deltamutex;
451 tbb::parallel_for(0, N, [&](int i){
452 std::vector<int> samples(m);
453 SampleWithoutReplacement(n, m, samples);
454 double hausdorff_dist = 0;
455 for (int j = 0; j < n; j++) {
456 double mj = distances[j][samples[0]];
457 for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
458 hausdorff_dist = (std::max)(hausdorff_dist, mj);
459 }
460 deltamutex.lock();
461 delta += hausdorff_dist / N;
462 deltamutex.unlock();
463 });
464 #else
465 for (int i = 0; i < N; i++) {
466 std::vector<int> samples(m);
467 SampleWithoutReplacement(n, m, samples);
468 double hausdorff_dist = 0;
469 for (int j = 0; j < n; j++) {
470 double mj = distances[j][samples[0]];
471 for (int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
472 hausdorff_dist = (std::max)(hausdorff_dist, mj);
473 }
474 delta += hausdorff_dist / N;
475 }
476 #endif
477
478 if (verbose) std::clog << "delta = " << delta << std::endl;
479 set_graph_from_rips(delta, distance);
480 return delta;
481 }
482
483 // *******************************************************************************************************************
484 // Functions.
485 // *******************************************************************************************************************
486
487 public: // Set function from file.
493 void set_function_from_file(const std::string& func_file_name) {
494 int i = 0;
495 std::ifstream input(func_file_name);
496 std::string line;
497 double f;
498 while (std::getline(input, line)) {
499 std::stringstream stream(line);
500 stream >> f;
501 func.push_back(f);
502 i++;
503 }
504 functional_cover = true;
505 cover_name = func_file_name;
506 }
507
508 public: // Set function from kth coordinate
515 if(point_cloud[0].size() > 0){
516 for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]);
517 functional_cover = true;
518 cover_name = "coordinate " + std::to_string(k);
519 }
520 else{
521 std::cerr << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl;
522 for (int i = 0; i < n; i++) func.push_back(0.0);
523 functional_cover = true;
524 cover_name = "null";
525 }
526 }
527
528 public: // Set function from vector.
534 template <class InputRange>
535 void set_function_from_range(InputRange const& function) {
536 for (int i = 0; i < n; i++) func.push_back(function[i]);
537 functional_cover = true;
538 }
539
540 // *******************************************************************************************************************
541 // Covers.
542 // *******************************************************************************************************************
543
544 public: // Automatic tuning of resolution.
553 if (!functional_cover) {
554 std::cerr << "Cover needs to come from the preimages of a function." << std::endl;
555 return 0;
556 }
557 if (type != "Nerve" && type != "GIC") {
558 std::cerr << "Type of complex needs to be specified." << std::endl;
559 return 0;
560 }
561
562 double reso = 0;
563 Index_map index = boost::get(boost::vertex_index, one_skeleton);
564
565 if (type == "GIC") {
566 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
567 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
568 reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
569 func[index[boost::target(*ei, one_skeleton)]]));
570 if (verbose) std::clog << "resolution = " << reso << std::endl;
571 resolution_double = reso;
572 }
573
574 if (type == "Nerve") {
575 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
576 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
577 reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
578 func[index[boost::target(*ei, one_skeleton)]]) /
579 gain);
580 if (verbose) std::clog << "resolution = " << reso << std::endl;
581 resolution_double = reso;
582 }
583
584 return reso;
585 }
586
587 public:
593 void set_resolution_with_interval_length(double reso) { resolution_double = reso; }
599 void set_resolution_with_interval_number(int reso) { resolution_int = reso; }
605 void set_gain(double g = 0.3) { gain = g; }
606
607 public: // Set cover with preimages of function.
612 if (resolution_double == -1 && resolution_int == -1) {
613 std::cerr << "Number and/or length of intervals not specified" << std::endl;
614 return;
615 }
616 if (gain == -1) {
617 std::cerr << "Gain not specified" << std::endl;
618 return;
619 }
620
621 // Read function values and compute min and max
622 double minf = (std::numeric_limits<float>::max)();
623 double maxf = std::numeric_limits<float>::lowest();
624 for (int i = 0; i < n; i++) {
625 minf = (std::min)(minf, func[i]);
626 maxf = (std::max)(maxf, func[i]);
627 }
628 if (verbose) std::clog << "Min function value = " << minf << " and Max function value = " << maxf << std::endl;
629
630 // Compute cover of im(f)
631 std::vector<std::pair<double, double> > intervals;
632 int res;
633
634 if (resolution_double == -1) { // Case we use an integer for the number of intervals.
635 double incr = (maxf - minf) / resolution_int;
636 double x = minf;
637 double alpha = (incr * gain) / (2 - 2 * gain);
638 double y = minf + incr + alpha;
639 std::pair<double, double> interm(x, y);
640 intervals.push_back(interm);
641 for (int i = 1; i < resolution_int - 1; i++) {
642 x = minf + i * incr - alpha;
643 y = minf + (i + 1) * incr + alpha;
644 std::pair<double, double> inter(x, y);
645 intervals.push_back(inter);
646 }
647 x = minf + (resolution_int - 1) * incr - alpha;
648 y = maxf;
649 std::pair<double, double> interM(x, y);
650 intervals.push_back(interM);
651 res = intervals.size();
652 if (verbose) {
653 for (int i = 0; i < res; i++)
654 std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
655 << std::endl;
656 }
657 } else {
658 if (resolution_int == -1) { // Case we use a double for the length of the intervals.
659 double x = minf;
660 double y = x + resolution_double;
661 while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
662 std::pair<double, double> inter(x, y);
663 intervals.push_back(inter);
664 x = y - gain * resolution_double;
665 y = x + resolution_double;
666 }
667 std::pair<double, double> interM(x, maxf);
668 intervals.push_back(interM);
669 res = intervals.size();
670 if (verbose) {
671 for (int i = 0; i < res; i++)
672 std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
673 << std::endl;
674 }
675 } else { // Case we use an integer and a double for the length of the intervals.
676 double x = minf;
677 double y = x + resolution_double;
678 int count = 0;
679 while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
680 std::pair<double, double> inter(x, y);
681 intervals.push_back(inter);
682 count++;
683 x = y - gain * resolution_double;
684 y = x + resolution_double;
685 }
686 res = intervals.size();
687 if (verbose) {
688 for (int i = 0; i < res; i++)
689 std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]"
690 << std::endl;
691 }
692 }
693 }
694
695 // Sort points according to function values
696 std::vector<int> points(n);
697 for (int i = 0; i < n; i++) points[i] = i;
698 std::sort(points.begin(), points.end(), [this](int p1, int p2){return (this->func[p1] < this->func[p2]);});
699
700 int id = 0;
701 int pos = 0;
702 Index_map index = boost::get(boost::vertex_index, one_skeleton); // int maxc = -1;
703 std::map<int, std::vector<int> > preimages;
704 std::map<int, double> funcstd;
705
706 if (verbose) std::clog << "Computing preimages..." << std::endl;
707 for (int i = 0; i < res; i++) {
708 // Find points in the preimage
709 std::pair<double, double> inter1 = intervals[i];
710 int tmp = pos;
711 double u, v;
712
713 if (i != res - 1) {
714 if (i != 0) {
715 std::pair<double, double> inter3 = intervals[i - 1];
716 while (func[points[tmp]] < inter3.second && tmp != n) {
717 preimages[i].push_back(points[tmp]);
718 tmp++;
719 }
720 u = inter3.second;
721 } else {
722 u = inter1.first;
723 }
724
725 std::pair<double, double> inter2 = intervals[i + 1];
726 while (func[points[tmp]] < inter2.first && tmp != n) {
727 preimages[i].push_back(points[tmp]);
728 tmp++;
729 }
730 v = inter2.first;
731 pos = tmp;
732 while (func[points[tmp]] < inter1.second && tmp != n) {
733 preimages[i].push_back(points[tmp]);
734 tmp++;
735 }
736
737 } else {
738 std::pair<double, double> inter3 = intervals[i - 1];
739 while (func[points[tmp]] < inter3.second && tmp != n) {
740 preimages[i].push_back(points[tmp]);
741 tmp++;
742 }
743 while (tmp != n) {
744 preimages[i].push_back(points[tmp]);
745 tmp++;
746 }
747 u = inter3.second;
748 v = inter1.second;
749 }
750
751 funcstd[i] = 0.5 * (u + v);
752 }
753
754 #ifdef GUDHI_USE_TBB
755 if (verbose) std::clog << "Computing connected components (parallelized)..." << std::endl;
756 std::mutex covermutex, idmutex;
757 tbb::parallel_for(0, res, [&](int i){
758 // Compute connected components
759 Graph G = one_skeleton.create_subgraph();
760 int num = preimages[i].size();
761 std::vector<int> component(num);
762 for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
763 boost::connected_components(G, &component[0]);
764 int max = 0;
765
766 // For each point in preimage
767 for (int j = 0; j < num; j++) {
768 // Update number of components in preimage
769 if (component[j] > max) max = component[j];
770
771 // Identify component with Cantor polynomial N^2 -> N
772 int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
773
774 // Update covers
775 covermutex.lock();
776 cover[preimages[i][j]].push_back(identifier);
777 cover_back[identifier].push_back(preimages[i][j]);
778 cover_fct[identifier] = i;
779 cover_std[identifier] = funcstd[i];
780 cover_color[identifier].second += func_color[preimages[i][j]];
781 cover_color[identifier].first += 1;
782 covermutex.unlock();
783 }
784
785 // Maximal dimension is total number of connected components
786 idmutex.lock();
787 id += max + 1;
788 idmutex.unlock();
789 });
790 #else
791 if (verbose) std::clog << "Computing connected components..." << std::endl;
792 for (int i = 0; i < res; i++) {
793 // Compute connected components
794 Graph G = one_skeleton.create_subgraph();
795 int num = preimages[i].size();
796 std::vector<int> component(num);
797 for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
798 boost::connected_components(G, &component[0]);
799 int max = 0;
800
801 // For each point in preimage
802 for (int j = 0; j < num; j++) {
803 // Update number of components in preimage
804 if (component[j] > max) max = component[j];
805
806 // Identify component with Cantor polynomial N^2 -> N
807 int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
808
809 // Update covers
810 cover[preimages[i][j]].push_back(identifier);
811 cover_back[identifier].push_back(preimages[i][j]);
812 cover_fct[identifier] = i;
813 cover_std[identifier] = funcstd[i];
814 cover_color[identifier].second += func_color[preimages[i][j]];
815 cover_color[identifier].first += 1;
816 }
817
818 // Maximal dimension is total number of connected components
819 id += max + 1;
820 }
821 #endif
822
823 maximal_dim = id - 1;
824 for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
825 iit->second.second /= iit->second.first;
826 }
827
828 public: // Set cover from file.
835 void set_cover_from_file(const std::string& cover_file_name) {
836 int i = 0;
837 int cov;
838 std::vector<int> cov_elts, cov_number;
839 std::ifstream input(cover_file_name);
840 std::string line;
841 while (std::getline(input, line)) {
842 cov_elts.clear();
843 std::stringstream stream(line);
844 while (stream >> cov) {
845 cov_elts.push_back(cov);
846 cov_number.push_back(cov);
847 cover_fct[cov] = cov;
848 cover_color[cov].second += func_color[i];
849 cover_color[cov].first++;
850 cover_back[cov].push_back(i);
851 }
852 cover[i] = cov_elts;
853 i++;
854 }
855
856 std::sort(cov_number.begin(), cov_number.end());
857 std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
858 cov_number.resize(std::distance(cov_number.begin(), it));
859
860 maximal_dim = cov_number.size() - 1;
861 for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
862 cover_name = cover_file_name;
863 }
864
865 public: // Set cover from Voronoi
872 template <typename Distance>
873 void set_cover_from_Voronoi(Distance distance, int m = 100) {
874 voronoi_subsamples.resize(m);
875 SampleWithoutReplacement(n, m, voronoi_subsamples);
876 if (distances.size() == 0) compute_pairwise_distances(distance);
877 set_graph_weights();
878 Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
879 Index_map index = boost::get(boost::vertex_index, one_skeleton);
880 std::vector<double> mindist(n);
881 for (int j = 0; j < n; j++) mindist[j] = (std::numeric_limits<double>::max)();
882
883 // Compute the geodesic distances to subsamples with Dijkstra
884 #ifdef GUDHI_USE_TBB
885 if (verbose) std::clog << "Computing geodesic distances (parallelized)..." << std::endl;
886 std::mutex coverMutex; std::mutex mindistMutex;
887 tbb::parallel_for(0, m, [&](int i){
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 coverMutex.lock(); mindistMutex.lock();
895 for (int j = 0; j < n; j++)
896 if (mindist[j] > dmap[j]) {
897 mindist[j] = dmap[j];
898 if (cover[j].size() == 0)
899 cover[j].push_back(i);
900 else
901 cover[j][0] = i;
902 }
903 coverMutex.unlock(); mindistMutex.unlock();
904 });
905 #else
906 for (int i = 0; i < m; i++) {
907 if (verbose) std::clog << "Computing geodesic distances to seed " << i << "..." << std::endl;
908 int seed = voronoi_subsamples[i];
909 std::vector<double> dmap(n);
910 boost::dijkstra_shortest_paths(
911 one_skeleton, vertices[seed],
912 boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
913
914 for (int j = 0; j < n; j++)
915 if (mindist[j] > dmap[j]) {
916 mindist[j] = dmap[j];
917 if (cover[j].size() == 0)
918 cover[j].push_back(i);
919 else
920 cover[j][0] = i;
921 }
922 }
923 #endif
924
925 for (int i = 0; i < n; i++) {
926 cover_back[cover[i][0]].push_back(i);
927 cover_color[cover[i][0]].second += func_color[i];
928 cover_color[cover[i][0]].first++;
929 }
930 for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first;
931 maximal_dim = m - 1;
932 cover_name = "Voronoi";
933 }
934
935 public: // return subset of data corresponding to a node
942 const std::vector<int>& subpopulation(int c) { return cover_back[name2idinv[c]]; }
943
944 // *******************************************************************************************************************
945 // Visualization.
946 // *******************************************************************************************************************
947
948 public: // Set color from file.
955 void set_color_from_file(const std::string& color_file_name) {
956 int i = 0;
957 std::ifstream input(color_file_name);
958 std::string line;
959 double f;
960 while (std::getline(input, line)) {
961 std::stringstream stream(line);
962 stream >> f;
963 func_color.push_back(f);
964 i++;
965 }
966 color_name = color_file_name;
967 }
968
969 public: // Set color from kth coordinate
976 if(point_cloud[0].size() > 0){
977 for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]);
978 color_name = "coordinate ";
979 color_name.append(std::to_string(k));
980 }
981 else{
982 std::cerr << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl;
983 for (int i = 0; i < n; i++) func.push_back(0.0);
984 functional_cover = true;
985 cover_name = "null";
986 }
987 }
988
989 public: // Set color from vector.
995 void set_color_from_range(std::vector<double> color) {
996 for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
997 }
998
999 public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
1004 void plot_DOT() {
1005 std::string mapp = point_cloud_name + "_sc.dot";
1006 std::ofstream graphic(mapp);
1007
1008 double maxv = std::numeric_limits<double>::lowest();
1009 double minv = (std::numeric_limits<double>::max)();
1010 for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1011 maxv = (std::max)(maxv, iit->second.second);
1012 minv = (std::min)(minv, iit->second.second);
1013 }
1014
1015 int k = 0;
1016 std::vector<int> nodes;
1017 nodes.clear();
1018
1019 graphic << "graph GIC {" << std::endl;
1020 int id = 0;
1021 for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1022 if (iit->second.first > mask) {
1023 nodes.push_back(iit->first);
1024 name2id[iit->first] = id;
1025 name2idinv[id] = iit->first;
1026 id++;
1027 graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
1028 << ":" << iit->second.first << "\" style=filled fillcolor=\""
1029 << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl;
1030 k++;
1031 }
1032 }
1033 int ke = 0;
1034 int num_simplices = simplices.size();
1035 for (int i = 0; i < num_simplices; i++)
1036 if (simplices[i].size() == 2) {
1037 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
1038 graphic << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];"
1039 << std::endl;
1040 ke++;
1041 }
1042 }
1043 graphic << "}";
1044 graphic.close();
1045 std::clog << mapp << " file generated. It can be visualized with e.g. neato." << std::endl;
1046 }
1047
1048 public: // Create a .txt file that can be compiled with KeplerMapper.
1052 void write_info() {
1053 int num_simplices = simplices.size();
1054 int num_edges = 0;
1055 std::string mapp = point_cloud_name + "_sc.txt";
1056 std::ofstream graphic(mapp);
1057
1058 for (int i = 0; i < num_simplices; i++)
1059 if (simplices[i].size() == 2)
1060 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++;
1061
1062 graphic << point_cloud_name << std::endl;
1063 graphic << cover_name << std::endl;
1064 graphic << color_name << std::endl;
1065 graphic << resolution_double << " " << gain << std::endl;
1066 graphic << cover_color.size() << " " << num_edges << std::endl;
1067
1068 int id = 0;
1069 for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1070 graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl;
1071 name2id[iit->first] = id;
1072 name2idinv[id] = iit->first;
1073 id++;
1074 }
1075
1076 for (int i = 0; i < num_simplices; i++)
1077 if (simplices[i].size() == 2)
1078 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
1079 graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl;
1080 graphic.close();
1081 std::clog << mapp
1082 << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
1083 << std::endl;
1084 }
1085
1086 public: // Create a .off file that can be visualized (e.g. with Geomview).
1091 void plot_OFF() {
1092 assert(cover_name == "Voronoi");
1093
1094 int m = voronoi_subsamples.size();
1095 int numedges = 0;
1096 int numfaces = 0;
1097 std::vector<std::vector<int> > edges, faces;
1098 int numsimplices = simplices.size();
1099
1100 std::string mapp = point_cloud_name + "_sc.off";
1101 std::ofstream graphic(mapp);
1102
1103 graphic << "OFF" << std::endl;
1104 for (int i = 0; i < numsimplices; i++) {
1105 if (simplices[i].size() == 2) {
1106 numedges++;
1107 edges.push_back(simplices[i]);
1108 }
1109 if (simplices[i].size() == 3) {
1110 numfaces++;
1111 faces.push_back(simplices[i]);
1112 }
1113 }
1114 graphic << m << " " << numedges + numfaces << std::endl;
1115 for (int i = 0; i < m; i++) {
1116 if (data_dimension <= 3) {
1117 for (int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1118 for (int j = data_dimension; j < 3; j++) graphic << 0 << " ";
1119 graphic << std::endl;
1120 } else {
1121 for (int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] << " ";
1122 }
1123 }
1124 for (int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl;
1125 for (int i = 0; i < numfaces; i++)
1126 graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
1127 graphic.close();
1128 std::clog << mapp << " generated. It can be visualized with e.g. geomview." << std::endl;
1129 }
1130
1131 // *******************************************************************************************************************
1132 // Extended Persistence Diagrams.
1133 // *******************************************************************************************************************
1134
1135 public:
1139 Persistence_diagram compute_PD() {
1140 Simplex_tree st;
1141
1142 // Compute max and min
1143 double maxf = std::numeric_limits<double>::lowest();
1144 double minf = (std::numeric_limits<double>::max)();
1145 for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1146 maxf = (std::max)(maxf, it->second);
1147 minf = (std::min)(minf, it->second);
1148 }
1149
1150 // Build filtration
1151 for (auto const& simplex : simplices) {
1152 std::vector<int> splx = simplex;
1153 splx.push_back(-2);
1154 st.insert_simplex_and_subfaces(splx, -3);
1155 }
1156
1157 for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1158 int vertex = it->first; float val = it->second;
1159 int vert[] = {vertex}; int edge[] = {vertex, -2};
1160 if(st.find(vert) != st.null_simplex()){
1161 st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf));
1162 st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf));
1163 }
1164 }
1166
1167 // Compute PD
1170
1171 // Output PD
1172 int max_dim = st.dimension();
1173 for (int i = 0; i < max_dim; i++) {
1174 std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i);
1175 int num_bars = bars.size(); if(i == 0) num_bars -= 1;
1176 if(verbose) std::clog << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
1177 for (int j = 0; j < num_bars; j++) {
1178 double birth = bars[j].first;
1179 double death = bars[j].second;
1180 if (i == 0 && std::isinf(death)) continue;
1181 if (birth < 0)
1182 birth = minf + (birth + 2) * (maxf - minf);
1183 else
1184 birth = minf + (2 - birth) * (maxf - minf);
1185 if (death < 0)
1186 death = minf + (death + 2) * (maxf - minf);
1187 else
1188 death = minf + (2 - death) * (maxf - minf);
1189 PD.push_back(std::pair<double, double>(birth, death));
1190 if (verbose) std::clog << " [" << birth << ", " << death << "]" << std::endl;
1191 }
1192 }
1193 return PD;
1194 }
1195
1196 public:
1202 void compute_distribution(unsigned int N = 100) {
1203 unsigned int sz = distribution.size();
1204 if (sz < N) {
1205 for (unsigned int i = 0; i < N - sz; i++) {
1206 if (verbose) std::clog << "Computing " << i << "th bootstrap, bottleneck distance = ";
1207
1208 Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true;
1209
1210 std::vector<int> boot(this->n);
1211 for (int j = 0; j < this->n; j++) {
1212 double u = GetUniform();
1213 int id = std::floor(u * (this->n)); boot[j] = id;
1214 Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]);
1215 boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
1216 }
1217 Cboot.set_color_from_range(Cboot.func);
1218
1219 for (int j = 0; j < n; j++) {
1220 std::vector<double> dist(n);
1221 for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]];
1222 Cboot.distances.push_back(dist);
1223 }
1224
1226 Cboot.set_gain();
1229 Cboot.find_simplices();
1230 Cboot.compute_PD();
1231 double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD);
1232 if (verbose) std::clog << db << std::endl;
1233 distribution.push_back(db);
1234 }
1235
1236 std::sort(distribution.begin(), distribution.end());
1237 }
1238 }
1239
1240 public:
1247 unsigned int N = distribution.size();
1248 double d = distribution[std::floor(alpha * N)];
1249 if (verbose) std::clog << "Distance corresponding to confidence " << alpha << " is " << d << std::endl;
1250 return d;
1251 }
1252
1253 public:
1260 unsigned int N = distribution.size();
1261 double level = 1;
1262 for (unsigned int i = 0; i < N; i++)
1263 if (distribution[i] >= d){ level = i * 1.0 / N; break; }
1264 if (verbose) std::clog << "Confidence level of distance " << d << " is " << level << std::endl;
1265 return level;
1266 }
1267
1268 public:
1274 double distancemin = (std::numeric_limits<double>::max)(); int N = PD.size();
1275 for (int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first));
1276 double p_value = 1 - compute_confidence_level_from_distance(distancemin);
1277 if (verbose) std::clog << "p value = " << p_value << std::endl;
1278 return p_value;
1279 }
1280
1281 // *******************************************************************************************************************
1282 // Computation of simplices.
1283 // *******************************************************************************************************************
1284
1285 public:
1291 template <typename SimplicialComplex>
1292 void create_complex(SimplicialComplex& complex) {
1293 unsigned int dimension = 0;
1294 for (auto const& simplex : simplices) {
1295 int numvert = simplex.size();
1296 double filt = std::numeric_limits<double>::lowest();
1297 for (int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt);
1298 complex.insert_simplex_and_subfaces(simplex, filt);
1299 if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
1300 }
1301 }
1302
1303 public:
1307 if (type != "Nerve" && type != "GIC") {
1308 std::cerr << "Type of complex needs to be specified." << std::endl;
1309 return;
1310 }
1311
1312 if (type == "Nerve") {
1313 for(int i = 0; i < n; i++) simplices.push_back(cover[i]);
1314 std::sort(simplices.begin(), simplices.end());
1315 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1316 simplices.resize(std::distance(simplices.begin(), it));
1317 }
1318
1319 if (type == "GIC") {
1320 Index_map index = boost::get(boost::vertex_index, one_skeleton);
1321
1322 if (functional_cover) {
1323 // Computes the simplices in the GIC by looking at all the edges of the graph and adding the
1324 // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals.
1325
1326 if (gain >= 0.5)
1327 throw std::invalid_argument(
1328 "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
1329
1330 // Loop on all edges.
1331 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1332 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
1333 int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
1334 for (int i = 0; i < nums; i++) {
1335 int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
1336 int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
1337 for (int j = 0; j < numt; j++) {
1338 int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
1339 if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
1340 std::vector<int> edge(2);
1341 edge[0] = (std::min)(vs, vt);
1342 edge[1] = (std::max)(vs, vt);
1343 simplices.push_back(edge);
1344 goto afterLoop;
1345 }
1346 }
1347 }
1348 afterLoop:;
1349 }
1350 std::sort(simplices.begin(), simplices.end());
1351 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1352 simplices.resize(std::distance(simplices.begin(), it));
1353
1354 } else {
1355 // Find edges to keep
1356 Simplex_tree st;
1357 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1358 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
1359 if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
1360 cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
1361 std::vector<int> edge(2);
1362 edge[0] = index[boost::source(*ei, one_skeleton)];
1363 edge[1] = index[boost::target(*ei, one_skeleton)];
1365 }
1366
1367 // st.insert_graph(one_skeleton);
1368
1369 // Build the Simplex Tree corresponding to the graph
1370 st.expansion(maximal_dim);
1371
1372 // Find simplices of GIC
1373 simplices.clear();
1374 for (auto simplex : st.complex_simplex_range()) {
1375 if (!st.has_children(simplex)) {
1376 std::vector<int> simplx;
1377 for (auto vertex : st.simplex_vertex_range(simplex)) {
1378 unsigned int sz = cover[vertex].size();
1379 for (unsigned int i = 0; i < sz; i++) {
1380 simplx.push_back(cover[vertex][i]);
1381 }
1382 }
1383 std::sort(simplx.begin(), simplx.end());
1384 std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
1385 simplx.resize(std::distance(simplx.begin(), it));
1386 simplices.push_back(simplx);
1387 }
1388 }
1389 std::sort(simplices.begin(), simplices.end());
1390 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1391 simplices.resize(std::distance(simplices.begin(), it));
1392 }
1393 }
1394 }
1395};
1396
1397} // namespace cover_complex
1398
1399} // namespace Gudhi
1400
1401#endif // GIC_H_
Compute the Euclidean distance between two Points given by a range of coordinates....
Definition: distance_functions.h:32
Options::Filtration_value Filtration_value
Type for the value of the filtration function.
Definition: Simplex_tree.h:86
void assign_filtration(Simplex_handle sh, Filtration_value fv)
Sets the filtration value of a simplex.
Definition: Simplex_tree.h:546
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:1386
Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const
Returns a range over the vertices of a simplex.
Definition: Simplex_tree.h:282
bool has_children(SimplexHandle sh) const
Returns true if the node in the simplex tree pointed by sh has children.
Definition: Simplex_tree.h:628
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:641
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1170
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Definition: Simplex_tree.h:600
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:807
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:556
Complex_simplex_range complex_simplex_range()
Returns a range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:237
Cover complex data structure.
Definition: GIC.h:87
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:493
double set_automatic_resolution()
Computes the optimal length of intervals (i.e. the smallest interval length avoiding discretization a...
Definition: GIC.h:552
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:873
void set_resolution_with_interval_number(int reso)
Sets a number of intervals from a value stored in memory.
Definition: GIC.h:599
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:207
Persistence_diagram compute_PD()
Computes the extended persistence diagram of the complex.
Definition: GIC.h:1139
double compute_distance_from_confidence_level(double alpha)
Computes the bottleneck distance threshold corresponding to a specific confidence level.
Definition: GIC.h:1246
void set_graph_from_rips(double threshold, Distance distance)
Creates a graph G from a Rips complex.
Definition: GIC.h:348
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:835
void set_graph_from_file(const std::string &graph_file_name)
Creates a graph G from a file containing the edges.
Definition: GIC.h:315
void create_complex(SimplicialComplex &complex)
Creates the simplicial complex.
Definition: GIC.h:1292
void set_type(const std::string &t)
Specifies whether the type of the output simplicial complex.
Definition: GIC.h:176
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:378
void find_simplices()
Computes the simplices of the simplicial complex.
Definition: GIC.h:1306
void set_function_from_range(InputRange const &function)
Creates the function f from a vector stored in memory.
Definition: GIC.h:535
void set_cover_from_function()
Creates a cover C from the preimages of the function f.
Definition: GIC.h:611
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:955
void write_info()
Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e....
Definition: GIC.h:1052
void plot_OFF()
Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC....
Definition: GIC.h:1091
void plot_DOT()
Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial comple...
Definition: GIC.h:1004
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:995
void set_subsampling(double constant, double power)
Sets the constants used to subsample the data set. These constants are explained in .
Definition: GIC.h:194
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:217
double set_graph_from_automatic_rips(Distance distance, int N=100)
Creates a graph G from a Rips complex whose threshold value is automatically tuned with subsampling—s...
Definition: GIC.h:439
void set_resolution_with_interval_length(double reso)
Sets a length of intervals from a value stored in memory.
Definition: GIC.h:593
const std::vector< int > & subpopulation(int c)
Returns the data subset corresponding to a specific node of the created complex.
Definition: GIC.h:942
void set_gain(double g=0.3)
Sets a gain from a value stored in memory (default value 0.3).
Definition: GIC.h:605
void set_function_from_coordinate(int k)
Creates the function f from the k-th coordinate of the point cloud P.
Definition: GIC.h:514
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:975
void compute_distribution(unsigned int N=100)
Computes bootstrapped distances distribution.
Definition: GIC.h:1202
double compute_confidence_level_from_distance(double d)
Computes the confidence level of a specific bottleneck distance threshold.
Definition: GIC.h:1259
void set_graph_from_OFF()
Creates a graph G from the triangulation given by the input .OFF file.
Definition: GIC.h:332
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:1273
bool read_point_cloud(const std::string &off_file_name)
Reads and stores the input point cloud from .(n)OFF file.
Definition: GIC.h:232
void set_verbose(bool verb=false)
Specifies whether the program should display information or not.
Definition: GIC.h:184
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:52
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:681
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:172
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:156
Rips complex data structure.
Definition: Rips_complex.h:45
Global distance functions.
Graph simplicial complex methods.
double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, double e=(std::numeric_limits< double >::min)())
Function to compute the Bottleneck distance between two persistence diagrams.
Definition: Bottleneck.h:116
This file includes common file reader for GUDHI.
Value type for a filtration function on a cell complex.
Definition: FiltrationValue.h:20