/*    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 SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_COMPLEX_H_
#define SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_COMPLEX_H_

#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 <iostream>
#include <fstream>
#include <sstream>
#include <memory>
#include <map>
#include <list>
#include <set>
#include <vector>
#include <string>
#include <algorithm>
#include <utility>

#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/Skeleton_blocker/internal/Trie.h"

#include "gudhi/Utils.h"

namespace Gudhi {

namespace skbl {

/**
 *@class Skeleton_blocker_complex
 *@brief Abstract Simplicial Complex represented with a skeleton/blockers pair.
 *@ingroup 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:
	/**
	 * @brief The type of stored vertex node, specified by the template SkeletonBlockerDS
	 */
	typedef typename SkeletonBlockerDS::Graph_vertex Graph_vertex;

	/**
	 * @brief The type of stored  edge node, specified by the template SkeletonBlockerDS
	 */
	typedef typename SkeletonBlockerDS::Graph_edge Graph_edge;

	typedef typename SkeletonBlockerDS::Root_vertex_handle Root_vertex_handle;

	/**
	 * @brief The type of an handle to a vertex of the complex.
	 */
	typedef typename SkeletonBlockerDS::Vertex_handle Vertex_handle;
	typedef typename Root_vertex_handle::boost_vertex_handle boost_vertex_handle;

	/**
	 * @brief A ordered set of integers that represents a simplex.
	 */
	typedef Skeleton_blocker_simplex<Vertex_handle> Simplex_handle;
	typedef Skeleton_blocker_simplex<Root_vertex_handle> Root_simplex_handle;

	/**
	 * @brief Handle to a blocker of the complex.
	 */
	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:
	/**
	 * @brief Handle to an edge of the complex.
	 */
	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;

	/**
	 * @details If 'x' is a Vertex_handle of a vertex in the complex then degree[x] = d is its degree.
	 *
	 * This quantity is updated when adding/removing edge.
	 *
	 * This is useful because the operation
	 * list.size() is done in linear time.
	 */
	std::vector<boost_vertex_handle> degree_;
	Graph skeleton; /** 1-skeleton of the simplicial complex. */

	/** Each vertex can access to the blockers passing through it. */
	BlockerMap blocker_map_;

public:
	/////////////////////////////////////////////////////////////////////////////
	/** @name Constructors, Destructors
	 */
	//@{
	/**
	 *@brief constructs a simplicial complex with a given number of vertices and a visitor.
	 */
	explicit Skeleton_blocker_complex(int num_vertices_ = 0, Visitor* visitor_ = NULL)
	: visitor(visitor_) {
		clear();
		for (int i = 0; i < num_vertices_; ++i) {
			add_vertex();
		}
	}

private:
	typedef Trie<Skeleton_blocker_complex<SkeletonBlockerDS>> STrie;
	template<typename SimpleHandleOutputIterator>
	/**
	 * return a vector of simplex trie for every vertex
	 */
	std::vector<STrie*> compute_cofaces(SimpleHandleOutputIterator simplex_begin,SimpleHandleOutputIterator simplex_end){
		std::vector<STrie*> cofaces(num_vertices(),0);
		for(auto i = 0u ; i < num_vertices(); ++i)
			cofaces[i] = new STrie(Vertex_handle(i));
		for(auto s_it = simplex_begin; s_it != simplex_end; ++s_it){
			if(s_it->dimension()>=1)
				cofaces[s_it->first_vertex()]->add_simplex(*s_it);
		}
		return cofaces;
	}


