Multiple fixes

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

Former-commit-id: 89e8ccf59283c5775db2f058bf5439bd487d79aa

2 files changed, 287 insertions, 215 deletions

diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 5eef283e..57e9b987 100644
--- a/src/Simplex_tree/include/gudhi/Simplex_tree.h
+++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h
@@ -20,8 +20,8 @@
  *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-#ifndef SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
-#define SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
+#ifndef SIMPLEX_TREE_H_
+#define SIMPLEX_TREE_H_
 
 #include <gudhi/Simplex_tree/Simplex_tree_node_explicit_storage.h>
 #include <gudhi/Simplex_tree/Simplex_tree_siblings.h>
@@ -37,9 +37,10 @@
 #include <algorithm>
 #include <utility>
 #include <vector>
+#include <functional>  // for greater<>
 
 namespace Gudhi {
-  
+
 /** \defgroup simplex_tree Filtered Complexes
  *
  * A simplicial complex \f$\mathbf{K}\f$
@@ -74,6 +75,7 @@ namespace Gudhi {
  * \copyright GNU General Public License v3.
  * @{
  */
+
 /**
  * \brief Simplex Tree data structure for representing simplicial complexes.
  *
@@ -86,9 +88,8 @@ namespace Gudhi {
  *
  */
 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
-//         , bool ContiguousVertexHandles = true   //true is Vertex_handles are exactly the set [0;n)
+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
 >
 class Simplex_tree {
  public:
@@ -111,20 +112,13 @@ class Simplex_tree {
   /* Type of dictionary Vertex_handle -> Node for traversing the simplex tree. */
   typedef typename boost::container::flat_map<Vertex_handle, Node> Dictionary;
 
-  friend class Simplex_tree_node_explicit_storage< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
-  friend class Simplex_tree_siblings< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle>, Dictionary>;
-  friend class Simplex_tree_simplex_vertex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
-  friend class Simplex_tree_boundary_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
-  friend class Simplex_tree_complex_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
-  friend class Simplex_tree_skeleton_simplex_iterator< Simplex_tree<FiltrationValue, SimplexKey, VertexHandle> >;
-
   /* \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;
 
  public:
   /** \brief Handle type to a simplex contained in the simplicial complex represented
-   * byt he simplex tree. */
+   * by the simplex tree. */
   typedef typename Dictionary::iterator Simplex_handle;
 
  private:
@@ -158,6 +152,8 @@ class Simplex_tree {
   typedef Simplex_tree_simplex_vertex_iterator<Simplex_tree> Simplex_vertex_iterator;
   /** \brief Range over the vertices of a simplex. */
   typedef boost::iterator_range<Simplex_vertex_iterator> Simplex_vertex_range;
+  /** \brief Range over the cofaces of a simplex. */
+  typedef std::vector<Simplex_handle> Cofaces_simplex_range;
   /** \brief Iterator over the simplices of the boundary of a simplex.
    *
    * 'value_type' is Simplex_handle. */
@@ -189,13 +185,12 @@ class Simplex_tree {
   /** \name Range and iterator methods
    * @{ */
 
-  /** \brief Returns a range over the vertices of the simplicial complex.
-   *
+  /** \brief Returns a range over the vertices of the simplicial complex. 
    * The order is increasing according to < on Vertex_handles.*/
   Complex_vertex_range complex_vertex_range() {
     return Complex_vertex_range(
-        boost::make_transform_iterator(root_.members_.begin(), return_first()),
-        boost::make_transform_iterator(root_.members_.end(), return_first()));
+                                boost::make_transform_iterator(root_.members_.begin(), return_first()),
+                                boost::make_transform_iterator(root_.members_.end(), return_first()));
   }
 
   /** \brief Returns a range over the simplices of the simplicial complex.
@@ -246,6 +241,7 @@ class Simplex_tree {
   Filtration_simplex_range filtration_simplex_range() {
     return filtration_simplex_range(Indexing_tag());
   }
+
   /** \brief Returns a range over the vertices of a simplex.
    *
    * The order in which the vertices are visited is the decreasing order for < on Vertex_handles,
@@ -253,6 +249,7 @@ class Simplex_tree {
    * equal to \f$(-1)^{\text{dim} \sigma}\f$ the canonical orientation on the simplex.
    */
   Simplex_vertex_range simplex_vertex_range(Simplex_handle sh) {
+    assert(sh != null_simplex());  // Empty simplex
     return Simplex_vertex_range(Simplex_vertex_iterator(this, sh),
                                 Simplex_vertex_iterator(this));
   }
@@ -283,12 +280,11 @@ class Simplex_tree {
   /** \brief Constructs an empty simplex tree. */
   Simplex_tree()
       : null_vertex_(-1),
-        threshold_(0),
-        num_simplices_(0),
-        root_(NULL, null_vertex_),
-        filtration_vect_(),
-        dimension_(-1) {
-  }
+      threshold_(0),
+      num_simplices_(0),
+      root_(NULL, 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_), num_simplices_(copy.num_simplices_), root_(NULL, -1, std::vector<std::pair<Vertex_handle, Node>> (copy.root_.members().begin(), copy.root_.members().end())), filtration_vect_(copy.filtration_vect_), dimension_(copy.dimension_) {
@@ -330,7 +326,7 @@ class Simplex_tree {
   }
   /** @} */  // end constructor/destructor
  private:
-  /** Recursive deletion. */
+  // Recursive deletion
   void rec_delete(Siblings * sib) {
     for (auto sh = sib->members().begin(); sh != sib->members().end(); ++sh) {
       if (has_children(sh)) {
@@ -379,20 +375,11 @@ class Simplex_tree {
   /** \brief Checks whether or not two simplex trees are equal. */
   bool operator ==(Simplex_tree st2)
   {
-	  if (root_.members().size() != st2.root()->members().size() || this->null_vertex_ != st2.null_vertex_ || this->threshold_ != st2.threshold_ || this->num_simplices_ != st2.num_simplices_ || this->dimension_ != st2.dimension_)
-		  return false;
-	  for (auto sh1 = root_.members().begin(), sh2 = st2.root()->members().begin(); sh1 != root_.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;
+	  return (this->null_vertex_ == st2.null_vertex_ 
+	  		&& this->threshold_ == st2.threshold_ 
+			&& this->num_simplices_ == st2.num_simplices_ 
+			&& this->dimension_ == st2.dimension_ 
+			&& rec_equal(&root_, &st2.root_));
   }
 
   /** rec_equal: Checks recursively whether or not two simplex trees are equal, using depth first search. */
@@ -419,100 +406,109 @@ class Simplex_tree {
   /** \brief Returns the key associated to a simplex.
    *
    * The filtration must be initialized. */
-  Simplex_key key(Simplex_handle sh) {
+  static Simplex_key key(Simplex_handle sh) {
 	  return sh->second.key();
   }
+
   /** \brief Returns the simplex associated to a key.
    *
    * The filtration must be initialized. */
-  Simplex_handle simplex(Simplex_key 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. */
-  Filtration_value filtration(Simplex_handle sh) {
-	  if (sh != null_simplex()) {
-		  return sh->second.filtration();
-	  } else {
-		  return INFINITY;
-	  }  // filtration(); }
-}
-/** \brief Returns an upper bound of the filtration values of the simplices. */
-Filtration_value filtration() {
-	return threshold_;
-}
-/** \brief Returns a Simplex_handle different from all Simplex_handles
- * associated to the simplices in the simplicial complex.
- *
- * One can call filtration(null_simplex()). */
-Simplex_handle null_simplex() {
-	return Dictionary_it(NULL);
-}
-/** \brief Returns a key different for all keys associated to the
- * simplices of the simplicial complex. */
-Simplex_key null_key() {
-	return -1;
-}
-/** \brief Returns a Vertex_handle different from all Vertex_handles associated
- * to the vertices of the simplicial complex. */
-Vertex_handle null_vertex() {
-	return null_vertex_;
-}
-/** \brief Returns the number of vertices in the complex. */
-size_t num_vertices() {
-	return root_.members_.size();
-}
-/** \brief Returns the number of simplices in the complex.
- *
- * Does not count the empty simplex. */
-const unsigned int& num_simplices() const {
-	return num_simplices_;
-}
+  static Filtration_value filtration(Simplex_handle sh) {
+    if (sh != null_simplex()) {
+      return sh->second.filtration();
+    } else {
+      return INFINITY;
+    }
+  }
 
-/** \brief Returns the dimension of a simplex.
- *
- * Must be different from null_simplex().*/
-int dimension(Simplex_handle sh) {
-	Siblings * curr_sib = self_siblings(sh);
-	int dim = 0;
-	while (curr_sib != NULL) {
-		++dim;
-		curr_sib = curr_sib->oncles();
-	}
-	return dim - 1;
-}
-/** \brief Returns an upper bound on the dimension of the simplicial complex. */
-int dimension() {
-	return dimension_;
-}
+  /** \brief Returns an upper bound of the filtration values of the simplices. */
+  Filtration_value filtration() const {
+    return threshold_;
+  }
 
-/** \brief Returns true iff the node in the simplex tree pointed by
- * sh has children.*/
-bool has_children(Simplex_handle sh) {
-	return (sh->second.children()->parent() == sh->first);
-}
+  /** \brief Returns a Simplex_handle different from all Simplex_handles
+   * associated to the simplices in the simplicial complex.
+   *
+   * One can call filtration(null_simplex()). */
+  static Simplex_handle null_simplex() {
+    return Dictionary_it(NULL);
+  }
 
-/** \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
- */
+  /** \brief Returns a key different for all keys associated to the
+   * simplices of the simplicial complex. */
+  static Simplex_key null_key() {
+    return -1;
+  }
+
+  /** \brief Returns a Vertex_handle different from all Vertex_handles associated
+   * to the vertices of the simplicial complex. */
+  Vertex_handle null_vertex() const {
+    return null_vertex_;
+  }
+
+  /** \brief Returns the number of vertices in the complex. */
+  size_t num_vertices() const {
+    return root_.members_.size();
+  }
+
+  /** \brief Returns the number of simplices in the complex.
+   *
+   * Does not count the empty simplex. */
+  unsigned int num_simplices() const {
+    return num_simplices_;
+  }
+
+  /** \brief Returns the dimension of a simplex.
+   *
+   * Must be different from null_simplex().*/
+  int dimension(Simplex_handle sh) {
+    Siblings * curr_sib = self_siblings(sh);
+    int dim = 0;
+    while (curr_sib != NULL) {
+      ++dim;
+      curr_sib = curr_sib->oncles();
+    }
+    return dim - 1;
+  }
+
+  /** \brief Returns an upper bound on the dimension of the simplicial complex. */
+  int dimension() const {
+    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
+   * 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>
   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_rec(copy);
+	return find_simplex(copy);
   }
 
  private:
-  /** recursive function of find */
-  Simplex_handle find_rec(std::vector<Vertex_handle> & simplex) {
+  /** 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];
@@ -534,13 +530,13 @@ template<class InputVertexRange>
   Simplex_handle find_vertex(Vertex_handle v) {
     return root_.members_.begin() + v;
   }
-//{ return root_.members_.find(v); }
+  //{ 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) {
+                                                 Filtration_value filtration) {
     if (simplex.empty()) {
       return std::pair<Simplex_handle, bool>(null_simplex(), true);
     }
@@ -555,12 +551,15 @@ template<class InputVertexRange>
         curr_sib = res_insert.first->second.children();
     }
     res_insert = curr_sib->members_.emplace(*vi, Node(curr_sib, filtration));
-    if (!res_insert.second) {  // if already in the complex
-      if (res_insert.first->second.filtration() > filtration) {  // if filtration value modified
+    if (!res_insert.second) {
+      // if already in the complex
+      if (res_insert.first->second.filtration() > filtration) {
+        // if filtration value modified
         res_insert.first->second.assign_filtration(filtration);
         return res_insert;
       }
-      return std::pair<Simplex_handle, bool>(null_simplex(), false);  // if filtration value unchanged
+      // if filtration value unchanged
+      return std::pair<Simplex_handle, bool>(null_simplex(), false);
     }
     // otherwise the insertion has succeeded
     return res_insert;
@@ -602,10 +601,10 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
    *
    * @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.
-  */
+   */
   template<class InputVertexRange>
   void insert_simplex_and_subfaces(const InputVertexRange& Nsimplex,
-                                         Filtration_value filtration = 0.0) {
+                                   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);
@@ -624,7 +623,7 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
         RandomAccessVertexRange NsimplexMinusOne;
         for (unsigned int NListIter = 0; NListIter < Nsimplex.size() - 1; NListIter++) {
           // (N-1)-Simplex creation
-          NsimplexMinusOne.push_back( Nsimplex[(NIndex + NListIter) % Nsimplex.size()]);
+          NsimplexMinusOne.push_back(Nsimplex[(NIndex + NListIter) % Nsimplex.size()]);
         }
         // (N-1)-Simplex recursive call
         insert_simplex_and_subfaces(NsimplexMinusOne, filtration);
@@ -658,9 +657,8 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
    * and edge. sh must point to a 1-dimensional simplex. This is an
    * optimized version of the boundary computation. */
   std::pair<Simplex_handle, Simplex_handle> endpoints(Simplex_handle sh) {
-    return std::pair<Simplex_handle, Simplex_handle>(
-        root_.members_.find(sh->first),
-        root_.members_.find(self_siblings(sh)->parent()));
+    assert(dimension(sh) == 1);
+    return { find_vertex(sh->first), find_vertex(self_siblings(sh)->parent()) };
   }
 
   /** Returns the Siblings containing a simplex.*/
@@ -671,12 +669,12 @@ 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 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 << " ";}
@@ -688,15 +686,16 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
     return &root_;
   }
 
- public:
   /** Set an upper bound for the filtration values. */
   void set_filtration(Filtration_value fil) {
     threshold_ = fil;
   }
+
   /** Set a number of simplices for the simplicial complex. */
-  void set_num_simplices(const unsigned int& num_simplices) {
+  void set_num_simplices(unsigned int num_simplices) {
     num_simplices_ = num_simplices;
   }
+
   /** Set a dimension for the simplicial complex. */
   void set_dimension(int dimension) {
     dimension_ = dimension;
@@ -721,17 +720,15 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
   void initialize_filtration() {
     filtration_vect_.clear();
     filtration_vect_.reserve(num_simplices());
-    for (auto cpx_it = complex_simplex_range().begin();
-        cpx_it != complex_simplex_range().end(); ++cpx_it) {
-      filtration_vect_.push_back(*cpx_it);
-    }
+    for (Simplex_handle sh : complex_simplex_range())
+      filtration_vect_.push_back(sh);
 
     // the stable sort ensures the ordering heuristic
     std::stable_sort(filtration_vect_.begin(), filtration_vect_.end(),
                      is_before_in_filtration(this));
   }
 
-/** \brief Contracts two edges : the contracted edge is erased from the three, and the remaining edge receives all the links of the contracted one.
+/** \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
 */
@@ -779,7 +776,6 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
 			std::pair<Vertex_handle, Node> member(sh->first, Node(sibTarget, sh->second.filtration()));
 			if (has_children(sh))
 			{
-				std::cout << "children" << std::endl;
 				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)
@@ -808,6 +804,92 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
 	}
 
  private:
+  /** Recursive search of cofaces
+   * This function uses DFS
+   *\param vertices contains a list of vertices, which represent the vertices of the simplex not found yet.
+   *\param curr_nbVertices represents the number of vertices of the simplex we reached by going through the tree.
+   *\param cofaces contains a list of Simplex_handle, representing all the cofaces asked.
+   *\param star true if we need the star of the simplex
+   *\param nbVertices number of vertices of the cofaces we search
+   * Prefix actions : When the bottom vertex matches with the current vertex in the tree, we remove the bottom vertex from vertices.
+   * Infix actions : Then we call or not the recursion.
+   * Postfix actions : Finally, we add back the removed vertex into vertices, and remove this vertex from curr_nbVertices so that we didn't change the parameters.
+   * If the vertices list is empty, we need to check if curr_nbVertices matches with the dimension of the cofaces asked.
+   */
+  void rec_coface(std::vector<Vertex_handle> &vertices, Siblings *curr_sib, int curr_nbVertices,
+                  std::vector<Simplex_handle>& cofaces, bool star, int nbVertices) {
+    if (!(star || curr_nbVertices <= nbVertices))  // dimension of actual simplex <= nbVertices
+      return;
+    for (Simplex_handle simplex = curr_sib->members().begin(); simplex != curr_sib->members().end(); ++simplex) {
+      if (vertices.empty()) {
+        // If we reached the end of the vertices, and the simplex has more vertices than the given simplex
+        // => we found a coface
+
+        // Add a coface if we wan't the star or if the number of vertices of the current simplex matches with nbVertices
+        bool addCoface = (star || curr_nbVertices == nbVertices);
+        if (addCoface)
+          cofaces.push_back(simplex);
+        if ((!addCoface || star) && has_children(simplex))  // Rec call
+          rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
+      } else {
+        if (simplex->first == vertices.back()) {
+          // If curr_sib matches with the top vertex
+          bool equalDim = (star || curr_nbVertices == nbVertices);  // dimension of actual simplex == nbVertices
+          bool addCoface = vertices.size() == 1 && equalDim;
+          if (addCoface)
+            cofaces.push_back(simplex);
+          if ((!addCoface || star) && has_children(simplex)) {
+            // Rec call
+            Vertex_handle tmp = vertices.back();
+            vertices.pop_back();
+            rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
+            vertices.push_back(tmp);
+          }
+        } else if (simplex->first > vertices.back()) {
+          return;
+        } else {
+          // (simplex->first < vertices.back()
+          if (has_children(simplex))
+            rec_coface(vertices, simplex->second.children(), curr_nbVertices + 1, cofaces, star, nbVertices);
+        }
+      }
+    }
+  }
+
+ public:
+  /** \brief Compute the star of a n simplex
+   * \param simplex represent the simplex of which we search the star
+   * \return Vector of Simplex_handle, empty vector if no cofaces found.
+   */
+
+  Cofaces_simplex_range star_simplex_range(const Simplex_handle simplex) {
+    return cofaces_simplex_range(simplex, 0);
+  }
+
+  /** \brief Compute the cofaces of a n simplex
+   * \param simplex represent the n-simplex of which we search the n+codimension cofaces
+   * \param codimension The function returns the n+codimension-cofaces of the n-simplex. If codimension = 0, 
+   * return all cofaces (equivalent of star function)
+   * \return Vector of Simplex_handle, empty vector if no cofaces found.
+   */
+
+  Cofaces_simplex_range cofaces_simplex_range(const Simplex_handle simplex, int codimension) {
+    Cofaces_simplex_range cofaces;
+    // codimension must be positive or null integer
+    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_
+      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()));
+    return cofaces;
+  }
+
+ private:
   /** \brief Returns true iff the list of vertices of sh1
    * is smaller than the list of vertices of sh2 w.r.t.
    * lexicographic order on the lists read in reverse.
@@ -830,6 +912,7 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
     }
     return ((it1 == rg1.end()) && (it2 != rg2.end()));
   }
+
   /** \brief StrictWeakOrdering, for the simplices, defined by the filtration.
    *
    * It corresponds to the partial order
@@ -838,15 +921,14 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
    * to be smaller. The filtration function must be monotonic. */
   struct is_before_in_filtration {
     explicit is_before_in_filtration(Simplex_tree * st)
-        : st_(st) {
-    }
+        : 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);
       }
-
-      return st_->reverse_lexicographic_order(sh1, sh2);  // is sh1 a proper subface of sh2
+      // is sh1 a proper subface of sh2
+      return st_->reverse_lexicographic_order(sh1, sh2);
     }
 
     Simplex_tree * st_;
@@ -873,7 +955,8 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
    *                                    must be undirected_tag. */
   template<class OneSkeletonGraph>
   void insert_graph(const OneSkeletonGraph& skel_graph) {
-    assert(num_simplices() == 0);  // the simplex tree must be empty
+    // the simplex tree must be empty
+    assert(num_simplices() == 0);
 
     if (boost::num_vertices(skel_graph) == 0) {
       return;
@@ -891,30 +974,32 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
     typename boost::graph_traits<OneSkeletonGraph>::vertex_iterator v_it,
         v_it_end;
     for (std::tie(v_it, v_it_end) = boost::vertices(skel_graph); v_it != v_it_end;
-        ++v_it) {
+         ++v_it) {
       root_.members_.emplace_hint(
-          root_.members_.end(), *v_it,
-          Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
+                                  root_.members_.end(), *v_it,
+                                  Node(&root_, boost::get(vertex_filtration_t(), skel_graph, *v_it)));
     }
     typename boost::graph_traits<OneSkeletonGraph>::edge_iterator e_it,
         e_it_end;
     for (std::tie(e_it, e_it_end) = boost::edges(skel_graph); e_it != e_it_end;
-        ++e_it) {
+         ++e_it) {
       auto u = source(*e_it, skel_graph);
       auto v = target(*e_it, skel_graph);
-      if (u < v) {  // count edges only once { std::swap(u,v); } // u < v
+      if (u < v) {
+        // count edges only once { std::swap(u,v); } // u < v
         auto sh = find_vertex(u);
         if (!has_children(sh)) {
           sh->second.assign_children(new Siblings(&root_, sh->first));
         }
 
         sh->second.children()->members().emplace(
-            v,
-            Node(sh->second.children(),
-                 boost::get(edge_filtration_t(), skel_graph, *e_it)));
+                                                 v,
+                                                 Node(sh->second.children(),
+                                                      boost::get(edge_filtration_t(), skel_graph, *e_it)));
       }
     }
   }
+
   /** \brief Expands the Simplex_tree containing only its one skeleton
    * until dimension max_dim.
    *
@@ -929,7 +1014,7 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
   void expansion(int max_dim) {
     dimension_ = max_dim;
     for (Dictionary_it root_it = root_.members_.begin();
-        root_it != root_.members_.end(); ++root_it) {
+         root_it != root_.members_.end(); ++root_it) {
       if (has_children(root_it)) {
         siblings_expansion(root_it->second.children(), max_dim - 1);
       }
@@ -940,7 +1025,7 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
  private:
   /** \brief Recursive expansion of the simplex tree.*/
   void siblings_expansion(Siblings * siblings,  // must contain elements
-      int k) {
+                          int k) {
     if (dimension_ > k) {
       dimension_ = k;
     }
@@ -951,68 +1036,56 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
 
     static std::vector<std::pair<Vertex_handle, Node> > inter;  // static, not thread-safe.
     for (Dictionary_it s_h = siblings->members().begin();
-        s_h != siblings->members().end(); ++s_h, ++next) {
+         s_h != siblings->members().end(); ++s_h, ++next) {
       Simplex_handle root_sh = find_vertex(s_h->first);
       if (has_children(root_sh)) {
         intersection(
-            inter,                      // output intersection
-            next,                       // begin
-            siblings->members().end(),  // end
-            root_sh->second.children()->members().begin(),
-            root_sh->second.children()->members().end(),
-            s_h->second.filtration());
+                     inter,  // output intersection
+                     next,  // begin
+                     siblings->members().end(),  // end
+                     root_sh->second.children()->members().begin(),
+                     root_sh->second.children()->members().end(),
+                     s_h->second.filtration());
         if (inter.size() != 0) {
           this->num_simplices_ += inter.size();
           Siblings * new_sib = new Siblings(siblings,  // oncles
-              s_h->first,  // parent
-              inter);      // boost::container::ordered_unique_range_t
+                                            s_h->first,  // parent
+                                            inter);  // boost::container::ordered_unique_range_t
           inter.clear();
           s_h->second.assign_children(new_sib);
           siblings_expansion(new_sib, k - 1);
         } else {
-          s_h->second.assign_children(siblings);  // ensure the children property
+          // ensure the children property
+          s_h->second.assign_children(siblings);
           inter.clear();
         }
       }
     }
   }
+
   /** \brief Intersects Dictionary 1 [begin1;end1) with Dictionary 2 [begin2,end2)
    * and assigns the maximal possible Filtration_value to the Nodes. */
-  void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
-                    Dictionary_it begin1, Dictionary_it end1,
-                    Dictionary_it begin2, Dictionary_it end2,
-                    Filtration_value filtration) {
+  static void intersection(std::vector<std::pair<Vertex_handle, Node> >& intersection,
+                           Dictionary_it begin1, Dictionary_it end1,
+                           Dictionary_it begin2, Dictionary_it end2,
+                           Filtration_value filtration_) {
     if (begin1 == end1 || begin2 == end2)
       return;  // ----->>
     while (true) {
       if (begin1->first == begin2->first) {
-        intersection.push_back(
-            std::pair<Vertex_handle, Node>(
-                begin1->first,
-                Node(NULL, maximum(begin1->second.filtration(), begin2->second.filtration(), filtration))));
-        ++begin1;
-        ++begin2;
-        if (begin1 == end1 || begin2 == end2)
+        Filtration_value filt = (std::max)({begin1->second.filtration(), begin2->second.filtration(), filtration_});
+        intersection.emplace_back(begin1->first, Node(NULL, filt));
+        if (++begin1 == end1 || ++begin2 == end2)
+          return;  // ----->>
+      } else if (begin1->first < begin2->first) {
+        if (++begin1 == end1)
+          return;
+      } else /* begin1->first > begin2->first */ {
+        if (++begin2 == end2)
           return;  // ----->>
-      } else {
-        if (begin1->first < begin2->first) {
-          ++begin1;
-          if (begin1 == end1)
-            return;
-        } else {
-          ++begin2;
-          if (begin2 == end2)
-            return;  // ----->>
-        }
       }
     }
   }
-  /** Maximum over 3 values.*/
-  Filtration_value maximum(Filtration_value a, Filtration_value b,
-                           Filtration_value c) {
-    Filtration_value max = (a < b) ? b : a;
-    return ((max < c) ? c : max);
-  }
 
  public:
   /** \brief Write the hasse diagram of the simplicial complex in os.
@@ -1047,8 +1120,9 @@ std::pair<Simplex_handle, bool> insert_simplex(const RandomAccessVertexRange & s
 };
 
 // Print a Simplex_tree in os.
-template<typename T1, typename T2, typename T3>
-std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) {
+
+template<typename...T>
+std::ostream& operator<<(std::ostream & os, Simplex_tree<T...> & st) {
   for (auto sh : st.filtration_simplex_range()) {
     os << st.dimension(sh) << " ";
     for (auto v : st.simplex_vertex_range(sh)) {
@@ -1058,25 +1132,30 @@ std::ostream& operator<<(std::ostream & os, Simplex_tree<T1, T2, T3> & st) {
   }
   return os;
 }
-template<typename T1, typename T2, typename T3>
-std::istream& operator>>(std::istream & is, Simplex_tree<T1, T2, T3> & st) {
+
+template<typename...T>
+std::istream& operator>>(std::istream & is, Simplex_tree<T...> & st) {
   // assert(st.num_simplices() == 0);
 
-  std::vector<typename Simplex_tree<T1, T2, T3>::Vertex_handle> simplex;
-  typename Simplex_tree<T1, T2, T3>::Filtration_value fil;
-  typename Simplex_tree<T1, T2, T3>::Filtration_value max_fil = 0;
+  typedef Simplex_tree<T...> ST;
+  std::vector<typename ST::Vertex_handle> simplex;
+  typename ST::Filtration_value fil;
+  typename ST::Filtration_value max_fil = 0;
   int max_dim = -1;
   size_t num_simplices = 0;
-  while (read_simplex(is, simplex, fil)) {  // read all simplices in the file as a list of vertices
+  while (read_simplex(is, simplex, fil)) {
+    // read all simplices in the file as a list of vertices
     ++num_simplices;
-    int dim = static_cast<int>(simplex.size() - 1);  // Warning : simplex_size needs to be casted in int - Can be 0
+    // Warning : simplex_size needs to be casted in int - Can be 0
+    int dim = static_cast<int> (simplex.size() - 1);
     if (max_dim < dim) {
       max_dim = dim;
     }
     if (max_fil < fil) {
       max_fil = fil;
     }
-    st.insert_simplex(simplex, fil);  // insert every simplex in the simplex tree
+    // insert every simplex in the simplex tree
+    st.insert_simplex(simplex, fil);
     simplex.clear();
   }
   st.set_num_simplices(num_simplices);
@@ -1085,8 +1164,8 @@ std::istream& operator>>(std::istream & is, Simplex_tree<T1, T2, T3> & st) {
 
   return is;
 }
-
 /** @} */  // end defgroup simplex_tree
+
 }  // namespace Gudhi
 
-#endif  // SRC_SIMPLEX_TREE_INCLUDE_GUDHI_SIMPLEX_TREE_H_
+#endif  // SIMPLEX_TREE_H_
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 977fafa1..de350f2d 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
@@ -77,8 +77,8 @@ class Simplex_tree_siblings {
         parent_(parent),
         members_(boost::container::ordered_unique_range, members.begin(),
                  members.end()) {
-    for (auto map_it = members_.begin(); map_it != members_.end(); map_it++) {
-      map_it->second.assign_children(this);
+    for (auto& map_el : members_) {
+      map_el.second.assign_children(this);
     }
   }
 
@@ -90,19 +90,12 @@ class Simplex_tree_siblings {
    * present in the node.
    */
   void insert(Vertex_handle v, Filtration_value filtration_value) {
-    typename Dictionary::iterator sh = members_.find(v);
-    if (sh != members_.end() && sh->second.filtration() > filtration_value) {
-      sh->second.assign_filtration(filtration_value);
-      return;
-    }
-    if (sh == members_.end()) {
-      members_.insert(
-          std::pair<Vertex_handle, Node>(v, Node(this, filtration_value)));
-      return;
-    }
+    auto ins = members_.emplace(v, Node(this, filtration_value));
+    if (!ins.second && filtration(ins.first) > filtration_value)
+      ins.first->second.assign_filtration(filtration_value);
   }
 
-  typename Dictionary::iterator find(Vertex_handle v) {
+  Dictionary_it find(Vertex_handle v) {
     return members_.find(v);
   }
 
@@ -110,7 +103,7 @@ class Simplex_tree_siblings {
     return oncles_;
   }
 
-  Vertex_handle parent() {
+  Vertex_handle parent() const {
     return parent_;
   }
 
@@ -118,7 +111,7 @@ class Simplex_tree_siblings {
     return members_;
   }
 
-  size_t size() {
+  size_t size() const {
     return members_.size();
   }