lemon/vf2pp.h
changeset 1415 959d78f3fe0e
parent 1413 e68f0ef37e77
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/vf2pp.h	Wed Oct 17 22:56:43 2018 +0200
     1.3 @@ -0,0 +1,860 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     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 + * provided that this copyright notice appears in all copies. For
    1.13 + * precise terms see the accompanying LICENSE file.
    1.14 + *
    1.15 + * This software is provided "AS IS" with no warranty of any kind,
    1.16 + * express or implied, and with no claim as to its suitability for any
    1.17 + * purpose.
    1.18 + *
    1.19 + */
    1.20 +
    1.21 +#ifndef LEMON_VF2PP_H
    1.22 +#define LEMON_VF2PP_H
    1.23 +
    1.24 +///\ingroup graph_properties
    1.25 +///\file
    1.26 +///\brief VF2 Plus Plus algorithm.
    1.27 +
    1.28 +#include <lemon/core.h>
    1.29 +#include <lemon/concepts/graph.h>
    1.30 +#include <lemon/bits/vf2_internals.h>
    1.31 +
    1.32 +#include <vector>
    1.33 +#include <algorithm>
    1.34 +#include <utility>
    1.35 +
    1.36 +namespace lemon {
    1.37 +
    1.38 +  ///%VF2 Plus Plus algorithm class \cite VF2PP.
    1.39 +
    1.40 +  ///\ingroup graph_isomorphism This class provides an efficient
    1.41 +  ///implementation of the %VF2 Plus Plus algorithm \cite VF2PP
    1.42 +  ///for variants of the (Sub)graph Isomorphism problem.
    1.43 +  ///
    1.44 +  ///There is also a \ref vf2pp() "function-type interface" called
    1.45 +  ///\ref vf2pp() for the %VF2 Plus Plus algorithm, which is probably
    1.46 +  ///more convenient in most use cases.
    1.47 +  ///
    1.48 +  ///\tparam G1 The type of the graph to be embedded.
    1.49 +  ///The default type is \ref ListGraph.
    1.50 +  ///\tparam G2 The type of the graph g1 will be embedded into.
    1.51 +  ///The default type is \ref ListGraph.
    1.52 +  ///\tparam M The type of the NodeMap storing the mapping.
    1.53 +  ///By default, it is G1::NodeMap<G2::Node>
    1.54 +  ///\tparam M1 The type of the NodeMap storing the integer node labels of G1.
    1.55 +  ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
    1.56 +  ///different labels. By default, it is G1::NodeMap<int>.
    1.57 +  ///\tparam M2 The type of the NodeMap storing the integer node labels of G2.
    1.58 +  ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
    1.59 +  ///different labels. By default, it is G2::NodeMap<int>.
    1.60 +  ///
    1.61 +  ///\sa vf2pp()
    1.62 +#ifdef DOXYGEN
    1.63 +  template<class G1, class G2, class M, class M1, class M2 >
    1.64 +#else
    1.65 +  template<class G1 = ListGraph,
    1.66 +           class G2 = ListGraph,
    1.67 +           class M = typename G1::template NodeMap<G2::Node>,
    1.68 +           class M1 = typename G1::template NodeMap<int>,
    1.69 +           class M2 = typename G2::template NodeMap<int> >
    1.70 +#endif
    1.71 +  class Vf2pp {
    1.72 +    //The graph to be embedded
    1.73 +    const G1 &_g1;
    1.74 +
    1.75 +    //The graph into which g1 is to be embedded
    1.76 +    const G2 &_g2;
    1.77 +
    1.78 +    //Current depth in the search tree
    1.79 +    int _depth;
    1.80 +
    1.81 +    //The current mapping. _mapping[v1]=v2 iff v1 has been mapped to v2,
    1.82 +    //where v1 is a node of G1 and v2 is a node of G2
    1.83 +    M &_mapping;
    1.84 +
    1.85 +    //_order[i] is a node of g1 for which a pair is searched if depth=i
    1.86 +    std::vector<typename G1::Node> _order;
    1.87 +
    1.88 +    //_conn[v2] = number of covered neighbours of v2
    1.89 +    typename G2::template NodeMap<int> _conn;
    1.90 +
    1.91 +    //_currEdgeIts[i] is the last used edge iterator in the i-th
    1.92 +    //depth to find a pair for node _order[i]
    1.93 +    std::vector<typename G2::IncEdgeIt> _currEdgeIts;
    1.94 + 
    1.95 +    //_rNewLabels1[v] is a pair of form
    1.96 +    //(label; num. of uncov. nodes with such label and no covered neighbours)
    1.97 +    typename G1::template NodeMap<std::vector<std::pair<int,int> > >
    1.98 +    _rNewLabels1;
    1.99 +
   1.100 +    //_rInOutLabels1[v] is the number of covered neighbours of v for each label
   1.101 +    //in form (label,number of such labels)
   1.102 +    typename G1::template NodeMap<std::vector<std::pair<int,int> > >
   1.103 +    _rInOutLabels1;
   1.104 +
   1.105 +    //_intLabels1[v]==i means that node v has label i in _g1
   1.106 +    //(i is in {0,1,2,..,K-1}, where K is the number of diff. labels)
   1.107 +    M1 &_intLabels1;
   1.108 +
   1.109 +    //_intLabels2[v]==i means that node v has label i in _g2
   1.110 +    //(i is in {0,1,2,..,K-1}, where K is the number of diff. labels)
   1.111 +    M2 &_intLabels2;
   1.112 +
   1.113 +    //largest label
   1.114 +    const int _maxLabel;
   1.115 +
   1.116 +    //lookup tables for manipulating with label class cardinalities
   1.117 +    //(after use they have to be reset to 0..0)
   1.118 +    std::vector<int> _labelTmp1,_labelTmp2;
   1.119 +
   1.120 +    MappingType _mapping_type;
   1.121 +
   1.122 +    //indicates whether the mapping or the labels must be deleted in the destructor
   1.123 +    bool _deallocMappingAfterUse,_deallocLabelsAfterUse;
   1.124 +
   1.125 +
   1.126 +    //improved cutting function
   1.127 +    template<MappingType MT>
   1.128 +    bool cutByLabels(const typename G1::Node n1,const typename G2::Node n2) {
   1.129 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.130 +        const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.131 +        if(_conn[currNode]>0)
   1.132 +          --_labelTmp1[_intLabels2[currNode]];
   1.133 +        else if(MT!=SUBGRAPH&&_conn[currNode]==0)
   1.134 +          --_labelTmp2[_intLabels2[currNode]];
   1.135 +      }
   1.136 +
   1.137 +      bool ret=1;
   1.138 +      if(ret) {
   1.139 +        for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
   1.140 +          _labelTmp1[_rInOutLabels1[n1][i].first]+=_rInOutLabels1[n1][i].second;
   1.141 +
   1.142 +        if(MT!=SUBGRAPH)
   1.143 +          for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
   1.144 +            _labelTmp2[_rNewLabels1[n1][i].first]+=_rNewLabels1[n1][i].second;
   1.145 +
   1.146 +        switch(MT) {
   1.147 +        case INDUCED:
   1.148 +          for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
   1.149 +            if(_labelTmp1[_rInOutLabels1[n1][i].first]>0) {
   1.150 +              ret=0;
   1.151 +              break;
   1.152 +            }
   1.153 +          if(ret)
   1.154 +            for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
   1.155 +              if(_labelTmp2[_rNewLabels1[n1][i].first]>0) {
   1.156 +                ret=0;
   1.157 +                break;
   1.158 +              }
   1.159 +          break;
   1.160 +        case SUBGRAPH:
   1.161 +          for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
   1.162 +            if(_labelTmp1[_rInOutLabels1[n1][i].first]>0) {
   1.163 +              ret=0;
   1.164 +              break;
   1.165 +            }
   1.166 +          break;
   1.167 +        case ISOMORPH:
   1.168 +          for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
   1.169 +            if(_labelTmp1[_rInOutLabels1[n1][i].first]!=0) {
   1.170 +              ret=0;
   1.171 +              break;
   1.172 +            }
   1.173 +          if(ret)
   1.174 +            for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
   1.175 +              if(_labelTmp2[_rNewLabels1[n1][i].first]!=0) {
   1.176 +                ret=0;
   1.177 +                break;
   1.178 +              }
   1.179 +          break;
   1.180 +        default:
   1.181 +          return false;
   1.182 +        }
   1.183 +        for(unsigned int i = 0; i < _rInOutLabels1[n1].size(); ++i)
   1.184 +          _labelTmp1[_rInOutLabels1[n1][i].first]=0;
   1.185 +
   1.186 +        if(MT!=SUBGRAPH)
   1.187 +          for(unsigned int i = 0; i < _rNewLabels1[n1].size(); ++i)
   1.188 +            _labelTmp2[_rNewLabels1[n1][i].first]=0;
   1.189 +      }
   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 +        _labelTmp1[_intLabels2[currNode]]=0;
   1.194 +        if(MT!=SUBGRAPH&&_conn[currNode]==0)
   1.195 +          _labelTmp2[_intLabels2[currNode]]=0;
   1.196 +      }
   1.197 +
   1.198 +      return ret;
   1.199 +    }
   1.200 +
   1.201 +
   1.202 +    //try to exclude the matching of n1 and n2
   1.203 +    template<MappingType MT>
   1.204 +    bool feas(const typename G1::Node n1,const typename G2::Node n2) {
   1.205 +      if(_intLabels1[n1]!=_intLabels2[n2])
   1.206 +        return 0;
   1.207 +
   1.208 +      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
   1.209 +        const typename G1::Node& currNode=_g1.oppositeNode(n1,e1);
   1.210 +        if(_mapping[currNode]!=INVALID)
   1.211 +          --_conn[_mapping[currNode]];
   1.212 +      }
   1.213 +
   1.214 +      bool isIso=1;
   1.215 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.216 +        int& connCurrNode = _conn[_g2.oppositeNode(n2,e2)];
   1.217 +        if(connCurrNode<-1)
   1.218 +          ++connCurrNode;
   1.219 +        else if(MT!=SUBGRAPH&&connCurrNode==-1) {
   1.220 +          isIso=0;
   1.221 +          break;
   1.222 +        }
   1.223 +      }
   1.224 +
   1.225 +      if(isIso)
   1.226 +        for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
   1.227 +          const typename G2::Node& currNodePair =
   1.228 +            _mapping[_g1.oppositeNode(n1,e1)];
   1.229 +          int& connCurrNodePair=_conn[currNodePair];
   1.230 +          if(currNodePair!=INVALID&&connCurrNodePair!=-1) {
   1.231 +            switch(MT){
   1.232 +            case INDUCED:
   1.233 +            case ISOMORPH:
   1.234 +              isIso=0;
   1.235 +              break;
   1.236 +            case SUBGRAPH:
   1.237 +              if(connCurrNodePair<-1)
   1.238 +                isIso=0;
   1.239 +              break;
   1.240 +            }
   1.241 +            connCurrNodePair=-1;
   1.242 +          }
   1.243 +        }
   1.244 +      else
   1.245 +        for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1) {
   1.246 +          const typename G2::Node currNode=_mapping[_g1.oppositeNode(n1,e1)];
   1.247 +          if(currNode!=INVALID/*&&_conn[currNode]!=-1*/)
   1.248 +            _conn[currNode]=-1;
   1.249 +        }
   1.250 +
   1.251 +      return isIso&&cutByLabels<MT>(n1,n2);
   1.252 +    }
   1.253 +
   1.254 +    //maps n1 to n2
   1.255 +    void addPair(const typename G1::Node n1,const typename G2::Node n2) {
   1.256 +      _conn[n2]=-1;
   1.257 +      _mapping.set(n1,n2);
   1.258 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2) {
   1.259 +        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
   1.260 +        if(currConn!=-1)
   1.261 +          ++currConn;
   1.262 +      }
   1.263 +    }
   1.264 +
   1.265 +    //removes mapping of n1 to n2
   1.266 +    void subPair(const typename G1::Node n1,const typename G2::Node n2) {
   1.267 +      _conn[n2]=0;
   1.268 +      _mapping.set(n1,INVALID);
   1.269 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2){
   1.270 +        int& currConn = _conn[_g2.oppositeNode(n2,e2)];
   1.271 +        if(currConn>0)
   1.272 +          --currConn;
   1.273 +        else if(currConn==-1)
   1.274 +          ++_conn[n2];
   1.275 +      }
   1.276 +    }
   1.277 +
   1.278 +    void processBfsTree(typename G1::Node source,unsigned int& orderIndex,
   1.279 +                         typename G1::template NodeMap<int>& dm1,
   1.280 +                         typename G1::template NodeMap<bool>& added) {
   1.281 +      _order[orderIndex]=source;
   1.282 +      added[source]=1;
   1.283 +
   1.284 +      unsigned int endPosOfLevel=orderIndex,
   1.285 +        startPosOfLevel=orderIndex,
   1.286 +        lastAdded=orderIndex;
   1.287 +
   1.288 +      typename G1::template NodeMap<int> currConn(_g1,0);
   1.289 +
   1.290 +      while(orderIndex<=lastAdded){
   1.291 +        typename G1::Node currNode = _order[orderIndex];
   1.292 +        for(typename G1::IncEdgeIt e(_g1,currNode); e!=INVALID; ++e) {
   1.293 +          typename G1::Node n = _g1.oppositeNode(currNode,e);
   1.294 +          if(!added[n]) {
   1.295 +            _order[++lastAdded]=n;
   1.296 +            added[n]=1;
   1.297 +          }
   1.298 +        }
   1.299 +        if(orderIndex>endPosOfLevel){
   1.300 +          for(unsigned int j = startPosOfLevel; j <= endPosOfLevel; ++j) {
   1.301 +            int minInd=j;
   1.302 +            for(unsigned int i = j+1; i <= endPosOfLevel; ++i)
   1.303 +              if(currConn[_order[i]]>currConn[_order[minInd]]||
   1.304 +                 (currConn[_order[i]]==currConn[_order[minInd]]&&
   1.305 +                  (dm1[_order[i]]>dm1[_order[minInd]]||
   1.306 +                   (dm1[_order[i]]==dm1[_order[minInd]]&&
   1.307 +                    _labelTmp1[_intLabels1[_order[minInd]]]>
   1.308 +                    _labelTmp1[_intLabels1[_order[i]]]))))
   1.309 +                minInd=i;
   1.310 +
   1.311 +            --_labelTmp1[_intLabels1[_order[minInd]]];
   1.312 +            for(typename G1::IncEdgeIt e(_g1,_order[minInd]); e!=INVALID; ++e)
   1.313 +              ++currConn[_g1.oppositeNode(_order[minInd],e)];
   1.314 +            std::swap(_order[j],_order[minInd]);
   1.315 +          }
   1.316 +          startPosOfLevel=endPosOfLevel+1;
   1.317 +          endPosOfLevel=lastAdded;
   1.318 +        }
   1.319 +        ++orderIndex;
   1.320 +      }
   1.321 +    }
   1.322 +
   1.323 +
   1.324 +    //we will find pairs for the nodes of g1 in this order
   1.325 +    void initOrder(){
   1.326 +      for(typename G2::NodeIt n2(_g2); n2!=INVALID; ++n2)
   1.327 +        ++_labelTmp1[_intLabels2[n2]];
   1.328 +
   1.329 +      typename G1::template NodeMap<int> dm1(_g1,0);
   1.330 +      for(typename G1::EdgeIt e(_g1); e!=INVALID; ++e) {
   1.331 +        ++dm1[_g1.u(e)];
   1.332 +        ++dm1[_g1.v(e)];
   1.333 +      }
   1.334 +
   1.335 +      typename G1::template NodeMap<bool> added(_g1,0);
   1.336 +      unsigned int orderIndex=0;
   1.337 +
   1.338 +      for(typename G1::NodeIt n(_g1); n!=INVALID;) {
   1.339 +        if(!added[n]){
   1.340 +          typename G1::Node minNode = n;
   1.341 +          for(typename G1::NodeIt n1(_g1,minNode); n1!=INVALID; ++n1)
   1.342 +            if(!added[n1] &&
   1.343 +               (_labelTmp1[_intLabels1[minNode]]>
   1.344 +                _labelTmp1[_intLabels1[n1]]||(dm1[minNode]<dm1[n1]&&
   1.345 +                                             _labelTmp1[_intLabels1[minNode]]==
   1.346 +                                             _labelTmp1[_intLabels1[n1]])))
   1.347 +              minNode=n1;
   1.348 +          processBfsTree(minNode,orderIndex,dm1,added);
   1.349 +        }
   1.350 +        else
   1.351 +          ++n;
   1.352 +      }
   1.353 +      for(unsigned int i = 0; i < _labelTmp1.size(); ++i)
   1.354 +        _labelTmp1[i]=0;
   1.355 +    }
   1.356 +
   1.357 +
   1.358 +    template<MappingType MT>
   1.359 +    bool extMatch(){
   1.360 +      while(_depth>=0) {
   1.361 +        if(_depth==static_cast<int>(_order.size())) {
   1.362 +          //all nodes of g1 are mapped to nodes of g2
   1.363 +          --_depth;
   1.364 +          return true;
   1.365 +        }
   1.366 +        typename G1::Node& nodeOfDepth = _order[_depth];
   1.367 +        const typename G2::Node& pairOfNodeOfDepth = _mapping[nodeOfDepth];
   1.368 +        typename G2::IncEdgeIt &edgeItOfDepth = _currEdgeIts[_depth];
   1.369 +        //the node of g2 whose neighbours are the candidates for
   1.370 +        //the pair of _order[_depth]
   1.371 +        typename G2::Node currPNode;
   1.372 +        if(edgeItOfDepth==INVALID){
   1.373 +          typename G1::IncEdgeIt fstMatchedE(_g1,nodeOfDepth);
   1.374 +          //if _mapping[_order[_depth]]!=INVALID, we don't need fstMatchedE
   1.375 +          if(pairOfNodeOfDepth==INVALID) {
   1.376 +            for(; fstMatchedE!=INVALID &&
   1.377 +                  _mapping[_g1.oppositeNode(nodeOfDepth,
   1.378 +                                            fstMatchedE)]==INVALID;
   1.379 +                ++fstMatchedE); //find fstMatchedE, it could be preprocessed
   1.380 +          }
   1.381 +          if(fstMatchedE==INVALID||pairOfNodeOfDepth!=INVALID) {
   1.382 +            //We found no covered neighbours, this means that
   1.383 +            //the graph is not connected (or _depth==0). Each
   1.384 +            //uncovered (and there are some other properties due
   1.385 +            //to the spec. problem types) node of g2 is
   1.386 +            //candidate. We can read the iterator of the last
   1.387 +            //tried node from the match if it is not the first
   1.388 +            //try (match[nodeOfDepth]!=INVALID)
   1.389 +            typename G2::NodeIt n2(_g2);
   1.390 +            //if it's not the first try
   1.391 +            if(pairOfNodeOfDepth!=INVALID) {
   1.392 +              n2=++typename G2::NodeIt(_g2,pairOfNodeOfDepth);
   1.393 +              subPair(nodeOfDepth,pairOfNodeOfDepth);
   1.394 +            }
   1.395 +            for(; n2!=INVALID; ++n2)
   1.396 +              if(MT!=SUBGRAPH) {
   1.397 +                if(_conn[n2]==0&&feas<MT>(nodeOfDepth,n2))
   1.398 +                  break;
   1.399 +              }
   1.400 +              else if(_conn[n2]>=0&&feas<MT>(nodeOfDepth,n2))
   1.401 +                break;
   1.402 +            // n2 is the next candidate
   1.403 +            if(n2!=INVALID) {
   1.404 +              addPair(nodeOfDepth,n2);
   1.405 +              ++_depth;
   1.406 +            }
   1.407 +            else // there are no more candidates
   1.408 +              --_depth;
   1.409 +            continue;
   1.410 +          }
   1.411 +          else {
   1.412 +            currPNode=_mapping[_g1.oppositeNode(nodeOfDepth,
   1.413 +                                                fstMatchedE)];
   1.414 +            edgeItOfDepth=typename G2::IncEdgeIt(_g2,currPNode);
   1.415 +          }
   1.416 +        }
   1.417 +        else {
   1.418 +          currPNode=_g2.oppositeNode(pairOfNodeOfDepth,
   1.419 +                                     edgeItOfDepth);
   1.420 +          subPair(nodeOfDepth,pairOfNodeOfDepth);
   1.421 +          ++edgeItOfDepth;
   1.422 +        }
   1.423 +        for(; edgeItOfDepth!=INVALID; ++edgeItOfDepth) {
   1.424 +          const typename G2::Node currNode =
   1.425 +            _g2.oppositeNode(currPNode, edgeItOfDepth);
   1.426 +          if(_conn[currNode]>0&&feas<MT>(nodeOfDepth,currNode)) {
   1.427 +            addPair(nodeOfDepth,currNode);
   1.428 +            break;
   1.429 +          }
   1.430 +        }
   1.431 +        edgeItOfDepth==INVALID?--_depth:++_depth;
   1.432 +      }
   1.433 +      return false;
   1.434 +    }
   1.435 +
   1.436 +    //calculate the lookup table for cutting the search tree
   1.437 +    void initRNew1tRInOut1t(){
   1.438 +      typename G1::template NodeMap<int> tmp(_g1,0);
   1.439 +      for(unsigned int i=0; i<_order.size(); ++i) {
   1.440 +        tmp[_order[i]]=-1;
   1.441 +        for(typename G1::IncEdgeIt e1(_g1,_order[i]); e1!=INVALID; ++e1) {
   1.442 +          const typename G1::Node currNode=_g1.oppositeNode(_order[i],e1);
   1.443 +          if(tmp[currNode]>0)
   1.444 +            ++_labelTmp1[_intLabels1[currNode]];
   1.445 +          else if(tmp[currNode]==0)
   1.446 +            ++_labelTmp2[_intLabels1[currNode]];
   1.447 +        }
   1.448 +        //_labelTmp1[i]=number of neightbours with label i in set rInOut
   1.449 +        //_labelTmp2[i]=number of neightbours with label i in set rNew
   1.450 +        for(typename G1::IncEdgeIt e1(_g1,_order[i]); e1!=INVALID; ++e1) {
   1.451 +          const int& currIntLabel = _intLabels1[_g1.oppositeNode(_order[i],e1)];
   1.452 +          if(_labelTmp1[currIntLabel]>0) {
   1.453 +            _rInOutLabels1[_order[i]]
   1.454 +              .push_back(std::make_pair(currIntLabel,
   1.455 +                                        _labelTmp1[currIntLabel]));
   1.456 +            _labelTmp1[currIntLabel]=0;
   1.457 +          }
   1.458 +          else if(_labelTmp2[currIntLabel]>0) {
   1.459 +            _rNewLabels1[_order[i]].
   1.460 +              push_back(std::make_pair(currIntLabel,_labelTmp2[currIntLabel]));
   1.461 +            _labelTmp2[currIntLabel]=0;
   1.462 +          }
   1.463 +        }
   1.464 +
   1.465 +        for(typename G1::IncEdgeIt e1(_g1,_order[i]); e1!=INVALID; ++e1) {
   1.466 +          int& tmpCurrNode=tmp[_g1.oppositeNode(_order[i],e1)];
   1.467 +          if(tmpCurrNode!=-1)
   1.468 +            ++tmpCurrNode;
   1.469 +        }
   1.470 +      }
   1.471 +    }
   1.472 +
   1.473 +    int getMaxLabel() const{
   1.474 +      int m=-1;
   1.475 +      for(typename G1::NodeIt n1(_g1); n1!=INVALID; ++n1) {
   1.476 +        const int& currIntLabel = _intLabels1[n1];
   1.477 +        if(currIntLabel>m)
   1.478 +          m=currIntLabel;
   1.479 +      }
   1.480 +      for(typename G2::NodeIt n2(_g2); n2!=INVALID; ++n2) {
   1.481 +        const int& currIntLabel = _intLabels2[n2];
   1.482 +        if(currIntLabel>m)
   1.483 +          m=currIntLabel;
   1.484 +      }
   1.485 +      return m;
   1.486 +    }
   1.487 +
   1.488 +  public:
   1.489 +    ///Constructor
   1.490 +
   1.491 +    ///Constructor.
   1.492 +    ///\param g1 The graph to be embedded.
   1.493 +    ///\param g2 The graph \e g1 will be embedded into.
   1.494 +    ///\param m The type of the NodeMap storing the mapping.
   1.495 +    ///By default, it is G1::NodeMap<G2::Node>
   1.496 +    ///\param intLabel1 The NodeMap storing the integer node labels of G1.
   1.497 +    ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
   1.498 +    ///different labels.
   1.499 +    ///\param intLabel1 The NodeMap storing the integer node labels of G2.
   1.500 +    ///The labels must be the numbers {0,1,2,..,K-1}, where K is the number of
   1.501 +    ///different labels.
   1.502 +    Vf2pp(const G1 &g1, const G2 &g2,M &m, M1 &intLabels1, M2 &intLabels2) :
   1.503 +      _g1(g1), _g2(g2), _depth(0), _mapping(m), _order(countNodes(g1),INVALID),
   1.504 +      _conn(g2,0), _currEdgeIts(countNodes(g1),INVALID), _rNewLabels1(_g1),
   1.505 +      _rInOutLabels1(_g1), _intLabels1(intLabels1) ,_intLabels2(intLabels2),
   1.506 +      _maxLabel(getMaxLabel()), _labelTmp1(_maxLabel+1),_labelTmp2(_maxLabel+1),
   1.507 +      _mapping_type(SUBGRAPH), _deallocMappingAfterUse(0),
   1.508 +      _deallocLabelsAfterUse(0)
   1.509 +    {
   1.510 +      initOrder();
   1.511 +      initRNew1tRInOut1t();
   1.512 +
   1.513 +      //reset mapping
   1.514 +      for(typename G1::NodeIt n(g1);n!=INVALID;++n)
   1.515 +        m[n]=INVALID;
   1.516 +    }
   1.517 +
   1.518 +    ///Destructor
   1.519 +
   1.520 +    ///Destructor.
   1.521 +    ///
   1.522 +    ~Vf2pp()
   1.523 +    {
   1.524 +      if(_deallocMappingAfterUse)
   1.525 +        delete &_mapping;
   1.526 +      if(_deallocLabelsAfterUse) {
   1.527 +        delete &_intLabels1;
   1.528 +        delete &_intLabels2;
   1.529 +      }
   1.530 +    }
   1.531 +
   1.532 +    ///Returns the current mapping type.
   1.533 +
   1.534 +    ///Returns the current mapping type.
   1.535 +    ///
   1.536 +    MappingType mappingType() const
   1.537 +    {
   1.538 +      return _mapping_type;
   1.539 +    }
   1.540 +
   1.541 +    ///Sets the mapping type
   1.542 +
   1.543 +    ///Sets the mapping type.
   1.544 +    ///
   1.545 +    ///The mapping type is set to \ref SUBGRAPH by default.
   1.546 +    ///
   1.547 +    ///\sa See \ref MappingType for the possible values.
   1.548 +    void mappingType(MappingType m_type)
   1.549 +    {
   1.550 +      _mapping_type = m_type;
   1.551 +    }
   1.552 +
   1.553 +    ///Finds a mapping.
   1.554 +
   1.555 +    ///This method finds a mapping from g1 into g2 according to the mapping
   1.556 +    ///type set by \ref mappingType(MappingType) "mappingType()".
   1.557 +    ///
   1.558 +    ///By subsequent calls, it returns all possible mappings one-by-one.
   1.559 +    ///
   1.560 +    ///\retval true if a mapping is found.
   1.561 +    ///\retval false if there is no (more) mapping.
   1.562 +    bool find()
   1.563 +    {
   1.564 +      switch(_mapping_type)
   1.565 +        {
   1.566 +        case SUBGRAPH:
   1.567 +          return extMatch<SUBGRAPH>();
   1.568 +        case INDUCED:
   1.569 +          return extMatch<INDUCED>();
   1.570 +        case ISOMORPH:
   1.571 +          return extMatch<ISOMORPH>();
   1.572 +        default:
   1.573 +          return false;
   1.574 +        }
   1.575 +    }
   1.576 +  };
   1.577 +
   1.578 +  template<typename G1, typename G2>
   1.579 +  class Vf2ppWizardBase {
   1.580 +  protected:
   1.581 +    typedef G1 Graph1;
   1.582 +    typedef G2 Graph2;
   1.583 +
   1.584 +    const G1 &_g1;
   1.585 +    const G2 &_g2;
   1.586 +
   1.587 +    MappingType _mapping_type;
   1.588 +
   1.589 +    typedef typename G1::template NodeMap<typename G2::Node> Mapping;
   1.590 +    bool _local_mapping;
   1.591 +    void *_mapping;
   1.592 +    void createMapping() {
   1.593 +      _mapping = new Mapping(_g1);
   1.594 +    }
   1.595 +
   1.596 +    bool _local_nodeLabels;
   1.597 +    typedef typename G1::template NodeMap<int> NodeLabels1;
   1.598 +    typedef typename G2::template NodeMap<int> NodeLabels2;
   1.599 +    void *_nodeLabels1, *_nodeLabels2;
   1.600 +    void createNodeLabels() {
   1.601 +      _nodeLabels1 = new NodeLabels1(_g1,0);
   1.602 +      _nodeLabels2 = new NodeLabels2(_g2,0);
   1.603 +    }
   1.604 +
   1.605 +    Vf2ppWizardBase(const G1 &g1,const G2 &g2)
   1.606 +      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH),
   1.607 +        _local_mapping(1), _local_nodeLabels(1) { }
   1.608 +  };
   1.609 +
   1.610 +
   1.611 +  /// \brief Auxiliary class for the function-type interface of %VF2
   1.612 +  /// Plus Plus algorithm.
   1.613 +  ///
   1.614 +  /// This auxiliary class implements the named parameters of
   1.615 +  /// \ref vf2pp() "function-type interface" of \ref Vf2pp algorithm.
   1.616 +  ///
   1.617 +  /// \warning This class is not to be used directly.
   1.618 +  ///
   1.619 +  /// \tparam TR The traits class that defines various types used by the
   1.620 +  /// algorithm.
   1.621 +  template<typename TR>
   1.622 +  class Vf2ppWizard : public TR {
   1.623 +    typedef TR Base;
   1.624 +    typedef typename TR::Graph1 Graph1;
   1.625 +    typedef typename TR::Graph2 Graph2;
   1.626 +    typedef typename TR::Mapping Mapping;
   1.627 +    typedef typename TR::NodeLabels1 NodeLabels1;
   1.628 +    typedef typename TR::NodeLabels2 NodeLabels2;
   1.629 +
   1.630 +    using TR::_g1;
   1.631 +    using TR::_g2;
   1.632 +    using TR::_mapping_type;
   1.633 +    using TR::_mapping;
   1.634 +    using TR::_nodeLabels1;
   1.635 +    using TR::_nodeLabels2;
   1.636 +
   1.637 +  public:
   1.638 +    ///Constructor
   1.639 +    Vf2ppWizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {}
   1.640 +
   1.641 +    ///Copy constructor
   1.642 +    Vf2ppWizard(const Base &b) : Base(b) {}
   1.643 +
   1.644 +
   1.645 +    template<typename T>
   1.646 +    struct SetMappingBase : public Base {
   1.647 +      typedef T Mapping;
   1.648 +      SetMappingBase(const Base &b) : Base(b) {}
   1.649 +    };
   1.650 +
   1.651 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.652 +    ///the mapping.
   1.653 +    ///
   1.654 +    ///\ref named-templ-param "Named parameter" function for setting
   1.655 +    ///the map that stores the found embedding.
   1.656 +    template<typename T>
   1.657 +    Vf2ppWizard< SetMappingBase<T> > mapping(const T &t) {
   1.658 +      Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
   1.659 +      Base::_local_mapping = 0;
   1.660 +      return Vf2ppWizard<SetMappingBase<T> >(*this);
   1.661 +    }
   1.662 +
   1.663 +    template<typename NL1, typename NL2>
   1.664 +    struct SetNodeLabelsBase : public Base {
   1.665 +      typedef NL1 NodeLabels1;
   1.666 +      typedef NL2 NodeLabels2;
   1.667 +      SetNodeLabelsBase(const Base &b) : Base(b) { }
   1.668 +    };
   1.669 +
   1.670 +    ///\brief \ref named-templ-param "Named parameter" for setting the
   1.671 +    ///node labels.
   1.672 +    ///
   1.673 +    ///\ref named-templ-param "Named parameter" function for setting
   1.674 +    ///the node labels.
   1.675 +    ///
   1.676 +    ///\param nodeLabels1 A \ref concepts::ReadMap "readable node map"
   1.677 +    ///of g1 with integer values. In case of K different labels, the labels
   1.678 +    ///must be the numbers {0,1,..,K-1}.
   1.679 +    ///\param nodeLabels2 A \ref concepts::ReadMap "readable node map"
   1.680 +    ///of g2 with integer values. In case of K different labels, the labels
   1.681 +    ///must be the numbers {0,1,..,K-1}.
   1.682 +    template<typename NL1, typename NL2>
   1.683 +    Vf2ppWizard< SetNodeLabelsBase<NL1,NL2> >
   1.684 +    nodeLabels(const NL1 &nodeLabels1, const NL2 &nodeLabels2) {
   1.685 +      Base::_local_nodeLabels = 0;
   1.686 +      Base::_nodeLabels1=
   1.687 +        reinterpret_cast<void*>(const_cast<NL1*>(&nodeLabels1));
   1.688 +      Base::_nodeLabels2=
   1.689 +        reinterpret_cast<void*>(const_cast<NL2*>(&nodeLabels2));
   1.690 +      return Vf2ppWizard<SetNodeLabelsBase<NL1,NL2> >
   1.691 +        (SetNodeLabelsBase<NL1,NL2>(*this));
   1.692 +    }
   1.693 +
   1.694 +
   1.695 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.696 +    ///the mapping type.
   1.697 +    ///
   1.698 +    ///\ref named-templ-param "Named parameter" for setting
   1.699 +    ///the mapping type.
   1.700 +    ///
   1.701 +    ///The mapping type is set to \ref SUBGRAPH by default.
   1.702 +    ///
   1.703 +    ///\sa See \ref MappingType for the possible values.
   1.704 +    Vf2ppWizard<Base> &mappingType(MappingType m_type) {
   1.705 +      _mapping_type = m_type;
   1.706 +      return *this;
   1.707 +    }
   1.708 +
   1.709 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.710 +    ///the mapping type to \ref INDUCED.
   1.711 +    ///
   1.712 +    ///\ref named-templ-param "Named parameter" for setting
   1.713 +    ///the mapping type to \ref INDUCED.
   1.714 +    Vf2ppWizard<Base> &induced() {
   1.715 +      _mapping_type = INDUCED;
   1.716 +      return *this;
   1.717 +    }
   1.718 +
   1.719 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.720 +    ///the mapping type to \ref ISOMORPH.
   1.721 +    ///
   1.722 +    ///\ref named-templ-param "Named parameter" for setting
   1.723 +    ///the mapping type to \ref ISOMORPH.
   1.724 +    Vf2ppWizard<Base> &iso() {
   1.725 +      _mapping_type = ISOMORPH;
   1.726 +      return *this;
   1.727 +    }
   1.728 +
   1.729 +    ///Runs the %VF2 Plus Plus algorithm.
   1.730 +
   1.731 +    ///This method runs the VF2 Plus Plus algorithm.
   1.732 +    ///
   1.733 +    ///\retval true if a mapping is found.
   1.734 +    ///\retval false if there is no mapping.
   1.735 +    bool run() {
   1.736 +      if(Base::_local_mapping)
   1.737 +        Base::createMapping();
   1.738 +      if(Base::_local_nodeLabels)
   1.739 +        Base::createNodeLabels();
   1.740 +
   1.741 +      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >
   1.742 +        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
   1.743 +            *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
   1.744 +            *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
   1.745 +
   1.746 +      alg.mappingType(_mapping_type);
   1.747 +
   1.748 +      const bool ret = alg.find();
   1.749 +
   1.750 +      if(Base::_local_nodeLabels) {
   1.751 +        delete reinterpret_cast<NodeLabels1*>(_nodeLabels1);
   1.752 +        delete reinterpret_cast<NodeLabels2*>(_nodeLabels2);
   1.753 +      }
   1.754 +      if(Base::_local_mapping)
   1.755 +        delete reinterpret_cast<Mapping*>(_mapping);
   1.756 +
   1.757 +      return ret;
   1.758 +    }
   1.759 +
   1.760 +    ///Get a pointer to the generated Vf2pp object.
   1.761 +
   1.762 +    ///Gives a pointer to the generated Vf2pp object.
   1.763 +    ///
   1.764 +    ///\return Pointer to the generated Vf2pp object.
   1.765 +    ///\warning Don't forget to delete the referred Vf2pp object after use.
   1.766 +    Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >*
   1.767 +    getPtrToVf2ppObject(){
   1.768 +      if(Base::_local_mapping)
   1.769 +        Base::createMapping();
   1.770 +      if(Base::_local_nodeLabels)
   1.771 +        Base::createNodeLabels();
   1.772 +
   1.773 +      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2 >* ptr =
   1.774 +        new Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2>
   1.775 +        (_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
   1.776 +         *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
   1.777 +         *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
   1.778 +      ptr->mappingType(_mapping_type);
   1.779 +      if(Base::_local_mapping)
   1.780 +        ptr->_deallocMappingAfterUse=true;
   1.781 +      if(Base::_local_nodeLabels)
   1.782 +        ptr->_deallocLabelMapsAfterUse=true;
   1.783 +
   1.784 +      return ptr;
   1.785 +    }
   1.786 +
   1.787 +    ///Counts the number of mappings.
   1.788 +
   1.789 +    ///This method counts the number of mappings.
   1.790 +    ///
   1.791 +    /// \return The number of mappings.
   1.792 +    int count() {
   1.793 +      if(Base::_local_mapping)
   1.794 +        Base::createMapping();
   1.795 +      if(Base::_local_nodeLabels)
   1.796 +        Base::createNodeLabels();
   1.797 +
   1.798 +      Vf2pp<Graph1, Graph2, Mapping, NodeLabels1, NodeLabels2>
   1.799 +        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping),
   1.800 +            *reinterpret_cast<NodeLabels1*>(_nodeLabels1),
   1.801 +            *reinterpret_cast<NodeLabels2*>(_nodeLabels2));
   1.802 +
   1.803 +      alg.mappingType(_mapping_type);
   1.804 +
   1.805 +      int ret = 0;
   1.806 +      while(alg.find())
   1.807 +        ++ret;
   1.808 +
   1.809 +      if(Base::_local_nodeLabels) {
   1.810 +        delete reinterpret_cast<NodeLabels1*>(_nodeLabels1);
   1.811 +        delete reinterpret_cast<NodeLabels2*>(_nodeLabels2);
   1.812 +      }
   1.813 +      if(Base::_local_mapping)
   1.814 +        delete reinterpret_cast<Mapping*>(_mapping);
   1.815 +
   1.816 +      return ret;
   1.817 +    }
   1.818 +  };
   1.819 +
   1.820 +
   1.821 +  ///Function-type interface for VF2 Plus Plus algorithm.
   1.822 +
   1.823 +  /// \ingroup graph_isomorphism
   1.824 +  ///Function-type interface for VF2 Plus Plus algorithm.
   1.825 +  ///
   1.826 +  ///This function has several \ref named-func-param "named parameters"
   1.827 +  ///declared as the members of class \ref Vf2ppWizard.
   1.828 +  ///The following examples show how to use these parameters.
   1.829 +  ///\code
   1.830 +  ///  ListGraph::NodeMap<ListGraph::Node> m(g);
   1.831 +  ///  // Find an embedding of graph g1 into graph g2
   1.832 +  ///  vf2pp(g1,g2).mapping(m).run();
   1.833 +  ///
   1.834 +  ///  // Check whether graphs g1 and g2 are isomorphic
   1.835 +  ///  bool is_iso = vf2pp(g1,g2).iso().run();
   1.836 +  ///
   1.837 +  ///  // Count the number of isomorphisms
   1.838 +  ///  int num_isos = vf2pp(g1,g2).iso().count();
   1.839 +  ///
   1.840 +  ///  // Iterate through all the induced subgraph mappings
   1.841 +  ///  // of graph g1 into g2 using the labels c1 and c2
   1.842 +  ///  auto* myVf2pp = vf2pp(g1,g2).mapping(m).nodeLabels(c1,c2)
   1.843 +  ///  .induced().getPtrToVf2Object();
   1.844 +  ///  while(myVf2pp->find()){
   1.845 +  ///    //process the current mapping m
   1.846 +  ///  }
   1.847 +  ///  delete myVf22pp;
   1.848 +  ///\endcode
   1.849 +  ///\warning Don't forget to put the \ref Vf2ppWizard::run() "run()",
   1.850 +  ///\ref Vf2ppWizard::count() "count()" or
   1.851 +  ///the \ref Vf2ppWizard::getPtrToVf2ppObject() "getPtrToVf2ppObject()"
   1.852 +  ///to the end of the expression.
   1.853 +  ///\sa Vf2ppWizard
   1.854 +  ///\sa Vf2pp
   1.855 +  template<class G1, class G2>
   1.856 +  Vf2ppWizard<Vf2ppWizardBase<G1,G2> > vf2pp(const G1 &g1, const G2 &g2) {
   1.857 +    return Vf2ppWizard<Vf2ppWizardBase<G1,G2> >(g1,g2);
   1.858 +  }
   1.859 +
   1.860 +}
   1.861 +
   1.862 +#endif
   1.863 +