public:
	/**
	 * @brief Constructor with a list of simplices.
	 * @details is_flag_complex indicates if the complex is a flag complex or not (to know if blockers have to be computed or not).
	 */
	template<typename SimpleHandleOutputIterator>
	Skeleton_blocker_complex(SimpleHandleOutputIterator simplex_begin,SimpleHandleOutputIterator simplex_end,bool is_flag_complex = false,Visitor* visitor_ = NULL)
	: num_vertices_(0),
	  num_blockers_(0),
	  visitor(visitor_){
		std::vector<std::pair<Vertex_handle,Vertex_handle>> edges;

		//first pass to add vertices and edges
		int num_vertex = -1;
		for(auto s_it = simplex_begin; s_it != simplex_end; ++s_it){
			if(s_it->dimension()==0) num_vertex = (std::max)(num_vertex,s_it->first_vertex().vertex);
			if(s_it->dimension()==1) edges.emplace_back(s_it->first_vertex(),s_it->last_vertex());
		}
		while(num_vertex-->=0) add_vertex();

		for(const auto& e : edges)
			add_edge(e.first,e.second);

		if(!is_flag_complex){
			//need to compute blockers
			std::vector<STrie*> cofaces(compute_cofaces(simplex_begin,simplex_end));
			std::vector<Simplex_handle> max_faces;

			for(auto i = 0u ; i < num_vertices(); ++i){
				auto max_faces_i = cofaces[i]->maximal_faces();
				max_faces.insert(max_faces.end(),max_faces_i.begin(),max_faces_i.end());
			}

			// all time here

			// for every maximal face s, it picks its first vertex v and check for all nv \in N(v)
			// if s union nv is in the complex, if not it is a blocker.
			for(auto max_face : max_faces){
				if(max_face.dimension()>0){
					auto first_v = max_face.first_vertex();
					for(auto nv : vertex_range(first_v)){
						if(! max_face.contains(nv) && max_face.first_vertex()<nv){
							//check that all edges in vertices(max_face)\cup nv are here
							//since max_face is a simplex, we only need to check that edges
							// (x nv) where x \in max_face are present
							bool all_edges_here = true;
							for(auto x : max_face)
								if(!contains_edge(x,nv)){
									all_edges_here = false;
									break;
								}
							if(all_edges_here){ //eg  this->contains(max_face)
								max_face.add_vertex(nv);
								if(!cofaces[first_v]->contains(max_face)){
									// if there exists a blocker included in max_face, we remove it
									// as it is not a minimum missing face
									// the other alternative would be to check to all properfaces
									// are in the complex before adding a blocker but that
									// would be more expensive if there are few blockers
									std::vector<Blocker_handle> blockers_to_remove;
									for(auto b : blocker_range(first_v))
										if(b->contains(max_face))
											blockers_to_remove.push_back(b);
									for(auto b : blockers_to_remove)
										this->delete_blocker(b);

									add_blocker(max_face);
								}
								max_face.remove_vertex(nv);
							}
						}
					}
				}
			}
			for(auto i = 0u ; i < num_vertices(); ++i)
				delete(cofaces[i]);
		}
	}

	// 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;
	}

	/**
	 * The destructor delete all blockers allocated.
	 */
	virtual ~Skeleton_blocker_complex() {
		clear();
	}

	/**
	 * @details Clears the simplicial complex. After a call to this function,
	 * blockers are destroyed. The 1-skeleton and the set of blockers
	 * are both empty.
	 */
	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();
	}

	/**
	 *@brief allows to change the visitor.
	 */
	void set_visitor(Visitor* other_visitor) {
		visitor = other_visitor;
	}

	//@}

	/////////////////////////////////////////////////////////////////////////////
	/** @name Vertices operations
	 */
	//@{
