diff options
author | Gard Spreemann <gspr@nonempty.org> | 2020-05-20 08:46:13 +0200 |
---|---|---|
committer | Gard Spreemann <gspr@nonempty.org> | 2020-05-20 08:46:13 +0200 |
commit | 2e879a40ee1bda90526d87f131b77ed365f3676b (patch) | |
tree | a3da06477fe2e8f3e0aa5c0be16f2eecb5bcaf9d /src/Nerve_GIC/include/gudhi | |
parent | b521123d85ae166e01c814292370d63856251326 (diff) | |
parent | 9b3079646ee3f6a494b83e864b3e10b8a93597d0 (diff) |
Merge branch 'dfsg/latest' into debian/sid
Diffstat (limited to 'src/Nerve_GIC/include/gudhi')
-rw-r--r-- | src/Nerve_GIC/include/gudhi/GIC.h | 90 |
1 files changed, 38 insertions, 52 deletions
diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 2a6d4788..1b1f9323 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -139,19 +139,9 @@ class Cover_complex { for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G); } - // Thread local is not available on XCode version < V.8 - // If not available, random engine is a class member. -#ifndef GUDHI_CAN_USE_CXX11_THREAD_LOCAL - std::default_random_engine re; -#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL - // Find random number in [0,1]. double GetUniform() { - // Thread local is not available on XCode version < V.8 - // If available, random engine is defined for each thread. -#ifdef GUDHI_CAN_USE_CXX11_THREAD_LOCAL thread_local std::default_random_engine re; -#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL std::uniform_real_distribution<double> Dist(0, 1); return Dist(re); } @@ -344,7 +334,7 @@ class Cover_complex { if (num_edges(one_skeleton_OFF)) one_skeleton = one_skeleton_OFF; else - std::cout << "No triangulation read in OFF file!" << std::endl; + std::cerr << "No triangulation read in OFF file!" << std::endl; } public: // Set graph from Rips complex. @@ -407,7 +397,7 @@ class Cover_complex { std::ifstream input(distance, std::ios::out | std::ios::binary); if (input.good()) { - if (verbose) std::cout << "Reading distances..." << std::endl; + if (verbose) std::clog << "Reading distances..." << std::endl; for (int i = 0; i < n; i++) { for (int j = i; j < n; j++) { input.read((char*)&d, 8); @@ -417,12 +407,12 @@ class Cover_complex { } input.close(); } else { - if (verbose) std::cout << "Computing distances..." << std::endl; + if (verbose) std::clog << "Computing distances..." << std::endl; input.close(); std::ofstream output(distance, std::ios::out | std::ios::binary); for (int i = 0; i < n; i++) { int state = (int)floor(100 * (i * 1.0 + 1) / n) % 10; - if (state == 0 && verbose) std::cout << "\r" << state << "%" << std::flush; + if (state == 0 && verbose) std::clog << "\r" << state << "%" << std::flush; for (int j = i; j < n; j++) { double dis = ref_distance(point_cloud[i], point_cloud[j]); distances[i][j] = dis; @@ -431,7 +421,7 @@ class Cover_complex { } } output.close(); - if (verbose) std::cout << std::endl; + if (verbose) std::clog << std::endl; } } @@ -451,14 +441,12 @@ class Cover_complex { m = (std::min)(m, n - 1); double delta = 0; - if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl; - if (verbose) std::cout << "Subsampling " << m << " points" << std::endl; + if (verbose) std::clog << n << " points in R^" << data_dimension << std::endl; + if (verbose) std::clog << "Subsampling " << m << " points" << std::endl; if (distances.size() == 0) compute_pairwise_distances(distance); - // This cannot be parallelized if thread_local is not defined - // thread_local is not defined for XCode < v.8 - #if defined(GUDHI_USE_TBB) && defined(GUDHI_CAN_USE_CXX11_THREAD_LOCAL) + #ifdef GUDHI_USE_TBB std::mutex deltamutex; tbb::parallel_for(0, N, [&](int i){ std::vector<int> samples(m); @@ -487,7 +475,7 @@ class Cover_complex { } #endif - if (verbose) std::cout << "delta = " << delta << std::endl; + if (verbose) std::clog << "delta = " << delta << std::endl; set_graph_from_rips(delta, distance); return delta; } @@ -530,7 +518,7 @@ class Cover_complex { cover_name = "coordinate " + std::to_string(k); } else{ - std::cout << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl; + std::cerr << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl; for (int i = 0; i < n; i++) func.push_back(0.0); functional_cover = true; cover_name = "null"; @@ -563,11 +551,11 @@ class Cover_complex { */ double set_automatic_resolution() { if (!functional_cover) { - std::cout << "Cover needs to come from the preimages of a function." << std::endl; + std::cerr << "Cover needs to come from the preimages of a function." << std::endl; return 0; } if (type != "Nerve" && type != "GIC") { - std::cout << "Type of complex needs to be specified." << std::endl; + std::cerr << "Type of complex needs to be specified." << std::endl; return 0; } @@ -579,7 +567,7 @@ class Cover_complex { for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - func[index[boost::target(*ei, one_skeleton)]])); - if (verbose) std::cout << "resolution = " << reso << std::endl; + if (verbose) std::clog << "resolution = " << reso << std::endl; resolution_double = reso; } @@ -589,7 +577,7 @@ class Cover_complex { reso = (std::max)(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] - func[index[boost::target(*ei, one_skeleton)]]) / gain); - if (verbose) std::cout << "resolution = " << reso << std::endl; + if (verbose) std::clog << "resolution = " << reso << std::endl; resolution_double = reso; } @@ -622,11 +610,11 @@ class Cover_complex { */ void set_cover_from_function() { if (resolution_double == -1 && resolution_int == -1) { - std::cout << "Number and/or length of intervals not specified" << std::endl; + std::cerr << "Number and/or length of intervals not specified" << std::endl; return; } if (gain == -1) { - std::cout << "Gain not specified" << std::endl; + std::cerr << "Gain not specified" << std::endl; return; } @@ -637,7 +625,7 @@ class Cover_complex { minf = (std::min)(minf, func[i]); maxf = (std::max)(maxf, func[i]); } - if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; + if (verbose) std::clog << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; // Compute cover of im(f) std::vector<std::pair<double, double> > intervals; @@ -663,7 +651,7 @@ class Cover_complex { res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; } } else { @@ -681,7 +669,7 @@ class Cover_complex { res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; } } else { // Case we use an integer and a double for the length of the intervals. @@ -698,7 +686,7 @@ class Cover_complex { res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" + std::clog << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; } } @@ -715,7 +703,7 @@ class Cover_complex { std::map<int, std::vector<int> > preimages; std::map<int, double> funcstd; - if (verbose) std::cout << "Computing preimages..." << std::endl; + if (verbose) std::clog << "Computing preimages..." << std::endl; for (int i = 0; i < res; i++) { // Find points in the preimage std::pair<double, double> inter1 = intervals[i]; @@ -764,7 +752,7 @@ class Cover_complex { } #ifdef GUDHI_USE_TBB - if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl; + if (verbose) std::clog << "Computing connected components (parallelized)..." << std::endl; std::mutex covermutex, idmutex; tbb::parallel_for(0, res, [&](int i){ // Compute connected components @@ -800,7 +788,7 @@ class Cover_complex { idmutex.unlock(); }); #else - if (verbose) std::cout << "Computing connected components..." << std::endl; + if (verbose) std::clog << "Computing connected components..." << std::endl; for (int i = 0; i < res; i++) { // Compute connected components Graph G = one_skeleton.create_subgraph(); @@ -894,7 +882,7 @@ class Cover_complex { // Compute the geodesic distances to subsamples with Dijkstra #ifdef GUDHI_USE_TBB - if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl; + if (verbose) std::clog << "Computing geodesic distances (parallelized)..." << std::endl; std::mutex coverMutex; std::mutex mindistMutex; tbb::parallel_for(0, m, [&](int i){ int seed = voronoi_subsamples[i]; @@ -916,7 +904,7 @@ class Cover_complex { }); #else for (int i = 0; i < m; i++) { - if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; + if (verbose) std::clog << "Computing geodesic distances to seed " << i << "..." << std::endl; int seed = voronoi_subsamples[i]; std::vector<double> dmap(n); boost::dijkstra_shortest_paths( @@ -991,7 +979,7 @@ class Cover_complex { color_name.append(std::to_string(k)); } else{ - std::cout << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl; + std::cerr << "Only pairwise distances provided---cannot access " << k << "th coordinate; returning null vector instead" << std::endl; for (int i = 0; i < n; i++) func.push_back(0.0); functional_cover = true; cover_name = "null"; @@ -1054,7 +1042,7 @@ class Cover_complex { } graphic << "}"; graphic.close(); - std::cout << mapp << " file generated. It can be visualized with e.g. neato." << std::endl; + std::clog << mapp << " file generated. It can be visualized with e.g. neato." << std::endl; } public: // Create a .txt file that can be compiled with KeplerMapper. @@ -1090,7 +1078,7 @@ class Cover_complex { if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl; graphic.close(); - std::cout << mapp + std::clog << mapp << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; } @@ -1137,7 +1125,7 @@ class Cover_complex { for (int i = 0; i < numfaces; i++) graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl; graphic.close(); - std::cout << mapp << " generated. It can be visualized with e.g. geomview." << std::endl; + std::clog << mapp << " generated. It can be visualized with e.g. geomview." << std::endl; } // ******************************************************************************************************************* @@ -1185,7 +1173,7 @@ class Cover_complex { for (int i = 0; i < max_dim; i++) { std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i); int num_bars = bars.size(); if(i == 0) num_bars -= 1; - if(verbose) std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl; + if(verbose) std::clog << num_bars << " interval(s) in dimension " << i << ":" << std::endl; for (int j = 0; j < num_bars; j++) { double birth = bars[j].first; double death = bars[j].second; @@ -1199,7 +1187,7 @@ class Cover_complex { else death = minf + (2 - death) * (maxf - minf); PD.push_back(std::pair<double, double>(birth, death)); - if (verbose) std::cout << " [" << birth << ", " << death << "]" << std::endl; + if (verbose) std::clog << " [" << birth << ", " << death << "]" << std::endl; } } return PD; @@ -1213,11 +1201,9 @@ class Cover_complex { */ void compute_distribution(unsigned int N = 100) { unsigned int sz = distribution.size(); - if (sz >= N) { - std::cout << "Already done!" << std::endl; - } else { + if (sz < N) { for (unsigned int i = 0; i < N - sz; i++) { - if (verbose) std::cout << "Computing " << i << "th bootstrap, bottleneck distance = "; + if (verbose) std::clog << "Computing " << i << "th bootstrap, bottleneck distance = "; Cover_complex Cboot; Cboot.n = this->n; Cboot.data_dimension = this->data_dimension; Cboot.type = this->type; Cboot.functional_cover = true; @@ -1243,7 +1229,7 @@ class Cover_complex { Cboot.find_simplices(); Cboot.compute_PD(); double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD); - if (verbose) std::cout << db << std::endl; + if (verbose) std::clog << db << std::endl; distribution.push_back(db); } @@ -1260,7 +1246,7 @@ class Cover_complex { double compute_distance_from_confidence_level(double alpha) { unsigned int N = distribution.size(); double d = distribution[std::floor(alpha * N)]; - if (verbose) std::cout << "Distance corresponding to confidence " << alpha << " is " << d << std::endl; + if (verbose) std::clog << "Distance corresponding to confidence " << alpha << " is " << d << std::endl; return d; } @@ -1275,7 +1261,7 @@ class Cover_complex { double level = 1; for (unsigned int i = 0; i < N; i++) if (distribution[i] >= d){ level = i * 1.0 / N; break; } - if (verbose) std::cout << "Confidence level of distance " << d << " is " << level << std::endl; + if (verbose) std::clog << "Confidence level of distance " << d << " is " << level << std::endl; return level; } @@ -1288,7 +1274,7 @@ class Cover_complex { double distancemin = (std::numeric_limits<double>::max)(); int N = PD.size(); for (int i = 0; i < N; i++) distancemin = (std::min)(distancemin, 0.5 * std::abs(PD[i].second - PD[i].first)); double p_value = 1 - compute_confidence_level_from_distance(distancemin); - if (verbose) std::cout << "p value = " << p_value << std::endl; + if (verbose) std::clog << "p value = " << p_value << std::endl; return p_value; } @@ -1319,7 +1305,7 @@ class Cover_complex { */ void find_simplices() { if (type != "Nerve" && type != "GIC") { - std::cout << "Type of complex needs to be specified." << std::endl; + std::cerr << "Type of complex needs to be specified." << std::endl; return; } |