Skeleton_blocker_complex.h

/*    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(s):       David Salinas
 *
 *    Copyright (C) 2014  INRIA Sophia Antipolis-Mediterranee (France)
 *
 *    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 GUDHI_SKELETON_BLOCKER_COMPLEX_H
#define GUDHI_SKELETON_BLOCKER_COMPLEX_H


#include <map>
#include <list>
#include <set>
#include <vector>
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
#include <memory>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/connected_components.hpp>
#include <boost/iterator/transform_iterator.hpp>
#include <boost/range/adaptor/map.hpp>

#include "gudhi/Skeleton_blocker/iterators/Skeleton_blockers_iterators.h"
#include "gudhi/Skeleton_blocker_link_complex.h"
#include "gudhi/Skeleton_blocker/Skeleton_blocker_link_superior.h"
#include "gudhi/Skeleton_blocker/Skeleton_blocker_sub_complex.h"
#include "gudhi/Skeleton_blocker/Skeleton_blocker_simplex.h"

#include "gudhi/Skeleton_blocker/Skeleton_blocker_complex_visitor.h"
#include "gudhi/Skeleton_blocker/internal/Top_faces.h"
#include "gudhi/Utils.h"


namespace Gudhi {

namespace skbl {


template<class SkeletonBlockerDS>
class Skeleton_blocker_complex
{
    template<class ComplexType> friend class Complex_vertex_iterator;
    template<class ComplexType> friend class Complex_neighbors_vertices_iterator;
    template<class ComplexType> friend class Complex_edge_iterator;
    template<class ComplexType> friend class Complex_edge_around_vertex_iterator;


    template<class ComplexType> friend class Skeleton_blocker_link_complex;
    template<class ComplexType> friend class Skeleton_blocker_link_superior;
    template<class ComplexType> friend class Skeleton_blocker_sub_complex;

public:

    typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;

    typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;

    typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;

    typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
    typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;


    typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
    typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;


    typedef Simplex_handle* Blocker_handle;



    typedef typename Root_simplex_handle::Simplex_vertex_const_iterator Root_simplex_iterator;
    typedef typename Simplex_handle::Simplex_vertex_const_iterator Simplex_handle_iterator;


protected:

    typedef typename boost::adjacency_list
            < boost::setS, //edges
            boost::vecS, // vertices
            boost::undirectedS,
            Graph_vertex,
            Graph_edge
            > Graph;
    //todo/remark : edges are not sorted, it heavily penalizes computation for SuperiorLink
    // (eg Link with greater vertices)
    // that burdens simplex iteration / complex initialization via list of simplices.
    // to avoid that, one should modify the graph by storing two lists of adjacency for every
    // vertex, the one with superior and the one with lower vertices, that way there is
    // no more extra cost for computation of SuperiorLink
    typedef typename boost::graph_traits<Graph>::vertex_iterator boost_vertex_iterator;
    typedef typename boost::graph_traits<Graph>::edge_iterator boost_edge_iterator;

protected:
    typedef typename boost::graph_traits<Graph>::adjacency_iterator boost_adjacency_iterator;

public:
    typedef typename boost::graph_traits<Graph>::edge_descriptor Edge_handle;


protected:
    typedef std::multimap<Vertex_handle,Simplex_handle *> BlockerMap;
    typedef typename std::multimap<Vertex_handle,Simplex_handle *>::value_type BlockerPair;
    typedef typename std::multimap<Vertex_handle,Simplex_handle *>::iterator BlockerMapIterator;
    typedef typename std::multimap<Vertex_handle,Simplex_handle *>::const_iterator BlockerMapConstIterator;

protected:
    int num_vertices_;
    int num_blockers_;

    typedef Skeleton_blocker_complex_visitor<Vertex_handle> Visitor;
    //  typedef Visitor* Visitor_ptr;
    Visitor* visitor;

    std::vector<boost_vertex_handle> degree_;
    Graph skeleton; 
    //todo remove!!!

    BlockerMap blocker_map_;




public:





    Skeleton_blocker_complex(int num_vertices_ = 0,Visitor* visitor_=NULL):visitor(visitor_){
        clear();
        for (int i=0; i<num_vertices_; ++i){
            add_vertex();
        }
    }

private:

    class Simplices_sets_from_list{
    private:
        typedef  std::set< Simplex_handle> Container_simplices;
        typedef typename Container_simplices::iterator Simplices_iterator;
        int dimension_;
        std::vector<Container_simplices > simplices_;

    public:
        Simplices_sets_from_list(std::list<Simplex_handle>& simplices):
            dimension_(-1){
            assert(!simplices.empty());

            for(auto simplex = simplices.begin() ; simplex != simplices.end(); ++simplex ){
                dimension_ = std::max(dimension_,(int)simplex->dimension());
            }
            simplices_ = std::vector<Container_simplices >(dimension_+1);

            // compute k-simplices
            for(auto simplex = simplices.begin() ; simplex != simplices.end(); ++simplex ){
                simplices_[simplex->dimension()].insert(*simplex);
            }
        }

        Simplices_iterator begin(int k){
            assert(0<= k && k<= dimension_);
            return simplices_[k].begin();
        }

        Simplices_iterator end(int k){
            assert(0<= k && k<= dimension_);
            return simplices_[k].end();
        }


        Container_simplices& simplices(int k){
            return simplices_[k];
        }

        int dimension(){
            return dimension_;
        }

        bool contains(const Simplex_handle& simplex) const{
            if(simplex.dimension()>dimension_)
                return false;
            else
                return simplices_[simplex.dimension()].find(simplex)!= simplices_[simplex.dimension()].end();
        }

        void print(){
            for(int i = 0; i < dimension_; ++i){
                std::cout << i<<"-simplices"<<std::endl;
                auto l = simplices_[i];
                for(auto s : l){
                    std::cout << s<<std::endl;
                }
            }
        }
    };


    void compute_next_expand(
            Simplices_sets_from_list& simplices,
            int dim,
            std::list<Simplex_handle>& next_expand)
    {
        next_expand.clear();

        //      high_resolution_clock::time_point tbeginfor = high_resolution_clock::now();
        //      auto durationlooplink = std::chrono::duration_cast<std::chrono::microseconds>( tbeginfor - tbeginfor ).count();

        for(auto sigma = simplices.begin(dim); sigma != simplices.end(dim); ++sigma){
            //          high_resolution_clock::time_point tbeg = high_resolution_clock::now();
            Simplex_handle t(*sigma);
            Skeleton_blocker_link_superior<Skeleton_blocker_complex> link(*this,t);
            // xxx all time here, most likely because accessing superior edges needs passing through lower one
            // currently

            //          durationlooplink += std::chrono::duration_cast<std::chrono::microseconds>( high_resolution_clock::now() - tbeg ).count();

            for(auto v : link.vertex_range()){
                Vertex_handle v_in_complex(*this->get_address(  link.get_id(v)) );
                t.add_vertex(v_in_complex);
                next_expand.push_back(t);
                t.remove_vertex(v_in_complex);
            }
        }
        //      high_resolution_clock::time_point t2 = high_resolution_clock::now();
        //      auto durationlooptotal = std::chrono::duration_cast<std::chrono::microseconds>( t2 - tbeginfor ).count();
        //      DBGVALUE(durationlooptotal);
        //      DBGVALUE(durationlooplink);
    }



public:
    Skeleton_blocker_complex(std::list<Simplex_handle>& simplices,Visitor* visitor_=NULL):
        num_vertices_(0),num_blockers_(0),
        visitor(visitor_){
        Simplices_sets_from_list set_simplices(simplices);

        int dim = set_simplices.dimension();


        // add 1-skeleton to the complex
        for(auto v_it = set_simplices.begin(0); v_it != set_simplices.end(0); ++v_it)
            add_vertex();

        for(auto e_it = set_simplices.begin(1); e_it != set_simplices.end(1); ++e_it){
            Vertex_handle a = e_it->first_vertex();
            Vertex_handle b = e_it->last_vertex();
            assert(contains_vertex(a) && contains_vertex(b));
            add_edge(a,b);
        }

        // then add blockers
        for(int current_dim = 1 ; current_dim <=dim ; ++current_dim){
            std::list<Simplex_handle> expansion_simplices;
            compute_next_expand(set_simplices,current_dim,expansion_simplices);

            for(const auto &simplex : expansion_simplices) {
                if(!set_simplices.contains(simplex)){
                    add_blocker(simplex);
                }
            }
        }
    }


    // We cannot use the default copy constructor since we need
    // to make a copy of each of the blockers
    Skeleton_blocker_complex(const Skeleton_blocker_complex& copy){
        visitor = NULL;
        degree_ = copy.degree_;
        skeleton = Graph(copy.skeleton);
        num_vertices_ = copy.num_vertices_;

        num_blockers_ = 0;
        // we copy the blockers
        for (auto blocker : copy.const_blocker_range()){
            add_blocker(*blocker);
        }
    }

    Skeleton_blocker_complex& operator=(const Skeleton_blocker_complex& copy){
        clear();
        visitor = NULL;
        degree_ = copy.degree_;
        skeleton = Graph(copy.skeleton);
        num_vertices_ = copy.num_vertices_;

        num_blockers_ = 0;
        // we copy the blockers
        for (auto blocker : copy.const_blocker_range())
            add_blocker(*blocker);
        return *this;
    }


    virtual ~Skeleton_blocker_complex(){
        clear();
    }

    virtual void clear(){
        // xxx for now the responsabilty of freeing the visitor is for
        // the user
        visitor = NULL;

        degree_.clear();
        num_vertices_ =0;

        remove_blockers();

        skeleton.clear();
    }

    void set_visitor(Visitor* other_visitor){
        visitor = other_visitor;
    }







public:

    Vertex_handle operator[](Root_vertex_handle global) const{
        auto local(get_address(global));
        assert(local);
        return *local;
    }

    Graph_vertex& operator[](Vertex_handle address){
        assert(0<=address.vertex && address.vertex< boost::num_vertices(skeleton));
        return skeleton[address.vertex];
    }

    const Graph_vertex& operator[](Vertex_handle address) const{
        assert(0<=address.vertex && address.vertex< boost::num_vertices(skeleton));
        return skeleton[address.vertex];
    }

    Vertex_handle add_vertex(){
        Vertex_handle address(boost::add_vertex(skeleton));
        num_vertices_++;
        (*this)[address].activate();
        // safe since we now that we are in the root complex and the field 'address' and 'id'
        // are identical for every vertices
        (*this)[address].set_id(Root_vertex_handle(address.vertex));
        degree_.push_back(0);
        if (visitor) visitor->on_add_vertex(address);
        return address;
    }


    void remove_vertex(Vertex_handle address){
        assert(contains_vertex(address));
        // We remove b
        boost::clear_vertex(address.vertex,skeleton);
        (*this)[address].deactivate();
        num_vertices_--;
        degree_[address.vertex]=-1;
        if (visitor) visitor->on_remove_vertex(address);
    }

    bool contains_vertex(Vertex_handle u) const{
        if (u.vertex<0 || u.vertex>=boost::num_vertices(skeleton)) return false;
        return (*this)[u].is_active();
    }

    bool contains_vertex(Root_vertex_handle u) const{
        boost::optional<Vertex_handle> address = get_address(u);
        return address &&  (*this)[*address].is_active();
    }


    bool contains_vertices(const Simplex_handle & sigma) const{
        for (auto vertex : sigma)
            if(!contains_vertex(vertex)) return false;
        return true;
    }

    virtual boost::optional<Vertex_handle> get_address(Root_vertex_handle id) const{
        boost::optional<Vertex_handle> res;
        if ( id.vertex< boost::num_vertices(skeleton) ) res = Vertex_handle(id.vertex);//xxx
        return res;
    }



    Root_vertex_handle get_id(Vertex_handle local) const{
        assert(0<=local.vertex && local.vertex< boost::num_vertices(skeleton));
        return (*this)[local].get_id();
    }


    Vertex_handle convert_handle_from_another_complex(
            const Skeleton_blocker_complex& other,Vertex_handle vh_in_other) const{
        auto vh_in_current_complex = get_address(other.get_id(vh_in_other));
        assert(vh_in_current_complex);
        return *vh_in_current_complex;
    }

    int degree(Vertex_handle local) const{
        assert(0<=local.vertex && local.vertex< boost::num_vertices(skeleton));
        return degree_[local.vertex];
    }




public:

    boost::optional<Edge_handle> operator[](const std::pair<Vertex_handle,Vertex_handle>& ab) const{
        boost::optional<Edge_handle> res;
        std::pair<Edge_handle,bool> edge_pair(boost::edge(ab.first.vertex,ab.second.vertex,skeleton));
        if (edge_pair.second)
            res = edge_pair.first;
        return res;
    }

    Graph_edge& operator[](Edge_handle edge_handle){
        return skeleton[edge_handle];
    }

    const Graph_edge& operator[](Edge_handle edge_handle) const{
        return skeleton[edge_handle];
    }

    Vertex_handle first_vertex(Edge_handle edge_handle) const{
        return source(edge_handle,skeleton);
    }

    Vertex_handle second_vertex(Edge_handle edge_handle) const{
        return target(edge_handle,skeleton);
    }

    Simplex_handle get_vertices(Edge_handle edge_handle) const{
        auto edge((*this)[edge_handle]);
        return Simplex_handle((*this)[edge.first()],(*this)[edge.second()]);
    }

    Edge_handle add_edge(Vertex_handle a, Vertex_handle b){
        assert(contains_vertex(a) && contains_vertex(b));
        assert(a!=b);

        auto edge_handle((*this)[std::make_pair(a,b)]);
        //      std::pair<Edge_handle,bool> pair_descr_bool = (*this)[std::make_pair(a,b)];
        //      Edge_handle edge_descr;
        //      bool edge_present = pair_descr_bool.second;
        if (!edge_handle)
        {
            edge_handle = boost::add_edge(a.vertex,b.vertex,skeleton).first;
            (*this)[*edge_handle].setId(get_id(a),get_id(b));
            degree_[a.vertex]++;
            degree_[b.vertex]++;
            if (visitor) visitor->on_add_edge(a,b);
        }
        return *edge_handle;
    }

    void add_edges(const Simplex_handle & sigma){
        Simplex_handle_iterator i, j;
        for (i = sigma.begin() ; i != sigma.end() ; ++i)
            for (j = i, j++ ; j != sigma.end() ; ++j)
                add_edge(*i,*j);
    }

    virtual Edge_handle remove_edge(Vertex_handle a, Vertex_handle b){
        bool found;
        Edge_handle edge;
        tie(edge,found) = boost::edge(a.vertex,b.vertex,skeleton);
        if (found)
        {
            if (visitor) visitor->on_remove_edge(a,b);
            //      if (heapCollapse.Contains(edge)) heapCollapse.Delete(edge);
            boost::remove_edge(a.vertex,b.vertex,skeleton);
            degree_[a.vertex]--;
            degree_[b.vertex]--;
        }
        return edge;
    }


    void remove_edge(Edge_handle edge){
        assert(contains_vertex(first_vertex(edge)));
        assert(contains_vertex(second_vertex(edge)));
        remove_edge(first_vertex(edge),second_vertex(edge));
    }


    void keep_only_vertices(){
        remove_blockers();

        for(auto u : vertex_range()){
            while (this->degree(u)> 0)
            {
                Vertex_handle v(*(adjacent_vertices(u.vertex, this->skeleton).first));
                this->remove_edge(u,v);
            }
        }
    }


    bool contains_edge(Vertex_handle a, Vertex_handle b) const{
        //if (a.vertex<0 || b.vertex <0) return false;
        return boost::edge(a.vertex,b.vertex,skeleton).second;
    }


    bool contains_edges(const Simplex_handle & sigma) const{
        for (auto i = sigma.begin() ; i != sigma.end() ; ++i){
            if(!contains_vertex(*i)) return false;
            for (auto j=i; ++j != sigma.end() ; ){
                if (!contains_edge(*i,*j))
                    return false;
            }
        }
        return true;
    }





    Blocker_handle add_blocker(const Simplex_handle& blocker){
        assert(blocker.dimension()>1);
        if (contains_blocker(blocker))
        {
            //std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
            return 0;
        }
        else{
            if (visitor) visitor->on_add_blocker(blocker);
            Blocker_handle blocker_pt = new Simplex_handle(blocker);
            num_blockers_++;
            auto vertex = blocker_pt->begin();
            while(vertex != blocker_pt->end())
            {
                blocker_map_.insert(BlockerPair(*vertex,blocker_pt));
                ++vertex;
            }
            return blocker_pt;
        }
    }

protected:
    void add_blocker(Blocker_handle blocker){
        if (contains_blocker(*blocker))
        {
            //std::cerr << "ATTEMPT TO ADD A BLOCKER ALREADY THERE ---> BLOCKER IGNORED" << endl;
            return;
        }
        else{
            if (visitor) visitor->on_add_blocker(*blocker);
            num_blockers_++;
            auto vertex = blocker->begin();
            while(vertex != blocker->end())
            {
                blocker_map_.insert(BlockerPair(*vertex,blocker));
                ++vertex;
            }
        }
    }


protected:
    void remove_blocker(const Blocker_handle sigma, Vertex_handle v){
        Complex_blocker_around_vertex_iterator blocker;
        for (blocker = blocker_range(v).begin();
                blocker != blocker_range(v).end();
                ++blocker
        ){
            if (*blocker == sigma) break;
        }
        if (*blocker != sigma){
            std::cerr << "bug ((*blocker).second == sigma) ie try to remove a blocker not present\n";
            assert(false);
        }
        else{
            blocker_map_.erase(blocker.current_position());
        }
    }

public:
    void remove_blocker(const Blocker_handle sigma){
        for (auto vertex : *sigma){
            remove_blocker(sigma,vertex);
        }
        num_blockers_--;
    }




    void remove_blockers(){
        // Desallocate the blockers
        while (!blocker_map_.empty()){
            delete_blocker(blocker_map_.begin()->second);
        }
        num_blockers_ = 0;
        blocker_map_.clear();
    }

protected:
    void remove_blocker(const Simplex_handle& sigma){
        assert(contains_blocker(sigma));
        for (auto vertex : sigma)
            remove_blocker(sigma,vertex);
        num_blockers_--;
    }




public:
    void delete_blocker(Blocker_handle sigma){
        if (visitor) visitor->on_delete_blocker(sigma);
        remove_blocker(sigma);
        delete sigma;
    }


    bool contains_blocker(const Blocker_handle s) const{
        if (s->dimension()<2)
            return false;

        Vertex_handle a = s->first_vertex();

        for (auto blocker : const_blocker_range(a)){
            if ( s == *blocker )
                return true;
        }
        return false;
    }

    bool contains_blocker(const Simplex_handle & s) const{
        if (s.dimension()<2)
            return false;

        Vertex_handle a = s.first_vertex();

        for (auto blocker : const_blocker_range(a)){
            if ( s == *blocker )
                return true;
        }
        return false;
    }


private:
    bool blocks(const Simplex_handle & sigma) const{

        for(auto blocker : const_blocker_range()){
            if ( sigma.contains(*blocker) )
                return true;
        }
        return false;
    }




protected:
    virtual void add_neighbours(Vertex_handle v, Simplex_handle & n,bool keep_only_superior=false) const{
        boost_adjacency_iterator ai, ai_end;
        for (tie(ai, ai_end) = adjacent_vertices(v.vertex, skeleton); ai != ai_end; ++ai){
            if (keep_only_superior){
                if (*ai>v.vertex)
                    n.add_vertex(Vertex_handle(*ai));
            }
            else
                n.add_vertex(Vertex_handle(*ai));
        }
    }

    virtual void add_neighbours(const Simplex_handle &alpha, Simplex_handle & res,bool keep_only_superior=false) const{
        res.clear();
        auto alpha_vertex = alpha.begin();
        add_neighbours(*alpha_vertex,res,keep_only_superior);
        for (alpha_vertex = (alpha.begin())++ ; alpha_vertex != alpha.end() ; ++alpha_vertex)
            keep_neighbours(*alpha_vertex,res,keep_only_superior);
    }

    virtual void keep_neighbours(Vertex_handle v, Simplex_handle& res,bool keep_only_superior=false) const{
        Simplex_handle nv;
        add_neighbours(v,nv,keep_only_superior);
        res.intersection(nv);
    }

    virtual void remove_neighbours(Vertex_handle v, Simplex_handle & res,bool keep_only_superior=false) const{
        Simplex_handle nv;
        add_neighbours(v,nv,keep_only_superior);
        res.difference(nv);
    }


public:

     //xxx rename get_address et place un using dans sub_complex
    boost::optional<Simplex_handle> get_simplex_address(const Root_simplex_handle& s) const
    {
        boost::optional<Simplex_handle> res;

        Simplex_handle s_address;
        //Root_simplex_const_iterator i;
        for (auto i = s.begin() ; i != s.end() ; ++i)
        {
            boost::optional<Vertex_handle> address = get_address(*i);
            if (!address)
                return res;
            else
                s_address.add_vertex(*address);
        }
        res = s_address;
        return res;
    }

    Root_simplex_handle get_id(const Simplex_handle& local_simplex) const{
        Root_simplex_handle global_simplex;
        for (auto x = local_simplex.begin(); x!= local_simplex.end();++x){
            global_simplex.add_vertex(get_id(*x));

        }
        return global_simplex;

    }


    virtual bool contains(const Simplex_handle & s) const{
        if (s.dimension() == -1 ) return false;
        else
            if (s.dimension() ==0 ){
                return contains_vertex(s.first_vertex());
            }
            else
                return ( contains_edges(s) && !blocks(s) );
    }

    /*
     * @brief returnrs true iff the complex is empty.
     */
    bool empty() const{
        return num_vertices()==0;
    }

    /*
     * @brief returns the number of vertices in the complex.
     */
    int num_vertices() const{
        //remark boost::num_vertices(skeleton) counts deactivated vertices
        return num_vertices_;
    }

    /*
     * @brief returns the number of edges in the complex.
     * @details currently in O(n)
     */
    //  todo cache the value
    int num_edges() const{
        return boost::num_edges(skeleton);
    }

    /*
     * @brief returns the number of blockers in the complex.
     */
    int num_blockers() const{
        return num_blockers_;
    }

    /*
     * @brief returns true iff the graph of the 1-skeleton of the complex is complete.
     */
    bool complete() const{
        return (num_vertices()*(num_vertices()-1))/2 == num_edges();
    }

    int num_connected_components() const{
        int num_vert_collapsed = skeleton.vertex_set().size() - num_vertices();
        std::vector<int> component(skeleton.vertex_set().size());
        return boost::connected_components(this->skeleton,&component[0]) - num_vert_collapsed;
    }

    bool is_cone() const{
        if (num_vertices()==0) return false;
        if (num_vertices()==1) return true;
        for(auto vi :  vertex_range()){
            //xxx todo faire une methode bool is_in_blocker(Vertex_handle)
            if (blocker_map_.find(vi)==blocker_map_.end()){
                // no blocker passes through the vertex, we just need to
                // check if the current vertex is linked to all others vertices of the complex
                if (degree_[vi.vertex] == num_vertices()-1)
                    return true;
            }
        }
        return false;
    }





    typedef Complex_vertex_iterator<Skeleton_blocker_complex> CVI; //todo rename

    //  /**
    //   * @brief Range over the vertices of the simplicial complex.
    //   * Methods .begin() and .end() return a Complex_vertex_iterator.
    //   */
    typedef boost::iterator_range     < Complex_vertex_iterator<Skeleton_blocker_complex> > Complex_vertex_range;

    Complex_vertex_range vertex_range() const
    {
        auto begin = Complex_vertex_iterator<Skeleton_blocker_complex>(this);
        auto end = Complex_vertex_iterator<Skeleton_blocker_complex>(this,0);
        return Complex_vertex_range(begin,end);
    }

    typedef boost::iterator_range     < Complex_neighbors_vertices_iterator<Skeleton_blocker_complex> > Complex_neighbors_vertices_range;

    Complex_neighbors_vertices_range vertex_range(Vertex_handle v) const
    {
        auto begin = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(this,v);
        auto end = Complex_neighbors_vertices_iterator<Skeleton_blocker_complex>(this,v,0);
        return Complex_neighbors_vertices_range(begin,end);
    }





    typedef boost::iterator_range  <Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
            Complex_edge_range;

    Complex_edge_range edge_range() const
    {
        auto begin = Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this);
        auto end = Complex_edge_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,0);
        return Complex_edge_range(begin,end);
    }



    typedef boost::iterator_range  <Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>>
            Complex_edge_around_vertex_range;
    Complex_edge_around_vertex_range edge_range(Vertex_handle v) const
    {
        auto begin = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,v);
        auto end = Complex_edge_around_vertex_iterator<Skeleton_blocker_complex<SkeletonBlockerDS>>(this,v,0);
        return Complex_edge_around_vertex_range(begin,end);
    }