public:
	/**
	 * @brief Return a local Vertex_handle of a vertex given a global one.
	 * @remark Assume that the vertex is present in the complex.
	 */
	Vertex_handle operator[](Root_vertex_handle global) const {
		auto local(get_address(global));
		assert(local);
		return *local;
	}

	/**
	 * @brief Return the vertex node associated to local Vertex_handle.
	 * @remark Assume that the vertex is present in the complex.
	 */
	Graph_vertex& operator[](Vertex_handle address) {
		assert(
				0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
		return skeleton[address.vertex];
	}

	/**
	 * @brief Return the vertex node associated to local Vertex_handle.
	 * @remark Assume that the vertex is present in the complex.
	 */
	const Graph_vertex& operator[](Vertex_handle address) const {
		assert(
				0 <= address.vertex && address.vertex < boost::num_vertices(skeleton));
		return skeleton[address.vertex];
	}

	/**
	 * @brief Adds a vertex to the simplicial complex and returns its Vertex_handle.
	 */
	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;
	}

	/**
	 * @brief Remove a vertex from the simplicial complex
	 * @remark It just deactivates the vertex with a boolean flag but does not
	 * remove it from vertices from complexity issues.
	 */
	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();
	}

	/**
	 * @return true iff the simplicial complex contains all vertices
	 * of simplex sigma
	 */
	bool contains_vertices(const Simplex_handle & sigma) const {
		for (auto vertex : sigma)
			if (!contains_vertex(vertex))
				return false;
		return true;
	}

	/**
	 * @brief Given an Id return the address of the vertex having this Id in the complex.
	 * @remark For a simplicial complex, the address is the id but it may not be the case for a SubComplex.
	 */
	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;
	}

	/**
	 * return the id of a vertex of adress local present in the graph
	 */
	Root_vertex_handle get_id(Vertex_handle local) const {
		assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
		return (*this)[local].get_id();
	}

	/**
	 * @brief Convert an address of a vertex of a complex to the address in
	 * the current complex.
	 * @details
	 * If the current complex is a sub (or sup) complex of 'other', it converts
	 * the address of a vertex v expressed in 'other' to the address of the vertex
	 * v in the current one.
	 * @remark this methods uses Root_vertex_handle to identify the vertex and
	 * assumes the vertex is present in the current complex.
	 */
	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;
	}

	/**
	 * @brief return the graph degree of a vertex.
	 */
	int degree(Vertex_handle local) const {
		assert(0 <= local.vertex && local.vertex < boost::num_vertices(skeleton));
		return degree_[local.vertex];
	}

	//@}

	/////////////////////////////////////////////////////////////////////////////
	/** @name Edges operations
	 */
	//@{
public:
	/**
	 * @brief return an edge handle if the two vertices forms
	 * an edge in the complex
	 */
	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;
	}

	/**
	 * @brief returns the stored node associated to an edge
	 */
	Graph_edge& operator[](Edge_handle edge_handle) {
		return skeleton[edge_handle];
	}

	/**
	 * @brief returns the stored node associated to an edge
	 */
	const Graph_edge& operator[](Edge_handle edge_handle) const {
		return skeleton[edge_handle];
	}

	/**
	 * @brief returns the first vertex of an edge
	 * @details it assumes that the edge is present in the complex
	 */
	Vertex_handle first_vertex(Edge_handle edge_handle) const {
		return static_cast<Vertex_handle>(source(edge_handle, skeleton));
	}

	/**
	 * @brief returns the first vertex of an edge
	 * @details it assumes that the edge is present in the complex
	 */
	Vertex_handle second_vertex(Edge_handle edge_handle) const {
		return static_cast<Vertex_handle>(target(edge_handle, skeleton));
	}

	/**
	 * @brief returns the simplex made with the two vertices of the edge
	 * @details it assumes that the edge is present in the complex

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

	/**
	 * @brief Adds an edge between vertices a and b and all its cofaces.
	 */
	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;
	}

	/**
	 * @brief Adds all edges and their cofaces of a simplex to the simplicial complex.
	 */
	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);
	}

	/**
	 * @brief Removes an edge from the simplicial complex and all its cofaces.
	 * @details returns the former Edge_handle representing the edge
	 */
	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;
	}

	/**
	 * @brief Removes edge and its cofaces from the simplicial complex.
	 */
	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));
	}

	/**
	 * @brief The complex is reduced to its set of vertices.
	 * All the edges and blockers are removed.
	 */
	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);
			}
		}
	}

	/**
	 * @return true iff the simplicial complex contains an edge between
	 * vertices a and b
	 */
	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;
	}

	/**
	 * @return true iff the simplicial complex contains all vertices
	 * and all edges of simplex sigma
	 */
	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;
	}
	//@}

	/////////////////////////////////////////////////////////////////////////////
	/** @name Blockers operations
	 */
	//@{
	/**
	 * @brief Adds the simplex to the set of blockers and
	 * returns a Blocker_handle toward it if was not present before and 0 otherwise.
	 */
	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:
	/**
	 * @brief Adds the simplex to the set of blockers
	 */
	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:
	/**
	 * Removes sigma from the blocker map of vertex v
	 */
	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:
	/**
	 * @brief Removes the simplex from the set of blockers.
	 * @remark sigma has to belongs to the set of blockers
	 */
	void remove_blocker(const Blocker_handle sigma) {
		for (auto vertex : *sigma)
			remove_blocker(sigma, vertex);
		num_blockers_--;
	}

	/**
	 * @brief Remove all blockers, in other words, it expand the simplicial
	 * complex to the smallest flag complex that contains it.
	 */
	void remove_blockers() {
		// Desallocate the blockers
		while (!blocker_map_.empty()) {
			delete_blocker(blocker_map_.begin()->second);
		}
		num_blockers_ = 0;
		blocker_map_.clear();
	}

