Merge from trunk.

hand-merge README that moved to *.md.


git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-glisse@3201 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: f1922a265c79937e507bbb963875d21fae88069e

20 files changed, 734 insertions, 603 deletions

diff --git a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
index 3a0d8154..2b648425 100644
--- a/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
+++ b/src/Nerve_GIC/doc/Intro_graph_induced_complex.h
@@ -70,15 +70,15 @@ namespace cover_complex {
  *
  * When launching:
  *
- * \code $> ./Nerve ../../../../data/points/human.off 2 10 0.3 --v
+ * \code $> ./Nerve ../../data/points/human.off 2 10 0.3 -v
  * \endcode
  *
  * the program output is:
  *
  * \include Nerve_GIC/Nerve.txt
  *
- * The program also writes a file SC.txt. The first three lines in this file are the location of the input point cloud
- * and the function used to compute the cover.
+ * The program also writes a file ../../data/points/human_sc.txt. The first three lines in this file are the location
+ * of the input point cloud and the function used to compute the cover.
  * 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.
@@ -113,12 +113,12 @@ namespace cover_complex {
  *
  * When launching:
  *
- * \code $> ./VoronoiGIC ../../../../data/points/human.off 700 --v
+ * \code $> ./VoronoiGIC ../../data/points/human.off 700 -v
  * \endcode
  *
  * the program outputs SC.off. Using e.g.
  *
- * \code $> geomview SC.off
+ * \code $> geomview ../../data/points/human_sc.off
  * \endcode
  *
  * one can obtain the following visualization:
@@ -146,7 +146,7 @@ namespace cover_complex {
  *
  * When launching:
  *
- * \code $> ./CoordGIC ../../../../data/points/KleinBottle5D.off 0 --v
+ * \code $> ./CoordGIC ../../data/points/KleinBottle5D.off 0 -v
  * \endcode
  *
  * the program outputs SC.dot. Using e.g.
@@ -169,15 +169,13 @@ namespace cover_complex {
  *
  * When launching:
  *
- * \code $> ./FuncGIC ../../data/points/COIL_database/lucky_cat.off ../../data/points/COIL_database/lucky_cat_PCA1 --v
+ * \code $> ./FuncGIC ../../data/points/COIL_database/lucky_cat.off ../../data/points/COIL_database/lucky_cat_PCA1 -v
  * \endcode
  *
  * the program outputs again SC.dot which gives the following visualization after using neato:
  *
  * \image html "funcGICvisu.jpg" "Visualization with neato"
  *
- * \copyright GNU General Public License v3.                         
- * \verbatim  Contact: gudhi-users@lists.gforge.inria.fr \endverbatim
  */
 /** @} */  // end defgroup cover_complex
 
@@ -186,31 +184,3 @@ namespace cover_complex {
 }  // namespace Gudhi
 
 #endif  // DOC_COVER_COMPLEX_INTRO_COVER_COMPLEX_H_
-
-
-/* * \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 outputs SC.txt, which can be visualized with python and firefox as before:
- *
- * \image html "gicvisu.jpg" "Visualization with KeplerMapper"
- * */
-
-
-/* * Using e.g.
- *
- * \code $> python KeplerMapperVisuFromTxtFile.py && firefox SC.html
- * \endcode */
diff --git a/src/Nerve_GIC/doc/funcGICvisu.jpg b/src/Nerve_GIC/doc/funcGICvisu.jpg
index f3da45ac..36b64dde 100644
--- a/src/Nerve_GIC/doc/funcGICvisu.jpg
+++ b/src/Nerve_GIC/doc/funcGICvisu.jpg
diff --git a/src/Nerve_GIC/doc/funcGICvisu.pdf b/src/Nerve_GIC/doc/funcGICvisu.pdf
new file mode 100644
index 00000000..d7456cd3
--- /dev/null
+++ b/src/Nerve_GIC/doc/funcGICvisu.pdf
diff --git a/src/Nerve_GIC/example/CMakeLists.txt b/src/Nerve_GIC/example/CMakeLists.txt
index 461b6db2..f2626927 100644
--- a/src/Nerve_GIC/example/CMakeLists.txt
+++ b/src/Nerve_GIC/example/CMakeLists.txt
@@ -1,29 +1,26 @@
 cmake_minimum_required(VERSION 2.6)
 project(Nerve_GIC_examples)
 
-add_executable ( Nerve Nerve.cpp )
-add_executable ( CoordGIC CoordGIC.cpp )
-add_executable ( FuncGIC FuncGIC.cpp )
-add_executable ( VoronoiGIC VoronoiGIC.cpp )
+if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
 
-if (TBB_FOUND)
-  target_link_libraries(Nerve ${TBB_LIBRARIES})
-  target_link_libraries(CoordGIC ${TBB_LIBRARIES})
-  target_link_libraries(FuncGIC ${TBB_LIBRARIES})
-  target_link_libraries(VoronoiGIC ${TBB_LIBRARIES})
-endif()
+  add_executable ( CoordGIC CoordGIC.cpp )
+  add_executable ( FuncGIC FuncGIC.cpp )
 
-file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+  if (TBB_FOUND)
+    target_link_libraries(CoordGIC ${TBB_LIBRARIES})
+    target_link_libraries(FuncGIC ${TBB_LIBRARIES})
+  endif()
 
-add_test(NAME Nerve_GIC_example_nerve COMMAND $<TARGET_FILE:Nerve>
-    "${CMAKE_SOURCE_DIR}/data/points/human.off" "2" "10" "0.3")
+  # Copy files for not to pollute sources when testing
+  file(COPY "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+  file(COPY "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+  file(COPY "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat_PCA1" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
 
-add_test(NAME Nerve_GIC_example_VoronoiGIC COMMAND $<TARGET_FILE:VoronoiGIC>
-    "${CMAKE_SOURCE_DIR}/data/points/human.off" "100")
+  add_test(NAME Nerve_GIC_example_CoordGIC COMMAND $<TARGET_FILE:CoordGIC>
+      "tore3D_1307.off" "0")
 
-add_test(NAME Nerve_GIC_example_CoordGIC COMMAND $<TARGET_FILE:CoordGIC>
-    "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "0")
+  add_test(NAME Nerve_GIC_example_FuncGIC COMMAND $<TARGET_FILE:FuncGIC>
+      "lucky_cat.off"
+      "lucky_cat_PCA1")
 
-add_test(NAME Nerve_GIC_example_FuncGIC COMMAND $<TARGET_FILE:FuncGIC>
-    "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat.off"
-    "${CMAKE_SOURCE_DIR}/data/points/COIL_database/lucky_cat_PCA1")
+endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
diff --git a/src/Nerve_GIC/example/CoordGIC.cpp b/src/Nerve_GIC/example/CoordGIC.cpp
index c03fcbb3..c92cf235 100644
--- a/src/Nerve_GIC/example/CoordGIC.cpp
+++ b/src/Nerve_GIC/example/CoordGIC.cpp
@@ -27,8 +27,8 @@
 
 void usage(int nbArgs, char *const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off coordinate [--v] \n";
-  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 2 --v \n";
+  std::cerr << "Usage: " << progName << " filename.off coordinate [-v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 2 -v \n";
   exit(-1);  // ----- >>
 }
 
diff --git a/src/Nerve_GIC/example/FuncGIC.cpp b/src/Nerve_GIC/example/FuncGIC.cpp
index 3762db4e..cb0f0d63 100644
--- a/src/Nerve_GIC/example/FuncGIC.cpp
+++ b/src/Nerve_GIC/example/FuncGIC.cpp
@@ -27,9 +27,9 @@
 
 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 << "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";
+                                              "../../data/points/COIL_database/lucky_cat_PCA1 -v \n";
   exit(-1);  // ----- >>
 }
 
diff --git a/src/Nerve_GIC/example/GIC.cpp b/src/Nerve_GIC/example/GIC.cpp
deleted file mode 100644
index 2bc24a4d..00000000
--- a/src/Nerve_GIC/example/GIC.cpp
+++ /dev/null
@@ -1,95 +0,0 @@
-/*    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
- *
- *    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/>.
- */
-
-#include <gudhi/GIC.h>
-
-#include <string>
-#include <vector>
-
-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]);
-
-  using Point = std::vector<float>;
-
-  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<Point> GIC;
-  GIC.set_verbose(verb);
-
-  bool check = GIC.read_point_cloud(off_file_name);
-
-  if (!check) {
-    std::cout << "Incorrect OFF file." << std::endl;
-  } else {
-    GIC.set_color_from_coordinate(coord);
-    GIC.set_function_from_coordinate(coord);
-
-    GIC.set_graph_from_rips(threshold, Gudhi::Euclidean_distance());
-
-    GIC.set_resolution_with_interval_length(resolution);
-    GIC.set_gain(gain);
-    GIC.set_cover_from_function();
-
-    GIC.find_GIC_simplices();
-
-    GIC.plot_TXT_for_KeplerMapper();
-
-    Gudhi::Simplex_tree<> stree;
-    GIC.create_complex(stree);
-
-    // ----------------------------------------------------------------------------
-    // Display information about the graph induced complex
-    // ----------------------------------------------------------------------------
-
-    if (verb) {
-      std::cout << "Graph induced complex is of dimension " << stree.dimension() << " - " << stree.num_simplices()
-                << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
-
-      std::cout << "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)) {
-          std::cout << vertex << " ";
-        }
-        std::cout << std::endl;
-      }
-    }
-  }
-
-  return 0;
-}
diff --git a/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py b/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py
deleted file mode 100755
index d2897774..00000000
--- a/src/Nerve_GIC/example/KeplerMapperVisuFromTxtFile.py
+++ /dev/null
@@ -1,72 +0,0 @@
-#!/usr/bin/env python
-
-import km
-import numpy as np
-from collections import defaultdict
-
-"""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 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/>.
-"""
-
-__author__ = "Mathieu Carriere"
-__copyright__ = "Copyright (C) 2017 INRIA"
-__license__ = "GPL v3"
-
-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)
diff --git a/src/Nerve_GIC/example/Nerve.txt b/src/Nerve_GIC/example/Nerve.txt
deleted file mode 100644
index 2a861c5f..00000000
--- a/src/Nerve_GIC/example/Nerve.txt
+++ /dev/null
@@ -1,43 +0,0 @@
-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/include/gudhi/GIC.h b/src/Nerve_GIC/include/gudhi/GIC.h
index 9f107a7e..40ff7a4a 100644
--- a/src/Nerve_GIC/include/gudhi/GIC.h
+++ b/src/Nerve_GIC/include/gudhi/GIC.h
@@ -30,16 +30,27 @@
 #include <gudhi/Rips_complex.h>
 #include <gudhi/Points_off_io.h>
 #include <gudhi/distance_functions.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Bottleneck.h>