private:
    typedef Skeleton_blocker_link_complex<Skeleton_blocker_complex<SkeletonBlockerDS> > Link;
    typedef Skeleton_blocker_link_superior<Skeleton_blocker_complex<SkeletonBlockerDS> > Superior_link;
public:
    typedef Triangle_around_vertex_iterator<Skeleton_blocker_complex,Superior_link> Superior_triangle_around_vertex_iterator;


    typedef boost::iterator_range     < Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link> > Complex_triangle_around_vertex_range;

    Complex_triangle_around_vertex_range triangle_range(Vertex_handle v) const
    {
        auto begin = Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link>(this,v);
        auto end = Triangle_around_vertex_iterator<Skeleton_blocker_complex,Link>(this,v,0);
        return Complex_triangle_around_vertex_range(begin,end);
    }


    typedef boost::iterator_range<Triangle_iterator<Skeleton_blocker_complex> > Complex_triangle_range;

    Complex_triangle_range triangle_range() const
    {
        auto end = Triangle_iterator<Skeleton_blocker_complex>(this,0);
        if(empty()){
            return Complex_triangle_range(end,end);
        }
        else{
            auto begin = Triangle_iterator<Skeleton_blocker_complex>(this);
            return Complex_triangle_range(begin,end);
        }

    }





    typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex,Link> Complex_simplex_around_vertex_iterator;

    typedef boost::iterator_range     < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;

    Complex_simplex_around_vertex_range simplex_range(Vertex_handle v) const
    {
        assert(contains_vertex(v));
        return Complex_simplex_around_vertex_range(
                Complex_simplex_around_vertex_iterator(this,v),
                Complex_simplex_around_vertex_iterator(this,v,true)
        );
    }

    //  typedef Simplex_iterator<Skeleton_blocker_complex,Superior_link> Complex_simplex_iterator;
    typedef Simplex_iterator<Skeleton_blocker_complex> Complex_simplex_iterator;

    typedef boost::iterator_range     < Complex_simplex_iterator > Complex_simplex_range;

    Complex_simplex_range simplex_range() const
    {
        Complex_simplex_iterator end(this,true);
        if(empty()){
            return Complex_simplex_range(end,end);
        }
        else{
            Complex_simplex_iterator begin(this);
            return Complex_simplex_range(begin,end);
        }
    }




