GUDHI 2.1.0 as released by upstream in a tarball.upstream/2.1.0

43 files changed, 1114 insertions, 2447 deletions

diff --git a/example/Bitmap_cubical_complex/Bitmap_cubical_complex.cpp b/example/Bitmap_cubical_complex/Bitmap_cubical_complex.cpp
deleted file mode 100644
index 67735ba1..00000000
--- a/example/Bitmap_cubical_complex/Bitmap_cubical_complex.cpp
+++ /dev/null
@@ -1,80 +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):       Pawel Dlotko
- *
- *    Copyright (C) 2015  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/>.
- */
-
-
-#include <gudhi/reader_utils.h>
-#include <gudhi/Bitmap_cubical_complex.h>
-#include <gudhi/Persistent_cohomology.h>
-
-// standard stuff
-#include <iostream>
-#include <string>
-#include <vector>
-#include <cstddef>
-
-int main(int argc, char** argv) {
-  std::cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " <<
-      "complexes provided in text files in Perseus style (the only numbered in the first line is a dimension D of a" <<
-      "bitmap. In the lines I between 2 and D+1 there are numbers of top dimensional cells in the direction I. Let " <<
-      "N denote product of the numbers in the lines between 2 and D. In the lines D+2 to D+2+N there are " <<
-      "filtrations of top dimensional cells. We assume that the cells are in the lexicographical order. See " <<
-      "CubicalOneSphere.txt or CubicalTwoSphere.txt for example.\n" << std::endl;
-
-  if (argc != 2) {
-    std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " <<
-        "the input. The program will now terminate.\n";
-    return 1;
-  }
-
-  typedef Gudhi::cubical_complex::Bitmap_cubical_complex_base<double> Bitmap_cubical_complex_base;
-  typedef Gudhi::cubical_complex::Bitmap_cubical_complex<Bitmap_cubical_complex_base> Bitmap_cubical_complex;
-  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
-  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
-
-  Bitmap_cubical_complex b(argv[1]);
-
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(b);
-  int p = 2;
-  double min_persistence = 0;
-
-  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  std::string output_file_name(argv[1]);
-  output_file_name += "_persistence";
-
-  std::size_t last_in_path = output_file_name.find_last_of("/\\");
-
-  if (last_in_path != std::string::npos) {
-    output_file_name = output_file_name.substr(last_in_path+1);
-  }
-
-  std::ofstream out(output_file_name.c_str());
-  pcoh.output_diagram(out);
-  out.close();
-
-  std::cout << "Result in file: " << output_file_name << "\n";
-
-  return 0;
-}
-
diff --git a/example/Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp b/example/Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
deleted file mode 100644
index 122160a2..00000000
--- a/example/Bitmap_cubical_complex/Bitmap_cubical_complex_periodic_boundary_conditions.cpp
+++ /dev/null
@@ -1,82 +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):       Pawel Dlotko
- *
- *    Copyright (C) 2015  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/>.
- */
-
-
-#include <gudhi/reader_utils.h>
-#include <gudhi/Bitmap_cubical_complex.h>
-#include <gudhi/Bitmap_cubical_complex_periodic_boundary_conditions_base.h>
-#include <gudhi/Persistent_cohomology.h>
-
-// standard stuff
-#include <iostream>
-#include <sstream>
-#include <vector>
-#include <string>
-
-int main(int argc, char** argv) {
-  std::cout << "This program computes persistent homology, by using " <<
-      "Bitmap_cubical_complex_periodic_boundary_conditions class, of cubical complexes provided in text files in " <<
-      "Perseus style (the only numbered in the first line is a dimension D of a bitmap. In the lines I between 2 " <<
-      "and D+1 there are numbers of top dimensional cells in the direction I. Let N denote product of the numbers " <<
-      "in the lines between 2 and D. In the lines D+2 to D+2+N there are filtrations of top dimensional cells. We " <<
-      "assume that the cells are in the lexicographical order. See CubicalOneSphere.txt or CubicalTwoSphere.txt for" <<
-      " example.\n" << std::endl;
-
-  if (argc != 2) {
-    std::cerr << "Wrong number of parameters. Please provide the name of a file with a Perseus style bitmap at " <<
-        "the input. The program will now terminate.\n";
-    return 1;
-  }
-
-  typedef Gudhi::cubical_complex::Bitmap_cubical_complex_periodic_boundary_conditions_base<double> Bitmap_base;
-  typedef Gudhi::cubical_complex::Bitmap_cubical_complex< Bitmap_base > Bitmap_cubical_complex;
-
-  Bitmap_cubical_complex b(argv[1]);
-
-  typedef Gudhi::persistent_cohomology::Field_Zp Field_Zp;
-  typedef Gudhi::persistent_cohomology::Persistent_cohomology<Bitmap_cubical_complex, Field_Zp> Persistent_cohomology;
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(b, true);
-
-  int p = 2;
-  double min_persistence = 0;
-  pcoh.init_coefficients(p);  // initializes the coefficient field for homology
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  std::string output_file_name(argv[1]);
-  output_file_name += "_persistence";
-
-  std::size_t last_in_path = output_file_name.find_last_of("/\\");
-
-  if (last_in_path != std::string::npos) {
-    output_file_name = output_file_name.substr(last_in_path+1);
-  }
-
-  std::ofstream out(output_file_name.c_str());
-  pcoh.output_diagram(out);
-  out.close();
-
-  std::cout << "Result in file: " << output_file_name << "\n";
-
-  return 0;
-}
-
diff --git a/example/Bitmap_cubical_complex/CMakeLists.txt b/example/Bitmap_cubical_complex/CMakeLists.txt
index a0401619..99304aa4 100644
--- a/example/Bitmap_cubical_complex/CMakeLists.txt
+++ b/example/Bitmap_cubical_complex/CMakeLists.txt
@@ -1,17 +1,6 @@
 cmake_minimum_required(VERSION 2.6)
 project(Bitmap_cubical_complex_examples)
 
-add_executable ( Bitmap_cubical_complex Bitmap_cubical_complex.cpp )
-if (TBB_FOUND)
-  target_link_libraries(Bitmap_cubical_complex ${TBB_LIBRARIES})
-endif()
-
-add_test(NAME Bitmap_cubical_complex_example_persistence_one_sphere COMMAND $<TARGET_FILE:Bitmap_cubical_complex>
-    "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalOneSphere.txt")
-
-add_test(NAME Bitmap_cubical_complex_example_persistence_two_sphere COMMAND $<TARGET_FILE:Bitmap_cubical_complex>
-    "${CMAKE_SOURCE_DIR}/data/bitmap/CubicalTwoSphere.txt")
-
 add_executable ( Random_bitmap_cubical_complex Random_bitmap_cubical_complex.cpp )
 if (TBB_FOUND)
   target_link_libraries(Random_bitmap_cubical_complex ${TBB_LIBRARIES})
@@ -19,19 +8,4 @@ endif()
 add_test(NAME Bitmap_cubical_complex_example_random COMMAND $<TARGET_FILE:Random_bitmap_cubical_complex>
     "2" "100" "100")
 
-add_executable ( Bitmap_cubical_complex_periodic_boundary_conditions Bitmap_cubical_complex_periodic_boundary_conditions.cpp )
-if (TBB_FOUND)
-  target_link_libraries(Bitmap_cubical_complex_periodic_boundary_conditions ${TBB_LIBRARIES})
-endif()
-
-add_test(NAME Bitmap_cubical_complex_example_periodic_boundary_conditions_2d_torus
-    COMMAND $<TARGET_FILE:Bitmap_cubical_complex_periodic_boundary_conditions>
-    "${CMAKE_SOURCE_DIR}/data/bitmap/2d_torus.txt")
-
-add_test(NAME Bitmap_cubical_complex_example_periodic_boundary_conditions_3d_torus
-    COMMAND $<TARGET_FILE:Bitmap_cubical_complex_periodic_boundary_conditions>
-    "${CMAKE_SOURCE_DIR}/data/bitmap/3d_torus.txt")
-
-install(TARGETS Bitmap_cubical_complex DESTINATION bin)
 install(TARGETS Random_bitmap_cubical_complex DESTINATION bin)
-install(TARGETS Bitmap_cubical_complex_periodic_boundary_conditions DESTINATION bin)
diff --git a/example/Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp b/example/Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp
index 16ad65a0..f70558f2 100644
--- a/example/Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp
+++ b/example/Bitmap_cubical_complex/Random_bitmap_cubical_complex.cpp
@@ -20,7 +20,6 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-
 // for persistence algorithm
 #include <gudhi/reader_utils.h>
 #include <gudhi/Bitmap_cubical_complex.h>
@@ -34,10 +33,11 @@
 int main(int argc, char** argv) {
   srand(time(0));
 
-  std::cout << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical " <<
-      "complexes. The first parameter of the program is the dimension D of the bitmap. The next D parameters are " <<
-      "number of top dimensional cubes in each dimension of the bitmap. The program will create random cubical " <<
-      "complex of that sizes and compute persistent homology of it." << std::endl;
+  std::cout
+      << "This program computes persistent homology, by using bitmap_cubical_complex class, of cubical "
+      << "complexes. The first parameter of the program is the dimension D of the bitmap. The next D parameters are "
+      << "number of top dimensional cubes in each dimension of the bitmap. The program will create random cubical "
+      << "complex of that sizes and compute persistent homology of it." << std::endl;
 
   int p = 2;
   double min_persistence = 0;
@@ -47,16 +47,16 @@ int main(int argc, char** argv) {
     return 1;
   }
 
-  size_t dimensionOfBitmap = (size_t) atoi(argv[1]);
-  std::vector< unsigned > sizes;
+  size_t dimensionOfBitmap = (size_t)atoi(argv[1]);
+  std::vector<unsigned> sizes;
   size_t multipliers = 1;
   for (size_t dim = 0; dim != dimensionOfBitmap; ++dim) {
-    unsigned sizeInThisDimension = (unsigned) atoi(argv[2 + dim]);
+    unsigned sizeInThisDimension = (unsigned)atoi(argv[2 + dim]);
     sizes.push_back(sizeInThisDimension);
     multipliers *= sizeInThisDimension;
   }
 
-  std::vector< double > data;
+  std::vector<double> data;
   for (size_t i = 0; i != multipliers; ++i) {
     data.push_back(rand() / static_cast<double>(RAND_MAX));
   }
@@ -80,4 +80,3 @@ int main(int argc, char** argv) {
 
   return 0;
 }
-
diff --git a/example/Bottleneck_distance/CMakeLists.txt b/example/Bottleneck_distance/CMakeLists.txt
index eac617db..6095d6eb 100644
--- a/example/Bottleneck_distance/CMakeLists.txt
+++ b/example/Bottleneck_distance/CMakeLists.txt
@@ -2,37 +2,21 @@ cmake_minimum_required(VERSION 2.6)
 project(Bottleneck_distance_examples)
 
 if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
-  add_executable (bottleneck_read_file_example bottleneck_read_file_example.cpp)
   add_executable (bottleneck_basic_example bottleneck_basic_example.cpp)
+  add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp)
+  target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_PROGRAM_OPTIONS_LIBRARY})
 
   if (TBB_FOUND)
-    target_link_libraries(bottleneck_read_file_example ${TBB_LIBRARIES})
+    target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES})
     target_link_libraries(bottleneck_basic_example ${TBB_LIBRARIES})
   endif(TBB_FOUND)
-
+  
   add_test(NAME Bottleneck_distance_example_basic COMMAND $<TARGET_FILE:bottleneck_basic_example>)
