Backmerge of trunk

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

Former-commit-id: 812e0b84d187dd5e30a6a18c12612ebc9bf9206a

10 files changed, 537 insertions, 269 deletions

diff --git a/src/Simplex_tree/concept/IndexingTag.h b/src/Simplex_tree/concept/IndexingTag.h
index d690da11..1dcdd756 100644
--- a/src/Simplex_tree/concept/IndexingTag.h
+++ b/src/Simplex_tree/concept/IndexingTag.h
@@ -25,6 +25,6 @@
   * continuous maps to a cell complex, and compute its persistent
   * homology.
   *
-  * Must be linear_indexing_tag. 
+  * Must be `Gudhi::linear_indexing_tag`. 
   */
 struct IndexingTag {};
diff --git a/src/Simplex_tree/concept/SimplexKey.h b/src/Simplex_tree/concept/SimplexKey.h
index ce5b2382..7fdbdd84 100644
--- a/src/Simplex_tree/concept/SimplexKey.h
+++ b/src/Simplex_tree/concept/SimplexKey.h
@@ -22,7 +22,7 @@
 
  /** \brief Key type used as simplex identifier.
   *
-  * Must be <CODE>int</CODE>
+  * Must be a signed integer type.
   */
  struct SimplexKey {};
- 
\ No newline at end of file
+ 
diff --git a/src/Simplex_tree/concept/SimplexTreeOptions.h b/src/Simplex_tree/concept/SimplexTreeOptions.h
new file mode 100644
index 00000000..a50a2bf1
--- /dev/null
+++ b/src/Simplex_tree/concept/SimplexTreeOptions.h
@@ -0,0 +1,41 @@
+ /*    This file is part of the Gudhi Library. The Gudhi library 
+  *    (Geometric Understanding in Higher Dimensions) is a generic C++ 
+  *    library for computational topology.
+  *
+  *    Author(s):       Marc Glisse
+  *
+  *    Copyright (C) 2015  INRIA Saclay - Ile-de-France (France)
+  *
+  *    This program is free software: you can redistribute it and/or modify
+  *    it under the terms of the GNU General Public License as published by
+  *    the Free Software Foundation, either version 3 of the License, or
+  *    (at your option) any later version.
+  *
+  *    This program is distributed in the hope that it will be useful,
+  *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+  *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  *    GNU General Public License for more details.
+  *
+  *    You should have received a copy of the GNU General Public License
+  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+  */
+  
+/** \brief Concept of the template parameter for the class `Gudhi::Simplex_tree<SimplexTreeOptions>`.
+ *
+ * One model for this is `Gudhi::Simplex_tree_options_full_featured`. If you want to provide your own, it is recommended that you derive from it and override some parts instead of writing a class from scratch.
+ */
+struct SimplexTreeOptions {
+  /// Forced for now.
+  typedef IndexingTag Indexing_tag;
+  /// Must be a signed integer type. It admits a total order <.
+  typedef VertexHandle Vertex_handle;
+  /// Must be comparable with operator<.
+  typedef FiltrationValue Filtration_value;
+  /// Must be a signed integer type.
+  typedef SimplexKey Simplex_key;
+  /// If true, each simplex has extra storage for one `Simplex_key`. Necessary for `Persistent_cohomology`.
+  static constexpr bool store_key;
+  /// If true, each simplex has extra storage for one `Filtration_value`, and this value is propagated by operations like `Gudhi::Simplex_tree<SimplexTreeOptions>::expansion`. Without it, `Persistent_cohomology` degenerates to computing usual (non-persistent) cohomology.
+  static constexpr bool store_filtration;
+};
+
diff --git a/src/Simplex_tree/concept/VertexHandle.h b/src/Simplex_tree/concept/VertexHandle.h
index 491f0f56..3efbba61 100644
--- a/src/Simplex_tree/concept/VertexHandle.h
+++ b/src/Simplex_tree/concept/VertexHandle.h
@@ -22,5 +22,6 @@
   
 /** \brief Handle type for the vertices of a cell complex.
   * 
-  * Must be int.*/
+  * Must be a signed integer type. <code>operator&lt;</code> defines a total order on it.
+  */
 struct VertexHandle {};
