module toplex_map added

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

Former-commit-id: 61f8f8a2d3920dbc30b86dd808a0e316383af137

10 files changed, 1060 insertions, 0 deletions

diff --git a/src/Toplex_map/doc/Intro_Toplex_map.h b/src/Toplex_map/doc/Intro_Toplex_map.h
new file mode 100644
index 00000000..da9562ec
--- /dev/null
+++ b/src/Toplex_map/doc/Intro_Toplex_map.h
@@ -0,0 +1,57 @@
+/*    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:       François Godi
+ *
+ *    Copyright (C) 2017  INRIA
+ *
+ *    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/>.
+ */
+
+#ifndef DOC_TOPLEX_MAP_H_
+#define DOC_TOPLEX_MAP_H_
+
+// needs namespace for Doxygen to link on classes
+namespace Gudhi {
+
+/**  \defgroup toplex_map Toplex Map
+ * 
+ * \author    Fran&ccedil;ois Godi
+ * @{
+ * 
+ * \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,
+ * 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$.
+ *
+ * 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
+ * containing this vertex.
+ *
+ * \image html map.png
+ *
+ */
+/** @} */  // end defgroup toplex_map
+
+}  // namespace Gudhi
+
+#endif  // DOC_TOPLEX_MAP_H_
diff --git a/src/Toplex_map/doc/map.png b/src/Toplex_map/doc/map.png
new file mode 100644
index 00000000..d1987043
--- /dev/null
+++ b/src/Toplex_map/doc/map.png
diff --git a/src/Toplex_map/example/CMakeLists.txt b/src/Toplex_map/example/CMakeLists.txt
new file mode 100644
index 00000000..2341fe06
--- /dev/null
+++ b/src/Toplex_map/example/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/example/chrono.cpp b/src/Toplex_map/example/chrono.cpp
new file mode 100644
index 00000000..d93d1e1f
--- /dev/null
+++ b/src/Toplex_map/example/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/include/gudhi/Fake_simplex_tree.h b/src/Toplex_map/include/gudhi/Fake_simplex_tree.h
new file mode 100644
index 00000000..5c7e7b12
--- /dev/null
+++ b/src/Toplex_map/include/gudhi/Fake_simplex_tree.h
@@ -0,0 +1,207 @@
+#ifndef FAKE_SIMPLEX_TREE_H
+#define FAKE_SIMPLEX_TREE_H
+
+#include <gudhi/Filtered_toplex_map.h>
+#include <boost/graph/adjacency_list.hpp>
+
+namespace Gudhi {
+
+class Fake_simplex_tree : public Filtered_toplex_map {
+
+public:
+
+    typedef Vertex Vertex_handle;
+
+    typedef Simplex_ptr Simplex_handle;
+
+    /** \brief Inserts 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. */
+    void expansion(int max_dim);
+
+    /** \brief Returns the number of vertices in the simplicial complex. */
+    std::size_t num_vertices();
+
+    Simplex_ptr_set candidates() const;
+
+    std::size_t dimension() const;
+
+    std::size_t num_simplices() const;
+
+    std::size_t num_vertices() const;
+
+    Simplex simplex_vertex_range(Simplex_ptr &sptr) const;
+
+    std::vector<Simplex_ptr> max_simplices() const;
+
+    std::unordered_set<Simplex_ptr> filtration_simplex_range() const;
+
+    std::unordered_set<Simplex_ptr> skeleton_simplex_range(int d=std::numeric_limits<int>::max()) const;
+
+    std::size_t dimension(Simplex_ptr& sptr) const;
+
+    void assign_filtration(Simplex_ptr& f_simplex, Filtration_value alpha_complex_filtration);
+
+    void make_filtration_non_decreasing();
+
+protected:
+
+    /** \internal Does all the facets of the given simplex belong to the complex ?
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    bool all_facets_inside(const Input_vertex_range &vertex_range) const;
+
+};
+
+template<class OneSkeletonGraph>
+void Fake_simplex_tree::insert_graph(const OneSkeletonGraph& skel_graph){
+    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) {
+        auto u = source(*e_it, skel_graph);
+        auto v = target(*e_it, skel_graph);
+        if(u<v){
+            Simplex s;
+            s.emplace(u);
+            s.emplace(v);
+            insert_simplex_and_subfaces(s);
+        }
+    }
+}
+
+void Fake_simplex_tree::expansion(int max_dim){
+    for(int d=3; d <= max_dim; d++){
+        auto cs = candidates();
+        if(cs.empty()) return;
+        for(const Simplex_ptr& sptr: cs){
+            Simplex sigma = *sptr;
+            insert_simplex_and_subfaces(sigma);
+        }
+    }
+}
+
+template <typename Input_vertex_range>
+bool Fake_simplex_tree::all_facets_inside(const Input_vertex_range &vertex_range) const{
+    Simplex sigma(vertex_range);
+    for(const Simplex& s : facets(sigma))
+        if(!filtrations.count(get_key(s))) return false;
+    return true;
+}
+
+Simplex_ptr_set Fake_simplex_tree::candidates() const{
+    Simplex_ptr_set c;
+    std::unordered_map<Simplex_ptr, std::vector<Vertex>, Toplex_map::Sptr_hash, Toplex_map::Sptr_equal> facets_to_max;
+    for(const auto& kv : filtrations){
+        Simplex sigma (*(kv.first));
+        for(Vertex v : sigma){
+            sigma.erase(v);
+            auto sptr = get_key(sigma);
+            if(!facets_to_max.count(sptr)) facets_to_max.emplace(sptr, std::vector<Vertex>());
+            facets_to_max.at(sptr).emplace_back(v);
+            sigma.insert(v);
+        }
+    }
+    for(const auto& kv : facets_to_max){
+        std::unordered_set<Vertex> facets(kv.second.begin(), kv.second.end());
+        for(Vertex v : kv.second){
+            facets.erase(v);
+            for(Vertex w : facets){
+                Simplex sigma(*(kv.first));
+                sigma.insert(v);
+                sigma.insert(w);
+                if(all_facets_inside(sigma))
+                    c.emplace(get_key(sigma));
+            }
+            facets.emplace(v);
+        }
+    }
+    return c;
+}
+
+std::size_t Fake_simplex_tree::dimension() const {
+    std::size_t max = 0;
+    for(auto kv : filtrations)
+        max = std::max(max, kv.first->size());
+    return max;
+}
+
+std::size_t Fake_simplex_tree::num_simplices() const {
+    return filtration_simplex_range().size();
+}
+
+std::size_t Fake_simplex_tree::num_vertices() const {
+    std::unordered_set<Vertex> vertices;
+    for(auto kv : filtrations)
+        for (Vertex v : *(kv.first))
+            vertices.emplace(v);
+    return vertices.size();
+}
+
+Simplex Fake_simplex_tree::simplex_vertex_range(Simplex_ptr& sptr) const {
+    return *sptr;
+}
+
+std::unordered_set<Simplex_ptr> Fake_simplex_tree::filtration_simplex_range() const{
+    std::vector<Simplex_ptr> m = max_simplices();
+    std::unordered_set<Simplex_ptr> seen;
+    while(m.begin()!=m.end()){
+        Simplex_ptr& sptr = m.back();
+        m.pop_back();
+        if(seen.find(sptr)!=seen.end()){
+            seen.emplace(sptr);
+            for(Simplex& sigma : facets(*sptr))
+                m.emplace_back(get_key(sigma));
+        }
+    }
+    return seen;
+}
+
+std::unordered_set<Simplex_ptr> Fake_simplex_tree::skeleton_simplex_range(int d) const{
+    std::unordered_set<Simplex_ptr> simplices;
+    for(auto sptr: filtration_simplex_range())
+        if(sptr->size()<=d)
+            simplices.emplace(sptr);
+    return simplices;
+}
+
+std::vector<Simplex_ptr> Fake_simplex_tree::max_simplices() const{
+    std::vector<Simplex_ptr> s;
+    for(auto kv : filtrations)
+        s.emplace_back(kv.first);
+    return s;
+}
+
+std::size_t Fake_simplex_tree::dimension(Simplex_ptr& sptr) const{
+    return sptr->size();
+}
+
+
+void Fake_simplex_tree::assign_filtration(Simplex_ptr& f_simplex, Filtration_value alpha_complex_filtration){
+    filtrations.emplace(f_simplex,alpha_complex_filtration);
+}
+
+void Fake_simplex_tree::make_filtration_non_decreasing(){
+    for(auto yt = filtrations.begin(); yt != filtrations.end(); ++yt)
+        for (auto it = toplex_maps.begin(); it != toplex_maps.end(); ++it){
+            if(it->first == yt -> second)
+                break;
+            if(it->second.membership(*(yt->first)))
+                for(const Simplex_ptr& sptr : it->second.maximal_cofaces(*(yt->first))){
+                    it->second.erase_maximal(sptr);
+                    toplex_maps.at(yt->second).insert_simplex(*sptr);
+                    filtrations.emplace(sptr,yt->second);
+                }
+        }
+
+}
+
+
+
+} //namespace Gudhi
+
+#endif /* FAKE_SIMPLEX_TREE_H */
+
diff --git a/src/Toplex_map/include/gudhi/Filtered_toplex_map.h b/src/Toplex_map/include/gudhi/Filtered_toplex_map.h
new file mode 100644
index 00000000..4b626f11
--- /dev/null
+++ b/src/Toplex_map/include/gudhi/Filtered_toplex_map.h
@@ -0,0 +1,45 @@
+#ifndef FILTERED_TOPLEX_MAP_H
+#define FILTERED_TOPLEX_MAP_H
+
+#include <gudhi/Toplex_map.h>
+#include <limits>
+
+#define filtration_upper_bound std::numeric_limits<Filtration_value>::max()
+
+namespace Gudhi {
+
+typedef double Filtration_value;
+
+class Filtered_toplex_map {
+
+public:
+    template <typename Input_vertex_range>
+    void insert_simplex_and_subfaces(const Input_vertex_range &vertex_range, Filtration_value f = filtration_upper_bound);
+
+    template <typename Input_vertex_range>
+    Filtration_value filtration(const Input_vertex_range &vertex_range) const;
+
+protected:
+    std::unordered_map<Filtration_value, Toplex_map> toplex_maps;
+    std::unordered_map<Simplex_ptr, Filtration_value> filtrations;
+
+};
+
+template <typename Input_vertex_range>
+void Filtered_toplex_map::insert_simplex_and_subfaces(const Input_vertex_range &vertex_range, Filtration_value f){
+    if(!toplex_maps.count(f)) toplex_maps.emplace(f,Toplex_map());
+    toplex_maps.at(f).insert_simplex(vertex_range);
+    filtrations.emplace(get_key(vertex_range),f);
+}
+
+template <typename Input_vertex_range>
+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;
+}
+
+} //namespace Gudhi
+
+#endif /* FILTERED_TOPLEX_MAP_H */
diff --git a/src/Toplex_map/include/gudhi/Lazy_toplex_map.h b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h
new file mode 100644
index 00000000..3ffe8214
--- /dev/null
+++ b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h
@@ -0,0 +1,218 @@
+#ifndef LAZY_TOPLEX_MAP_H
+#define LAZY_TOPLEX_MAP_H
+
+#include <gudhi/Toplex_map.h>
+#include <boost/heap/fibonacci_heap.hpp>
+
+namespace Gudhi {
+
+class Lazy_Toplex_map {
+
+public:
+    template <typename Input_vertex_range>
+    void insert_max_simplex(const Input_vertex_range &vertex_range);
+    template <typename Input_vertex_range>
+    bool insert_simplex(const Input_vertex_range &vertex_range);
+    template <typename Input_vertex_range>
+    void remove_simplex(const Input_vertex_range &vertex_range);
+
+    template <typename Input_vertex_range>
+    bool membership(const Input_vertex_range &vertex_range);
+    template <typename Input_vertex_range>
+    bool all_facets_inside(const Input_vertex_range &vertex_range);
+
+    Vertex contraction(const Vertex x, const Vertex y);
+
+    std::size_t num_simplices() const;
+
+private:
+    template <typename Input_vertex_range>
+    void erase_max(const Input_vertex_range &vertex_range);
+    template <typename Input_vertex_range>
+    Vertex best_index(const Input_vertex_range &vertex_range);
+    void clean(const Vertex v);
+
+    std::unordered_map<Vertex, Simplex_ptr_set> t0;
+    bool empty_toplex; // Is the empty simplex a toplex ?
+
+    typedef boost::heap::fibonacci_heap<std::pair<std::size_t,Vertex>> PriorityQueue;
+    PriorityQueue cleaning_priority;
+    std::unordered_map<Vertex, PriorityQueue::handle_type> cp_handles;
+
+    std::unordered_map<Vertex, std::size_t> gamma0_lbounds;
+    std::size_t get_gamma0_lbound(const Vertex v) const;
+
+    std::size_t size_lbound = 0;
+    std::size_t size = 0;
+
+    const double alpha = 2; //time
+    const double betta = 3; //memory
+};
+
+template <typename Input_vertex_range>
+void Lazy_Toplex_map::insert_max_simplex(const Input_vertex_range &vertex_range){
+    for(const Vertex& v : vertex_range)
+        if(!gamma0_lbounds.count(v)) gamma0_lbounds.emplace(v,1);
+        else gamma0_lbounds[v]++;
+    size_lbound++;
+    insert_simplex(vertex_range);
+}
+
+template <typename Input_vertex_range>
+bool Lazy_Toplex_map::insert_simplex(const Input_vertex_range &vertex_range){
+    Simplex sigma(vertex_range.begin(),vertex_range.end());
+    empty_toplex = (sigma.size()==0); //vérifier la gestion de empty face
+    Simplex_ptr sptr = std::make_shared<Simplex>(sigma);
+    bool inserted = false;
+    for(const Vertex& v : sigma){
+        if(!t0.count(v)){
+            t0.emplace(v, Simplex_ptr_set());
+            auto v_handle = cleaning_priority.push(std::make_pair(0, v));
+            cp_handles.emplace(v, v_handle);
+        }
+        inserted = t0.at(v).emplace(sptr).second;
+        cleaning_priority.update(cp_handles.at(v), std::make_pair(t0.at(v).size() - get_gamma0_lbound(v),v));
+    }
+    if(inserted)
+        size++;
+    if(size > size_lbound * betta)
+        clean(cleaning_priority.top().second);
+    return inserted;
+}
+
+template <typename Input_vertex_range>
+void Lazy_Toplex_map::remove_simplex(const Input_vertex_range &vertex_range){
+    if(vertex_range.begin()==vertex_range.end()){
+        t0.clear();
+        gamma0_lbounds.clear();
+        cleaning_priority.clear();
+        size_lbound = 0;
+        size = 0;
+        empty_toplex = false;
+    }
+    else {
+        const Vertex& v = best_index(vertex_range);
+        //Copy constructor needed because the set is modified
+        if(t0.count(v)) for(const Simplex_ptr& sptr : Simplex_ptr_set(t0.at(v)))
+            if(included(vertex_range, *sptr)){
+                erase_max(*sptr);
+                for(const Simplex& f : facets(vertex_range))
+                    insert_max_simplex(f);
+            }
+    }
+}
+
+template <typename Input_vertex_range>
+bool Lazy_Toplex_map::membership(const Input_vertex_range &vertex_range){
+    if(t0.size()==0 && !empty_toplex) return false; //empty complex
+    if(vertex_range.begin()==vertex_range.end()) return true; //empty query simplex
+    Vertex v = best_index(vertex_range);
+    if(!t0.count(v))  return false;
+    for(const Simplex_ptr& sptr : t0.at(v))
+        if(included(vertex_range, *sptr)) return true;
+    return false;
+}
+
+template <typename Input_vertex_range>
+bool Lazy_Toplex_map::all_facets_inside(const Input_vertex_range &vertex_range){
+    Simplex sigma(vertex_range.begin(),vertex_range.end());
+    Vertex v = best_index(sigma);
+    if(!t0.count(v))  return false;
+    Simplex f = sigma; f.erase(v);
+    if(!membership(f)) return false;
+    std::unordered_set<Vertex> facets_inside;
+    for(const Simplex_ptr& sptr : t0.at(v))
+        for(const Vertex& w : sigma){
+            f = sigma; f.erase(w);
+            if(included(f, *sptr)) facets_inside.insert(w);
+        }
+    return facets_inside.size() == sigma.size() - 1;
+}
+
+/* Returns the remaining vertex */
+Vertex Lazy_Toplex_map::contraction(const Vertex x, const Vertex y){
+    if(!t0.count(x)) return y;
+    if(!t0.count(y)) return x;
+    Vertex k, d;
+    if(t0.at(x).size() > t0.at(y).size())
+        k=x, d=y;
+    else
+        k=y, d=x;
+    //Copy constructor needed because the set is modified
+    for(const Simplex_ptr& sptr : Simplex_ptr_set(t0.at(d))){
+        Simplex sigma(*sptr);
+        erase_max(sigma);
+        sigma.erase(d);
+        sigma.insert(k);
+        insert_simplex(sigma);
+    }
+    t0.erase(d);
+    return k;
+}
+
+/* No facets insert_simplexed */
+template <typename Input_vertex_range>
+inline void Lazy_Toplex_map::erase_max(const Input_vertex_range &vertex_range){
+    Simplex sigma(vertex_range.begin(),vertex_range.end());
+    empty_toplex = false;
+    Simplex_ptr sptr = std::make_shared<Simplex>(sigma);
+    bool erased;
+    for(const Vertex& v : sigma){
+        erased = t0.at(v).erase(sptr) > 0;
+        if(t0.at(v).size()==0)
+            t0.erase(v);
+    }
+    if (erased)
+        size--;
+}
+
+template <typename Input_vertex_range>
+Vertex Lazy_Toplex_map::best_index(const Input_vertex_range &vertex_range){
+    Simplex tau(vertex_range.begin(),vertex_range.end());
+    std::size_t min = std::numeric_limits<size_t>::max(); Vertex arg_min = -1;
+    for(const Vertex& v : tau)
+        if(!t0.count(v)) return v;
+        else if(t0.at(v).size() < min)
+            min = t0.at(v).size(), arg_min = v;
+    if(min > alpha * get_gamma0_lbound(arg_min))
+        clean(arg_min);
+    return arg_min;
+}
+
+std::size_t Lazy_Toplex_map::get_gamma0_lbound(const Vertex v) const{
+    return gamma0_lbounds.count(v) ? gamma0_lbounds.at(v) : 0;
+}
+
+
+void Lazy_Toplex_map::clean(const Vertex v){
+    Toplex_map toplices;
+    std::unordered_map<int, std::vector<Simplex>> dsorted_simplices;
+    int max_dim = 0;
+    for(const Simplex_ptr& sptr : Simplex_ptr_set(t0.at(v))){
+        if(sptr->size() > max_dim){
+            for(int d = max_dim+1; d<=sptr->size(); d++)
+                dsorted_simplices.emplace(d, std::vector<Simplex>());
+            max_dim = sptr->size();
+        }
+        dsorted_simplices[sptr->size()].emplace_back(*sptr);
+        erase_max(*sptr);
+    }
+    for(int d = max_dim; d>=1; d--)
+        for(const Simplex &s : dsorted_simplices.at(d))
+            if(!toplices.membership(s))
+                toplices.insert_independent_simplex(s);
+    Simplex sv; sv.insert(v);
+    auto clean_cofaces = toplices.maximal_cofaces(sv);
+    size_lbound = size_lbound - get_gamma0_lbound(v) + clean_cofaces.size();
+    gamma0_lbounds[v] = clean_cofaces.size();
+    for(const Simplex_ptr& sptr : clean_cofaces)
+        insert_simplex(*sptr);
+}
+
+std::size_t Lazy_Toplex_map::num_simplices() const{
+    return size;
+}
+
+} //namespace Gudhi
+
+#endif /* LAZY_TOPLEX_MAP_H */
diff --git a/src/Toplex_map/include/gudhi/Toplex_map.h b/src/Toplex_map/include/gudhi/Toplex_map.h
new file mode 100644
index 00000000..0b6cad37
--- /dev/null
+++ b/src/Toplex_map/include/gudhi/Toplex_map.h
@@ -0,0 +1,307 @@
+#ifndef TOPLEX_MAP_H
+#define TOPLEX_MAP_H
+
+#include <vector>
+#include <unordered_set>
+#include <unordered_map>
+#include <memory>
+
+#define vertex_upper_bound std::numeric_limits<Vertex>::max()
+
+namespace Gudhi {
+
+/** Vertex is the type of vertices.
+ * \ingroup toplex_map   */
+typedef std::size_t Vertex;
+
+/** Simplex is the type of simplices.
+ * \ingroup toplex_map   */
+typedef std::unordered_set<Vertex> Simplex;
+
+/** A Toplex_map represents the simplicial complex.
+ * A "toplex" is a maximal simplex.
+ * \ingroup toplex_map   */
+class Toplex_map {
+
+public:
+    /** The type of the pointers to maximal simplices.
+     * \ingroup toplex_map   */
+    typedef std::shared_ptr<Simplex> Simplex_ptr;
+
+    struct Sptr_hash{ std::size_t operator()(const Simplex_ptr& s) const; };
+    struct Sptr_equal{ std::size_t operator()(const Simplex_ptr& a, const Simplex_ptr& b) const; };
+    /** The type of the sets of Simplex_ptr.
+     * \ingroup toplex_map   */
+    typedef std::unordered_set<Simplex_ptr, Sptr_hash, Sptr_equal> Simplex_ptr_set;
+
+    /** \brief Adds the given simplex to the complex.
+     * Nothing happens if the simplex has a coface in the complex.
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    void insert_simplex(const Input_vertex_range &vertex_range);
+
+    /** \brief Removes the given simplex and its cofaces from the complex.
+     * Its faces are kept inside.
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    void remove_simplex(const Input_vertex_range &vertex_range);
+
+    /** Does a simplex belong to the complex ?
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    bool membership(const Input_vertex_range &vertex_range) const;
+
+    /** Does a simplex is a toplex ?
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    bool maximality(const Input_vertex_range &vertex_range) const;
+
+    /** Gives a set of pointers to the maximal cofaces of a simplex.
+     * Gives the toplices if given the empty simplex.
+     * Gives not more than max_number maximal cofaces if max_number is strictly positive.
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    Simplex_ptr_set maximal_cofaces(const Input_vertex_range &vertex_range, const std::size_t max_number = 0) const;
+
+    /** Contracts one edge in the complex.
+     * The edge has to verify the link condition if you want to preserve topology.
+     * Returns the remaining vertex.
+     * \ingroup toplex_map   */
+    Vertex contraction(const Vertex x, const Vertex y);
+
+    /** Adds the given simplex to the complex.
+     * The simplex must not have neither maximal face nor coface in the complex.
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    void insert_independent_simplex(const Input_vertex_range &vertex_range);
+
+
+    /** \internal Removes a toplex without adding facets after.
+     * \ingroup toplex_map   */
+    void erase_maximal(const Simplex_ptr& sptr);
+
+    /** Removes a vertex from any simplex containing it.
+     * \ingroup toplex_map   */
+    void remove_vertex(const Vertex x);
+
+    /** \brief Number of maximal simplices.
+         * /!\ Not efficient !
+         * \ingroup toplex_map   */
+    std::size_t num_simplices() const;
+
+protected:
+    /** \internal Gives an index in order to look for a simplex quickly.
+     * \ingroup toplex_map   */
+    template <typename Input_vertex_range>
+    Vertex best_index(const Input_vertex_range &vertex_range) const;
+    
+    /** \internal The map from vertices to toplices
+     * \ingroup toplex_map   */
+    std::unordered_map<Vertex, Simplex_ptr_set> t0;
+
+};
+
+typedef Toplex_map::Simplex_ptr Simplex_ptr;
+typedef Toplex_map::Simplex_ptr_set Simplex_ptr_set;
+
+// Pointers are also used as key in the hash sets.
+template <typename Input_vertex_range>
+Simplex_ptr get_key(const Input_vertex_range &vertex_range);
+
+// Is the first simplex a face of the second ?
+template <typename Input_vertex_range1, typename Input_vertex_range2>
+bool included(const Input_vertex_range1 &vertex_range1, const Input_vertex_range2 &vertex_range2);
+
+// All the facets of the given simplex.
+template <typename Input_vertex_range>
+std::vector<Simplex> facets(const Input_vertex_range &vertex_range);
+
+template <typename Input_vertex_range>
+void Toplex_map::insert_simplex(const Input_vertex_range &vertex_range){
+    if(membership(vertex_range)) return;
+    bool replace_facets = true;
+    for(const Simplex& facet : facets(vertex_range))
+        if(!maximality(facet))
+        {
+            replace_facets=false;
+            break;
+        }
+    if(replace_facets)
+        for(const Simplex& facet : facets(vertex_range))
+            erase_maximal(get_key(facet));
+    else
+        for(const Vertex& v : vertex_range)
+            if(t0.count(v))  for(const Simplex_ptr& fptr : Simplex_ptr_set(t0.at(v)))
+                //Copy constructor needed because the set is modified
+                if(included(*fptr,vertex_range)) erase_maximal(fptr);
+    // We erase all the maximal faces of the simplex
+    insert_independent_simplex(vertex_range);
+}
+
+template <typename Input_vertex_range>
+void Toplex_map::remove_simplex(const Input_vertex_range &vertex_range){
+    if(vertex_range.begin()==vertex_range.end())
+        t0.clear();
+    // Removal of the empty simplex means cleaning everything
+    else {
+        const Vertex& v = best_index(vertex_range);
+        if(t0.count(v)) for(const Simplex_ptr& sptr : Simplex_ptr_set(t0.at(v)))
+            //Copy constructor needed because the set is modified
+            if(included(vertex_range, *sptr)){
+                erase_maximal(sptr);
+                for(const Simplex& f : facets(vertex_range))
+                    if(!membership(f)) insert_independent_simplex(f);
+                // We add the facets which are new maximal simplices
+            }
+    }
+}
+
+template <typename Input_vertex_range>
+bool Toplex_map::membership(const Input_vertex_range &vertex_range) const{
+    if(t0.size()==0) return false;
+    const Vertex& v = best_index(vertex_range);
+    if(!t0.count(v))  return false;
+    if(maximality(vertex_range)) return true;
+    for(const Simplex_ptr& sptr : t0.at(v))
+        if(included(vertex_range, *sptr))
+            return true;
+    return false;
+}
+
+template <typename Input_vertex_range>
+bool Toplex_map::maximality(const Input_vertex_range &vertex_range) const{
+    const Vertex& v =  best_index(vertex_range);
+    if(!t0.count(v)) return false;
+    return t0.at(v).count(get_key(vertex_range));
+}
+
+template <typename Input_vertex_range>
+Simplex_ptr_set Toplex_map::maximal_cofaces(const Input_vertex_range &vertex_range, const std::size_t max_number) const{
+    Simplex_ptr_set cofaces;
+    if(maximality(vertex_range))
+        cofaces.emplace(get_key(vertex_range));
+    else if(vertex_range.begin()==vertex_range.end())
+        for(const auto& kv : t0)
+            for(const Simplex_ptr& sptr : kv.second){
+                //kv.second is a Simplex_ptr_set
+                cofaces.emplace(sptr);
+                if(cofaces.size()==max_number)
+                    return cofaces;
+            }
+    else {
+        const Vertex& v = best_index(vertex_range);
+        if(t0.count(v)) for(const Simplex_ptr& sptr : t0.at(v))
+            if(included(vertex_range, *sptr)){
+                cofaces.emplace(sptr);
+                if(cofaces.size()==max_number)
+                    return cofaces;
+            }
+    }
+    return cofaces;
+}
+
+Vertex Toplex_map::contraction(const Vertex x, const Vertex y){
+    if(!t0.count(x)) return y;
+    if(!t0.count(y)) return x;
+    int k, d;
+    if(t0.at(x).size() > t0.at(y).size())
+        k=x, d=y;
+    else
+        k=y, d=x;
+    for(const Simplex_ptr& sptr : Simplex_ptr_set(t0.at(d))){
+        //Copy constructor needed because the set is modified
+        Simplex sigma(*sptr);
+        erase_maximal(sptr);
+        sigma.erase(d);
+        sigma.insert(k);
+        insert_simplex(sigma);
+    }
+    return k;
+}
+
+template <typename Input_vertex_range>
+void Toplex_map::insert_independent_simplex(const Input_vertex_range &vertex_range){
+    for(const Vertex& v : vertex_range){
+        if(!t0.count(v)) t0.emplace(v, Simplex_ptr_set());
+        t0.at(v).emplace(get_key(vertex_range));
+    }
+}
+
+void Toplex_map::remove_vertex(const Vertex x){
+    for(const Simplex_ptr& sptr : Simplex_ptr_set(t0.at(x))){
+        Simplex sigma(*sptr);
+        erase_maximal(sptr);
+        sigma.erase(x);
+        insert_simplex(sigma);
+    }
+}
+
+std::size_t Toplex_map::num_simplices() const{
+    return maximal_cofaces(Simplex()).size();
+}
+
+inline void Toplex_map::erase_maximal(const Simplex_ptr& sptr){
+    Simplex sigma(*sptr);
+    if (sptr->size()==0)
+        sigma.insert(vertex_upper_bound);
+    for(const Vertex& v : sigma){
+        t0.at(v).erase(sptr);
+        if(t0.at(v).size()==0) t0.erase(v);
+    }
+}
+
+template <typename Input_vertex_range>
+Vertex Toplex_map::best_index(const Input_vertex_range &vertex_range) const{
+    std::size_t min = std::numeric_limits<size_t>::max();
+    Vertex arg_min = vertex_upper_bound;
+    for(const Vertex& v : vertex_range)
+        if(!t0.count(v)) return v;
+        else if(t0.at(v).size() < min)
+            min = t0.at(v).size(), arg_min = v;
+    return arg_min;
+}
+
+std::size_t Toplex_map::Sptr_equal::operator()(const Simplex_ptr& s1, const Simplex_ptr& s2) const {
+    if (s1->size() != s2->size()) return false;
+    return included(*s1,*s2);
+    // inclusion tests equality for same size simplices
+}
+
+std::size_t Toplex_map::Sptr_hash::operator()(const Simplex_ptr& s) const {
+    std::hash<double> h_f;
+    //double hash works better than int hash
+    size_t h = 0;
+    for(const Vertex& v : *s)
+        h += h_f(static_cast<double>(v));
+    return h;
+}
+
+template <typename Input_vertex_range>
+Simplex_ptr get_key(const Input_vertex_range &vertex_range){
+    Simplex s(vertex_range.begin(), vertex_range.end());
+    return std::make_shared<Simplex>(s);
+}
+
+template <typename Input_vertex_range1, typename Input_vertex_range2>
+bool included(const Input_vertex_range1 &vertex_range1, const Input_vertex_range2 &vertex_range2){
+    Simplex s2(vertex_range2.begin(), vertex_range2.end());
+    for(const Vertex& v : vertex_range1)
+        if(!s2.count(v)) return false;
+    return true;
+}
+
+template <typename Input_vertex_range>
+std::vector<Simplex> facets(const Input_vertex_range &vertex_range){
+    std::vector<Simplex> facets;
+    Simplex f(vertex_range.begin(), vertex_range.end());
+    for(const Vertex& v : vertex_range){
+        f.erase(v);
+        facets.emplace_back(f);
+        f.insert(v);
+    }
+    return facets;
+}
+
+} //namespace Gudhi
+
+#endif /* TOPLEX_MAP_H */
diff --git a/src/Toplex_map/test/CMakeLists.txt b/src/Toplex_map/test/CMakeLists.txt
new file mode 100644
index 00000000..223ebccb
--- /dev/null
+++ b/src/Toplex_map/test/CMakeLists.txt
@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 2.6)
+project(Toplex_map_tests)
+
+add_executable ( ToplexMapUT test.cpp )
+target_link_libraries(ToplexMapUT ${Boost_SYSTEM_LIBRARY} ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+
+
+# Unitary tests
+add_test(NAME SalUT 
+  COMMAND ${CMAKE_CURRENT_BINARY_DIR}/ToplexMapUT
+    ${CMAKE_SOURCE_DIR}/src/Toplex_map/test/test.txt
+    # XML format for Jenkins xUnit plugin 
+    --log_format=XML --log_sink=${CMAKE_SOURCE_DIR}/ToplexMapUT.xml --log_level=test_suite --report_level=no)
+
diff --git a/src/Toplex_map/test/test.cpp b/src/Toplex_map/test/test.cpp
new file mode 100644
index 00000000..3f4d96c2
--- /dev/null
+++ b/src/Toplex_map/test/test.cpp
@@ -0,0 +1,71 @@
+#include <iostream>
+#include <gudhi/Lazy_toplex_map.h>
+#include <gudhi/Filtered_toplex_map.h>
+#include <gudhi/Fake_simplex_tree.h>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "toplex map"
+#include <boost/test/unit_test.hpp>
+
+using namespace Gudhi;
+
+std::vector<Vertex> sigma1 = {1, 2, 3, 4};
+std::vector<Vertex> sigma2 = {5, 2, 3, 6};
+std::vector<Vertex> sigma3 = {5};
+std::vector<Vertex> sigma4 = {5, 2, 3};
+std::vector<Vertex> sigma5 = {5, 2, 7};
+std::vector<Vertex> sigma6 = {4, 5, 3};
+std::vector<Vertex> sigma7 = {4, 5, 9};
+std::vector<Vertex> sigma8 = {1, 2, 3, 6};
+
+
+BOOST_AUTO_TEST_CASE(toplexmap) {
+    Toplex_map K;
+    K.insert_simplex(sigma1);
+    K.insert_simplex(sigma2);
+    K.insert_simplex(sigma3);
+    K.insert_simplex(sigma6);
+    K.insert_simplex(sigma7);
+    BOOST_CHECK(K.membership(sigma4));
+    BOOST_CHECK(!K.maximality(sigma5));
+    BOOST_CHECK(!K.membership(sigma5));
+    K.contraction(4,5);
+    BOOST_CHECK(!K.membership(sigma6));
+}
+
+BOOST_AUTO_TEST_CASE(ltoplexmap) {
+    Lazy_Toplex_map K;
+    K.insert_simplex(sigma1);
+    K.insert_simplex(sigma2);
+    K.insert_simplex(sigma3);
+    K.insert_simplex(sigma6);
+    K.insert_simplex(sigma7);
+    BOOST_CHECK(K.membership(sigma4));
+    BOOST_CHECK(!K.membership(sigma5));
+    K.contraction(4,5);
+    BOOST_CHECK(!K.membership(sigma6));
+}
+
+BOOST_AUTO_TEST_CASE(ftoplexmap) {
+    Filtered_toplex_map K;
+    K.insert_simplex_and_subfaces(sigma1, 2.);
+    K.insert_simplex_and_subfaces(sigma2, 2.);
+    K.insert_simplex_and_subfaces(sigma6, 1.);
+    K.insert_simplex_and_subfaces(sigma7, 1.);
+    BOOST_CHECK(K.filtration(sigma4)==2.);
+    BOOST_CHECK(K.filtration(sigma3)==1.);
+}
+
+/*
+BOOST_AUTO_TEST_CASE(toplexmap_candidates) {
+    Toplex_map K;
+    K.insert_simplex(sigma1);
+    K.insert_simplex(sigma2);
+    K.remove_simplex(sigma1);
+    K.remove_simplex(sigma2);
+    auto c = K.candidates();
+    BOOST_CHECK(c.count(get_key(sigma1)));
+    BOOST_CHECK(c.count(get_key(sigma2)));
+    BOOST_CHECK(c.size()==2);
+}
+*/