Changes piled up for months

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

Former-commit-id: 0447901e608890eb607456fd12f3ea53547b8f10

11 files changed, 3851 insertions, 820 deletions

diff --git a/src/Witness_complex/example/CMakeLists.txt b/src/Witness_complex/example/CMakeLists.txt
index 23919b4a..ff372d16 100644
--- a/src/Witness_complex/example/CMakeLists.txt
+++ b/src/Witness_complex/example/CMakeLists.txt
@@ -1,16 +1,7 @@
 cmake_minimum_required(VERSION 2.6)
 project(GUDHIWitnessComplex)
-#cmake -DCGAL_DIR=~/GitDrive/CGAL/ ../../..
 
-#if (CGAL_FOUND)
-   #message(STATUS "CGAL version: ${CGAL_VERSION}.")
-   #include( ${CGAL_USE_FILE} )
-
-   #find_package(Eigen3 3.1.0)
-   #include( ${EIGEN3_USE_FILE} )
-
-   #INCLUDE_DIRECTORIES(${EIGEN3_INCLUDE_DIR})
-   #INCLUDE_DIRECTORIES(${CGAL_INCLUDE_DIRS})
+# A simple example
    add_executable ( simple_witness_complex simple_witness_complex.cpp )
    add_test(simple_witness_complex ${CMAKE_CURRENT_BINARY_DIR}/simple_witness_complex)
 
@@ -21,7 +12,7 @@ project(GUDHIWitnessComplex)
    add_executable( witness_complex_from_off witness_complex_from_off.cpp )
    
    add_executable( witness_complex_from_wl_matrix witness_complex_from_wl_matrix.cpp )
-#endif()
+
 
 # An example with Simplex-tree using CGAL alpha_shapes_3
 
@@ -51,6 +42,7 @@ if(CGAL_FOUND)
     if (EIGEN3_FOUND)
       message(STATUS "Eigen3 version: ${EIGEN3_VERSION}.")
       include( ${EIGEN3_USE_FILE} )
+      message(STATUS "Eigen3 use file: ${EIGEN3_USE_FILE}.")
       include_directories (BEFORE "../../include")
 
       add_executable ( witness_complex_knn_landmarks witness_complex_knn_landmarks.cpp )
