//
// This file is part of MorphoGraphX - http://www.MorphoGraphX.org
// Copyright (C) 2012-2015 Richard S. Smith and collaborators.
//
// If you use MorphoGraphX in your work, please cite:
//   http://dx.doi.org/10.7554/eLife.05864
//
// MorphoGraphX is free software, and is licensed under under the terms of the 
// GNU General (GPL) Public License version 2.0, http://www.gnu.org/licenses.
// 
#ifndef VV_GRAPH_HPP
#define VV_GRAPH_HPP

/**
 * \file VVGraph.hpp
 *
 * Contain the definition of the VVGraph template class.
 */

#include <Config.hpp>

#include <CircIterator.hpp>
#include <Edge.hpp>
#include <Forall.hpp>
#include <MemberIterator.hpp>
#include <Tie.hpp>
#include <UnorderedMap.hpp>
#include <Vertex.hpp>

#include <algorithm>
#include <iostream>
#include <iterator>
#include <list>
#include <memory>
#include <utility>
#include <vector>

namespace std 
{
  //using namespace tr1;
}

namespace mgx
{
  /**
   * Contains all the classes related to the graphs.
   *
   * For now, the only graph available is the VVGraph.
   */

  // This was used by the forall macro, now used only for VV
  namespace ForAll 
	{
    struct BaseForwardIter 
		{
      virtual ~BaseForwardIter() {}
    };

    template <typename iterator> struct ForwardIter : public BaseForwardIter 
  	{
      typedef typename std::iterator_traits<iterator>::value_type value_type;
      typedef typename std::iterator_traits<iterator>::reference reference;
  
      ForwardIter(const iterator& fst, const iterator& lst) : it(fst), end(lst), brk(0) {}
  
      const reference value() const
      {
        ++brk;
        return *it;
      }
      bool is_end() const { return (it == end) || brk; }
      mutable iterator it;
      iterator end;
      mutable int brk;
    };
  
    template <typename Pair, typename iterator> 
  	inline ForwardIter<iterator> forwardIter(const Pair& range, const iterator*)
    {
      return ForwardIter<iterator>(range.first, range.second);
    }
  
    template <typename iterator>
    inline const ForwardIter<iterator>* castForwardIter(const BaseForwardIter* base, const iterator*)
    {
      return static_cast<const ForwardIter<iterator>*>(base);
    }
  
    template <typename T> inline T* pointer(const T&) 
    {
      return 0;
    }
  
    template <typename Container>
    std::pair<typename Container::iterator, typename Container::iterator> make_range(Container& cont)
    {
      return std::make_pair(cont.begin(), cont.end());
    }
  
    template <typename Container>
    std::pair<typename Container::const_iterator, 
  					        typename Container::const_iterator> make_range(const Container& cont)
    {
      return std::make_pair(cont.begin(), cont.end());
    }
    
    template <typename Iterator> std::pair<Iterator, Iterator> 
  	make_range(const std::pair<Iterator, Iterator>& cont)
    {
      return cont;
    }
    
    template <typename Container>
    std::pair<typename Container::reverse_iterator, typename Container::reverse_iterator>
    make_reverse_range(Container& cont)
    {
      return std::make_pair(cont.rbegin(), cont.rend());
    }
    
    template <typename Container>
    std::pair<typename Container::const_reverse_iterator, typename Container::const_reverse_iterator>
    make_reverse_range(const Container& cont)
    {
      return std::make_pair(cont.rbegin(), cont.rend());
    }
    
    template <typename Iterator>
    std::pair<std::reverse_iterator<Iterator>, std::reverse_iterator<Iterator> >
    make_reverse_range(const std::pair<Iterator, Iterator>& cont)
    {
      typedef std::reverse_iterator<Iterator> reverse_it;
      return make_pair(reverse_it(cont.second), reverse_it(cont.first));
    }
  }
  
  /**
   * \class _EmptyEdgeContent VVGraph.hpp <VVGraph.hpp>
   *
   * Empty class used as default for edge content
   */
  struct _EmptyEdgeContent {
    _EmptyEdgeContent& operator=(const _EmptyEdgeContent& /*other*/) 
    {
      return *this;
    }
  };
  
  /**
   * \class VVGraph VVGraph.hpp <VVGraph.hpp>
   * Class representing a VV graph.
   *
   * A VV graph is a graph rotation system, that is, a graph where the neighbors
   * of vertices are stored in a circular order. 
   *
   * \c VertexContent is the type used to hold the vertex data.
   *
   * \c EdgeContent is the type used to hold the edge data. If no type is
   * specified for \c EdgeContent, then an empty structure is used.
   *
   */
  template <typename VertexContent, typename EdgeContent = _EmptyEdgeContent>
  class VVGraph 
  {
  public:
    ///\name Smart pointer types
    //@{
    /**
     * Smart pointer on a vertex
     */
    typedef Vertex<VertexContent> vertex_t;
    /**
     * Weak pointer on an edge
     */
    typedef Edge<EdgeContent> edge_t;
    /**
     * Weak pointer on a constant edge
     */
    typedef Edge<const EdgeContent> const_edge_t;
    //@}
  
    /**
     * Type of the value of the graph as standard container.
     */
    typedef vertex_t value_type;
  
    /**
     * Default constructor
     *
     * Creates an empty VV graph
     */
    VVGraph() : modified(false) {}
  
    /**
     * Copy constructor.
     */
    VVGraph(const VVGraph& copy);

    /**
     * Destructor.
     */
    ~VVGraph()
     {
       if(neighborhood.size() > 0) {
         typename neighborhood_t::iterator it;
         for(it = neighborhood.begin(); it != neighborhood.end(); ++it) {
           delete *it;
           *it = 0;
         }
       }
     }
  
    /**
     * Construct from an existing list of vertices
     *
     * \param checkUnique If false, there must be not be duplicated vertex in the sequence
     */
    template <typename Iterator> VVGraph(Iterator begin, Iterator end, bool checkUnique = true)
    {
      insert(begin, end, checkUnique);
    }
  
    /**
     * Construct from an existing list of vertices
     *
     * \param checkUnique If false, there must be not be duplicated vertex in the sequence
     */
    template <typename Container> VVGraph(const Container& c, bool checkUnique = true)
    {
      insert(ForAll::make_range(c), checkUnique);
    }

  protected:
    struct single_neighborhood_t;
  
    /**
     * Structure maintaining the data for a single neighbor
     */
    struct neighbor_t : public EdgeContent {
      /**
       * Type of the list of outgoing neighbors
       */
      //RSS typedef std::list<neighbor_t> edge_list_t;
      // swaped list for vector, does not seem to make much difference performance-wise, i
      // and should be smaller
      typedef std::vector<neighbor_t> edge_list_t;
  
      /**
       * Constructor
       */
      neighbor_t(const vertex_t& tgt, const EdgeContent& e) : EdgeContent(e), target(tgt) {}
  
      neighbor_t(const neighbor_t& copy) : EdgeContent(copy), target(copy.target) {}
  
      void clear_edge() 
      {
        static_cast<EdgeContent&>(*this) = EdgeContent();
      }
  
      ~neighbor_t() {}
  
