VF2 algorithm added (#597)
authorPeter Madarasi <madarasip@caesar.elte.hu>
Mon, 30 Mar 2015 17:42:30 +0200
changeset 1141a037254714b3
parent 1140 f8ec64f78b5f
child 1142 2f479109a71d
VF2 algorithm added (#597)

The implementation of this feature was sponsored by QuantumBio Inc.
lemon/vf2.h
test/CMakeLists.txt
test/vf2_test.cc
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/vf2.h	Mon Mar 30 17:42:30 2015 +0200
     1.3 @@ -0,0 +1,523 @@
     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
     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_VF2_H
    1.22 +#define LEMON_VF2_H
    1.23 +
    1.24 +#include <lemon/core.h>
    1.25 +#include <lemon/concepts/graph.h>
    1.26 +#include <lemon/dfs.h>
    1.27 +#include <lemon/bfs.h>
    1.28 +#include <test/test_tools.h>
    1.29 +
    1.30 +#include <vector>
    1.31 +#include <set>
    1.32 +
    1.33 +namespace lemon
    1.34 +{
    1.35 +  namespace bits
    1.36 +  {
    1.37 +    namespace vf2
    1.38 +    {
    1.39 +      class AlwaysEq
    1.40 +      {
    1.41 +      public:
    1.42 +        template<class T1, class T2>
    1.43 +        bool operator()(T1, T2) const
    1.44 +        {
    1.45 +          return true;
    1.46 +        }
    1.47 +      };
    1.48 +
    1.49 +      template<class M1, class M2>
    1.50 +      class MapEq
    1.51 +      {
    1.52 +        const M1 &_m1;
    1.53 +        const M2 &_m2;
    1.54 +      public:
    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 +        {
    1.58 +          return _m1[k1] == _m2[k2];
    1.59 +        }
    1.60 +      };
    1.61 +
    1.62 +      template <class G>
    1.63 +      class DfsLeaveOrder : public DfsVisitor<G>
    1.64 +      {
    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 +          _order[--i]=node;
    1.75 +        }
    1.76 +      };
    1.77 +
    1.78 +      template <class G>
    1.79 +      class BfsLeaveOrder : public BfsVisitor<G>
    1.80 +      {
    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 +          _order[i++]=node;
    1.91 +        }
    1.92 +      };
    1.93 +    }
    1.94 +  }
    1.95 +
    1.96 +  enum MappingType {
    1.97 +    SUBGRAPH = 0,
    1.98 +    INDUCED = 1,
    1.99 +    ISOMORPH = 2
   1.100 +  };
   1.101 +
   1.102 +  template<class G1, class G2, class I, class NEq = bits::vf2::AlwaysEq >
   1.103 +  class Vf2
   1.104 +  {
   1.105 +    //Current depth in the DFS tree.
   1.106 +    int _depth;
   1.107 +    //Functor with bool operator()(G1::Node,G2::Node), which returns 1
   1.108 +    //if and only if the 2 nodes are equivalent.
   1.109 +    NEq _nEq;
   1.110 +
   1.111 +    typename G2::template NodeMap<int> _conn;
   1.112 +    //Current matching. We index it by the nodes of g1, and match[v] is
   1.113 +    //a node of g2.
   1.114 +    I &_match;
   1.115 +    //order[i] is the node of g1, for which we find a pair if depth=i
   1.116 +    std::vector<typename G1::Node> order;
   1.117 +    //currEdgeIts[i] is an edge iterator, witch is last used in the ith
   1.118 +    //depth to find a pair for order[i].
   1.119 +    std::vector<typename G2::IncEdgeIt> currEdgeIts;
   1.120 +    //The small graph.
   1.121 +    const G1 &_g1;
   1.122 +    //The big graph.
   1.123 +    const G2 &_g2;
   1.124 +    //lookup tables for cut the searchtree
   1.125 +    typename G1::template NodeMap<int> rNew1t,rInOut1t;
   1.126 +
   1.127 +    MappingType _mapping_type;
   1.128 +
   1.129 +    //cut the search tree
   1.130 +    template<MappingType MT>
   1.131 +    bool cut(const typename G1::Node n1,const typename G2::Node n2) const
   1.132 +    {
   1.133 +      int rNew2=0,rInOut2=0;
   1.134 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.135 +        {
   1.136 +          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.137 +          if(_conn[currNode]>0)
   1.138 +            ++rInOut2;
   1.139 +          else if(MT!=SUBGRAPH&&_conn[currNode]==0)
   1.140 +            ++rNew2;
   1.141 +        }
   1.142 +      switch(MT)
   1.143 +        {
   1.144 +        case INDUCED:
   1.145 +          return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
   1.146 +        case SUBGRAPH:
   1.147 +          return rInOut1t[n1]<=rInOut2;
   1.148 +        case ISOMORPH:
   1.149 +          return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
   1.150 +        default:
   1.151 +          return false;
   1.152 +        }
   1.153 +    }
   1.154 +
   1.155 +    template<MappingType MT>
   1.156 +    bool feas(const typename G1::Node n1,const typename G2::Node n2)
   1.157 +    {
   1.158 +      if(!_nEq(n1,n2))
   1.159 +        return 0;
   1.160 +
   1.161 +      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
   1.162 +        {
   1.163 +          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
   1.164 +          if(_match[currNode]!=INVALID)
   1.165 +            --_conn[_match[currNode]];
   1.166 +        }
   1.167 +      bool isIso=1;
   1.168 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.169 +        {
   1.170 +          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.171 +          if(_conn[currNode]<-1)
   1.172 +            ++_conn[currNode];
   1.173 +          else if(MT!=SUBGRAPH&&_conn[currNode]==-1)
   1.174 +            {
   1.175 +              isIso=0;
   1.176 +              break;
   1.177 +            }
   1.178 +        }
   1.179 +
   1.180 +      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
   1.181 +        {
   1.182 +          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
   1.183 +          if(_match[currNode]!=INVALID&&_conn[_match[currNode]]!=-1)
   1.184 +            {
   1.185 +              switch(MT)
   1.186 +                {
   1.187 +                case INDUCED:
   1.188 +                case ISOMORPH:
   1.189 +                  isIso=0;
   1.190 +                  break;
   1.191 +                case SUBGRAPH:
   1.192 +                  if(_conn[_match[currNode]]<-1)
   1.193 +                    isIso=0;
   1.194 +                  break;
   1.195 +                }
   1.196 +              _conn[_match[currNode]]=-1;
   1.197 +            }
   1.198 +        }
   1.199 +      return isIso&&cut<MT>(n1,n2);
   1.200 +    }
   1.201 +
   1.202 +    void addPair(const typename G1::Node n1,const typename G2::Node n2)
   1.203 +    {
   1.204 +      _conn[n2]=-1;
   1.205 +      _match.set(n1,n2);
   1.206 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.207 +        if(_conn[_g2.oppositeNode(n2,e2)]!=-1)
   1.208 +          ++_conn[_g2.oppositeNode(n2,e2)];
   1.209 +    }
   1.210 +
   1.211 +    void subPair(const typename G1::Node n1,const typename G2::Node n2)
   1.212 +    {
   1.213 +      _conn[n2]=0;
   1.214 +      _match.set(n1,INVALID);
   1.215 +      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
   1.216 +        {
   1.217 +          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
   1.218 +          if(_conn[currNode]>0)
   1.219 +            --_conn[currNode];
   1.220 +          else if(_conn[currNode]==-1)
   1.221 +            ++_conn[n2];
   1.222 +        }
   1.223 +    }
   1.224 +
   1.225 +    void setOrder()//we will find pairs for the nodes of g1 in this order
   1.226 +    {
   1.227 +      // bits::vf2::DfsLeaveOrder<G1> v(_g1,order);
   1.228 +      //   DfsVisit<G1,bits::vf2::DfsLeaveOrder<G1> >dfs(_g1, v);
   1.229 +      //   dfs.run();
   1.230 +
   1.231 +      //it is more efficient in practice than DFS
   1.232 +      bits::vf2::BfsLeaveOrder<G1> v(_g1,order);
   1.233 +      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> >bfs(_g1, v);
   1.234 +      bfs.run();
   1.235 +    }
   1.236 +
   1.237 +  public:
   1.238 +
   1.239 +    template<MappingType MT>
   1.240 +    bool extMatch()
   1.241 +    {
   1.242 +      while(_depth>=0)
   1.243 +        {
   1.244 +          //there are not nodes in g1, which has not pair in g2.
   1.245 +          if(_depth==static_cast<int>(order.size()))
   1.246 +            {
   1.247 +              --_depth;
   1.248 +              return true;
   1.249 +            }
   1.250 +          //the node of g2, which neighbours are the candidates for
   1.251 +          //the pair of order[_depth]
   1.252 +          typename G2::Node currPNode;
   1.253 +          if(currEdgeIts[_depth]==INVALID)
   1.254 +            {
   1.255 +              typename G1::IncEdgeIt fstMatchedE(_g1,order[_depth]);
   1.256 +              //if _match[order[_depth]]!=INVALID, we dont use
   1.257 +              //fstMatchedE
   1.258 +              if(_match[order[_depth]]==INVALID)
   1.259 +                for(; fstMatchedE!=INVALID &&
   1.260 +                      _match[_g1.oppositeNode(order[_depth],
   1.261 +                                              fstMatchedE)]==INVALID;
   1.262 +                    ++fstMatchedE) ; //find fstMatchedE
   1.263 +              if(fstMatchedE==INVALID||_match[order[_depth]]!=INVALID)
   1.264 +                {
   1.265 +                  //We did not find an covered neighbour, this means
   1.266 +                  //the graph is not connected(or _depth==0).  Every
   1.267 +                  //uncovered(and there are some other properties due
   1.268 +                  //to the spec. problem types) node of g2 is
   1.269 +                  //candidate.  We can read the iterator of the last
   1.270 +                  //tryed node from the match if it is not the first
   1.271 +                  //try(match[order[_depth]]!=INVALID)
   1.272 +                  typename G2::NodeIt n2(_g2);
   1.273 +                  //if its not the first try
   1.274 +                  if(_match[order[_depth]]!=INVALID)
   1.275 +                    {
   1.276 +                      n2=++typename G2::NodeIt(_g2,_match[order[_depth]]);
   1.277 +                      subPair(order[_depth],_match[order[_depth]]);
   1.278 +                    }
   1.279 +                  for(; n2!=INVALID; ++n2)
   1.280 +                    if(MT!=SUBGRAPH&&_conn[n2]==0)
   1.281 +                      {
   1.282 +                        if(feas<MT>(order[_depth],n2))
   1.283 +                          break;
   1.284 +                      }
   1.285 +                    else if(MT==SUBGRAPH&&_conn[n2]>=0)
   1.286 +                      if(feas<MT>(order[_depth],n2))
   1.287 +                        break;
   1.288 +                  // n2 is the next candidate
   1.289 +                  if(n2!=INVALID)
   1.290 +                    {
   1.291 +                      addPair(order[_depth],n2);
   1.292 +                      ++_depth;
   1.293 +                    }
   1.294 +                  else // there is no more candidate
   1.295 +                    --_depth;
   1.296 +                  continue;
   1.297 +                }
   1.298 +              else
   1.299 +                {
   1.300 +                  currPNode=_match[_g1.oppositeNode(order[_depth],fstMatchedE)];
   1.301 +                  currEdgeIts[_depth]=typename G2::IncEdgeIt(_g2,currPNode);
   1.302 +                }
   1.303 +            }
   1.304 +          else
   1.305 +            {
   1.306 +              currPNode=_g2.oppositeNode(_match[order[_depth]],
   1.307 +                                         currEdgeIts[_depth]);
   1.308 +              subPair(order[_depth],_match[order[_depth]]);
   1.309 +              ++currEdgeIts[_depth];
   1.310 +            }
   1.311 +          for(; currEdgeIts[_depth]!=INVALID; ++currEdgeIts[_depth])
   1.312 +            {
   1.313 +              const typename G2::Node currNode =
   1.314 +                _g2.oppositeNode(currPNode, currEdgeIts[_depth]);
   1.315 +              //if currNode is uncovered
   1.316 +              if(_conn[currNode]>0&&feas<MT>(order[_depth],currNode))
   1.317 +                {
   1.318 +                  addPair(order[_depth],currNode);
   1.319 +                  break;
   1.320 +                }
   1.321 +            }
   1.322 +          currEdgeIts[_depth]==INVALID?--_depth:++_depth;
   1.323 +        }
   1.324 +      return false;
   1.325 +    }
   1.326 +
   1.327 +    //calc. the lookup table for cut the searchtree
   1.328 +    void setRNew1tRInOut1t()
   1.329 +    {
   1.330 +      typename G1::template NodeMap<int> tmp(_g1,0);
   1.331 +      for(unsigned int i=0; i<order.size(); ++i)
   1.332 +        {
   1.333 +          tmp[order[i]]=-1;
   1.334 +          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
   1.335 +            {
   1.336 +              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
   1.337 +              if(tmp[currNode]>0)
   1.338 +                ++rInOut1t[order[i]];
   1.339 +              else if(tmp[currNode]==0)
   1.340 +                ++rNew1t[order[i]];
   1.341 +            }
   1.342 +          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
   1.343 +            {
   1.344 +              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
   1.345 +              if(tmp[currNode]!=-1)
   1.346 +                ++tmp[currNode];
   1.347 +            }
   1.348 +        }
   1.349 +    }
   1.350 +  public:
   1.351 +    Vf2(const G1 &g1, const G2 &g2,I &i, const NEq &nEq = NEq() ) :
   1.352 +      _nEq(nEq), _conn(g2,0), _match(i), order(countNodes(g1)),
   1.353 +      currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNew1t(g1,0),
   1.354 +      rInOut1t(g1,0), _mapping_type(SUBGRAPH)
   1.355 +    {
   1.356 +      _depth=0;
   1.357 +      setOrder();
   1.358 +      setRNew1tRInOut1t();
   1.359 +    }
   1.360 +
   1.361 +    MappingType mappingType() const { return _mapping_type; }
   1.362 +    void mappingType(MappingType m_type) { _mapping_type = m_type; }
   1.363 +
   1.364 +    bool find()
   1.365 +    {
   1.366 +      switch(_mapping_type)
   1.367 +        {
   1.368 +        case SUBGRAPH:
   1.369 +          return extMatch<SUBGRAPH>();
   1.370 +        case INDUCED:
   1.371 +          return extMatch<INDUCED>();
   1.372 +        case ISOMORPH:
   1.373 +          return extMatch<ISOMORPH>();
   1.374 +        default:
   1.375 +          return false;
   1.376 +        }
   1.377 +    }
   1.378 +  };
   1.379 +
   1.380 +  template<class G1, class G2>
   1.381 +  class Vf2WizardBase
   1.382 +  {
   1.383 +  protected:
   1.384 +    typedef G1 Graph1;
   1.385 +    typedef G2 Graph2;
   1.386 +
   1.387 +    const G1 &_g1;
   1.388 +    const G2 &_g2;
   1.389 +
   1.390 +    MappingType _mapping_type;
   1.391 +
   1.392 +    typedef typename G1::template NodeMap<typename G2::Node> Mapping;
   1.393 +    bool _local_mapping;
   1.394 +    void *_mapping;
   1.395 +    void createMapping()
   1.396 +    {
   1.397 +      _mapping = new Mapping(_g1);
   1.398 +    }
   1.399 +
   1.400 +    typedef bits::vf2::AlwaysEq NodeEq;
   1.401 +    NodeEq _node_eq;
   1.402 +
   1.403 +    Vf2WizardBase(const G1 &g1,const G2 &g2)
   1.404 +      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) {}
   1.405 +  };
   1.406 +
   1.407 +  template<class TR>
   1.408 +  class Vf2Wizard : public TR
   1.409 +  {
   1.410 +    typedef TR Base;
   1.411 +    typedef typename TR::Graph1 Graph1;
   1.412 +    typedef typename TR::Graph2 Graph2;
   1.413 +
   1.414 +    typedef typename TR::Mapping Mapping;
   1.415 +    typedef typename TR::NodeEq NodeEq;
   1.416 +
   1.417 +    using TR::_g1;
   1.418 +    using TR::_g2;
   1.419 +    using TR::_mapping_type;
   1.420 +    using TR::_mapping;
   1.421 +    using TR::_node_eq;
   1.422 +
   1.423 +  public:
   1.424 +    ///Copy constructor
   1.425 +    Vf2Wizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {}
   1.426 +
   1.427 +    ///Copy constructor
   1.428 +    Vf2Wizard(const Base &b) : Base(b) {}
   1.429 +
   1.430 +
   1.431 +    template<class T>
   1.432 +    struct SetMappingBase : public Base {
   1.433 +      typedef T Mapping;
   1.434 +      SetMappingBase(const Base &b) : Base(b) {}
   1.435 +    };
   1.436 +
   1.437 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.438 +    ///the mapping.
   1.439 +    ///
   1.440 +    ///\ref named-templ-param "Named parameter" function for setting
   1.441 +    ///the map that stores the found embedding.
   1.442 +    template<class T>
   1.443 +    Vf2Wizard< SetMappingBase<T> > mapping(const T &t)
   1.444 +    {
   1.445 +      Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
   1.446 +      Base::_local_mapping = false;
   1.447 +      return Vf2Wizard<SetMappingBase<T> >(*this);
   1.448 +    }
   1.449 +
   1.450 +    template<class NE>
   1.451 +    struct SetNodeEqBase : public Base {
   1.452 +      typedef NE NodeEq;
   1.453 +      NodeEq _node_eq;
   1.454 +      SetNodeEqBase(const Base &b, const NE &node_eq)
   1.455 +        : Base(b), _node_eq(node_eq) {}
   1.456 +    };
   1.457 +
   1.458 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.459 +    ///the node equivalence relation.
   1.460 +    ///
   1.461 +    ///\ref named-templ-param "Named parameter" function for setting
   1.462 +    ///the equivalence relation between the nodes.
   1.463 +    template<class T>
   1.464 +    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq)
   1.465 +    {
   1.466 +      return Vf2Wizard<SetNodeEqBase<T> >(SetNodeEqBase<T>(*this,node_eq));
   1.467 +    }
   1.468 +
   1.469 +    ///\brief \ref named-templ-param "Named parameter" for setting
   1.470 +    ///the node labels.
   1.471 +    ///
   1.472 +    ///\ref named-templ-param "Named parameter" function for setting
   1.473 +    ///the node labels defining equivalence relation between them.
   1.474 +    template<class M1, class M2>
   1.475 +    Vf2Wizard< SetNodeEqBase<bits::vf2::MapEq<M1,M2> > >
   1.476 +    nodeLabels(const M1 &m1,const M2 &m2)
   1.477 +    {
   1.478 +      return nodeEq(bits::vf2::MapEq<M1,M2>(m1,m2));
   1.479 +    }
   1.480 +
   1.481 +    Vf2Wizard<Base> &mappingType(MappingType m_type)
   1.482 +    {
   1.483 +      _mapping_type = m_type;
   1.484 +      return *this;
   1.485 +    }
   1.486 +
   1.487 +    Vf2Wizard<Base> &induced()
   1.488 +    {
   1.489 +      _mapping_type = INDUCED;
   1.490 +      return *this;
   1.491 +    }
   1.492 +
   1.493 +    Vf2Wizard<Base> &iso()
   1.494 +    {
   1.495 +      _mapping_type = ISOMORPH;
   1.496 +      return *this;
   1.497 +    }
   1.498 +
   1.499 +    bool run()
   1.500 +    {
   1.501 +      if(Base::_local_mapping)
   1.502 +        Base::createMapping();
   1.503 +
   1.504 +      Vf2<Graph1, Graph2, Mapping, NodeEq >
   1.505 +        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
   1.506 +
   1.507 +      alg.mappingType(_mapping_type);
   1.508 +
   1.509 +      bool ret = alg.find();
   1.510 +
   1.511 +      if(Base::_local_mapping)
   1.512 +        delete reinterpret_cast<Mapping*>(_mapping);
   1.513 +
   1.514 +      return ret;
   1.515 +    }
   1.516 +  };
   1.517 +
   1.518 +  template<class G1, class G2>
   1.519 +  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2)
   1.520 +  {
   1.521 +    return Vf2Wizard<Vf2WizardBase<G1,G2> >(g1,g2);
   1.522 +  }
   1.523 +
   1.524 +}
   1.525 +
   1.526 +#endif
     2.1 --- a/test/CMakeLists.txt	Thu May 07 11:42:19 2015 +0200
     2.2 +++ b/test/CMakeLists.txt	Mon Mar 30 17:42:30 2015 +0200
     2.3 @@ -54,6 +54,7 @@
     2.4    time_measure_test
     2.5    tsp_test
     2.6    unionfind_test
     2.7 +  vf2_test
     2.8  )
     2.9  
    2.10  IF(LEMON_HAVE_LP)
     3.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     3.2 +++ b/test/vf2_test.cc	Mon Mar 30 17:42:30 2015 +0200
     3.3 @@ -0,0 +1,361 @@
     3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     3.5 + *
     3.6 + * This file is a part of LEMON, a generic C++ optimization library.
     3.7 + *
     3.8 + * Copyright (C) 2015
     3.9 + * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
    3.10 + *
    3.11 + * Permission to use, modify and distribute this software is granted
    3.12 + * provided that this copyright notice appears in all copies. For
    3.13 + * precise terms see the accompanying LICENSE file.
    3.14 + *
    3.15 + * This software is provided "AS IS" with no warranty of any kind,
    3.16 + * express or implied, and with no claim as to its suitability for any
    3.17 + * purpose.
    3.18 + *
    3.19 + */
    3.20 +
    3.21 +#include <lemon/vf2.h>
    3.22 +#include <lemon/concepts/digraph.h>
    3.23 +#include <lemon/smart_graph.h>
    3.24 +#include <lemon/lgf_reader.h>
    3.25 +#include <lemon/concepts/maps.h>
    3.26 +
    3.27 +#include <test/test_tools.h>
    3.28 +#include <sstream>
    3.29 +
    3.30 +using namespace lemon;
    3.31 +
    3.32 +char petersen_lgf[] =
    3.33 +  "@nodes\n"
    3.34 +  "label col1 col2\n"
    3.35 +  "0 1 1\n"
    3.36 +  "1 1 2\n"
    3.37 +  "2 1 3\n"
    3.38 +  "3 1 4\n"
    3.39 +  "4 2 5\n"
    3.40 +  "5 2 1\n"
    3.41 +  "6 2 2\n"
    3.42 +  "7 2 3\n"
    3.43 +  "8 2 4\n"
    3.44 +  "9 2 5\n"
    3.45 +  "@arcs\n"
    3.46 +  "     -\n"
    3.47 +  "0 1\n"
    3.48 +  "1 2\n"
    3.49 +  "2 3\n"
    3.50 +  "3 4\n"
    3.51 +  "4 0\n"
    3.52 +  "0 5\n"
    3.53 +  "1 6\n"
    3.54 +  "2 7\n"
    3.55 +  "3 8\n"
    3.56 +  "4 9\n"
    3.57 +  "5 8\n"
    3.58 +  "5 7\n"
    3.59 +  "9 6\n"
    3.60 +  "9 7\n"
    3.61 +  "6 8\n";
    3.62 +
    3.63 +char c5_lgf[] =
    3.64 +  "@nodes\n"
    3.65 +  "label col\n"
    3.66 +  "0 1\n"
    3.67 +  "1 2\n"
    3.68 +  "2 3\n"
    3.69 +  "3 4\n"
    3.70 +  "4 5\n"
    3.71 +  "@arcs\n"
    3.72 +  "     -\n"
    3.73 +  "0 1\n"
    3.74 +  "1 2\n"
    3.75 +  "2 3\n"
    3.76 +  "3 4\n"
    3.77 +  "4 0\n";
    3.78 +
    3.79 +char c7_lgf[] =
    3.80 +  "@nodes\n"
    3.81 +  "label\n"
    3.82 +  "0\n"
    3.83 +  "1\n"
    3.84 +  "2\n"
    3.85 +  "3\n"
    3.86 +  "4\n"
    3.87 +  "5\n"
    3.88 +  "6\n"
    3.89 +  "@arcs\n"
    3.90 +  "     -\n"
    3.91 +  "0 1\n"
    3.92 +  "1 2\n"
    3.93 +  "2 3\n"
    3.94 +  "3 4\n"
    3.95 +  "4 5\n"
    3.96 +  "5 6\n"
    3.97 +  "6 0\n";
    3.98 +
    3.99 +char c10_lgf[] =
   3.100 +  "@nodes\n"
   3.101 +  "label\n"
   3.102 +  "0\n"
   3.103 +  "1\n"
   3.104 +  "2\n"
   3.105 +  "3\n"
   3.106 +  "4\n"
   3.107 +  "5\n"
   3.108 +  "6\n"
   3.109 +  "7\n"
   3.110 +  "8\n"
   3.111 +  "9\n"
   3.112 +  "@arcs\n"
   3.113 +  "     -\n"
   3.114 +  "0 1\n"
   3.115 +  "1 2\n"
   3.116 +  "2 3\n"
   3.117 +  "3 4\n"
   3.118 +  "4 5\n"
   3.119 +  "5 6\n"
   3.120 +  "6 7\n"
   3.121 +  "7 8\n"
   3.122 +  "8 9\n"
   3.123 +  "9 0\n";
   3.124 +
   3.125 +char p10_lgf[] =
   3.126 +  "@nodes\n"
   3.127 +  "label\n"
   3.128 +  "0\n"
   3.129 +  "1\n"
   3.130 +  "2\n"
   3.131 +  "3\n"
   3.132 +  "4\n"
   3.133 +  "5\n"
   3.134 +  "6\n"
   3.135 +  "7\n"
   3.136 +  "8\n"
   3.137 +  "9\n"
   3.138 +  "@arcs\n"
   3.139 +  "     -\n"
   3.140 +  "0 1\n"
   3.141 +  "1 2\n"
   3.142 +  "2 3\n"
   3.143 +  "3 4\n"
   3.144 +  "4 5\n"
   3.145 +  "5 6\n"
   3.146 +  "6 7\n"
   3.147 +  "7 8\n"
   3.148 +  "8 9\n";
   3.149 +
   3.150 +SmartGraph petersen, c5, c7, c10, p10;
   3.151 +SmartGraph::NodeMap<int> petersen_col1(petersen);
   3.152 +SmartGraph::NodeMap<int> petersen_col2(petersen);
   3.153 +SmartGraph::NodeMap<int> c5_col(c5);
   3.154 +
   3.155 +void make_graphs() {
   3.156 +  std::stringstream ss(petersen_lgf);
   3.157 +  graphReader(petersen, ss)
   3.158 +    .nodeMap("col1",petersen_col1)
   3.159 +    .nodeMap("col2",petersen_col2)
   3.160 +    .run();
   3.161 +
   3.162 +  ss.clear();
   3.163 +  ss.str("");
   3.164 +  ss<<c5_lgf;
   3.165 +  //std::stringstream ss2(c5_lgf);
   3.166 +  graphReader(c5, ss)
   3.167 +    .nodeMap("col",c5_col)
   3.168 +    .run();
   3.169 +
   3.170 +  ss.clear();
   3.171 +  ss.str("");
   3.172 +  ss<<c7_lgf;
   3.173 +  graphReader(c7, ss).run();
   3.174 +
   3.175 +  ss.clear();
   3.176 +  ss.str("");
   3.177 +  ss<<c10_lgf;
   3.178 +  graphReader(c10, ss).run();
   3.179 +
   3.180 +  ss.clear();
   3.181 +  ss.str("");
   3.182 +  ss<<p10_lgf;
   3.183 +  graphReader(p10, ss).run();
   3.184 +
   3.185 +}
   3.186 +
   3.187 +class EqComparable {
   3.188 +public:
   3.189 +  bool operator==(const EqComparable&) { return false; }
   3.190 +};
   3.191 +
   3.192 +template<class A, class B>
   3.193 +class EqClass {
   3.194 +public:
   3.195 +  bool operator()(A, B) { return false; }
   3.196 +};
   3.197 +
   3.198 +template<class G1,class G2>
   3.199 +void checkVf2Compile()
   3.200 +{
   3.201 +  G1 g;
   3.202 +  G2 h;
   3.203 +  concepts::ReadWriteMap<typename G1::Node, typename G2::Node> r;
   3.204 +  bool succ;
   3.205 +  ::lemon::ignore_unused_variable_warning(succ);
   3.206 +
   3.207 +  succ = vf2(g,h).run();
   3.208 +  succ = vf2(g,h).induced().run();
   3.209 +  succ = vf2(g,h).iso().run();
   3.210 +  succ = vf2(g,h).mapping(r).run();
   3.211 +  succ = vf2(g,h).induced().mapping(r).run();
   3.212 +  succ = vf2(g,h).iso().mapping(r).run();
   3.213 +  concepts::ReadMap<typename G1::Node, EqComparable> l1;
   3.214 +  concepts::ReadMap<typename G2::Node, EqComparable> l2;
   3.215 +  succ = vf2(g,h).nodeLabels(l1,l2).mapping(r).run();
   3.216 +  succ = vf2(g,h).nodeEq(EqClass<typename G1::Node,typename G2::Node>())
   3.217 +    .mapping(r).run();
   3.218 +}
   3.219 +
   3.220 +void justCompile()
   3.221 +{
   3.222 +  checkVf2Compile<concepts::Graph,concepts::Graph>();
   3.223 +  checkVf2Compile<concepts::Graph,SmartGraph>();
   3.224 +  checkVf2Compile<SmartGraph,concepts::Graph>();
   3.225 +}
   3.226 +
   3.227 +template<class G1, class G2, class I>
   3.228 +void checkSub(const G1 &g1, const G2 &g2, const I &i)
   3.229 +{
   3.230 +  {
   3.231 +    std::set<typename G2::Node> image;
   3.232 +    for(typename G1::NodeIt n(g1);n!=INVALID;++n)
   3.233 +      {
   3.234 +          check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
   3.235 +          check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
   3.236 +          image.insert(i[n]);
   3.237 +      }
   3.238 +  }
   3.239 +  for(typename G1::EdgeIt e(g1);e!=INVALID;++e)
   3.240 +    check(findEdge(g2,i[g1.u(e)],i[g1.v(e)])!=INVALID,
   3.241 +          "Wrong isomorphism: missing edge(checkSub).");
   3.242 +}
   3.243 +
   3.244 +template<class G1, class G2, class I>
   3.245 +void checkInd(const G1 &g1, const G2 &g2, const I &i)
   3.246 +  {
   3.247 +  std::set<typename G2::Node> image;
   3.248 +  for(typename G1::NodeIt n(g1);n!=INVALID;++n)
   3.249 +    {
   3.250 +    check(i[n]!=INVALID, "Wrong isomorphism: incomplete mapping.");
   3.251 +    check(image.count(i[n])==0,"Wrong isomorphism: not injective.");
   3.252 +    image.insert(i[n]);
   3.253 +    }
   3.254 +  for(typename G1::NodeIt n(g1); n!=INVALID; ++n)
   3.255 +    for(typename G1::NodeIt m(g1); m!=INVALID; ++m)
   3.256 +      if((findEdge(g1,n,m)==INVALID) != (findEdge(g2,i[n],i[m])==INVALID))
   3.257 +        {
   3.258 +        std::cout << "Wrong isomorphism: edge mismatch";
   3.259 +        exit(1);
   3.260 +        }
   3.261 +  }
   3.262 +
   3.263 +template<class G1,class G2>
   3.264 +int checkSub(const G1 &g1, const G2 &g2)
   3.265 +{
   3.266 +  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
   3.267 +  if(vf2(g1,g2).mapping(iso).run())
   3.268 +    {
   3.269 +      checkSub(g1,g2,iso);
   3.270 +      return true;
   3.271 +    }
   3.272 +  else return false;
   3.273 +}
   3.274 +
   3.275 +template<class G1,class G2>
   3.276 +int checkInd(const G1 &g1, const G2 &g2)
   3.277 +{
   3.278 +  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
   3.279 +  if(vf2(g1,g2).induced().mapping(iso).run())
   3.280 +    {
   3.281 +      checkInd(g1,g2,iso);
   3.282 +      return true;
   3.283 +    }
   3.284 +  else return false;
   3.285 +}
   3.286 +
   3.287 +template<class G1,class G2>
   3.288 +int checkIso(const G1 &g1, const G2 &g2)
   3.289 +{
   3.290 +  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
   3.291 +  if(vf2(g1,g2).iso().mapping(iso).run())
   3.292 +    {
   3.293 +      check(countNodes(g1)==countNodes(g2),
   3.294 +            "Wrong iso alg.: they are not isomophic.");
   3.295 +      checkInd(g1,g2,iso);
   3.296 +      return true;
   3.297 +    }
   3.298 +  else return false;
   3.299 +}
   3.300 +
   3.301 +template<class G1, class G2, class L1, class L2, class I>
   3.302 +void checkLabel(const G1 &g1, const G2 &,
   3.303 +                const L1 &l1, const L2 &l2,const I &i)
   3.304 +{
   3.305 +  for(typename G1::NodeIt n(g1);n!=INVALID;++n)
   3.306 +    {
   3.307 +      check(l1[n]==l2[i[n]],"Wrong isomorphism: label mismatch.");
   3.308 +    }
   3.309 +}
   3.310 +
   3.311 +template<class G1,class G2,class L1,class L2>
   3.312 +int checkSub(const G1 &g1, const G2 &g2, const L1 &l1, const L2 &l2)
   3.313 +{
   3.314 +  typename G1:: template NodeMap<typename G2::Node> iso(g1,INVALID);
   3.315 +  if(vf2(g1,g2).nodeLabels(l1,l2).mapping(iso).run())
   3.316 +    {
   3.317 +      checkSub(g1,g2,iso);
   3.318 +      checkLabel(g1,g2,l1,l2,iso);
   3.319 +      return true;
   3.320 +    }
   3.321 +  else return false;
   3.322 +}
   3.323 +
   3.324 +int main() {
   3.325 +  make_graphs();
   3.326 +  check(checkSub(c5,petersen), "There should exist a C5->Petersen mapping.");
   3.327 +  check(!checkSub(c7,petersen),
   3.328 +        "There should not exist a C7->Petersen mapping.");
   3.329 +  check(checkSub(p10,petersen), "There should exist a P10->Petersen mapping.");
   3.330 +  check(!checkSub(c10,petersen),
   3.331 +        "There should not exist a C10->Petersen mapping.");
   3.332 +  check(checkSub(petersen,petersen),
   3.333 +        "There should exist a Petersen->Petersen mapping.");
   3.334 +
   3.335 +  check(checkInd(c5,petersen),
   3.336 +        "There should exist a C5->Petersen spanned mapping.");
   3.337 +  check(!checkInd(c7,petersen),
   3.338 +        "There should exist a C7->Petersen spanned mapping.");
   3.339 +  check(!checkInd(p10,petersen),
   3.340 +        "There should not exist a P10->Petersen spanned mapping.");
   3.341 +  check(!checkInd(c10,petersen),
   3.342 +        "There should not exist a C10->Petersen spanned mapping.");
   3.343 +  check(checkInd(petersen,petersen),
   3.344 +        "There should exist a Petersen->Petersen spanned mapping.");
   3.345 +
   3.346 +  check(!checkSub(petersen,c10),
   3.347 +        "There should not exist a Petersen->C10 mapping.");
   3.348 +  check(checkSub(p10,c10),
   3.349 +        "There should exist a P10->C10 mapping.");
   3.350 +  check(!checkInd(p10,c10),
   3.351 +        "There should not exist a P10->C10 spanned mapping.");
   3.352 +  check(!checkSub(c10,p10),
   3.353 +        "There should not exist a C10->P10 mapping.");
   3.354 +
   3.355 +  check(!checkIso(p10,c10),
   3.356 +        "P10 and C10 are not isomorphic.");
   3.357 +  check(checkIso(c10,c10),
   3.358 +        "C10 and C10 are not isomorphic.");
   3.359 +
   3.360 +  check(!checkSub(c5,petersen,c5_col,petersen_col1),
   3.361 +        "There should exist a C5->Petersen mapping.");
   3.362 +  check(checkSub(c5,petersen,c5_col,petersen_col2),
   3.363 +        "There should exist a C5->Petersen mapping.");
   3.364 +}