diff --git a/src/Witness_complex/example/protected_sets/output_tikz.h b/src/Witness_complex/example/protected_sets/output_tikz.h
new file mode 100644
index 00000000..edfd9a5f
--- /dev/null
+++ b/src/Witness_complex/example/protected_sets/output_tikz.h
@@ -0,0 +1,67 @@
+#ifndef OUTPUT_TIKZ_H
+#define OUTPUT_TIKZ_H
+
+#include <vector>
+#include <string>
+#include <algorithm>
+#include <fstream>
+#include <cmath>
+
+void write_tikz_plot(std::vector<FT> data, std::string filename)
+{
+  int n = data.size();
+  FT vmax = *(std::max_element(data.begin(), data.end()));
+  //std::cout << std::log10(vmax) << " " << std::floor(std::log10(vmax));
+  
+  FT order10 = pow(10,std::floor(std::log10(vmax)));
+  int digit = std::floor( vmax / order10) + 1;
+  if (digit == 4 || digit == 6) digit = 5;
+  if (digit > 6) digit = 10;
+  FT plot_max = digit*order10;
+  std::cout << plot_max << " " << vmax;
+  FT hstep = 10.0/(n-1);
+  FT wstep = 10.0 / plot_max;
+
+  std::cout << "(eps_max-eps_min)/(N-48) = " << (vmax-*data.begin())/(data.size()-48) << "\n";
+  std::ofstream ofs(filename, std::ofstream::out);
+
+  ofs <<
+    "\\documentclass{standalone}\n" <<
+    "\\usepackage[utf8]{inputenc}\n" <<
+    "\\usepackage{amsmath}\n" <<
+    "\\usepackage{tikz}\n\n" <<
+    "\\begin{document}\n" <<
+    "\\begin{tikzpicture}\n";
+  
+  ofs <<  "\\draw[->] (0,0) -- (0,11);" << std::endl <<
+    "\\draw[->] (0,0) -- (11,0);" << std::endl <<
+    "\\foreach \\i in {1,...,10}" << std::endl <<
+    "\\draw (0,\\i) -- (-0.05,\\i);" << std::endl <<
+    "\\foreach \\i in {1,...,10}" << std::endl <<
+    "\\draw (\\i,0) -- (\\i,-0.05);" << std::endl << std::endl <<
+
+    "\\foreach \\i in {1,...,10}" << std::endl <<
+    "\\draw[dashed] (-0.05,\\i) -- (11,\\i);" << std::endl << std::endl <<
+    
+    "\\node at (-0.5,11) {$*$}; " << std::endl <<
+    "\\node at (11,-0.5) {$*$}; " << std::endl <<
+    "\\node at (-0.5,-0.5) {0}; " << std::endl <<
+    "\\node at (-0.5,10) {" << plot_max << "}; " << std::endl <<
+    "%\\node at (10,-0.5) {2}; " << std::endl;
+
+  ofs << "\\draw[red] (0," << wstep*data[0] << ")";
+  for (int i = 1; i < n; ++i)
+    ofs << " -- (" << hstep*i << "," << wstep*data[i] << ")";
+  ofs << ";\n";
+
+  ofs <<
+    "\\end{tikzpicture}\n" <<
+    "\\end{document}";
+  
+  ofs.close();
+
+
+  
+}
+
+#endif
diff --git a/src/Witness_complex/example/protected_sets/protected_sets.h b/src/Witness_complex/example/protected_sets/protected_sets.h
new file mode 100644
index 00000000..ec627808
--- /dev/null
+++ b/src/Witness_complex/example/protected_sets/protected_sets.h
@@ -0,0 +1,597 @@
+#ifndef PROTECTED_SETS_H
+#define PROTECTED_SETS_H
+
+#include <algorithm>
+#include <CGAL/Cartesian_d.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/Euclidean_distance.h>
+#include <CGAL/Kernel_d/Sphere_d.h>
+#include <CGAL/Kernel_d/Hyperplane_d.h>
+#include <CGAL/Kernel_d/Vector_d.h>
+
+#include <CGAL/point_generators_d.h>
+#include <CGAL/constructions_d.h>
+
+
+typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
+typedef K::Point_d Point_d;
+typedef K::Vector_d Vector_d;
+typedef K::Oriented_side_d Oriented_side_d;
+typedef K::Has_on_positive_side_d Has_on_positive_side_d;
+typedef K::Sphere_d Sphere_d;
+typedef K::Hyperplane_d Hyperplane_d;
+
+typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation;
+typedef Delaunay_triangulation::Facet Facet;
+typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex;
+typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle;
+
+typedef std::vector<Point_d> Point_Vector;
+typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
+
+FT  _sfty = pow(10,-14);
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+// AUXILLARY FUNCTIONS
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well 
+ *  
+ *  W is the initial point vector
+ *  chosen_landmark is the index of the chosen point in W
+ *  landmarks_ind is the vector of indices of already chosen points in W
+ *  delaunay is the Delaunay triangulation
+ *  landmark_count is the current number of chosen vertices
+ *  torus is true iff you are working on a flat torus [-1,1]^d
+ *  OUT: Vertex handle to the newly inserted point
+ */
+Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus)
+{
+  if (!torus)
+    {
+      Delaunay_vertex v =delaunay.insert(W[chosen_landmark]);
+      landmarks_ind.push_back(chosen_landmark);
+      landmark_count++;
+      return v;
+    }
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      Delaunay_vertex v;
+      for (int i = 0; i < nb_cells; ++i)
+        {
+          std::vector<FT> point;
+          int cell_i = i;
+          for (int l = 0; l < D; ++l)
+            {
+              point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1));
+              cell_i /= 3;
+            }
+          v = delaunay.insert(point);
+        }
+      landmarks_ind.push_back(chosen_landmark);
+      landmark_count++;
+      return v;
+    }
+}
+
+/** Small check if the vertex v is in the full cell fc
+ */
+
+bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc)
+{
+  for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it)
+    if (*v_it == v)
+      return true;
+  return false;
+}
+
+/** Fill chosen point vector from indices with copies if you are working on a flat torus
+ *  
+ *  IN:  W is the point vector
+ *  OUT: landmarks is the output vector
+ *  IN:  landmarks_ind is the vector of indices
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */
+
+void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus)
+{
+  if (!torus)
+    for (unsigned j = 0; j < landmarks_ind.size(); ++j)
+      landmarks.push_back(W[landmarks_ind[j]]);
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      int nbL = landmarks_ind.size();
+      // Fill landmarks
+      for (int i = 0; i < nb_cells-1; ++i)
+        for (int j = 0; j < nbL; ++j)
+          {
+            int cell_i = i;
+            Point_d point;
+            for (int l = 0; l < D; ++l)
+              {
+                point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1));
+                cell_i /= 3;
+              }
+            landmarks.push_back(point);
+          }
+    }
+}
+
+/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices
+ *
+ *  IN: t is the Delaunay triangulation
+ *  OUT: full_cells is the output vector
+ */
+
+void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells)
+{
+  // Store vertex indices in a map
+  int ind = 0; //index of a vertex
+  std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex;
+  for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+    if (t.is_infinite(v_it))
+      continue;
+    else
+      index_of_vertex[v_it] = ind++;
+  // Write full cells as vectors in full_cells
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      std::vector<int> cell;
+      for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it)
+        cell.push_back(index_of_vertex[*v_it]);
+      full_cells.push_back(cell);
+    }
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// IS VIOLATED TEST
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Check if a newly created cell is protected from old vertices
+ *  
+ *  t is the Delaunay triangulation
+ *  vertices is the vector containing the point to insert and a facet f in t
+ *  v1 is the vertex of t, such that f and v1 form a simplex
+ *  v2 is the vertex of t, such that f and v2 form another simplex
+ *  delta is the protection constant
+ *  power_protection is true iff the delta-power protection is used
+ */
+
+bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0)
+{
+  assert(vertices.size() == vertices[0].size() ||
+         vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1
+  assert(v1 != v2);
+  if (vertices.size() == vertices[0].size() + 1)
+  // FINITE CASE
+    {
+      Sphere_d cs(vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0]));
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it)) 
+          {
+            //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); 
+            if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+              {
+                FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point());
+                if (!power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+                    return true;
+                if (power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+                    return true;                 
+              }
+          }
+      */
+      // Check if the simplex is thick enough
+      Hyperplane_d tau_h(vertices.begin()+1, vertices.end());
+      Vector_d orth_tau = tau_h.orthogonal_vector();
+      /*
+               p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1));
+      */
+      //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n";
+      FT orth_length = sqrt(orth_tau.squared_length());
+      K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it;
+      // Compute the altitude
+      FT h = 0;
+      for (o_it = orth_tau.cartesian_begin(),
+           p_it = vertices.begin()->cartesian_begin(),
+           s_it = (vertices.begin()+1)->cartesian_begin();
+           o_it != orth_tau.cartesian_end();
+           ++o_it, ++p_it, ++s_it)
+        h += (*o_it)*(*p_it - *s_it)/orth_length;
+      h = fabs(h);
+      // Is the center inside the box?
+      bool inside_the_box = true;
+      for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it)
+        if (*c_it > 1.0 || *c_it < -1.0)
+          {
+            inside_the_box = false; break;
+          }
+      if (inside_the_box && h/r < theta0)
+        return true;
+      if (!t.is_infinite(v1))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+              return true;
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+              return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+              return true;
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+              return true;
+        }
+    }
+  else
+  // INFINITE CASE
+    {
+      Delaunay_triangulation::Vertex_iterator v = t.vertices_begin();
+      while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end())
+        v++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE);
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it))
+          if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+            {
+              std::vector<FT> coords;
+              Point_d p = v_it->point();
+              auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+              for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+                coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+              Point_d p_delta = Point_d(coords);
+              bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+              bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+              if (!p_is_inside && p_delta_is_inside)
+                return true;
+            }
+      */
+      if (!t.is_infinite(v1))
+        {
+          std::vector<FT> coords;
+          Point_d p = v1->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          std::vector<FT> coords;
+          Point_d p = v2->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+    }
+  return false;
+}
+  
+/** Auxillary recursive function to check if the point p violates the protection of the cell c and
+ *  if there is a violation of an eventual new cell
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  c is the current cell (simplex)
+ *  parent_cell is the parent cell (simplex)
+ *  index is the index of the facet between c and parent_cell from parent_cell's point of view
+ *  D is the dimension of the triangulation
+ *  delta is the protection constant
+ *  marked_cells is the vector of all visited cells containing p in their circumscribed ball
+ *  power_protection is true iff you are working with delta-power protection
+ *
+ *  OUT: true iff inserting p hasn't produced any violation so far
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0)
+{
+  Euclidean_distance ed;
+  std::vector<Point_d> vertices;
+  if (!t.is_infinite(c))
+    {
+      // if the cell is finite, we look if the protection is violated
+      for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it)
+        vertices.push_back((*v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(ed.transformed_distance(center_cs, vertices[0]));
+      FT dist2 = ed.transformed_distance(center_cs, p);
+      // if the new point is inside the protection ball of a non conflicting simplex
+      if (!power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+          return true;
+      if (power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+          return true;    
+      // if the new point is inside the circumscribing ball : continue violation searching on neighbours
+      //if (dist2 < r*r)
+      //if (dist2 < (5*r+delta)*(5*r+delta))
+      if (dist2 < r*r)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0))
+                return true;
+            }
+        }
+      // if the new point is outside the protection sphere
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is guaranteed to be finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;            
+        }
+    } 
+  else
+    {
+      // Inside of the convex hull is + side. Outside is - side.
+      for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+        if (!t.is_infinite(*vh_it))
+          vertices.push_back((*vh_it)->point());
+      Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin();
+      while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c))
+        v_it++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE);
+      //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point());
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      std::vector<FT> coords;
+      auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+      for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+        coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+      Point_d p_delta = Point_d(coords);
+      bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO);
+      bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+
+      // If we work with power protection, we just ignore any conflicts
+      if (!power_protection && !p_is_inside && p_delta_is_inside)
+        return true;
+      //if the cell is infinite we look at the neighbours regardless
+      if (p_is_inside)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);    
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0))
+                return true;
+            }
+        }
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is finite if the parent cell is finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              if (!t.is_infinite(parent_cell->vertex(i)))
+                vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;
+        }
+    }
+  //c->tds_data().clear_visited();
+  //marked_cells.pop_back();
+  return false;
+}
+
+/** Checks if inserting the point p in t will make conflicts
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  D is the dimension of triangulation
+ *  delta is the protection constant
+ *  power_protection is true iff you are working with delta-power protection
+ *  OUT: true iff inserting p produces a violation of delta-protection.
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0)
+{
+  Euclidean_distance ed;
+  Delaunay_triangulation::Vertex_handle v;
+  Delaunay_triangulation::Face f(t.current_dimension()); 
+  Delaunay_triangulation::Facet ft; 
+  Delaunay_triangulation::Full_cell_handle c; 
+  Delaunay_triangulation::Locate_type lt;
+  std::vector<Full_cell_handle> marked_cells;
+  c = t.locate(p, lt, f, ft, v);
+  bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0);
+  for (Full_cell_handle fc : marked_cells)
+    fc->tds_data().clear();
+  return violation_existing_cells;
+}
+
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!//
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////
+// LANDMARK CHOICE PROCEDURE
+///////////////////////////////////////////////////////////////////////
+
+/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call.
+ *
+ *  IN:  W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d
+ *  IN:  nbP is the size of W
+ *  OUT: landmarks is the output vector for the points
+ *  OUT: landmarks_ind is the output vector for the indices of the selected points in W
+ *  IN:  delta is the constant of protection
+ *  OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */ 
+
+void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT delta, std::vector<std::vector<int>>& full_cells, bool torus, bool power_protection, FT theta0)
+{
+  bool return_ = true;
+  unsigned D = W[0].size();
+  Torus_distance td;
+  Euclidean_distance ed;
+  Delaunay_triangulation t(D);
+  CGAL::Random rand;
+  int landmark_count = 0;
+  std::list<int> index_list;
+  //      shuffle the list of indexes (via a vector)
+  {
+    std::vector<int> temp_vector;
+    for (int i = 0; i < nbP; ++i)
+      temp_vector.push_back(i);
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed));
+    //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end());
+    for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it)
+      index_list.push_front(*it);
+  }
+  if (!torus)
+    for (unsigned pos1 = 0; pos1 < D+1; ++pos1)
+      {
+        std::vector<FT> point;
+        for (unsigned i = 0; i < pos1; ++i)
+          point.push_back(-1);
+        if (pos1 != D)
+          point.push_back(1);
+        for (unsigned i = pos1+1; i < D; ++i)
+          point.push_back(0);
+        assert(point.size() == D);
+        W[index_list.front()] = Point_d(point);
+        insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+        index_list.pop_front();
+      }
+  else if (D == 2)
+    {
+      for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 2; ++j)
+          {
+            W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0});
+            insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+            index_list.pop_front();
+            W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j});
+            insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+            index_list.pop_front();
+          }
+    }
+  else
+    std::cout << "No torus starter available for dim>2\n";
+  std::list<int>::iterator list_it = index_list.begin();
+  while (list_it != index_list.end())
+    {
+      if (!is_violating_protection(W[*list_it], t, D, delta, power_protection, theta0))
+        {
+          // If no conflicts then insert in every copy of T^3
+
+          insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus);
+          if (return_)
+            {
+              index_list.erase(list_it);
+              list_it = index_list.begin();
+            }
+          else
+            index_list.erase(list_it++);
+          /*
+          // PIECE OF CODE FOR DEBUGGING PURPOSES
+          
+          Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count);
+          if (triangulation_is_protected(t, delta))
+            {
+              index_list.erase(list_it);
+              list_it = index_list.begin();
+            }
+          else
+            { //THAT'S WHERE SOMETHING'S WRONG
+              t.remove(inserted_v);
+              landmarks_ind.pop_back();
+              landmark_count--;
+              write_delaunay_mesh(t, W[*list_it], is2d);
+              is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore
+            }
+          */
+          //std::cout << "index_list_size() = " << index_list.size() << "\n";
+        }
+      else
+        {
+          list_it++;
+          //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n";
+        }
+      //if (list_it != index_list.end())
+      //  write_delaunay_mesh(t, W[*list_it], is2d);
+    }
+  fill_landmarks(W, landmarks, landmarks_ind, torus);
+  fill_full_cell_vector(t, full_cells);
+  /*
+  if (triangulation_is_protected(t, delta))
+    std::cout << "Triangulation is ok\n";
+  else
+    {
+      std::cout << "Triangulation is BAD!! T_T しくしく!\n";
+    }
+  */
+  //write_delaunay_mesh(t, W[0], is2d);
+  //std::cout << t << std::endl;
+}
+
+#endif
diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp b/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp
new file mode 100644
index 00000000..f3df3f1e
--- /dev/null
+++ b/src/Witness_complex/example/protected_sets/protected_sets_paper.cpp
@@ -0,0 +1,610 @@
+#ifndef PROTECTED_SETS_H
+#define PROTECTED_SETS_H
+
+#include <algorithm>
+#include <CGAL/Cartesian_d.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/Euclidean_distance.h>
+#include <CGAL/Kernel_d/Sphere_d.h>
+#include <CGAL/Kernel_d/Hyperplane_d.h>
+#include <CGAL/Kernel_d/Vector_d.h>
+
+typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
+typedef K::Point_d Point_d;
+typedef K::Vector_d Vector_d;
+typedef K::Oriented_side_d Oriented_side_d;
+typedef K::Has_on_positive_side_d Has_on_positive_side_d;
+typedef K::Sphere_d Sphere_d;
+typedef K::Hyperplane_d Hyperplane_d;
+
+typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation;
+typedef Delaunay_triangulation::Facet Facet;
+typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex;
+typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle;
+
+typedef std::vector<Point_d> Point_Vector;
+typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
+
+FT  _sfty = pow(10,-14);
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+// AUXILLARY FUNCTIONS
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well 
+ *  
+ *  W is the initial point vector
+ *  chosen_landmark is the index of the chosen point in W
+ *  landmarks_ind is the vector of indices of already chosen points in W
+ *  delaunay is the Delaunay triangulation
+ *  landmark_count is the current number of chosen vertices
+ *  torus is true iff you are working on a flat torus [-1,1]^d
+ *  OUT: Vertex handle to the newly inserted point
+ */
+Delaunay_vertex insert_delaunay_landmark_with_copies(Point_d& p, Delaunay_triangulation& delaunay, int& landmark_count, bool torus)
+{
+  if (!torus)
+    {
+      Delaunay_vertex v =delaunay.insert(p);
+      landmark_count++;
+      return v;
+    }
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      Delaunay_vertex v;
+      for (int i = 0; i < nb_cells; ++i)
+        {
+          std::vector<FT> point;
+          int cell_i = i;
+          for (int l = 0; l < D; ++l)
+            {
+              point.push_back(p[l] + 2.0*(cell_i%3-1));
+              cell_i /= 3;
+            }
+          v = delaunay.insert(point);
+        }
+      landmark_count++;
+      return v;
+    }
+}
+
+/** Small check if the vertex v is in the full cell fc
+ */
+
+bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc)
+{
+  for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it)
+    if (*v_it == v)
+      return true;
+  return false;
+}
+
+/** Fill chosen point vector from indices with copies if you are working on a flat torus
+ *  
+ *  IN:  W is the point vector
+ *  OUT: landmarks is the output vector
+ *  IN:  landmarks_ind is the vector of indices
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */
+
+void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus)
+{
+  if (!torus)
+    for (unsigned j = 0; j < landmarks_ind.size(); ++j)
+      landmarks.push_back(W[landmarks_ind[j]]);
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      int nbL = landmarks_ind.size();
+      // Fill landmarks
+      for (int i = 0; i < nb_cells-1; ++i)
+        for (int j = 0; j < nbL; ++j)
+          {
+            int cell_i = i;
+            Point_d point;
+            for (int l = 0; l < D; ++l)
+              {
+                point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1));
+                cell_i /= 3;
+              }
+            landmarks.push_back(point);
+          }
+    }
+}
+
+/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices
+ *
+ *  IN: t is the Delaunay triangulation
+ *  OUT: full_cells is the output vector
+ */
+
+void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells)
+{
+  // Store vertex indices in a map
+  int ind = 0; //index of a vertex
+  std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex;
+  for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+    if (t.is_infinite(v_it))
+      continue;
+    else
+      index_of_vertex[v_it] = ind++;
+  // Write full cells as vectors in full_cells
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      std::vector<int> cell;
+      for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it)
+        cell.push_back(index_of_vertex[*v_it]);
+      full_cells.push_back(cell);
+    }
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// IS VIOLATED TEST
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Check if a newly created cell is protected from old vertices
+ *  
+ *  t is the Delaunay triangulation
+ *  vertices is the vector containing the point to insert and a facet f in t
+ *  v1 is the vertex of t, such that f and v1 form a simplex
+ *  v2 is the vertex of t, such that f and v2 form another simplex
+ *  delta is the protection constant
+ *  power_protection is true iff the delta-power protection is used
+ */
+
+bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0)
+{
+  assert(vertices.size() == vertices[0].size() ||
+         vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1
+  assert(v1 != v2);
+  if (vertices.size() == vertices[0].size() + 1)
+  // FINITE CASE
+    {
+      Sphere_d cs(vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0]));
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it)) 
+          {
+            //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); 
+            if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+              {
+                FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point());
+                if (!power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+                    return true;
+                if (power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+                    return true;                 
+              }
+          }
+      */
+      // Check if the simplex is thick enough
+      Hyperplane_d tau_h(vertices.begin()+1, vertices.end());
+      Vector_d orth_tau = tau_h.orthogonal_vector();
+      /*
+               p_s1 = Vector_d(*(vertices.begin()), *(vertices.begin()+1));
+      */
+      //std::cout << "||orth_tau|| = " << sqrt(orth_tau.squared_length()) << "\n";
+      FT orth_length = sqrt(orth_tau.squared_length());
+      K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it;
+      // Compute the altitude
+      FT h = 0;
+      for (o_it = orth_tau.cartesian_begin(),
+           p_it = vertices.begin()->cartesian_begin(),
+           s_it = (vertices.begin()+1)->cartesian_begin();
+           o_it != orth_tau.cartesian_end();
+           ++o_it, ++p_it, ++s_it)
+        h += (*o_it)*(*p_it - *s_it)/orth_length;
+      h = fabs(h);
+      // Is the center inside the box?
+      bool inside_the_box = true;
+      for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it)
+        if (*c_it > 1.0 || *c_it < -1.0)
+          {
+            inside_the_box = false; break;
+          }
+      if (inside_the_box && h/r < theta0)
+        return true;
+      if (!t.is_infinite(v1))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+              return true;
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+              return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+              return true;
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+              return true;
+        }
+    }
+  else
+  // INFINITE CASE
+    {
+      Delaunay_triangulation::Vertex_iterator v = t.vertices_begin();
+      while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end())
+        v++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE);
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it))
+          if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+            {
+              std::vector<FT> coords;
+              Point_d p = v_it->point();
+              auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+              for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+                coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+              Point_d p_delta = Point_d(coords);
+              bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+              bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+              if (!p_is_inside && p_delta_is_inside)
+                return true;
+            }
+      */
+      if (!t.is_infinite(v1))
+        {
+          std::vector<FT> coords;
+          Point_d p = v1->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          std::vector<FT> coords;
+          Point_d p = v2->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+    }
+  return false;
+}
+  
+/** Auxillary recursive function to check if the point p violates the protection of the cell c and
+ *  if there is a violation of an eventual new cell
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  c is the current cell (simplex)
+ *  parent_cell is the parent cell (simplex)
+ *  index is the index of the facet between c and parent_cell from parent_cell's point of view
+ *  D is the dimension of the triangulation
+ *  delta is the protection constant
+ *  marked_cells is the vector of all visited cells containing p in their circumscribed ball
+ *  power_protection is true iff you are working with delta-power protection
+ *
+ *  OUT: true iff inserting p hasn't produced any violation so far
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0)
+{
+  Euclidean_distance ed;
+  std::vector<Point_d> vertices;
+  if (!t.is_infinite(c))
+    {
+      // if the cell is finite, we look if the protection is violated
+      for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it)
+        vertices.push_back((*v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(ed.transformed_distance(center_cs, vertices[0]));
+      FT dist2 = ed.transformed_distance(center_cs, p);
+      // if the new point is inside the protection ball of a non conflicting simplex
+      if (!power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+          return true;
+      if (power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+          return true;    
+      // if the new point is inside the circumscribing ball : continue violation searching on neighbours
+      //if (dist2 < r*r)
+      //if (dist2 < (5*r+delta)*(5*r+delta))
+      if (dist2 < r*r)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0))
+                return true;
+            }
+        }
+      // if the new point is outside the protection sphere
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is guaranteed to be finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;            
+        }
+    } 
+  else
+    {
+      // Inside of the convex hull is + side. Outside is - side.
+      for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+        if (!t.is_infinite(*vh_it))
+          vertices.push_back((*vh_it)->point());
+      Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin();
+      while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c))
+        v_it++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE);
+      //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point());
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      std::vector<FT> coords;
+      auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+      for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+        coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+      Point_d p_delta = Point_d(coords);
+      bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO);
+      bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+
+      // If we work with power protection, we just ignore any conflicts
+      if (!power_protection && !p_is_inside && p_delta_is_inside)
+        return true;
+      //if the cell is infinite we look at the neighbours regardless
+      if (p_is_inside)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);    
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0))
+                return true;
+            }
+        }
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is finite if the parent cell is finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              if (!t.is_infinite(parent_cell->vertex(i)))
+                vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;
+        }
+    }
+  //c->tds_data().clear_visited();
+  //marked_cells.pop_back();
+  return false;
+}
+
+/** Checks if inserting the point p in t will make conflicts
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  D is the dimension of triangulation
+ *  delta is the protection constant
+ *  power_protection is true iff you are working with delta-power protection
+ *  OUT: true iff inserting p produces a violation of delta-protection.
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0)
+{
+  Euclidean_distance ed;
+  Delaunay_triangulation::Vertex_handle v;
+  Delaunay_triangulation::Face f(t.current_dimension()); 
+  Delaunay_triangulation::Facet ft; 
+  Delaunay_triangulation::Full_cell_handle c; 
+  Delaunay_triangulation::Locate_type lt;
+  std::vector<Full_cell_handle> marked_cells;
+  c = t.locate(p, lt, f, ft, v);
+  bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0);
+  for (Full_cell_handle fc : marked_cells)
+    fc->tds_data().clear();
+  return violation_existing_cells;
+}
+
+//////////////////////////////////////////////////////////////////////
+// INITIALIZATION
+//////////////////////////////////////////////////////////////////////
+
+void initialize(Search_Tree& W, Delaunay& t, int D, int width, bool torus)
+{
+  if (!torus)
+    std::cout << "Non-toric case is not supported\n";
+  else
+    {
+      if (D == 2)
+        {
+          FT stepx = 2.0/width;
+          FT stepy = sqrt(3)/width;
+          for (int i = 0; i < width; ++i)
+            for (int j = 0; j < floor(2*width/sqrt(3)); ++j)
+              {
+                insert_delaunay_landmark_with_copies(Point_d(step*i,))
+              }
+        }
+      else (D == 3)
+        {
+          
+        }
+      else std::cout << "T^d with d>3 not supported";
+    }
+}
+  
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!//
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////
+// LANDMARK CHOICE PROCEDURE AS IN PAPER
+///////////////////////////////////////////////////////////////////////
+
+/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call.
+ *
+ *  IN:  W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d
+ *  IN:  nbP is the size of W
+ *  OUT: landmarks is the output vector for the points
+ *  OUT: landmarks_ind is the output vector for the indices of the selected points in W
+ *  IN:  delta is the constant of protection
+ *  OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */ 
+
+template<class Search_Tree>
+void protected_delaunay_refinement(Search_Tree& W, int nbP, Point_Vector& landmarks, FT delta, bool torus, bool power_protection, FT theta0)
+{
+  bool return_ = true;
+  unsigned D = W[0].size();
+  Torus_distance td;
+  Euclidean_distance ed;
+  Delaunay_triangulation t(D);
+  CGAL::Random rand;
+  int landmark_count = 0;
+  //std::list<int> index_list;
+  //      shuffle the list of indexes (via a vector)
+  // {
+  //   std::vector<int> temp_vector;
+  //   for (int i = 0; i < nbP; ++i)
+  //     temp_vector.push_back(i);
+  //   unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+  //   std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed));
+  //   //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end());
+  //   for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it)
+  //     index_list.push_front(*it);
+  // }
+  if (torus)
+  if (D == 2)
+    // \T^2
+    {
+      for (int i = 0; i < 4; ++i)
+        for (int j = 0; j < 2; ++j)
+          {
+            W[index_list.front()] = Point_d(std::vector<FT>{i*0.5, j*1.0});
+            insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+            index_list.pop_front();
+            W[index_list.front()] = Point_d(std::vector<FT>{0.25+i*0.5, 0.5+j});
+            insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+            index_list.pop_front();
+          }
+    }
+  else if (D == 3)
+    {
+      
+    }
+    //std::cout << "No torus starter available for dim>2\n";
+  std::list<int>::iterator list_it = index_list.begin();
+  while (list_it != index_list.end())
+    {
+      if (!is_violating_protection(W[*list_it], t, D, delta, power_protection, theta0))
+        {
+          // If no conflicts then insert in every copy of T^3
+
+          insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count, torus);
+          if (return_)
+            {
+              index_list.erase(list_it);
+              list_it = index_list.begin();
+            }
+          else
+            index_list.erase(list_it++);
+          /*
+          // PIECE OF CODE FOR DEBUGGING PURPOSES
+          
+          Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count);
+          if (triangulation_is_protected(t, delta))
+            {
+              index_list.erase(list_it);
+              list_it = index_list.begin();
+            }
+          else
+            { //THAT'S WHERE SOMETHING'S WRONG
+              t.remove(inserted_v);
+              landmarks_ind.pop_back();
+              landmark_count--;
+              write_delaunay_mesh(t, W[*list_it], is2d);
+              is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore
+            }
+          */
+          //std::cout << "index_list_size() = " << index_list.size() << "\n";
+        }
+      else
+        {
+          list_it++;
+          //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n";
+        }
+      //if (list_it != index_list.end())
+      //  write_delaunay_mesh(t, W[*list_it], is2d);
+    }
+  fill_landmarks(W, landmarks, landmarks_ind, torus);
+  fill_full_cell_vector(t, full_cells);
+  /*
+  if (triangulation_is_protected(t, delta))
+    std::cout << "Triangulation is ok\n";
+  else
+    {
+      std::cout << "Triangulation is BAD!! T_T しくしく!\n";
+    }
+  */
+  //write_delaunay_mesh(t, W[0], is2d);
+  //std::cout << t << std::endl;
+}
+
+#endif
diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper.h b/src/Witness_complex/example/protected_sets/protected_sets_paper.h
new file mode 100644
index 00000000..61fcc75b
--- /dev/null
+++ b/src/Witness_complex/example/protected_sets/protected_sets_paper.h
@@ -0,0 +1,917 @@
+#ifndef PROTECTED_SETS_H
+#define PROTECTED_SETS_H
+
+#include <algorithm>
+#include <CGAL/Cartesian_d.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/Euclidean_distance.h>
+#include <CGAL/Kernel_d/Sphere_d.h>
+#include <CGAL/Kernel_d/Hyperplane_d.h>
+#include <CGAL/Kernel_d/Vector_d.h>
+
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Kd_tree.h>
+#include <CGAL/Fuzzy_sphere.h>
+
+#include <boost/heap/fibonacci_heap.hpp>
+#include <boost/heap/policies.hpp>
+
+#include "output_tikz.h"
+
+typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
+typedef K::Point_d Point_d;
+typedef K::Line_d Line_d;
+typedef K::Vector_d Vector_d;
+typedef K::Oriented_side_d Oriented_side_d;
+typedef K::Has_on_positive_side_d Has_on_positive_side_d;
+typedef K::Sphere_d Sphere_d;
+typedef K::Hyperplane_d Hyperplane_d;
+
+typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation;
+typedef Delaunay_triangulation::Facet Facet;
+typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex;
+typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle;
+
+typedef std::vector<Point_d> Point_Vector;
+typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
+
+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_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> K_neighbor_search;
+typedef K_neighbor_search::Tree Tree;
+typedef K_neighbor_search::Distance Distance;
+typedef K_neighbor_search::iterator KNS_iterator;
+typedef K_neighbor_search::iterator KNS_range;
+typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
+
+
+FT  _sfty = pow(10,-14);
+
+bool experiment1, experiment2 = false;
+
+/* Experiment 1: epsilon as function on time **********************/
+std::vector<FT> eps_vector;
+
+/* Experiment 2: R/epsilon on delta *******************************/
+std::vector<FT> epsratio_vector;
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+// AUXILLARY FUNCTIONS
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well 
+ *  
+ *  W is the initial point vector
+ *  chosen_landmark is the index of the chosen point in W
+ *  landmarks_ind is the vector of indices of already chosen points in W
+ *  delaunay is the Delaunay triangulation
+ *  landmark_count is the current number of chosen vertices
+ *  torus is true iff you are working on a flat torus [-1,1]^d
+ *  OUT: Vertex handle to the newly inserted point
+ */
+Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus)
+{
+  if (!torus)
+    {
+      Delaunay_vertex v =delaunay.insert(W[chosen_landmark]);
+      landmarks_ind.push_back(chosen_landmark);
+      landmark_count++;
+      return v;
+    }
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      Delaunay_vertex v;
+      for (int i = 0; i < nb_cells; ++i)
+        {
+          std::vector<FT> point;
+          int cell_i = i;
+          for (int l = 0; l < D; ++l)
+            {
+              point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1));
+              cell_i /= 3;
+            }
+          if (i == nb_cells/2)
+            v = delaunay.insert(point); //v = center point
+          else
+            delaunay.insert(point);
+        }
+      landmarks_ind.push_back(chosen_landmark);
+      landmark_count++;
+      return v;
+    }
+}
+
+/** Small check if the vertex v is in the full cell fc
+ */
+
+bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc)
+{
+  for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it)
+    if (*v_it == v)
+      return true;
+  return false;
+}
+
+/** Fill chosen point vector from indices with copies if you are working on a flat torus
+ *  
+ *  IN:  W is the point vector
+ *  OUT: landmarks is the output vector
+ *  IN:  landmarks_ind is the vector of indices
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */
+
+void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus)
+{
+  if (!torus)
+    for (unsigned j = 0; j < landmarks_ind.size(); ++j)
+      landmarks.push_back(W[landmarks_ind[j]]);
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      int nbL = landmarks_ind.size();
+      // Fill landmarks
+      for (int i = 0; i < nb_cells-1; ++i)
+        for (int j = 0; j < nbL; ++j)
+          {
+            int cell_i = i;
+            Point_d point;
+            for (int l = 0; l < D; ++l)
+              {
+                point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1));
+                cell_i /= 3;
+              }
+            landmarks.push_back(point);
+          }
+    }
+}
+
+/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices
+ *
+ *  IN: t is the Delaunay triangulation
+ *  OUT: full_cells is the output vector
+ */
+
+void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells)
+{
+  // Store vertex indices in a map
+  int ind = 0; //index of a vertex
+  std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex;
+  for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+    if (t.is_infinite(v_it))
+      continue;
+    else
+      index_of_vertex[v_it] = ind++;
+  // Write full cells as vectors in full_cells
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      std::vector<int> cell;
+      for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it)
+        cell.push_back(index_of_vertex[*v_it]);
+      full_cells.push_back(cell);
+    }
+}
+
+bool sphere_intersects_cube(Point_d& c, FT r)
+{
+  bool in_cube = true;
+  //  int i = 0, D = p.size();
+  for (auto xi = c.cartesian_begin(); xi != c.cartesian_end(); ++xi)
+    // if ((*xi < 1.0 || *xi > -1.0) &&
+    //     (*xi-r < 1.0 || *xi-r > -1.0) &&
+    //     (*xi+r < 1.0 || *xi+r > -1.0))
+
+    if ((*xi-r < -1.0 && *xi+r < -1.0) ||
+        (*xi-r > 1.0  && *xi+r >  1.0 ))
+      {
+        in_cube = false; break;
+      }
+  return in_cube;
+}
+
+/** Recursive function for checking if the simplex is good,
+ *  meaning it does not contain a k-face, which is not theta0^(k-1) thick 
+ */
+
+bool is_theta0_good(std::vector<Point_d>& vertices, FT theta0)
+{
+  if (theta0 > 1)
+    {
+      std::cout << "Warning! theta0 is set > 1\n";
+      return false;
+    }
+  int D = vertices.size()-1;
+  if (D <= 1)
+    return true; // Edges are always good
+  //******** Circumscribed sphere
+  Euclidean_distance ed;
+  Sphere_d cs(vertices.begin(), vertices.end());
+  FT r = sqrt(cs.squared_radius());
+  for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it)
+    {
+      std::vector<Point_d> facet;
+      for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it)
+        if (f_it != v_it)
+          facet.push_back(*f_it);
+      // Compute the altitude
+
+      if (vertices[0].size() == 3 && D == 2)
+        {
+          //Vector_d l = facet[0] - facet[1];
+          FT orth_length2 = ed.transformed_distance(facet[0],facet[1]);
+          K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it;
+          FT h = 0;
+          // Scalar product = <sp,l>
+          FT scalar = 0;
+          for (p_it = v_it->cartesian_begin(),
+               s_it = facet[0].cartesian_begin(),
+               l_it = facet[1].cartesian_begin();
+               p_it != v_it->cartesian_end();
+               ++l_it, ++p_it, ++s_it)
+            scalar += (*l_it - *s_it)*(*p_it - *s_it);
+          // Gram-Schmidt for one vector
+          for (p_it = v_it->cartesian_begin(),
+               s_it = facet[0].cartesian_begin(),
+               l_it = facet[1].cartesian_begin();
+               p_it != v_it->cartesian_end();
+               ++l_it, ++p_it, ++s_it)
+            {
+              FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2;
+              h += hx*hx;
+            }
+          h = sqrt(h);
+      
+          if (h/(2*r) < pow(theta0, D-1))
+            return false;
+          if (!is_theta0_good(facet, theta0))
+            return false;
+        }
+      else
+        {
+      Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it);
+      Vector_d orth_tau = tau_h.orthogonal_vector();
+      FT orth_length = sqrt(orth_tau.squared_length());
+      K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it;
+      FT h = 0;
+      for (o_it = orth_tau.cartesian_begin(),
+           p_it = v_it->cartesian_begin(),
+           s_it = (facet.begin())->cartesian_begin();
+           o_it != orth_tau.cartesian_end();
+           ++o_it, ++p_it, ++s_it)
+        h += (*o_it)*(*p_it - *s_it)/orth_length;
+      h = fabs(h);
+      if (h/(2*r) < pow(theta0, D-1))
+        return false;
+      if (!is_theta0_good(facet, theta0))
+        return false;
+        }
+    }
+  return true;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// IS VIOLATED TEST
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Check if a newly created cell is protected from old vertices
+ *  
+ *  t is the Delaunay triangulation
+ *  vertices is the vector containing the point to insert and a facet f in t
+ *  v1 is the vertex of t, such that f and v1 form a simplex
+ *  v2 is the vertex of t, such that f and v2 form another simplex
+ *  delta is the protection constant
+ *  power_protection is true iff the delta-power protection is used
+ */
+
+bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta, bool power_protection, FT theta0)
+{
+  assert(vertices.size() == vertices[0].size() ||
+         vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1
+  assert(v1 != v2);
+  if (vertices.size() == vertices[0].size() + 1)
+  // FINITE CASE
+    {
+      Sphere_d cs(vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0]));
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it)) 
+          {
+            //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); 
+            if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+              {
+                FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point());
+                if (!power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+                    return true;
+                if (power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+                    return true;                 
+              }
+          }
+      */
+      // Check if the simplex is theta0-good
+      if (!is_theta0_good(vertices, theta0))
+        return true;
+      // Is the center inside the box? (only Euclidean case)
+      // if (!torus)
+      //   {
+      //     bool inside_the_box = true;
+      //     for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it)
+      //       if (*c_it > 1.0 || *c_it < -1.0)
+      //         {
+      //           inside_the_box = false; break;
+      //         }
+      //     if (inside_the_box && h/r < theta0)
+      //       return true;
+      //   }
+      // Check the two vertices (if not infinite)
+      if (!t.is_infinite(v1))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+              return true;
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+              return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+              return true;
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+              return true;
+        }
+    }
+  else
+  // INFINITE CASE
+    {
+      Delaunay_triangulation::Vertex_iterator v = t.vertices_begin();
+      while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end())
+        v++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE);
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it))
+          if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+            {
+              std::vector<FT> coords;
+              Point_d p = v_it->point();
+              auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+              for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+                coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+              Point_d p_delta = Point_d(coords);
+              bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+              bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+              if (!p_is_inside && p_delta_is_inside)
+                return true;
+            }
+      */
+      if (!t.is_infinite(v1))
+        {
+          std::vector<FT> coords;
+          Point_d p = v1->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          std::vector<FT> coords;
+          Point_d p = v2->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+    }
+  return false;
+}
+  
+/** Auxillary recursive function to check if the point p violates the protection of the cell c and
+ *  if there is a violation of an eventual new cell
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  c is the current cell (simplex)
+ *  parent_cell is the parent cell (simplex)
+ *  index is the index of the facet between c and parent_cell from parent_cell's point of view
+ *  D is the dimension of the triangulation
+ *  delta is the protection constant
+ *  marked_cells is the vector of all visited cells containing p in their circumscribed ball
+ *  power_protection is true iff you are working with delta-power protection
+ *
+ *  OUT: true iff inserting p hasn't produced any violation so far
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0)
+{
+  Euclidean_distance ed;
+  std::vector<Point_d> vertices;
+  if (!t.is_infinite(c))
+    {
+      // if the cell is finite, we look if the protection is violated
+      for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it)
+        vertices.push_back((*v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(ed.transformed_distance(center_cs, vertices[0]));
+      FT dist2 = ed.transformed_distance(center_cs, p);
+      // if the new point is inside the protection ball of a non conflicting simplex
+      if (!power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+          return true;
+      if (power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+          return true;    
+      // if the new point is inside the circumscribing ball : continue violation searching on neighbours
+      //if (dist2 < r*r)
+      //if (dist2 < (5*r+delta)*(5*r+delta))
+      if (dist2 < r*r)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0))
+                return true;
+            }
+        }
+      // if the new point is outside the protection sphere
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is guaranteed to be finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;            
+        }
+    } 
+  else
+    {
+      // Inside of the convex hull is + side. Outside is - side.
+      for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+        if (!t.is_infinite(*vh_it))
+          vertices.push_back((*vh_it)->point());
+      Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin();
+      while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c))
+        v_it++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE);
+      //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point());
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      std::vector<FT> coords;
+      auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+      for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+        coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+      Point_d p_delta = Point_d(coords);
+      bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO);
+      bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+
+      // If we work with power protection, we just ignore any conflicts
+      if (!power_protection && !p_is_inside && p_delta_is_inside)
+        return true;
+      //if the cell is infinite we look at the neighbours regardless
+      if (p_is_inside)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);    
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells, power_protection, theta0))
+                return true;
+            }
+        }
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is finite if the parent cell is finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              if (!t.is_infinite(parent_cell->vertex(i)))
+                vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta, power_protection, theta0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;
+        }
+    }
+  //c->tds_data().clear_visited();
+  //marked_cells.pop_back();
+  return false;
+}
+
+/** Checks if inserting the point p in t will make conflicts
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  D is the dimension of triangulation
+ *  delta is the protection constant
+ *  power_protection is true iff you are working with delta-power protection
+ *  OUT: true iff inserting p produces a violation of delta-protection.
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta, bool power_protection, FT theta0)
+{
+  Euclidean_distance ed;
+  Delaunay_triangulation::Vertex_handle v;
+  Delaunay_triangulation::Face f(t.current_dimension()); 
+  Delaunay_triangulation::Facet ft; 
+  Delaunay_triangulation::Full_cell_handle c; 
+  Delaunay_triangulation::Locate_type lt;
+  std::vector<Full_cell_handle> marked_cells;
+  //c = t.locate(p, lt, f, ft, v);
+  c = t.locate(p);
+  bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells, power_protection, theta0);
+  for (Full_cell_handle fc : marked_cells)
+    fc->tds_data().clear();
+  return violation_existing_cells;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// INITIALIZATION
+////////////////////////////////////////////////////////////////////////
+
+// Query for a sphere near a cite in all copies of a torus
+// OUT points_inside
+void torus_search(Tree& treeW, int D, Point_d cite, FT r, std::vector<int>& points_inside)
+{
+  int nb_cells = pow(3, D);
+  Delaunay_vertex v;
+  for (int i = 0; i < nb_cells; ++i)
+    {
+      std::vector<FT> cite_copy;
+      int cell_i = i;
+      for (int l = 0; l < D; ++l)
+        {
+          cite_copy.push_back(cite[l] + 2.0*(cell_i%3-1));
+          cell_i /= 3;
+        }
+      Fuzzy_sphere fs(cite_copy, r, 0, treeW.traits());
+      treeW.search(std::insert_iterator<std::vector<int>>(points_inside, points_inside.end()), fs);
+    }
+}
+
+  
+void initialize_torus(Point_Vector& W, Tree& treeW, Delaunay_triangulation& t, FT epsilon, std::vector<int>& landmarks_ind, int& landmark_count)
+{
+  int D = W[0].size();
+  if (D == 2)
+    {
+      int xw = 6, yw = 4;
+      // Triangular lattice close to regular triangles h=0.866a ~ 0.875a   : 48p
+      for (int i = 0; i < xw; ++i)
+        for (int j = 0; j < yw; ++j)
+          {
+            Point_d cite1(std::vector<FT>{2.0/xw*i, 1.0/yw*j});
+            std::vector<int> points_inside;
+            torus_search(treeW, D, cite1, epsilon, points_inside);
+            assert(points_inside.size() > 0);
+            insert_delaunay_landmark_with_copies(W, *(points_inside.begin()),
+                                                 landmarks_ind, t, landmark_count, true);
+            Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 1.0/yw*(j+0.5)});
+            points_inside.clear();
+            torus_search(treeW, D, cite2, epsilon, points_inside);
+            assert(points_inside.size() > 0);
+            insert_delaunay_landmark_with_copies(W, *(points_inside.begin()),
+                                                 landmarks_ind, t, landmark_count, true);            
+          }
+    }
+  else if (D == 3)
+    {
+      int wd = 3;
+      // Body-centered cubic lattice : 54p
+      for (int i = 0; i < wd; ++i)
+        for (int j = 0; j < wd; ++j)
+          for (int k = 0; k < wd; ++k)
+            {
+              Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k});
+              std::vector<int> points_inside;
+              torus_search(treeW, D, cite1, epsilon, points_inside);
+              assert(points_inside.size() > 0);
+              insert_delaunay_landmark_with_copies(W, *(points_inside.begin()),
+                                                   landmarks_ind, t, landmark_count, true);
+              Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)});
+              points_inside.clear();
+              torus_search(treeW, D, cite2, epsilon, points_inside);
+              assert(points_inside.size() > 0);
+              insert_delaunay_landmark_with_copies(W, *(points_inside.begin()),
+                                                   landmarks_ind, t, landmark_count, true);
+            }
+    }
+}
+  
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!//
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+
+// Struct for R_max_heap elements
+
+struct R_max_handle
+{
+  FT value;
+  Point_d center;
+  
+  R_max_handle(FT value_, Point_d c): value(value_), center(c)
+  {}
+};
+
+struct R_max_compare
+{
+  bool operator()(const R_max_handle& rmh1, const R_max_handle& rmh2) const
+  {
+    return rmh1.value < rmh2.value;
+  }
+};
+
+// typedef boost::heap::fibonacci_heap<R_max_handle, boost::heap::compare<R_max_compare>> Heap;
+
+// void make_heap(Delaunay_triangulation& t, Heap& R_max_heap)
+// {
+//   R_max_heap.clear();
+//   for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) 
+//     {
+//       if (t.is_infinite(fc_it))
+//         continue;
+//       Point_Vector vertices;
+//       for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+//         vertices.push_back((*fc_v_it)->point());
+//       Sphere_d cs( vertices.begin(), vertices.end());
+//       Point_d csc = cs.center();
+//       FT r = sqrt(cs.squared_radius());
+//       // A ball is in the heap, if it intersects the cube
+//       bool accepted = sphere_intersects_cube(csc, sqrt(r)); 
+//       if (!accepted)
+//         continue;
+//       R_max_heap.push(R_max_handle(r, fc_it, csc));
+//     }
+// }
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+// SAMPLING RADIUS 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+R_max_handle sampling_radius(Delaunay_triangulation& t)
+{
+  FT epsilon2 = 0;
+  Point_d final_center;
+  Point_d control_point;
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin()));
+      if (epsilon2 < r2)
+        {
+          epsilon2 = r2;
+          final_center = csc;         
+          control_point = (*vertices.begin());
+        }
+    }
+  return R_max_handle(sqrt(epsilon2), final_center);
+}
+
+///////////////////////////////////////////////////////////////////////
+// LANDMARK CHOICE PROCEDURE
+///////////////////////////////////////////////////////////////////////
+
+/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call.
+ *
+ *  IN:  W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d
+ *  IN:  nbP is the size of W
+ *  OUT: landmarks is the output vector for the points
+ *  OUT: landmarks_ind is the output vector for the indices of the selected points in W
+ *  IN:  delta is the constant of protection
+ *  OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */ 
+
+void protected_delaunay(Point_Vector& W,
+                        //Point_Vector& landmarks,
+                        std::vector<int>& landmarks_ind,
+                        FT delta,
+                        FT epsilon,
+                        FT alpha,
+                        FT theta0,
+                        //std::vector<std::vector<int>>& full_cells,
+                        bool torus,
+                        bool power_protection
+                        )
+{
+  //bool return_ = true;
+  unsigned D = W[0].size();
+  int nbP = W.size();
+  Torus_distance td;
+  Euclidean_distance ed;
+  Delaunay_triangulation t(D);
+  CGAL::Random rand;
+  int landmark_count = 0;
+  std::list<int> index_list;
+  //****************** Kd Tree W
+  STraits traits(&(W[0]));
+  Tree treeW(boost::counting_iterator<std::ptrdiff_t>(0),
+         boost::counting_iterator<std::ptrdiff_t>(nbP),
+         typename Tree::Splitter(),
+         traits);
+  //      shuffle the list of indexes (via a vector)
+  {
+    std::vector<int> temp_vector;
+    for (int i = 0; i < nbP; ++i)
+      temp_vector.push_back(i);
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed));
+    //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end());
+    for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it)
+      index_list.push_front(*it);
+  }
+  //******************** Initialize point set
+  if (!torus)
+    for (unsigned pos1 = 0; pos1 < D+1; ++pos1)
+      {
+        std::vector<FT> point;
+        for (unsigned i = 0; i < pos1; ++i)
+          point.push_back(-1);
+        if (pos1 != D)
+          point.push_back(1);
+        for (unsigned i = pos1+1; i < D; ++i)
+          point.push_back(0);
+        assert(point.size() == D);
+        W[index_list.front()] = Point_d(point);
+        insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+        index_list.pop_front();
+      }
+  else
+    initialize_torus(W, treeW, t, epsilon, landmarks_ind, landmark_count);
+  //std::cout << "Size of treeW: " << treeW.size() << "\n";
+  //std::cout << "Size of t: " << t.number_of_vertices() << "\n";
+  //******************* Initialize heap for R_max
+  //Heap R_max_heap;
+  //make_heap(t, R_max_heap);
+
+  
+  R_max_handle rh = sampling_radius(t);
+  FT epsilon0 = rh.value; 
+  if (experiment1) eps_vector.push_back(pow(1/rh.value,D));
+  //******************** Iterative algorithm
+  std::vector<int> candidate_points;
+  torus_search(treeW, D,
+               rh.center,
+               alpha*rh.value,
+               candidate_points);
+  std::list<int>::iterator list_it;
+  std::vector<int>::iterator cp_it = candidate_points.begin();
+  while (cp_it != candidate_points.end())
+    {
+      if (!is_violating_protection(W[*cp_it], t, D, delta, power_protection, theta0))
+        {
+          insert_delaunay_landmark_with_copies(W, *cp_it, landmarks_ind, t, landmark_count, torus);
+          //make_heap(t, R_max_heap);
+          rh = sampling_radius(t);
+          if (experiment1) eps_vector.push_back(pow(1/rh.value,D));
+          //std::cout << "rhvalue = " << rh.value << "\n";
+          //std::cout << "D = " << 
+          candidate_points.clear();
+          torus_search(treeW, D,
+                       rh.center,
+                       alpha*rh.value,
+                       candidate_points);
+          /*
+          // PIECE OF CODE FOR DEBUGGING PURPOSES
+          
+          Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count);
+          if (triangulation_is_protected(t, delta))
+            {
+              index_list.erase(list_it);
+              list_it = index_list.begin();
+            }
+          else
+            { //THAT'S WHERE SOMETHING'S WRONG
+              t.remove(inserted_v);
+              landmarks_ind.pop_back();
+              landmark_count--;
+              write_delaunay_mesh(t, W[*list_it], is2d);
+              is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore
+            }
+          */
+          //std::cout << "index_list_size() = " << index_list.size() << "\n";
+        }
+      else
+        {
+          cp_it++;
+          //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n";
+        }
+      //if (list_it != index_list.end())
+      //  write_delaunay_mesh(t, W[*list_it], is2d);
+    }
+  if (experiment2) epsratio_vector.push_back(rh.value/epsilon0);
+  std::cout << "The iteration ended when cp_count = " << candidate_points.size() << "\n";
+  std::cout << "alphaRmax = " << alpha*rh.value << "\n";
+  std::cout << "epsilon' = " << rh.value << "\n";
+  std::cout << "nbL = " << landmarks_ind.size() << "\n";
+  //fill_landmarks(W, landmarks, landmarks_ind, torus);
+  //fill_full_cell_vector(t, full_cells);
+  /*
+  if (triangulation_is_protected(t, delta))
+    std::cout << "Triangulation is ok\n";
+  else
+    {
+      std::cout << "Triangulation is BAD!! T_T しくしく!\n";
+    }
+  */
+  //write_delaunay_mesh(t, W[0], is2d);
+  //std::cout << t << std::endl;
+}
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Series of experiments
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void start_experiments(Point_Vector& W, FT theta0, std::vector<int>& landmarks_ind, FT epsilon)
+{
+  // Experiment 1
+  experiment1 = true;
+  protected_delaunay(W, landmarks_ind, 0.1*epsilon, epsilon, 0.5, 0, true, true);
+  write_tikz_plot(eps_vector,"epstime.tikz");
+  experiment1 = false;
+
+  // Experiment 2
+  // experiment2 = true;
+  // for (FT delta = 0; delta < epsilon; delta += 0.1*epsilon)
+  //   protected_delaunay(W, landmarks_ind, delta, epsilon, 0.5, 0, true, true);
+  // write_tikz_plot(epsratio_vector,"epsratio_delta.tikz");
+  // experiment2 = false;
+
+}
+  
+#endif
diff --git a/src/Witness_complex/example/protected_sets/protected_sets_paper2.h b/src/Witness_complex/example/protected_sets/protected_sets_paper2.h
new file mode 100644
index 00000000..04b5e3bc
--- /dev/null
+++ b/src/Witness_complex/example/protected_sets/protected_sets_paper2.h
@@ -0,0 +1,1384 @@
+#ifndef PROTECTED_SETS_H
+#define PROTECTED_SETS_H
+
+#include <algorithm>
+#include <CGAL/Cartesian_d.h>
+#include <CGAL/Epick_d.h>
+#include <CGAL/Euclidean_distance.h>
+#include <CGAL/Kernel_d/Sphere_d.h>
+#include <CGAL/Kernel_d/Hyperplane_d.h>
+#include <CGAL/Kernel_d/Vector_d.h>
+
+#include <CGAL/Orthogonal_k_neighbor_search.h>
+#include <CGAL/Kd_tree.h>
+#include <CGAL/Fuzzy_sphere.h>
+
+#include <boost/heap/fibonacci_heap.hpp>
+#include <boost/heap/policies.hpp>
+
+#include "output_tikz.h"
+#include "../output.h"
+#include "../generators.h"
+
+#include <CGAL/point_generators_d.h>
+
+
+typedef CGAL::Epick_d<CGAL::Dynamic_dimension_tag> K;
+typedef K::Point_d Point_d;
+typedef K::Line_d Line_d;
+typedef K::Vector_d Vector_d;
+typedef K::Oriented_side_d Oriented_side_d;
+typedef K::Has_on_positive_side_d Has_on_positive_side_d;
+typedef K::Sphere_d Sphere_d;
+typedef K::Hyperplane_d Hyperplane_d;
+
+typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation;
+typedef Delaunay_triangulation::Facet Facet;
+typedef Delaunay_triangulation::Vertex_handle Delaunay_vertex;
+typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle;
+
+typedef std::vector<Point_d> Point_Vector;
+typedef CGAL::Euclidean_distance<Traits_base> Euclidean_distance;
+
+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_k_neighbor_search<STraits, CGAL::Distance_adapter<std::ptrdiff_t,Point_d*,Euclidean_distance>> K_neighbor_search;
+typedef K_neighbor_search::Tree Tree;
+typedef K_neighbor_search::Distance Distance;
+typedef K_neighbor_search::iterator KNS_iterator;
+typedef K_neighbor_search::iterator KNS_range;
+typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
+
+typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator;
+
+
+FT  _sfty = pow(10,-14);
+
+bool experiment1, experiment2, experiment3, experiment5 = false;
+
+/* Experiment 1: epsilon as function on time **********************/
+std::vector<FT> eps_vector;
+
+/* Experiment 2: R/epsilon on alpha *******************************/
+std::vector<FT> epsratio_vector;
+std::vector<FT> epsslope_vector;
+
+/* Experiment 3: theta on delta ***********************************/
+std::vector<FT> thetamin_vector; FT curr_theta;
+std::vector<FT> gammamin_vector;
+
+/* Statistical data ***********************************************/
+int refused_case1, refused_case2, refused_bad, refused_centers1, refused_centers2;
+
+void initialize_statistics()
+{
+  refused_case1 = 0;
+  refused_case2 = 0;
+  refused_bad   = 0;
+  refused_centers1 = 0;
+  refused_centers2 = 0;
+}
+
+void print_statistics()
+{
+  std::cout << " * Old simplex not protected: " << refused_case1 << "\n";
+  std::cout << " * New simplex not protected: " << refused_case2 << "\n";
+  std::cout << " * New simplex not good:      " << refused_bad << "\n";
+  std::cout << " * New-old centers too close:      " << refused_centers1 << "\n";
+  std::cout << " * New-new centers too close:      " << refused_centers2 << "\n";
+}
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+// AUXILLARY FUNCTIONS
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Insert a point in Delaunay triangulation. If you are working in a flat torus, the procedure adds all the 3^d copies in adjacent cubes as well 
+ *  
+ *  W is the initial point vector
+ *  chosen_landmark is the index of the chosen point in W
+ *  landmarks_ind is the vector of indices of already chosen points in W
+ *  delaunay is the Delaunay triangulation
+ *  landmark_count is the current number of chosen vertices
+ *  torus is true iff you are working on a flat torus [-1,1]^d
+ *  OUT: Vertex handle to the newly inserted point
+ */
+Delaunay_vertex insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count, bool torus)
+{
+  if (!torus)
+    {
+      Delaunay_vertex v =delaunay.insert(W[chosen_landmark]);
+      landmarks_ind.push_back(chosen_landmark);
+      landmark_count++;
+      return v;
+    }
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      Delaunay_vertex v;
+      for (int i = 0; i < nb_cells; ++i)
+        {
+          std::vector<FT> point;
+          int cell_i = i;
+          for (int l = 0; l < D; ++l)
+            {
+              point.push_back(W[chosen_landmark][l] + 2.0*(cell_i%3-1));
+              cell_i /= 3;
+            }
+          if (i == nb_cells/2)
+            v = delaunay.insert(point); //v = center point
+          else
+            delaunay.insert(point);
+        }
+      landmarks_ind.push_back(chosen_landmark);
+      landmark_count++;
+      return v;
+    }
+}
+
+/** Small check if the vertex v is in the full cell fc
+ */
+
+bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc)
+{
+  for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it)
+    if (*v_it == v)
+      return true;
+  return false;
+}
+
+/** Fill chosen point vector from indices with copies if you are working on a flat torus
+ *  
+ *  IN:  W is the point vector
+ *  OUT: landmarks is the output vector
+ *  IN:  landmarks_ind is the vector of indices
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */
+
+void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind, bool torus)
+{
+  if (!torus)
+    for (unsigned j = 0; j < landmarks_ind.size(); ++j)
+      landmarks.push_back(W[landmarks_ind[j]]);
+  else
+    {
+      int D = W[0].size();
+      int nb_cells = pow(3, D);
+      int nbL = landmarks_ind.size();
+      // Fill landmarks
+      for (int i = 0; i < nb_cells-1; ++i)
+        for (int j = 0; j < nbL; ++j)
+          {
+            int cell_i = i;
+            Point_d point;
+            for (int l = 0; l < D; ++l)
+              {
+                point.push_back(W[landmarks_ind[j]][l] + 2.0*(cell_i-1));
+                cell_i /= 3;
+              }
+            landmarks.push_back(point);
+          }
+    }
+}
+
+/** Fill a vector of all simplices in the Delaunay triangulation giving integer indices to vertices
+ *
+ *  IN: t is the Delaunay triangulation
+ *  OUT: full_cells is the output vector
+ */
+
+void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells)
+{
+  // Store vertex indices in a map
+  int ind = 0; //index of a vertex
+  std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex;
+  for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+    if (t.is_infinite(v_it))
+      continue;
+    else
+      index_of_vertex[v_it] = ind++;
+  // Write full cells as vectors in full_cells
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      std::vector<int> cell;
+      for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it)
+        cell.push_back(index_of_vertex[*v_it]);
+      full_cells.push_back(cell);
+    }
+}
+
+bool sphere_intersects_cube(Point_d& c, FT r)
+{
+  bool in_cube = true;
+  //  int i = 0, D = p.size();
+  for (auto xi = c.cartesian_begin(); xi != c.cartesian_end(); ++xi)
+    // if ((*xi < 1.0 || *xi > -1.0) &&
+    //     (*xi-r < 1.0 || *xi-r > -1.0) &&
+    //     (*xi+r < 1.0 || *xi+r > -1.0))
+
+    if ((*xi-r < -1.0 && *xi+r < -1.0) ||
+        (*xi-r > 1.0  && *xi+r >  1.0 ))
+      {
+        in_cube = false; break;
+      }
+  return in_cube;
+}
+
+/** Recursive function for checking if the simplex is good,
+ *  meaning it does not contain a k-face, which is not theta0^(k-1) thick 
+ */
+
+bool is_theta0_good(std::vector<Point_d>& vertices, FT theta0)
+{
+  if (theta0 > 1)
+    {
+      std::cout << "Warning! theta0 is set > 1\n";
+      return false;
+    }
+  int D = vertices.size()-1;
+  if (D <= 1)
+    return true; // Edges are always good
+  //******** Circumscribed sphere
+  Euclidean_distance ed;
+  Sphere_d cs(vertices.begin(), vertices.end());
+  FT r = sqrt(cs.squared_radius());
+  for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it)
+    {
+      std::vector<Point_d> facet;
+      for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it)
+        if (f_it != v_it)
+          facet.push_back(*f_it);
+      // Compute the altitude
+
+      if (vertices[0].size() == 3 && D == 2)
+        {
+          //Vector_d l = facet[0] - facet[1];
+          FT orth_length2 = ed.transformed_distance(facet[0],facet[1]);
+          K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it;
+          FT h = 0;
+          // Scalar product = <sp,l>
+          FT scalar = 0;
+          for (p_it = v_it->cartesian_begin(),
+               s_it = facet[0].cartesian_begin(),
+               l_it = facet[1].cartesian_begin();
+               p_it != v_it->cartesian_end();
+               ++l_it, ++p_it, ++s_it)
+            scalar += (*l_it - *s_it)*(*p_it - *s_it);
+          // Gram-Schmidt for one vector
+          for (p_it = v_it->cartesian_begin(),
+               s_it = facet[0].cartesian_begin(),
+               l_it = facet[1].cartesian_begin();
+               p_it != v_it->cartesian_end();
+               ++l_it, ++p_it, ++s_it)
+            {
+              FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2;
+              h += hx*hx;
+            }
+          h = sqrt(h);
+      
+          if (h/(2*r) < pow(theta0, D-1))
+            return false;
+          if (!is_theta0_good(facet, theta0))
+            return false;
+        }
+      else
+        {
+          Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it);
+          Vector_d orth_tau = tau_h.orthogonal_vector();
+          FT orth_length = sqrt(orth_tau.squared_length());
+          K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it;
+          FT h = 0;
+          for (o_it = orth_tau.cartesian_begin(),
+                 p_it = v_it->cartesian_begin(),
+                 s_it = (facet.begin())->cartesian_begin();
+               o_it != orth_tau.cartesian_end();
+               ++o_it, ++p_it, ++s_it)
+            h += (*o_it)*(*p_it - *s_it)/orth_length;
+          h = fabs(h);
+          if (experiment3 && thetamin_vector[thetamin_vector.size()-1] > pow(h/(2*r), 1.0/(D-1)))
+            {
+              thetamin_vector[thetamin_vector.size()-1] = pow(h/(2*r), 1.0/(D-1));
+              //std::cout << "theta=" << h/(2*r) << ", ";
+            }
+          if (h/(2*r) < pow(theta0, D-1))
+              return false;
+          if (!is_theta0_good(facet, theta0))
+            return false;
+        }
+    }
+  return true;
+}
+
+/** Recursive function for checking the goodness of a simplex,
+ *  meaning it does not contain a k-face, which is not theta0^(k-1) thick 
+ */
+
+FT theta(std::vector<Point_d>& vertices)
+{
+  FT curr_value = 1.0;
+  int D = vertices.size()-1;
+  if (D <= 1)
+    return 1; // Edges are always good
+  //******** Circumscribed sphere
+  Euclidean_distance ed;
+  Sphere_d cs(vertices.begin(), vertices.end());
+  FT r = sqrt(cs.squared_radius());
+  for (std::vector<Point_d>::iterator v_it = vertices.begin(); v_it != vertices.end(); ++v_it)
+    {
+      std::vector<Point_d> facet;
+      for (std::vector<Point_d>::iterator f_it = vertices.begin(); f_it != vertices.end(); ++f_it)
+        if (f_it != v_it)
+          facet.push_back(*f_it);
+      // Compute the altitude
+      curr_value = std::min(curr_value, theta(facet)); // Check the corresponding facet
+      if (vertices[0].size() == 3 && D == 2)
+        {
+          //Vector_d l = facet[0] - facet[1];
+          FT orth_length2 = ed.transformed_distance(facet[0],facet[1]);
+          K::Cartesian_const_iterator_d l_it, p_it, s_it, c_it;
+          FT h = 0;
+          // Scalar product = <sp,l>
+          FT scalar = 0;
+          for (p_it = v_it->cartesian_begin(),
+               s_it = facet[0].cartesian_begin(),
+               l_it = facet[1].cartesian_begin();
+               p_it != v_it->cartesian_end();
+               ++l_it, ++p_it, ++s_it)
+            scalar += (*l_it - *s_it)*(*p_it - *s_it);
+          // Gram-Schmidt for one vector
+          for (p_it = v_it->cartesian_begin(),
+               s_it = facet[0].cartesian_begin(),
+               l_it = facet[1].cartesian_begin();
+               p_it != v_it->cartesian_end();
+               ++l_it, ++p_it, ++s_it)
+            {
+              FT hx = (*p_it - *s_it) - scalar*(*l_it - *s_it)/orth_length2;
+              h += hx*hx;
+            }
+          h = sqrt(h);
+          curr_value = std::min(curr_value, std::pow(h/(2*r), 1.0/(D-1)));
+        }
+      else
+        {
+          Hyperplane_d tau_h(facet.begin(), facet.end(), *v_it);
+          Vector_d orth_tau = tau_h.orthogonal_vector();
+          FT orth_length = sqrt(orth_tau.squared_length());
+          K::Cartesian_const_iterator_d o_it, p_it, s_it, c_it;
+          FT h = 0;
+          for (o_it = orth_tau.cartesian_begin(),
+                 p_it = v_it->cartesian_begin(),
+                 s_it = (facet.begin())->cartesian_begin();
+               o_it != orth_tau.cartesian_end();
+               ++o_it, ++p_it, ++s_it)
+            h += (*o_it)*(*p_it - *s_it)/orth_length;
+          h = fabs(h);
+          curr_value = std::min(curr_value, pow(h/(2*r), 1.0/(D-1)));
+        }
+    }
+  return curr_value;
+}
+
+// Doubling in a way 1->2->5->10
+void double_round(int& i)
+{
+  FT order10 = pow(10,std::floor(std::log10(i)));
+  int digit = std::floor( i / order10);
+  std::cout << digit;
+  if (digit == 1)
+    i *= 2;
+  else if (digit == 2)
+    i = 5*i/2;
+  else if (digit == 5)
+    i *= 2;
+  else
+    std::cout << "digit not correct. digit = " << digit << std::endl;
+}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// IS VIOLATED TEST
+////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+/** Check if a newly created cell is protected from old vertices
+ *  
+ *  t is the Delaunay triangulation
+ *  vertices is the vector containing the point to insert and a facet f in t
+ *  v1 is the vertex of t, such that f and v1 form a simplex
+ *  v2 is the vertex of t, such that f and v2 form another simplex
+ *  delta is the protection constant
+ *  power_protection is true iff the delta-power protection is used
+ */
+
+bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, const Delaunay_vertex& v1, const Delaunay_vertex v2, FT delta0, bool power_protection, FT theta0, FT gamma0)
+{
+  assert(vertices.size() == vertices[0].size() ||
+         vertices.size() == vertices[0].size() + 1); //simplex size = d | d+1
+  assert(v1 != v2);
+  if (vertices.size() == vertices[0].size() + 1)
+  // FINITE CASE
+    {
+      Sphere_d cs(vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0]));
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it)) 
+          {
+            //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); 
+            if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+              {
+                FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point());
+                if (!power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= (r+delta)*(r+delta))
+                    return true;
+                if (power_protection)
+                  if (dist2 >= r*r-_sfty && dist2 <= r*r+delta*delta) 
+                    return true;                 
+              }
+          }
+      */
+      // Is the center inside the box? (only Euclidean case)
+      // if (!torus)
+      //   {
+      //     bool inside_the_box = true;
+      //     for (c_it = center_cs.cartesian_begin(); c_it != center_cs.cartesian_end(); ++c_it)
+      //       if (*c_it > 1.0 || *c_it < -1.0)
+      //         {
+      //           inside_the_box = false; break;
+      //         }
+      //     if (inside_the_box && h/r < theta0)
+      //       return true;
+      //   }
+      // Check the two vertices (if not infinite)
+      if (!t.is_infinite(v1))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v1->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0))
+              { refused_case2++; return true;}
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+r*r*delta0*delta0) 
+              { refused_case2++; return true;}
+          // Check if the centers are not too close
+          std::vector<Point_d> sigma(vertices);
+          sigma[0] = v1->point();
+          Sphere_d cs_sigma(sigma.begin(), sigma.end());
+          Point_d csc_sigma = cs_sigma.center();
+          FT r_sigma = sqrt(cs_sigma.squared_radius());
+          FT dcc = sqrt(Euclidean_distance().transformed_distance(center_cs, csc_sigma));
+          if (experiment3 && dcc/r < gammamin_vector[gammamin_vector.size()-1])
+            gammamin_vector[gammamin_vector.size()-1] = dcc/r;
+          if (experiment3 && dcc/r_sigma < gammamin_vector[gammamin_vector.size()-1])
+            gammamin_vector[gammamin_vector.size()-1] = dcc/r_sigma;
+          if (dcc < r*gamma0 || dcc < r_sigma*gamma0)
+            { refused_centers1++; return true; }
+        }
+      if (!t.is_infinite(v2))
+        {
+          FT dist2 = Euclidean_distance().transformed_distance(center_cs, v2->point());
+          if (!power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0))
+              { refused_case2++; return true;}
+          if (power_protection)
+            if (dist2 >= r*r-_sfty && dist2 <= r*r+r*r*delta0*delta0) 
+              { refused_case2++; return true;}
+          // Check if the centers are not too close
+          std::vector<Point_d> sigma(vertices);
+          sigma[0] = v2->point();
+          Sphere_d cs_sigma(sigma.begin(), sigma.end());
+          Point_d csc_sigma = cs_sigma.center();
+          FT r_sigma = sqrt(cs_sigma.squared_radius());
+          FT dcc = sqrt(Euclidean_distance().transformed_distance(center_cs, csc_sigma));
+          if (experiment3 && dcc/r < gammamin_vector[gammamin_vector.size()-1])
+            gammamin_vector[gammamin_vector.size()-1] = dcc/r;
+          if (experiment3 && dcc/r_sigma < gammamin_vector[gammamin_vector.size()-1])
+            gammamin_vector[gammamin_vector.size()-1] = dcc/r_sigma;
+          if (dcc < r*gamma0 || dcc < r_sigma*gamma0)
+            { refused_centers1++; return true; }
+        }
+      // Check if the simplex is theta0-good
+      if (!is_theta0_good(vertices, theta0))
+        { refused_bad++; return true;}
+      
+    }
+  else
+  // INFINITE CASE
+    {
+      Delaunay_triangulation::Vertex_iterator v = t.vertices_begin();
+      while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end())
+        v++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE);
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      /*
+      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        if (!t.is_infinite(v_it))
+          if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
+            {
+              std::vector<FT> coords;
+              Point_d p = v_it->point();
+              auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+              for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+                coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
+              Point_d p_delta = Point_d(coords);
+              bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+              bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+              if (!p_is_inside && p_delta_is_inside)
+                return true;
+            }
+      */
+      if (!t.is_infinite(v1))
+        {
+          std::vector<FT> coords;
+          Point_d p = v1->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta0 / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+      if (!t.is_infinite(v2))
+        {
+          std::vector<FT> coords;
+          Point_d p = v2->point();
+          auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+          for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+            coords.push_back((*p_i) - (*orth_i) * delta0 / sqrt(orth_v.squared_length()));
+          Point_d p_delta = Point_d(coords);
+          bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
+          bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+          if (!power_protection && !p_is_inside && p_delta_is_inside)
+            return true;
+        }
+    }
+  return false;
+}
+  
+/** Auxillary recursive function to check if the point p violates the protection of the cell c and
+ *  if there is a violation of an eventual new cell
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  c is the current cell (simplex)
+ *  parent_cell is the parent cell (simplex)
+ *  index is the index of the facet between c and parent_cell from parent_cell's point of view
+ *  D is the dimension of the triangulation
+ *  delta is the protection constant
+ *  marked_cells is the vector of all visited cells containing p in their circumscribed ball
+ *  power_protection is true iff you are working with delta-power protection
+ *
+ *  OUT: true iff inserting p hasn't produced any violation so far
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta0, std::vector<Full_cell_handle>& marked_cells, bool power_protection, FT theta0, FT gamma0)
+{
+  Euclidean_distance ed;
+  std::vector<Point_d> vertices;
+  if (!t.is_infinite(c))
+    {
+      // if the cell is finite, we look if the protection is violated
+      for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it)
+        vertices.push_back((*v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d center_cs = cs.center();
+      FT r = sqrt(ed.transformed_distance(center_cs, vertices[0]));
+      FT dist2 = ed.transformed_distance(center_cs, p);
+      // if the new point is inside the protection ball of a non conflicting simplex
+      if (!power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= (r+r*delta0)*(r+r*delta0))
+          { refused_case1++; return true;}
+      if (power_protection)
+        if (dist2 >= r*r-_sfty && dist2 <= r*r+delta0*delta0*r*r) 
+          { refused_case1++; return true;}    
+      // if the new point is inside the circumscribing ball : continue violation searching on neighbours
+      //if (dist2 < r*r)
+      //if (dist2 < (5*r+delta)*(5*r+delta))
+      if (dist2 < r*r)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta0, marked_cells, power_protection, theta0, gamma0))
+                return true;
+            }
+        }
+      // if the new point is outside the protection sphere
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is guaranteed to be finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta0, power_protection, theta0, gamma0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;            
+        }
+    } 
+  else
+    {
+      // Inside of the convex hull is + side. Outside is - side.
+      for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+        if (!t.is_infinite(*vh_it))
+          vertices.push_back((*vh_it)->point());
+      Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin();
+      while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c))
+        v_it++;
+      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE);
+      //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point());
+      Vector_d orth_v = facet_plane.orthogonal_vector();
+      std::vector<FT> coords;
+      auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
+      for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
+        coords.push_back((*p_i) - (*orth_i) * delta0 / sqrt(orth_v.squared_length()));
+      Point_d p_delta = Point_d(coords);
+      bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p) && (Oriented_side_d()(facet_plane, p) != CGAL::ZERO);
+      bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
+
+      // If we work with power protection, we just ignore any conflicts
+      if (!power_protection && !p_is_inside && p_delta_is_inside)
+        return true;
+      //if the cell is infinite we look at the neighbours regardless
+      if (p_is_inside)
+        {
+          c->tds_data().mark_visited();
+          marked_cells.push_back(c);    
+          for (int i = 0; i < D+1; ++i)
+            {
+              Full_cell_handle next_c = c->neighbor(i);
+              if (next_c->tds_data().is_clear() &&
+                  is_violating_protection(p, t, next_c, c, i, D, delta0, marked_cells, power_protection, theta0, gamma0))
+                return true;
+            }
+        }
+      else
+        {
+          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
+          // the new simplex is finite if the parent cell is finite
+          vertices.clear(); vertices.push_back(p);
+          for (int i = 0; i < D+1; ++i)
+            if (i != index)
+              if (!t.is_infinite(parent_cell->vertex(i)))
+                vertices.push_back(parent_cell->vertex(i)->point());
+          Delaunay_vertex vertex_to_check = t.infinite_vertex();
+          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
+            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
+              {
+                vertex_to_check = *vh_it; break;
+              }
+          if (new_cell_is_violated(t, vertices, vertex_to_check, parent_cell->vertex(index), delta0, power_protection, theta0, gamma0)) 
+          //if (new_cell_is_violated(t, vertices, vertex_to_check->point(), delta)) 
+            return true;
+        }
+    }
+  //c->tds_data().clear_visited();
+  //marked_cells.pop_back();
+  return false;
+}
+
+/** Checks if inserting the point p in t will make conflicts
+ *
+ *  p is the point to insert
+ *  t is the current triangulation
+ *  D is the dimension of triangulation
+ *  delta is the protection constant
+ *  power_protection is true iff you are working with delta-power protection
+ *  OUT: true iff inserting p produces a violation of delta-protection.
+ */
+
+bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta0, bool power_protection, FT theta0, FT gamma0)
+{
+  Euclidean_distance ed;
+  Delaunay_triangulation::Vertex_handle v;
+  Delaunay_triangulation::Face f(t.current_dimension()); 
+  Delaunay_triangulation::Facet ft; 
+  Delaunay_triangulation::Full_cell_handle c; 
+  Delaunay_triangulation::Locate_type lt;
+  std::vector<Full_cell_handle> marked_cells;
+  //c = t.locate(p, lt, f, ft, v);
+  c = t.locate(p);
+  bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta0, marked_cells, power_protection, theta0, gamma0);
+  for (Full_cell_handle fc : marked_cells)
+    fc->tds_data().clear();
+  return violation_existing_cells;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// INITIALIZATION
+////////////////////////////////////////////////////////////////////////
+
+// Query for a sphere near a cite in all copies of a torus
+// OUT points_inside
+void torus_search(Tree& treeW, int D, Point_d cite, FT r, std::vector<int>& points_inside)
+{
+  int nb_cells = pow(3, D);
+  Delaunay_vertex v;
+  for (int i = 0; i < nb_cells; ++i)
+    {
+      std::vector<FT> cite_copy;
+      int cell_i = i;
+      for (int l = 0; l < D; ++l)
+        {
+          cite_copy.push_back(cite[l] + 2.0*(cell_i%3-1));
+          cell_i /= 3;
+        }
+      Fuzzy_sphere fs(cite_copy, r, 0, treeW.traits());
+      treeW.search(std::insert_iterator<std::vector<int>>(points_inside, points_inside.end()), fs);
+    }
+}
+
+  
+void initialize_torus(Point_Vector& W, Tree& treeW, Delaunay_triangulation& t, FT epsilon, std::vector<int>& landmarks_ind, int& landmark_count, std::vector<bool>& point_taken)
+{
+  initialize_statistics();
+  int D = W[0].size();
+  if (D == 2)
+    {
+      int xw = 6, yw = 4;
+      // Triangular lattice close to regular triangles h=0.866a ~ 0.875a   : 48p
+      for (int i = 0; i < xw; ++i)
+        for (int j = 0; j < yw; ++j)
+          {
+            Point_d cite1(std::vector<FT>{2.0/xw*i, 2.0/yw*j});
+            std::vector<int> points_inside;
+            torus_search(treeW, D, cite1, epsilon, points_inside);
+            //std::cout << "i=" << i << ", j=" << j << " "; print_vector(points_inside); std::cout << "\n";
+            std::vector<int>::iterator p_it = points_inside.begin();
+            while (p_it != points_inside.end() && point_taken[*p_it])
+              ++p_it;
+            assert(p_it != points_inside.end());
+            //W[*p_it] = cite1;  // debug purpose
+            insert_delaunay_landmark_with_copies(W, *p_it,
+                                                 landmarks_ind, t, landmark_count, true);
+            point_taken[*p_it] = true;
+
+            Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 2.0/yw*(j+0.5)});
+            points_inside.clear();
+            torus_search(treeW, D, cite2, epsilon, points_inside);
+            //std::cout << "i=" << i << ", j=" << j << " "; print_vector(points_inside); std::cout << "\n";
+            p_it = points_inside.begin();
+            while (p_it != points_inside.end() && point_taken[*p_it])
+              ++p_it;
+            assert(p_it != points_inside.end());
+            //W[*p_it] = cite2;  // debug purpose
+            insert_delaunay_landmark_with_copies(W, *p_it,
+                                                 landmarks_ind, t, landmark_count, true);
+            point_taken[*p_it] = true;            
+          }
+    }
+  else if (D == 3)
+    {
+      int wd = 3;
+      // Body-centered cubic lattice : 54p
+      for (int i = 0; i < wd; ++i)
+        for (int j = 0; j < wd; ++j)
+          for (int k = 0; k < wd; ++k)
+            {
+              Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k});
+              std::vector<int> points_inside;
+              torus_search(treeW, D, cite1, epsilon, points_inside);
+              std::vector<int>::iterator p_it = points_inside.begin();
+              while (p_it != points_inside.end() && point_taken[*p_it])
+                ++p_it;
+              assert(p_it != points_inside.end());
+              insert_delaunay_landmark_with_copies(W, *(points_inside.begin()),
+                                                   landmarks_ind, t, landmark_count, true);
+              point_taken[*p_it] = true;
+              
+              Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)});
+              points_inside.clear();
+              torus_search(treeW, D, cite2, epsilon, points_inside);
+              p_it = points_inside.begin();
+              while (p_it != points_inside.end() && point_taken[*p_it])
+                ++p_it;
+              assert(p_it != points_inside.end());
+              insert_delaunay_landmark_with_copies(W, *(points_inside.begin()),
+                                                   landmarks_ind, t, landmark_count, true);
+              point_taken[*p_it] = true;
+            }
+    }
+  //write_mesh
+}
+  
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+//!!!!!!!!!!!!! THE INTERFACE FOR LANDMARK CHOICE IS BELOW !!!!!!!!!!//
+///////////////////////////////////////////////////////////////////////
+///////////////////////////////////////////////////////////////////////
+
+// Struct for R_max_heap elements
+
+struct R_max_handle
+{
+  FT value;
+  Point_d center;
+  
+  R_max_handle(FT value_, Point_d c): value(value_), center(c)
+  {}
+};
+
+struct R_max_compare
+{
+  bool operator()(const R_max_handle& rmh1, const R_max_handle& rmh2) const
+  {
+    return rmh1.value < rmh2.value;
+  }
+};
+
+// typedef boost::heap::fibonacci_heap<R_max_handle, boost::heap::compare<R_max_compare>> Heap;
+
+// void make_heap(Delaunay_triangulation& t, Heap& R_max_heap)
+// {
+//   R_max_heap.clear();
+//   for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it) 
+//     {
+//       if (t.is_infinite(fc_it))
+//         continue;
+//       Point_Vector vertices;
+//       for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+//         vertices.push_back((*fc_v_it)->point());
+//       Sphere_d cs( vertices.begin(), vertices.end());
+//       Point_d csc = cs.center();
+//       FT r = sqrt(cs.squared_radius());
+//       // A ball is in the heap, if it intersects the cube
+//       bool accepted = sphere_intersects_cube(csc, sqrt(r)); 
+//       if (!accepted)
+//         continue;
+//       R_max_heap.push(R_max_handle(r, fc_it, csc));
+//     }
+// }
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+// SAMPLING RADIUS 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+R_max_handle sampling_radius(Delaunay_triangulation& t)
+{
+  FT epsilon2 = 0;
+  Point_d final_center;
+  Point_d control_point;
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin()));
+      if (epsilon2 < r2)
+        {
+          epsilon2 = r2;
+          final_center = csc;         
+          control_point = (*vertices.begin());
+        }
+    }
+  return R_max_handle(sqrt(epsilon2), final_center);
+}
+
+FT sampling_fatness(Delaunay_triangulation& t)
+{
+  FT curr_theta = 1.0;
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
+    {
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      FT theta_f = theta(vertices);
+      curr_theta = std::min(curr_theta, theta_f);
+      //std::cout << "theta(sigma) = " << theta_f << "\n";
+    }
+  return curr_theta;
+}
+
+// Generate an epsilon sample for a given epsilon
+void generate_epsilon_sample_torus(Point_Vector& W, FT epsilon, int dim, Delaunay_triangulation& t)
+{
+  W.clear();
+  t.clear();
+  int point_count = 0;
+  std::vector<int> point_ind;
+  // std::vector<FT> coords;
+  FT curr_eps = 2*dim;
+  // Initialize
+  // for (int i = 0; i < dim; ++i)
+  //   coords.push_back(-1);
+  // R_max_handle rmh(2*sqrt(dim), Point_d(coords));
+  // int N = dim; std::floor(std::pow(1/epsilon,dim));
+  // std::cout << N << "\n";
+  typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator;
+  Random_cube_iterator rp(dim, 1.0);
+  W.push_back(*rp++);
+  insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true);
+  curr_eps = sampling_radius(t).value;      
+  while (curr_eps > epsilon)
+    {
+      
+      W.push_back(*rp++);
+      insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true);
+      
+      Point_d c = sampling_radius(t).center;
+      W.push_back(c);
+      insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true);
+      curr_eps = sampling_radius(t).value;
+      
+      std::cout << "curr_eps = " << curr_eps << "\n";
+    }
+  // Iterate and insert in a torus
+  // while (rmh.value > epsilon)
+  //   {
+  //     W.push_back(rmh.center);
+  //     insert_delaunay_landmark_with_copies(W, W.size()-1, point_ind, t, point_count, true);
+  //     rmh = sampling_radius(t);
+  //     //std::cout << rmh.value;
+  //   }
+}
+
+///////////////////////////////////////////////////////////////////////
+// LANDMARK CHOICE PROCEDURE
+///////////////////////////////////////////////////////////////////////
+
+/** Procedure to compute a maximal protected subset from a point cloud. All OUTs should be empty at call.
+ *
+ *  IN:  W is the initial point cloud having type Epick_d<Dynamic_dimension_tag>::Point_d
+ *  IN:  nbP is the size of W
+ *  OUT: landmarks is the output vector for the points
+ *  OUT: landmarks_ind is the output vector for the indices of the selected points in W
+ *  IN:  delta is the constant of protection
+ *  OUT: full_cells is the output vector of the simplices in the final Delaunay triangulation
+ *  IN:  torus is true iff you are working on a flat torus [-1,1]^d
+ */ 
+
+void protected_delaunay(Point_Vector& W,
+                        //Point_Vector& landmarks,
+                        std::vector<int>& landmarks_ind,
+                        FT alpha,
+                        FT epsilon,
+                        FT delta0,
+                        FT theta0,
+                        FT gamma0,
+                        //std::vector<std::vector<int>>& full_cells,
+                        bool torus,
+                        bool power_protection
+                        )
+{
+  //bool return_ = true;
+  unsigned D = W[0].size();
+  int nbP = W.size();
+  //FT beta = 1/(1-alpha);
+  //FT Ad = pow((4*alpha + 8*beta)/alpha, D);
+  //FT theta0 = 1/Ad;
+  //FT delta0 = pow(1/Ad,D);
+  Torus_distance td;
+  Euclidean_distance ed;
+  Delaunay_triangulation t(D);
+  std::vector<bool> point_taken(nbP,false);
+  CGAL::Random rand;
+  int landmark_count = 0;
+  std::list<int> index_list;
+  //****************** Kd Tree W
+  STraits traits(&(W[0]));
+  Tree treeW(boost::counting_iterator<std::ptrdiff_t>(0),
+         boost::counting_iterator<std::ptrdiff_t>(nbP),
+         typename Tree::Splitter(),
+         traits);
+  //      shuffle the list of indexes (via a vector)
+  {
+    std::vector<int> temp_vector;
+    for (int i = 0; i < nbP; ++i)
+      temp_vector.push_back(i);
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed));
+    //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end());
+    for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it)
+      index_list.push_front(*it);
+  }
+  //******************** Initialize point set
+  if (!torus)
+    for (unsigned pos1 = 0; pos1 < D+1; ++pos1)
+      {
+        std::vector<FT> point;
+        for (unsigned i = 0; i < pos1; ++i)
+          point.push_back(-1);
+        if (pos1 != D)
+          point.push_back(1);
+        for (unsigned i = pos1+1; i < D; ++i)
+          point.push_back(0);
+        assert(point.size() == D);
+        W[index_list.front()] = Point_d(point);
+        insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count, torus);
+        index_list.pop_front();
+      }
+  else
+    initialize_torus(W, treeW, t, epsilon, landmarks_ind, landmark_count, point_taken);
+  //std::cout << "Size of treeW: " << treeW.size() << "\n";
+  //std::cout << "Size of t: " << t.number_of_vertices() << "\n";
+  //******************* Initialize heap for R_max
+  //Heap R_max_heap;
+  //make_heap(t, R_max_heap);
+
+  
+  R_max_handle rh = sampling_radius(t);
+  FT epsilon0 = rh.value; 
+  if (experiment1) eps_vector.push_back(pow(1/rh.value,D));
+  //******************** Iterative algorithm
+  std::vector<int> candidate_points;
+  torus_search(treeW, D,
+               rh.center,
+               alpha*rh.value,
+               candidate_points);
+  std::list<int>::iterator list_it;
+  std::vector<int>::iterator cp_it = candidate_points.begin();
+  while (cp_it != candidate_points.end())
+    {
+      if (!point_taken[*cp_it] && !is_violating_protection(W[*cp_it], t, D, delta0, power_protection, theta0, gamma0))
+        {
+          Delaunay_vertex v = insert_delaunay_landmark_with_copies(W, *cp_it, landmarks_ind, t, landmark_count, torus);
+          {
+            // Simple check if the new cells don't have centers too close one to another
+            std::vector<Full_cell_handle> inc_cells;
+            std::back_insert_iterator<std::vector<Full_cell_handle>> out(inc_cells);
+            t.tds().incident_full_cells(v, out);
+
+            std::vector<Sphere_d> spheres;
+            for (auto i_it = inc_cells.begin(); i_it != inc_cells.end(); ++i_it)
+              {
+                std::vector<Point_d> vertices;
+                for (auto v_it = (*i_it)->vertices_begin(); v_it != (*i_it)->vertices_end(); ++v_it)
+                  vertices.push_back((*v_it)->point());
+                spheres.push_back(Sphere_d(vertices.begin(), vertices.end()));
+              }
+            for (auto s_it = spheres.begin(); s_it != spheres.end(); ++s_it)
+              for (auto t_it = s_it+1; t_it != spheres.end(); ++t_it)
+                {
+                  FT ddc2 = ed.transformed_distance(s_it->center(),t_it->center());
+                  if (ddc2 < gamma0*gamma0*s_it->squared_radius() ||
+                      ddc2 < gamma0*gamma0*t_it->squared_radius())
+                    { refused_centers2++; }
+                }
+          }
+            
+          //std::cout << *cp_it << ",\n";
+          //make_heap(t, R_max_heap);
+          point_taken[*cp_it] = true;
+          rh = sampling_radius(t);
+          if (experiment1) eps_vector.push_back(pow(1/rh.value,D));
+          //std::cout << "rhvalue = " << rh.value << "\n";
+          //std::cout << "D = " << 
+          candidate_points.clear();
+          torus_search(treeW, D,
+                       rh.center,
+                       alpha*rh.value,
+                       candidate_points);
+          cp_it = candidate_points.begin();
+          /*
+          // PIECE OF CODE FOR DEBUGGING PURPOSES
+          
+          Delaunay_vertex inserted_v = insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count);
+          if (triangulation_is_protected(t, delta))
+            {
+              index_list.erase(list_it);
+              list_it = index_list.begin();
+            }
+          else
+            { //THAT'S WHERE SOMETHING'S WRONG
+              t.remove(inserted_v);
+              landmarks_ind.pop_back();
+              landmark_count--;
+              write_delaunay_mesh(t, W[*list_it], is2d);
+              is_violating_protection(W[*list_it], t_old, D, delta); //Called for encore
+            }
+          */
+          //std::cout << "index_list_size() = " << index_list.size() << "\n";
+        }
+      else
+        {
+          cp_it++;
+          //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n";
+        }
+      //if (list_it != index_list.end())
+      //  write_delaunay_mesh(t, W[*list_it], is2d);
+    }
+  
+  if (experiment2) epsratio_vector.push_back(rh.value/epsilon0);
+  if (experiment2) epsslope_vector.push_back( (pow(1/rh.value,D)-pow(1/epsilon0,D))/(landmarks_ind.size() - 48) );
+  std::cout << "The iteration ended when cp_count = " << candidate_points.size() << "\n";
+  std::cout << "alphaRmax = " << alpha*rh.value << "\n";
+  std::cout << "epsilon' = " << rh.value << "\n";
+  std::cout << "nbL = " << landmarks_ind.size() << "\n";
+  print_statistics();
+  //print_vector(landmarks_ind); std::cout << std::endl;
+  //std::sort(landmarks_ind.begin(), landmarks_ind.end());
+  print_vector(landmarks_ind); std::cout << std::endl;
+  if (experiment3) thetamin_vector[thetamin_vector.size()-1] = sampling_fatness(t);
+  std::cout << "theta = " << sampling_fatness(t) << "\n";
+  //fill_landmarks(W, landmarks, landmarks_ind, torus);
+  //fill_full_cell_vector(t, full_cells);
+  /*
+  if (triangulation_is_protected(t, delta))
+    std::cout << "Triangulation is ok\n";
+  else
+    {
+      std::cout << "Triangulation is BAD!! T_T しくしく!\n";
+    }
+  */
+  write_delaunay_mesh(t, W[0], true);
+  //std::cout << t << std::endl;
+}
+
+void run_experiment5(Point_Vector& W,
+                     int D,
+                     FT alpha,
+                     FT epsilon,
+                     FT delta0,
+                     FT theta0,
+                     FT gamma0,
+                     //std::vector<std::vector<int>>& full_cells,
+                     bool torus,
+                     bool power_protection
+                     )
+{
+  // INITIALIZATION
+  Delaunay_triangulation t(D);
+  std::vector<int> landmarks_ind;
+  int landmark_count = 0;
+  initialize_statistics();
+  if (D == 2)
+    {
+      int xw = 6, yw = 4;
+      // Triangular lattice close to regular triangles h=0.866a ~ 0.875a   : 48p
+      for (int i = 0; i < xw; ++i)
+        for (int j = 0; j < yw; ++j)
+          {
+            Point_d cite1(std::vector<FT>{2.0/xw*i, 2.0/yw*j});
+            W.push_back(cite1);  // debug purpose
+            insert_delaunay_landmark_with_copies(W, W.size()-1,
+                                                 landmarks_ind, t, landmark_count, true);
+
+            Point_d cite2(std::vector<FT>{2.0/xw*(i+0.5), 2.0/yw*(j+0.5)});
+            W.push_back(cite2);  // debug purpose
+            insert_delaunay_landmark_with_copies(W, W.size()-1,
+                                                 landmarks_ind, t, landmark_count, true);
+          }
+    }
+  else if (D == 3)
+    {
+      int wd = 3;
+      // Body-centered cubic lattice : 54p
+      for (int i = 0; i < wd; ++i)
+        for (int j = 0; j < wd; ++j)
+          for (int k = 0; k < wd; ++k)
+            {
+              Point_d cite1(std::vector<FT>{2.0/wd*i, 2.0/wd*j, 2.0/wd*k});
+              W.push_back(cite1);  // debug purpose
+              insert_delaunay_landmark_with_copies(W, W.size()-1,
+                                                   landmarks_ind, t, landmark_count, true);
+
+              Point_d cite2(std::vector<FT>{2.0/wd*(i+0.5), 2.0/wd*(j+0.5), 2.0/wd*(k+0.5)});
+              W.push_back(cite2);  // debug purpose
+              insert_delaunay_landmark_with_copies(W, W.size()-1,
+                                                   landmarks_ind, t, landmark_count, true);
+            }
+    }
+
+  // ITERATIONS
+  R_max_handle rh = sampling_radius(t);
+  Point_d rp = *(Random_point_iterator(D, alpha*rh.value));
+  int death_count = 0;
+  std::cout << "death count " << death_count << " rp = " << rp << "\n";
+  while (death_count < 100)
+    {
+      std::vector<FT> coords;
+      for (auto c_it = rh.center.cartesian_begin(),
+           r_it = rp.cartesian_begin();
+           c_it != rh.center.cartesian_end();
+           ++c_it, ++r_it)
+        coords.push_back(*c_it + *r_it);
+      Point_d new_p(coords);
+      if (!is_violating_protection(new_p, t, D, delta0, power_protection, theta0, gamma0))
+        {
+          W.push_back(new_p);
+          insert_delaunay_landmark_with_copies(W, W.size()-1, landmarks_ind, t, landmark_count, torus);
+          rh = sampling_radius(t);
+          rp = *(Random_point_iterator(D, alpha*rh.value));
+          death_count = 0;
+          std::cout << "death count " << death_count << " rp = " << rp << "\n";
+        }
+      else
+        {
+          rp = *(Random_point_iterator(D, alpha*rh.value));
+          death_count++;
+          std::cout << "death count " << death_count << " rp = " << rp << "\n";
+        }
+      //Point_d new_p = (*rp++) + Vector_d;
+    }
+}
+
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Series of experiments
+///////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+void start_experiments(Point_Vector& W, FT alpha, std::vector<int>& landmarks_ind, FT epsilon)
+{
+  int experiment_no = 1;
+  FT delta0 = 0.1;
+  FT theta0 = 0.1;
+  FT gamma0 = 0.01;
+  std::string suffix;
+  //std::cout << "ようこそジプシー我が神秘の部屋へ:\n";
+  while (experiment_no != 0)
+    {
+      std::cout << "Enter experiment no (0 to exit): ";
+      std::cin >> experiment_no;
+      switch (experiment_no)
+        {
+        case 1:
+          // Experiment 1
+          experiment1 = true;         
+          eps_vector = {};
+          std::cout << "Enter delta0: "; std::cin >> delta0;
+          std::cout << "Enter theta0: "; std::cin >> theta0;
+          std::cout << "Enter gamma0: "; std::cin >> gamma0;
+          protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true);
+          write_tikz_plot(eps_vector,"epstime.tikz");
+          experiment1 = false;
+          break;
+          
+        case 2:
+          // Experiment 2
+          suffix = "";
+          experiment2 = true;
+          epsratio_vector = {0};
+          epsslope_vector = {0};
+          std::cout << "File name suffix: ";
+          std::cin >> suffix;
+          for (FT alpha = 0.01; alpha < 0.999; alpha += 0.01)
+            {
+              landmarks_ind.clear();
+              std::cout << "Test for alpha = " << alpha << "\n";
+              protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true);
+            }
+          write_tikz_plot(epsratio_vector,"epsratio_alpha." + suffix + ".tex");
+          write_tikz_plot(epsslope_vector,"epsslope_alpha." + suffix + ".tex");
+          experiment2 = false;
+          break;
+
+        case 3:
+          // Experiment 3
+          experiment3 = true;
+          thetamin_vector = {};
+          gammamin_vector = {};
+          theta0 = 0;
+          gamma0 = 0;
+          for (FT delta0 = 0; delta0 < 0.999; delta0 += 0.05)
+            {
+              landmarks_ind.clear();
+              thetamin_vector.push_back(1.0); //0.7489 fatness of the initialization
+              gammamin_vector.push_back(10);
+              std::cout << "Test for delta0 = " << delta0 << "\n";
+              protected_delaunay(W, landmarks_ind, alpha, epsilon, delta0, theta0, gamma0, true, true);
+            }
+          write_tikz_plot(thetamin_vector,"thetamin_delta.tex");
+          write_tikz_plot(gammamin_vector,"gammamin_delta.tex");
+          experiment3 = false;
+          break;
+
+        // case 4:
+        //   // Experiment 4
+        //   {
+        //     int dim;
+        //     std::cout << "Enter dimension: ";
+        //     std::cin >> dim;
+        //     Delaunay_triangulation t(dim);
+        //     // for (FT eps = 0.7; eps < 1.1; eps += 0.1)
+        //   //     {
+        //   //       generate_epsilon_sample_torus(W, eps, dim, t);
+        //   //       for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+        //   //         {
+        //   //           if (t.is_infinite(v_it))
+        //   //             continue;
+        //   //           bool in_cube = true; 
+        //   //           for (auto xi = v_it->cartesian_begin(); xi != v_it->cartesian_end(); ++xi)
+        //   //             if (*xi > 1.0 || *xi < -1.0)
+        //   //               {
+        //   //                 in_cube = false; break;
+        //   //               }
+        //   //           if (!in_cube)
+        //   //             continue;
+        //   //           for (auto t.tds().incident_full_cells())
+        //   //         }
+        //   //       std::cout << "eps = " << eps << ", real epsilon = " << sampling_radius(t).value << "\n";
+        //   //     }
+        //   // }
+        //   break;
+
+
+        case 5:
+          // Experiment 5
+          experiment5 = true;
+          //     std::cout << "Enter dimension: ";
+          //     std::cin >> dim;
+
+          landmarks_ind.clear();
+          W.clear();
+          run_experiment5(W,  alpha, epsilon, delta0, theta0, gamma0, true, true);
+          experiment5 = false;
+          break;
+        }
+
+    }
+
+}
+  
+#endif
diff --git a/src/Witness_complex/example/witness_complex_cube.cpp b/src/Witness_complex/example/witness_complex_cube.cpp
index a9a2959b..e448c55d 100644
--- a/src/Witness_complex/example/witness_complex_cube.cpp
+++ b/src/Witness_complex/example/witness_complex_cube.cpp
@@ -20,6 +20,11 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
+// Avoiding the max arity issue with CGAL
+#ifndef BOOST_PARAMETER_MAX_ARITY
+# define BOOST_PARAMETER_MAX_ARITY 12
+#endif
+
 #include <iostream>
 #include <fstream>
 #include <ctime>