      neighbor_t& operator=(const neighbor_t& copy)
      {
        target = copy.target;
        static_cast<EdgeContent&>(*this) = static_cast<const EdgeContent&>(copy);
        return *this;
      }
  
      /**
       * Target of the edge
       *
       * i.e. neighbor itself
       */
      vertex_t target;
  
      /**
       * Equality of two neighbors/edge
       *
       * A single neighbor is equal to another if the vertex and the content are
       * the same
       */
      bool operator==(const neighbor_t& other) const
      {
        return (target == other.target)
             && (static_cast<const EdgeContent&>(*this) == static_cast<const EdgeContent&>(other));
      }
    };
  
    /**
     * Type of the list of outgoing neighbors
     */
    typedef typename neighbor_t::edge_list_t edge_list_t;
  
    /**
     * Type of the neighborhood of a vertex
     */
    struct single_neighborhood_t 
    {
  
      single_neighborhood_t() {}
  
      single_neighborhood_t(const single_neighborhood_t& copy) : edges(copy.edges) {}
  
      /**
       * List of the outgoing edges
       */
      edge_list_t edges;
  
      single_neighborhood_t& operator=(const single_neighborhood_t& other)
      {
        return *this;
      }
  
      /**
       * Equality of two neighborhood
       *
       * Only the outgoing edges are compared!
       */
      bool operator==(const single_neighborhood_t& other) const 
      {
        return edges == other.edges;
      }
    };
  
  public:
    struct NeighborhoodPair 
    {
      NeighborhoodPair(const vertex_t& v) : vertex(v) {}
  
      NeighborhoodPair(const NeighborhoodPair& copy)
        : vertex(copy.vertex), single_neighborhood(copy.single_neighborhood) {}
  
      NeighborhoodPair* copy() const
      {
        NeighborhoodPair* np = new NeighborhoodPair(*this);
        return np;
      }
  
      vertex_t vertex;
      single_neighborhood_t single_neighborhood;
    };
  
    /**
     * Type of the list of vertexes, together with the neighborhood
     */
    typedef std::vector<NeighborhoodPair*> neighborhood_t;
  
    /**
     * Type of the fast-search map
     */
    typedef std::unordered_map<vertex_t, NeighborhoodPair*> lookup_t;
  
    /**
     * Shortcut for the value_type of the neighborhood
     */
    typedef typename neighborhood_t::value_type neighborhood_value_type;
  
    /**
     * Type of a size
     */
    typedef typename neighborhood_t::size_type size_type;
  
  ///\name Iterators
  //@{
  /**
   * Iterator on the vertexes
   */
    typedef SelectMemberPointerIterator<typename neighborhood_t::iterator, 
                  const vertex_t, &NeighborhoodPair::vertex> iterator;
    /**
     * Constant iterator on the vertexes
     */
    typedef SelectMemberPointerIterator<typename neighborhood_t::const_iterator, 
                    const vertex_t, &NeighborhoodPair::vertex> const_iterator;
    /**
     * Iterator on the neighbors of a vertex
     */
    typedef SelectMemberIterator<typename edge_list_t::iterator, vertex_t, 
                                                &neighbor_t::target> neighbor_iterator;
    /**
     * Range of the neighbors of a vertex. The \c first element of the pair is an
     * iterator on the beginning of the range, the \c second element is the
     * past-the-end iterator.
     */
    typedef std::pair<neighbor_iterator, neighbor_iterator> neighbor_iterator_pair;
  
    /**
     * Iterator used to iterate over the neighbors, but specifying the
     * starting point
     */
    typedef CircIterator<neighbor_iterator> circ_neighbor_iterator;
  
    /**
     * Range of circular iterators
     */
    typedef std::pair<circ_neighbor_iterator, circ_neighbor_iterator> circ_neighbor_iterator_pair;
  
    /**
     * Constant iterator on the neighbors of a vertex
     */
    typedef SelectMemberIterator<typename edge_list_t::const_iterator, 
                                                  const vertex_t, &neighbor_t::target>
      const_neighbor_iterator;
    /**
     * Constant range of the neighbors of a vertex. The \c first element of the
     * pair is a constant iterator on the beginning of the range, the \c second
     * element is the past-the-end constant iterator.
     */
    typedef std::pair<const_neighbor_iterator, const_neighbor_iterator> 
                                                         const_neighbor_iterator_pair;
  
    /**
     * Iterator used to iterate over the neighbors, but spcifying the
     * starting point
     */
    typedef CircIterator<const_neighbor_iterator> const_circ_neighbor_iterator;
  
    /**
     * Range of circular iterators
     */
    typedef std::pair<const_circ_neighbor_iterator, const_circ_neighbor_iterator> 
                                                      const_circ_neighbor_iterator_pair;
  
    //@}
  
    ///\name Vertex set editing methods
    //@{
  
    /**
     * Remove a vertex from the graph.
     *
     * \param v vertex to erase
     *
     * \returns 1 if the vertex was erased, 0 else.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Add vertices to S
     * vertex v = S.any();
     * S.erase(v); // Remove a vertex of S
     * \endcode
     */
    size_type erase(const vertex_t& v);
  
    /**
     * Insert a new vertex in the graph.
     *
     * \returns An true if successful, false if vertex was already there
     *
     * Example:
     * \code
     * vvGraph S;
     * vertex v; // Create a vertex
     * S.insert(v); // Insert it in the graph
     * \endcode
     */
    bool insert(const vertex_t& v, bool checkUnique = true);
  
    //@}
  
    ///\name Vertex set lookup methods
    //@{
  
    /**
     * Return a vertex from the graph.
     *
     * \returns A vertex of the graph, or a null vertex if the graph is empty.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Add vertices in S
     * const vertex& v = S.any(); // Get a vertex in S
     * \endcode
     */
    const vertex_t& any() const;
  
    /**
     * Return the element of index \c idx.
     *
     * \note This is a convenience function whose complexity if O(idx)
     *
     * \returns The element \c idx position after the first one (using
     * iterators)
     *
     * Example:
     * \code
     * vvGraph S;
     * // Fill in S
     * size_t i = 0, selected = 0;
     * forall(const vertex& v, S)
     * {
     *   if(condition(v))
     *   {
     *     selected = i;
     *     break;
     *   }
     *   ++i;
     * }
     * // ...
     * const vertex& v = S[selected]; // Work with the vertex previously selected
     * \endcode
     * This construct is useful is case you need to refer to your vertices by
     * 32bits numbers. It is true if you use selection with OpenGL for
     * example.
     */
    const vertex_t& operator[](size_type idx) const;
  
    /**
     * Test if \c v is in the graph.
     *
     * \returns \c true is \c v is in the graph.
     *
     * Example:
     * \code
     * vvGraph S;
     * vertex v;
     * S.insert(v);
     * assert(S.contains(v));
     * \endcode
     */
    bool contains(const vertex_t& v) const;
  