diff --git a/src/Simplex_tree/example/CMakeLists.txt b/src/Simplex_tree/example/CMakeLists.txt
index 1a3cdfbf..2f924490 100644
--- a/src/Simplex_tree/example/CMakeLists.txt
+++ b/src/Simplex_tree/example/CMakeLists.txt
@@ -7,15 +7,18 @@ add_test(simplex_tree_from_file_3 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_
 
 add_executable ( simple_simplex_tree simple_simplex_tree.cpp )
 add_test(simple_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/simple_simplex_tree)
-   
+
+add_executable ( mini_simplex_tree mini_simplex_tree.cpp )
+add_test(mini_simplex_tree ${CMAKE_CURRENT_BINARY_DIR}/mini_simplex_tree)
+
 
 # An example with Simplex-tree using CGAL alpha_shapes_3
 if(GMP_FOUND AND CGAL_FOUND)
    message("CGAL_lib = ${CGAL_LIBRARIES_DIR}")
-   message("GMP_LIBRARIES = ${GMP_LIBRARIES}")   
+   message("GMP_LIBRARIES = ${GMP_LIBRARIES}")
    INCLUDE_DIRECTORIES(${GMP_INCLUDE_DIR})
    INCLUDE_DIRECTORIES(${CGAL_INCLUDE_DIRS})
    add_executable ( simplex_tree_from_alpha_shapes_3 simplex_tree_from_alpha_shapes_3.cpp )
    target_link_libraries(simplex_tree_from_alpha_shapes_3 ${GMP_LIBRARIES} ${CGAL_LIBRARY})
    add_test(simplex_tree_from_alpha_shapes_3 ${CMAKE_CURRENT_BINARY_DIR}/simplex_tree_from_alpha_shapes_3 ${CMAKE_SOURCE_DIR}/data/points/bunny_5000)
-endif()  
+endif()
diff --git a/src/Simplex_tree/example/mini_simplex_tree.cpp b/src/Simplex_tree/example/mini_simplex_tree.cpp
new file mode 100644
index 00000000..08d626d3
--- /dev/null
+++ b/src/Simplex_tree/example/mini_simplex_tree.cpp
@@ -0,0 +1,68 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s):       Marc Glisse
+ *
+ *    Copyright (C) 2015  INRIA Saclay - Ile-de-France (France)
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/Simplex_tree.h>
+#include <iostream>
+#include <initializer_list>
+
+using namespace Gudhi;
+
+struct MyOptions : Simplex_tree_options_full_featured {
+  // Not doing persistence, so we don't need those
+  static const bool store_key = false;
+  static const bool store_filtration = false;
+  // I have few vertices
+  typedef short Vertex_handle;
+};
+typedef Simplex_tree<MyOptions> ST;
+
+// Dictionary should be private, but for now this is the easiest way.
+static_assert(sizeof(ST::Dictionary::value_type) < sizeof(Simplex_tree<>::Dictionary::value_type),
+    "Not storing the filtration and key should save some space");
+
+int main() {
+  ST st;
+
+  /* Complex to build. */
+  /*    1              */
+  /*    o              */
+  /*   /X\             */
+  /*  o---o---o        */
+  /*  2   0   3        */
+
+  auto triangle012 = {0, 1, 2};
+  auto edge03 = {0, 3};
+  st.insert_simplex_and_subfaces(triangle012);
+  st.insert_simplex_and_subfaces(edge03);
+  // FIXME: Remove this line
+  st.set_dimension(2);
+
+  auto edge02 = {0, 2};
+  ST::Simplex_handle e = st.find(edge02);
+  assert(st.filtration(e) == 0); // We are not using filtrations so everything has value 0
+  for(ST::Simplex_handle t : st.cofaces_simplex_range(e, 1)) // Only coface is 012
+    {
+      for(ST::Vertex_handle v : st.simplex_vertex_range(t)) // v in { 0, 1, 2 }
+        std::cout << v;
+      std::cout << '\n';
+    }
+}
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 279327f7..6a47083c 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -28,6 +28,9 @@
 #include <gudhi/Simplex_tree/Simplex_tree_iterators.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>