+
+#include <boost/config.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/graph/connected_components.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+#include <boost/graph/subgraph.hpp>
+#include <boost/graph/graph_utility.hpp>
 
 #include <iostream>
 #include <vector>
 #include <map>
 #include <string>
 #include <limits>     // for numeric_limits
-#include <utility>    // for pair<>
+#include <utility>    // for std::pair<>
 #include <algorithm>  // for std::max
 #include <random>
 #include <cassert>
+#include <cmath>
 
 namespace Gudhi {
 
@@ -48,6 +59,13 @@ namespace cover_complex {
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
 using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
+using Persistence_diagram = std::vector<std::pair<double, double> >;
+using Graph = boost::subgraph<
+    boost::adjacency_list<boost::setS, boost::vecS, boost::undirectedS, boost::no_property,
+                          boost::property<boost::edge_index_t, int, boost::property<boost::edge_weight_t, double> > > >;
+using Vertex_t = boost::graph_traits<Graph>::vertex_descriptor;
+using Index_map = boost::property_map<Graph, boost::vertex_index_t>::type;
+using Weight_map = boost::property_map<Graph, boost::edge_weight_t>::type;
 
 /**
  * \class Cover_complex
@@ -72,25 +90,40 @@ using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
 template <typename Point>
 class Cover_complex {
  private:
-  // Graph_induced_complex(std::map<int, double> fun){func = fun;}
   bool verbose = false;  // whether to display information.
-  std::vector<Point> point_cloud;
-  std::vector<std::vector<int> > one_skeleton;
-  typedef int Cover_t;  // elements of cover C are indexed by integers.
-  std::vector<std::vector<Cover_t> > simplices;
-  std::map<int, std::vector<Cover_t> > cover;
-  std::map<Cover_t, std::vector<int> > cover_back;
-  int maximal_dim;     // maximal dimension of output simplicial complex.
-  int data_dimension;  // dimension of input data.
-  int n;               // number of points.
-  std::map<Cover_t, int>
+  std::string type;      // Nerve or GIC
+
+  std::vector<Point> point_cloud;               // input point cloud.
+  std::vector<std::vector<double> > distances;  // all pairwise distances.
+  int maximal_dim;                              // maximal dimension of output simplicial complex.
+  int data_dimension;                           // dimension of input data.
+  int n;                                        // number of points.
+
+  std::map<int, double> func;  // function used to compute the output simplicial complex.
+  std::map<int, double>
+      func_color;                 // function used to compute the colors of the nodes of the output simplicial complex.
+  bool functional_cover = false;  // whether we use a cover with preimages of a function or not.
+
+  Graph one_skeleton_OFF;          // one-skeleton given by the input OFF file (if it exists).
+  Graph one_skeleton;              // one-skeleton used to compute the connected components.
+  std::vector<Vertex_t> vertices;  // vertices of one_skeleton.
+
+  std::vector<std::vector<int> > simplices;  // simplices of output simplicial complex.
+  std::vector<int> voronoi_subsamples;       // Voronoi germs (in case of Voronoi cover).
+
+  Persistence_diagram PD;
+  std::vector<double> distribution;
+
+  std::map<int, std::vector<int> >
+      cover;  // function associating to each data point its vectors of cover elements to which it belongs.
+  std::map<int, std::vector<int> >
+      cover_back;  // inverse of cover, in order to get the data points associated to a specific cover element.
+  std::map<int, double> cover_std;  // standard function (induced by func) used to compute the extended persistence
+                                    // diagram of the output simplicial complex.
+  std::map<int, int>
       cover_fct;  // integer-valued function that allows to state if two elements of the cover are consecutive or not.
-  std::map<Cover_t, std::pair<int, double> >
-      cover_color;  // size and coloring of the vertices of the output simplicial complex.
-  Simplex_tree st;
-
-  std::map<int, std::vector<int> > adjacency_matrix;
-  std::vector<std::vector<double> > distances;
+  std::map<int, std::pair<int, double> >
+      cover_color;  // size and coloring (induced by func_color) of the vertices of the output simplicial complex.
 
   int resolution_int = -1;
   double resolution_double = -1;
@@ -99,14 +132,11 @@ class Cover_complex {
   double rate_power = 0.001;  // Power in the subsampling.
   int mask = 0;               // Ignore nodes containing less than mask points.
 
-  std::map<int, double> func;
-  std::map<int, double> func_color;
-  std::vector<int> voronoi_subsamples;
+  std::map<int, int> name2id, name2idinv;
+
   std::string cover_name;
   std::string point_cloud_name;
   std::string color_name;
-  std::string type;               // Nerve or GIC
-  bool functional_cover = false;  // whether we use a cover with preimages of a function or not
 
   // Point comparator
   struct Less {
@@ -120,21 +150,16 @@ class Cover_complex {
     }
   };
 
-  // DFS
- private:
-  void dfs(std::map<int, std::vector<int> >& G, int p, std::vector<int>& cc, std::map<int, bool>& visit) {
-    cc.push_back(p);
-    visit[p] = true;
-    int neighb = G[p].size();
-    for (int i = 0; i < neighb; i++)
-      if (visit.find(G[p][i]) != visit.end())
-        if (!(visit[G[p][i]])) dfs(G, G[p][i], cc, visit);
+  // Remove all edges of a graph.
+  void remove_edges(Graph& G) {
+    boost::graph_traits<Graph>::edge_iterator ei, ei_end;
+    for (boost::tie(ei, ei_end) = boost::edges(G); ei != ei_end; ++ei) boost::remove_edge(*ei, G);
   }
 
   // Find random number in [0,1].
   double GetUniform() {
-    static std::default_random_engine re;
-    static std::uniform_real_distribution<double> Dist(0, 1);
+    thread_local std::default_random_engine re;
+    thread_local std::uniform_real_distribution<double> Dist(0, 1);
     return Dist(re);
   }
 
@@ -155,24 +180,14 @@ class Cover_complex {
     }
   }
 
- private:
-  void fill_adjacency_matrix_from_st() {
-    std::vector<int> empty;
-    for (int i = 0; i < n; i++) adjacency_matrix[i] = empty;
-    for (auto simplex : st.complex_simplex_range()) {
-      if (st.dimension(simplex) == 1) {
-        std::vector<int> vertices;
-        for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex);
-        adjacency_matrix[vertices[0]].push_back(vertices[1]);
-        adjacency_matrix[vertices[1]].push_back(vertices[0]);
-      }
-    }
-  }
+  // *******************************************************************************************************************
+  // Utils.
+  // *******************************************************************************************************************
 
  public:
   /** \brief Specifies whether the type of the output simplicial complex.
    *
-   * @param[in] t string (either "GIC" or "Nerve").
+   * @param[in] t std::string (either "GIC" or "Nerve").
    *
    */
   void set_type(const std::string& t) { type = t; }
@@ -221,14 +236,15 @@ class Cover_complex {
 
     char comment = '#';
     while (comment == '#') {
-      getline(input, line);
-      if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) comment = line[line.find_first_not_of(' ')];
+      std::getline(input, line);
+      if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
+        comment = line[line.find_first_not_of(' ')];
     }
-    if (std::strcmp((char*)line.c_str(), "nOFF") == 0) {
+    if (strcmp((char*)line.c_str(), "nOFF") == 0) {
       comment = '#';
       while (comment == '#') {
-        getline(input, line);
-        if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace))
+        std::getline(input, line);
+        if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
           comment = line[line.find_first_not_of(' ')];
       }
       std::stringstream stream(line);
@@ -238,10 +254,11 @@ class Cover_complex {
     }
 
     comment = '#';
-    int numedges, numfaces, i, num;
+    int numedges, numfaces, i, dim;
     while (comment == '#') {
-      getline(input, line);
-      if (!line.empty() && !std::all_of(line.begin(), line.end(), isspace)) comment = line[line.find_first_not_of(' ')];
+      std::getline(input, line);
+      if (!line.empty() && !all_of(line.begin(), line.end(), (int (*)(int))isspace))
+        comment = line[line.find_first_not_of(' ')];
     }
     std::stringstream stream(line);
     stream >> n;
@@ -250,34 +267,31 @@ class Cover_complex {
 
     i = 0;
     while (i < n) {
-      getline(input, line);
+      std::getline(input, line);
       if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
-          !std::all_of(line.begin(), line.end(), isspace)) {
+          !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
+        std::stringstream iss(line);
         std::vector<double> point;
-        std::istringstream iss(line);
         point.assign(std::istream_iterator<double>(iss), std::istream_iterator<double>());
         point_cloud.emplace_back(point.begin(), point.begin() + data_dimension);
+        boost::add_vertex(one_skeleton_OFF);
+        vertices.push_back(boost::add_vertex(one_skeleton));
         i++;
       }
     }
 
     i = 0;
     while (i < numfaces) {
-      getline(input, line);
+      std::getline(input, line);
       if (!line.empty() && line[line.find_first_not_of(' ')] != '#' &&
-          !std::all_of(line.begin(), line.end(), isspace)) {
+          !all_of(line.begin(), line.end(), (int (*)(int))isspace)) {
         std::vector<int> simplex;
-        std::istringstream iss(line);
+        std::stringstream iss(line);
         simplex.assign(std::istream_iterator<int>(iss), std::istream_iterator<int>());
-        num = simplex[0];
-        std::vector<int> edge(2);
-        for (int j = 1; j <= num; j++) {
-          for (int k = j + 1; k <= num; k++) {
-            edge[0] = simplex[j];
-            edge[1] = simplex[k];
-            one_skeleton.push_back(edge);
-          }
-        }
+        dim = simplex[0];
+        for (int j = 1; j <= dim; j++)
+          for (int k = j + 1; k <= dim; k++)
+            boost::add_edge(vertices[simplex[j]], vertices[simplex[k]], one_skeleton_OFF);
         i++;
       }
     }
@@ -298,22 +312,16 @@ class Cover_complex {
            *
            */
   void set_graph_from_file(const std::string& graph_file_name) {
+    remove_edges(one_skeleton);
     int neighb;
     std::ifstream input(graph_file_name);
     std::string line;
-    int edge[2];
-    int n = 0;
+    int source;
     while (std::getline(input, line)) {
       std::stringstream stream(line);
-      stream >> edge[0];
-      while (stream >> neighb) {
-        edge[1] = neighb;
-        st.insert_simplex_and_subfaces(edge);
-      }
-      n++;
+      stream >> source;
+      while (stream >> neighb) boost::add_edge(vertices[source], vertices[neighb], one_skeleton);
     }
-
-    fill_adjacency_matrix_from_st();
   }
 
  public:  // Set graph from OFF file.
@@ -321,13 +329,11 @@ class Cover_complex {
            *
            */
   void set_graph_from_OFF() {
-    int num_edges = one_skeleton.size();
-    if (num_edges > 0) {
-      for (int i = 0; i < num_edges; i++) st.insert_simplex_and_subfaces(one_skeleton[i]);
-      fill_adjacency_matrix_from_st();
-    } else {
+    remove_edges(one_skeleton);
+    if (num_edges(one_skeleton_OFF))
+      one_skeleton = one_skeleton_OFF;
+    else
       std::cout << "No triangulation read in OFF file!" << std::endl;
-    }
   }
 
  public:  // Set graph from Rips complex.
@@ -339,9 +345,27 @@ class Cover_complex {
            */
   template <typename Distance>
   void set_graph_from_rips(double threshold, Distance distance) {
-    Rips_complex rips_complex_from_points(point_cloud, threshold, distance);
-    rips_complex_from_points.create_complex(st, 1);
-    fill_adjacency_matrix_from_st();
+    remove_edges(one_skeleton);
+    if (distances.size() == 0) compute_pairwise_distances(distance);
+    for (int i = 0; i < n; i++) {
+      for (int j = i + 1; j < n; j++) {
+        if (distances[i][j] <= threshold) {
+          boost::add_edge(vertices[i], vertices[j], one_skeleton);
+          boost::put(boost::edge_weight, one_skeleton, boost::edge(vertices[i], vertices[j], one_skeleton).first,
+                     distances[i][j]);
+        }
+      }
+    }
+  }
+
+ public:
+  void set_graph_weights() {
+    Index_map index = boost::get(boost::vertex_index, one_skeleton);
+    Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
+    boost::graph_traits<Graph>::edge_iterator ei, ei_end;
+    for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
+      boost::put(weight, *ei,
+                 distances[index[boost::source(*ei, one_skeleton)]][index[boost::target(*ei, one_skeleton)]]);
   }
 
  public:  // Pairwise distances.
@@ -352,9 +376,8 @@ class Cover_complex {
     double d;
     std::vector<double> zeros(n);
     for (int i = 0; i < n; i++) distances.push_back(zeros);
-    std::string distance = point_cloud_name;
-    distance.append("_dist");
-    std::ifstream input(distance.c_str(), std::ios::out | std::ios::binary);
+    std::string distance = point_cloud_name + "_dist";
+    std::ifstream input(distance, std::ios::out | std::ios::binary);
 
     if (input.good()) {
       if (verbose) std::cout << "Reading distances..." << std::endl;
@@ -420,10 +443,7 @@ class Cover_complex {
     }
 
     if (verbose) std::cout << "delta = " << delta << std::endl;
-    Rips_complex rips_complex_from_points(point_cloud, delta, distance);
-    rips_complex_from_points.create_complex(st, 1);
-    fill_adjacency_matrix_from_st();
-
+    set_graph_from_rips(delta, distance);
     return delta;
   }
 
@@ -438,15 +458,15 @@ class Cover_complex {
            *
            */
   void set_function_from_file(const std::string& func_file_name) {
-    int vertex_id = 0;
+    int i = 0;
     std::ifstream input(func_file_name);
     std::string line;
     double f;
     while (std::getline(input, line)) {
       std::stringstream stream(line);
       stream >> f;
-      func.emplace(vertex_id, f);
-      vertex_id++;
+      func.emplace(i, f);
+      i++;
     }
     functional_cover = true;
     cover_name = func_file_name;
@@ -460,10 +480,8 @@ class Cover_complex {
            */
   void set_function_from_coordinate(int k) {
     for (int i = 0; i < n; i++) func.emplace(i, point_cloud[i][k]);
-    char coordinate[100];
-    sprintf(coordinate, "coordinate %d", k);
     functional_cover = true;
-    cover_name = coordinate;
+    cover_name = "coordinate " + std::to_string(k);
   }
 
  public:  // Set function from vector.
@@ -474,12 +492,8 @@ class Cover_complex {
            */
   template <class InputRange>
   void set_function_from_range(InputRange const& function) {
+    for (int i = 0; i < n; i++) func.emplace(i, function[i]);
     functional_cover = true;
-    int index = 0;
-    for (auto v : function) {
-      func.emplace(index, v);
-      index++;
-    }
   }
 
   // *******************************************************************************************************************
@@ -505,27 +519,23 @@ class Cover_complex {
     }
 
     double reso = 0;
+    Index_map index = boost::get(boost::vertex_index, one_skeleton);
 
     if (type == "GIC") {
-      for (auto simplex : st.complex_simplex_range()) {
-        if (st.dimension(simplex) == 1) {
-          std::vector<int> vertices;
-          for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex);
-          reso = std::max(reso, std::abs(func[vertices[0]] - func[vertices[1]]));
-        }
-      }
+      boost::graph_traits<Graph>::edge_iterator ei, ei_end;
+      for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
+        reso = std::max(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
+                                       func[index[boost::target(*ei, one_skeleton)]]));
       if (verbose) std::cout << "resolution = " << reso << std::endl;
       resolution_double = reso;
     }
 
     if (type == "Nerve") {
-      for (auto simplex : st.complex_simplex_range()) {
-        if (st.dimension(simplex) == 1) {
-          std::vector<int> vertices;
-          for (auto vertex : st.simplex_vertex_range(simplex)) vertices.push_back(vertex);
-          reso = std::max(reso, (std::abs(func[vertices[0]] - func[vertices[1]])) / gain);
-        }
-      }
+      boost::graph_traits<Graph>::edge_iterator ei, ei_end;
+      for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
+        reso = std::max(reso, std::abs(func[index[boost::source(*ei, one_skeleton)]] -
+                                       func[index[boost::target(*ei, one_skeleton)]]) /
+                                  gain);
       if (verbose) std::cout << "resolution = " << reso << std::endl;
       resolution_double = reso;
     }
@@ -568,15 +578,12 @@ class Cover_complex {
     }
 
     // Read function values and compute min and max
-    std::map<int, double>::iterator it;
-    double maxf, minf;
-    minf = std::numeric_limits<float>::max();
-    maxf = std::numeric_limits<float>::min();
-    for (it = func.begin(); it != func.end(); it++) {
-      minf = std::min(minf, it->second);
-      maxf = std::max(maxf, it->second);
+    double minf = std::numeric_limits<float>::max();
+    double maxf = std::numeric_limits<float>::lowest();
+    for (int i = 0; i < n; i++) {
+      minf = std::min(minf, func[i]);
+      maxf = std::max(maxf, func[i]);
     }
-    int n = func.size();
     if (verbose) std::cout << "Min function value = " << minf << " and Max function value = " << maxf << std::endl;
 
     // Compute cover of im(f)
@@ -648,78 +655,96 @@ class Cover_complex {
     std::vector<int> points(n);
     for (int i = 0; i < n; i++) points[i] = i;
     std::sort(points.begin(), points.end(), Less(this->func));
+
     int id = 0;
     int pos = 0;
+    Index_map index = boost::get(boost::vertex_index, one_skeleton);  // int maxc = -1;
+    std::map<int, std::vector<int> > preimages;
+    std::map<int, double> funcstd;
 
+    if (verbose) std::cout << "Computing preimages..." << std::endl;
     for (int i = 0; i < res; i++) {
       // Find points in the preimage
-      std::map<int, std::vector<int> > prop;
       std::pair<double, double> inter1 = intervals[i];
       int tmp = pos;
+      double u, v;
 
       if (i != res - 1) {
         if (i != 0) {
           std::pair<double, double> inter3 = intervals[i - 1];
           while (func[points[tmp]] < inter3.second && tmp != n) {
-            prop[points[tmp]] = adjacency_matrix[points[tmp]];
+            preimages[i].push_back(points[tmp]);
             tmp++;
           }
+          u = inter3.second;
+        } else {
+          u = inter1.first;
         }
 
         std::pair<double, double> inter2 = intervals[i + 1];
         while (func[points[tmp]] < inter2.first && tmp != n) {
-          prop[points[tmp]] = adjacency_matrix[points[tmp]];
+          preimages[i].push_back(points[tmp]);
           tmp++;
         }
-
+        v = inter2.first;
         pos = tmp;
         while (func[points[tmp]] < inter1.second && tmp != n) {
-          prop[points[tmp]] = adjacency_matrix[points[tmp]];
+          preimages[i].push_back(points[tmp]);
           tmp++;
         }
 
       } else {
         std::pair<double, double> inter3 = intervals[i - 1];
         while (func[points[tmp]] < inter3.second && tmp != n) {
-          prop[points[tmp]] = adjacency_matrix[points[tmp]];
+          preimages[i].push_back(points[tmp]);
           tmp++;
         }
-
         while (tmp != n) {
-          prop[points[tmp]] = adjacency_matrix[points[tmp]];
+          preimages[i].push_back(points[tmp]);
           tmp++;
         }
+        u = inter3.second;
+        v = inter1.second;
       }
 
-      // Compute the connected components with DFS
-      std::map<int, bool> visit;
-      if (verbose) std::cout << "Preimage of interval " << i << std::endl;
-      for (std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++)
-        visit[it->first] = false;
-      if (!(prop.empty())) {
-        for (std::map<int, std::vector<int> >::iterator it = prop.begin(); it != prop.end(); it++) {
-          if (!(visit[it->first])) {
-            std::vector<int> cc;
-            cc.clear();
-            dfs(prop, it->first, cc, visit);
-            int cci = cc.size();
-            if (verbose) std::cout << "one CC with " << cci << " points, ";
-            double average_col = 0;
-            for (int j = 0; j < cci; j++) {
-              cover[cc[j]].push_back(id);
-              cover_back[id].push_back(cc[j]);
-              average_col += func_color[cc[j]] / cci;
-            }
-            cover_fct[id] = i;
-            cover_color[id] = std::pair<int, double>(cci, average_col);
-            id++;
-          }
-        }
+      funcstd[i] = 0.5 * (u + v);
+    }
+
+    if (verbose) std::cout << "Computing connected components..." << std::endl;
+    // #pragma omp parallel for
+    for (int i = 0; i < res; i++) {
+      // Compute connected components
+      Graph G = one_skeleton.create_subgraph();
+      int num = preimages[i].size();
+      std::vector<int> component(num);
+      for (int j = 0; j < num; j++) boost::add_vertex(index[vertices[preimages[i][j]]], G);
+      boost::connected_components(G, &component[0]);
+      int max = 0;
+
+      // For each point in preimage
+      for (int j = 0; j < num; j++) {
+        // Update number of components in preimage
+        if (component[j] > max) max = component[j];
+
+        // Identify component with Cantor polynomial N^2 -> N
+        int identifier = (std::pow(i + component[j], 2) + 3 * i + component[j]) / 2;
+
+        // Update covers
+        cover[preimages[i][j]].push_back(identifier);
+        cover_back[identifier].push_back(preimages[i][j]);
+        cover_fct[identifier] = i;
+        cover_std[identifier] = funcstd[i];
+        cover_color[identifier].second += func_color[preimages[i][j]];
+        cover_color[identifier].first += 1;
       }
-      if (verbose) std::cout << std::endl;
+
+      // Maximal dimension is total number of connected components
+      id += max + 1;
     }
 
     maximal_dim = id - 1;
+    for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++)
+      iit->second.second /= iit->second.first;
   }
 
  public:  // Set cover from file.
@@ -730,9 +755,9 @@ class Cover_complex {
   *
   */
   void set_cover_from_file(const std::string& cover_file_name) {
-    int vertex_id = 0;
-    Cover_t cov;
-    std::vector<Cover_t> cov_elts, cov_number;
+    int i = 0;
+    int cov;
+    std::vector<int> cov_elts, cov_number;
     std::ifstream input(cover_file_name);
     std::string line;
     while (std::getline(input, line)) {
@@ -742,17 +767,18 @@ class Cover_complex {
         cov_elts.push_back(cov);
         cov_number.push_back(cov);
         cover_fct[cov] = cov;
-        cover_color[cov].second += func_color[vertex_id];
+        cover_color[cov].second += func_color[i];
         cover_color[cov].first++;
-        cover_back[cov].push_back(vertex_id);
+        cover_back[cov].push_back(i);
       }
-      cover[vertex_id] = cov_elts;
-      vertex_id++;
+      cover[i] = cov_elts;
+      i++;
     }
-    std::vector<Cover_t>::iterator it;
+
     std::sort(cov_number.begin(), cov_number.end());
-    it = std::unique(cov_number.begin(), cov_number.end());
+    std::vector<int>::iterator it = std::unique(cov_number.begin(), cov_number.end());
     cov_number.resize(std::distance(cov_number.begin(), it));
+
     maximal_dim = cov_number.size() - 1;
     for (int i = 0; i <= maximal_dim; i++) cover_color[i].second /= cover_color[i].first;
     cover_name = cover_file_name;
@@ -770,52 +796,25 @@ class Cover_complex {
     voronoi_subsamples.resize(m);
     SampleWithoutReplacement(n, m, voronoi_subsamples);
     if (distances.size() == 0) compute_pairwise_distances(distance);
+    set_graph_weights();
+    Weight_map weight = boost::get(boost::edge_weight, one_skeleton);
+    Index_map index = boost::get(boost::vertex_index, one_skeleton);
     std::vector<double> mindist(n);
     for (int j = 0; j < n; j++) mindist[j] = std::numeric_limits<double>::max();
 
     // Compute the geodesic distances to subsamples with Dijkstra
+    // #pragma omp parallel for
     for (int i = 0; i < m; i++) {
       if (verbose) std::cout << "Computing geodesic distances to seed " << i << "..." << std::endl;
       int seed = voronoi_subsamples[i];
-      std::vector<double> dist(n);
-      std::vector<int> process(n);
-      for (int j = 0; j < n; j++) {
-        dist[j] = std::numeric_limits<double>::max();
-        process[j] = j;
-      }
-      dist[seed] = 0;
-      int curr_size = process.size();
-      int min_point, min_index;
-      double min_dist;
-      std::vector<int> neighbors;
-      int num_neighbors;
-
-      while (curr_size > 0) {
-        min_dist = std::numeric_limits<double>::max();
-        min_index = -1;
-        min_point = -1;
-        for (int j = 0; j < curr_size; j++) {
-          if (dist[process[j]] < min_dist) {
-            min_point = process[j];
-            min_dist = dist[process[j]];
-            min_index = j;
-          }
-        }
-        assert(min_index != -1);
-        process.erase(process.begin() + min_index);
-        assert(min_point != -1);
-        neighbors = adjacency_matrix[min_point];
-        num_neighbors = neighbors.size();
-        for (int j = 0; j < num_neighbors; j++) {
-          double d = dist[min_point] + distances[min_point][neighbors[j]];
-          dist[neighbors[j]] = std::min(dist[neighbors[j]], d);
-        }
-        curr_size = process.size();
-      }
+      std::vector<double> dmap(n);
+      boost::dijkstra_shortest_paths(
+          one_skeleton, vertices[seed],
+          boost::weight_map(weight).distance_map(boost::make_iterator_property_map(dmap.begin(), index)));
 
       for (int j = 0; j < n; j++)
-        if (mindist[j] > dist[j]) {
-          mindist[j] = dist[j];
+        if (mindist[j] > dmap[j]) {
+          mindist[j] = dmap[j];
           if (cover[j].size() == 0)
             cover[j].push_back(i);
           else
@@ -840,7 +839,7 @@ class Cover_complex {
           * @result cover_back(c) vector of IDs of data points.
           *
           */
-  const std::vector<int>& subpopulation(Cover_t c) { return cover_back[c]; }
+  const std::vector<int>& subpopulation(int c) { return cover_back[name2idinv[c]]; }
 
   // *******************************************************************************************************************
   // Visualization.
@@ -854,15 +853,15 @@ class Cover_complex {
    *
    */
   void set_color_from_file(const std::string& color_file_name) {
-    int vertex_id = 0;
+    int i = 0;
     std::ifstream input(color_file_name);
     std::string line;
     double f;
     while (std::getline(input, line)) {
       std::stringstream stream(line);
       stream >> f;
-      func_color.emplace(vertex_id, f);
-      vertex_id++;
+      func_color.emplace(i, f);
+      i++;
     }
     color_name = color_file_name;
   }
@@ -874,7 +873,7 @@ class Cover_complex {
            *
            */
   void set_color_from_coordinate(int k = 0) {
-    for (int i = 0; i < n; i++) func_color.emplace(i, point_cloud[i][k]);
+    for (int i = 0; i < n; i++) func_color[i] = point_cloud[i][k];
     color_name = "coordinate ";
     color_name.append(std::to_string(k));
   }
@@ -886,7 +885,7 @@ class Cover_complex {
    *
    */
   void set_color_from_vector(std::vector<double> color) {
-    for (unsigned int i = 0; i < color.size(); i++) func_color.emplace(i, color[i]);
+    for (unsigned int i = 0; i < color.size(); i++) func_color[i] = color[i];
   }
 
  public:  // Create a .dot file that can be compiled with neato to produce a .pdf file.
@@ -895,26 +894,30 @@ class Cover_complex {
    * of its 1-skeleton in a .pdf file.
    */
   void plot_DOT() {
-    char mapp[11] = "SC.dot";
+    std::string mapp = point_cloud_name + "_sc.dot";
     std::ofstream graphic(mapp);
-    graphic << "graph GIC {" << std::endl;
-    double maxv, minv;
-    maxv = std::numeric_limits<double>::min();
-    minv = std::numeric_limits<double>::max();
-    for (std::map<Cover_t, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end();
-         iit++) {
+
+    double maxv = std::numeric_limits<double>::lowest();
+    double minv = std::numeric_limits<double>::max();
+    for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
       maxv = std::max(maxv, iit->second.second);
       minv = std::min(minv, iit->second.second);
     }
+
     int k = 0;
     std::vector<int> nodes;
     nodes.clear();
-    for (std::map<Cover_t, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end();
-         iit++) {
+
+    graphic << "graph GIC {" << std::endl;
+    int id = 0;
+    for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
       if (iit->second.first > mask) {
         nodes.push_back(iit->first);
-        graphic << iit->first << "[shape=circle fontcolor=black color=black label=\"" << iit->first << ":"
-                << iit->second.first << "\" style=filled fillcolor=\""
+        name2id[iit->first] = id;
+        name2idinv[id] = iit->first;
+        id++;
+        graphic << name2id[iit->first] << "[shape=circle fontcolor=black color=black label=\"" << name2id[iit->first]
+                << ":" << iit->second.first << "\" style=filled fillcolor=\""
                 << (1 - (maxv - iit->second.second) / (maxv - minv)) * 0.6 << ", 1, 1\"]" << std::endl;
         k++;
       }
@@ -924,13 +927,14 @@ class Cover_complex {
     for (int i = 0; i < num_simplices; i++)
       if (simplices[i].size() == 2) {
         if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask) {
-          graphic << "  " << simplices[i][0] << " -- " << simplices[i][1] << " [weight=15];" << std::endl;
+          graphic << "  " << name2id[simplices[i][0]] << " -- " << name2id[simplices[i][1]] << " [weight=15];"
+                  << std::endl;
           ke++;
         }
       }
     graphic << "}";
     graphic.close();
-    std::cout << "SC.dot generated. It can be visualized with e.g. neato." << std::endl;
+    std::cout << mapp << " file generated. It can be visualized with e.g. neato." << std::endl;
   }
 
  public:  // Create a .txt file that can be compiled with KeplerMapper.
@@ -940,8 +944,9 @@ class Cover_complex {
   void write_info() {
     int num_simplices = simplices.size();
     int num_edges = 0;
-    char mapp[11] = "SC.txt";
+    std::string mapp = point_cloud_name + "_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++;
@@ -952,16 +957,21 @@ class Cover_complex {
     graphic << resolution_double << " " << gain << std::endl;
     graphic << cover_color.size() << " " << num_edges << std::endl;
 
-    for (std::map<Cover_t, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end();
-         iit++)
-      graphic << iit->first << " " << iit->second.second << " " << iit->second.first << std::endl;
+    int id = 0;
+    for (std::map<int, std::pair<int, double> >::iterator iit = cover_color.begin(); iit != cover_color.end(); iit++) {
+      graphic << id << " " << iit->second.second << " " << iit->second.first << std::endl;
+      name2id[iit->first] = id;
+      name2idinv[id] = iit->first;
+      id++;
+    }
 
     for (int i = 0; i < num_simplices; i++)
       if (simplices[i].size() == 2)
         if (cover_color[simplices[i][0]].first > mask && cover_color[simplices[i][1]].first > mask)
-          graphic << simplices[i][0] << " " << simplices[i][1] << std::endl;
+          graphic << name2id[simplices[i][0]] << " " << name2id[simplices[i][1]] << std::endl;
     graphic.close();
-    std::cout << "SC.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
+    std::cout << mapp
+              << " generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox."
               << std::endl;
   }
 
@@ -972,14 +982,17 @@ class Cover_complex {
            */
   void plot_OFF() {
     assert(cover_name == "Voronoi");
-    char gic[11] = "SC.off";
-    std::ofstream graphic(gic);
-    graphic << "OFF" << std::endl;
+
     int m = voronoi_subsamples.size();
     int numedges = 0;
     int numfaces = 0;
     std::vector<std::vector<int> > edges, faces;
     int numsimplices = simplices.size();
+
+    std::string mapp = point_cloud_name + "_sc.off";
+    std::ofstream graphic(mapp);
+
+    graphic << "OFF" << std::endl;
     for (int i = 0; i < numsimplices; i++) {
       if (simplices[i].size() == 2) {
         numedges++;
@@ -1004,26 +1017,183 @@ class Cover_complex {
     for (int i = 0; i < numfaces; i++)
       graphic << 3 << " " << faces[i][0] << " " << faces[i][1] << " " << faces[i][2] << std::endl;
     graphic.close();
-    std::cout << "SC.off generated. It can be visualized with e.g. geomview." << std::endl;
+    std::cout << mapp << " generated. It can be visualized with e.g. geomview." << std::endl;
   }
 
   // *******************************************************************************************************************
+  // Extended Persistence Diagrams.
+  // *******************************************************************************************************************
+
+ public:
+  /** \brief Computes the extended persistence diagram of the complex.
+   *
+   */
+  void compute_PD() {
+    Simplex_tree st;
+
+    // Compute max and min
+    double maxf = std::numeric_limits<double>::lowest();
+    double minf = std::numeric_limits<double>::max();
+    for (std::map<int, double>::iterator it = cover_std.begin(); it != cover_std.end(); it++) {
+      maxf = std::max(maxf, it->second);
+      minf = std::min(minf, it->second);
+    }
+
+    for (auto const& simplex : simplices) {
+      // Add a simplex and a cone on it
+      std::vector<int> splx = simplex;
+      splx.push_back(-2);
+      st.insert_simplex_and_subfaces(splx);
+    }
+
+    // Build filtration
+    for (auto simplex : st.complex_simplex_range()) {
+      double filta = std::numeric_limits<double>::lowest();
+      double filts = filta;
+      bool ascending = true;
+      for (auto vertex : st.simplex_vertex_range(simplex)) {
+        if (vertex == -2) {
+          ascending = false;
+          continue;
+        }
+        filta = std::max(-2 + (cover_std[vertex] - minf) / (maxf - minf), filta);
+        filts = std::max(2 - (cover_std[vertex] - minf) / (maxf - minf), filts);
+      }
+      if (ascending)
+        st.assign_filtration(simplex, filta);
+      else
+        st.assign_filtration(simplex, filts);
+    }
+    int magic[] = {-2};
+    st.assign_filtration(st.find(magic), -3);
+
+    // Compute PD
+    st.initialize_filtration();
+    Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Gudhi::persistent_cohomology::Field_Zp> pcoh(st);
+    pcoh.init_coefficients(2);
+    pcoh.compute_persistent_cohomology();
+
+    // Output PD
+    int max_dim = st.dimension();
+    for (int i = 0; i < max_dim; i++) {
+      std::vector<std::pair<double, double> > bars = pcoh.intervals_in_dimension(i);
+      int num_bars = bars.size();
+      if(verbose)  std::cout << num_bars << " interval(s) in dimension " << i << ":" << std::endl;
+      for (int j = 0; j < num_bars; j++) {
+        double birth = bars[j].first;
+        double death = bars[j].second;
+        if (i == 0 && std::isinf(death)) continue;
+        if (birth < 0)
+          birth = minf + (birth + 2) * (maxf - minf);
+        else
+          birth = minf + (2 - birth) * (maxf - minf);
+        if (death < 0)
+          death = minf + (death + 2) * (maxf - minf);
+        else
+          death = minf + (2 - death) * (maxf - minf);
+        PD.push_back(std::pair<double, double>(birth, death));
+        if (verbose) std::cout << "  [" << birth << ", " << death << "]" << std::endl;
+      }
+    }
+  }
+
+ public:
+  /** \brief Computes bootstrapped distances distribution.
+   *
+   * @param[in] N number of bootstrap iterations.
+   *
+   */
+  template <typename SimplicialComplex>
+  void compute_distribution(int N = 100) {
+    if (distribution.size() >= N) {
+      std::cout << "Already done!" << std::endl;
+    } else {
+      for (int i = 0; i < N - distribution.size(); i++) {
+        Cover_complex Cboot;
+        Cboot.n = this->n;
+        std::vector<int> boot(this->n);
+        for (int j = 0; j < this->n; j++) {
+          double u = GetUniform();
+          int id = std::floor(u * (this->n));
+          boot[j] = id;
+          Cboot.point_cloud[j] = this->point_cloud[id];
+          Cboot.func.emplace(j, this->func[id]);
+        }
+        for (int j = 0; j < n; j++) {
+          std::vector<double> dist(n);
+          for (int k = 0; k < n; k++) dist[k] = distances[boot[j]][boot[k]];
+          Cboot.distances.push_back(dist);
+        }
+
+        Cboot.set_graph_from_automatic_rips(Gudhi::Euclidean_distance());
+        Cboot.set_automatic_resolution();
+        Cboot.set_gain();
+        Cboot.set_cover_from_function();
+        Cboot.find_simplices();
+        Cboot.compute_PD();
+
+        distribution.push_back(Gudhi::persistence_diagram::bottleneck_distance(this->PD, Cboot.PD));
+      }
+
+      std::sort(distribution.begin(), distribution.end());
+    }
+  }
+
+ public:
+  /** \brief Computes the bottleneck distance threshold corresponding to a specific confidence level.
+   *
+   * @param[in] alpha Confidence level.
+   *
+   */
+  double compute_distance_from_confidence_level(double alpha) {
+    int N = distribution.size();
+    return distribution[std::floor(alpha * N)];
+  }
+
+ public:
+  /** \brief Computes the confidence level of a specific bottleneck distance threshold.
+   *
+   * @param[in] d Bottleneck distance.
+   *
+   */
+  double compute_confidence_level_from_distance(double d) {
+    int N = distribution.size();
+    for (int i = 0; i < N; i++)
+      if (distribution[i] > d) return i * 1.0 / N;
+  }
+
+ public:
+  /** \brief Computes the p-value, i.e. the opposite of the confidence level of the largest bottleneck
+   * distance preserving the points in the persistence diagram of the output simplicial complex.
+   *
+   */
+  double compute_p_value() {
+    double distancemin = -std::numeric_limits<double>::lowest();
+    int N = PD.size();
+    for (int i = 0; i < N; i++) distancemin = std::min(distancemin, 0.5 * (PD[i].second - PD[i].first));
+    return 1 - compute_confidence_level_from_distance(distancemin);
+  }
+
+  // *******************************************************************************************************************
+  // Computation of simplices.
   // *******************************************************************************************************************
 
  public:
   /** \brief Creates the simplicial complex.
    *
-   * @param[in] complex SimplicialComplexForRips to be created.
+   * @param[in] complex SimplicialComplex to be created.
    *
    */
-  template <typename SimplicialComplexForRips>
-  void create_complex(SimplicialComplexForRips& complex) {
+  template <typename SimplicialComplex>
+  void create_complex(SimplicialComplex& complex) {
     unsigned int dimension = 0;
     for (auto const& simplex : simplices) {
-      complex.insert_simplex_and_subfaces(simplex);
+      int numvert = simplex.size();
+      double filt = std::numeric_limits<double>::lowest();
+      for (int i = 0; i < numvert; i++) filt = std::max(cover_color[simplex[i]].second, filt);
+      complex.insert_simplex_and_subfaces(simplex, filt);
       if (dimension < simplex.size() - 1) dimension = simplex.size() - 1;
     }
-    complex.set_dimension(dimension);
   }
 
  public:
@@ -1036,15 +1206,16 @@ class Cover_complex {
     }
 
     if (type == "Nerve") {
-      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;
+      for(auto& simplex : cover)
+        simplices.push_back(simplex.second);
       std::sort(simplices.begin(), simplices.end());
-      it = std::unique(simplices.begin(), simplices.end());
+      std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
       simplices.resize(std::distance(simplices.begin(), it));
     }
 
     if (type == "GIC") {
+      Index_map index = boost::get(boost::vertex_index, one_skeleton);
+
       if (functional_cover) {
         // Computes the simplices in the GIC by looking at all the edges of the graph and adding the
         // corresponding edges in the GIC if the images of the endpoints belong to consecutive intervals.
@@ -1053,81 +1224,44 @@ class Cover_complex {
           throw std::invalid_argument(
               "the output of this function is correct ONLY if the cover is minimal, i.e. the gain is less than 0.5.");
 
-        int v1, v2;
-
-        // Loop on all points.
-        for (std::map<int, std::vector<Cover_t> >::iterator it = cover.begin(); it != cover.end(); it++) {
-          int vid = it->first;
-          std::vector<int> neighbors = adjacency_matrix[vid];
-          int num_neighb = neighbors.size();
-
-          // Find cover of current point (vid).
-          if (cover[vid].size() == 2)
-            v1 = std::min(cover[vid][0], cover[vid][1]);
-          else
-            v1 = cover[vid][0];
-          std::vector<int> node(1);
-          node[0] = v1;
-          simplices.push_back(node);
-
-          // Loop on neighbors.
-          for (int i = 0; i < num_neighb; i++) {
-            int neighb = neighbors[i];
-
-            // Find cover of neighbor (neighb).
-            if (cover[neighb].size() == 2)
-              v2 = std::max(cover[neighb][0], cover[neighb][1]);
-            else
-              v2 = cover[neighb][0];
-
-            // If neighbor is in next interval, add edge.
-            if (cover_fct[v2] == cover_fct[v1] + 1) {
-              std::vector<int> edge(2);
-              edge[0] = v1;
-              edge[1] = v2;
-              simplices.push_back(edge);
+        // Loop on all edges.
+        boost::graph_traits<Graph>::edge_iterator ei, ei_end;
+        for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei) {
+          int nums = cover[index[boost::source(*ei, one_skeleton)]].size();
+          for (int i = 0; i < nums; i++) {
+            int vs = cover[index[boost::source(*ei, one_skeleton)]][i];
+            int numt = cover[index[boost::target(*ei, one_skeleton)]].size();
+            for (int j = 0; j < numt; j++) {
+              int vt = cover[index[boost::target(*ei, one_skeleton)]][j];
+              if (cover_fct[vs] == cover_fct[vt] + 1 || cover_fct[vt] == cover_fct[vs] + 1) {
+                std::vector<int> edge(2);
+                edge[0] = std::min(vs, vt);
+                edge[1] = std::max(vs, vt);
+                simplices.push_back(edge);
+                goto afterLoop;
+              }
             }
           }
+        afterLoop:;
         }
-        std::vector<std::vector<Cover_t> >::iterator it;
         std::sort(simplices.begin(), simplices.end());
-        it = std::unique(simplices.begin(), simplices.end());
+        std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
         simplices.resize(std::distance(simplices.begin(), it));
 
       } else {
-        // Find IDs of edges to remove
-        std::vector<int> simplex_to_remove;
-        int simplex_id = 0;
-        for (auto simplex : st.complex_simplex_range()) {
-          if (st.dimension(simplex) == 1) {
-            std::vector<std::vector<Cover_t> > comp;
-            for (auto vertex : st.simplex_vertex_range(simplex)) comp.push_back(cover[vertex]);
-            if (comp[0].size() == 1 && comp[0] == comp[1]) simplex_to_remove.push_back(simplex_id);
+        // Find edges to keep
+        Simplex_tree st;
+        boost::graph_traits<Graph>::edge_iterator ei, ei_end;
+        for (boost::tie(ei, ei_end) = boost::edges(one_skeleton); ei != ei_end; ++ei)
+          if (!(cover[index[boost::target(*ei, one_skeleton)]].size() == 1 &&
+                cover[index[boost::target(*ei, one_skeleton)]] == cover[index[boost::source(*ei, one_skeleton)]])) {
+            std::vector<int> edge(2);
+            edge[0] = index[boost::source(*ei, one_skeleton)];
+            edge[1] = index[boost::target(*ei, one_skeleton)];
+            st.insert_simplex_and_subfaces(edge);
           }
-          simplex_id++;
-        }
 
-        // Remove edges
-        if (simplex_to_remove.size() > 1) {
-          int current_id = 1;
-          auto simplex = st.complex_simplex_range().begin();
-          int num_rem = 0;
-          for (int i = 0; i < simplex_id - 1; i++) {
-            int j = i + 1;
-            auto simplex_tmp = simplex;
-            simplex_tmp++;
-            if (j == simplex_to_remove[current_id]) {
-              st.remove_maximal_simplex(*simplex_tmp);
-              current_id++;
-              num_rem++;
-            } else {
-              simplex++;
-            }
-          }
-          simplex = st.complex_simplex_range().begin();
-          for (int i = 0; i < simplex_to_remove[0]; i++) simplex++;
-          st.remove_maximal_simplex(*simplex);
-        }
+        // st.insert_graph(one_skeleton);
 
         // Build the Simplex Tree corresponding to the graph
         st.expansion(maximal_dim);
@@ -1136,23 +1270,21 @@ class Cover_complex {
         simplices.clear();
         for (auto simplex : st.complex_simplex_range()) {
           if (!st.has_children(simplex)) {
-            std::vector<Cover_t> simplx;
+            std::vector<int> simplx;
             for (auto vertex : st.simplex_vertex_range(simplex)) {
               unsigned int sz = cover[vertex].size();
               for (unsigned int i = 0; i < sz; i++) {
                 simplx.push_back(cover[vertex][i]);
               }
             }
-
             std::sort(simplx.begin(), simplx.end());
-            std::vector<Cover_t>::iterator it = std::unique(simplx.begin(), simplx.end());
+            std::vector<int>::iterator it = std::unique(simplx.begin(), simplx.end());
             simplx.resize(std::distance(simplx.begin(), it));
             simplices.push_back(simplx);
           }
         }
-        std::vector<std::vector<Cover_t> >::iterator it;
         std::sort(simplices.begin(), simplices.end());
-        it = std::unique(simplices.begin(), simplices.end());
+        std::vector<std::vector<int> >::iterator it = std::unique(simplices.begin(), simplices.end());
         simplices.resize(std::distance(simplices.begin(), it));
       }
     }
diff --git a/src/Nerve_GIC/test/CMakeLists.txt b/src/Nerve_GIC/test/CMakeLists.txt
index 03fe47ca..c35cdff7 100644
--- a/src/Nerve_GIC/test/CMakeLists.txt
+++ b/src/Nerve_GIC/test/CMakeLists.txt
@@ -1,14 +1,17 @@
 cmake_minimum_required(VERSION 2.6)
 project(Graph_induced_complex_tests)
 
-include(GUDHI_test_coverage)
+if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
+  include(GUDHI_test_coverage)
 
-add_executable ( Nerve_GIC_test_unit test_GIC.cpp )
-target_link_libraries(Nerve_GIC_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
-if (TBB_FOUND)
-  target_link_libraries(Nerve_GIC_test_unit ${TBB_LIBRARIES})
-endif()
+  add_executable ( Nerve_GIC_test_unit test_GIC.cpp )
+  target_link_libraries(Nerve_GIC_test_unit ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(Nerve_GIC_test_unit ${TBB_LIBRARIES})
+  endif()
 
-file(COPY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+  file(COPY data DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
 
-gudhi_add_coverage_test(Nerve_GIC_test_unit)
+  gudhi_add_coverage_test(Nerve_GIC_test_unit)
+
+endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
diff --git a/src/Nerve_GIC/test/test_GIC.cpp b/src/Nerve_GIC/test/test_GIC.cpp
index a8b1e7f7..d633753c 100644
--- a/src/Nerve_GIC/test/test_GIC.cpp
+++ b/src/Nerve_GIC/test/test_GIC.cpp
@@ -24,11 +24,11 @@
 #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>
diff --git a/src/Nerve_GIC/utilities/CMakeLists.txt b/src/Nerve_GIC/utilities/CMakeLists.txt
new file mode 100644
index 00000000..7762c8a0
--- /dev/null
+++ b/src/Nerve_GIC/utilities/CMakeLists.txt
@@ -0,0 +1,24 @@
+cmake_minimum_required(VERSION 2.6)
+project(Nerve_GIC_examples)
+
+if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
+
+  add_executable ( Nerve Nerve.cpp )
+  add_executable ( VoronoiGIC VoronoiGIC.cpp )
+
+  if (TBB_FOUND)
+    target_link_libraries(Nerve ${TBB_LIBRARIES})
+    target_link_libraries(VoronoiGIC ${TBB_LIBRARIES})
+  endif()
+
+  file(COPY KeplerMapperVisuFromTxtFile.py km.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+  # Copy files for not to pollute sources when testing
+  file(COPY "${CMAKE_SOURCE_DIR}/data/points/human.off" DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
+
+  add_test(NAME Nerve_GIC_utilities_nerve COMMAND $<TARGET_FILE:Nerve>
+      "human.off" "2" "10" "0.3")
+
+  add_test(NAME Nerve_GIC_utilities_VoronoiGIC COMMAND $<TARGET_FILE:VoronoiGIC>
+      "human.off" "100")
+
+endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
diff --git a/src/Nerve_GIC/utilities/KeplerMapperVisuFromTxtFile.py b/src/Nerve_GIC/utilities/KeplerMapperVisuFromTxtFile.py
new file mode 100755
index 00000000..c811f610
--- /dev/null
+++ b/src/Nerve_GIC/utilities/KeplerMapperVisuFromTxtFile.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+
+import km
+import numpy as np
+from collections import defaultdict
+import argparse
+
+"""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 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/>.
+"""
+
+__author__ = "Mathieu Carriere"
+__copyright__ = "Copyright (C) 2017 INRIA"
+__license__ = "GPL v3"
+
+parser = argparse.ArgumentParser(description='Creates an html Keppler Mapper '
+                                 'file to visualize a SC.txt file.',
+                                 epilog='Example: '
+                                 './KeplerMapperVisuFromTxtFile.py '
+                                 '-f ../../data/points/human.off_sc.txt'
+                                 '- Constructs an human.off_sc.html file.')
+parser.add_argument("-f", "--file", type=str, required=True)
+
+args = parser.parse_args()
+
+with open(args.file, 'r') as f:
+    network = {}
+    mapper = km.KeplerMapper(verbose=0)
+    data = np.zeros((3,3))
+    projected_data = mapper.fit_transform( data, projection="sum", scaler=None )
+
+    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
+
+    html_output_filename = args.file.rsplit('.', 1)[0] + '.html'
+    mapper.visualize(network, color_function = color, path_html=html_output_filename, 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)
+    message = repr(html_output_filename) + " is generated. You can now use your favorite web browser to visualize it."
+    print(message)
+
+
+    f.close()
diff --git a/src/Nerve_GIC/example/Nerve.cpp b/src/Nerve_GIC/utilities/Nerve.cpp
index 4d5b009b..aefc3874 100644
--- a/src/Nerve_GIC/example/Nerve.cpp
+++ b/src/Nerve_GIC/utilities/Nerve.cpp
@@ -27,8 +27,8 @@
 
 void usage(int nbArgs, char *const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [--v] \n";
-  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 --v \n";
+  std::cerr << "Usage: " << progName << " filename.off coordinate resolution gain [-v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 2 10 0.3 -v \n";
   exit(-1);  // ----- >>
 }
 
@@ -72,6 +72,7 @@ int main(int argc, char **argv) {
 
     Gudhi::Simplex_tree<> stree;
     SC.create_complex(stree);
+    SC.compute_PD();
 
     // ----------------------------------------------------------------------------
     // Display information about the graph induced complex
diff --git a/src/Nerve_GIC/utilities/Nerve.txt b/src/Nerve_GIC/utilities/Nerve.txt
new file mode 100644
index 00000000..839ff45e
--- /dev/null
+++ b/src/Nerve_GIC/utilities/Nerve.txt
@@ -0,0 +1,63 @@
+Min function value = -0.979672 and Max function value = 0.816414
+Interval 0 = [-0.979672, -0.761576]
+Interval 1 = [-0.838551, -0.581967]
+Interval 2 = [-0.658942, -0.402359]
+Interval 3 = [-0.479334, -0.22275]
+Interval 4 = [-0.299725, -0.0431415]
+Interval 5 = [-0.120117, 0.136467]
+Interval 6 = [0.059492, 0.316076]
+Interval 7 = [0.239101, 0.495684]
+Interval 8 = [0.418709, 0.675293]
+Interval 9 = [0.598318, 0.816414]
+Computing preimages...
+Computing connected components...
+.txt generated. It can be visualized with e.g. python KeplerMapperVisuFromTxtFile.py and firefox.
+5 interval(s) in dimension 0:
+  [-0.909111, 0.00817529]
+  [-0.171433, 0.367392]
+  [-0.171433, 0.367392]
+  [-0.909111, 0.745853]
+0 interval(s) in dimension 1:
+Nerve is of dimension 1 - 41 simplices - 21 vertices.
+Iterator on Nerve simplices
+1 
+0 
+4 
+4 0 
+2 
+2 1 
+8 
+8 2 
+5 
+5 4 
+9 
+9 8 
+13 
+13 5 
+14 
+14 9 
+19 
+19 13 
+25 
+32 
+20 
+32 20 
+33 
+33 25 
+26 
+26 14 
+26 19 
+42 
+42 26 
+34 
+34 33 
+27 
+27 20 
+35 
+35 27 
+35 34 
+42 35 
+44 
+44 35 
+54 
+54 44 
\ No newline at end of file
diff --git a/src/Nerve_GIC/example/VoronoiGIC.cpp b/src/Nerve_GIC/utilities/VoronoiGIC.cpp
index 32431cc2..54bb871e 100644
--- a/src/Nerve_GIC/example/VoronoiGIC.cpp
+++ b/src/Nerve_GIC/utilities/VoronoiGIC.cpp
@@ -27,8 +27,8 @@
 
 void usage(int nbArgs, char *const progName) {
   std::cerr << "Error: Number of arguments (" << nbArgs << ") is not correct\n";
-  std::cerr << "Usage: " << progName << " filename.off N [--v] \n";
-  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 100 --v \n";
+  std::cerr << "Usage: " << progName << " filename.off N [-v] \n";
+  std::cerr << "       i.e.: " << progName << " ../../data/points/human.off 100 -v \n";
   exit(-1);  // ----- >>
 }
 
diff --git a/src/Nerve_GIC/utilities/covercomplex.md b/src/Nerve_GIC/utilities/covercomplex.md
new file mode 100644
index 00000000..f33cb2e0
--- /dev/null
+++ b/src/Nerve_GIC/utilities/covercomplex.md
@@ -0,0 +1,62 @@
+
+
+# Cover complex #
+
+
+## Nerve ##
+This program builds the Nerve of a point cloud sampled on an OFF file.
+The cover C comes from the preimages of intervals covering a coordinate function,
+which are then refined into their connected components using the triangulation of the .OFF file.
+
+The program also writes a file SC.txt.
+The first three lines in this file are the location of the input point cloud and the function used to compute the cover.
+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.
+
+**Usage**
+
+`Nerve <OFF input file> coordinate resolution gain [-v]`
+
+where
+
+* `coordinate` is the coordinate function to cover
+* `resolution` is the number of the intervals
+* `gain` is the gain for each interval
+* `-v` is optional, it activates verbose mode.
+
+**Example**
+
+`Nerve ../../data/points/human.off 2 10 0.3`
+
+* 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).
+
+`python KeplerMapperVisuFromTxtFile.py -f ../../data/points/human.off_sc.txt`
+
+* Constructs `human.off_sc.html` file. You can now use your favorite web browser to visualize it.
+
+## VoronoiGIC ##
+
+This util builds the Graph Induced Complex (GIC) of a point cloud.
+It subsamples *N* 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 program also writes a file `*_sc.off`, that is an OFF file that can be visualized with GeomView.
+
+**Usage**
+
+`VoroniGIC <OFF input file> samples_number [-v]`
+
+where
+
+* `samples_number` is the number of samples to take from the point cloud
+* `-v` is optional, it activates verbose mode.
+
+**Example**
+
+`VoroniGIC ../../data/points/human.off 700`
+
+* Builds the Voronoi Graph Induced Complex with 700 subsamples from `human.off` file.
+`../../data/points/human_sc.off` can be visualized with GeomView.
+
diff --git a/src/Nerve_GIC/example/km.py b/src/Nerve_GIC/utilities/km.py
index 53024aab..53024aab 100755
--- a/src/Nerve_GIC/example/km.py
+++ b/src/Nerve_GIC/utilities/km.py
diff --git a/src/Nerve_GIC/example/km.py.COPYRIGHT b/src/Nerve_GIC/utilities/km.py.COPYRIGHT
index bef7b121..bef7b121 100644
--- a/src/Nerve_GIC/example/km.py.COPYRIGHT
+++ b/src/Nerve_GIC/utilities/km.py.COPYRIGHT