From f4bacd6ca6db4ef85a030cd505715174e4db6f6d Mon Sep 17 00:00:00 2001 From: mcarrier Date: Tue, 19 Dec 2017 14:06:48 +0000 Subject: changed data structure to use boost graphs git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3085 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 614b6362ceeb8237c224d12b523677dcbc82fba8 --- src/Nerve_GIC/doc/Intro_graph_induced_complex.h | 2 +- src/Nerve_GIC/example/FuncGIC.cpp | 1 + src/Nerve_GIC/example/Nerve.cpp | 15 +- src/Nerve_GIC/include/gudhi/GIC.h | 1077 +++++++++++++---------- 4 files changed, 627 insertions(+), 468 deletions(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h index 3a0d8154..6668d850 100644 --- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h +++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h @@ -169,7 +169,7 @@ namespace cover_complex { * * When launching: * - * \code $> ./FuncGIC ../../data/points/COIL_database/lucky_cat.off ../../data/points/COIL_database/lucky_cat_PCA1 --v + * \code $> ./FuncGIC ../../../../data/points/COIL_database/lucky_cat.off ../../../../data/points/COIL_database/lucky_cat_PCA1 --v * \endcode * * the program outputs again SC.dot which gives the following visualization after using neato: diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp index 3762db4e..3583c66f 100644 --- a/src/Nerve_GIC/example/FuncGIC.cpp +++ b/src/Nerve_GIC/example/FuncGIC.cpp @@ -71,6 +71,7 @@ int main(int argc, char **argv) { Gudhi::Simplex_tree<> stree; GIC.create_complex(stree); + GIC.compute_PD >(); // -------------------------------------------- // Display information about the functional GIC diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp index 4d5b009b..7634c6f4 100644 --- a/src/Nerve_GIC/example/Nerve.cpp +++ b/src/Nerve_GIC/example/Nerve.cpp @@ -25,19 +25,21 @@ #include #include +using namespace std; + void usage(int nbArgs, char *const progName) { - std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; - std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [--v] \n"; - std::cerr << " i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 --v \n"; + cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; + cerr << "Usage: " << progName << " filename.off coordinate resolution gain [--v] \n"; + cerr << " i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 --v \n"; exit(-1); // ----- >> } int main(int argc, char **argv) { if ((argc != 5) && (argc != 6)) usage(argc, argv[0]); - using Point = std::vector; + using Point = vector; - std::string off_file_name(argv[1]); + string off_file_name(argv[1]); int coord = atoi(argv[2]); int resolution = atoi(argv[3]); double gain = atof(argv[4]); @@ -54,7 +56,7 @@ int main(int argc, char **argv) { bool check = SC.read_point_cloud(off_file_name); if (!check) { - std::cout << "Incorrect OFF file." << std::endl; + cout << "Incorrect OFF file." << endl; } else { SC.set_type("Nerve"); @@ -72,6 +74,7 @@ int main(int argc, char **argv) { Gudhi::Simplex_tree<> stree; SC.create_complex(stree); + SC.compute_PD >(); // ---------------------------------------------------------------------------- // Display information about the graph induced complex diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 9f107a7e..d8c6abd1 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -30,6 +30,16 @@ #include #include #include +#include +#include + +#include +#include +#include +#include +#include +#include +#include #include #include @@ -41,13 +51,21 @@ #include #include +using namespace boost; +using namespace std; + namespace Gudhi { namespace cover_complex { -using Simplex_tree = Gudhi::Simplex_tree<>; -using Filtration_value = Simplex_tree::Filtration_value; -using Rips_complex = Gudhi::rips_complex::Rips_complex; +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex; +using PersistenceDiagram = vector >; +using Graph = subgraph > > >; +using vertex_t = graph_traits::vertex_descriptor; +using IndexMap = property_map::type; +using WeightMap = property_map::type; /** * \class Cover_complex @@ -72,25 +90,35 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex; template class Cover_complex { private: - // Graph_induced_complex(std::map fun){func = fun;} + bool verbose = false; // whether to display information. - std::vector point_cloud; - std::vector > one_skeleton; - typedef int Cover_t; // elements of cover C are indexed by integers. - std::vector > simplices; - std::map > cover; - std::map > cover_back; + + vector point_cloud; int maximal_dim; // maximal dimension of output simplicial complex. int data_dimension; // dimension of input data. int n; // number of points. - std::map - cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. - std::map > - cover_color; // size and coloring of the vertices of the output simplicial complex. - Simplex_tree st; - std::map > adjacency_matrix; - std::vector > distances; + vector > distances; + + map func; // function used to compute the output simplicial complex. + map 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). + Graph one_skeleton; // one-skeleton used to compute the connected components. + vector vertices; + vector > simplices; + + vector voronoi_subsamples; + + PersistenceDiagram PD; + vector distribution; + + map > cover; + map > cover_back; + map cover_std; // standard function (induced by func) used to compute the extended persistence diagram of the output simplicial complex. + map cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. + map > cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex. int resolution_int = -1; double resolution_double = -1; @@ -99,19 +127,15 @@ class Cover_complex { double rate_power = 0.001; // Power in the subsampling. int mask = 0; // Ignore nodes containing less than mask points. - std::map func; - std::map func_color; - std::vector voronoi_subsamples; - std::string cover_name; - std::string point_cloud_name; - std::string color_name; - std::string type; // Nerve or GIC - bool functional_cover = false; // whether we use a cover with preimages of a function or not + string cover_name; + string point_cloud_name; + string color_name; + string type; // Nerve or GIC // Point comparator struct Less { - Less(std::map func) { Fct = func; } - std::map Fct; + Less(map func) { Fct = func; } + map Fct; bool operator()(int a, int b) { if (Fct[a] == Fct[b]) return a < b; @@ -120,54 +144,23 @@ class Cover_complex { } }; - // DFS - private: - void dfs(std::map >& G, int p, std::vector& cc, std::map& visit) { - cc.push_back(p); - visit[p] = true; - int neighb = G[p].size(); - for (int i = 0; i < neighb; i++) - if (visit.find(G[p][i]) != visit.end()) - if (!(visit[G[p][i]])) dfs(G, G[p][i], cc, visit); - } - // Find random number in [0,1]. double GetUniform() { - static std::default_random_engine re; - static std::uniform_real_distribution Dist(0, 1); + thread_local default_random_engine re; + thread_local uniform_real_distribution Dist(0, 1); return Dist(re); } // Subsample points. - void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector& samples) { - int t = 0; - int m = 0; - double u; - while (m < sampleSize) { + void SampleWithoutReplacement(int populationSize, int sampleSize, vector & samples) { + int t = 0; int m = 0; double u; + while (m < sampleSize){ u = GetUniform(); - if ((populationSize - t) * u >= sampleSize - m) { - t++; - } else { - samples[m] = t; - t++; - m++; - } + if ((populationSize - t) * u >= sampleSize - m) t++; + else{ samples[m] = t; t++; m++;} } } - private: - void fill_adjacency_matrix_from_st() { - std::vector empty; - for (int i = 0; i < n; i++) adjacency_matrix[i] = empty; - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector vertices; - for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); - adjacency_matrix[vertices[0]].push_back(vertices[1]); - adjacency_matrix[vertices[1]].push_back(vertices[0]); - } - } - } public: /** \brief Specifies whether the type of the output simplicial complex. @@ -214,36 +207,36 @@ class Cover_complex { * @param[in] off_file_name name of the input .OFF or .nOFF file. * */ - bool read_point_cloud(const std::string& off_file_name) { + bool read_point_cloud(const string & off_file_name) { point_cloud_name = off_file_name; - std::ifstream input(off_file_name); - std::string line; + ifstream input(off_file_name); + string line; char comment = '#'; while (comment == '#') { getline(input, line); - if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) comment = line[line.find_first_not_of(' ')]; + if (!line.empty() && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } - if (std::strcmp((char*)line.c_str(), "nOFF") == 0) { + if (strcmp((char*)line.c_str(), "nOFF") == 0) { comment = '#'; while (comment == '#') { getline(input, line); - if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) + if (!line.empty() && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } - std::stringstream stream(line); + stringstream stream(line); stream >> data_dimension; } else { data_dimension = 3; } comment = '#'; - int numedges, numfaces, i, num; + int numedges, numfaces, i, dim; while (comment == '#') { getline(input, line); - if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) comment = line[line.find_first_not_of(' ')]; + if (!line.empty() && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } - std::stringstream stream(line); + stringstream stream(line); stream >> n; stream >> numfaces; stream >> numedges; @@ -251,12 +244,10 @@ class Cover_complex { i = 0; while (i < n) { getline(input, line); - if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && - !std::all_of(line.begin(), line.end(), isspace)) { - std::vector point; - std::istringstream iss(line); - point.assign(std::istream_iterator(iss), std::istream_iterator()); + if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) { + istringstream iss(line); vector point; point.assign(istream_iterator(iss), istream_iterator()); point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); + add_vertex(one_skeleton_OFF); vertices.push_back(add_vertex(one_skeleton)); i++; } } @@ -264,20 +255,12 @@ class Cover_complex { i = 0; while (i < numfaces) { getline(input, line); - if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && - !std::all_of(line.begin(), line.end(), isspace)) { - std::vector simplex; - std::istringstream iss(line); - simplex.assign(std::istream_iterator(iss), std::istream_iterator()); - num = simplex[0]; - std::vector edge(2); - for (int j = 1; j <= num; j++) { - for (int k = j + 1; k <= num; k++) { - edge[0] = simplex[j]; - edge[1] = simplex[k]; - one_skeleton.push_back(edge); - } - } + if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) { + vector simplex; istringstream iss(line); + simplex.assign(istream_iterator(iss), istream_iterator()); dim = simplex[0]; + for (int j = 1; j <= dim; j++) + for (int k = j + 1; k <= dim; k++) + add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF); i++; } } @@ -297,23 +280,12 @@ class Cover_complex { * each edge being represented by the IDs of its two nodes. * */ - void set_graph_from_file(const std::string& graph_file_name) { - int neighb; - std::ifstream input(graph_file_name); - std::string line; - int edge[2]; - int n = 0; - while (std::getline(input, line)) { - std::stringstream stream(line); - stream >> edge[0]; - while (stream >> neighb) { - edge[1] = neighb; - st.insert_simplex_and_subfaces(edge); - } - n++; + void set_graph_from_file(const string & graph_file_name){ + int neighb; ifstream input(graph_file_name); string line; int source; + while (getline(input, line)){ + stringstream stream(line); stream >> source; + while (stream >> neighb) add_edge(vertices[source], vertices[neighb], one_skeleton); } - - fill_adjacency_matrix_from_st(); } public: // Set graph from OFF file. @@ -321,13 +293,8 @@ class Cover_complex { * */ void set_graph_from_OFF() { - int num_edges = one_skeleton.size(); - if (num_edges > 0) { - for (int i = 0; i < num_edges; i++) st.insert_simplex_and_subfaces(one_skeleton[i]); - fill_adjacency_matrix_from_st(); - } else { - std::cout << "No triangulation read in OFF file!" << std::endl; - } + if(num_edges(one_skeleton_OFF)) one_skeleton = one_skeleton_OFF; + else cout << "No triangulation read in OFF file!" << endl; } public: // Set graph from Rips complex. @@ -339,9 +306,23 @@ class Cover_complex { */ template void set_graph_from_rips(double threshold, Distance distance) { - Rips_complex rips_complex_from_points(point_cloud, threshold, distance); - rips_complex_from_points.create_complex(st, 1); - fill_adjacency_matrix_from_st(); + if(distances.size() == 0) compute_pairwise_distances(distance); + for(int i = 0; i < n; i++){ + for(int j = i+1; j < n; j++){ + if(distances[i][j] <= threshold){ + add_edge(vertices[i], vertices[j], one_skeleton); + put(edge_weight, one_skeleton, edge(vertices[i], vertices[j], one_skeleton).first, distances[i][j]); + } + } + } + } + + public: + void set_graph_weights(){ + IndexMap index = get(vertex_index, one_skeleton); WeightMap weight = get(edge_weight, one_skeleton); + graph_traits::edge_iterator ei, ei_end; + for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) + put(weight, *ei, distances[index[source(*ei, one_skeleton)]][index[target(*ei, one_skeleton)]]); } public: // Pairwise distances. @@ -349,39 +330,34 @@ class Cover_complex { */ template void compute_pairwise_distances(Distance ref_distance) { - double d; - std::vector zeros(n); - for (int i = 0; i < n; i++) distances.push_back(zeros); - std::string distance = point_cloud_name; + double d; vector zeros(n); for (int i = 0; i < n; i++) distances.push_back(zeros); + string distance = point_cloud_name; distance.append("_dist"); - std::ifstream input(distance.c_str(), std::ios::out | std::ios::binary); + ifstream input(distance.c_str(), ios::out | ios::binary); if (input.good()) { if (verbose) std::cout << "Reading distances..." << std::endl; for (int i = 0; i < n; i++) { for (int j = i; j < n; j++) { input.read((char*)&d, 8); - distances[i][j] = d; - distances[j][i] = d; + distances[i][j] = d; distances[j][i] = d; } } input.close(); } else { - if (verbose) std::cout << "Computing distances..." << std::endl; - input.close(); - std::ofstream output(distance, std::ios::out | std::ios::binary); + if (verbose) cout << "Computing distances..." << endl; + input.close(); ofstream output(distance, ios::out | 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) cout << "\r" << state << "%" << flush; for (int j = i; j < n; j++) { double dis = ref_distance(point_cloud[i], point_cloud[j]); - distances[i][j] = dis; - distances[j][i] = dis; + distances[i][j] = dis; distances[j][i] = dis; output.write((char*)&dis, 8); } } output.close(); - if (verbose) std::cout << std::endl; + if (verbose) cout << endl; } } @@ -397,13 +373,13 @@ class Cover_complex { */ template 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); + int m = floor(n / exp((1 + rate_power) * log(log(n) / log(rate_constant)))); + m = min(m, n - 1); + vector samples(m); 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) cout << n << " points in R^" << data_dimension << endl; + if (verbose) cout << "Subsampling " << m << " points" << endl; if (distances.size() == 0) compute_pairwise_distances(distance); @@ -413,17 +389,14 @@ class Cover_complex { 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); + for (int k = 1; k < m; k++) mj = min(mj, distances[j][samples[k]]); + hausdorff_dist = max(hausdorff_dist, mj); } delta += hausdorff_dist / N; } - if (verbose) std::cout << "delta = " << delta << std::endl; - Rips_complex rips_complex_from_points(point_cloud, delta, distance); - rips_complex_from_points.create_complex(st, 1); - fill_adjacency_matrix_from_st(); - + if (verbose) cout << "delta = " << delta << endl; + set_graph_from_rips(delta, distance); return delta; } @@ -438,15 +411,9 @@ class Cover_complex { * */ void set_function_from_file(const std::string& func_file_name) { - int vertex_id = 0; - std::ifstream input(func_file_name); - std::string line; - double f; - while (std::getline(input, line)) { - std::stringstream stream(line); - stream >> f; - func.emplace(vertex_id, f); - vertex_id++; + int i = 0; ifstream input(func_file_name); string line; double f; + while (getline(input, line)) { + stringstream stream(line); stream >> f; func.emplace(i, f); i++; } functional_cover = true; cover_name = func_file_name; @@ -473,13 +440,9 @@ class Cover_complex { * */ template - void set_function_from_range(InputRange const& function) { + void set_function_from_range(InputRange const& f) { + for (int i = 0; i < n; i++) func.emplace(i, f[i]); functional_cover = true; - int index = 0; - for (auto v : function) { - func.emplace(index, v); - index++; - } } // ******************************************************************************************************************* @@ -504,29 +467,21 @@ class Cover_complex { return 0; } - double reso = 0; + double reso = 0; IndexMap index = get(vertex_index, one_skeleton); if (type == "GIC") { - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector vertices; - for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); - reso = std::max(reso, std::abs(func[vertices[0]] - func[vertices[1]])); - } - } - if (verbose) std::cout << "resolution = " << reso << std::endl; + graph_traits::edge_iterator ei, ei_end; + for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) + reso = max(reso, abs(func[index[source(*ei, one_skeleton)]] - func[index[target(*ei, one_skeleton)]])); + if (verbose) cout << "resolution = " << reso << endl; resolution_double = reso; } if (type == "Nerve") { - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector vertices; - for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); - reso = std::max(reso, (std::abs(func[vertices[0]] - func[vertices[1]])) / gain); - } - } - if (verbose) std::cout << "resolution = " << reso << std::endl; + graph_traits::edge_iterator ei, ei_end; + for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) + reso = max(reso, abs(func[index[source(*ei, one_skeleton)]] - func[index[target(*ei, one_skeleton)]]) / gain); + if (verbose) cout << "resolution = " << reso << endl; resolution_double = reso; } @@ -559,77 +514,69 @@ 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; + cout << "Number and/or length of intervals not specified" << endl; return; } if (gain == -1) { - std::cout << "Gain not specified" << std::endl; + cout << "Gain not specified" << endl; return; } // Read function values and compute min and max - std::map::iterator it; - double maxf, minf; - minf = std::numeric_limits::max(); - maxf = std::numeric_limits::min(); - for (it = func.begin(); it != func.end(); it++) { - minf = std::min(minf, it->second); - maxf = std::max(maxf, it->second); + double minf = numeric_limits::max(); double maxf = numeric_limits::lowest(); + for (int i = 0; i < n; i++) { + minf = min(minf, func[i]); maxf = max(maxf, func[i]); } - int n = func.size(); - if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; + if (verbose) cout << "Min function value = " << minf << " and Max function value = " << maxf << endl; // Compute cover of im(f) - std::vector > intervals; - int res; + vector > intervals; int res; if (resolution_double == -1) { // Case we use an integer for the number of intervals. double incr = (maxf - minf) / resolution_int; double x = minf; double alpha = (incr * gain) / (2 - 2 * gain); double y = minf + incr + alpha; - std::pair interm(x, y); + pair interm(x, y); intervals.push_back(interm); for (int i = 1; i < resolution_int - 1; i++) { x = minf + i * incr - alpha; y = minf + (i + 1) * incr + alpha; - std::pair inter(x, y); + pair inter(x, y); intervals.push_back(inter); } x = minf + (resolution_int - 1) * incr - alpha; y = maxf; - std::pair interM(x, y); + pair interM(x, y); intervals.push_back(interM); res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" - << std::endl; + cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << endl; } } else { if (resolution_int == -1) { // Case we use a double for the length of the intervals. double x = minf; double y = x + resolution_double; while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { - std::pair inter(x, y); + pair inter(x, y); intervals.push_back(inter); x = y - gain * resolution_double; y = x + resolution_double; } - std::pair interM(x, maxf); + pair interM(x, maxf); intervals.push_back(interM); res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" - << std::endl; + cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << endl; } } else { // Case we use an integer and a double for the length of the intervals. double x = minf; double y = x + resolution_double; int count = 0; while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { - std::pair inter(x, y); + pair inter(x, y); intervals.push_back(inter); count++; x = y - gain * resolution_double; @@ -638,88 +585,80 @@ 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::endl; + cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << endl; } } } // 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)); - int id = 0; - int pos = 0; + vector points(n); for (int i = 0; i < n; i++) points[i] = i; + sort(points.begin(), points.end(), Less(this->func)); + + int id = 0; int pos = 0; int maxc = -1; IndexMap index = get(vertex_index, one_skeleton); for (int i = 0; i < res; i++) { + // Find points in the preimage - std::map > prop; - std::pair inter1 = intervals[i]; - int tmp = pos; + vector indices; pair inter1 = intervals[i]; + int tmp = pos; double u, v; Graph G = one_skeleton.create_subgraph(); if (i != res - 1) { + if (i != 0) { - std::pair inter3 = intervals[i - 1]; - while (func[points[tmp]] < inter3.second && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; + pair inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n){ + add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } + u = inter3.second; } + else u = inter1.first; - std::pair inter2 = intervals[i + 1]; - while (func[points[tmp]] < inter2.first && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; + pair inter2 = intervals[i + 1]; + while (func[points[tmp]] < inter2.first && tmp != n){ + add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } + v = inter2.first; + pos = tmp; - while (func[points[tmp]] < inter1.second && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; + while (func[points[tmp]] < inter1.second && tmp != n){ + add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } } else { - std::pair inter3 = intervals[i - 1]; - while (func[points[tmp]] < inter3.second && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; + pair inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n){ + add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } - while (tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; - tmp++; + while (tmp != n){ + add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } + + u = inter3.second; v = inter1.second; + } - // Compute the connected components with DFS - std::map visit; - if (verbose) std::cout << "Preimage of interval " << i << std::endl; - for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) - visit[it->first] = false; - if (!(prop.empty())) { - for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) { - if (!(visit[it->first])) { - std::vector cc; - cc.clear(); - dfs(prop, it->first, cc, visit); - int cci = cc.size(); - if (verbose) std::cout << "one CC with " << cci << " points, "; - double average_col = 0; - for (int j = 0; j < cci; j++) { - cover[cc[j]].push_back(id); - cover_back[id].push_back(cc[j]); - average_col += func_color[cc[j]] / cci; - } - cover_fct[id] = i; - cover_color[id] = std::pair(cci, average_col); - id++; - } - } + int num = num_vertices(G); vector component(num); + + // Compute connected components + connected_components(G, &component[0]); int maxct = maxc + 1; + + // Update covers + for(int j = 0; j < num; j++){ + maxc = max(maxc, maxct + component[j]); + cover [indices[j]] .push_back(maxct + component[j]); + cover_back [maxct + component[j]] .push_back(indices[j]); + cover_fct [maxct + component[j]] = i; + cover_std [maxct + component[j]] = 0.5*(u+v); + cover_color [maxct + component[j]] .second += func_color[indices[j]]; //= pair(cci, average_col); + cover_color [maxct + component[j]] .first += 1; } - if (verbose) std::cout << std::endl; } maximal_dim = id - 1; + for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) + iit->second.second /= iit->second.first; } public: // Set cover from file. @@ -729,30 +668,27 @@ class Cover_complex { * @param[in] cover_file_name name of the input cover file. * */ - void set_cover_from_file(const std::string& cover_file_name) { - int vertex_id = 0; - Cover_t cov; - std::vector cov_elts, cov_number; - std::ifstream input(cover_file_name); - std::string line; - while (std::getline(input, line)) { + void set_cover_from_file(const string & cover_file_name) { + int i = 0; int cov; vector cov_elts, cov_number; + ifstream input(cover_file_name); string line; + while (getline(input, line)) { cov_elts.clear(); - std::stringstream stream(line); + stringstream stream(line); while (stream >> cov) { cov_elts.push_back(cov); cov_number.push_back(cov); - cover_fct[cov] = cov; - cover_color[cov].second += func_color[vertex_id]; - cover_color[cov].first++; - cover_back[cov].push_back(vertex_id); + cover_fct [cov] = cov; + cover_color [cov] .second += func_color[i]; + cover_color [cov] .first++; + cover_back [cov] .push_back(i); } - cover[vertex_id] = cov_elts; - vertex_id++; + cover[i] = cov_elts; i++; } - std::vector::iterator it; - std::sort(cov_number.begin(), cov_number.end()); - it = std::unique(cov_number.begin(), cov_number.end()); - cov_number.resize(std::distance(cov_number.begin(), it)); + + sort(cov_number.begin(), cov_number.end()); + vector::iterator it = unique(cov_number.begin(), cov_number.end()); + cov_number.resize(distance(cov_number.begin(), it)); + maximal_dim = cov_number.size() - 1; for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first; cover_name = cover_file_name; @@ -766,60 +702,24 @@ class Cover_complex { * */ template - void set_cover_from_Voronoi(Distance distance, int m = 100) { - voronoi_subsamples.resize(m); - SampleWithoutReplacement(n, m, voronoi_subsamples); - if (distances.size() == 0) compute_pairwise_distances(distance); - std::vector mindist(n); - for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits::max(); + void set_cover_from_Voronoi(Distance distance, int m = 100){ + + voronoi_subsamples.resize(m); SampleWithoutReplacement(n, m, voronoi_subsamples); + if (distances.size() == 0) compute_pairwise_distances(distance); set_graph_weights(); + WeightMap weight = get(edge_weight, one_skeleton); IndexMap index = get(vertex_index, one_skeleton); + vector mindist(n); for (int j = 0; j < n; j++) mindist[j] = numeric_limits::max(); // Compute the geodesic distances to subsamples with Dijkstra 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 dist(n); - std::vector process(n); - for (int j = 0; j < n; j++) { - dist[j] = std::numeric_limits::max(); - process[j] = j; - } - dist[seed] = 0; - int curr_size = process.size(); - int min_point, min_index; - double min_dist; - std::vector neighbors; - int num_neighbors; - - while (curr_size > 0) { - min_dist = std::numeric_limits::max(); - min_index = -1; - min_point = -1; - for (int j = 0; j < curr_size; j++) { - if (dist[process[j]] < min_dist) { - min_point = process[j]; - min_dist = dist[process[j]]; - min_index = j; - } - } - assert(min_index != -1); - process.erase(process.begin() + min_index); - assert(min_point != -1); - neighbors = adjacency_matrix[min_point]; - num_neighbors = neighbors.size(); - for (int j = 0; j < num_neighbors; j++) { - double d = dist[min_point] + distances[min_point][neighbors[j]]; - dist[neighbors[j]] = std::min(dist[neighbors[j]], d); - } - curr_size = process.size(); - } + + if (verbose) cout << "Computing geodesic distances to seed " << i << "..." << endl; + int seed = voronoi_subsamples[i]; vector dmap(n); + dijkstra_shortest_paths(one_skeleton, vertices[seed], weight_map(weight).distance_map(make_iterator_property_map(dmap.begin(), index))); for (int j = 0; j < n; j++) - if (mindist[j] > dist[j]) { - mindist[j] = dist[j]; - if (cover[j].size() == 0) - cover[j].push_back(i); - else - cover[j][0] = i; + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; + if (cover[j].size() == 0) cover[j].push_back(i); else cover[j][0] = i; } } @@ -831,6 +731,7 @@ class Cover_complex { for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first; maximal_dim = m - 1; cover_name = "Voronoi"; + } public: // return subset of data corresponding to a node @@ -840,7 +741,7 @@ class Cover_complex { * @result cover_back(c) vector of IDs of data points. * */ - const std::vector& subpopulation(Cover_t c) { return cover_back[c]; } + const vector & subpopulation(int c) { return cover_back[c]; } // ******************************************************************************************************************* // Visualization. @@ -853,16 +754,17 @@ class Cover_complex { * @param[in] color_file_name name of the input color file. * */ - void set_color_from_file(const std::string& color_file_name) { - int vertex_id = 0; - std::ifstream input(color_file_name); - std::string line; + void set_color_from_file(const string & color_file_name) { + int i = 0; + ifstream input(color_file_name); + string line; double f; - while (std::getline(input, line)) { - std::stringstream stream(line); + while (getline(input, line)) { + stringstream stream(line); + //stream >> one_skeleton[vertices[i]].color; stream >> f; - func_color.emplace(vertex_id, f); - vertex_id++; + func_color.emplace(i, f); + i++; } color_name = color_file_name; } @@ -874,9 +776,9 @@ class Cover_complex { * */ void set_color_from_coordinate(int k = 0) { - for (int i = 0; i < n; i++) func_color.emplace(i, point_cloud[i][k]); + for (int i = 0; i < n; i++) func_color[i] = point_cloud[i][k]; color_name = "coordinate "; - color_name.append(std::to_string(k)); + color_name.append(to_string(k)); } public: // Set color from vector. @@ -885,8 +787,8 @@ class Cover_complex { * @param[in] color input vector of values. * */ - void set_color_from_vector(std::vector color) { - for (unsigned int i = 0; i < color.size(); i++) func_color.emplace(i, color[i]); + void set_color_from_vector(vector c) { + for (unsigned int i = 0; i < c.size(); i++) func_color[i] = c[i]; } public: // Create a .dot file that can be compiled with neato to produce a .pdf file. @@ -895,22 +797,18 @@ class Cover_complex { * of its 1-skeleton in a .pdf file. */ void plot_DOT() { - char mapp[11] = "SC.dot"; - std::ofstream graphic(mapp); - graphic << "graph GIC {" << std::endl; - double maxv, minv; - maxv = std::numeric_limits::min(); - minv = std::numeric_limits::max(); - for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); - iit++) { - maxv = std::max(maxv, iit->second.second); - minv = std::min(minv, iit->second.second); + + char mapp[100]; sprintf(mapp, "%s_sc.dot",point_cloud_name.c_str()); ofstream graphic(mapp); + + double maxv = numeric_limits::lowest(); double minv = numeric_limits::max(); + for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { + maxv = max(maxv, iit->second.second); minv = min(minv, iit->second.second); } - int k = 0; - std::vector nodes; - nodes.clear(); - for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); - iit++) { + + int k = 0; vector nodes; nodes.clear(); + + graphic << "graph GIC {" << endl; + for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { if (iit->second.first > mask) { nodes.push_back(iit->first); graphic << iit->first << "[shape=circle fontcolor=black color=black label=\"" << iit->first << ":" @@ -930,7 +828,7 @@ class Cover_complex { } graphic << "}"; graphic.close(); - std::cout << "SC.dot generated. It can be visualized with e.g. neato." << std::endl; + cout << ".dot file generated. It can be visualized with e.g. neato." << endl; } public: // Create a .txt file that can be compiled with KeplerMapper. @@ -938,22 +836,21 @@ class Cover_complex { * KeplerMapper. */ void write_info() { - int num_simplices = simplices.size(); - int num_edges = 0; - char mapp[11] = "SC.txt"; - std::ofstream graphic(mapp); + + int num_simplices = simplices.size(); int num_edges = 0; + char mapp[100]; sprintf(mapp, "%s_sc.txt",point_cloud_name.c_str()); ofstream graphic(mapp); + for (int i = 0; i < num_simplices; i++) if (simplices[i].size() == 2) if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++; - graphic << point_cloud_name << std::endl; - graphic << cover_name << std::endl; - graphic << color_name << std::endl; - graphic << resolution_double << " " << gain << std::endl; - graphic << cover_color.size() << " " << num_edges << std::endl; + graphic << point_cloud_name << endl; + graphic << cover_name << endl; + graphic << color_name << endl; + graphic << resolution_double << " " << gain << endl; + graphic << cover_color.size() << " " << num_edges << endl; - for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); - iit++) + for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) graphic << iit->first << " " << iit->second.second << " " << iit->second.first << std::endl; for (int i = 0; i < num_simplices; i++) @@ -961,8 +858,8 @@ class Cover_complex { if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) graphic << simplices[i][0] << " " << simplices[i][1] << std::endl; graphic.close(); - std::cout << "SC.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." - << std::endl; + cout << ".txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << endl; + } public: // Create a .off file that can be visualized (e.g. with Geomview). @@ -971,15 +868,15 @@ class Cover_complex { * the remaining coordinates of the points embedded in 3D are set to 0. */ void plot_OFF() { + assert(cover_name == "Voronoi"); - char gic[11] = "SC.off"; - std::ofstream graphic(gic); - graphic << "OFF" << std::endl; - int m = voronoi_subsamples.size(); - int numedges = 0; - int numfaces = 0; - std::vector > edges, faces; + + int m = voronoi_subsamples.size(); int numedges = 0; int numfaces = 0; vector > edges, faces; int numsimplices = simplices.size(); + + char gic[100]; sprintf(gic, "%s_sc.off",point_cloud_name.c_str()); ofstream graphic(gic); + + graphic << "OFF" << std::endl; for (int i = 0; i < numsimplices; i++) { if (simplices[i].size() == 2) { numedges++; @@ -1004,7 +901,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 << "SC.off generated. It can be visualized with e.g. geomview." << std::endl; + cout << ".off generated. It can be visualized with e.g. geomview." << endl; } // ******************************************************************************************************************* @@ -1013,121 +910,180 @@ class Cover_complex { public: /** \brief Creates the simplicial complex. * - * @param[in] complex SimplicialComplexForRips to be created. + * @param[in] complex SimplicialComplex to be created. * */ - template - void create_complex(SimplicialComplexForRips& complex) { + template + void create_complex(SimplicialComplex& complex){ unsigned int dimension = 0; for (auto const& simplex : simplices) { - complex.insert_simplex_and_subfaces(simplex); + int numvert = simplex.size(); double filt = std::numeric_limits::lowest(); + for(int i = 0; i < numvert; i++) filt = std::max(cover_color[simplex[i]].second, filt); + complex.insert_simplex_and_subfaces(simplex, filt); if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; } complex.set_dimension(dimension); } + public: + /** \brief Computes the extended persistence diagram of the complex. + * + */ + template + void compute_PD(){ + + SimplicialComplex streef, streeb; unsigned int dimension = 0; + for (auto const & simplex : simplices) { + int numvert = simplex.size(); double filtM = numeric_limits::lowest(); double filtm = filtM; + for(int i = 0; i < numvert; i++){filtM = max(cover_std[simplex[i]], filtM); filtm = max(-cover_std[simplex[i]], filtm);} + streef.insert_simplex_and_subfaces(simplex, filtM); streeb.insert_simplex_and_subfaces(simplex, filtm); + if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; + } streef.set_dimension(dimension); streeb.set_dimension(dimension); + + streef.initialize_filtration(); + Gudhi::persistent_cohomology::Persistent_cohomology pcohf(streef); + pcohf.init_coefficients(2); pcohf.compute_persistent_cohomology(); + pcohf.output_diagram(); + + streeb.initialize_filtration(); + Gudhi::persistent_cohomology::Persistent_cohomology pcohb(streeb); + pcohb.init_coefficients(2); pcohb.compute_persistent_cohomology(); + pcohb.output_diagram(); + + //PD = pcohf.get_persistent_pairs(); + + } + + public: + /** \brief Computes bootstrapped distances distribution. + * + * @param[in] N number of bootstrap iterations. + * + */ + template + void compute_distribution(int N = 100){ + + if(distribution.size() >= N) std::cout << "Already done!" << std::endl; + else{ + for(int i = 0; i < N-distribution.size(); i++){ + + Cover_complex Cboot; Cboot.n = this->n; 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]); + } + for(int j = 0; j < n; j++){ + vector dist(n); + for(int k = 0; k < n; k++) + dist[k] = distances[boot[j]][boot[k]]; + Cboot.distances.push_back(dist); + } + + Cboot.set_graph_from_automatic_rips(Gudhi::Euclidean_distance()); + Cboot.set_automatic_resolution(); Cboot.set_gain(); Cboot.set_cover_from_function(); + Cboot.find_simplices(); Cboot.compute_PD >(); + + distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); + + } + + std::sort(distribution.begin(), distribution.end()); + + } + } + + public: + /** \brief Computes the bottleneck distance corresponding to a specific confidence level. + * + * @param[in] alpha Confidence level. + * + */ + double compute_distance_from_confidence_level(double alpha){ + int N = distribution.size(); + return distribution[std::floor(alpha*N)]; + } + + public: + /** \brief Computes the confidence level of a specific bottleneck distance. + * + * @param[in] d Bottleneck distance. + * + */ + 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; + } + + public: + /** \brief Computes the p-value, i.e. the opposite of the confidence level of the largest bottleneck + * distance preserving the points in the persistence diagram of the output simplicial complex. + * + */ + double compute_p_value(){ + 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); + } + public: /** \brief Computes the simplices of the simplicial complex. */ void find_simplices() { if (type != "Nerve" && type != "GIC") { - std::cout << "Type of complex needs to be specified." << std::endl; + cout << "Type of complex needs to be specified." << endl; return; } if (type == "Nerve") { - for (std::map >::iterator it = cover.begin(); it != cover.end(); it++) + for (map >::iterator it = cover.begin(); it != cover.end(); it++) simplices.push_back(it->second); - std::vector >::iterator it; - std::sort(simplices.begin(), simplices.end()); - it = std::unique(simplices.begin(), simplices.end()); - simplices.resize(std::distance(simplices.begin(), it)); + sort(simplices.begin(), simplices.end()); + vector >::iterator it = unique(simplices.begin(), simplices.end()); + simplices.resize(distance(simplices.begin(), it)); } if (type == "GIC") { + + IndexMap index = get(vertex_index, one_skeleton); + if (functional_cover) { // Computes the simplices in the GIC by looking at all the edges of the graph and adding the // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals. if (gain >= 0.5) - throw std::invalid_argument( - "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); - - int v1, v2; - - // Loop on all points. - for (std::map >::iterator it = cover.begin(); it != cover.end(); it++) { - int vid = it->first; - std::vector neighbors = adjacency_matrix[vid]; - int num_neighb = neighbors.size(); - - // Find cover of current point (vid). - if (cover[vid].size() == 2) - v1 = std::min(cover[vid][0], cover[vid][1]); - else - v1 = cover[vid][0]; - std::vector node(1); - node[0] = v1; - simplices.push_back(node); - - // Loop on neighbors. - for (int i = 0; i < num_neighb; i++) { - int neighb = neighbors[i]; - - // Find cover of neighbor (neighb). - if (cover[neighb].size() == 2) - v2 = std::max(cover[neighb][0], cover[neighb][1]); - else - v2 = cover[neighb][0]; - - // If neighbor is in next interval, add edge. - if (cover_fct[v2] == cover_fct[v1] + 1) { - std::vector edge(2); - edge[0] = v1; - edge[1] = v2; - simplices.push_back(edge); + throw invalid_argument("the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); + + // Loop on all edges. + graph_traits::edge_iterator ei, ei_end; + for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei){ + int nums = cover[index[source(*ei, one_skeleton)]].size(); + for(int i = 0; i < nums; i++){ + int vs = cover[index[source(*ei, one_skeleton)]][i]; + int numt = cover[index[target(*ei, one_skeleton)]].size(); + for(int j = 0; j < numt; j++){ + int vt = cover[index[target(*ei, one_skeleton)]][j]; + if(cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1){ + vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); + } } } } - std::vector >::iterator it; - std::sort(simplices.begin(), simplices.end()); - it = std::unique(simplices.begin(), simplices.end()); - simplices.resize(std::distance(simplices.begin(), it)); + sort(simplices.begin(), simplices.end()); + vector >::iterator it = unique(simplices.begin(), simplices.end()); + simplices.resize(distance(simplices.begin(), it)); } else { - // Find IDs of edges to remove - std::vector simplex_to_remove; - int simplex_id = 0; - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector > comp; - for (auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); - if (comp[0].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id); - } - simplex_id++; - } - // Remove edges - if (simplex_to_remove.size() > 1) { - int current_id = 1; - auto simplex = st.complex_simplex_range().begin(); - int num_rem = 0; - for (int i = 0; i < simplex_id - 1; i++) { - int j = i + 1; - auto simplex_tmp = simplex; - simplex_tmp++; - if (j == simplex_to_remove[current_id]) { - st.remove_maximal_simplex(*simplex_tmp); - current_id++; - num_rem++; - } else { - simplex++; - } + // Find edges to keep + Simplex_tree st; graph_traits::edge_iterator ei, ei_end; + for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) + if( !( cover[index[target(*ei, one_skeleton)]].size() == 1 && + cover[index[target(*ei, one_skeleton)]] == cover[index[source(*ei, one_skeleton)]]) ){ + vector edge(2); edge[0] = index[source(*ei, one_skeleton)]; edge[1] = index[target(*ei, one_skeleton)]; + st.insert_simplex_and_subfaces(edge); } - simplex = st.complex_simplex_range().begin(); - for (int i = 0; i < simplex_to_remove[0]; i++) simplex++; - st.remove_maximal_simplex(*simplex); - } + + //st.insert_graph(one_skeleton); // Build the Simplex Tree corresponding to the graph st.expansion(maximal_dim); @@ -1136,24 +1092,23 @@ class Cover_complex { simplices.clear(); for (auto simplex : st.complex_simplex_range()) { if (!st.has_children(simplex)) { - std::vector simplx; + vector simplx; for (auto vertex : st.simplex_vertex_range(simplex)) { unsigned int sz = cover[vertex].size(); for (unsigned int i = 0; i < sz; i++) { simplx.push_back(cover[vertex][i]); } } - - std::sort(simplx.begin(), simplx.end()); - std::vector::iterator it = std::unique(simplx.begin(), simplx.end()); - simplx.resize(std::distance(simplx.begin(), it)); + sort(simplx.begin(), simplx.end()); + vector::iterator it = unique(simplx.begin(), simplx.end()); + simplx.resize(distance(simplx.begin(), it)); simplices.push_back(simplx); } } - std::vector >::iterator it; - std::sort(simplices.begin(), simplices.end()); - it = std::unique(simplices.begin(), simplices.end()); - simplices.resize(std::distance(simplices.begin(), it)); + sort(simplices.begin(), simplices.end()); + vector >::iterator it = unique(simplices.begin(), simplices.end()); + simplices.resize(distance(simplices.begin(), it)); + } } } @@ -1164,3 +1119,203 @@ class Cover_complex { } // namespace Gudhi #endif // GIC_H_ + + + + + + + + + + + + + + + + + +/*Old code. + + private: + void fill_adjacency_matrix_from_st() { + std::vector empty; + for (int i = 0; i < n; i++) adjacency_matrix[i] = empty; + for (auto simplex : st.complex_simplex_range()) { + if (st.dimension(simplex) == 1) { + std::vector vertices; + for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); + adjacency_matrix[vertices[0]].push_back(vertices[1]); + adjacency_matrix[vertices[1]].push_back(vertices[0]); + } + } + } + +std::vector simplex_to_remove; +int simplex_id = 0; +for (auto simplex : st.complex_simplex_range()) { + if (st.dimension(simplex) == 1) { + std::vector > comp; + for (auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); + if (comp[0].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id); + } + simplex_id++; +} + +// Remove edges +if (simplex_to_remove.size() > 1) { + int current_id = 1; + auto simplex = st.complex_simplex_range().begin(); + int num_rem = 0; + for (int i = 0; i < simplex_id - 1; i++) { + int j = i + 1; + auto simplex_tmp = simplex; + simplex_tmp++; + if (j == simplex_to_remove[current_id]) { + st.remove_maximal_simplex(*simplex_tmp); + current_id++; + num_rem++; + } else { + simplex++; + } + } + simplex = st.complex_simplex_range().begin(); + for (int i = 0; i < simplex_to_remove[0]; i++) simplex++; + st.remove_maximal_simplex(*simplex); +} + + +if(cover[index[source(*ei, one_skeleton)]].size() == 1){ + vs = cover[index[source(*ei, one_skeleton)]][0]; + vm = cover_fct[vs]; +} +else{ + vs0 = cover[index[source(*ei, one_skeleton)]][0]; + vs1 = cover[index[source(*ei, one_skeleton)]][1]; + vm = min(cover_fct[vs0], cover_fct[vs1]); + if(vm == cover_fct[vs0]) vs = vs0; else vs = vs1; +} + +if(cover[index[target(*ei, one_skeleton)]].size() == 1){ + vt = cover[index[target(*ei, one_skeleton)]][0]; + vM = cover_fct[vt]; +} +else{ + vt0 = cover[index[target(*ei, one_skeleton)]][0]; + vt1 = cover[index[target(*ei, one_skeleton)]][1]; + vM = max(cover_fct[vt0], cover_fct[vt1]); + if(vM == cover_fct[vt0]) vt = vt0; else vt = vt1; +} + +if(vM == vm + 1){ + //if(max(cover_fct[cover[index[target(*ei, one_skeleton)]][0]], cover_fct[cover[index[target(*ei, one_skeleton)]][1]])== min(cover_fct[index[source(*ei, one_skeleton)]][0], cover_fct[index[source(*ei, one_skeleton)]][1]) + 1){ + vector edge(2); edge[0] = vs; edge[1] = vt; + simplices.push_back(edge); +} + + +for (map >::iterator it = cover.begin(); it != cover.end(); it++) { +int vid = it->first; +vector neighbors = adjacency_matrix[vid]; +int num_neighb = neighbors.size(); + +// Find cover of current point (vid). +if (cover[vid].size() == 2) v1 = std::min(cover[vid][0], cover[vid][1]); +else v1 = cover[vid][0]; +vector node(1); node[0] = v1; +simplices.push_back(node); + +// Loop on neighbors. +for (int i = 0; i < num_neighb; i++) { + int neighb = neighbors[i]; + + // Find cover of neighbor (neighb). + if (cover[neighb].size() == 2) v2 = std::max(cover[neighb][0], cover[neighb][1]); + else v2 = cover[neighb][0]; + + // If neighbor is in next interval, add edge. + if (cover_fct[v2] == cover_fct[v1] + 1) { + vector edge(2); edge[0] = v1; edge[1] = v2; + simplices.push_back(edge); break; + } +} +} + + + std::vector dist(n); + std::vector process(n); + for (int j = 0; j < n; j++) { + dist[j] = std::numeric_limits::max(); + process[j] = j; + } + dist[seed] = 0; + int curr_size = process.size(); + int min_point, min_index; + double min_dist; + std::vector neighbors; + int num_neighbors; + + while (curr_size > 0) { + min_dist = std::numeric_limits::max(); + min_index = -1; + min_point = -1; + for (int j = 0; j < curr_size; j++) { + if (dist[process[j]] < min_dist) { + min_point = process[j]; + min_dist = dist[process[j]]; + min_index = j; + } + } + assert(min_index != -1); + process.erase(process.begin() + min_index); + assert(min_point != -1); + neighbors = adjacency_matrix[min_point]; + num_neighbors = neighbors.size(); + for (int j = 0; j < num_neighbors; j++) { + double d = dist[min_point] + distances[min_point][neighbors[j]]; + dist[neighbors[j]] = std::min(dist[neighbors[j]], d); + } + curr_size = process.size(); + } + + + // Compute the connected components with DFS + std::map visit; + if (verbose) std::cout << "Preimage of interval " << i << std::endl; + for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) + visit[it->first] = false; + if (!(prop.empty())) { + for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) { + if (!(visit[it->first])) { + std::vector cc; + cc.clear(); + dfs(prop, it->first, cc, visit); + int cci = cc.size(); + if (verbose) std::cout << "one CC with " << cci << " points, "; + double average_col = 0; + for (int j = 0; j < cci; j++) { + cover[cc[j]].push_back(id); + cover_back[id].push_back(cc[j]); + average_col += func_color[cc[j]] / cci; + } + cover_fct[id] = i; + cover_std[id] = std::pair(cci, 0.5*(u+v)); + cover_color[id] = std::pair(cci, average_col); + id++; + } + } + } + + + // DFS + private: + void dfs(std::map >& G, int p, std::vector& cc, std::map& visit) { + cc.push_back(p); + visit[p] = true; + int neighb = G[p].size(); + for (int i = 0; i < neighb; i++) + if (visit.find(G[p][i]) != visit.end()) + if (!(visit[G[p][i]])) dfs(G, G[p][i], cc, visit); + } +*/ -- cgit v1.2.3 From f679fa8ac13ec529ac4569bec2124b7ec1e7938b Mon Sep 17 00:00:00 2001 From: mcarrier Date: Tue, 19 Dec 2017 14:50:25 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3087 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 0d5e8163d4ea01c9a34ec6b6c6d4bd4fbedffb68 --- src/Nerve_GIC/include/gudhi/GIC.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index d8c6abd1..1a5f891b 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -1063,10 +1063,11 @@ class Cover_complex { for(int j = 0; j < numt; j++){ int vt = cover[index[target(*ei, one_skeleton)]][j]; if(cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1){ - vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); + vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); goto afterLoop; } } } + afterLoop: ; } sort(simplices.begin(), simplices.end()); vector >::iterator it = unique(simplices.begin(), simplices.end()); -- cgit v1.2.3 From 5d1909dfc9a81743ced07bc3dfad8fb28d3c29e2 Mon Sep 17 00:00:00 2001 From: mcarrier Date: Tue, 19 Dec 2017 22:41:09 +0000 Subject: removed namespace boost and std git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3091 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: a1c367d57f6d2874c4008e08be68a39d6f821706 --- src/Nerve_GIC/include/gudhi/GIC.h | 434 +++++++++++++++++++------------------- 1 file changed, 220 insertions(+), 214 deletions(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 1a5f891b..3cd8a92a 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -46,14 +46,11 @@ #include #include #include // for numeric_limits -#include // for pair<> +#include // for std::pair<> #include // for std::max #include #include -using namespace boost; -using namespace std; - namespace Gudhi { namespace cover_complex { @@ -61,11 +58,11 @@ namespace cover_complex { using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; using Rips_complex = Gudhi::rips_complex::Rips_complex; -using PersistenceDiagram = vector >; -using Graph = subgraph > > >; -using vertex_t = graph_traits::vertex_descriptor; -using IndexMap = property_map::type; -using WeightMap = property_map::type; +using PersistenceDiagram = std::vector >; +using Graph = boost::subgraph > > >; +using vertex_t = boost::graph_traits::vertex_descriptor; +using IndexMap = boost::property_map::type; +using WeightMap = boost::property_map::type; /** * \class Cover_complex @@ -93,32 +90,32 @@ class Cover_complex { bool verbose = false; // whether to display information. - vector point_cloud; + std::vector point_cloud; int maximal_dim; // maximal dimension of output simplicial complex. int data_dimension; // dimension of input data. int n; // number of points. - vector > distances; + std::vector > distances; - map func; // function used to compute the output simplicial complex. - map func_color; // function used to compute the colors of the nodes of the output simplicial complex. + 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. 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). Graph one_skeleton; // one-skeleton used to compute the connected components. - vector vertices; - vector > simplices; + std::vector vertices; + std::vector > simplices; - vector voronoi_subsamples; + std::vector voronoi_subsamples; PersistenceDiagram PD; - vector distribution; + std::vector distribution; - map > cover; - map > cover_back; - map cover_std; // standard function (induced by func) used to compute the extended persistence diagram of the output simplicial complex. - map cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. - map > cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex. + std::map > cover; + std::map > cover_back; + std::map cover_std; // standard function (induced by func) used to compute the extended persistence diagram of the output simplicial complex. + std::map cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. + std::map > cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex. int resolution_int = -1; double resolution_double = -1; @@ -127,15 +124,15 @@ class Cover_complex { double rate_power = 0.001; // Power in the subsampling. int mask = 0; // Ignore nodes containing less than mask points. - string cover_name; - string point_cloud_name; - string color_name; - string type; // Nerve or GIC + std::string cover_name; + std::string point_cloud_name; + std::string color_name; + std::string type; // Nerve or GIC // Point comparator struct Less { - Less(map func) { Fct = func; } - map Fct; + Less(std::map func) { Fct = func; } + std::map Fct; bool operator()(int a, int b) { if (Fct[a] == Fct[b]) return a < b; @@ -144,15 +141,21 @@ class Cover_complex { } }; + // 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); + } + // Find random number in [0,1]. double GetUniform() { - thread_local default_random_engine re; - thread_local uniform_real_distribution Dist(0, 1); + thread_local std::default_random_engine re; + thread_local std::uniform_real_distribution Dist(0, 1); return Dist(re); } // Subsample points. - void SampleWithoutReplacement(int populationSize, int sampleSize, vector & samples) { + void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector & samples) { int t = 0; int m = 0; double u; while (m < sampleSize){ u = GetUniform(); @@ -165,10 +168,10 @@ class Cover_complex { public: /** \brief Specifies whether the type of the output simplicial complex. * - * @param[in] t string (either "GIC" or "Nerve"). + * @param[in] t std::string (either "GIC" or "Nerve"). * */ - void set_type(const std::string& t) { type = t; } + void set_type(const std::string & t) { type = t; } public: /** \brief Specifies whether the program should display information or not. @@ -207,24 +210,24 @@ class Cover_complex { * @param[in] off_file_name name of the input .OFF or .nOFF file. * */ - bool read_point_cloud(const string & off_file_name) { + bool read_point_cloud(const std::string & off_file_name) { point_cloud_name = off_file_name; - ifstream input(off_file_name); - string line; + std::ifstream input(off_file_name); + std::string line; char comment = '#'; while (comment == '#') { - getline(input, line); + std::getline(input, line); if (!line.empty() && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } if (strcmp((char*)line.c_str(), "nOFF") == 0) { comment = '#'; while (comment == '#') { - getline(input, line); + std::getline(input, line); if (!line.empty() && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } - stringstream stream(line); + std::stringstream stream(line); stream >> data_dimension; } else { data_dimension = 3; @@ -233,34 +236,34 @@ class Cover_complex { comment = '#'; int numedges, numfaces, i, dim; while (comment == '#') { - getline(input, line); + std::getline(input, line); if (!line.empty() && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } - stringstream stream(line); + std::stringstream stream(line); stream >> n; stream >> numfaces; stream >> numedges; i = 0; while (i < n) { - getline(input, line); + std::getline(input, line); if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) { - istringstream iss(line); vector point; point.assign(istream_iterator(iss), istream_iterator()); + std::stringstream iss(line); std::vector point; point.assign(std::istream_iterator(iss), std::istream_iterator()); point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); - add_vertex(one_skeleton_OFF); vertices.push_back(add_vertex(one_skeleton)); + boost::add_vertex(one_skeleton_OFF); vertices.push_back(boost::add_vertex(one_skeleton)); i++; } } i = 0; while (i < numfaces) { - getline(input, line); + std::getline(input, line); if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && !all_of(line.begin(), line.end(), (int(*)(int))isspace)) { - vector simplex; istringstream iss(line); - simplex.assign(istream_iterator(iss), istream_iterator()); dim = simplex[0]; + std::vector simplex; std::stringstream iss(line); + simplex.assign(std::istream_iterator(iss), std::istream_iterator()); dim = simplex[0]; for (int j = 1; j <= dim; j++) for (int k = j + 1; k <= dim; k++) - add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF); + boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF); i++; } } @@ -280,11 +283,12 @@ class Cover_complex { * each edge being represented by the IDs of its two nodes. * */ - void set_graph_from_file(const string & graph_file_name){ - int neighb; ifstream input(graph_file_name); string line; int source; - while (getline(input, line)){ - stringstream stream(line); stream >> source; - while (stream >> neighb) add_edge(vertices[source], vertices[neighb], one_skeleton); + void set_graph_from_file(const std::string & graph_file_name){ + remove_edges(one_skeleton); + int neighb; std::ifstream input(graph_file_name); std::string line; int source; + while (std::getline(input, line)){ + std::stringstream stream(line); stream >> source; + while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton); } } @@ -293,8 +297,9 @@ class Cover_complex { * */ void set_graph_from_OFF() { + remove_edges(one_skeleton); if(num_edges(one_skeleton_OFF)) one_skeleton = one_skeleton_OFF; - else cout << "No triangulation read in OFF file!" << endl; + else std::cout << "No triangulation read in OFF file!" << std::endl; } public: // Set graph from Rips complex. @@ -306,12 +311,13 @@ class Cover_complex { */ template void set_graph_from_rips(double threshold, Distance distance) { + remove_edges(one_skeleton); if(distances.size() == 0) compute_pairwise_distances(distance); for(int i = 0; i < n; i++){ for(int j = i+1; j < n; j++){ if(distances[i][j] <= threshold){ - add_edge(vertices[i], vertices[j], one_skeleton); - put(edge_weight, one_skeleton, edge(vertices[i], vertices[j], one_skeleton).first, distances[i][j]); + boost::add_edge(vertices[i], vertices[j], one_skeleton); + boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first, distances[i][j]); } } } @@ -319,10 +325,10 @@ class Cover_complex { public: void set_graph_weights(){ - IndexMap index = get(vertex_index, one_skeleton); WeightMap weight = get(edge_weight, one_skeleton); - graph_traits::edge_iterator ei, ei_end; - for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) - put(weight, *ei, distances[index[source(*ei, one_skeleton)]][index[target(*ei, one_skeleton)]]); + IndexMap index = boost::get(boost::vertex_index, one_skeleton); WeightMap weight = boost::get(boost::edge_weight, one_skeleton); + boost::graph_traits::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) + boost::put(weight, *ei, distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]); } public: // Pairwise distances. @@ -330,10 +336,10 @@ class Cover_complex { */ template void compute_pairwise_distances(Distance ref_distance) { - double d; vector zeros(n); for (int i = 0; i < n; i++) distances.push_back(zeros); - string distance = point_cloud_name; + double d; std::vector zeros(n); for (int i = 0; i < n; i++) distances.push_back(zeros); + std::string distance = point_cloud_name; distance.append("_dist"); - ifstream input(distance.c_str(), ios::out | ios::binary); + std::ifstream input(distance.c_str(), std::ios::out | std::ios::binary); if (input.good()) { if (verbose) std::cout << "Reading distances..." << std::endl; @@ -345,11 +351,11 @@ class Cover_complex { } input.close(); } else { - if (verbose) cout << "Computing distances..." << endl; - input.close(); ofstream output(distance, ios::out | ios::binary); + if (verbose) std::cout << "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) cout << "\r" << state << "%" << flush; + if (state == 0 && verbose) std::cout << "\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; distances[j][i] = dis; @@ -357,7 +363,7 @@ class Cover_complex { } } output.close(); - if (verbose) cout << endl; + if (verbose) std::cout << std::endl; } } @@ -374,12 +380,12 @@ class Cover_complex { template double set_graph_from_automatic_rips(Distance distance, int N = 100) { int m = floor(n / exp((1 + rate_power) * log(log(n) / log(rate_constant)))); - m = min(m, n - 1); - vector samples(m); + m = std::min(m, n - 1); + std::vector samples(m); double delta = 0; - if (verbose) cout << n << " points in R^" << data_dimension << endl; - if (verbose) cout << "Subsampling " << m << " points" << endl; + if (verbose) std::cout << n << " points in R^" << data_dimension << std::endl; + if (verbose) std::cout << "Subsampling " << m << " points" << std::endl; if (distances.size() == 0) compute_pairwise_distances(distance); @@ -389,13 +395,13 @@ class Cover_complex { 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 = min(mj, distances[j][samples[k]]); - hausdorff_dist = max(hausdorff_dist, mj); + 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; } - if (verbose) cout << "delta = " << delta << endl; + if (verbose) std::cout << "delta = " << delta << std::endl; set_graph_from_rips(delta, distance); return delta; } @@ -411,9 +417,9 @@ class Cover_complex { * */ void set_function_from_file(const std::string& func_file_name) { - int i = 0; ifstream input(func_file_name); string line; double f; - while (getline(input, line)) { - stringstream stream(line); stream >> f; func.emplace(i, f); i++; + int i = 0; std::ifstream input(func_file_name); std::string line; double f; + while (std::getline(input, line)) { + std::stringstream stream(line); stream >> f; func.emplace(i, f); i++; } functional_cover = true; cover_name = func_file_name; @@ -467,21 +473,21 @@ class Cover_complex { return 0; } - double reso = 0; IndexMap index = get(vertex_index, one_skeleton); + double reso = 0; IndexMap index = boost::get(boost::vertex_index, one_skeleton); if (type == "GIC") { - graph_traits::edge_iterator ei, ei_end; - for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) - reso = max(reso, abs(func[index[source(*ei, one_skeleton)]] - func[index[target(*ei, one_skeleton)]])); - if (verbose) cout << "resolution = " << reso << endl; + boost::graph_traits::edge_iterator ei, ei_end; + 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; resolution_double = reso; } if (type == "Nerve") { - graph_traits::edge_iterator ei, ei_end; - for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) - reso = max(reso, abs(func[index[source(*ei, one_skeleton)]] - func[index[target(*ei, one_skeleton)]]) / gain); - if (verbose) cout << "resolution = " << reso << endl; + boost::graph_traits::edge_iterator ei, ei_end; + 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)]]) / gain); + if (verbose) std::cout << "resolution = " << reso << std::endl; resolution_double = reso; } @@ -514,69 +520,69 @@ class Cover_complex { */ void set_cover_from_function() { if (resolution_double == -1 && resolution_int == -1) { - cout << "Number and/or length of intervals not specified" << endl; + std::cout << "Number and/or length of intervals not specified" << std::endl; return; } if (gain == -1) { - cout << "Gain not specified" << endl; + std::cout << "Gain not specified" << std::endl; return; } // Read function values and compute min and max - double minf = numeric_limits::max(); double maxf = numeric_limits::lowest(); + double minf = std::numeric_limits::max(); double maxf = std::numeric_limits::lowest(); for (int i = 0; i < n; i++) { - minf = min(minf, func[i]); maxf = max(maxf, func[i]); + minf = std::min(minf, func[i]); maxf = std::max(maxf, func[i]); } - if (verbose) cout << "Min function value = " << minf << " and Max function value = " << maxf << endl; + if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; // Compute cover of im(f) - vector > intervals; int res; + std::vector > intervals; int res; if (resolution_double == -1) { // Case we use an integer for the number of intervals. double incr = (maxf - minf) / resolution_int; double x = minf; double alpha = (incr * gain) / (2 - 2 * gain); double y = minf + incr + alpha; - pair interm(x, y); + std::pair interm(x, y); intervals.push_back(interm); for (int i = 1; i < resolution_int - 1; i++) { x = minf + i * incr - alpha; y = minf + (i + 1) * incr + alpha; - pair inter(x, y); + std::pair inter(x, y); intervals.push_back(inter); } x = minf + (resolution_int - 1) * incr - alpha; y = maxf; - pair interM(x, y); + std::pair interM(x, y); intervals.push_back(interM); res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << endl; + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; } } else { if (resolution_int == -1) { // Case we use a double for the length of the intervals. double x = minf; double y = x + resolution_double; while (y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { - pair inter(x, y); + std::pair inter(x, y); intervals.push_back(inter); x = y - gain * resolution_double; y = x + resolution_double; } - pair interM(x, maxf); + std::pair interM(x, maxf); intervals.push_back(interM); res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << endl; + std::cout << "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. double x = minf; double y = x + resolution_double; int count = 0; while (count < resolution_int && y <= maxf && maxf - (y - gain * resolution_double) >= resolution_double) { - pair inter(x, y); + std::pair inter(x, y); intervals.push_back(inter); count++; x = y - gain * resolution_double; @@ -585,79 +591,79 @@ class Cover_complex { res = intervals.size(); if (verbose) { for (int i = 0; i < res; i++) - cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << endl; + std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; } } } // Sort points according to function values - vector points(n); for (int i = 0; i < n; i++) points[i] = i; - sort(points.begin(), points.end(), Less(this->func)); + std::vector points(n); for (int i = 0; i < n; i++) points[i] = i; + std::sort(points.begin(), points.end(), Less(this->func)); - int id = 0; int pos = 0; int maxc = -1; IndexMap index = get(vertex_index, one_skeleton); + int id = 0; int pos = 0; int maxc = -1; IndexMap index = boost::get(boost::vertex_index, one_skeleton); for (int i = 0; i < res; i++) { // Find points in the preimage - vector indices; pair inter1 = intervals[i]; + std::vector indices; std::pair inter1 = intervals[i]; int tmp = pos; double u, v; Graph G = one_skeleton.create_subgraph(); if (i != res - 1) { if (i != 0) { - pair inter3 = intervals[i - 1]; + std::pair inter3 = intervals[i - 1]; while (func[points[tmp]] < inter3.second && tmp != n){ - add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; + boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } u = inter3.second; } else u = inter1.first; - pair inter2 = intervals[i + 1]; + std::pair inter2 = intervals[i + 1]; while (func[points[tmp]] < inter2.first && tmp != n){ - add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; + boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } v = inter2.first; pos = tmp; while (func[points[tmp]] < inter1.second && tmp != n){ - add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; + boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } } else { - pair inter3 = intervals[i - 1]; + std::pair inter3 = intervals[i - 1]; while (func[points[tmp]] < inter3.second && tmp != n){ - add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; + boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } while (tmp != n){ - add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; + boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; } u = inter3.second; v = inter1.second; } - int num = num_vertices(G); vector component(num); + int num = num_vertices(G); std::vector component(num); // Compute connected components - connected_components(G, &component[0]); int maxct = maxc + 1; + boost::connected_components(G, &component[0]); int maxct = maxc + 1; // Update covers for(int j = 0; j < num; j++){ - maxc = max(maxc, maxct + component[j]); + maxc = std::max(maxc, maxct + component[j]); cover [indices[j]] .push_back(maxct + component[j]); cover_back [maxct + component[j]] .push_back(indices[j]); cover_fct [maxct + component[j]] = i; cover_std [maxct + component[j]] = 0.5*(u+v); - cover_color [maxct + component[j]] .second += func_color[indices[j]]; //= pair(cci, average_col); + cover_color [maxct + component[j]] .second += func_color[indices[j]]; cover_color [maxct + component[j]] .first += 1; } } maximal_dim = id - 1; - for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) + for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) iit->second.second /= iit->second.first; } @@ -668,12 +674,12 @@ class Cover_complex { * @param[in] cover_file_name name of the input cover file. * */ - void set_cover_from_file(const string & cover_file_name) { - int i = 0; int cov; vector cov_elts, cov_number; - ifstream input(cover_file_name); string line; - while (getline(input, line)) { + void set_cover_from_file(const std::string & cover_file_name) { + int i = 0; int cov; std::vector cov_elts, cov_number; + std::ifstream input(cover_file_name); std::string line; + while (std::getline(input, line)) { cov_elts.clear(); - stringstream stream(line); + std::stringstream stream(line); while (stream >> cov) { cov_elts.push_back(cov); cov_number.push_back(cov); @@ -685,9 +691,9 @@ class Cover_complex { cover[i] = cov_elts; i++; } - sort(cov_number.begin(), cov_number.end()); - vector::iterator it = unique(cov_number.begin(), cov_number.end()); - cov_number.resize(distance(cov_number.begin(), it)); + std::sort(cov_number.begin(), cov_number.end()); + std::vector::iterator it = std::unique(cov_number.begin(), cov_number.end()); + cov_number.resize(std::distance(cov_number.begin(), it)); maximal_dim = cov_number.size() - 1; for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first; @@ -706,15 +712,15 @@ class Cover_complex { voronoi_subsamples.resize(m); SampleWithoutReplacement(n, m, voronoi_subsamples); if (distances.size() == 0) compute_pairwise_distances(distance); set_graph_weights(); - WeightMap weight = get(edge_weight, one_skeleton); IndexMap index = get(vertex_index, one_skeleton); - vector mindist(n); for (int j = 0; j < n; j++) mindist[j] = numeric_limits::max(); + WeightMap weight = boost::get(boost::edge_weight, one_skeleton); IndexMap index = boost::get(boost::vertex_index, one_skeleton); + std::vector mindist(n); for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits::max(); // Compute the geodesic distances to subsamples with Dijkstra for (int i = 0; i < m; i++) { - if (verbose) cout << "Computing geodesic distances to seed " << i << "..." << endl; - int seed = voronoi_subsamples[i]; vector dmap(n); - dijkstra_shortest_paths(one_skeleton, vertices[seed], weight_map(weight).distance_map(make_iterator_property_map(dmap.begin(), index))); + 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]) { @@ -724,9 +730,9 @@ class Cover_complex { } for (int i = 0; i < n; i++) { - cover_back[cover[i][0]].push_back(i); - cover_color[cover[i][0]].second += func_color[i]; - cover_color[cover[i][0]].first++; + cover_back [cover[i][0]] .push_back(i); + cover_color [cover[i][0]] .second += func_color[i]; + cover_color [cover[i][0]] .first++; } for (int i = 0; i < m; i++) cover_color[i].second /= cover_color[i].first; maximal_dim = m - 1; @@ -741,7 +747,7 @@ class Cover_complex { * @result cover_back(c) vector of IDs of data points. * */ - const vector & subpopulation(int c) { return cover_back[c]; } + const std::vector & subpopulation(int c) { return cover_back[c]; } // ******************************************************************************************************************* // Visualization. @@ -754,13 +760,13 @@ class Cover_complex { * @param[in] color_file_name name of the input color file. * */ - void set_color_from_file(const string & color_file_name) { + void set_color_from_file(const std::string & color_file_name) { int i = 0; - ifstream input(color_file_name); - string line; + std::ifstream input(color_file_name); + std::string line; double f; - while (getline(input, line)) { - stringstream stream(line); + while (std::getline(input, line)) { + std::stringstream stream(line); //stream >> one_skeleton[vertices[i]].color; stream >> f; func_color.emplace(i, f); @@ -778,7 +784,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]; color_name = "coordinate "; - color_name.append(to_string(k)); + color_name.append(std::to_string(k)); } public: // Set color from vector. @@ -787,7 +793,7 @@ class Cover_complex { * @param[in] color input vector of values. * */ - void set_color_from_vector(vector c) { + void set_color_from_vector(std::vector c) { for (unsigned int i = 0; i < c.size(); i++) func_color[i] = c[i]; } @@ -798,17 +804,17 @@ class Cover_complex { */ void plot_DOT() { - char mapp[100]; sprintf(mapp, "%s_sc.dot",point_cloud_name.c_str()); ofstream graphic(mapp); + char mapp[100]; sprintf(mapp, "%s_sc.dot",point_cloud_name.c_str()); std::ofstream graphic(mapp); - double maxv = numeric_limits::lowest(); double minv = numeric_limits::max(); - for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { - maxv = max(maxv, iit->second.second); minv = min(minv, iit->second.second); + double maxv = std::numeric_limits::lowest(); double minv = std::numeric_limits::max(); + for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { + maxv = std::max(maxv, iit->second.second); minv = std::min(minv, iit->second.second); } - int k = 0; vector nodes; nodes.clear(); + int k = 0; std::vector nodes; nodes.clear(); - graphic << "graph GIC {" << endl; - for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { + graphic << "graph GIC {" << std::endl; + for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { if (iit->second.first > mask) { nodes.push_back(iit->first); graphic << iit->first << "[shape=circle fontcolor=black color=black label=\"" << iit->first << ":" @@ -828,7 +834,7 @@ class Cover_complex { } graphic << "}"; graphic.close(); - cout << ".dot file generated. It can be visualized with e.g. neato." << endl; + std::cout << ".dot file generated. It can be visualized with e.g. neato." << std::endl; } public: // Create a .txt file that can be compiled with KeplerMapper. @@ -838,19 +844,19 @@ class Cover_complex { void write_info() { int num_simplices = simplices.size(); int num_edges = 0; - char mapp[100]; sprintf(mapp, "%s_sc.txt",point_cloud_name.c_str()); ofstream graphic(mapp); + char mapp[100]; sprintf(mapp, "%s_sc.txt",point_cloud_name.c_str()); std::ofstream graphic(mapp); for (int i = 0; i < num_simplices; i++) if (simplices[i].size() == 2) if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) num_edges++; - graphic << point_cloud_name << endl; - graphic << cover_name << endl; - graphic << color_name << endl; - graphic << resolution_double << " " << gain << endl; - graphic << cover_color.size() << " " << num_edges << endl; + graphic << point_cloud_name << std::endl; + graphic << cover_name << std::endl; + graphic << color_name << std::endl; + graphic << resolution_double << " " << gain << std::endl; + graphic << cover_color.size() << " " << num_edges << std::endl; - for (map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) + for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) graphic << iit->first << " " << iit->second.second << " " << iit->second.first << std::endl; for (int i = 0; i < num_simplices; i++) @@ -858,7 +864,7 @@ class Cover_complex { if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) graphic << simplices[i][0] << " " << simplices[i][1] << std::endl; graphic.close(); - cout << ".txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << endl; + std::cout << ".txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; } @@ -871,10 +877,10 @@ class Cover_complex { assert(cover_name == "Voronoi"); - int m = voronoi_subsamples.size(); int numedges = 0; int numfaces = 0; vector > edges, faces; + int m = voronoi_subsamples.size(); int numedges = 0; int numfaces = 0; std::vector > edges, faces; int numsimplices = simplices.size(); - char gic[100]; sprintf(gic, "%s_sc.off",point_cloud_name.c_str()); ofstream graphic(gic); + char gic[100]; sprintf(gic, "%s_sc.off",point_cloud_name.c_str()); std::ofstream graphic(gic); graphic << "OFF" << std::endl; for (int i = 0; i < numsimplices; i++) { @@ -901,7 +907,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(); - cout << ".off generated. It can be visualized with e.g. geomview." << endl; + std::cout << ".off generated. It can be visualized with e.g. geomview." << std::endl; } // ******************************************************************************************************************* @@ -934,8 +940,8 @@ class Cover_complex { SimplicialComplex streef, streeb; unsigned int dimension = 0; for (auto const & simplex : simplices) { - int numvert = simplex.size(); double filtM = numeric_limits::lowest(); double filtm = filtM; - for(int i = 0; i < numvert; i++){filtM = max(cover_std[simplex[i]], filtM); filtm = max(-cover_std[simplex[i]], filtm);} + int numvert = simplex.size(); double filtM = std::numeric_limits::lowest(); double filtm = filtM; + for(int i = 0; i < numvert; i++){filtM = std::max(cover_std[simplex[i]], filtM); filtm = std::max(-cover_std[simplex[i]], filtm);} streef.insert_simplex_and_subfaces(simplex, filtM); streeb.insert_simplex_and_subfaces(simplex, filtm); if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; } streef.set_dimension(dimension); streeb.set_dimension(dimension); @@ -973,7 +979,7 @@ class Cover_complex { Cboot.point_cloud[j] = this->point_cloud[id]; Cboot.func.emplace(j,this->func[id]); } for(int j = 0; j < n; j++){ - vector dist(n); + std::vector dist(n); for(int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]]; Cboot.distances.push_back(dist); @@ -1030,57 +1036,57 @@ class Cover_complex { */ void find_simplices() { if (type != "Nerve" && type != "GIC") { - cout << "Type of complex needs to be specified." << endl; + std::cout << "Type of complex needs to be specified." << std::endl; return; } if (type == "Nerve") { - for (map >::iterator it = cover.begin(); it != cover.end(); it++) + for (std::map >::iterator it = cover.begin(); it != cover.end(); it++) simplices.push_back(it->second); - sort(simplices.begin(), simplices.end()); - vector >::iterator it = unique(simplices.begin(), simplices.end()); - simplices.resize(distance(simplices.begin(), it)); + std::sort(simplices.begin(), simplices.end()); + std::vector >::iterator it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); } if (type == "GIC") { - IndexMap index = get(vertex_index, one_skeleton); + IndexMap index = boost::get(boost::vertex_index, one_skeleton); if (functional_cover) { // Computes the simplices in the GIC by looking at all the edges of the graph and adding the // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals. if (gain >= 0.5) - throw invalid_argument("the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); + throw std::invalid_argument("the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5."); // Loop on all edges. - graph_traits::edge_iterator ei, ei_end; - for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei){ - int nums = cover[index[source(*ei, one_skeleton)]].size(); + boost::graph_traits::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei){ + int nums = cover[index[boost::source(*ei, one_skeleton)]].size(); for(int i = 0; i < nums; i++){ - int vs = cover[index[source(*ei, one_skeleton)]][i]; - int numt = cover[index[target(*ei, one_skeleton)]].size(); + int vs = cover[index[boost::source(*ei, one_skeleton)]][i]; + int numt = cover[index[boost::target(*ei, one_skeleton)]].size(); for(int j = 0; j < numt; j++){ - int vt = cover[index[target(*ei, one_skeleton)]][j]; + int vt = cover[index[boost::target(*ei, one_skeleton)]][j]; if(cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1){ - vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); goto afterLoop; + std::vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); goto afterLoop; } } } afterLoop: ; } - sort(simplices.begin(), simplices.end()); - vector >::iterator it = unique(simplices.begin(), simplices.end()); - simplices.resize(distance(simplices.begin(), it)); + std::sort(simplices.begin(), simplices.end()); + std::vector >::iterator it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); } else { // Find edges to keep - Simplex_tree st; graph_traits::edge_iterator ei, ei_end; - for (tie(ei, ei_end) = edges(one_skeleton); ei != ei_end; ++ei) - if( !( cover[index[target(*ei, one_skeleton)]].size() == 1 && - cover[index[target(*ei, one_skeleton)]] == cover[index[source(*ei, one_skeleton)]]) ){ - vector edge(2); edge[0] = index[source(*ei, one_skeleton)]; edge[1] = index[target(*ei, one_skeleton)]; + Simplex_tree st; boost::graph_traits::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) + if( !( cover[index[boost::target(*ei, one_skeleton)]].size() == 1 && + cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]]) ){ + std::vector edge(2); edge[0] = index[boost::source(*ei, one_skeleton)]; edge[1] = index[boost::target(*ei, one_skeleton)]; st.insert_simplex_and_subfaces(edge); } @@ -1093,22 +1099,22 @@ class Cover_complex { simplices.clear(); for (auto simplex : st.complex_simplex_range()) { if (!st.has_children(simplex)) { - vector simplx; + std::vector simplx; for (auto vertex : st.simplex_vertex_range(simplex)) { unsigned int sz = cover[vertex].size(); for (unsigned int i = 0; i < sz; i++) { simplx.push_back(cover[vertex][i]); } } - sort(simplx.begin(), simplx.end()); - vector::iterator it = unique(simplx.begin(), simplx.end()); - simplx.resize(distance(simplx.begin(), it)); + std::sort(simplx.begin(), simplx.end()); + std::vector::iterator it = std::unique(simplx.begin(), simplx.end()); + simplx.resize(std::distance(simplx.begin(), it)); simplices.push_back(simplx); } } - sort(simplices.begin(), simplices.end()); - vector >::iterator it = unique(simplices.begin(), simplices.end()); - simplices.resize(distance(simplices.begin(), it)); + std::sort(simplices.begin(), simplices.end()); + std::vector >::iterator it = std::unique(simplices.begin(), simplices.end()); + simplices.resize(std::distance(simplices.begin(), it)); } } @@ -1141,11 +1147,11 @@ class Cover_complex { private: void fill_adjacency_matrix_from_st() { - std::vector empty; + std::std::vector empty; for (int i = 0; i < n; i++) adjacency_matrix[i] = empty; for (auto simplex : st.complex_simplex_range()) { if (st.dimension(simplex) == 1) { - std::vector vertices; + std::std::vector vertices; for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex); adjacency_matrix[vertices[0]].push_back(vertices[1]); adjacency_matrix[vertices[1]].push_back(vertices[0]); @@ -1153,11 +1159,11 @@ class Cover_complex { } } -std::vector simplex_to_remove; +std::std::vector simplex_to_remove; int simplex_id = 0; for (auto simplex : st.complex_simplex_range()) { if (st.dimension(simplex) == 1) { - std::vector > comp; + std::std::vector > comp; for (auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]); if (comp[0].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id); } @@ -1211,20 +1217,20 @@ else{ if(vM == vm + 1){ //if(max(cover_fct[cover[index[target(*ei, one_skeleton)]][0]], cover_fct[cover[index[target(*ei, one_skeleton)]][1]])== min(cover_fct[index[source(*ei, one_skeleton)]][0], cover_fct[index[source(*ei, one_skeleton)]][1]) + 1){ - vector edge(2); edge[0] = vs; edge[1] = vt; + std::vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); } -for (map >::iterator it = cover.begin(); it != cover.end(); it++) { +for (std::map >::iterator it = cover.begin(); it != cover.end(); it++) { int vid = it->first; -vector neighbors = adjacency_matrix[vid]; +std::vector neighbors = adjacency_matrix[vid]; int num_neighb = neighbors.size(); // Find cover of current point (vid). if (cover[vid].size() == 2) v1 = std::min(cover[vid][0], cover[vid][1]); else v1 = cover[vid][0]; -vector node(1); node[0] = v1; +std::vector node(1); node[0] = v1; simplices.push_back(node); // Loop on neighbors. @@ -1237,15 +1243,15 @@ for (int i = 0; i < num_neighb; i++) { // If neighbor is in next interval, add edge. if (cover_fct[v2] == cover_fct[v1] + 1) { - vector edge(2); edge[0] = v1; edge[1] = v2; + std::vector edge(2); edge[0] = v1; edge[1] = v2; simplices.push_back(edge); break; } } } - std::vector dist(n); - std::vector process(n); + std::std::vector dist(n); + std::std::vector process(n); for (int j = 0; j < n; j++) { dist[j] = std::numeric_limits::max(); process[j] = j; @@ -1254,7 +1260,7 @@ for (int i = 0; i < num_neighb; i++) { int curr_size = process.size(); int min_point, min_index; double min_dist; - std::vector neighbors; + std::std::vector neighbors; int num_neighbors; while (curr_size > 0) { @@ -1282,18 +1288,18 @@ for (int i = 0; i < num_neighb; i++) { // Compute the connected components with DFS - std::map visit; - if (verbose) std::cout << "Preimage of interval " << i << std::endl; - for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) + std::std::map visit; + if (verbose) std::std::cout << "Preimage of interval " << i << std::std::endl; + for (std::std::map >::iterator it = prop.begin(); it != prop.end(); it++) visit[it->first] = false; if (!(prop.empty())) { - for (std::map >::iterator it = prop.begin(); it != prop.end(); it++) { + for (std::std::map >::iterator it = prop.begin(); it != prop.end(); it++) { if (!(visit[it->first])) { - std::vector cc; + std::std::vector cc; cc.clear(); dfs(prop, it->first, cc, visit); int cci = cc.size(); - if (verbose) std::cout << "one CC with " << cci << " points, "; + if (verbose) std::std::cout << "one CC with " << cci << " points, "; double average_col = 0; for (int j = 0; j < cci; j++) { cover[cc[j]].push_back(id); @@ -1301,8 +1307,8 @@ for (int i = 0; i < num_neighb; i++) { average_col += func_color[cc[j]] / cci; } cover_fct[id] = i; - cover_std[id] = std::pair(cci, 0.5*(u+v)); - cover_color[id] = std::pair(cci, average_col); + cover_std[id] = std::std::pair(cci, 0.5*(u+v)); + cover_color[id] = std::std::pair(cci, average_col); id++; } } @@ -1311,7 +1317,7 @@ for (int i = 0; i < num_neighb; i++) { // DFS private: - void dfs(std::map >& G, int p, std::vector& cc, std::map& visit) { + void dfs(std::std::map >& G, int p, std::std::vector& cc, std::std::map& visit) { cc.push_back(p); visit[p] = true; int neighb = G[p].size(); -- cgit v1.2.3 From ea53f69c13512465be3371f6378b5367f62164ab Mon Sep 17 00:00:00 2001 From: mcarrier Date: Wed, 20 Dec 2017 18:36:10 +0000 Subject: added extended persistence git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3092 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: a20ecf24407ce88654ad01257d8e44f8395c6b62 --- src/Nerve_GIC/example/FuncGIC.cpp | 2 +- src/Nerve_GIC/example/Nerve.cpp | 2 +- src/Nerve_GIC/include/gudhi/GIC.h | 168 +++++++++++++++++++++++--------------- 3 files changed, 104 insertions(+), 68 deletions(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp index 3583c66f..027ff525 100644 --- a/src/Nerve_GIC/example/FuncGIC.cpp +++ b/src/Nerve_GIC/example/FuncGIC.cpp @@ -71,7 +71,7 @@ int main(int argc, char **argv) { Gudhi::Simplex_tree<> stree; GIC.create_complex(stree); - GIC.compute_PD >(); + GIC.compute_PD(); // -------------------------------------------- // Display information about the functional GIC diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp index 7634c6f4..aaf33cff 100644 --- a/src/Nerve_GIC/example/Nerve.cpp +++ b/src/Nerve_GIC/example/Nerve.cpp @@ -74,7 +74,7 @@ int main(int argc, char **argv) { Gudhi::Simplex_tree<> stree; SC.create_complex(stree); - SC.compute_PD >(); + SC.compute_PD(); // ---------------------------------------------------------------------------- // Display information about the graph induced complex diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 3cd8a92a..e72480b7 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -31,7 +31,7 @@ #include #include #include -#include +//#include #include #include @@ -124,6 +124,8 @@ class Cover_complex { double rate_power = 0.001; // Power in the subsampling. int mask = 0; // Ignore nodes containing less than mask points. + std::map name2id, name2idinv; + std::string cover_name; std::string point_cloud_name; std::string color_name; @@ -600,66 +602,72 @@ class Cover_complex { std::vector points(n); for (int i = 0; i < n; i++) points[i] = i; std::sort(points.begin(), points.end(), Less(this->func)); - int id = 0; int pos = 0; int maxc = -1; IndexMap index = boost::get(boost::vertex_index, one_skeleton); + int id = 0; int pos = 0; IndexMap index = boost::get(boost::vertex_index, one_skeleton); //int maxc = -1; + std::map > preimages; std::map funcstd; - for (int i = 0; i < res; i++) { + if(verbose) std::cout << "Computing preimages..." << std::endl; + for (int i = 0; i < res; i++){ // Find points in the preimage - std::vector indices; std::pair inter1 = intervals[i]; - int tmp = pos; double u, v; Graph G = one_skeleton.create_subgraph(); + std::pair inter1 = intervals[i]; int tmp = pos; double u, v; if (i != res - 1) { if (i != 0) { std::pair inter3 = intervals[i - 1]; - while (func[points[tmp]] < inter3.second && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } + while (func[points[tmp]] < inter3.second && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} u = inter3.second; } else u = inter1.first; std::pair inter2 = intervals[i + 1]; - while (func[points[tmp]] < inter2.first && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } - + while (func[points[tmp]] < inter2.first && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} v = inter2.first; - pos = tmp; - while (func[points[tmp]] < inter1.second && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } + while (func[points[tmp]] < inter1.second && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} } else { - std::pair inter3 = intervals[i - 1]; - while (func[points[tmp]] < inter3.second && tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } - - while (tmp != n){ - boost::add_vertex(index[vertices[points[tmp]]], G); indices.push_back(points[tmp]); tmp++; - } + std::pair inter3 = intervals[i - 1]; + while (func[points[tmp]] < inter3.second && tmp != n){preimages[i].push_back(points[tmp]); tmp++;} + while (tmp != n){preimages[i].push_back(points[tmp]); tmp++;} u = inter3.second; v = inter1.second; } - int num = num_vertices(G); std::vector component(num); + funcstd[i] = 0.5*(u+v); + + } + + if(verbose) std::cout << "Computing connected components..." << std::endl; + for (int i = 0; i < res; i++){ // Compute connected components - boost::connected_components(G, &component[0]); int maxct = maxc + 1; + 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; - // Update covers + // For each point in preimage for(int j = 0; j < num; j++){ - maxc = std::max(maxc, maxct + component[j]); - cover [indices[j]] .push_back(maxct + component[j]); - cover_back [maxct + component[j]] .push_back(indices[j]); - cover_fct [maxct + component[j]] = i; - cover_std [maxct + component[j]] = 0.5*(u+v); - cover_color [maxct + component[j]] .second += func_color[indices[j]]; - cover_color [maxct + component[j]] .first += 1; + + // 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; + } maximal_dim = id - 1; @@ -747,7 +755,7 @@ class Cover_complex { * @result cover_back(c) vector of IDs of data points. * */ - const std::vector & subpopulation(int c) { return cover_back[c]; } + const std::vector & subpopulation(int c) { return cover_back[name2idinv[c]]; } // ******************************************************************************************************************* // Visualization. @@ -767,9 +775,7 @@ class Cover_complex { double f; while (std::getline(input, line)) { std::stringstream stream(line); - //stream >> one_skeleton[vertices[i]].color; - stream >> f; - func_color.emplace(i, f); + stream >> f; func_color.emplace(i, f); i++; } color_name = color_file_name; @@ -813,11 +819,11 @@ class Cover_complex { int k = 0; std::vector nodes; nodes.clear(); - graphic << "graph GIC {" << std::endl; + graphic << "graph GIC {" << std::endl; int id = 0; for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) { if (iit->second.first > mask) { - nodes.push_back(iit->first); - graphic << iit->first << "[shape=circle fontcolor=black color=black label=\"" << iit->first << ":" + nodes.push_back(iit->first); name2id[iit->first] = id; name2idinv[id] = iit->first; id++; + graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first] << ":" << iit->second.first << "\" style=filled fillcolor=\"" << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl; k++; @@ -828,7 +834,7 @@ class Cover_complex { for (int i = 0; i < num_simplices; i++) if (simplices[i].size() == 2) { if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) { - graphic << " " << simplices[i][0] << " -- " << simplices[i][1] << " [weight=15];" << std::endl; + graphic << " " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];" << std::endl; ke++; } } @@ -856,13 +862,14 @@ class Cover_complex { graphic << resolution_double << " " << gain << std::endl; graphic << cover_color.size() << " " << num_edges << std::endl; - for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) - graphic << iit->first << " " << iit->second.second << " " << iit->second.first << std::endl; + int id = 0; + for (std::map >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++){ + graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl; name2id[iit->first] = id; name2idinv[id] = iit->first; id++;} for (int i = 0; i < num_simplices; i++) if (simplices[i].size() == 2) if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) - graphic << simplices[i][0] << " " << simplices[i][1] << std::endl; + graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl; graphic.close(); std::cout << ".txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; @@ -928,35 +935,64 @@ class Cover_complex { complex.insert_simplex_and_subfaces(simplex, filt); if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; } - complex.set_dimension(dimension); + //complex.set_dimension(dimension); } public: /** \brief Computes the extended persistence diagram of the complex. * */ - template + //template void compute_PD(){ - SimplicialComplex streef, streeb; unsigned int dimension = 0; - for (auto const & simplex : simplices) { - int numvert = simplex.size(); double filtM = std::numeric_limits::lowest(); double filtm = filtM; - for(int i = 0; i < numvert; i++){filtM = std::max(cover_std[simplex[i]], filtM); filtm = std::max(-cover_std[simplex[i]], filtm);} - streef.insert_simplex_and_subfaces(simplex, filtM); streeb.insert_simplex_and_subfaces(simplex, filtm); - if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; - } streef.set_dimension(dimension); streeb.set_dimension(dimension); + Simplex_tree sc; - streef.initialize_filtration(); - Gudhi::persistent_cohomology::Persistent_cohomology pcohf(streef); - pcohf.init_coefficients(2); pcohf.compute_persistent_cohomology(); - pcohf.output_diagram(); + // Add magic point + std::vector node = {-2}; sc.insert_simplex(node); + + // Compute max and min + double maxf = std::numeric_limits::lowest(); double minf = std::numeric_limits::max(); + for(std::map::iterator it = cover_std.begin(); it != cover_std.end(); it++){ + maxf = std::max(maxf, it->second); minf = std::min(minf, it->second); + } - streeb.initialize_filtration(); - Gudhi::persistent_cohomology::Persistent_cohomology pcohb(streeb); - pcohb.init_coefficients(2); pcohb.compute_persistent_cohomology(); - pcohb.output_diagram(); + // Build filtration + for (auto const & simplex : simplices) { + int numvert = simplex.size(); double filta = std::numeric_limits::lowest(); double filts = filta; + for(int i = 0; i < numvert; i++) filta = std::max( -2 + (cover_std[simplex[i]] - minf)/(maxf - minf), filta ); + for(int i = 0; i < numvert; i++) filts = std::max( 2 - (cover_std[simplex[i]] - minf)/(maxf - minf), filts ); + std::vector splx = simplex; splx.push_back(node[0]); + sc.insert_simplex_and_subfaces(simplex, filta); sc.insert_simplex_and_subfaces(splx, filts); + //for(int i = 0; i < numvert; i++) std::cout << simplex[i] << ", "; std::cout << std::endl << " with filter value " << filta << std::endl; + //for(int i = 0; i < numvert + 1; i++) std::cout << splx[i] << ", "; std::cout << std::endl << " with filter value " << filts << std::endl; + } - //PD = pcohf.get_persistent_pairs(); + // Compute PD + sc.initialize_filtration(); Gudhi::persistent_cohomology::Persistent_cohomology pcoh(sc); + pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); + + // Output PD + std::vector > bars = pcoh.intervals_in_dimension(0); double birth_cc, death_cc; + int num_bars = bars.size(); std::cout << num_bars - 1 << " interval(s) in dimension 0:" << std::endl; + bool found_birth = false; bool found_death = false; + for(int j = 0; j < num_bars; j++){ + double birth = bars[j].first; double death = bars[j].second; + if(birth == 0){ std::cout << birth << " " << death << std::endl; if(found_birth) birth = birth_cc; else{death_cc = death; found_death = true; continue;} } + if(std::isinf(death)){ std::cout << birth << " " << death << std::endl; if(found_death) death = death_cc; else{birth_cc = birth; found_birth = true; continue;} } + if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); + if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); + std::cout << " " << birth << ", " << death << std::endl; + } + + for(int i = 1; i < sc.dimension(); i++){ + bars = pcoh.intervals_in_dimension(i); num_bars = bars.size(); std::cout << 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; + if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); + if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); + std::cout << " " << birth << ", " << death << std::endl; + } + } } @@ -987,9 +1023,9 @@ class Cover_complex { Cboot.set_graph_from_automatic_rips(Gudhi::Euclidean_distance()); Cboot.set_automatic_resolution(); Cboot.set_gain(); Cboot.set_cover_from_function(); - Cboot.find_simplices(); Cboot.compute_PD >(); + Cboot.find_simplices(); Cboot.compute_PD(); - distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); + //distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); } @@ -1069,7 +1105,7 @@ class Cover_complex { for(int j = 0; j < numt; j++){ int vt = cover[index[boost::target(*ei, one_skeleton)]][j]; if(cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1){ - std::vector edge(2); edge[0] = vs; edge[1] = vt; simplices.push_back(edge); goto afterLoop; + std::vector edge(2); edge[0] = std::min(vs,vt); edge[1] = std::max(vs,vt); simplices.push_back(edge); goto afterLoop; } } } -- cgit v1.2.3 From caef9f19009ddb82d8aefee55af345db79c1363f Mon Sep 17 00:00:00 2001 From: mcarrier Date: Thu, 21 Dec 2017 14:51:15 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3094 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 277086858a1634a27fc271f18b3d85f2b669f9bd --- src/Nerve_GIC/include/gudhi/GIC.h | 101 ++++++++++++++++++++++++-------------- 1 file changed, 63 insertions(+), 38 deletions(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index e72480b7..0441735d 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -942,58 +942,83 @@ class Cover_complex { /** \brief Computes the extended persistence diagram of the complex. * */ - //template void compute_PD(){ - Simplex_tree sc; + Graph G; - // Add magic point - std::vector node = {-2}; sc.insert_simplex(node); - // Compute max and min - double maxf = std::numeric_limits::lowest(); double minf = std::numeric_limits::max(); + std::map transf, transf_inv; + double maxf = std::numeric_limits::lowest(); double minf = std::numeric_limits::max(); int id = 0; for(std::map::iterator it = cover_std.begin(); it != cover_std.end(); it++){ - maxf = std::max(maxf, it->second); minf = std::min(minf, it->second); + maxf = std::max(maxf, it->second); minf = std::min(minf, it->second); boost::add_vertex(id, G); transf[id] = it->first; transf_inv[it->first] = id; id++; } - // Build filtration - for (auto const & simplex : simplices) { - int numvert = simplex.size(); double filta = std::numeric_limits::lowest(); double filts = filta; - for(int i = 0; i < numvert; i++) filta = std::max( -2 + (cover_std[simplex[i]] - minf)/(maxf - minf), filta ); - for(int i = 0; i < numvert; i++) filts = std::max( 2 - (cover_std[simplex[i]] - minf)/(maxf - minf), filts ); - std::vector splx = simplex; splx.push_back(node[0]); - sc.insert_simplex_and_subfaces(simplex, filta); sc.insert_simplex_and_subfaces(splx, filts); - //for(int i = 0; i < numvert; i++) std::cout << simplex[i] << ", "; std::cout << std::endl << " with filter value " << filta << std::endl; - //for(int i = 0; i < numvert + 1; i++) std::cout << splx[i] << ", "; std::cout << std::endl << " with filter value " << filts << std::endl; + // Build graph + for(auto const & simplex : simplices){ + int numvert = simplex.size(); + for(int i = 0; i < numvert; i++) for(int j = 0; j < numvert; j++) boost::add_edge(transf_inv[simplex[i]], transf_inv[simplex[j]], G); } - // Compute PD - sc.initialize_filtration(); Gudhi::persistent_cohomology::Persistent_cohomology pcoh(sc); - pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); - - // Output PD - std::vector > bars = pcoh.intervals_in_dimension(0); double birth_cc, death_cc; - int num_bars = bars.size(); std::cout << num_bars - 1 << " interval(s) in dimension 0:" << std::endl; - bool found_birth = false; bool found_death = false; - for(int j = 0; j < num_bars; j++){ - double birth = bars[j].first; double death = bars[j].second; - if(birth == 0){ std::cout << birth << " " << death << std::endl; if(found_birth) birth = birth_cc; else{death_cc = death; found_death = true; continue;} } - if(std::isinf(death)){ std::cout << birth << " " << death << std::endl; if(found_death) death = death_cc; else{birth_cc = birth; found_birth = true; continue;} } - if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); - if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); - std::cout << " " << birth << ", " << death << std::endl; + // Find connected components of G + int numvert = cover_std.size(); std::vector component(numvert); boost::connected_components(G, &component[0]); + int numcomp = 0; for(int i = 0; i < numvert; i++) numcomp = std::max(numcomp, component[i]); numcomp += 1; + + std::map complexes; + for(auto const & simplex : simplices){ + // Identify connected component + int compid = component[transf_inv[simplex[0]]]; + // Add simplices + complexes[compid].insert_simplex_and_subfaces(simplex); + // Add cone on simplices + std::vector splx = simplex; splx.push_back(-2); complexes[compid].insert_simplex_and_subfaces(splx); } + + + for(std::map::iterator it = complexes.begin(); it != complexes.end(); it++){ + + if(verbose) std::cout << "PD of connected component " << it->first << ":" << std::endl; + Simplex_tree sc = it->second; + + // Build filtration + for(auto simplex : sc.complex_simplex_range()){ + double filta = std::numeric_limits::lowest(); double filts = filta; bool ascending = true; + for(auto vertex : sc.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) sc.assign_filtration(simplex, filta); else sc.assign_filtration(simplex, filts); std::vector magic = {-2}; sc.assign_filtration(sc.find(magic), 0); + } + + // Compute PD + sc.initialize_filtration(); Gudhi::persistent_cohomology::Persistent_cohomology pcoh(sc); + pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); - for(int i = 1; i < sc.dimension(); i++){ - bars = pcoh.intervals_in_dimension(i); num_bars = bars.size(); std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl; + // Output PD + std::vector > bars = pcoh.intervals_in_dimension(0); double birth_cc, death_cc; + int num_bars = bars.size(); std::cout << num_bars - 1 << " interval(s) in dimension 0:" << std::endl; + bool found_birth = false; bool found_death = false; for(int j = 0; j < num_bars; j++){ double birth = bars[j].first; double death = bars[j].second; - if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); - if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); - std::cout << " " << birth << ", " << death << std::endl; + if(birth == 0){ if(found_birth) birth = birth_cc; else{death_cc = death; found_death = true; continue;} } + if(std::isinf(death)){ if(found_death) death = death_cc; else{birth_cc = birth; found_birth = true; continue;} } + if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); + if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); + PD.push_back(std::pair(birth, death)); + if(verbose) std::cout << " " << birth << ", " << death << std::endl; } - } + for(int i = 1; i < sc.dimension(); i++){ + bars = pcoh.intervals_in_dimension(i); num_bars = bars.size(); std::cout << 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; + if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); + if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); + PD.push_back(std::pair(birth, death)); + if(verbose) std::cout << " " << birth << ", " << death << std::endl; + } + } + } } public: @@ -1025,7 +1050,7 @@ class Cover_complex { Cboot.set_automatic_resolution(); Cboot.set_gain(); Cboot.set_cover_from_function(); Cboot.find_simplices(); Cboot.compute_PD(); - //distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); + distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); } -- cgit v1.2.3 From 2e1d09e3bfd7203c92a11471eab9e6bbcdd36944 Mon Sep 17 00:00:00 2001 From: mcarrier Date: Fri, 22 Dec 2017 11:03:28 +0000 Subject: cleaned the code git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3096 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 0a24fb8ff0300dc06c48fdb163d053387ae21c29 --- biblio/how_to_cite_gudhi.bib | 9 ++ src/Nerve_GIC/doc/funcGICvisu.jpg | Bin 71647 -> 68388 bytes src/Nerve_GIC/doc/funcGICvisu.pdf | Bin 0 -> 11347 bytes src/Nerve_GIC/example/FuncGIC.cpp | 1 - src/Nerve_GIC/example/Nerve.txt | 86 ++++++---- src/Nerve_GIC/include/gudhi/GIC.h | 290 +++++++++++++++++++++++----------- src/cmake/modules/GUDHI_modules.cmake | 6 +- 7 files changed, 261 insertions(+), 131 deletions(-) create mode 100644 src/Nerve_GIC/doc/funcGICvisu.pdf (limited to 'src/Nerve_GIC') diff --git a/biblio/how_to_cite_gudhi.bib b/biblio/how_to_cite_gudhi.bib index 59c05a5b..28014ade 100644 --- a/biblio/how_to_cite_gudhi.bib +++ b/biblio/how_to_cite_gudhi.bib @@ -122,3 +122,12 @@ , url = "http://gudhi.gforge.inria.fr/python/latest/" , year = 2016 } + +@incollection{gudhi:CoverComplex +, author = "Mathieu Carri\`ere" +, title = "Cover complex" +, publisher = "{GUDHI Editorial Board}" +, booktitle = "{GUDHI} User and Reference Manual" +, url = "http://gudhi.gforge.inria.fr/doc/latest/group__bottleneck__distance.html" +, year = 2017 +} diff --git a/src/Nerve_GIC/doc/funcGICvisu.jpg b/src/Nerve_GIC/doc/funcGICvisu.jpg index f3da45ac..36b64dde 100644 Binary files a/src/Nerve_GIC/doc/funcGICvisu.jpg and b/src/Nerve_GIC/doc/funcGICvisu.jpg differ diff --git a/src/Nerve_GIC/doc/funcGICvisu.pdf b/src/Nerve_GIC/doc/funcGICvisu.pdf new file mode 100644 index 00000000..d7456cd3 Binary files /dev/null and b/src/Nerve_GIC/doc/funcGICvisu.pdf differ diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp index 027ff525..3762db4e 100644 --- a/src/Nerve_GIC/example/FuncGIC.cpp +++ b/src/Nerve_GIC/example/FuncGIC.cpp @@ -71,7 +71,6 @@ int main(int argc, char **argv) { Gudhi::Simplex_tree<> stree; GIC.create_complex(stree); - GIC.compute_PD(); // -------------------------------------------- // Display information about the functional GIC diff --git a/src/Nerve_GIC/example/Nerve.txt b/src/Nerve_GIC/example/Nerve.txt index 2a861c5f..839ff45e 100644 --- a/src/Nerve_GIC/example/Nerve.txt +++ b/src/Nerve_GIC/example/Nerve.txt @@ -1,43 +1,63 @@ +Min function value = -0.979672 and Max function value = 0.816414 +Interval 0 = [-0.979672, -0.761576] +Interval 1 = [-0.838551, -0.581967] +Interval 2 = [-0.658942, -0.402359] +Interval 3 = [-0.479334, -0.22275] +Interval 4 = [-0.299725, -0.0431415] +Interval 5 = [-0.120117, 0.136467] +Interval 6 = [0.059492, 0.316076] +Interval 7 = [0.239101, 0.495684] +Interval 8 = [0.418709, 0.675293] +Interval 9 = [0.598318, 0.816414] +Computing preimages... +Computing connected components... +.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox. +5 interval(s) in dimension 0: + [-0.909111, 0.00817529] + [-0.171433, 0.367392] + [-0.171433, 0.367392] + [-0.909111, 0.745853] +0 interval(s) in dimension 1: Nerve is of dimension 1 - 41 simplices - 21 vertices. Iterator on Nerve simplices -0 1 -2 -2 0 -3 -3 1 +0 4 -4 3 -5 -5 2 -6 -6 4 -7 -7 5 +4 0 +2 +2 1 8 +8 2 +5 +5 4 9 -9 6 -10 -10 7 -11 -12 -12 8 +9 8 13 -13 9 -13 10 +13 5 14 -14 11 -15 -15 13 -16 -16 12 -17 -17 14 -18 -18 15 -18 16 -18 17 +14 9 19 -19 18 +19 13 +25 +32 20 -20 19 +32 20 +33 +33 25 +26 +26 14 +26 19 +42 +42 26 +34 +34 33 +27 +27 20 +35 +35 27 +35 34 +42 35 +44 +44 35 +54 +54 44 \ No newline at end of file diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 0441735d..3dd28841 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -89,32 +89,32 @@ class Cover_complex { private: bool verbose = false; // whether to display information. + std::string type; // Nerve or GIC - std::vector point_cloud; - int maximal_dim; // maximal dimension of output simplicial complex. - int data_dimension; // dimension of input data. - int n; // number of points. - - std::vector > distances; + std::vector point_cloud; // input point cloud. + std::vector > distances; // all pairwise distances. + int maximal_dim; // maximal dimension of output simplicial 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. - bool functional_cover = false; // whether we use a cover with preimages of a function or not. + 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). - Graph one_skeleton; // one-skeleton used to compute the connected components. - std::vector vertices; - std::vector > simplices; + Graph one_skeleton_OFF; // one-skeleton given by the input OFF file (if it exists). + Graph one_skeleton; // one-skeleton used to compute the connected components. + std::vector vertices; // vertices of one_skeleton. - std::vector voronoi_subsamples; + std::vector > simplices; // simplices of output simplicial complex. + std::vector voronoi_subsamples; // Voronoi germs (in case of Voronoi cover). PersistenceDiagram PD; std::vector distribution; - std::map > cover; - std::map > cover_back; - std::map cover_std; // standard function (induced by func) used to compute the extended persistence diagram of the output simplicial complex. - std::map cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. + std::map > cover; // function associating to each data point its vectors 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 diagram of the output simplicial complex. + std::map cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. std::map > cover_color; // size and coloring (induced by func_color) of the vertices of the output simplicial complex. int resolution_int = -1; @@ -129,7 +129,15 @@ class Cover_complex { std::string cover_name; std::string point_cloud_name; std::string color_name; - std::string type; // Nerve or GIC + + + + + + + + + // Point comparator struct Less { @@ -167,6 +175,14 @@ class Cover_complex { } + + + + + // ******************************************************************************************************************* + // Utils. + // ******************************************************************************************************************* + public: /** \brief Specifies whether the type of the output simplicial complex. * @@ -273,6 +289,23 @@ class Cover_complex { return input.is_open(); } + + + + + + + + + + + + + + + + + // ******************************************************************************************************************* // Graphs. // ******************************************************************************************************************* @@ -408,6 +441,29 @@ class Cover_complex { return delta; } + + + + + + + + + + + + + + + + + + + + + + + // ******************************************************************************************************************* // Functions. // ******************************************************************************************************************* @@ -453,6 +509,27 @@ class Cover_complex { functional_cover = true; } + + + + + + + + + + + + + + + + + + + + + // ******************************************************************************************************************* // Covers. // ******************************************************************************************************************* @@ -640,6 +717,7 @@ class Cover_complex { } if(verbose) std::cout << "Computing connected components..." << std::endl; + // #pragma omp parallel for for (int i = 0; i < res; i++){ // Compute connected components @@ -724,6 +802,7 @@ class Cover_complex { std::vector mindist(n); 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; @@ -757,6 +836,18 @@ class Cover_complex { */ const std::vector & subpopulation(int c) { return cover_back[name2idinv[c]]; } + + + + + + + + + + + + // ******************************************************************************************************************* // Visualization. // ******************************************************************************************************************* @@ -917,110 +1008,85 @@ class Cover_complex { std::cout << ".off generated. It can be visualized with e.g. geomview." << std::endl; } + + + + + + + + + + + + + + + + + + + + + // ******************************************************************************************************************* + // Extended Persistence Diagrams. // ******************************************************************************************************************* - public: - /** \brief Creates the simplicial complex. - * - * @param[in] complex SimplicialComplex to be created. - * - */ - template - void create_complex(SimplicialComplex& complex){ - unsigned int dimension = 0; - for (auto const& simplex : simplices) { - int numvert = simplex.size(); double filt = std::numeric_limits::lowest(); - for(int i = 0; i < numvert; i++) filt = std::max(cover_color[simplex[i]].second, filt); - complex.insert_simplex_and_subfaces(simplex, filt); - if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; - } - //complex.set_dimension(dimension); - } - public: /** \brief Computes the extended persistence diagram of the complex. * */ void compute_PD(){ - Graph G; + Simplex_tree st; // Compute max and min - std::map transf, transf_inv; - double maxf = std::numeric_limits::lowest(); double minf = std::numeric_limits::max(); int id = 0; + double maxf = std::numeric_limits::lowest(); double minf = std::numeric_limits::max(); for(std::map::iterator it = cover_std.begin(); it != cover_std.end(); it++){ - maxf = std::max(maxf, it->second); minf = std::min(minf, it->second); boost::add_vertex(id, G); transf[id] = it->first; transf_inv[it->first] = id; id++; + maxf = std::max(maxf, it->second); minf = std::min(minf, it->second); } - - // Build graph - for(auto const & simplex : simplices){ - int numvert = simplex.size(); - for(int i = 0; i < numvert; i++) for(int j = 0; j < numvert; j++) boost::add_edge(transf_inv[simplex[i]], transf_inv[simplex[j]], G); - } - - // Find connected components of G - int numvert = cover_std.size(); std::vector component(numvert); boost::connected_components(G, &component[0]); - int numcomp = 0; for(int i = 0; i < numvert; i++) numcomp = std::max(numcomp, component[i]); numcomp += 1; - - std::map complexes; + for(auto const & simplex : simplices){ - // Identify connected component - int compid = component[transf_inv[simplex[0]]]; // Add simplices - complexes[compid].insert_simplex_and_subfaces(simplex); + st.insert_simplex_and_subfaces(simplex); // Add cone on simplices - std::vector splx = simplex; splx.push_back(-2); complexes[compid].insert_simplex_and_subfaces(splx); + std::vector splx = simplex; splx.push_back(-2); st.insert_simplex_and_subfaces(splx); } - - - for(std::map::iterator it = complexes.begin(); it != complexes.end(); it++){ - - if(verbose) std::cout << "PD of connected component " << it->first << ":" << std::endl; - Simplex_tree sc = it->second; - // Build filtration - for(auto simplex : sc.complex_simplex_range()){ - double filta = std::numeric_limits::lowest(); double filts = filta; bool ascending = true; - for(auto vertex : sc.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) sc.assign_filtration(simplex, filta); else sc.assign_filtration(simplex, filts); std::vector magic = {-2}; sc.assign_filtration(sc.find(magic), 0); + // 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); + } + std::vector magic = {-2}; st.assign_filtration(st.find(magic), -3); - // Compute PD - sc.initialize_filtration(); Gudhi::persistent_cohomology::Persistent_cohomology pcoh(sc); - pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); + // Compute PD + st.initialize_filtration(); Gudhi::persistent_cohomology::Persistent_cohomology pcoh(st); + pcoh.init_coefficients(2); pcoh.compute_persistent_cohomology(); - // Output PD - std::vector > bars = pcoh.intervals_in_dimension(0); double birth_cc, death_cc; - int num_bars = bars.size(); std::cout << num_bars - 1 << " interval(s) in dimension 0:" << std::endl; - bool found_birth = false; bool found_death = false; + // 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(); std::cout << 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; - if(birth == 0){ if(found_birth) birth = birth_cc; else{death_cc = death; found_death = true; continue;} } - if(std::isinf(death)){ if(found_death) death = death_cc; else{birth_cc = birth; found_birth = true; continue;} } + if(i == 0 && std::isinf(death)) continue; if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); PD.push_back(std::pair(birth, death)); - if(verbose) std::cout << " " << birth << ", " << death << std::endl; - } - - for(int i = 1; i < sc.dimension(); i++){ - bars = pcoh.intervals_in_dimension(i); num_bars = bars.size(); std::cout << 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; - if(birth < 0) birth = minf + (birth + 2)*(maxf - minf); else birth = minf + (2 - birth)*(maxf - minf); - if(death < 0) death = minf + (death + 2)*(maxf - minf); else death = minf + (2 - death)*(maxf - minf); - PD.push_back(std::pair(birth, death)); - if(verbose) std::cout << " " << birth << ", " << death << std::endl; - } + if(verbose) std::cout << " [" << birth << ", " << death << "]" << std::endl; } } + } + public: /** \brief Computes bootstrapped distances distribution. * @@ -1050,7 +1116,7 @@ class Cover_complex { Cboot.set_automatic_resolution(); Cboot.set_gain(); Cboot.set_cover_from_function(); Cboot.find_simplices(); Cboot.compute_PD(); - distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); + //distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); } @@ -1060,7 +1126,7 @@ class Cover_complex { } public: - /** \brief Computes the bottleneck distance corresponding to a specific confidence level. + /** \brief Computes the bottleneck distance threshold corresponding to a specific confidence level. * * @param[in] alpha Confidence level. * @@ -1071,7 +1137,7 @@ class Cover_complex { } public: - /** \brief Computes the confidence level of a specific bottleneck distance. + /** \brief Computes the confidence level of a specific bottleneck distance threshold. * * @param[in] d Bottleneck distance. * @@ -1092,6 +1158,42 @@ class Cover_complex { return 1-compute_confidence_level_from_distance(distancemin); } + + + + + + + + + + + + + + + + // ******************************************************************************************************************* + // Computation of simplices. + // ******************************************************************************************************************* + + public: + /** \brief Creates the simplicial complex. + * + * @param[in] complex SimplicialComplex to be created. + * + */ + template + void create_complex(SimplicialComplex& complex){ + unsigned int dimension = 0; + for (auto const& simplex : simplices) { + int numvert = simplex.size(); double filt = std::numeric_limits::lowest(); + for(int i = 0; i < numvert; i++) filt = std::max(cover_color[simplex[i]].second, filt); + complex.insert_simplex_and_subfaces(simplex, filt); + if (dimension < simplex.size() - 1) dimension = simplex.size() - 1; + } + } + public: /** \brief Computes the simplices of the simplicial complex. */ diff --git a/src/cmake/modules/GUDHI_modules.cmake b/src/cmake/modules/GUDHI_modules.cmake index 74c5fb55..205ee8a1 100644 --- a/src/cmake/modules/GUDHI_modules.cmake +++ b/src/cmake/modules/GUDHI_modules.cmake @@ -16,10 +16,10 @@ function(add_gudhi_module file_path) endfunction(add_gudhi_module) -option(WITH_GUDHI_BENCHMARK "Activate/desactivate benchmark compilation" OFF) +option(WITH_GUDHI_BENCHMARK "Activate/desactivate benchmark compilation" ON) option(WITH_GUDHI_EXAMPLE "Activate/desactivate examples compilation and installation" ON) -option(WITH_GUDHI_PYTHON "Activate/desactivate python module compilation and installation" OFF) -option(WITH_GUDHI_TEST "Activate/desactivate examples compilation and installation" OFF) +option(WITH_GUDHI_PYTHON "Activate/desactivate python module compilation and installation" ON) +option(WITH_GUDHI_TEST "Activate/desactivate examples compilation and installation" ON) option(WITH_GUDHI_UTILITIES "Activate/desactivate utilities compilation and installation" ON) if (WITH_GUDHI_BENCHMARK) -- cgit v1.2.3 From b2fc8e3851008880595653fe50aee25a8af92e58 Mon Sep 17 00:00:00 2001 From: mcarrier Date: Fri, 22 Dec 2017 11:37:18 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3099 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: fa1505962b7b6140611935cd960edcdb0c48dd5c --- src/Nerve_GIC/include/gudhi/GIC.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h index 3dd28841..2bec27db 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -31,7 +31,7 @@ #include #include #include -//#include +#include #include #include @@ -1116,7 +1116,7 @@ class Cover_complex { Cboot.set_automatic_resolution(); Cboot.set_gain(); Cboot.set_cover_from_function(); Cboot.find_simplices(); Cboot.compute_PD(); - //distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); + distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD,Cboot.PD)); } -- cgit v1.2.3 From 3eb6aac7d217ff5004c3d50caac63fc1cb463e1c Mon Sep 17 00:00:00 2001 From: mcarrier Date: Wed, 3 Jan 2018 15:31:30 +0000 Subject: git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/Nerve_GIC@3112 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 0e3300cbc4d9c1375a7a162d175a982c3145c42b --- src/Nerve_GIC/example/CoordGIC.cpp | 2 +- src/Nerve_GIC/example/FuncGIC.cpp | 4 ++-- src/Nerve_GIC/example/GIC.cpp | 2 +- src/Nerve_GIC/example/Nerve.cpp | 2 +- src/Nerve_GIC/example/VoronoiGIC.cpp | 2 +- 5 files changed, 6 insertions(+), 6 deletions(-) (limited to 'src/Nerve_GIC') diff --git a/src/Nerve_GIC/example/CoordGIC.cpp b/src/Nerve_GIC/example/CoordGIC.cpp index c03fcbb3..0d06c2ba 100644 --- a/src/Nerve_GIC/example/CoordGIC.cpp +++ b/src/Nerve_GIC/example/CoordGIC.cpp @@ -28,7 +28,7 @@ void usage(int nbArgs, char *const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; std::cerr << "Usage: " << progName << " filename.off coordinate [--v] \n"; - std::cerr << " i.e.: " << progName << " ../../data/points/human.off 2 --v \n"; + std::cerr << " i.e.: " << progName << " ../../../../data/points/human.off 2 --v \n"; exit(-1); // ----- >> } diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp index 3762db4e..081b3a2f 100644 --- a/src/Nerve_GIC/example/FuncGIC.cpp +++ b/src/Nerve_GIC/example/FuncGIC.cpp @@ -28,8 +28,8 @@ void usage(int nbArgs, char *const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; std::cerr << "Usage: " << progName << " filename.off function [--v] \n"; - std::cerr << " i.e.: " << progName << " ../../data/points/COIL_database/lucky_cat.off " - "../../data/points/COIL_database/lucky_cat_PCA1 --v \n"; + std::cerr << " i.e.: " << progName << " ../../../../data/points/COIL_database/lucky_cat.off " + "../../../../data/points/COIL_database/lucky_cat_PCA1 --v \n"; exit(-1); // ----- >> } diff --git a/src/Nerve_GIC/example/GIC.cpp b/src/Nerve_GIC/example/GIC.cpp index 2bc24a4d..2f10fd17 100644 --- a/src/Nerve_GIC/example/GIC.cpp +++ b/src/Nerve_GIC/example/GIC.cpp @@ -28,7 +28,7 @@ void usage(int nbArgs, char *const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; std::cerr << "Usage: " << progName << " filename.off threshold coordinate resolution gain [--v] \n"; - std::cerr << " i.e.: " << progName << " ../../data/points/human.off 0.075 2 0.075 0 --v \n"; + std::cerr << " i.e.: " << progName << " ../../../../data/points/human.off 0.075 2 0.075 0 --v \n"; exit(-1); // ----- >> } diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp index aaf33cff..6460836d 100644 --- a/src/Nerve_GIC/example/Nerve.cpp +++ b/src/Nerve_GIC/example/Nerve.cpp @@ -30,7 +30,7 @@ using namespace std; void usage(int nbArgs, char *const progName) { cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; cerr << "Usage: " << progName << " filename.off coordinate resolution gain [--v] \n"; - cerr << " i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 --v \n"; + cerr << " i.e.: " << progName << " ../../../../data/points/human.off 2 10 0.3 --v \n"; exit(-1); // ----- >> } diff --git a/src/Nerve_GIC/example/VoronoiGIC.cpp b/src/Nerve_GIC/example/VoronoiGIC.cpp index 32431cc2..41ac9e5f 100644 --- a/src/Nerve_GIC/example/VoronoiGIC.cpp +++ b/src/Nerve_GIC/example/VoronoiGIC.cpp @@ -28,7 +28,7 @@ void usage(int nbArgs, char *const progName) { std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n"; std::cerr << "Usage: " << progName << " filename.off N [--v] \n"; - std::cerr << " i.e.: " << progName << " ../../data/points/human.off 100 --v \n"; + std::cerr << " i.e.: " << progName << " ../../../../data/points/human.off 100 --v \n"; exit(-1); // ----- >> } -- cgit v1.2.3