@@ -37,6 +42,10 @@
 #include "gudhi/Witness_complex.h"
 #include "gudhi/reader_utils.h"
 #include "Torus_distance.h" 
+#include "generators.h"
+#include "output.h"
+//#include "protected_sets/protected_sets.h"
+#include "protected_sets/protected_sets_paper2.h"
 
 #include <CGAL/Cartesian_d.h>
 #include <CGAL/Search_traits.h>
@@ -106,8 +115,6 @@ typedef std::vector<Point_d> Point_Vector;
 
 //typedef K::Equal_d Equal_d;
 //typedef CGAL::Random_points_in_cube_d<CGAL::Point_d<CGAL::Cartesian_d<FT> > > Random_cube_iterator;
-typedef CGAL::Random_points_in_cube_d<Point_d> Random_cube_iterator;
-typedef CGAL::Random_points_in_ball_d<Point_d> Random_point_iterator;
 
 typedef CGAL::Delaunay_triangulation<K> Delaunay_triangulation;
 typedef Delaunay_triangulation::Facet Facet;
@@ -117,449 +124,84 @@ typedef Delaunay_triangulation::Full_cell_handle Full_cell_handle;
 typedef K::Sphere_d Sphere_d;
 typedef K::Hyperplane_d Hyperplane_d;
 