    /**
     * Get a reference on the vertex in the graph
     *
     * If the vertex does not exist in the graph, the function returns a null
     * vertex.
     *
     * Example:
     * \code
     * vvGraph S1, S2;
     * // Initialize S1 and S2 with the same vertices but different neighborhood
     * forall(const vertex& v1, S1)
     * {
     *   const vertex& v2 = S2.reference(v1);
     *   // Use v2 in S2 and v1 in v1 for speed optimizatioN
     * }
     * \endcode
     */
    const vertex_t& reference(vertex_t v) const
    {
      const NeighborhoodPair* found = this->findVertex(v);
      if(found)
        return found->vertex;
      return vertex_t::null;     // vertex_t(0);
    }
  
    /**
     * Return the number of vertexes on the graph
     */
    size_type size() const 
    {
      return neighborhood.size();
    }
  
    /**
     * Test if there is any vertex in the graph
     */
    bool empty() const 
   {
      return neighborhood.empty();
    }

    /**
     * Test if vertex list has been modified
     */
    bool isModified() const 
    {
      return modified;
    }

    /**
     * Set modified status
     */
    void setModified(bool _modified) 
    {
       modified = _modified;
    }
  
  
    //@}
  
    ///\name Neighborhood lookup methods
    //@{
  
    /**
     * Return a vertex in the neighborhood of \c v.
     *
     * \returns A vertex in the neighborhood of \c v, or a null vertex if \c v is
     * not in the graph or \c v has no neighborhood.
     *
     * Example:
     * \code
     * vvGraph S;
     * forall(const vertex& v, S)
     * {
     *   const vertex& n = S.anyIn(v);
     *   if(n)
     *   {
     *     edge e = S.edge(v,n); // e exist !
     *   }
     * }
     * \endcode
     */
    const vertex_t& anyIn(const vertex_t& v) const;
  
    /**
     * Returns the number of neighbors of \c v.
     *
     * If \c v is not in the graph, the behavior is undefined.
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& v, S)
     * {
     *   size_t nb_neighbors = S.valence(v); // Number of neighbors of v in S
     *   // ...
     * }
     * \endcode
     */
    size_type valence(const vertex_t& v) const;
  
    /**
     * Test if a vertex has a neighbor.
     *
     * Returns false if it does. If a vertex is not in the graph, it is
     * considered as having no neighbor. In the same way, the null vertex has
     * no neighbors.
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(cnst vertex& v, S)
     * {
     *   if(!S.empty(v))
     *   {
     *      const vertex& n = S.anyIn(v);
     *      // Working with n a neighbor of v
     *   }
     * }
     * \endcode
     */
    bool empty(const vertex_t& v) const;
  
    /**
     * Returns the nth vertex after \c neighbor in the neighborhood of \c v.
     *
     * \returns The vertex found or a null vertex if there is any problem.
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     const vertex& m, S.nextTo(v);
     *     // m is the vertex after n in v
     *   }
     * }
     * \endcode
     */
    const vertex_t& nextTo(const vertex_t& v, const vertex_t& neighbor, unsigned int n = 1) const;
  
    /**
     * Returns the nth vertex before ref in the neighborhood of v.
     *
     * If ref is not in the neighborhood of v, return a null vertex.
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     const vertex& m, S.prevTo(v);
     *     // m is the vertex before n in v
     *   }
     * }
     * \endcode
     */
    const vertex_t& prevTo(const vertex_t& v, const vertex_t& ref, unsigned int n = 1) const;
  
    /**
     * Returns the edge from \c src to \c tgt.
     *
     * If the edge does not exists, returns a null edge.
     *
     * Example:
     * \code
     * struct VertexContent {};
     * struct EdgeContent { int a; }
     * typedef VVGraph<VertexContent,EdgeContent> vvGraph;
     * typedef vvGraph::vertex_t vertex;
     * typedef vvGraph::edge_t edge;
     * // ...
     *
     * vvGraph S;
     * // initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     edge e = S.edge(v,n);
     *     e->a = 10;
     *   }
     * }
     * \endcode
     */
    edge_t edge(const vertex_t& src, const vertex_t& tgt);
  
    /**
     * Returns the edge from \c src to \c tgt.
     *
     * If the edge does not exists, returns a null edge.
     *
     * Example:
     * \code
     * struct VertexContent {};
     * struct EdgeContent { int a; }
     * typedef VVGraph<VertexContent,EdgeContent> vvGraph;
     * typedef vvGraph::vertex_t vertex;
     * typedef vvGraph::edge_t edge;
     * typedef vvGraph::const_edge_t const_edge;
     * // ...
     *
     * void fct(const vvGraph& S)
     * {
     *   forall(const vertex& v, S)
     *   {
     *     forall(const vertex& n, S.neighbors(v))
     *     {
     *       const_edge e = S.edge(v,n);
     *       cout << e->a << endl; // Here e->a cannot be modified
     *     }
     *   }
     * }
     * \endcode
     */
    const_edge_t edge(const vertex_t& src, const vertex_t& tgt) const;
  
    /**
     * Return the constant range of neighbors of \c v
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     // n is a neighbor of v
     *   }
     * }
     * \endcode
     */
    const_neighbor_iterator_pair neighbors(const vertex_t& v) const;
    /**
     * Return the range of neighbors of \c v
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     // n is a neighbor of v
     *   }
     * }
     * \endcode
     */
    neighbor_iterator_pair neighbors(const vertex_t& v);
  
    /**
     * Return the range of neighbors of \c v, starting at \c n
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& m, S.neighbors(v,n))
     * {
     *   // m is a neighbor of v
     *   // on the first loop, m == n
     * }
     * \endcode
     */
    circ_neighbor_iterator_pair neighbors(const vertex_t& v, const vertex_t& n);
  
    /**
     * Return the range of neighbors of \c v, starting at \c n
     *
     * Example:
     * \code
     * vvGraph S;
     * // initialize S
     * forall(const vertex& m, S.neighbors(v,n))
     * {
     *   // m is a neighbor of v
     *   // on the first loop, m == n
     * }
     * \endcode
     */
    const_circ_neighbor_iterator_pair neighbors(const vertex_t& v, const vertex_t& n) const;
  
    /**
     * Return the source vertex of the edge
     *
     * Example:
     * \code
     * void fct(edge e, const vvGraph& S)
     * {
     *   const vertex& src = S.source(e);
     *   const vertex& tgt = S.target(e);
     *   // Work with the vertices
     * }
     * \endcode
     */
    const vertex_t& source(const edge_t& edge) const
    {
      const NeighborhoodPair* found = this->findVertex(vertex_t(edge.source()));
      if(found)
        return found->vertex;
      return vertex_t::null;     // vertex_t(0);
    }
  
    /**
     * Return the target vertex of the edge
     *
     * Example:
     * \code
     * void fct(edge e, const vvGraph& S)
     * {
     *   const vertex& src = S.source(e);
     *   const vertex& tgt = S.target(e);
     *   // Work with the vertices
     * }
     * \endcode
     */
    const vertex_t& target(const edge_t& edge) const
    {
      const NeighborhoodPair* found = this->findVertex(vertex_t(edge.target()));
      if(found)
        return found->vertex;
      return vertex_t::null;     // vertex_t(0);
    }
  
    //@}
  
