summaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/CMakeLists.txt1
-rw-r--r--src/Doxyfile3
-rw-r--r--src/Nerve_GIC/doc/COPYRIGHT19
-rw-r--r--src/Nerve_GIC/doc/Intro_graph_induced_complex.h187
-rw-r--r--src/Nerve_GIC/doc/gic_complex.pngbin0 -> 56761 bytes
-rw-r--r--src/Nerve_GIC/doc/nerve.pngbin0 -> 45129 bytes
-rw-r--r--src/Nerve_GIC/example/CMakeLists.txt27
-rw-r--r--src/Nerve_GIC/example/GIC.cpp65
-rw-r--r--src/Nerve_GIC/example/GIC.txt89
-rw-r--r--src/Nerve_GIC/example/GICvoronoi.cpp60
-rw-r--r--src/Nerve_GIC/example/GICvoronoi.txt563
-rw-r--r--src/Nerve_GIC/example/MapperDeltaCoord.cpp62
-rw-r--r--src/Nerve_GIC/example/MapperDeltaCoord.txt143
-rw-r--r--src/Nerve_GIC/example/MapperDeltaFunc.cpp62
-rw-r--r--src/Nerve_GIC/example/MapperDeltaFunc.txt10
-rw-r--r--src/Nerve_GIC/example/Nerve.cpp64
-rw-r--r--src/Nerve_GIC/example/Nerve.txt43
-rwxr-xr-xsrc/Nerve_GIC/example/km.py390
-rwxr-xr-xsrc/Nerve_GIC/example/visu.py44
-rw-r--r--src/Nerve_GIC/include/gudhi/GIC.h911
-rw-r--r--src/Nerve_GIC/test/CMakeLists.txt23
-rw-r--r--src/Nerve_GIC/test/data/cloud5
-rw-r--r--src/Nerve_GIC/test/data/cover3
-rw-r--r--src/Nerve_GIC/test/data/graph3
-rw-r--r--src/Nerve_GIC/test/test_GIC.cpp93
-rw-r--r--src/cmake/modules/GUDHI_user_version_target.txt2
-rw-r--r--src/common/doc/main_page.h18
27 files changed, 2888 insertions, 2 deletions
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ebcb6888..545a606b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -52,6 +52,7 @@ add_subdirectory(example/Spatial_searching)
add_subdirectory(example/Subsampling)
add_subdirectory(example/Tangential_complex)
add_subdirectory(example/Bottleneck_distance)
+add_subdirectory(example/Nerve_GIC)
# data points generator
add_subdirectory(data/points/generator)
diff --git a/src/Doxyfile b/src/Doxyfile
index 9cd3894c..904a9a99 100644
--- a/src/Doxyfile
+++ b/src/Doxyfile
@@ -851,7 +851,8 @@ IMAGE_PATH = doc/Skeleton_blocker/ \
doc/Subsampling/ \
doc/Spatial_searching/ \
doc/Tangential_complex/ \
- doc/Bottleneck_distance/
+ doc/Bottleneck_distance/ \
+ doc/Nerve_GIC/
# The INPUT_FILTER tag can be used to specify a program that doxygen should
# invoke to filter for each input file. Doxygen will invoke the filter program
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
new file mode 100644
index 00000000..fb4b20ad
--- /dev/null
+++ b/src/Nerve_GIC/doc/gic_complex.png
Binary files differ
diff --git a/src/Nerve_GIC/doc/nerve.png b/src/Nerve_GIC/doc/nerve.png
new file mode 100644
index 00000000..b66da4a4
--- /dev/null
+++ b/src/Nerve_GIC/doc/nerve.png
Binary files differ
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);
+}
+
diff --git a/src/cmake/modules/GUDHI_user_version_target.txt b/src/cmake/modules/GUDHI_user_version_target.txt
index ca6bcd34..0d409cc3 100644
--- a/src/cmake/modules/GUDHI_user_version_target.txt
+++ b/src/cmake/modules/GUDHI_user_version_target.txt
@@ -47,7 +47,7 @@ if (NOT CMAKE_VERSION VERSION_LESS 2.8.11)
add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
copy_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI)
- set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Bottleneck_distance;Contraction;Hasse_complex;Persistent_cohomology;Rips_complex;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex")
+ set(GUDHI_MODULES "common;Alpha_complex;Bitmap_cubical_complex;Bottleneck_distance;Contraction;Hasse_complex;Nerve_GIC;Persistent_cohomology;Rips_complex;Simplex_tree;Skeleton_blocker;Spatial_searching;Subsampling;Tangential_complex;Witness_complex")
set(GUDHI_DIRECTORIES "doc;example;concept")
set(GUDHI_INCLUDE_DIRECTORIES "include/gudhi;include/gudhi_patches")
diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h
index bd4615f5..8d282fca 100644
--- a/src/common/doc/main_page.h
+++ b/src/common/doc/main_page.h
@@ -93,6 +93,24 @@
</td>
</tr>
</table>
+ \subsection GICDataStructure Nerves and Graph Induced Complexes
+ \image html "gic_complex.png" "Graph Induced Complex of a point cloud."
+<table border="0">
+ <tr>
+ <td width="25%">
+ <b>Author:</b> Mathieu Carri&egrave;re<br>
+ <b>Introduced in:</b> GUDHI 2.0.1<br>
+ <b>Copyright:</b> GPL v3<br>
+ </td>
+ <td width="75%">
+ Nerves and Graph Induced Complexes are simplicial complexes that provably contain
+ topological information about the input data. They can be computed with a cover of the
+ data, that often comes from the preimage of a family of intervals covering the image
+ of a scalar-valued function defined on the data. <br>
+ <b>User manual:</b> \ref graph_induced_complex - <b>Reference manual:</b> Gudhi::graph_induced_complex::Graph_induced_complex
+ </td>
+ </tr>
+</table>
\subsection SkeletonBlockerDataStructure Skeleton blocker
\image html "ds_representation.png" "Skeleton blocker representation"
<table border="0">