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