    ///\name Neighborhood editing methods
    //@{
    /**
     * Erase the edge from \c src to \c tgt.
     *
     * \returns True if successful.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     if(condition)
     *     {
     *       S.eraseEdge(v,n);
     *       break; // The iteration must not continue after the neighborhood has been altered
     *     }
     *   }
     * }
     * \endcode
     */
    size_type eraseEdge(const vertex_t& src, const vertex_t& tgt);
  
    /**
     * Replace a vertex by another in a neighborhood.
     *
     * \param[in] v Vertex whose neighborhood is changed.
     * \param[in] neighbor Vertex to replace
     * \param[in] new_neighbor Vertex replacing \c neighbor
     *
     * If \c new_neighbor is already in the neighborhood of \c v, then the
     * operation fails and nothing is changed.
     *
     * \returns The edge between \c v and \c new_neighbor or a null edge if
     * anything goes wrong.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     if(condition)
     *     {
     *       vertex n1;
     *       S.replace(v,n,n1); // Now n1 is where n was in v
     *       break; // The iteration must not continue after the neighborhood has been altered
     *     }
     *   }
     * }
     * \endcode
     */
    edge_t replace(const vertex_t& v, const vertex_t& neighbor, const vertex_t& new_neighbor);
  
    /**
     * Insert neighbor in the neighborhood of v after reference.
     *
     * If \c new_neighbor is already in the neighborhood of \c v, the insertion
     * fails.
     *
     * \returns The just created edge if everything succeed, or a null edge.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     if(condition)
     *     {
     *       vertex n1;
     *       S.spliceAfter(v,n,n1); // Now n1 is after n in v
     *       break; // The iteration must not continue after the neighborhood has been altered
     *     }
     *   }
     * }
     * \endcode
     */
    edge_t spliceAfter(const vertex_t& v, const vertex_t& neighbor, const vertex_t& new_neighbor);
  
    /**
     * Insert neighbor in the neighborhood of v before reference.
     *
     * If \c new_neighbor is already in the neighborhood of \c v, the , checkUnique = trueinsertion
     * fails.
     *
     * \returns The just created edge if everything succeed, or a null edge.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * forall(const vertex& v, S)
     * {
     *   forall(const vertex& n, S.neighbors(v))
     *   {
     *     if(condition)
     *     {
     *       vertex n1;
     *       S.spliceBefore(v,n,n1); // Now n1 is before n in v
     *       break; // The iteration must not continue after the neighborhood has been altered
     *     }
     *   }
     * }
     * \endcode
     */
    edge_t spliceBefore(const vertex_t& v, const vertex_t& neighbor, const vertex_t& new_neighbor);
  
    /**
     * Insert a new edge in the graph, without ordering.
     *
     * If \c new_neighbor is already in the neighborhood of \c v, the insertion
     * fails.
     *
     * \returns The just created edge if everything succeed, or a null edge.
     *
     * Example:
     * \code
     * vvGraph S;
     * vertex v1, v2;
     * S.insert(v1);
     * S.insert(v2);
     * S.insertEdge(v1, v2); // Insert the edge between v1 and v2
     * S.insertEdge(v2, v1); // Insert the edge between v2 and v1
     * \endcode
     */
    edge_t insertEdge(const vertex_t& src, const vertex_t& tgt);
  
    /**
     * Clear the neighborhood of a vertex.
     *
     * All edges, to and from \c v will be erased.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * forall(const vertex& v, S)
     * {
     *   if(condition)
     *   {
     *     S.eraseAllEdges(v);
     *     assert(S.empty(v)); // No outgoing edges
     *     assert(S.iEmpty(v)); // No incoming edges
     *   }
     * }
     * \endcode
     */
    bool eraseAllEdges(const vertex_t& v);
  
    /**
     * Clear the outgoing edges of a vertex
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * forall(const vertex& v, S)
     * {
     *   if(condition)
     *   {
     *     S.clear(v);
     *     assert(S.empty(v));
     *   }
     * }
     * \endcode
     */
    bool clear(const vertex_t& v);

    /**
     * Test for unique triangle
     *
     * Only true for one of (v,n,m) (n,m,v) or (m,v,n). Only for use with
     * triangle meshes.
     */
    bool uniqueTri(const vertex_t& v, const vertex_t& n, const vertex_t& m) const;

    /**
     * Test for unique line
     *
     * Only true for one of (v,n) or (n,v). Only use in graphs that always have
     * edges in both directions.
     */
    bool uniqueLine(const vertex_t& v, const vertex_t& n) const;
    //@}
  
    ///\name Whole structure methods
    //@{
    /**
     * Copy the graph into another graph
     */
    VVGraph& operator=(const VVGraph& other);
  
    /**
     * Reverse the order of all the neighbors
     */
    VVGraph& reverse();
  
    /**
     * Swap the content of both graphs
     */
    void swap(VVGraph& other);
  
    /**
     * Test equality for the graphs.
     *
     * Two graphs are equal if they shared the sames vertices and their
     * neighborhood are equivalent.
     *
     * \note Current implementation do not consider invariance of
     * neighborhood by rotation.
     */
    bool operator==(const VVGraph& other) const;
  
    /**
     * Clear the graph.
     */
    void clear();
  
    /**
     * Clear all the edges in the graph
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * vvGraph::size_type s1 = S.size(); // Save the number of vertices
     * S.eraseAllEdges();
     * assert(s1 == S.size()); // The number of vertices didn't change
     * forall(const vertex& v, S)
     * {
     *   assert(S.empty(v)); // No neighbor
     * }
     * \endcode
     */
    void eraseAllEdges();
  
    /**
     * Extract the subgraph containing the \p vertexes
     *
     * The subgraph is the set of vertexes and the edges whose source and target
     * are in the extracted vertexes.
     *
     * Example:
     * \code
     * vvGraph S;
     * // Initialize S
     * std::vector<vertex> to_keep;
     * // Fill in to_keep with a subset of the vertices of S
     * vvGraph S1 = S.subgraph(to_keep);
     * \endcode
     */
    template <typename VertexContainer> VVGraph subgraph(const VertexContainer& vertexes) const;
    //@}
  
    ///\name STL compatibility methods
    //@{
  
    /**
     * Erase element at position \p pos
     */
    void erase(iterator pos);
  
    /**
     * Erase the edge pointed to by the iterator.
     */
    size_type eraseEdge(const vertex_t& v, neighbor_iterator pos);
  
    /**
     * Erase the edge pointed to by the iterator.
     */
    size_type eraseEdge(const vertex_t& v, circ_neighbor_iterator pos);
  
    /**
     * Insert all the vertexes from first to last-1 in the graph.
     */
    template <typename Iterator> 
    void insert(Iterator first, Iterator last, bool checkUnique = true);
  
    /**
     * Insert all the vertexes from range.first to range.last-1 in the graph.
     */
    template <typename Iterator> 
    void insert(const std::pair<Iterator, Iterator>& range, bool checkUnique = true);
  
    /**
     * Returns an iterator on the beginning of the set of vertexes of the graph.
     */
    iterator begin() {
      return neighborhood.begin();
    }
    /**
     * Returns an iterator on the end of the set of vertexes of the graph.
     */
    iterator end() {
      return neighborhood.end();
    }
  