+/*//////////////////////////////////////
+ * GLOBAL VARIABLES ********************
+ *//////////////////////////////////////
 
-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_grid(Point_Vector& W, int width, int D)
-{
-  int nb_points = 1;
-  for (int i = 0; i < D; ++i)
-    nb_points *= width;
-  for (int i = 0; i < nb_points; ++i)
-    {
-      std::vector<double> point;
-      int cell_i = i;
-      for (int l = 0; l < D; ++l)
-        {
-          point.push_back(0.01*(cell_i%width));
-          cell_i /= width;
-        }
-      W.push_back(point);
-    }
-}
-
-void generate_points_random_box(Point_Vector& W, int nbP, int dim)
-{
-  /*
-  Random_cube_iterator rp(dim, 1);
-  for (int i = 0; i < nbP; i++)
-    {
-      std::vector<double> point;
-      for (auto it = rp->cartesian_begin(); it != rp->cartesian_end(); ++it)
-        point.push_back(*it);
-      W.push_back(Point_d(point));
-      rp++;
-    }
-  */
-  Random_cube_iterator rp(dim, 1.0);
-  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();
-}
+//NA bool toric=false;
+bool power_protection = true;
+bool grid_points = true;
+bool is2d = true;
+//FT  _sfty = pow(10,-14);
+bool torus = false;
 
 
-void write_points( std::string file_name, std::vector< Point_d > & points)
-{
-  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(std::string file_name, Witness_complex<>& witness_complex, Point_Vector& landmarks)
-{
-  std::ofstream ofs (file_name, std::ofstream::out);
-  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();
-}
-
-
-void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark, std::vector<int>& landmarks_ind, Delaunay_triangulation& delaunay, int& landmark_count)
-{
-  delaunay.insert(W[chosen_landmark]);
-  landmarks_ind.push_back(chosen_landmark);
-  landmark_count++;
-}
-
-bool vertex_is_in_full_cell(Delaunay_triangulation::Vertex_handle v, Full_cell_handle fc)
-{
-  for (auto v_it = fc->vertices_begin(); v_it != fc->vertices_end(); ++v_it)
-    if (*v_it == v)
-      return true;
-  return false;
-}
-
-bool new_cell_is_violated(Delaunay_triangulation& t, std::vector<Point_d>& vertices, bool is_infinite, const Point_d& p, FT delta)
-{
-  if (!is_infinite)
-  // FINITE CASE
-    {
-      Sphere_d cs(vertices.begin(), vertices.end());
-      Point_d center_cs = cs.center();
-      FT r = sqrt(Euclidean_distance().transformed_distance(center_cs, vertices[0]));
-      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
-        if (!t.is_infinite(v_it)) 
-          {
-            //CGAL::Oriented_side side = Oriented_side_d()(cs, (v_it)->point()); 
-            if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
-              {
-                FT dist2 = Euclidean_distance().transformed_distance(center_cs, (v_it)->point());
-                //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta))
-                if (dist2 >= r*r && dist2 <= r*r+delta*delta) 
-                  return true;                 
-              }
-          }
-    }
-  else
-  // INFINITE CASE
-    {
-      Delaunay_triangulation::Vertex_iterator v = t.vertices_begin();
-      while (t.is_infinite(v) || std::find(vertices.begin(), vertices.end(), v->point()) == vertices.end())
-        v++;
-      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v->point(), CGAL::ON_POSITIVE_SIDE);
-      Vector_d orth_v = facet_plane.orthogonal_vector();
-      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
-        if (!t.is_infinite(v_it))
-          if (std::find(vertices.begin(), vertices.end(), v_it->point()) == vertices.end())
-            {
-              std::vector<FT> coords;
-              Point_d p = v_it->point();
-              auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
-              for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
-                coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
-              Point_d p_delta = Point_d(coords);
-              bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
-              bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
-              if (!p_is_inside && p_delta_is_inside)
-                return true;
-            }
-    }
-  return false;
-}
-  
-
-bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, Full_cell_handle c, Full_cell_handle parent_cell, int index, int D, FT delta, std::vector<Full_cell_handle>& marked_cells)
-{
-  Euclidean_distance ed;
-  std::vector<Point_d> vertices;
-  if (!t.is_infinite(c))
-    {
-      // if the cell is finite, we look if the protection is violated
-      for (auto v_it = c->vertices_begin(); v_it != c->vertices_end(); ++v_it)
-        vertices.push_back((*v_it)->point());
-      Sphere_d cs( vertices.begin(), vertices.end());
-      Point_d center_cs = cs.center();
-      FT r = sqrt(ed.transformed_distance(center_cs, vertices[0]));
-      FT dist2 = ed.transformed_distance(center_cs, p);
-      // if the new point is inside the protection ball of a non conflicting simplex
-      //if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta))
-      if (dist2 >= r*r && dist2 <= r*r+delta*delta) 
-        return true;    
-      c->tds_data().mark_visited();
-      marked_cells.push_back(c);
-      // if the new point is inside the circumscribing ball : continue violation searching on neughbours
-      if (dist2 < r*r)
-        for (int i = 0; i < D+1; ++i)
-          {
-            Full_cell_handle next_c = c->neighbor(i);
-            if (next_c->tds_data().is_clear() &&
-                is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells))
-              return true;
-          }
-      // if the new point is outside the protection sphere
-      else
-        {
-          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
-          // the new simplex is guaranteed to be finite
-          vertices.clear(); vertices.push_back(p);
-          for (int i = 0; i < D+1; ++i)
-            if (i != index)
-              vertices.push_back(parent_cell->vertex(i)->point());
-          Delaunay_vertex vertex_to_check;
-          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
-            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
-              {
-                vertex_to_check = *vh_it; break;
-              }
-          if (new_cell_is_violated(t, vertices, false, vertex_to_check->point(), delta)) 
-            return true;            
-        }
-    } 
-  else
-    {
-      // Inside of the convex hull is + side. Outside is - side.
-      for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
-        if (!t.is_infinite(*vh_it))
-          vertices.push_back((*vh_it)->point());
-      Delaunay_triangulation::Vertex_iterator v_it = t.vertices_begin();
-      while (t.is_infinite(v_it) || vertex_is_in_full_cell(v_it, c))
-        v_it++;
-      Hyperplane_d facet_plane(vertices.begin(), vertices.end(), v_it->point(), CGAL::ON_POSITIVE_SIDE);
-      //CGAL::Oriented_side outside = Oriented_side_d()(facet_plane, v_it->point());
-      Vector_d orth_v = facet_plane.orthogonal_vector();
-      std::vector<FT> coords;
-      auto orth_i = orth_v.cartesian_begin(), p_i = p.cartesian_begin();
-      for (; orth_i != orth_v.cartesian_end(); ++orth_i, ++p_i)
-        coords.push_back((*p_i) - (*orth_i) * delta / sqrt(orth_v.squared_length()));
-      Point_d p_delta = Point_d(coords);
-      bool p_is_inside = !Has_on_positive_side_d()(facet_plane, p);
-      bool p_delta_is_inside = !Has_on_positive_side_d()(facet_plane, p_delta);
-
-      if (!p_is_inside && p_delta_is_inside)
-        return true;
-      //if the cell is infinite we look at the neighbours regardless
-      c->tds_data().mark_visited();
-      marked_cells.push_back(c);
-      if (p_is_inside)
-        for (int i = 0; i < D+1; ++i)
-          {
-            Full_cell_handle next_c = c->neighbor(i);
-            if (next_c->tds_data().is_clear() &&
-                is_violating_protection(p, t, next_c, c, i, D, delta, marked_cells))
-              return true;
-          }
-      else
-        {
-          // facet f is on the border of the conflict zone : check protection of simplex {p,f}
-          // the new simplex is finite if the parent cell is finite
-          vertices.clear(); vertices.push_back(p);
-          bool new_simplex_is_finite = false;
-          for (int i = 0; i < D+1; ++i)
-            if (i != index)
-              {
-                if (t.is_infinite(parent_cell->vertex(i)))
-                  new_simplex_is_finite = true;
-                else
-                  vertices.push_back(parent_cell->vertex(i)->point());
-              }
-          Delaunay_vertex vertex_to_check;
-          for (auto vh_it = c->vertices_begin(); vh_it != c->vertices_end(); ++vh_it)
-            if (!vertex_is_in_full_cell(*vh_it, parent_cell))
-              {
-                vertex_to_check = *vh_it; break;
-              }
-          if (new_cell_is_violated(t, vertices, new_simplex_is_finite, vertex_to_check->point(), delta)) 
-            return true;
-        }
-    }
-  return false;
-}
-
-bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta)
-{
-  Euclidean_distance ed;
-  Delaunay_triangulation::Vertex_handle v;
-  Delaunay_triangulation::Face f(t.current_dimension()); 
-  Delaunay_triangulation::Facet ft; 
-  Delaunay_triangulation::Full_cell_handle c; 
-  Delaunay_triangulation::Locate_type lt;
-  std::vector<Full_cell_handle> marked_cells;
-  c = t.locate(p, lt, f, ft, v);
-  bool violation_existing_cells = is_violating_protection(p, t, c, c, 0, D, delta, marked_cells);
-  for (Full_cell_handle fc : marked_cells)
-    fc->tds_data().clear();
-  return violation_existing_cells;
-}
-
-bool old_is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta)
+bool triangulation_is_protected(Delaunay_triangulation& t, FT delta)
 {
+  std::cout << "Start protection verification\n";
   Euclidean_distance ed;
-  Delaunay_triangulation::Vertex_handle v;
-  Delaunay_triangulation::Face f(t.current_dimension()); 
-  Delaunay_triangulation::Facet ft; 
-  Delaunay_triangulation::Full_cell_handle c; 
-  Delaunay_triangulation::Locate_type lt;
-  c = t.locate(p, lt, f, ft, v);
-  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
-    if (!t.is_infinite(fc_it))
-      {
-        std::vector<Point_d> vertices;
-        for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it)
-          vertices.push_back((*v_it)->point());
-        Sphere_d cs( vertices.begin(), vertices.end());
-        Point_d center_cs = cs.center();
-        FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(1)->point()));
-        FT dist2 = ed.transformed_distance(center_cs, p);
-        //if the new point is inside the protection ball of a non conflicting simplex
-        if (dist2 >= r*r && dist2 <= (r+delta)*(r+delta)) 
-          return true;
-      }
-  t.insert(p, c);
-  return false;
-}
-
-void write_delaunay_mesh(Delaunay_triangulation& t, const Point_d& p)
-{
-  std::ofstream ofs ("delaunay.mesh", std::ofstream::out);
-  int nbV = t.number_of_vertices()+1;
-  ofs << "MeshVersionFormatted 1\nDimension 2\n";
-  ofs << "Vertices\n" << nbV << "\n";
-  int ind = 1; //index of a vertex
+  // Fill the map Vertices -> Numbers
   std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex;
+  int ind = 0;
   for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
     {
       if (t.is_infinite(v_it))
         continue;
-      for (auto coord = v_it->point().cartesian_begin(); coord != v_it->point().cartesian_end(); ++coord)
-        ofs << *coord << " ";
-      ofs << "508\n";
       index_of_vertex[v_it] = ind++;
     }
-  for (auto coord = p.cartesian_begin(); coord != p.cartesian_end(); ++coord)
-    ofs << *coord << " ";
-  ofs << "208\n";
-  /*
-  int nbFacets = 0;
-  for (auto ft_it = t.finite_facets_begin(); ft_it != t.finite_facets_end(); ++ft_it)
-    nbFacets++;
-  ofs << "\nEdges\n" << nbFacets << "\n\n";
-  for (auto ft_it = t.facets_begin(); ft_it != t.facets_end(); ++ft_it)
-    {
-      if (t.is_infinite(ft_it))
-        continue;
-      for (auto vh_it = ft_it->vertices_begin(); vh_it != ft_it->vertices_end(); ++vh_it)
-        ofs << index_of_vertex[*vh_it] << " "; 
-    }
-  */
-  ofs << "Triangles " << t.number_of_finite_full_cells()+1 << "\n";
-  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
-    {
-      if (t.is_infinite(fc_it))
-        continue;
-      for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it)
-        ofs << index_of_vertex[*vh_it] << " ";
-      ofs << "508\n";
-    }
-  ofs << nbV << " " << nbV << " " << nbV << " " << 208 << "\n";
-  ofs << "End\n";
-  ofs.close();
-}
-
-bool triangulation_is_protected(Delaunay_triangulation& t, FT delta)
-{
-  // Verification part
-  Euclidean_distance ed;
   for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
     if (!t.is_infinite(fc_it))
-      for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
-        if (!t.is_infinite(v_it))
+      {
+        std::vector<Point_d> vertices;
+        for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+          vertices.push_back((*fc_v_it)->point());                
+        Sphere_d cs( vertices.begin(), vertices.end());
+        Point_d center_cs = cs.center();
+        FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(0)->point()));
+        for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
+          if (!t.is_infinite(v_it))
           //check if vertex belongs to the face
