Merge simplex tree copy/move constructor branch back into the trunk

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

Former-commit-id: ceacd34f54b14a2d2c4e0a68b234a91912b215b8

2 files changed, 225 insertions, 85 deletions

diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h
index 03e15a7d..b911552d 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>
@@ -266,7 +269,8 @@ 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). If the complex has changed since the last time the filtration was initialized, please call `initialize_filtration()` to recompute it. */
+   * 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 filtration_simplex_range(Indexing_tag) {
     if (filtration_vect_.empty()) {
       initialize_filtration();
@@ -318,10 +322,49 @@ 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_(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. */
   ~Simplex_tree() {
     for (auto sh = root_.members().begin(); sh != root_.members().end(); ++sh) {
@@ -343,6 +386,40 @@ 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.
@@ -384,7 +461,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
@@ -431,7 +508,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();
     }
@@ -449,26 +526,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();
@@ -488,6 +573,37 @@ class Simplex_tree {
   }
   //{ 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;
+    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();
+    }
+    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
+        res_insert.first->second.assign_filtration(filtration);
+        return res_insert;
+      }
+      // if filtration value unchanged
+      return std::pair<Simplex_handle, bool>(null_simplex(), false);
+    }
+    // otherwise the insertion has succeeded
+    return res_insert;
+  }
+
+ 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
@@ -497,7 +613,7 @@ class Simplex_tree {
    * 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.
    *
@@ -509,72 +625,107 @@ class Simplex_tree {
    * 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,
+   * 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) {
-    if (simplex.empty()) {
-      return std::pair<Simplex_handle, bool>(null_simplex(), true);
-    }
-    // must be sorted in increasing order
-    sort(simplex.begin(), simplex.end());
+    auto first = std::begin(simplex);
+    auto last = std::end(simplex);
 
-    Siblings * curr_sib = &root_;
-    std::pair<Simplex_handle, bool> res_insert;
-    typename RandomAccessVertexRange::iterator vi;
-    for (vi = simplex.begin(); 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();
-    }
-    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
-        res_insert.first->second.assign_filtration(filtration);
-        return res_insert;
-      }
-      // if filtration value unchanged
-      return std::pair<Simplex_handle, bool>(null_simplex(), false);
-    }
-    // otherwise the insertion has succeeded
-    return res_insert;
+    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 RandomAccessVertexRange>
-  void insert_simplex_and_subfaces(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);
+  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(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(Nsimplex, filtration);
-    } else if (Nsimplex.size() == 1) {
-      // 1-Simplex insert - End of recursivity
-      std::pair<Simplex_handle, bool> returned = insert_simplex(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) {
@@ -598,17 +749,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() {
@@ -920,7 +1060,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) {
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);
-    }
+	}
   }
 
   /*