From c524232f734de875d69e2f190f01a6c976024368 Mon Sep 17 00:00:00 2001 From: Gard Spreemann Date: Thu, 14 Jun 2018 20:39:01 +0200 Subject: GUDHI 2.2.0 as released by upstream in a tarball. --- include/gudhi/GIC.h | 333 +++++++++++++++++++++++++++++++--------------------- 1 file changed, 202 insertions(+), 131 deletions(-) (limited to 'include/gudhi/GIC.h') 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 +#include +#endif + #include #include #include @@ -99,9 +104,8 @@ class Cover_complex { int data_dimension; // dimension of input data. int n; // number of points. - std::map func; // function used to compute the output simplicial complex. - std::map - func_color; // function used to compute the colors of the nodes of the output simplicial complex. + std::vector func; // function used to compute the output simplicial complex. + std::vector 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 distribution; - std::map > - cover; // function associating to each data point its vectors of cover elements to which it belongs. + std::vector > + cover; // function associating to each data point the vector of cover elements to which it belongs. std::map > cover_back; // inverse of cover, in order to get the data points associated to a specific cover element. std::map 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 func) { Fct = func; } - std::map 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::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 Dist(0, 1); +#endif // GUDHI_CAN_USE_CXX11_THREAD_LOCAL + std::uniform_real_distribution 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 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 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 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 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 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 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 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 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 >::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::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 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 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 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 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 splx = simplex; - splx.push_back(-2); - st.insert_simplex_and_subfaces(splx); + std::vector 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::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::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 pcoh(st); - pcoh.init_coefficients(2); + Gudhi::persistent_cohomology::Persistent_cohomology 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 > 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 - 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 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 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::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 >::iterator it = std::unique(simplices.begin(), simplices.end()); simplices.resize(std::distance(simplices.begin(), it)); -- cgit v1.2.3