-          if (!vertex_is_in_full_cell(v_it, fc_it))
-            {
-              std::vector<Point_d> vertices;
-              for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
-                vertices.push_back((*fc_v_it)->point());
-              Sphere_d cs( vertices.begin(), vertices.end());
-              Point_d center_cs = cs.center();
-              FT r = sqrt(ed.transformed_distance(center_cs, fc_it->vertex(0)->point()));
-              FT dist2 = ed.transformed_distance(center_cs, v_it->point());
-              //if the new point is inside the protection ball of a non conflicting simplex
-              //std::cout << "Dist^2 = " << dist2 << " (r+delta)*(r+delta) = " << (r+delta)*(r+delta) << " r^2 = " << r*r <<"\n";
-              //if (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r)
-              if (dist2 <= r*r+delta*delta && dist2 >= r*r)
-                {
-                  write_delaunay_mesh(t, v_it->point());
-                  std::cout << "Problematic vertex " << *v_it << " ";
-                  std::cout << "Problematic cell " << *fc_it << "\n";
-                  std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", r^2+delta^2 = " << r*r+delta*delta << "\n";
-                  return false;
-                }
-            }
-        
+            if (!vertex_is_in_full_cell(v_it, fc_it))
+              {
+                FT dist2 = ed.transformed_distance(center_cs, v_it->point());
+                //if the new point is inside the protection ball of a non conflicting simplex
+                //std::cout << "Dist^2 = " << dist2 << " (r+delta)*(r+delta) = " << (r+delta)*(r+delta) << " r^2 = " << r*r <<"\n";
+                if (!power_protection)
+                  if (dist2 <= (r+delta)*(r+delta) && dist2 >= r*r)
+                    {
+                      write_delaunay_mesh(t, v_it->point(), is2d);
+                      // Output the problems
+                      std::cout << "Problematic vertex " << index_of_vertex[v_it] << " ";
+                      std::cout << "Problematic cell ";
+                      for (auto vh_it = fc_it->vertices_begin(); vh_it != fc_it->vertices_end(); ++vh_it)
+                        if (!t.is_infinite(*vh_it))
+                          std::cout << index_of_vertex[*vh_it] << " ";
+                      std::cout << "\n";
+                      std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", (r+delta)^2 = " << (r+delta)*(r+delta) << "\n";
+                      return false;
+                    }
+                if (power_protection)
+                  if (dist2 <= r*r+delta*delta && dist2 >= r*r)
+                    {
+                      write_delaunay_mesh(t, v_it->point(), is2d);
+                      std::cout << "Problematic vertex " << *v_it << " ";
+                      std::cout << "Problematic cell " << *fc_it << "\n";
+                      std::cout << "r^2 = " << r*r << ", d^2 = " << dist2 << ", r^2+delta^2 = " << r*r+delta*delta << "\n";
+                      return false;
+                    }
+              }
+      } 
   return true;
 }
 
