diff options
Diffstat (limited to 'src/Nerve_GIC')
23 files changed, 2866 insertions, 0 deletions
diff --git a/src/Nerve_GIC/doc/COPYRIGHT b/src/Nerve_GIC/doc/COPYRIGHT new file mode 100644 index 00000000..0c36a526 --- /dev/null +++ b/src/Nerve_GIC/doc/COPYRIGHT @@ -0,0 +1,19 @@ +The files of this directory are part of the Gudhi Library. The Gudhi library +(Geometric Understanding in Higher Dimensions) is a generic C++ library for +computational topology. + +Author(s): Mathieu Carrière + +Copyright (C) 2017 INRIA + +This program is free software: you can redistribute it and/or modify it under +the terms of the GNU General Public License as published by the Free Software +Foundation, either version 3 of the License, or (at your option) any later +version. + +This program is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS +FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along with +this program. If not, see <http://www.gnu.org/licenses/>. diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h new file mode 100644 index 00000000..3a6c6c85 --- /dev/null +++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h @@ -0,0 +1,187 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Clément Maria, Pawel Dlotko, Vincent Rouvreau + * + * Copyright (C) 2016 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef DOC_GRAPH_INDUCED_COMPLEX_INTRO_GRAPH_INDUCED_COMPLEX_H_ +#define DOC_GRAPH_INDUCED_COMPLEX_INTRO_GRAPH_INDUCED_COMPLEX_H_ + +namespace Gudhi { + +namespace graph_induced_complex { + +/** \defgroup graph_induced_complex Graph induced complex + * + * \author Mathieu Carrière + * + * @{ + * + * Visualizations of the simplicial complexes require neato, python and firefox!! + * + * \section covers Covers + * + * Nerves and Graph Induced Complexes require a cover C of the input point cloud P, + * that is a set of subsets of P whose union is P itself. + * Very often, this cover is obtained from the preimage of a family of intervals covering + * the image of some scalar-valued function f defined on P. This family is parameterized + * by its resolution, which can be either the number or the length of the intervals, + * and its gain, which is the overlap percentage between consecutive intervals (ordered by their first values). + * + * \section nerves Nerves + * + * \subsection nervedefinition Nerve definition + * + * Assume you are given a cover C of your point cloud P. Then, the Nerve of this cover + * is the simplicial complex that has one k-simplex per k-fold intersection of cover elements. + * See also <a target="_blank" href="https://en.wikipedia.org/wiki/Nerve_of_a_covering"> Wikipedia </a>. + * + * \image html "nerve.png" "Nerve of a double torus" + * + * \subsection nerveexample Example + * + * This example builds the Nerve of a point cloud sampled on a 3D human shape (human.off). + * The cover C comes from the preimages of intervals (10 intervals with gain 0.3) + * covering the height function (coordinate 2), + * which are then refined into their connected components using the triangulation of the .OFF file. + * + * \include Nerve_GIC/Nerve.cpp + * + * When launching: + * + * \code $> ./Nerve ../../../data/points/human.off 2 10 0.3 --v + * \endcode + * + * the program output is: + * + * \include Nerve_GIC/Nerve.txt + * + * The first three lines are requirements for visualization with Kepler-Mapper. + * The fourth line contains the number of vertices nv and edges ne of the Nerve. + * The next nv lines represent the vertices. Each line contains the vertex ID, + * the number of data points it contains, and their average color function value. + * Finally, the next ne lines represent the edges, characterized by the ID of their vertices. + * + * + * \section gic Graph Induced Complexes (GIC) + * + * \subsection gicdefinition GIC definition + * + * Again, assume you are given a cover C of your point cloud P. Moreover, assume + * you are also given a graph G built on top of P. Then, for any clique in G + * whose nodes all belong to different elements of C, the GIC includes a corresponding + * simplex, whose dimension is the number of nodes in the clique minus one. + * See <a target="_blank" href="https://arxiv.org/abs/1304.0662"> this article </a> + * for more details. + * + * \image html "gic_complex.png" "GIC of a point cloud." + * + * \subsection gicexample Example with cover from function + * + * This example builds the GIC of a point cloud sampled on a 3D human shape (human.off). + * The cover C comes from the preimages of intervals (with length 0.075 and gain 0) + * covering the height function (coordinate 2), + * and the graph G comes from a Rips complex built with threshold 0.075. + * Note that if the gain is too big, the number of cliques increases a lot, + * which make the computation time much larger. + * + * \include Nerve_GIC/GIC.cpp + * + * When launching: + * + * \code $> ./GIC ../../../data/points/human.off 0.075 2 0.075 0 --v + * \endcode + * + * the program output is: + * + * \include Nerve_GIC/GIC.txt + * + * \subsection gicexamplevor Example with cover from Voronoï + * + * This example builds the GIC of a point cloud sampled on a 3D human shape (human.off). + * We randomly subsampled 100 points in the point cloud, which act as seeds of + * a geodesic Voronoï diagram. Each cell of the diagram is then an element of C. + * The graph G (used to compute both the geodesics for Voronoï and the GIC) + * comes from the triangulation of the human shape. + * + * \include Nerve_GIC/GICvoronoi.cpp + * + * When launching: + * + * \code $> ./GICvoronoi ../../../data/points/human.off 100 --v + * \endcode + * + * the program output is: + * + * \include Nerve_GIC/GICvoronoi.txt + * + * \subsection mapperdeltadefinition Mapper Delta + * + * If one restricts to the cliques in G whose nodes all belong to preimages of consecutive + * intervals (assuming the cover of the height function is minimal, i.e. no more than + * two intervals can intersect at a time), the GIC is of dimension one, i.e. a graph. + * We call this graph the Mapper Delta, since it is related to the usual Mapper. See + * <a target="_blank" href="https://arxiv.org/abs/1511.05823"> this article </a> for more details. + * + * \subsection mapperdeltaexample Example + * + * Mapper Delta comes with optimal selection for the Rips threshold, + * the resolution and the gain of the function cover. In this example, + * we compute the Mapper Delta of a point cloud sampled on a 3D human shape (human.off), + * where the graph G comes from a Rips complex with optimal threshold, + * and the cover C comes from the preimages of intervals covering the height function (coordinate 2), + * with optimal resolution and gain. Note that optimal threshold, resolution and gain + * also exist for the Nerve of this cover. + * + * \include Nerve_GIC/MapperDeltaCoord.cpp + * + * When launching: + * + * \code $> ./MapperDeltaCoord ../../../data/points/human.off 2 --v + * \endcode + * + * the program output is: + * + * \include MapperDeltaCoord.txt + * + * We also provide an example on a set of 72 pictures taken around the same object (lucky_cat.off). + * The function is now the first eigenfunction given by PCA, whose values + * are written in a file (lucky_cat_PCA1). Threshold, resolution and gain are automatically selected as before. + * + * \include Nerve_GIC/MapperDeltaFunc.cpp + * + * When launching: + * + * \code $> ./MapperDeltaFunc ../../../data/points/COIL_database/lucky_cat.off ../../../data/points/COIL_database/lucky_cat_PCA1 --v + * \endcode + * + * the program output is: + * + * \include MapperDeltaFunc.txt + * + * \copyright GNU General Public License v3. + * \verbatim Contact: gudhi-users@lists.gforge.inria.fr \endverbatim + */ +/** @} */ // end defgroup graph_induced_complex + +} // namespace graph_induced_complex + +} // namespace Gudhi + +#endif // DOC_GRAPH_INDUCED_COMPLEX_INTRO_GRAPH_INDUCED_COMPLEX_H_ diff --git a/src/Nerve_GIC/doc/gic_complex.png b/src/Nerve_GIC/doc/gic_complex.png Binary files differnew file mode 100644 index 00000000..fb4b20ad --- /dev/null +++ b/src/Nerve_GIC/doc/gic_complex.png diff --git a/src/Nerve_GIC/doc/nerve.png b/src/Nerve_GIC/doc/nerve.png Binary files differnew file mode 100644 index 00000000..b66da4a4 --- /dev/null +++ b/src/Nerve_GIC/doc/nerve.png diff --git a/src/Nerve_GIC/example/CMakeLists.txt b/src/Nerve_GIC/example/CMakeLists.txt new file mode 100644 index 00000000..b2c501c3 --- /dev/null +++ b/src/Nerve_GIC/example/CMakeLists.txt @@ -0,0 +1,27 @@ +cmake_minimum_required(VERSION 2.6) +project(Nerve_GIC_examples) + +add_executable ( Nerve Nerve.cpp ) +target_link_libraries(Nerve ${Boost_SYSTEM_LIBRARY}) + +add_executable ( GIC GIC.cpp ) +target_link_libraries(GIC ${Boost_SYSTEM_LIBRARY}) + +add_executable ( MapperDeltaCoord MapperDeltaCoord.cpp ) +target_link_libraries(MapperDeltaCoord ${Boost_SYSTEM_LIBRARY}) + +add_executable ( MapperDeltaFunc MapperDeltaFunc.cpp ) +target_link_libraries(MapperDeltaFunc ${Boost_SYSTEM_LIBRARY}) + +add_executable ( GICvoronoi GICvoronoi.cpp ) +target_link_libraries(MapperDeltaFunc ${Boost_SYSTEM_LIBRARY}) + +file(COPY visu.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +if (TBB_FOUND) + target_link_libraries(Nerve ${TBB_LIBRARIES}) + target_link_libraries(GIC ${TBB_LIBRARIES}) + target_link_libraries(MapperDeltaCoord ${TBB_LIBRARIES}) + target_link_libraries(MapperDeltaFunc ${TBB_LIBRARIES}) + target_link_libraries(GICvoronoi ${TBB_LIBRARIES}) +endif()
\ No newline at end of file diff --git a/src/Nerve_GIC/example/GIC.cpp b/src/Nerve_GIC/example/GIC.cpp new file mode 100644 index 00000000..1889cb33 --- /dev/null +++ b/src/Nerve_GIC/example/GIC.cpp @@ -0,0 +1,65 @@ +#include <gudhi/GIC.h> + +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"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 6) && (argc != 7)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + double threshold = atof(argv[2]); + int coord = atoi(argv[3]); + double resolution = atof(argv[4]); + double gain = atof(argv[5]); + bool verb = 0; if(argc == 7) verb = 1; + + // ---------------------------------------------------------------------------- + // Init of a graph induced complex from an OFF file + // ---------------------------------------------------------------------------- + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + GIC.set_verbose(verb); + + GIC.read_point_cloud(off_file_name); + + GIC.set_color_from_coordinate(coord); + GIC.set_function_from_coordinate(coord); + + GIC.set_graph_from_rips(threshold); + + GIC.set_resolution_double(resolution); GIC.set_gain(gain); + GIC.set_cover_from_function(1); + + GIC.find_GIC_simplices(); + + GIC.plot_with_KeplerMapper(); + + Simplex_tree stree; GIC.create_complex(stree); + + std::streambuf* streambufffer = std::cout.rdbuf(); + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if(verb){ + output_stream << "Graph induced complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on graph induced complex simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << std::endl; + } + } + + return 0; +} diff --git a/src/Nerve_GIC/example/GIC.txt b/src/Nerve_GIC/example/GIC.txt new file mode 100644 index 00000000..79f61e92 --- /dev/null +++ b/src/Nerve_GIC/example/GIC.txt @@ -0,0 +1,89 @@ +Graph induced complex is of dimension 1 - 87 simplices - 44 vertices. +Iterator on graph induced complex simplices +0 +1 +2 +2 1 +3 +3 0 +4 +4 2 +5 +5 3 +6 +6 5 +7 +7 4 +8 +8 7 +9 +9 6 +10 +10 8 +11 +11 9 +12 +12 11 +13 +13 10 +14 +14 12 +15 +15 13 +16 +16 14 +17 +17 15 +18 +18 16 +18 17 +19 +19 18 +20 +21 +22 +22 20 +23 +23 19 +24 +24 21 +25 +25 23 +26 +26 24 +27 +27 22 +28 +28 25 +29 +29 26 +30 +30 27 +31 +31 28 +32 +32 30 +33 +33 29 +34 +34 32 +35 +35 31 +36 +36 33 +37 +37 34 +37 35 +37 36 +38 +38 37 +39 +39 38 +40 +40 39 +41 +41 40 +42 +42 41 +43 +43 42 diff --git a/src/Nerve_GIC/example/GICvoronoi.cpp b/src/Nerve_GIC/example/GICvoronoi.cpp new file mode 100644 index 00000000..9bf3de5e --- /dev/null +++ b/src/Nerve_GIC/example/GICvoronoi.cpp @@ -0,0 +1,60 @@ +#include <gudhi/GIC.h> + +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"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + int m = atoi(argv[2]); + bool verb = 0; if(argc == 4) verb = 1; + + // ---------------------------------------------------------------------------- + // Init of a graph induced complex from an OFF file + // ---------------------------------------------------------------------------- + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + GIC.set_verbose(verb); + + GIC.read_point_cloud(off_file_name); + + GIC.set_color_from_coordinate(); + + GIC.set_graph_from_OFF(off_file_name); + + GIC.set_cover_from_Voronoi(m); + + GIC.find_GIC_simplices(); + + GIC.plot_with_Geomview(); + + Simplex_tree stree; GIC.create_complex(stree); + + std::streambuf* streambufffer = std::cout.rdbuf(); + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if(verb){ + output_stream << "Graph induced complex is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on graph induced complex simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << std::endl; + } + } + + return 0; +} diff --git a/src/Nerve_GIC/example/GICvoronoi.txt b/src/Nerve_GIC/example/GICvoronoi.txt new file mode 100644 index 00000000..34517933 --- /dev/null +++ b/src/Nerve_GIC/example/GICvoronoi.txt @@ -0,0 +1,563 @@ +Graph induced complex is of dimension 2 - 561 simplices - 100 vertices. +Iterator on graph induced complex simplices +0 +1 +2 +3 +3 0 +4 +4 1 +5 +6 +7 +8 +9 +9 1 +10 +11 +11 1 +11 9 +11 9 1 +12 +12 1 +12 11 +12 11 1 +13 +13 2 +14 +14 11 +15 +15 4 +16 +16 2 +16 15 +17 +17 3 +18 +19 +19 8 +20 +20 11 +20 12 +20 12 11 +21 +22 +22 12 +23 +23 22 +24 +24 19 +25 +25 0 +25 3 +25 18 +26 +26 1 +26 4 +26 4 1 +27 +27 26 +28 +28 27 +29 +30 +31 +31 30 +32 +32 30 +32 31 +32 31 30 +33 +33 1 +33 4 +33 4 1 +33 24 +33 29 +34 +34 5 +34 21 +35 +35 2 +35 16 +35 16 2 +36 +36 17 +37 +37 24 +37 29 +38 +38 3 +38 25 +38 25 3 +39 +39 8 +39 25 +39 38 +39 38 25 +40 +40 7 +40 19 +40 24 +40 24 19 +40 37 +40 37 24 +40 39 +41 +41 5 +41 34 +41 34 5 +42 +42 21 +42 41 +43 +43 5 +43 41 +43 41 5 +43 42 +43 42 41 +44 +44 14 +44 27 +45 +45 21 +46 +46 8 +46 19 +46 19 8 +46 39 +46 39 8 +47 +47 0 +47 3 +47 3 0 +47 19 +47 38 +47 38 3 +47 39 +47 39 38 +47 40 +47 40 19 +47 40 39 +47 46 +47 46 19 +47 46 39 +48 +48 37 +48 45 +49 +49 10 +49 37 +49 48 +49 48 37 +50 +50 10 +50 37 +50 45 +50 49 +50 49 10 +50 49 37 +51 +51 45 +51 48 +51 48 45 +51 50 +51 50 45 +52 +52 27 +52 30 +52 31 +52 31 30 +53 +53 10 +53 45 +53 48 +53 48 45 +53 50 +53 50 10 +54 +54 6 +55 +55 54 +56 +56 44 +57 +57 13 +58 +58 24 +58 29 +58 33 +58 33 24 +58 33 29 +58 37 +58 37 24 +58 37 29 +59 +59 7 +59 29 +60 +60 1 +60 12 +60 12 1 +60 22 +60 22 12 +60 26 +60 26 1 +60 27 +60 27 26 +60 44 +60 44 27 +60 56 +60 56 44 +61 +61 5 +61 21 +61 34 +61 34 5 +61 34 21 +61 42 +61 42 21 +61 43 +61 43 5 +61 43 42 +62 +62 6 +62 54 +62 54 6 +62 55 +62 55 54 +63 +63 6 +63 62 +63 62 6 +64 +64 10 +64 48 +64 49 +64 49 10 +64 49 48 +64 53 +64 53 10 +64 53 48 +65 +65 21 +65 45 +65 45 21 +65 50 +65 50 45 +65 53 +65 53 50 +66 +66 14 +66 44 +66 44 14 +66 56 +66 56 44 +67 +67 13 +67 57 +67 57 13 +68 +68 2 +68 13 +68 13 2 +68 35 +68 35 2 +69 +69 21 +69 34 +69 34 21 +69 41 +69 41 34 +69 42 +69 42 21 +69 42 41 +69 65 +69 65 21 +70 +70 11 +70 14 +70 14 11 +70 20 +70 20 11 +70 56 +70 66 +70 66 14 +70 66 56 +71 +71 6 +71 54 +71 54 6 +71 55 +71 55 54 +71 62 +71 62 55 +71 63 +71 63 6 +71 63 62 +72 +72 37 +72 48 +72 48 37 +72 51 +72 51 48 +73 +73 3 +73 17 +73 17 3 +73 36 +73 36 17 +74 +74 7 +74 8 +74 39 +74 39 8 +74 40 +74 40 7 +74 40 39 +74 52 +74 59 +74 59 7 +75 +75 4 +75 15 +75 15 4 +76 +76 15 +76 75 +76 75 15 +77 +77 4 +77 75 +77 75 4 +78 +78 4 +78 26 +78 26 4 +78 27 +78 27 26 +78 28 +78 28 27 +78 75 +78 77 +78 77 4 +78 77 75 +79 +79 9 +79 11 +79 11 9 +79 14 +79 14 11 +79 27 +79 30 +79 32 +79 32 30 +79 44 +79 44 14 +79 44 27 +79 52 +79 52 27 +79 52 30 +80 +80 0 +80 25 +80 25 0 +80 39 +80 39 25 +80 47 +80 47 0 +80 47 39 +81 +81 2 +81 4 +81 15 +81 15 4 +81 16 +81 16 2 +81 16 15 +81 27 +81 29 +81 33 +81 33 4 +81 33 29 +81 52 +81 52 27 +81 59 +81 59 29 +81 74 +81 74 52 +81 74 59 +81 75 +81 78 +81 78 27 +81 78 75 +82 +82 2 +82 75 +82 81 +82 81 2 +82 81 75 +83 +83 0 +83 3 +83 3 0 +83 25 +83 25 0 +83 73 +83 73 3 +84 +84 1 +84 8 +84 19 +84 19 8 +84 24 +84 24 19 +84 33 +84 33 1 +84 33 24 +84 74 +84 74 8 +85 +85 17 +85 36 +85 36 17 +85 73 +85 73 36 +86 +86 6 +86 63 +86 63 6 +87 +87 3 +87 17 +87 17 3 +87 18 +87 25 +87 25 3 +87 25 18 +88 +88 17 +88 18 +88 85 +88 85 17 +88 87 +88 87 17 +88 87 18 +89 +89 63 +89 86 +89 86 63 +90 +90 7 +90 37 +90 40 +90 40 7 +90 40 37 +90 59 +90 59 7 +91 +91 21 +91 45 +91 45 21 +91 53 +91 53 45 +91 65 +91 65 21 +91 65 53 +92 +92 9 +92 55 +92 71 +92 71 55 +93 +93 2 +93 15 +93 16 +93 16 15 +93 35 +93 35 2 +93 35 16 +93 75 +93 76 +93 76 15 +93 76 75 +93 82 +93 82 2 +93 82 75 +94 +94 1 +94 9 +94 9 1 +94 31 +94 32 +94 32 31 +94 52 +94 52 31 +94 55 +94 74 +94 74 52 +94 79 +94 79 9 +94 79 32 +94 84 +94 84 1 +94 84 74 +94 92 +94 92 9 +94 92 55 +95 +95 18 +95 25 +95 25 18 +95 83 +95 83 25 +95 85 +95 88 +95 88 18 +95 88 85 +96 +96 63 +96 86 +96 86 63 +96 89 +96 89 63 +96 89 86 +97 +97 73 +97 83 +97 83 73 +97 85 +97 85 73 +97 95 +97 95 83 +97 95 85 +98 +98 12 +98 20 +98 20 12 +98 22 +98 22 12 +98 23 +98 23 22 +98 56 +98 60 +98 60 22 +98 60 56 +98 70 +98 70 20 +98 70 56 +99 +99 29 +99 37 +99 37 29 +99 50 +99 50 37 +99 51 +99 51 50 +99 59 +99 59 29 +99 72 +99 72 37 +99 72 51 +99 90 +99 90 37 +99 90 59 diff --git a/src/Nerve_GIC/example/MapperDeltaCoord.cpp b/src/Nerve_GIC/example/MapperDeltaCoord.cpp new file mode 100644 index 00000000..9445e988 --- /dev/null +++ b/src/Nerve_GIC/example/MapperDeltaCoord.cpp @@ -0,0 +1,62 @@ +#include <gudhi/GIC.h> + +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"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + int coord = atoi(argv[2]); + bool verb = 0; if(argc == 4) verb = 1; + + // --------------------------------------- + // Init of a Mapper Delta from an OFF file + // --------------------------------------- + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + GIC.set_verbose(verb); + + GIC.read_point_cloud(off_file_name); + + GIC.set_color_from_coordinate(coord); + GIC.set_function_from_coordinate(coord); + + GIC.set_graph_from_automatic_rips(); + + GIC.set_automatic_resolution_for_GICMAP(); GIC.set_gain(); + GIC.set_cover_from_function(1); + + GIC.find_GICMAP_simplices_with_functional_minimal_cover(); + + GIC.plot_with_KeplerMapper(); + + Simplex_tree stree; GIC.create_complex(stree); + + std::streambuf* streambufffer = std::cout.rdbuf(); + std::ostream output_stream(streambufffer); + + // ------------------------------------------ + // Display information about the Mapper Delta + // ------------------------------------------ + + if(verb){ + output_stream << "Mapper Delta is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on Mapper Delta simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << std::endl; + } + } + + return 0; +} diff --git a/src/Nerve_GIC/example/MapperDeltaCoord.txt b/src/Nerve_GIC/example/MapperDeltaCoord.txt new file mode 100644 index 00000000..faf832c7 --- /dev/null +++ b/src/Nerve_GIC/example/MapperDeltaCoord.txt @@ -0,0 +1,143 @@ +Mapper Delta is of dimension 1 - 141 simplices - 71 vertices. +Iterator on Mapper Delta simplices +0 +1 +2 +2 1 +3 +3 0 +4 +4 2 +5 +5 3 +6 +6 4 +7 +7 5 +8 +8 6 +9 +9 7 +10 +10 9 +11 +11 8 +12 +12 11 +13 +13 10 +14 +14 12 +15 +15 13 +16 +16 14 +17 +17 15 +18 +18 17 +19 +19 16 +20 +20 18 +21 +21 19 +22 +22 20 +23 +23 21 +24 +24 22 +25 +25 23 +26 +26 24 +27 +27 25 +28 +28 26 +28 27 +29 +29 28 +30 +30 29 +31 +32 +33 +33 30 +34 +34 31 +35 +35 32 +36 +36 34 +37 +37 33 +38 +38 35 +39 +39 37 +40 +40 38 +41 +41 36 +42 +42 39 +43 +43 40 +44 +44 41 +45 +45 42 +46 +46 44 +47 +47 43 +48 +48 45 +49 +49 46 +50 +50 47 +51 +51 48 +52 +52 49 +53 +53 50 +54 +54 51 +55 +55 52 +56 +56 53 +57 +57 54 +58 +58 55 +59 +59 56 +60 +60 57 +60 58 +60 59 +61 +61 60 +62 +62 61 +63 +63 62 +64 +64 63 +65 +65 64 +66 +66 65 +67 +67 66 +68 +68 67 +69 +69 68 +70 +70 69 diff --git a/src/Nerve_GIC/example/MapperDeltaFunc.cpp b/src/Nerve_GIC/example/MapperDeltaFunc.cpp new file mode 100644 index 00000000..448323b1 --- /dev/null +++ b/src/Nerve_GIC/example/MapperDeltaFunc.cpp @@ -0,0 +1,62 @@ +#include <gudhi/GIC.h> + +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"; + exit(-1); // ----- >> +} + +int main(int argc, char **argv) { + if ((argc != 3) && (argc != 4)) usage(argc, (argv[0] - 1)); + + std::string off_file_name(argv[1]); + std::string func_file_name = argv[2]; + bool verb = 0; if(argc == 4) verb = 1; + + // --------------------------------------- + // Init of a Mapper Delta from an OFF file + // --------------------------------------- + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + GIC.set_verbose(verb); + + GIC.read_point_cloud(off_file_name); + + GIC.set_color_from_file(func_file_name); + GIC.set_function_from_file(func_file_name); + + GIC.set_graph_from_automatic_rips(); + + GIC.set_automatic_resolution_for_GICMAP(); GIC.set_gain(); + GIC.set_cover_from_function(1); + + GIC.find_GICMAP_simplices_with_functional_minimal_cover(); + + GIC.plot_with_KeplerMapper(); + + Simplex_tree stree; GIC.create_complex(stree); + + std::streambuf* streambufffer = std::cout.rdbuf(); + std::ostream output_stream(streambufffer); + + // ------------------------------------------ + // Display information about the Mapper Delta + // ------------------------------------------ + + if(verb){ + output_stream << "Mapper Delta is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on Mapper Delta simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << std::endl; + } + } + + return 0; +} diff --git a/src/Nerve_GIC/example/MapperDeltaFunc.txt b/src/Nerve_GIC/example/MapperDeltaFunc.txt new file mode 100644 index 00000000..d13d0695 --- /dev/null +++ b/src/Nerve_GIC/example/MapperDeltaFunc.txt @@ -0,0 +1,10 @@ +Mapper Delta is of dimension 1 - 8 simplices - 4 vertices. +Iterator on Mapper Delta simplices +0 +1 +1 0 +2 +2 0 +3 +3 1 +3 2 diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/example/Nerve.cpp new file mode 100644 index 00000000..84f74625 --- /dev/null +++ b/src/Nerve_GIC/example/Nerve.cpp @@ -0,0 +1,64 @@ +#include <gudhi/GIC.h> + +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/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] - 1)); + + std::string off_file_name(argv[1]); + int coord = atoi(argv[2]); + int resolution = atoi(argv[3]); + double gain = atof(argv[4]); + bool verb = 0; if(argc == 6) verb = 1; + + // -------------------------------- + // Init of a Nerve from an OFF file + // -------------------------------- + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + GIC.set_verbose(verb); + + GIC.read_point_cloud(off_file_name); + + GIC.set_color_from_coordinate(coord); + GIC.set_function_from_coordinate(coord); + + GIC.set_graph_from_OFF(off_file_name); + + GIC.set_resolution_int(resolution); GIC.set_gain(gain); + GIC.set_cover_from_function(0); + + GIC.find_Nerve_simplices(); + + GIC.plot_with_KeplerMapper(); + + Simplex_tree stree; GIC.create_complex(stree); + + std::streambuf* streambufffer = std::cout.rdbuf(); + std::ostream output_stream(streambufffer); + + // ---------------------------------------------------------------------------- + // Display information about the graph induced complex + // ---------------------------------------------------------------------------- + + if(verb){ + output_stream << "Nerve is of dimension " << stree.dimension() << + " - " << stree.num_simplices() << " simplices - " << + stree.num_vertices() << " vertices." << std::endl; + + output_stream << "Iterator on Nerve simplices" << std::endl; + for (auto f_simplex : stree.filtration_simplex_range()) { + for (auto vertex : stree.simplex_vertex_range(f_simplex)) { + output_stream << vertex << " "; + } + output_stream << std::endl; + } + } + + return 0; +} diff --git a/src/Nerve_GIC/example/Nerve.txt b/src/Nerve_GIC/example/Nerve.txt new file mode 100644 index 00000000..2a861c5f --- /dev/null +++ b/src/Nerve_GIC/example/Nerve.txt @@ -0,0 +1,43 @@ +Nerve is of dimension 1 - 41 simplices - 21 vertices. +Iterator on Nerve simplices +0 +1 +2 +2 0 +3 +3 1 +4 +4 3 +5 +5 2 +6 +6 4 +7 +7 5 +8 +9 +9 6 +10 +10 7 +11 +12 +12 8 +13 +13 9 +13 10 +14 +14 11 +15 +15 13 +16 +16 12 +17 +17 14 +18 +18 15 +18 16 +18 17 +19 +19 18 +20 +20 19 diff --git a/src/Nerve_GIC/example/km.py b/src/Nerve_GIC/example/km.py new file mode 100755 index 00000000..a1a08156 --- /dev/null +++ b/src/Nerve_GIC/example/km.py @@ -0,0 +1,390 @@ +from __future__ import division +import numpy as np +from collections import defaultdict +import json +import itertools +from sklearn import cluster, preprocessing, manifold +from datetime import datetime +import sys + +class KeplerMapper(object): + # With this class you can build topological networks from (high-dimensional) data. + # + # 1) Fit a projection/lens/function to a dataset and transform it. + # For instance "mean_of_row(x) for x in X" + # 2) Map this projection with overlapping intervals/hypercubes. + # Cluster the points inside the interval + # (Note: we cluster on the inverse image/original data to lessen projection loss). + # If two clusters/nodes have the same members (due to the overlap), then: + # connect these with an edge. + # 3) Visualize the network using HTML and D3.js. + # + # functions + # --------- + # fit_transform: Create a projection (lens) from a dataset + # map: Apply Mapper algorithm on this projection and build a simplicial complex + # visualize: Turns the complex dictionary into a HTML/D3.js visualization + + def __init__(self, verbose=2): + self.verbose = verbose + + self.chunk_dist = [] + self.overlap_dist = [] + self.d = [] + self.nr_cubes = 0 + self.overlap_perc = 0 + self.clusterer = False + + def fit_transform(self, X, projection="sum", scaler=preprocessing.MinMaxScaler()): + # Creates the projection/lens from X. + # + # Input: X. Input features as a numpy array. + # Output: projected_X. original data transformed to a projection (lens). + # + # parameters + # ---------- + # projection: Projection parameter is either a string, + # a scikit class with fit_transform, like manifold.TSNE(), + # or a list of dimension indices. + # scaler: if None, do no scaling, else apply scaling to the projection + # Default: Min-Max scaling + + self.scaler = scaler + self.projection = str(projection) + + # Detect if projection is a class (for scikit-learn) + if str(type(projection))[1:6] == "class": #TODO: de-ugly-fy + reducer = projection + if self.verbose > 0: + try: + projection.set_params(**{"verbose":self.verbose}) + except: + pass + print("\n..Projecting data using: \n\t%s\n"%str(projection)) + X = reducer.fit_transform(X) + + # Detect if projection is a string (for standard functions) + if isinstance(projection, str): + if self.verbose > 0: + print("\n..Projecting data using: %s"%(projection)) + # Stats lenses + if projection == "sum": # sum of row + X = np.sum(X, axis=1).reshape((X.shape[0],1)) + if projection == "mean": # mean of row + X = np.mean(X, axis=1).reshape((X.shape[0],1)) + if projection == "median": # mean of row + X = np.median(X, axis=1).reshape((X.shape[0],1)) + if projection == "max": # max of row + X = np.max(X, axis=1).reshape((X.shape[0],1)) + if projection == "min": # min of row + X = np.min(X, axis=1).reshape((X.shape[0],1)) + if projection == "std": # std of row + X = np.std(X, axis=1).reshape((X.shape[0],1)) + + if projection == "dist_mean": # Distance of x to mean of X + X_mean = np.mean(X, axis=0) + X = np.sum(np.sqrt((X - X_mean)**2), axis=1).reshape((X.shape[0],1)) + + # Detect if projection is a list (with dimension indices) + if isinstance(projection, list): + if self.verbose > 0: + print("\n..Projecting data using: %s"%(str(projection))) + X = X[:,np.array(projection)] + + # Scaling + if scaler is not None: + if self.verbose > 0: + print("\n..Scaling with: %s\n"%str(scaler)) + X = scaler.fit_transform(X) + + return X + + def map(self, projected_X, inverse_X=None, clusterer=cluster.DBSCAN(eps=0.5,min_samples=3), nr_cubes=10, overlap_perc=0.1): + # This maps the data to a simplicial complex. Returns a dictionary with nodes and links. + # + # Input: projected_X. A Numpy array with the projection/lens. + # Output: complex. A dictionary with "nodes", "links" and "meta information" + # + # parameters + # ---------- + # projected_X projected_X. A Numpy array with the projection/lens. Required. + # inverse_X Numpy array or None. If None then the projection itself is used for clustering. + # clusterer Scikit-learn API compatible clustering algorithm. Default: DBSCAN + # nr_cubes Int. The number of intervals/hypercubes to create. + # overlap_perc Float. The percentage of overlap "between" the intervals/hypercubes. + + start = datetime.now() + + # Helper function + def cube_coordinates_all(nr_cubes, nr_dimensions): + # Helper function to get origin coordinates for our intervals/hypercubes + # Useful for looping no matter the number of cubes or dimensions + # Example: if there are 4 cubes per dimension and 3 dimensions + # return the bottom left (origin) coordinates of 64 hypercubes, + # as a sorted list of Numpy arrays + # TODO: elegance-ify... + l = [] + for x in range(nr_cubes): + l += [x] * nr_dimensions + return [np.array(list(f)) for f in sorted(set(itertools.permutations(l,nr_dimensions)))] + + nodes = defaultdict(list) + links = defaultdict(list) + complex = {} + self.nr_cubes = nr_cubes + self.clusterer = clusterer + self.overlap_perc = overlap_perc + + if self.verbose > 0: + print("Mapping on data shaped %s using dimensions\n"%(str(projected_X.shape))) + + # If inverse image is not provided, we use the projection as the inverse image (suffer projection loss) + if inverse_X is None: + inverse_X = projected_X + + # We chop up the min-max column ranges into 'nr_cubes' parts + self.chunk_dist = (np.max(projected_X, axis=0) - np.min(projected_X, axis=0))/nr_cubes + + # We calculate the overlapping windows distance + self.overlap_dist = self.overlap_perc * self.chunk_dist + + # We find our starting point + self.d = np.min(projected_X, axis=0) + + # Use a dimension index array on the projected X + # (For now this uses the entire dimensionality, but we keep for experimentation) + di = np.array([x for x in range(projected_X.shape[1])]) + + # Prefix'ing the data with ID's + ids = np.array([x for x in range(projected_X.shape[0])]) + projected_X = np.c_[ids,projected_X] + inverse_X = np.c_[ids,inverse_X] + + # Subdivide the projected data X in intervals/hypercubes with overlap + if self.verbose > 0: + total_cubes = len(cube_coordinates_all(nr_cubes,projected_X.shape[1])) + print("Creating %s hypercubes."%total_cubes) + + for i, coor in enumerate(cube_coordinates_all(nr_cubes,di.shape[0])): + # Slice the hypercube + hypercube = projected_X[ np.invert(np.any((projected_X[:,di+1] >= self.d[di] + (coor * self.chunk_dist[di])) & + (projected_X[:,di+1] < self.d[di] + (coor * self.chunk_dist[di]) + self.chunk_dist[di] + self.overlap_dist[di]) == False, axis=1 )) ] + + if self.verbose > 1: + print("There are %s points in cube_%s / %s with starting range %s"% + (hypercube.shape[0],i,total_cubes,self.d[di] + (coor * self.chunk_dist[di]))) + + # If at least one sample inside the hypercube + if hypercube.shape[0] > 0: + # Cluster the data point(s) in the cube, skipping the id-column + # Note that we apply clustering on the inverse image (original data samples) that fall inside the cube. + inverse_x = inverse_X[[int(nn) for nn in hypercube[:,0]]] + + clusterer.fit(inverse_x[:,1:]) + + if self.verbose > 1: + print("Found %s clusters in cube_%s\n"%(np.unique(clusterer.labels_[clusterer.labels_ > -1]).shape[0],i)) + + #Now for every (sample id in cube, predicted cluster label) + for a in np.c_[hypercube[:,0],clusterer.labels_]: + if a[1] != -1: #if not predicted as noise + cluster_id = str(coor[0])+"_"+str(i)+"_"+str(a[1])+"_"+str(coor)+"_"+str(self.d[di] + (coor * self.chunk_dist[di])) # TODO: de-rudimentary-ify + nodes[cluster_id].append( int(a[0]) ) # Append the member id's as integers + else: + if self.verbose > 1: + print("Cube_%s is empty.\n"%(i)) + + # Create links when clusters from different hypercubes have members with the same sample id. + candidates = itertools.combinations(nodes.keys(),2) + for candidate in candidates: + # if there are non-unique members in the union + if len(nodes[candidate[0]]+nodes[candidate[1]]) != len(set(nodes[candidate[0]]+nodes[candidate[1]])): + links[candidate[0]].append( candidate[1] ) + + # Reporting + if self.verbose > 0: + nr_links = 0 + for k in links: + nr_links += len(links[k]) + print("\ncreated %s edges and %s nodes in %s."%(nr_links,len(nodes),str(datetime.now()-start))) + + complex["nodes"] = nodes + complex["links"] = links + complex["meta"] = self.projection + + return complex + + def visualize(self, complex, color_function="", path_html="mapper_visualization_output.html", title="My Data", + graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=None, width_html=0, + height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=0,gain=0,minimum=0,maximum=0): + # Turns the dictionary 'complex' in a html file with d3.js + # + # Input: complex. Dictionary (output from calling .map()) + # Output: a HTML page saved as a file in 'path_html'. + # + # parameters + # ---------- + # color_function string. Not fully implemented. Default: "" (distance to origin) + # path_html file path as string. Where to save the HTML page. + # title string. HTML page document title and first heading. + # graph_link_distance int. Edge length. + # graph_gravity float. "Gravity" to center of layout. + # graph_charge int. charge between nodes. + # custom_tooltips None or Numpy Array. You could use "y"-label array for this. + # width_html int. Width of canvas. Default: 0 (full width) + # height_html int. Height of canvas. Default: 0 (full height) + # show_tooltips bool. default:True + # show_title bool. default:True + # show_meta bool. default:True + + # Format JSON for D3 graph + json_s = {} + json_s["nodes"] = [] + json_s["links"] = [] + k2e = {} # a key to incremental int dict, used for id's when linking + + for e, k in enumerate(complex["nodes"]): + # Tooltip and node color formatting, TODO: de-mess-ify + if custom_tooltips is not None: + tooltip_s = "<h2>Cluster %s</h2>"%k + " ".join(str(custom_tooltips[complex["nodes"][k][0]]).split(" ")) + if maximum == minimum: + tooltip_i = 0 + else: + tooltip_i = int(30*(custom_tooltips[complex["nodes"][k][0]]-minimum)/(maximum-minimum)) + json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(complex["nodes"][k][2])), "color": tooltip_i}) + else: + tooltip_s = "<h2>Cluster %s</h2>Contains %s members."%(k,len(complex["nodes"][k])) + json_s["nodes"].append({"name": str(k), "tooltip": tooltip_s, "group": 2 * int(np.log(len(complex["nodes"][k]))), "color": str(k.split("_")[0])}) + k2e[k] = e + for k in complex["links"]: + for link in complex["links"][k]: + json_s["links"].append({"source": k2e[k], "target":k2e[link],"value":1}) + + # Width and height of graph in HTML output + if width_html == 0: + width_css = "100%" + width_js = 'document.getElementById("holder").offsetWidth-20' + else: + width_css = "%spx" % width_html + width_js = "%s" % width_html + if height_html == 0: + height_css = "100%" + height_js = 'document.getElementById("holder").offsetHeight-20' + else: + height_css = "%spx" % height_html + height_js = "%s" % height_html + + # Whether to show certain UI elements or not + if show_tooltips == False: + tooltips_display = "display: none;" + else: + tooltips_display = "" + + if show_meta == False: + meta_display = "display: none;" + else: + meta_display = "" + + if show_title == False: + title_display = "display: none;" + else: + title_display = "" + + with open(path_html,"wb") as outfile: + html = """<!DOCTYPE html> + <meta charset="utf-8"> + <meta name="generator" content="KeplerMapper"> + <title>%s | KeplerMapper</title> + <link href='https://fonts.googleapis.com/css?family=Roboto:700,300' rel='stylesheet' type='text/css'> + <style> + * {margin: 0; padding: 0;} + html { height: 100%%;} + body {background: #111; height: 100%%; font: 100 16px Roboto, Sans-serif;} + .link { stroke: #999; stroke-opacity: .333; } + .divs div { border-radius: 50%%; background: red; position: absolute; } + .divs { position: absolute; top: 0; left: 0; } + #holder { position: relative; width: %s; height: %s; background: #111; display: block;} + h1 { %s padding: 20px; color: #fafafa; text-shadow: 0px 1px #000,0px -1px #000; position: absolute; font: 300 30px Roboto, Sans-serif;} + h2 { text-shadow: 0px 1px #000,0px -1px #000; font: 700 16px Roboto, Sans-serif;} + .meta { position: absolute; opacity: 0.9; width: 220px; top: 80px; left: 20px; display: block; %s background: #000; line-height: 25px; color: #fafafa; border: 20px solid #000; font: 100 16px Roboto, Sans-serif;} + div.tooltip { position: absolute; width: 380px; display: block; %s padding: 20px; background: #000; border: 0px; border-radius: 3px; pointer-events: none; z-index: 999; color: #FAFAFA;} + } + </style> + <body> + <div id="holder"> + <h1>%s</h1> + <p class="meta"> + <b>Lens</b><br>%s<br><br> + <b>Length of intervals</b><br>%s<br><br> + <b>Overlap percentage</b><br>%s%%<br><br> + <b>Color Function</b><br>%s + </p> + </div> + <script src="https://cdnjs.cloudflare.com/ajax/libs/d3/3.5.5/d3.min.js"></script> + <script> + var width = %s, + height = %s; + var color = d3.scale.ordinal() + .domain(["0","1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13","14","15","16","17","18","19","20","21","22","23","24","25","26","27","28","29","30"]) + .range(["#FF0000","#FF1400","#FF2800","#FF3c00","#FF5000","#FF6400","#FF7800","#FF8c00","#FFa000","#FFb400","#FFc800","#FFdc00","#FFf000","#fdff00","#b0ff00","#65ff00","#17ff00","#00ff36","#00ff83","#00ffd0","#00e4ff","#00c4ff","#00a4ff","#00a4ff","#0084ff","#0064ff","#0044ff","#0022ff","#0002ff","#0100ff","#0300ff","#0500ff"]); + var force = d3.layout.force() + .charge(%s) + .linkDistance(%s) + .gravity(%s) + .size([width, height]); + var svg = d3.select("#holder").append("svg") + .attr("width", width) + .attr("height", height); + + var div = d3.select("#holder").append("div") + .attr("class", "tooltip") + .style("opacity", 0.0); + + var divs = d3.select('#holder').append('div') + .attr('class', 'divs') + .attr('style', function(d) { return 'overflow: hidden; width: ' + width + 'px; height: ' + height + 'px;'; }); + + graph = %s; + force + .nodes(graph.nodes) + .links(graph.links) + .start(); + var link = svg.selectAll(".link") + .data(graph.links) + .enter().append("line") + .attr("class", "link") + .style("stroke-width", function(d) { return Math.sqrt(d.value); }); + var node = divs.selectAll('div') + .data(graph.nodes) + .enter().append('div') + .on("mouseover", function(d) { + div.transition() + .duration(200) + .style("opacity", .9); + div .html(d.tooltip + "<br/>") + .style("left", (d3.event.pageX + 100) + "px") + .style("top", (d3.event.pageY - 28) + "px"); + }) + .on("mouseout", function(d) { + div.transition() + .duration(500) + .style("opacity", 0); + }) + .call(force.drag); + + node.append("title") + .text(function(d) { return d.name; }); + force.on("tick", function() { + link.attr("x1", function(d) { return d.source.x; }) + .attr("y1", function(d) { return d.source.y; }) + .attr("x2", function(d) { return d.target.x; }) + .attr("y2", function(d) { return d.target.y; }); + node.attr("cx", function(d) { return d.x; }) + .attr("cy", function(d) { return d.y; }) + .attr('style', function(d) { return 'width: ' + (d.group * 2) + 'px; height: ' + (d.group * 2) + 'px; ' + 'left: '+(d.x-(d.group))+'px; ' + 'top: '+(d.y-(d.group))+'px; background: '+color(d.color)+'; box-shadow: 0px 0px 3px #111; box-shadow: 0px 0px 33px '+color(d.color)+', inset 0px 0px 5px rgba(0, 0, 0, 0.2);'}) + ; + }); + </script>"""%(title,width_css, height_css, title_display, meta_display, tooltips_display, title,complex["meta"],res,gain*100,color_function,width_js,height_js,graph_charge,graph_link_distance,graph_gravity,json.dumps(json_s)) + outfile.write(html.encode("utf-8")) + if self.verbose > 0: + print("\nWrote d3.js graph to '%s'"%path_html)
\ No newline at end of file diff --git a/src/Nerve_GIC/example/visu.py b/src/Nerve_GIC/example/visu.py new file mode 100755 index 00000000..06d22abd --- /dev/null +++ b/src/Nerve_GIC/example/visu.py @@ -0,0 +1,44 @@ +import km +import numpy as np +from collections import defaultdict + +network = {} +mapper = km.KeplerMapper(verbose=0) +data = np.zeros((3,3)) +projected_data = mapper.fit_transform( data, projection="sum", scaler=None ) + +f = open('SC.txt','r') +nodes = defaultdict(list) +links = defaultdict(list) +custom = defaultdict(list) + +dat = f.readline() +lens = f.readline() +color = f.readline(); +param = [float(i) for i in f.readline().split(" ")] + +nums = [int(i) for i in f.readline().split(" ")] +num_nodes = nums[0] +num_edges = nums[1] + +for i in range(0,num_nodes): + point = [float(j) for j in f.readline().split(" ")] + nodes[ str(int(point[0])) ] = [ int(point[0]), point[1], int(point[2]) ] + links[ str(int(point[0])) ] = [] + custom[ int(point[0]) ] = point[1] + +m = min([custom[i] for i in range(0,num_nodes)]) +M = max([custom[i] for i in range(0,num_nodes)]) + +for i in range(0,num_edges): + edge = [int(j) for j in f.readline().split(" ")] + links[ str(edge[0]) ].append( str(edge[1]) ) + links[ str(edge[1]) ].append( str(edge[0]) ) + +network["nodes"] = nodes +network["links"] = links +network["meta"] = lens + +mapper.visualize(network, color_function = color, path_html="SC.html", title=dat, +graph_link_distance=30, graph_gravity=0.1, graph_charge=-120, custom_tooltips=custom, width_html=0, +height_html=0, show_tooltips=True, show_title=True, show_meta=True, res=param[0],gain=param[1], minimum=m,maximum=M)
\ No newline at end of file diff --git a/src/Nerve_GIC/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h new file mode 100644 index 00000000..72e08513 --- /dev/null +++ b/src/Nerve_GIC/include/gudhi/GIC.h @@ -0,0 +1,911 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author: Mathieu Carriere + * + * Copyright (C) 2017 INRIA + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#ifndef GIC_H_ +#define GIC_H_ + +#include <gudhi/Debug_utils.h> +#include <gudhi/graph_simplicial_complex.h> +#include <gudhi/reader_utils.h> +#include <gudhi/Simplex_tree.h> +#include <gudhi/Rips_complex.h> +#include <gudhi/Points_off_io.h> +#include <gudhi/distance_functions.h> +#include <gudhi/Kd_tree_search.h> +#include <CGAL/Epick_d.h> + +#include <boost/graph/adjacency_list.hpp> + +#include <iostream> +#include <vector> +#include <map> +#include <string> +#include <limits> // for numeric_limits +#include <utility> // for pair<> +#include <functional> // for greater<> +#include <stdexcept> +#include <initializer_list> +#include <algorithm> // for std::max +#include <cstdint> // for std::uint32_t +#include <random> +#include <cassert> + +#define CONSTANT 10 +#define ETA 0.001 +#define MASK 0 + +using Simplex_tree = Gudhi::Simplex_tree<>; +using Filtration_value = Simplex_tree::Filtration_value; +using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>; +using Point = std::vector<float>; + +std::map<int, double> func; +std::map<int, double> func_color; + +namespace Gudhi { + +namespace graph_induced_complex { + + +/** + * \class Graph_induced_complex + * \brief Graph induced complex data structure. + * + * \ingroup graph_induced_complex + * + * \details + * The data structure is a simplicial complex, representing a + * Graph Induced simplicial Complex (GIC) or a Nerve, + * and whose simplices are computed with a cover C of a point + * cloud P, which often comes from the preimages of intervals + * covering the image of a function f defined on P. + * These intervals are parameterized by their resolution + * (either their length or their number) + * and their gain (percentage of overlap). + * To compute a GIC, one also needs a graph G built on top of P, + * whose cliques with vertices belonging to different elements of C + * correspond to the simplices of the GIC. + * + */ + +class Graph_induced_complex { + + private: + bool verbose; // whether to display information. + std::vector<Point> point_cloud; + 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; + 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> 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; + int resolution_int; + double resolution_double; + double gain; + std::vector<int> voronoi_subsamples; + std::string cover_name; + std::string point_cloud_name; + std::string color_name; + + // Simplex comparator + private: + bool simplex_comp(const std::vector<Cover_t>& s1, const std::vector<Cover_t>& s2){ + if(s1.size() == s2.size()){ + return s1[0] < s2[0]; + } + else return s1.size() < s2.size(); + } + + // Point comparator + private: + static bool functional_comp(const int& a, const int& b){ + if(func[a] == func[b]) return a < b; + else return func[a] < func[b]; + } + + // DFS + private: + void dfs(std::map<int,std::vector<int> >& G, const 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); + } + + // Find random number in [0,1]. + double GetUniform(){ + static std::default_random_engine re; + static std::uniform_real_distribution<double> Dist(0,1); + return Dist(re); + } + + // Subsample points. + void SampleWithoutReplacement(int populationSize, int sampleSize, std::vector<int> & samples){ + int& n = sampleSize; int& N = populationSize; + int t = 0; int m = 0; double u; + while (m < n){ + u = GetUniform(); + if ( (N - t)*u >= n - m ) + t++; + else{samples[m] = t; t++; m++;} + } + } + + public: + void set_verbose(bool verb = 0){verbose = verb;} + + public: + void read_point_cloud(const std::string& off_file_name){ + Gudhi::Points_off_reader<Point> off_reader = Points_off_reader<Point>(off_file_name); + point_cloud = off_reader.get_point_cloud(); + point_cloud_name = off_file_name; + n = point_cloud.size(); + data_dimension = point_cloud[0].size(); + } + + // ******************************************************************************************************************* + // Graphs. + // ******************************************************************************************************************* + + public: // Set graph from file. + /** \brief Creates the graph G from a file containing the edges. + * + * @param[in] graph_file_name name of the input graph file. + * + */ + void set_graph_from_file(const std::string& graph_file_name){ + int neighb; int vid; std::ifstream input(graph_file_name); std::string line; std::vector<int> edge(2); int n = 0; + while(std::getline(input,line)){ + std::stringstream stream(line); stream >> vid; edge[0] = vid; + while(stream >> neighb){edge[1] = neighb; st.insert_simplex_and_subfaces(edge);} + n++; + } + + std::vector<int> empty; + for(int i = 0; i < n; i++) adjacency_matrix.insert(std::pair<int,std::vector<int> >(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]); + } + } + } + + public: // Set graph from OFF file. + /** \brief Creates the graph G from the triangulation given by an .OFF file. + * + * @param[in] off_file_name name of the input .OFF file. + * + */ + void set_graph_from_OFF(const std::string& off_file_name){ + int numedges, numfaces, i; std::vector<int> edge(2); double x; int num; std::vector<int> simplex; + std::ifstream input(off_file_name); std::string line; getline(input, line); + input >> n; input >> numfaces; input >> numedges; + i = 0; while(i < n){input >> x; input >> x; input >> x; i++;} + i = 0; while(i < numfaces){ + simplex.clear(); input >> num; + for(int j = 0; j < num; j++){int k; input >> k; simplex.push_back(k);} + for(int j = 0; j < num; j++){ + for(int k = j+1; k < num; k++){ + edge[0] = simplex[j]; edge[1] = simplex[k]; + st.insert_simplex_and_subfaces(edge); + } + } + i++; + } + + std::vector<int> empty; + for(int i = 0; i < n; i++) adjacency_matrix.insert(std::pair<int,std::vector<int> >(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]); + } + } + } + + public: // Set graph from Rips complex. + /** \brief Creates the graph G from a Rips complex. + * + * @param[in] threshold threshold value for the Rips complex. + * @param[in] off_file_name name of the input .OFF file. + * + */ + void set_graph_from_rips(const double& threshold){ + Rips_complex rips_complex_from_points(point_cloud, threshold, Euclidean_distance()); + rips_complex_from_points.create_complex(st, 1); + + std::vector<int> empty; + for(int i = 0; i < n; i++) adjacency_matrix.insert(std::pair<int,std::vector<int> >(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]); + } + } + + } + + public: // Pairwise distances. + /** \brief Computes all pairwise distances (Euclidean norm). + */ + void compute_pairwise_distances(){ + + double d; std::vector<double> dumb(n); for(int i = 0; i < n; i++) distances.push_back(dumb); + char distance[100]; sprintf(distance,"%s_dist",(char*) point_cloud_name.c_str()); + std::ifstream input(distance, std::ios::out | std::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; + } + } + input.close(); + } + + else{ + 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++){ + if( (int) floor( 100*(i*1.0+1)/(n*1.0) ) %10 == 0 ) + if(verbose) std::cout << "\r" << floor( 100*(i*1.0+1)/(n*1.0) ) << "%" << std::flush; + for (int j = i; j < n; j++){ + double dis = 0; for(int k = 0; k < data_dimension; k++) + dis += pow(point_cloud[i][k]-point_cloud[j][k],2); + dis = std::sqrt(dis); distances[i][j] = dis; distances[j][i] = dis; + output.write((char*) &dis,8); + } + } + output.close(); if(verbose) std::cout << std::endl; + } + + } + + public: // Automatic tuning of Rips complex. + /** \brief Creates the graph G from a Rips complex whose threshold value is automatically tuned with subsampling. + * + * @param[in] off_file_name name of the input .OFF file. + * @param[in] N number of subsampling iteration (default value 100). + * + */ + void set_graph_from_automatic_rips(int N = 100){ + + int m = floor(n/pow(log(n)/log(CONSTANT),1+ETA)); m = std::min(m,n-1); + std::vector<int> samples(m); double delta = 0; + + if(verbose) std::cout << n << " points in R^" << data_dimension << std::endl; + if(verbose) std::cout << "Subsampling " << m << " points" << std::endl; + + if(distances.size() == 0) compute_pairwise_distances(); + + //#pragma omp parallel for + for (int i = 0; i < N; i++){ + + SampleWithoutReplacement(n,m,samples); + double hausdorff_dist = 0; + for (int j = 0; j < n; j++){ + double mj = distances[j][samples[0]]; for (int k = 1; k < m; k++) mj = std::min(mj, distances[j][samples[k]]); + hausdorff_dist = std::max(hausdorff_dist, mj); + } + delta += hausdorff_dist/N; + + } + + if(verbose) std::cout << "delta = " << delta << std::endl; + Rips_complex rips_complex_from_points(point_cloud, delta, Euclidean_distance()); + rips_complex_from_points.create_complex(st, 1); + + std::vector<int> empty; + for(int i = 0; i < n; i++) adjacency_matrix.insert(std::pair<int,std::vector<int> >(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]); + } + } + + } + + + // ******************************************************************************************************************* + // Functions. + // ******************************************************************************************************************* + + public: // Set function from file. + /** \brief Creates the function f from a file containing the function values. + * + * @param[in] func_file_name name of the input function file. + * + */ + 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.insert(std::pair<int,double>(vertex_id, f)); vertex_id++; + } + cover_name = func_file_name; + } + + public: // Set function from kth coordinate + /** \brief Creates the function f from the k-th coordinate of the point cloud P. + * + * @param[in] k coordinate to use (start at 0). + * + */ + void set_function_from_coordinate(const int& k){ + for(int i = 0; i < n; i++) func.insert(std::pair<int,double>(i,point_cloud[i][k])); + char coordinate[100]; sprintf(coordinate, "coordinate %d", k); + cover_name = coordinate; + } + + public: // Set function from vector. + /** \brief Creates the function f from a vector stored in memory. + * + * @param[in] function input vector of values. + * + */ + void set_function_from_vector(const std::vector<double>& function){ + for(int i = 0; i < n; i++) func.insert(std::pair<int,double>(i, function[i])); + } + + // ******************************************************************************************************************* + // Covers. + // ******************************************************************************************************************* + + public: // Set cover from file. + /** \brief Creates the cover C from a file containing the cover elements of each point (the order has to be the same + * as in the input file!). + * + * @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<Cover_t> cov_elts, cov_number; + std::ifstream input(cover_file_name); std::string line; + while(std::getline(input,line)){ + cov_elts.clear(); std::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[vertex_id] = cov_elts; vertex_id++; + } + std::vector<Cover_t>::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)); 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; + } + + public: // Set cover from Voronoi + /** \brief Creates the cover C from the Voronoï cells of a subsampling of the point cloud. + * + * @param[in] m number of points in the subsample. + * + */ + void set_cover_from_Voronoi(const int& m){ + + voronoi_subsamples.resize(m); SampleWithoutReplacement(n,m,voronoi_subsamples); + if(distances.size() == 0) compute_pairwise_distances(); + 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 + 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(); + } + + 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; + } + } + + for(int i = 0; i < n; 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; cover_name = "Voronoi"; + + } + + public: // Automatic tuning of resolution for Mapper Delta. + /** \brief Computes the optimal length of intervals for a Mapper Delta. + */ + void set_automatic_resolution_for_GICMAP(){ + double reso = 0; + 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]])); + } + } + if(verbose) std::cout << "resolution = " << reso << std::endl; + resolution_double = reso; + } + + public: + /** \brief Sets a length of intervals from a value stored in memory. + * + * @param[in] reso length of intervals. + * + */ + void set_resolution_double(double reso){resolution_double = reso;} + /** \brief Sets a number of intervals from a value stored in memory. + * + * @param[in] reso number of intervals. + * + */ + void set_resolution_int(int reso){resolution_int = reso;} + /** \brief Sets a gain from a value stored in memory (default value 0.3). + * + * @param[in] g gain. + * + */ + void set_gain(double g = 0.3){gain = g;} + + public: // Automatic tuning of resolution for Mapper Point. + /** \brief Computes the optimal length of intervals for a standard Mapper. + */ + void set_automatic_resolution_for_MAP(const double& gain){ + double reso = 0; + 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); + } + } + if(verbose) std::cout << "resolution = " << reso << std::endl; + resolution_double = reso; + } + + public: // Set cover with preimages of function. + /** \brief Creates a cover C from the preimages of the function f. + * + * @param[in] token boolean specifying whether we use the length or the number of intervals for the cover of im(f). + * + */ + void set_cover_from_function(const bool& token){ + + // 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);} + int n = func.size(); if(verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl; + + // Compute cover of im(f) + std::vector<std::pair<double,double> > intervals; int res; + + if(!token){ // 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<double, double> 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<double, double> inter(x,y); intervals.push_back(inter); + } + x = minf + (resolution_int-1)*incr - alpha; y = maxf; + std::pair<double, double> interM(x,y); intervals.push_back(interM); res = intervals.size(); + for(int i = 0; i < res; i++) + if(verbose) std::cout << "Interval " << i << " = [" << intervals[i].first << ", " << intervals[i].second << "]" << std::endl; + } + + else{ // 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<double, double> inter(x,y); intervals.push_back(inter); + x = y - gain*resolution_double; + y = x + resolution_double; + } + std::pair<double, double> interM(x,maxf); intervals.push_back(interM); res = intervals.size(); + for(int i = 0; i < res; i++) if(verbose) std::cout << "Interval " << i << " = [" << \ + intervals[i].first << ", " << intervals[i].second << "]" << std::endl; + } + + // Sort points according to function values + std::vector<int> points(n); for(int i = 0; i < n; i++) points[i] = i; + std::sort(points.begin(),points.end(),functional_comp); + + int id = 0; int pos = 0; + + for(int i = 0; i < res; i++){ + + // Find points in the preimage + std::map<int,std::vector<int> > prop; prop.clear(); + std::pair<double, double> inter1 = intervals[i]; + int tmp = pos; + + if(i != res-1){ + + if(i != 0){ + std::pair<double, double> inter3 = intervals[i-1]; + while(func[points[tmp]] < inter3.second && tmp != n){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + } + + std::pair<double, double> inter2 = intervals[i+1]; + while(func[points[tmp]] < inter2.first && tmp != n){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + + pos = tmp; + while(func[points[tmp]] < inter1.second && tmp != n){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + + } + + else{ + + std::pair<double, double> inter3 = intervals[i-1]; + while(func[points[tmp]] < inter3.second && tmp != n){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + + while(tmp != n){ + prop.insert(std::pair<int,std::vector<int> >(points[tmp],adjacency_matrix[points[tmp]])); + tmp++; + } + + } + + // 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.insert(std::pair<int,bool>(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); average_col += func_color[cc[j]]/cci;} + cover_fct[id] = i; cover_color[id] = std::pair<int,double>(cci,average_col); + id++; + } + } + } + if(verbose) std::cout << std::endl; + } + + maximal_dim = id-1; + + } + + // ******************************************************************************************************************* + // Visualization. + // ******************************************************************************************************************* + + public: // Set color from file. + /** \brief Computes the function used to color the nodes of the simplicial complex from a file containing the function values. + * + * @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; double f; + while(std::getline(input,line)){ + std::stringstream stream(line); stream >> f; + func_color.insert(std::pair<int,double>(vertex_id, f)); vertex_id++; + } + color_name = color_file_name; + } + + public: // Set color from kth coordinate + /** \brief Computes the function used to color the nodes of the simplicial complex from the k-th coordinate. + * + * @param[in] k coordinate to use (start at 0). + * @param[in] off_file_name name of the input .OFF file. + * + */ + void set_color_from_coordinate(int k = 0){ + for(int i = 0; i < n; i++) func_color.insert(std::pair<int,double>(i, point_cloud[i][k])); + char coordinate[100]; sprintf(coordinate, "coordinate %d", k); + color_name = coordinate; + } + + public: // Set color from vector. + /** \brief Computes the function used to color the nodes of the simplicial complex from a vector stored in memory. + * + * @param[in] color input vector of values. + * + */ + void set_color_from_vector(const std::vector<double>& color){ + for(unsigned int i = 0; i < color.size(); i++) func_color.insert(std::pair<int,double>(i, color[i])); + } + + public: // Create a .dot file that can be compiled with neato to produce a .pdf file. + /** \brief Creates a .dot file for neato once the simplicial complex is computed to get a nice visualization + * of its 1-skeleton in a .pdf file. + */ + void plot_with_neato(){ + char mapp[11] = "SC.dot"; std::ofstream graphic(mapp); graphic << "graph Mapper {" << 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++){ + 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++){ + 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=\"" \ + << (1-(maxv-iit->second.second)/(maxv-minv))*0.6 << ", 1, 1\"]" << std::endl; + k++; + } + } + int ke = 0; int num_simplices = simplices.size(); + 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; ke++;} + graphic << "}"; graphic.close(); + char command[100]; sprintf(command, "neato SC.dot -Tpdf -o SC_visu.pdf && rm SC.dot"); + int systemRet = system(command); + if(systemRet == -1) std::cout << "Visualization failed. Do you have neato?" << std::endl; + } + + public: // Create a .txt file that can be compiled with KeplerMapper to produce a .html file. + /** \brief Creates a .txt file for KeplerMapper once the simplicial complex is computed to get a nice visualization + * of its 1-skeleton in browser. + */ + void plot_with_KeplerMapper(){ + + int num_simplices = simplices.size(); int num_edges = 0; + char mapp[11] = "SC.txt"; 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 << 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 (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; + + 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.close(); + char command[100]; sprintf(command, "python visu.py && firefox SC.html"); + int systemRet = system(command); + if(systemRet == -1) std::cout << "Visualization failed. Do you have python and firefox?" << std::endl; + } + + + public: // Create a .off file that can be visualized with Geomview. + /** \brief Creates a .off file for visualization with Geomview. + * For GIC computed with Voronoi only. + */ + void plot_with_Geomview(){ + + assert(data_dimension <= 3); + 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(); + for (int i = 0; i < numsimplices; i++) { + if(simplices[i].size() == 2){ numedges++; + edges.push_back(simplices[i]); + } + if(simplices[i].size() == 3){ numfaces++; + faces.push_back(simplices[i]); + } + } + graphic << m << " " << numedges + numfaces << std::endl; + for(int i = 0; i < m; i++) graphic << point_cloud[voronoi_subsamples[i]][0] << " " \ + << point_cloud[voronoi_subsamples[i]][1] << " " \ + << point_cloud[voronoi_subsamples[i]][2] << std::endl; + for(int i = 0; i < numedges; i++) graphic << 2 << " " << edges[i][0] << " " << edges[i][1] << std::endl; + for(int i = 0; i < numfaces; i++) graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl; + graphic.close(); + + char command[100]; sprintf(command, "geomview SC.off"); + int systemRet = system(command); + if(systemRet == -1) std::cout << "Visualization failed. Do you have geomview?" << std::endl; + + } + + // ******************************************************************************************************************* + // ******************************************************************************************************************* + + + public: + /** \brief Creates the simplicial complex. + * + * @param[in] complex SimplicialComplexForGIC to be created. + * + */ + template<typename SimplicialComplexForGIC> + void create_complex(SimplicialComplexForGIC & complex) { + size_t sz = simplices.size(); unsigned int dimension = 0; + for(unsigned int i = 0; i < sz; i++){ + complex.insert_simplex_and_subfaces(simplices[i]); + if(dimension < simplices[i].size()-1) dimension = simplices[i].size()-1; + } + complex.set_dimension(dimension); + } + + public: + /** \brief Finds the maximal clique formed by different elements of the cover in a set of points. + * + * @param[in] cover_elts vector of points represented by vectors of cover elements (the ones to which they belong). + * + */ + void find_all_simplices(const std::vector<std::vector<Cover_t> > & cover_elts){ + int num_nodes = cover_elts.size(); + std::vector<Cover_t> simplex; + for(int i = 0; i < num_nodes; i++) + for(unsigned int j = 0; j < cover_elts[i].size(); j++) + simplex.push_back(cover_elts[i][j]); + std::sort(simplex.begin(),simplex.end()); std::vector<Cover_t>::iterator it = std::unique(simplex.begin(),simplex.end()); + simplex.resize(std::distance(simplex.begin(),it)); + simplices.push_back(simplex); + } + + public: + /** \brief Computes the simplices in the Nerve of the cover C. + */ + void find_Nerve_simplices(){ + for(std::map<int,std::vector<Cover_t> >::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()); + simplices.resize(std::distance(simplices.begin(),it)); + } + + public: + /** \brief Computes the simplices in the GIC of the graph G and the cover C. + */ + void find_GIC_simplices() { + + // 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[1].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); + } + + // Build the Simplex Tree corresponding to the graph + st.expansion(maximal_dim); + + // Find simplices of GIC + simplices.clear(); + for (auto simplex : st.complex_simplex_range()) { + if(!st.has_children(simplex)){ + std::vector<std::vector<Cover_t> > cover_elts; + for (auto vertex : st.simplex_vertex_range(simplex)) cover_elts.push_back(cover[vertex]); + find_all_simplices(cover_elts); + } + } + std::vector<std::vector<Cover_t> >::iterator it; + std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); + simplices.resize(std::distance(simplices.begin(),it)); + } + + public: + /** \brief Computes the simplices in the Mapper Delta by looking at all the edges of the graph + * and adding the corresponding edges in the Mapper Delta if the images of the endpoints belong + * to consecutive intervals. + * + * \remark WARNING: the output of this function is correct ONLY if the cover is minimal, i.e. + * the gain is less than 0.5!!! + * + */ + void find_GICMAP_simplices_with_functional_minimal_cover(){ + + int v1, v2; assert(gain < 0.5); + + // 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); + } + } + } + std::vector<std::vector<Cover_t> >::iterator it; + std::sort(simplices.begin(),simplices.end()); it = std::unique(simplices.begin(),simplices.end()); + simplices.resize(std::distance(simplices.begin(),it)); + } + +}; + +} // namespace graph_induced_complex + +} // namespace Gudhi + +#endif // GIC_H_ diff --git a/src/Nerve_GIC/test/CMakeLists.txt b/src/Nerve_GIC/test/CMakeLists.txt new file mode 100644 index 00000000..627778fa --- /dev/null +++ b/src/Nerve_GIC/test/CMakeLists.txt @@ -0,0 +1,23 @@ +cmake_minimum_required(VERSION 2.6) +project(Graph_induced_complex_tests) + +if (GCOVR_PATH) + # for gcovr to make coverage reports - Corbera Jenkins plugin + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fprofile-arcs -ftest-coverage") +endif() +if (GPROF_PATH) + # for gprof to make coverage reports - Jenkins + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pg") +endif() + +add_executable ( graph_induced_complex_UT test_GIC.cpp ) +target_link_libraries(graph_induced_complex_UT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY}) +if (TBB_FOUND) + target_link_libraries(graph_induced_complex_UT ${TBB_LIBRARIES}) +endif() + +file(COPY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/) + +add_test(graph_induced_complex_UT ${CMAKE_CURRENT_BINARY_DIR}/graph_induced_complex_UT + # XML format for Jenkins xUnit plugin + --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/graph_induced_complex_UT.xml --log_level=test_suite --report_level=no) diff --git a/src/Nerve_GIC/test/data/cloud b/src/Nerve_GIC/test/data/cloud new file mode 100644 index 00000000..82fc5c79 --- /dev/null +++ b/src/Nerve_GIC/test/data/cloud @@ -0,0 +1,5 @@ +OFF +3 0 0 +0 0 0 +2 1 0 +4 0 0
\ No newline at end of file diff --git a/src/Nerve_GIC/test/data/cover b/src/Nerve_GIC/test/data/cover new file mode 100644 index 00000000..5f5fbe75 --- /dev/null +++ b/src/Nerve_GIC/test/data/cover @@ -0,0 +1,3 @@ +1 +2 +3
\ No newline at end of file diff --git a/src/Nerve_GIC/test/data/graph b/src/Nerve_GIC/test/data/graph new file mode 100644 index 00000000..37142800 --- /dev/null +++ b/src/Nerve_GIC/test/data/graph @@ -0,0 +1,3 @@ +0 1 +0 2 +1 2
\ No newline at end of file diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp new file mode 100644 index 00000000..e8051cf0 --- /dev/null +++ b/src/Nerve_GIC/test/test_GIC.cpp @@ -0,0 +1,93 @@ +/* This file is part of the Gudhi Library. The Gudhi library + * (Geometric Understanding in Higher Dimensions) is a generic C++ + * library for computational topology. + * + * Author(s): Mathieu Carrière + * + * Copyright (C) 2017 INRIA Saclay (France) + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE "graph_induced_complex" + +#include <boost/test/unit_test.hpp> +#include <cmath> // float comparison +#include <limits> +#include <string> +#include <vector> +#include <algorithm> // std::max +#include <gudhi/GIC.h> +#include <gudhi/distance_functions.h> +#include <gudhi/reader_utils.h> + +bool are_almost_the_same(float a, float b) { + return std::fabs(a - b) < std::numeric_limits<float>::epsilon(); +} + +BOOST_AUTO_TEST_CASE(check_nerve) { + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); + std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name); + std::string cover_file_name("data/cover"); GIC.set_cover_from_file(cover_file_name); + GIC.set_color_from_coordinate(); + GIC.find_Nerve_simplices(); Simplex_tree stree; GIC.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 3); + BOOST_CHECK((stree.num_simplices()-stree.num_vertices()) == 0); + BOOST_CHECK(stree.dimension() == 0); +} + +BOOST_AUTO_TEST_CASE(check_GICMAP) { + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); GIC.set_color_from_coordinate(); + std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name); + std::string cover_file_name("data/cover"); GIC.set_cover_from_file(cover_file_name); + GIC.find_GICMAP_simplices_with_functional_minimal_cover(); Simplex_tree stree; GIC.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 3); + BOOST_CHECK((stree.num_simplices()-stree.num_vertices()) == 2); + BOOST_CHECK(stree.dimension() == 1); +} + +BOOST_AUTO_TEST_CASE(check_GICcover) { + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); GIC.set_color_from_coordinate(); + std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name); + std::string cover_file_name("data/cover"); GIC.set_cover_from_file(cover_file_name); + GIC.find_GIC_simplices(); Simplex_tree stree; GIC.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 3); + BOOST_CHECK((stree.num_simplices()-stree.num_vertices()) == 4); + BOOST_CHECK(stree.dimension() == 2); +} + + +BOOST_AUTO_TEST_CASE(check_GICvoronoi) { + + Gudhi::graph_induced_complex::Graph_induced_complex GIC; + std::string cloud_file_name("data/cloud"); GIC.read_point_cloud(cloud_file_name); GIC.set_color_from_coordinate(); + std::string graph_file_name("data/graph"); GIC.set_graph_from_file(graph_file_name); + GIC.set_cover_from_Voronoi(2); + GIC.find_GIC_simplices(); Simplex_tree stree; GIC.create_complex(stree); + + BOOST_CHECK(stree.num_vertices() == 2); + BOOST_CHECK((stree.num_simplices()-stree.num_vertices()) == 1); + BOOST_CHECK(stree.dimension() == 1); +} + |