protected:
	/**
	 * Removes the simplex sigma from the set of blockers.
	 * sigma has to belongs to the set of blockers
	 *
	 * @remark contrarily to delete_blockers does not call the destructor
	 */
	void remove_blocker(const Simplex_handle& sigma) {
		assert(contains_blocker(sigma));
		for (auto vertex : sigma)
			remove_blocker(sigma, vertex);
		num_blockers_--;
	}

public:
	/**
	 * Removes the simplex s from the set of blockers
	 * and desallocate s.
	 */
	void delete_blocker(Blocker_handle sigma) {
		if (visitor)
			visitor->on_delete_blocker(sigma);
		remove_blocker(sigma);
		delete sigma;
	}

	/**
	 * @return true iff s is a blocker of the simplicial complex
	 */
	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;
	}

	/**
	 * @return true iff s is a blocker of the simplicial complex
	 */
	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:
	/**
	 * @return true iff a blocker of the simplicial complex
	 * is a face of sigma.
	 */
	bool blocks(const Simplex_handle & sigma) const {
		for (auto blocker : const_blocker_range()) {
			if (sigma.contains(*blocker))
				return true;
		}
		return false;
	}

	//@}

protected:
	/**
	 * @details Adds to simplex the neighbours of v e.g. \f$ n \leftarrow n \cup N(v) \f$.
	 * If keep_only_superior is true then only vertices that are greater than v are added.
	 */
	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));
			}
		}
	}

	/**
	 * @details Add to simplex res all vertices which are
	 * neighbours of alpha: ie \f$ res \leftarrow res \cup N(alpha) \f$.
	 *
	 * If 'keep_only_superior' is true then only vertices that are greater than alpha are added.
	 * todo revoir
	 *
	 */
	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);
	}

	/**
	 * @details Remove from simplex n all vertices which are
	 * not neighbours of v e.g. \f$ res \leftarrow res \cap N(v) \f$.
	 * If 'keep_only_superior' is true then only vertices that are greater than v are keeped.
	 */
	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);
	}

	/**
	 * @details Remove from simplex all vertices which are
	 * neighbours of v eg \f$ res \leftarrow res \setminus N(v) \f$.
	 * If 'keep_only_superior' is true then only vertices that are greater than v are added.
	 */
	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:
	/**
	 * @brief Compute the local vertices of 's' in the current complex
	 * If one of them is not present in the complex then the return value is uninitialized.
	 *
	 *
	 */
	// 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;
	}

	/**
	 * @brief returns a simplex with vertices which are the id of vertices of the
	 * argument.
	 */
	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;
	}

	/**
	 * @brief returns true iff the simplex s belongs to the simplicial
	 * complex.
	 */
	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();
	}

	/**
	 * @brief returns the number of connected components in the graph of the 1-skeleton.
	 */
	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;
	}

	/**
	 * @brief %Test if the complex is a cone.
	 * @details Runs in O(n) where n is the number of vertices.
	 */
	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;
	}

	//@}

	/////////////////////////////////////////////////////////////////////////////
	/** @name Vertex iterators
	 */
	//@{
	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;

	/**
	 * @brief Returns a Complex_vertex_range over all vertices of the complex
	 */
	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;

	/**
	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes trough v
	 */
	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);
	}

	//@}

	/** @name Edge iterators
	 */
	//@{

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

	/**
	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex
	 */
	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;
	/**
	 * @brief Returns a Complex_edge_range over all edges of the simplicial complex that passes
	 * through 'v'
	 */
	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);
	}

	//@}

	/** @name Triangles iterators
	 */
	//@{
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;

	/**
	 * @brief Range over triangles around a vertex of the simplicial complex.
	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
	 *
	 */
	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;

	/**
	 * @brief Range over triangles of the simplicial complex.
	 * Methods .begin() and .end() return a Triangle_around_vertex_iterator.
	 *
	 */
	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);
		}
	}

	//@}

	/** @name Simplices iterators
	 */
	//@{
	typedef Simplex_around_vertex_iterator<Skeleton_blocker_complex, Link> Complex_simplex_around_vertex_iterator;

	/**
	 * @brief Range over the simplices of the simplicial complex around a vertex.
	 * Methods .begin() and .end() return a Complex_simplex_around_vertex_iterator.
	 */
	typedef boost::iterator_range < Complex_simplex_around_vertex_iterator > Complex_simplex_around_vertex_range;

	/**
	 * @brief Returns a Complex_simplex_around_vertex_range over all the simplices around a vertex of the complex
	 */
	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;

	/**
	 * @brief Returns a Complex_simplex_range over all the simplices of the complex
	 */
	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);
		}
	}

	//@}

	/** @name Blockers iterators
	 */
	//@{
