summaryrefslogtreecommitdiff
path: root/src/Nerve_GIC/include/gudhi/GIC.h
diff options
context:
space:
mode:
Diffstat (limited to 'src/Nerve_GIC/include/gudhi/GIC.h')
-rw-r--r--src/Nerve_GIC/include/gudhi/GIC.h108
1 files changed, 53 insertions, 55 deletions
diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 2a6d4788..047fba61 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -17,6 +17,14 @@
#include <mutex>
#endif
+#if __has_include(<CGAL/version.h>)
+# define GUDHI_GIC_USE_CGAL 1
+# include <gudhi/Bottleneck.h>
+#elif __has_include(<hera/bottleneck.h>)
+# define GUDHI_GIC_USE_HERA 1
+# include <hera/bottleneck.h>
+#endif
+
#include <gudhi/Debug_utils.h>
#include <gudhi/graph_simplicial_complex.h>
#include <gudhi/reader_utils.h>
@@ -25,7 +33,6 @@
#include <gudhi/Points_off_io.h>
#include <gudhi/distance_functions.h>
#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Bottleneck.h>
#include <boost/config.hpp>
#include <boost/graph/graph_traits.hpp>
@@ -35,8 +42,6 @@
#include <boost/graph/subgraph.hpp>
#include <boost/graph/graph_utility.hpp>
-#include <CGAL/version.h> // for CGAL_VERSION_NR
-
#include <iostream>
#include <vector>
#include <map>
@@ -139,19 +144,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 +339,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 +402,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 +412,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 +426,7 @@ class Cover_complex {
}
}
output.close();
- if (verbose) std::cout << std::endl;
+ if (verbose) std::clog << std::endl;
}
}
@@ -451,14 +446,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 +480,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 +523,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 +556,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 +572,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 +582,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 +615,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 +630,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 +656,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 +674,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 +691,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 +708,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 +757,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 +793,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 +887,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 +909,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 +984,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 +1047,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 +1083,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 +1130,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 +1178,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 +1192,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 +1206,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;
@@ -1242,8 +1233,15 @@ class Cover_complex {
Cboot.set_cover_from_function();
Cboot.find_simplices();
Cboot.compute_PD();
+#ifdef GUDHI_GIC_USE_CGAL
double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD);
- if (verbose) std::cout << db << std::endl;
+#elif defined GUDHI_GIC_USE_HERA
+ double db = hera::bottleneckDistExact(this->PD, Cboot.PD);
+#else
+ double db;
+ throw std::logic_error("This function requires CGAL or Hera for the bottleneck distance.");
+#endif
+ if (verbose) std::clog << db << std::endl;
distribution.push_back(db);
}
@@ -1260,7 +1258,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 +1273,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 +1286,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 +1317,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;
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-06 22:49:11 +0000