    /**
     * Returns an iterator on the beginning of the set of vertexes of the graph.
     */
    const_iterator begin() const {
      return neighborhood.begin();
    }
    /**
     * Returns an iterator on the end of the set of vertexes of the graph.
     */
    const_iterator end() const {
      return neighborhood.end();
    }
  
    /**
     * Count the number of vertexes \p v in the graph (can be 0 or 1 only)
     */
    size_type count(const vertex_t& v) const;
  
    /**
     * Find the vertex \c n in the neighborhood of \c v.
     *
     * \returns
     *  - An iterator on \c n in the neighborhood of \c v, if found
     *  - An iterator on the end of the neighborhood if \c n is not found in \c
     *  v
     *  - Otherwise, the result is undefined
     */
    neighbor_iterator findIn(const vertex_t& v, const vertex_t& n);
    /**
     * Find the vertex \c n in the neighborhood of \c v.
     *
     * \returns
     *  - An iterator on \c n in the neighborhood of \c v, if found
     *  - An iterator on the end of the neighborhood if \c n is not found in \c
     *  v
     *  - Otherwise, the result is undefined
     */
    const_neighbor_iterator findIn(const vertex_t& v, const vertex_t& n) const;
  
    /**
     * Clear the neighborhood of a vertex.
     *
     * All edges, to and from \c *it will be erased.
     */
    void eraseAllEdges(iterator it);
  
    /**
     * Clear the outgoing edges of a vertex.
     */
    void clear(iterator it);
  
    /**
     * Reserve space for more vertices
     */
    void reserve(size_type s) {
      neighborhood.reserve(s);
    }
  
    //@}
  
  protected:
    void eraseAllEdges(NeighborhoodPair* it);
    /**
     * Type of the result of the search for a vertex in a neighborhood
     *
     * This data structure is made such that all edit or lookup operations can
     * be done from that data structure.
     */
    template <class Neighborhood, class Iterator> struct search_result_t 
    {
      /**
       * Default constructor
       *
       * By default, the search is unsuccessful
       */
      search_result_t() : found(false), neighborhood(0) {}
  
      /**
       * Successful constructor
       */
      search_result_t(Iterator i, Neighborhood* n, bool ok = true)
        : found(ok) , it(i), neighborhood(n) {}
  
      /**
       * Copy constructor
       */
      search_result_t(const search_result_t& copy)
        : found(copy.found), it(copy.it), neighborhood(copy.neighborhood) {}
  
      /**
       * Convert the object to true if the search is successful
       */
      operator bool() 
      {
        return found;
      }
  
      /**
       * True if the search was completely successful
       */
      bool found;
      /**
       * Iterator pointing in the edge list
       */
      Iterator it;
      /**
       * Neighborhood containing the element
       */
      Neighborhood* neighborhood;
    };
  
    /**
     * Iterator on the outgoing edges
     */
    typedef typename edge_list_t::iterator int_neighbor_iterator;
    /**
     * Constant iterator on the outgoing edges
     */
    typedef typename edge_list_t::const_iterator int_const_neighbor_iterator;
  
    /**
     * Result of a search in the neighborhood
     */
    typedef search_result_t<single_neighborhood_t, int_neighbor_iterator> neighbor_found_t;
    /**
     * Constant result of a search in the neighborhood
     */
    typedef search_result_t<const single_neighborhood_t, int_const_neighbor_iterator> 
                                                                      const_neighbor_found_t;
  
    /**
     * Find a vertex in the graph and returns the iterator on it.
     */
    NeighborhoodPair* findVertex(const vertex_t& v);
  
    /**
     * Constant version of findVertex(const vertex_t&)
     */
    const NeighborhoodPair* findVertex(const vertex_t& v) const;
  
    /**
     * Find the iterator for a given neighborhood pair
     */
    typename neighborhood_t::iterator findIterator(const NeighborhoodPair* np);
    typename neighborhood_t::const_iterator findIterator(const NeighborhoodPair* np) const;
  
    /**
     * Find a vertex neighbor in the neighborhood of v.
     *
     * If 'v' is in the graph and 'neighbor' is in its neighborhood, the
     * result is convertible to true, 'neighborhood' points toward the
     * neighborhood structure and 'it' points toward the neighbor.
     *
     * If 'v' is in the graph, but 'neighbor' is not in the neighborhood, the
     * result is convertible to false, 'neighborhood' points toward the
     * neighborhood structure and 'it' is neighborhood->edges.end().
     *
     * If v is not in the graph, the result is convertible to false and the
     * address of the stored neighborhood is 0.
     */
    neighbor_found_t findInVertex(const vertex_t& v, const vertex_t& neighbor);
  
    /**
     * Constant version of findInVertex(const vertex_t&, const vertex_t&)
     */
    const_neighbor_found_t findInVertex(const vertex_t& v, const vertex_t& neighbor) const;
  
    /**
     * Main data structure containing everything
     */
    neighborhood_t neighborhood;
  
    /**
     * Hash map to optimise look-up
     */
    lookup_t lookup;
  
    /**
     * If true, the vertex list has been modified
     */
    bool modified;
  };
  
  // Clear the VVGraph
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::clear()
  {
    if(neighborhood.size() == 0)
      return;
    typename neighborhood_t::iterator it;
    for(it = neighborhood.begin(); it != neighborhood.end(); ++it) {
      delete *it;
      *it = 0;
    }
    neighborhood.clear();
    lookup.clear();
    modified = true;
  }
  
  // Erase vertex at iterator, RSS can we check it is valid?
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::erase(iterator it)
  {
    eraseAllEdges(it);
    lookup.erase(*it);
    delete *it.base();
    *it.base() = neighborhood.back();
    neighborhood.pop_back();
    modified = true;
  }
  
  // Erase a vertex, return 1 is successful
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::size_type
  VVGraph<VertexContent, EdgeContent>::erase(const vertex_t& v)
  {
    if(eraseAllEdges(v)) {
      NeighborhoodPair* np = lookup[v];
      lookup.erase(v);
      typename neighborhood_t::iterator it = this->findIterator(np);
      delete np;
      *it = neighborhood.back();
      neighborhood.pop_back();
      modified = true;
      return 1;
    }
    return 0;
  }

  // Erase all edges of a vertex
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::eraseAllEdges()
  {
    if(neighborhood.size() == 0)
      return;
    for(typename neighborhood_t::iterator it = neighborhood.begin(); it != neighborhood.end(); ++it)
      it->single_neighborhood.edges.clear();
    modified = true;
  }
 
  // Erase all edges of a vertex, return true if sucessful
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::eraseAllEdges(const vertex_t& v)
  {
    NeighborhoodPair* it_found = this->findVertex(v);
    if(it_found) {
      eraseAllEdges(it_found);
      return true;
    }
    return false;
  }
  
  // Erase all edges of a neighborhood pair
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::eraseAllEdges(NeighborhoodPair* np)
  {
    const vertex_t& v = np->vertex;

    NeighborhoodPair* it_found = this->findVertex(v);
    edge_list_t& lst = it_found->single_neighborhood.edges;
    for(typename edge_list_t::iterator it = lst.begin(); it != lst.end(); ++it)
      eraseEdge(it->target, v);

    np->single_neighborhood.edges.clear();
    modified = true;
  }
  
