Junk out

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

Former-commit-id: ce2dfd7033a6a8366f9861ced4ed64ac41dfb82e

23 files changed, 21 insertions, 5670 deletions

diff --git a/src/Witness_complex/example/Landmark_choice_random_knn.h b/src/Witness_complex/example/Landmark_choice_random_knn.h
deleted file mode 100644
index fb91e116..00000000
--- a/src/Witness_complex/example/Landmark_choice_random_knn.h
+++ /dev/null
@@ -1,142 +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) 2015  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/>.
- */
-
-#ifndef LANDMARK_CHOICE_BY_RANDOM_KNN_H_
-#define LANDMARK_CHOICE_BY_RANDOM_KNN_H_
-
-#include <utility>  // for pair<>
-#include <vector>
-#include <cstddef> // for ptrdiff_t type
-
-//#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-//#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_incremental_neighbor_search.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-#include <CGAL/Kd_tree.h>
-#include <CGAL/Euclidean_distance.h>
-//#include <CGAL/Kernel_d/Vector_d.h>
-#include <CGAL/Random.h>
-
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
-typedef K::FT FT;
-typedef K::Point_d Point_d;
-typedef CGAL::Search_traits< FT,
-                             Point_d,
-                             typename K::Cartesian_const_iterator_d,
-                             typename K::Construct_cartesian_const_iterator_d > Traits_base;
-typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
-typedef CGAL::Search_traits_adapter< std::ptrdiff_t,
-                                     Point_d*,
-                                     Traits_base> STraits;
-typedef CGAL::Distance_adapter< std::ptrdiff_t,
-                                Point_d*,
-                                Euclidean_distance > Distance_adapter;
-typedef CGAL::Orthogonal_incremental_neighbor_search< STraits,
-                                                      Distance_adapter > Neighbor_search;
-  typedef CGAL::Orthogonal_k_neighbor_search< STraits > Neighbor_search2;
-typedef Neighbor_search::Tree Tree;
-
-  
-  /** \brief Landmark choice strategy by taking random vertices for landmarks.
-   *  \details It chooses nbL distinct landmarks from a random access range `points`
-   *  and outputs a matrix {witness}*{closest landmarks} in knn.
-   */
-  template <typename KNearestNeighbours,
-            typename Point_random_access_range,
-            typename Distance_matrix>
-  void landmark_choice_by_random_knn(Point_random_access_range const & points,
-                                     int nbL,
-                                     FT alpha,
-                                     unsigned limD,
-                                     KNearestNeighbours & knn,
-                                     Distance_matrix & distances) {
-    int nbP = points.end() - points.begin();
-    assert(nbP >= nbL);
-    std::vector<Point_d> landmarks;
-    std::vector<int> landmarks_ind;
-    Point_d p;
-    int chosen_landmark;
-    CGAL::Random rand;
-    // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety
-    int current_number_of_landmarks = 0;  // counter for landmarks
-    for (; current_number_of_landmarks != nbL; current_number_of_landmarks++) {
-      do chosen_landmark = rand.get_int(0,nbP);
-      while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end());
-      p = points[chosen_landmark];
-      landmarks.push_back(p);
-      landmarks_ind.push_back(chosen_landmark);
-    }
-    // std::cout << "Choice finished!" << std::endl;
-    
-    //int dim = points.begin()->size();
-    knn = KNearestNeighbours(nbP);
-    distances = Distance_matrix(nbP);
-    STraits traits(&(landmarks[0]));
-    CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance> adapter(&(landmarks[0]));
-    Euclidean_distance ed;
-    Tree landmark_tree(boost::counting_iterator<std::ptrdiff_t>(0),
-                       boost::counting_iterator<std::ptrdiff_t>(nbL),
-                       typename Tree::Splitter(),
-                       traits);
-    for (int points_i = 0; points_i < nbP; points_i++) {
-      Point_d const & w = points[points_i];
-      Neighbor_search search(landmark_tree,
-                             w,
-                             FT(0),
-                             true,
-                             adapter);
-      Neighbor_search::iterator search_it = search.begin();
-      // Neighbor_search2 search(landmark_tree,
-      //                         w, limD+1,
-      //                         FT(0),
-      //                         true,
-      //                         adapter);
-      // Neighbor_search2::iterator search_it = search.begin();
-      
-      while (knn[points_i].size() < limD) {
-        distances[points_i].push_back(sqrt(search_it->second));
-        knn[points_i].push_back((search_it++)->first);
-      }
-      FT dtow = distances[points_i][limD-1];
-      
-      if (alpha != 0)
-        while (search_it != search.end() && search_it->second < dtow + alpha) {
-          distances[points_i].push_back(sqrt(search_it->second));
-          knn[points_i].push_back((search_it++)->first);
-        }
-      //std::cout << "k = " << knn[points_i].size() << std::endl;
-    }
-  }
-
-}  // namespace witness_complex
-
-}  // namespace Gudhi
-
-#endif  // LANDMARK_CHOICE_BY_RANDOM_POINT_H_
diff --git a/src/Witness_complex/example/Landmark_choice_sparsification.h b/src/Witness_complex/example/Landmark_choice_sparsification.h
deleted file mode 100644
index d5d5fba1..00000000
--- a/src/Witness_complex/example/Landmark_choice_sparsification.h
+++ /dev/null
@@ -1,230 +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) 2015  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/>.
- */
-
-#ifndef LANDMARK_CHOICE_BY_SPARSIFICATION_H_
-#define LANDMARK_CHOICE_BY_SPARSIFICATION_H_
-
-#include <utility>  // for pair<>
-#include <vector>
-#include <cstddef> // for ptrdiff_t type
-#include <algorithm>
-
-//#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-//#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_incremental_neighbor_search.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-#include <CGAL/Kd_tree.h>
-#include <CGAL/Euclidean_distance.h>
-//#include <CGAL/Kernel_d/Vector_d.h>
-#include <CGAL/Random.h>
-#include <CGAL/Fuzzy_sphere.h>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
-  typedef K::FT FT;
-  typedef K::Point_d Point_d;
-  typedef CGAL::Search_traits< FT,
-                               Point_d,
-                               typename K::Cartesian_const_iterator_d,
-                               typename K::Construct_cartesian_const_iterator_d > Traits_base;
-  typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
-  typedef CGAL::Search_traits_adapter< std::ptrdiff_t,
-                                       Point_d*,
-                                       Traits_base> STraits;
-  typedef CGAL::Distance_adapter< std::ptrdiff_t,
-                                  Point_d*,
-                                  Euclidean_distance > Distance_adapter;
-  typedef CGAL::Orthogonal_incremental_neighbor_search< STraits,
-                                                        Distance_adapter > Neighbor_search;
-  typedef CGAL::Orthogonal_k_neighbor_search< STraits > Neighbor_search2;
-  typedef Neighbor_search::Tree Tree;
-  typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
-
-  /** \brief Landmark selection function by as a sub-epsilon-net of the
-   *  given set of points.
-   */
-  template <typename Point_random_access_range>
-  void landmark_choice_by_sparsification(Point_random_access_range & points,
-                                         unsigned nbL,
-                                         FT mu_epsilon,
-                                         Point_random_access_range & landmarks)
-  {
-    unsigned nbP = points.end() - points.begin();
-    assert(nbP >= nbL);
-    CGAL::Random rand;
-    // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety
-    STraits points_traits(&(points[0]));
-    CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance> points_adapter(&(points[0]));
-    std::vector<bool> dropped_points(nbP, false);
-    
-    Tree witness_tree(boost::counting_iterator<std::ptrdiff_t>(0),
-                      boost::counting_iterator<std::ptrdiff_t>(nbP),
-                      typename Tree::Splitter(),
-                      points_traits);
-    
-    for (unsigned points_i = 0; points_i < nbP; points_i++) {
-      if (dropped_points[points_i])
-        continue;
-      Point_d & w = points[points_i];
-      Fuzzy_sphere fs(w, mu_epsilon, 0, points_traits);
-      std::vector<int> close_neighbors;
-      witness_tree.search(std::insert_iterator<std::vector<int>>(close_neighbors,close_neighbors.begin()),fs);
-      for (int i: close_neighbors)
-        dropped_points[i] = true;
-    }
-
-    for (unsigned points_i = 0; points_i < nbP; points_i++) {
-      if (dropped_points[points_i])
-        landmarks.push_back(points[points_i]);
-    }
-
-    if (nbL < landmarks.size()) {
-      std::random_shuffle(landmarks.begin(), landmarks.end());
-      landmarks.resize(nbL);
-    }
-  }
-    
-
-
-  
-  /** \brief Landmark choice strategy by taking random vertices for landmarks.
-   *  \details It chooses nbL distinct landmarks from a random access range `points`
-   *  and outputs a matrix {witness}*{closest landmarks} in knn.
-   */
-  template <typename KNearestNeighbours,
-            typename Point_random_access_range,
-            typename Distance_matrix>
-  void build_distance_matrix(Point_random_access_range const & points,
-                             Point_random_access_range & landmarks,
-                             FT alpha,
-                             unsigned limD,
-                             KNearestNeighbours & knn,
-                             Distance_matrix & distances)
-  {
-    int nbP = points.end() - points.begin();
-    knn = KNearestNeighbours(nbP);
-    distances = Distance_matrix(nbP);
-    STraits traits(&(landmarks[0]));
-    CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance> adapter(&(landmarks[0]));
-    Euclidean_distance ed;
-    Tree landmark_tree(boost::counting_iterator<std::ptrdiff_t>(0),
-                       boost::counting_iterator<std::ptrdiff_t>(landmarks.size()),
-                       typename Tree::Splitter(),
-                       traits);
-    for (int points_i = 0; points_i < nbP; points_i++) {
-      Point_d const & w = points[points_i];
-      Neighbor_search search(landmark_tree,
-                             w,
-                             FT(0),
-                             true,
-                             adapter);
-      Neighbor_search::iterator search_it = search.begin();
-      // Neighbor_search2 search(landmark_tree,
-      //                         w, limD+1,
-      //                         FT(0),
-      //                         true,
-      //                         adapter);
-      // Neighbor_search2::iterator search_it = search.begin();
-      
-      while (knn[points_i].size() <= limD) {
-        distances[points_i].push_back(search_it->second); //!sq_dist
-        knn[points_i].push_back((search_it++)->first);
-      }
-      FT dtow = distances[points_i][limD];
-      
-      while (search_it != search.end() && search_it->second < dtow + alpha) {
-        distances[points_i].push_back(search_it->second);
-        knn[points_i].push_back((search_it++)->first);
-      }
-      //std::cout << "k = " << knn[points_i].size() << std::endl;
-    }
-  }
- 
-  /*
-  template <typename Kernel, typename Point_container>
-  std::vector<typename Point_container::value_type>
-  sparsify_point_set(const Kernel &k,
-                     Point_container const& input_pts,
-                     typename Kernel::FT min_squared_dist)
-  {
-    typedef typename CGAL::Tangential_complex_::Point_cloud_data_structure<Kernel, Point_container> Points_ds;
-    typedef typename Points_ds::INS_iterator      INS_iterator;
-    typedef typename Points_ds::INS_range         INS_range;
-
-  typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object();
-
-  // Create the output container
-  std::vector<typename Point_container::value_type> output;
-
-  Points_ds points_ds(input_pts);
-
-  std::vector<bool> dropped_points(input_pts.size(), false);
-
-  // Parse the input points, and add them if they are not too close to
-  // the other points
-  std::size_t pt_idx = 0;
-  for (typename Point_container::const_iterator it_pt = input_pts.begin() ;
-       it_pt != input_pts.end();
-       ++it_pt, ++pt_idx)
-  {
-    if (dropped_points[pt_idx])
-      continue;
-
-    output.push_back(*it_pt);
-
-    INS_range ins_range = points_ds.query_incremental_ANN(*it_pt);
-
-    // If another point Q is closer that min_squared_dist, mark Q to be dropped
-    for (INS_iterator nn_it = ins_range.begin() ;
-        nn_it != ins_range.end() ;
-        ++nn_it)
-    {
-      std::size_t neighbor_point_idx = nn_it->first;
-      // If the neighbor is too close, we drop the neighbor
-      if (nn_it->second < min_squared_dist)
-      {
-        // N.B.: If neighbor_point_idx < pt_idx, 
-        // dropped_points[neighbor_point_idx] is already true but adding a
-        // test doesn't make things faster, so why bother?
-        dropped_points[neighbor_point_idx] = true;
-      }
-      else
-        break;
-    }
-  }
-
-  return output;
-}
-  */  
-
-
-}  // namespace witness_complex
-
-}  // namespace Gudhi
-
-#endif  // LANDMARK_CHOICE_BY_RANDOM_POINT_H_
diff --git a/src/Witness_complex/example/a0_persistence.cpp b/src/Witness_complex/example/a0_persistence.cpp
deleted file mode 100644
index 422185ca..00000000
--- a/src/Witness_complex/example/a0_persistence.cpp
+++ /dev/null
@@ -1,638 +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 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/Simplex_tree.h>
-#include <gudhi/A0_complex.h>
-#include <gudhi/Relaxed_witness_complex.h>
-#include <gudhi/Dim_lists.h>
-#include <gudhi/reader_utils.h>
-#include <gudhi/Persistent_cohomology.h>
-#include "Landmark_choice_random_knn.h"
-#include "Landmark_choice_sparsification.h"
-
-#include <iostream>
-#include <fstream>
-#include <ctime>
-#include <utility>
-#include <algorithm>
-#include <set>
-#include <queue>
-#include <iterator>
-#include <string>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-#include "generators.h"
-#include "output.h"
-#include "output_tikz.h"
-
-using namespace Gudhi;
-using namespace Gudhi::witness_complex;
-using namespace Gudhi::persistent_cohomology;
-
-typedef std::vector<Point_d> Point_Vector;
-typedef A0_complex< Simplex_tree<> > A0Complex;
-typedef Simplex_tree<>::Simplex_handle Simplex_handle;
-
-typedef A0_complex< Simplex_tree<> > SRWit;
-typedef Relaxed_witness_complex< Simplex_tree<> > WRWit;
-
-/**
- * \brief Customized version of read_points
- * which takes into account a possible nbP first line
- *
- */
-inline void
-read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) {
-  std::ifstream in_file(file_name.c_str(), std::ios::in);
-  if (!in_file.is_open()) {
-    std::cerr << "Unable to open file " << file_name << std::endl;
-    return;
-  }
-  std::string line;
-  double x;
-  while (getline(in_file, line)) {
-    std::vector< double > point;
-    std::istringstream iss(line);
-    while (iss >> x) {
-      point.push_back(x);
-    }
-    if (point.size() != 1)
-      points.push_back(point);
-  }
-  in_file.close();
-}
-
-void output_experiment_information(char * const file_name)
-{
-    std::cout << "Enter a valid experiment number. Usage: "
-              << file_name << " exp_no options\n";
-    std::cout << "Experiment description:\n"
-              << "0 nbP nbL dim alpha limD mu_epsilon: "
-              << "Build persistence diagram on relaxed witness complex "
-              << "built from a point cloud on (dim-1)-dimensional sphere "
-              << "consisting of nbP witnesses and nbL landmarks. "
-              << "The maximal relaxation is alpha and the limit on simplicial complex "
-              << "dimension is limD.\n";
-    std::cout << "1 file_name nbL alpha limD: "
-              << "Build persistence diagram on relaxed witness complex "
-              << "build from a point cloud stored in a file and nbL landmarks. "
-              << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n";
-}
-
-void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT mu_epsilon = 0.1)
-{
-  clock_t start, end;
-  Simplex_tree<> simplex_tree;
-
-  // Choose landmarks
-  std::vector<std::vector< int > > knn;
-  std::vector<std::vector< FT > > distances;
-  start = clock();
-  //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances);
-
-  std::vector<Point_d> landmarks;
-  Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                landmarks,  // aka landmarks
-                                                alpha2,
-                                                limD,
-                                                knn,
-                                                distances);
-  end = clock();
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Choice of " << nbL << " landmarks took "
-            << time << " s. \n";
-  // Compute witness complex
-  start = clock();
-  A0Complex rw(distances,
-                knn,
-                simplex_tree,
-                nbL,
-                alpha2,
-                limD);
-  end = clock();
-  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Witness complex for " << nbL << " landmarks took "
-            << time << " s. \n";
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  int streedim = 0;
-  for (auto s: simplex_tree.complex_simplex_range())
-    if (simplex_tree.dimension(s) > streedim)
-      streedim = simplex_tree.dimension(s);
-  std::cout << "Dimension of the simplicial complex is " << streedim << std::endl;
-
-  //std::cout << simplex_tree << "\n";
-  
-  // Compute the persistence diagram of the complex
-  simplex_tree.set_dimension(limD);
-  persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, true);
-  int p = 3;
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  start = clock();
-  pcoh.compute_persistent_cohomology( alpha2/10 );
-  end = clock();
-  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Persistence diagram took "
-            << time << " s. \n";
-  pcoh.output_diagram();
-
-  int chi = 0;
-  for (auto sh: simplex_tree.complex_simplex_range())
-    chi += 1-2*(simplex_tree.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-
-  // Gudhi::witness_complex::Dim_lists<Simplex_tree<>> simplices(simplex_tree, limD); 
-  
-  // simplices.collapse();
-  // simplices.output_simplices();
-
-  // Simplex_tree<> collapsed_tree;
-  // for (auto sh: simplices) {
-  //   std::vector<int> vertices;
-  //   for (int v: collapsed_tree.simplex_vertex_range(sh))
-  //     vertices.push_back(v);
-  //   collapsed_tree.insert_simplex(vertices);
-  // } 
-  std::vector<int> landmarks_ind(nbL); 
-  for (unsigned i = 0; i != distances.size(); ++i) {
-    if (distances[i][0] == 0)
-      landmarks_ind[knn[i][0]] = i;
-  }
-  //write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplices, false, true);
-  write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, "witness_before_collapse.mesh");
-
-  // collapsed_tree.set_dimension(limD);
-  // persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh2(collapsed_tree, true);
-  // pcoh2.init_coefficients( p ); //initilizes the coefficient field for homology
-  // pcoh2.compute_persistent_cohomology( alpha2/10 );
-  // pcoh2.output_diagram();
-
-  // chi = 0;
-  // for (auto sh: simplices)
-  //   chi += 1-2*(simplex_tree.dimension(sh)%2);
-  // std::cout << "Euler characteristic is " << chi << std::endl;
-  // write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, "witness_after_collapse.mesh");
-}
-
-void rips_experiment(Point_Vector & points, double threshold, int dim_max)
-{
-  typedef std::vector<double> Point_t;
-  typedef Simplex_tree<Simplex_tree_options_fast_persistence> ST;
-  clock_t start, end;
-  ST st;
-
-  // Compute the proximity graph of the points
-  start = clock();
-  Graph_t prox_graph = compute_proximity_graph(points, threshold
-                                               , euclidean_distance<Point_t>);
-  // Construct the Rips complex in a Simplex Tree
-  // insert the proximity graph in the simplex tree
-  st.insert_graph(prox_graph);
-  // expand the graph until dimension dim_max
-  st.expansion(dim_max);
-  end = clock();
-  
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Rips complex took "
-            << time << " s. \n";
-  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();
-
-  // Compute the persistence diagram of the complex
-  persistent_cohomology::Persistent_cohomology<ST, Field_Zp > pcoh(st);
-  // initializes the coefficient field for homology
-  int p = 3;
-  double min_persistence = -1; //threshold/5;
-  pcoh.init_coefficients(p);
-  pcoh.compute_persistent_cohomology(min_persistence);
-  pcoh.output_diagram();
-}
-
-
-int experiment0 (int argc, char * const argv[])
-{
-  if (argc != 8) {
-    std::cerr << "Usage: " << argv[0]
-              << " 0 nbP nbL dim alpha limD mu_epsilon\n";
-    return 0;
-  }
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  int nbP       = atoi(argv[2]);
-  int nbL       = atoi(argv[3]);
-  int dim       = atoi(argv[4]);
-  double alpha  = atof(argv[5]);
-  int limD      = atoi(argv[6]);
-  double mu_epsilon = atof(argv[7]);
-
-  // Read the point file
-  Point_Vector point_vector;
-  generate_points_sphere(point_vector, nbP, dim);
-  std::cout << "Successfully generated " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-  rw_experiment(point_vector, nbL, alpha, limD);
-  return 0;
-}
-
-// int experiment1 (int argc, char * const argv[])
-// {
-//   if (argc != 3) {
-//     std::cerr << "Usage: " << argv[0]
-//               << " 1 file_name\n";
-//     return 0;
-//   }
-//   /*
-//     boost::filesystem::path p;
-//     for (; argc > 2; --argc, ++argv)
-//     p /= argv[1];
-//   */
-  
-//   std::string file_name = argv[2];
-
-//   // Read the point file
-//   Point_Vector point_vector;
-//   read_points_cust(file_name, point_vector);
-//   std::cout << "The file contains " << point_vector.size() << " points.\n";
-//   std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-//   bool ok = false;
-//   int nbL, limD;
-//   double alpha;
-//   while (!ok) {
-//     std::cout << "Relaxed witness complex: parameters nbL, alpha, limD.\n";
-//     std::cout << "Enter nbL: ";
-//     std::cin >> nbL;
-//     std::cout << "Enter alpha: ";
-//     std::cin >> alpha;
-//     std::cout << "Enter limD: ";
-//     std::cin >> limD;
-//     std::cout << "Start relaxed witness complex...\n";
-//     rw_experiment(point_vector, nbL, alpha, limD);
-//     std::cout << "Is the result correct? [y/n]: ";
-//     char answer;
-//     std::cin >> answer;
-//     switch (answer) {
-//     case 'n':
-//       ok = false; break;
-//     default :
-//       ok = true; break;
-//     }
-//   }
-//   // ok = false;
-//   // while (!ok) {
-//   //   std::cout << "Rips complex: parameters threshold, limD.\n";
-//   //   std::cout << "Enter threshold: ";
-//   //   std::cin >> alpha;
-//   //   std::cout << "Enter limD: ";
-//   //   std::cin >> limD;
-//   //   std::cout << "Start Rips complex...\n";
-//   //   rips_experiment(point_vector, alpha, limD);
-//   //   std::cout << "Is the result correct? [y/n]: ";
-//   //   char answer;
-//   //   std::cin >> answer;
-//   //   switch (answer) {
-//   //   case 'n':
-//   //     ok = false; break;
-//   //   default :
-//   //     ok = true; break;
-//   //   }
-//   // }
-//   return 0;
-// }
-
-int experiment1 (int argc, char * const argv[])
-{
-  if (argc != 8) {
-    std::cerr << "Usage: " << argv[0]
-              << " 1 file_name nbL alpha mu_epsilon limD experiment_name\n";
-    return 0;
-  }
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  std::string file_name = argv[2];
-  int nbL = atoi(argv[3]), limD = atoi(argv[6]);
-  double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]);
-  std::string experiment_name = argv[7];
-
-  // Read the point file
-  Point_Vector point_vector;
-  read_points_cust(file_name, point_vector);
-  std::cout << "The file contains " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-  
-  Simplex_tree<> simplex_tree;
-  std::vector<std::vector< int > > knn;
-  std::vector<std::vector< FT > > distances;
-  std::vector<Point_d> landmarks;
-  Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                landmarks,  // aka landmarks
-                                                alpha2,
-                                                limD,
-                                                knn,
-                                                distances);
-
-  rw_experiment(point_vector, nbL, alpha2, limD, mu_epsilon);
-  
-  // ok = false;
-  // while (!ok) {
-  //   std::cout << "Rips complex: parameters threshold, limD.\n";
-  //   std::cout << "Enter threshold: ";
-  //   std::cin >> alpha;
-  //   std::cout << "Enter limD: ";
-  //   std::cin >> limD;
-  //   std::cout << "Start Rips complex...\n";
-  //   rips_experiment(point_vector, alpha, limD);
-  //   std::cout << "Is the result correct? [y/n]: ";
-  //   char answer;
-  //   std::cin >> answer;
-  //   switch (answer) {
-  //   case 'n':
-  //     ok = false; break;
-  //   default :
-  //     ok = true; break;
-  //   }
-  // }
-  return 0;
-}
-
-
-/********************************************************************************************
- * Length of the good interval experiment
- *******************************************************************************************/
-
-struct Pers_endpoint {
-  double alpha;
-  bool start;
-  int dim;
-  Pers_endpoint(double alpha_, bool start_, int dim_)
-    : alpha(alpha_), start(start_), dim(dim_)
-  {}
-};
-
-/*
-struct less_than_key {
-  inline bool operator() (const MyStruct& struct1, const MyStruct& struct2) {
-    return (struct1.key < struct2.key);
-  }
-};
-*/
-
-double good_interval_length(const std::vector<int> & desired_homology, Simplex_tree<> & simplex_tree, double alpha2)
-{
-  int nbL = simplex_tree.num_vertices();
-  int p = 3;
-  persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, true);
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  pcoh.compute_persistent_cohomology( -1 );
-  std::ofstream out_stream("pers_diag.tmp");
-  pcoh.output_diagram(out_stream);
-  out_stream.close();
-  std::ifstream in_stream("pers_diag.tmp", std::ios::in);
-  std::string line;
-  std::vector<Pers_endpoint> pers_endpoints;
-  while (getline(in_stream, line)) {
-    int p, dim;
-    double alpha_start, alpha_end;
-    std::istringstream iss(line);
-    iss >> p >> dim >> alpha_start >> alpha_end;
-    if (alpha_start != alpha_end) {
-      if (alpha_end < alpha_start)
-        alpha_end = alpha2;
-      pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim));
-      pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim));
-    }
-  }
-  std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl;
-  in_stream.close();
-  std::sort(pers_endpoints.begin(),
-            pers_endpoints.end(),
-            [](const Pers_endpoint & p1, const Pers_endpoint & p2){
-              return p1.alpha < p2.alpha;}
-            );
-  write_barcodes("pers_diag.tmp", alpha2);
-  /*
-  for (auto p: pers_endpoints) {
-    std::cout << p.alpha << " " << p.dim << " " << p.start << "\n";
-  }
-  */
-  std::vector<int> current_homology(nbL-1,0);
-  current_homology[0] = 1; // for the compulsary "0 0 inf" entry
-  double good_start = 0, good_end = 0;
-  double sum_intervals = 0;
-  int num_pieces = 0;
-  bool interval_in_process = (desired_homology == current_homology);
-  for (auto p: pers_endpoints) {
-    /*
-    std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start
-              << " [";
-    for (int v: current_homology)
-      std::cout << v << " ";
-    std::cout << "]\n";
-    */
-    if (p.start)
-      current_homology[p.dim]++;
-    else
-      current_homology[p.dim]--;
-    if (interval_in_process) {
-      good_end = p.alpha;
-      sum_intervals += good_end - good_start;
-      std::cout << "good_start = " << good_start
-                << ", good_end = " << good_end << "\n";
-                
-      Gudhi::witness_complex::Dim_lists<Simplex_tree<>> simplices(simplex_tree, nbL-1, (good_end - good_start)/2);
-      simplices.collapse();
-      simplices.output_simplices();
-      interval_in_process = false;
-      //break;
-    }
-    else if (desired_homology == current_homology) {
-      interval_in_process = true;
-      good_start = p.alpha;
-      num_pieces++;
-    }
-  }
-  std::cout << "Number of good homology intervals: " << num_pieces << "\n";
-  return sum_intervals;
-}
-
-void run_comparison(std::vector<std::vector< int > > const & knn,
-                    std::vector<std::vector< FT > > const & distances,
-                    unsigned nbL,
-                    double alpha2,
-                    std::vector<int>& desired_homology)
-{
-  clock_t start, end;
-  Simplex_tree<> simplex_tree;
-
-  start = clock();
-  SRWit srwit(distances,
-              knn,
-              simplex_tree,
-              nbL,
-              alpha2,
-              nbL-1);
-  end = clock();
-  std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl;
-  simplex_tree.set_dimension(nbL-1);
-
-  std::cout << "Good homology interval length for SRWit is "
-                << good_interval_length(desired_homology, simplex_tree, alpha2) << "\n";
-  std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
-
-  /*
-  Simplex_tree<> simplex_tree2;
-  start = clock();
-  WRWit wrwit(distances,
-              knn,
-              simplex_tree2,
-              nbL,
-              alpha2,
-              nbL-1);
-  end = clock();
-  std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl;
-  simplex_tree.set_dimension(nbL-1);
-  
-  std::cout << "Good homology interval length for WRWit is "
-            << good_interval_length(desired_homology, simplex_tree2, alpha2) << "\n";
-  std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
-  */
-}
-
-int experiment2(int argc, char * const argv[])
-{
-  for (unsigned d = 2; d < 2; d++) {
-    // Sphere S^d 
-    Point_Vector point_vector;
-    unsigned N = 1;
-    double alpha2 = 2.4 - 0.4*d; 
-    switch (d) {
-    case 1: alpha2 = 2.2; break;
-    case 2: alpha2 = 1.8; break;
-    case 3: alpha2 = 1.5; break;
-    case 4: alpha2 = 1.4; break;
-    default: alpha2 = 1.4; break;
-    }
-    std::cout << "alpha2 = " << alpha2 << "\n";
-    unsigned nbL = 20;
-    std::vector<int> desired_homology(nbL-1,0);
-    desired_homology[0] = 1; desired_homology[d] = 1;
-    
-    
-    for (unsigned i = 1; i <= N; ++i) {
-      unsigned nbW = 1000*i;//, nbL = 20;
-      double mu_epsilon = 1/sqrt(nbL);
-      std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl;
-      generate_points_sphere(point_vector, i*1000, d+1);
-      std::vector<Point_d> landmarks;
-      Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-
-      std::vector<std::vector< int > > knn;
-      std::vector<std::vector< FT > > distances;
-    
-      std::cout << "|L| after sparsification: " << landmarks.size() << "\n";
-      Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                    landmarks,  // aka landmarks
-                                                    alpha2,
-                                                    nbL-1,
-                                                    knn,
-                                                    distances);
-      run_comparison(knn, distances, nbL, alpha2, desired_homology);
-    }
-  }
-  {
-    // SO(3)
-    Point_Vector point_vector;
-    double alpha2 = 0.6; 
-    std::cout << "alpha2 = " << alpha2 << "\n";
-    unsigned nbL = 150;
-    std::vector<int> desired_homology(nbL-1,0);
-    desired_homology[0] = 1; desired_homology[1] = 1; desired_homology[2] = 1; //Kl
-    // desired_homology[0] = 1; desired_homology[3] = 1; //SO3
-
-    double mu_epsilon = 1/sqrt(nbL);
-    if (argc < 3) std::cerr << "No file name indicated!\n";
-    read_points_cust(argv[2], point_vector);
-    int nbW = point_vector.size();
-    std::cout << "Running test SO(3), |W|=" << nbW << ", |L|=" << nbL << std::endl;
-    std::vector<Point_d> landmarks;
-    Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-    
-    std::vector<std::vector< int > > knn;
-    std::vector<std::vector< FT > > distances;
-    
-    std::cout << "|L| after sparsification: " << landmarks.size() << "\n";
-    Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                  landmarks,  // aka landmarks
-                                                  alpha2,
-                                                  nbL-1,
-                                                  knn,
-                                                  distances);
-    run_comparison(knn, distances, nbL, alpha2, desired_homology);
-  }
-  return 0;
-}
-
-int experiment3(int argc, char * const argv[])
-{
-  // COLLAPSES EXPERIMENT
-  
-  return 0;
-}
-
-int main (int argc, char * const argv[])
-{
-  if (argc == 1) {
-    output_experiment_information(argv[0]);
-    return 1;
-  }
-  switch (atoi(argv[1])) {
-  case 0 :
-    return experiment0(argc, argv);
-    break;
-  case 1 :
-    return experiment1(argc, argv);
-    break;
-  case 2 :
-    return experiment2(argc, argv);
-    break;
-  default :
-    output_experiment_information(argv[0]);
-    return 1;
-  }
-}
diff --git a/src/Witness_complex/example/bench_rwit.cpp b/src/Witness_complex/example/bench_rwit.cpp
deleted file mode 100644
index 2d3a009c..00000000
--- a/src/Witness_complex/example/bench_rwit.cpp
+++ /dev/null
@@ -1,709 +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 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/Simplex_tree.h>
-#include <gudhi/A0_complex.h>
-#include <gudhi/Relaxed_witness_complex.h>
-#include <gudhi/Dim_lists.h>
-#include <gudhi/reader_utils.h>
-#include <gudhi/Persistent_cohomology.h>
-#include <gudhi/Good_links.h>
-#include "Landmark_choice_random_knn.h"
-#include "Landmark_choice_sparsification.h"
-
-#include <iostream>
-#include <fstream>
-#include <ctime>
-#include <utility>
-#include <algorithm>
-#include <set>
-#include <queue>
-#include <iterator>
-#include <string>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-#include <boost/program_options.hpp>
-
-#include "generators.h"
-#include "output.h"
-#include "output_tikz.h"
-
-using namespace Gudhi;
-using namespace Gudhi::witness_complex;
-using namespace Gudhi::persistent_cohomology;
-
-typedef std::vector<Point_d> Point_Vector;
-//typedef Simplex_tree<Simplex_tree_options_fast_persistence> STree;
-typedef Simplex_tree<> STree;
-typedef A0_complex< STree> A0Complex;
-typedef STree::Simplex_handle Simplex_handle;
-
-typedef A0_complex< STree > SRWit;
-typedef Relaxed_witness_complex< STree > WRWit;
-
-/** Program options ***************************************************************
-***********************************************************************************
-***********************************************************************************
-***********************************************************************************
-**********************************************************************************/
-
-
-void program_options(int argc, char * const argv[]
-                     , int & experiment_number
-                     , std::string & filepoints
-                     , std::string & landmark_file
-                     , std::string & experiment_name
-                     , int & nbL
-                     , double & alpha2_s
-                     , double & alpha2_w
-                     , double & mu_epsilon
-                     , int & dim_max
-                     , std::vector<int> & desired_homology
-                     , double & min_persistence) {
-  namespace po = boost::program_options;
-  po::options_description hidden("Hidden options");
-  hidden.add_options()
-    ("option", po::value<int>(& experiment_number),
-     "Experiment id.")
-    ("input-file", po::value<std::string>(&filepoints),
-       "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>(&experiment_name)->default_value("witness"),
-       "The prefix of all the output files. Default is 'witness'")
-      ("landmarks,L", po::value<int>(&nbL)->default_value(0),
-       "Number of landmarks.")
-      ( "landmark-file,l", po::value<std::string>(&landmark_file),
-        "Name of a fike containing landmarks")
-      ("alpha2_s,A", po::value<double>(&alpha2_s)->default_value(0),
-       "Relaxation parameter for the strong complex.")
-      ("alpha2_w,a", po::value<double>(&alpha2_w)->default_value(0),
-       "Relaxation parameter for the weak complex.")
-      ("mu_epsilon,e", po::value<double>(&mu_epsilon)->default_value(0),
-       "Sparsification parameter.")
-      ("cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
-       "Maximal dimension of the Witness complex we want to compute.")
-    
-    ("homology,H", po::value<std::vector<int>>(&desired_homology)->multitoken(),
-    "The desired Betti numbers.")
-      ("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("option", 1);
-  pos.add("input-file", 2);
-
-  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/3Z \n";
-    std::cout << "of a Strong relaxed 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();
-  }
-}
-
-
-
-
-
-/**
- * \brief Customized version of read_points
- * which takes into account a possible nbP first line
- *
- */
-inline void
-read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) {
-  std::ifstream in_file(file_name.c_str(), std::ios::in);
-  if (!in_file.is_open()) {
-    std::cerr << "Unable to open file " << file_name << std::endl;
-    return;
-  }
-  std::string line;
-  double x;
-  while (getline(in_file, line)) {
-    std::vector< double > point;
-    std::istringstream iss(line);
-    while (iss >> x) {
-      point.push_back(x);
-    }
-    if (point.size() != 1)
-      points.push_back(point);
-  }
-  in_file.close();
-}
-
-void rips(Point_Vector & points, double alpha2, int dim_max, STree& st)
-{
-  Graph_t prox_graph = compute_proximity_graph(points, sqrt(alpha2)
-                                               , euclidean_distance<std::vector<FT> >);
-  // Construct the Rips complex in a Simplex Tree
-  // insert the proximity graph in the simplex tree
-  st.insert_graph(prox_graph);
-  // expand the graph until dimension dim_max
-  st.expansion(dim_max);
-}
-
-void output_experiment_information(char * const file_name)
-{
-    std::cout << "Enter a valid experiment number. Usage: "
-              << file_name << " exp_no options\n";
-    std::cout << "Experiment description:\n"
-              << "0 nbP nbL dim alpha limD mu_epsilon: "
-              << "Build persistence diagram on relaxed witness complex "
-              << "built from a point cloud on (dim-1)-dimensional sphere "
-              << "consisting of nbP witnesses and nbL landmarks. "
-              << "The maximal relaxation is alpha and the limit on simplicial complex "
-              << "dimension is limD.\n";
-    std::cout << "1 file_name nbL alpha limD: "
-              << "Build persistence diagram on relaxed witness complex "
-              << "build from a point cloud stored in a file and nbL landmarks. "
-              << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n";
-}
-
-void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha2, int limD, FT mu_epsilon = 0.1,
-                   std::string mesh_filename = "witness")
-{
-  clock_t start, end;
-  STree simplex_tree;
-
-  // Choose landmarks
-  std::vector<std::vector< int > > knn;
-  std::vector<std::vector< FT > > distances;
-  start = clock();
-  //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances);
-
-  std::vector<Point_d> landmarks;
-  Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                landmarks,  // aka landmarks
-                                                alpha2,
-                                                limD,
-                                                knn,
-                                                distances);
-  end = clock();
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Choice of " << nbL << " landmarks took "
-            << time << " s. \n";
-  // Compute witness complex
-  start = clock();
-  A0Complex rw(distances,
-                knn,
-                simplex_tree,
-                nbL,
-                alpha2,
-                limD);
-  end = clock();
-  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Witness complex for " << nbL << " landmarks took "
-            << time << " s. \n";
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  //std::cout << simplex_tree << "\n";
-  
-  // Compute the persistence diagram of the complex
-  simplex_tree.set_dimension(limD);
-  persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh(simplex_tree, true);
-  int p = 3;
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  start = clock();
-  pcoh.compute_persistent_cohomology( alpha2/10 );
-  end = clock();
-  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Persistence diagram took "
-            << time << " s. \n";
-  pcoh.output_diagram();
-
-  int chi = 0;
-  for (auto sh: simplex_tree.complex_simplex_range())
-    chi += 1-2*(simplex_tree.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-
-  Gudhi::witness_complex::Dim_lists<STree> simplices(simplex_tree, limD); 
-  
-  // std::vector<Simplex_handle> simplices;
-  std::cout << "Starting collapses...\n";
-  simplices.collapse();
-  simplices.output_simplices();
-
-  STree collapsed_tree;
-  for (auto sh: simplices) {
-    std::vector<int> vertices;
-    for (int v: collapsed_tree.simplex_vertex_range(sh))
-      vertices.push_back(v);
-    collapsed_tree.insert_simplex(vertices);
-  } 
-  std::vector<int> landmarks_ind(nbL); 
-  for (unsigned i = 0; i != distances.size(); ++i) {
-    if (distances[i][0] == 0)
-      landmarks_ind[knn[i][0]] = i;
-  }
-  //write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplices, false, true);
-  write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), false, true, mesh_filename+"_before_collapse.mesh");
-
-  collapsed_tree.set_dimension(limD);
-  persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh2(collapsed_tree, true);
-  pcoh2.init_coefficients( p ); //initilizes the coefficient field for homology
-  pcoh2.compute_persistent_cohomology( alpha2/10 );
-  pcoh2.output_diagram();
-
-  chi = 0;
-  for (auto sh: simplices)
-    chi += 1-2*(simplex_tree.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-  write_witness_mesh(point_vector, landmarks_ind, collapsed_tree, collapsed_tree.complex_simplex_range(), false, true, mesh_filename+"_after_collapse.mesh");
-  Gudhi::Good_links<STree> gl(collapsed_tree);
-  if (gl.complex_is_pseudomanifold())
-    std::cout << "Collapsed complex is a pseudomanifold.\n";
-  else
-    std::cout << "Collapsed complex is NOT a pseudomanifold.\n";
-  bool good = true;
-  for (auto v: collapsed_tree.complex_vertex_range())
-    if (!gl.has_good_link(v)) {
-      std::cout << "Bad link around " << v << std::endl;
-      good = false;
-    }
-  if (good)
-    std::cout << "All links are good.\n";
-  else
-    std::cout << "There are bad links.\n";
-}
-
-void rips_experiment(Point_Vector & points, double threshold, int dim_max)
-{
-  typedef STree ST;
-  clock_t start, end;
-  ST st;
-
-  // Compute the proximity graph of the points
-  start = clock();
-  rips(points, threshold, dim_max, st);
-  end = clock();
-  
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Rips complex took "
-            << time << " s. \n";
-  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();
-
-  // Compute the persistence diagram of the complex
-  persistent_cohomology::Persistent_cohomology<ST, Field_Zp > pcoh(st);
-  // initializes the coefficient field for homology
-  int p = 3;
-  double min_persistence = -1; //threshold/5;
-  pcoh.init_coefficients(p);
-  pcoh.compute_persistent_cohomology(min_persistence);
-  pcoh.output_diagram();
-}
-
-
-int experiment0 (int argc, char * const argv[])
-{
-  if (argc != 8) {
-    std::cerr << "Usage: " << argv[0]
-              << " 0 nbP nbL dim alpha limD mu_epsilon\n";
-    return 0;
-  }
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  int nbP       = atoi(argv[2]);
-  int nbL       = atoi(argv[3]);
-  int dim       = atoi(argv[4]);
-  double alpha  = atof(argv[5]);
-  int limD      = atoi(argv[6]);
-  double mu_epsilon = atof(argv[7]);
-
-  // Read the point file
-  Point_Vector point_vector;
-  generate_points_sphere(point_vector, nbP, dim);
-  std::cout << "Successfully generated " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-  rw_experiment(point_vector, nbL, alpha, limD);
-  return 0;
-}
-
-
-/********************************************************************************************
- * Length of the good interval experiment
- *******************************************************************************************/
-
-struct Pers_endpoint {
-  double alpha;
-  bool start;
-  int dim;
-  Pers_endpoint(double alpha_, bool start_, int dim_)
-    : alpha(alpha_), start(start_), dim(dim_)
-  {}
-};
-
-/*
-struct less_than_key {
-  inline bool operator() (const MyStruct& struct1, const MyStruct& struct2) {
-    return (struct1.key < struct2.key);
-  }
-};
-*/
-
-double good_interval_length(const std::vector<int> & desired_homology, STree & simplex_tree, double alpha2)
-{
-  int nbL = simplex_tree.num_vertices();
-  int p = 3;
-  persistent_cohomology::Persistent_cohomology< STree, Field_Zp > pcoh(simplex_tree, true);
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  pcoh.compute_persistent_cohomology( -1 );
-  std::ofstream out_stream("pers_diag.tmp");
-  pcoh.output_diagram(out_stream);
-  out_stream.close();
-  std::ifstream in_stream("pers_diag.tmp", std::ios::in);
-  std::string line;
-  std::vector<Pers_endpoint> pers_endpoints;
-  while (getline(in_stream, line)) {
-    unsigned p, dim;
-    double alpha_start, alpha_end;
-    std::istringstream iss(line);
-    iss >> p >> dim >> alpha_start >> alpha_end;
-    if (iss.fail())
-      alpha_end = alpha2;
-    //std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n";
-    //if (dim < desired_homology.size()+1)
-    if (alpha_start != alpha_end) {
-    //   if (alpha_end < alpha_start)
-    //     alpha_end = alpha2;
-      pers_endpoints.push_back(Pers_endpoint(alpha_start, true, dim));
-      pers_endpoints.push_back(Pers_endpoint(alpha_end, false, dim));
-      std::cout << p << " " << dim << " " << alpha_start << " " << alpha_end << "\n";
-    }
-  }
-  std::cout << "desired_homology.size() = " << desired_homology.size() << "\n";
-  for (auto nd: desired_homology)
-    std::cout << nd << std::endl;
-  std::cout << "Pers_endpoints.size = " << pers_endpoints.size() << std::endl;
-  in_stream.close();
-  std::sort(pers_endpoints.begin(),
-            pers_endpoints.end(),
-            [](const Pers_endpoint & p1, const Pers_endpoint & p2){
-              return p1.alpha < p2.alpha;}
-            );
-  write_barcodes("pers_diag.tmp", alpha2);
-  /*
-  for (auto p: pers_endpoints) {
-    std::cout << p.alpha << " " << p.dim << " " << p.start << "\n";
-  }
-  */
-  std::vector<int> current_homology(desired_homology.size(),0);
-  //current_homology[0] = 1; // for the compulsary "0 0 inf" entry
-  double good_start = 0, good_end = 0;
-  double sum_intervals = 0;
-  int num_pieces = 0;
-  bool interval_in_process = (desired_homology == current_homology);
-  for (auto p: pers_endpoints) {
-    /*
-      std::cout << p.alpha << " " << p.dim << " ";
-    if (p.start)
-      std::cout << "s\n";
-    else
-      std::cout << "e\n";
-    */
-    /*
-    std::cout << "Treating " << p.alpha << " " << p.dim << " " << p.start
-              << " [";
-    for (int v: current_homology)
-      std::cout << v << " ";
-    std::cout << "]\n";
-    */
-    if (p.start)
-      current_homology[p.dim]++;
-    else
-      current_homology[p.dim]--;
-    if (interval_in_process) {
-      good_end = p.alpha;
-      sum_intervals += good_end - good_start;
-      std::cout << "good_start = " << good_start
-                << ", good_end = " << good_end << "\n";
-                
-      Gudhi::witness_complex::Dim_lists<STree> simplices(simplex_tree, nbL-1, (good_end - good_start)/2);
-      //simplices.collapse();
-      //simplices.output_simplices();
-      interval_in_process = false;
-      //break;
-    }
-    else if (desired_homology == current_homology) {
-      interval_in_process = true;
-      good_start = p.alpha;
-      num_pieces++;
-    }
-  }
-  std::cout << "Number of good homology intervals: " << num_pieces << "\n";
-  return sum_intervals;
-}
-
-
-
-void run_comparison(std::vector<std::vector< int > > const & knn,
-                    std::vector<std::vector< FT > > const & distances,
-                    Point_Vector & points,
-                    unsigned nbL,
-                    unsigned limD,
-                    double alpha2_s,
-                    double alpha2_w,
-                    std::vector<int>& desired_homology)
-{
-  clock_t start, end;
-  STree simplex_tree;
-
-  //std::cout << "alpha2 = " << alpha2_s << "\n";
-  start = clock();
-  SRWit srwit(distances,
-              knn,
-              simplex_tree,
-              nbL,
-              alpha2_s,
-              limD);
-  end = clock();
-  std::cout << "SRWit.size = " << simplex_tree.num_simplices() << std::endl;
-  simplex_tree.set_dimension(desired_homology.size());
-
-  std::cout << "Good homology interval length for SRWit is "
-                << good_interval_length(desired_homology, simplex_tree, alpha2_s) << "\n";
-  std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
-  int chi = 0;
-  for (auto sh: simplex_tree.complex_simplex_range())
-    chi += 1-2*(simplex_tree.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-
-
-  STree simplex_tree2;
-  std::cout << "alpha2 = " << alpha2_w << "\n";
-  start = clock();
-  WRWit wrwit(distances,
-              knn,
-              simplex_tree2,
-              nbL,
-              alpha2_w,
-              limD);
-  end = clock();
-  std::cout << "WRWit.size = " << simplex_tree2.num_simplices() << std::endl;
-  simplex_tree2.set_dimension(nbL-1);
-  
-  std::cout << "Good homology interval length for WRWit is "
-            << good_interval_length(desired_homology, simplex_tree2, alpha2_w) << "\n";
-  std::cout << "Time: " << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
-  chi = 0;
-  for (auto sh: simplex_tree2.complex_simplex_range())
-    chi += 1-2*(simplex_tree2.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-
-  //write_witness_mesh(points, landmarks_ind, simplex_tree2, simplex_tree2.complex_simplex_range(), false, true, "wrwit.mesh");
-
-  
-}
-
-int experiment1 (int argc, char * const argv[])
-{
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  // std::string file_name = argv[2];
-  // int nbL = atoi(argv[3]), limD = atoi(argv[6]);
-  // double alpha2 = atof(argv[4]), mu_epsilon = atof(argv[5]);
-  // std::string experiment_name = argv[7];
-
-  int option = 1;
-  std::string file_name, landmark_file;
-  int nbL = 0, limD;
-  double alpha2_s, alpha2_w, mu_epsilon, min_pers;
-  std::string experiment_name;
-  std::vector<int> desired_homology = {1};
-  std::vector<Point_d> landmarks;
-  
-  program_options(argc, argv, option, file_name, landmark_file, experiment_name, nbL, alpha2_s, alpha2_w, mu_epsilon, limD, desired_homology, min_pers);
-
-  // Read the point file
-  Point_Vector point_vector;
-  read_points_cust(file_name, point_vector);
-  //std::cout << "The file contains " << point_vector.size() << " points.\n";
-  //std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-  //std::cout << "Limit dimension for the complexes is " << limD << ".\n";
-
-  if (landmark_file == "")
-    Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-    //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha2_s, limD, knn, distances);
-  else
-    read_points_cust(landmark_file, landmarks);
-  nbL = landmarks.size();
-  STree simplex_tree;
-  std::vector<std::vector< int > > knn;
-  std::vector<std::vector< FT > > distances;
-
-  //Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                landmarks,  // aka landmarks
-                                                alpha2_s,
-                                                limD,
-                                                knn,
-                                                distances);
-
-  run_comparison(knn, distances, point_vector, nbL, limD, alpha2_s, alpha2_w, desired_homology);
-  return 0;
-}
-
-
-int experiment2(int argc, char * const argv[])
-{
-  for (unsigned d = 3; d < 4; d++) {
-    // Sphere S^d 
-    Point_Vector point_vector;
-    unsigned N = 1;
-    double alpha2 = 2.4 - 0.4*d; 
-    switch (d) {
-    case 1: alpha2 = 2.2; break;
-    case 2: alpha2 = 1.7; break;
-    case 3: alpha2 = 1.5; break;
-    case 4: alpha2 = 1.4; break;
-    default: alpha2 = 1.4; break;
-    }
-    unsigned nbL = 20;
-    std::vector<int> desired_homology(nbL-1,0);
-    desired_homology[0] = 1; desired_homology[d] = 1;
-    
-    
-    for (unsigned i = 1; i <= N; ++i) {
-      unsigned nbW = 1000*i;//, nbL = 20;
-      double mu_epsilon = 1/sqrt(nbL);
-      std::cout << "Running test S"<< d <<", |W|=" << nbW << ", |L|=" << nbL << std::endl;
-      generate_points_sphere(point_vector, i*1000, d+1);
-      std::vector<Point_d> landmarks;
-
-      Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-      
-      std::vector<std::vector< int > > knn;
-      std::vector<std::vector< FT > > distances;
-
-      
-      std::cout << "|L| after sparsification: " << landmarks.size() << "\n";
-      
-      Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                    landmarks,  // aka landmarks
-                                                    alpha2,
-                                                    nbL-1,
-                                                    knn,
-                                                    distances);
-      run_comparison(knn, distances, point_vector, nbL, nbL-1, alpha2, alpha2, desired_homology);
-    }
-  }
-  /*
-  {
-    // SO(3)
-    Point_Vector point_vector;
-    double alpha2 = 0.6; 
-    std::cout << "alpha2 = " << alpha2 << "\n";
-    unsigned nbL = 150;
-    std::vector<int> desired_homology(nbL-1,0);
-    desired_homology[0] = 1; desired_homology[1] = 1; desired_homology[2] = 1; //Kl
-    // desired_homology[0] = 1; desired_homology[3] = 1; //SO3
-
-    double mu_epsilon = 1/sqrt(nbL);
-    if (argc < 3) std::cerr << "No file name indicated!\n";
-    read_points_cust(argv[2], point_vector);
-    int nbW = point_vector.size();
-    std::cout << "Running test SO(3), |W|=" << nbW << ", |L|=" << nbL << std::endl;
-    std::vector<Point_d> landmarks;
-    Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-    
-    std::vector<std::vector< int > > knn;
-    std::vector<std::vector< FT > > distances;
-    
-    std::cout << "|L| after sparsification: " << landmarks.size() << "\n";
-    Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                  landmarks,  // aka landmarks
-                                                  alpha2,
-                                                  nbL-1,
-                                                  knn,
-                                                  distances);
-    run_comparison(knn, distances, point_vector, nbL, alpha2, desired_homology);
-  }
-  */
-  return 0;
-}
-
-int experiment3(int argc, char * const argv[])
-{
-  // Both witnesses and landmarks are given as input
-  
-  
-  return 0;
-}
-
-int main (int argc, char * const argv[])
-{
-  if (argc == 1) {
-    output_experiment_information(argv[0]);
-    return 1;
-  }
-  switch (atoi(argv[1])) {
-  case 0 :
-    return experiment0(argc, argv);
-    break;
-  case 1 :
-    return experiment1(argc, argv);
-    break;
-  case 2 :
-    return experiment2(argc, argv);
-    break;
-  case 3 :
-    return experiment3(argc, argv);
-    break;
-  default :
-    output_experiment_information(argv[0]);
-    return 1;
-  }
-}
diff --git a/src/Witness_complex/example/color-picker.cpp b/src/Witness_complex/example/color-picker.cpp
deleted file mode 100644
index fcea033e..00000000
--- a/src/Witness_complex/example/color-picker.cpp
+++ /dev/null
@@ -1,37 +0,0 @@
-#include <iostream>
-#include <assert.h>
-
-struct Color {
-  double r, g, b;
-  Color() : r(0), g(0), b(0) { }
-  Color(double r_, double g_, double b_)
-    : r(r_), g(g_), b(b_) { } 
-};
-
-// Convert a double in [0,1] to a color
-template< typename ColorType >
-void double_to_color(double t, ColorType & color) {
-  assert(t >= 0 && t <= 1);
-  double s = 1.0/6;
-  if (t >= 0 && t <= s)
-    color = ColorType(1, 6*t, 0);
-  else if (t > s && t <= 2*s)
-    color = ColorType(1-6*(t-s), 1, 0);
-  else if (t > 2*s && t <= 3*s)
-    color = ColorType(0, 1, 6*(t-2*s));
-  else if (t > 3*s && t <= 4*s)
-    color = ColorType(0, 1-6*(t-3*s), 1);
-  else if (t > 4*s && t <= 5*s)
-    color = ColorType(6*(t-4*s), 0, 1);
-  else
-    color = ColorType(1, 0, 1-6*(t-5*s));
-}
-
-int main (int argc, char * const argv[]) {
-  double number_to_convert = 0;
-  std::cout << "Enter a real number in [0,1]: ";
-  std::cin >> number_to_convert;
-  Color color;
-  double_to_color(number_to_convert, color);
-  std::cout << "The corresponding color in rgb is: (" << color.r << ", " << color.g << ", " << color.b << ")\n";
-}
diff --git a/src/Witness_complex/example/off_writer.h b/src/Witness_complex/example/off_writer.h
deleted file mode 100644
index 8cafe9e2..00000000
--- a/src/Witness_complex/example/off_writer.h
+++ /dev/null
@@ -1,11 +0,0 @@
-#ifndef OFF_WRITER_H_
-#define OFF_WRITER_H_
-
-
-class Off_Writer {
-
-private:
-
-}
-
-#endif
diff --git a/src/Witness_complex/example/relaxed_delaunay.cpp b/src/Witness_complex/example/relaxed_delaunay.cpp
deleted file mode 100644
index b0190446..00000000
--- a/src/Witness_complex/example/relaxed_delaunay.cpp
+++ /dev/null
@@ -1,180 +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 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/Simplex_tree.h>
-#include <gudhi/Relaxed_witness_complex.h>
-#include <gudhi/Dim_lists.h>
-#include <gudhi/reader_utils.h>
-#include <gudhi/Persistent_cohomology.h>
-#include "Landmark_choice_random_knn.h"
-#include "Landmark_choice_sparsification.h"
-
-#include <iostream>
-#include <fstream>
-#include <ctime>
-#include <utility>
-#include <algorithm>
-#include <set>
-#include <queue>
-#include <iterator>
-#include <string>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-#include <CGAL/Epick_d.h>
-#include <CGAL/Delaunay_triangulation.h>
-//#include <CGAL/Sphere_d.h>
-
-#include "generators.h"
-#include "output.h"
-
-using namespace Gudhi;
-using namespace Gudhi::witness_complex;
-using namespace Gudhi::persistent_cohomology;
-
-typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
-typedef K::Point_d Point_d;
-typedef K::Sphere_d Sphere_d;
-typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation;
-
-typedef std::vector<Point_d> Point_Vector;
-typedef Relaxed_witness_complex< Simplex_tree<> > RelaxedWitnessComplex;
-typedef Simplex_tree<>::Simplex_handle Simplex_handle;
-
-
-
-/**
- * \brief Customized version of read_points
- * which takes into account a possible nbP first line
- *
- */
-inline void
-read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) {
-  std::ifstream in_file(file_name.c_str(), std::ios::in);
-  if (!in_file.is_open()) {
-    std::cerr << "Unable to open file " << file_name << std::endl;
-    return;
-  }
-  std::string line;
-  double x;
-  while (getline(in_file, line)) {
-    std::vector< double > point;
-    std::istringstream iss(line);
-    while (iss >> x) {
-      point.push_back(x);
-    }
-    if (point.size() != 1)
-      points.push_back(point);
-  }
-  in_file.close();
-}
-
-int main (int argc, char * const argv[])
-{
-  if (argc != 4) {
-    std::cerr << "Usage: " << argv[0]
-              << " 1 file_name alpha limD\n";
-    return 0;
-  }
-  std::string file_name = argv[1];
-  double alpha2  = atof(argv[2]);
-  int limD = atoi(argv[3]);
-  
-  // Read points
-  Point_Vector point_vector;
-  read_points_cust(file_name, point_vector);
-  generate_points_random_box(point_vector, 200, 2);
-  write_points(file_name, point_vector);
-  
-  std::cout << "The file contains " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-  // 1. Compute Delaunay centers
-  Delaunay_triangulation delaunay(point_vector[0].size());
-  delaunay.insert(point_vector.begin(), point_vector.end());
-  Point_Vector del_centers;
-  for (auto f_it = delaunay.full_cells_begin(); f_it != delaunay.full_cells_end(); ++f_it) {
-    if (delaunay.is_infinite(f_it))
-      continue;
-    Point_Vector vertices;
-    for (auto v_it = f_it->vertices_begin(); v_it != f_it->vertices_end(); ++v_it)
-      vertices.push_back((*v_it)->point());
-    Sphere_d sphere(vertices.begin(), vertices.end());
-    del_centers.push_back(sphere.center());
-  }
-  std::cout << "Delaunay center count: " << del_centers.size() << ".\n";
-
-  // 2. Build Relaxed Witness Complex
-  std::vector<std::vector<int>> knn;
-  std::vector<std::vector<double>> distances;
-  Gudhi::witness_complex::build_distance_matrix(del_centers,   // aka witnesses
-                                                point_vector,  // aka landmarks
-                                                alpha2,
-                                                limD,
-                                                knn,
-                                                distances);
-  
-  write_wl("wl_distances.txt", distances);
-  Simplex_tree<> simplex_tree;
-  Gudhi::witness_complex::Relaxed_witness_complex<Simplex_tree<>> rwc(distances,
-                                                                      knn,
-                                                                      simplex_tree,
-                                                                      point_vector.size(),
-                                                                      alpha2,
-                                                                      limD);
-  std::vector<int> dim_simplices(limD+1);
-  for (auto sh: simplex_tree.complex_simplex_range()) {
-    dim_simplices[simplex_tree.dimension(sh)]++;
-  }
-  for (unsigned i =0; i != dim_simplices.size(); ++i)
-    std::cout << "dim[" << i << "]: " << dim_simplices[i] << " simplices.\n";
-
-  std::vector<int> landmarks_ind;
-  for (unsigned i = 0; i < point_vector.size(); ++i)
-    landmarks_ind.push_back(i);
-  write_witness_mesh(point_vector, landmarks_ind, simplex_tree, simplex_tree.complex_simplex_range(), true, true, "relaxed_delaunay.mesh");
-  
-  // 3. Check if the thing is Relaxed Delaunay
-  for (auto sh: simplex_tree.complex_simplex_range()) {
-    Point_Vector vertices;
-    for (auto v: simplex_tree.simplex_vertex_range(sh))
-      vertices.push_back(point_vector[v]);
-    Sphere_d sphere(vertices.begin(), vertices.end());
-    Point_d center = sphere.center();
-    double r2 = sphere.squared_radius();
-    typename K::Squared_distance_d dist2;
-    std::vector<int> v_inds;
-    for (auto v: simplex_tree.simplex_vertex_range(sh))
-      v_inds.push_back(v);
-    auto range_begin = std::begin(v_inds);
-    auto range_end = std::end(v_inds);
-    if (simplex_tree.dimension(sh) == (int)point_vector[0].size())
-    for (auto v: simplex_tree.complex_vertex_range())
-      if (std::find(range_begin, range_end, v) == range_end) {
-        if (dist2(point_vector[v], center) < r2 - alpha2)
-          std::cout << "WARNING! The vertex " << point_vector[v] << " is inside the (r2-alpha2)-ball (" << center << ", " << r2 << ") distance is " << dist2(point_vector[v], center) << "\n";
-      }
-  }
-}
diff --git a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp b/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp
deleted file mode 100644
index 067321ce..00000000
--- a/src/Witness_complex/example/relaxed_witness_complex_sphere.cpp
+++ /dev/null
@@ -1,461 +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) 2015  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 <iostream>
-#include <fstream>
-#include <ctime>
-#include <utility>
-#include <algorithm>
-#include <set>
-#include <queue>
-#include <iterator>
-
-#include <sys/types.h>
-#include <sys/stat.h>
-//#include <stdlib.h>
-
-//#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/Relaxed_witness_complex.h"
-#include "gudhi/reader_utils.h"
-#include "gudhi/Collapse/Collapse.h"
-//#include <boost/filesystem.hpp> 
-
-//#include <CGAL/Delaunay_triangulation.h>
-#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_incremental_neighbor_search.h>
-#include <CGAL/Kd_tree.h>
-#include <CGAL/Euclidean_distance.h>
-
-#include <CGAL/Kernel_d/Vector_d.h>
-#include <CGAL/point_generators_d.h>
-#include <CGAL/constructions_d.h>
-#include <CGAL/Fuzzy_sphere.h>
-#include <CGAL/Random.h>
-
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-using namespace Gudhi;
-//using namespace boost::filesystem;
-
-typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
-typedef K::FT FT;
-typedef K::Point_d Point_d;
-typedef CGAL::Search_traits<
-  FT, Point_d,
-  typename K::Cartesian_const_iterator_d,
-  typename K::Construct_cartesian_const_iterator_d> Traits_base;
-typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
-
-typedef std::vector< Vertex_handle > typeVectorVertex;
-
-//typedef std::pair<typeVectorVertex, Filtration_value> typeSimplex;
-//typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool;
-
-typedef CGAL::Search_traits_adapter<
-  std::ptrdiff_t, Point_d*, Traits_base> STraits;
-//typedef K TreeTraits;
-//typedef CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance > Euclidean_adapter;
-//typedef CGAL::Kd_tree<STraits> Kd_tree;
-typedef CGAL::Orthogonal_incremental_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> Neighbor_search;
-typedef Neighbor_search::Tree Tree;
-typedef Neighbor_search::Distance Distance;
-typedef Neighbor_search::iterator KNS_iterator;
-typedef Neighbor_search::iterator KNS_range;
-typedef boost::container::flat_map<int, int> Point_etiquette_map;
-typedef CGAL::Kd_tree<STraits> Tree2;
-
-typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
-
-typedef std::vector<Point_d> Point_Vector;
-
-//typedef K::Equal_d Equal_d;
-typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator;
-typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator;
-
-bool toric=false;
-
-/**
- * \brief Customized version of read_points
- * which takes into account a possible nbP first line
- *
- */
-inline void
-read_points_cust ( std::string file_name , Point_Vector & points)
-{  
-  std::ifstream in_file (file_name.c_str(),std::ios::in);
-  if(!in_file.is_open())
-    {
-      std::cerr << "Unable to open file " << file_name << std::endl;
-      return;
-    }
-  std::string line;
-  double x;
-  while( getline ( in_file , line ) )
-    {
-      std::vector< double > point;
-      std::istringstream iss( line );
-      while(iss >> x) { point.push_back(x); }
-      Point_d p(point.begin(), point.end());
-      if (point.size() != 1)
-        points.push_back(p);
-    }
-  in_file.close();
-}
-
-
-void generate_points_sphere(Point_Vector& W, int nbP, int dim)
-{
-  CGAL::Random_points_on_sphere_d<Point_d> rp(dim,1);
-  for (int i = 0; i < nbP; i++)
-    W.push_back(*rp++);
-}
-
-
-void write_wl( std::string file_name, std::vector< std::vector <int> > & WL)
-{
-  std::ofstream ofs (file_name, std::ofstream::out);
-  for (auto w : WL)
-    {
-      for (auto l: w)
-        ofs << l << " ";
-      ofs << "\n";
-    }
-  ofs.close();
-}
-
-void write_rl( std::string file_name, std::vector< std::vector <std::vector<int>::iterator> > & rl)
-{
-  std::ofstream ofs (file_name, std::ofstream::out);
-  for (auto w : rl)
-    {
-      for (auto l: w)
-        ofs << *l << " ";
-      ofs << "\n";
-    }
-  ofs.close();
-}
-
-std::vector<Point_d> convert_to_torus(std::vector< Point_d>& points)
-{
-  std::vector< Point_d > points_torus;
-  for (auto p: points)
-    {
-      FT theta = M_PI*p[0];
-      FT phi = M_PI*p[1];
-      std::vector<FT> p_torus;
-      p_torus.push_back((1+0.2*cos(theta))*cos(phi));
-      p_torus.push_back((1+0.2*cos(theta))*sin(phi));
-      p_torus.push_back(0.2*sin(theta));
-      points_torus.push_back(Point_d(p_torus));
-    }
-  return points_torus;
-}
-
-
-void write_points_torus( std::string file_name, std::vector< Point_d > & points)
-{
-  std::ofstream ofs (file_name, std::ofstream::out);
-  std::vector<Point_d> points_torus = convert_to_torus(points);
-  for (auto w : points_torus)
-    {
-      for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it)
-        ofs << *it << " ";
-      ofs << "\n";
-    }
-  ofs.close();
-}
-
-
-void write_points( std::string file_name, std::vector< Point_d > & points)
-{
-  if (toric) write_points_torus(file_name, points);
-  else
-    {
-      std::ofstream ofs (file_name, std::ofstream::out);
-      for (auto w : points)
-        {
-          for (auto it = w.cartesian_begin(); it != w.cartesian_end(); ++it)
-            ofs << *it << " ";
-          ofs << "\n";
-        }
-      ofs.close();
-    }
-}
-
-
-void write_edges_torus(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks)
-{
-  std::ofstream ofs (file_name, std::ofstream::out);
-  Point_Vector l_torus = convert_to_torus(landmarks);
-  for (auto u: witness_complex.complex_vertex_range())
-    for (auto v: witness_complex.complex_vertex_range())
-      {
-        typeVectorVertex edge = {u,v};
-        if (u < v && witness_complex.find(edge) != witness_complex.null_simplex())   
-          {
-            for (auto it = l_torus[u].cartesian_begin(); it != l_torus[u].cartesian_end(); ++it)
-              ofs << *it << " ";
-            ofs << "\n";
-            for (auto it = l_torus[v].cartesian_begin(); it != l_torus[v].cartesian_end(); ++it)
-              ofs << *it << " ";
-            ofs << "\n\n\n";
-          }
-    }
-  ofs.close();
-}
-
-void write_edges(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks)
-{
-  std::ofstream ofs (file_name, std::ofstream::out);
-  if (toric) write_edges_torus(file_name, witness_complex, landmarks);
-  else
-    {
-      for (auto u: witness_complex.complex_vertex_range())
-        for (auto v: witness_complex.complex_vertex_range())
-          {
-            typeVectorVertex edge = {u,v};
-            if (u < v && witness_complex.find(edge) != witness_complex.null_simplex())   
-              {
-                for (auto it = landmarks[u].cartesian_begin(); it != landmarks[u].cartesian_end(); ++it)
-                  ofs << *it << " ";
-                ofs << "\n";
-                for (auto it = landmarks[v].cartesian_begin(); it != landmarks[v].cartesian_end(); ++it)
-                  ofs << *it << " ";
-                ofs << "\n\n\n";
-              }
-          }
-      ofs.close();
-    }
-}
-
-
-/** Function that chooses landmarks from W and place it in the kd-tree L.
- *  Note: nbL hould be removed if the code moves to Witness_complex
- */
-void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind)
-{
-  std::cout << "Enter landmark choice to kd tree\n";
-  //std::vector<Point_d> landmarks;
-  int chosen_landmark;
-  //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false);
-  Point_d* p;
-  CGAL::Random rand;
-  for (int i = 0; i < nbL; i++)
-    {
-      //      while (!res.second)
-      //  {
-      do chosen_landmark = rand.get_int(0,nbP);
-      while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end());
-      //rand++;
-      //std::cout << "Chose " << chosen_landmark << std::endl;
-      p = &W[chosen_landmark];
-      //L_i.emplace(chosen_landmark,i);
-      //  }
-      landmarks.push_back(*p);
-      landmarks_ind.push_back(chosen_landmark);
-      //std::cout << "Added landmark " << chosen_landmark << std::endl;
-    }
- }
-
-
-void landmarks_to_witness_complex(Point_Vector &W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT alpha)
-{
-  //********************Preface: origin point
-  unsigned D = W[0].size();
-  std::vector<FT> orig_vector;
-  for (unsigned i = 0; i < D; i++)
-    orig_vector.push_back(0);
-  Point_d origin(orig_vector);
-  //Distance dist;
-  //dist.transformed_distance(0,1);
-  //******************** Constructing a WL matrix
-  int nbP = W.size();
-  int nbL = landmarks.size();
-  STraits traits(&(landmarks[0]));
-  Euclidean_distance ed;
-  std::vector< std::vector <int> > WL(nbP);
-  std::vector< std::vector< typename std::vector<int>::iterator > > ope_limits(nbP);
-  Tree L(boost::counting_iterator<std::ptrdiff_t>(0),
-         boost::counting_iterator<std::ptrdiff_t>(nbL),
-         typename Tree::Splitter(),
-         traits);
-
-  std::cout << "Enter (D+1) nearest landmarks\n";
-  //std::cout << "Size of the tree is " << L.size() << std::endl;
-  for (int i = 0; i < nbP; i++)
-    {
-      //std::cout << "Entered witness number " << i << std::endl;
-      Point_d& w = W[i];
-      std::queue< typename std::vector<int>::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark
-      Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])));
-      Neighbor_search::iterator search_it = search.begin();
-      
-      //Incremental search and filling WL
-      while (WL[i].size() < D)
-        WL[i].push_back((search_it++)->first);
-      FT dtow = ed.transformed_distance(w, landmarks[WL[i][D-1]]);
-      while (search_it->second < dtow + alpha)
-        WL[i].push_back((search_it++)->first);
-
-      //Filling the (1+epsilon)-limits table
-      for (std::vector<int>::iterator wl_it = WL[i].begin(); wl_it != WL[i].end(); ++wl_it)
-        {
-          ope_queue.push(wl_it);
-          FT d_to_curr_l = ed.transformed_distance(w, landmarks[*wl_it]);
-          //std::cout << "d_to_curr_l=" << d_to_curr_l << std::endl;
-          //std::cout << "d_to_front+alpha=" << d_to_curr_l << std::endl;
-          while (d_to_curr_l > alpha + ed.transformed_distance(w, landmarks[*(ope_queue.front())]))
-            {
-              ope_limits[i].push_back(wl_it);
-              ope_queue.pop();
-            }
-        }
-      while (ope_queue.size() > 0)
-        {
-          ope_limits[i].push_back(WL[i].end());
-          ope_queue.pop();
-        }
-      //std::cout << "Safely constructed a point\n";
-      ////Search D+1 nearest neighbours from the tree of landmarks L
-      /*
-      if (w[0]>0.95)
-        std::cout << i << std::endl;
-      */
-      //K_neighbor_search search(L, w, D, FT(0), true,
-      //                         CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])) );
-      //std::cout << "Safely found nearest landmarks\n";
-      /* 
-      for(K_neighbor_search::iterator it = search.begin(); it != search.end(); ++it)
-        {
-          //std::cout << "Entered KNN_it with point at distance " << it->second << "\n";
-          //Point_etiquette_map::iterator itm = L_i.find(it->first);
-          //assert(itm != L_i.end());
-          //std::cout << "Entered KNN_it with point at distance " << it->second << "\n";
-          WL[i].push_back(it->first);
-          //std::cout << "ITFIRST " << it->first << std::endl;
-          //std::cout << i << " " << it->first << ": " << it->second << std::endl; 
-        }
-      */
-    }
-  //std::cout << "\n";
-  
-  //std::string out_file = "wl_result";
-  write_wl("wl_result",WL);
-  write_rl("rl_result",ope_limits);
-  
-  //******************** Constructng a witness complex
-  std::cout << "Entered witness complex construction\n";
-  Witness_complex<> witnessComplex;
-  witnessComplex.setNbL(nbL);
-  witnessComplex.relaxed_witness_complex(WL, ope_limits);
-  char buffer[100];
-  int i = sprintf(buffer,"stree_result.txt");
-  
-  if (i >= 0)
-    {
-      std::string out_file = (std::string)buffer;
-      std::ofstream ofs (out_file, std::ofstream::out);
-      witnessComplex.st_to_file(ofs);
-      ofs.close();
-    }
-  write_edges("landmarks/edges", witnessComplex, landmarks);
-  std::cout << Distance().transformed_distance(Point_d(std::vector<double>({0.1,0.1})), Point_d(std::vector<double>({1.9,1.9}))) << std::endl;
-}
-
-
-int main (int argc, char * const argv[])
-{
-  
-  if (argc != 5)
-    {
-      std::cerr << "Usage: " << argv[0]
-                << " nbP nbL dim alpha\n";
-      return 0;
-    }
-  /*
-  boost::filesystem::path p;
-  for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  int nbP       = atoi(argv[1]);
-  int nbL       = atoi(argv[2]);
-  int dim       = atoi(argv[3]);
-  double alpha  = atof(argv[4]);
-  //clock_t start, end;
-  //Construct the Simplex Tree
-  Witness_complex<> witnessComplex;
- 
-  std::cout << "Let the carnage begin!\n";
-  Point_Vector point_vector;
-  //read_points_cust(file_name, point_vector);
-  generate_points_sphere(point_vector, nbP, dim);
-  /*
-  for (auto &p: point_vector)
-    {
-      assert(std::count(point_vector.begin(),point_vector.end(),p) == 1);
-    }
-  */
-  //std::cout << "Successfully read the points\n";
-  //witnessComplex.setNbL(nbL);
-  Point_Vector L;
-  std::vector<int> chosen_landmarks;
-  landmark_choice(point_vector, nbP, nbL, L, chosen_landmarks);
-  //start = clock();
-  
-  write_points("landmarks/initial_pointset",point_vector);
-  write_points("landmarks/initial_landmarks",L);
-  
-  landmarks_to_witness_complex(point_vector, L, chosen_landmarks, alpha);
-  //end = clock();
-  
-  /*
-  std::cout << "Landmark choice took "
-            << (double)(end-start)/CLOCKS_PER_SEC << " s. \n";
-  start = clock();
-  witnessComplex.witness_complex(WL);
-  //
-  end = clock();
-  std::cout << "Howdy world! The process took "
-       << (double)(end-start)/CLOCKS_PER_SEC << " s. \n";
-  */
-  
-  /*
-  out_file = "output/"+file_name+"_"+argv[2]+".stree";
-  std::ofstream ofs (out_file, std::ofstream::out);
-  witnessComplex.st_to_file(ofs);
-  ofs.close();
-
-  out_file = "output/"+file_name+"_"+argv[2]+".badlinks";
-  std::ofstream ofs2(out_file, std::ofstream::out);
-  witnessComplex.write_bad_links(ofs2);
-  ofs2.close();
-  */
-}
diff --git a/src/Witness_complex/example/relaxed_witness_persistence.cpp b/src/Witness_complex/example/relaxed_witness_persistence.cpp
deleted file mode 100644
index b4837594..00000000
--- a/src/Witness_complex/example/relaxed_witness_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):       Siargey Kachanovich
- *
- *    Copyright (C) 2016  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/Simplex_tree.h>
-#include <gudhi/Relaxed_witness_complex.h>
-#include <gudhi/Dim_lists.h>
-#include <gudhi/reader_utils.h>
-#include <gudhi/Persistent_cohomology.h>
-#include "Landmark_choice_random_knn.h"
-#include "Landmark_choice_sparsification.h"
-
-#include <iostream>
-#include <fstream>
-#include <ctime>
-#include <utility>
-#include <algorithm>
-#include <set>
-#include <queue>
-#include <iterator>
-#include <string>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-#include "generators.h"
-#include "output.h"
-
-using namespace Gudhi;
-using namespace Gudhi::witness_complex;
-using namespace Gudhi::persistent_cohomology;
-
-typedef std::vector<Point_d> Point_Vector;
-typedef Relaxed_witness_complex< Simplex_tree<> > RelaxedWitnessComplex;
-typedef Simplex_tree<>::Simplex_handle Simplex_handle;
-
-/**
- * \brief Customized version of read_points
- * which takes into account a possible nbP first line
- *
- */
-inline void
-read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) {
-  std::ifstream in_file(file_name.c_str(), std::ios::in);
-  if (!in_file.is_open()) {
-    std::cerr << "Unable to open file " << file_name << std::endl;
-    return;
-  }
-  std::string line;
-  double x;
-  while (getline(in_file, line)) {
-    std::vector< double > point;
-    std::istringstream iss(line);
-    while (iss >> x) {
-      point.push_back(x);
-    }
-    if (point.size() != 1)
-      points.push_back(point);
-  }
-  in_file.close();
-}
-
-void output_experiment_information(char * const file_name)
-{
-    std::cout << "Enter a valid experiment number. Usage: "
-              << file_name << " exp_no options\n";
-    std::cout << "Experiment description:\n"
-              << "0 nbP nbL dim alpha limD mu_epsilon: "
-              << "Build persistence diagram on relaxed witness complex "
-              << "built from a point cloud on (dim-1)-dimensional sphere "
-              << "consisting of nbP witnesses and nbL landmarks. "
-              << "The maximal relaxation is alpha and the limit on simplicial complex "
-              << "dimension is limD.\n";
-    std::cout << "1 file_name nbL alpha limD: "
-              << "Build persistence diagram on relaxed witness complex "
-              << "build from a point cloud stored in a file and nbL landmarks. "
-              << "The maximal relaxation is alpha and the limit on simplicial complex dimension is limD\n";
-}
-
-void rw_experiment(Point_Vector & point_vector, int nbL, FT alpha, int limD, FT mu_epsilon = 0.1)
-{
-  clock_t start, end;
-  Simplex_tree<> simplex_tree;
-
-  // Choose landmarks
-  std::vector<std::vector< int > > knn;
-  std::vector<std::vector< FT > > distances;
-  start = clock();
-  //Gudhi::witness_complex::landmark_choice_by_random_knn(point_vector, nbL, alpha, limD, knn, distances);
-  std::vector<Point_d> landmarks;
-  Gudhi::witness_complex::landmark_choice_by_sparsification(point_vector, nbL, mu_epsilon, landmarks);
-  Gudhi::witness_complex::build_distance_matrix(point_vector,   // aka witnesses
-                                                landmarks,  // aka landmarks
-                                                alpha,
-                                                limD,
-                                                knn,
-                                                distances);
-  end = clock();
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Choice of " << nbL << " landmarks took "
-            << time << " s. \n";
-  // Compute witness complex
-  start = clock();
-  RelaxedWitnessComplex rw(distances,
-                           knn,
-                           simplex_tree,
-                           nbL,
-                           alpha,
-                           limD);
-  end = clock();
-  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Relaxed witness complex for " << nbL << " landmarks took "
-            << time << " s. \n";
-  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
-  // std::cout << simplex_tree << "\n";
-  
-  // Compute the persistence diagram of the complex
-  simplex_tree.set_dimension(limD);
-  persistent_cohomology::Persistent_cohomology< Simplex_tree<>, Field_Zp > pcoh(simplex_tree, true);
-  int p = 3;
-  pcoh.init_coefficients( p ); //initilizes the coefficient field for homology
-  start = clock();
-  pcoh.compute_persistent_cohomology( alpha/1000 );
-  end = clock();
-  time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Persistence diagram took "
-            << time << " s. \n";
-  pcoh.output_diagram();
-
-  int chi = 0;
-  for (auto sh: simplex_tree.complex_simplex_range())
-    chi += 1-2*(simplex_tree.dimension(sh)%2);
-  std::cout << "Euler characteristic is " << chi << std::endl;
-
-}
-
-
-int experiment0 (int argc, char * const argv[])
-{
-  if (argc != 8) {
-    std::cerr << "Usage: " << argv[0]
-              << " 0 nbP nbL dim alpha limD mu_epsilon\n";
-    return 0;
-  }
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  int nbP       = atoi(argv[2]);
-  int nbL       = atoi(argv[3]);
-  int dim       = atoi(argv[4]);
-  double alpha  = atof(argv[5]);
-  int limD      = atoi(argv[6]);
-  double mu_epsilon = atof(argv[7]);
-
-  // Read the point file
-  Point_Vector point_vector;
-  generate_points_sphere(point_vector, nbP, dim);
-  std::cout << "Successfully generated " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-  rw_experiment(point_vector, nbL, alpha, limD);
-}
-
-int experiment1 (int argc, char * const argv[])
-{
-  if (argc != 3) {
-    std::cerr << "Usage: " << argv[0]
-              << " 1 file_name\n";
-    return 0;
-  }
-  /*
-    boost::filesystem::path p;
-    for (; argc > 2; --argc, ++argv)
-    p /= argv[1];
-  */
-  
-  std::string file_name = argv[2];
-
-  // Read the point file
-  Point_Vector point_vector;
-  read_points_cust(file_name, point_vector);
-  std::cout << "The file contains " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-  bool ok = false;
-  int nbL, limD;
-  double alpha;
-  while (!ok) {
-    std::cout << "Relaxed witness complex: parameters nbL, alpha, limD.\n";
-    std::cout << "Enter nbL: ";
-    std::cin >> nbL;
-    std::cout << "Enter alpha: ";
-    std::cin >> alpha;
-    std::cout << "Enter limD: ";
-    std::cin >> limD;
-    std::cout << "Start relaxed witness complex...\n";
-    rw_experiment(point_vector, nbL, alpha, limD);
-    std::cout << "Is the result correct? [y/n]: ";
-    char answer;
-    std::cin >> answer;
-    switch (answer) {
-    case 'n':
-      ok = false; break;
-    default :
-      ok = true; break;
-    }
-  }
-  // ok = false;
-  // while (!ok) {
-  //   std::cout << "Rips complex: parameters threshold, limD.\n";
-  //   std::cout << "Enter threshold: ";
-  //   std::cin >> alpha;
-  //   std::cout << "Enter limD: ";
-  //   std::cin >> limD;
-  //   std::cout << "Start Rips complex...\n";
-  //   rips_experiment(point_vector, alpha, limD);
-  //   std::cout << "Is the result correct? [y/n]: ";
-  //   char answer;
-  //   std::cin >> answer;
-  //   switch (answer) {
-  //   case 'n':
-  //     ok = false; break;
-  //   default :
-  //     ok = true; break;
-  //   }
-  // }
-}
-
-int main (int argc, char * const argv[])
-{
-  if (argc == 1) {
-    output_experiment_information(argv[0]);
-    return 1;
-  }
-  switch (atoi(argv[1])) {
-  case 0 :
-    return experiment0(argc, argv);
-  case 1 :
-    return experiment1(argc, argv);
-  default :
-    output_experiment_information(argv[0]);
-    return 1;
-  }
-}
diff --git a/src/Witness_complex/example/strange_sphere.cpp b/src/Witness_complex/example/strange_sphere.cpp
deleted file mode 100644
index 9188bd8e..00000000
--- a/src/Witness_complex/example/strange_sphere.cpp
+++ /dev/null
@@ -1,219 +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) 2015  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/>.
- */
-#define BOOST_PARAMETER_MAX_ARITY 12
-
-
-#include <sys/types.h>
-#include <sys/stat.h>
-
-#include <gudhi/Simplex_tree.h>
-#include <gudhi/Strange_witness_complex.h>
-#include <gudhi/Landmark_choice_by_random_point.h>
-#include <gudhi/reader_utils.h>
-#include "output.h"
-
-#include <iostream>
-#include <fstream>
-#include <ctime>
-#include <utility>
-#include <string>
-#include <vector>
-
-#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_incremental_neighbor_search.h>
-#include <CGAL/Kd_tree.h>
-#include <CGAL/Euclidean_distance.h>
-
-#include <CGAL/Kernel_d/Vector_d.h>
-#include <CGAL/point_generators_d.h>
-#include <CGAL/constructions_d.h>
-#include <CGAL/Fuzzy_sphere.h>
-#include <CGAL/Random.h>
-
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-
-#include "generators.h"
-
-using namespace Gudhi;
-using namespace Gudhi::witness_complex;
-
-typedef std::vector< Vertex_handle > typeVectorVertex;
-typedef Strange_witness_complex< Simplex_tree<> > WitnessComplex;
-
-typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
-typedef K::FT FT;
-typedef K::Point_d Point_d;
-typedef CGAL::Search_traits<
-  FT, Point_d,
-  typename K::Cartesian_const_iterator_d,
-  typename K::Construct_cartesian_const_iterator_d> Traits_base;
-typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
-
-typedef CGAL::Search_traits_adapter<
-  std::ptrdiff_t, Point_d*, Traits_base> STraits;
-typedef CGAL::Orthogonal_incremental_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> Neighbor_search;
-typedef Neighbor_search::Tree Tree;
-typedef Neighbor_search::Distance Distance;
-typedef Neighbor_search::iterator KNS_iterator;
-typedef Neighbor_search::iterator KNS_range;
-typedef boost::container::flat_map<int, int> Point_etiquette_map;
-typedef CGAL::Kd_tree<STraits> Tree2;
-typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
-typedef std::vector<Point_d> Point_Vector;
-
-/** Function that chooses landmarks from W and place it in the kd-tree L.
- *  Note: nbL hould be removed if the code moves to Witness_complex
- */
-void landmark_choice(Point_Vector &W, int nbP, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind)
-{
-  std::cout << "Enter landmark choice to kd tree\n";
-  //std::vector<Point_d> landmarks;
-  int chosen_landmark;
-  //std::pair<Point_etiquette_map::iterator,bool> res = std::make_pair(L_i.begin(),false);
-  Point_d* p;
-  CGAL::Random rand;
-  for (int i = 0; i < nbL; i++)
-    {
-      //      while (!res.second)
-      //  {
-      do chosen_landmark = rand.get_int(0,nbP);
-      while (std::find(landmarks_ind.begin(), landmarks_ind.end(), chosen_landmark) != landmarks_ind.end());
-      //rand++;
-      //std::cout << "Chose " << chosen_landmark << std::endl;
-      p = &W[chosen_landmark];
-      //L_i.emplace(chosen_landmark,i);
-      //  }
-      landmarks.push_back(*p);
-      landmarks_ind.push_back(chosen_landmark);
-      //std::cout << "Added landmark " << chosen_landmark << std::endl;
-    }
- }
-
-void landmarks_to_witness_complex(Point_Vector &W, int nbL)
-{
-  //********************Preface: origin point
-  unsigned D = W[0].size();
-  std::vector<FT> orig_vector;
-  for (unsigned i = 0; i < D; i++)
-    orig_vector.push_back(0);
-  Point_d origin(orig_vector);
-  //Distance dist;
-  //dist.transformed_distance(0,1);
-  //******************** Constructing a WL matrix
-  int nbP = W.size();
-  Point_Vector landmarks;
-  std::vector<int> landmarks_ind;
-  landmark_choice(W, nbP, nbL, landmarks, landmarks_ind);
-  
-  
-  STraits traits(&(landmarks[0]));
-  Euclidean_distance ed;
-  std::vector< std::vector <int> > WL(nbP);
-  std::vector< std::vector< typename std::vector<int>::iterator > > ope_limits(nbP);
-  Tree L(boost::counting_iterator<std::ptrdiff_t>(0),
-         boost::counting_iterator<std::ptrdiff_t>(nbL),
-         typename Tree::Splitter(),
-         traits);
-
-  std::cout << "Enter (D+1) nearest landmarks\n";
-  //std::cout << "Size of the tree is " << L.size() << std::endl;
-  for (int i = 0; i < nbP; i++) {
-    //std::cout << "Entered witness number " << i << std::endl;
-    Point_d& w = W[i];
-    std::queue< typename std::vector<int>::iterator > ope_queue; // queue of points at (1+epsilon) distance to current landmark
-    Neighbor_search search(L, w, FT(0), true, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>(&(landmarks[0])));
-    Neighbor_search::iterator search_it = search.begin();
-      
-    //Incremental search and filling WL
-    while (WL[i].size() < D)
-      WL[i].push_back((search_it++)->first);
-  }
-  write_wl("WL.txt", WL); //!
-  //******************** Constructng a witness complex
-  std::cout << "Entered witness complex construction\n";
-
-  Simplex_tree<> simplex_tree;
-  WitnessComplex(WL, simplex_tree, nbL, W[0].size()-1);
-  //std::cout << simplex_tree << "\n"; //!
-  write_witness_mesh(W, landmarks_ind, simplex_tree, false, true);
-}
-
-
-
-/** 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) {
-  std::ofstream ofs(filename, std::ofstream::out);
-  for (auto entry : data)
-    ofs << entry.first << ", " << entry.second << "\n";
-  ofs.close();
-}
-
-
-
-int main(int argc, char * const argv[]) {
-  if (argc != 4) {
-    std::cerr << "Usage: " << argv[0]
-        << " nbP nbL dim\n";
-    return 0;
-  }
-
-  int nbP = atoi(argv[1]);
-  int nbL = atoi(argv[2]);
-  int dim = atoi(argv[3]);
-  clock_t start, end;
-
-  // Construct the Simplex Tree
-  Simplex_tree<> simplex_tree;
-
-  std::vector< std::pair<int, double> > l_time;
-
-  // Read the point file
-  Point_Vector point_vector;
-  generate_points_sphere(point_vector, nbP, dim);
-  std::cout << "Successfully generated " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
-
-  // Choose landmarks
-  start = clock();
-  std::vector<std::vector< int > > knn;
-  
-  landmarks_to_witness_complex(point_vector, nbL);
-
-  end = clock();
-  double time = static_cast<double>(end - start) / CLOCKS_PER_SEC;
-  std::cout << "Witness complex for " << nbL << " landmarks took "
-            << time << " s. \n";
-  l_time.push_back(std::make_pair(nbP, time));
-  //write_data(l_time, "w_time.dat");
-}
diff --git a/src/Witness_complex/example/witness_complex_from_file.cpp b/src/Witness_complex/example/witness_complex_from_file.cpp
index 5e9f0e81..e5859b2a 100644
--- a/src/Witness_complex/example/witness_complex_from_file.cpp
+++ b/src/Witness_complex/example/witness_complex_from_file.cpp
@@ -25,18 +25,22 @@
 
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Witness_complex.h>
-#include <gudhi/Construct_closest_landmark_table.h>
 #include <gudhi/pick_n_random_points.h>
 #include <gudhi/reader_utils.h>
 
+#include <CGAL/Epick_d.h>
+
 #include <iostream>
 #include <fstream>
 #include <ctime>
 #include <string>
 #include <vector>
 
+typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
+typedef typename K::Point_d Point_d;
+typedef typename Gudhi::witness_complex::Witness_complex<K> Witness_complex;
 typedef std::vector< Vertex_handle > typeVectorVertex;
-typedef std::vector< std::vector <double> > Point_Vector;
+typedef std::vector< Point_d > Point_vector;
 
 /**
  * \brief Customized version of read_points
@@ -44,7 +48,7 @@ typedef std::vector< std::vector <double> > Point_Vector;
  *
  */
 inline void