private:
	/**
	 * @brief Iterator over the blockers adjacent to a vertex
	 */
	typedef Blocker_iterator_around_vertex_internal<
			typename std::multimap<Vertex_handle, Simplex_handle *>::iterator,
			Blocker_handle>
	Complex_blocker_around_vertex_iterator;

	/**
	 * @brief Iterator over (constant) blockers adjacent to a vertex
	 */
	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:
	/**
	 * @brief Returns a range of the blockers of the complex passing through a vertex
	 */
	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);
	}

	/**
	 * @brief Returns a range of the blockers of the complex passing through a vertex
	 */
	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:
	/**
	 * @brief Iterator over the blockers.
	 */
	typedef Blocker_iterator_internal<
			typename std::multimap<Vertex_handle,Simplex_handle *>::iterator,
			Blocker_handle>
	Complex_blocker_iterator;

	/**
	 * @brief Iterator over the (constant) blockers.
	 */
	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:
	/**
	 * @brief Returns a range of the blockers of the complex
	 */
	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);
	}

	/**
	 * @brief Returns a range of the blockers of the complex
	 */
	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);
	}

	//@}

	/////////////////////////////////////////////////////////////////////////////
	/** @name Print and IO methods
	 */
	//@{
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() << "} ";
		stream << std::endl;
		return stream.str();
	}

	std::string blockers_to_string() const {
		std::ostringstream stream;

		for(auto b : const_blocker_range())
			stream<<*b<<std::endl;
		return stream.str();
	}

	//@}
};

/**
 * build a simplicial complex from a collection
 * of top faces.
 * return the total number of simplices
 */