  // Erase alls edges of a vertex at iterator
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::eraseAllEdges(iterator it_found)
  {
    eraseAllEdges(*it_found.base());
  }
  
  // Clear a vertex, erase all edges
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::clear(const vertex_t& v)
  {
    typename neighborhood_t::iterator it_found = this->findVertex(v);
    if(it_found != neighborhood.end()) {
      clear(it_found);
      return true;
    }
    return false;
  }
  
  // Clear a vertex?
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::clear(iterator it)
  {
    const vertex_t& v = *it;
    it.base()->single_neighborhood.edges.clear();
    modified = true;
  }
  
  // Return 1 if vertex exists
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::size_type
  VVGraph<VertexContent, EdgeContent>::count(const vertex_t& v) const
  {
    return lookup.count(v);
  }
  
  // Inset a vertex in the graph and return an iterator to it
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::insert(const vertex_t& v, bool checkUnique)
  {
    typename lookup_t::iterator it_found = lookup.end();
    if(checkUnique)
      it_found = lookup.find(v);
    if(!checkUnique or it_found == lookup.end()) {
      NeighborhoodPair* np = new NeighborhoodPair(v);
      neighborhood.push_back(np);
      lookup[v] = np;
      iterator it = neighborhood.end();
      --it;
      modified = true;
      return true;
    }
    return false;
  }
  
  // Insert a range of vertices into the graph 
  template <typename VertexContent, typename EdgeContent>
  template <typename Iterator>
  void VVGraph<VertexContent, EdgeContent>::insert(const std::pair<Iterator, Iterator>& range, bool checkUnique)
  {
    insert(range.first, range.second, checkUnique);
  }
  
  // Insert a range of vertices into the graph 
  template <typename VertexContent, typename EdgeContent>
  template <typename Iterator>
  void VVGraph<VertexContent, EdgeContent>::insert(Iterator first, Iterator last, bool checkUnique)
  {
    if(checkUnique) {
      for(Iterator it = first; it != last; ++it)
        insert(*it);
    } else {
      for(Iterator it = first; it != last; ++it) {
        NeighborhoodPair* np = new NeighborhoodPair(*it);
        neighborhood.push_back(np);
        lookup[*it] = np;
      }
    }
    modified = true;
  }
  
  // Return any vertex in the graph
  template <typename VertexContent, typename EdgeContent>
  const typename VVGraph<VertexContent, EdgeContent>::vertex_t&
  VVGraph<VertexContent, EdgeContent>::any() const
  {
    if(neighborhood.empty())
      return vertex_t::null;
    return *begin();
  }
  
  // Return vertex at []
  template <typename VertexContent, typename EdgeContent>
  const typename VVGraph<VertexContent, EdgeContent>::vertex_t& VVGraph<VertexContent, EdgeContent>
  ::operator[](size_type idx) const
  {
    typename neighborhood_t::const_iterator it = neighborhood.begin();
    std::advance(it, idx);
    return (*it)->vertex;
  }
  
