lemon/vf2.h
changeset 1415 959d78f3fe0e
parent 1410 d9f79b81ef6c
     1.1 --- a/lemon/vf2.h	Wed Oct 17 19:22:52 2018 +0200
     1.2 +++ b/lemon/vf2.h	Wed Oct 17 22:56:43 2018 +0200
     1.3 @@ -2,7 +2,7 @@
     1.4   *
     1.5   * This file is a part of LEMON, a generic C++ optimization library.
     1.6   *
     1.7 - * Copyright (C) 2015
     1.8 + * Copyright (C) 2015-2017
     1.9   * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
    1.10   *
    1.11   * Permission to use, modify and distribute this software is granted
    1.12 @@ -24,97 +24,50 @@
    1.13  
    1.14  #include <lemon/core.h>
    1.15  #include <lemon/concepts/graph.h>
    1.16 -#include <lemon/dfs.h>
    1.17  #include <lemon/bfs.h>
    1.18 +#include <lemon/bits/vf2_internals.h>
    1.19  
    1.20  #include <vector>
    1.21 -#include <set>
    1.22  
    1.23 -namespace lemon
    1.24 -{
    1.25 -  namespace bits
    1.26 -  {
    1.27 -    namespace vf2
    1.28 -    {
    1.29 -      class AlwaysEq
    1.30 -      {
    1.31 +namespace lemon {
    1.32 +  namespace bits {
    1.33 +    namespace vf2 {
    1.34 +
    1.35 +      class AlwaysEq {
    1.36        public:
    1.37          template<class T1, class T2>
    1.38 -        bool operator()(T1, T2) const
    1.39 -        {
    1.40 +        bool operator()(T1&, T2&) const {
    1.41            return true;
    1.42          }
    1.43        };
    1.44  
    1.45        template<class M1, class M2>
    1.46 -      class MapEq
    1.47 -      {
    1.48 +      class MapEq{
    1.49          const M1 &_m1;
    1.50          const M2 &_m2;
    1.51        public:
    1.52 -        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) {}
    1.53 -        bool operator()(typename M1::Key k1, typename M2::Key k2) const
    1.54 -        {
    1.55 +        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) { }
    1.56 +        bool operator()(typename M1::Key k1, typename M2::Key k2) const {
    1.57            return _m1[k1] == _m2[k2];
    1.58          }
    1.59        };
    1.60  
    1.61        template <class G>
    1.62 -      class DfsLeaveOrder : public DfsVisitor<G>
    1.63 -      {
    1.64 -        const G &_g;
    1.65 -        std::vector<typename G::Node> &_order;
    1.66 -        int i;
    1.67 -      public:
    1.68 -        DfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
    1.69 -          : i(countNodes(g)), _g(g), _order(order)
    1.70 -        {}
    1.71 -        void leave(const typename G::Node &node)
    1.72 -        {
    1.73 -          _order[--i]=node;
    1.74 -        }
    1.75 -      };
    1.76 -
    1.77 -      template <class G>
    1.78 -      class BfsLeaveOrder : public BfsVisitor<G>
    1.79 -      {
    1.80 +      class BfsLeaveOrder : public BfsVisitor<G> {
    1.81          int i;
    1.82          const G &_g;
    1.83          std::vector<typename G::Node> &_order;
    1.84        public:
    1.85          BfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
    1.86 -          : i(0), _g(g), _order(order)
    1.87 -        {}
    1.88 -        void process(const typename G::Node &node)
    1.89 -        {
    1.90 +          : i(0), _g(g), _order(order){
    1.91 +        }
    1.92 +        void process(const typename G::Node &node) {
    1.93            _order[i++]=node;
    1.94          }
    1.95        };
    1.96      }
    1.97    }
    1.98  
    1.99 -  ///Graph mapping types.
   1.100 -
   1.101 -  ///\ingroup graph_isomorphism
   1.102 -  ///The \ref Vf2 "VF2" algorithm is capable of finding different kind of
   1.103 -  ///embeddings, this enum specifies its type.
   1.104 -  ///
   1.105 -  ///See \ref graph_isomorphism for a more detailed description.
   1.106 -  enum Vf2MappingType {
   1.107 -    /// Subgraph isomorphism
   1.108 -    SUBGRAPH = 0,
   1.109 -    /// Induced subgraph isomorphism
   1.110 -    INDUCED = 1,
   1.111 -    /// Graph isomorphism
   1.112 -
   1.113 -    /// If the two graph has the same number of nodes, than it is
   1.114 -    /// equivalent to \ref INDUCED, and if they also have the same
   1.115 -    /// number of edges, then it is also equivalent to \ref SUBGRAPH.
   1.116 -    ///
   1.117 -    /// However, using this setting is faster than the other two
   1.118 -    /// options.
   1.119 -    ISOMORPH = 2
   1.120 -  };
   1.121  
   1.122    ///%VF2 algorithm class.
   1.123  
   1.124 @@ -124,271 +77,252 @@
   1.125    ///
   1.126    ///There is also a \ref vf2() "function-type interface" called \ref vf2()
   1.127    ///for the %VF2 algorithm, which is probably more convenient in most
   1.128 -  ///use-cases.
   1.129 +  ///use cases.
   1.130    ///
   1.131    ///\tparam G1 The type of the graph to be embedded.
   1.132 -  ///The default type is \ref ListDigraph.
   1.133 +  ///The default type is \ref ListGraph.
   1.134    ///\tparam G2 The type of the graph g1 will be embedded into.
   1.135 -  ///The default type is \ref ListDigraph.
   1.136 +  ///The default type is \ref ListGraph.
   1.137    ///\tparam M The type of the NodeMap storing the mapping.
   1.138    ///By default, it is G1::NodeMap<G2::Node>
   1.139    ///\tparam NEQ A bool-valued binary functor determinining whether a node is
   1.140 -  ///mappable to another. By default it is an always true operator.
   1.141 +  ///mappable to another. By default, it is an always-true operator.
   1.142    ///
   1.143    ///\sa vf2()
   1.144  #ifdef DOXYGEN
   1.145    template<class G1, class G2, class M, class NEQ >
   1.146  #else
   1.147 -  template<class G1=ListDigraph,
   1.148 -           class G2=ListDigraph,
   1.149 +  template<class G1 = ListGraph,
   1.150 +           class G2 = ListGraph,
   1.151             class M = typename G1::template NodeMap<G2::Node>,
   1.152             class NEQ = bits::vf2::AlwaysEq >
   1.153  #endif
   1.154 -  class Vf2
   1.155 -  {
   1.156 -    //Current depth in the DFS tree.
   1.157 -    int _depth;
   1.158 +  class Vf2 {
   1.159 +    //The graph to be embedded
   1.160 +    const G1 &_g1;
   1.161 +
   1.162 +    //The graph into which g1 is to be embedded
   1.163 +    const G2 &_g2;
   1.164 +
   1.165      //Functor with bool operator()(G1::Node,G2::Node), which returns 1
   1.166 -    //if and only if the 2 nodes are equivalent.
   1.167 +    //if and only if the two nodes are equivalent
   1.168      NEQ _nEq;
   1.169  
   1.170 +    //Current depth in the search tree
   1.171 +    int _depth;
   1.172 +
   1.173 +    //The current mapping. _mapping[v1]=v2 iff v1 has been mapped to v2,
   1.174 +    //where v1 is a node of G1 and v2 is a node of G2
   1.175 +    M &_mapping;
   1.176 +
   1.177 +    //_order[i] is the node of g1 for which a pair is searched if depth=i
   1.178 +    std::vector<typename G1::Node> _order;
   1.179 + 
   1.180 +    //_conn[v2] = number of covered neighbours of v2
   1.181      typename G2::template NodeMap<int> _conn;
   1.182 -    //Current mapping. We index it by the nodes of g1, and match[v] is
   1.183 -    //a node of g2.
   1.184 -    M &_mapping;
   1.185 -    //order[i] is the node of g1, for which we find a pair if depth=i
   1.186 -    std::vector<typename G1::Node> order;
   1.187 -    //currEdgeIts[i] is an edge iterator, witch is last used in the ith
   1.188 -    //depth to find a pair for order[i].
   1.189 -    std::vector<typename G2::IncEdgeIt> currEdgeIts;
   1.190 -    //The small graph.
   1.191 -    const G1 &_g1;
   1.192 -    //The big graph.
   1.193 -    const G2 &_g2;
   1.194 -    //lookup tables for cut the searchtree
   1.195 -    typename G1::template NodeMap<int> rNew1t,rInOut1t;
   1.196  
   1.197 -    Vf2MappingType _mapping_type;
   1.198 +    //_currEdgeIts[i] is the last used edge iterator in the i-th
   1.199 +    //depth to find a pair for node _order[i]
   1.200 +    std::vector<typename G2::IncEdgeIt> _currEdgeIts;
   1.201 +
   1.202 +    //lookup tables for cutting the searchtree
   1.203 +    typename G1::template NodeMap<int> _rNew1t, _rInOut1t;
   1.204 +
   1.205 +    MappingType _mapping_type;
   1.206 +
   1.207 +    bool _deallocMappingAfterUse;
   1.208  
   1.209      //cut the search tree
   1.210 -    template<Vf2MappingType MT>
   1.211 -    bool cut(const typename G1::Node n1,const typename G2::Node n2) const
   1.212 -    {
   1.213 +    template<MappingType MT>
   1.214 +    bool cut(const typename G1::Node n1,const typename G2::Node n2) const {
   1.215        int rNew2=0,rInOut2=0;
   1.216 -      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.217 -        {
   1.218 -          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.219 -          if(_conn[currNode]>0)
   1.220 -            ++rInOut2;
   1.221 -          else if(MT!=SUBGRAPH&&_conn[currNode]==0)
   1.222 -            ++rNew2;
   1.223 -        }
   1.224 -      switch(MT)
   1.225 -        {
   1.226 -        case INDUCED:
   1.227 -          return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
   1.228 -        case SUBGRAPH:
   1.229 -          return rInOut1t[n1]<=rInOut2;
   1.230 -        case ISOMORPH:
   1.231 -          return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
   1.232 -        default:
   1.233 -          return false;
   1.234 -        }
   1.235 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.236 +        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.237 +        if(_conn[currNode]>0)
   1.238 +          ++rInOut2;
   1.239 +        else if(MT!=SUBGRAPH&&_conn[currNode]==0)
   1.240 +          ++rNew2;
   1.241 +      }
   1.242 +      switch(MT) {
   1.243 +      case INDUCED:
   1.244 +        return _rInOut1t[n1]<=rInOut2&&_rNew1t[n1]<=rNew2;
   1.245 +      case SUBGRAPH:
   1.246 +        return _rInOut1t[n1]<=rInOut2;
   1.247 +      case ISOMORPH:
   1.248 +        return _rInOut1t[n1]==rInOut2&&_rNew1t[n1]==rNew2;
   1.249 +      default:
   1.250 +        return false;
   1.251 +      }
   1.252      }
   1.253  
   1.254 -    template<Vf2MappingType MT>
   1.255 -    bool feas(const typename G1::Node n1,const typename G2::Node n2)
   1.256 -    {
   1.257 +    template<MappingType MT>
   1.258 +    bool feas(const typename G1::Node n1,const typename G2::Node n2) {
   1.259        if(!_nEq(n1,n2))
   1.260          return 0;
   1.261  
   1.262 -      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
   1.263 -        {
   1.264 -          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
   1.265 -          if(_mapping[currNode]!=INVALID)
   1.266 -            --_conn[_mapping[currNode]];
   1.267 +      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
   1.268 +        const typename G1::Node& currNode=_g1.oppositeNode(n1,e1);
   1.269 +        if(_mapping[currNode]!=INVALID)
   1.270 +          --_conn[_mapping[currNode]];
   1.271 +      }
   1.272 +      bool isIso=1;
   1.273 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.274 +        int& connCurrNode = _conn[_g2.oppositeNode(n2,e2)];
   1.275 +        if(connCurrNode<-1)
   1.276 +          ++connCurrNode;
   1.277 +        else if(MT!=SUBGRAPH&&connCurrNode==-1) {
   1.278 +          isIso=0;
   1.279 +          break;
   1.280          }
   1.281 -      bool isIso=1;
   1.282 -      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.283 -        {
   1.284 -          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.285 -          if(_conn[currNode]<-1)
   1.286 -            ++_conn[currNode];
   1.287 -          else if(MT!=SUBGRAPH&&_conn[currNode]==-1)
   1.288 -            {
   1.289 +      }
   1.290 +
   1.291 +      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
   1.292 +        const typename G2::Node& currNodePair=_mapping[_g1.oppositeNode(n1,e1)];
   1.293 +        int& connCurrNodePair=_conn[currNodePair];
   1.294 +        if(currNodePair!=INVALID&&connCurrNodePair!=-1) {
   1.295 +          switch(MT) {
   1.296 +          case INDUCED:
   1.297 +          case ISOMORPH:
   1.298 +            isIso=0;
   1.299 +            break;
   1.300 +          case SUBGRAPH:
   1.301 +            if(connCurrNodePair<-1)
   1.302                isIso=0;
   1.303 -              break;
   1.304 -            }
   1.305 +            break;
   1.306 +          }
   1.307 +          connCurrNodePair=-1;
   1.308          }
   1.309 -
   1.310 -      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
   1.311 -        {
   1.312 -          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
   1.313 -          if(_mapping[currNode]!=INVALID&&_conn[_mapping[currNode]]!=-1)
   1.314 -            {
   1.315 -              switch(MT)
   1.316 -                {
   1.317 -                case INDUCED:
   1.318 -                case ISOMORPH:
   1.319 -                  isIso=0;
   1.320 -                  break;
   1.321 -                case SUBGRAPH:
   1.322 -                  if(_conn[_mapping[currNode]]<-1)
   1.323 -                    isIso=0;
   1.324 -                  break;
   1.325 -                }
   1.326 -              _conn[_mapping[currNode]]=-1;
   1.327 -            }
   1.328 -        }
   1.329 +      }
   1.330        return isIso&&cut<MT>(n1,n2);
   1.331      }
   1.332  
   1.333 -    void addPair(const typename G1::Node n1,const typename G2::Node n2)
   1.334 -    {
   1.335 +    void addPair(const typename G1::Node n1,const typename G2::Node n2) {
   1.336        _conn[n2]=-1;
   1.337        _mapping.set(n1,n2);
   1.338 -      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.339 -        if(_conn[_g2.oppositeNode(n2,e2)]!=-1)
   1.340 -          ++_conn[_g2.oppositeNode(n2,e2)];
   1.341 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.342 +        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
   1.343 +        if(currConn!=-1)
   1.344 +          ++currConn;
   1.345 +      }
   1.346      }
   1.347  
   1.348 -    void subPair(const typename G1::Node n1,const typename G2::Node n2)
   1.349 -    {
   1.350 +    void subPair(const typename G1::Node n1,const typename G2::Node n2) {
   1.351        _conn[n2]=0;
   1.352        _mapping.set(n1,INVALID);
   1.353 -      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.354 -        {
   1.355 -          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.356 -          if(_conn[currNode]>0)
   1.357 -            --_conn[currNode];
   1.358 -          else if(_conn[currNode]==-1)
   1.359 -            ++_conn[n2];
   1.360 -        }
   1.361 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.362 +        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
   1.363 +        if(currConn>0)
   1.364 +          --currConn;
   1.365 +        else if(currConn==-1)
   1.366 +          ++_conn[n2];
   1.367 +      }
   1.368      }
   1.369  
   1.370 -    void setOrder()//we will find pairs for the nodes of g1 in this order
   1.371 -    {
   1.372 -      // bits::vf2::DfsLeaveOrder<G1> v(_g1,order);
   1.373 -      //   DfsVisit<G1,bits::vf2::DfsLeaveOrder<G1> >dfs(_g1, v);
   1.374 -      //   dfs.run();
   1.375 -
   1.376 -      //it is more efficient in practice than DFS
   1.377 -      bits::vf2::BfsLeaveOrder<G1> v(_g1,order);
   1.378 -      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> >bfs(_g1, v);
   1.379 +    void initOrder() {
   1.380 +      //determine the order in which we will find pairs for the nodes of g1
   1.381 +      //BFS order is more efficient in practice than DFS
   1.382 +      bits::vf2::BfsLeaveOrder<G1> v(_g1,_order);
   1.383 +      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> > bfs(_g1, v);
   1.384        bfs.run();
   1.385      }
   1.386  
   1.387 -    template<Vf2MappingType MT>
   1.388 -    bool extMatch()
   1.389 -    {
   1.390 -      while(_depth>=0)
   1.391 -        {
   1.392 -          //there are not nodes in g1, which has not pair in g2.
   1.393 -          if(_depth==static_cast<int>(order.size()))
   1.394 -            {
   1.395 +    template<MappingType MT>
   1.396 +    bool extMatch() {
   1.397 +      while(_depth>=0) {
   1.398 +        if(_depth==static_cast<int>(_order.size())) {
   1.399 +          //all nodes of g1 are mapped to nodes of g2
   1.400 +          --_depth;
   1.401 +          return true;
   1.402 +        }
   1.403 +        typename G1::Node& nodeOfDepth = _order[_depth];
   1.404 +        const typename G2::Node& pairOfNodeOfDepth = _mapping[nodeOfDepth];
   1.405 +        typename G2::IncEdgeIt &edgeItOfDepth = _currEdgeIts[_depth];
   1.406 +        //the node of g2 whose neighbours are the candidates for
   1.407 +        //the pair of nodeOfDepth
   1.408 +        typename G2::Node currPNode;
   1.409 +        if(edgeItOfDepth==INVALID) {
   1.410 +          typename G1::IncEdgeIt fstMatchedE(_g1,nodeOfDepth);
   1.411 +          //if pairOfNodeOfDepth!=INVALID, we don't use fstMatchedE
   1.412 +          if(pairOfNodeOfDepth==INVALID) {
   1.413 +            for(; fstMatchedE!=INVALID &&
   1.414 +                  _mapping[_g1.oppositeNode(nodeOfDepth,
   1.415 +                                            fstMatchedE)]==INVALID;
   1.416 +                ++fstMatchedE) ; //find fstMatchedE
   1.417 +          }
   1.418 +          if(fstMatchedE==INVALID||pairOfNodeOfDepth!=INVALID) {
   1.419 +            //We found no covered neighbours, this means that
   1.420 +            //the graph is not connected (or _depth==0). Each
   1.421 +            //uncovered (and there are some other properties due
   1.422 +            //to the spec. problem types) node of g2 is
   1.423 +            //candidate. We can read the iterator of the last
   1.424 +            //tried node from the match if it is not the first
   1.425 +            //try (match[nodeOfDepth]!=INVALID)
   1.426 +            typename G2::NodeIt n2(_g2);
   1.427 +            //if it's not the first try
   1.428 +            if(pairOfNodeOfDepth!=INVALID) {
   1.429 +              n2=++typename G2::NodeIt(_g2,pairOfNodeOfDepth);
   1.430 +              subPair(nodeOfDepth,pairOfNodeOfDepth);
   1.431 +            }
   1.432 +            for(; n2!=INVALID; ++n2)
   1.433 +              if(MT!=SUBGRAPH) {
   1.434 +                if(_conn[n2]==0&&feas<MT>(nodeOfDepth,n2))
   1.435 +                  break;
   1.436 +              }
   1.437 +              else if(_conn[n2]>=0&&feas<MT>(nodeOfDepth,n2))
   1.438 +                break;
   1.439 +            // n2 is the next candidate
   1.440 +            if(n2!=INVALID){
   1.441 +              addPair(nodeOfDepth,n2);
   1.442 +              ++_depth;
   1.443 +            }
   1.444 +            else // there are no more candidates
   1.445                --_depth;
   1.446 -              return true;
   1.447 -            }
   1.448 -          //the node of g2, which neighbours are the candidates for
   1.449 -          //the pair of order[_depth]
   1.450 -          typename G2::Node currPNode;
   1.451 -          if(currEdgeIts[_depth]==INVALID)
   1.452 -            {
   1.453 -              typename G1::IncEdgeIt fstMatchedE(_g1,order[_depth]);
   1.454 -              //if _mapping[order[_depth]]!=INVALID, we dont use
   1.455 -              //fstMatchedE
   1.456 -              if(_mapping[order[_depth]]==INVALID)
   1.457 -                for(; fstMatchedE!=INVALID &&
   1.458 -                      _mapping[_g1.oppositeNode(order[_depth],
   1.459 -                                              fstMatchedE)]==INVALID;
   1.460 -                    ++fstMatchedE) ; //find fstMatchedE
   1.461 -              if(fstMatchedE==INVALID||_mapping[order[_depth]]!=INVALID)
   1.462 -                {
   1.463 -                  //We did not find an covered neighbour, this means
   1.464 -                  //the graph is not connected(or _depth==0).  Every
   1.465 -                  //uncovered(and there are some other properties due
   1.466 -                  //to the spec. problem types) node of g2 is
   1.467 -                  //candidate.  We can read the iterator of the last
   1.468 -                  //tryed node from the match if it is not the first
   1.469 -                  //try(match[order[_depth]]!=INVALID)
   1.470 -                  typename G2::NodeIt n2(_g2);
   1.471 -                  //if its not the first try
   1.472 -                  if(_mapping[order[_depth]]!=INVALID)
   1.473 -                    {
   1.474 -                      n2=++typename G2::NodeIt(_g2,_mapping[order[_depth]]);
   1.475 -                      subPair(order[_depth],_mapping[order[_depth]]);
   1.476 -                    }
   1.477 -                  for(; n2!=INVALID; ++n2)
   1.478 -                    if(MT!=SUBGRAPH&&_conn[n2]==0)
   1.479 -                      {
   1.480 -                        if(feas<MT>(order[_depth],n2))
   1.481 -                          break;
   1.482 -                      }
   1.483 -                    else if(MT==SUBGRAPH&&_conn[n2]>=0)
   1.484 -                      if(feas<MT>(order[_depth],n2))
   1.485 -                        break;
   1.486 -                  // n2 is the next candidate
   1.487 -                  if(n2!=INVALID)
   1.488 -                    {
   1.489 -                      addPair(order[_depth],n2);
   1.490 -                      ++_depth;
   1.491 -                    }
   1.492 -                  else // there is no more candidate
   1.493 -                    --_depth;
   1.494 -                  continue;
   1.495 -                }
   1.496 -              else
   1.497 -                {
   1.498 -                  currPNode=_mapping[_g1.oppositeNode(order[_depth],
   1.499 -                                                      fstMatchedE)];
   1.500 -                  currEdgeIts[_depth]=typename G2::IncEdgeIt(_g2,currPNode);
   1.501 -                }
   1.502 -            }
   1.503 -          else
   1.504 -            {
   1.505 -              currPNode=_g2.oppositeNode(_mapping[order[_depth]],
   1.506 -                                         currEdgeIts[_depth]);
   1.507 -              subPair(order[_depth],_mapping[order[_depth]]);
   1.508 -              ++currEdgeIts[_depth];
   1.509 -            }
   1.510 -          for(; currEdgeIts[_depth]!=INVALID; ++currEdgeIts[_depth])
   1.511 -            {
   1.512 -              const typename G2::Node currNode =
   1.513 -                _g2.oppositeNode(currPNode, currEdgeIts[_depth]);
   1.514 -              //if currNode is uncovered
   1.515 -              if(_conn[currNode]>0&&feas<MT>(order[_depth],currNode))
   1.516 -                {
   1.517 -                  addPair(order[_depth],currNode);
   1.518 -                  break;
   1.519 -                }
   1.520 -            }
   1.521 -          currEdgeIts[_depth]==INVALID?--_depth:++_depth;
   1.522 +            continue;
   1.523 +          }
   1.524 +          else {
   1.525 +            currPNode=_mapping[_g1.oppositeNode(nodeOfDepth,
   1.526 +                                                fstMatchedE)];
   1.527 +            edgeItOfDepth=typename G2::IncEdgeIt(_g2,currPNode);
   1.528 +          }
   1.529          }
   1.530 +        else {
   1.531 +          currPNode=_g2.oppositeNode(pairOfNodeOfDepth,
   1.532 +                                     edgeItOfDepth);
   1.533 +          subPair(nodeOfDepth,pairOfNodeOfDepth);
   1.534 +          ++edgeItOfDepth;
   1.535 +        }
   1.536 +        for(; edgeItOfDepth!=INVALID; ++edgeItOfDepth) {
   1.537 +          const typename G2::Node currNode =
   1.538 +            _g2.oppositeNode(currPNode, edgeItOfDepth);
   1.539 +          if(_conn[currNode]>0&&feas<MT>(nodeOfDepth,currNode)) {
   1.540 +            addPair(nodeOfDepth,currNode);
   1.541 +            break;
   1.542 +          }
   1.543 +        }
   1.544 +        edgeItOfDepth==INVALID?--_depth:++_depth;
   1.545 +      }
   1.546        return false;
   1.547      }
   1.548  
   1.549 -    //calc. the lookup table for cut the searchtree
   1.550 -    void setRNew1tRInOut1t()
   1.551 -    {
   1.552 +    //calculate the lookup table for cutting the search tree
   1.553 +    void initRNew1tRInOut1t() {
   1.554        typename G1::template NodeMap<int> tmp(_g1,0);
   1.555 -      for(unsigned int i=0; i<order.size(); ++i)
   1.556 -        {
   1.557 -          tmp[order[i]]=-1;
   1.558 -          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
   1.559 -            {
   1.560 -              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
   1.561 -              if(tmp[currNode]>0)
   1.562 -                ++rInOut1t[order[i]];
   1.563 -              else if(tmp[currNode]==0)
   1.564 -                ++rNew1t[order[i]];
   1.565 -            }
   1.566 -          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
   1.567 -            {
   1.568 -              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
   1.569 -              if(tmp[currNode]!=-1)
   1.570 -                ++tmp[currNode];
   1.571 -            }
   1.572 +      for(unsigned int i=0; i<_order.size(); ++i) {
   1.573 +        const typename G1::Node& orderI = _order[i];
   1.574 +        tmp[orderI]=-1;
   1.575 +        for(typename G1::IncEdgeIt e1(_g1,orderI); e1!=INVALID; ++e1) {
   1.576 +          const typename G1::Node currNode=_g1.oppositeNode(orderI,e1);
   1.577 +          if(tmp[currNode]>0)
   1.578 +            ++_rInOut1t[orderI];
   1.579 +          else if(tmp[currNode]==0)
   1.580 +            ++_rNew1t[orderI];
   1.581          }
   1.582 +        for(typename G1::IncEdgeIt e1(_g1,orderI); e1!=INVALID; ++e1) {
   1.583 +          const typename G1::Node currNode=_g1.oppositeNode(orderI,e1);
   1.584 +          if(tmp[currNode]!=-1)
   1.585 +            ++tmp[currNode];
   1.586 +        }
   1.587 +      }
   1.588      }
   1.589    public:
   1.590      ///Constructor
   1.591 @@ -399,62 +333,73 @@
   1.592      ///\param g2 The graph \e g1 will be embedded into.
   1.593      ///\param m \ref concepts::ReadWriteMap "read-write" NodeMap
   1.594      ///storing the found mapping.
   1.595 -    ///\param neq A bool-valued binary functor determinining whether a node is
   1.596 +    ///\param neq A bool-valued binary functor determining whether a node is
   1.597      ///mappable to another. By default it is an always true operator.
   1.598 -    Vf2(const G1 &g1, const G2 &g2,M &m, const NEQ &neq = NEQ() ) :
   1.599 -      _nEq(neq), _conn(g2,0), _mapping(m), order(countNodes(g1)),
   1.600 -      currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNew1t(g1,0),
   1.601 -      rInOut1t(g1,0), _mapping_type(SUBGRAPH)
   1.602 +    Vf2(const G1 &g1, const G2 &g2, M &m, const NEQ &neq = NEQ() ) :
   1.603 +      _g1(g1), _g2(g2), _nEq(neq), _depth(0), _mapping(m),
   1.604 +      _order(countNodes(g1)), _conn(g2,0),
   1.605 +      _currEdgeIts(countNodes(g1),INVALID), _rNew1t(g1,0), _rInOut1t(g1,0),
   1.606 +      _mapping_type(SUBGRAPH), _deallocMappingAfterUse(0)
   1.607      {
   1.608 -      _depth=0;
   1.609 -      setOrder();
   1.610 -      setRNew1tRInOut1t();
   1.611 +      initOrder();
   1.612 +      initRNew1tRInOut1t();
   1.613        for(typename G1::NodeIt n(g1);n!=INVALID;++n)
   1.614          m[n]=INVALID;
   1.615      }
   1.616  
   1.617 +    ///Destructor
   1.618 +
   1.619 +    ///Destructor.
   1.620 +    ///
   1.621 +
   1.622 +    ~Vf2(){
   1.623 +      if(_deallocMappingAfterUse)
   1.624 +        delete &_mapping;
   1.625 +    }
   1.626 +
   1.627      ///Returns the current mapping type
   1.628  
   1.629      ///Returns the current mapping type
   1.630      ///
   1.631 -    Vf2MappingType mappingType() const { return _mapping_type; }
   1.632 +    MappingType mappingType() const {
   1.633 +      return _mapping_type;
   1.634 +    }
   1.635      ///Sets mapping type
   1.636  
   1.637      ///Sets mapping type.
   1.638      ///
   1.639      ///The mapping type is set to \ref SUBGRAPH by default.
   1.640      ///
   1.641 -    ///\sa See \ref Vf2MappingType for the possible values.
   1.642 -    void mappingType(Vf2MappingType m_type) { _mapping_type = m_type; }
   1.643 +    ///\sa See \ref MappingType for the possible values.
   1.644 +    void mappingType(MappingType m_type) {
   1.645 +      _mapping_type = m_type;
   1.646 +    }
   1.647  
   1.648      ///Finds a mapping
   1.649  
   1.650 -    ///It finds a mapping between from g1 into g2 according to the mapping
   1.651 -    ///type set by \ref mappingType(Vf2MappingType) "mappingType()".
   1.652 +    ///It finds a mapping from g1 into g2 according to the mapping
   1.653 +    ///type set by \ref mappingType(MappingType) "mappingType()".
   1.654      ///
   1.655      ///By subsequent calls, it returns all possible mappings one-by-one.
   1.656      ///
   1.657      ///\retval true if a mapping is found.
   1.658      ///\retval false if there is no (more) mapping.
   1.659 -    bool find()
   1.660 -    {
   1.661 -      switch(_mapping_type)
   1.662 -        {
   1.663 -        case SUBGRAPH:
   1.664 -          return extMatch<SUBGRAPH>();
   1.665 -        case INDUCED:
   1.666 -          return extMatch<INDUCED>();
   1.667 -        case ISOMORPH:
   1.668 -          return extMatch<ISOMORPH>();
   1.669 -        default:
   1.670 -          return false;
   1.671 -        }
   1.672 +    bool find() {
   1.673 +      switch(_mapping_type) {
   1.674 +      case SUBGRAPH:
   1.675 +        return extMatch<SUBGRAPH>();
   1.676 +      case INDUCED:
   1.677 +        return extMatch<INDUCED>();
   1.678 +      case ISOMORPH:
   1.679 +        return extMatch<ISOMORPH>();
   1.680 +      default:
   1.681 +        return false;
   1.682 +      }
   1.683      }
   1.684    };
   1.685  
   1.686    template<class G1, class G2>
   1.687 -  class Vf2WizardBase
   1.688 -  {
   1.689 +  class Vf2WizardBase {
   1.690    protected:
   1.691      typedef G1 Graph1;
   1.692      typedef G2 Graph2;
   1.693 @@ -462,35 +407,37 @@
   1.694      const G1 &_g1;
   1.695      const G2 &_g2;
   1.696  
   1.697 -    Vf2MappingType _mapping_type;
   1.698 +    MappingType _mapping_type;
   1.699  
   1.700      typedef typename G1::template NodeMap<typename G2::Node> Mapping;
   1.701      bool _local_mapping;
   1.702      void *_mapping;
   1.703 -    void createMapping()
   1.704 -    {
   1.705 +    void createMapping() {
   1.706        _mapping = new Mapping(_g1);
   1.707      }
   1.708  
   1.709 +    void *myVf2; //used in Vf2Wizard::find
   1.710 +
   1.711 +
   1.712      typedef bits::vf2::AlwaysEq NodeEq;
   1.713      NodeEq _node_eq;
   1.714  
   1.715      Vf2WizardBase(const G1 &g1,const G2 &g2)
   1.716 -      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) {}
   1.717 +      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) { }
   1.718    };
   1.719  
   1.720 +
   1.721    /// Auxiliary class for the function-type interface of %VF2 algorithm.
   1.722 -
   1.723 +  ///
   1.724    /// This auxiliary class implements the named parameters of
   1.725    /// \ref vf2() "function-type interface" of \ref Vf2 algorithm.
   1.726    ///
   1.727 -  /// \warning This class should only be used through the function \ref vf2().
   1.728 +  /// \warning This class is not to be used directly.
   1.729    ///
   1.730    /// \tparam TR The traits class that defines various types used by the
   1.731    /// algorithm.
   1.732    template<class TR>
   1.733 -  class Vf2Wizard : public TR
   1.734 -  {
   1.735 +  class Vf2Wizard : public TR {
   1.736      typedef TR Base;
   1.737      typedef typename TR::Graph1 Graph1;
   1.738      typedef typename TR::Graph2 Graph2;
   1.739 @@ -511,9 +458,12 @@
   1.740      ///Copy constructor
   1.741      Vf2Wizard(const Base &b) : Base(b) {}
   1.742  
   1.743 +    ///Copy constructor
   1.744 +    Vf2Wizard(const Vf2Wizard &b) : Base(b) {}
   1.745 +
   1.746  
   1.747      template<class T>
   1.748 -    struct SetMappingBase : public Base {
   1.749 +    struct SetMappingBase : public Base{
   1.750        typedef T Mapping;
   1.751        SetMappingBase(const Base &b) : Base(b) {}
   1.752      };
   1.753 @@ -524,8 +474,7 @@
   1.754      ///\ref named-templ-param "Named parameter" function for setting
   1.755      ///the map that stores the found embedding.
   1.756      template<class T>
   1.757 -    Vf2Wizard< SetMappingBase<T> > mapping(const T &t)
   1.758 -    {
   1.759 +    Vf2Wizard< SetMappingBase<T> > mapping(const T &t) {
   1.760        Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
   1.761        Base::_local_mapping = false;
   1.762        return Vf2Wizard<SetMappingBase<T> >(*this);
   1.763 @@ -536,7 +485,8 @@
   1.764        typedef NE NodeEq;
   1.765        NodeEq _node_eq;
   1.766        SetNodeEqBase(const Base &b, const NE &node_eq)
   1.767 -        : Base(b), _node_eq(node_eq) {}
   1.768 +        : Base(b), _node_eq(node_eq){
   1.769 +      }
   1.770      };
   1.771  
   1.772      ///\brief \ref named-templ-param "Named parameter" for setting
   1.773 @@ -549,8 +499,7 @@
   1.774      ///whether a node is mappable to another. By default it is an
   1.775      ///always true operator.
   1.776      template<class T>
   1.777 -    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq)
   1.778 -    {
   1.779 +    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq) {
   1.780        return Vf2Wizard<SetNodeEqBase<T> >(SetNodeEqBase<T>(*this,node_eq));
   1.781      }
   1.782  
   1.783 @@ -560,16 +509,15 @@
   1.784      ///\ref named-templ-param "Named parameter" function for setting
   1.785      ///the node labels defining equivalence relation between them.
   1.786      ///
   1.787 -    ///\param m1 It is arbitrary \ref concepts::ReadMap "readable node map"
   1.788 +    ///\param m1 An arbitrary \ref concepts::ReadMap "readable node map"
   1.789      ///of g1.
   1.790 -    ///\param m2 It is arbitrary \ref concepts::ReadMap "readable node map"
   1.791 +    ///\param m2 An arbitrary \ref concepts::ReadMap "readable node map"
   1.792      ///of g2.
   1.793      ///
   1.794      ///The value type of these maps must be equal comparable.
   1.795      template<class M1, class M2>
   1.796      Vf2Wizard< SetNodeEqBase<bits::vf2::MapEq<M1,M2> > >
   1.797 -    nodeLabels(const M1 &m1,const M2 &m2)
   1.798 -    {
   1.799 +    nodeLabels(const M1 &m1,const M2 &m2){
   1.800        return nodeEq(bits::vf2::MapEq<M1,M2>(m1,m2));
   1.801      }
   1.802  
   1.803 @@ -581,9 +529,8 @@
   1.804      ///
   1.805      ///The mapping type is set to \ref SUBGRAPH by default.
   1.806      ///
   1.807 -    ///\sa See \ref Vf2MappingType for the possible values.
   1.808 -    Vf2Wizard<Base> &mappingType(Vf2MappingType m_type)
   1.809 -    {
   1.810 +    ///\sa See \ref MappingType for the possible values.
   1.811 +    Vf2Wizard<Base> &mappingType(MappingType m_type) {
   1.812        _mapping_type = m_type;
   1.813        return *this;
   1.814      }
   1.815 @@ -593,8 +540,7 @@
   1.816      ///
   1.817      ///\ref named-templ-param "Named parameter" for setting
   1.818      ///the mapping type to \ref INDUCED.
   1.819 -    Vf2Wizard<Base> &induced()
   1.820 -    {
   1.821 +    Vf2Wizard<Base> &induced() {
   1.822        _mapping_type = INDUCED;
   1.823        return *this;
   1.824      }
   1.825 @@ -604,20 +550,19 @@
   1.826      ///
   1.827      ///\ref named-templ-param "Named parameter" for setting
   1.828      ///the mapping type to \ref ISOMORPH.
   1.829 -    Vf2Wizard<Base> &iso()
   1.830 -    {
   1.831 +    Vf2Wizard<Base> &iso() {
   1.832        _mapping_type = ISOMORPH;
   1.833        return *this;
   1.834      }
   1.835  
   1.836 +
   1.837      ///Runs VF2 algorithm.
   1.838  
   1.839      ///This method runs VF2 algorithm.
   1.840      ///
   1.841      ///\retval true if a mapping is found.
   1.842 -    ///\retval false if there is no (more) mapping.
   1.843 -    bool run()
   1.844 -    {
   1.845 +    ///\retval false if there is no mapping.
   1.846 +    bool run(){
   1.847        if(Base::_local_mapping)
   1.848          Base::createMapping();
   1.849  
   1.850 @@ -633,6 +578,46 @@
   1.851  
   1.852        return ret;
   1.853      }
   1.854 +
   1.855 +    ///Get a pointer to the generated Vf2 object.
   1.856 +
   1.857 +    ///Gives a pointer to the generated Vf2 object.
   1.858 +    ///
   1.859 +    ///\return Pointer to the generated Vf2 object.
   1.860 +    ///\warning Don't forget to delete the referred Vf2 object after use.
   1.861 +    Vf2<Graph1, Graph2, Mapping, NodeEq >* getPtrToVf2Object() {
   1.862 +      if(Base::_local_mapping)
   1.863 +        Base::createMapping();
   1.864 +      Vf2<Graph1, Graph2, Mapping, NodeEq >* ptr =
   1.865 +        new Vf2<Graph1, Graph2, Mapping, NodeEq>
   1.866 +        (_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
   1.867 +      ptr->mappingType(_mapping_type);
   1.868 +      if(Base::_local_mapping)
   1.869 +        ptr->_deallocMappingAfterUse = true;
   1.870 +      return ptr;
   1.871 +    }
   1.872 +
   1.873 +    ///Counts the number of mappings.
   1.874 +
   1.875 +    ///This method counts the number of mappings.
   1.876 +    ///
   1.877 +    /// \return The number of mappings.
   1.878 +    int count() {
   1.879 +      if(Base::_local_mapping)
   1.880 +        Base::createMapping();
   1.881 +
   1.882 +      Vf2<Graph1, Graph2, Mapping, NodeEq>
   1.883 +        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
   1.884 +      if(Base::_local_mapping)
   1.885 +        alg._deallocMappingAfterUse = true;
   1.886 +      alg.mappingType(_mapping_type);
   1.887 +
   1.888 +      int ret = 0;
   1.889 +      while(alg.find())
   1.890 +        ++ret;
   1.891 +
   1.892 +      return ret;
   1.893 +    }
   1.894    };
   1.895  
   1.896    ///Function-type interface for VF2 algorithm.
   1.897 @@ -644,20 +629,32 @@
   1.898    ///declared as the members of class \ref Vf2Wizard.
   1.899    ///The following examples show how to use these parameters.
   1.900    ///\code
   1.901 -  ///  // Find an embedding of graph g into graph h
   1.902 +  ///  // Find an embedding of graph g1 into graph g2
   1.903    ///  ListGraph::NodeMap<ListGraph::Node> m(g);
   1.904 -  ///  vf2(g,h).mapping(m).run();
   1.905 +  ///  vf2(g1,g2).mapping(m).run();
   1.906    ///
   1.907 -  ///  // Check whether graphs g and h are isomorphic
   1.908 -  ///  bool is_iso = vf2(g,h).iso().run();
   1.909 +  ///  // Check whether graphs g1 and g2 are isomorphic
   1.910 +  ///  bool is_iso = vf2(g1,g2).iso().run();
   1.911 +  ///
   1.912 +  ///  // Count the number of isomorphisms
   1.913 +  ///  int num_isos = vf2(g1,g2).iso().count();
   1.914 +  ///
   1.915 +  ///  // Iterate through all the induced subgraph mappings of graph g1 into g2
   1.916 +  ///  auto* myVf2 = vf2(g1,g2).mapping(m).nodeLabels(c1,c2)
   1.917 +  ///  .induced().getPtrToVf2Object();
   1.918 +  ///  while(myVf2->find()){
   1.919 +  ///    //process the current mapping m
   1.920 +  ///  }
   1.921 +  ///  delete myVf22;
   1.922    ///\endcode
   1.923 -  ///\warning Don't forget to put the \ref Vf2Wizard::run() "run()"
   1.924 +  ///\warning Don't forget to put the \ref Vf2Wizard::run() "run()",
   1.925 +  ///\ref Vf2Wizard::count() "count()" or
   1.926 +  ///the \ref Vf2Wizard::getPtrToVf2Object() "getPtrToVf2Object()"
   1.927    ///to the end of the expression.
   1.928    ///\sa Vf2Wizard
   1.929    ///\sa Vf2
   1.930    template<class G1, class G2>
   1.931 -  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2)
   1.932 -  {
   1.933 +  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2) {
   1.934      return Vf2Wizard<Vf2WizardBase<G1,G2> >(g1,g2);
   1.935    }
   1.936