-read_points_cust(std::string file_name, std::vector< std::vector< double > > & points) {
+read_points_cust(std::string file_name, Point_vector & points) {
   std::ifstream in_file(file_name.c_str(), std::ios::in);
   if (!in_file.is_open()) {
     std::cerr << "Unable to open file " << file_name << std::endl;
@@ -59,44 +63,44 @@ read_points_cust(std::string file_name, std::vector< std::vector< double > > & p
       point.push_back(x);
     }
     if (point.size() != 1)
-      points.push_back(point);
+      points.push_back(Point_d(point));
   }
   in_file.close();
 }
 
 int main(int argc, char * const argv[]) {
-  if (argc != 3) {
+  if (argc != 4) {
     std::cerr << "Usage: " << argv[0]
-        << " path_to_point_file nbL \n";
+        << " path_to_point_file nbL alpha^2\n";
     return 0;
   }
 
   std::string file_name = argv[1];
   int nbL = atoi(argv[2]);
+  double alpha2 = atof(argv[3]);
   clock_t start, end;
 
   // Construct the Simplex Tree
   Gudhi::Simplex_tree<> simplex_tree;
 
   // Read the point file
-  Point_Vector point_vector, landmarks;
+  Point_vector point_vector, landmarks;
   read_points_cust(file_name, point_vector);
   std::cout << "Successfully read " << point_vector.size() << " points.\n";
-  std::cout << "Ambient dimension is " << point_vector[0].size() << ".\n";
+  std::cout << "Ambient dimension is " << point_vector[0].dimension() << ".\n";
 
   // Choose landmarks
-  start = clock();
-  std::vector<std::vector< int > > knn;
-  Gudhi::subsampling::pick_n_random_points(point_vector, 100, std::back_inserter(landmarks));
-  Gudhi::witness_complex::construct_closest_landmark_table(point_vector, landmarks, knn);
-  end = clock();
-  std::cout << "Landmark choice for " << nbL << " landmarks took "
-      << static_cast<double>(end - start) / CLOCKS_PER_SEC << " s. \n";
+  Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks));
 
   // Compute witness complex
   start = clock();
-  Gudhi::witness_complex::witness_complex(knn, nbL, point_vector[0].size(), simplex_tree);
+  Witness_complex witness_complex(landmarks.begin(),
+                                  landmarks.end(),
+                                  point_vector.begin(),
+                                  point_vector.end());
+  witness_complex.create_complex(simplex_tree, alpha2);
   end = clock();
   std::cout << "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/src/Witness_complex/include/gudhi/A0_complex.h b/src/Witness_complex/include/gudhi/A0_complex.h
deleted file mode 100644
index 2018e1c8..00000000
--- a/src/Witness_complex/include/gudhi/A0_complex.h
+++ /dev/null
@@ -1,337 +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) 2015  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/>.
- */
-
-#ifndef A0_COMPLEX_H_
-#define A0_COMPLEX_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-// Needed for nearest neighbours
-#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /** \addtogroup simplex_tree
-   *  Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$:
-   *  \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$
-   *  and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that
-   *  w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. 
-   */
-template< class Simplicial_complex >
-class A0_complex {
-
-private:    
-  struct Active_witness {
-    int witness_id;
-    int landmark_id;
-    
-    Active_witness(int witness_id_, int landmark_id_)
-      : witness_id(witness_id_),
-        landmark_id(landmark_id_) { }
-    };
-
-private:
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Simplicial_complex::Vertex_handle Vertex_handle;
-  typedef typename Simplicial_complex::Filtration_value FT;
-
-  typedef std::vector< double > Point_t;
-  typedef std::vector< Point_t > Point_Vector;
-
-  // typedef typename Simplicial_complex::Filtration_value Filtration_value;
-  typedef std::vector< Vertex_handle > typeVectorVertex;
-  typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex;
-  typedef std::pair< Simplex_handle, bool > typePairSimplexBool;
-
-  typedef int Witness_id;
-  typedef int Landmark_id;
-  typedef std::list< Vertex_handle > ActiveWitnessList;
-
- private:
-  int nbL;  // Number of landmarks
-  Simplicial_complex& sc;  // Simplicial complex
-  
-  public:
-    /** @name Constructor
-     */
-
-    //@{
-  
-    /**
-     *  \brief Iterative construction of the relaxed witness complex.
-     *  \details The witness complex is written in sc_ basing on a matrix knn
-     *  of k nearest neighbours of the form {witnesses}x{landmarks} and 
-     *  and a matrix distances of distances to these landmarks from witnesses.
-     *  The parameter alpha defines relaxation and
-     *  limD defines the 
-     *
-     *  The line lengths in one given matrix can differ, 
-     *  however both matrices have the same corresponding line lengths.
-     *
-     *  The type KNearestNeighbors can be seen as 
-     *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-     *  Witness_range and Closest_landmark_range are random access ranges.
-     *
-     *  Constructor takes into account at most (dim+1) 
-     *  first landmarks from each landmark range to construct simplices.
-     *
-     *  Landmarks are supposed to be in [0,nbL_-1]
-     */
-    
-    template< typename KNearestNeighbours >
-    A0_complex(std::vector< std::vector<double> > const & distances,
-                            KNearestNeighbours const & knn,
-                            Simplicial_complex & sc_, 
-                            int nbL_,
-                            double alpha2,
-                            unsigned limD) : nbL(nbL_), sc(sc_) {
-      int nbW = knn.size();
-      typeVectorVertex vv;
-      //int counter = 0;
-      /* The list of still useful witnesses
-       * it will diminuish in the course of iterations
-       */
-      ActiveWitnessList active_w;// = new ActiveWitnessList();
-      for (int i = 0; i != nbL; ++i) {
-        // initial fill of 0-dimensional simplices
-        // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
-        //counter++;
-        vv = {i};
-        sc.insert_simplex(vv, Filtration_value(0.0));
-        /* TODO Error if not inserted : normally no need here though*/
-      }
-      for (int i=0; i != nbW; ++i) {
-        // int i_end = limD+1;
-        // if (knn[i].size() < limD+1)
-        //   i_end = knn[i].size();
-        // double dist_wL = *(distances[i].begin());
-        // while (distances[i][i_end] > dist_wL + alpha2)
-        //   i_end--;
-        // add_all_witnessed_faces(distances[i].begin(),
-        //                         knn[i].begin(),
-        //                         knn[i].begin() + i_end + 1);
-        unsigned j_end = 0;
-        while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) {
-          std::vector<int> simplex;
-          for (unsigned j = 0; j <= j_end; ++j)
-            simplex.push_back(knn[i][j]);
-          assert(distances[i][j_end] - distances[i][0] >= 0);
-          sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]);
-          j_end++; 
-        }
-      }
-      sc.set_dimension(limD);
-    }
-
-  //@}
-  
-private:
-  /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness
-   * Iterator is needed to know until how far we can take landmarks to form simplexes
-   * simplex is the prefix of the simplexes to insert
-   * The output value indicates if the witness rests active or not
-   */
-  void add_all_witnessed_faces(std::vector<double>::const_iterator curr_d,
-                               std::vector<int>::const_iterator curr_l,
-                               std::vector<int>::const_iterator end)
-  {
-    std::vector<int> simplex;
-    std::vector<int>::const_iterator l_end = curr_l;
-    for (; l_end != end; ++l_end) {
-      std::vector<int>::const_iterator l_it = curr_l;
-      std::vector<double>::const_iterator d_it = curr_d;
-      simplex = {};
-      for (; l_it != l_end; ++l_it, ++d_it)
-        simplex.push_back(*l_it);
-      sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d);
-    }
-  }
-  
-  /** \brief Check if the facets of the k-dimensional simplex witnessed 
-   *  by witness witness_id are already in the complex.
-   *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
-   */
-  bool all_faces_in(std::vector<int>& simplex, double* filtration_value)
-  {
-    std::vector< int > facet;
-    for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it)
-      {
-        facet.clear();
-        for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it)
-          if (it != not_it)
-            facet.push_back(*it);
-        Simplex_handle facet_sh = sc.find(facet);
-        if (facet_sh == sc.null_simplex())
-          return false;
-        else if (sc.filtration(facet_sh) > *filtration_value)
-          *filtration_value = sc.filtration(facet_sh);
-      }
-    return true;
-  }
-
-  bool is_face(Simplex_handle face, Simplex_handle coface)
-  {
-    // vertex range is sorted in decreasing order
-    auto fvr = sc.simplex_vertex_range(face);
-    auto cfvr = sc.simplex_vertex_range(coface);
-    auto fv_it = fvr.begin();
-    auto cfv_it = cfvr.begin();
-    while (fv_it != fvr.end() && cfv_it != cfvr.end()) {
-      if (*fv_it < *cfv_it)
-        ++cfv_it;
-      else if (*fv_it == *cfv_it) {
-        ++cfv_it;
-        ++fv_it;
-      }
-      else
-        return false;
-      
-    }
-    return (fv_it == fvr.end());
-  }
-
-  // void erase_simplex(Simplex_handle sh)
-  // {
-  //   auto siblings = sc.self_siblings(sh);
-  //   auto oncles = siblings->oncles();
-  //   int prev_vertex = siblings->parent();
-  //   siblings->members().erase(sh->first);
-  //   if (siblings->members().empty()) {
-  //     typename typedef Simplicial_complex::Siblings Siblings;
-  //     oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex));
-  //     assert(!sc.has_children(oncles->members().find(prev_vertex)));
-  //     //delete siblings;
-  //   }
-
-  // }
-  
-  void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh)
-  {
-    erase_simplex(coface_sh);
-    erase_simplex(face_sh);
-  }
-
-public:
-  // void collapse(std::vector<Simplex_handle>& simplices)
-  // {
-  //   // Get a vector of simplex handles ordered by filtration value
-  //   std::cout << sc << std::endl;
-  //   //std::vector<Simplex_handle> simplices;
-  //   for (Simplex_handle sh: sc.filtration_simplex_range())
-  //     simplices.push_back(sh);
-  //   // std::sort(simplices.begin(),
-  //   //           simplices.end(),
-  //   //           [&](Simplex_handle sh1, Simplex_handle sh2)
-  //   //           { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2);
-  //   //             return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); });
-  //   // Double iteration
-  //   auto face_it = simplices.rbegin();
-  //   while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) {
-  //     int coface_count = 0;
-  //     auto reduced_coface = simplices.rbegin();
-  //     for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it)
-  //       if (face_it != coface_it && is_face(*face_it, *coface_it)) {
-  //         coface_count++;
-  //         if (coface_count == 1)
-  //           reduced_coface = coface_it;
-  //         else
-  //           break;
-  //       }
-  //     if (coface_count == 1) {
-  //       std::cout << "Erase ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base())))
-  //         std::cout << v << " ";
-        
-  //       simplices.erase(--(reduced_coface.base()));
-  //       //elementary_collapse(*face_it, *reduced_coface);
-  //       std::cout << ") and then ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--face_it.base())))
-  //         std::cout << v << " ";
-  //       std::cout << ")\n";
-  //       simplices.erase(--((face_it++).base()));
-  //       //face_it = simplices.rbegin();
-  //       //std::cout << "Size of vector: " << simplices.size() << "\n";
-  //     } 
-  //     else
-  //       face_it++;
-  //   }
-  //   sc.initialize_filtration();
-  //   //std::cout << sc << std::endl;
-  // }
-
-  template <class Dim_lists>
-  void collapse(Dim_lists& dim_lists)
-  {
-    dim_lists.collapse();
-  }
-  
-private:
-  
-  /** Collapse recursively boundary faces of the given simplex
-   *  if its filtration is bigger than alpha_lim.  
-   */
-  void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim)
-  {
-    for (Simplex_handle face_it : sc.boundary_simplex_range()) {
-      
-    }
-      
-  }
-  
-}; //class Relaxed_witness_complex
-
-} // namespace witness_complex
-  
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h b/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h
deleted file mode 100644
index ef711c34..00000000
--- a/src/Witness_complex/include/gudhi/Construct_closest_landmark_table.h
+++ /dev/null
@@ -1,90 +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) 2015  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/>.
- */
-
-#ifndef CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_
-#define CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_
-
-#include <boost/range/size.hpp>
-
-#include <queue>  // for priority_queue<>
-#include <utility>  // for pair<>
-#include <iterator>
-#include <vector>
-#include <set>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /**
-   *  \ingroup witness_complex
-   * \brief Construct the closest landmark tables for all witnesses.
-   *  \details Output a table 'knn', each line of which represents a witness and
-   *   consists of landmarks sorted by
-   *   euclidean distance from the corresponding witness.
-   *
-   *  The type WitnessContainer is a random access range and
-   *  the type LandmarkContainer is a range.
-   *  The type KNearestNeighbors can be seen as 
-   *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-   *  Witness_range and Closest_landmark_range are random access ranges and
-   *  Vertex_handle is the label type of a vertex in a simplicial complex.
-   *  Closest_landmark_range needs to have push_back operation.
-   */
-
-  template <typename WitnessContainer,
-            typename LandmarkContainer,
-            typename KNearestNeighbours>
-  void construct_closest_landmark_table(WitnessContainer const &points,
-                                        LandmarkContainer const &landmarks,
-                                        KNearestNeighbours &knn) {
-    int nbP = boost::size(points);
-    assert(nbP >= boost::size(landmarks));
-
-    int dim = boost::size(*std::begin(points));
-    typedef std::pair<double, int> dist_i;
-    typedef bool (*comp)(dist_i, dist_i);
-    knn = KNearestNeighbours(nbP);
-    for (int points_i = 0; points_i < nbP; points_i++) {
-      std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) {
-          return j1.first > j2.first;
-        });
-      typename LandmarkContainer::const_iterator landmarks_it;
-      int landmarks_i = 0;
-      for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end();
-           ++landmarks_it, landmarks_i++) {
-        dist_i dist = std::make_pair(euclidean_distance(points[points_i], *landmarks_it), landmarks_i);
-        l_heap.push(dist);
-      }
-      for (int i = 0; i < dim + 1; i++) {
-        dist_i dist = l_heap.top();
-        knn[points_i].push_back(dist.second);
-        l_heap.pop();
-      }
-    }
-  }
-
-}  // namespace witness_complex
-
-}  // namespace Gudhi
-
-#endif  // CONSTRUCT_CLOSEST_LANDMARK_TABLE_H_
diff --git a/src/Witness_complex/include/gudhi/Dim_list_iterator.h b/src/Witness_complex/include/gudhi/Dim_list_iterator.h
deleted file mode 100644
index 3f23e8c9..00000000
--- a/src/Witness_complex/include/gudhi/Dim_list_iterator.h
+++ /dev/null
@@ -1,155 +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) 2015  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/>.
- */
-
-#ifndef DIM_LISTS_ITERATOR_H_
-#define DIM_LISTS_ITERATOR_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /** \addtogroup simplex_tree
-   *  Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$:
-   *  \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$
-   *  and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that
-   *  w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. 
-   */
-template< class Simplex_handle >
-class Dim_lists_iterator {
-
-private:
-
-  typedef typename std::list<Simplex_handle>::iterator List_iterator;
-  typedef typename std::vector<std::list<Simplex_handle>>::reverse_iterator Vector_iterator;
-  typedef typename Gudhi::witness_complex::Dim_lists_iterator<Simplex_handle> Iterator;
-  
-  
-  List_iterator sh_;
-  Vector_iterator curr_list_;
-  typename std::vector<std::list<Simplex_handle>>& table_;
-
-public:
-
-  Dim_lists_iterator(List_iterator sh,
-                     Vector_iterator curr_list,
-                     typename std::vector<std::list<Simplex_handle>>& table)
-    : sh_(sh), curr_list_(curr_list), table_(table)
-  {
-  }
-
-  Simplex_handle operator*()
-  {
-    return *sh_;
-  }
-  
-  Iterator operator++()
-  {
-    increment();
-    return (*this);
-  }
-
-  Iterator operator++(int)
-  {
-    Iterator prev_it(sh_, curr_list_, table_);
-    increment();
-    return prev_it;
-  }
-
-  Iterator dim_begin()
-  {
-    return Iterator(curr_list_->begin(), curr_list_, table_);
-  }  
-
-  Iterator dim_end()
-  {
-    return Iterator(curr_list_->end(), curr_list_, table_);
-  }
-  
-  Iterator dimp1_begin()
-  {
-    return Iterator((curr_list_-1)->begin(), curr_list_-1, table_); 
-  }
-
-  Iterator dimp1_end()
-  {
-    return Iterator((curr_list_-1)->end(), curr_list_-1, table_); 
-  }
-  
-  bool operator==(const Iterator& it2) const
-  {
-    return (sh_ == it2.sh_);
-  }
-
-  bool operator!=(const Iterator& it2) const
-  {
-    return (sh_ != it2.sh_);
-  }
-
-  void remove_incr()
-  {
-    
-  }
-  
-private:
-  
-  void increment()
-  {
-    if (++sh_ == curr_list_->end())
-      if (++curr_list_ != table_.rend())
-        sh_ = curr_list_->begin();
-    // The iterator of the end of the table is the end of the last list
-  }
-
-  
-}; //class Dim_lists_iterator
-
-} // namespace witness_complex
-  
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Dim_lists.h b/src/Witness_complex/include/gudhi/Dim_lists.h
deleted file mode 100644
index 1d1b03c3..00000000
--- a/src/Witness_complex/include/gudhi/Dim_lists.h
+++ /dev/null
@@ -1,195 +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) 2015  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/>.
- */
-
-#ifndef DIM_LISTS_H_
-#define DIM_LISTS_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include <gudhi/Dim_list_iterator.h>
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /** \addtogroup simplex_tree
-   *  Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$:
-   *  \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$
-   *  and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that
-   *  w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. 
-   */
-template< class Simplicial_complex >
-class Dim_lists {
-
-private:
-
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Gudhi::witness_complex::Dim_lists_iterator<Simplex_handle> Iterator;
-  
-  std::vector<std::list<Simplex_handle>> table_;
-  Simplicial_complex& sc_;
-
-public:
-
-  Dim_lists(Simplicial_complex & sc, int dim, double alpha_max = 100)
-    : sc_(sc)
-  {
-    table_ = std::vector<std::list<Simplex_handle>>(dim+1);
-    for (auto sh: sc.filtration_simplex_range()) {
-      if (sc_.filtration(sh) < alpha_max)
-        table_[sc.dimension(sh)].push_front(sh);
-    }
-    auto t_it = table_.rbegin();
-    while (t_it->empty()) {
-      t_it++;
-      table_.pop_back();
-    }
-  }
-
-  Iterator begin()
-  {
-    return Iterator(table_.rbegin()->begin(), table_.rbegin(), table_);
-  }
-
-  Iterator end()
-  {
-    return Iterator(table_[0].end(), table_.rend(), table_);
-  }
-  
-  unsigned size()
-  {
-    unsigned curr_size = 0;
-    for (auto l: table_)
-      curr_size += l.size();
-    return curr_size;
-  }
-  
-  bool is_face(Simplex_handle face, Simplex_handle coface)
-  {
-    // vertex range is sorted in decreasing order
-    auto fvr = sc_.simplex_vertex_range(face);
-    auto cfvr = sc_.simplex_vertex_range(coface);
-    auto fv_it = fvr.begin();
-    auto cfv_it = cfvr.begin();
-    while (fv_it != fvr.end() && cfv_it != cfvr.end()) {
-      if (*fv_it < *cfv_it)
-        ++cfv_it;
-      else if (*fv_it == *cfv_it) {
-        ++cfv_it;
-        ++fv_it;
-      }
-      else
-        return false;
-      
-    }
-    return (fv_it == fvr.end());
-  }
-
-  void output_simplices() {
-    std::cout << "Size of vector: " << size() << std::endl;
-    for (auto line_it = table_.rbegin(); line_it != table_.rend(); ++line_it)
-      for (auto sh: *line_it) {
-        std::cout << sc_.dimension(sh) << " ";
-        for (auto v : sc_.simplex_vertex_range(sh))
-          std::cout << v << " ";
-        std::cout << sc_.filtration(sh) << "\n";
-      }
-  }
-  
-  void collapse()
-  {
-    auto coface_list_it = table_.rbegin();
-    auto face_list_it = table_.rbegin()+1;
-    for ( ;
-          face_list_it != table_.rend();
-          ++face_list_it) {
-      auto face_it = face_list_it->begin();
-      while (face_it != face_list_it->end() && sc_.filtration(*face_it) != 0) {
-        int coface_count = 0;
-        auto reduced_coface = coface_list_it->begin();
-        for (auto coface_it = coface_list_it->begin(); coface_it != coface_list_it->end() && sc_.filtration(*coface_it) != 0; ++coface_it)
-          if (is_face(*face_it, *coface_it)) {
-            coface_count++;
-            if (coface_count == 1)
-              reduced_coface = coface_it;
-            else
-              break;
-          }
-        if (coface_count == 1) {
-          /*
-          std::cout << "Erase ( ";
-          for (auto v: sc_.simplex_vertex_range(*reduced_coface))
-            std::cout << v << " ";        
-          */
-          coface_list_it->erase(reduced_coface);
-          /*
-          std::cout << ") and then ( ";
-          for (auto v: sc_.simplex_vertex_range(*face_it))
-            std::cout << v << " ";
-          std::cout << ")\n";
-          */
-          face_list_it->erase(face_it);
-          face_it = face_list_it->begin();
-        } 
-        else
-          face_it++;
-      }
-      if ((coface_list_it++)->empty())
-        table_.pop_back();
-    }
-  }
-  
-  bool is_pseudomanifold()
-  {
-
-    return true;
-  }
-  
-}; //class Dim_lists
-
-} // namespace witness_complex
-  
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Good_links.h b/src/Witness_complex/include/gudhi/Good_links.h
deleted file mode 100644
index c5e993e5..00000000
--- a/src/Witness_complex/include/gudhi/Good_links.h
+++ /dev/null
@@ -1,449 +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) 2015  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/>.
- */
-
-#ifndef GOOD_LINKS_H_
-#define GOOD_LINKS_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include <gudhi/Dim_list_iterator.h>
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-  
-template < typename Simplicial_complex >
-class Good_links {
-  
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Simplicial_complex::Vertex_handle Vertex_handle;
-  typedef std::vector<Vertex_handle> Vertex_vector;
-  
-  // graph is an adjacency list
-  typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Adj_graph;
-  // map that gives to a certain simplex its node in graph and its dimension
-  //typedef std::pair<boost::vecS,Color> Reference;
-  typedef boost::graph_traits<Adj_graph>::vertex_descriptor Vertex_t;
-  typedef boost::graph_traits<Adj_graph>::edge_descriptor Edge_t;
-  typedef boost::graph_traits<Adj_graph>::adjacency_iterator Adj_it;
-  typedef std::pair<Adj_it, Adj_it> Out_edge_it;
-  
-  typedef boost::container::flat_map<Simplex_handle, Vertex_t> Graph_map;
-  typedef boost::container::flat_map<Vertex_t, Simplex_handle> Inv_graph_map;
-
-public:
-  Good_links(Simplicial_complex& sc): sc_(sc)
-  {
-    int dim = 0;
-    for (auto sh: sc_.complex_simplex_range())
-      if (sc_.dimension(sh) > dim)
-        dim = sc_.dimension(sh);
-    dimension = dim;
-    count_good = Vertex_vector(dim);
-    count_bad = Vertex_vector(dim);
-  }
-  
-private:
-  
-  Simplicial_complex& sc_;
-  unsigned dimension;
-  std::vector<int> count_good;
-  std::vector<int> count_bad;
-  
-  void add_vertices_to_link_graph(Vertex_vector& star_vertices,
-                                  typename Vertex_vector::iterator curr_v,
-                                  Adj_graph& adj_graph,
-                                  Graph_map& d_map,
-                                  Graph_map& f_map,
-                                  Vertex_vector& curr_simplex,
-                                  int curr_d,
-                                  int link_dimension)
-  {
-    Simplex_handle sh;
-    Vertex_t vert;
-    typename Vertex_vector::iterator it;
-    //std::pair<typename Graph_map::iterator,bool> resPair;
-    //typename Graph_map::iterator resPair;
-    //Add vertices
-    //std::cout << "Entered add vertices\n";
-    for (it = curr_v; it != star_vertices.end(); ++it) {
-      curr_simplex.push_back(*it);               //push next vertex in question
-      curr_simplex.push_back(star_vertices[0]);  //push the center of the star
-      /*
-        std::cout << "Searching for ";
-        print_vector(curr_simplex);
-        std::cout << " curr_dim " << curr_d << " d " << dimension << "";
-      */
-      Vertex_vector curr_simplex_copy(curr_simplex);
-      sh = sc_.find(curr_simplex_copy);                   //a simplex of the star
-      curr_simplex.pop_back();                   //pop the center of the star
-      curr_simplex_copy = Vertex_vector(curr_simplex);
-      if (sh != sc_.null_simplex()) {
-        //std::cout << " added\n";
-        if (curr_d == link_dimension) {
-          sh = sc_.find(curr_simplex_copy);               //a simplex of the link
-          assert(sh != sc_.null_simplex());          //ASSERT!
-          vert = boost::add_vertex(adj_graph);
-          d_map.emplace(sh,vert);
-        }
-        else {        
-          if (curr_d == link_dimension-1) {
-            sh = sc_.find(curr_simplex_copy);               //a simplex of the link
-            assert(sh != sc_.null_simplex());
-            vert = boost::add_vertex(adj_graph);
-            f_map.emplace(sh,vert);
-          }
-          
-          //delete (&curr_simplex_copy); //Just so it doesn't stack
-          add_vertices_to_link_graph(star_vertices,
-                                     it+1,
-                                     adj_graph,
-                                     d_map,
-                                     f_map,
-                                     curr_simplex,
-                                     curr_d+1, link_dimension);
-        }
-      }
-      /*
-        else
-        std::cout << "\n";
-      */
-      curr_simplex.pop_back();                           //pop the vertex in question
-    }
-  }
-
-  void add_edges_to_link_graph(Adj_graph& adj_graph,
-                               Graph_map& d_map,
-                               Graph_map& f_map)
-    {
-      Simplex_handle sh;
-      // Add edges
-      //std::cout << "Entered add edges:\n";
-      typename Graph_map::iterator map_it;
-      for (auto d_map_pair : d_map) {
-        //std::cout << "*";
-        sh = d_map_pair.first;
-        Vertex_t d_vert = d_map_pair.second;
-        for (auto facet_sh : sc_.boundary_simplex_range(sh))
-          //for (auto f_map_it : f_map)
-          {
-            //std::cout << "'"; 
-            map_it = f_map.find(facet_sh);
-            //We must have all the facets in the graph at this point
-            assert(map_it != f_map.end());
-            Vertex_t f_vert = map_it->second;
-            //std::cout << "Added edge " << sh->first << "-" << map_it->first->first << "\n";
-            boost::add_edge(d_vert,f_vert,adj_graph);
-          }
-      }
-    }
-
-
-  
-  /* \brief Verifies if the simplices formed by vertices given by link_vertices 
-   * form a pseudomanifold.
-   * The idea is to make a bipartite graph, where vertices are the d- and (d-1)-dimensional
-   * faces and edges represent adjacency between them.
-   */
-  bool link_is_pseudomanifold(Vertex_vector& star_vertices,
-                              int dimension)
-  {
-    Adj_graph adj_graph;
-    Graph_map d_map, f_map;
-    // d_map = map for d-dimensional simplices
-    // f_map = map for its facets
-    Vertex_vector init_vector = {};
-    add_vertices_to_link_graph(star_vertices,
-                               star_vertices.begin()+1,
-                               adj_graph,
-                               d_map,
-                               f_map,
-                               init_vector,
-                               0, dimension);
-    //std::cout << "DMAP_SIZE: " << d_map.size() << "\n";
-    //std::cout << "FMAP_SIZE: " << f_map.size() << "\n";
-    add_edges_to_link_graph(adj_graph,
-                            d_map,
-                            f_map);
-    for (auto f_map_it : f_map) {
-      //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n";
-      if (boost::out_degree(f_map_it.second, adj_graph) != 2){
-        count_bad[dimension]++;
-        return false;
-      }
-    }
-    // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices
-    // What is left is to check the connexity
-    //std::vector<int> components(boost::num_vertices(adj_graph));
-    return true; //Forget the connexity
-    //return (boost::connected_components(adj_graph, &components[0]) == 1);
-  }
-
-  int star_dim(Vertex_vector& star_vertices,
-               typename Vertex_vector::iterator curr_v,
-               int curr_d,
-               Vertex_vector& curr_simplex,
-               typename std::vector< int >::iterator curr_dc)
-  {
-    //std::cout << "Entered star_dim for " << *(curr_v-1) << "\n";
-    Simplex_handle sh;
-    int final_d = curr_d;
-      typename Vertex_vector::iterator it;
-      typename Vertex_vector::iterator dc_it;
-      //std::cout << "Current vertex is " <<  
-      for (it = curr_v, dc_it = curr_dc; it != star_vertices.end(); ++it, ++dc_it)
-        {
-          curr_simplex.push_back(*it);
-          Vertex_vector curr_simplex_copy(curr_simplex);
-          /*
-          std::cout << "Searching for ";
-          print_vector(curr_simplex);
-          std::cout << " curr_dim " << curr_d << " final_dim " << final_d;
-          */
-          sh = sc_.find(curr_simplex_copy); //Need a copy because find sorts the vector and I want star center to be the first
-          if (sh != sc_.null_simplex())
-            {
-              //std::cout << " -> " << *it << "\n";
-              int d = star_dim(star_vertices,
-                               it+1,
-                               curr_d+1,
-                               curr_simplex,
-                               dc_it);
-              if (d >= final_d)
-                {
-                  final_d = d;
-                  //std::cout << d << " ";
-                  //print_vector(curr_simplex);
-		  //std::cout << std::endl;
-                }
-              if (d >= *dc_it)
-                *dc_it = d;
-            }
-          /*
-          else
-            std::cout << "\n";
-          */
-          curr_simplex.pop_back();
-        }
-      return final_d;
-    }
-
-public:
-  
-  /** \brief Returns true if the link is good
-   */
-  bool has_good_link(Vertex_handle v)
-  {
-    typedef Vertex_vector typeVectorVertex;
-    typeVectorVertex star_vertices;
-    // Fill star_vertices
-    star_vertices.push_back(v);
-    for (auto u: sc_.complex_vertex_range())
-      {
-        typeVectorVertex edge = {u,v};
-        if (u != v && sc_.find(edge) != sc_.null_simplex())   
-          star_vertices.push_back(u);
-      }
-    // Find the dimension
-    typeVectorVertex init_simplex = {star_vertices[0]};
-    bool is_pure = true;
-    std::vector<int> dim_coface(star_vertices.size(), 1);
-    int d = star_dim(star_vertices,
-                     star_vertices.begin()+1,
-                     0,
-                     init_simplex,
-                     dim_coface.begin()+1) - 1; //link_dim = star_dim - 1
-    
-    assert(init_simplex.size() == 1);
-    // if (!is_pure)
-    //   std::cout << "Found an impure star around " << v << "\n";
-    for (int dc: dim_coface)
-      is_pure = (dc == dim_coface[0]);
-    /*
-      if (d == count_good.size())
-      {
-      std::cout << "Found a star of dimension " << (d+1) << " around " << v << "\nThe star is ";
-      print_vector(star_vertices); std::cout << std::endl;
-      }
-    */
-    //if (d == -1) count_bad[0]++;
-    bool b= (is_pure && link_is_pseudomanifold(star_vertices,d));
-    if (d != -1) {if (b) count_good[d]++; else count_bad[d]++;}
-    if (!is_pure) count_bad[0]++;
-    return (d != -1 && b && is_pure);
-  }
-
-private:
-  void add_max_simplices_to_graph(Vertex_vector& star_vertices,
-                                  typename Vertex_vector::iterator curr_v,
-                                  Adj_graph& adj_graph,
-                                  Graph_map& d_map,
-                                  Graph_map& f_map,
-                                  Inv_graph_map& inv_d_map,
-                                  Vertex_vector& curr_simplex,
-                                  int curr_d,
-                                  int link_dimension)
-  {
-    Simplex_handle sh;
-    Vertex_t vert;
-    typename Vertex_vector::iterator it;
-    //std::pair<typename Graph_map::iterator,bool> resPair;
-    //typename Graph_map::iterator resPair;
-    //Add vertices
-    //std::cout << "Entered add vertices\n";
-    for (it = curr_v; it != star_vertices.end(); ++it) {
-      curr_simplex.push_back(*it);               //push next vertex in question
-      //curr_simplex.push_back(star_vertices[0]);  //push the center of the star
-      /*
-        std::cout << "Searching for ";
-        print_vector(curr_simplex);
-        std::cout << " curr_dim " << curr_d << " d " << dimension << "";
-      */
-      Vertex_vector curr_simplex_copy(curr_simplex);
-      sh = sc_.find(curr_simplex_copy);                   //a simplex of the star
-      //curr_simplex.pop_back();                   //pop the center of the star
-      curr_simplex_copy = Vertex_vector(curr_simplex);
-      if (sh != sc_.null_simplex()) {
-        //std::cout << " added\n";
-        if (curr_d == link_dimension) {
-          sh = sc_.find(curr_simplex_copy);               //a simplex of the link
-          assert(sh != sc_.null_simplex()); //ASSERT!
-          vert = boost::add_vertex(adj_graph);
-          d_map.emplace(sh,vert);
-          inv_d_map.emplace(vert,sh);
-        }
-        else {      
-          if (curr_d == link_dimension-1) {
-              sh = sc_.find(curr_simplex_copy);               //a simplex of the link
-              assert(sh != sc_.null_simplex());
-              vert = boost::add_vertex(adj_graph);
-              f_map.emplace(sh,vert);
-          }        
-          //delete (&curr_simplex_copy); //Just so it doesn't stack
-          add_max_simplices_to_graph(star_vertices,
-                                     it+1,
-                                     adj_graph,
-                                     d_map,
-                                     f_map,
-                                     inv_d_map,
-                                     curr_simplex,
-                                     curr_d+1,
-                                     link_dimension);
-        }
-      }
-      /*
-        else
-        std::cout << "\n";
-      */
-      curr_simplex.pop_back();                           //pop the vertex in question
-    }
-  }
-
-public:
-  bool complex_is_pseudomanifold()
-  {
-    Adj_graph adj_graph;
-    Graph_map d_map, f_map;
-    // d_map = map for d-dimensional simplices
-    // f_map = map for its facets
-    Inv_graph_map inv_d_map;
-    Vertex_vector init_vector = {};
-    std::vector<int> star_vertices;
-    for (int v: sc_.complex_vertex_range())
-      star_vertices.push_back(v);
-    add_max_simplices_to_graph(star_vertices,
-                               star_vertices.begin(),
-                               adj_graph,
-                               d_map,
-                               f_map,
-                               inv_d_map,
-                               init_vector,
-                               0, dimension);
-    //std::cout << "DMAP_SIZE: " << d_map.size() << "\n";
-    //std::cout << "FMAP_SIZE: " << f_map.size() << "\n";
-    add_edges_to_link_graph(adj_graph,
-                            d_map,
-                            f_map);
-    for (auto f_map_it : f_map) {
-      //std::cout << "Degree of " << f_map_it.first->first << " is " << boost::out_degree(f_map_it.second, adj_graph) << "\n";
-      if (boost::out_degree(f_map_it.second, adj_graph) != 2) {
-        // if (boost::out_degree(f_map_it.second, adj_graph) >= 3) {		 
-        //   // std::cout << "This simplex has 3+ cofaces: ";
-        //   // for(auto v : sc_.simplex_vertex_range(f_map_it.first))
-        //   //   std::cout << v << " ";
-        //   // std::cout << std::endl;
-        //   Adj_it ai, ai_end; 
-        //   for (std::tie(ai, ai_end) = boost::adjacent_vertices(f_map_it.second, adj_graph); ai != ai_end; ++ai) {
-        //     auto it = inv_d_map.find(*ai);
-        //     assert (it != inv_d_map.end());
-        //     typename Simplicial_complex::Simplex_handle sh = it->second;
-        //     for(auto v : sc_.simplex_vertex_range(sh))
-        //       std::cout << v << " ";
-        //     std::cout << std::endl;
-        //   }
-        // }
-        count_bad[dimension]++;
-        return false;
-      }
-    }
-    // At this point I know that all (d-1)-simplices are adjacent to exactly 2 d-simplices
-    // What is left is to check the connexity
-    //std::vector<int> components(boost::num_vertices(adj_graph));
-    return true; //Forget the connexity
-    //return (boost::connected_components(adj_graph, &components[0]) == 1);
-  }
-
-  int number_good_links(int dim)
-  {
-    return count_good[dim];
-  }
-  
-  int number_bad_links(int dim)
-  {
-    return count_bad[dim];
-  }
-  
-};
-    
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h
deleted file mode 100644
index df93155b..00000000
--- a/src/Witness_complex/include/gudhi/Landmark_choice_by_furthest_point.h
+++ /dev/null
@@ -1,105 +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) 2015  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/>.
- */
-
-#ifndef LANDMARK_CHOICE_BY_FURTHEST_POINT_H_
-#define LANDMARK_CHOICE_BY_FURTHEST_POINT_H_
-
-#include <boost/range/size.hpp>
-
-#include <limits>  // for numeric_limits<>
-#include <iterator>
-#include <algorithm>  // for sort
-#include <vector>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  typedef std::vector<int> typeVectorVertex;
-
-  /** 
-   *  \ingroup witness_complex
-   *  \brief Landmark choice strategy by iteratively adding the furthest witness from the
-   *  current landmark set as the new landmark. 
-   *  \details It chooses nbL landmarks from a random access range `points` and
-   *  writes {witness}*{closest landmarks} matrix in `knn`.
-   *
-   *  The type KNearestNeighbors can be seen as 
-   *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-   *  Witness_range and Closest_landmark_range are random access ranges 
-   *  
-   */
-
-  template <typename KNearestNeighbours,
-  typename Point_random_access_range>
-  void landmark_choice_by_furthest_point(Point_random_access_range const &points,
-                                         int nbL,
-                                         KNearestNeighbours &knn) {
-    int nb_points = boost::size(points);
-    assert(nb_points >= nbL);
-    // distance matrix witness x landmarks
-    std::vector<std::vector<double>> wit_land_dist(nb_points, std::vector<double>());
-    // landmark list
-    typeVectorVertex chosen_landmarks;
-
-    knn = KNearestNeighbours(nb_points, std::vector<int>());
-    int current_number_of_landmarks = 0;  // counter for landmarks
-    double curr_max_dist = 0;  // used for defining the furhest point from L
-    const double infty = std::numeric_limits<double>::infinity();  // infinity (see next entry)
-    std::vector< double > dist_to_L(nb_points, infty);  // vector of current distances to L from points
-
-    // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety
-    // or better yet std::uniform_int_distribution
-    int rand_int = rand() % nb_points;
-    int curr_max_w = rand_int;  // For testing purposes a pseudo-random number is used here
-
-    for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) {
-      // curr_max_w at this point is the next landmark
-      chosen_landmarks.push_back(curr_max_w);
-      unsigned i = 0;
-      for (auto& p : points) {
-        double curr_dist = euclidean_distance(p, *(std::begin(points) + chosen_landmarks[current_number_of_landmarks]));
-        wit_land_dist[i].push_back(curr_dist);
-        knn[i].push_back(current_number_of_landmarks);
-        if (curr_dist < dist_to_L[i])
-          dist_to_L[i] = curr_dist;
-        ++i;
-      }
-      curr_max_dist = 0;
-      for (i = 0; i < dist_to_L.size(); i++)
-        if (dist_to_L[i] > curr_max_dist) {
-          curr_max_dist = dist_to_L[i];
-          curr_max_w = i;
-        }
-    }
-    for (int i = 0; i < nb_points; ++i)
-      std::sort(std::begin(knn[i]),
-                std::end(knn[i]),
-                [&wit_land_dist, i](int a, int b) {
-                  return wit_land_dist[i][a] < wit_land_dist[i][b]; });
-  }
-
-}  // namespace witness_complex
-
-}  // namespace Gudhi
-
-#endif  // LANDMARK_CHOICE_BY_FURTHEST_POINT_H_
diff --git a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h b/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h
deleted file mode 100644
index ebf6aad1..00000000
--- a/src/Witness_complex/include/gudhi/Landmark_choice_by_random_point.h
+++ /dev/null
@@ -1,96 +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) 2015  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/>.
- */
-
-#ifndef LANDMARK_CHOICE_BY_RANDOM_POINT_H_
-#define LANDMARK_CHOICE_BY_RANDOM_POINT_H_
-
-#include <boost/range/size.hpp>
-
-#include <queue>  // for priority_queue<>
-#include <utility>  // for pair<>
-#include <iterator>
-#include <vector>
-#include <set>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /**
-   *  \ingroup witness_complex
-   * \brief Landmark choice strategy by taking random vertices for landmarks.
-   *  \details It chooses nbL distinct landmarks from a random access range `points`
-   *  and outputs a matrix {witness}*{closest landmarks} in knn.
-   *
-   *  The type KNearestNeighbors can be seen as 
-   *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-   *  Witness_range and Closest_landmark_range are random access ranges and
-   *  Vertex_handle is the label type of a vertex in a simplicial complex.
-   *  Closest_landmark_range needs to have push_back operation.
-   */
-
-  template <typename KNearestNeighbours,
-            typename Point_random_access_range>
-  void landmark_choice_by_random_point(Point_random_access_range const &points,
-                                       int nbL,
-                                       KNearestNeighbours &knn) {
-    int nbP = boost::size(points);
-    assert(nbP >= nbL);
-    std::set<int> landmarks;
-    int current_number_of_landmarks = 0;  // counter for landmarks
-
-    // TODO(SK) Consider using rand_r(...) instead of rand(...) for improved thread safety
-    int chosen_landmark = rand() % nbP;
-    for (current_number_of_landmarks = 0; current_number_of_landmarks != nbL; current_number_of_landmarks++) {
-      while (landmarks.find(chosen_landmark) != landmarks.end())
-        chosen_landmark = rand() % nbP;
-      landmarks.insert(chosen_landmark);
-    }
-
-    int dim = boost::size(*std::begin(points));
-    typedef std::pair<double, int> dist_i;
-    typedef bool (*comp)(dist_i, dist_i);
-    knn = KNearestNeighbours(nbP);
-    for (int points_i = 0; points_i < nbP; points_i++) {
-      std::priority_queue<dist_i, std::vector<dist_i>, comp> l_heap([](dist_i j1, dist_i j2) {
-          return j1.first > j2.first;
-        });
-      std::set<int>::iterator landmarks_it;
-      int landmarks_i = 0;
-      for (landmarks_it = landmarks.begin(), landmarks_i = 0; landmarks_it != landmarks.end();
-           ++landmarks_it, landmarks_i++) {
-        dist_i dist = std::make_pair(euclidean_distance(points[points_i], points[*landmarks_it]), landmarks_i);
-        l_heap.push(dist);
-      }
-      for (int i = 0; i < dim + 1; i++) {
-        dist_i dist = l_heap.top();
-        knn[points_i].push_back(dist.second);
-        l_heap.pop();
-      }
-    }
-  }
-
-}  // namespace witness_complex
-
-}  // namespace Gudhi
-
-#endif  // LANDMARK_CHOICE_BY_RANDOM_POINT_H_
diff --git a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h
deleted file mode 100644
index 2971149a..00000000
--- a/src/Witness_complex/include/gudhi/Relaxed_witness_complex.h
+++ /dev/null
@@ -1,379 +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) 2015  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/>.
- */
-
-#ifndef RELAXED_WITNESS_COMPLEX_H_
-#define RELAXED_WITNESS_COMPLEX_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-// Needed for nearest neighbours
-#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /** \addtogroup simplex_tree
-   *  Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$:
-   *  \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$
-   *  and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that
-   *  w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. 
-   */
-template< class Simplicial_complex >
-class Relaxed_witness_complex {
-
-private:    
-  struct Active_witness {
-    int witness_id;
-    int landmark_id;
-    
-    Active_witness(int witness_id_, int landmark_id_)
-      : witness_id(witness_id_),
-        landmark_id(landmark_id_) { }
-    };
-
-private:
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Simplicial_complex::Vertex_handle Vertex_handle;
-  typedef typename Simplicial_complex::Filtration_value FT;
-
-  typedef std::vector< double > Point_t;
-  typedef std::vector< Point_t > Point_Vector;
-
-  // typedef typename Simplicial_complex::Filtration_value Filtration_value;
-  typedef std::vector< Vertex_handle > typeVectorVertex;
-  typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex;
-  typedef std::pair< Simplex_handle, bool > typePairSimplexBool;
-
-  typedef int Witness_id;
-  typedef int Landmark_id;
-  typedef std::list< Vertex_handle > ActiveWitnessList;
-
- private:
-  int nbL;  // Number of landmarks
-  Simplicial_complex& sc;  // Simplicial complex
-  
-  public:
-    /** @name Constructor
-     */
-
-    //@{
-  
-    /**
-     *  \brief Iterative construction of the relaxed witness complex.
-     *  \details The witness complex is written in sc_ basing on a matrix knn
-     *  of k nearest neighbours of the form {witnesses}x{landmarks} and 
-     *  and a matrix distances of distances to these landmarks from witnesses.
-     *  The parameter alpha defines relaxation and
-     *  limD defines the 
-     *
-     *  The line lengths in one given matrix can differ, 
-     *  however both matrices have the same corresponding line lengths.
-     *
-     *  The type KNearestNeighbors can be seen as 
-     *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-     *  Witness_range and Closest_landmark_range are random access ranges.
-     *
-     *  Constructor takes into account at most (dim+1) 
-     *  first landmarks from each landmark range to construct simplices.
-     *
-     *  Landmarks are supposed to be in [0,nbL_-1]
-     */
-    
-    template< typename KNearestNeighbours >
-    Relaxed_witness_complex(std::vector< std::vector<double> > const & distances,
-                            KNearestNeighbours const & knn,
-                            Simplicial_complex & sc_, 
-                            int nbL_,
-                            double alpha,
-                            int limD) : nbL(nbL_), sc(sc_) {
-      int nbW = knn.size();
-      typeVectorVertex vv;
-      //int counter = 0;
-      /* The list of still useful witnesses
-       * it will diminuish in the course of iterations
-       */
-      ActiveWitnessList active_w;// = new ActiveWitnessList();
-      for (int i = 0; i != nbL; ++i) {
-        // initial fill of 0-dimensional simplices
-        // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
-        //counter++;
-        vv = {i};
-        sc.insert_simplex(vv, Filtration_value(0.0));
-        /* TODO Error if not inserted : normally no need here though*/
-      }
-      int k = 1; /* current dimension in iterative construction */
-      for (int i=0; i != nbW; ++i)
-        active_w.push_back(i);
-      while (!active_w.empty() && k <= limD && k < nbL )
-        {
-          typename ActiveWitnessList::iterator aw_it = active_w.begin();
-          while (aw_it != active_w.end())
-            {
-              std::vector<int> simplex;
-              bool ok = add_all_faces_of_dimension(k,
-                                                   alpha,
-                                                   std::numeric_limits<double>::infinity(),
-                                                   distances[*aw_it].begin(),
-                                                   knn[*aw_it].begin(),
-                                                   simplex,
-                                                   knn[*aw_it].end());
-              if (!ok)
-                active_w.erase(aw_it++); //First increase the iterator and then erase the previous element
-              else
-                aw_it++;
-            }
-          k++;
-        }
-      sc.set_dimension(limD);
-    }
-
-  //@}
-  
-private:
-  /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness
-   * Iterator is needed to know until how far we can take landmarks to form simplexes
-   * simplex is the prefix of the simplexes to insert
-   * The output value indicates if the witness rests active or not
-   */
-  bool add_all_faces_of_dimension(int dim,
-                                  double alpha,
-                                  double norelax_dist,
-                                  std::vector<double>::const_iterator curr_d,
-                                  std::vector<int>::const_iterator curr_l,
-                                  std::vector<int>& simplex,
-                                  std::vector<int>::const_iterator end)
-  {
-    if (curr_l == end)
-      return false;
-    bool will_be_active = false;
-    std::vector<int>::const_iterator l_it = curr_l;
-    std::vector<double>::const_iterator d_it = curr_d;
-    if (dim > 0)
-      for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) {
-        simplex.push_back(*l_it);
-        if (sc.find(simplex) != sc.null_simplex())
-          will_be_active = add_all_faces_of_dimension(dim-1,
-                                                      alpha,
-                                                      norelax_dist,
-                                                      d_it+1,
-                                                      l_it+1,
-                                                      simplex,
-                                                      end) || will_be_active;
-        simplex.pop_back();
-        // If norelax_dist is infinity, change to first omitted distance
-        if (*d_it <= norelax_dist)
-          norelax_dist = *d_it;
-        will_be_active = add_all_faces_of_dimension(dim-1,
-                                                    alpha,
-                                                    norelax_dist,
-                                                    d_it+1,
-                                                    l_it+1,
-                                                    simplex,
-                                                    end) || will_be_active;
-      } 
-    else if (dim == 0)
-      for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) {
-        simplex.push_back(*l_it);
-        double filtration_value = 0;
-        // if norelax_dist is infinite, relaxation is 0.
-        if (*d_it > norelax_dist) 
-          filtration_value = *d_it - norelax_dist; 
-        if (all_faces_in(simplex, &filtration_value)) {
-          will_be_active = true;
-          sc.insert_simplex(simplex, filtration_value);
-        }
-        simplex.pop_back();
-        // If norelax_dist is infinity, change to first omitted distance
-        if (*d_it < norelax_dist)
-          norelax_dist = *d_it;
-      } 
-    return will_be_active;
-  }
-  
-  /** \brief Check if the facets of the k-dimensional simplex witnessed 
-   *  by witness witness_id are already in the complex.
-   *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
-   */
-  bool all_faces_in(std::vector<int>& simplex, double* filtration_value)
-  {
-    std::vector< int > facet;
-    for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it)
-      {
-        facet.clear();
-        for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it)
-          if (it != not_it)
-            facet.push_back(*it);
-        Simplex_handle facet_sh = sc.find(facet);
-        if (facet_sh == sc.null_simplex())
-          return false;
-        else if (sc.filtration(facet_sh) > *filtration_value)
-          *filtration_value = sc.filtration(facet_sh);
-      }
-    return true;
-  }
-
-  bool is_face(Simplex_handle face, Simplex_handle coface)
-  {
-    // vertex range is sorted in decreasing order
-    auto fvr = sc.simplex_vertex_range(face);
-    auto cfvr = sc.simplex_vertex_range(coface);
-    auto fv_it = fvr.begin();
-    auto cfv_it = cfvr.begin();
-    while (fv_it != fvr.end() && cfv_it != cfvr.end()) {
-      if (*fv_it < *cfv_it)
-        ++cfv_it;
-      else if (*fv_it == *cfv_it) {
-        ++cfv_it;
-        ++fv_it;
-      }
-      else
-        return false;
-      
-    }
-    return (fv_it == fvr.end());
-  }
-
-  // void erase_simplex(Simplex_handle sh)
-  // {
-  //   auto siblings = sc.self_siblings(sh);
-  //   auto oncles = siblings->oncles();
-  //   int prev_vertex = siblings->parent();
-  //   siblings->members().erase(sh->first);
-  //   if (siblings->members().empty()) {
-  //     typename typedef Simplicial_complex::Siblings Siblings;
-  //     oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex));
-  //     assert(!sc.has_children(oncles->members().find(prev_vertex)));
-  //     //delete siblings;
-  //   }
-
-  // }
-  
-  void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh)
-  {
-    erase_simplex(coface_sh);
-    erase_simplex(face_sh);
-  }
-
-public:
-  // void collapse(std::vector<Simplex_handle>& simplices)
-  // {
-  //   // Get a vector of simplex handles ordered by filtration value
-  //   std::cout << sc << std::endl;
-  //   //std::vector<Simplex_handle> simplices;
-  //   for (Simplex_handle sh: sc.filtration_simplex_range())
-  //     simplices.push_back(sh);
-  //   // std::sort(simplices.begin(),
-  //   //           simplices.end(),
-  //   //           [&](Simplex_handle sh1, Simplex_handle sh2)
-  //   //           { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2);
-  //   //             return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); });
-  //   // Double iteration
-  //   auto face_it = simplices.rbegin();
-  //   while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) {
-  //     int coface_count = 0;
-  //     auto reduced_coface = simplices.rbegin();
-  //     for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it)
-  //       if (face_it != coface_it && is_face(*face_it, *coface_it)) {
-  //         coface_count++;
-  //         if (coface_count == 1)
-  //           reduced_coface = coface_it;
-  //         else
-  //           break;
-  //       }
-  //     if (coface_count == 1) {
-  //       std::cout << "Erase ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base())))
-  //         std::cout << v << " ";
-        
-  //       simplices.erase(--(reduced_coface.base()));
-  //       //elementary_collapse(*face_it, *reduced_coface);
-  //       std::cout << ") and then ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--face_it.base())))
-  //         std::cout << v << " ";
-  //       std::cout << ")\n";
-  //       simplices.erase(--((face_it++).base()));
-  //       //face_it = simplices.rbegin();
-  //       //std::cout << "Size of vector: " << simplices.size() << "\n";
-  //     } 
-  //     else
-  //       face_it++;
-  //   }
-  //   sc.initialize_filtration();
-  //   //std::cout << sc << std::endl;
-  // }
-
-  template <class Dim_lists>
-  void collapse(Dim_lists& dim_lists)
-  {
-    dim_lists.collapse();
-  }
-  
-private:
-  
-  /** Collapse recursively boundary faces of the given simplex
-   *  if its filtration is bigger than alpha_lim.  
-   */
-  void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim)
-  {
-    for (Simplex_handle face_it : sc.boundary_simplex_range()) {
-      
-    }
-      
-  }
-  
-}; //class Relaxed_witness_complex
-
-} // namespace witness_complex
-  
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Strange_witness_complex.h b/src/Witness_complex/include/gudhi/Strange_witness_complex.h
deleted file mode 100644
index c1c2912b..00000000
--- a/src/Witness_complex/include/gudhi/Strange_witness_complex.h
+++ /dev/null
@@ -1,232 +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) 2015  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/>.
- */
-
-#ifndef WITNESS_COMPLEX_H_
-#define WITNESS_COMPLEX_H_
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-
-#include <gudhi/distance_functions.h>
-
-#include <algorithm>
-#include <utility>
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <ctime>
-#include <iostream>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-/** 
-    \class Witness_complex
-    \brief Constructs the witness complex for the given set of witnesses and landmarks.
-    \ingroup witness_complex
- */
-template< class Simplicial_complex>
-class Strange_witness_complex {
- private:
-  struct Active_witness {
-    int witness_id;
-    int landmark_id;
-
-    Active_witness(int witness_id_, int landmark_id_)
-        : witness_id(witness_id_),
-        landmark_id(landmark_id_) { }
-  };
-
- private:
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Simplicial_complex::Vertex_handle Vertex_handle;
-
-  typedef std::vector< double > Point_t;
-  typedef std::vector< Point_t > Point_Vector;
-
-  // typedef typename Simplicial_complex::Filtration_value Filtration_value;
-  typedef std::vector< Vertex_handle > typeVectorVertex;
-  typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex;
-  typedef std::pair< Simplex_handle, bool > typePairSimplexBool;
-
-  typedef int Witness_id;
-  typedef int Landmark_id;
-  typedef std::list< Vertex_handle > ActiveWitnessList;
-
- private:
-  int nbL;  // Number of landmarks
-  Simplicial_complex& sc;  // Simplicial complex
-
- public:
-  /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-  /** @name Constructor
-   */
-
-  //@{
-
-  // Witness_range<Closest_landmark_range<Vertex_handle>>
-
-  /**
-   *  \brief Iterative construction of the witness complex.
-   *  \details The witness complex is written in sc_ basing on a matrix knn of
-   *  nearest neighbours of the form {witnesses}x{landmarks}.
-   *
-   *  The type KNearestNeighbors can be seen as 
-   *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-   *  Witness_range and Closest_landmark_range are random access ranges.
-   *
-   *  Constructor takes into account at most (dim+1) 
-   *  first landmarks from each landmark range to construct simplices.
-   *
-   *  Landmarks are supposed to be in [0,nbL_-1]
-   */
-  template< typename KNearestNeighbors >
-  Strange_witness_complex(KNearestNeighbors const & knn,
-                  Simplicial_complex & sc_,
-                  int nbL_,
-                  int dim) : nbL(nbL_), sc(sc_) {
-    // Construction of the active witness list
-    int nbW = knn.size();
-    typeVectorVertex vv;
-    int counter = 0;
-    /* The list of still useful witnesses
-     * it will diminuish in the course of iterations
-     */
-    for (int i = 0; i != nbL; ++i) {
-      // initial fill of 0-dimensional simplices
-      // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
-      counter++;
-      vv = {i};
-      sc.insert_simplex(vv);
-      // TODO(SK) Error if not inserted : normally no need here though
-    }
-    for (int i = 0; i != nbW; ++i) {
-      std::vector<int> simplex;
-      simplex.reserve(dim+1);
-      for (int j = 0; j <= dim; j++)
-        simplex.push_back(knn[i][j]);
-      sc.insert_simplex_and_subfaces(simplex);
-    }
-  }
-
-  //@}
-
- private:
-  /** \brief Check if the facets of the k-dimensional simplex witnessed 
-   *  by witness witness_id are already in the complex.
-   *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
-   */
-  template <typename KNearestNeighbors>
-  bool all_faces_in(KNearestNeighbors const &knn, int witness_id, int k) {
-    // std::cout << "All face in with the landmark " << inserted_vertex << std::endl;
-    std::vector< Vertex_handle > facet;
-    // Vertex_handle curr_vh = curr_sh->first;
-    // CHECK ALL THE FACETS
-    for (int i = 0; i != k + 1; ++i) {
-      facet = {};
-      for (int j = 0; j != k + 1; ++j) {
-        if (j != i) {
-          facet.push_back(knn[witness_id][j]);
-        }
-      }  // endfor
-      if (sc.find(facet) == sc.null_simplex())
-        return false;
-      // std::cout << "++++ finished loop safely\n";
-    }  // endfor
-    return true;
-  }
-
-  template <typename T>
-  static void print_vector(const std::vector<T>& v) {
-    std::cout << "[";
-    if (!v.empty()) {
-      std::cout << *(v.begin());
-      for (auto it = v.begin() + 1; it != v.end(); ++it) {
-        std::cout << ",";
-        std::cout << *it;
-      }
-    }
-    std::cout << "]";
-  }
-
- public:
-  /**
-   *  \brief Verification if every simplex in the complex is witnessed by witnesses in knn.
-   *  \param print_output =true will print the witnesses for each simplex
-   *  \remark Added for debugging purposes.
-   */
-  template< class KNearestNeighbors >
-  bool is_witness_complex(KNearestNeighbors const & knn, bool print_output) {
-    // bool final_result = true;
-    for (Simplex_handle sh : sc.complex_simplex_range()) {
-      bool is_witnessed = false;
-      typeVectorVertex simplex;
-      int nbV = 0;  // number of verticed in the simplex
-      for (int v : sc.simplex_vertex_range(sh))
-        simplex.push_back(v);
-      nbV = simplex.size();
-      for (typeVectorVertex w : knn) {
-        bool has_vertices = true;
-        for (int v : simplex)
-          if (std::find(w.begin(), w.begin() + nbV, v) == w.begin() + nbV) {
-            has_vertices = false;
-            // break;
-          }
-        if (has_vertices) {
-          is_witnessed = true;
-          if (print_output) {
-            std::cout << "The simplex ";
-            print_vector(simplex);
-            std::cout << " is witnessed by the witness ";
-            print_vector(w);
-            std::cout << std::endl;
-          }
-          break;
-        }
-      }
-      if (!is_witnessed) {
-        if (print_output) {
-          std::cout << "The following simplex is not witnessed ";
-          print_vector(simplex);
-          std::cout << std::endl;
-        }
-        assert(is_witnessed);
-        return false;
-      }
-    }
-    return true;
-  }
-};
-
-}  // namespace witness_complex
-
-}  // namespace Gudhi
-
-#endif  // WITNESS_COMPLEX_H_
diff --git a/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h
deleted file mode 100644
index ee863a42..00000000
--- a/src/Witness_complex/include/gudhi/Strong_relaxed_witness_complex.h
+++ /dev/null
@@ -1,337 +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) 2015  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/>.
- */
-
-#ifndef STRONG_RELAXED_WITNESS_COMPLEX_H_
-#define STRONG_RELAXED_WITNESS_COMPLEX_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-// Needed for nearest neighbours
-#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /** \addtogroup simplex_tree
-   *  Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$:
-   *  \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$
-   *  and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that
-   *  w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. 
-   */
-template< class Simplicial_complex >
-class Strong_relaxed_witness_complex {
-
-private:    
-  struct Active_witness {
-    int witness_id;
-    int landmark_id;
-    
-    Active_witness(int witness_id_, int landmark_id_)
-      : witness_id(witness_id_),
-        landmark_id(landmark_id_) { }
-    };
-
-private:
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Simplicial_complex::Vertex_handle Vertex_handle;
-  typedef typename Simplicial_complex::Filtration_value FT;
-
-  typedef std::vector< double > Point_t;
-  typedef std::vector< Point_t > Point_Vector;
-
-  // typedef typename Simplicial_complex::Filtration_value Filtration_value;
-  typedef std::vector< Vertex_handle > typeVectorVertex;
-  typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex;
-  typedef std::pair< Simplex_handle, bool > typePairSimplexBool;
-
-  typedef int Witness_id;
-  typedef int Landmark_id;
-  typedef std::list< Vertex_handle > ActiveWitnessList;
-
- private:
-  int nbL;  // Number of landmarks
-  Simplicial_complex& sc;  // Simplicial complex
-  
-  public:
-    /** @name Constructor
-     */
-
-    //@{
-  
-    /**
-     *  \brief Iterative construction of the relaxed witness complex.
-     *  \details The witness complex is written in sc_ basing on a matrix knn
-     *  of k nearest neighbours of the form {witnesses}x{landmarks} and 
-     *  and a matrix distances of distances to these landmarks from witnesses.
-     *  The parameter alpha defines relaxation and
-     *  limD defines the 
-     *
-     *  The line lengths in one given matrix can differ, 
-     *  however both matrices have the same corresponding line lengths.
-     *
-     *  The type KNearestNeighbors can be seen as 
-     *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-     *  Witness_range and Closest_landmark_range are random access ranges.
-     *
-     *  Constructor takes into account at most (dim+1) 
-     *  first landmarks from each landmark range to construct simplices.
-     *
-     *  Landmarks are supposed to be in [0,nbL_-1]
-     */
-    
-    template< typename KNearestNeighbours >
-    Strong_relaxed_witness_complex(std::vector< std::vector<double> > const & distances,
-                            KNearestNeighbours const & knn,
-                            Simplicial_complex & sc_, 
-                            int nbL_,
-                            double alpha2,
-                            unsigned limD) : nbL(nbL_), sc(sc_) {
-      int nbW = knn.size();
-      typeVectorVertex vv;
-      //int counter = 0;
-      /* The list of still useful witnesses
-       * it will diminuish in the course of iterations
-       */
-      ActiveWitnessList active_w;// = new ActiveWitnessList();
-      for (int i = 0; i != nbL; ++i) {
-        // initial fill of 0-dimensional simplices
-        // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
-        //counter++;
-        vv = {i};
-        sc.insert_simplex(vv, Filtration_value(0.0));
-        /* TODO Error if not inserted : normally no need here though*/
-      }
-      for (int i=0; i != nbW; ++i) {
-        // int i_end = limD+1;
-        // if (knn[i].size() < limD+1)
-        //   i_end = knn[i].size();
-        // double dist_wL = *(distances[i].begin());
-        // while (distances[i][i_end] > dist_wL + alpha2)
-        //   i_end--;
-        // add_all_witnessed_faces(distances[i].begin(),
-        //                         knn[i].begin(),
-        //                         knn[i].begin() + i_end + 1);
-        unsigned j_end = 0;
-        while (j_end < distances[i].size() && j_end <= limD && distances[i][j_end] <= distances[i][0] + alpha2) {
-          std::vector<int> simplex;
-          for (unsigned j = 0; j <= j_end; ++j)
-            simplex.push_back(knn[i][j]);
-          assert(distances[i][j_end] - distances[i][0] >= 0);
-          sc.insert_simplex_and_subfaces(simplex, distances[i][j_end] - distances[i][0]);
-          j_end++; 
-        }
-      }
-      sc.set_dimension(limD);
-    }
-
-  //@}
-  
-private:
-  /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness
-   * Iterator is needed to know until how far we can take landmarks to form simplexes
-   * simplex is the prefix of the simplexes to insert
-   * The output value indicates if the witness rests active or not
-   */
-  void add_all_witnessed_faces(std::vector<double>::const_iterator curr_d,
-                               std::vector<int>::const_iterator curr_l,
-                               std::vector<int>::const_iterator end)
-  {
-    std::vector<int> simplex;
-    std::vector<int>::const_iterator l_end = curr_l;
-    for (; l_end != end; ++l_end) {
-      std::vector<int>::const_iterator l_it = curr_l;
-      std::vector<double>::const_iterator d_it = curr_d;
-      simplex = {};
-      for (; l_it != l_end; ++l_it, ++d_it)
-        simplex.push_back(*l_it);
-      sc.insert_simplex_and_subfaces(simplex, *(d_it--) - *curr_d);
-    }
-  }
-  
-  /** \brief Check if the facets of the k-dimensional simplex witnessed 
-   *  by witness witness_id are already in the complex.
-   *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
-   */
-  bool all_faces_in(std::vector<int>& simplex, double* filtration_value)
-  {
-    std::vector< int > facet;
-    for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it)
-      {
-        facet.clear();
-        for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it)
-          if (it != not_it)
-            facet.push_back(*it);
-        Simplex_handle facet_sh = sc.find(facet);
-        if (facet_sh == sc.null_simplex())
-          return false;
-        else if (sc.filtration(facet_sh) > *filtration_value)
-          *filtration_value = sc.filtration(facet_sh);
-      }
-    return true;
-  }
-
-  bool is_face(Simplex_handle face, Simplex_handle coface)
-  {
-    // vertex range is sorted in decreasing order
-    auto fvr = sc.simplex_vertex_range(face);
-    auto cfvr = sc.simplex_vertex_range(coface);
-    auto fv_it = fvr.begin();
-    auto cfv_it = cfvr.begin();
-    while (fv_it != fvr.end() && cfv_it != cfvr.end()) {
-      if (*fv_it < *cfv_it)
-        ++cfv_it;
-      else if (*fv_it == *cfv_it) {
-        ++cfv_it;
-        ++fv_it;
-      }
-      else
-        return false;
-      
-    }
-    return (fv_it == fvr.end());
-  }
-
-  // void erase_simplex(Simplex_handle sh)
-  // {
-  //   auto siblings = sc.self_siblings(sh);
-  //   auto oncles = siblings->oncles();
-  //   int prev_vertex = siblings->parent();
-  //   siblings->members().erase(sh->first);
-  //   if (siblings->members().empty()) {
-  //     typename typedef Simplicial_complex::Siblings Siblings;
-  //     oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex));
-  //     assert(!sc.has_children(oncles->members().find(prev_vertex)));
-  //     //delete siblings;
-  //   }
-
-  // }
-  
-  void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh)
-  {
-    erase_simplex(coface_sh);
-    erase_simplex(face_sh);
-  }
-
-public:
-  // void collapse(std::vector<Simplex_handle>& simplices)
-  // {
-  //   // Get a vector of simplex handles ordered by filtration value
-  //   std::cout << sc << std::endl;
-  //   //std::vector<Simplex_handle> simplices;
-  //   for (Simplex_handle sh: sc.filtration_simplex_range())
-  //     simplices.push_back(sh);
-  //   // std::sort(simplices.begin(),
-  //   //           simplices.end(),
-  //   //           [&](Simplex_handle sh1, Simplex_handle sh2)
-  //   //           { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2);
-  //   //             return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); });
-  //   // Double iteration
-  //   auto face_it = simplices.rbegin();
-  //   while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) {
-  //     int coface_count = 0;
-  //     auto reduced_coface = simplices.rbegin();
-  //     for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it)
-  //       if (face_it != coface_it && is_face(*face_it, *coface_it)) {
-  //         coface_count++;
-  //         if (coface_count == 1)
-  //           reduced_coface = coface_it;
-  //         else
-  //           break;
-  //       }
-  //     if (coface_count == 1) {
-  //       std::cout << "Erase ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base())))
-  //         std::cout << v << " ";
-        
-  //       simplices.erase(--(reduced_coface.base()));
-  //       //elementary_collapse(*face_it, *reduced_coface);
-  //       std::cout << ") and then ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--face_it.base())))
-  //         std::cout << v << " ";
-  //       std::cout << ")\n";
-  //       simplices.erase(--((face_it++).base()));
-  //       //face_it = simplices.rbegin();
-  //       //std::cout << "Size of vector: " << simplices.size() << "\n";
-  //     } 
-  //     else
-  //       face_it++;
-  //   }
-  //   sc.initialize_filtration();
-  //   //std::cout << sc << std::endl;
-  // }
-
-  template <class Dim_lists>
-  void collapse(Dim_lists& dim_lists)
-  {
-    dim_lists.collapse();
-  }
-  
-private:
-  
-  /** Collapse recursively boundary faces of the given simplex
-   *  if its filtration is bigger than alpha_lim.  
-   */
-  void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim)
-  {
-    for (Simplex_handle face_it : sc.boundary_simplex_range()) {
-      
-    }
-      
-  }
-  
-}; //class Strong_relaxed_witness_complex
-
-} // namespace witness_complex
-  
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h b/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h
deleted file mode 100644
index a1aedd8e..00000000
--- a/src/Witness_complex/include/gudhi/Weak_relaxed_witness_complex.h
+++ /dev/null
@@ -1,379 +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) 2015  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/>.
- */
-
-#ifndef WEAK_RELAXED_WITNESS_COMPLEX_H_
-#define WEAK_RELAXED_WITNESS_COMPLEX_H_
-
-#include <boost/container/flat_map.hpp>
-#include <boost/iterator/transform_iterator.hpp>
-#include <algorithm>
-#include <utility>
-#include "gudhi/reader_utils.h"
-#include "gudhi/distance_functions.h"
-#include "gudhi/Simplex_tree.h"
-#include <vector>
-#include <list>
-#include <set>
-#include <queue>
-#include <limits>
-#include <math.h>
-#include <ctime>
-#include <iostream>
-
-// Needed for nearest neighbours
-#include <CGAL/Cartesian_d.h>
-#include <CGAL/Search_traits.h>
-#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/property_map.h>
-#include <CGAL/Epick_d.h>
-#include <CGAL/Orthogonal_k_neighbor_search.h>
-
-#include <boost/tuple/tuple.hpp>
-#include <boost/iterator/zip_iterator.hpp>
-#include <boost/iterator/counting_iterator.hpp>
-#include <boost/range/iterator_range.hpp>
-
-// Needed for the adjacency graph in bad link search
-#include <boost/graph/graph_traits.hpp>
-#include <boost/graph/adjacency_list.hpp>
-#include <boost/graph/connected_components.hpp>
-
-namespace Gudhi {
-
-namespace witness_complex {
-
-  /** \addtogroup simplex_tree
-   *  Witness complex is a simplicial complex defined on two sets of points in \f$\mathbf{R}^D\f$:
-   *  \f$W\f$ set of witnesses and \f$L \subseteq W\f$ set of landmarks. The simplices are based on points in \f$L\f$
-   *  and a simplex belongs to the witness complex if and only if it is witnessed (there exists a point \f$w \in W\f$ such that
-   *  w is closer to the vertices of this simplex than others) and all of its faces are witnessed as well. 
-   */
-template< class Simplicial_complex >
-class Weak_relaxed_witness_complex {
-
-private:    
-  struct Active_witness {
-    int witness_id;
-    int landmark_id;
-    
-    Active_witness(int witness_id_, int landmark_id_)
-      : witness_id(witness_id_),
-        landmark_id(landmark_id_) { }
-    };
-
-private:
-  typedef typename Simplicial_complex::Simplex_handle Simplex_handle;
-  typedef typename Simplicial_complex::Vertex_handle Vertex_handle;
-  typedef typename Simplicial_complex::Filtration_value FT;
-
-  typedef std::vector< double > Point_t;
-  typedef std::vector< Point_t > Point_Vector;
-
-  // typedef typename Simplicial_complex::Filtration_value Filtration_value;
-  typedef std::vector< Vertex_handle > typeVectorVertex;
-  typedef std::pair< typeVectorVertex, Filtration_value> typeSimplex;
-  typedef std::pair< Simplex_handle, bool > typePairSimplexBool;
-
-  typedef int Witness_id;
-  typedef int Landmark_id;
-  typedef std::list< Vertex_handle > ActiveWitnessList;
-
- private:
-  int nbL;  // Number of landmarks
-  Simplicial_complex& sc;  // Simplicial complex
-  
-  public:
-    /** @name Constructor
-     */
-
-    //@{
-  
-    /**
-     *  \brief Iterative construction of the relaxed witness complex.
-     *  \details The witness complex is written in sc_ basing on a matrix knn
-     *  of k nearest neighbours of the form {witnesses}x{landmarks} and 
-     *  and a matrix distances of distances to these landmarks from witnesses.
-     *  The parameter alpha defines relaxation and
-     *  limD defines the 
-     *
-     *  The line lengths in one given matrix can differ, 
-     *  however both matrices have the same corresponding line lengths.
-     *
-     *  The type KNearestNeighbors can be seen as 
-     *  Witness_range<Closest_landmark_range<Vertex_handle>>, where
-     *  Witness_range and Closest_landmark_range are random access ranges.
-     *
-     *  Constructor takes into account at most (dim+1) 
-     *  first landmarks from each landmark range to construct simplices.
-     *
-     *  Landmarks are supposed to be in [0,nbL_-1]
-     */
-    
-    template< typename KNearestNeighbours >
-    Weak_relaxed_witness_complex(std::vector< std::vector<double> > const & distances,
-                            KNearestNeighbours const & knn,
-                            Simplicial_complex & sc_, 
-                            int nbL_,
-                            double alpha,
-                            int limD) : nbL(nbL_), sc(sc_) {
-      int nbW = knn.size();
-      typeVectorVertex vv;
-      //int counter = 0;
-      /* The list of still useful witnesses
-       * it will diminuish in the course of iterations
-       */
-      ActiveWitnessList active_w;// = new ActiveWitnessList();
-      for (int i = 0; i != nbL; ++i) {
-        // initial fill of 0-dimensional simplices
-        // by doing it we don't assume that landmarks are necessarily witnesses themselves anymore
-        //counter++;
-        vv = {i};
-        sc.insert_simplex(vv, Filtration_value(0.0));
-        /* TODO Error if not inserted : normally no need here though*/
-      }
-      int k = 1; /* current dimension in iterative construction */
-      for (int i=0; i != nbW; ++i)
-        active_w.push_back(i);
-      while (!active_w.empty() && k <= limD && k < nbL )
-        {
-          typename ActiveWitnessList::iterator aw_it = active_w.begin();
-          while (aw_it != active_w.end())
-            {
-              std::vector<int> simplex;
-              bool ok = add_all_faces_of_dimension(k,
-                                                   alpha,
-                                                   std::numeric_limits<double>::infinity(),
-                                                   distances[*aw_it].begin(),
-                                                   knn[*aw_it].begin(),
-                                                   simplex,
-                                                   knn[*aw_it].end());
-              if (!ok)
-                active_w.erase(aw_it++); //First increase the iterator and then erase the previous element
-              else
-                aw_it++;
-            }
-          k++;
-        }
-      sc.set_dimension(limD);
-    }
-
-  //@}
-  
-private:
-  /* \brief Adds recursively all the faces of a certain dimension dim witnessed by the same witness
-   * Iterator is needed to know until how far we can take landmarks to form simplexes
-   * simplex is the prefix of the simplexes to insert
-   * The output value indicates if the witness rests active or not
-   */
-  bool add_all_faces_of_dimension(int dim,
-                                  double alpha,
-                                  double norelax_dist,
-                                  std::vector<double>::const_iterator curr_d,
-                                  std::vector<int>::const_iterator curr_l,
-                                  std::vector<int>& simplex,
-                                  std::vector<int>::const_iterator end)
-  {
-    if (curr_l == end)
-      return false;
-    bool will_be_active = false;
-    std::vector<int>::const_iterator l_it = curr_l;
-    std::vector<double>::const_iterator d_it = curr_d;
-    if (dim > 0)
-      for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) {
-        simplex.push_back(*l_it);
-        if (sc.find(simplex) != sc.null_simplex())
-          will_be_active = add_all_faces_of_dimension(dim-1,
-                                                      alpha,
-                                                      norelax_dist,
-                                                      d_it+1,
-                                                      l_it+1,
-                                                      simplex,
-                                                      end) || will_be_active;
-        simplex.pop_back();
-        // If norelax_dist is infinity, change to first omitted distance
-        if (*d_it <= norelax_dist)
-          norelax_dist = *d_it;
-        will_be_active = add_all_faces_of_dimension(dim-1,
-                                                    alpha,
-                                                    norelax_dist,
-                                                    d_it+1,
-                                                    l_it+1,
-                                                    simplex,
-                                                    end) || will_be_active;
-      } 
-    else if (dim == 0)
-      for (; *d_it - alpha <= norelax_dist && l_it != end; ++l_it, ++d_it) {
-        simplex.push_back(*l_it);
-        double filtration_value = 0;
-        // if norelax_dist is infinite, relaxation is 0.
-        if (*d_it > norelax_dist) 
-          filtration_value = *d_it - norelax_dist; 
-        if (all_faces_in(simplex, &filtration_value)) {
-          will_be_active = true;
-          sc.insert_simplex(simplex, filtration_value);
-        }
-        simplex.pop_back();
-        // If norelax_dist is infinity, change to first omitted distance
-        if (*d_it < norelax_dist)
-          norelax_dist = *d_it;
-      } 
-    return will_be_active;
-  }
-  
-  /** \brief Check if the facets of the k-dimensional simplex witnessed 
-   *  by witness witness_id are already in the complex.
-   *  inserted_vertex is the handle of the (k+1)-th vertex witnessed by witness_id
-   */
-  bool all_faces_in(std::vector<int>& simplex, double* filtration_value)
-  {
-    std::vector< int > facet;
-    for (std::vector<int>::iterator not_it = simplex.begin(); not_it != simplex.end(); ++not_it)
-      {
-        facet.clear();
-        for (std::vector<int>::iterator it = simplex.begin(); it != simplex.end(); ++it)
-          if (it != not_it)
-            facet.push_back(*it);
-        Simplex_handle facet_sh = sc.find(facet);
-        if (facet_sh == sc.null_simplex())
-          return false;
-        else if (sc.filtration(facet_sh) > *filtration_value)
-          *filtration_value = sc.filtration(facet_sh);
-      }
-    return true;
-  }
-
-  bool is_face(Simplex_handle face, Simplex_handle coface)
-  {
-    // vertex range is sorted in decreasing order
-    auto fvr = sc.simplex_vertex_range(face);
-    auto cfvr = sc.simplex_vertex_range(coface);
-    auto fv_it = fvr.begin();
-    auto cfv_it = cfvr.begin();
-    while (fv_it != fvr.end() && cfv_it != cfvr.end()) {
-      if (*fv_it < *cfv_it)
-        ++cfv_it;
-      else if (*fv_it == *cfv_it) {
-        ++cfv_it;
-        ++fv_it;
-      }
-      else
-        return false;
-      
-    }
-    return (fv_it == fvr.end());
-  }
-
-  // void erase_simplex(Simplex_handle sh)
-  // {
-  //   auto siblings = sc.self_siblings(sh);
-  //   auto oncles = siblings->oncles();
-  //   int prev_vertex = siblings->parent();
-  //   siblings->members().erase(sh->first);
-  //   if (siblings->members().empty()) {
-  //     typename typedef Simplicial_complex::Siblings Siblings;
-  //     oncles->members().find(prev_vertex)->second.assign_children(new Siblings(oncles, prev_vertex));
-  //     assert(!sc.has_children(oncles->members().find(prev_vertex)));
-  //     //delete siblings;
-  //   }
-
-  // }
-  
-  void elementary_collapse(Simplex_handle face_sh, Simplex_handle coface_sh)
-  {
-    erase_simplex(coface_sh);
-    erase_simplex(face_sh);
-  }
-
-public:
-  // void collapse(std::vector<Simplex_handle>& simplices)
-  // {
-  //   // Get a vector of simplex handles ordered by filtration value
-  //   std::cout << sc << std::endl;
-  //   //std::vector<Simplex_handle> simplices;
-  //   for (Simplex_handle sh: sc.filtration_simplex_range())
-  //     simplices.push_back(sh);
-  //   // std::sort(simplices.begin(),
-  //   //           simplices.end(),
-  //   //           [&](Simplex_handle sh1, Simplex_handle sh2)
-  //   //           { double f1 = sc.filtration(sh1), f2 = sc.filtration(sh2);
-  //   //             return (f1 > f2) || (f1 >= f2 && sc.dimension(sh1) > sc.dimension(sh2)); });
-  //   // Double iteration
-  //   auto face_it = simplices.rbegin();
-  //   while (face_it != simplices.rend() && sc.filtration(*face_it) != 0) {
-  //     int coface_count = 0;
-  //     auto reduced_coface = simplices.rbegin();
-  //     for (auto coface_it = simplices.rbegin(); coface_it != simplices.rend() && sc.filtration(*coface_it) != 0; ++coface_it)
-  //       if (face_it != coface_it && is_face(*face_it, *coface_it)) {
-  //         coface_count++;
-  //         if (coface_count == 1)
-  //           reduced_coface = coface_it;
-  //         else
-  //           break;
-  //       }
-  //     if (coface_count == 1) {
-  //       std::cout << "Erase ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--reduced_coface.base())))
-  //         std::cout << v << " ";
-        
-  //       simplices.erase(--(reduced_coface.base()));
-  //       //elementary_collapse(*face_it, *reduced_coface);
-  //       std::cout << ") and then ( ";
-  //       for (auto v: sc.simplex_vertex_range(*(--face_it.base())))
-  //         std::cout << v << " ";
-  //       std::cout << ")\n";
-  //       simplices.erase(--((face_it++).base()));
-  //       //face_it = simplices.rbegin();
-  //       //std::cout << "Size of vector: " << simplices.size() << "\n";
-  //     } 
-  //     else
-  //       face_it++;
-  //   }
-  //   sc.initialize_filtration();
-  //   //std::cout << sc << std::endl;
-  // }
-
-  template <class Dim_lists>
-  void collapse(Dim_lists& dim_lists)
-  {
-    dim_lists.collapse();
-  }
-  
-private:
-  
-  /** Collapse recursively boundary faces of the given simplex
-   *  if its filtration is bigger than alpha_lim.  
-   */
-  void rec_boundary_collapse(Simplex_handle sh, FT alpha_lim)
-  {
-    for (Simplex_handle face_it : sc.boundary_simplex_range()) {
-      
-    }
-      
-  }
-  
-}; //class Weak_relaxed_witness_complex
-
-} // namespace witness_complex
-  
-} // namespace Gudhi
-
-#endif
diff --git a/src/Witness_complex/include/gudhi/Witness_complex.h b/src/Witness_complex/include/gudhi/Witness_complex.h
index e732cb18..a16f9270 100644
--- a/src/Witness_complex/include/gudhi/Witness_complex.h
+++ b/src/Witness_complex/include/gudhi/Witness_complex.h
@@ -173,8 +173,6 @@ private:
         else
           aw_it++;
       } 
-      std::cout << "Active witnesses after dim=" << k << " is finished: " << active_witnesses.size() << "\n";
-      std::cout << complex << "\n";
       k++;
     }
     return true;
@@ -233,11 +231,8 @@ private:
         simplex.push_back(l_it->first);
         double filtration_value = 0;
         // if norelax_dist is infinite, relaxation is 0.
-        //std::cout << "landmark_id=" << l_it->first << " distance=" << l_it->second << "\n";
-        // std::size_t landmark_id = l_it->first;
-        // double distance = l_it->second;
         if (l_it->second > norelax_dist2) 
-          filtration_value = l_it->second - norelax_dist2; 
+          filtration_value = l_it->second - norelax_dist2;
         if (all_faces_in(simplex, &filtration_value, sc)) {
           will_be_active = true;
           sc.insert_simplex(simplex, filtration_value);