17 #include <tbb/parallel_for.h>
22 #include <gudhi/Debug_utils.h>
23 #include <gudhi/graph_simplicial_complex.h>
25 #include <gudhi/Simplex_tree.h>
26 #include <gudhi/Rips_complex.h>
27 #include <gudhi/Points_off_io.h>
29 #include <gudhi/Persistent_cohomology.h>
30 #include <gudhi/Bottleneck.h>
32 #include <boost/config.hpp>
33 #include <boost/graph/graph_traits.hpp>
34 #include <boost/graph/adjacency_list.hpp>
35 #include <boost/graph/connected_components.hpp>
36 #include <boost/graph/dijkstra_shortest_paths.hpp>
37 #include <boost/graph/subgraph.hpp>
38 #include <boost/graph/graph_utility.hpp>
40 #include <CGAL/version.h>
55 namespace cover_complex {
60 using Persistence_diagram = std::vector<std::pair<double, double> >;
61 using Graph = boost::subgraph<
62 boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
63 boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
64 using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
65 using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
66 using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
88 template <
typename Po
int>
94 std::vector<Point> point_cloud;
95 std::vector<std::vector<double> > distances;
100 std::vector<double> func;
101 std::vector<double> func_color;
102 bool functional_cover =
false;
104 Graph one_skeleton_OFF;
106 std::vector<Vertex_t> vertices;
108 std::vector<std::vector<int> > simplices;
109 std::vector<int> voronoi_subsamples;
111 Persistence_diagram PD;
112 std::vector<double> distribution;
114 std::vector<std::vector<int> >
116 std::map<int, std::vector<int> >
118 std::map<int, double> cover_std;
122 std::map<int, std::pair<int, double> >
125 int resolution_int = -1;
126 double resolution_double = -1;
128 double rate_constant = 10;
129 double rate_power = 0.001;
132 std::map<int, int> name2id, name2idinv;
134 std::string cover_name;
135 std::string point_cloud_name;
136 std::string color_name;
139 void remove_edges(Graph& G) {
140 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
141 for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
145 double GetUniform() {
146 thread_local std::default_random_engine re;
147 std::uniform_real_distribution<double> Dist(0, 1);
152 void SampleWithoutReplacement(
int populationSize,
int sampleSize, std::vector<int>& samples) {
156 while (m < sampleSize) {
158 if ((populationSize - t) * u >= sampleSize - m) {
197 rate_constant = constant;
220 n = point_cloud.size(); data_dimension = point_cloud[0].size();
221 point_cloud_name =
"matrix"; cover.resize(n);
222 for(
int i = 0; i < n; i++){
223 boost::add_vertex(one_skeleton_OFF);
224 vertices.push_back(boost::add_vertex(one_skeleton));
226 this->point_cloud = point_cloud;
235 point_cloud_name = off_file_name;
236 std::ifstream input(off_file_name);
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(
' ')];
245 if (strcmp((
char*)line.c_str(),
"nOFF") == 0) {
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(
' ')];
252 std::stringstream stream(line);
253 stream >> data_dimension;
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(
' ')];
265 std::stringstream stream(line);
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();
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>());
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);
302 return input.is_open();
318 remove_edges(one_skeleton);
320 std::ifstream input(graph_file_name);
323 while (std::getline(input, line)) {
324 std::stringstream stream(line);
326 while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
335 remove_edges(one_skeleton);
336 if (num_edges(one_skeleton_OFF))
337 one_skeleton = one_skeleton_OFF;
339 std::cerr <<
"No triangulation read in OFF file!" << std::endl;
349 template <
typename 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,
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)]]);
381 n = distance_matrix.size(); data_dimension = 0; point_cloud_name =
"matrix";
382 cover.resize(n); point_cloud.resize(n);
383 for(
int i = 0; i < n; i++){
384 boost::add_vertex(one_skeleton_OFF);
385 vertices.push_back(boost::add_vertex(one_skeleton));
387 distances = distance_matrix;
393 template <
typename Distance>
394 void compute_pairwise_distances(Distance ref_distance) {
396 std::vector<double> zeros(n);
397 for (
int i = 0; i < n; i++) distances.push_back(zeros);
398 std::string distance = point_cloud_name +
"_dist";
399 std::ifstream input(distance, std::ios::out | std::ios::binary);
402 if (verbose) std::clog <<
"Reading distances..." << std::endl;
403 for (
int i = 0; i < n; i++) {
404 for (
int j = i; j < n; j++) {
405 input.read((
char*)&d, 8);
412 if (verbose) std::clog <<
"Computing distances..." << std::endl;
414 std::ofstream output(distance, std::ios::out | std::ios::binary);
415 for (
int i = 0; i < n; i++) {
416 int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10;
417 if (state == 0 && verbose) std::clog <<
"\r" << state <<
"%" << std::flush;
418 for (
int j = i; j < n; j++) {
419 double dis = ref_distance(point_cloud[i], point_cloud[j]);
420 distances[i][j] = dis;
421 distances[j][i] = dis;
422 output.write((
char*)&dis, 8);
426 if (verbose) std::clog << std::endl;
440 template <
typename Distance>
442 int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant))));
443 m = (std::min)(m, n - 1);
446 if (verbose) std::clog << n <<
" points in R^" << data_dimension << std::endl;
447 if (verbose) std::clog <<
"Subsampling " << m <<
" points" << std::endl;
449 if (distances.size() == 0) compute_pairwise_distances(distance);
452 std::mutex deltamutex;
453 tbb::parallel_for(0, N, [&](
int 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);
463 delta += hausdorff_dist / N;
467 for (
int i = 0; i < N; i++) {
468 std::vector<int> samples(m);
469 SampleWithoutReplacement(n, m, samples);
470 double hausdorff_dist = 0;
471 for (
int j = 0; j < n; j++) {
472 double mj = distances[j][samples[0]];
473 for (
int k = 1; k < m; k++) mj = (std::min)(mj, distances[j][samples[k]]);
474 hausdorff_dist = (std::max)(hausdorff_dist, mj);
476 delta += hausdorff_dist / N;
480 if (verbose) std::clog <<
"delta = " << delta << std::endl;
497 std::ifstream input(func_file_name);
500 while (std::getline(input, line)) {
501 std::stringstream stream(line);
506 functional_cover =
true;
507 cover_name = func_file_name;
517 if(point_cloud[0].size() > 0){
518 for (
int i = 0; i < n; i++) func.push_back(point_cloud[i][k]);
519 functional_cover =
true;
520 cover_name =
"coordinate " + std::to_string(k);
523 std::cerr <<
"Only pairwise distances provided---cannot access " << k <<
"th coordinate; returning null vector instead" << std::endl;
524 for (
int i = 0; i < n; i++) func.push_back(0.0);
525 functional_cover =
true;
536 template <
class InputRange>
538 for (
int i = 0; i < n; i++) func.push_back(
function[i]);
539 functional_cover =
true;
555 if (!functional_cover) {
556 std::cerr <<
"Cover needs to come from the preimages of a function." << std::endl;
559 if (type !=
"Nerve" && type !=
"GIC") {
560 std::cerr <<
"Type of complex needs to be specified." << std::endl;
565 Index_map index = boost::get(boost::vertex_index, one_skeleton);
568 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
569 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
570 reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
571 func[index[boost::target(*ei, one_skeleton)]]));
572 if (verbose) std::clog <<
"resolution = " << reso << std::endl;
573 resolution_double = reso;
576 if (type ==
"Nerve") {
577 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
578 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
579 reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
580 func[index[boost::target(*ei, one_skeleton)]]) /
582 if (verbose) std::clog <<
"resolution = " << reso << std::endl;
583 resolution_double = reso;
614 if (resolution_double == -1 && resolution_int == -1) {
615 std::cerr <<
"Number and/or length of intervals not specified" << std::endl;
619 std::cerr <<
"Gain not specified" << std::endl;
624 double minf = (std::numeric_limits<float>::max)();
625 double maxf = std::numeric_limits<float>::lowest();
626 for (
int i = 0; i < n; i++) {
627 minf = (std::min)(minf, func[i]);
628 maxf = (std::max)(maxf, func[i]);
630 if (verbose) std::clog <<
"Min function value = " << minf <<
" and Max function value = " << maxf << std::endl;
633 std::vector<std::pair<double, double> > intervals;
636 if (resolution_double == -1) {
637 double incr = (maxf - minf) / resolution_int;
639 double alpha = (incr * gain) / (2 - 2 * gain);
640 double y = minf + incr + alpha;
641 std::pair<double, double> interm(x, y);
642 intervals.push_back(interm);
643 for (
int i = 1; i < resolution_int - 1; i++) {
644 x = minf + i * incr - alpha;
645 y = minf + (i + 1) * incr + alpha;
646 std::pair<double, double> inter(x, y);
647 intervals.push_back(inter);
649 x = minf + (resolution_int - 1) * incr - alpha;
651 std::pair<double, double> interM(x, y);
652 intervals.push_back(interM);
653 res = intervals.size();
655 for (
int i = 0; i < res; i++)
656 std::clog <<
"Interval " << i <<
" = [" << intervals[i].first <<
", " << intervals[i].second <<
"]"
660 if (resolution_int == -1) {
662 double y = x + resolution_double;
663 while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
664 std::pair<double, double> inter(x, y);
665 intervals.push_back(inter);
666 x = y - gain * resolution_double;
667 y = x + resolution_double;
669 std::pair<double, double> interM(x, maxf);
670 intervals.push_back(interM);
671 res = intervals.size();
673 for (
int i = 0; i < res; i++)
674 std::clog <<
"Interval " << i <<
" = [" << intervals[i].first <<
", " << intervals[i].second <<
"]"
679 double y = x + resolution_double;
681 while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) {
682 std::pair<double, double> inter(x, y);
683 intervals.push_back(inter);
685 x = y - gain * resolution_double;
686 y = x + resolution_double;
688 res = intervals.size();
690 for (
int i = 0; i < res; i++)
691 std::clog <<
"Interval " << i <<
" = [" << intervals[i].first <<
", " << intervals[i].second <<
"]"
698 std::vector<int> points(n);
699 for (
int i = 0; i < n; i++) points[i] = i;
700 std::sort(points.begin(), points.end(), [
this](
int p1,
int p2){return (this->func[p1] < this->func[p2]);});
704 Index_map index = boost::get(boost::vertex_index, one_skeleton);
705 std::map<int, std::vector<int> > preimages;
706 std::map<int, double> funcstd;
708 if (verbose) std::clog <<
"Computing preimages..." << std::endl;
709 for (
int i = 0; i < res; i++) {
711 std::pair<double, double> inter1 = intervals[i];
717 std::pair<double, double> inter3 = intervals[i - 1];
718 while (func[points[tmp]] < inter3.second && tmp != n) {
719 preimages[i].push_back(points[tmp]);
727 std::pair<double, double> inter2 = intervals[i + 1];
728 while (func[points[tmp]] < inter2.first && tmp != n) {
729 preimages[i].push_back(points[tmp]);
734 while (func[points[tmp]] < inter1.second && tmp != n) {
735 preimages[i].push_back(points[tmp]);
740 std::pair<double, double> inter3 = intervals[i - 1];
741 while (func[points[tmp]] < inter3.second && tmp != n) {
742 preimages[i].push_back(points[tmp]);
746 preimages[i].push_back(points[tmp]);
753 funcstd[i] = 0.5 * (u + v);
757 if (verbose) std::clog <<
"Computing connected components (parallelized)..." << std::endl;
758 std::mutex covermutex, idmutex;
759 tbb::parallel_for(0, res, [&](
int i){
761 Graph G = one_skeleton.create_subgraph();
762 int num = preimages[i].size();
763 std::vector<int> component(num);
764 for (
int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
765 boost::connected_components(G, &component[0]);
769 for (
int j = 0; j < num; j++) {
771 if (component[j] > max) max = component[j];
774 int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
778 cover[preimages[i][j]].push_back(identifier);
779 cover_back[identifier].push_back(preimages[i][j]);
780 cover_fct[identifier] = i;
781 cover_std[identifier] = funcstd[i];
782 cover_color[identifier].second += func_color[preimages[i][j]];
783 cover_color[identifier].first += 1;
793 if (verbose) std::clog <<
"Computing connected components..." << std::endl;
794 for (
int i = 0; i < res; i++) {
796 Graph G = one_skeleton.create_subgraph();
797 int num = preimages[i].size();
798 std::vector<int> component(num);
799 for (
int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
800 boost::connected_components(G, &component[0]);
804 for (
int j = 0; j < num; j++) {
806 if (component[j] > max) max = component[j];
809 int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
812 cover[preimages[i][j]].push_back(identifier);
813 cover_back[identifier].push_back(preimages[i][j]);
814 cover_fct[identifier] = i;
815 cover_std[identifier] = funcstd[i];
816 cover_color[identifier].second += func_color[preimages[i][j]];
817 cover_color[identifier].first += 1;
825 maximal_dim =
id - 1;
826 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
827 iit->second.second /= iit->second.first;
840 std::vector<int> cov_elts, cov_number;
841 std::ifstream input(cover_file_name);
843 while (std::getline(input, line)) {
845 std::stringstream stream(line);
846 while (stream >> cov) {
847 cov_elts.push_back(cov);
848 cov_number.push_back(cov);
849 cover_fct[cov] = cov;
850 cover_color[cov].second += func_color[i];
851 cover_color[cov].first++;
852 cover_back[cov].push_back(i);
858 std::sort(cov_number.begin(), cov_number.end());
859 std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
860 cov_number.resize(std::distance(cov_number.begin(), it));
862 maximal_dim = cov_number.size() - 1;
863 for (
int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
864 cover_name = cover_file_name;
874 template <
typename Distance>
876 voronoi_subsamples.resize(m);
877 SampleWithoutReplacement(n, m, voronoi_subsamples);
878 if (distances.size() == 0) compute_pairwise_distances(distance);
880 Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
881 Index_map index = boost::get(boost::vertex_index, one_skeleton);
882 std::vector<double> mindist(n);
883 for (
int j = 0; j < n; j++) mindist[j] = (std::numeric_limits<double>::max)();
887 if (verbose) std::clog <<
"Computing geodesic distances (parallelized)..." << std::endl;
888 std::mutex coverMutex; std::mutex mindistMutex;
889 tbb::parallel_for(0, m, [&](
int i){
890 int seed = voronoi_subsamples[i];
891 std::vector<double> dmap(n);
892 boost::dijkstra_shortest_paths(
893 one_skeleton, vertices[seed],
894 boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
896 coverMutex.lock(); mindistMutex.lock();
897 for (
int j = 0; j < n; j++)
898 if (mindist[j] > dmap[j]) {
899 mindist[j] = dmap[j];
900 if (cover[j].size() == 0)
901 cover[j].push_back(i);
905 coverMutex.unlock(); mindistMutex.unlock();
908 for (
int i = 0; i < m; i++) {
909 if (verbose) std::clog <<
"Computing geodesic distances to seed " << i <<
"..." << std::endl;
910 int seed = voronoi_subsamples[i];
911 std::vector<double> dmap(n);
912 boost::dijkstra_shortest_paths(
913 one_skeleton, vertices[seed],
914 boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
916 for (
int j = 0; j < n; j++)
917 if (mindist[j] > dmap[j]) {
918 mindist[j] = dmap[j];
919 if (cover[j].size() == 0)
920 cover[j].push_back(i);
927 for (
int i = 0; i < n; i++) {
928 cover_back[cover[i][0]].push_back(i);
929 cover_color[cover[i][0]].second += func_color[i];
930 cover_color[cover[i][0]].first++;
932 for (
int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first;
934 cover_name =
"Voronoi";
944 const std::vector<int>&
subpopulation(
int c) {
return cover_back[name2idinv[c]]; }
959 std::ifstream input(color_file_name);
962 while (std::getline(input, line)) {
963 std::stringstream stream(line);
965 func_color.push_back(f);
968 color_name = color_file_name;
978 if(point_cloud[0].size() > 0){
979 for (
int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]);
980 color_name =
"coordinate ";
981 color_name.append(std::to_string(k));
984 std::cerr <<
"Only pairwise distances provided---cannot access " << k <<
"th coordinate; returning null vector instead" << std::endl;
985 for (
int i = 0; i < n; i++) func.push_back(0.0);
986 functional_cover =
true;
998 for (
unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
1007 std::string mapp = point_cloud_name +
"_sc.dot";
1008 std::ofstream graphic(mapp);
1010 double maxv = std::numeric_limits<double>::lowest();
1011 double minv = (std::numeric_limits<double>::max)();
1012 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1013 maxv = (std::max)(maxv, iit->second.second);
1014 minv = (std::min)(minv, iit->second.second);
1018 std::vector<int> nodes;
1021 graphic <<
"graph GIC {" << std::endl;
1023 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1024 if (iit->second.first > mask) {
1025 nodes.push_back(iit->first);
1026 name2id[iit->first] = id;
1027 name2idinv[id] = iit->first;
1029 graphic << name2id[iit->first] <<
"[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
1030 <<
":" << iit->second.first <<
"\" style=filled fillcolor=\""
1031 << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 <<
", 1, 1\"]" << std::endl;
1036 int num_simplices = simplices.size();
1037 for (
int i = 0; i < num_simplices; i++)
1038 if (simplices[i].size() == 2) {
1039 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
1040 graphic <<
" " << name2id[simplices[i][0]] <<
" -- " << name2id[simplices[i][1]] <<
" [weight=15];"
1047 std::clog << mapp <<
" file generated. It can be visualized with e.g. neato." << std::endl;
1055 int num_simplices = simplices.size();
1057 std::string mapp = point_cloud_name +
"_sc.txt";
1058 std::ofstream graphic(mapp);
1060 for (
int i = 0; i < num_simplices; i++)
1061 if (simplices[i].size() == 2)
1062 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++;
1064 graphic << point_cloud_name << std::endl;
1065 graphic << cover_name << std::endl;
1066 graphic << color_name << std::endl;
1067 graphic << resolution_double <<
" " << gain << std::endl;
1068 graphic << cover_color.size() <<
" " << num_edges << std::endl;
1071 for (std::map<
int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
1072 graphic <<
id <<
" " << iit->second.second <<
" " << iit->second.first << std::endl;
1073 name2id[iit->first] = id;
1074 name2idinv[id] = iit->first;
1078 for (
int i = 0; i < num_simplices; i++)
1079 if (simplices[i].size() == 2)
1080 if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
1081 graphic << name2id[simplices[i][0]] <<
" " << name2id[simplices[i][1]] << std::endl;
1084 <<
" generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
1094 assert(cover_name ==
"Voronoi");
1096 int m = voronoi_subsamples.size();
1099 std::vector<std::vector<int> > edges, faces;
1100 int numsimplices = simplices.size();
1102 std::string mapp = point_cloud_name +
"_sc.off";
1103 std::ofstream graphic(mapp);
1105 graphic <<
"OFF" << std::endl;
1106 for (
int i = 0; i < numsimplices; i++) {
1107 if (simplices[i].size() == 2) {
1109 edges.push_back(simplices[i]);
1111 if (simplices[i].size() == 3) {
1113 faces.push_back(simplices[i]);
1116 graphic << m <<
" " << numedges + numfaces << std::endl;
1117 for (
int i = 0; i < m; i++) {
1118 if (data_dimension <= 3) {
1119 for (
int j = 0; j < data_dimension; j++) graphic << point_cloud[voronoi_subsamples[i]][j] <<
" ";
1120 for (
int j = data_dimension; j < 3; j++) graphic << 0 <<
" ";
1121 graphic << std::endl;
1123 for (
int j = 0; j < 3; j++) graphic << point_cloud[voronoi_subsamples[i]][j] <<
" ";
1126 for (
int i = 0; i < numedges; i++) graphic << 2 <<
" " << edges[i][0] <<
" " << edges[i][1] << std::endl;
1127 for (
int i = 0; i < numfaces; i++)
1128 graphic << 3 <<
" " << faces[i][0] <<
" " << faces[i][1] <<
" " << faces[i][2] << std::endl;
1130 std::clog << mapp <<
" generated. It can be visualized with e.g. geomview." << std::endl;
1145 double maxf = std::numeric_limits<double>::lowest();
1146 double minf = (std::numeric_limits<double>::max)();
1147 for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1148 maxf = (std::max)(maxf, it->second);
1149 minf = (std::min)(minf, it->second);
1153 for (
auto const& simplex : simplices) {
1154 std::vector<int> splx = simplex;
1159 for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
1160 int vertex = it->first;
float val = it->second;
1161 int vert[] = {vertex};
int edge[] = {vertex, -2};
1175 for (
int i = 0; i < max_dim; i++) {
1177 int num_bars = bars.size();
if(i == 0) num_bars -= 1;
1178 if(verbose) std::clog << num_bars <<
" interval(s) in dimension " << i <<
":" << std::endl;
1179 for (
int j = 0; j < num_bars; j++) {
1180 double birth = bars[j].first;
1181 double death = bars[j].second;
1182 if (i == 0 && std::isinf(death))
continue;
1184 birth = minf + (birth + 2) * (maxf - minf);
1186 birth = minf + (2 - birth) * (maxf - minf);
1188 death = minf + (death + 2) * (maxf - minf);
1190 death = minf + (2 - death) * (maxf - minf);
1191 PD.push_back(std::pair<double, double>(birth, death));
1192 if (verbose) std::clog <<
" [" << birth <<
", " << death <<
"]" << std::endl;
1205 unsigned int sz = distribution.size();
1207 for (
unsigned int i = 0; i < N - sz; i++) {
1208 if (verbose) std::clog <<
"Computing " << i <<
"th bootstrap, bottleneck distance = ";
1210 Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover =
true;
1212 std::vector<int> boot(this->n);
1213 for (
int j = 0; j < this->n; j++) {
1214 double u = GetUniform();
1215 int id = std::floor(u * (this->n)); boot[j] = id;
1216 Cboot.point_cloud.push_back(this->point_cloud[
id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[
id]);
1217 boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
1221 for (
int j = 0; j < n; j++) {
1222 std::vector<double> dist(n);
1223 for (
int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]];
1224 Cboot.distances.push_back(dist);
1234 if (verbose) std::clog << db << std::endl;
1235 distribution.push_back(db);
1238 std::sort(distribution.begin(), distribution.end());
1249 unsigned int N = distribution.size();
1250 double d = distribution[std::floor(alpha * N)];
1251 if (verbose) std::clog <<
"Distance corresponding to confidence " << alpha <<
" is " << d << std::endl;
1262 unsigned int N = distribution.size();
1264 for (
unsigned int i = 0; i < N; i++)
1265 if (distribution[i] >= d){ level = i * 1.0 / N;
break; }
1266 if (verbose) std::clog <<
"Confidence level of distance " << d <<
" is " << level << std::endl;
1276 double distancemin = (std::numeric_limits<double>::max)();
int N = PD.size();
1277 for (
int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first));
1279 if (verbose) std::clog <<
"p value = " << p_value << std::endl;
1293 template <
typename SimplicialComplex>
1295 unsigned int dimension = 0;
1296 for (
auto const& simplex : simplices) {
1297 int numvert = simplex.size();
1298 double filt = std::numeric_limits<double>::lowest();
1299 for (
int i = 0; i < numvert; i++) filt = (std::max)(cover_color[simplex[i]].second, filt);
1300 complex.insert_simplex_and_subfaces(simplex, filt);
1301 if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
1309 if (type !=
"Nerve" && type !=
"GIC") {
1310 std::cerr <<
"Type of complex needs to be specified." << std::endl;
1314 if (type ==
"Nerve") {
1315 for(
int i = 0; i < n; i++) simplices.push_back(cover[i]);
1316 std::sort(simplices.begin(), simplices.end());
1317 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1318 simplices.resize(std::distance(simplices.begin(), it));
1321 if (type ==
"GIC") {
1322 Index_map index = boost::get(boost::vertex_index, one_skeleton);
1324 if (functional_cover) {
1329 throw std::invalid_argument(
1330 "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
1333 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1334 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
1335 int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
1336 for (
int i = 0; i < nums; i++) {
1337 int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
1338 int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
1339 for (
int j = 0; j < numt; j++) {
1340 int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
1341 if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
1342 std::vector<int> edge(2);
1343 edge[0] = (std::min)(vs, vt);
1344 edge[1] = (std::max)(vs, vt);
1345 simplices.push_back(edge);
1352 std::sort(simplices.begin(), simplices.end());
1353 std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
1354 simplices.resize(std::distance(simplices.begin(), it));
1359 boost::graph_traits<Graph>::edge_iterator ei, ei_end;
1360 for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
1361 if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
1362 cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
1363 std::vector<int> edge(2);
1364 edge[0] = index[boost::source(*ei, one_skeleton)];
1365 edge[1] = index[boost::target(*ei, one_skeleton)];
1378 std::vector<int> simplx;
1380 unsigned int sz = cover[vertex].size();
1381 for (
unsigned int i = 0; i < sz; i++) {
1382 simplx.push_back(cover[vertex][i]);
1385 std::sort(simplx.begin(), simplx.end());
1386 std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
1387 simplx.resize(std::distance(simplx.begin(), it));
1388 simplices.push_back(simplx);
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));