private:
    typedef Blocker_iterator_around_vertex_internal<
            typename std::multimap<Vertex_handle,Simplex_handle *>::iterator,
            Blocker_handle>
    Complex_blocker_around_vertex_iterator;

    typedef Blocker_iterator_around_vertex_internal<
            typename std::multimap<Vertex_handle,Simplex_handle *>::const_iterator,
            const Blocker_handle>
    Const_complex_blocker_around_vertex_iterator;

    typedef boost::iterator_range  <Complex_blocker_around_vertex_iterator> Complex_blocker_around_vertex_range;
    typedef boost::iterator_range  <Const_complex_blocker_around_vertex_iterator>   Const_complex_blocker_around_vertex_range;

public:


    Complex_blocker_around_vertex_range blocker_range(Vertex_handle v)
    {
        auto begin = Complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
        auto end = Complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
        return Complex_blocker_around_vertex_range(begin,end);
    }

    Const_complex_blocker_around_vertex_range const_blocker_range(Vertex_handle v) const
    {
        auto begin = Const_complex_blocker_around_vertex_iterator(blocker_map_.lower_bound(v));
        auto end = Const_complex_blocker_around_vertex_iterator(blocker_map_.upper_bound(v));
        return Const_complex_blocker_around_vertex_range(begin,end);
    }




