diff options
-rw-r--r-- | biblio/how_to_cite_gudhi.bib | 9 | ||||
-rw-r--r-- | src/Nerve_GIC/doc/Intro_graph_induced_complex.h | 8 | ||||
-rw-r--r-- | src/Nerve_GIC/doc/funcGICvisu.jpg | bin | 71647 -> 68388 bytes | |||
-rw-r--r-- | src/Nerve_GIC/doc/funcGICvisu.pdf | bin | 0 -> 11347 bytes | |||
-rw-r--r-- | src/Nerve_GIC/example/Nerve.cpp | 1 | ||||
-rw-r--r-- | src/Nerve_GIC/example/Nerve.txt | 86 | ||||
-rw-r--r-- | src/Nerve_GIC/include/gudhi/GIC.h | 767 |
7 files changed, 519 insertions, 352 deletions
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/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h index 3a0d8154..7578cc53 100644 --- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h +++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h @@ -70,7 +70,7 @@ namespace cover_complex { * * When launching: * - * \code $> ./Nerve ../../../../data/points/human.off 2 10 0.3 --v + * \code $> ./Nerve ../../data/points/human.off 2 10 0.3 --v * \endcode * * the program output is: @@ -113,7 +113,7 @@ namespace cover_complex { * * When launching: * - * \code $> ./VoronoiGIC ../../../../data/points/human.off 700 --v + * \code $> ./VoronoiGIC ../../data/points/human.off 700 --v * \endcode * * the program outputs SC.off. Using e.g. @@ -146,7 +146,7 @@ namespace cover_complex { * * When launching: * - * \code $> ./CoordGIC ../../../../data/points/KleinBottle5D.off 0 --v + * \code $> ./CoordGIC ../../data/points/KleinBottle5D.off 0 --v * \endcode * * the program outputs SC.dot. Using e.g. @@ -201,7 +201,7 @@ namespace cover_complex { * * When launching: * - * \code $> ./GIC ../../../../data/points/human.off 0.075 2 0.075 0 --v + * \code $> ./GIC ../../data/points/human.off 0.075 2 0.075 0 --v * \endcode * * the program outputs SC.txt, which can be visualized with python and firefox as before: diff --git a/src/Nerve_GIC/doc/funcGICvisu.jpg b/src/Nerve_GIC/doc/funcGICvisu.jpg Binary files differindex f3da45ac..36b64dde 100644 --- a/src/Nerve_GIC/doc/funcGICvisu.jpg +++ b/src/Nerve_GIC/doc/funcGICvisu.jpg diff --git a/src/Nerve_GIC/doc/funcGICvisu.pdf b/src/Nerve_GIC/doc/funcGICvisu.pdf Binary files differnew file mode 100644 index 00000000..d7456cd3 --- /dev/null +++ b/src/Nerve_GIC/doc/funcGICvisu.pdf diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp index 4d5b009b..6abdedc7 100644 --- a/src/Nerve_GIC/example/Nerve.cpp +++ b/src/Nerve_GIC/example/Nerve.cpp @@ -72,6 +72,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/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 9f107a7e..9dc754f4 100644 --- a/src/Nerve_GIC/include/gudhi/GIC.h +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -30,13 +30,23 @@ #include <gudhi/Rips_complex.h> #include <gudhi/Points_off_io.h> #include <gudhi/distance_functions.h> +#include <gudhi/Persistent_cohomology.h> +#include <gudhi/Bottleneck.h> + +#include <boost/config.hpp> +#include <boost/graph/graph_traits.hpp> +#include <boost/graph/adjacency_list.hpp> +#include <boost/graph/connected_components.hpp> +#include <boost/graph/dijkstra_shortest_paths.hpp> +#include <boost/graph/subgraph.hpp> +#include <boost/graph/graph_utility.hpp> #include <iostream> #include <vector> #include <map> #include <string> #include <limits> // for numeric_limits -#include <utility> // for pair<> +#include <utility> // for std::pair<> #include <algorithm> // for std::max #include <random> #include <cassert> @@ -48,6 +58,13 @@ namespace cover_complex { using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using PersistenceDiagram = std::vector<std::pair<double, double> >; +using Graph = boost::subgraph< + boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property, + boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >; +using vertex_t = boost::graph_traits<Graph>::vertex_descriptor; +using IndexMap = boost::property_map<Graph, boost::vertex_index_t>::type; +using WeightMap = boost::property_map<Graph, boost::edge_weight_t>::type; /** * \class Cover_complex @@ -72,25 +89,40 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; template <typename Point> class Cover_complex { private: - // Graph_induced_complex(std::map<int, double> fun){func = fun;} bool verbose = false; // whether to display information. - std::vector<Point> point_cloud; - std::vector<std::vector<int> > one_skeleton; - typedef int Cover_t; // elements of cover C are indexed by integers. - std::vector<std::vector<Cover_t> > simplices; - std::map<int, std::vector<Cover_t> > cover; - std::map<Cover_t, std::vector<int> > cover_back; - int maximal_dim; // maximal dimension of output simplicial complex. - int data_dimension; // dimension of input data. - int n; // number of points. - std::map<Cover_t, int> + std::string type; // Nerve or GIC + + std::vector<Point> point_cloud; // input point cloud. + std::vector<std::vector<double> > 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<int, double> func; // function used to compute the output simplicial complex. + std::map<int, double> + func_color; // function used to compute the colors of the nodes of the output simplicial complex. + 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<vertex_t> vertices; // vertices of one_skeleton. + + std::vector<std::vector<int> > simplices; // simplices of output simplicial complex. + std::vector<int> voronoi_subsamples; // Voronoi germs (in case of Voronoi cover). + + PersistenceDiagram PD; + std::vector<double> distribution; + + std::map<int, std::vector<int> > + cover; // function associating to each data point its vectors of cover elements to which it belongs. + std::map<int, std::vector<int> > + cover_back; // inverse of cover, in order to get the data points associated to a specific cover element. + std::map<int, double> cover_std; // standard function (induced by func) used to compute the extended persistence + // diagram of the output simplicial complex. + std::map<int, int> cover_fct; // integer-valued function that allows to state if two elements of the cover are consecutive or not. - std::map<Cover_t, std::pair<int, double> > - cover_color; // size and coloring of the vertices of the output simplicial complex. - Simplex_tree st; - - std::map<int, std::vector<int> > adjacency_matrix; - std::vector<std::vector<double> > distances; + std::map<int, std::pair<int, double> > + 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,14 +131,11 @@ class Cover_complex { double rate_power = 0.001; // Power in the subsampling. int mask = 0; // Ignore nodes containing less than mask points. - std::map<int, double> func; - std::map<int, double> func_color; - std::vector<int> voronoi_subsamples; + std::map<int, int> name2id, name2idinv; + 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 // Point comparator struct Less { @@ -120,21 +149,16 @@ class Cover_complex { } }; - // DFS - private: - void dfs(std::map<int, std::vector<int> >& G, int p, std::vector<int>& cc, std::map<int, bool>& 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); + // Remove all edges of a graph. + void remove_edges(Graph& G) { + boost::graph_traits<Graph>::edge_iterator ei, ei_end; + for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G); } // Find random number in [0,1]. double GetUniform() { - static std::default_random_engine re; - static std::uniform_real_distribution<double> Dist(0, 1); + thread_local std::default_random_engine re; + thread_local std::uniform_real_distribution<double> Dist(0, 1); return Dist(re); } @@ -145,9 +169,9 @@ class Cover_complex { double u; while (m < sampleSize) { u = GetUniform(); - if ((populationSize - t) * u >= sampleSize - m) { + if ((populationSize - t) * u >= sampleSize - m) t++; - } else { + else { samples[m] = t; t++; m++; @@ -155,24 +179,14 @@ class Cover_complex { } } - private: - void fill_adjacency_matrix_from_st() { - std::vector<int> 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<int> 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]); - } - } - } + // ******************************************************************************************************************* + // Utils. + // ******************************************************************************************************************* 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; } @@ -221,14 +235,15 @@ class Cover_complex { 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(' ')]; + std::getline(input, line); + 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)) + std::getline(input, line); + if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace)) comment = line[line.find_first_not_of(' ')]; } std::stringstream stream(line); @@ -238,10 +253,11 @@ class Cover_complex { } 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(' ')]; + std::getline(input, line); + if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace)) + comment = line[line.find_first_not_of(' ')]; } std::stringstream stream(line); stream >> n; @@ -250,34 +266,31 @@ class Cover_complex { i = 0; while (i < n) { - getline(input, line); + std::getline(input, line); if (!line.empty() && line[line.find_first_not_of(' ')] != '#' && - !std::all_of(line.begin(), line.end(), isspace)) { + !all_of(line.begin(), line.end(), (int (*)(int))isspace)) { + std::stringstream iss(line); std::vector<double> point; - std::istringstream iss(line); point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>()); point_cloud.emplace_back(point.begin(), point.begin() + data_dimension); + 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(' ')] != '#' && - !std::all_of(line.begin(), line.end(), isspace)) { + !all_of(line.begin(), line.end(), (int (*)(int))isspace)) { std::vector<int> simplex; - std::istringstream iss(line); + std::stringstream iss(line); simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>()); - num = simplex[0]; - std::vector<int> 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); - } - } + dim = simplex[0]; + for (int j = 1; j <= dim; j++) + for (int k = j + 1; k <= dim; k++) + boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF); i++; } } @@ -298,22 +311,16 @@ class Cover_complex { * */ 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 edge[2]; - int n = 0; + int source; 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++; + stream >> source; + while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton); } - - fill_adjacency_matrix_from_st(); } public: // Set graph from OFF file. @@ -321,13 +328,11 @@ 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 { + remove_edges(one_skeleton); + if (num_edges(one_skeleton_OFF)) + one_skeleton = one_skeleton_OFF; + else std::cout << "No triangulation read in OFF file!" << std::endl; - } } public: // Set graph from Rips complex. @@ -339,9 +344,27 @@ class Cover_complex { */ template <typename Distance> 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(); + 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) { + 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]); + } + } + } + } + + public: + void set_graph_weights() { + IndexMap index = boost::get(boost::vertex_index, one_skeleton); + WeightMap weight = boost::get(boost::edge_weight, one_skeleton); + boost::graph_traits<Graph>::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. @@ -397,7 +420,7 @@ class Cover_complex { */ template <typename Distance> 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)))); + int m = floor(n / exp((1 + rate_power) * log(log(n) / log(rate_constant)))); m = std::min(m, n - 1); std::vector<int> samples(m); double delta = 0; @@ -420,10 +443,7 @@ class Cover_complex { } 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(); - + set_graph_from_rips(delta, distance); return delta; } @@ -438,15 +458,15 @@ class Cover_complex { * */ void set_function_from_file(const std::string& func_file_name) { - int vertex_id = 0; + 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(vertex_id, f); - vertex_id++; + func.emplace(i, f); + i++; } functional_cover = true; cover_name = func_file_name; @@ -474,12 +494,8 @@ class Cover_complex { */ template <class InputRange> void set_function_from_range(InputRange const& function) { + for (int i = 0; i < n; i++) func.emplace(i, function[i]); functional_cover = true; - int index = 0; - for (auto v : function) { - func.emplace(index, v); - index++; - } } // ******************************************************************************************************************* @@ -505,27 +521,23 @@ class Cover_complex { } double reso = 0; + IndexMap index = boost::get(boost::vertex_index, one_skeleton); if (type == "GIC") { - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector<int> 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]])); - } - } + boost::graph_traits<Graph>::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") { - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector<int> 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); - } - } + boost::graph_traits<Graph>::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; } @@ -568,15 +580,12 @@ class Cover_complex { } // Read function values and compute min and max - std::map<int, double>::iterator it; - double maxf, minf; - minf = std::numeric_limits<float>::max(); - maxf = std::numeric_limits<float>::min(); - for (it = func.begin(); it != func.end(); it++) { - minf = std::min(minf, it->second); - maxf = std::max(maxf, it->second); + double minf = std::numeric_limits<float>::max(); + double maxf = std::numeric_limits<float>::lowest(); + for (int i = 0; i < n; i++) { + minf = std::min(minf, func[i]); + maxf = std::max(maxf, func[i]); } - int n = func.size(); if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; // Compute cover of im(f) @@ -648,78 +657,95 @@ class Cover_complex { std::vector<int> 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; + IndexMap index = boost::get(boost::vertex_index, one_skeleton); // int maxc = -1; + std::map<int, std::vector<int> > preimages; + std::map<int, double> funcstd; + if (verbose) std::cout << "Computing preimages..." << std::endl; for (int i = 0; i < res; i++) { // Find points in the preimage - std::map<int, std::vector<int> > prop; std::pair<double, double> inter1 = intervals[i]; int tmp = pos; + double u, v; if (i != res - 1) { if (i != 0) { std::pair<double, double> inter3 = intervals[i - 1]; while (func[points[tmp]] < inter3.second && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; + preimages[i].push_back(points[tmp]); tmp++; } - } + u = inter3.second; + } else + u = inter1.first; std::pair<double, double> inter2 = intervals[i + 1]; while (func[points[tmp]] < inter2.first && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; + preimages[i].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]]; + preimages[i].push_back(points[tmp]); tmp++; } } else { std::pair<double, double> inter3 = intervals[i - 1]; while (func[points[tmp]] < inter3.second && tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; + preimages[i].push_back(points[tmp]); tmp++; } - while (tmp != n) { - prop[points[tmp]] = adjacency_matrix[points[tmp]]; + preimages[i].push_back(points[tmp]); tmp++; } + u = inter3.second; + v = inter1.second; } - // Compute the connected components with DFS - std::map<int, bool> visit; - if (verbose) std::cout << "Preimage of interval " << i << std::endl; - for (std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++) - visit[it->first] = false; - if (!(prop.empty())) { - for (std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++) { - if (!(visit[it->first])) { - std::vector<int> 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<int, double>(cci, average_col); - id++; - } - } + funcstd[i] = 0.5 * (u + v); + } + + if (verbose) std::cout << "Computing connected components..." << std::endl; + // #pragma omp parallel for + for (int i = 0; i < res; i++) { + // Compute connected components + Graph G = one_skeleton.create_subgraph(); + int num = preimages[i].size(); + std::vector<int> component(num); + for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G); + boost::connected_components(G, &component[0]); + int max = 0; + + // For each point in preimage + for (int j = 0; j < num; j++) { + // Update number of components in preimage + if (component[j] > max) max = component[j]; + + // Identify component with Cantor polynomial N^2 -> N + int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2; + + // Update covers + cover[preimages[i][j]].push_back(identifier); + cover_back[identifier].push_back(preimages[i][j]); + cover_fct[identifier] = i; + cover_std[identifier] = funcstd[i]; + cover_color[identifier].second += func_color[preimages[i][j]]; + cover_color[identifier].first += 1; } - if (verbose) std::cout << std::endl; + + // Maximal dimension is total number of connected components + id += max + 1; } maximal_dim = id - 1; + for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) + iit->second.second /= iit->second.first; } public: // Set cover from file. @@ -730,9 +756,9 @@ class Cover_complex { * */ void set_cover_from_file(const std::string& cover_file_name) { - int vertex_id = 0; - Cover_t cov; - std::vector<Cover_t> cov_elts, cov_number; + int i = 0; + int cov; + std::vector<int> cov_elts, cov_number; std::ifstream input(cover_file_name); std::string line; while (std::getline(input, line)) { @@ -742,17 +768,18 @@ class Cover_complex { 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].second += func_color[i]; cover_color[cov].first++; - cover_back[cov].push_back(vertex_id); + cover_back[cov].push_back(i); } - cover[vertex_id] = cov_elts; - vertex_id++; + cover[i] = cov_elts; + i++; } - std::vector<Cover_t>::iterator it; + std::sort(cov_number.begin(), cov_number.end()); - it = std::unique(cov_number.begin(), cov_number.end()); + std::vector<int>::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; cover_name = cover_file_name; @@ -770,52 +797,25 @@ 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 = boost::get(boost::edge_weight, one_skeleton); + IndexMap index = boost::get(boost::vertex_index, one_skeleton); std::vector<double> mindist(n); for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max(); // Compute the geodesic distances to subsamples with Dijkstra + // #pragma omp parallel for for (int i = 0; i < m; i++) { if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl; int seed = voronoi_subsamples[i]; - std::vector<double> dist(n); - std::vector<int> process(n); - for (int j = 0; j < n; j++) { - dist[j] = std::numeric_limits<double>::max(); - process[j] = j; - } - dist[seed] = 0; - int curr_size = process.size(); - int min_point, min_index; - double min_dist; - std::vector<int> neighbors; - int num_neighbors; - - while (curr_size > 0) { - min_dist = std::numeric_limits<double>::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(); - } + std::vector<double> dmap(n); + boost::dijkstra_shortest_paths( + one_skeleton, vertices[seed], + boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index))); for (int j = 0; j < n; j++) - if (mindist[j] > dist[j]) { - mindist[j] = dist[j]; + if (mindist[j] > dmap[j]) { + mindist[j] = dmap[j]; if (cover[j].size() == 0) cover[j].push_back(i); else @@ -840,7 +840,7 @@ class Cover_complex { * @result cover_back(c) vector of IDs of data points. * */ - const std::vector<int>& subpopulation(Cover_t c) { return cover_back[c]; } + const std::vector<int>& subpopulation(int c) { return cover_back[name2idinv[c]]; } // ******************************************************************************************************************* // Visualization. @@ -854,15 +854,15 @@ class Cover_complex { * */ void set_color_from_file(const std::string& color_file_name) { - int vertex_id = 0; + int i = 0; std::ifstream input(color_file_name); std::string line; double f; while (std::getline(input, line)) { std::stringstream stream(line); stream >> f; - func_color.emplace(vertex_id, f); - vertex_id++; + func_color.emplace(i, f); + i++; } color_name = color_file_name; } @@ -874,7 +874,7 @@ 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)); } @@ -886,7 +886,7 @@ class Cover_complex { * */ void set_color_from_vector(std::vector<double> color) { - for (unsigned int i = 0; i < color.size(); i++) func_color.emplace(i, color[i]); + for (unsigned int i = 0; i < color.size(); i++) func_color[i] = color[i]; } public: // Create a .dot file that can be compiled with neato to produce a .pdf file. @@ -895,26 +895,31 @@ class Cover_complex { * of its 1-skeleton in a .pdf file. */ void plot_DOT() { - char mapp[11] = "SC.dot"; + char mapp[100]; + sprintf(mapp, "%s_sc.dot", point_cloud_name.c_str()); std::ofstream graphic(mapp); - graphic << "graph GIC {" << std::endl; - double maxv, minv; - maxv = std::numeric_limits<double>::min(); - minv = std::numeric_limits<double>::max(); - for (std::map<Cover_t, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); - iit++) { + + double maxv = std::numeric_limits<double>::lowest(); + double minv = std::numeric_limits<double>::max(); + for (std::map<int, std::pair<int, double> >::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; std::vector<int> nodes; nodes.clear(); - for (std::map<Cover_t, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); - iit++) { + + graphic << "graph GIC {" << std::endl; + int id = 0; + for (std::map<int, std::pair<int, double> >::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 << ":" - << iit->second.first << "\" style=filled fillcolor=\"" + 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++; } @@ -924,13 +929,14 @@ 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++; } } graphic << "}"; graphic.close(); - std::cout << "SC.dot generated. It can be visualized with e.g. neato." << std::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. @@ -940,8 +946,10 @@ class Cover_complex { void write_info() { int num_simplices = simplices.size(); int num_edges = 0; - char mapp[11] = "SC.txt"; + 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++; @@ -952,16 +960,20 @@ class Cover_complex { graphic << resolution_double << " " << gain << std::endl; graphic << cover_color.size() << " " << num_edges << std::endl; - for (std::map<Cover_t, std::pair<int, double> >::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<int, std::pair<int, double> >::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 << "SC.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." + std::cout << ".txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox." << std::endl; } @@ -972,14 +984,18 @@ class Cover_complex { */ 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<std::vector<int> > edges, faces; int numsimplices = simplices.size(); + + 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++) { if (simplices[i].size() == 2) { numedges++; @@ -1004,26 +1020,185 @@ 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; + std::cout << ".off generated. It can be visualized with e.g. geomview." << std::endl; + } + + // ******************************************************************************************************************* + // Extended Persistence Diagrams. + // ******************************************************************************************************************* + + public: + /** \brief Computes the extended persistence diagram of the complex. + * + */ + void compute_PD() { + Simplex_tree st; + + // Compute max and min + double maxf = std::numeric_limits<double>::lowest(); + double minf = std::numeric_limits<double>::max(); + for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) { + maxf = std::max(maxf, it->second); + minf = std::min(minf, it->second); + } + + for (auto const& simplex : simplices) { + // Add simplices + st.insert_simplex_and_subfaces(simplex); + // Add cone on simplices + std::vector<int> splx = simplex; + splx.push_back(-2); + st.insert_simplex_and_subfaces(splx); + } + + // Build filtration + for (auto simplex : st.complex_simplex_range()) { + double filta = std::numeric_limits<double>::lowest(); + double filts = filta; + bool ascending = true; + for (auto vertex : st.simplex_vertex_range(simplex)) { + if (vertex == -2) { + ascending = false; + continue; + } + filta = std::max(-2 + (cover_std[vertex] - minf) / (maxf - minf), filta); + filts = std::max(2 - (cover_std[vertex] - minf) / (maxf - minf), filts); + } + if (ascending) + st.assign_filtration(simplex, filta); + else + st.assign_filtration(simplex, filts); + } + std::vector<int> magic = {-2}; + st.assign_filtration(st.find(magic), -3); + + // Compute PD + st.initialize_filtration(); + Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st); + pcoh.init_coefficients(2); + pcoh.compute_persistent_cohomology(); + + // Output PD + int max_dim = st.dimension(); + for (int i = 0; i < max_dim; i++) { + std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i); + int num_bars = bars.size(); + 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 (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<double, double>(birth, death)); + if (verbose) std::cout << " [" << birth << ", " << death << "]" << std::endl; + } + } + } + + public: + /** \brief Computes bootstrapped distances distribution. + * + * @param[in] N number of bootstrap iterations. + * + */ + template <typename SimplicialComplex> + 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<int> boot(this->n); + for (int j = 0; j < this->n; j++) { + double u = GetUniform(); + int id = std::floor(u * (this->n)); + boot[j] = id; + Cboot.point_cloud[j] = this->point_cloud[id]; + Cboot.func.emplace(j, this->func[id]); + } + for (int j = 0; j < n; j++) { + std::vector<double> dist(n); + for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]]; + 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 threshold 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 threshold. + * + * @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<double>::lowest(); + int N = PD.size(); + for (int i = 0; i < N; i++) distancemin = std::min(distancemin, 0.5 * (PD[i].second - PD[i].first)); + return 1 - compute_confidence_level_from_distance(distancemin); } // ******************************************************************************************************************* + // Computation of simplices. // ******************************************************************************************************************* public: /** \brief Creates the simplicial complex. * - * @param[in] complex SimplicialComplexForRips to be created. + * @param[in] complex SimplicialComplex to be created. * */ - template <typename SimplicialComplexForRips> - void create_complex(SimplicialComplexForRips& complex) { + template <typename SimplicialComplex> + 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<double>::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: @@ -1036,15 +1211,16 @@ class Cover_complex { } if (type == "Nerve") { - for (std::map<int, std::vector<Cover_t> >::iterator it = cover.begin(); it != cover.end(); it++) + for (std::map<int, std::vector<int> >::iterator it = cover.begin(); it != cover.end(); it++) simplices.push_back(it->second); - std::vector<std::vector<Cover_t> >::iterator it; std::sort(simplices.begin(), simplices.end()); - it = std::unique(simplices.begin(), simplices.end()); + std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); simplices.resize(std::distance(simplices.begin(), it)); } if (type == "GIC") { + 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. @@ -1053,81 +1229,44 @@ class Cover_complex { 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<int, std::vector<Cover_t> >::iterator it = cover.begin(); it != cover.end(); it++) { - int vid = it->first; - std::vector<int> 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<int> 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<int> edge(2); - edge[0] = v1; - edge[1] = v2; - simplices.push_back(edge); + // Loop on all edges. + boost::graph_traits<Graph>::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[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[boost::target(*ei, one_skeleton)]][j]; + if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) { + std::vector<int> edge(2); + edge[0] = std::min(vs, vt); + edge[1] = std::max(vs, vt); + simplices.push_back(edge); + goto afterLoop; + } } } + afterLoop:; } - std::vector<std::vector<Cover_t> >::iterator it; std::sort(simplices.begin(), simplices.end()); - it = std::unique(simplices.begin(), simplices.end()); + std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); simplices.resize(std::distance(simplices.begin(), it)); } else { - // Find IDs of edges to remove - std::vector<int> simplex_to_remove; - int simplex_id = 0; - for (auto simplex : st.complex_simplex_range()) { - if (st.dimension(simplex) == 1) { - std::vector<std::vector<Cover_t> > 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); + // Find edges to keep + Simplex_tree st; + boost::graph_traits<Graph>::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<int> edge(2); + edge[0] = index[boost::source(*ei, one_skeleton)]; + edge[1] = index[boost::target(*ei, one_skeleton)]; + st.insert_simplex_and_subfaces(edge); } - 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); - } + // st.insert_graph(one_skeleton); // Build the Simplex Tree corresponding to the graph st.expansion(maximal_dim); @@ -1136,23 +1275,21 @@ class Cover_complex { simplices.clear(); for (auto simplex : st.complex_simplex_range()) { if (!st.has_children(simplex)) { - std::vector<Cover_t> simplx; + std::vector<int> 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<Cover_t>::iterator it = std::unique(simplx.begin(), simplx.end()); + std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end()); simplx.resize(std::distance(simplx.begin(), it)); simplices.push_back(simplx); } } - std::vector<std::vector<Cover_t> >::iterator it; std::sort(simplices.begin(), simplices.end()); - it = std::unique(simplices.begin(), simplices.end()); + std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end()); simplices.resize(std::distance(simplices.begin(), it)); } } |