  // Return any neighbor
  template <typename VertexContent, typename EdgeContent>
  const typename VVGraph<VertexContent, EdgeContent>::vertex_t&
  VVGraph<VertexContent, EdgeContent>::anyIn(const vertex_t &v) const
  {
    if(v.isNull())
      return vertex_t::null;     // vertex_t(0);
    const NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found || it_found->single_neighborhood.edges.empty())
      return vertex_t::null;     // vertex_t(0);
    const edge_list_t& lst = it_found->single_neighborhood.edges;
    return lst.front().target;
  }
  
  // Return the valence
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::size_type
  VVGraph<VertexContent, EdgeContent>::valence(const vertex_t& v) const
  {
    if(v.isNull())
      return 0;
    const NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return 0;
    return it_found->single_neighborhood.edges.size();
  }
  
  // Is graph empty
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::empty(const vertex_t& v) const
  {
    if(v.isNull())
      return true;
    const NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return true;
    return it_found->single_neighborhood.edges.empty();
  }
  
  // Find a vertex
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::neighborhood_t::iterator
  VVGraph<VertexContent, EdgeContent>::findIterator(const NeighborhoodPair* np)
  {
    return std::find(neighborhood.begin(), neighborhood.end(), np);
  }
  
  // Find a vertex
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::neighborhood_t::const_iterator
  VVGraph<VertexContent, EdgeContent>::findIterator(const NeighborhoodPair* np) const
  {
    return std::find(neighborhood.begin(), neighborhood.end(), np);
  }
  
  // Find a vertex
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::NeighborhoodPair*
  VVGraph<VertexContent, EdgeContent>::findVertex(const vertex_t &v)
  {
    typename lookup_t::const_iterator it_found = lookup.find(v);
    if(it_found != lookup.end())
      return it_found->second;
    return 0;
  }
  
  // Find a vertex
  template <typename VertexContent, typename EdgeContent>
  const typename VVGraph<VertexContent, EdgeContent>::NeighborhoodPair*
  VVGraph<VertexContent, EdgeContent>::findVertex(const vertex_t &v) const
  {
    typename lookup_t::const_iterator it_found = lookup.find(v);
    if(it_found != lookup.end())
      return it_found->second;
    return 0;
  }
  
  // Find a vertex in a neighborhood
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::neighbor_found_t
  VVGraph<VertexContent, EdgeContent>::findInVertex(const vertex_t& v, const vertex_t& n)
  {
    if(v.isNull() || n.isNull() || v == n)
      return neighbor_found_t();
    NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return neighbor_found_t();
    edge_list_t& lst = it_found->single_neighborhood.edges;
    single_neighborhood_t* neighborhood = &it_found->single_neighborhood;
    for(typename edge_list_t::iterator it = lst.begin(); it != lst.end(); ++it) {
      if(it->target == n) {
        return neighbor_found_t(it, neighborhood);
      }
    }
    return neighbor_found_t(lst.end(), neighborhood, false);
  }
  
  // Find a vertex in a neighborhood
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::const_neighbor_found_t
  VVGraph<VertexContent, EdgeContent>::findInVertex(const vertex_t& v, const vertex_t& n) const
  {
    if(v.isNull() || n.isNull() || v == n)
      return const_neighbor_found_t();
    const NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return const_neighbor_found_t();
    const edge_list_t& lst = it_found->single_neighborhood.edges;
    const single_neighborhood_t* neighborhood = &it_found->single_neighborhood;
    for(typename edge_list_t::const_iterator it = lst.begin(); it != lst.end(); ++it) {
      if(it->target == n) {
        return const_neighbor_found_t(it, neighborhood);
      }
    }
    return const_neighbor_found_t(lst.end(), neighborhood, false);
  }
  
  // Find a vertex in a neighborhood
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::neighbor_iterator
  VVGraph<VertexContent, EdgeContent>::findIn(const vertex_t& v, const vertex_t& n)
  {
    neighbor_found_t found = this->findInVertex(v, n);
    if(found.neighborhood)
      return neighbor_iterator(found.it);
    return neighbor_iterator();
  }
  
  // Find a vertex in a neighborhood
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::const_neighbor_iterator
  VVGraph<VertexContent, EdgeContent>::findIn(const vertex_t& v, const vertex_t& n) const
  {
    const_neighbor_found_t found = this->findInVertex(v, n);
    if(found.neighborhood)
      return const_neighbor_iterator(found.it);
    return const_neighbor_iterator(found.it);
  }
  
  // Erase an edge
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::size_type
  VVGraph<VertexContent, EdgeContent>::eraseEdge(const vertex_t& v, neighbor_iterator pos)
  {
    return eraseEdge(v, *pos);
  }
  
  // Erase an edge
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::size_type
  VVGraph<VertexContent, EdgeContent>::eraseEdge(const vertex_t& v, circ_neighbor_iterator pos)
  {
    return eraseEdge(v, pos.base());
  }
  
  // Erase an edge
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::size_type
  VVGraph<VertexContent, EdgeContent>::eraseEdge(const vertex_t& v, const vertex_t& n)
  {
    neighbor_found_t found = this->findInVertex(v, n);
    if(!found)
      return 0;
    found.neighborhood->edges.erase(found.it);
    modified = true;
    return 1;
  }
  
  // See if graph contains a vertex
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::contains(const vertex_t& v) const
  {
    return this->findVertex(v);
  }
  
  // Replace a vertex in a neighborhood with another
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::edge_t
  VVGraph<VertexContent, EdgeContent>::replace(const vertex_t& v, const vertex_t& neighbor,
                                                      const vertex_t& new_neighbor)
  {
    if(new_neighbor.isNull() or (v == new_neighbor) or (neighbor == new_neighbor))
      return edge_t();
    neighbor_found_t found = this->findInVertex(v, neighbor);
    if(!found)
      return edge_t();
    //NeighborhoodPair* n_found = this->findVertex(neighbor);

    found.it->target = new_neighbor;
    found.it->clear_edge();
    modified = true;
    return edge_t(v.id(), new_neighbor.id(), &*found.it);
  }
  
  // Get next vertex in a neighborhood
  template <typename VertexContent, typename EdgeContent>
  const typename VVGraph<VertexContent, EdgeContent>::vertex_t&
  VVGraph<VertexContent, EdgeContent>::nextTo(const vertex_t &v, const vertex_t &ref, unsigned int n) const
  {
    const_neighbor_found_t found = this->findInVertex(v, ref);
    if(!found)
      return vertex_t::null;     // vertex_t(0);
    typename edge_list_t::const_iterator it = found.it;
    const edge_list_t& lst = found.neighborhood->edges;
    for(unsigned int i = 0; i < n; ++i) {
      ++it;
      if(it == lst.end())
        it = lst.begin();
    }
    return it->target;
  }
  
  // Get previous vertex in a neighborhood
  template <typename VertexContent, typename EdgeContent>
  const typename VVGraph<VertexContent, EdgeContent>::vertex_t&
  VVGraph<VertexContent, EdgeContent>::prevTo(const vertex_t &v, const vertex_t &ref, unsigned int n) const
  {
    const_neighbor_found_t found = this->findInVertex(v, ref);
    if(!found)
      return vertex_t::null;     // vertex_t(0);
    typename edge_list_t::const_iterator it = found.it;
    const edge_list_t& lst = found.neighborhood->edges;
    for(unsigned int i = 0; i < n; ++i) {
      if(it == lst.begin())
        it = lst.end();
      --it;
    }
    return it->target;
  }
  
  // Get an edge
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::edge_t
  VVGraph<VertexContent, EdgeContent>::edge(const vertex_t& src, const vertex_t& tgt)
  {
    neighbor_found_t found = this->findInVertex(src, tgt);
    if(!found)
      return edge_t();
    return edge_t(src.id(), tgt.id(), &*found.it);
  }
  
  // Get an edge
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::const_edge_t
  VVGraph<VertexContent, EdgeContent>::edge(const vertex_t& src, const vertex_t& tgt) const
  {
    const_neighbor_found_t found = this->findInVertex(src, tgt);
    if(!found)
      return const_edge_t();
    return const_edge_t(src.id(), tgt.id(), &*found.it);
  }
  
  // Insert a vertex in a neighborhood after another
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::edge_t
  VVGraph<VertexContent, EdgeContent>::spliceAfter(const vertex_t& v, const vertex_t& neighbor,
                                                          const vertex_t& new_neighbor)
  {
    if(new_neighbor.isNull() || v == new_neighbor)
      return edge_t();
    neighbor_found_t found = this->findInVertex(v, neighbor);
    if(!found)
      return edge_t();
    ++found.it;
    typename edge_list_t::iterator new_edge_it
      = found.neighborhood->edges.insert(found.it, neighbor_t(new_neighbor, EdgeContent()));
    modified = true;
    return edge_t(v.id(), new_neighbor.id(), &*new_edge_it);
  }
  
  // Insert a vertex in a neighborhood before another
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::edge_t
  VVGraph<VertexContent, EdgeContent>::spliceBefore(const vertex_t& v, const vertex_t& neighbor,
                                                           const vertex_t& new_neighbor)
  {
    if(new_neighbor.isNull() || v == new_neighbor)
      return edge_t();
    neighbor_found_t found = this->findInVertex(v, neighbor);
    if(!found)
      return edge_t();
    typename edge_list_t::iterator new_edge_it
      = found.neighborhood->edges.insert(found.it, neighbor_t(new_neighbor, EdgeContent()));
    modified = true;
    return edge_t(v.id(), new_neighbor.id(), &*new_edge_it);
  }
  
  // Insert an edge
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::edge_t
  VVGraph<VertexContent, EdgeContent>::insertEdge(const vertex_t& v, const vertex_t& new_neighbor)
  {
    if(v.isNull() || new_neighbor.isNull() || v == new_neighbor)
      return edge_t();
    NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return edge_t();
    edge_list_t& lst = it_found->single_neighborhood.edges;
    //RSS change from list to vec
    //RSS lst.push_front(neighbor_t(new_neighbor, EdgeContent()));
    //RSS typename edge_list_t::iterator it = lst.begin();
    lst.push_back(neighbor_t(new_neighbor, EdgeContent()));
    modified = true;
    return edge_t(v.id(), new_neighbor.id(), &lst.back());
  }
  
  // Return the neighborhood
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::const_neighbor_iterator_pair
  VVGraph<VertexContent, EdgeContent>::neighbors(const vertex_t& v) const
  {
    const_neighbor_iterator_pair result;
    if(v.isNull())
      return result;
    const NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return result;
    const edge_list_t& lst = it_found->single_neighborhood.edges;
    result.first = const_neighbor_iterator(lst.begin());
    result.second = const_neighbor_iterator(lst.end());
    return result;
  }

  // Return the neighborhood
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::neighbor_iterator_pair
  VVGraph<VertexContent, EdgeContent>::neighbors(const vertex_t& v)
  {
    neighbor_iterator_pair result;
    if(v.isNull())
      return result;
    NeighborhoodPair* it_found = this->findVertex(v);
    if(not it_found)
      return result;
    edge_list_t& lst = it_found->single_neighborhood.edges;
    result.first = neighbor_iterator(lst.begin());
    result.second = neighbor_iterator(lst.end());
    return result;
  }
  
  // Return the neighborhood starting at n
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::const_circ_neighbor_iterator_pair
  VVGraph<VertexContent, EdgeContent>::neighbors(const vertex_t& v, const vertex_t& n) const
  {
    const_circ_neighbor_iterator_pair result;
    if(v.isNull())
      return result;
    const_neighbor_found_t found = this->findInVertex(v, n);
    if(!found)
      return result;
    const edge_list_t& lst = found.neighborhood->edges;
    result.first = const_circ_neighbor_iterator(lst.begin(), lst.end(), found.it);
    result.second = const_circ_neighbor_iterator(lst.begin(), lst.end());
    return result;
  }
  
  // Return the neighborhood starting at n
  template <typename VertexContent, typename EdgeContent>
  typename VVGraph<VertexContent, EdgeContent>::circ_neighbor_iterator_pair
  VVGraph<VertexContent, EdgeContent>::neighbors(const vertex_t& v, const vertex_t& n)
  {
    circ_neighbor_iterator_pair result;
    if(v.isNull())
      return result;
    neighbor_found_t found = this->findInVertex(v, n);
    if(!found)
      return result;
    edge_list_t& lst = found.neighborhood->edges;
    result.first = circ_neighbor_iterator(lst.begin(), lst.end(), found.it);
    result.second = circ_neighbor_iterator(lst.begin(), lst.end());
    return result;
  }
  
  // Get a subgraph of the graph
  template <typename VertexContent, typename EdgeContent>
  template <typename VertexContainer>
  VVGraph<VertexContent, EdgeContent>
  VVGraph<VertexContent, EdgeContent>::subgraph(const VertexContainer& verts) const
  {
    typename VertexContainer::const_iterator it;
    VVGraph result;
    // First, insert the vertexes in the result graph
    for(it = verts.begin(); it != verts.end(); ++it) {
      result.insert(*it);
    }
    const NeighborhoodPair* it_orign;
    typename neighborhood_t::iterator it_n;
    //RSS typename edge_list_t::const_reverse_iterator it_sn;
    typename edge_list_t::const_iterator it_sn;
    // Then, creates the in and out edges
    for(it_n = result.neighborhood.begin(); it_n != result.neighborhood.end(); ++it_n) {
      const vertex_t& v = (*it_n)->vertex;
      // For each vertex, find out which edges to add
      it_orign = this->findVertex(v);
      const edge_list_t& orig_lst = it_orign->single_neighborhood.edges;
      edge_list_t& lst = (*it_n)->single_neighborhood.edges;
      //RSS change from list to vec
      // Going backward because it is easier to retrieve the iterator ...
      //RSS for(it_sn = orig_lst.rbegin(); it_sn != orig_lst.rend(); ++it_sn) {
      for(it_sn = orig_lst.begin(); it_sn != orig_lst.end(); ++it_sn) {
        // For each existing edge from the vertex, keep only the ones whose
        // target are on the subgraph
        if(result.contains(it_sn->target)) {
          // Insert the outgoing edge
          // RSSlst.push_front(neighbor_t(it_sn->target, *it_sn));
          lst.push_back(neighbor_t(it_sn->target, *it_sn));
        }
      }
    }
    return result;
  }
  
  // Copy constructor
  template <typename VertexContent, typename EdgeContent>
  VVGraph<VertexContent, EdgeContent>::VVGraph(const VVGraph& copy)
  {
    *this = copy;
  }
  
  // Operator ==
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::operator==(const VVGraph& other) const
  {
    if(neighborhood.size() != other.neighborhood.size())
      return false;
    for(typename neighborhood_t::const_iterator it = neighborhood.begin(); it != neighborhood.end(); ++it) {
      const vertex_t& v1 = it->vertex;
      typename lookup_t::const_iterator it_found = other.lookup.find(v1);
      if(it_found == other.lookup.end())
        return false;
      if(it->single_neighborhood.edges.empty()) {
        if(!it_found->single_neighborhood.single_neighborhood.edges.empty())
          return false;
      } else {
        const edge_list_t& lst = it->single_neighborhood.edges;
        const edge_list_t& olst = it_found->second->single_neighborhood.edges;
        if(lst.size() != olst.size())
          return false;
        const vertex_t& v2 = lst.begin()->target;
        bool found = false;
        for(typename edge_list_t::const_iterator it_olst = olst.begin(); it_olst != olst.end(); ++it_olst) {
          if(it_olst->target == v2) {
            found = true;
            for(typename edge_list_t::const_iterator it_lst = lst.begin(); it_lst != lst.end(); ++it_lst) {
              if(it_lst->target != it_olst->target)
                return false;
              ++it_olst;
              if(it_olst == olst.end())
                it_olst = olst.begin();
            }
            break;
          }
        }
        if(!found)
          return false;
      }
    }
    return true;
  }
  
  // Return graph in reverse
  template <typename VertexContent, typename EdgeContent>
  VVGraph<VertexContent, EdgeContent>& VVGraph<VertexContent, EdgeContent>::reverse()
  {
    typename neighborhood_t::iterator it;
    for(it = neighborhood.begin(); it != neighborhood.end(); ++it) {
      std::reverse((*it)->single_neighborhood.edges.begin(), (*it)->single_neighborhood.edges.end());
    }
    return *this;
  }
  
  // Operator =
  template <typename VertexContent, typename EdgeContent>
  VVGraph<VertexContent, EdgeContent>& VVGraph<VertexContent, EdgeContent>::operator=(const VVGraph& other)
  {
    typename neighborhood_t::const_iterator it_o;
    typename neighborhood_t::iterator it;
    for(it = neighborhood.begin(); it != neighborhood.end(); ++it) {
      delete *it;
      *it = 0;
    }
    neighborhood.clear();
    lookup.clear();
    it = neighborhood.end();
    for(it_o = other.neighborhood.begin(); it_o != other.neighborhood.end(); ++it_o) {
      NeighborhoodPair* np = (*it_o)->copy();
      neighborhood.push_back(np);
      lookup[(*it_o)->vertex] = np;
    }
    return *this;
  }
  
  // Swap
  template <typename VertexContent, typename EdgeContent>
  void VVGraph<VertexContent, EdgeContent>::swap(VVGraph& other)
  {
    neighborhood.swap(other.neighborhood);
    lookup.swap(other.lookup);
  }

  // Test for unique triangle, only use with triangle meshes
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::uniqueTri(const vertex_t& v, const vertex_t& n, const vertex_t& m) const
  {
    if(v <= n or v <= m or n == m or !edge(v, n) or !edge(v, m) or !edge(n, m))
      return false;
    else
      return true;
  }

  // Test for unique line, only use with graphs that always have edges in both directions
  template <typename VertexContent, typename EdgeContent>
  bool VVGraph<VertexContent, EdgeContent>::uniqueLine(const vertex_t& v, const vertex_t& n) const
  {
    if(v <= n or !edge(v, n))
      return false;
    else
      return true;
  }
}

#endif