-void fill_landmarks(Point_Vector& W, Point_Vector& landmarks, std::vector<int>& landmarks_ind)
-{
-  for (unsigned j = 0; j < landmarks_ind.size(); ++j)
-    landmarks.push_back(W[landmarks_ind[j]]);
-}
-
-void fill_full_cell_vector(Delaunay_triangulation& t, std::vector<std::vector<int>>& full_cells)
-{
-  // Store vertex indices in a map
-  int ind = 0; //index of a vertex
-  std::map<Delaunay_triangulation::Vertex_handle, int> index_of_vertex;
-  for (auto v_it = t.vertices_begin(); v_it != t.vertices_end(); ++v_it)
-    if (t.is_infinite(v_it))
-      continue;
-    else
-      index_of_vertex[v_it] = ind++;
-  // Write full cells as vectors in full_cells
-  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
-    {
-      if (t.is_infinite(fc_it))
-        continue;
-      std::vector<int> cell;
-      for (auto v_it = fc_it->vertices_begin(); v_it != fc_it->vertices_end(); ++v_it)
-        cell.push_back(index_of_vertex[*v_it]);
-      full_cells.push_back(cell);
-    }
-}
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
+// SAMPLING RADIUS 
+//////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-FT sampling_radius(Delaunay_triangulation& t)
+FT sampling_radius(Delaunay_triangulation& t, FT epsilon0)
 {
-  FT epsilon2 = 4.0;
+  FT epsilon2 = 0;
+  Point_d control_point;
   for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
     {
       if (t.is_infinite(fc_it))
@@ -578,134 +220,106 @@ FT sampling_radius(Delaunay_triangulation& t)
       if (!in_cube)
         continue;
       FT r2 = Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin()));
