# HG changeset patch
# User Alpar Juttner <alpar@cs.elte.hu>
# Date 1539809803 -7200
# Node ID 959d78f3fe0e5358e3ae072c0885f0d3f204cfe5
# Parent  c8d0179a32a22dd4d8fc33951fa38baec8ecf5fb# Parent  73e29215aaa434a3ffb638710b0c873101d5a600
Merge #597

diff -r c8d0179a32a2 -r 959d78f3fe0e doc/references.bib
--- a/doc/references.bib	Wed Oct 17 19:22:52 2018 +0200
+++ b/doc/references.bib	Wed Oct 17 22:56:43 2018 +0200
@@ -43,6 +43,17 @@
   pages =        {67--118}
 }
 
+@article{VF2PP,
+  author =       {Alp\'ar J\"uttner and  P\'eter Madarasi},
+  title =        {{VF2++} — An improved subgraph isomorphism algorithm},
+  journal =      {Discrete Applied Mathematics},
+  year =         {2018},
+  volume =       {242},
+  pages =        {69--81},
+  url =          {https://doi.org/10.1016/j.dam.2018.02.018}
+}
+
+
 
 %%%%% Other libraries %%%%%%
 
diff -r c8d0179a32a2 -r 959d78f3fe0e lemon/bits/vf2_internals.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/bits/vf2_internals.h	Wed Oct 17 22:56:43 2018 +0200
@@ -0,0 +1,47 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2015-2017
+ * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef VF2_INTERNALS_H
+#define VF2_INTERNALS_H
+
+
+///\ingroup graph_properties
+///\file
+///\brief Mapping types for graph matching algorithms.
+
+namespace lemon {
+  ///\ingroup graph_isomorphism
+  ///The \ref Vf2 "VF2" algorithm is capable of finding different kind of
+  ///graph embeddings, this enum specifies their types.
+  ///
+  ///See \ref graph_isomorphism for a more detailed description.
+  enum MappingType {
+    /// Subgraph isomorphism
+    SUBGRAPH = 0,
+    /// Induced subgraph isomorphism
+    INDUCED = 1,
+    /// Graph isomorphism
+    ///
+    /// If the two graphs have the same number of nodes, than it is
+    /// equivalent to \ref INDUCED, and if they also have the same
+    /// number of edges, then it is also equivalent to \ref SUBGRAPH.
+    ///
+    /// However, using this setting is faster than the other two options.
+    ISOMORPH = 2
+  };
+}
+#endif
diff -r c8d0179a32a2 -r 959d78f3fe0e lemon/vf2.h
--- a/lemon/vf2.h	Wed Oct 17 19:22:52 2018 +0200
+++ b/lemon/vf2.h	Wed Oct 17 22:56:43 2018 +0200
@@ -2,7 +2,7 @@
  *
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2015
+ * Copyright (C) 2015-2017
  * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
  *
  * Permission to use, modify and distribute this software is granted
@@ -24,97 +24,50 @@
 
 #include <lemon/core.h>
 #include <lemon/concepts/graph.h>
-#include <lemon/dfs.h>
 #include <lemon/bfs.h>
+#include <lemon/bits/vf2_internals.h>
 
 #include <vector>
-#include <set>
 
-namespace lemon
-{
-  namespace bits
-  {
-    namespace vf2
-    {
-      class AlwaysEq
-      {
+namespace lemon {
+  namespace bits {
+    namespace vf2 {
+
+      class AlwaysEq {
       public:
         template<class T1, class T2>
-        bool operator()(T1, T2) const
-        {
+        bool operator()(T1&, T2&) const {
           return true;
         }
       };
 
       template<class M1, class M2>
-      class MapEq
-      {
+      class MapEq{
         const M1 &_m1;
         const M2 &_m2;
       public:
-        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) {}
-        bool operator()(typename M1::Key k1, typename M2::Key k2) const
-        {
+        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) { }
+        bool operator()(typename M1::Key k1, typename M2::Key k2) const {
           return _m1[k1] == _m2[k2];
         }
       };
 
       template <class G>
-      class DfsLeaveOrder : public DfsVisitor<G>
-      {
-        const G &_g;
-        std::vector<typename G::Node> &_order;
-        int i;
-      public:
-        DfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
-          : i(countNodes(g)), _g(g), _order(order)
-        {}
-        void leave(const typename G::Node &node)
-        {
-          _order[--i]=node;
-        }
-      };
-
-      template <class G>
-      class BfsLeaveOrder : public BfsVisitor<G>
-      {
+      class BfsLeaveOrder : public BfsVisitor<G> {
         int i;
         const G &_g;
         std::vector<typename G::Node> &_order;
       public:
         BfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
-          : i(0), _g(g), _order(order)
-        {}
-        void process(const typename G::Node &node)
-        {
+          : i(0), _g(g), _order(order){
+        }
+        void process(const typename G::Node &node) {
           _order[i++]=node;
         }
       };
     }
   }
 
-  ///Graph mapping types.
-
-  ///\ingroup graph_isomorphism
-  ///The \ref Vf2 "VF2" algorithm is capable of finding different kind of
-  ///embeddings, this enum specifies its type.
-  ///
-  ///See \ref graph_isomorphism for a more detailed description.
-  enum Vf2MappingType {
-    /// Subgraph isomorphism
-    SUBGRAPH = 0,
-    /// Induced subgraph isomorphism
-    INDUCED = 1,
-    /// Graph isomorphism
-
-    /// If the two graph has the same number of nodes, than it is
-    /// equivalent to \ref INDUCED, and if they also have the same
-    /// number of edges, then it is also equivalent to \ref SUBGRAPH.
-    ///
-    /// However, using this setting is faster than the other two
-    /// options.
-    ISOMORPH = 2
-  };
 
   ///%VF2 algorithm class.
 
@@ -124,271 +77,252 @@
   ///
   ///There is also a \ref vf2() "function-type interface" called \ref vf2()
   ///for the %VF2 algorithm, which is probably more convenient in most
-  ///use-cases.
+  ///use cases.
   ///
   ///\tparam G1 The type of the graph to be embedded.
-  ///The default type is \ref ListDigraph.
+  ///The default type is \ref ListGraph.
   ///\tparam G2 The type of the graph g1 will be embedded into.
-  ///The default type is \ref ListDigraph.
+  ///The default type is \ref ListGraph.
   ///\tparam M The type of the NodeMap storing the mapping.
   ///By default, it is G1::NodeMap<G2::Node>
   ///\tparam NEQ A bool-valued binary functor determinining whether a node is
-  ///mappable to another. By default it is an always true operator.
+  ///mappable to another. By default, it is an always-true operator.
   ///
   ///\sa vf2()
 #ifdef DOXYGEN
   template<class G1, class G2, class M, class NEQ >
 #else
-  template<class G1=ListDigraph,
-           class G2=ListDigraph,
+  template<class G1 = ListGraph,
+           class G2 = ListGraph,
            class M = typename G1::template NodeMap<G2::Node>,
            class NEQ = bits::vf2::AlwaysEq >
 #endif
-  class Vf2
-  {
-    //Current depth in the DFS tree.
-    int _depth;
+  class Vf2 {
+    //The graph to be embedded
+    const G1 &_g1;
+
+    //The graph into which g1 is to be embedded
+    const G2 &_g2;
+
     //Functor with bool operator()(G1::Node,G2::Node), which returns 1
-    //if and only if the 2 nodes are equivalent.
+    //if and only if the two nodes are equivalent
     NEQ _nEq;
 
+    //Current depth in the search tree
+    int _depth;
+
+    //The current mapping. _mapping[v1]=v2 iff v1 has been mapped to v2,
+    //where v1 is a node of G1 and v2 is a node of G2
+    M &_mapping;
+
+    //_order[i] is the node of g1 for which a pair is searched if depth=i
+    std::vector<typename G1::Node> _order;
+ 
+    //_conn[v2] = number of covered neighbours of v2
     typename G2::template NodeMap<int> _conn;
-    //Current mapping. We index it by the nodes of g1, and match[v] is
-    //a node of g2.
-    M &_mapping;
-    //order[i] is the node of g1, for which we find a pair if depth=i
-    std::vector<typename G1::Node> order;
-    //currEdgeIts[i] is an edge iterator, witch is last used in the ith
-    //depth to find a pair for order[i].
-    std::vector<typename G2::IncEdgeIt> currEdgeIts;
-    //The small graph.
-    const G1 &_g1;
-    //The big graph.
-    const G2 &_g2;
-    //lookup tables for cut the searchtree
-    typename G1::template NodeMap<int> rNew1t,rInOut1t;
 
-    Vf2MappingType _mapping_type;
+    //_currEdgeIts[i] is the last used edge iterator in the i-th
+    //depth to find a pair for node _order[i]
+    std::vector<typename G2::IncEdgeIt> _currEdgeIts;
+
+    //lookup tables for cutting the searchtree
+    typename G1::template NodeMap<int> _rNew1t, _rInOut1t;
+
+    MappingType _mapping_type;
+
+    bool _deallocMappingAfterUse;
 
     //cut the search tree
-    template<Vf2MappingType MT>
-    bool cut(const typename G1::Node n1,const typename G2::Node n2) const
-    {
+    template<MappingType MT>
+    bool cut(const typename G1::Node n1,const typename G2::Node n2) const {
       int rNew2=0,rInOut2=0;
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        {
-          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
-          if(_conn[currNode]>0)
-            ++rInOut2;
-          else if(MT!=SUBGRAPH&&_conn[currNode]==0)
-            ++rNew2;
-        }
-      switch(MT)
-        {
-        case INDUCED:
-          return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
-        case SUBGRAPH:
-          return rInOut1t[n1]<=rInOut2;
-        case ISOMORPH:
-          return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
-        default:
-          return false;
-        }
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
+        if(_conn[currNode]>0)
+          ++rInOut2;
+        else if(MT!=SUBGRAPH&&_conn[currNode]==0)
+          ++rNew2;
+      }
+      switch(MT) {
+      case INDUCED:
+        return _rInOut1t[n1]<=rInOut2&&_rNew1t[n1]<=rNew2;
+      case SUBGRAPH:
+        return _rInOut1t[n1]<=rInOut2;
+      case ISOMORPH:
+        return _rInOut1t[n1]==rInOut2&&_rNew1t[n1]==rNew2;
+      default:
+        return false;
+      }
     }
 
-    template<Vf2MappingType MT>
-    bool feas(const typename G1::Node n1,const typename G2::Node n2)
-    {
+    template<MappingType MT>
+    bool feas(const typename G1::Node n1,const typename G2::Node n2) {
       if(!_nEq(n1,n2))
         return 0;
 
-      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
-        {
-          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
-          if(_mapping[currNode]!=INVALID)
-            --_conn[_mapping[currNode]];
+      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+        const typename G1::Node& currNode=_g1.oppositeNode(n1,e1);
+        if(_mapping[currNode]!=INVALID)
+          --_conn[_mapping[currNode]];
+      }
+      bool isIso=1;
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& connCurrNode = _conn[_g2.oppositeNode(n2,e2)];
+        if(connCurrNode<-1)
+          ++connCurrNode;
+        else if(MT!=SUBGRAPH&&connCurrNode==-1) {
+          isIso=0;
+          break;
         }
-      bool isIso=1;
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        {
-          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
-          if(_conn[currNode]<-1)
-            ++_conn[currNode];
-          else if(MT!=SUBGRAPH&&_conn[currNode]==-1)
-            {
+      }
+
+      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+        const typename G2::Node& currNodePair=_mapping[_g1.oppositeNode(n1,e1)];
+        int& connCurrNodePair=_conn[currNodePair];
+        if(currNodePair!=INVALID&&connCurrNodePair!=-1) {
+          switch(MT) {
+          case INDUCED:
+          case ISOMORPH:
+            isIso=0;
+            break;
+          case SUBGRAPH:
+            if(connCurrNodePair<-1)
               isIso=0;
-              break;
-            }
+            break;
+          }
+          connCurrNodePair=-1;
         }
-
-      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
-        {
-          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
-          if(_mapping[currNode]!=INVALID&&_conn[_mapping[currNode]]!=-1)
-            {
-              switch(MT)
-                {
-                case INDUCED:
-                case ISOMORPH:
-                  isIso=0;
-                  break;
-                case SUBGRAPH:
-                  if(_conn[_mapping[currNode]]<-1)
-                    isIso=0;
-                  break;
-                }
-              _conn[_mapping[currNode]]=-1;
-            }
-        }
+      }
       return isIso&&cut<MT>(n1,n2);
     }
 
-    void addPair(const typename G1::Node n1,const typename G2::Node n2)
-    {
+    void addPair(const typename G1::Node n1,const typename G2::Node n2) {
       _conn[n2]=-1;
       _mapping.set(n1,n2);
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        if(_conn[_g2.oppositeNode(n2,e2)]!=-1)
-          ++_conn[_g2.oppositeNode(n2,e2)];
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn!=-1)
+          ++currConn;
+      }
     }
 
-    void subPair(const typename G1::Node n1,const typename G2::Node n2)
-    {
+    void subPair(const typename G1::Node n1,const typename G2::Node n2) {
       _conn[n2]=0;
       _mapping.set(n1,INVALID);
-      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
-        {
-          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
-          if(_conn[currNode]>0)
-            --_conn[currNode];
-          else if(_conn[currNode]==-1)
-            ++_conn[n2];
-        }
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn>0)
+          --currConn;
+        else if(currConn==-1)
+          ++_conn[n2];
+      }
     }
 
-    void setOrder()//we will find pairs for the nodes of g1 in this order
-    {
-      // bits::vf2::DfsLeaveOrder<G1> v(_g1,order);
-      //   DfsVisit<G1,bits::vf2::DfsLeaveOrder<G1> >dfs(_g1, v);
-      //   dfs.run();
-
-      //it is more efficient in practice than DFS
-      bits::vf2::BfsLeaveOrder<G1> v(_g1,order);
-      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> >bfs(_g1, v);
+    void initOrder() {
+      //determine the order in which we will find pairs for the nodes of g1
+      //BFS order is more efficient in practice than DFS
+      bits::vf2::BfsLeaveOrder<G1> v(_g1,_order);
+      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> > bfs(_g1, v);
       bfs.run();
     }
 
-    template<Vf2MappingType MT>
-    bool extMatch()
-    {
-      while(_depth>=0)
-        {
-          //there are not nodes in g1, which has not pair in g2.
-          if(_depth==static_cast<int>(order.size()))
-            {
+    template<MappingType MT>
+    bool extMatch() {
+      while(_depth>=0) {
+        if(_depth==static_cast<int>(_order.size())) {
+          //all nodes of g1 are mapped to nodes of g2
+          --_depth;
+          return true;
+        }
+        typename G1::Node& nodeOfDepth = _order[_depth];
+        const typename G2::Node& pairOfNodeOfDepth = _mapping[nodeOfDepth];
+        typename G2::IncEdgeIt &edgeItOfDepth = _currEdgeIts[_depth];
+        //the node of g2 whose neighbours are the candidates for
+        //the pair of nodeOfDepth
+        typename G2::Node currPNode;
+        if(edgeItOfDepth==INVALID) {
+          typename G1::IncEdgeIt fstMatchedE(_g1,nodeOfDepth);
+          //if pairOfNodeOfDepth!=INVALID, we don't use fstMatchedE
+          if(pairOfNodeOfDepth==INVALID) {
+            for(; fstMatchedE!=INVALID &&
+                  _mapping[_g1.oppositeNode(nodeOfDepth,
+                                            fstMatchedE)]==INVALID;
+                ++fstMatchedE) ; //find fstMatchedE
+          }
+          if(fstMatchedE==INVALID||pairOfNodeOfDepth!=INVALID) {
+            //We found no covered neighbours, this means that
+            //the graph is not connected (or _depth==0). Each
+            //uncovered (and there are some other properties due
+            //to the spec. problem types) node of g2 is
+            //candidate. We can read the iterator of the last
+            //tried node from the match if it is not the first
+            //try (match[nodeOfDepth]!=INVALID)
+            typename G2::NodeIt n2(_g2);
+            //if it's not the first try
+            if(pairOfNodeOfDepth!=INVALID) {
+              n2=++typename G2::NodeIt(_g2,pairOfNodeOfDepth);
+              subPair(nodeOfDepth,pairOfNodeOfDepth);
+            }
+            for(; n2!=INVALID; ++n2)
+              if(MT!=SUBGRAPH) {
+                if(_conn[n2]==0&&feas<MT>(nodeOfDepth,n2))
+                  break;
+              }
+              else if(_conn[n2]>=0&&feas<MT>(nodeOfDepth,n2))
+                break;
+            // n2 is the next candidate
+            if(n2!=INVALID){
+              addPair(nodeOfDepth,n2);
+              ++_depth;
+            }
+            else // there are no more candidates
               --_depth;
-              return true;
-            }
-          //the node of g2, which neighbours are the candidates for
-          //the pair of order[_depth]
-          typename G2::Node currPNode;
-          if(currEdgeIts[_depth]==INVALID)
-            {
-              typename G1::IncEdgeIt fstMatchedE(_g1,order[_depth]);
-              //if _mapping[order[_depth]]!=INVALID, we dont use
-              //fstMatchedE
-              if(_mapping[order[_depth]]==INVALID)
-                for(; fstMatchedE!=INVALID &&
-                      _mapping[_g1.oppositeNode(order[_depth],
-                                              fstMatchedE)]==INVALID;
-                    ++fstMatchedE) ; //find fstMatchedE
-              if(fstMatchedE==INVALID||_mapping[order[_depth]]!=INVALID)
-                {
-                  //We did not find an covered neighbour, this means
-                  //the graph is not connected(or _depth==0).  Every
-                  //uncovered(and there are some other properties due
-                  //to the spec. problem types) node of g2 is
-                  //candidate.  We can read the iterator of the last
-                  //tryed node from the match if it is not the first
-                  //try(match[order[_depth]]!=INVALID)
-                  typename G2::NodeIt n2(_g2);
-                  //if its not the first try
-                  if(_mapping[order[_depth]]!=INVALID)
-                    {
-                      n2=++typename G2::NodeIt(_g2,_mapping[order[_depth]]);
-                      subPair(order[_depth],_mapping[order[_depth]]);
-                    }
-                  for(; n2!=INVALID; ++n2)
-                    if(MT!=SUBGRAPH&&_conn[n2]==0)
-                      {
-                        if(feas<MT>(order[_depth],n2))
-                          break;
-                      }
-                    else if(MT==SUBGRAPH&&_conn[n2]>=0)
-                      if(feas<MT>(order[_depth],n2))
-                        break;
-                  // n2 is the next candidate
-                  if(n2!=INVALID)
-                    {
-                      addPair(order[_depth],n2);
-                      ++_depth;
-                    }
-                  else // there is no more candidate
-                    --_depth;
-                  continue;
-                }
-              else
-                {
-                  currPNode=_mapping[_g1.oppositeNode(order[_depth],
-                                                      fstMatchedE)];
-                  currEdgeIts[_depth]=typename G2::IncEdgeIt(_g2,currPNode);
-                }
-            }
-          else
-            {
-              currPNode=_g2.oppositeNode(_mapping[order[_depth]],
-                                         currEdgeIts[_depth]);
-              subPair(order[_depth],_mapping[order[_depth]]);
-              ++currEdgeIts[_depth];
-            }
-          for(; currEdgeIts[_depth]!=INVALID; ++currEdgeIts[_depth])
-            {
-              const typename G2::Node currNode =
-                _g2.oppositeNode(currPNode, currEdgeIts[_depth]);
-              //if currNode is uncovered
-              if(_conn[currNode]>0&&feas<MT>(order[_depth],currNode))
-                {
-                  addPair(order[_depth],currNode);
-                  break;
-                }
-            }
-          currEdgeIts[_depth]==INVALID?--_depth:++_depth;
+            continue;
+          }
+          else {
+            currPNode=_mapping[_g1.oppositeNode(nodeOfDepth,
+                                                fstMatchedE)];
+            edgeItOfDepth=typename G2::IncEdgeIt(_g2,currPNode);
+          }
         }
+        else {
+          currPNode=_g2.oppositeNode(pairOfNodeOfDepth,
+                                     edgeItOfDepth);
+          subPair(nodeOfDepth,pairOfNodeOfDepth);
+          ++edgeItOfDepth;
+        }
+        for(; edgeItOfDepth!=INVALID; ++edgeItOfDepth) {
+          const typename G2::Node currNode =
+            _g2.oppositeNode(currPNode, edgeItOfDepth);
+          if(_conn[currNode]>0&&feas<MT>(nodeOfDepth,currNode)) {
+            addPair(nodeOfDepth,currNode);
+            break;
+          }
+        }
+        edgeItOfDepth==INVALID?--_depth:++_depth;
+      }
       return false;
     }
 
-    //calc. the lookup table for cut the searchtree
-    void setRNew1tRInOut1t()
-    {
+    //calculate the lookup table for cutting the search tree
+    void initRNew1tRInOut1t() {
       typename G1::template NodeMap<int> tmp(_g1,0);
-      for(unsigned int i=0; i<order.size(); ++i)
-        {
-          tmp[order[i]]=-1;
-          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
-            {
-              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
-              if(tmp[currNode]>0)
-                ++rInOut1t[order[i]];
-              else if(tmp[currNode]==0)
-                ++rNew1t[order[i]];
-            }
-          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
-            {
-              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
-              if(tmp[currNode]!=-1)
-                ++tmp[currNode];
-            }
+      for(unsigned int i=0; i<_order.size(); ++i) {
+        const typename G1::Node& orderI = _order[i];
+        tmp[orderI]=-1;
+        for(typename G1::IncEdgeIt e1(_g1,orderI); e1!=INVALID; ++e1) {
+          const typename G1::Node currNode=_g1.oppositeNode(orderI,e1);
+          if(tmp[currNode]>0)
+            ++_rInOut1t[orderI];
+          else if(tmp[currNode]==0)
+            ++_rNew1t[orderI];
         }
+        for(typename G1::IncEdgeIt e1(_g1,orderI); e1!=INVALID; ++e1) {
+          const typename G1::Node currNode=_g1.oppositeNode(orderI,e1);
+          if(tmp[currNode]!=-1)
+            ++tmp[currNode];
+        }
+      }
     }
   public:
     ///Constructor
@@ -399,62 +333,73 @@
     ///\param g2 The graph \e g1 will be embedded into.
     ///\param m \ref concepts::ReadWriteMap "read-write" NodeMap
     ///storing the found mapping.
-    ///\param neq A bool-valued binary functor determinining whether a node is
+    ///\param neq A bool-valued binary functor determining whether a node is
     ///mappable to another. By default it is an always true operator.
-    Vf2(const G1 &g1, const G2 &g2,M &m, const NEQ &neq = NEQ() ) :
-      _nEq(neq), _conn(g2,0), _mapping(m), order(countNodes(g1)),
-      currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNew1t(g1,0),
-      rInOut1t(g1,0), _mapping_type(SUBGRAPH)
+    Vf2(const G1 &g1, const G2 &g2, M &m, const NEQ &neq = NEQ() ) :
+      _g1(g1), _g2(g2), _nEq(neq), _depth(0), _mapping(m),
+      _order(countNodes(g1)), _conn(g2,0),
+      _currEdgeIts(countNodes(g1),INVALID), _rNew1t(g1,0), _rInOut1t(g1,0),
+      _mapping_type(SUBGRAPH), _deallocMappingAfterUse(0)
     {
-      _depth=0;
-      setOrder();
-      setRNew1tRInOut1t();
+      initOrder();
+      initRNew1tRInOut1t();
       for(typename G1::NodeIt n(g1);n!=INVALID;++n)
         m[n]=INVALID;
     }
 
+    ///Destructor
+
+    ///Destructor.
+    ///
+
+    ~Vf2(){
+      if(_deallocMappingAfterUse)
+        delete &_mapping;
+    }
+
     ///Returns the current mapping type
 
     ///Returns the current mapping type
     ///
-    Vf2MappingType mappingType() const { return _mapping_type; }
+    MappingType mappingType() const {
+      return _mapping_type;
+    }
     ///Sets mapping type
 
     ///Sets mapping type.
     ///
     ///The mapping type is set to \ref SUBGRAPH by default.
     ///
-    ///\sa See \ref Vf2MappingType for the possible values.
-    void mappingType(Vf2MappingType m_type) { _mapping_type = m_type; }
+    ///\sa See \ref MappingType for the possible values.
+    void mappingType(MappingType m_type) {
+      _mapping_type = m_type;
+    }
 
     ///Finds a mapping
 
-    ///It finds a mapping between from g1 into g2 according to the mapping
-    ///type set by \ref mappingType(Vf2MappingType) "mappingType()".
+    ///It finds a mapping from g1 into g2 according to the mapping
+    ///type set by \ref mappingType(MappingType) "mappingType()".
     ///
     ///By subsequent calls, it returns all possible mappings one-by-one.
     ///
     ///\retval true if a mapping is found.
     ///\retval false if there is no (more) mapping.
-    bool find()
-    {
-      switch(_mapping_type)
-        {
-        case SUBGRAPH:
-          return extMatch<SUBGRAPH>();
-        case INDUCED:
-          return extMatch<INDUCED>();
-        case ISOMORPH:
-          return extMatch<ISOMORPH>();
-        default:
-          return false;
-        }
+    bool find() {
+      switch(_mapping_type) {
+      case SUBGRAPH:
+        return extMatch<SUBGRAPH>();
+      case INDUCED:
+        return extMatch<INDUCED>();
+      case ISOMORPH:
+        return extMatch<ISOMORPH>();
+      default:
+        return false;
+      }
     }
   };
 
   template<class G1, class G2>
-  class Vf2WizardBase
-  {
+  class Vf2WizardBase {
   protected:
     typedef G1 Graph1;
     typedef G2 Graph2;
@@ -462,35 +407,37 @@
     const G1 &_g1;
     const G2 &_g2;
 
-    Vf2MappingType _mapping_type;
+    MappingType _mapping_type;
 
     typedef typename G1::template NodeMap<typename G2::Node> Mapping;
     bool _local_mapping;
     void *_mapping;
-    void createMapping()
-    {
+    void createMapping() {
       _mapping = new Mapping(_g1);
     }
 
+    void *myVf2; //used in Vf2Wizard::find
+
+
     typedef bits::vf2::AlwaysEq NodeEq;
     NodeEq _node_eq;
 
     Vf2WizardBase(const G1 &g1,const G2 &g2)
-      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) {}
+      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) { }
   };
 
+
   /// Auxiliary class for the function-type interface of %VF2 algorithm.
-
+  ///
   /// This auxiliary class implements the named parameters of
   /// \ref vf2() "function-type interface" of \ref Vf2 algorithm.
   ///
-  /// \warning This class should only be used through the function \ref vf2().
+  /// \warning This class is not to be used directly.
   ///
   /// \tparam TR The traits class that defines various types used by the
   /// algorithm.
   template<class TR>
-  class Vf2Wizard : public TR
-  {
+  class Vf2Wizard : public TR {
     typedef TR Base;
     typedef typename TR::Graph1 Graph1;
     typedef typename TR::Graph2 Graph2;
@@ -511,9 +458,12 @@
     ///Copy constructor
     Vf2Wizard(const Base &b) : Base(b) {}
 
+    ///Copy constructor
+    Vf2Wizard(const Vf2Wizard &b) : Base(b) {}
+
 
     template<class T>
-    struct SetMappingBase : public Base {
+    struct SetMappingBase : public Base{
       typedef T Mapping;
       SetMappingBase(const Base &b) : Base(b) {}
     };
@@ -524,8 +474,7 @@
     ///\ref named-templ-param "Named parameter" function for setting
     ///the map that stores the found embedding.
     template<class T>
-    Vf2Wizard< SetMappingBase<T> > mapping(const T &t)
-    {
+    Vf2Wizard< SetMappingBase<T> > mapping(const T &t) {
       Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
       Base::_local_mapping = false;
       return Vf2Wizard<SetMappingBase<T> >(*this);
@@ -536,7 +485,8 @@
       typedef NE NodeEq;
       NodeEq _node_eq;
       SetNodeEqBase(const Base &b, const NE &node_eq)
-        : Base(b), _node_eq(node_eq) {}
+        : Base(b), _node_eq(node_eq){
+      }
     };
 
     ///\brief \ref named-templ-param "Named parameter" for setting
@@ -549,8 +499,7 @@
     ///whether a node is mappable to another. By default it is an
     ///always true operator.
     template<class T>
-    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq)
-    {
+    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq) {
       return Vf2Wizard<SetNodeEqBase<T> >(SetNodeEqBase<T>(*this,node_eq));
     }
 
@@ -560,16 +509,15 @@
     ///\ref named-templ-param "Named parameter" function for setting
     ///the node labels defining equivalence relation between them.
     ///
-    ///\param m1 It is arbitrary \ref concepts::ReadMap "readable node map"
+    ///\param m1 An arbitrary \ref concepts::ReadMap "readable node map"
     ///of g1.
-    ///\param m2 It is arbitrary \ref concepts::ReadMap "readable node map"
+    ///\param m2 An arbitrary \ref concepts::ReadMap "readable node map"
     ///of g2.
     ///
     ///The value type of these maps must be equal comparable.
     template<class M1, class M2>
     Vf2Wizard< SetNodeEqBase<bits::vf2::MapEq<M1,M2> > >
-    nodeLabels(const M1 &m1,const M2 &m2)
-    {
+    nodeLabels(const M1 &m1,const M2 &m2){
       return nodeEq(bits::vf2::MapEq<M1,M2>(m1,m2));
     }
 
@@ -581,9 +529,8 @@
     ///
     ///The mapping type is set to \ref SUBGRAPH by default.
     ///
-    ///\sa See \ref Vf2MappingType for the possible values.
-    Vf2Wizard<Base> &mappingType(Vf2MappingType m_type)
-    {
+    ///\sa See \ref MappingType for the possible values.
+    Vf2Wizard<Base> &mappingType(MappingType m_type) {
       _mapping_type = m_type;
       return *this;
     }
@@ -593,8 +540,7 @@
     ///
     ///\ref named-templ-param "Named parameter" for setting
     ///the mapping type to \ref INDUCED.
-    Vf2Wizard<Base> &induced()
-    {
+    Vf2Wizard<Base> &induced() {
       _mapping_type = INDUCED;
       return *this;
     }
@@ -604,20 +550,19 @@
     ///
     ///\ref named-templ-param "Named parameter" for setting
     ///the mapping type to \ref ISOMORPH.
-    Vf2Wizard<Base> &iso()
-    {
+    Vf2Wizard<Base> &iso() {
       _mapping_type = ISOMORPH;
       return *this;
     }
 
+
     ///Runs VF2 algorithm.
 
     ///This method runs VF2 algorithm.
     ///
     ///\retval true if a mapping is found.
-    ///\retval false if there is no (more) mapping.
-    bool run()
-    {
+    ///\retval false if there is no mapping.
+    bool run(){
       if(Base::_local_mapping)
         Base::createMapping();
 
@@ -633,6 +578,46 @@
 
       return ret;
     }
+
+    ///Get a pointer to the generated Vf2 object.
+
+    ///Gives a pointer to the generated Vf2 object.
+    ///
+    ///\return Pointer to the generated Vf2 object.
+    ///\warning Don't forget to delete the referred Vf2 object after use.
+    Vf2<Graph1, Graph2, Mapping, NodeEq >* getPtrToVf2Object() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+      Vf2<Graph1, Graph2, Mapping, NodeEq >* ptr =
+        new Vf2<Graph1, Graph2, Mapping, NodeEq>
+        (_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
+      ptr->mappingType(_mapping_type);
+      if(Base::_local_mapping)
+        ptr->_deallocMappingAfterUse = true;
+      return ptr;
+    }
+
+    ///Counts the number of mappings.
+
+    ///This method counts the number of mappings.
+    ///
+    /// \return The number of mappings.
+    int count() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+
+      Vf2<Graph1, Graph2, Mapping, NodeEq>
+        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
+      if(Base::_local_mapping)
+        alg._deallocMappingAfterUse = true;
+      alg.mappingType(_mapping_type);
+
+      int ret = 0;
+      while(alg.find())
+        ++ret;
+
+      return ret;
+    }
   };
 
   ///Function-type interface for VF2 algorithm.
@@ -644,20 +629,32 @@
   ///declared as the members of class \ref Vf2Wizard.
   ///The following examples show how to use these parameters.
   ///\code
-  ///  // Find an embedding of graph g into graph h
+  ///  // Find an embedding of graph g1 into graph g2
   ///  ListGraph::NodeMap<ListGraph::Node> m(g);
-  ///  vf2(g,h).mapping(m).run();
+  ///  vf2(g1,g2).mapping(m).run();
   ///
-  ///  // Check whether graphs g and h are isomorphic
-  ///  bool is_iso = vf2(g,h).iso().run();
+  ///  // Check whether graphs g1 and g2 are isomorphic
+  ///  bool is_iso = vf2(g1,g2).iso().run();
+  ///
+  ///  // Count the number of isomorphisms
+  ///  int num_isos = vf2(g1,g2).iso().count();
+  ///
+  ///  // Iterate through all the induced subgraph mappings of graph g1 into g2
+  ///  auto* myVf2 = vf2(g1,g2).mapping(m).nodeLabels(c1,c2)
+  ///  .induced().getPtrToVf2Object();
+  ///  while(myVf2->find()){
+  ///    //process the current mapping m
+  ///  }
+  ///  delete myVf22;
   ///\endcode
-  ///\warning Don't forget to put the \ref Vf2Wizard::run() "run()"
+  ///\warning Don't forget to put the \ref Vf2Wizard::run() "run()",
+  ///\ref Vf2Wizard::count() "count()" or
+  ///the \ref Vf2Wizard::getPtrToVf2Object() "getPtrToVf2Object()"
   ///to the end of the expression.
   ///\sa Vf2Wizard
   ///\sa Vf2
   template<class G1, class G2>
-  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2)
-  {
+  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2) {
     return Vf2Wizard<Vf2WizardBase<G1,G2> >(g1,g2);
   }
 
diff -r c8d0179a32a2 -r 959d78f3fe0e lemon/vf2pp.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/lemon/vf2pp.h	Wed Oct 17 22:56:43 2018 +0200
@@ -0,0 +1,860 @@
+/* -*- mode: C++; indent-tabs-mode: nil; -*-
+ *
+ * This file is a part of LEMON, a generic C++ optimization library.
+ *
+ * Copyright (C) 2015-2017
+ * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
+ *
+ * Permission to use, modify and distribute this software is granted
+ * provided that this copyright notice appears in all copies. For
+ * precise terms see the accompanying LICENSE file.
+ *
+ * This software is provided "AS IS" with no warranty of any kind,
+ * express or implied, and with no claim as to its suitability for any
+ * purpose.
+ *
+ */
+
+#ifndef LEMON_VF2PP_H
+#define LEMON_VF2PP_H
+
+///\ingroup graph_properties
+///\file
+///\brief VF2 Plus Plus algorithm.
+
+#include <lemon/core.h>
+#include <lemon/concepts/graph.h>
+#include <lemon/bits/vf2_internals.h>
+
+#include <vector>
+#include <algorithm>
+#include <utility>
+
+namespace lemon {
+
+  ///%VF2 Plus Plus algorithm class \cite VF2PP.
+
+  ///\ingroup graph_isomorphism This class provides an efficient
+  ///implementation of the %VF2 Plus Plus algorithm \cite VF2PP
+  ///for variants of the (Sub)graph Isomorphism problem.
+  ///
+  ///There is also a \ref vf2pp() "function-type interface" called
+  ///\ref vf2pp() for the %VF2 Plus Plus algorithm, which is probably
+  ///more convenient in most use cases.
+  ///
+  ///\tparam G1 The type of the graph to be embedded.
+  ///The default type is \ref ListGraph.
+  ///\tparam G2 The type of the graph g1 will be embedded into.
+  ///The default type is \ref ListGraph.
+  ///\tparam M The type of the NodeMap storing the mapping.
+  ///By default, it is G1::NodeMap<G2::Node>
+  ///\tparam M1 The type of the NodeMap storing the integer node labels of G1.
+  ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+  ///different labels. By default, it is G1::NodeMap<int>.
+  ///\tparam M2 The type of the NodeMap storing the integer node labels of G2.
+  ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+  ///different labels. By default, it is G2::NodeMap<int>.
+  ///
+  ///\sa vf2pp()
+#ifdef DOXYGEN
+  template<class G1, class G2, class M, class M1, class M2 >
+#else
+  template<class G1 = ListGraph,
+           class G2 = ListGraph,
+           class M = typename G1::template NodeMap<G2::Node>,
+           class M1 = typename G1::template NodeMap<int>,
+           class M2 = typename G2::template NodeMap<int> >
+#endif
+  class Vf2pp {
+    //The graph to be embedded
+    const G1 &_g1;
+
+    //The graph into which g1 is to be embedded
+    const G2 &_g2;
+
+    //Current depth in the search tree
+    int _depth;
+
+    //The current mapping. _mapping[v1]=v2 iff v1 has been mapped to v2,
+    //where v1 is a node of G1 and v2 is a node of G2
+    M &_mapping;
+
+    //_order[i] is a node of g1 for which a pair is searched if depth=i
+    std::vector<typename G1::Node> _order;
+
+    //_conn[v2] = number of covered neighbours of v2
+    typename G2::template NodeMap<int> _conn;
+
+    //_currEdgeIts[i] is the last used edge iterator in the i-th
+    //depth to find a pair for node _order[i]
+    std::vector<typename G2::IncEdgeIt> _currEdgeIts;
+ 
+    //_rNewLabels1[v] is a pair of form
+    //(label; num. of uncov. nodes with such label and no covered neighbours)
+    typename G1::template NodeMap<std::vector<std::pair<int,int> > >
+    _rNewLabels1;
+
+    //_rInOutLabels1[v] is the number of covered neighbours of v for each label
+    //in form (label,number of such labels)
+    typename G1::template NodeMap<std::vector<std::pair<int,int> > >
+    _rInOutLabels1;
+
+    //_intLabels1[v]==i means that node v has label i in _g1
+    //(i is in {0,1,2,..,K-1}, where K is the number of diff. labels)
+    M1 &_intLabels1;
+
+    //_intLabels2[v]==i means that node v has label i in _g2
+    //(i is in {0,1,2,..,K-1}, where K is the number of diff. labels)
+    M2 &_intLabels2;
+
+    //largest label
+    const int _maxLabel;
+
+    //lookup tables for manipulating with label class cardinalities
+    //(after use they have to be reset to 0..0)
+    std::vector<int> _labelTmp1,_labelTmp2;
+
+    MappingType _mapping_type;
+
+    //indicates whether the mapping or the labels must be deleted in the destructor
+    bool _deallocMappingAfterUse,_deallocLabelsAfterUse;
+
+
+    //improved cutting function
+    template<MappingType MT>
+    bool cutByLabels(const typename G1::Node n1,const typename G2::Node n2) {
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
+        if(_conn[currNode]>0)
+          --_labelTmp1[_intLabels2[currNode]];
+        else if(MT!=SUBGRAPH&&_conn[currNode]==0)
+          --_labelTmp2[_intLabels2[currNode]];
+      }
+
+      bool ret=1;
+      if(ret) {
+        for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
+          _labelTmp1[_rInOutLabels1[n1][i].first]+=_rInOutLabels1[n1][i].second;
+
+        if(MT!=SUBGRAPH)
+          for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
+            _labelTmp2[_rNewLabels1[n1][i].first]+=_rNewLabels1[n1][i].second;
+
+        switch(MT) {
+        case INDUCED:
+          for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
+            if(_labelTmp1[_rInOutLabels1[n1][i].first]>0) {
+              ret=0;
+              break;
+            }
+          if(ret)
+            for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
+              if(_labelTmp2[_rNewLabels1[n1][i].first]>0) {
+                ret=0;
+                break;
+              }
+          break;
+        case SUBGRAPH:
+          for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
+            if(_labelTmp1[_rInOutLabels1[n1][i].first]>0) {
+              ret=0;
+              break;
+            }
+          break;
+        case ISOMORPH:
+          for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
+            if(_labelTmp1[_rInOutLabels1[n1][i].first]!=0) {
+              ret=0;
+              break;
+            }
+          if(ret)
+            for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
+              if(_labelTmp2[_rNewLabels1[n1][i].first]!=0) {
+                ret=0;
+                break;
+              }
+          break;
+        default:
+          return false;
+        }
+        for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
+          _labelTmp1[_rInOutLabels1[n1][i].first]=0;
+
+        if(MT!=SUBGRAPH)
+          for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
+            _labelTmp2[_rNewLabels1[n1][i].first]=0;
+      }
+
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
+        _labelTmp1[_intLabels2[currNode]]=0;
+        if(MT!=SUBGRAPH&&_conn[currNode]==0)
+          _labelTmp2[_intLabels2[currNode]]=0;
+      }
+
+      return ret;
+    }
+
+
+    //try to exclude the matching of n1 and n2
+    template<MappingType MT>
+    bool feas(const typename G1::Node n1,const typename G2::Node n2) {
+      if(_intLabels1[n1]!=_intLabels2[n2])
+        return 0;
+
+      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+        const typename G1::Node& currNode=_g1.oppositeNode(n1,e1);
+        if(_mapping[currNode]!=INVALID)
+          --_conn[_mapping[currNode]];
+      }
+
+      bool isIso=1;
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& connCurrNode = _conn[_g2.oppositeNode(n2,e2)];
+        if(connCurrNode<-1)
+          ++connCurrNode;
+        else if(MT!=SUBGRAPH&&connCurrNode==-1) {
+          isIso=0;
+          break;
+        }
+      }
+
+      if(isIso)
+        for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+          const typename G2::Node& currNodePair =
+            _mapping[_g1.oppositeNode(n1,e1)];
+          int& connCurrNodePair=_conn[currNodePair];
+          if(currNodePair!=INVALID&&connCurrNodePair!=-1) {
+            switch(MT){
+            case INDUCED:
+            case ISOMORPH:
+              isIso=0;
+              break;
+            case SUBGRAPH:
+              if(connCurrNodePair<-1)
+                isIso=0;
+              break;
+            }
+            connCurrNodePair=-1;
+          }
+        }
+      else
+        for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
+          const typename G2::Node currNode=_mapping[_g1.oppositeNode(n1,e1)];
+          if(currNode!=INVALID/*&&_conn[currNode]!=-1*/)
+            _conn[currNode]=-1;
+        }
+
+      return isIso&&cutByLabels<MT>(n1,n2);
+    }
+
+    //maps n1 to n2
+    void addPair(const typename G1::Node n1,const typename G2::Node n2) {
+      _conn[n2]=-1;
+      _mapping.set(n1,n2);
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn!=-1)
+          ++currConn;
+      }
+    }
+
+    //removes mapping of n1 to n2
+    void subPair(const typename G1::Node n1,const typename G2::Node n2) {
+      _conn[n2]=0;
+      _mapping.set(n1,INVALID);
+      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2){
+        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
+        if(currConn>0)
+          --currConn;
+        else if(currConn==-1)
+          ++_conn[n2];
+      }
+    }
+
+    void processBfsTree(typename G1::Node source,unsigned int& orderIndex,
+                         typename G1::template NodeMap<int>& dm1,
+                         typename G1::template NodeMap<bool>& added) {
+      _order[orderIndex]=source;
+      added[source]=1;
+
+      unsigned int endPosOfLevel=orderIndex,
+        startPosOfLevel=orderIndex,
+        lastAdded=orderIndex;
+
+      typename G1::template NodeMap<int> currConn(_g1,0);
+
+      while(orderIndex<=lastAdded){
+        typename G1::Node currNode = _order[orderIndex];
+        for(typename G1::IncEdgeIt e(_g1,currNode); e!=INVALID; ++e) {
+          typename G1::Node n = _g1.oppositeNode(currNode,e);
+          if(!added[n]) {
+            _order[++lastAdded]=n;
+            added[n]=1;
+          }
+        }
+        if(orderIndex>endPosOfLevel){
+          for(unsigned int j = startPosOfLevel; j <= endPosOfLevel; ++j) {
+            int minInd=j;
+            for(unsigned int i = j+1; i <= endPosOfLevel; ++i)
+              if(currConn[_order[i]]>currConn[_order[minInd]]||
+                 (currConn[_order[i]]==currConn[_order[minInd]]&&
+                  (dm1[_order[i]]>dm1[_order[minInd]]||
+                   (dm1[_order[i]]==dm1[_order[minInd]]&&
+                    _labelTmp1[_intLabels1[_order[minInd]]]>
+                    _labelTmp1[_intLabels1[_order[i]]]))))
+                minInd=i;
+
+            --_labelTmp1[_intLabels1[_order[minInd]]];
+            for(typename G1::IncEdgeIt e(_g1,_order[minInd]); e!=INVALID; ++e)
+              ++currConn[_g1.oppositeNode(_order[minInd],e)];
+            std::swap(_order[j],_order[minInd]);
+          }
+          startPosOfLevel=endPosOfLevel+1;
+          endPosOfLevel=lastAdded;
+        }
+        ++orderIndex;
+      }
+    }
+
+
+    //we will find pairs for the nodes of g1 in this order
+    void initOrder(){
+      for(typename G2::NodeIt n2(_g2); n2!=INVALID; ++n2)
+        ++_labelTmp1[_intLabels2[n2]];
+
+      typename G1::template NodeMap<int> dm1(_g1,0);
+      for(typename G1::EdgeIt e(_g1); e!=INVALID; ++e) {
+        ++dm1[_g1.u(e)];
+        ++dm1[_g1.v(e)];
+      }
+
+      typename G1::template NodeMap<bool> added(_g1,0);
+      unsigned int orderIndex=0;
+
+      for(typename G1::NodeIt n(_g1); n!=INVALID;) {
+        if(!added[n]){
+          typename G1::Node minNode = n;
+          for(typename G1::NodeIt n1(_g1,minNode); n1!=INVALID; ++n1)
+            if(!added[n1] &&
+               (_labelTmp1[_intLabels1[minNode]]>
+                _labelTmp1[_intLabels1[n1]]||(dm1[minNode]<dm1[n1]&&
+                                             _labelTmp1[_intLabels1[minNode]]==
+                                             _labelTmp1[_intLabels1[n1]])))
+              minNode=n1;
+          processBfsTree(minNode,orderIndex,dm1,added);
+        }
+        else
+          ++n;
+      }
+      for(unsigned int i = 0; i < _labelTmp1.size(); ++i)
+        _labelTmp1[i]=0;
+    }
+
+
+    template<MappingType MT>
+    bool extMatch(){
+      while(_depth>=0) {
+        if(_depth==static_cast<int>(_order.size())) {
+          //all nodes of g1 are mapped to nodes of g2
+          --_depth;
+          return true;
+        }
+        typename G1::Node& nodeOfDepth = _order[_depth];
+        const typename G2::Node& pairOfNodeOfDepth = _mapping[nodeOfDepth];
+        typename G2::IncEdgeIt &edgeItOfDepth = _currEdgeIts[_depth];
+        //the node of g2 whose neighbours are the candidates for
+        //the pair of _order[_depth]
+        typename G2::Node currPNode;
+        if(edgeItOfDepth==INVALID){
+          typename G1::IncEdgeIt fstMatchedE(_g1,nodeOfDepth);
+          //if _mapping[_order[_depth]]!=INVALID, we don't need fstMatchedE
+          if(pairOfNodeOfDepth==INVALID) {
+            for(; fstMatchedE!=INVALID &&
+                  _mapping[_g1.oppositeNode(nodeOfDepth,
+                                            fstMatchedE)]==INVALID;
+                ++fstMatchedE); //find fstMatchedE, it could be preprocessed
+          }
+          if(fstMatchedE==INVALID||pairOfNodeOfDepth!=INVALID) {
+            //We found no covered neighbours, this means that
+            //the graph is not connected (or _depth==0). Each
+            //uncovered (and there are some other properties due
+            //to the spec. problem types) node of g2 is
+            //candidate. We can read the iterator of the last
+            //tried node from the match if it is not the first
+            //try (match[nodeOfDepth]!=INVALID)
+            typename G2::NodeIt n2(_g2);
+            //if it's not the first try
+            if(pairOfNodeOfDepth!=INVALID) {
+              n2=++typename G2::NodeIt(_g2,pairOfNodeOfDepth);
+              subPair(nodeOfDepth,pairOfNodeOfDepth);
+            }
+            for(; n2!=INVALID; ++n2)
+              if(MT!=SUBGRAPH) {
+                if(_conn[n2]==0&&feas<MT>(nodeOfDepth,n2))
+                  break;
+              }
+              else if(_conn[n2]>=0&&feas<MT>(nodeOfDepth,n2))
+                break;
+            // n2 is the next candidate
+            if(n2!=INVALID) {
+              addPair(nodeOfDepth,n2);
+              ++_depth;
+            }
+            else // there are no more candidates
+              --_depth;
+            continue;
+          }
+          else {
+            currPNode=_mapping[_g1.oppositeNode(nodeOfDepth,
+                                                fstMatchedE)];
+            edgeItOfDepth=typename G2::IncEdgeIt(_g2,currPNode);
+          }
+        }
+        else {
+          currPNode=_g2.oppositeNode(pairOfNodeOfDepth,
+                                     edgeItOfDepth);
+          subPair(nodeOfDepth,pairOfNodeOfDepth);
+          ++edgeItOfDepth;
+        }
+        for(; edgeItOfDepth!=INVALID; ++edgeItOfDepth) {
+          const typename G2::Node currNode =
+            _g2.oppositeNode(currPNode, edgeItOfDepth);
+          if(_conn[currNode]>0&&feas<MT>(nodeOfDepth,currNode)) {
+            addPair(nodeOfDepth,currNode);
+            break;
+          }
+        }
+        edgeItOfDepth==INVALID?--_depth:++_depth;
+      }
+      return false;
+    }
+
+    //calculate the lookup table for cutting the search tree
+    void initRNew1tRInOut1t(){
+      typename G1::template NodeMap<int> tmp(_g1,0);
+      for(unsigned int i=0; i<_order.size(); ++i) {
+        tmp[_order[i]]=-1;
+        for(typename G1::IncEdgeIt e1(_g1,_order[i]); e1!=INVALID; ++e1) {
+          const typename G1::Node currNode=_g1.oppositeNode(_order[i],e1);
+          if(tmp[currNode]>0)
+            ++_labelTmp1[_intLabels1[currNode]];
+          else if(tmp[currNode]==0)
+            ++_labelTmp2[_intLabels1[currNode]];
+        }
+        //_labelTmp1[i]=number of neightbours with label i in set rInOut
+        //_labelTmp2[i]=number of neightbours with label i in set rNew
+        for(typename G1::IncEdgeIt e1(_g1,_order[i]); e1!=INVALID; ++e1) {
+          const int& currIntLabel = _intLabels1[_g1.oppositeNode(_order[i],e1)];
+          if(_labelTmp1[currIntLabel]>0) {
+            _rInOutLabels1[_order[i]]
+              .push_back(std::make_pair(currIntLabel,
+                                        _labelTmp1[currIntLabel]));
+            _labelTmp1[currIntLabel]=0;
+          }
+          else if(_labelTmp2[currIntLabel]>0) {
+            _rNewLabels1[_order[i]].
+              push_back(std::make_pair(currIntLabel,_labelTmp2[currIntLabel]));
+            _labelTmp2[currIntLabel]=0;
+          }
+        }
+
+        for(typename G1::IncEdgeIt e1(_g1,_order[i]); e1!=INVALID; ++e1) {
+          int& tmpCurrNode=tmp[_g1.oppositeNode(_order[i],e1)];
+          if(tmpCurrNode!=-1)
+            ++tmpCurrNode;
+        }
+      }
+    }
+
+    int getMaxLabel() const{
+      int m=-1;
+      for(typename G1::NodeIt n1(_g1); n1!=INVALID; ++n1) {
+        const int& currIntLabel = _intLabels1[n1];
+        if(currIntLabel>m)
+          m=currIntLabel;
+      }
+      for(typename G2::NodeIt n2(_g2); n2!=INVALID; ++n2) {
+        const int& currIntLabel = _intLabels2[n2];
+        if(currIntLabel>m)
+          m=currIntLabel;
+      }
+      return m;
+    }
+
+  public:
+    ///Constructor
+
+    ///Constructor.
+    ///\param g1 The graph to be embedded.
+    ///\param g2 The graph \e g1 will be embedded into.
+    ///\param m The type of the NodeMap storing the mapping.
+    ///By default, it is G1::NodeMap<G2::Node>
+    ///\param intLabel1 The NodeMap storing the integer node labels of G1.
+    ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+    ///different labels.
+    ///\param intLabel1 The NodeMap storing the integer node labels of G2.
+    ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
+    ///different labels.
+    Vf2pp(const G1 &g1, const G2 &g2,M &m, M1 &intLabels1, M2 &intLabels2) :
+      _g1(g1), _g2(g2), _depth(0), _mapping(m), _order(countNodes(g1),INVALID),
+      _conn(g2,0), _currEdgeIts(countNodes(g1),INVALID), _rNewLabels1(_g1),
+      _rInOutLabels1(_g1), _intLabels1(intLabels1) ,_intLabels2(intLabels2),
+      _maxLabel(getMaxLabel()), _labelTmp1(_maxLabel+1),_labelTmp2(_maxLabel+1),
+      _mapping_type(SUBGRAPH), _deallocMappingAfterUse(0),
+      _deallocLabelsAfterUse(0)
+    {
+      initOrder();
+      initRNew1tRInOut1t();
+
+      //reset mapping
+      for(typename G1::NodeIt n(g1);n!=INVALID;++n)
+        m[n]=INVALID;
+    }
+
+    ///Destructor
+
+    ///Destructor.
+    ///
+    ~Vf2pp()
+    {
+      if(_deallocMappingAfterUse)
+        delete &_mapping;
+      if(_deallocLabelsAfterUse) {
+        delete &_intLabels1;
+        delete &_intLabels2;
+      }
+    }
+
+    ///Returns the current mapping type.
+
+    ///Returns the current mapping type.
+    ///
+    MappingType mappingType() const
+    {
+      return _mapping_type;
+    }
+
+    ///Sets the mapping type
+
+    ///Sets the mapping type.
+    ///
+    ///The mapping type is set to \ref SUBGRAPH by default.
+    ///
+    ///\sa See \ref MappingType for the possible values.
+    void mappingType(MappingType m_type)
+    {
+      _mapping_type = m_type;
+    }
+
+    ///Finds a mapping.
+
+    ///This method finds a mapping from g1 into g2 according to the mapping
+    ///type set by \ref mappingType(MappingType) "mappingType()".
+    ///
+    ///By subsequent calls, it returns all possible mappings one-by-one.
+    ///
+    ///\retval true if a mapping is found.
+    ///\retval false if there is no (more) mapping.
+    bool find()
+    {
+      switch(_mapping_type)
+        {
+        case SUBGRAPH:
+          return extMatch<SUBGRAPH>();
+        case INDUCED:
+          return extMatch<INDUCED>();
+        case ISOMORPH:
+          return extMatch<ISOMORPH>();
+        default:
+          return false;
+        }
+    }
+  };
+
+  template<typename G1, typename G2>
+  class Vf2ppWizardBase {
+  protected:
+    typedef G1 Graph1;
+    typedef G2 Graph2;
+
+    const G1 &_g1;
+    const G2 &_g2;
+
+    MappingType _mapping_type;
+
+    typedef typename G1::template NodeMap<typename G2::Node> Mapping;
+    bool _local_mapping;
+    void *_mapping;
+    void createMapping() {
+      _mapping = new Mapping(_g1);
+    }
+
+    bool _local_nodeLabels;
+    typedef typename G1::template NodeMap<int> NodeLabels1;
+    typedef typename G2::template NodeMap<int> NodeLabels2;
+    void *_nodeLabels1, *_nodeLabels2;
+    void createNodeLabels() {
+      _nodeLabels1 = new NodeLabels1(_g1,0);
+      _nodeLabels2 = new NodeLabels2(_g2,0);
+    }
+
+    Vf2ppWizardBase(const G1 &g1,const G2 &g2)
+      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH),
+        _local_mapping(1), _local_nodeLabels(1) { }
+  };
+
+
+  /// \brief Auxiliary class for the function-type interface of %VF2
+  /// Plus Plus algorithm.
+  ///
+  /// This auxiliary class implements the named parameters of
+  /// \ref vf2pp() "function-type interface" of \ref Vf2pp algorithm.
+  ///
+  /// \warning This class is not to be used directly.
+  ///
+  /// \tparam TR The traits class that defines various types used by the
+  /// algorithm.
+  template<typename TR>
+  class Vf2ppWizard : public TR {
+    typedef TR Base;
+    typedef typename TR::Graph1 Graph1;
+    typedef typename TR::Graph2 Graph2;
+    typedef typename TR::Mapping Mapping;
+    typedef typename TR::NodeLabels1 NodeLabels1;
+    typedef typename TR::NodeLabels2 NodeLabels2;
+
+    using TR::_g1;
+    using TR::_g2;
+    using TR::_mapping_type;
+    using TR::_mapping;
+    using TR::_nodeLabels1;
+    using TR::_nodeLabels2;
+
+  public:
+    ///Constructor
+    Vf2ppWizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {}
+
+    ///Copy constructor
+    Vf2ppWizard(const Base &b) : Base(b) {}
+
+
+    template<typename T>
+    struct SetMappingBase : public Base {
+      typedef T Mapping;
+      SetMappingBase(const Base &b) : Base(b) {}
+    };
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping.
+    ///
+    ///\ref named-templ-param "Named parameter" function for setting
+    ///the map that stores the found embedding.
+    template<typename T>
+    Vf2ppWizard< SetMappingBase<T> > mapping(const T &t) {
+      Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
+      Base::_local_mapping = 0;
+      return Vf2ppWizard<SetMappingBase<T> >(*this);
+    }
+
+    template<typename NL1, typename NL2>
+    struct SetNodeLabelsBase : public Base {
+      typedef NL1 NodeLabels1;
+      typedef NL2 NodeLabels2;
+      SetNodeLabelsBase(const Base &b) : Base(b) { }
+    };
+
+    ///\brief \ref named-templ-param "Named parameter" for setting the
+    ///node labels.
+    ///
+    ///\ref named-templ-param "Named parameter" function for setting
+    ///the node labels.
+    ///
+    ///\param nodeLabels1 A \ref concepts::ReadMap "readable node map"
+    ///of g1 with integer values. In case of K different labels, the labels
+    ///must be the numbers {0,1,..,K-1}.
+    ///\param nodeLabels2 A \ref concepts::ReadMap "readable node map"
+    ///of g2 with integer values. In case of K different labels, the labels
+    ///must be the numbers {0,1,..,K-1}.
+    template<typename NL1, typename NL2>
+    Vf2ppWizard< SetNodeLabelsBase<NL1,NL2> >
+    nodeLabels(const NL1 &nodeLabels1, const NL2 &nodeLabels2) {
+      Base::_local_nodeLabels = 0;
+      Base::_nodeLabels1=
+        reinterpret_cast<void*>(const_cast<NL1*>(&nodeLabels1));
+      Base::_nodeLabels2=
+        reinterpret_cast<void*>(const_cast<NL2*>(&nodeLabels2));
+      return Vf2ppWizard<SetNodeLabelsBase<NL1,NL2> >
+        (SetNodeLabelsBase<NL1,NL2>(*this));
+    }
+
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping type.
+    ///
+    ///\ref named-templ-param "Named parameter" for setting
+    ///the mapping type.
+    ///
+    ///The mapping type is set to \ref SUBGRAPH by default.
+    ///
+    ///\sa See \ref MappingType for the possible values.
+    Vf2ppWizard<Base> &mappingType(MappingType m_type) {
+      _mapping_type = m_type;
+      return *this;
+    }
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref INDUCED.
+    ///
+    ///\ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref INDUCED.
+    Vf2ppWizard<Base> &induced() {
+      _mapping_type = INDUCED;
+      return *this;
+    }
+
+    ///\brief \ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref ISOMORPH.
+    ///
+    ///\ref named-templ-param "Named parameter" for setting
+    ///the mapping type to \ref ISOMORPH.
+    Vf2ppWizard<Base> &iso() {
+      _mapping_type = ISOMORPH;
+      return *this;
+    }
+
+    ///Runs the %VF2 Plus Plus algorithm.
+
+    ///This method runs the VF2 Plus Plus algorithm.
+    ///
+    ///\retval true if a mapping is found.
+    ///\retval false if there is no mapping.
+    bool run() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+      if(Base::_local_nodeLabels)
+        Base::createNodeLabels();
+
+      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >
+        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
+            *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
+            *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
+
+      alg.mappingType(_mapping_type);
+
+      const bool ret = alg.find();
+
+      if(Base::_local_nodeLabels) {
+        delete reinterpret_cast<NodeLabels1*>(_nodeLabels1);
+        delete reinterpret_cast<NodeLabels2*>(_nodeLabels2);
+      }
+      if(Base::_local_mapping)
+        delete reinterpret_cast<Mapping*>(_mapping);
+
+      return ret;
+    }
+
+    ///Get a pointer to the generated Vf2pp object.
+
+    ///Gives a pointer to the generated Vf2pp object.
+    ///
+    ///\return Pointer to the generated Vf2pp object.
+    ///\warning Don't forget to delete the referred Vf2pp object after use.
+    Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >*
+    getPtrToVf2ppObject(){
+      if(Base::_local_mapping)
+        Base::createMapping();
+      if(Base::_local_nodeLabels)
+        Base::createNodeLabels();
+
+      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >* ptr =
+        new Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2>
+        (_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
+         *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
+         *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
+      ptr->mappingType(_mapping_type);
+      if(Base::_local_mapping)
+        ptr->_deallocMappingAfterUse=true;
+      if(Base::_local_nodeLabels)
+        ptr->_deallocLabelMapsAfterUse=true;
+
+      return ptr;
+    }
+
+    ///Counts the number of mappings.
+
+    ///This method counts the number of mappings.
+    ///
+    /// \return The number of mappings.
+    int count() {
+      if(Base::_local_mapping)
+        Base::createMapping();
+      if(Base::_local_nodeLabels)
+        Base::createNodeLabels();
+
+      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2>
+        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
+            *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
+            *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
+
+      alg.mappingType(_mapping_type);
+
+      int ret = 0;
+      while(alg.find())
+        ++ret;
+
+      if(Base::_local_nodeLabels) {
+        delete reinterpret_cast<NodeLabels1*>(_nodeLabels1);
+        delete reinterpret_cast<NodeLabels2*>(_nodeLabels2);
+      }
+      if(Base::_local_mapping)
+        delete reinterpret_cast<Mapping*>(_mapping);
+
+      return ret;
+    }
+  };
+
+
+  ///Function-type interface for VF2 Plus Plus algorithm.
+
+  /// \ingroup graph_isomorphism
+  ///Function-type interface for VF2 Plus Plus algorithm.
+  ///
+  ///This function has several \ref named-func-param "named parameters"
+  ///declared as the members of class \ref Vf2ppWizard.
+  ///The following examples show how to use these parameters.
+  ///\code
+  ///  ListGraph::NodeMap<ListGraph::Node> m(g);
+  ///  // Find an embedding of graph g1 into graph g2
+  ///  vf2pp(g1,g2).mapping(m).run();
+  ///
+  ///  // Check whether graphs g1 and g2 are isomorphic
+  ///  bool is_iso = vf2pp(g1,g2).iso().run();
+  ///
+  ///  // Count the number of isomorphisms
+  ///  int num_isos = vf2pp(g1,g2).iso().count();
+  ///
+  ///  // Iterate through all the induced subgraph mappings
+  ///  // of graph g1 into g2 using the labels c1 and c2
+  ///  auto* myVf2pp = vf2pp(g1,g2).mapping(m).nodeLabels(c1,c2)
+  ///  .induced().getPtrToVf2Object();
+  ///  while(myVf2pp->find()){
+  ///    //process the current mapping m
+  ///  }
+  ///  delete myVf22pp;
+  ///\endcode
+  ///\warning Don't forget to put the \ref Vf2ppWizard::run() "run()",
+  ///\ref Vf2ppWizard::count() "count()" or
+  ///the \ref Vf2ppWizard::getPtrToVf2ppObject() "getPtrToVf2ppObject()"
+  ///to the end of the expression.
+  ///\sa Vf2ppWizard
+  ///\sa Vf2pp
+  template<class G1, class G2>
+  Vf2ppWizard<Vf2ppWizardBase<G1,G2> > vf2pp(const G1 &g1, const G2 &g2) {
+    return Vf2ppWizard<Vf2ppWizardBase<G1,G2> >(g1,g2);
+  }
+
+}
+
+#endif
+
diff -r c8d0179a32a2 -r 959d78f3fe0e test/vf2_test.cc
--- a/test/vf2_test.cc	Wed Oct 17 19:22:52 2018 +0200
+++ b/test/vf2_test.cc	Wed Oct 17 22:56:43 2018 +0200
@@ -2,7 +2,7 @@
  *
  * This file is a part of LEMON, a generic C++ optimization library.
  *
- * Copyright (C) 2015
+ * Copyright (C) 2015-2017
  * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
  *
  * Permission to use, modify and distribute this software is granted
@@ -16,6 +16,7 @@
  */
 
 #include <lemon/vf2.h>
+#include <lemon/vf2pp.h>
 #include <lemon/concepts/digraph.h>
 #include <lemon/smart_graph.h>
 #include <lemon/lgf_reader.h>
@@ -152,49 +153,66 @@
 void make_graphs() {
   std::stringstream ss(petersen_lgf);
   graphReader(petersen, ss)
-    .nodeMap("col1",petersen_col1)
-    .nodeMap("col2",petersen_col2)
+    .nodeMap("col1", petersen_col1)
+    .nodeMap("col2", petersen_col2)
     .run();
 
   ss.clear();
   ss.str("");
-  ss<<c5_lgf;
+  ss << c5_lgf;
   //std::stringstream ss2(c5_lgf);
   graphReader(c5, ss)
-    .nodeMap("col",c5_col)
+    .nodeMap("col", c5_col)
     .run();
 
   ss.clear();
   ss.str("");
-  ss<<c7_lgf;
+  ss << c7_lgf;
   graphReader(c7, ss).run();
 
   ss.clear();
   ss.str("");
-  ss<<c10_lgf;
+  ss << c10_lgf;
   graphReader(c10, ss).run();
 
   ss.clear();
   ss.str("");
-  ss<<p10_lgf;
+  ss << p10_lgf;
   graphReader(p10, ss).run();
 
 }
 
 class EqComparable {
 public:
-  bool operator==(const EqComparable&) { return false; }
+  bool operator==(const EqComparable&) {
+    return false;
+  }
 };
 
 template<class A, class B>
 class EqClass {
 public:
-  bool operator()(A, B) { return false; }
+  bool operator()(A, B){
+    return false;
+  }
+};
+
+class IntConvertible1 {
+public:
+  operator int() {
+    return 0;
+  }
+};
+
+class IntConvertible2 {
+public:
+  operator int() {
+    return 0;
+  }
 };
 
 template<class G1,class G2>
-void checkVf2Compile()
-{
+void checkVf2Compile() {
   G1 g;
   G2 h;
   concepts::ReadWriteMap<typename G1::Node, typename G2::Node> r;
@@ -205,8 +223,15 @@
   succ = vf2(g,h).induced().run();
   succ = vf2(g,h).iso().run();
   succ = vf2(g,h).mapping(r).run();
+
+  Vf2<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+      EqClass<typename G1::Node,typename G2::Node> >
+    myVf2(g,h,r,EqClass<typename G1::Node,typename G2::Node>());
+  myVf2.find();
+
   succ = vf2(g,h).induced().mapping(r).run();
   succ = vf2(g,h).iso().mapping(r).run();
+
   concepts::ReadMap<typename G1::Node, EqComparable> l1;
   concepts::ReadMap<typename G2::Node, EqComparable> l2;
   succ = vf2(g,h).nodeLabels(l1,l2).mapping(r).run();
@@ -214,153 +239,221 @@
     .mapping(r).run();
 }
 
-void justCompile()
-{
+void vf2Compile() {
   checkVf2Compile<concepts::Graph,concepts::Graph>();
   checkVf2Compile<concepts::Graph,SmartGraph>();
   checkVf2Compile<SmartGraph,concepts::Graph>();
 }
 
-template<class G1, class G2, class I>
-void checkSub(const G1 &g1, const G2 &g2, const I &i)
-{
-  {
-    std::set<typename G2::Node> image;
-    for(typename G1::NodeIt n(g1);n!=INVALID;++n)
-      {
-          check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
-          check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
-          image.insert(i[n]);
-      }
-  }
-  for(typename G1::EdgeIt e(g1);e!=INVALID;++e)
-    check(findEdge(g2,i[g1.u(e)],i[g1.v(e)])!=INVALID,
-          "Wrong isomorphism: missing edge(checkSub).");
+template<class G1,class G2>
+void checkVf2ppCompile() {
+  G1 g;
+  G2 h;
+  concepts::ReadWriteMap<typename G1::Node, typename G2::Node> r;
+  bool succ;
+  ::lemon::ignore_unused_variable_warning(succ);
+
+  succ = vf2pp(g,h).run();
+  succ = vf2pp(g,h).induced().run();
+  succ = vf2pp(g,h).iso().run();
+  succ = vf2pp(g,h).mapping(r).run();
+  succ = vf2pp(g,h).induced().mapping(r).run();
+  succ = vf2pp(g,h).iso().mapping(r).run();
+
+  concepts::ReadMap<typename G1::Node, int> c1;
+  concepts::ReadMap<typename G2::Node, int> c2;
+  Vf2pp<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+        concepts::ReadMap<typename G1::Node, int>,
+        concepts::ReadMap<typename G2::Node, int> >
+    myVf2pp(g,h,r,c1,c2);
+  myVf2pp.find();
+
+  succ = vf2pp(g,h).nodeLabels(c1,c2).mapping(r).run();
+  succ = vf2pp(g,h).nodeLabels(c1,c2).mapping(r).run();
+
+  concepts::ReadMap<typename G1::Node, char> c1_c;
+  concepts::ReadMap<typename G2::Node, char> c2_c;
+  Vf2pp<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+        concepts::ReadMap<typename G1::Node, char>,
+        concepts::ReadMap<typename G2::Node, char> >
+    myVf2pp_c(g,h,r,c1_c,c2_c);
+  myVf2pp_c.find();
+
+  succ = vf2pp(g,h).nodeLabels(c1_c,c2_c).mapping(r).run();
+  succ = vf2pp(g,h).nodeLabels(c1_c,c2_c).mapping(r).run();
+
+  concepts::ReadMap<typename G1::Node, IntConvertible1> c1_IntConv;
+  concepts::ReadMap<typename G2::Node, IntConvertible2> c2_IntConv;
+  Vf2pp<G1,G2,concepts::ReadWriteMap<typename G1::Node, typename G2::Node>,
+        concepts::ReadMap<typename G1::Node, IntConvertible1>,
+        concepts::ReadMap<typename G2::Node, IntConvertible2> >
+    myVf2pp_IntConv(g,h,r,c1_IntConv,c2_IntConv);
+  myVf2pp_IntConv.find();
+
+  succ = vf2pp(g,h).nodeLabels(c1_IntConv,c2_IntConv).mapping(r).run();
+  succ = vf2pp(g,h).nodeLabels(c1_IntConv,c2_IntConv).mapping(r).run();
+}
+
+void vf2ppCompile() {
+  checkVf2ppCompile<concepts::Graph,concepts::Graph>();
+  checkVf2ppCompile<concepts::Graph,SmartGraph>();
+  checkVf2ppCompile<SmartGraph,concepts::Graph>();
 }
 
 template<class G1, class G2, class I>
-void checkInd(const G1 &g1, const G2 &g2, const I &i)
-  {
+void checkSubIso(const G1 &g1, const G2 &g2, const I &i) {
   std::set<typename G2::Node> image;
-  for(typename G1::NodeIt n(g1);n!=INVALID;++n)
-    {
+  for (typename G1::NodeIt n(g1);n!=INVALID;++n){
     check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
     check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
     image.insert(i[n]);
-    }
-  for(typename G1::NodeIt n(g1); n!=INVALID; ++n)
-    for(typename G1::NodeIt m(g1); m!=INVALID; ++m)
-      if((findEdge(g1,n,m)==INVALID) != (findEdge(g2,i[n],i[m])==INVALID))
-        {
+  }
+  for (typename G1::EdgeIt e(g1);e!=INVALID;++e) {
+    check(findEdge(g2,i[g1.u(e)],i[g1.v(e)])!=INVALID,
+          "Wrong isomorphism: missing edge(checkSub).");
+  }
+}
+
+template<class G1, class G2, class I>
+void checkIndIso(const G1 &g1, const G2 &g2, const I &i) {
+  std::set<typename G2::Node> image;
+  for (typename G1::NodeIt n(g1);n!=INVALID;++n) {
+    check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
+    check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
+    image.insert(i[n]);
+  }
+  for (typename G1::NodeIt n(g1); n!=INVALID; ++n) {
+    for (typename G1::NodeIt m(g1); m!=INVALID; ++m) {
+      if((findEdge(g1,n,m)==INVALID) != (findEdge(g2,i[n],i[m])==INVALID)) {
         std::cout << "Wrong isomorphism: edge mismatch";
         exit(1);
-        }
+      }
+    }
   }
-
-template<class G1,class G2>
-int checkSub(const G1 &g1, const G2 &g2)
-{
-  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).mapping(iso).run())
-    {
-      checkSub(g1,g2,iso);
-      return true;
-    }
-  else return false;
 }
 
-template<class G1,class G2>
-int checkInd(const G1 &g1, const G2 &g2)
-{
+template<class G1,class G2,class T>
+bool checkSub(const G1 &g1, const G2 &g2, const T &vf2) {
   typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).induced().mapping(iso).run())
-    {
-      checkInd(g1,g2,iso);
-      return true;
-    }
-  else return false;
+  if (const_cast<T&>(vf2).mapping(iso).run()) {
+    checkSubIso(g1,g2,iso);
+    return true;
+  }
+  return false;
 }
 
-template<class G1,class G2>
-int checkIso(const G1 &g1, const G2 &g2)
-{
+template<class G1,class G2,class T>
+bool checkInd(const G1 &g1, const G2 &g2, const T &vf2) {
   typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).iso().mapping(iso).run())
-    {
-      check(countNodes(g1)==countNodes(g2),
-            "Wrong iso alg.: they are not isomophic.");
-      checkInd(g1,g2,iso);
-      return true;
-    }
-  else return false;
+  if (const_cast<T&>(vf2).induced().mapping(iso).run()) {
+    checkIndIso(g1,g2,iso);
+    return true;
+  }
+  return false;
+}
+
+template<class G1,class G2,class T>
+bool checkIso(const G1 &g1, const G2 &g2, const T &vf2) {
+  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
+  if (const_cast<T&>(vf2).iso().mapping(iso).run()) {
+    check(countNodes(g1)==countNodes(g2),
+          "Wrong iso alg.: they are not isomophic.");
+    checkIndIso(g1,g2,iso);
+    return true;
+  }
+  return false;
 }
 
 template<class G1, class G2, class L1, class L2, class I>
 void checkLabel(const G1 &g1, const G2 &,
-                const L1 &l1, const L2 &l2,const I &i)
-{
-  for(typename G1::NodeIt n(g1);n!=INVALID;++n)
-    {
-      check(l1[n]==l2[i[n]],"Wrong isomorphism: label mismatch.");
-    }
+                const L1 &l1, const L2 &l2, const I &i) {
+  for (typename G1::NodeIt n(g1);n!=INVALID;++n) {
+    check(l1[n]==l2[i[n]],"Wrong isomorphism: label mismatch.");
+  }
+}
+
+template<class G1,class G2,class L1,class L2,class T>
+bool checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2, const T &vf2) {
+  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
+  if (const_cast<T&>(vf2).nodeLabels(l1,l2).mapping(iso).run()){
+    checkSubIso(g1,g2,iso);
+    checkLabel(g1,g2,l1,l2,iso);
+    return true;
+  }
+  return false;
+}
+
+template<class G1,class G2>
+void checkSub(const G1 &g1, const G2 &g2, bool expected, const char* msg) {
+  check(checkSub(g1,g2,vf2(g1,g2)) == expected, msg);
+  check(checkSub(g1,g2,vf2pp(g1,g2)) == expected, msg);
+}
+
+template<class G1,class G2>
+void checkInd(const G1 &g1, const G2 &g2, bool expected, const char* msg) {
+  check(checkInd(g1,g2,vf2(g1,g2)) == expected, msg);
+  check(checkInd(g1,g2,vf2pp(g1,g2)) == expected, msg);
+}
+
+template<class G1,class G2>
+void checkIso(const G1 &g1, const G2 &g2, bool expected, const char* msg) {
+  check(checkIso(g1,g2,vf2(g1,g2)) == expected, msg);
+  check(checkIso(g1,g2,vf2pp(g1,g2)) == expected, msg);
 }
 
 template<class G1,class G2,class L1,class L2>
-int checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2)
-{
-  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
-  if(vf2(g1,g2).nodeLabels(l1,l2).mapping(iso).run())
-    {
-      checkSub(g1,g2,iso);
-      checkLabel(g1,g2,l1,l2,iso);
-      return true;
-    }
-  else return false;
+void checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2, bool expected, const char* msg) {
+  check(checkSub(g1,g2,l1,l2,vf2(g1,g2)) == expected, msg);
+  check(checkSub(g1,g2,l1,l2,vf2pp(g1,g2)) == expected, msg);
 }
 
 int main() {
   make_graphs();
-  check(checkSub(c5,petersen), "There should exist a C5->Petersen mapping.");
-  check(!checkSub(c7,petersen),
-        "There should not exist a C7->Petersen mapping.");
-  check(checkSub(p10,petersen), "There should exist a P10->Petersen mapping.");
-  check(!checkSub(c10,petersen),
-        "There should not exist a C10->Petersen mapping.");
-  check(checkSub(petersen,petersen),
-        "There should exist a Petersen->Petersen mapping.");
 
-  check(checkInd(c5,petersen),
-        "There should exist a C5->Petersen spanned mapping.");
-  check(!checkInd(c7,petersen),
-        "There should exist a C7->Petersen spanned mapping.");
-  check(!checkInd(p10,petersen),
-        "There should not exist a P10->Petersen spanned mapping.");
-  check(!checkInd(c10,petersen),
-        "There should not exist a C10->Petersen spanned mapping.");
-  check(checkInd(petersen,petersen),
+  checkSub(c5,petersen,true,
+      "There should exist a C5->Petersen mapping.");
+  checkSub(c7,petersen,false,
+      "There should not exist a C7->Petersen mapping.");
+  checkSub(p10,petersen,true,
+      "There should exist a P10->Petersen mapping.");
+  checkSub(c10,petersen,false,
+      "There should not exist a C10->Petersen mapping.");
+  checkSub(petersen,petersen,true,
+      "There should exist a Petersen->Petersen mapping.");
+
+  checkInd(c5,petersen,true,
+      "There should exist a C5->Petersen spanned mapping.");
+  checkInd(c7,petersen,false,
+      "There should exist a C7->Petersen spanned mapping.");
+  checkInd(p10,petersen,false,
+      "There should not exist a P10->Petersen spanned mapping.");
+  checkInd(c10,petersen,false,
+      "There should not exist a C10->Petersen spanned mapping.");
+  checkInd(petersen,petersen,true,
         "There should exist a Petersen->Petersen spanned mapping.");
 
-  check(!checkSub(petersen,c10),
-        "There should not exist a Petersen->C10 mapping.");
-  check(checkSub(p10,c10),
-        "There should exist a P10->C10 mapping.");
-  check(!checkInd(p10,c10),
-        "There should not exist a P10->C10 spanned mapping.");
-  check(!checkSub(c10,p10),
-        "There should not exist a C10->P10 mapping.");
+  checkSub(petersen,c10,false,
+      "There should not exist a Petersen->C10 mapping.");
+  checkSub(p10,c10,true,
+      "There should exist a P10->C10 mapping.");
+  checkInd(p10,c10,false,
+      "There should not exist a P10->C10 spanned mapping.");
+  checkSub(c10,p10,false,
+      "There should not exist a C10->P10 mapping.");
 
-  check(!checkIso(p10,c10),
-        "P10 and C10 are not isomorphic.");
-  check(checkIso(c10,c10),
-        "C10 and C10 are isomorphic.");
+  checkIso(p10,c10,false,
+      "P10 and C10 are not isomorphic.");
+  checkIso(c10,c10,true,
+      "C10 and C10 are isomorphic.");
 
-  check(!vf2(p10,c10).iso().run(),
-        "P10 and C10 are not isomorphic.");
-  check(vf2(c10,c10).iso().run(),
-        "C10 and C10 are isomorphic.");
+  checkSub(c5,petersen,c5_col,petersen_col1,false,
+      "There should exist a C5->Petersen mapping.");
+  checkSub(c5,petersen,c5_col,petersen_col2,true,
+      "There should exist a C5->Petersen mapping.");
 
-  check(!checkSub(c5,petersen,c5_col,petersen_col1),
-        "There should exist a C5->Petersen mapping.");
-  check(checkSub(c5,petersen,c5_col,petersen_col2),
-        "There should exist a C5->Petersen mapping.");
+  check(!vf2(p10,c10).iso().run(), "P10 and C10 are not isomorphic.");
+  check(!vf2pp(p10,c10).iso().run(), "P10 and C10 are not isomorphic.");
+
+  check(vf2(c10,c10).iso().run(), "C10 and C10 are isomorphic.");
+  check(vf2pp(c10,c10).iso().run(), "C10 and C10 are isomorphic.");
 }