16#include <tbb/parallel_for.h>
20#if __has_include(<CGAL/version.h>)
21# define GUDHI_GIC_USE_CGAL 1
22# include <gudhi/Bottleneck.h>
23#elif __has_include(<hera/bottleneck.h>)
24# define GUDHI_GIC_USE_HERA 1
28# include <type_traits>
30using ssize_t = std::make_signed_t<std::size_t>;
32# include <hera/bottleneck.h>
35#include <gudhi/Debug_utils.h>
38#include <gudhi/Simplex_tree.h>
39#include <gudhi/Rips_complex.h>
40#include <gudhi/Points_off_io.h>
42#include <gudhi/Persistent_cohomology.h>
44#include <boost/config.hpp>
45#include <boost/graph/graph_traits.hpp>
46#include <boost/graph/adjacency_list.hpp>
47#include <boost/graph/connected_components.hpp>
48#include <boost/graph/dijkstra_shortest_paths.hpp>
49#include <boost/graph/subgraph.hpp>
50#include <boost/graph/graph_utility.hpp>
65namespace cover_complex {
70using Persistence_diagram = std::vector<std::pair<double, double> >;
71using Graph = boost::subgraph<
72 boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
73 boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
74using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
75using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
76using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
98template <
typename Po
int>
101 bool verbose =
false;
104 std::vector<Point> point_cloud;
105 std::vector<std::vector<double> > distances;
110 std::vector<double> func;
111 std::vector<double> func_color;
112 bool functional_cover =
false;
114 Graph one_skeleton_OFF;
116 std::vector<Vertex_t> vertices;
118 std::vector<std::vector<int> > simplices;
119 std::vector<int> voronoi_subsamples;
121 Persistence_diagram PD;
122 std::vector<double> distribution;
124 std::vector<std::vector<int> >
126 std::map<int, std::vector<int> >
128 std::map<int, double> cover_std;
132 std::map<int, std::pair<int, double> >
135 int resolution_int = -1;
136 double resolution_double = -1;
138 double rate_constant = 10;
139 double rate_power = 0.001;
142 std::map<int, int> name2id, name2idinv;
144 std::string cover_name;
145 std::string point_cloud_name;
146 std::string color_name;
149 void remove_edges(Graph& G) {
150 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
151 for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
155 double GetUniform() {
156 thread_local std::default_random_engine re;
157 std::uniform_real_distribution<double> Dist(0, 1);
162 void SampleWithoutReplacement(
int populationSize,
int sampleSize, std::vector<int>& samples) {
166 while (m < sampleSize) {
168 if ((populationSize - t) * u >= sampleSize - m) {
207 rate_constant = constant;
230 this->num_points = point_cloud.size(); data_dimension = point_cloud[0].size();
231 point_cloud_name =
"cloud"; cover.resize(this->num_points);
232 for(
int i = 0; i < this->num_points; i++){
233 boost::add_vertex(one_skeleton_OFF);
234 vertices.push_back(boost::add_vertex(one_skeleton));
236 this->point_cloud = point_cloud;
245 point_cloud_name = off_file_name;
246 std::ifstream input(off_file_name);
250 while (comment ==
'#') {
251 std::getline(input, line);
252 if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(
int))isspace))
253 comment = line[line.find_first_not_of(
' ')];
255 if (strcmp((
char*)line.c_str(),
"nOFF") == 0) {
257 while (comment ==
'#') {
258 std::getline(input, line);
259 if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(
int))isspace))
260 comment = line[line.find_first_not_of(
' ')];
262 std::stringstream stream(line);
263 stream >> data_dimension;
269 int numedges, numfaces, i, dim;
270 while (comment ==
'#') {
271 std::getline(input, line);
272 if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(
int))isspace))
273 comment = line[line.find_first_not_of(
' ')];
275 std::stringstream stream(line);
276 stream >> this->num_points;
281 while (i < this->num_points) {
282 std::getline(input, line);
283 if (!line.empty() && line[line.find_first_not_of(
' ')] !=
'#' &&
284 !all_of(line.begin(), line.end(), (int (*)(
int))isspace)) {
285 std::stringstream iss(line);
286 std::vector<double> point;
287 point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
288 point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
289 boost::add_vertex(one_skeleton_OFF);
290 vertices.push_back(boost::add_vertex(one_skeleton));
291 cover.emplace_back();
297 while (i < numfaces) {
298 std::getline(input, line);
299 if (!line.empty() && line[line.find_first_not_of(
' ')] !=
'#' &&
300 !all_of(line.begin(), line.end(), (int (*)(
int))isspace)) {
301 std::vector<int> simplex;
302 std::stringstream iss(line);
303 simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>());
305 for (
int j = 1; j <= dim; j++)
306 for (
int k = j + 1; k <= dim; k++)
307 boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF);
312 return input.is_open();
328 remove_edges(one_skeleton);
330 std::ifstream input(graph_file_name);
333 while (std::getline(input, line)) {
334 std::stringstream stream(line);
336 while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
345 remove_edges(one_skeleton);
346 if (num_edges(one_skeleton_OFF))
347 one_skeleton = one_skeleton_OFF;
349 std::cerr <<
"No triangulation read in OFF file!" << std::endl;
359 template <
typename Distance>
361 remove_edges(one_skeleton);
362 if (distances.size() == 0) compute_pairwise_distances(distance);
363 for (
int i = 0; i < this->num_points; i++) {
364 for (
int j = i + 1; j < this->num_points; j++) {
365 if (distances[i][j] <= threshold) {
366 boost::add_edge(vertices[i], vertices[j], one_skeleton);
367 boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first,
375 void set_graph_weights() {
376 Index_map index = boost::get(boost::vertex_index, one_skeleton);
377 Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
378 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
379 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
380 boost::put(weight, *ei,
381 distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]);
391 this->num_points = distance_matrix.size(); data_dimension = 0; point_cloud_name =
"matrix";
392 cover.resize(this->num_points); point_cloud.resize(this->num_points);
393 for(
int i = 0; i < this->num_points; i++){
394 boost::add_vertex(one_skeleton_OFF);
395 vertices.push_back(boost::add_vertex(one_skeleton));
397 distances = distance_matrix;
403 template <
typename Distance>
404 void compute_pairwise_distances(Distance ref_distance) {
405 std::vector<double> zeros(this->num_points);
406 for (
int i = 0; i < this->num_points; i++) distances.push_back(zeros);
407 if (verbose) std::clog <<
"Computing distances..." << std::endl;
408 for (
int i = 0; i < this->num_points; i++) {
409 int state = 100 * (i + 1) / this->num_points;
410 if (verbose && state % 10 == 0) std::clog <<
"\r" << state <<
"%" << std::flush;
411 for (
int j = i; j < this->num_points; j++) {
412 double dis = ref_distance(point_cloud[i], point_cloud[j]);
413 distances[i][j] = dis;
414 distances[j][i] = dis;
417 if (verbose) std::clog << std::endl;
430 template <
typename Distance>
432 int m = floor(this->num_points / std::exp((1 + rate_power) * std::log(std::log(this->num_points) / std::log(rate_constant))));
433 m = (std::min)(m, this->num_points - 1);
436 if (verbose) std::clog << this->num_points <<
" points in R^" << data_dimension << std::endl;
437 if (verbose) std::clog <<
"Subsampling " << m <<
" points" << std::endl;
439 if (distances.size() == 0) compute_pairwise_distances(distance);
442 std::mutex deltamutex;
443 tbb::parallel_for(0, N, [&](
int i){
444 std::vector<int> samples(m);
445 SampleWithoutReplacement(this->num_points, m, samples);
446 double hausdorff_dist = 0;
447 for (
int j = 0; j < this->num_points; j++) {
448 double mj = distances[j][samples[0]];
449 for (
int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
450 hausdorff_dist = (std::max)(hausdorff_dist, mj);
453 delta += hausdorff_dist / N;
457 for (
int i = 0; i < N; i++) {
458 std::vector<int> samples(m);
459 SampleWithoutReplacement(this->num_points, m, samples);
460 double hausdorff_dist = 0;
461 for (
int j = 0; j < this->num_points; j++) {
462 double mj = distances[j][samples[0]];
463 for (
int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
464 hausdorff_dist = (std::max)(hausdorff_dist, mj);
466 delta += hausdorff_dist / N;
470 if (verbose) std::clog <<
"delta = " << delta << std::endl;
487 std::ifstream input(func_file_name);
490 while (std::getline(input, line)) {
491 std::stringstream stream(line);
496 functional_cover =
true;
497 cover_name = func_file_name;
507 if(point_cloud[0].size() > 0){
508 for (
int i = 0; i < this->num_points; i++) func.push_back(point_cloud[i][k]);
509 functional_cover =
true;
510 cover_name =
"coordinate " + std::to_string(k);
513 std::cerr <<
"Only pairwise distances provided---cannot access " << k <<
"th coordinate; returning null vector instead" << std::endl;
514 for (
int i = 0; i < this->num_points; i++) func.push_back(0.0);
515 functional_cover =
true;
526 template <
class InputRange>
528 for (
int i = 0; i < this->num_points; i++) func.push_back(function[i]);
529 functional_cover =
true;
545 if (!functional_cover) {
546 std::cerr <<
"Cover needs to come from the preimages of a function." << std::endl;
549 if (type !=
"Nerve" && type !=
"GIC") {
550 std::cerr <<
"Type of complex needs to be specified." << std::endl;
555 Index_map index = boost::get(boost::vertex_index, one_skeleton);
558 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
559 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
560 reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
561 func[index[boost::target(*ei, one_skeleton)]]));
562 if (verbose) std::clog <<
"resolution = " << reso << std::endl;
563 resolution_double = reso;
566 if (type ==
"Nerve") {
567 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
568 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
569 reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
570 func[index[boost::target(*ei, one_skeleton)]]) /
572 if (verbose) std::clog <<
"resolution = " << reso << std::endl;
573 resolution_double = reso;
604 if (resolution_double == -1 && resolution_int == -1) {
605 std::cerr <<
"Number and/or length of intervals not specified" << std::endl;
609 std::cerr <<
"Gain not specified" << std::endl;
614 double minf = (std::numeric_limits<float>::max)();
615 double maxf = std::numeric_limits<float>::lowest();
616 for (
int i = 0; i < this->num_points; i++) {
617 minf = (std::min)(minf, func[i]);
618 maxf = (std::max)(maxf, func[i]);
620 if (verbose) std::clog <<
"Min function value = " << minf <<
" and Max function value = " << maxf << std::endl;
623 std::vector<std::pair<double, double> > intervals;
626 if (resolution_double == -1) {
627 double incr = (maxf - minf) / resolution_int;
629 double alpha = (incr * gain) / (2 - 2 * gain);
630 double y = minf + incr + alpha;
631 std::pair<double, double> interm(x, y);
632 intervals.push_back(interm);
633 for (
int i = 1; i < resolution_int - 1; i++) {
634 x = minf + i * incr - alpha;
635 y = minf + (i + 1) * incr + alpha;
636 std::pair<double, double> inter(x, y);
637 intervals.push_back(inter);
639 x = minf + (resolution_int - 1) * incr - alpha;
641 std::pair<double, double> interM(x, y);
642 intervals.push_back(interM);
643 res = intervals.size();
645 for (
int i = 0; i < res; i++)
646 std::clog <<
"Interval " << i <<
" = [" << intervals[i].first <<
", " << intervals[i].second <<
"]"
650 if (resolution_int == -1) {
652 double y = x + resolution_double;
653 while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
654 std::pair<double, double> inter(x, y);
655 intervals.push_back(inter);
656 x = y - gain * resolution_double;
657 y = x + resolution_double;
659 std::pair<double, double> interM(x, maxf);
660 intervals.push_back(interM);
661 res = intervals.size();
663 for (
int i = 0; i < res; i++)
664 std::clog <<
"Interval " << i <<
" = [" << intervals[i].first <<
", " << intervals[i].second <<
"]"
669 double y = x + resolution_double;
671 while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
672 std::pair<double, double> inter(x, y);
673 intervals.push_back(inter);
675 x = y - gain * resolution_double;
676 y = x + resolution_double;
678 res = intervals.size();
680 for (
int i = 0; i < res; i++)
681 std::clog <<
"Interval " << i <<
" = [" << intervals[i].first <<
", " << intervals[i].second <<
"]"
688 std::vector<int> points(this->num_points);
689 for (
int i = 0; i < this->num_points; i++) points[i] = i;
690 std::sort(points.begin(), points.end(), [
this](
int p1,
int p2){return (this->func[p1] < this->func[p2]);});
694 Index_map index = boost::get(boost::vertex_index, one_skeleton);
695 std::map<int, std::vector<int> > preimages;
696 std::map<int, double> funcstd;
698 if (verbose) std::clog <<
"Computing preimages..." << std::endl;
699 for (
int i = 0; i < res; i++) {
701 std::pair<double, double> inter1 = intervals[i];
707 std::pair<double, double> inter3 = intervals[i - 1];
708 while (func[points[tmp]] < inter3.second && tmp != this->num_points) {
709 preimages[i].push_back(points[tmp]);
717 std::pair<double, double> inter2 = intervals[i + 1];
718 while (func[points[tmp]] < inter2.first && tmp != this->num_points) {
719 preimages[i].push_back(points[tmp]);
724 while (func[points[tmp]] < inter1.second && tmp != this->num_points) {
725 preimages[i].push_back(points[tmp]);
730 std::pair<double, double> inter3 = intervals[i - 1];
731 while (func[points[tmp]] < inter3.second && tmp != this->num_points) {
732 preimages[i].push_back(points[tmp]);
735 while (tmp != this->num_points) {
736 preimages[i].push_back(points[tmp]);
743 funcstd[i] = 0.5 * (u + v);
747 if (verbose) std::clog <<
"Computing connected components (parallelized)..." << std::endl;
748 std::mutex covermutex, idmutex;
749 tbb::parallel_for(0, res, [&](
int i){
751 Graph G = one_skeleton.create_subgraph();
752 int num = preimages[i].size();
753 std::vector<int> component(num);
754 for (
int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
755 boost::connected_components(G, &component[0]);
759 for (
int j = 0; j < num; j++) {
761 if (component[j] > max) max = component[j];
764 int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
768 cover[preimages[i][j]].push_back(identifier);
769 cover_back[identifier].push_back(preimages[i][j]);
770 cover_fct[identifier] = i;
771 cover_std[identifier] = funcstd[i];
772 cover_color[identifier].second += func_color[preimages[i][j]];
773 cover_color[identifier].first += 1;
783 if (verbose) std::clog <<
"Computing connected components..." << std::endl;
784 for (
int i = 0; i < res; i++) {
786 Graph G = one_skeleton.create_subgraph();
787 int num = preimages[i].size();
788 std::vector<int> component(num);
789 for (
int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
790 boost::connected_components(G, &component[0]);
794 for (
int j = 0; j < num; j++) {
796 if (component[j] > max) max = component[j];
799 int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
802 cover[preimages[i][j]].push_back(identifier);
803 cover_back[identifier].push_back(preimages[i][j]);
804 cover_fct[identifier] = i;
805 cover_std[identifier] = funcstd[i];
806 cover_color[identifier].second += func_color[preimages[i][j]];
807 cover_color[identifier].first += 1;
815 maximal_dim =
id - 1;
816 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
817 iit->second.second /= iit->second.first;
830 std::vector<int> cov_elts, cov_number;
831 std::ifstream input(cover_file_name);
833 while (std::getline(input, line)) {
835 std::stringstream stream(line);
836 while (stream >> cov) {
837 cov_elts.push_back(cov);
838 cov_number.push_back(cov);
839 cover_fct[cov] = cov;
840 cover_color[cov].second += func_color[i];
841 cover_color[cov].first++;
842 cover_back[cov].push_back(i);
848 std::sort(cov_number.begin(), cov_number.end());
849 std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
850 cov_number.resize(std::distance(cov_number.begin(), it));
852 maximal_dim = cov_number.size() - 1;
853 for (
int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
854 cover_name = cover_file_name;
864 template <
class AssignmentRange>
866 std::vector<int> cov_elts, cov_number;
867 for(
int i=0; i < static_cast<int>(assignments.size()); i++){
869 for (
int cov : assignments[i]){
870 cov_elts.push_back(cov);
871 cov_number.push_back(cov);
872 cover_fct[cov] = cov;
873 auto& cc = cover_color[cov];
874 cc.second += func_color[i];
876 cover_back[cov].push_back(i);
881 std::sort(cov_number.begin(), cov_number.end());
882 std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
883 cov_number.resize(std::distance(cov_number.begin(), it));
885 maximal_dim = cov_number.size() - 1;
886 for (
int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
896 template <
typename Distance>
898 voronoi_subsamples.resize(m);
899 SampleWithoutReplacement(this->num_points, m, voronoi_subsamples);
900 if (distances.size() == 0) compute_pairwise_distances(distance);
902 Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
903 Index_map index = boost::get(boost::vertex_index, one_skeleton);
904 std::vector<double> mindist(this->num_points);
905 for (
int j = 0; j < this->num_points; j++) mindist[j] = (std::numeric_limits<double>::max)();
909 if (verbose) std::clog <<
"Computing geodesic distances (parallelized)..." << std::endl;
910 std::mutex coverMutex; std::mutex mindistMutex;
911 tbb::parallel_for(0, m, [&](
int i){
912 int seed = voronoi_subsamples[i];
913 std::vector<double> dmap(this->num_points);
914 boost::dijkstra_shortest_paths(
915 one_skeleton, vertices[seed],
916 boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
918 coverMutex.lock(); mindistMutex.lock();
919 for (
int j = 0; j < this->num_points; j++)
920 if (mindist[j] > dmap[j]) {
921 mindist[j] = dmap[j];
922 if (cover[j].size() == 0)
923 cover[j].push_back(i);
927 coverMutex.unlock(); mindistMutex.unlock();
930 for (
int i = 0; i < m; i++) {
931 if (verbose) std::clog <<
"Computing geodesic distances to seed " << i <<
"..." << std::endl;
932 int seed = voronoi_subsamples[i];
933 std::vector<double> dmap(this->num_points);
934 boost::dijkstra_shortest_paths(
935 one_skeleton, vertices[seed],
936 boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
938 for (
int j = 0; j < this->num_points; j++)
939 if (mindist[j] > dmap[j]) {
940 mindist[j] = dmap[j];
941 if (cover[j].size() == 0)
942 cover[j].push_back(i);
949 for (
int i = 0; i < this->num_points; i++) {
950 cover_back[cover[i][0]].push_back(i);
951 cover_color[cover[i][0]].second += func_color[i];
952 cover_color[cover[i][0]].first++;
954 for (
int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first;
956 cover_name =
"Voronoi";
976 double subcolor(
int c) {
return cover_color[c].second; }
991 std::ifstream input(color_file_name);
994 while (std::getline(input, line)) {
995 std::stringstream stream(line);
997 func_color.push_back(f);
1000 color_name = color_file_name;
1010 if(point_cloud[0].size() > 0){
1011 for (
int i = 0; i < this->num_points; i++) func_color.push_back(point_cloud[i][k]);
1012 color_name =
"coordinate ";
1013 color_name.append(std::to_string(k));
1016 std::cerr <<
"Only pairwise distances provided---cannot access " << k <<
"th coordinate; returning null vector instead" << std::endl;
1017 for (
int i = 0; i < this->num_points; i++) func.push_back(0.0);
1018 functional_cover =
true;
1019 cover_name =
"null";
1030 for (
unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
1039 std::string mapp = point_cloud_name +
"_sc.dot";
1040 std::ofstream graphic(mapp);
1042 double maxv = std::numeric_limits<double>::lowest();
1043 double minv = (std::numeric_limits<double>::max)();
1044 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1045 maxv = (std::max)(maxv, iit->second.second);
1046 minv = (std::min)(minv, iit->second.second);
1050 std::vector<int> nodes;
1053 graphic <<
"graph GIC {" << std::endl;
1055 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1056 if (iit->second.first > mask) {
1057 nodes.push_back(iit->first);
1058 name2id[iit->first] = id;
1059 name2idinv[id] = iit->first;
1061 graphic << name2id[iit->first] <<
"[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
1062 <<
":" << iit->second.first <<
"\" style=filled fillcolor=\""
1063 << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 <<
", 1, 1\"]" << std::endl;
1068 int num_simplices = simplices.size();
1069 for (
int i = 0; i < num_simplices; i++)
1070 if (simplices[i].size() == 2) {
1071 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
1072 graphic <<
" " << name2id[simplices[i][0]] <<
" -- " << name2id[simplices[i][1]] <<
" [weight=15];"
1079 std::clog << mapp <<
" file generated. It can be visualized with e.g. neato." << std::endl;
1087 int num_simplices = simplices.size();
1089 std::string mapp = point_cloud_name +
"_sc.txt";
1090 std::ofstream graphic(mapp);
1092 for (
int i = 0; i < num_simplices; i++)
1093 if (simplices[i].size() == 2)
1094 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++;
1096 graphic << point_cloud_name << std::endl;
1097 graphic << cover_name << std::endl;
1098 graphic << color_name << std::endl;
1099 graphic << resolution_double <<
" " << gain << std::endl;
1100 graphic << cover_color.size() <<
" " << num_edges << std::endl;
1103 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1104 graphic <<
id <<
" " << iit->second.second <<
" " << iit->second.first << std::endl;
1105 name2id[iit->first] = id;
1106 name2idinv[id] = iit->first;
1110 for (
int i = 0; i < num_simplices; i++)
1111 if (simplices[i].size() == 2)
1112 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
1113 graphic << name2id[simplices[i][0]] <<
" " << name2id[simplices[i][1]] << std::endl;
1116 <<
" generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
1126 assert(cover_name ==
"Voronoi");
1128 int m = voronoi_subsamples.size();
1131 std::vector<std::vector<int> > edges, faces;
1132 int numsimplices = simplices.size();
1134 std::string mapp = point_cloud_name +
"_sc.off";
1135 std::ofstream graphic(mapp);
1137 graphic <<
"OFF" << std::endl;
1138 for (
int i = 0; i < numsimplices; i++) {
1139 if (simplices[i].size() == 2) {
1141 edges.push_back(simplices[i]);
1143 if (simplices[i].size() == 3) {
1145 faces.push_back(simplices[i]);
1148 graphic << m <<
" " << numedges + numfaces << std::endl;
1149 for (
int i = 0; i < m; i++) {
1150 if (data_dimension <= 3) {
1151 for (
int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] <<
" ";
1152 for (
int j = data_dimension; j < 3; j++) graphic << 0 <<
" ";
1153 graphic << std::endl;
1155 for (
int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] <<
" ";
1158 for (
int i = 0; i < numedges; i++) graphic << 2 <<
" " << edges[i][0] <<
" " << edges[i][1] << std::endl;
1159 for (
int i = 0; i < numfaces; i++)
1160 graphic << 3 <<
" " << faces[i][0] <<
" " << faces[i][1] <<
" " << faces[i][2] << std::endl;
1162 std::clog << mapp <<
" generated. It can be visualized with e.g. geomview." << std::endl;
1177 double maxf = std::numeric_limits<double>::lowest();
1178 double minf = (std::numeric_limits<double>::max)();
1179 for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1180 maxf = (std::max)(maxf, it->second);
1181 minf = (std::min)(minf, it->second);
1185 for (
auto const& simplex : simplices) {
1186 std::vector<int> splx = simplex;
1191 for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1192 int vertex = it->first;
float val = it->second;
1193 int vert[] = {vertex};
int edge[] = {vertex, -2};
1207 for (
int i = 0; i < max_dim; i++) {
1209 int num_bars = bars.size();
if(i == 0) num_bars -= 1;
1210 if(verbose) std::clog << num_bars <<
" interval(s) in dimension " << i <<
":" << std::endl;
1211 for (
int j = 0; j < num_bars; j++) {
1212 double birth = bars[j].first;
1213 double death = bars[j].second;
1214 if (i == 0 && std::isinf(death))
continue;
1216 birth = minf + (birth + 2) * (maxf - minf);
1218 birth = minf + (2 - birth) * (maxf - minf);
1220 death = minf + (death + 2) * (maxf - minf);
1222 death = minf + (2 - death) * (maxf - minf);
1223 PD.push_back(std::pair<double, double>(birth, death));
1224 if (verbose) std::clog <<
" [" << birth <<
", " << death <<
"]" << std::endl;
1237 unsigned int sz = distribution.size();
1239 for (
unsigned int i = 0; i < N - sz; i++) {
1240 if (verbose) std::clog <<
"Computing " << i <<
"th bootstrap, bottleneck distance = ";
1242 Cover_complex Cboot; Cboot.num_points = this->num_points; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover =
true;
1244 std::vector<int> boot(this->num_points);
1245 for (
int j = 0; j < this->num_points; j++) {
1246 double u = GetUniform();
1247 int id = std::floor(u * (this->num_points)); boot[j] = id;
1248 Cboot.point_cloud.push_back(this->point_cloud[
id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[
id]);
1249 boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
1253 for (
int j = 0; j < this->num_points; j++) {
1254 std::vector<double> dist(this->num_points);
1255 for (
int k = 0; k < this->num_points; k++) dist[k] = distances[boot[j]][boot[k]];
1256 Cboot.distances.push_back(dist);
1265#ifdef GUDHI_GIC_USE_CGAL
1267#elif defined GUDHI_GIC_USE_HERA
1268 double db = hera::bottleneckDistExact(this->PD, Cboot.PD);
1271 throw std::logic_error(
"This function requires CGAL or Hera for the bottleneck distance.");
1273 if (verbose) std::clog << db << std::endl;
1274 distribution.push_back(db);
1277 std::sort(distribution.begin(), distribution.end());
1288 unsigned int N = distribution.size();
1289 double d = distribution[std::floor(alpha * N)];
1290 if (verbose) std::clog <<
"Distance corresponding to confidence " << alpha <<
" is " << d << std::endl;
1301 unsigned int N = distribution.size();
1303 for (
unsigned int i = 0; i < N; i++)
1304 if (distribution[i] >= d){ level = i * 1.0 / N;
break; }
1305 if (verbose) std::clog <<
"Confidence level of distance " << d <<
" is " << level << std::endl;
1315 double distancemin = (std::numeric_limits<double>::max)();
int N = PD.size();
1316 for (
int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first));
1318 if (verbose) std::clog <<
"p value = " << p_value << std::endl;
1332 template <
typename SimplicialComplex>
1334 unsigned int dimension = 0;
1335 for (
auto const& simplex : simplices) {
1336 int numvert = simplex.size();
1337 double filt = std::numeric_limits<double>::lowest();
1338 for (
int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt);
1339 complex.insert_simplex_and_subfaces(simplex, filt);
1340 if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
1348 if (type !=
"Nerve" && type !=
"GIC") {
1349 std::cerr <<
"Type of complex needs to be specified." << std::endl;
1353 if (type ==
"Nerve") {
1355 std::sort(simplices.begin(), simplices.end());
1356 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1357 simplices.resize(std::distance(simplices.begin(), it));
1360 if (type ==
"GIC") {
1361 Index_map index = boost::get(boost::vertex_index, one_skeleton);
1363 if (functional_cover) {
1368 throw std::invalid_argument(
1369 "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
1372 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1373 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
1374 int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
1375 for (
int i = 0; i < nums; i++) {
1376 int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
1377 int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
1378 for (
int j = 0; j < numt; j++) {
1379 int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
1380 if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
1381 std::vector<int> edge(2);
1382 edge[0] = (std::min)(vs, vt);
1383 edge[1] = (std::max)(vs, vt);
1384 simplices.push_back(edge);
1391 std::sort(simplices.begin(), simplices.end());
1392 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1393 simplices.resize(std::distance(simplices.begin(), it));
1398 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1399 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
1400 if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
1401 cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
1402 std::vector<int> edge(2);
1403 edge[0] = index[boost::source(*ei, one_skeleton)];
1404 edge[1] = index[boost::target(*ei, one_skeleton)];
1417 std::vector<int> simplx;
1419 unsigned int sz = cover[vertex].size();
1420 for (
unsigned int i = 0; i < sz; i++) {
1421 simplx.push_back(cover[vertex][i]);
1424 std::sort(simplx.begin(), simplx.end());
1425 std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
1426 simplx.resize(std::distance(simplx.begin(), it));
1427 simplices.push_back(simplx);
1430 std::sort(simplices.begin(), simplices.end());
1431 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1432 simplices.resize(std::distance(simplices.begin(), it));
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:90
void assign_filtration(Simplex_handle sh, Filtration_value fv)
Sets the filtration value of a simplex.
Definition: Simplex_tree.h:550
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:1410
Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) const
Returns a range over the vertices of a simplex.
Definition: Simplex_tree.h:286
bool has_children(SimplexHandle sh) const
Returns true if the node in the simplex tree pointed by sh has children.
Definition: Simplex_tree.h:637
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:650
void expansion(int max_dim)
Expands the Simplex_tree containing only its one skeleton until dimension max_dim.
Definition: Simplex_tree.h:1194
int dimension(Simplex_handle sh)
Returns the dimension of a simplex.
Definition: Simplex_tree.h:609
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:818
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:560
Complex_simplex_range complex_simplex_range()
Returns a range over the simplices of the simplicial complex.
Definition: Simplex_tree.h:241
Cover complex data structure.
Definition: GIC.h:99
void set_function_from_file(const std::string &func_file_name)
Creates the function f from a file containing the function values.
Definition: GIC.h:485
double set_automatic_resolution()
Computes the optimal length of intervals (i.e. the smallest interval length avoiding discretization a...
Definition: GIC.h:544
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:897
void set_resolution_with_interval_number(int reso)
Sets a number of intervals from a value stored in memory.
Definition: GIC.h:591
void set_mask(int nodemask)
Sets the mask, which is a threshold integer such that nodes in the complex that contain a number of d...
Definition: GIC.h:219
Persistence_diagram compute_PD()
Computes the extended persistence diagram of the complex.
Definition: GIC.h:1173
double compute_distance_from_confidence_level(double alpha)
Computes the bottleneck distance threshold corresponding to a specific confidence level.
Definition: GIC.h:1287
void set_graph_from_rips(double threshold, Distance distance)
Creates a graph G from a Rips complex.
Definition: GIC.h:360
void set_cover_from_file(const std::string &cover_file_name)
Creates the cover C from a file containing the cover elements of each point (the order has to be the ...
Definition: GIC.h:827
void set_graph_from_file(const std::string &graph_file_name)
Creates a graph G from a file containing the edges.
Definition: GIC.h:327
void set_cover_from_range(AssignmentRange const &assignments)
Creates the cover C from a vector of assignments stored in memory. The assignments,...
Definition: GIC.h:865
void create_complex(SimplicialComplex &complex)
Creates the simplicial complex.
Definition: GIC.h:1333
void set_type(const std::string &t)
Specifies whether the type of the output simplicial complex.
Definition: GIC.h:188
void set_distances_from_range(const std::vector< std::vector< double > > &distance_matrix)
Reads and stores the distance matrices from vector stored in memory.
Definition: GIC.h:390
void find_simplices()
Computes the simplices of the simplicial complex.
Definition: GIC.h:1347
void set_function_from_range(InputRange const &function)
Creates the function f from a vector stored in memory.
Definition: GIC.h:527
void set_cover_from_function()
Creates a cover C from the preimages of the function f.
Definition: GIC.h:603
double subcolor(int c)
Returns the mean color corresponding to a specific node of the created complex.
Definition: GIC.h:976
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:989
void write_info()
Creates a .txt file called SC.txt describing the 1-skeleton, which can then be plotted with e....
Definition: GIC.h:1086
void plot_OFF()
Creates a .off file called SC.off for 3D visualization, which contains the 2-skeleton of the GIC....
Definition: GIC.h:1125
void plot_DOT()
Creates a .dot file called SC.dot for neato (part of the graphviz package) once the simplicial comple...
Definition: GIC.h:1038
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:1029
void set_subsampling(double constant, double power)
Sets the constants used to subsample the data set. These constants are explained in .
Definition: GIC.h:206
void set_point_cloud_from_range(const std::vector< std::vector< double > > &point_cloud)
Reads and stores the input point cloud from vector stored in memory.
Definition: GIC.h:229
double set_graph_from_automatic_rips(Distance distance, int N=100)
Creates a graph G from a Rips complex whose threshold value is automatically tuned with subsampling—s...
Definition: GIC.h:431
void set_resolution_with_interval_length(double reso)
Sets a length of intervals from a value stored in memory.
Definition: GIC.h:585
const std::vector< int > & subpopulation(int c)
Returns the data subset corresponding to a specific node of the created complex.
Definition: GIC.h:966
void set_gain(double g=0.3)
Sets a gain from a value stored in memory (default value 0.3).
Definition: GIC.h:597
void set_function_from_coordinate(int k)
Creates the function f from the k-th coordinate of the point cloud P.
Definition: GIC.h:506
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:1009
void compute_distribution(unsigned int N=100)
Computes bootstrapped distances distribution.
Definition: GIC.h:1236
double compute_confidence_level_from_distance(double d)
Computes the confidence level of a specific bottleneck distance threshold.
Definition: GIC.h:1300
void set_graph_from_OFF()
Creates a graph G from the triangulation given by the input .OFF file.
Definition: GIC.h:344
double compute_p_value()
Computes the p-value, i.e. the opposite of the confidence level of the largest bottleneck distance pr...
Definition: GIC.h:1314
bool read_point_cloud(const std::string &off_file_name)
Reads and stores the input point cloud from .(n)OFF file.
Definition: GIC.h:244
void set_verbose(bool verb=false)
Specifies whether the program should display information or not.
Definition: GIC.h:196
Computes the persistent cohomology of a filtered complex.
Definition: Persistent_cohomology.h:53
std::vector< std::pair< Filtration_value, Filtration_value > > intervals_in_dimension(int dimension)
Returns persistence intervals for a given dimension.
Definition: Persistent_cohomology.h:682
void compute_persistent_cohomology(Filtration_value min_interval_length=0)
Compute the persistent homology of the filtered simplicial complex.
Definition: Persistent_cohomology.h:173
void init_coefficients(int charac)
Initializes the coefficient field.
Definition: Persistent_cohomology.h:157
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