private:

    typedef Blocker_iterator_internal<
            typename std::multimap<Vertex_handle,Simplex_handle *>::iterator,
            Blocker_handle>
    Complex_blocker_iterator;

    typedef Blocker_iterator_internal<
            typename std::multimap<Vertex_handle,Simplex_handle *>::const_iterator,
            const Blocker_handle>
    Const_complex_blocker_iterator;

    typedef boost::iterator_range  <Complex_blocker_iterator>   Complex_blocker_range;
    typedef boost::iterator_range  <Const_complex_blocker_iterator> Const_complex_blocker_range;


public:

    Complex_blocker_range blocker_range()
    {
        auto begin = Complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end() );
        auto end = Complex_blocker_iterator(blocker_map_.end() , blocker_map_.end() );
        return Complex_blocker_range(begin,end);
    }

    Const_complex_blocker_range const_blocker_range() const
    {
        auto begin = Const_complex_blocker_iterator(blocker_map_.begin(), blocker_map_.end() );
        auto end = Const_complex_blocker_iterator(blocker_map_.end() , blocker_map_.end() );
        return Const_complex_blocker_range(begin,end);
    }










public:

    std::string to_string() const{
        std::ostringstream stream;
        stream<<num_vertices()<<" vertices:\n"<<vertices_to_string()<<std::endl;
        stream<<num_edges()<<" edges:\n"<<edges_to_string()<<std::endl;
        stream<<num_blockers()<<" blockers:\n"<<blockers_to_string()<<std::endl;
        return stream.str();
    }

    std::string vertices_to_string() const{
        std::ostringstream stream;
        for(auto vertex : vertex_range())
            stream << "("<<(*this)[vertex].get_id()<<"),";
        stream<< std::endl;
        return stream.str();
    }

    std::string edges_to_string() const{
        std::ostringstream stream;
        for(auto edge : edge_range())
            stream << "("<< (*this)[edge].first()<<","<< (*this)[edge].second() << ")"<<" id = "<< (*this)[edge].index()<< std::endl;
        stream<< std::endl;
        return stream.str();
    }


    std::string blockers_to_string() const{
        std::ostringstream stream;
        for (auto bl:blocker_map_){
            stream << bl.first << " => " << bl.second << ":"<<*bl.second <<"\n";
        }
        return stream.str();
    }


};



template<typename Complex,typename SimplexHandleIterator>
unsigned make_complex_from_top_faces(Complex& complex,SimplexHandleIterator begin,SimplexHandleIterator end){
    typedef typename Complex::Simplex_handle Simplex_handle;

    int dimension = 0;
    for(auto top_face = begin; top_face != end; ++top_face)
        dimension = std::max(dimension,top_face->dimension());

    std::vector< std::set<Simplex_handle> > simplices_per_dimension(dimension+1);


    for(auto top_face = begin; top_face != end; ++top_face){
        register_faces(simplices_per_dimension,*top_face);
    }

    // compute list of simplices
    std::list<Simplex_handle> simplices;
    for(int dim = 0 ; dim <= dimension ; ++dim){
        std::copy(
                simplices_per_dimension[dim].begin(),
                simplices_per_dimension[dim].end(),
                std::back_inserter(simplices)
        );
    }
    complex = Complex(simplices);
    return simplices.size();
}

} // namespace skbl

}  // namespace GUDHI



#endif