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 --- 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 ++++++++++++++++++++++++++------------ 5 files changed, 249 insertions(+), 128 deletions(-) create mode 100644 src/Nerve_GIC/doc/funcGICvisu.pdf (limited to 'src/Nerve_GIC') 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. */ -- cgit v1.2.3