-
-  add_test(NAME Bottleneck_read_file_example
-      COMMAND $<TARGET_FILE:bottleneck_read_file_example>
-      "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers" "${CMAKE_SOURCE_DIR}/data/persistence_diagram/second.pers")
-
-  install(TARGETS bottleneck_read_file_example DESTINATION bin)
-  install(TARGETS bottleneck_basic_example DESTINATION bin)
-
-endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
-
-# Eigen3 and CGAL > 4.7.0 is required for alpha complex
-# CGAL > 4.8.1 is required for bottleneck distance =>
-if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
-  add_executable (alpha_rips_persistence_bottleneck_distance alpha_rips_persistence_bottleneck_distance.cpp)
-  target_link_libraries(alpha_rips_persistence_bottleneck_distance ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
   add_test(NAME Bottleneck_distance_example_alpha_rips_persistence_bottleneck
-           COMMAND $<TARGET_FILE:alpha_rips_persistence_bottleneck_distance>
-           "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.15" "-m" "0.12" "-d" "3" "-p" "3")
+      COMMAND $<TARGET_FILE:alpha_rips_persistence_bottleneck_distance>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.15" "-m" "0.12" "-d" "3" "-p" "3")
 
+  install(TARGETS bottleneck_basic_example DESTINATION bin)
   install(TARGETS alpha_rips_persistence_bottleneck_distance DESTINATION bin)
-  if (TBB_FOUND)
-    target_link_libraries(alpha_rips_persistence_bottleneck_distance ${TBB_LIBRARIES})
-  endif(TBB_FOUND)
-endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
+  
+endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
diff --git a/example/Bottleneck_distance/README b/example/Bottleneck_distance/README
new file mode 100644
index 00000000..01bcd74a
--- /dev/null
+++ b/example/Bottleneck_distance/README
@@ -0,0 +1,19 @@
+# Bottleneck_distance #

+

+## `alpha_rips_persistence_bottleneck_distance` ##

+This program computes the persistent homology with coefficient field Z/pZ of a Rips complex defined on a set of input points. The output diagram contains one bar per line, written with the convention:

+

+`p dim birth death`

+

+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients.

+

+Usage:

+`alpha_rips_persistence_bottleneck_distance [options] <OFF input file>`

+

+Allowed options:

+

+* `-h [ --help ]` Produce help message

+* `-r [ --max-edge-length ]` (default = inf) Maximal length of an edge for the Rips complex construction.`

+* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute.`

+* `-p [ --field-charac ]` (default = 11)     Characteristic p of the coefficient field Z/pZ for computing homology.

+* `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals.

diff --git a/example/Bottleneck_distance/bottleneck_read_file_example.cpp b/example/Bottleneck_distance/bottleneck_read_file_example.cpp
deleted file mode 100644
index 24d73c57..00000000
--- a/example/Bottleneck_distance/bottleneck_read_file_example.cpp
+++ /dev/null
@@ -1,50 +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.
- *
- *    Authors:       Francois Godi, small modifications by Pawel Dlotko
- *
- *    Copyright (C) 2015  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/Bottleneck.h>
-#include <gudhi/reader_utils.h>
-#include <iostream>
-#include <vector>
-#include <utility>  // for pair
-#include <string>
-#include <limits>  // for numeric_limits
-
-int main(int argc, char** argv) {
-  if (argc < 3) {
-    std::cout << "To run this program please provide as an input two files with persistence diagrams. Each file" <<
-        " should contain a birth-death pair per line. Third, optional parameter is an error bound on a bottleneck" <<
-        " distance (set by default to the smallest positive double value). If you set the error bound to 0, be" <<
-        " aware this version is exact but expensive. The program will now terminate \n";
-    return -1;
-  }
-  std::vector<std::pair<double, double>> diag1 = Gudhi::read_persistence_intervals_in_dimension(argv[1]);
-  std::vector<std::pair<double, double>> diag2 = Gudhi::read_persistence_intervals_in_dimension(argv[2]);
-
-  double tolerance = std::numeric_limits<double>::min();
-  if (argc == 4) {
-    tolerance = atof(argv[3]);
-  }
-  double b = Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, tolerance);
-  std::cout << "The distance between the diagrams is : " << b << ". The tolerance is : " << tolerance << std::endl;
-
-  return 0;
-}
diff --git a/example/Contraction/CMakeLists.txt b/example/Contraction/CMakeLists.txt
index 83594c0e..a92d1685 100644
--- a/example/Contraction/CMakeLists.txt
+++ b/example/Contraction/CMakeLists.txt
@@ -1,9 +1,10 @@
 cmake_minimum_required(VERSION 2.6)
 project(Contraction_examples)
 
-
 add_executable(RipsContraction Rips_contraction.cpp)
+
 add_executable(GarlandHeckbert Garland_heckbert.cpp)
+target_link_libraries(GarlandHeckbert ${Boost_TIMER_LIBRARY})
 
 add_test(NAME Contraction_example_tore3D_0.2 COMMAND $<TARGET_FILE:RipsContraction>
     "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "0.2")
diff --git a/example/Contraction/Garland_heckbert.cpp b/example/Contraction/Garland_heckbert.cpp
index f0cde95e..2b0dc973 100644
--- a/example/Contraction/Garland_heckbert.cpp
+++ b/example/Contraction/Garland_heckbert.cpp
@@ -29,7 +29,6 @@
 #include <gudhi/Edge_contraction.h>
 #include <gudhi/Skeleton_blocker.h>
 #include <gudhi/Off_reader.h>
-#include <gudhi/Clock.h>
 
 #include <iostream>
 
@@ -165,8 +164,6 @@ int main(int argc, char *argv[]) {
 
   int num_contractions = atoi(argv[3]);
 
-  Gudhi::Clock contraction_chrono("Time to simplify and enumerate simplices");
-
   // constructs the contractor object with Garland Heckbert policies.
   Complex_contractor contractor(complex,
                                 new GH_cost(complex),
@@ -182,8 +179,6 @@ int main(int argc, char *argv[]) {
       complex.num_edges() << " edges and " <<
       complex.num_triangles() << " triangles." << std::endl;
 
-  std::cout << contraction_chrono;
-
   // write simplified complex
   Gudhi::skeleton_blocker::Skeleton_blocker_off_writer<Complex> off_writer(argv[2], complex);
 
diff --git a/example/Nerve_GIC/CMakeLists.txt b/example/Nerve_GIC/CMakeLists.txt
new file mode 100644
index 00000000..f2626927
--- /dev/null
+++ b/example/Nerve_GIC/CMakeLists.txt
@@ -0,0 +1,26 @@
+cmake_minimum_required(VERSION 2.6)
+project(Nerve_GIC_examples)
+
+if (NOT CGAL_VERSION VERSION_LESS 4.8.1)
+
+  add_executable ( CoordGIC CoordGIC.cpp )
+  add_executable ( FuncGIC FuncGIC.cpp )
+
+  if (TBB_FOUND)
+    target_link_libraries(CoordGIC ${TBB_LIBRARIES})
+    target_link_libraries(FuncGIC ${TBB_LIBRARIES})
+  endif()
+
+  # 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_CoordGIC COMMAND $<TARGET_FILE:CoordGIC>
+      "tore3D_1307.off" "0")
+
+  add_test(NAME Nerve_GIC_example_FuncGIC COMMAND $<TARGET_FILE:FuncGIC>
+      "lucky_cat.off"
+      "lucky_cat_PCA1")
+
+endif (NOT CGAL_VERSION VERSION_LESS 4.8.1)
diff --git a/example/Nerve_GIC/CoordGIC.cpp b/example/Nerve_GIC/CoordGIC.cpp
new file mode 100644
index 00000000..c92cf235
--- /dev/null
+++ b/example/Nerve_GIC/CoordGIC.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
+ *
+ *    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 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]);
+
+  using Point = std::vector<float>;
+
+  std::string off_file_name(argv[1]);
+  int coord = atoi(argv[2]);
+  bool verb = 0;
+  if (argc == 4) verb = 1;
+
+  // -----------------------------------------
+  // Init of a functional GIC from an OFF file
+  // -----------------------------------------
+
+  Gudhi::cover_complex::Cover_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_type("GIC");
+
+    GIC.set_color_from_coordinate(coord);
+    GIC.set_function_from_coordinate(coord);
+
+    GIC.set_graph_from_automatic_rips(Gudhi::Euclidean_distance());
+    GIC.set_automatic_resolution();
+    GIC.set_gain();
+    GIC.set_cover_from_function();
+
+    GIC.find_simplices();
+
+    GIC.plot_DOT();
+
+    Gudhi::Simplex_tree<> stree;
+    GIC.create_complex(stree);
+
+    // --------------------------------------------
+    // Display information about the functional GIC
+    // --------------------------------------------
+
+    if (verb) {
+      std::cout << "Functional GIC is of dimension " << stree.dimension() << " - " << stree.num_simplices()
+                << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
+
+      std::cout << "Iterator on functional GIC 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/example/Nerve_GIC/FuncGIC.cpp b/example/Nerve_GIC/FuncGIC.cpp
new file mode 100644
index 00000000..cb0f0d63
--- /dev/null
+++ b/example/Nerve_GIC/FuncGIC.cpp
@@ -0,0 +1,94 @@
+/*    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 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]);
+
+  using Point = std::vector<float>;
+
+  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 functional GIC from an OFF file
+  // -----------------------------------------
+
+  Gudhi::cover_complex::Cover_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_type("GIC");
+
+    GIC.set_color_from_file(func_file_name);
+    GIC.set_function_from_file(func_file_name);
+
+    GIC.set_graph_from_automatic_rips(Gudhi::Euclidean_distance());
+    GIC.set_automatic_resolution();
+    GIC.set_gain();
+    GIC.set_cover_from_function();
+
+    GIC.find_simplices();
+
+    GIC.plot_DOT();
+
+    Gudhi::Simplex_tree<> stree;
+    GIC.create_complex(stree);
+
+    // --------------------------------------------
+    // Display information about the functional GIC
+    // --------------------------------------------
+
+    if (verb) {
+      std::cout << "Functional GIC is of dimension " << stree.dimension() << " - " << stree.num_simplices()
+                << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
+
+      std::cout << "Iterator on functional GIC 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/example/Persistence_representations/CMakeLists.txt b/example/Persistence_representations/CMakeLists.txt
new file mode 100644
index 00000000..eb3258f8
--- /dev/null
+++ b/example/Persistence_representations/CMakeLists.txt
@@ -0,0 +1,27 @@
+cmake_minimum_required(VERSION 2.6)
+project(Persistence_representations_example)
+
+add_executable ( Persistence_representations_example_landscape_on_grid persistence_landscape_on_grid.cpp )
+add_test(NAME Persistence_representations_example_landscape_on_grid
+    COMMAND $<TARGET_FILE:Persistence_representations_example_landscape_on_grid>)
+
+add_executable ( Persistence_representations_example_landscape persistence_landscape.cpp )
+add_test(NAME Persistence_representations_example_landscape
+    COMMAND $<TARGET_FILE:Persistence_representations_example_landscape>)
+
+add_executable ( Persistence_representations_example_intervals persistence_intervals.cpp )
+add_test(NAME Persistence_representations_example_intervals
+    COMMAND $<TARGET_FILE:Persistence_representations_example_intervals>
+    "${CMAKE_SOURCE_DIR}/data/persistence_diagram/first.pers")
+
+add_executable ( Persistence_representations_example_vectors persistence_vectors.cpp )
+add_test(NAME Persistence_representations_example_vectors
+    COMMAND $<TARGET_FILE:Persistence_representations_example_vectors>)
+
+add_executable ( Persistence_representations_example_heat_maps persistence_heat_maps.cpp )
+add_test(NAME Persistence_representations_example_heat_maps
+    COMMAND $<TARGET_FILE:Persistence_representations_example_heat_maps>)
+
+
+
+
diff --git a/example/Persistence_representations/persistence_heat_maps.cpp b/example/Persistence_representations/persistence_heat_maps.cpp
new file mode 100644
index 00000000..2a472ac6
--- /dev/null
+++ b/example/Persistence_representations/persistence_heat_maps.cpp
@@ -0,0 +1,80 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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/>.
+ */
+
+#include <gudhi/Persistence_heat_maps.h>
+
+#include <iostream>
+#include <vector>
+#include <utility>
+
+using constant_scaling_function = Gudhi::Persistence_representations::constant_scaling_function;
+using Persistence_heat_maps = Gudhi::Persistence_representations::Persistence_heat_maps<constant_scaling_function>;
+
+int main(int argc, char** argv) {
+  // create two simple vectors with birth--death pairs:
+
+  std::vector<std::pair<double, double> > persistence1;
+  std::vector<std::pair<double, double> > persistence2;
+
+  persistence1.push_back(std::make_pair(1, 2));
+  persistence1.push_back(std::make_pair(6, 8));
+  persistence1.push_back(std::make_pair(0, 4));
+  persistence1.push_back(std::make_pair(3, 8));
+
+  persistence2.push_back(std::make_pair(2, 9));
+  persistence2.push_back(std::make_pair(1, 6));
+  persistence2.push_back(std::make_pair(3, 5));
+  persistence2.push_back(std::make_pair(6, 10));
+
+  // over here we define a function we sill put on a top on every birth--death pair in the persistence interval. It can
+  // be anything. Over here we will use standard Gaussian
+  std::vector<std::vector<double> > filter = Gudhi::Persistence_representations::create_Gaussian_filter(5, 1);
+
+  // creating two heat maps.
+  Persistence_heat_maps hm1(persistence1, filter, false, 20, 0, 11);
+  Persistence_heat_maps hm2(persistence2, filter, false, 20, 0, 11);
+
+  std::vector<Persistence_heat_maps*> vector_of_maps;
+  vector_of_maps.push_back(&hm1);
+  vector_of_maps.push_back(&hm2);
+
+  // compute median/mean of a vector of heat maps:
+  Persistence_heat_maps mean;
+  mean.compute_mean(vector_of_maps);
+  Persistence_heat_maps median;
+  median.compute_median(vector_of_maps);
+
+  // to compute L^1 distance between hm1 and hm2:
+  std::cout << "The L^1 distance is : " << hm1.distance(hm2, 1) << std::endl;
+
+  // to average of hm1 and hm2:
+  std::vector<Persistence_heat_maps*> to_average;
+  to_average.push_back(&hm1);
+  to_average.push_back(&hm2);
+  Persistence_heat_maps av;
+  av.compute_average(to_average);
+
+  // to compute scalar product of hm1 and hm2:
+  std::cout << "Scalar product is : " << hm1.compute_scalar_product(hm2) << std::endl;
+
+  return 0;
+}
diff --git a/example/Persistence_representations/persistence_intervals.cpp b/example/Persistence_representations/persistence_intervals.cpp
new file mode 100644
index 00000000..c1ceb458
--- /dev/null
+++ b/example/Persistence_representations/persistence_intervals.cpp
@@ -0,0 +1,89 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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/>.
+ */
+
+#include <gudhi/Persistence_intervals.h>
+
+#include <iostream>
+#include <utility>
+#include <vector>
+
+using Persistence_intervals = Gudhi::Persistence_representations::Persistence_intervals;
+
+int main(int argc, char** argv) {
+  if (argc != 2) {
+    std::cout << "To run this program, please provide the name of a file with persistence diagram \n";
+    return 1;
+  }
+
+  Persistence_intervals p(argv[1]);
+  std::pair<double, double> min_max_ = p.get_x_range();
+  std::cout << "Birth-death range : " << min_max_.first << " " << min_max_.second << std::endl;
+
+  std::vector<double> dominant_ten_intervals_length = p.length_of_dominant_intervals(10);
+  std::cout << "Length of ten dominant intervals : " << std::endl;
+  for (size_t i = 0; i != dominant_ten_intervals_length.size(); ++i) {
+    std::cout << dominant_ten_intervals_length[i] << std::endl;
+  }
+
+  std::vector<std::pair<double, double> > ten_dominant_intervals = p.dominant_intervals(10);
+  std::cout << "Here are the dominant intervals : " << std::endl;
+  for (size_t i = 0; i != ten_dominant_intervals.size(); ++i) {
+    std::cout << "( " << ten_dominant_intervals[i].first << "," << ten_dominant_intervals[i].second << std::endl;
+  }
+
+  std::vector<size_t> histogram = p.histogram_of_lengths(10);
+  std::cout << "Here is the histogram of barcode's length : " << std::endl;
+  for (size_t i = 0; i != histogram.size(); ++i) {
+    std::cout << histogram[i] << " ";
+  }
+  std::cout << std::endl;
+
+  std::vector<size_t> cumulative_histogram = p.cumulative_histogram_of_lengths(10);
+  std::cout << "Cumulative histogram : " << std::endl;
+  for (size_t i = 0; i != cumulative_histogram.size(); ++i) {
+    std::cout << cumulative_histogram[i] << " ";
+  }
+  std::cout << std::endl;
+
+  std::vector<double> char_funct_diag = p.characteristic_function_of_diagram(min_max_.first, min_max_.second);
+  std::cout << "Characteristic function of diagram : " << std::endl;
+  for (size_t i = 0; i != char_funct_diag.size(); ++i) {
+    std::cout << char_funct_diag[i] << " ";
+  }
+  std::cout << std::endl;
+
+  std::vector<double> cumul_char_funct_diag =
+      p.cumulative_characteristic_function_of_diagram(min_max_.first, min_max_.second);
+  std::cout << "Cumulative characteristic function of diagram : " << std::endl;
+  for (size_t i = 0; i != cumul_char_funct_diag.size(); ++i) {
+    std::cout << cumul_char_funct_diag[i] << " ";
+  }
+  std::cout << std::endl;
+
+  std::cout << "Persistence Betti numbers \n";
+  std::vector<std::pair<double, size_t> > pbns = p.compute_persistent_betti_numbers();
+  for (size_t i = 0; i != pbns.size(); ++i) {
+    std::cout << pbns[i].first << " " << pbns[i].second << std::endl;
+  }
+
+  return 0;
+}
diff --git a/example/Persistence_representations/persistence_landscape.cpp b/example/Persistence_representations/persistence_landscape.cpp
new file mode 100644
index 00000000..400a9ae1
--- /dev/null
+++ b/example/Persistence_representations/persistence_landscape.cpp
@@ -0,0 +1,86 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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/>.
+ */
+
+#include <gudhi/Persistence_landscape.h>
+
+#include <iostream>
+#include <vector>
+#include <utility>
+
+using Persistence_landscape = Gudhi::Persistence_representations::Persistence_landscape;
+
+int main(int argc, char** argv) {
+  // create two simple vectors with birth--death pairs:
+
+  std::vector<std::pair<double, double> > persistence1;
+  std::vector<std::pair<double, double> > persistence2;
+
+  persistence1.push_back(std::make_pair(1, 2));
+  persistence1.push_back(std::make_pair(6, 8));
+  persistence1.push_back(std::make_pair(0, 4));
+  persistence1.push_back(std::make_pair(3, 8));
+
+  persistence2.push_back(std::make_pair(2, 9));
+  persistence2.push_back(std::make_pair(1, 6));
+  persistence2.push_back(std::make_pair(3, 5));
+  persistence2.push_back(std::make_pair(6, 10));
+
+  // create two persistence landscapes based on persistence1 and persistence2:
+  Persistence_landscape l1(persistence1);
+  Persistence_landscape l2(persistence2);
+
+  // This is how to compute integral of landscapes:
+  std::cout << "Integral of the first landscape : " << l1.compute_integral_of_landscape() << std::endl;
+  std::cout << "Integral of the second landscape : " << l2.compute_integral_of_landscape() << std::endl;
+
+  // And here how to write landscapes to stream:
+  std::cout << "l1 : " << l1 << std::endl;
+  std::cout << "l2 : " << l2 << std::endl;
+
+  // Arithmetic operations on landscapes:
+  Persistence_landscape sum = l1 + l2;
+  std::cout << "sum : " << sum << std::endl;
+
+  // here are the maxima of the functions:
+  std::cout << "Maximum of l1 : " << l1.compute_maximum() << std::endl;
+  std::cout << "Maximum of l2 : " << l2.compute_maximum() << std::endl;
+
+  // here are the norms of landscapes:
+  std::cout << "L^1 Norm of l1 : " << l1.compute_norm_of_landscape(1.) << std::endl;
+  std::cout << "L^1 Norm of l2 : " << l2.compute_norm_of_landscape(1.) << std::endl;
+
+  // here is the average of landscapes:
+  Persistence_landscape average;
+  average.compute_average({&l1, &l2});
+  std::cout << "average : " << average << std::endl;
+
+  // here is the distance of landscapes:
+  std::cout << "Distance : " << l1.distance(l2) << std::endl;
+
+  // here is the scalar product of landscapes:
+  std::cout << "Scalar product : " << l1.compute_scalar_product(l2) << std::endl;
+
+  // here is how to create a file which is suitable for visualization via gnuplot:
+  average.plot("average_landscape");
+
+  return 0;
+}
diff --git a/example/Persistence_representations/persistence_landscape_on_grid.cpp b/example/Persistence_representations/persistence_landscape_on_grid.cpp
new file mode 100644
index 00000000..b201b397
--- /dev/null
+++ b/example/Persistence_representations/persistence_landscape_on_grid.cpp
@@ -0,0 +1,82 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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/>.
+ */
+
+#include <gudhi/Persistence_landscape_on_grid.h>
+
+#include <iostream>
+#include <utility>
+#include <vector>
+
+using Persistence_landscape_on_grid = Gudhi::Persistence_representations::Persistence_landscape_on_grid;
+
+int main(int argc, char** argv) {
+  // create two simple vectors with birth--death pairs:
+
+  std::vector<std::pair<double, double> > persistence1;
+  std::vector<std::pair<double, double> > persistence2;
+
+  persistence1.push_back(std::make_pair(1, 2));
+  persistence1.push_back(std::make_pair(6, 8));
+  persistence1.push_back(std::make_pair(0, 4));
+  persistence1.push_back(std::make_pair(3, 8));
+
+  persistence2.push_back(std::make_pair(2, 9));
+  persistence2.push_back(std::make_pair(1, 6));
+  persistence2.push_back(std::make_pair(3, 5));
+  persistence2.push_back(std::make_pair(6, 10));
+
+  // create two persistence landscapes based on persistence1 and persistence2:
+  Persistence_landscape_on_grid l1(persistence1, 0, 11, 20);
+  Persistence_landscape_on_grid l2(persistence2, 0, 11, 20);
+
+  // This is how to compute integral of landscapes:
+  std::cout << "Integral of the first landscape : " << l1.compute_integral_of_landscape() << std::endl;
+  std::cout << "Integral of the second landscape : " << l2.compute_integral_of_landscape() << std::endl;
+
+  // And here how to write landscapes to stream:
+  std::cout << "l1 : " << l1 << std::endl;
+  std::cout << "l2 : " << l2 << std::endl;
+
+  // here are the maxima of the functions:
+  std::cout << "Maximum of l1 : " << l1.compute_maximum() << std::endl;
+  std::cout << "Maximum of l2 : " << l2.compute_maximum() << std::endl;
+
+  // here are the norms of landscapes:
+  std::cout << "L^1 Norm of l1 : " << l1.compute_norm_of_landscape(1.) << std::endl;
+  std::cout << "L^1 Norm of l2 : " << l2.compute_norm_of_landscape(1.) << std::endl;
+
+  // here is the average of landscapes:
+  Persistence_landscape_on_grid average;
+  average.compute_average({&l1, &l2});
+  std::cout << "average : " << average << std::endl;
+
+  // here is the distance of landscapes:
+  std::cout << "Distance : " << l1.distance(l2) << std::endl;
+
+  // here is the scalar product of landscapes:
+  std::cout << "Scalar product : " << l1.compute_scalar_product(l2) << std::endl;
+
+  // here is how to create a file which is suitable for visualization via gnuplot:
+  average.plot("average_landscape");
+
+  return 0;
+}
diff --git a/example/Persistence_representations/persistence_vectors.cpp b/example/Persistence_representations/persistence_vectors.cpp
new file mode 100644
index 00000000..834ae644
--- /dev/null
+++ b/example/Persistence_representations/persistence_vectors.cpp
@@ -0,0 +1,74 @@
+/*    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):       Pawel Dlotko
+ *
+ *    Copyright (C) 2016  INRIA (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/>.
+ */
+
+#include <gudhi/Persistence_vectors.h>
+
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <iomanip>
+#include <limits>
+#include <utility>
+
+using Vector_distances_in_diagram =
+    Gudhi::Persistence_representations::Vector_distances_in_diagram<Gudhi::Euclidean_distance>;
+
+int main(int argc, char** argv) {
+  // create two simple vectors with birth--death pairs:
+
+  std::vector<std::pair<double, double> > persistence1;
+  std::vector<std::pair<double, double> > persistence2;
+
+  persistence1.push_back(std::make_pair(1, 2));
+  persistence1.push_back(std::make_pair(6, 8));
+  persistence1.push_back(std::make_pair(0, 4));
+  persistence1.push_back(std::make_pair(3, 8));
+
+  persistence2.push_back(std::make_pair(2, 9));
+  persistence2.push_back(std::make_pair(1, 6));
+  persistence2.push_back(std::make_pair(3, 5));
+  persistence2.push_back(std::make_pair(6, 10));
+
+  // create two persistence vectors based on persistence1 and persistence2:
+  Vector_distances_in_diagram v1(persistence1, std::numeric_limits<size_t>::max());
+  Vector_distances_in_diagram v2(persistence2, std::numeric_limits<size_t>::max());
+
+  // writing to a stream:
+  std::cout << "v1 : " << v1 << std::endl;
+  std::cout << "v2 : " << v2 << std::endl;
+
+  // averages:
+  Vector_distances_in_diagram average;
+  average.compute_average({&v1, &v2});
+  std::cout << "Average : " << average << std::endl;
+
+  // computations of distances:
+  std::cout << "l^1 distance : " << v1.distance(v2) << std::endl;
+
+  // computations of scalar product:
+  std::cout << "Scalar product of l1 and l2 : " << v1.compute_scalar_product(v2) << std::endl;
+
+  // create a file with a gnuplot script:
+  v1.plot("plot_of_vector_representation");
+
+  return 0;
+}
diff --git a/example/Persistent_cohomology/CMakeLists.txt b/example/Persistent_cohomology/CMakeLists.txt
index f47de4c3..18e2913b 100644
--- a/example/Persistent_cohomology/CMakeLists.txt
+++ b/example/Persistent_cohomology/CMakeLists.txt
@@ -5,12 +5,6 @@ add_executable(plain_homology plain_homology.cpp)
 
 add_executable(persistence_from_simple_simplex_tree persistence_from_simple_simplex_tree.cpp)
 
-add_executable(rips_distance_matrix_persistence rips_distance_matrix_persistence.cpp)
-target_link_libraries(rips_distance_matrix_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
-add_executable(rips_persistence rips_persistence.cpp)
-target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
 add_executable(rips_persistence_step_by_step rips_persistence_step_by_step.cpp)
 target_link_libraries(rips_persistence_step_by_step ${Boost_PROGRAM_OPTIONS_LIBRARY})
 
@@ -23,8 +17,6 @@ target_link_libraries(persistence_from_file ${Boost_PROGRAM_OPTIONS_LIBRARY})
 if (TBB_FOUND)
   target_link_libraries(plain_homology ${TBB_LIBRARIES})
   target_link_libraries(persistence_from_simple_simplex_tree ${TBB_LIBRARIES})
-  target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES})
-  target_link_libraries(rips_persistence ${TBB_LIBRARIES})
   target_link_libraries(rips_persistence_step_by_step ${TBB_LIBRARIES})
   target_link_libraries(rips_persistence_via_boundary_matrix ${TBB_LIBRARIES})
   target_link_libraries(persistence_from_file ${TBB_LIBRARIES})
@@ -33,10 +25,6 @@ endif()
 add_test(NAME Persistent_cohomology_example_plain_homology COMMAND $<TARGET_FILE:plain_homology>)
 add_test(NAME Persistent_cohomology_example_from_simple_simplex_tree COMMAND $<TARGET_FILE:persistence_from_simple_simplex_tree>
     "1" "0")
-add_test(NAME Persistent_cohomology_example_from_rips_distance_matrix COMMAND $<TARGET_FILE:rips_distance_matrix_persistence>
-    "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0")
-add_test(NAME Persistent_cohomology_example_from_rips_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence>
-    "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3")
 add_test(NAME Persistent_cohomology_example_from_rips_step_by_step_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence_step_by_step>
     "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3")
 add_test(NAME Persistent_cohomology_example_via_boundary_matrix COMMAND $<TARGET_FILE:rips_persistence_via_boundary_matrix>
@@ -48,8 +36,6 @@ add_test(NAME Persistent_cohomology_example_from_file_3_3_100 COMMAND $<TARGET_F
 
 install(TARGETS plain_homology DESTINATION bin)
 install(TARGETS persistence_from_simple_simplex_tree DESTINATION bin)
-install(TARGETS rips_distance_matrix_persistence DESTINATION bin)
-install(TARGETS rips_persistence DESTINATION bin)
 install(TARGETS rips_persistence_step_by_step DESTINATION bin)
 install(TARGETS rips_persistence_via_boundary_matrix DESTINATION bin)
 install(TARGETS persistence_from_file DESTINATION bin)
@@ -69,53 +55,15 @@ if(GMP_FOUND)
 endif(GMP_FOUND)
 
 if(CGAL_FOUND)
-  add_executable(alpha_complex_3d_persistence alpha_complex_3d_persistence.cpp)
-  target_link_libraries(alpha_complex_3d_persistence ${CGAL_LIBRARY})
-  add_executable(exact_alpha_complex_3d_persistence exact_alpha_complex_3d_persistence.cpp)
-  target_link_libraries(exact_alpha_complex_3d_persistence ${CGAL_LIBRARY})
-  add_executable(weighted_alpha_complex_3d_persistence weighted_alpha_complex_3d_persistence.cpp)
-  target_link_libraries(weighted_alpha_complex_3d_persistence ${CGAL_LIBRARY})
-
-  if (TBB_FOUND)
-    target_link_libraries(alpha_complex_3d_persistence ${TBB_LIBRARIES})
-    target_link_libraries(exact_alpha_complex_3d_persistence ${TBB_LIBRARIES})
-    target_link_libraries(weighted_alpha_complex_3d_persistence ${TBB_LIBRARIES})
-  endif(TBB_FOUND)
-  add_test(NAME Persistent_cohomology_example_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
-      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45")
-  add_test(NAME Persistent_cohomology_example_exact_alpha_complex_3d COMMAND $<TARGET_FILE:exact_alpha_complex_3d_persistence>
-      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "2" "0.45")
-  add_test(NAME Persistent_cohomology_example_weighted_alpha_complex_3d COMMAND $<TARGET_FILE:weighted_alpha_complex_3d_persistence>
-      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.weights" "2" "0.45")
-
-  install(TARGETS alpha_complex_3d_persistence DESTINATION bin)
-  install(TARGETS exact_alpha_complex_3d_persistence DESTINATION bin)
-  install(TARGETS weighted_alpha_complex_3d_persistence DESTINATION bin)
-
   if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
-    add_executable (alpha_complex_persistence alpha_complex_persistence.cpp)
-    target_link_libraries(alpha_complex_persistence
-        ${CGAL_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
-    add_executable(periodic_alpha_complex_3d_persistence periodic_alpha_complex_3d_persistence.cpp)
-    target_link_libraries(periodic_alpha_complex_3d_persistence ${CGAL_LIBRARY})
-
     add_executable(custom_persistence_sort custom_persistence_sort.cpp)
     target_link_libraries(custom_persistence_sort ${CGAL_LIBRARY})
 
     if (TBB_FOUND)
-      target_link_libraries(alpha_complex_persistence ${TBB_LIBRARIES})
-      target_link_libraries(periodic_alpha_complex_3d_persistence ${TBB_LIBRARIES})
       target_link_libraries(custom_persistence_sort ${TBB_LIBRARIES})
     endif(TBB_FOUND)
-    add_test(NAME Persistent_cohomology_example_alpha_complex COMMAND $<TARGET_FILE:alpha_complex_persistence>
-        "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-p" "2" "-m" "0.45")
-    add_test(NAME Persistent_cohomology_example_periodic_alpha_complex_3d COMMAND $<TARGET_FILE:periodic_alpha_complex_3d_persistence>
-        "${CMAKE_SOURCE_DIR}/data/points/grid_10_10_10_in_0_1.off" "${CMAKE_SOURCE_DIR}/data/points/iso_cuboid_3_in_0_1.txt" "2" "0")
     add_test(NAME Persistent_cohomology_example_custom_persistence_sort COMMAND $<TARGET_FILE:custom_persistence_sort>)
 
-    install(TARGETS alpha_complex_persistence DESTINATION bin)
-    install(TARGETS periodic_alpha_complex_3d_persistence DESTINATION bin)
     install(TARGETS custom_persistence_sort DESTINATION bin)
 
   endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.7.0)
diff --git a/example/Persistent_cohomology/README b/example/Persistent_cohomology/README
index 794b94ae..f39d9584 100644
--- a/example/Persistent_cohomology/README
+++ b/example/Persistent_cohomology/README
@@ -1,43 +1,14 @@
-To build the example, run in a Terminal:
+To build the examples, run in a Terminal:
 
-cd /path-to-example/
+cd /path-to-examples/
 cmake .
 make
 
 ***********************************************************************************************************************
 Example of use of RIPS:
 
-Computation of the persistent homology with Z/2Z coefficients of the Rips complex on points 
-sampling a Klein bottle:
-
-./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 2
-
-output:
-2  0 0 inf 
-2  1 0.0983494 inf 
-2  1 0.104347 inf 
-2  2 0.138335 inf 
-
-
-Every line is of this format: p1*...*pr   dim b d
-where
- 	p1*...*pr is the product of prime numbers pi such that the homology feature exists in homology with Z/piZ coefficients. 
-	dim is the dimension of the homological feature,
-    b and d are respectively the birth and death of the feature and
-   
-
-
-with Z/3Z coefficients:
-
-./rips_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.5 -d 3 -p 3
-
-output:
-3  0 0 inf 
-3  1 0.0983494 inf 
-3  1 0.104347 inf 
-3  2 0.138335 inf 
-
-and the computation with Z/2Z and Z/3Z coefficients simultaneously:
+Computation of the persistent homology with Z/2Z and Z/3Z coefficients simultaneously of the Rips complex 
+on points sampling a 3D torus:
 
 ./rips_multifield_persistence ../../data/points/tore3D_1307.off -r 0.25 -m 0.12 -d 3 -p 2 -q 3
 
@@ -53,7 +24,13 @@ output:
 6  0 0 0.12047 
 6  0 0 0.120414 
 
-and finally the computation with all Z/pZ for 2 <= p <= 71 (20 first prime numbers):
+Every line is of this format: p1*...*pr   dim b d
+where
+ 	p1*...*pr is the product of prime numbers pi such that the homology feature exists in homology with Z/piZ coefficients. 
+	dim is the dimension of the homological feature,
+    b and d are respectively the birth and death of the feature and
+
+and the computation with all Z/pZ for 2 <= p <= 71 (20 first prime numbers):
 
  ./rips_multifield_persistence ../../data/points/Kl.off -r 0.25 -m 0.5 -d 3 -p 2 -q 71
 
@@ -70,82 +47,6 @@ output:
 557940830126698960967415390  0 0 0.120414 
 
 ***********************************************************************************************************************
-Example of use of ALPHA:
-
-For a more verbose mode, please run cmake with option "DEBUG_TRACES=TRUE" and recompile the programs.
-
-1) 3D special case
-------------------
-Computation of the persistent homology with Z/2Z coefficients of the alpha complex on points 
-sampling a torus 3D:
-
-./alpha_complex_3d_persistence ../../data/points/tore3D_300.off 2 0.45
-
-output:
-Simplex_tree dim: 3
-2  0 0 inf 
-2  1 0.0682162 1.0001 
-2  1 0.0934117 1.00003 
-2  2 0.56444 1.03938 
-
-Here we retrieve expected Betti numbers on a tore 3D:
-Betti numbers[0] = 1
-Betti numbers[1] = 2
-Betti numbers[2] = 1
-
-N.B.: - alpha_complex_3d_persistence accepts only OFF files in 3D dimension.
-      - filtration values are alpha square values
-
-2) d-Dimension case
--------------------
-Computation of the persistent homology with Z/2Z coefficients of the alpha complex on points 
-sampling a torus 3D:
-
-./alpha_complex_persistence -r 32 -p 2 -m 0.45 ../../data/points/tore3D_300.off
-
-output:
-Alpha complex is of dimension 3 - 9273 simplices - 300 vertices.
-Simplex_tree dim: 3
-2  0 0 inf 
-2  1 0.0682162 1.0001 
-2  1 0.0934117 1.00003 
-2  2 0.56444 1.03938 
-
-Here we retrieve expected Betti numbers on a tore 3D:
-Betti numbers[0] = 1
-Betti numbers[1] = 2
-Betti numbers[2] = 1
-
-N.B.: - alpha_complex_persistence accepts OFF files in d-Dimension.
-      - filtration values are alpha square values
-
-3) 3D periodic special case
----------------------------
-./periodic_alpha_complex_3d_persistence ../../data/points/grid_10_10_10_in_0_1.off ../../data/points/iso_cuboid_3_in_0_1.txt 3 1.0
-
-output:
-Periodic Delaunay computed.
-Simplex_tree dim: 3
-3  0 0 inf 
-3  1 0.0025 inf 
-3  1 0.0025 inf 
-3  1 0.0025 inf 
-3  2 0.005 inf 
-3  2 0.005 inf 
-3  2 0.005 inf 
-3  3 0.0075 inf 
-
-Here we retrieve expected Betti numbers on a tore 3D:
-Betti numbers[0] = 1
-Betti numbers[1] = 3
-Betti numbers[2] = 3
-Betti numbers[3] = 1
-
-N.B.: - periodic_alpha_complex_3d_persistence accepts only OFF files in 3D dimension. In this example, the periodic cube
-is hard coded to { x = [0,1]; y = [0,1]; z = [0,1] }
-      - filtration values are alpha square values
-
-***********************************************************************************************************************
 Example of use of PLAIN HOMOLOGY:
 
 This example computes the plain homology of the following simplicial complex without filtration values:
diff --git a/example/Persistent_cohomology/alpha_complex_3d_helper.h b/example/Persistent_cohomology/alpha_complex_3d_helper.h
deleted file mode 100644
index 7865e4ec..00000000
--- a/example/Persistent_cohomology/alpha_complex_3d_helper.h
+++ /dev/null
@@ -1,76 +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):       Vincent Rouvreau
- *
- *    Copyright (C) 2014  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/>.
- */
-
-#ifndef ALPHA_COMPLEX_3D_HELPER_H_
-#define ALPHA_COMPLEX_3D_HELPER_H_
-
-template<class Vertex_list, class Cell_handle>
-Vertex_list from_cell(const Cell_handle& ch) {
-  Vertex_list the_list;
-  for (auto i = 0; i < 4; i++) {
-#ifdef DEBUG_TRACES
-    std::cout << "from cell[" << i << "]=" << ch->vertex(i)->point() << std::endl;
-#endif  // DEBUG_TRACES
-    the_list.push_back(ch->vertex(i));
-  }
-  return the_list;
-}
-
-template<class Vertex_list, class Facet>
-Vertex_list from_facet(const Facet& fct) {
-  Vertex_list the_list;
-  for (auto i = 0; i < 4; i++) {
-    if (fct.second != i) {
-#ifdef DEBUG_TRACES
-      std::cout << "from facet=[" << i << "]" << fct.first->vertex(i)->point() << std::endl;
-#endif  // DEBUG_TRACES
-      the_list.push_back(fct.first->vertex(i));
-    }
-  }
-  return the_list;
-}
-
-template<class Vertex_list, class Edge_3>
-Vertex_list from_edge(const Edge_3& edg) {
-  Vertex_list the_list;
-  for (auto i = 0; i < 4; i++) {
-    if ((edg.second == i) || (edg.third == i)) {
-#ifdef DEBUG_TRACES
-      std::cout << "from edge[" << i << "]=" << edg.first->vertex(i)->point() << std::endl;
-#endif  // DEBUG_TRACES
-      the_list.push_back(edg.first->vertex(i));
-    }
-  }
-  return the_list;
-}
-
-template<class Vertex_list, class Vertex_handle>
-Vertex_list from_vertex(const Vertex_handle& vh) {
-  Vertex_list the_list;
-#ifdef DEBUG_TRACES
-  std::cout << "from vertex=" << vh->point() << std::endl;
-#endif  // DEBUG_TRACES
-  the_list.push_back(vh);
-  return the_list;
-}
-
-#endif  // ALPHA_COMPLEX_3D_HELPER_H_
diff --git a/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp
deleted file mode 100644
index f63ff0f6..00000000
--- a/example/Persistent_cohomology/alpha_complex_3d_persistence.cpp
+++ /dev/null
@@ -1,242 +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):       Vincent Rouvreau
- *
- *    Copyright (C) 2014  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 <boost/variant.hpp>
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_3D_off_io.h>
-
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Delaunay_triangulation_3.h>
-#include <CGAL/Alpha_shape_3.h>
-#include <CGAL/iterator.h>
-
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <tuple>
-#include <map>
-#include <utility>
-#include <list>
-#include <vector>
-
-#include "alpha_complex_3d_helper.h"
-
-// Alpha_shape_3 templates type definitions
-using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel>;
-using Fb = CGAL::Alpha_shape_cell_base_3<Kernel>;
-using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
-using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
-using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>;
-
-// From file type definition
-using Point_3 = Kernel::Point_3;
-
-// filtration with alpha values needed type definition
-using Alpha_value_type = Alpha_shape_3::FT;
-using Object = CGAL::Object;
-using Dispatch = CGAL::Dispatch_output_iterator<
-          CGAL::cpp11::tuple<Object, Alpha_value_type>,
-          CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >,
-          std::back_insert_iterator< std::vector<Alpha_value_type> > > >;
-using Cell_handle = Alpha_shape_3::Cell_handle;
-using Facet = Alpha_shape_3::Facet;
-using Edge_3 = Alpha_shape_3::Edge;
-using Vertex_handle = Alpha_shape_3::Vertex_handle;
-using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>;
-
-// gudhi type definition
-using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = ST::Filtration_value;
-using Simplex_tree_vertex = ST::Vertex_handle;
-using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >;
-using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >;
-using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >;
-
-void usage(const std::string& progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " <<
-               "> 0] min_persistence[float >= -1.0]\n";
-  std::cerr << "  path_to_OFF_file is the path to your points cloud in OFF format.\n";
-  exit(-1);
-}
-
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 4) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[0]);
-  }
-
-  int coeff_field_characteristic = atoi(argv[2]);
-
-  Filtration_value min_persistence = 0.0;
-  int returnedScanValue = sscanf(argv[3], "%f", &min_persistence);
-  if ((returnedScanValue == EOF) || (min_persistence < -1.0)) {
-    std::cerr << "Error: " << argv[3] << " is not correct\n";
-    usage(argv[0]);
-  }
-
-  // Read points from file
-  std::string offInputFile(argv[1]);
-  // Read the OFF file (input file name given as parameter) and triangulate points
-  Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile);
-  // Check the read operation was correct
-  if (!off_reader.is_valid()) {
-    std::cerr << "Unable to read file " << offInputFile << std::endl;
-    usage(argv[0]);
-  }
-
-  // Retrieve the triangulation
-  std::vector<Point_3> lp = off_reader.get_point_cloud();
-
-  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
-  Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL);
-#ifdef DEBUG_TRACES
-  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
-#endif  // DEBUG_TRACES
-
-  // filtration with alpha values from alpha shape
-  std::vector<Object> the_objects;
-  std::vector<Alpha_value_type> the_alpha_values;
-
-  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                  std::back_inserter(the_alpha_values));
-
-  as.filtration_with_alpha_values(disp);
-#ifdef DEBUG_TRACES
-  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
-#endif  // DEBUG_TRACES
-
-  Alpha_shape_3::size_type count_vertices = 0;
-  Alpha_shape_3::size_type count_edges = 0;
-  Alpha_shape_3::size_type count_facets = 0;
-  Alpha_shape_3::size_type count_cells = 0;
-
-  // Loop on objects vector
-  Vertex_list vertex_list;
-  ST simplex_tree;
-  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
-  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-  int dim_max = 0;
-  Filtration_value filtration_max = 0.0;
-  for (auto object_iterator : the_objects) {
-    // Retrieve Alpha shape vertex list from object
-    if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
-      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
-      count_cells++;
-      if (dim_max < 3) {
-        // Cell is of dim 3
-        dim_max = 3;
-      }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
-      vertex_list = from_facet<Vertex_list, Facet>(*facet);
-      count_facets++;
-      if (dim_max < 2) {
-        // Facet is of dim 2
-        dim_max = 2;
-      }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
-      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
-      count_edges++;
-      if (dim_max < 1) {
-        // Edge_3 is of dim 1
-        dim_max = 1;
-      }
-    } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
-      count_vertices++;
-      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
-    }
-    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
-    Simplex_tree_vector_vertex the_simplex_tree;
-    for (auto the_alpha_shape_vertex : vertex_list) {
-      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
-      if (the_map_iterator == map_cgal_simplex_tree.end()) {
-        // alpha shape not found
-        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-        map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex));
-      } else {
-        // alpha shape found
-        Simplex_tree_vertex vertex = the_map_iterator->second;
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-      }
-    }
-    // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator);
-#ifdef DEBUG_TRACES
-    std::cout << "filtration = " << filtr << std::endl;
-#endif  // DEBUG_TRACES
-    if (filtr > filtration_max) {
-      filtration_max = filtr;
-    }
-    simplex_tree.insert_simplex(the_simplex_tree, filtr);
-    if (the_alpha_value_iterator != the_alpha_values.end())
-      ++the_alpha_value_iterator;
-    else
-      std::cout << "This shall not happen" << std::endl;
-  }
-  simplex_tree.set_dimension(dim_max);
-
-#ifdef DEBUG_TRACES
-  std::cout << "vertices \t\t" << count_vertices << std::endl;
-  std::cout << "edges \t\t" << count_edges << std::endl;
-  std::cout << "facets \t\t" << count_facets << std::endl;
-  std::cout << "cells \t\t" << count_cells << std::endl;
-
-
-  std::cout << "Information of the Simplex Tree: " << std::endl;
-  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
-  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
-  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
-#endif  // DEBUG_TRACES
-
-#ifdef DEBUG_TRACES
-  std::cout << "Iterator on vertices: " << std::endl;
-  for (auto vertex : simplex_tree.complex_vertex_range()) {
-    std::cout << vertex << " ";
-  }
-#endif  // DEBUG_TRACES
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
-  // Compute the persistence diagram of the complex
-  PCOH pcoh(simplex_tree);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(coeff_field_characteristic);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  pcoh.output_diagram();
-
-  return 0;
-}
diff --git a/example/Persistent_cohomology/alpha_complex_persistence.cpp b/example/Persistent_cohomology/alpha_complex_persistence.cpp
deleted file mode 100644
index 9e84e91f..00000000
--- a/example/Persistent_cohomology/alpha_complex_persistence.cpp
+++ /dev/null
@@ -1,125 +0,0 @@
-#include <boost/program_options.hpp>
-
-#include <CGAL/Epick_d.h>
-
-#include <gudhi/Alpha_complex.h>
-#include <gudhi/Persistent_cohomology.h>
-// to construct a simplex_tree from alpha complex
-#include <gudhi/Simplex_tree.h>
-
-#include <iostream>
-#include <string>
-#include <limits>  // for numeric_limits
-
-using Simplex_tree = Gudhi::Simplex_tree<>;
-using Filtration_value = Simplex_tree::Filtration_value;
-
-void program_options(int argc, char * argv[]
-                     , std::string & off_file_points
-                     , std::string & output_file_diag
-                     , Filtration_value & alpha_square_max_value
-                     , int & coeff_field_characteristic
-                     , Filtration_value & min_persistence);
-
-int main(int argc, char **argv) {
-  std::string off_file_points;
-  std::string output_file_diag;
-  Filtration_value alpha_square_max_value;
-  int coeff_field_characteristic;
-  Filtration_value min_persistence;
-
-  program_options(argc, argv, off_file_points, output_file_diag, alpha_square_max_value,
-                  coeff_field_characteristic, min_persistence);
-
-  // ----------------------------------------------------------------------------
-  // Init of an alpha complex from an OFF file
-  // ----------------------------------------------------------------------------
-  using Kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >;
-  Gudhi::alpha_complex::Alpha_complex<Kernel> alpha_complex_from_file(off_file_points);
-
-  Simplex_tree simplex;
-  if (alpha_complex_from_file.create_complex(simplex, alpha_square_max_value)) {
-    // ----------------------------------------------------------------------------
-    // Display information about the alpha complex
-    // ----------------------------------------------------------------------------
-    std::cout << "Simplicial complex is of dimension " << simplex.dimension() <<
-        " - " << simplex.num_simplices() << " simplices - " <<
-        simplex.num_vertices() << " vertices." << std::endl;
-
-    // Sort the simplices in the order of the filtration
-    simplex.initialize_filtration();
-
-    std::cout << "Simplex_tree dim: " << simplex.dimension() << std::endl;
-    // Compute the persistence diagram of the complex
-    Gudhi::persistent_cohomology::Persistent_cohomology< Simplex_tree,
-        Gudhi::persistent_cohomology::Field_Zp > pcoh(simplex);
-    // initializes the coefficient field for homology
-    pcoh.init_coefficients(coeff_field_characteristic);
-
-    pcoh.compute_persistent_cohomology(min_persistence);
-
-    // Output the diagram in filediag
-    if (output_file_diag.empty()) {
-      pcoh.output_diagram();
-    } else {
-      std::cout << "Result in file: " << output_file_diag << std::endl;
-      std::ofstream out(output_file_diag);
-      pcoh.output_diagram(out);
-      out.close();
-    }
-  }
-
-  return 0;
-}
-
-void program_options(int argc, char * argv[]
-                     , std::string & off_file_points
-                     , std::string & output_file_diag
-                     , Filtration_value & alpha_square_max_value
-                     , int & coeff_field_characteristic
-                     , Filtration_value & min_persistence) {
-  namespace po = boost::program_options;
-  po::options_description hidden("Hidden options");
-  hidden.add_options()
-      ("input-file", po::value<std::string>(&off_file_points),
-       "Name of file containing a point set. Format is one point per line:   X1 ... Xd ");
-
-  po::options_description visible("Allowed options", 100);
-  visible.add_options()
-      ("help,h", "produce help message")
-      ("output-file,o", po::value<std::string>(&output_file_diag)->default_value(std::string()),
-       "Name of file in which the persistence diagram is written. Default print in std::cout")
-      ("max-alpha-square-value,r",
-      po::value<Filtration_value>(&alpha_square_max_value)->default_value(std::numeric_limits<Filtration_value>::infinity()),
-       "Maximal alpha square value for the Alpha complex construction.")
-      ("field-charac,p", po::value<int>(&coeff_field_characteristic)->default_value(11),
-       "Characteristic p of the coefficient field Z/pZ for computing homology.")
-      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
-
-  po::positional_options_description pos;
-  pos.add("input-file", 1);
-
-  po::options_description all;
-  all.add(visible).add(hidden);
-
-  po::variables_map vm;
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), vm);
-  po::notify(vm);
-
-  if (vm.count("help") || !vm.count("input-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of an Alpha complex defined on a set of input points.\n \n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-    std::cout << visible << std::endl;
-    std::abort();
-  }
-}
diff --git a/example/Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp
deleted file mode 100644
index 09561d03..00000000
--- a/example/Persistent_cohomology/exact_alpha_complex_3d_persistence.cpp
+++ /dev/null
@@ -1,244 +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):       Vincent Rouvreau
- *
- *    Copyright (C) 2014  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/>.
- */
-
-#include <boost/variant.hpp>
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_3D_off_io.h>
-
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Delaunay_triangulation_3.h>
-#include <CGAL/Alpha_shape_3.h>
-#include <CGAL/iterator.h>
-
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <tuple>
-#include <map>
-#include <utility>
-#include <list>
-#include <vector>
-
-#include "alpha_complex_3d_helper.h"
-
-// Alpha_shape_3 templates type definitions
-using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-using Exact_tag = CGAL::Tag_true;
-using Vb = CGAL::Alpha_shape_vertex_base_3<Kernel, CGAL::Default, Exact_tag>;
-using Fb = CGAL::Alpha_shape_cell_base_3<Kernel, CGAL::Default, Exact_tag>;
-using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
-using Triangulation_3 = CGAL::Delaunay_triangulation_3<Kernel, Tds>;
-using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3, Exact_tag>;
-
-// From file type definition
-using Point_3 = Kernel::Point_3;
-
-// filtration with alpha values needed type definition
-using Alpha_value_type = Alpha_shape_3::FT;
-using Object = CGAL::Object;
-using Dispatch = CGAL::Dispatch_output_iterator<
-          CGAL::cpp11::tuple<Object, Alpha_value_type>,
-          CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >,
-          std::back_insert_iterator< std::vector<Alpha_value_type> > > >;
-using Cell_handle = Alpha_shape_3::Cell_handle;
-using Facet = Alpha_shape_3::Facet;
-using Edge_3 = Alpha_shape_3::Edge;
-using Vertex_handle = Alpha_shape_3::Vertex_handle;
-using Vertex_list = std::list<Vertex_handle>;
-
-// gudhi type definition
-using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = ST::Filtration_value;
-using Simplex_tree_vertex = ST::Vertex_handle;
-using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >;
-using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >;
-using PCOH = Gudhi::persistent_cohomology::Persistent_cohomology< ST, Gudhi::persistent_cohomology::Field_Zp >;
-
-void usage(char * const progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file coeff_field_characteristic[integer " <<
-               "> 0] min_persistence[float >= -1.0]\n";
-  std::cerr << "  path_to_OFF_file is the path to your points cloud in OFF format.\n";
-  exit(-1);
-}
-
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 4) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[0]);
-  }
-
-  int coeff_field_characteristic = atoi(argv[2]);
-
-  Filtration_value min_persistence = 0.0;
-  int returnedScanValue = sscanf(argv[3], "%f", &min_persistence);
-  if ((returnedScanValue == EOF) || (min_persistence < -1.0)) {
-    std::cerr << "Error: " << argv[3] << " is not correct\n";
-    usage(argv[0]);
-  }
-
-  // Read points from file
-  std::string offInputFile(argv[1]);
-  // Read the OFF file (input file name given as parameter) and triangulate points
-  Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile);
-  // Check the read operation was correct
-  if (!off_reader.is_valid()) {
-    std::cerr << "Unable to read file " << offInputFile << std::endl;
-    usage(argv[0]);
-  }
-
-  // Retrieve the triangulation
-  std::vector<Point_3> lp = off_reader.get_point_cloud();
-
-  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
-  Alpha_shape_3 as(lp.begin(), lp.end(), 0, Alpha_shape_3::GENERAL);
-#ifdef DEBUG_TRACES
-  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
-#endif  // DEBUG_TRACES
-
-  // filtration with alpha values from alpha shape
-  std::vector<Object> the_objects;
-  std::vector<Alpha_value_type> the_alpha_values;
-
-  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                  std::back_inserter(the_alpha_values));
-
-  as.filtration_with_alpha_values(disp);
-#ifdef DEBUG_TRACES
-  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
-#endif  // DEBUG_TRACES
-
-  Alpha_shape_3::size_type count_vertices = 0;
-  Alpha_shape_3::size_type count_edges = 0;
-  Alpha_shape_3::size_type count_facets = 0;
-  Alpha_shape_3::size_type count_cells = 0;
-
-  // Loop on objects vector
-  Vertex_list vertex_list;
-  ST simplex_tree;
-  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
-  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-  int dim_max = 0;
-  Filtration_value filtration_max = 0.0;
-  for (auto object_iterator : the_objects) {
-    // Retrieve Alpha shape vertex list from object
-    if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
-      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
-      count_cells++;
-      if (dim_max < 3) {
-        // Cell is of dim 3
-        dim_max = 3;
-      }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
-      vertex_list = from_facet<Vertex_list, Facet>(*facet);
-      count_facets++;
-      if (dim_max < 2) {
-        // Facet is of dim 2
-        dim_max = 2;
-      }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
-      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
-      count_edges++;
-      if (dim_max < 1) {
-        // Edge_3 is of dim 1
-        dim_max = 1;
-      }
-    } else if (const Vertex_handle * vertex = CGAL::object_cast<Vertex_handle>(&object_iterator)) {
-      count_vertices++;
-      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
-    }
-    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
-    Simplex_tree_vector_vertex the_simplex_tree;
-    for (auto the_alpha_shape_vertex : vertex_list) {
-      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
-      if (the_map_iterator == map_cgal_simplex_tree.end()) {
-        // alpha shape not found
-        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-        map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex));
-      } else {
-        // alpha shape found
-        Simplex_tree_vertex vertex = the_map_iterator->second;
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-      }
-    }
-    // Construction of the simplex_tree
-    // you can also use the_alpha_value_iterator->exact()
-    Filtration_value filtr = /*std::sqrt*/CGAL::to_double(the_alpha_value_iterator->exact());
-#ifdef DEBUG_TRACES
-    std::cout << "filtration = " << filtr << std::endl;
-#endif  // DEBUG_TRACES
-    if (filtr > filtration_max) {
-      filtration_max = filtr;
-    }
-    simplex_tree.insert_simplex(the_simplex_tree, filtr);
-    if (the_alpha_value_iterator != the_alpha_values.end())
-      ++the_alpha_value_iterator;
-    else
-      std::cout << "This shall not happen" << std::endl;
-  }
-  simplex_tree.set_dimension(dim_max);
-
-#ifdef DEBUG_TRACES
-  std::cout << "vertices \t\t" << count_vertices << std::endl;
-  std::cout << "edges \t\t" << count_edges << std::endl;
-  std::cout << "facets \t\t" << count_facets << std::endl;
-  std::cout << "cells \t\t" << count_cells << std::endl;
-
-
-  std::cout << "Information of the Simplex Tree: " << std::endl;
-  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
-  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
-  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
-#endif  // DEBUG_TRACES
-
-#ifdef DEBUG_TRACES
-  std::cout << "Iterator on vertices: " << std::endl;
-  for (auto vertex : simplex_tree.complex_vertex_range()) {
-    std::cout << vertex << " ";
-  }
-#endif  // DEBUG_TRACES
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
-  // Compute the persistence diagram of the complex
-  PCOH pcoh(simplex_tree);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(coeff_field_characteristic);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  pcoh.output_diagram();
-
-  return 0;
-}
diff --git a/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp
deleted file mode 100644
index 8140a3c5..00000000
--- a/example/Persistent_cohomology/periodic_alpha_complex_3d_persistence.cpp
+++ /dev/null
@@ -1,268 +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):       Vincent Rouvreau
- *
- *    Copyright (C) 2014  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 <boost/variant.hpp>
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_3D_off_io.h>
-
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Periodic_3_Delaunay_triangulation_traits_3.h>
-#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
-#include <CGAL/Alpha_shape_3.h>
-#include <CGAL/iterator.h>
-
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <tuple>
-#include <map>
-#include <utility>
-#include <list>
-#include <vector>
-#include <cstdlib>
-
-#include "alpha_complex_3d_helper.h"
-
-// Traits
-using K = CGAL::Exact_predicates_inexact_constructions_kernel;
-using PK = CGAL::Periodic_3_Delaunay_triangulation_traits_3<K>;
-// Vertex type
-using DsVb = CGAL::Periodic_3_triangulation_ds_vertex_base_3<>;
-using Vb = CGAL::Triangulation_vertex_base_3<PK, DsVb>;
-using AsVb = CGAL::Alpha_shape_vertex_base_3<PK, Vb>;
-// Cell type
-using DsCb = CGAL::Periodic_3_triangulation_ds_cell_base_3<>;
-using Cb = CGAL::Triangulation_cell_base_3<PK, DsCb>;
-using AsCb = CGAL::Alpha_shape_cell_base_3<PK, Cb>;
-using Tds = CGAL::Triangulation_data_structure_3<AsVb, AsCb>;
-using P3DT3 = CGAL::Periodic_3_Delaunay_triangulation_3<PK, Tds>;
-using Alpha_shape_3 = CGAL::Alpha_shape_3<P3DT3>;
-using Point_3 = PK::Point_3;
-
-// filtration with alpha values needed type definition
-using Alpha_value_type = Alpha_shape_3::FT;
-using Object = CGAL::Object;
-using Dispatch = CGAL::Dispatch_output_iterator<
-    CGAL::cpp11::tuple<Object, Alpha_value_type>,
-    CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >,
-    std::back_insert_iterator< std::vector<Alpha_value_type> > > >;
-using Cell_handle = Alpha_shape_3::Cell_handle;
-using Facet = Alpha_shape_3::Facet;
-using Edge_3 = Alpha_shape_3::Edge;
-using Vertex_handle = Alpha_shape_3::Vertex_handle;
-using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>;
-
-// gudhi type definition
-using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = ST::Filtration_value;
-using Simplex_tree_vertex = ST::Vertex_handle;
-using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >;
-using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >;
-using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<
-    ST, Gudhi::persistent_cohomology::Field_Zp >;
-
-void usage(char * const progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file path_to_iso_cuboid_3_file coeff_field_characteristic[" <<
-               "integer > 0] min_persistence[float >= -1.0]\n" <<
-               "  path_to_OFF_file is the path to your points cloud in OFF format.\n" <<
-               "  path_to_iso_cuboid_3_file is the path to the iso cuboid file with the following format :\n" <<
-               "    x_min y_min z_min x_max y_max z_max\n" <<
-               "  In this example, the periodic cube will be " <<
-               "{ x = [x_min,x_max]; y = [y_min,y_max]; z = [z_min,z_max] }.\n" <<
-               "  For more information, please refer to\n" <<
-               "    https://doc.cgal.org/latest/Kernel_23/classCGAL_1_1Iso__cuboid__3.html\n";
-
-  exit(-1);
-}
-
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 5) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[0]);
-  }
-
-  int coeff_field_characteristic = atoi(argv[3]);
-  Filtration_value min_persistence = strtof(argv[4], nullptr);
-
-  // Read points from file
-  std::string offInputFile(argv[1]);
-  // Read the OFF file (input file name given as parameter) and triangulate points
-  Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile);
-  // Check the read operation was correct
-  if (!off_reader.is_valid()) {
-    std::cerr << "Unable to read file " << offInputFile << std::endl;
-    usage(argv[0]);
-  }
-
-  // Read iso_cuboid_3 information from file
-  std::ifstream iso_cuboid_str(argv[2]);
-  double x_min, y_min, z_min, x_max, y_max, z_max;
-  if (iso_cuboid_str.good()) {
-    iso_cuboid_str >> x_min >> y_min >> z_min >> x_max >> y_max >> z_max;
-  } else {
-    std::cerr << "Unable to read file " << argv[2] << std::endl;
-    usage(argv[0]);
-  }
-
-  // Retrieve the triangulation
-  std::vector<Point_3> lp = off_reader.get_point_cloud();
-
-  // Define the periodic cube
-  P3DT3 pdt(PK::Iso_cuboid_3(x_min, y_min, z_min, x_max, y_max, z_max));
-  // Heuristic for inserting large point sets (if pts is reasonably large)
-  pdt.insert(lp.begin(), lp.end(), true);
-  // As pdt won't be modified anymore switch to 1-sheeted cover if possible
-  if (pdt.is_triangulation_in_1_sheet()) pdt.convert_to_1_sheeted_covering();
-  std::cout << "Periodic Delaunay computed." << std::endl;
-
-  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode. This is the default mode
-  // Maybe need to set it to GENERAL mode
-  Alpha_shape_3 as(pdt, 0, Alpha_shape_3::GENERAL);
-
-  // filtration with alpha values from alpha shape
-  std::vector<Object> the_objects;
-  std::vector<Alpha_value_type> the_alpha_values;
-
-  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                  std::back_inserter(the_alpha_values));
-
-  as.filtration_with_alpha_values(disp);
-#ifdef DEBUG_TRACES
-  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
-#endif  // DEBUG_TRACES
-
-  Alpha_shape_3::size_type count_vertices = 0;
-  Alpha_shape_3::size_type count_edges = 0;
-  Alpha_shape_3::size_type count_facets = 0;
-  Alpha_shape_3::size_type count_cells = 0;
-
-  // Loop on objects vector
-  Vertex_list vertex_list;
-  ST simplex_tree;
-  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
-  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-  int dim_max = 0;
-  Filtration_value filtration_max = 0.0;
-  for (auto object_iterator : the_objects) {
-    // Retrieve Alpha shape vertex list from object
-    if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
-      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
-      count_cells++;
-      if (dim_max < 3) {
-        // Cell is of dim 3
-        dim_max = 3;
-      }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
-      vertex_list = from_facet<Vertex_list, Facet>(*facet);
-      count_facets++;
-      if (dim_max < 2) {
-        // Facet is of dim 2
-        dim_max = 2;
-      }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
-      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
-      count_edges++;
-      if (dim_max < 1) {
-        // Edge_3 is of dim 1
-        dim_max = 1;
-      }
-    } else if (const Alpha_shape_3::Vertex_handle * vertex =
-               CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
-      count_vertices++;
-      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
-    }
-    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
-    Simplex_tree_vector_vertex the_simplex_tree;
-    for (auto the_alpha_shape_vertex : vertex_list) {
-      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
-      if (the_map_iterator == map_cgal_simplex_tree.end()) {
-        // alpha shape not found
-        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-        map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex));
-      } else {
-        // alpha shape found
-        Simplex_tree_vertex vertex = the_map_iterator->second;
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-      }
-    }
-    // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator);
-#ifdef DEBUG_TRACES
-    std::cout << "filtration = " << filtr << std::endl;
-#endif  // DEBUG_TRACES
-    if (filtr > filtration_max) {
-      filtration_max = filtr;
-    }
-    simplex_tree.insert_simplex(the_simplex_tree, filtr);
-    if (the_alpha_value_iterator != the_alpha_values.end())
-      ++the_alpha_value_iterator;
-    else
-      std::cout << "This shall not happen" << std::endl;
-  }
-  simplex_tree.set_dimension(dim_max);
-
-#ifdef DEBUG_TRACES
-  std::cout << "vertices \t\t" << count_vertices << std::endl;
-  std::cout << "edges \t\t" << count_edges << std::endl;
-  std::cout << "facets \t\t" << count_facets << std::endl;
-  std::cout << "cells \t\t" << count_cells << std::endl;
-
-
-  std::cout << "Information of the Simplex Tree: " << std::endl;
-  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
-  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
-  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
-#endif  // DEBUG_TRACES
-
-#ifdef DEBUG_TRACES
-  std::cout << "Iterator on vertices: " << std::endl;
-  for (auto vertex : simplex_tree.complex_vertex_range()) {
-    std::cout << vertex << " ";
-  }
-#endif  // DEBUG_TRACES
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(simplex_tree, true);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(coeff_field_characteristic);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  pcoh.output_diagram();
-
-  return 0;
-}
diff --git a/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp b/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp
index 8214d66a..8ef479d4 100644
--- a/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp
+++ b/example/Persistent_cohomology/persistence_from_simple_simplex_tree.cpp
@@ -142,7 +142,6 @@ int main(int argc, char * const argv[]) {
   /*   An edge [11,6]         */
   /*   An edge [10,12,2]      */
 
-  st.set_dimension(2);
 
   std::cout << "The complex contains " << st.num_simplices() << " simplices - " << st.num_vertices() << " vertices "
       << std::endl;
diff --git a/example/Persistent_cohomology/plain_homology.cpp b/example/Persistent_cohomology/plain_homology.cpp
index 50f692f2..a5ae09c8 100644
--- a/example/Persistent_cohomology/plain_homology.cpp
+++ b/example/Persistent_cohomology/plain_homology.cpp
@@ -64,8 +64,6 @@ int main() {
   st.insert_simplex_and_subfaces(edge03);
   st.insert_simplex(edge13);
   st.insert_simplex(vertex4);
-  // FIXME: Remove this line
-  st.set_dimension(2);
 
   // Sort the simplices in the order of the filtration
   st.initialize_filtration();
diff --git a/example/Persistent_cohomology/rips_distance_matrix_persistence.cpp b/example/Persistent_cohomology/rips_distance_matrix_persistence.cpp
deleted file mode 100644
index d38808c7..00000000
--- a/example/Persistent_cohomology/rips_distance_matrix_persistence.cpp
+++ /dev/null
@@ -1,144 +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):       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/>.
- */
-
-#include <gudhi/Rips_complex.h>
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/reader_utils.h>
-
-#include <boost/program_options.hpp>
-
-#include <string>
-#include <vector>
-#include <limits>  // infinity
-
-// Types definition
-using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = Simplex_tree::Filtration_value;
-using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
-using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
-using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >;
-using Distance_matrix = std::vector<std::vector<Filtration_value>>;
-
-void program_options(int argc, char * argv[]
-                     , std::string & csv_matrix_file
-                     , std::string & filediag
-                     , Filtration_value & threshold
-                     , int & dim_max
-                     , int & p
-                     , Filtration_value & min_persistence);
-
-int main(int argc, char * argv[]) {
-  std::string csv_matrix_file;
-  std::string filediag;
-  Filtration_value threshold;
-  int dim_max;
-  int p;
-  Filtration_value min_persistence;
-
-  program_options(argc, argv, csv_matrix_file, filediag, threshold, dim_max, p, min_persistence);
-
-  Distance_matrix distances = Gudhi::read_lower_triangular_matrix_from_csv_file<Filtration_value>(csv_matrix_file);
-  Rips_complex rips_complex_from_file(distances, threshold);
-
-  // Construct the Rips complex in a Simplex Tree
-  Simplex_tree simplex_tree;
-
-  rips_complex_from_file.create_complex(simplex_tree, dim_max);
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(simplex_tree);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(p);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  // Output the diagram in filediag
-  if (filediag.empty()) {
-    pcoh.output_diagram();
-  } else {
-    std::ofstream out(filediag);
-    pcoh.output_diagram(out);
-    out.close();
-  }
-  return 0;
-}
-
-void program_options(int argc, char * argv[]
-                     , std::string & csv_matrix_file
-                     , std::string & filediag
-                     , Filtration_value & threshold
-                     , int & dim_max
-                     , int & p
-                     , Filtration_value & min_persistence) {
-  namespace po = boost::program_options;
-  po::options_description hidden("Hidden options");
-  hidden.add_options()
-      ("input-file", po::value<std::string>(&csv_matrix_file),
-       "Name of file containing a distance matrix. Can be square or lower triangular matrix. Separator is ';'.");
-
-  po::options_description visible("Allowed options", 100);
-  visible.add_options()
-      ("help,h", "produce help message")
-      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
-       "Name of file in which the persistence diagram is written. Default print in std::cout")
-      ("max-edge-length,r",
-       po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
-       "Maximal length of an edge for the Rips complex construction.")
-      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
-       "Maximal dimension of the Rips complex we want to compute.")
-      ("field-charac,p", po::value<int>(&p)->default_value(11),
-       "Characteristic p of the coefficient field Z/pZ for computing homology.")
-      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
-
-  po::positional_options_description pos;
-  pos.add("input-file", 1);
-
-  po::options_description all;
-  all.add(visible).add(hidden);
-
-  po::variables_map vm;
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), vm);
-  po::notify(vm);
-
-  if (vm.count("help") || !vm.count("input-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of a Rips complex defined on a set of distance matrix.\n \n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-    std::cout << visible << std::endl;
-    std::abort();
-  }
-}
diff --git a/example/Persistent_cohomology/rips_persistence.cpp b/example/Persistent_cohomology/rips_persistence.cpp
deleted file mode 100644
index d504798b..00000000
--- a/example/Persistent_cohomology/rips_persistence.cpp
+++ /dev/null
@@ -1,147 +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):       Clément Maria
- *
- *    Copyright (C) 2014  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/Rips_complex.h>
-#include <gudhi/distance_functions.h>
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_off_io.h>
-
-#include <boost/program_options.hpp>
-
-#include <string>
-#include <vector>
-#include <limits>  // infinity
-
-// Types definition
-using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = Simplex_tree::Filtration_value;
-using Rips_complex = Gudhi::rips_complex::Rips_complex<Filtration_value>;
-using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
-using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >;
-using Point = std::vector<double>;
-using Points_off_reader = Gudhi::Points_off_reader<Point>;
-
-void program_options(int argc, char * argv[]
-                     , std::string & off_file_points
-                     , std::string & filediag
-                     , Filtration_value & threshold
-                     , int & dim_max
-                     , int & p
-                     , Filtration_value & min_persistence);
-
-int main(int argc, char * argv[]) {
-  std::string off_file_points;
-  std::string filediag;
-  Filtration_value threshold;
-  int dim_max;
-  int p;
-  Filtration_value min_persistence;
-
-  program_options(argc, argv, off_file_points, filediag, threshold, dim_max, p, min_persistence);
-
-  Points_off_reader off_reader(off_file_points);
-  Rips_complex rips_complex_from_file(off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance());
-
-  // Construct the Rips complex in a Simplex Tree
-  Simplex_tree simplex_tree;
-
-  rips_complex_from_file.create_complex(simplex_tree, dim_max);
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(simplex_tree);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(p);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  // Output the diagram in filediag
-  if (filediag.empty()) {
-    pcoh.output_diagram();
-  } else {
-    std::ofstream out(filediag);
-    pcoh.output_diagram(out);
-    out.close();
-  }
-
-  return 0;
-}
-
-void program_options(int argc, char * argv[]
-                     , std::string & off_file_points
-                     , std::string & filediag
-                     , Filtration_value & threshold
-                     , int & dim_max
-                     , int & p
-                     , Filtration_value & min_persistence) {
-  namespace po = boost::program_options;
-  po::options_description hidden("Hidden options");
-  hidden.add_options()
-      ("input-file", po::value<std::string>(&off_file_points),
-       "Name of an OFF file containing a point set.\n");
-
-  po::options_description visible("Allowed options", 100);
-  visible.add_options()
-      ("help,h", "produce help message")
-      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
-       "Name of file in which the persistence diagram is written. Default print in std::cout")
-      ("max-edge-length,r",
-       po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
-       "Maximal length of an edge for the Rips complex construction.")
-      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
-       "Maximal dimension of the Rips complex we want to compute.")
-      ("field-charac,p", po::value<int>(&p)->default_value(11),
-       "Characteristic p of the coefficient field Z/pZ for computing homology.")
-      ("min-persistence,m", po::value<Filtration_value>(&min_persistence),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals");
-
-  po::positional_options_description pos;
-  pos.add("input-file", 1);
-
-  po::options_description all;
-  all.add(visible).add(hidden);
-
-  po::variables_map vm;
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), vm);
-  po::notify(vm);
-
-  if (vm.count("help") || !vm.count("input-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of a Rips complex defined on a set of input points.\n \n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-    std::cout << visible << std::endl;
-    std::abort();
-  }
-}
diff --git a/example/Persistent_cohomology/rips_persistence_step_by_step.cpp b/example/Persistent_cohomology/rips_persistence_step_by_step.cpp
index 554eeba6..c1de0ef8 100644
--- a/example/Persistent_cohomology/rips_persistence_step_by_step.cpp
+++ b/example/Persistent_cohomology/rips_persistence_step_by_step.cpp
@@ -45,14 +45,7 @@
 using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
 using Vertex_handle = Simplex_tree::Vertex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
-using Graph_t = boost::adjacency_list < boost::vecS, boost::vecS, boost::undirectedS
-, boost::property < vertex_filtration_t, Filtration_value >
-, boost::property < edge_filtration_t, Filtration_value >
->;
-using Edge_t = std::pair< Vertex_handle, Vertex_handle >;
-
-template< typename InputPointRange, typename Distance >
-Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance);
+using Proximity_graph = Gudhi::Proximity_graph<Simplex_tree>;
 
 using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
 using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp >;
@@ -81,8 +74,9 @@ int main(int argc, char * argv[]) {
   Points_off_reader off_reader(off_file_points);
 
   // Compute the proximity graph of the points
-  Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold
-                                               , Gudhi::Euclidean_distance());
+  Proximity_graph prox_graph = Gudhi::compute_proximity_graph<Simplex_tree>(off_reader.get_point_cloud(),
+                                                                            threshold,
+                                                                            Gudhi::Euclidean_distance());
 
   // Construct the Rips complex in a Simplex Tree
   Simplex_tree st;
@@ -170,48 +164,3 @@ void program_options(int argc, char * argv[]
     std::abort();
   }
 }
-
-/** Output the proximity graph of the points.
- *
- * If points contains n elements, the proximity graph is the graph
- * with n vertices, and an edge [u,v] iff the distance function between
- * points u and v is smaller than threshold.
- *
- * The type PointCloud furnishes .begin() and .end() methods, that return
- * iterators with value_type Point.
- */
-template< typename InputPointRange, typename Distance >
-Graph_t compute_proximity_graph(InputPointRange &points, Filtration_value threshold, Distance distance) {
-  std::vector< Edge_t > edges;
-  std::vector< Filtration_value > edges_fil;
-
-  Vertex_handle idx_u, idx_v;
-  Filtration_value fil;
-  idx_u = 0;
-  for (auto it_u = points.begin(); it_u != points.end(); ++it_u) {
-    idx_v = idx_u + 1;
-    for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) {
-      fil = distance(*it_u, *it_v);
-      if (fil <= threshold) {
-        edges.emplace_back(idx_u, idx_v);
-        edges_fil.push_back(fil);
-      }
-    }
-    ++idx_u;
-  }
-
-  Graph_t skel_graph(edges.begin()
-                     , edges.end()
-                     , edges_fil.begin()
-                     , idx_u);  // number of points labeled from 0 to idx_u-1
-
-  auto vertex_prop = boost::get(vertex_filtration_t(), skel_graph);
-
-  boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end;
-  for (std::tie(vi, vi_end) = boost::vertices(skel_graph);
-       vi != vi_end; ++vi) {
-    boost::put(vertex_prop, *vi, 0.);
-  }
-
-  return skel_graph;
-}
diff --git a/example/Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp b/example/Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp
deleted file mode 100644
index 72268511..00000000
--- a/example/Persistent_cohomology/weighted_alpha_complex_3d_persistence.cpp
+++ /dev/null
@@ -1,267 +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):       Vincent Rouvreau
- *
- *    Copyright (C) 2014  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 <boost/variant.hpp>
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_3D_off_io.h>
-
-#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
-#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
-#include <CGAL/Regular_triangulation_3.h>
-#include <CGAL/Alpha_shape_3.h>
-#include <CGAL/iterator.h>
-
-#include <fstream>
-#include <cmath>
-#include <string>
-#include <tuple>
-#include <map>
-#include <utility>
-#include <list>
-#include <vector>
-#include <cstdlib>
-
-#include "alpha_complex_3d_helper.h"
-
-// Traits
-using Kernel = CGAL::Exact_predicates_inexact_constructions_kernel;
-using Gt = CGAL::Regular_triangulation_euclidean_traits_3<Kernel>;
-using Vb = CGAL::Alpha_shape_vertex_base_3<Gt>;
-using Fb = CGAL::Alpha_shape_cell_base_3<Gt>;
-using Tds = CGAL::Triangulation_data_structure_3<Vb, Fb>;
-using Triangulation_3 = CGAL::Regular_triangulation_3<Gt, Tds>;
-using Alpha_shape_3 = CGAL::Alpha_shape_3<Triangulation_3>;
-
-// From file type definition
-using Point_3 = Gt::Bare_point;
-using Weighted_point_3 = Gt::Weighted_point;
-
-// filtration with alpha values needed type definition
-using Alpha_value_type = Alpha_shape_3::FT;
-using Object = CGAL::Object;
-using Dispatch = CGAL::Dispatch_output_iterator<
-          CGAL::cpp11::tuple<Object, Alpha_value_type>,
-          CGAL::cpp11::tuple<std::back_insert_iterator< std::vector<Object> >,
-          std::back_insert_iterator< std::vector<Alpha_value_type> > > >;
-using Cell_handle = Alpha_shape_3::Cell_handle;
-using Facet = Alpha_shape_3::Facet;
-using Edge_3 = Alpha_shape_3::Edge;
-using Vertex_handle = Alpha_shape_3::Vertex_handle;
-using Vertex_list = std::list<Alpha_shape_3::Vertex_handle>;
-
-// gudhi type definition
-using ST = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
-using Filtration_value = ST::Filtration_value;
-using Simplex_tree_vertex = ST::Vertex_handle;
-using Alpha_shape_simplex_tree_map = std::map<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex >;
-using Alpha_shape_simplex_tree_pair = std::pair<Alpha_shape_3::Vertex_handle, Simplex_tree_vertex>;
-using Simplex_tree_vector_vertex = std::vector< Simplex_tree_vertex >;
-using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<
-    ST, Gudhi::persistent_cohomology::Field_Zp >;
-
-void usage(char * const progName) {
-  std::cerr << "Usage:\n" << progName << " path_to_OFF_file path_to_weight_file coeff_field_characteristic[integer " <<
-               "> 0] min_persistence[float >= -1.0]\n";
-  std::cerr << "  path_to_OFF_file is the path to your points cloud in OFF format.\n";
-  std::cerr << "  path_to_weight_file is the path to the weights of your points cloud (one value per line.)\n";
-  std::cerr << "  Weights values are explained on CGAL documentation:\n";
-  std::cerr << "    https://doc.cgal.org/latest/Alpha_shapes_3/index.html#title0\n";
-  std::cerr << "    https://doc.cgal.org/latest/Triangulation_3/index.html#Triangulation3secclassRegulartriangulation\n";
-  exit(-1);
-}
-
-int main(int argc, char * const argv[]) {
-  // program args management
-  if (argc != 5) {
-    std::cerr << "Error: Number of arguments (" << argc << ") is not correct\n";
-    usage(argv[0]);
-  }
-
-  int coeff_field_characteristic = atoi(argv[3]);
-  Filtration_value min_persistence = strtof(argv[4], nullptr);
-
-  // Read points from file
-  std::string offInputFile(argv[1]);
-  // Read the OFF file (input file name given as parameter) and triangulate points
-  Gudhi::Points_3D_off_reader<Point_3> off_reader(offInputFile);
-  // Check the read operation was correct
-  if (!off_reader.is_valid()) {
-    std::cerr << "Unable to read file " << offInputFile << std::endl;
-    usage(argv[0]);
-  }
-
-  // Retrieve the triangulation
-  std::vector<Point_3> lp = off_reader.get_point_cloud();
-
-  // Read weights information from file
-  std::ifstream weights_ifstr(argv[2]);
-  std::vector<Weighted_point_3> wp;
-  if (weights_ifstr.good()) {
-    double weight = 0.0;
-    std::size_t index = 0;
-    wp.reserve(lp.size());
-    // Attempt read the weight in a double format, return false if it fails
-    while ((weights_ifstr >> weight) && (index < lp.size())) {
-      wp.push_back(Weighted_point_3(lp[index], weight));
-      index++;
-    }
-    if (index != lp.size()) {
-      std::cerr << "Bad number of weights in file " << argv[2] << std::endl;
-      usage(argv[0]);
-    }
-  } else {
-    std::cerr << "Unable to read file " << argv[2] << std::endl;
-    usage(argv[0]);
-  }
-
-  // alpha shape construction from points. CGAL has a strange behavior in REGULARIZED mode.
-  Alpha_shape_3 as(wp.begin(), wp.end(), 0, Alpha_shape_3::GENERAL);
-#ifdef DEBUG_TRACES
-  std::cout << "Alpha shape computed in GENERAL mode" << std::endl;
-#endif  // DEBUG_TRACES
-
-  // filtration with alpha values from alpha shape
-  std::vector<Object> the_objects;
-  std::vector<Alpha_value_type> the_alpha_values;
-
-  Dispatch disp = CGAL::dispatch_output<Object, Alpha_value_type>(std::back_inserter(the_objects),
-                                                                  std::back_inserter(the_alpha_values));
-
-  as.filtration_with_alpha_values(disp);
-#ifdef DEBUG_TRACES
-  std::cout << "filtration_with_alpha_values returns : " << the_objects.size() << " objects" << std::endl;
-#endif  // DEBUG_TRACES
-
-  Alpha_shape_3::size_type count_vertices = 0;
-  Alpha_shape_3::size_type count_edges = 0;
-  Alpha_shape_3::size_type count_facets = 0;
-  Alpha_shape_3::size_type count_cells = 0;
-
-  // Loop on objects vector
-  Vertex_list vertex_list;
-  ST simplex_tree;
-  Alpha_shape_simplex_tree_map map_cgal_simplex_tree;
-  std::vector<Alpha_value_type>::iterator the_alpha_value_iterator = the_alpha_values.begin();
-  int dim_max = 0;
-  Filtration_value filtration_max = 0.0;
-  for (auto object_iterator : the_objects) {
-    // Retrieve Alpha shape vertex list from object
-    if (const Cell_handle * cell = CGAL::object_cast<Cell_handle>(&object_iterator)) {
-      vertex_list = from_cell<Vertex_list, Cell_handle>(*cell);
-      count_cells++;
-      if (dim_max < 3) {
-        // Cell is of dim 3
-        dim_max = 3;
-      }
-    } else if (const Facet * facet = CGAL::object_cast<Facet>(&object_iterator)) {
-      vertex_list = from_facet<Vertex_list, Facet>(*facet);
-      count_facets++;
-      if (dim_max < 2) {
-        // Facet is of dim 2
-        dim_max = 2;
-      }
-    } else if (const Edge_3 * edge = CGAL::object_cast<Edge_3>(&object_iterator)) {
-      vertex_list = from_edge<Vertex_list, Edge_3>(*edge);
-      count_edges++;
-      if (dim_max < 1) {
-        // Edge_3 is of dim 1
-        dim_max = 1;
-      }
-    } else if (const Alpha_shape_3::Vertex_handle * vertex =
-               CGAL::object_cast<Alpha_shape_3::Vertex_handle>(&object_iterator)) {
-      count_vertices++;
-      vertex_list = from_vertex<Vertex_list, Vertex_handle>(*vertex);
-    }
-    // Construction of the vector of simplex_tree vertex from list of alpha_shapes vertex
-    Simplex_tree_vector_vertex the_simplex_tree;
-    for (auto the_alpha_shape_vertex : vertex_list) {
-      Alpha_shape_simplex_tree_map::iterator the_map_iterator = map_cgal_simplex_tree.find(the_alpha_shape_vertex);
-      if (the_map_iterator == map_cgal_simplex_tree.end()) {
-        // alpha shape not found
-        Simplex_tree_vertex vertex = map_cgal_simplex_tree.size();
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] not found - insert " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-        map_cgal_simplex_tree.insert(Alpha_shape_simplex_tree_pair(the_alpha_shape_vertex, vertex));
-      } else {
-        // alpha shape found
-        Simplex_tree_vertex vertex = the_map_iterator->second;
-#ifdef DEBUG_TRACES
-        std::cout << "vertex [" << the_alpha_shape_vertex->point() << "] found in " << vertex << std::endl;
-#endif  // DEBUG_TRACES
-        the_simplex_tree.push_back(vertex);
-      }
-    }
-    // Construction of the simplex_tree
-    Filtration_value filtr = /*std::sqrt*/(*the_alpha_value_iterator);
-#ifdef DEBUG_TRACES
-    std::cout << "filtration = " << filtr << std::endl;
-#endif  // DEBUG_TRACES
-    if (filtr > filtration_max) {
-      filtration_max = filtr;
-    }
-    simplex_tree.insert_simplex(the_simplex_tree, filtr);
-    if (the_alpha_value_iterator != the_alpha_values.end())
-      ++the_alpha_value_iterator;
-    else
-      std::cout << "This shall not happen" << std::endl;
-  }
-  simplex_tree.set_dimension(dim_max);
-
-#ifdef DEBUG_TRACES
-  std::cout << "vertices \t\t" << count_vertices << std::endl;
-  std::cout << "edges \t\t" << count_edges << std::endl;
-  std::cout << "facets \t\t" << count_facets << std::endl;
-  std::cout << "cells \t\t" << count_cells << std::endl;
-
-
-  std::cout << "Information of the Simplex Tree: " << std::endl;
-  std::cout << "  Number of vertices = " << simplex_tree.num_vertices() << " ";
-  std::cout << "  Number of simplices = " << simplex_tree.num_simplices() << std::endl << std::endl;
-  std::cout << "  Dimension = " << simplex_tree.dimension() << " ";
-#endif  // DEBUG_TRACES
-
-#ifdef DEBUG_TRACES
-  std::cout << "Iterator on vertices: " << std::endl;
-  for (auto vertex : simplex_tree.complex_vertex_range()) {
-    std::cout << vertex << " ";
-  }
-#endif  // DEBUG_TRACES
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  std::cout << "Simplex_tree dim: " << simplex_tree.dimension() << std::endl;
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(simplex_tree, true);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(coeff_field_characteristic);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  pcoh.output_diagram();
-
-  return 0;
-}
diff --git a/example/Simplex_tree/CMakeLists.txt b/example/Simplex_tree/CMakeLists.txt
index e22cc92c..b33b2d05 100644
--- a/example/Simplex_tree/CMakeLists.txt
+++ b/example/Simplex_tree/CMakeLists.txt
@@ -35,4 +35,21 @@ if(GMP_FOUND AND CGAL_FOUND)
 
   install(TARGETS Simplex_tree_example_alpha_shapes_3_from_off DESTINATION bin)
 
+  add_executable ( Simplex_tree_example_cech_complex_cgal_mini_sphere_3d cech_complex_cgal_mini_sphere_3d.cpp )
+  target_link_libraries(Simplex_tree_example_cech_complex_cgal_mini_sphere_3d ${Boost_PROGRAM_OPTIONS_LIBRARY} ${CGAL_LIBRARY})
+  if (TBB_FOUND)
+    target_link_libraries(Simplex_tree_example_cech_complex_cgal_mini_sphere_3d ${TBB_LIBRARIES})
+  endif()
+  add_test(NAME Simplex_tree_example_cech_complex_cgal_mini_sphere_3d COMMAND $<TARGET_FILE:Simplex_tree_example_cech_complex_cgal_mini_sphere_3d>
+      "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" -r 0.3 -d 3)
+
+  install(TARGETS Simplex_tree_example_alpha_shapes_3_from_off DESTINATION bin)
 endif()
+
+add_executable ( Simplex_tree_example_graph_expansion_with_blocker graph_expansion_with_blocker.cpp )
+if (TBB_FOUND)
+  target_link_libraries(Simplex_tree_example_graph_expansion_with_blocker ${TBB_LIBRARIES})
+endif()
+add_test(NAME Simplex_tree_example_graph_expansion_with_blocker COMMAND $<TARGET_FILE:Simplex_tree_example_graph_expansion_with_blocker>)
+
+install(TARGETS Simplex_tree_example_graph_expansion_with_blocker DESTINATION bin)
diff --git a/example/Simplex_tree/cech_complex_cgal_mini_sphere_3d.cpp b/example/Simplex_tree/cech_complex_cgal_mini_sphere_3d.cpp
new file mode 100644
index 00000000..9bd51106
--- /dev/null
+++ b/example/Simplex_tree/cech_complex_cgal_mini_sphere_3d.cpp
@@ -0,0 +1,221 @@
+/*    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
+ *
+ *    Copyright (C) 2014  INRIA Sophia Antipolis-Méditerranée (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/>.
+ */
+
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Points_off_io.h>
+
+#include <CGAL/Epick_d.h>
+#include <CGAL/Min_sphere_of_spheres_d.h>
+#include <CGAL/Min_sphere_of_points_d_traits_d.h>
+
+#include <boost/program_options.hpp>
+
+#include <string>
+#include <vector>
+#include <limits>   // infinity
+#include <utility>  // for pair
+#include <map>
+
+// -------------------------------------------------------------------------------
+// cech_complex_cgal_mini_sphere_3d is an example of each step that is required to
+// build a Cech over a Simplex_tree. Please refer to cech_persistence to see
+// how to do the same thing with the Cech_complex wrapper for less detailed
+// steps.
+// -------------------------------------------------------------------------------
+
+// Types definition
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using Vertex_handle = Simplex_tree::Vertex_handle;
+using Simplex_handle = Simplex_tree::Simplex_handle;
+using Filtration_value = Simplex_tree::Filtration_value;
+using Siblings = Simplex_tree::Siblings;
+using Graph_t = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+                                      boost::property<Gudhi::vertex_filtration_t, Filtration_value>,
+                                      boost::property<Gudhi::edge_filtration_t, Filtration_value> >;
+using Edge_t = std::pair<Vertex_handle, Vertex_handle>;
+
+using Kernel = CGAL::Epick_d<CGAL::Dimension_tag<3> >;
+using Point = Kernel::Point_d;
+using Traits = CGAL::Min_sphere_of_points_d_traits_d<Kernel, Filtration_value, 3>;
+using Min_sphere = CGAL::Min_sphere_of_spheres_d<Traits>;
+
+using Points_off_reader = Gudhi::Points_off_reader<Point>;
+
+class Cech_blocker {
+ public:
+  bool operator()(Simplex_handle sh) {
+    std::vector<Point> points;
+#if DEBUG_TRACES
+    std::cout << "Cech_blocker on [";
+#endif  // DEBUG_TRACES
+    for (auto vertex : simplex_tree_.simplex_vertex_range(sh)) {
+      points.push_back(point_cloud_[vertex]);
+#if DEBUG_TRACES
+      std::cout << vertex << ", ";
+#endif  // DEBUG_TRACES
+    }
+    Min_sphere ms(points.begin(), points.end());
+    Filtration_value radius = ms.radius();
+#if DEBUG_TRACES
+    std::cout << "] - radius = " << radius << " - returns " << (radius > threshold_) << std::endl;
+#endif  // DEBUG_TRACES
+    simplex_tree_.assign_filtration(sh, radius);
+    return (radius > threshold_);
+  }
+  Cech_blocker(Simplex_tree& simplex_tree, Filtration_value threshold, const std::vector<Point>& point_cloud)
+      : simplex_tree_(simplex_tree), threshold_(threshold), point_cloud_(point_cloud) {}
+
+ private:
+  Simplex_tree simplex_tree_;
+  Filtration_value threshold_;
+  std::vector<Point> point_cloud_;
+};
+
+template <typename InputPointRange>
+Graph_t compute_proximity_graph(InputPointRange& points, Filtration_value threshold);
+
+void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& threshold, int& dim_max);
+
+int main(int argc, char* argv[]) {
+  std::string off_file_points;
+  Filtration_value threshold;
+  int dim_max;
+
+  program_options(argc, argv, off_file_points, threshold, dim_max);
+
+  // Extract the points from the file filepoints
+  Points_off_reader off_reader(off_file_points);
+
+  // Compute the proximity graph of the points
+  Graph_t prox_graph = compute_proximity_graph(off_reader.get_point_cloud(), threshold);
+
+  // Min_sphere sph1(off_reader.get_point_cloud()[0], off_reader.get_point_cloud()[1], off_reader.get_point_cloud()[2]);
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree st;
+  // insert the proximity graph in the simplex tree
+  st.insert_graph(prox_graph);
+  // expand the graph until dimension dim_max
+  st.expansion_with_blockers(dim_max, Cech_blocker(st, threshold, off_reader.get_point_cloud()));
+
+  std::cout << "The complex contains " << st.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << st.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  st.initialize_filtration();
+
+#if DEBUG_TRACES
+  std::cout << "********************************************************************\n";
+  // Display the Simplex_tree - Can not be done in the middle of 2 inserts
+  std::cout << "* The complex contains " << st.num_simplices() << " simplices - dimension=" << st.dimension() << "\n";
+  std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
+  for (auto f_simplex : st.filtration_simplex_range()) {
+    std::cout << "   "
+              << "[" << st.filtration(f_simplex) << "] ";
+    for (auto vertex : st.simplex_vertex_range(f_simplex)) {
+      std::cout << static_cast<int>(vertex) << " ";
+    }
+    std::cout << std::endl;
+  }
+#endif  // DEBUG_TRACES
+  return 0;
+}
+
+void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& threshold, int& dim_max) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of an OFF file containing a 3d point set.\n");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "max-edge-length,r",
+      po::value<Filtration_value>(&threshold)->default_value(std::numeric_limits<Filtration_value>::infinity()),
+      "Maximal length of an edge for the Cech complex construction.")(
+      "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+      "Maximal dimension of the Cech complex we want to compute.");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  if (vm.count("help") || !vm.count("input-file")) {
+    std::cout << std::endl;
+    std::cout << "Construct a Cech complex defined on a set of input points.\n \n";
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
+
+/** Output the proximity graph of the points.
+ *
+ * If points contains n elements, the proximity graph is the graph
+ * with n vertices, and an edge [u,v] iff the distance function between
+ * points u and v is smaller than threshold.
+ *
+ * The type PointCloud furnishes .begin() and .end() methods, that return
+ * iterators with value_type Point.
+ */
+template <typename InputPointRange>
+Graph_t compute_proximity_graph(InputPointRange& points, Filtration_value threshold) {
+  std::vector<Edge_t> edges;
+  std::vector<Filtration_value> edges_fil;
+
+  Kernel k;
+  Vertex_handle idx_u, idx_v;
+  Filtration_value fil;
+  idx_u = 0;
+  for (auto it_u = points.begin(); it_u != points.end(); ++it_u) {
+    idx_v = idx_u + 1;
+    for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) {
+      fil = k.squared_distance_d_object()(*it_u, *it_v);
+      // For Cech Complex, threshold is a radius (distance /2)
+      fil = std::sqrt(fil) / 2.;
+      if (fil <= threshold) {
+        edges.emplace_back(idx_u, idx_v);
+        edges_fil.push_back(fil);
+      }
+    }
+    ++idx_u;
+  }
+
+  Graph_t skel_graph(edges.begin(), edges.end(), edges_fil.begin(),
+                     idx_u);  // number of points labeled from 0 to idx_u-1
+
+  auto vertex_prop = boost::get(Gudhi::vertex_filtration_t(), skel_graph);
+
+  boost::graph_traits<Graph_t>::vertex_iterator vi, vi_end;
+  for (std::tie(vi, vi_end) = boost::vertices(skel_graph); vi != vi_end; ++vi) {
+    boost::put(vertex_prop, *vi, 0.);
+  }
+
+  return skel_graph;
+}
diff --git a/example/Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp b/example/Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp
index ff2eebcb..d8289ba9 100644
--- a/example/Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp
+++ b/example/Simplex_tree/example_alpha_shapes_3_simplex_tree_from_off_file.cpp
@@ -28,6 +28,8 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Delaunay_triangulation_3.h>
 #include <CGAL/Alpha_shape_3.h>
+#include <CGAL/Alpha_shape_vertex_base_3.h>
+#include <CGAL/Alpha_shape_cell_base_3.h>
 #include <CGAL/iterator.h>
 
 #include <fstream>
diff --git a/example/Simplex_tree/graph_expansion_with_blocker.cpp b/example/Simplex_tree/graph_expansion_with_blocker.cpp
new file mode 100644
index 00000000..0d458cbd
--- /dev/null
+++ b/example/Simplex_tree/graph_expansion_with_blocker.cpp
@@ -0,0 +1,77 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2014
+ *
+ *    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/Simplex_tree.h>
+
+#include <iostream>
+
+using Simplex_tree = Gudhi::Simplex_tree<>;
+using Simplex_handle = Simplex_tree::Simplex_handle;
+
+int main(int argc, char* const argv[]) {
+  // Construct the Simplex Tree with a 1-skeleton graph example
+  Simplex_tree simplexTree;
+
+  simplexTree.insert_simplex({0, 1}, 0.);
+  simplexTree.insert_simplex({0, 2}, 1.);
+  simplexTree.insert_simplex({0, 3}, 2.);
+  simplexTree.insert_simplex({1, 2}, 3.);
+  simplexTree.insert_simplex({1, 3}, 4.);
+  simplexTree.insert_simplex({2, 3}, 5.);
+  simplexTree.insert_simplex({2, 4}, 6.);
+  simplexTree.insert_simplex({3, 6}, 7.);
+  simplexTree.insert_simplex({4, 5}, 8.);
+  simplexTree.insert_simplex({4, 6}, 9.);
+  simplexTree.insert_simplex({5, 6}, 10.);
+  simplexTree.insert_simplex({6}, 10.);
+
+  simplexTree.expansion_with_blockers(3, [&](Simplex_handle sh) {
+    bool result = false;
+    std::cout << "Blocker on [";
+    // User can loop on the vertices from the given simplex_handle i.e.
+    for (auto vertex : simplexTree.simplex_vertex_range(sh)) {
+      // We block the expansion, if the vertex '6' is in the given list of vertices
+      if (vertex == 6) result = true;
+      std::cout << vertex << ", ";
+    }
+    std::cout << "] ( " << simplexTree.filtration(sh);
+    // User can re-assign a new filtration value directly in the blocker (default is the maximal value of boudaries)
+    simplexTree.assign_filtration(sh, simplexTree.filtration(sh) + 1.);
+
+    std::cout << " + 1. ) = " << result << std::endl;
+
+    return result;
+  });
+
+  std::cout << "********************************************************************\n";
+  std::cout << "* The complex contains " << simplexTree.num_simplices() << " simplices";
+  std::cout << "   - dimension " << simplexTree.dimension() << "\n";
+  std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
+  for (auto f_simplex : simplexTree.filtration_simplex_range()) {
+    std::cout << "   "
+              << "[" << simplexTree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")";
+    std::cout << std::endl;
+  }
+
+  return 0;
+}
diff --git a/example/Simplex_tree/mini_simplex_tree.cpp b/example/Simplex_tree/mini_simplex_tree.cpp
index ad99df23..19e45361 100644
--- a/example/Simplex_tree/mini_simplex_tree.cpp
+++ b/example/Simplex_tree/mini_simplex_tree.cpp
@@ -52,8 +52,6 @@ int main() {
   auto edge03 = {0, 3};
   st.insert_simplex_and_subfaces(triangle012);
   st.insert_simplex_and_subfaces(edge03);
-  // FIXME: Remove this line
-  st.set_dimension(2);
 
   auto edge02 = {0, 2};
   ST::Simplex_handle e = st.find(edge02);
diff --git a/example/Simplex_tree/simple_simplex_tree.cpp b/example/Simplex_tree/simple_simplex_tree.cpp
index d8318f03..828977c2 100644
--- a/example/Simplex_tree/simple_simplex_tree.cpp
+++ b/example/Simplex_tree/simple_simplex_tree.cpp
@@ -30,10 +30,10 @@
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Vertex_handle = Simplex_tree::Vertex_handle;
 using Filtration_value = Simplex_tree::Filtration_value;
-using typeVectorVertex = std::vector< Vertex_handle >;
-using typePairSimplexBool = std::pair< Simplex_tree::Simplex_handle, bool >;
+using typeVectorVertex = std::vector<Vertex_handle>;
+using typePairSimplexBool = std::pair<Simplex_tree::Simplex_handle, bool>;
 
-int main(int argc, char * const argv[]) {
+int main(int argc, char* const argv[]) {
   const Filtration_value FIRST_FILTRATION_VALUE = 0.1;
   const Filtration_value SECOND_FILTRATION_VALUE = 0.2;
   const Filtration_value THIRD_FILTRATION_VALUE = 0.3;
@@ -54,7 +54,7 @@ int main(int argc, char * const argv[]) {
 
   // ++ FIRST
   std::cout << "   * INSERT 0" << std::endl;
-  typeVectorVertex firstSimplexVector = { 0 };
+  typeVectorVertex firstSimplexVector = {0};
   typePairSimplexBool returnValue =
       simplexTree.insert_simplex(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
@@ -66,9 +66,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ SECOND
   std::cout << "   * INSERT 1" << std::endl;
-  typeVectorVertex secondSimplexVector = { 1 };
-  returnValue =
-      simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typeVectorVertex secondSimplexVector = {1};
+  returnValue = simplexTree.insert_simplex(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + 1 INSERTED" << std::endl;
@@ -78,9 +77,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ THIRD
   std::cout << "   * INSERT (0,1)" << std::endl;
-  typeVectorVertex thirdSimplexVector = { 0, 1 };
-  returnValue =
-      simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex thirdSimplexVector = {0, 1};
+  returnValue = simplexTree.insert_simplex(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (0,1) INSERTED" << std::endl;
@@ -90,9 +88,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ FOURTH
   std::cout << "   * INSERT 2" << std::endl;
-  typeVectorVertex fourthSimplexVector = { 2 };
-  returnValue =
-      simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typeVectorVertex fourthSimplexVector = {2};
+  returnValue = simplexTree.insert_simplex(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + 2 INSERTED" << std::endl;
@@ -102,9 +99,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ FIFTH
   std::cout << "   * INSERT (2,0)" << std::endl;
-  typeVectorVertex fifthSimplexVector = { 2, 0 };
-  returnValue =
-      simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex fifthSimplexVector = {2, 0};
+  returnValue = simplexTree.insert_simplex(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,0) INSERTED" << std::endl;
@@ -114,9 +110,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ SIXTH
   std::cout << "   * INSERT (2,1)" << std::endl;
-  typeVectorVertex sixthSimplexVector = { 2, 1 };
-  returnValue =
-      simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex sixthSimplexVector = {2, 1};
+  returnValue = simplexTree.insert_simplex(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,1) INSERTED" << std::endl;
@@ -126,9 +121,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ SEVENTH
   std::cout << "   * INSERT (2,1,0)" << std::endl;
-  typeVectorVertex seventhSimplexVector = { 2, 1, 0 };
-  returnValue =
-      simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
+  typeVectorVertex seventhSimplexVector = {2, 1, 0};
+  returnValue = simplexTree.insert_simplex(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,1,0) INSERTED" << std::endl;
@@ -138,9 +132,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ EIGHTH
   std::cout << "   * INSERT 3" << std::endl;
-  typeVectorVertex eighthSimplexVector = { 3 };
-  returnValue =
-      simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
+  typeVectorVertex eighthSimplexVector = {3};
+  returnValue = simplexTree.insert_simplex(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + 3 INSERTED" << std::endl;
@@ -150,9 +143,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ NINETH
   std::cout << "   * INSERT (3,0)" << std::endl;
-  typeVectorVertex ninethSimplexVector = { 3, 0 };
-  returnValue =
-      simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
+  typeVectorVertex ninethSimplexVector = {3, 0};
+  returnValue = simplexTree.insert_simplex(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (3,0) INSERTED" << std::endl;
@@ -162,7 +154,7 @@ int main(int argc, char * const argv[]) {
 
   // ++ TENTH
   std::cout << "   * INSERT 0 (already inserted)" << std::endl;
-  typeVectorVertex tenthSimplexVector = { 0 };
+  typeVectorVertex tenthSimplexVector = {0};
   // With a different filtration value
   returnValue = simplexTree.insert_simplex(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
 
@@ -174,9 +166,8 @@ int main(int argc, char * const argv[]) {
 
   // ++ ELEVENTH
   std::cout << "   * INSERT (2,1,0) (already inserted)" << std::endl;
-  typeVectorVertex eleventhSimplexVector = { 2, 1, 0 };
-  returnValue =
-      simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
+  typeVectorVertex eleventhSimplexVector = {2, 1, 0};
+  returnValue = simplexTree.insert_simplex(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
 
   if (returnValue.second == true) {
     std::cout << "   + (2,1,0) INSERTED" << std::endl;
@@ -185,7 +176,6 @@ int main(int argc, char * const argv[]) {
   }
 
   // ++ GENERAL VARIABLE SET
-  simplexTree.set_dimension(2);  // Max dimension = 2 -> (2,1,0)
 
   std::cout << "********************************************************************\n";
   // Display the Simplex_tree - Can not be done in the middle of 2 inserts
@@ -193,10 +183,9 @@ int main(int argc, char * const argv[]) {
   std::cout << "   - dimension " << simplexTree.dimension() << "\n";
   std::cout << "* Iterator on Simplices in the filtration, with [filtration value]:\n";
   for (auto f_simplex : simplexTree.filtration_simplex_range()) {
-    std::cout << "   " << "[" << simplexTree.filtration(f_simplex) << "] ";
-    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) {
-      std::cout << static_cast<int>(vertex) << " ";
-    }
+    std::cout << "   "
+              << "[" << simplexTree.filtration(f_simplex) << "] ";
+    for (auto vertex : simplexTree.simplex_vertex_range(f_simplex)) std::cout << "(" << vertex << ")";
     std::cout << std::endl;
   }
   //   [0.1] 0
@@ -219,7 +208,7 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  typeVectorVertex unknownSimplexVector = { 15 };
+  typeVectorVertex unknownSimplexVector = {15};
   simplexFound = simplexTree.find(unknownSimplexVector);
   std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
@@ -234,7 +223,7 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  typeVectorVertex otherSimplexVector = { 1, 15 };
+  typeVectorVertex otherSimplexVector = {1, 15};
   simplexFound = simplexTree.find(otherSimplexVector);
   std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
@@ -242,12 +231,38 @@ int main(int argc, char * const argv[]) {
   else
     std::cout << "***- NO IT ISN'T\n";
 
-  typeVectorVertex invSimplexVector = { 1, 2, 0 };
+  typeVectorVertex invSimplexVector = {1, 2, 0};
   simplexFound = simplexTree.find(invSimplexVector);
   std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != simplexTree.null_simplex())
     std::cout << "***+ YES IT IS!\n";
   else
     std::cout << "***- NO IT ISN'T\n";
+
+  simplexFound = simplexTree.find({0, 1});
+  std::cout << "**************IS THE SIMPLEX {0,1} IN THE SIMPLEX TREE ?\n";
+  if (simplexFound != simplexTree.null_simplex())
+    std::cout << "***+ YES IT IS!\n";
+  else
+    std::cout << "***- NO IT ISN'T\n";
+
+  std::cout << "**************COFACES OF {0,1} IN CODIMENSION 1 ARE\n";
+  for (auto& simplex : simplexTree.cofaces_simplex_range(simplexTree.find({0, 1}), 1)) {
+    for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")";
+    std::cout << std::endl;
+  }
+
+  std::cout << "**************STARS OF {0,1} ARE\n";
+  for (auto& simplex : simplexTree.star_simplex_range(simplexTree.find({0, 1}))) {
+    for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")";
+    std::cout << std::endl;
+  }
+
+  std::cout << "**************BOUNDARIES OF {0,1,2} ARE\n";
+  for (auto& simplex : simplexTree.boundary_simplex_range(simplexTree.find({0, 1, 2}))) {
+    for (auto vertex : simplexTree.simplex_vertex_range(simplex)) std::cout << "(" << vertex << ")";
+    std::cout << std::endl;
+  }
+
   return 0;
 }
diff --git a/example/Witness_complex/CMakeLists.txt b/example/Witness_complex/CMakeLists.txt
index cbc53902..a8231392 100644
--- a/example/Witness_complex/CMakeLists.txt
+++ b/example/Witness_complex/CMakeLists.txt
@@ -13,39 +13,23 @@ install(TARGETS Witness_complex_example_nearest_landmark_table DESTINATION bin)
 # CGAL and Eigen3 are required for Euclidean version of Witness
 if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0)
   add_executable( Witness_complex_example_off example_witness_complex_off.cpp )
-  add_executable( Witness_complex_example_strong_off example_strong_witness_complex_off.cpp )
   add_executable ( Witness_complex_example_sphere example_witness_complex_sphere.cpp )
-
-  add_executable ( Witness_complex_example_witness_persistence example_witness_complex_persistence.cpp )
-  target_link_libraries(Witness_complex_example_witness_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
-  add_executable ( Witness_complex_example_strong_witness_persistence example_strong_witness_persistence.cpp )
-  target_link_libraries(Witness_complex_example_strong_witness_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
-
-  if (TBB_FOUND)
-    target_link_libraries(Witness_complex_example_witness_persistence ${TBB_LIBRARIES})
-    target_link_libraries(Witness_complex_example_strong_witness_persistence ${TBB_LIBRARIES})
-  endif()
+  
+  add_executable( Witness_complex_example_strong_off example_strong_witness_complex_off.cpp )
+  target_link_libraries(Witness_complex_example_strong_off)
 
   add_test(NAME Witness_complex_example_off_test_torus
       COMMAND $<TARGET_FILE:Witness_complex_example_off>
       "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "20" "1.0" "3")
+  add_test(NAME Witness_complex_example_test_sphere_10
+      COMMAND $<TARGET_FILE:Witness_complex_example_sphere> "10")
   add_test(NAME Witness_complex_example_strong_off_test_torus
       COMMAND $<TARGET_FILE:Witness_complex_example_strong_off>
       "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "20" "1.0" "3")
-  add_test(NAME Witness_complex_example_test_sphere_10
-      COMMAND $<TARGET_FILE:Witness_complex_example_sphere> "10")
-  add_test(NAME Witness_complex_example_test_torus_persistence
-      COMMAND $<TARGET_FILE:Witness_complex_example_witness_persistence>
-      "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5")
-  add_test(NAME Witness_complex_example_strong_test_torus_persistence
-      COMMAND $<TARGET_FILE:Witness_complex_example_strong_witness_persistence>
-      "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-l" "20" "-a" "0.5")
-
+      
   install(TARGETS Witness_complex_example_off DESTINATION bin)
-  install(TARGETS Witness_complex_example_strong_off DESTINATION bin)
   install(TARGETS Witness_complex_example_sphere DESTINATION bin)
-  install(TARGETS Witness_complex_example_witness_persistence DESTINATION bin)
-  install(TARGETS Witness_complex_example_strong_witness_persistence DESTINATION bin)
+  install(TARGETS Witness_complex_example_strong_off DESTINATION bin)
+
 
 endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.6.0)
diff --git a/example/Witness_complex/example_strong_witness_complex_off.cpp b/example/Witness_complex/example_strong_witness_complex_off.cpp
index 0ee9ee90..346bef6d 100644
--- a/example/Witness_complex/example_strong_witness_complex_off.cpp
+++ b/example/Witness_complex/example_strong_witness_complex_off.cpp
@@ -23,6 +23,7 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Euclidean_strong_witness_complex.h>
 #include <gudhi/pick_n_random_points.h>
+#include <gudhi/choose_n_farthest_points.h>
 #include <gudhi/Points_off_io.h>
 
 #include <CGAL/Epick_d.h>
@@ -38,10 +39,9 @@ using Point_d = typename K::Point_d;
 using Witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex<K>;
 using Point_vector = std::vector<Point_d>;
 
-int main(int argc, char * const argv[]) {
+int main(int argc, char* const argv[]) {
   if (argc != 5) {
-    std::cerr << "Usage: " << argv[0]
-        << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n";
+    std::cerr << "Usage: " << argv[0] << " path_to_point_file number_of_landmarks max_squared_alpha limit_dimension\n";
     return 0;
   }
 
@@ -55,25 +55,25 @@ int main(int argc, char * const argv[]) {
   Point_vector point_vector, landmarks;
   Gudhi::Points_off_reader<Point_d> off_reader(file_name);
   if (!off_reader.is_valid()) {
-      std::cerr << "Strong witness complex - Unable to read file " << file_name << "\n";
-      exit(-1);  // ----- >>
-    }
+    std::cerr << "Strong witness complex - Unable to read file " << file_name << "\n";
+    exit(-1);  // ----- >>
+  }
   point_vector = Point_vector(off_reader.get_point_cloud());
 
   std::cout << "Successfully read " << point_vector.size() << " points.\n";
   std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n";
 
-  // Choose landmarks
-  Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
+  // Choose landmarks (decomment one of the following two lines)
+  // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point,
+                                               std::back_inserter(landmarks));
 
   // Compute witness complex
   start = clock();
-  Witness_complex witness_complex(landmarks,
-                                  point_vector);
+  Witness_complex witness_complex(landmarks, point_vector);
 
   witness_complex.create_complex(simplex_tree, alpha2, lim_dim);
   end = clock();
-  std::cout << "Strong witness complex took "
-      << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
+  std::cout << "Strong witness complex took " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
   std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n";
 }
diff --git a/example/Witness_complex/example_strong_witness_persistence.cpp b/example/Witness_complex/example_strong_witness_persistence.cpp
deleted file mode 100644
index f786fe7b..00000000
--- a/example/Witness_complex/example_strong_witness_persistence.cpp
+++ /dev/null
@@ -1,171 +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):       Siargey Kachanovich
- *
- *    Copyright (C) 2016  INRIA (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/>.
- */
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Euclidean_strong_witness_complex.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_off_io.h>
-#include <gudhi/pick_n_random_points.h>
-
-#include <boost/program_options.hpp>
-
-#include <CGAL/Epick_d.h>
-
-#include <string>
-#include <vector>
-#include <limits>  // infinity
-
-using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
-using Point_d = K::Point_d;
-
-using Point_vector = std::vector<Point_d>;
-using Strong_witness_complex = Gudhi::witness_complex::Euclidean_strong_witness_complex<K>;
-using SimplexTree = Gudhi::Simplex_tree<>;
-
-using Filtration_value = SimplexTree::Filtration_value;
-
-using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
-using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<SimplexTree, Field_Zp>;
-
-void program_options(int argc, char * argv[]
-                     , int & nbL
-                     , std::string & file_name
-                     , std::string & filediag
-                     , Filtration_value & max_squared_alpha
-                     , int & p
-                     , int & dim_max
-                     , Filtration_value & min_persistence);
-
-int main(int argc, char * argv[]) {
-  std::string file_name;
-  std::string filediag;
-  Filtration_value max_squared_alpha;
-  int p, nbL, lim_d;
-  Filtration_value min_persistence;
-  SimplexTree simplex_tree;
-
-  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence);
-
-  // Extract the points from the file file_name
-  Point_vector witnesses, landmarks;
-  Gudhi::Points_off_reader<Point_d> off_reader(file_name);
-  if (!off_reader.is_valid()) {
-      std::cerr << "Witness complex - Unable to read file " << file_name << "\n";
-      exit(-1);  // ----- >>
-  }
-  witnesses = Point_vector(off_reader.get_point_cloud());
-  std::cout << "Successfully read " << witnesses.size() << " points.\n";
-  std::cout << "Ambient dimension is " << witnesses[0].dimension() << ".\n";
-
-  // Choose landmarks from witnesses
-  Gudhi::subsampling::pick_n_random_points(witnesses, nbL, std::back_inserter(landmarks));
-
-  // Compute witness complex
-  Strong_witness_complex strong_witness_complex(landmarks,
-                                                witnesses);
-
-  strong_witness_complex.create_complex(simplex_tree, max_squared_alpha, lim_d);
-
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(simplex_tree);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(p);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  // Output the diagram in filediag
-  if (filediag.empty()) {
-    pcoh.output_diagram();
-  } else {
-    std::ofstream out(filediag);
-    pcoh.output_diagram(out);
-    out.close();
-  }
-
-  return 0;
-}
-
-void program_options(int argc, char * argv[]
-                     , int & nbL
-                     , std::string & file_name
-                     , std::string & filediag
-                     , Filtration_value & max_squared_alpha
-                     , int & p
-                     , int & dim_max
-                     , Filtration_value & min_persistence) {
-  namespace po = boost::program_options;
-
-  po::options_description hidden("Hidden options");
-  hidden.add_options()
-      ("input-file", po::value<std::string>(&file_name),
-      "Name of file containing a point set in off format.");
-
-  po::options_description visible("Allowed options", 100);
-  Filtration_value default_alpha = std::numeric_limits<Filtration_value>::infinity();
-  visible.add_options()
-      ("help,h", "produce help message")
-      ("landmarks,l", po::value<int>(&nbL),
-       "Number of landmarks to choose from the point cloud.")
-      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
-       "Name of file in which the persistence diagram is written. Default print in std::cout")
-      ("max-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha),
-       "Maximal squared relaxation parameter.")
-      ("field-charac,p", po::value<int>(&p)->default_value(11),
-       "Characteristic p of the coefficient field Z/pZ for computing homology.")
-      ("min-persistence,m", po::value<Filtration_value>(&min_persistence)->default_value(0),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals")
-      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
-       "Maximal dimension of the strong witness complex we want to compute.");
-
-  po::positional_options_description pos;
-  pos.add("input-file", 1);
-
-  po::options_description all;
-  all.add(visible).add(hidden);
-  po::variables_map vm;
-
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), vm);
-  po::notify(vm);
-
-  if (vm.count("help") || !vm.count("input-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of a Strong witness complex defined on a set of input points.\n \n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-    std::cout << visible << std::endl;
-    std::abort();
-  }
-}
-
diff --git a/example/Witness_complex/example_witness_complex_off.cpp b/example/Witness_complex/example_witness_complex_off.cpp
index b36dac0d..be11c955 100644
--- a/example/Witness_complex/example_witness_complex_off.cpp
+++ b/example/Witness_complex/example_witness_complex_off.cpp
@@ -4,6 +4,7 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Euclidean_witness_complex.h>
 #include <gudhi/pick_n_random_points.h>
+#include <gudhi/choose_n_farthest_points.h>
 #include <gudhi/Points_off_io.h>
 
 #include <CGAL/Epick_d.h>
@@ -44,8 +45,9 @@ int main(int argc, char * const argv[]) {
   std::cout << "Successfully read " << point_vector.size() << " points.\n";
   std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n";
 
-  // Choose landmarks
-  Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
+  // Choose landmarks (decomment one of the following two lines)
+  // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
+  Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks));
 
   // Compute witness complex
   start = clock();
diff --git a/example/Witness_complex/example_witness_complex_persistence.cpp b/example/Witness_complex/example_witness_complex_persistence.cpp
deleted file mode 100644
index a1146922..00000000
--- a/example/Witness_complex/example_witness_complex_persistence.cpp
+++ /dev/null
@@ -1,171 +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):       Siargey Kachanovich
- *
- *    Copyright (C) 2016  INRIA (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/>.
- */
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Euclidean_witness_complex.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Points_off_io.h>
-#include <gudhi/pick_n_random_points.h>
-
-#include <boost/program_options.hpp>
-
-#include <CGAL/Epick_d.h>
-
-#include <string>
-#include <vector>
-#include <limits>  // infinity
-
-using K = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
-using Point_d = K::Point_d;
-
-using Point_vector = std::vector<Point_d>;
-using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<K>;
-using SimplexTree = Gudhi::Simplex_tree<>;
-
-using Filtration_value = SimplexTree::Filtration_value;
-
-using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
-using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<SimplexTree, Field_Zp>;
-
-void program_options(int argc, char * argv[]
-                     , int & nbL
-                     , std::string & file_name
-                     , std::string & filediag
-                     , Filtration_value & max_squared_alpha
-                     , int & p
-                     , int & dim_max
-                     , Filtration_value & min_persistence);
-
-int main(int argc, char * argv[]) {
-  std::string file_name;
-  std::string filediag;
-  Filtration_value max_squared_alpha;
-  int p, nbL, lim_d;
-  Filtration_value min_persistence;
-  SimplexTree simplex_tree;
-
-  program_options(argc, argv, nbL, file_name, filediag, max_squared_alpha, p, lim_d, min_persistence);
-
-  // Extract the points from the file file_name
-  Point_vector witnesses, landmarks;
-  Gudhi::Points_off_reader<Point_d> off_reader(file_name);
-  if (!off_reader.is_valid()) {
-      std::cerr << "Witness complex - Unable to read file " << file_name << "\n";
-      exit(-1);  // ----- >>
-  }
-  witnesses = Point_vector(off_reader.get_point_cloud());
-  std::cout << "Successfully read " << witnesses.size() << " points.\n";
-  std::cout << "Ambient dimension is " << witnesses[0].dimension() << ".\n";
-
-  // Choose landmarks from witnesses
-  Gudhi::subsampling::pick_n_random_points(witnesses, nbL, std::back_inserter(landmarks));
-
-  // Compute witness complex
-  Witness_complex witness_complex(landmarks,
-                                  witnesses);
-
-  witness_complex.create_complex(simplex_tree, max_squared_alpha, lim_d);
-
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
-
-  // Sort the simplices in the order of the filtration
-  simplex_tree.initialize_filtration();
-
-  // Compute the persistence diagram of the complex
-  Persistent_cohomology pcoh(simplex_tree);
-  // initializes the coefficient field for homology
-  pcoh.init_coefficients(p);
-
-  pcoh.compute_persistent_cohomology(min_persistence);
-
-  // Output the diagram in filediag
-  if (filediag.empty()) {
-    pcoh.output_diagram();
-  } else {
-    std::ofstream out(filediag);
-    pcoh.output_diagram(out);
-    out.close();
-  }
-
-  return 0;
-}
-
-
-void program_options(int argc, char * argv[]
-                     , int & nbL
-                     , std::string & file_name
-                     , std::string & filediag
-                     , Filtration_value & max_squared_alpha
-                     , int & p
-                     , int & dim_max
-                     , Filtration_value & min_persistence) {
-  namespace po = boost::program_options;
-
-  po::options_description hidden("Hidden options");
-  hidden.add_options()
-      ("input-file", po::value<std::string>(&file_name),
-      "Name of file containing a point set in off format.");
-
-  Filtration_value default_alpha = std::numeric_limits<Filtration_value>::infinity();
-  po::options_description visible("Allowed options", 100);
-  visible.add_options()
-      ("help,h", "produce help message")
-      ("landmarks,l", po::value<int>(&nbL),
-       "Number of landmarks to choose from the point cloud.")
-      ("output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
-       "Name of file in which the persistence diagram is written. Default print in std::cout")
-      ("max-sq-alpha,a", po::value<Filtration_value>(&max_squared_alpha)->default_value(default_alpha),
-       "Maximal squared relaxation parameter.")
-      ("field-charac,p", po::value<int>(&p)->default_value(11),
-       "Characteristic p of the coefficient field Z/pZ for computing homology.")
-      ("min-persistence,m", po::value<Filtration_value>(&min_persistence)->default_value(0),
-       "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length intervals")
-      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
-       "Maximal dimension of the weak witness complex we want to compute.");
-
-  po::positional_options_description pos;
-  pos.add("input-file", 1);
-
-  po::options_description all;
-  all.add(visible).add(hidden);
-  po::variables_map vm;
-
-  po::store(po::command_line_parser(argc, argv).
-            options(all).positional(pos).run(), vm);
-  po::notify(vm);
-
-  if (vm.count("help") || !vm.count("input-file")) {
-    std::cout << std::endl;
-    std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
-    std::cout << "of a Weak witness complex defined on a set of input points.\n \n";
-    std::cout << "The output diagram contains one bar per line, written with the convention: \n";
-    std::cout << "   p   dim b d \n";
-    std::cout << "where dim is the dimension of the homological feature,\n";
-    std::cout << "b and d are respectively the birth and death of the feature and \n";
-    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
-
-    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
-    std::cout << visible << std::endl;
-    std::abort();
-  }
-}
diff --git a/example/Witness_complex/example_witness_complex_sphere.cpp b/example/Witness_complex/example_witness_complex_sphere.cpp
index 124fd99b..a6e9b11a 100644
--- a/example/Witness_complex/example_witness_complex_sphere.cpp
+++ b/example/Witness_complex/example_witness_complex_sphere.cpp
@@ -25,6 +25,7 @@
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Euclidean_witness_complex.h>
 #include <gudhi/pick_n_random_points.h>
+#include <gudhi/choose_n_farthest_points.h>
 #include <gudhi/reader_utils.h>
 
 #include <CGAL/Epick_d.h>
@@ -41,27 +42,25 @@
 /** Write a gnuplot readable file.
  *  Data range is a random access range of pairs (arg, value)
  */
-template < typename Data_range >
-void write_data(Data_range & data, std::string filename) {
+template <typename Data_range>
+void write_data(Data_range& data, std::string filename) {
   std::ofstream ofs(filename, std::ofstream::out);
-  for (auto entry : data)
-    ofs << entry.first << ", " << entry.second << "\n";
+  for (auto entry : data) ofs << entry.first << ", " << entry.second << "\n";
   ofs.close();
 }
 
-int main(int argc, char * const argv[]) {
+int main(int argc, char* const argv[]) {
   using Kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
   using Witness_complex = Gudhi::witness_complex::Euclidean_witness_complex<Kernel>;
 
   if (argc != 2) {
-    std::cerr << "Usage: " << argv[0]
-        << " number_of_landmarks \n";
+    std::cerr << "Usage: " << argv[0] << " number_of_landmarks \n";
     return 0;
   }
 
   int number_of_landmarks = atoi(argv[1]);
 
-  std::vector< std::pair<int, double> > l_time;
+  std::vector<std::pair<int, double> > l_time;
 
   // Generate points
   for (int nbP = 500; nbP < 10000; nbP += 500) {
@@ -75,16 +74,17 @@ int main(int argc, char * const argv[]) {
 
     // Choose landmarks
     start = clock();
-    Gudhi::subsampling::pick_n_random_points(point_vector, number_of_landmarks, std::back_inserter(landmarks));
+    // Gudhi::subsampling::pick_n_random_points(point_vector, number_of_landmarks, std::back_inserter(landmarks));
+    Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, number_of_landmarks,
+                                                 Gudhi::subsampling::random_starting_point,
+                                                 std::back_inserter(landmarks));
 
     // Compute witness complex
-    Witness_complex witness_complex(landmarks,
-                                    point_vector);
+    Witness_complex witness_complex(landmarks, point_vector);
     witness_complex.create_complex(simplex_tree, 0);
     end = clock();
     double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-    std::cout << "Witness complex for " << number_of_landmarks << " landmarks took "
-        << time << " s. \n";
+    std::cout << "Witness complex for " << number_of_landmarks << " landmarks took " << time << " s. \n";
     std::cout << "Number of simplices is: " << simplex_tree.num_simplices() << "\n";
     l_time.push_back(std::make_pair(nbP, time));
   }