@@ -72,6 +75,17 @@ namespace Gudhi {
  * \copyright GNU General Public License v3.
  * @{
  */
+
+/// Model of SimplexTreeOptions.
+struct Simplex_tree_options_full_featured {
+  typedef linear_indexing_tag Indexing_tag;
+  typedef int Vertex_handle;
+  typedef double Filtration_value;
+  typedef int Simplex_key;
+  static const bool store_key = true;
+  static const bool store_filtration = true;
+};
+
 /**
  * \brief Simplex Tree data structure for representing simplicial complexes.
  *
@@ -83,35 +97,63 @@ namespace Gudhi {
  * \implements FilteredComplex
  *
  */
-template<typename IndexingTag = linear_indexing_tag,
-typename FiltrationValue = double, typename SimplexKey = int  // must be a signed integer type
-, typename VertexHandle = int  // must be a signed integer type, int convertible to it
->
+
+template<typename SimplexTreeOptions = Simplex_tree_options_full_featured>
 class Simplex_tree {
  public:
-  typedef IndexingTag Indexing_tag;
+  typedef SimplexTreeOptions Options;
+  typedef typename Options::Indexing_tag Indexing_tag;
   /** \brief Type for the value of the filtration function.
    *
    * Must be comparable with <. */
-  typedef FiltrationValue Filtration_value;
+  typedef typename Options::Filtration_value Filtration_value;
   /** \brief Key associated to each simplex.
    *
    * Must be a signed integer type. */
-  typedef SimplexKey Simplex_key;
+  typedef typename Options::Simplex_key Simplex_key;
   /** \brief Type for the vertex handle.
    *
    * Must be a signed integer type. It admits a total order <. */
-  typedef VertexHandle Vertex_handle;
+  typedef typename Options::Vertex_handle Vertex_handle;
 
   /* Type of node in the simplex tree. */
   typedef Simplex_tree_node_explicit_storage<Simplex_tree> Node;
   /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */
+  // Note: this wastes space when Vertex_handle is 32 bits and Node is aligned on 64 bits. It would be better to use a flat_set (with our own comparator) where we can control the layout of the struct (put Vertex_handle and Simplex_key next to each other).
   typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;
 
   /* \brief Set of nodes sharing a same parent in the simplex tree. */
   /* \brief Set of nodes sharing a same parent in the simplex tree. */
   typedef Simplex_tree_siblings<Simplex_tree, Dictionary> Siblings;
 
+  struct Key_simplex_base_real {
+    Key_simplex_base_real() : key_(-1) {}
+    void assign_key(Simplex_key k) { key_ = k; }
+    Simplex_key key() const { return key_; }
+    private:
+    Simplex_key key_;
+  };
+  struct Key_simplex_base_dummy {
+    Key_simplex_base_dummy() {}
+    void assign_key(Simplex_key) { }
+    Simplex_key key() const { assert(false); return -1; }
+  };
+  typedef typename std::conditional<Options::store_key, Key_simplex_base_real, Key_simplex_base_dummy>::type Key_simplex_base;
+
+  struct Filtration_simplex_base_real {
+    Filtration_simplex_base_real() : filt_(0) {}
+    void assign_filtration(Filtration_value f) { filt_ = f; }
+    Filtration_value filtration() const { return filt_; }
+    private:
+    Filtration_value filt_;
+  };
+  struct Filtration_simplex_base_dummy {
+    Filtration_simplex_base_dummy() {}
+    void assign_filtration(Filtration_value f) { assert(f == 0); }
+    Filtration_value filtration() const { return 0; }
+  };
+  typedef typename std::conditional<Options::store_filtration, Filtration_simplex_base_real, Filtration_simplex_base_dummy>::type Filtration_simplex_base;
+
  public:
   /** \brief Handle type to a simplex contained in the simplicial complex represented
    * by the simplex tree. */
@@ -170,12 +212,12 @@ class Simplex_tree {
   /** \brief Range over the simplices of the skeleton of the simplicial complex, for a given 
    * dimension. */
   typedef boost::iterator_range<Skeleton_simplex_iterator> Skeleton_simplex_range;
+  /** \brief Range over the simplices of the simplicial complex, ordered by the filtration. */
+  typedef std::vector<Simplex_handle> Filtration_simplex_range;
   /** \brief Iterator over the simplices of the simplicial complex, ordered by the filtration.
    *
    * 'value_type' is Simplex_handle. */
-  typedef typename std::vector<Simplex_handle>::iterator Filtration_simplex_iterator;
-  /** \brief Range over the simplices of the simplicial complex, ordered by the filtration. */
-  typedef boost::iterator_range<Filtration_simplex_iterator> Filtration_simplex_range;
+  typedef typename Filtration_simplex_range::const_iterator Filtration_simplex_iterator;
 
   /* @} */  // end name range and iterator types
   /** \name Range and iterator methods
@@ -226,17 +268,13 @@ class Simplex_tree {
    * order is used.
    *
    * The filtration must be valid. If the filtration has not been initialized yet, the
-   * method initializes it (i.e. order the simplices). */
-  Filtration_simplex_range filtration_simplex_range(Indexing_tag) {
+   * method initializes it (i.e. order the simplices). If the complex has changed since the last time the filtration
+   * was initialized, please call `initialize_filtration()` to recompute it. */
+  Filtration_simplex_range const& filtration_simplex_range(Indexing_tag=Indexing_tag()) {
     if (filtration_vect_.empty()) {
       initialize_filtration();
     }
-    return Filtration_simplex_range(filtration_vect_.begin(),
-                                    filtration_vect_.end());
-  }
-
-  Filtration_simplex_range filtration_simplex_range() {
-    return filtration_simplex_range(Indexing_tag());
+    return filtration_vect_;
   }
 
   /** \brief Returns a range over the vertices of a simplex.
@@ -278,9 +316,47 @@ class Simplex_tree {
   Simplex_tree()
       : null_vertex_(-1),
       threshold_(0),
-      root_(NULL, null_vertex_),
+      root_(nullptr, null_vertex_),
+      filtration_vect_(),
+      dimension_(-1) { }
+
+  /** \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_(),
-      dimension_(-1) {
+      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_source) {
+    for (auto sh = sib->members().begin(), sh_source = sib_source->members().begin();
+         sh != sib->members().end(); ++sh, ++sh_source) {
+      if (has_children(sh_source)) {
+        Siblings * newsib = new Siblings(sib, sh_source->first);
+        newsib->members_.reserve(sh_source->second.children()->members().size());
+        for (auto & child : sh_source->second.children()->members())
+          newsib->members_.emplace_hint(newsib->members_.end(), child.first, Node(sib, child.second.filtration()));
+        rec_copy(newsib, sh_source->second.children());
+        sh->second.assign_children(newsib);
+      }
+    }
+  }
+
+  /** \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. */
@@ -304,23 +380,63 @@ class Simplex_tree {
   }
 
  public:
+  /** \brief Checks if two simplex trees are equal. */
+  bool operator==(Simplex_tree& st2) {
+    if ((null_vertex_ != st2.null_vertex_) ||
+        (threshold_ != st2.threshold_) ||
+        (dimension_ != st2.dimension_))
+      return false;
+    return rec_equal(&root_, &st2.root_);
+  }
+
+  /** \brief Checks if two simplex trees are different. */
+  bool operator!=(Simplex_tree& st2) {
+    return (!(*this == st2));
+  }
+
+ private:
+  /** rec_equal: Checks recursively whether or not two simplex trees are equal, using depth first search. */
+  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() && sh2 != s2->members().end()); ++sh1, ++sh2) {
+      if (sh1->first != sh2->first || sh1->second.filtration() != sh2->second.filtration())
+        return false;
+      if (has_children(sh1) != has_children(sh2))
+        return false;
+      // Recursivity on children only if both have children
+      else if (has_children(sh1))
+        if (!rec_equal(sh1->second.children(), sh2->second.children()))
+          return false;
+    }
+    return true;
+  }
+
+ public:
   /** \brief Returns the key associated to a simplex.
    *
-   * The filtration must be initialized. */
+   * The filtration must be initialized.
+   * \pre SimplexTreeOptions::store_key
+   */
   static Simplex_key key(Simplex_handle sh) {
     return sh->second.key();
   }
 
   /** \brief Returns the simplex associated to a key.
    *
-   * The filtration must be initialized. */
+   * The filtration must be initialized.
+   * \pre SimplexTreeOptions::store_key
+   */
   Simplex_handle simplex(Simplex_key key) const {
     return filtration_vect_[key];
   }
 
   /** \brief Returns the filtration value of a simplex.
    *
-   * Called on the null_simplex, returns INFINITY. */
+   * Called on the null_simplex, returns INFINITY.
+   * If SimplexTreeOptions::store_filtration is false, returns 0.
+   */
   static Filtration_value filtration(Simplex_handle sh) {
     if (sh != null_simplex()) {
       return sh->second.filtration();
@@ -348,7 +464,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
@@ -395,7 +511,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();
     }
@@ -413,26 +529,34 @@ class Simplex_tree {
     return (sh->second.children()->parent() == sh->first);
   }
 
- public:
-  /** \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 RandomAccessVertexRange must be a range for which .begin() and
-   * .end() return random access iterators, with <CODE>value_type</CODE>
-   * <CODE>Vertex_handle</CODE>.
+   * The type InputVertexRange must be a range of <CODE>Vertex_handle</CODE>
+   * on which we can call std::begin() function
    */
-  template<class RandomAccessVertexRange>
-  Simplex_handle find(RandomAccessVertexRange & s) {
-    if (s.begin() == s.end())  // Empty simplex
-      return null_simplex();
-
-    sort(s.begin(), s.end());
+  template<class InputVertexRange>
+  Simplex_handle find(const InputVertexRange & s) {
+    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(const std::vector<Vertex_handle> & simplex) {
     Siblings * tmp_sib = &root_;
     Dictionary_it tmp_dit;
-    Vertex_handle last = s[s.size() - 1];
-    for (auto v : s) {
+    Vertex_handle last = simplex.back();
+    for (auto v : simplex) {
       tmp_dit = tmp_sib->members_.find(v);
       if (tmp_dit == tmp_sib->members_.end()) {
         return null_simplex();
@@ -450,43 +574,17 @@ class Simplex_tree {
   Simplex_handle find_vertex(Vertex_handle v) {
     return root_.members_.begin() + v;
   }
-  
-  /** \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
-   * @param[in]  filtration the filtration value assigned to the new 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 filed of the
-   * output pair to the Simplex_handle of the simplex. Otherwise, we set the Simplex_handle part to
-   * null_simplex.
-   *
-   * All subsimplices do not necessary need to be already in the simplex tree to proceed to an
-   * insertion. However, the property of being a simplicial complex will be violated. This allows
-   * us to insert a stream of simplices contained in a simplicial complex without considering any
-   * order on them.
-   *
-   * The filtration value
-   * assigned to the new simplex must preserve the monotonicity of the filtration.
-   *
-   * 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(RandomAccessVertexRange & simplex,
-                                                 Filtration_value filtration) {
-    if (simplex.empty()) {
-      return std::pair<Simplex_handle, bool>(null_simplex(), true);
-    }
-    // must be sorted in increasing order
-    sort(simplex.begin(), simplex.end());
+  //{ return root_.members_.find(v); }
 
+ private:
+  /** \brief Inserts 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));
@@ -508,24 +606,77 @@ class Simplex_tree {
     return res_insert;
   }
 
-  /** \brief Insert a N-simplex and all his subfaces, from a N-simplex represented by a range of
+ 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
+   * @param[in]  filtration the filtration value assigned to the new 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.
+   *
+   * All subsimplices do not necessary need to be already in the simplex tree to proceed to an
+   * insertion. However, the property of being a simplicial complex will be violated. This allows
+   * us to insert a stream of simplices contained in a simplicial complex without considering any
+   * order on them.
+   *
+   * The filtration value
+   * assigned to the new simplex must preserve the monotonicity of the filtration.
+   *
+   * The type InputVertexRange must be a range for which .begin() and
+   * .end() return input iterators, with 'value_type' Vertex_handle. */
+  template<class InputVertexRange>
+  std::pair<Simplex_handle, bool> insert_simplex(const InputVertexRange & simplex,
+                                                 Filtration_value filtration = 0) {
+    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
    * 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 RandomAccessVertexRange>
-  std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const RandomAccessVertexRange& Nsimplex,
-                                                              Filtration_value filtration = 0.0) {
-    // Simplex copy
-    std::vector<Vertex_handle> the_simplex(Nsimplex.begin(), Nsimplex.end());
-    // must be sorted in increasing order
-    std::sort(the_simplex.begin(), the_simplex.end());
+  template<class InputVertexRange>
+  std::pair<Simplex_handle, bool> insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
+                                   Filtration_value filtration = 0) {
+    auto first = std::begin(Nsimplex);
+    auto last = std::end(Nsimplex);
+
+    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));
+
     std::vector<std::vector<Vertex_handle>> to_be_inserted;
     std::vector<std::vector<Vertex_handle>> to_be_propagated;
-    return rec_insert_simplex_and_subfaces(the_simplex, to_be_inserted, to_be_propagated, filtration);
+    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,
@@ -534,15 +685,15 @@ class Simplex_tree {
     std::pair<Simplex_handle, bool> insert_result;
     if (the_simplex.size() > 1) {
       // Get and remove last vertex
-      Vertex_handle last_vertex= the_simplex.back();
+      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);
@@ -557,26 +708,23 @@ class Simplex_tree {
 
       // insert all to_be_inserted
       for (auto& simplex_tbi : to_be_inserted) {
-        insert_result = insert_simplex(simplex_tbi, filtration);
+        insert_result = insert_vertex_vector(simplex_tbi, 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)){
+      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_simplex(first_simplex, filtration);
+
+      insert_result = insert_vertex_vector(first_simplex, filtration);
     } else {
         std::cerr << "Simplex_tree::rec_insert_simplex_and_subfaces - Recursivity error" << std::endl;
         exit(-1);
     }
-    
     return insert_result;
   }
 
@@ -603,17 +751,6 @@ class Simplex_tree {
       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() {
@@ -778,8 +915,9 @@ class Simplex_tree {
         : st_(st) { }
 
     bool operator()(const Simplex_handle sh1, const Simplex_handle sh2) const {
-      if (st_->filtration(sh1) != st_->filtration(sh2)) {
-        return st_->filtration(sh1) < st_->filtration(sh2);
+      // Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
+      if (sh1->second.filtration() != sh2->second.filtration()) {
+	return sh1->second.filtration() < sh2->second.filtration();
       }
       // is sh1 a proper subface of sh2
       return st_->reverse_lexicographic_order(sh1, sh2);
@@ -925,7 +1063,7 @@ class Simplex_tree {
     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) {
@@ -947,7 +1085,7 @@ class Simplex_tree {
    * of the simplex, and fil is its filtration value. */
   void print_hasse(std::ostream& os) {
     os << num_simplices() << " " << std::endl;
-    for (auto sh : filtration_simplex_range(Indexing_tag())) {
+    for (auto sh : filtration_simplex_range()) {
       os << dimension(sh) << " ";
       for (auto b_sh : boundary_simplex_range(sh)) {
         os << key(b_sh) << " ";
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h
index 1f1a34cc..c49e30b9 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h
@@ -39,62 +39,31 @@ namespace Gudhi {
  * It stores explicitely its own filtration value and its own Simplex_key.
  */
 template<class SimplexTree>
-class Simplex_tree_node_explicit_storage {
- public:
+struct Simplex_tree_node_explicit_storage : SimplexTree::Filtration_simplex_base, SimplexTree::Key_simplex_base {
   typedef typename SimplexTree::Siblings Siblings;
   typedef typename SimplexTree::Filtration_value Filtration_value;
   typedef typename SimplexTree::Simplex_key Simplex_key;
 
-  // Default constructor.
-  Simplex_tree_node_explicit_storage()
-      : children_(NULL),
-        simplex_key_(-1),
-        filtration_(0) {
-  }
-
-  Simplex_tree_node_explicit_storage(Siblings * sib,
-                                     Filtration_value filtration)
-      : children_(sib),
-        simplex_key_(-1),
-        filtration_(filtration) {
-  }
-
-  void assign_key(Simplex_key key) {
-    simplex_key_ = key;
+  Simplex_tree_node_explicit_storage(Siblings * sib = nullptr,
+                                     Filtration_value filtration = 0)
+      : children_(sib) {
+	  this->assign_filtration(filtration);
   }
 
   /*
-   * Assign a children to the node
+   * Assign children to the node
    */
   void assign_children(Siblings * children) {
     children_ = children;
   }
-  /*
-   *
-   */
-  void assign_filtration(double filtration_value) {
-    filtration_ = filtration_value;
-  }
-
-  Filtration_value filtration() {
-    return filtration_;
-  }
 
   /* Careful -> children_ can be NULL*/
   Siblings * children() {
     return children_;
   }
 
-  Simplex_key key() {
-    return simplex_key_;
-  }
-
  private:
   Siblings * children_;
-
-  // Data attached to simplex, explicit storage
-  Simplex_key simplex_key_;
-  Filtration_value filtration_;   // value in the filtration
 };
 
 /* @} */  // end addtogroup simplex_tree
diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h
index de350f2d..d20a91d7 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree/Simplex_tree_siblings.h
@@ -71,15 +71,15 @@ class Simplex_tree_siblings {
   /* \brief Constructor with initialized set of members.
    *
    * 'members' must be sorted and unique.*/
-  Simplex_tree_siblings(Simplex_tree_siblings * oncles, Vertex_handle parent,
-                        const std::vector<std::pair<Vertex_handle, Node> > & members)
+  template<typename RandomAccessVertexRange>
+  Simplex_tree_siblings(Simplex_tree_siblings * oncles, Vertex_handle parent, const RandomAccessVertexRange & members)
       : oncles_(oncles),
         parent_(parent),
         members_(boost::container::ordered_unique_range, members.begin(),
                  members.end()) {
     for (auto& map_el : members_) {
       map_el.second.assign_children(this);
-    }
+	}
   }
 
   /*
diff --git a/src/Simplex_tree/test/simplex_tree_unit_test.cpp b/src/Simplex_tree/test/simplex_tree_unit_test.cpp
index 9340aaa3..a4871cfd 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/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;
@@ -126,6 +125,7 @@ void test_simplex_tree_contains(typeST& simplexTree, typeSimplex& simplex, int p
   BOOST_CHECK(AreAlmostTheSame(simplexTree.filtration(*f_simplex), simplex.second));
 
   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)) {
     std::cout << "test_simplex_tree_contains - vertex=" << vertex << "||" << simplex.first.at(simplexIndex) << std::endl;
     BOOST_CHECK(vertex == simplex.first.at(simplexIndex));
@@ -170,22 +170,6 @@ void set_and_test_simplex_tree_dim_fil(typeST& simplexTree, int vectorSize, cons
   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_AUTO_TEST_CASE(simplex_tree_insertion) {
   const Filtration_value FIRST_FILTRATION_VALUE = 0.1;
   const Filtration_value SECOND_FILTRATION_VALUE = 0.2;
@@ -378,9 +362,8 @@ BOOST_AUTO_TEST_CASE(simplex_tree_insertion) {
 bool sort_in_decr_order (Vertex_handle i,Vertex_handle j) { return (i>j); }
 
 BOOST_AUTO_TEST_CASE(NSimplexAndSubfaces_tree_insertion) {
-  // TEST OF INSERTION WITH SUBFACES
   std::cout << "********************************************************************" << std::endl;
-  std::cout << "TEST OF INSERTION WITH SUBFACES" << std::endl;
+  std::cout << "TEST OF RECURSIVE INSERTION" << std::endl;
   typeST st;
   typePairSimplexBool returnValue;
   int position = 0;
@@ -569,114 +552,179 @@ 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);
-
-  /*
-     // TEST Off read
-     std::cout << "********************************************************************" << std::endl;
-     typeST st2;
-     st2.tree_from_off("test.off");
-     std::cout << st2;
-   */
+  //    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;
+
+  // Copy constructor  
+  typeST st_copy = st;
+  std::cout << "Printing a copy of st - address = " << &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;
+
+  // 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;
 }