template<typename Complex, typename SimplexHandleIterator>
Complex make_complex_from_top_faces(SimplexHandleIterator simplex_begin,SimplexHandleIterator simplex_end,bool is_flag_complex=false) {
	typedef typename Complex::Simplex_handle Simplex_handle;
	typedef typename Complex::Blocker_handle Blocker_handle;
	typedef typename Complex::Vertex_handle Vertex_handle;
	Complex complex;

	complex.clear();

	std::vector<std::pair<Vertex_handle,Vertex_handle>> edges;
	//first pass to add vertices and edges
	for(auto s_it = simplex_begin; s_it != simplex_end; ++s_it){
		// if meet simplex 9 12 15, need to have at least 16 vertices
		int max_vertex = 0;
		for(auto vh : *s_it )
			max_vertex=(std::max)(max_vertex,(int)vh);
		while( complex.num_vertices() <= max_vertex )
			complex.add_vertex();

		//for all pairs in s, add an edge
		for(auto u_it = s_it->begin(); u_it != s_it->end(); ++u_it)
			for(auto v_it = u_it; ++v_it != s_it->end(); /**/)
				complex.add_edge(*u_it,*v_it);
	}

	if(!is_flag_complex){
		//need to compute blockers

		//store a structure to decide faster if a simplex is in the complex defined by the max faces

		std::vector<std::set<Simplex_handle>> vertex_to_maxfaces(complex.num_vertices());
		for(auto s_it = simplex_begin; s_it != simplex_end; ++s_it)
			vertex_to_maxfaces[s_it->first_vertex()].insert(*s_it);

		// for every maximal face s, it picks its first vertex v and check for all nv \in N(v)
		// if s union nv is in the complex, if not it is a blocker.
		for(auto max_face = simplex_begin; max_face != simplex_end; ++max_face){
			if(max_face->dimension()>0){
				auto first_v = max_face->first_vertex();
				for(auto nv : complex.vertex_range(first_v)){
					if(! max_face->contains(nv) && max_face->first_vertex()<nv){
						//check that all edges in vertices(max_face)\cup nv are here
						//since max_face is a simplex, we only need to check that edges
						// (x nv) where x \in max_face are present
						bool all_edges_here = true;
						for(auto x : *max_face)
							if(!complex.contains_edge(x,nv)){
								all_edges_here = false;
								break;
							}
						if(all_edges_here){ //eg  this->contains(max_face)
							max_face->add_vertex(nv);
							if(vertex_to_maxfaces[first_v].find(*max_face)==vertex_to_maxfaces[first_v].end()){ //xxxx
								// if there exists a blocker included in max_face, we remove it
								// as it is not a minimum missing face
								// the other alternative would be to check to all properfaces
								// are in the complex before adding a blocker but that
								// would be more expensive if there are few blockers
								std::vector<Blocker_handle> blockers_to_remove;
								for(auto b : complex.blocker_range(first_v))
									if(b->contains(*max_face))
										blockers_to_remove.push_back(b);
								for(auto b : blockers_to_remove)
									complex.delete_blocker(b);
								complex.add_blocker(*max_face);
							}
							max_face->remove_vertex(nv);
						}
					}
				}
			}
		}
	}
	return complex;
	//	std::vector<Simplex_handle> simplices;
	//
	//	for (auto top_face = simplex_begin; top_face != simplex_end; ++top_face) {
	//		auto subfaces_topface = subfaces(*top_face);
	//		simplices.insert(simplices.end(),subfaces_topface.begin(),subfaces_topface.end());
	//	}
	//
	//	complex = Complex(simplices.begin(),simplices.end());
	//	return simplices.size();
}

}  // namespace skbl

}  // namespace Gudhi

#endif  // SRC_SKELETON_BLOCKER_INCLUDE_GUDHI_SKELETON_BLOCKER_COMPLEX_H_