Peer review fix

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

Former-commit-id: e2896c0a1dc03cb8f87771910b6c756a755f258c

2 files changed, 367 insertions, 414 deletions

diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index b9b8f0ea..b38e5813 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -26,13 +26,14 @@
 #include <gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h>
 #include <gudhi/Simplex_tree/Simplex_tree_siblings.h>
 #include <gudhi/Simplex_tree/Simplex_tree_iterators.h>
-#include <gudhi/reader_utils.h>
 #include <gudhi/Simplex_tree/indexing_tag.h>
 
+#include <gudhi/reader_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+
 #include <boost/container/flat_map.hpp>
 #include <boost/iterator/transform_iterator.hpp>
 #include <boost/graph/adjacency_list.hpp>
-#include <gudhi/graph_simplicial_complex.h>
 
 #include <algorithm>
 #include <utility>
@@ -234,8 +235,8 @@ class Simplex_tree {
   Filtration_simplex_range filtration_simplex_range(Indexing_tag) {
     if (filtration_vect_.empty()) {
       initialize_filtration();
-    } return Filtration_simplex_range(filtration_vect_.begin(),
-                                    filtration_vect_.end());
+    }
+    return Filtration_simplex_range(filtration_vect_.begin(), filtration_vect_.end());
   }
 
   Filtration_simplex_range filtration_simplex_range() {
@@ -281,45 +282,48 @@ class Simplex_tree {
   Simplex_tree()
       : null_vertex_(-1),
       threshold_(0),
-      root_(NULL, null_vertex_),
+      root_(nullptr, null_vertex_),
       filtration_vect_(),
       dimension_(-1) { }
 
-  /** \brief Copy; copy the whole tree structure. */
-  Simplex_tree(Simplex_tree& copy) : 	null_vertex_(copy.null_vertex_),
-  										threshold_(copy.threshold_),
-										root_(NULL, -1, copy.root_.members()),
-										filtration_vect_(copy.filtration_vect_),
-										dimension_(copy.dimension_) {
-		rec_copy(&root_, &copy.root_);
-	}
-
-   /** rec_copy: depth first search, inserts simplices when reaching a leaf.*/
-  void rec_copy(Siblings *sib, Siblings *sib_copy)
-  {
-	for (auto sh = sib->members().begin(), sh_copy = sib_copy->members().begin(); sh != sib->members().end(); ++sh, ++sh_copy) {
-		if (has_children(sh_copy)) {
-//			boost::container::flat_map<Vertex_handle, Node> copy(sh_copy->second.children()->members());
-			Siblings * newsib =  new Siblings (sib, sh_copy->first);
-			for (auto it = sh_copy->second.children()->members().begin(); it != sh_copy->second.children()->members().end(); ++it)
-				newsib->members_.emplace(it->first, Node(sib, it->second.filtration()));
-			rec_copy(newsib, sh_copy->second.children());
-			sh->second.assign_children(newsib);
-		}
-	}
+  /** \brief User-defined copy constructor reproduces the whole tree structure. */
+  Simplex_tree(const Simplex_tree& simplex_source)
+      : null_vertex_(simplex_source.null_vertex_),
+      threshold_(simplex_source.threshold_),
+      filtration_vect_(simplex_source.filtration_vect_),
+      dimension_(simplex_source.dimension_) {
+    auto root_source = simplex_source.root_;
+    auto memb_source = root_source.members();
+    root_ = Siblings(nullptr, null_vertex_, memb_source);
+    rec_copy(&root_, &root_source);
   }
 
+  /** \brief depth first search, inserts simplices when reaching a leaf. */
+  void rec_copy(Siblings *sib, Siblings *sib_copy) {
+    for (auto sh = sib->members().begin(), sh_copy = sib_copy->members().begin();
+         sh != sib->members().end(); ++sh, ++sh_copy) {
+      if (has_children(sh_copy)) {
+        Siblings * newsib = new Siblings(sib, sh_copy->first);
+        newsib->members_.reserve(sh_copy->second.children()->members().size());
+        for (auto it = sh_copy->second.children()->members().begin();
+             it != sh_copy->second.children()->members().end(); ++it)
+          newsib->members_.emplace_hint(newsib->members_.end(), it->first, Node(sib, it->second.filtration()));
+        rec_copy(newsib, sh_copy->second.children());
+        sh->second.assign_children(newsib);
+      }
+    }
+  }
 
-
-
-  /** \brief Move; moves the whole tree structure. */
-  Simplex_tree(Simplex_tree&& old) : 	null_vertex_(std::move(old.null_vertex_)),
-  										threshold_(std::move(old.threshold_)),
-										root_(std::move(old.root_)), filtration_vect_(std::move(old.filtration_vect_)),
-										dimension_(std::move(old.dimension_)) {
-		old.dimension_ = -1;
-		old.threshold_ = 0;
-	}
+  /** \brief User-defined move constructor moves the whole tree structure. */
+  Simplex_tree(Simplex_tree && old)
+      : null_vertex_(std::move(old.null_vertex_)),
+      threshold_(std::move(old.threshold_)),
+      root_(std::move(old.root_)), filtration_vect_(std::move(old.filtration_vect_)),
+      dimension_(std::move(old.dimension_)) {
+    old.dimension_ = -1;
+    old.threshold_ = 0;
+    old.root_ = Siblings(nullptr, null_vertex_);
+  }
 
   /** \brief Destructor; deallocates the whole tree structure. */
   ~Simplex_tree() {
@@ -341,71 +345,36 @@ class Simplex_tree {
     delete sib;
   }
 
- public: 
-  /** \brief Prints the simplex_tree hierarchically. 
-   * Since it prints the vertices recursively, one can watch its tree shape.
-   */
-  void print_tree()
-  {
-	  for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh)
-	  {
-		  std::cout << sh->first << " ";
-		  if (has_children(sh))
-		  {
-			  std::cout << "(";
-			  rec_print(sh->second.children());
-			  std::cout << ")";
-		  }
-		  std::cout << std::endl;
-	  }
-  }
-
-
-  /** \brief Recursively prints the simplex_tree, using depth first search. */
- private:
-  void rec_print(Siblings * sib)
-  {
-	  for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh)
-	  {
-		  std::cout << " " << sh->first << " ";
-		  if (has_children(sh))
-		  {
-			  std::cout << "(";
-			  rec_print(sh->second.children());
-			  std::cout << ")";
-		  }
-	  }
-  }
-
  public:
+  /** \brief Checks if two simplex trees are equal. */
+  bool operator==(Simplex_tree st2) {
+    return (this->null_vertex_ == st2.null_vertex_
+            && this->threshold_ == st2.threshold_
+            && this->dimension_ == st2.dimension_
+            && rec_equal(&root_, &st2.root_));
+  }
 
-  /** \brief Checks whether or not two simplex trees are equal. */
-  bool operator ==(Simplex_tree st2)
-  {
-	  return (this->null_vertex_ == st2.null_vertex_ 
-	  		&& this->threshold_ == st2.threshold_ 
-			&& this->dimension_ == st2.dimension_ 
-			&& rec_equal(&root_, &st2.root_));
+  /** \brief Checks if two simplex trees are different. */
+  bool operator!=(Simplex_tree st2) {
+    return (!(*this == st2));
   }
 
   /** rec_equal: Checks recursively whether or not two simplex trees are equal, using depth first search. */
  private:
-  bool rec_equal(Siblings * s1, Siblings * s2)
-  {
-	  if (s1->members().size() != s2->members().size())
-		  return false;
-	  for (auto sh1 = s1->members().begin(), sh2 = s2->members().begin(); sh1 != s1->members().end(); ++sh1, ++sh2)
-	  {
-		  if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
-			  return false;
-		  if (has_children(sh1))
-			  return rec_equal(sh1->second.children(), sh2->second.children());
-		  else if (has_children(sh2))
-		  	return false;
-		  else
-			  return true;
-	  }
-	  return true;
+  bool rec_equal(Siblings * s1, Siblings * s2) {
+    if (s1->members().size() != s2->members().size())
+      return false;
+    for (auto sh1 = s1->members().begin(), sh2 = s2->members().begin(); sh1 != s1->members().end(); ++sh1, ++sh2) {
+      if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
+        return false;
+      if (has_children(sh1))
+        return rec_equal(sh1->second.children(), sh2->second.children());
+      else if (has_children(sh2))
+        return false;
+      else
+        return true;
+    }
+    return true;
   }
 
  public:
@@ -444,7 +413,7 @@ class Simplex_tree {
    *
    * One can call filtration(null_simplex()). */
   static Simplex_handle null_simplex() {
-    return Dictionary_it(NULL);
+    return Dictionary_it(nullptr);
   }
 
   /** \brief Returns a key different for all keys associated to the
@@ -491,7 +460,7 @@ class Simplex_tree {
   int dimension(Simplex_handle sh) {
     Siblings * curr_sib = self_siblings(sh);
     int dim = 0;
-    while (curr_sib != NULL) {
+    while (curr_sib != nullptr) {
       ++dim;
       curr_sib = curr_sib->oncles();
     }
@@ -503,32 +472,36 @@ class Simplex_tree {
     return dimension_;
   }
 
-
   /** \brief Returns true iff the node in the simplex tree pointed by
    * sh has children.*/
   bool has_children(Simplex_handle sh) const {
     return (sh->second.children()->parent() == sh->first);
   }
 
-  /** \brief Given a range of Vertex_handles, returns the Simplex_handle
+    /** \brief Given a range of Vertex_handles, returns the Simplex_handle
    * of the simplex in the simplicial complex containing the corresponding
    * vertices. Return null_simplex() if the simplex is not in the complex.
    *
    * The type InputVertexRange must be a range of <CODE>Vertex_handle</CODE>
    * on which we can call std::begin() function
    */
-template<class InputVertexRange>
+  template<class InputVertexRange>
   Simplex_handle find(const InputVertexRange & s) {
-  	std::vector<Vertex_handle> copy(std::begin(s), std::end(s));
-  	std::sort(std::begin(copy), std::end(copy));
-	return find_simplex(copy);
+    auto first = std::begin(s);
+    auto last = std::end(s);
+
+    if (first == last)
+      return null_simplex();  // ----->>
+
+    // Copy before sorting
+    std::vector<Vertex_handle> copy(first, last);
+    std::sort(std::begin(copy), std::end(copy));
+    return find_simplex(copy);
   }
 
  private:
   /** Find function, with a sorted range of vertices. */
   Simplex_handle find_simplex(std::vector<Vertex_handle> & simplex) {
-    if (simplex.begin() == simplex.end())
-	  return null_simplex();
     Siblings * tmp_sib = &root_;
     Dictionary_it tmp_dit;
     Vertex_handle last = simplex[simplex.size() - 1];
@@ -553,22 +526,19 @@ template<class InputVertexRange>
   //{ return root_.members_.find(v); }
 
  private:
-  /** Recursively insert a simplex */
-  template <typename RandomAccessVertexRange>
-  std::pair<Simplex_handle, bool> insert_simplex_rec(RandomAccessVertexRange & simplex,
-                                                 Filtration_value filtration) {
-    if (simplex.empty()) {
-      return std::pair<Simplex_handle, bool>(null_simplex(), true);
-    }
+  /** \brief Recursively insert a simplex represented by a vector of vertex.
+   \warning the vector must be sorted by increasing vertex handle order */
+  std::pair<Simplex_handle, bool> insert_vertex_vector(const std::vector<Vertex_handle>& simplex,
+                                                     Filtration_value filtration) {
     Siblings * curr_sib = &root_;
     std::pair<Simplex_handle, bool> res_insert;
-    typename RandomAccessVertexRange::iterator vi;
-    for (vi = simplex.begin(); vi != simplex.end() - 1; ++vi) {
+    auto vi = simplex.begin();
+    for (; vi != simplex.end() - 1; ++vi) {
       res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
       if (!(has_children(res_insert.first))) {
         res_insert.first->second.assign_children(new Siblings(curr_sib, *vi));
       }
-        curr_sib = res_insert.first->second.children();
+      curr_sib = res_insert.first->second.children();
     }
     res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
     if (!res_insert.second) {
@@ -584,7 +554,8 @@ template<class InputVertexRange>
     // otherwise the insertion has succeeded
     return res_insert;
   }
-public:
+
+ public:
   /** \brief Insert a simplex, represented by a range of Vertex_handles, in the simplicial complex.
    *
    * @param[in]  simplex    range of Vertex_handles, representing the vertices of the new simplex
@@ -594,7 +565,7 @@ public:
    * to the new simplex.
    * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion
    * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration',
-   * we assign this simplex with the new value 'filtration', and set the Simplex_handle filed of the
+   * we assign this simplex with the new value 'filtration', and set the Simplex_handle field of the
    * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to
    * null_simplex.
    *
@@ -608,58 +579,105 @@ public:
    *
    * The type RandomAccessVertexRange must be a range for which .begin() and
    * .end() return random access iterators, with 'value_type' Vertex_handle. */
-template<class RandomAccessVertexRange>
-std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & simplex,
-		Filtration_value filtration) {
-	std::vector<Vertex_handle> copy(std::begin(simplex), std::end(simplex));;
-  	std::sort(std::begin(copy), std::end(copy));
-	return insert_simplex_rec(copy, filtration);
-}
+  template<class RandomAccessVertexRange>
+  std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & simplex,
+                                                 Filtration_value filtration) {
+    auto first = std::begin(simplex);
+    auto last = std::end(simplex);
+
+    if (first == last)
+      return std::pair<Simplex_handle, bool>(null_simplex(), true);  // ----->>
+
+    // Copy before sorting
+    std::vector<Vertex_handle> copy(first, last);
+    std::sort(std::begin(copy), std::end(copy));
+    return insert_vertex_vector(copy, filtration);
+  }
 
-  /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
+    /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
    * Vertex_handles, in the simplicial complex.
    *
    * @param[in]  Nsimplex   range of Vertex_handles, representing the vertices of the new N-simplex
    * @param[in]  filtration the filtration value assigned to the new N-simplex.
+   * The return type is a pair. If the new simplex is inserted successfully (i.e. it was not in the
+   * simplicial complex yet) the bool is set to true and the Simplex_handle is the handle assigned
+   * to the new simplex.
+   * If the insertion fails (the simplex is already there), the bool is set to false. If the insertion
+   * fails and the simplex already in the complex has a filtration value strictly bigger than 'filtration',
+   * we assign this simplex with the new value 'filtration', and set the Simplex_handle field of the
+   * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to
+   * null_simplex.
    */
   template<class InputVertexRange>
-  void insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
+  std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
                                    Filtration_value filtration = 0.0) {
-	std::vector<Vertex_handle> copy(std::begin(Nsimplex), std::end(Nsimplex));;
-  	std::sort(std::begin(copy), std::end(copy));
-	insert_simplex_and_subfaces_rec(copy, filtration);
-	}
+    auto first = std::begin(Nsimplex);
+    auto last = std::end(Nsimplex);
 
- private:
+    if (first == last)
+      return std::pair<Simplex_handle, bool>(null_simplex(), true);  // ----->>
 
-  /** Recursively insert simplex and all of its subfaces */
-  template <typename RandomAccessVertexRange>
-  void insert_simplex_and_subfaces_rec(RandomAccessVertexRange & Nsimplex,
-                                         Filtration_value filtration = 0.0) {
-
-    if (Nsimplex.size() > 1) {
-      for (unsigned int NIndex = 0; NIndex < Nsimplex.size(); NIndex++) {
-        // insert N (N-1)-Simplex
-        RandomAccessVertexRange NsimplexMinusOne;
-        for (unsigned int NListIter = 0; NListIter < Nsimplex.size() - 1; NListIter++) {
-          // (N-1)-Simplex creation
-          NsimplexMinusOne.push_back(Nsimplex[(NIndex + NListIter) % Nsimplex.size()]);
-        }
-        // (N-1)-Simplex recursive call
-        insert_simplex_and_subfaces(NsimplexMinusOne, filtration);
+    // Copy before sorting
+    std::vector<Vertex_handle> copy(first, last);
+    std::sort(std::begin(copy), std::end(copy));
+
+    std::vector<std::vector<Vertex_handle>> to_be_inserted;
+    std::vector<std::vector<Vertex_handle>> to_be_propagated;
+    return rec_insert_simplex_and_subfaces(copy, to_be_inserted, to_be_propagated, filtration);
+  }
+
+ private:
+  std::pair<Simplex_handle, bool> rec_insert_simplex_and_subfaces(std::vector<Vertex_handle>& the_simplex,
+                                                             std::vector<std::vector<Vertex_handle>>& to_be_inserted,
+                                                             std::vector<std::vector<Vertex_handle>>& to_be_propagated,
+                                                             Filtration_value filtration = 0.0) {
+    std::pair<Simplex_handle, bool> insert_result;
+    if (the_simplex.size() > 1) {
+      // Get and remove last vertex
+      Vertex_handle last_vertex = the_simplex.back();
+      the_simplex.pop_back();
+      // Recursive call after last vertex removal
+      insert_result = rec_insert_simplex_and_subfaces(the_simplex, to_be_inserted, to_be_propagated, filtration);
+
+      // Concatenation of to_be_inserted and to_be_propagated
+      to_be_inserted.insert(to_be_inserted.begin(), to_be_propagated.begin(), to_be_propagated.end());
+      to_be_propagated = to_be_inserted;
+
+      // to_be_inserted treatment
+      for (auto& simplex_tbi : to_be_inserted) {
+        simplex_tbi.push_back(last_vertex);
+      }
+      std::vector<Vertex_handle> last_simplex(1, last_vertex);
+      to_be_inserted.insert(to_be_inserted.begin(), last_simplex);
+      // i.e. (0,1,2) =>
+      // [to_be_inserted | to_be_propagated] = [(1) (0,1) | (0)]
+      // [to_be_inserted | to_be_propagated] = [(2) (0,2) (1,2) (0,1,2) | (0) (1) (0,1)]
+      // N.B. : it is important the last inserted to be the highest in dimension
+      // in order to return the "last" insert_simplex result
+
+      // insert all to_be_inserted
+      for (auto& simplex_tbi : to_be_inserted) {
+        insert_result = insert_vertex_vector(simplex_tbi, filtration);
       }
-      // N-Simplex insert
-      std::pair<Simplex_handle, bool> returned = insert_simplex_rec(Nsimplex, filtration);
-    } else if (Nsimplex.size() == 1) {
-      // 1-Simplex insert - End of recursivity
-      std::pair<Simplex_handle, bool> returned = insert_simplex_rec(Nsimplex, filtration);
+    } else if (the_simplex.size() == 1) {
+      // When reaching the end of recursivity, vector of simplices shall be empty and filled on back recursive
+      if ((to_be_inserted.size() != 0) || (to_be_propagated.size() != 0)) {
+        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Error vector not empty" << std::endl;
+        exit(-1);
+      }
+      std::vector<Vertex_handle> first_simplex(1, the_simplex.back());
+      // i.e. (0,1,2) => [to_be_inserted | to_be_propagated] = [(0) | ]
+      to_be_inserted.push_back(first_simplex);
+
+      insert_result = insert_vertex_vector(first_simplex, filtration);
     } else {
-      // Nothing to insert - empty vector
+        std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error" << std::endl;
+        exit(-1);
     }
+    return insert_result;
   }
 
  public:
-
   /** \brief Assign a value 'key' to the key of the simplex
    * represented by the Simplex_handle 'sh'. */
   void assign_key(Simplex_handle sh, Simplex_key key) {
@@ -683,17 +701,6 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
       return sh->second.children();
   }
 
-  // void display_simplex(Simplex_handle sh)
-  // {
-  //   std::cout << "   " << "[" << filtration(sh) << "] ";
-  //   for( auto vertex : simplex_vertex_range(sh) )
-  //     { std::cout << vertex << " "; }
-  // }
-
-  // void print(Simplex_handle sh, std::ostream& os = std::cout)
-  // {  for(auto v : simplex_vertex_range(sh)) {os << v << " ";}
-  //    os << std::endl; }
-
  public:
   /** Returns a pointer to the root nodes of the simplex tree. */
   Siblings * root() {
@@ -737,81 +744,6 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
                      is_before_in_filtration(this));
   }
 
-/** \brief Contracts two vertices : the contracted vertex is erased from the three, and the remaining vertex receives all the links of the contracted one.
-	\param remaining The value of the vertex within the other is contracted
-	\param deleted The value of the vertex to be contracted
-*/
-	void edge_contraction(Vertex_handle remaining, Vertex_handle deleted)
-	{
-		std::vector<Vertex_handle> accessRemaining, accessDeleted;
-		accessRemaining.push_back(remaining);
-		accessDeleted.push_back(deleted);
-		Simplex_handle shr = find(accessRemaining), shd = find(accessDeleted);
-		if (has_children(shd))
-		{
-			Siblings * sibDeleted, * sibRemaining; // We need the siblings
-			sibDeleted = shd->second.children();
-			sibRemaining = shr->second.children();
-			rec_insert(sibDeleted, sibRemaining, shd->first);
-		}
-		rec_delete(&root_, shd->first);
-	}
-	
-
-/** \brief recursively insert the members of a Sibling into another, any time the replacedVertex appears into the target Sibling
-	\param sibInserted The sibling to be inserted
-	\param sibTarget The target sibling on which the other sibling is copied
-	\param replacedVertex The replacedVertex (Vertex_handle) needed to insert the sibling
-*/
-	void rec_insert(Siblings * sibInserted, Siblings * sibTarget, Vertex_handle replacedVertex)
-	{
-		for (auto sh = sibTarget->members().begin(); sh != sibTarget->members().end(); ++sh)
-		{
-			if (has_children(sh))
-				rec_insert(sibInserted, sh->second.children(), replacedVertex);
-			if (sh->first == replacedVertex)
-				insert_members(sibTarget, sibInserted);
-		}
-	}
-
-/** \brief Copy a Sibling into another, which is possibly not empty
-	\param sibInserted The sibling to be copied
-	\param sibTarget The sibling on which we wan't to copy the other sibling
-*/
-	void insert_members(Siblings * sibTarget, Siblings * sibInserted)
-	{
-		for (auto sh = sibInserted->members().begin(); sh != sibInserted->members().end(); ++sh)
-		{
-			std::pair<Vertex_handle, Node> member(sh->first, Node(sibTarget, sh->second.filtration()));
-			if (has_children(sh))
-			{
-				std::vector<std::pair<Vertex_handle, Node>> v(sh->second.children()->members().begin(), sh->second.children()->members().end());
-				Siblings * newsib =  new Siblings (sibTarget, sh->first);
-				for (auto it = v.begin(); it != v.end(); ++it)
-					newsib->members_.emplace(it->first, Node(sibTarget, it->second.filtration()));
-				insert_members(newsib, sh->second.children());
-				member.second.assign_children(newsib);
-
-			}
-			sibTarget->members_.emplace(member);
-		}
-	}
-
-/** \brief Erase every occurencies of a vertex in a Sibling
-	\param sib The sibling on which we wan't to erase the elements
-	\param deletedVertex The value of the members we wan't to erase
-*/
-	void rec_delete(Siblings * sib, Vertex_handle deletedVertex)
-	{
-		for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh)
-		{
-			if (has_children(sh))
-				rec_delete(sh->second.children(), deletedVertex);
-			if (sh->first == deletedVertex)
-				sib->members_.erase(sh);
-		}
-	}
-
  private:
   /** Recursive search of cofaces
    * This function uses DFS
@@ -888,13 +820,13 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
     assert(codimension >= 0);
     Simplex_vertex_range rg = simplex_vertex_range(simplex);
     std::vector<Vertex_handle> copy(rg.begin(), rg.end());
-    if (codimension + static_cast<int>(copy.size()) > dimension_ + 1 ||
-        (codimension == 0 && static_cast<int>(copy.size()) > dimension_))  // n+codimension greater than dimension_
+    if (codimension + static_cast<int> (copy.size()) > dimension_ + 1 ||
+        (codimension == 0 && static_cast<int> (copy.size()) > dimension_))  // n+codimension greater than dimension_
       return cofaces;
     // must be sorted in decreasing order
     assert(std::is_sorted(copy.begin(), copy.end(), std::greater<Vertex_handle>()));
     bool star = codimension == 0;
-    rec_coface(copy, &root_, 1, cofaces, star, codimension + static_cast<int>(copy.size()));
+    rec_coface(copy, &root_, 1, cofaces, star, codimension + static_cast<int> (copy.size()));
     return cofaces;
   }
 
@@ -1080,7 +1012,7 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
     while (true) {
       if (begin1->first == begin2->first) {
         Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
-        intersection.emplace_back(begin1->first, Node(NULL, filt));
+        intersection.emplace_back(begin1->first, Node(nullptr, filt));
         if (++begin1 == end1 || ++begin2 == end2)
           return;  // ----->>
       } else if (begin1->first < begin2->first) {
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index f09a11ea..d9666c52 100644
--- a/src/Simplex_tree/test/simplex_tree_unit_test.cpp
+++ b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
@@ -9,8 +9,8 @@
 #define BOOST_TEST_MODULE "simplex_tree"
 #include <boost/test/included/unit_test.hpp>
 
-#include "gudhi/graph_simplicial_complex.h"
-#include "gudhi/reader_utils.h"
+//  ^
+// /!\ Nothing else from Simplex_tree shall be included to test includes are well defined.
 #include "gudhi/Simplex_tree.h"
 
 using namespace Gudhi;
@@ -59,7 +59,6 @@ void test_iterators_on_empty_simplex_tree(typeST& tst) {
 BOOST_AUTO_TEST_CASE(simplex_tree_when_empty) {
   const Filtration_value DEFAULT_FILTRATION_VALUE = 0;
 
-  // TEST OF DEFAULT CONSTRUCTOR
   std::cout << "********************************************************************" << std::endl;
   std::cout << "TEST OF DEFAULT CONSTRUCTOR" << std::endl;
   typeST st;
@@ -125,10 +124,9 @@ void test_simplex_tree_contains(typeST& simplexTree, typeSimplex& simplex, int p
   std::cout << "test_simplex_tree_contains - filtration=" << simplexTree.filtration(*f_simplex) << "||" << simplex.second << std::endl;
   BOOST_CHECK(AreAlmostTheSame(simplexTree.filtration(*f_simplex), simplex.second));
 
-  int simplexIndex=simplex.first.size()-1;
+  int simplexIndex = simplex.first.size() - 1;
   std::sort(simplex.first.begin(), simplex.first.end()); // if the simplex wasn't sorted, the next test could fail
-  for( auto vertex : simplexTree.simplex_vertex_range(*f_simplex) )
-  {
+  for (auto vertex : simplexTree.simplex_vertex_range(*f_simplex)) {
     std::cout << "test_simplex_tree_contains - vertex=" << vertex << "||" << simplex.first.at(simplexIndex) << std::endl;
     BOOST_CHECK(vertex == simplex.first.at(simplexIndex));
     BOOST_CHECK(simplexIndex >= 0);
@@ -159,33 +157,17 @@ void set_and_test_simplex_tree_dim_fil(typeST& simplexTree, int vectorSize, cons
     std::cout << "   set_and_test_simplex_tree_dim_fil - max_fil=" << max_fil
         << std::endl;
   }
-  
+
   BOOST_CHECK(simplexTree.dimension() == dim_max);
   BOOST_CHECK(AreAlmostTheSame(simplexTree.filtration(), max_fil));
 
   // Another way to count simplices:
-  long long int num_simp = 0;
+  size_t num_simp = 0;
   for (auto f_simplex : simplexTree.complex_simplex_range()) {
     num_simp++;
   }
-  
-  BOOST_CHECK(simplexTree.num_simplices() == num_simp);
-}
 
-void test_cofaces(typeST& st, std::vector<Vertex_handle> v, int dim, std::vector<typeST::Simplex_handle> res) {
-  typeST::Cofaces_simplex_range cofaces;
-  if (dim == 0)
-    cofaces = st.star_simplex_range(st.find(v));
-  else
-    cofaces = st.cofaces_simplex_range(st.find(v), dim);
-  for (auto simplex = cofaces.begin(); simplex != cofaces.end(); ++simplex) {
-    typeST::Simplex_vertex_range rg = st.simplex_vertex_range(*simplex);
-    for (auto vertex = rg.begin(); vertex != rg.end(); ++vertex) {
-      std::cout << "(" << *vertex << ")";
-    }
-    std::cout << std::endl;
-    BOOST_CHECK(std::find(res.begin(), res.end(), *simplex) != res.end());
-  }
+  BOOST_CHECK(simplexTree.num_simplices() == num_simp);
 }
 
 BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
@@ -201,7 +183,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ FIRST
   std::cout << "   - INSERT 0" << std::endl;
-  typeVectorVertex firstSimplexVector { 0 };
+  typeVectorVertex firstSimplexVector{0};
   BOOST_CHECK(firstSimplexVector.size() == 1);
   typeSimplex firstSimplex = std::make_pair(firstSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
   typePairSimplexBool returnValue = st.insert_simplex(firstSimplex.first, firstSimplex.second);
@@ -212,7 +194,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ SECOND
   std::cout << "   - INSERT 1" << std::endl;
-  typeVectorVertex secondSimplexVector { 1 };
+  typeVectorVertex secondSimplexVector{1};
   BOOST_CHECK(secondSimplexVector.size() == 1);
   typeSimplex secondSimplex = std::make_pair(secondSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
   returnValue = st.insert_simplex(secondSimplex.first, secondSimplex.second);
@@ -223,7 +205,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ THIRD
   std::cout << "   - INSERT (0,1)" << std::endl;
-  typeVectorVertex thirdSimplexVector { 0, 1 };
+  typeVectorVertex thirdSimplexVector{0, 1};
   BOOST_CHECK(thirdSimplexVector.size() == 2);
   typeSimplex thirdSimplex = std::make_pair(thirdSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
   returnValue = st.insert_simplex(thirdSimplex.first, thirdSimplex.second);
@@ -234,7 +216,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ FOURTH
   std::cout << "   - INSERT 2" << std::endl;
-  typeVectorVertex fourthSimplexVector { 2 };
+  typeVectorVertex fourthSimplexVector{2};
   BOOST_CHECK(fourthSimplexVector.size() == 1);
   typeSimplex fourthSimplex = std::make_pair(fourthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
   returnValue = st.insert_simplex(fourthSimplex.first, fourthSimplex.second);
@@ -245,7 +227,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ FIFTH
   std::cout << "   - INSERT (2,0)" << std::endl;
-  typeVectorVertex fifthSimplexVector { 2, 0 };
+  typeVectorVertex fifthSimplexVector{2, 0};
   BOOST_CHECK(fifthSimplexVector.size() == 2);
   typeSimplex fifthSimplex = std::make_pair(fifthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
   returnValue = st.insert_simplex(fifthSimplex.first, fifthSimplex.second);
@@ -256,7 +238,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ SIXTH
   std::cout << "   - INSERT (2,1)" << std::endl;
-  typeVectorVertex sixthSimplexVector { 2, 1 };
+  typeVectorVertex sixthSimplexVector{2, 1};
   BOOST_CHECK(sixthSimplexVector.size() == 2);
   typeSimplex sixthSimplex = std::make_pair(sixthSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
   returnValue = st.insert_simplex(sixthSimplex.first, sixthSimplex.second);
@@ -267,7 +249,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ SEVENTH
   std::cout << "   - INSERT (2,1,0)" << std::endl;
-  typeVectorVertex seventhSimplexVector { 2, 1, 0 };
+  typeVectorVertex seventhSimplexVector{2, 1, 0};
   BOOST_CHECK(seventhSimplexVector.size() == 3);
   typeSimplex seventhSimplex = std::make_pair(seventhSimplexVector, Filtration_value(THIRD_FILTRATION_VALUE));
   returnValue = st.insert_simplex(seventhSimplex.first, seventhSimplex.second);
@@ -278,7 +260,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ EIGHTH
   std::cout << "   - INSERT 3" << std::endl;
-  typeVectorVertex eighthSimplexVector { 3 };
+  typeVectorVertex eighthSimplexVector{3};
   BOOST_CHECK(eighthSimplexVector.size() == 1);
   typeSimplex eighthSimplex = std::make_pair(eighthSimplexVector, Filtration_value(FIRST_FILTRATION_VALUE));
   returnValue = st.insert_simplex(eighthSimplex.first, eighthSimplex.second);
@@ -289,7 +271,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ NINETH
   std::cout << "   - INSERT (3,0)" << std::endl;
-  typeVectorVertex ninethSimplexVector { 3, 0 };
+  typeVectorVertex ninethSimplexVector{3, 0};
   BOOST_CHECK(ninethSimplexVector.size() == 2);
   typeSimplex ninethSimplex = std::make_pair(ninethSimplexVector, Filtration_value(SECOND_FILTRATION_VALUE));
   returnValue = st.insert_simplex(ninethSimplex.first, ninethSimplex.second);
@@ -300,7 +282,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ TENTH
   std::cout << "   - INSERT 0 (already inserted)" << std::endl;
-  typeVectorVertex tenthSimplexVector { 0 };
+  typeVectorVertex tenthSimplexVector{0};
   BOOST_CHECK(tenthSimplexVector.size() == 1);
   // With a different filtration value
   typeSimplex tenthSimplex = std::make_pair(tenthSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
@@ -315,7 +297,7 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 
   // ++ ELEVENTH
   std::cout << "   - INSERT (2,1,0) (already inserted)" << std::endl;
-  typeVectorVertex eleventhSimplexVector { 2, 1, 0 };
+  typeVectorVertex eleventhSimplexVector{2, 1, 0};
   BOOST_CHECK(eleventhSimplexVector.size() == 3);
   typeSimplex eleventhSimplex = std::make_pair(eleventhSimplexVector, Filtration_value(FOURTH_FILTRATION_VALUE));
   returnValue = st.insert_simplex(eleventhSimplex.first, eleventhSimplex.second);
@@ -378,14 +360,13 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 }
 
 BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
-  // TEST OF INSERTION
   std::cout << "********************************************************************" << std::endl;
-  std::cout << "TEST OF INSERTION" << std::endl;
+  std::cout << "TEST OF RECURSIVE INSERTION" << std::endl;
   typeST st;
 
   // ++ FIRST
   std::cout << "   - INSERT (2,1,0)" << std::endl;
-  typeVectorVertex SimplexVector1 { 2, 1, 0 };
+  typeVectorVertex SimplexVector1{2, 1, 0};
   BOOST_CHECK(SimplexVector1.size() == 3);
   st.insert_simplex_and_subfaces(SimplexVector1);
 
@@ -393,7 +374,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
 
   // ++ SECOND
   std::cout << "   - INSERT 3" << std::endl;
-  typeVectorVertex SimplexVector2 { 3 };
+  typeVectorVertex SimplexVector2{3};
   BOOST_CHECK(SimplexVector2.size() == 1);
   st.insert_simplex_and_subfaces(SimplexVector2);
 
@@ -401,7 +382,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
 
   // ++ THIRD
   std::cout << "   - INSERT (0,3)" << std::endl;
-  typeVectorVertex SimplexVector3 { 3, 0 };
+  typeVectorVertex SimplexVector3{3, 0};
   BOOST_CHECK(SimplexVector3.size() == 2);
   st.insert_simplex_and_subfaces(SimplexVector3);
 
@@ -409,7 +390,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
 
   // ++ FOURTH
   std::cout << "   - INSERT (1,0) (already inserted)" << std::endl;
-  typeVectorVertex SimplexVector4 { 1, 0 };
+  typeVectorVertex SimplexVector4{1, 0};
   BOOST_CHECK(SimplexVector4.size() == 2);
   st.insert_simplex_and_subfaces(SimplexVector4);
 
@@ -417,7 +398,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
 
   // ++ FIFTH
   std::cout << "   - INSERT (3,4,5)" << std::endl;
-  typeVectorVertex SimplexVector5 { 3, 4, 5 };
+  typeVectorVertex SimplexVector5{3, 4, 5};
   BOOST_CHECK(SimplexVector5.size() == 3);
   st.insert_simplex_and_subfaces(SimplexVector5);
 
@@ -425,7 +406,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
 
   // ++ SIXTH
   std::cout << "   - INSERT (0,1,6,7)" << std::endl;
-  typeVectorVertex SimplexVector6 { 0, 1, 6, 7 };
+  typeVectorVertex SimplexVector6{0, 1, 6, 7};
   BOOST_CHECK(SimplexVector6.size() == 4);
   st.insert_simplex_and_subfaces(SimplexVector6);
 
@@ -462,7 +443,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
   // ------------------------------------------------------------------------------------------------------------------
   // Find in the simplex_tree
   // ------------------------------------------------------------------------------------------------------------------
-  typeVectorVertex simpleSimplexVector { 1 };
+  typeVectorVertex simpleSimplexVector{1};
   Simplex_tree<>::Simplex_handle simplexFound = st.find(simpleSimplexVector);
   std::cout << "**************IS THE SIMPLEX {1} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != st.null_simplex())
@@ -472,7 +453,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
   // Check it is found
   BOOST_CHECK(simplexFound != st.null_simplex());
 
-  typeVectorVertex unknownSimplexVector { 15 };
+  typeVectorVertex unknownSimplexVector{15};
   simplexFound = st.find(unknownSimplexVector);
   std::cout << "**************IS THE SIMPLEX {15} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != st.null_simplex())
@@ -491,7 +472,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
   // Check it is found
   BOOST_CHECK(simplexFound != st.null_simplex());
 
-  typeVectorVertex otherSimplexVector { 1, 15 };
+  typeVectorVertex otherSimplexVector{1, 15};
   simplexFound = st.find(otherSimplexVector);
   std::cout << "**************IS THE SIMPLEX {15,1} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != st.null_simplex())
@@ -501,7 +482,7 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
   // Check it is NOT found
   BOOST_CHECK(simplexFound == st.null_simplex());
 
-  typeVectorVertex invSimplexVector { 1, 2, 0 };
+  typeVectorVertex invSimplexVector{1, 2, 0};
   simplexFound = st.find(invSimplexVector);
   std::cout << "**************IS THE SIMPLEX {1,2,0} IN THE SIMPLEX TREE ?\n";
   if (simplexFound != st.null_simplex())
@@ -522,143 +503,183 @@ BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
     }
     std::cout << std::endl;
   }
+}
+
+void test_cofaces(typeST& st, std::vector<Vertex_handle> expected, int dim, std::vector<typeST::Simplex_handle> res) {
+  typeST::Cofaces_simplex_range cofaces;
+  if (dim == 0)
+    cofaces = st.star_simplex_range(st.find(expected));
+  else
+    cofaces = st.cofaces_simplex_range(st.find(expected), dim);
+  for (auto simplex = cofaces.begin(); simplex != cofaces.end(); ++simplex) {
+    typeST::Simplex_vertex_range rg = st.simplex_vertex_range(*simplex);
+    for (auto vertex = rg.begin(); vertex != rg.end(); ++vertex) {
+      std::cout << "(" << *vertex << ")";
+    }
+    std::cout << std::endl;
+    BOOST_CHECK(std::find(res.begin(), res.end(), *simplex) != res.end());
+  }
+}
 
+BOOST_AUTO_TEST_CASE(coface_on_simplex_tree) {
   std::cout << "********************************************************************" << std::endl;
-  // TEST COFACE ALGORITHM
+  std::cout << "TEST COFACE ALGORITHM" << std::endl;
+  typeST st;
+
+  typeVectorVertex SimplexVector{2, 1, 0};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  SimplexVector = {3, 0};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  SimplexVector = {3, 4, 5};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  SimplexVector = {0, 1, 6, 7};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  /* Inserted simplex:        */
+  /*    1   6                 */
+  /*    o---o                 */
+  /*   /X\7/                  */
+  /*  o---o---o---o           */
+  /*  2   0   3\X/4           */
+  /*            o             */
+  /*            5             */
+
+  // FIXME
   st.set_dimension(3);
-  std::cout << "COFACE ALGORITHM" << std::endl;
-  std::vector<Vertex_handle> v;
-  std::vector<Vertex_handle> simplex;
+
+  std::vector<Vertex_handle> simplex_result;
   std::vector<typeST::Simplex_handle> result;
-  v.push_back(3);
-  std::cout << "First test : " << std::endl;
-  std::cout << "Star of (3):" << std::endl;
-
-  simplex.push_back(3);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(3);
-  simplex.push_back(0);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(4);
-  simplex.push_back(3);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(5);
-  simplex.push_back(4);
-  simplex.push_back(3);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(5);
-  simplex.push_back(3);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  test_cofaces(st, v, 0, result);
-  v.clear();
+  std::cout << "First test - Star of (3):" << std::endl;
+
+  simplex_result = {3};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {3, 0};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {4, 3};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {5, 4, 3};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {5, 3};
+  result.push_back(st.find(simplex_result));
+  simplex_result.clear();
+
+  std::vector<Vertex_handle> vertex = {3};
+  test_cofaces(st, vertex, 0, result);
+  vertex.clear();
   result.clear();
 
-  v.push_back(1);
-  v.push_back(7);
-  std::cout << "Second test : " << std::endl;
-  std::cout << "Star of (1,7): " << std::endl;
-
-  simplex.push_back(7);
-  simplex.push_back(1);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(7);
-  simplex.push_back(6);
-  simplex.push_back(1);
-  simplex.push_back(0);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(7);
-  simplex.push_back(1);
-  simplex.push_back(0);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  simplex.push_back(7);
-  simplex.push_back(6);
-  simplex.push_back(1);
-  result.push_back(st.find(simplex));
-  simplex.clear();
-
-  test_cofaces(st, v, 0, result);
+  vertex.push_back(1);
+  vertex.push_back(7);
+  std::cout << "Second test - Star of (1,7): " << std::endl;
+
+  simplex_result = {7, 1};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {7, 6, 1, 0};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {7, 1, 0};
+  result.push_back(st.find(simplex_result));
+
+  simplex_result = {7, 6, 1};
+  result.push_back(st.find(simplex_result));
+
+  test_cofaces(st, vertex, 0, result);
   result.clear();
 
-  std::cout << "Third test : " << std::endl;
-  std::cout << "2-dimension Cofaces of simplex(1,7) : " << std::endl;
+  std::cout << "Third test - 2-dimension Cofaces of simplex(1,7) : " << std::endl;
 
-  simplex.push_back(7);
-  simplex.push_back(1);
-  simplex.push_back(0);
-  result.push_back(st.find(simplex));
-  simplex.clear();
+  simplex_result = {7, 1, 0};
+  result.push_back(st.find(simplex_result));
 
-  simplex.push_back(7);
-  simplex.push_back(6);
-  simplex.push_back(1);
-  result.push_back(st.find(simplex));
-  simplex.clear();
+  simplex_result = {7, 6, 1};
+  result.push_back(st.find(simplex_result));
 
-  test_cofaces(st, v, 1, result);
+  test_cofaces(st, vertex, 1, result);
   result.clear();
 
   std::cout << "Cofaces with a codimension too high (codimension + vetices > tree.dimension) :" << std::endl;
-  test_cofaces(st, v, 5, result);
-  //    std::cout << "Cofaces with an empty codimension" << std::endl;
-  //    test_cofaces(st, v, -1, result);
+  test_cofaces(st, vertex, 5, result);
+
+  //std::cout << "Cofaces with an empty codimension" << std::endl;
+  //test_cofaces(st, vertex, -1, result);
   //    std::cout << "Cofaces in an empty simplex tree" << std::endl;
   //   typeST empty_tree;
-  //    test_cofaces(empty_tree, v, 1, result);
-  //    std::cout << "Cofaces of an empty simplex" << std::endl;
-  //    v.clear();
-  //    test_cofaces(st, v, 1, result);
-
-	std::cout << "********************************************************************" << std::endl;
-	// TEST Copy constructor / Move 
-	std::cout << "Printing st" << std::endl;
-	std::cout << &st << std::endl;
-	st.print_tree();
-	typeST st_copy_1 = st;
-	BOOST_CHECK(st == st_copy_1);
-	typeST st_move = std::move(st);
-	std::cout << "Printing a copy of st" << std::endl;
-	std::cout << &st_copy_1 << std::endl;
-	st_copy_1.print_tree();
-	BOOST_CHECK(st_move == st_copy_1);
-	std::cout << "Printing a move of st" << std::endl;
-	std::cout << &st_move << std::endl;
-	st_move.print_tree();
-	typeST st_empty;
-	BOOST_CHECK(st == st_empty);
-	std::cout << "Printing st again" << std::endl;
-	std::cout << &st << std::endl;
-	st.print_tree();
-
-	std::cout << "********************************************************************" << std::endl;
-	// TEST Edge_contraction
-	typeST st_copy_2 = st_copy_1, st_copy_3 = st_copy_1, st_copy_4 = st_copy_1;
-	st_copy_1.edge_contraction(0, 3);
-	std::cout << "Printing a copy of st, with the edge (0, 3) contracted, 3 being contracted in 0" << std::endl;
-	st_copy_1.print_tree();
-	st_copy_2.edge_contraction(1, 3);
-	std::cout << "Printing a copy of st, with the edge (1, 3) contracted, 3 being contracted in 1" << std::endl;
-	st_copy_2.print_tree();
-	st_copy_3.edge_contraction(3, 4);
-	std::cout << "Printing a copy of st, with the edge (3, 4) contracted, 4 being contracted in 3" << std::endl;
-	st_copy_3.print_tree();
-	st_copy_4.edge_contraction(1, 6);
-	std::cout << "Printing a copy of st, with the edge (1, 6) contracted, 6 being contracted in 1" << std::endl;
-	st_copy_4.print_tree();
+  //    test_cofaces(empty_tree, vertex, 1, result);
+  //std::cout << "Cofaces of an empty simplex" << std::endl;
+  //vertex.clear();
+  // test_cofaces(st, vertex, 1, result);
+
+}
 
+BOOST_AUTO_TEST_CASE(copy_move_on_simplex_tree) {
+  std::cout << "********************************************************************" << std::endl;
+  std::cout << "TEST COPY MOVE CONSTRUCTORS" << std::endl;
+  typeST st;
+
+  typeVectorVertex SimplexVector{2, 1, 0};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  SimplexVector = {3, 0};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  SimplexVector = {3, 4, 5};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  SimplexVector = {0, 1, 6, 7};
+  st.insert_simplex_and_subfaces(SimplexVector);
+
+  /* Inserted simplex:        */
+  /*    1   6                 */
+  /*    o---o                 */
+  /*   /X\7/                  */
+  /*  o---o---o---o           */
+  /*  2   0   3\X/4           */
+  /*            o             */
+  /*            5             */
+
+  // FIXME
+  st.set_dimension(3);
+
+  std::cout << "Printing st - address = " << &st << std::endl;
+  //std::cout << st << std::endl;
+
+  // Copy constructor  
+  typeST st_copy = st;
+  std::cout << "Printing a copy of st - address = " << &st_copy << std::endl;
+  //std::cout << st_copy << std::endl;
+
+  // Check the data are the same
+  BOOST_CHECK(st == st_copy);
+  // Check there is a new simplex tree reference
+  BOOST_CHECK(&st != &st_copy);
+
+  // Move constructor  
+  typeST st_move = std::move(st);
+  std::cout << "Printing a move of st - address = " << &st_move << std::endl;
+  //std::cout << st_move << std::endl;
+
+  // Check the data are the same
+  BOOST_CHECK(st_move == st_copy);
+  // Check there is a new simplex tree reference
+  BOOST_CHECK(&st_move != &st_copy);
+  BOOST_CHECK(&st_move != &st);
+  
+  typeST st_empty;
+  // Check st has been emptied by the move
+  BOOST_CHECK(st == st_empty);
+  BOOST_CHECK(st.filtration() == 0);
+  BOOST_CHECK(st.dimension() == -1);
+  BOOST_CHECK(st.num_simplices() == 0);
+  BOOST_CHECK(st.num_vertices() == (size_t)0);
+  
+  std::cout << "Printing st once again- address = " << &st << std::endl;
+  //std::cout << st << std::endl;
 }