...

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

Former-commit-id: d872eef9342046002b81f55cc12760e799df0c1d

5 files changed, 182 insertions, 47 deletions

diff --git a/src/Toplex_map/benchmarks/CMakeLists.txt b/src/Toplex_map/benchmarks/CMakeLists.txt
new file mode 100644
index 00000000..2341fe06
--- /dev/null
+++ b/src/Toplex_map/benchmarks/CMakeLists.txt
@@ -0,0 +1,4 @@
+cmake_minimum_required(VERSION 2.6)
+project(Toplex_map_examples)
+
+add_executable(chrono chrono.cpp)
diff --git a/src/Toplex_map/benchmarks/chrono.cpp b/src/Toplex_map/benchmarks/chrono.cpp
new file mode 100644
index 00000000..d93d1e1f
--- /dev/null
+++ b/src/Toplex_map/benchmarks/chrono.cpp
@@ -0,0 +1,137 @@
+#include <iostream>
+#include <random>
+#include <chrono>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Lazy_toplex_map.h>
+
+using namespace Gudhi;
+
+typedef Simplex typeVectorVertex;
+typedef std::pair< Simplex_tree<>::Simplex_handle, bool > typePairSimplexBool;
+
+class ST_wrapper {
+
+public:
+    void insert_simplex(const Simplex& tau);
+    bool membership(const Simplex& tau);
+    Vertex contraction(const Vertex x, const Vertex y);
+    std::size_t num_simplices();
+
+private:
+    Simplex_tree<> simplexTree;
+    void erase_max(const Simplex& sigma);
+};
+
+void ST_wrapper::insert_simplex(const Simplex& tau){
+    simplexTree.insert_simplex_and_subfaces(tau);
+}
+
+bool ST_wrapper::membership(const Simplex& tau) {
+    return simplexTree.find(tau) != simplexTree.null_simplex();
+}
+
+void ST_wrapper::erase_max(const Simplex& sigma){
+    if(membership(sigma))
+        simplexTree.remove_maximal_simplex(simplexTree.find(sigma));
+}
+
+Vertex ST_wrapper::contraction(const Vertex x, const Vertex y){
+    Simplex sx; sx.insert(x);
+    auto hx = simplexTree.find(sx);
+    if(hx != simplexTree.null_simplex())
+        for(auto h : simplexTree.cofaces_simplex_range(hx,0)){
+            auto sr = simplexTree.simplex_vertex_range(h);
+            Simplex sigma(sr.begin(),sr.end());
+            erase_max(sigma);
+            sigma.erase(x);
+            sigma.insert(y);
+            insert_simplex(sigma);
+        }
+    return y;
+}
+
+std::size_t ST_wrapper::num_simplices(){
+    return simplexTree.num_simplices();
+}
+
+
+
+int n = 300;
+
+int nb_insert_simplex1 = 3000;
+int nb_membership1 = 4000;
+int nb_contraction = 300;
+int nb_insert_simplex2 = 3000;
+int nb_membership2 = 400000;
+
+Simplex random_simplex(int n, int d){
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dis(1, n);
+    Simplex s;
+    while(s.size()!=d)
+        s.insert(dis(gen));
+    return s;
+}
+
+std::vector<Simplex> r_vector_simplices(int n, int max_d, int m){
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dis(1, max_d);
+    std::vector<Simplex> v;
+    for(int i=0; i<m; i++)
+        v.push_back(random_simplex(n,dis(gen)));
+    return v;
+}
+
+template<typename complex_type>
+void chrono(int n, int d){
+    complex_type K;
+    std::vector<Simplex> simplices_insert_simplex1 = r_vector_simplices(n,d,nb_insert_simplex1);
+    std::vector<Simplex> simplices_membership1 = r_vector_simplices(n,d,nb_membership1);
+    std::vector<Simplex> simplices_insert_simplex2 = r_vector_simplices(n - 2*nb_contraction,d,nb_insert_simplex2);
+    std::vector<Simplex> simplices_membership2 = r_vector_simplices(n - 2*nb_contraction,d,nb_membership2);
+    std::chrono::time_point<std::chrono::system_clock> start, end;
+
+    for(const Simplex& s : simplices_insert_simplex1)
+        K.insert_simplex(s);
+
+    for(const Simplex& s : simplices_membership1)
+        K.membership(s);
+
+    start = std::chrono::system_clock::now();
+    for(int i = 0; i<=nb_contraction; i++)
+        K.contraction(n-2*i,n-2*i-1);
+    end = std::chrono::system_clock::now();
+    auto c3 = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count();
+
+    start = std::chrono::system_clock::now();
+    for(const Simplex& s : simplices_insert_simplex2)
+        K.insert_simplex(s);
+    end = std::chrono::system_clock::now();
+    auto c1 = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count();
+
+    start = std::chrono::system_clock::now();
+    for(const Simplex& s : simplices_membership2)
+        K.membership(s);
+    end = std::chrono::system_clock::now();
+    auto c2 = std::chrono::duration_cast<std::chrono::milliseconds>(end-start).count();
+
+    std::cout << c1 << "\t \t" << c2 << "\t \t" << c3 << "\t \t" << K.num_simplices() << std::endl;
+}
+
+int main(){
+    for(int d=5;d<=40;d+=5){
+        std::cout << "d=" << d << " \t  Insertions \t   Membership \t Contractions \t        Size" << std::endl;
+        std::cout << "T Map \t \t";
+        chrono<Toplex_map>(n,d);
+        std::cout << "Lazy \t \t";
+        chrono<Lazy_Toplex_map>(n,d);
+        if(d<=15){
+            std::cout << "ST \t \t";
+            chrono<ST_wrapper>(n,d);
+        }
+        std::cout << std::endl;
+    }
+}
diff --git a/src/Toplex_map/doc/Intro_Toplex_map.h b/src/Toplex_map/doc/Intro_Toplex_map.h
index da9562ec..6f4c1a1b 100644
--- a/src/Toplex_map/doc/Intro_Toplex_map.h
+++ b/src/Toplex_map/doc/Intro_Toplex_map.h
@@ -33,15 +33,15 @@ namespace Gudhi {
  * 
  * \section toplexmapdefinition Definition
  * 
- * Let's consider a simplicial complex, denote by $d$ its dimension
- * and by $k$ its number of maximal simplices.
- * Furthermore, denote by $\gamma_0$ the maximal number of toplices, i.e. maximal simplices,
+ * Let's consider a simplicial complex, denote by \f$d\f$ its dimension
+ * and by \f$k\f$ its number of maximal simplices.
+ * Furthermore, denote by \f$\gamma_0\f$ the maximal number of toplices, i.e. maximal simplices,
  * that contain a same vertex.
  *
  * The goal of the Toplex Map is both to represent the complex in optimal
  * O(kd) space and to provide fast standard operations such as : insertion, removal
  * and membership of a simplex, contraction of an edge, collapses. The time needed
- * for these operation is linear or quadratic in $\gamma_0$ and $d$.
+ * for these operation is linear or quadratic in \f$\gamma_0\f$ and \f$d\f$.
  *
  * Toplex map is composed firstly of a raw storage of toplices and secondly of a
  * map which associate any vertex to a set of pointers toward all toplices
diff --git a/src/Toplex_map/include/gudhi/Fake_simplex_tree.h b/src/Toplex_map/include/gudhi/Fake_simplex_tree.h
index 6a0782ea..3de962af 100644
--- a/src/Toplex_map/include/gudhi/Fake_simplex_tree.h
+++ b/src/Toplex_map/include/gudhi/Fake_simplex_tree.h
@@ -1,14 +1,14 @@
 #ifndef FAKE_SIMPLEX_TREE_H
 #define FAKE_SIMPLEX_TREE_H
 
+#include <cmath>
+
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Filtered_toplex_map.h>
 
 #include <boost/graph/adjacency_list.hpp>
 #include <boost/graph/bron_kerbosch_all_cliques.hpp>
 
-#define filtration_upper_bound std::numeric_limits<Filtration_value>::max()
-
 namespace Gudhi {
 
 struct Visitor {
@@ -35,33 +35,31 @@ public:
 
     typedef void Insertion_result_type;
 
-    /** \brief Inserts a given range `Gudhi::rips_complex::Rips_complex::OneSkeletonGraph` in the simplicial
+    /** \brief Inserts the flag complex of a given range `Gudhi::rips_complex::Rips_complex::OneSkeletonGraph` in the simplicial
      * complex. */
     template<class OneSkeletonGraph>
     void insert_graph(const OneSkeletonGraph& skel_graph);
 
-    /** \brief Expands the simplicial complex containing only its one skeleton until a given maximal dimension as
-     * explained in \ref ripsdefinition. */
+    /** \brief Do nothing */
     void expansion(int max_dim);
 
-    /** \brief Returns the number of vertices in the simplicial complex. */
+    /** \brief Returns the number of vertices stored i.e. the number of max simplices */
     std::size_t num_vertices() const;
 
-    Simplex_ptr_set candidates(int min_dim=-1) const;
-
     std::size_t dimension() const;
 
+    std::size_t dimension(Simplex_ptr& sptr) const;
+
     std::size_t num_simplices() const;
 
     Simplex simplex_vertex_range(const Simplex& s) const;
 
     std::vector<Simplex> max_simplices() const;
 
-    std::vector<Simplex> filtration_simplex_range() const;
+    std::vector<Simplex> filtration_simplex_range(int d=std::numeric_limits<int>::max()) const;
 
-    std::vector<Simplex> skeleton_simplex_range(int d=std::numeric_limits<int>::max()) const;
+    std::vector<Simplex> skeleton_simplex_range(int d) const;
 
-    std::size_t dimension(Simplex_ptr& sptr) const;
 
 protected:
 
@@ -74,8 +72,8 @@ protected:
 
 template<class OneSkeletonGraph>
 void Fake_simplex_tree::insert_graph(const OneSkeletonGraph& skel_graph){
-    toplex_maps.emplace(filtration_upper_bound,Toplex_map());
-    bron_kerbosch_all_cliques(skel_graph, Visitor(&(this->toplex_maps.at(filtration_upper_bound))));
+    toplex_maps.emplace(nan(""),Toplex_map());
+    bron_kerbosch_all_cliques(skel_graph, Visitor(&(this->toplex_maps.at(nan("")))));
 }
 
 void Fake_simplex_tree::expansion(int max_dim){}
@@ -95,6 +93,10 @@ std::size_t Fake_simplex_tree::dimension() const {
     return max-1;
 }
 
+std::size_t Fake_simplex_tree::dimension(Simplex_ptr& sptr) const{
+    return sptr->size();
+}
+
 std::size_t Fake_simplex_tree::num_simplices() const {
     //return filtration_simplex_range().size();
     return max_simplices().size();
@@ -112,42 +114,35 @@ Simplex Fake_simplex_tree::simplex_vertex_range(const Simplex& s) const {
     return s;
 }
 
-std::vector<Simplex> Fake_simplex_tree::filtration_simplex_range() const{
+std::vector<Simplex> Fake_simplex_tree::max_simplices() const{
+    std::vector<Simplex> max_s;
+    for(auto kv : toplex_maps)
+        for(const Simplex_ptr& sptr : kv.second.maximal_cofaces(Simplex()))
+            max_s.emplace_back(*sptr);
+    return max_s;
+}
+
+std::vector<Simplex> Fake_simplex_tree::filtration_simplex_range(int d) const{
     std::vector<Simplex> m = max_simplices();
-    std::vector<Simplex> seen1;
-    Simplex_ptr_set seen2;
+    std::vector<Simplex> range;
+    Simplex_ptr_set seen;
     while(m.begin()!=m.end()){
         Simplex s(m.back());
         m.pop_back();
-        if(seen2.find(get_key(s))==seen2.end()){
-            seen1.emplace_back(s);
-            seen2.emplace(get_key(s));
+        if(seen.find(get_key(s))==seen.end()){
+            if(s.size()-1<=d)
+                range.emplace_back(s);
+            seen.emplace(get_key(s));
             if(s.size()>0)
                 for(Simplex& sigma : facets(s))
                     m.emplace_back(sigma);
         }
     }
-    return seen1;
+    return range;
 }
 
 std::vector<Simplex> Fake_simplex_tree::skeleton_simplex_range(int d) const{
-    std::vector<Simplex> simplices;
-    for(const Simplex& s : max_simplices())
-        if(s.size()<=d)
-            simplices.emplace_back(s);
-    return simplices;
-}
-
-std::vector<Simplex> Fake_simplex_tree::max_simplices() const{
-    std::vector<Simplex> max_s;
-    for(auto kv : toplex_maps)
-        for(const Simplex_ptr& sptr : kv.second.maximal_cofaces(Simplex()))
-            max_s.emplace_back(*sptr);
-    return max_s;
-}
-
-std::size_t Fake_simplex_tree::dimension(Simplex_ptr& sptr) const{
-    return sptr->size();
+    return filtration_simplex_range(d);
 }
 
 } //namespace Gudhi
diff --git a/src/Toplex_map/include/gudhi/Filtered_toplex_map.h b/src/Toplex_map/include/gudhi/Filtered_toplex_map.h
index 3a0064dc..a0c24304 100644
--- a/src/Toplex_map/include/gudhi/Filtered_toplex_map.h
+++ b/src/Toplex_map/include/gudhi/Filtered_toplex_map.h
@@ -2,10 +2,9 @@
 #define FILTERED_TOPLEX_MAP_H
 
 #include <gudhi/Toplex_map.h>
+#include <map>
 #include <limits>
 
-#define filtration_upper_bound std::numeric_limits<Filtration_value>::max()
-
 namespace Gudhi {
 
 class Filtered_toplex_map {
@@ -14,7 +13,7 @@ public:
     typedef double Filtration_value;
 
     template <typename Input_vertex_range>
-    std::pair<Simplex, bool> insert_simplex_and_subfaces(const Input_vertex_range &vertex_range, Filtration_value f = filtration_upper_bound);
+    std::pair<Simplex, bool> insert_simplex_and_subfaces(const Input_vertex_range &vertex_range, Filtration_value f = nan(""));
 
     template <typename Input_vertex_range>
     Filtration_value filtration(const Input_vertex_range &vertex_range) const;
@@ -23,7 +22,7 @@ public:
     bool membership(const Input_vertex_range &vertex_range) const;
 
 protected:
-    std::unordered_map<Filtration_value, Toplex_map> toplex_maps;
+    std::map<Filtration_value, Toplex_map> toplex_maps;
 };
 
 template <typename Input_vertex_range>
@@ -40,8 +39,8 @@ template <typename Input_vertex_range>
 Filtered_toplex_map::Filtration_value Filtered_toplex_map::filtration(const Input_vertex_range &vertex_range) const{
     for(auto kv : toplex_maps)
         if(kv.second.membership(vertex_range))
-            return kv.first;
-    return filtration_upper_bound;
+            return kv.first; //min only because a map is ordered
+    return nan("");
 }
 
 template <typename Input_vertex_range>