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

4 files changed, 777 insertions, 0 deletions

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 */