summaryrefslogtreecommitdiff
path: root/include/gudhi/GIC.h
diff options
context:
space:
mode:
Diffstat (limited to 'include/gudhi/GIC.h')
-rw-r--r--include/gudhi/GIC.h333
1 files changed, 202 insertions, 131 deletions
diff --git a/include/gudhi/GIC.h b/include/gudhi/GIC.h
index 40ff7a4a..7aa95210 100644
--- a/include/gudhi/GIC.h
+++ b/include/gudhi/GIC.h
@@ -4,7 +4,7 @@
*
* Author: Mathieu Carriere
*
- * Copyright (C) 2017 INRIA
+ * Copyright (C) 2017 Inria
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
@@ -23,6 +23,11 @@
#ifndef GIC_H_
#define GIC_H_
+#ifdef GUDHI_USE_TBB
+#include <tbb/parallel_for.h>
+#include <tbb/mutex.h>
+#endif
+
#include <gudhi/Debug_utils.h>
#include <gudhi/graph_simplicial_complex.h>
#include <gudhi/reader_utils.h>
@@ -99,9 +104,8 @@ class Cover_complex {
int data_dimension; // dimension of input data.
int n; // number of points.
- std::map<int, double> func; // function used to compute the output simplicial complex.
- std::map<int, double>
- func_color; // function used to compute the colors of the nodes of the output simplicial complex.
+ std::vector<double> func; // function used to compute the output simplicial complex.
+ std::vector<double> func_color; // function used to compute the colors of the nodes of the output simplicial complex.
bool functional_cover = false; // whether we use a cover with preimages of a function or not.
Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists).
@@ -114,8 +118,8 @@ class Cover_complex {
Persistence_diagram PD;
std::vector<double> distribution;
- std::map<int, std::vector<int> >
- cover; // function associating to each data point its vectors of cover elements to which it belongs.
+ std::vector<std::vector<int> >
+ cover; // function associating to each data point the vector of cover elements to which it belongs.
std::map<int, std::vector<int> >
cover_back; // inverse of cover, in order to get the data points associated to a specific cover element.
std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence
@@ -138,28 +142,26 @@ class Cover_complex {
std::string point_cloud_name;
std::string color_name;
- // Point comparator
- struct Less {
- Less(std::map<int, double> func) { Fct = func; }
- std::map<int, double> Fct;
- bool operator()(int a, int b) {
- if (Fct[a] == Fct[b])
- return a < b;
- else
- return Fct[a] < Fct[b];
- }
- };
-
// Remove all edges of a graph.
void remove_edges(Graph& G) {
boost::graph_traits<Graph>::edge_iterator ei, ei_end;
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;
- thread_local std::uniform_real_distribution<double> Dist(0, 1);
+#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL
+ std::uniform_real_distribution<double> Dist(0, 1);
return Dist(re);
}
@@ -276,6 +278,7 @@ class Cover_complex {
point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
boost::add_vertex(one_skeleton_OFF);
vertices.push_back(boost::add_vertex(one_skeleton));
+ cover.emplace_back();
i++;
}
}
@@ -422,7 +425,6 @@ class Cover_complex {
double set_graph_from_automatic_rips(Distance distance, int N = 100) {
int m = floor(n / std::exp((1 + rate_power) * std::log(std::log(n) / std::log(rate_constant))));
m = std::min(m, n - 1);
- std::vector<int> samples(m);
double delta = 0;
if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl;
@@ -430,17 +432,36 @@ class Cover_complex {
if (distances.size() == 0) compute_pairwise_distances(distance);
- // #pragma omp parallel for
- for (int i = 0; i < N; i++) {
- SampleWithoutReplacement(n, m, samples);
- double hausdorff_dist = 0;
- for (int j = 0; j < n; j++) {
- double mj = distances[j][samples[0]];
- for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
- hausdorff_dist = std::max(hausdorff_dist, mj);
+ // 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)
+ tbb::mutex deltamutex;
+ tbb::parallel_for(0, N, [&](int i){
+ std::vector<int> samples(m);
+ SampleWithoutReplacement(n, m, samples);
+ double hausdorff_dist = 0;
+ for (int j = 0; j < n; j++) {
+ double mj = distances[j][samples[0]];
+ for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
+ hausdorff_dist = std::max(hausdorff_dist, mj);
+ }
+ deltamutex.lock();
+ delta += hausdorff_dist / N;
+ deltamutex.unlock();
+ });
+ #else
+ for (int i = 0; i < N; i++) {
+ std::vector<int> samples(m);
+ SampleWithoutReplacement(n, m, samples);
+ double hausdorff_dist = 0;
+ for (int j = 0; j < n; j++) {
+ double mj = distances[j][samples[0]];
+ for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]);
+ hausdorff_dist = std::max(hausdorff_dist, mj);
+ }
+ delta += hausdorff_dist / N;
}
- delta += hausdorff_dist / N;
- }
+ #endif
if (verbose) std::cout << "delta = " << delta << std::endl;
set_graph_from_rips(delta, distance);
@@ -465,7 +486,7 @@ class Cover_complex {
while (std::getline(input, line)) {
std::stringstream stream(line);
stream >> f;
- func.emplace(i, f);
+ func.push_back(f);
i++;
}
functional_cover = true;
@@ -479,7 +500,7 @@ class Cover_complex {
*
*/
void set_function_from_coordinate(int k) {
- for (int i = 0; i < n; i++) func.emplace(i, point_cloud[i][k]);
+ for (int i = 0; i < n; i++) func.push_back(point_cloud[i][k]);
functional_cover = true;
cover_name = "coordinate " + std::to_string(k);
}
@@ -492,7 +513,7 @@ class Cover_complex {
*/
template <class InputRange>
void set_function_from_range(InputRange const& function) {
- for (int i = 0; i < n; i++) func.emplace(i, function[i]);
+ for (int i = 0; i < n; i++) func.push_back(function[i]);
functional_cover = true;
}
@@ -654,7 +675,7 @@ class Cover_complex {
// Sort points according to function values
std::vector<int> points(n);
for (int i = 0; i < n; i++) points[i] = i;
- std::sort(points.begin(), points.end(), Less(this->func));
+ std::sort(points.begin(), points.end(), [=](const int & p1, const int & p2){return (this->func[p1] < this->func[p2]);});
int id = 0;
int pos = 0;
@@ -710,37 +731,74 @@ class Cover_complex {
funcstd[i] = 0.5 * (u + v);
}
- if (verbose) std::cout << "Computing connected components..." << std::endl;
- // #pragma omp parallel for
- for (int i = 0; i < res; i++) {
- // Compute connected components
- Graph G = one_skeleton.create_subgraph();
- int num = preimages[i].size();
- std::vector<int> component(num);
- for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
- boost::connected_components(G, &component[0]);
- int max = 0;
-
- // For each point in preimage
- for (int j = 0; j < num; j++) {
- // Update number of components in preimage
- if (component[j] > max) max = component[j];
-
- // Identify component with Cantor polynomial N^2 -> N
- int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
-
- // Update covers
- cover[preimages[i][j]].push_back(identifier);
- cover_back[identifier].push_back(preimages[i][j]);
- cover_fct[identifier] = i;
- cover_std[identifier] = funcstd[i];
- cover_color[identifier].second += func_color[preimages[i][j]];
- cover_color[identifier].first += 1;
- }
+ #ifdef GUDHI_USE_TBB
+ if (verbose) std::cout << "Computing connected components (parallelized)..." << std::endl;
+ tbb::mutex covermutex, idmutex;
+ tbb::parallel_for(0, res, [&](int i){
+ // Compute connected components
+ Graph G = one_skeleton.create_subgraph();
+ int num = preimages[i].size();
+ std::vector<int> component(num);
+ for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
+ boost::connected_components(G, &component[0]);
+ int max = 0;
+
+ // For each point in preimage
+ for (int j = 0; j < num; j++) {
+ // Update number of components in preimage
+ if (component[j] > max) max = component[j];
+
+ // Identify component with Cantor polynomial N^2 -> N
+ int identifier = ((i + component[j])*(i + component[j]) + 3 * i + component[j]) / 2;
+
+ // Update covers
+ covermutex.lock();
+ cover[preimages[i][j]].push_back(identifier);
+ cover_back[identifier].push_back(preimages[i][j]);
+ cover_fct[identifier] = i;
+ cover_std[identifier] = funcstd[i];
+ cover_color[identifier].second += func_color[preimages[i][j]];
+ cover_color[identifier].first += 1;
+ covermutex.unlock();
+ }
- // Maximal dimension is total number of connected components
- id += max + 1;
- }
+ // Maximal dimension is total number of connected components
+ idmutex.lock();
+ id += max + 1;
+ idmutex.unlock();
+ });
+ #else
+ if (verbose) std::cout << "Computing connected components..." << std::endl;
+ for (int i = 0; i < res; i++) {
+ // Compute connected components
+ Graph G = one_skeleton.create_subgraph();
+ int num = preimages[i].size();
+ std::vector<int> component(num);
+ for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
+ boost::connected_components(G, &component[0]);
+ int max = 0;
+
+ // For each point in preimage
+ for (int j = 0; j < num; j++) {
+ // Update number of components in preimage
+ if (component[j] > max) max = component[j];
+
+ // Identify component with Cantor polynomial N^2 -> N
+ int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
+
+ // Update covers
+ cover[preimages[i][j]].push_back(identifier);
+ cover_back[identifier].push_back(preimages[i][j]);
+ cover_fct[identifier] = i;
+ cover_std[identifier] = funcstd[i];
+ cover_color[identifier].second += func_color[preimages[i][j]];
+ cover_color[identifier].first += 1;
+ }
+
+ // Maximal dimension is total number of connected components
+ id += max + 1;
+ }
+ #endif
maximal_dim = id - 1;
for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
@@ -803,24 +861,46 @@ class Cover_complex {
for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max();
// Compute the geodesic distances to subsamples with Dijkstra
- // #pragma omp parallel for
- for (int i = 0; i < m; i++) {
- if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
- int seed = voronoi_subsamples[i];
- std::vector<double> dmap(n);
- boost::dijkstra_shortest_paths(
- one_skeleton, vertices[seed],
- boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
-
- for (int j = 0; j < n; j++)
- if (mindist[j] > dmap[j]) {
- mindist[j] = dmap[j];
- if (cover[j].size() == 0)
- cover[j].push_back(i);
- else
- cover[j][0] = i;
- }
- }
+ #ifdef GUDHI_USE_TBB
+ if (verbose) std::cout << "Computing geodesic distances (parallelized)..." << std::endl;
+ tbb::mutex coverMutex; tbb::mutex mindistMutex;
+ tbb::parallel_for(0, m, [&](int i){
+ int seed = voronoi_subsamples[i];
+ std::vector<double> dmap(n);
+ boost::dijkstra_shortest_paths(
+ one_skeleton, vertices[seed],
+ boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
+
+ coverMutex.lock(); mindistMutex.lock();
+ for (int j = 0; j < n; j++)
+ if (mindist[j] > dmap[j]) {
+ mindist[j] = dmap[j];
+ if (cover[j].size() == 0)
+ cover[j].push_back(i);
+ else
+ cover[j][0] = i;
+ }
+ coverMutex.unlock(); mindistMutex.unlock();
+ });
+ #else
+ for (int i = 0; i < m; i++) {
+ if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
+ int seed = voronoi_subsamples[i];
+ std::vector<double> dmap(n);
+ boost::dijkstra_shortest_paths(
+ one_skeleton, vertices[seed],
+ boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
+
+ for (int j = 0; j < n; j++)
+ if (mindist[j] > dmap[j]) {
+ mindist[j] = dmap[j];
+ if (cover[j].size() == 0)
+ cover[j].push_back(i);
+ else
+ cover[j][0] = i;
+ }
+ }
+ #endif
for (int i = 0; i < n; i++) {
cover_back[cover[i][0]].push_back(i);
@@ -860,7 +940,7 @@ class Cover_complex {
while (std::getline(input, line)) {
std::stringstream stream(line);
stream >> f;
- func_color.emplace(i, f);
+ func_color.push_back(f);
i++;
}
color_name = color_file_name;
@@ -873,7 +953,7 @@ class Cover_complex {
*
*/
void set_color_from_coordinate(int k = 0) {
- for (int i = 0; i < n; i++) func_color[i] = point_cloud[i][k];
+ for (int i = 0; i < n; i++) func_color.push_back(point_cloud[i][k]);
color_name = "coordinate ";
color_name.append(std::to_string(k));
}
@@ -885,7 +965,7 @@ class Cover_complex {
*
*/
void set_color_from_vector(std::vector<double> color) {
- for (unsigned int i = 0; i < color.size(); i++) func_color[i] = color[i];
+ for (unsigned int i = 0; i < color.size(); i++) func_color.push_back(color[i]);
}
public: // Create a .dot file that can be compiled with neato to produce a .pdf file.
@@ -1039,45 +1119,29 @@ class Cover_complex {
minf = std::min(minf, it->second);
}
+ // Build filtration
for (auto const& simplex : simplices) {
- // Add a simplex and a cone on it
- std::vector<int> splx = simplex;
- splx.push_back(-2);
- st.insert_simplex_and_subfaces(splx);
+ std::vector<int> splx = simplex; splx.push_back(-2);
+ st.insert_simplex_and_subfaces(splx, -3);
}
- // Build filtration
- for (auto simplex : st.complex_simplex_range()) {
- double filta = std::numeric_limits<double>::lowest();
- double filts = filta;
- bool ascending = true;
- for (auto vertex : st.simplex_vertex_range(simplex)) {
- if (vertex == -2) {
- ascending = false;
- continue;
- }
- filta = std::max(-2 + (cover_std[vertex] - minf) / (maxf - minf), filta);
- filts = std::max(2 - (cover_std[vertex] - minf) / (maxf - minf), filts);
- }
- if (ascending)
- st.assign_filtration(simplex, filta);
- else
- st.assign_filtration(simplex, filts);
+ for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
+ int vertex = it->first; float val = it->second;
+ int vert[] = {vertex}; int edge[] = {vertex, -2};
+ st.assign_filtration(st.find(vert), -2 + (val - minf)/(maxf - minf));
+ st.assign_filtration(st.find(edge), 2 - (val - minf)/(maxf - minf));
}
- int magic[] = {-2};
- st.assign_filtration(st.find(magic), -3);
+ st.make_filtration_non_decreasing();
// Compute PD
- st.initialize_filtration();
- Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st);
- pcoh.init_coefficients(2);
+ Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); pcoh.init_coefficients(2);
pcoh.compute_persistent_cohomology();
// Output PD
int max_dim = st.dimension();
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();
+ int num_bars = bars.size(); if(i == 0) num_bars -= 1;
if(verbose) std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
for (int j = 0; j < num_bars; j++) {
double birth = bars[j].first;
@@ -1103,22 +1167,25 @@ class Cover_complex {
* @param[in] N number of bootstrap iterations.
*
*/
- template <typename SimplicialComplex>
- void compute_distribution(int N = 100) {
- if (distribution.size() >= N) {
+ void compute_distribution(unsigned int N = 100) {
+ unsigned int sz = distribution.size();
+ if (sz >= N) {
std::cout << "Already done!" << std::endl;
} else {
- for (int i = 0; i < N - distribution.size(); i++) {
- Cover_complex Cboot;
- Cboot.n = this->n;
+ for (unsigned int i = 0; i < N - sz; i++) {
+ if (verbose) std::cout << "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;
+
std::vector<int> boot(this->n);
for (int j = 0; j < this->n; j++) {
double u = GetUniform();
- int id = std::floor(u * (this->n));
- boot[j] = id;
- Cboot.point_cloud[j] = this->point_cloud[id];
- Cboot.func.emplace(j, this->func[id]);
+ int id = std::floor(u * (this->n)); boot[j] = id;
+ Cboot.point_cloud.push_back(this->point_cloud[id]); Cboot.cover.emplace_back(); Cboot.func.push_back(this->func[id]);
+ boost::add_vertex(Cboot.one_skeleton_OFF); Cboot.vertices.push_back(boost::add_vertex(Cboot.one_skeleton));
}
+ Cboot.set_color_from_vector(Cboot.func);
+
for (int j = 0; j < n; j++) {
std::vector<double> dist(n);
for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]];
@@ -1131,8 +1198,9 @@ class Cover_complex {
Cboot.set_cover_from_function();
Cboot.find_simplices();
Cboot.compute_PD();
-
- distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD));
+ double db = Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD);
+ if (verbose) std::cout << db << std::endl;
+ distribution.push_back(db);
}
std::sort(distribution.begin(), distribution.end());
@@ -1146,7 +1214,7 @@ class Cover_complex {
*
*/
double compute_distance_from_confidence_level(double alpha) {
- int N = distribution.size();
+ unsigned int N = distribution.size();
return distribution[std::floor(alpha * N)];
}
@@ -1157,9 +1225,11 @@ class Cover_complex {
*
*/
double compute_confidence_level_from_distance(double d) {
- int N = distribution.size();
- for (int i = 0; i < N; i++)
- if (distribution[i] > d) return i * 1.0 / N;
+ unsigned int N = distribution.size();
+ double level = 1;
+ for (unsigned int i = 0; i < N; i++)
+ if (distribution[i] > d){ level = i * 1.0 / N; break; }
+ return level;
}
public:
@@ -1171,7 +1241,9 @@ class Cover_complex {
double distancemin = -std::numeric_limits<double>::lowest();
int N = PD.size();
for (int i = 0; i < N; i++) distancemin = std::min(distancemin, 0.5 * (PD[i].second - PD[i].first));
- return 1 - compute_confidence_level_from_distance(distancemin);
+ double p_value = 1 - compute_confidence_level_from_distance(distancemin);
+ if (verbose) std::cout << "p value = " << p_value << std::endl;
+ return p_value;
}
// *******************************************************************************************************************
@@ -1206,8 +1278,7 @@ class Cover_complex {
}
if (type == "Nerve") {
- for(auto& simplex : cover)
- simplices.push_back(simplex.second);
+ for(int i = 0; i < n; i++) simplices.push_back(cover[i]);
std::sort(simplices.begin(), simplices.end());
std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
simplices.resize(std::distance(simplices.begin(), it));
generated by cgit v1.2.3 (git 2.39.1) at 2024-04-27 19:03:09 +0000