-      if (epsilon2 > r2)
-        epsilon2 = r2;
+      if (epsilon2 < r2)
+        {
+          epsilon2 = r2;
+          control_point = (*vertices.begin());
+        }
+    }
+  if (epsilon2 < epsilon0*epsilon0)
+    {
+      std::cout << "ACHTUNG! E' < E\n";
+      std::cout << "eps = " << epsilon0 << " eps' = " << sqrt(epsilon2) << "\n";
+      write_delaunay_mesh(t, control_point, is2d);
     }
   return sqrt(epsilon2);
 }
 
-FT point_sampling_radius_by_delaunay(Point_Vector& points)
+FT point_sampling_radius_by_delaunay(Point_Vector& points, FT epsilon0)
 {
   Delaunay_triangulation t(points[0].size());
   t.insert(points.begin(), points.end());
-  return sampling_radius(t);
+  return sampling_radius(t, epsilon0);
 }
 
-void landmark_choice_protected_delaunay(Point_Vector& W, int nbP, Point_Vector& landmarks, std::vector<int>& landmarks_ind, FT delta, std::vector<std::vector<int>>& full_cells)
+// A little script to make a tikz histogram of epsilon distribution
+// Returns the average epsilon
+FT epsilon_histogram(Delaunay_triangulation& t, int n)
 {
-  unsigned D = W[0].size();
-  Torus_distance td;
-  Euclidean_distance ed;
-  Delaunay_triangulation t(D);
-  CGAL::Random rand;
-  int landmark_count = 0;
-  std::list<int> index_list;
-  //      shuffle the list of indexes (via a vector)
-  {
-    std::vector<int> temp_vector;
-    for (int i = 0; i < nbP; ++i)
-      temp_vector.push_back(i);
-    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
-    std::shuffle(temp_vector.begin(), temp_vector.end(), std::default_random_engine(seed));
-    //CGAL::spatial_sort(temp_vector.begin(), temp_vector.end());
-    for (std::vector<int>::iterator it = temp_vector.begin(); it != temp_vector.end(); ++it)
-      index_list.push_front(*it);
-  }
-  for (unsigned pos1 = 0; pos1 < D+1; ++pos1)
-    {
-      std::vector<FT> point;
-      for (unsigned i = 0; i < pos1; ++i)
-        point.push_back(-1);
-      if (pos1 != D)
-        point.push_back(1);
-      for (unsigned i = pos1+1; i < D; ++i)
-        point.push_back(0);
-      assert(point.size() == D);
-      W[index_list.front()] = Point_d(point);
-      insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count);
-      index_list.pop_front();
-    }
-  //      add the first D+1 vertices to form one finite cell
-  /*
-  for (int i = 0; i <= D+1; ++i)
-    {
-      t.insert(W[index_list.front()]);
-      insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count);
-      index_list.pop_front();
-    }
-  */
-  /*
-  {   
-    std::vector<FT> coords;
-    for (int i = 0; i < D; ++i)
-      coords.push_back(-1);
-    W[index_list.front()] = Point_d(coords);
-    insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count);
-    index_list.pop_front();
-    for (int i = 0; i < D; ++i)
-      {
-        coords.clear();
-        for (int j = 0; j < D; ++j)
-          if (i == j)
-            coords.push_back(1);
-          else
-            coords.push_back(-1);
-        W[index_list.front()] = Point_d(coords);
-        insert_delaunay_landmark_with_copies(W, index_list.front(), landmarks_ind, t, landmark_count);
-        index_list.pop_front();
-      }
-  }
-  */
-  //std::cout << t;
-  //assert(t.number_of_vertices() == D+1);
-  //assert(landmarks_ind.size() == D+1);
-  //assert(W[landmarks_ind[0]][0] == 0);
-  //      add other vertices if they don't violate protection
-  std::list<int>::iterator list_it = index_list.begin();
-  while (list_it != index_list.end())
+  FT epsilon_max = 0; //sampling_radius(t,0);
+  FT sum_epsilon = 0;
+  int count_simplices = 0;
+  std::vector<int> histo(n+1, 0);
+  for (auto fc_it = t.full_cells_begin(); fc_it != t.full_cells_end(); ++fc_it)
     {
-      if (!is_violating_protection(W[*list_it], t, D, delta))
-        {
-          // If no conflicts then insert in every copy of T^3      
-          is_violating_protection(W[*list_it], t, D, delta);
-          insert_delaunay_landmark_with_copies(W, *list_it, landmarks_ind, t, landmark_count);
-          index_list.erase(list_it);
-          list_it = index_list.begin();
-          //std::cout << "index_list_size() = " << index_list.size() << "\n";
-        }
-      else
-        {
-          list_it++;
-          //std::cout << "!!!!!WARNING!!!!! A POINT HAS BEEN OMITTED!!!\n";
-        }
-      //write_delaunay_mesh(t, W[*list_it]);
+      if (t.is_infinite(fc_it))
+        continue;
+      Point_Vector vertices;
+      for (auto fc_v_it = fc_it->vertices_begin(); fc_v_it != fc_it->vertices_end(); ++fc_v_it)
+        vertices.push_back((*fc_v_it)->point());
+      Sphere_d cs( vertices.begin(), vertices.end());
+      Point_d csc = cs.center();
+      bool in_cube = true; 
+      for (auto xi = csc.cartesian_begin(); xi != csc.cartesian_end(); ++xi)
+        if (*xi > 1.0 || *xi < -1.0)
+          {
+            in_cube = false; break;
+          }
+      if (!in_cube)
+        continue;
+      FT r = sqrt(Euclidean_distance().transformed_distance(cs.center(), *(vertices.begin())));
+      if (r > epsilon_max)
+        epsilon_max = r;
+      sum_epsilon += r;
+      count_simplices++;
+      histo[floor(r/epsilon_max*n)]++;
     }
-  fill_landmarks(W, landmarks, landmarks_ind);
-  fill_full_cell_vector(t, full_cells);
-  if (triangulation_is_protected(t, delta))
-    std::cout << "Triangulation is ok\n";
-  else
-    std::cout << "Triangulation is BAD!! T_T しくしく!\n";
-  write_delaunay_mesh(t, Point_d(std::vector<FT>({0,0})));
-  //std::cout << t << std::endl;
+  std::ofstream ofs ("histogram.tikz", std::ofstream::out);
+  FT barwidth = 20.0/n;
+  int max_value = *(std::max_element(histo.begin(), histo.end()));
+  std::cout << max_value << std::endl;
+  FT ten_power = pow(10, ceil(log10(max_value)));
+  FT max_histo = ten_power;
+  if (max_value/ten_power < 2)
+    max_histo = 0.2*ten_power;
+  if (max_value/ten_power < 5)
+    max_histo = 0.5*ten_power;  
+  std::cout << ceil(log10(max_value)) << std::endl << max_histo << std::endl;
+  FT unitht = max_histo/10.0;
+
+  ofs << "\\draw[->] (0,0) -- (0,11);\n" <<
+    "\\draw[->] (0,0) -- (21,0);\n" <<
+    "\\foreach \\i in {1,...,10}\n" <<
+    "\\draw (0,\\i) -- (-0.1,\\i);\n" <<
+    "\\foreach \\i in {1,...,20}\n" <<
+    "\\draw (\\i,0) -- (\\i,-0.1);\n" <<
+ 
+    "\\node at (-1,11) {$\\epsilon$};\n" << 
+    "\\node at (22,-1) {$\\epsilon/\\epsilon_{max}$};\n" << 
+    "\\node at (-0.5,-0.5) {0};\n" << 
+    "\\node at (-0.5,10) {" << max_histo << "};\n" << 
+    "\\node at (20,-0.5) {1};\n";
+    
+
+  for (int i = 0; i < n; ++i)
+    ofs << "\\draw (" << barwidth*i << "," << histo[i]/unitht << ") -- ("
+        << barwidth*(i+1) << "," << histo[i]/unitht << ") -- ("
+  << barwidth*(i+1) << ",0) -- (" << barwidth*i << ",0) -- cycle;\n";
+  
+  ofs.close();
+
+  //return sum_epsilon/count_simplices;
+  return epsilon_max;
 }
 
-template <typename T>
-void print_vector(std::vector<T> v)
+FT epsilon_histogram_by_delaunay(Point_Vector& points, int n)
 {
-  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 << "]";
+  Delaunay_triangulation t(points[0].size());
+  t.insert(points.begin(), points.end());
+  return epsilon_histogram(t, n);
 }
 
+
 int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std::vector<int>& landmarks_ind, std::vector<std::vector<int>>& full_cells)
 {
   //******************** Preface: origin point
@@ -764,7 +378,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std
         }
     }
   std::string out_file = "wl_result";
-  write_wl(out_file,WL);
+  //write_wl(out_file,WL);
   
   //******************** Constructng a witness complex
   std::cout << "Entered witness complex construction\n";
@@ -787,7 +401,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std
     not_in << " are not.\n";
   
   //******************** Making a set of bad link landmarks
-  /*
+  
   std::cout << "Entered bad links\n";
   std::set< int > perturbL;
   int count_badlinks = 0;
@@ -814,7 +428,7 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std
     if (count_bad[i] != 0)
       std::cout << "count_bad[" << i << "] = " << count_bad[i] << std::endl;
   std::cout << "\nBad links total: " << count_badlinks << " Points to perturb: " << perturbL.size() << std::endl;
-  */
+  
   //*********************** Perturb bad link landmarks  
   /*
   for (auto u: perturbL)
@@ -848,16 +462,19 @@ int landmark_perturbation(Point_Vector &W, int nbL, Point_Vector& landmarks, std
       ofs.close();
     }
   
-  write_edges("landmarks/edges", witnessComplex, landmarks);
+  //write_edges("landmarks/edges", witnessComplex, landmarks);
   /*
   return count_badlinks;
   */
   return 0;
 }
 
-
 int main (int argc, char * const argv[])
 {
+  power_protection = true;//false;
+  grid_points = false;//true;
+  torus = true;
+  
   if (argc != 4)
     {
       std::cerr << "Usage: " << argv[0]
@@ -866,40 +483,98 @@ int main (int argc, char * const argv[])
     }
   int nbP       = atoi(argv[1]);
   int dim       = atoi(argv[2]);
-  double delta  = atof(argv[3]);
+  double theta0 = atof(argv[3]);
+  //double delta  = atof(argv[3]);
+
+  is2d = (dim == 2);
   
   std::cout << "Let the carnage begin!\n";
   Point_Vector point_vector;
-  generate_points_random_box(point_vector, nbP, dim);
-  FT epsilon = point_sampling_radius_by_delaunay(point_vector);
+  if (grid_points)
+    {
+      generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim, torus);
+      nbP = (int)pow((int)pow(nbP, 1.0/dim), dim);
+    }
+  else
+    generate_points_random_box(point_vector, nbP, dim);
+  FT epsilon = point_sampling_radius_by_delaunay(point_vector, 0);
+  //FT epsilon = epsilon_histogram_by_delaunay(point_vector,50);
   std::cout << "Initial epsilon = " << epsilon << std::endl;
   Point_Vector L;
   std::vector<int> chosen_landmarks;
   //write_points("landmarks/initial_pointset",point_vector);
   //write_points("landmarks/initial_landmarks",L);
   CGAL::Timer timer;
+
+  int n = 1;
+  std::vector<FT> values(n,0);
+  std::vector<FT> time(n,0);
+
+  //FT step = 0.001;
+  //FT delta = 0.01*epsilon;
+  //FT alpha = 0.5;
+  //FT step = atof(argv[3]);
+
+  start_experiments(point_vector, theta0, chosen_landmarks, epsilon);
+    
+  // for (int i = 0; i < n; i++)
+  // //for (int i = 0; bl > 0; i++)
+  //   {
+  //     //std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n";
+  //     //double delta = pow(10, -(1.0*i)/2);
+  //     //delta = step*i*epsilon;
+  //     //theta0 = step*i;
+  //     std::cout << "delta/epsilon = " << delta/epsilon << std::endl;
+  //     std::cout << "theta0 = " << theta0 << std::endl;
+  //     // Averaging the result
+  //     int sum_values = 0;
+  //     int nb_iterations = 1;
+  //     std::vector<std::vector<int>> full_cells;
+  //     for (int i = 0; i < nb_iterations; ++i)
+  //       {
+  //         //L = {};
+  //         chosen_landmarks = {};
+  //         //full_cells = {};
+  //         //timer.start();
+  //         //protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, epsilon, alpha, theta0, full_cells, torus, power_protection);
+  //         protected_delaunay(point_vector, chosen_landmarks, delta, epsilon, alpha, theta0, torus, power_protection);
+  //         //timer.stop();
+  //         sum_values += chosen_landmarks.size();
+  //       }
+  //     //FT epsilon2 = point_sampling_radius_by_delaunay(L, epsilon);
+  //     //std::cout << "Final epsilon = " << epsilon2  << ". Ratio = " << epsilon2/epsilon << std::endl;
+  //     //write_points("landmarks/initial_landmarks",L);
+  //     //std::cout << "delta/epsilon' = " << delta/epsilon2 << std::endl;
+  //     FT nbL = (sum_values*1.0)/nb_iterations;
+  //     //values[i] = pow((1.0*nbL)/nbP, -1.0/dim);
+  //     values[i] = (1.0*nbL)/nbP;
+  //     std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl;
+  //     //landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells);
+  //     time[i] = timer.time();
+  //     timer.reset();
+  //     //write_points("landmarks/landmarks0",L);
+  //   }
+  
+  // // OUTPUT A PLOT
+  // FT hstep = 20.0/(n-1);
+  // FT wstep = 10.0;
+
+  // std::ofstream ofs("N'Nplot.tikz", std::ofstream::out);
+  // ofs << "\\draw[red] (0," << wstep*values[0] << ")";
+  // for (int i = 1; i < n; ++i)
+  //   ofs << " -- (" << hstep*i << "," << wstep*values[i] << ")";
+  // ofs << ";\n";
+  // ofs.close();
   /*
-  for (int i = 0; i < 11; i++)
-  //for (int i = 0; bl > 0; i++)
-    {
-      //std::cout << "========== Start iteration " << i << "== curr_min(" << curr_min << ")========\n";
-      double delta = pow(10, -(1.0*i)/2);
-      std::cout << "delta = " << delta << std::endl;
-      L = {}; chosen_landmarks = {};
-      std::vector<std::vector<int>> full_cells;
-      timer.start();
-      landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, full_cells);
-      timer.stop();
-      FT epsilon2 = point_sampling_radius_by_delaunay(L);
-      std::cout << "Final epsilon = " << epsilon2  << ". Ratio = " << epsilon/epsilon2 << std::endl;
-      write_points("landmarks/initial_landmarks",L);
-      int nbL = chosen_landmarks.size();
-      std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl;
-      landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells);
-      timer.reset();
-      //write_points("landmarks/landmarks0",L);
-    }
-  */
+  wstep = 0.1;
+  ofs = std::ofstream("time.tikz", std::ofstream::out);
+  ofs << "\\draw[red] (0," << wstep*time[0] << ")";
+  for (int i = 1; i < n; ++i)
+    ofs << " -- (" << hstep*i << "," << wstep*time[i] << ")";
+  ofs << ";\n";
+  ofs.close();
+  
+  
   std::vector<std::vector<int>> full_cells;
   timer.start();
   landmark_choice_protected_delaunay(point_vector, nbP, L, chosen_landmarks, delta, full_cells);
@@ -909,6 +584,7 @@ int main (int argc, char * const argv[])
   write_points("landmarks/initial_landmarks",L);
   int nbL = chosen_landmarks.size();
   std::cout << "Number of landmarks = " << nbL << ", time= " << timer.time() << "s"<< std::endl;
-  landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells);
+  //landmark_perturbation(point_vector, nbL, L, chosen_landmarks, full_cells);
   timer.reset();
+  */
 }
diff --git a/src/Witness_complex/example/witness_complex_flat_torus.cpp b/src/Witness_complex/example/witness_complex_flat_torus.cpp
index 69ef5fbf..49383154 100644
--- a/src/Witness_complex/example/witness_complex_flat_torus.cpp
+++ b/src/Witness_complex/example/witness_complex_flat_torus.cpp
@@ -776,8 +776,8 @@ int main (int argc, char * const argv[])
   std::cout << "Let the carnage begin!\n";
   Point_Vector point_vector;
   //read_points_cust(file_name, point_vector);
-  generate_points_random_box(point_vector, nbP, dim);
-  //generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim);
+  //generate_points_random_box(point_vector, nbP, dim);
+  generate_points_grid(point_vector, (int)pow(nbP, 1.0/dim), dim);
   //nbP = (int)(pow((int)pow(nbP, 1.0/dim), dim));
   /*
   for (auto &p: point_vector)
diff --git a/src/Witness_complex/example/witness_complex_knn_landmarks.cpp b/src/Witness_complex/example/witness_complex_knn_landmarks.cpp
index e4a1c324..c45bc0c1 100644
--- a/src/Witness_complex/example/witness_complex_knn_landmarks.cpp
+++ b/src/Witness_complex/example/witness_complex_knn_landmarks.cpp
@@ -32,6 +32,8 @@
 //#include "gudhi/graph_simplicial_complex.h"
 #include "gudhi/Witness_complex.h"
 #include "gudhi/reader_utils.h"
+#include "generators.h"
+#include "output.h"
 //#include <boost/filesystem.hpp>
 
 //#include <CGAL/Delaunay_triangulation.h>
@@ -73,60 +75,6 @@ typedef K_neighbor_search::iterator KNS_range;
 typedef boost::container::flat_map<int, int> Point_etiquette_map;
 
 typedef std::vector<Point_d> Point_Vector;
-/**
- * \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 read_points_to_tree (std::string file_name, Tree& tree)
-{
-  //I assume here that tree is empty
-  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> coords;
-      std::istringstream iss( line );
-      while(iss >> x) { coords.push_back(x); }
-      if (coords.size() != 1)
-        {
-          Point_d point(coords.begin(), coords.end());
-          tree.insert(point);
-        }
-    }
-  in_file.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
@@ -184,19 +132,6 @@ void d_nearest_landmarks(Point_Vector &W, Tree &L, Point_etiquette_map &L_i, std
     }
 }
 
-
-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();
-}
-
 int main (int argc, char * const argv[])
 {
   if (argc != 3)
@@ -270,6 +205,6 @@ int main (int argc, char * const argv[])
 
   out_file = "output/"+file_name+"_"+argv[2]+".badlinks";
   std::ofstream ofs2(out_file, std::ofstream::out);
-  witnessComplex.write_bad_links(ofs2);
+  //witnessComplex.write_bad_links(ofs2);
   ofs2.close();
 }
diff --git a/src/Witness_complex/example/witness_complex_protected_delaunay.cpp b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp
index 2f795a5f..77a167a5 100644
--- a/src/Witness_complex/example/witness_complex_protected_delaunay.cpp
+++ b/src/Witness_complex/example/witness_complex_protected_delaunay.cpp
@@ -268,6 +268,15 @@ void insert_delaunay_landmark_with_copies(Point_Vector& W, int chosen_landmark,
   landmark_count++;
 }
 
+
+
+
+////////////////////////////////////////////////////////////////////////
+// OLD CODE VVVVVVVV
+////////////////////////////////////////////////////////////////////////
+
+
+/*
 bool is_violating_protection(Point_d& p, Delaunay_triangulation& t, int D, FT delta)
 {
   Euclidean_distance ed;
@@ -592,3 +601,4 @@ int main (int argc, char * const argv[])
     }
   
 }
+*/
diff --git a/src/Witness_complex/example/witness_complex_sphere.cpp b/src/Witness_complex/example/witness_complex_sphere.cpp
index 550c9392..bf3015fa 100644
--- a/src/Witness_complex/example/witness_complex_sphere.cpp
+++ b/src/Witness_complex/example/witness_complex_sphere.cpp
@@ -35,6 +35,8 @@
 //#include "gudhi/graph_simplicial_complex.h"
 #include "gudhi/Witness_complex.h"
 #include "gudhi/reader_utils.h"
+#include "generators.h"
+#include "output.h"
 //#include <boost/filesystem.hpp> 
 
 //#include <CGAL/Delaunay_triangulation.h>
@@ -94,101 +96,9 @@ 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_grid(Point_Vector& W, int width, int D)
-{
- 
-}
-
-void generate_points_random_box(Point_Vector& W, int nbP, int dim)
-{
-  Random_cube_iterator rp(dim, 1);
-  for (int i = 0; i < nbP; i++)
-    {
-      W.push_back(*rp++);
-    }
-}
-
-/* NOT TORUS RELATED
- */
-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 read_points_to_tree (std::string file_name, Tree& tree)
-{
-  //I assume here that tree is empty
-  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> coords;
-      std::istringstream iss( line );
-      while(iss >> x) { coords.push_back(x); }
-      if (coords.size() != 1)
-        {
-          Point_d point(coords.begin(), coords.end());
-          tree.insert(point);
-        }
-    }
-  in_file.close();
-}
-*/
-
-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();
-}
-
-
 std::vector<Point_d> convert_to_torus(std::vector< Point_d>& points)
 {
   std::vector< Point_d > points_torus;
@@ -205,82 +115,6 @@ std::vector<Point_d> convert_to_torus(std::vector< Point_d>& points)
   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
  */
@@ -356,6 +190,7 @@ void landmark_choice_600cell(Point_Vector&W, int nbP, int nbL, Point_Vector& lan
 int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector<int>& landmarks_ind)
 {
   //********************Preface: origin point
+  clock_t start, end;
   int D = W[0].size();
   std::vector<FT> orig_vector;
   for (int i=0; i<D; i++)
@@ -383,6 +218,7 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector<
   */
   std::cout << "Enter (D+1) nearest landmarks\n";
   //std::cout << "Size of the tree is " << L.size() << std::endl;
+  start = clock();
   for (int i = 0; i < nbP; i++)
     {
       //std::cout << "Entered witness number " << i << std::endl;
@@ -416,7 +252,9 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector<
         }
     }
   //std::cout << "\n";
-  
+  end = clock();
+  std::cout << "Landmark choice for " << nbL << " landmarks took "
+            << (double)(end-start)/CLOCKS_PER_SEC << " s. \n";
   std::string out_file = "wl_result";
   write_wl(out_file,WL);
   
@@ -424,14 +262,19 @@ int landmark_perturbation(Point_Vector &W, Point_Vector& landmarks, std::vector<
   std::cout << "Entered witness complex construction\n";
   Witness_complex<> witnessComplex;
   witnessComplex.setNbL(nbL);
+  start = clock();
   witnessComplex.witness_complex(WL);
+  //
+  end = clock();
+  std::cout << "Howdy world! The process took "
+       << (double)(end-start)/CLOCKS_PER_SEC << " s. \n";
+  //witnessComplex.witness_complex(WL);
   /*
   if (witnessComplex.is_witness_complex(WL))
     std::cout << "!!YES. IT IS A WITNESS COMPLEX!!\n";
   else
     std::cout << "??NO. IT IS NOT A WITNESS COMPLEX??\n";   
   */ 
-  */
   //******************** Making a set of bad link landmarks
   std::cout << "Entered bad links\n";
   std::set< int > perturbL;
@@ -575,8 +418,8 @@ int main (int argc, char * const argv[])
       */
     }
   int bl = nbL, curr_min = bl;
-  write_points("landmarks/initial_pointset",point_vector);
-  write_points("landmarks/initial_landmarks",L);
+  //write_points("landmarks/initial_pointset",point_vector);
+  //write_points("landmarks/initial_landmarks",L);
   
   for (int i = 0; bl > 0; i++)
     //for (int i = 0; i < 1; i++)
@@ -585,7 +428,7 @@ int main (int argc, char * const argv[])
       bl=landmark_perturbation(point_vector, L, chosen_landmarks);
       if (bl < curr_min)
         curr_min=bl;
-      write_points("landmarks/landmarks0",L);
+      //write_points("landmarks/landmarks0",L);
     }
   //end = clock();