COIN-OR::LEMON - Graph Library

source: lemon/lemon/vf2.h @ 1350:a037254714b3

Last change on this file since 1350:a037254714b3 was 1350:a037254714b3, checked in by Peter Madarasi <madarasip@…>, 4 years ago

VF2 algorithm added (#597)

The implementation of this feature was sponsored by QuantumBio? Inc.

File size: 14.8 KB
Line 
1/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 *
3 * This file is a part of LEMON, a generic C++ optimization library.
4 *
5 * Copyright (C) 2015
6 * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
7 *
8 * Permission to use, modify and distribute this software is granted
9 * provided that this copyright notice appears in all copies. For
10 * precise terms see the accompanying LICENSE file.
11 *
12 * This software is provided "AS IS" with no warranty of any kind,
13 * express or implied, and with no claim as to its suitability for any
14 * purpose.
15 *
16 */
17
18#ifndef LEMON_VF2_H
19#define LEMON_VF2_H
20
21#include <lemon/core.h>
22#include <lemon/concepts/graph.h>
23#include <lemon/dfs.h>
24#include <lemon/bfs.h>
25#include <test/test_tools.h>
26
27#include <vector>
28#include <set>
29
30namespace lemon
31{
32  namespace bits
33  {
34    namespace vf2
35    {
36      class AlwaysEq
37      {
38      public:
39        template<class T1, class T2>
40        bool operator()(T1, T2) const
41        {
42          return true;
43        }
44      };
45
46      template<class M1, class M2>
47      class MapEq
48      {
49        const M1 &_m1;
50        const M2 &_m2;
51      public:
52        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) {}
53        bool operator()(typename M1::Key k1, typename M2::Key k2) const
54        {
55          return _m1[k1] == _m2[k2];
56        }
57      };
58
59      template <class G>
60      class DfsLeaveOrder : public DfsVisitor<G>
61      {
62        const G &_g;
63        std::vector<typename G::Node> &_order;
64        int i;
65      public:
66        DfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
67          : i(countNodes(g)), _g(g), _order(order)
68        {}
69        void leave(const typename G::Node &node)
70        {
71          _order[--i]=node;
72        }
73      };
74
75      template <class G>
76      class BfsLeaveOrder : public BfsVisitor<G>
77      {
78        int i;
79        const G &_g;
80        std::vector<typename G::Node> &_order;
81      public:
82        BfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
83          : i(0), _g(g), _order(order)
84        {}
85        void process(const typename G::Node &node)
86        {
87          _order[i++]=node;
88        }
89      };
90    }
91  }
92
93  enum MappingType {
94    SUBGRAPH = 0,
95    INDUCED = 1,
96    ISOMORPH = 2
97  };
98
99  template<class G1, class G2, class I, class NEq = bits::vf2::AlwaysEq >
100  class Vf2
101  {
102    //Current depth in the DFS tree.
103    int _depth;
104    //Functor with bool operator()(G1::Node,G2::Node), which returns 1
105    //if and only if the 2 nodes are equivalent.
106    NEq _nEq;
107
108    typename G2::template NodeMap<int> _conn;
109    //Current matching. We index it by the nodes of g1, and match[v] is
110    //a node of g2.
111    I &_match;
112    //order[i] is the node of g1, for which we find a pair if depth=i
113    std::vector<typename G1::Node> order;
114    //currEdgeIts[i] is an edge iterator, witch is last used in the ith
115    //depth to find a pair for order[i].
116    std::vector<typename G2::IncEdgeIt> currEdgeIts;
117    //The small graph.
118    const G1 &_g1;
119    //The big graph.
120    const G2 &_g2;
121    //lookup tables for cut the searchtree
122    typename G1::template NodeMap<int> rNew1t,rInOut1t;
123
124    MappingType _mapping_type;
125
126    //cut the search tree
127    template<MappingType MT>
128    bool cut(const typename G1::Node n1,const typename G2::Node n2) const
129    {
130      int rNew2=0,rInOut2=0;
131      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
132        {
133          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
134          if(_conn[currNode]>0)
135            ++rInOut2;
136          else if(MT!=SUBGRAPH&&_conn[currNode]==0)
137            ++rNew2;
138        }
139      switch(MT)
140        {
141        case INDUCED:
142          return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
143        case SUBGRAPH:
144          return rInOut1t[n1]<=rInOut2;
145        case ISOMORPH:
146          return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
147        default:
148          return false;
149        }
150    }
151
152    template<MappingType MT>
153    bool feas(const typename G1::Node n1,const typename G2::Node n2)
154    {
155      if(!_nEq(n1,n2))
156        return 0;
157
158      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
159        {
160          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
161          if(_match[currNode]!=INVALID)
162            --_conn[_match[currNode]];
163        }
164      bool isIso=1;
165      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
166        {
167          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
168          if(_conn[currNode]<-1)
169            ++_conn[currNode];
170          else if(MT!=SUBGRAPH&&_conn[currNode]==-1)
171            {
172              isIso=0;
173              break;
174            }
175        }
176
177      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
178        {
179          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
180          if(_match[currNode]!=INVALID&&_conn[_match[currNode]]!=-1)
181            {
182              switch(MT)
183                {
184                case INDUCED:
185                case ISOMORPH:
186                  isIso=0;
187                  break;
188                case SUBGRAPH:
189                  if(_conn[_match[currNode]]<-1)
190                    isIso=0;
191                  break;
192                }
193              _conn[_match[currNode]]=-1;
194            }
195        }
196      return isIso&&cut<MT>(n1,n2);
197    }
198
199    void addPair(const typename G1::Node n1,const typename G2::Node n2)
200    {
201      _conn[n2]=-1;
202      _match.set(n1,n2);
203      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
204        if(_conn[_g2.oppositeNode(n2,e2)]!=-1)
205          ++_conn[_g2.oppositeNode(n2,e2)];
206    }
207
208    void subPair(const typename G1::Node n1,const typename G2::Node n2)
209    {
210      _conn[n2]=0;
211      _match.set(n1,INVALID);
212      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
213        {
214          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
215          if(_conn[currNode]>0)
216            --_conn[currNode];
217          else if(_conn[currNode]==-1)
218            ++_conn[n2];
219        }
220    }
221
222    void setOrder()//we will find pairs for the nodes of g1 in this order
223    {
224      // bits::vf2::DfsLeaveOrder<G1> v(_g1,order);
225      //   DfsVisit<G1,bits::vf2::DfsLeaveOrder<G1> >dfs(_g1, v);
226      //   dfs.run();
227
228      //it is more efficient in practice than DFS
229      bits::vf2::BfsLeaveOrder<G1> v(_g1,order);
230      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> >bfs(_g1, v);
231      bfs.run();
232    }
233
234  public:
235
236    template<MappingType MT>
237    bool extMatch()
238    {
239      while(_depth>=0)
240        {
241          //there are not nodes in g1, which has not pair in g2.
242          if(_depth==static_cast<int>(order.size()))
243            {
244              --_depth;
245              return true;
246            }
247          //the node of g2, which neighbours are the candidates for
248          //the pair of order[_depth]
249          typename G2::Node currPNode;
250          if(currEdgeIts[_depth]==INVALID)
251            {
252              typename G1::IncEdgeIt fstMatchedE(_g1,order[_depth]);
253              //if _match[order[_depth]]!=INVALID, we dont use
254              //fstMatchedE
255              if(_match[order[_depth]]==INVALID)
256                for(; fstMatchedE!=INVALID &&
257                      _match[_g1.oppositeNode(order[_depth],
258                                              fstMatchedE)]==INVALID;
259                    ++fstMatchedE) ; //find fstMatchedE
260              if(fstMatchedE==INVALID||_match[order[_depth]]!=INVALID)
261                {
262                  //We did not find an covered neighbour, this means
263                  //the graph is not connected(or _depth==0).  Every
264                  //uncovered(and there are some other properties due
265                  //to the spec. problem types) node of g2 is
266                  //candidate.  We can read the iterator of the last
267                  //tryed node from the match if it is not the first
268                  //try(match[order[_depth]]!=INVALID)
269                  typename G2::NodeIt n2(_g2);
270                  //if its not the first try
271                  if(_match[order[_depth]]!=INVALID)
272                    {
273                      n2=++typename G2::NodeIt(_g2,_match[order[_depth]]);
274                      subPair(order[_depth],_match[order[_depth]]);
275                    }
276                  for(; n2!=INVALID; ++n2)
277                    if(MT!=SUBGRAPH&&_conn[n2]==0)
278                      {
279                        if(feas<MT>(order[_depth],n2))
280                          break;
281                      }
282                    else if(MT==SUBGRAPH&&_conn[n2]>=0)
283                      if(feas<MT>(order[_depth],n2))
284                        break;
285                  // n2 is the next candidate
286                  if(n2!=INVALID)
287                    {
288                      addPair(order[_depth],n2);
289                      ++_depth;
290                    }
291                  else // there is no more candidate
292                    --_depth;
293                  continue;
294                }
295              else
296                {
297                  currPNode=_match[_g1.oppositeNode(order[_depth],fstMatchedE)];
298                  currEdgeIts[_depth]=typename G2::IncEdgeIt(_g2,currPNode);
299                }
300            }
301          else
302            {
303              currPNode=_g2.oppositeNode(_match[order[_depth]],
304                                         currEdgeIts[_depth]);
305              subPair(order[_depth],_match[order[_depth]]);
306              ++currEdgeIts[_depth];
307            }
308          for(; currEdgeIts[_depth]!=INVALID; ++currEdgeIts[_depth])
309            {
310              const typename G2::Node currNode =
311                _g2.oppositeNode(currPNode, currEdgeIts[_depth]);
312              //if currNode is uncovered
313              if(_conn[currNode]>0&&feas<MT>(order[_depth],currNode))
314                {
315                  addPair(order[_depth],currNode);
316                  break;
317                }
318            }
319          currEdgeIts[_depth]==INVALID?--_depth:++_depth;
320        }
321      return false;
322    }
323
324    //calc. the lookup table for cut the searchtree
325    void setRNew1tRInOut1t()
326    {
327      typename G1::template NodeMap<int> tmp(_g1,0);
328      for(unsigned int i=0; i<order.size(); ++i)
329        {
330          tmp[order[i]]=-1;
331          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
332            {
333              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
334              if(tmp[currNode]>0)
335                ++rInOut1t[order[i]];
336              else if(tmp[currNode]==0)
337                ++rNew1t[order[i]];
338            }
339          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
340            {
341              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
342              if(tmp[currNode]!=-1)
343                ++tmp[currNode];
344            }
345        }
346    }
347  public:
348    Vf2(const G1 &g1, const G2 &g2,I &i, const NEq &nEq = NEq() ) :
349      _nEq(nEq), _conn(g2,0), _match(i), order(countNodes(g1)),
350      currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNew1t(g1,0),
351      rInOut1t(g1,0), _mapping_type(SUBGRAPH)
352    {
353      _depth=0;
354      setOrder();
355      setRNew1tRInOut1t();
356    }
357
358    MappingType mappingType() const { return _mapping_type; }
359    void mappingType(MappingType m_type) { _mapping_type = m_type; }
360
361    bool find()
362    {
363      switch(_mapping_type)
364        {
365        case SUBGRAPH:
366          return extMatch<SUBGRAPH>();
367        case INDUCED:
368          return extMatch<INDUCED>();
369        case ISOMORPH:
370          return extMatch<ISOMORPH>();
371        default:
372          return false;
373        }
374    }
375  };
376
377  template<class G1, class G2>
378  class Vf2WizardBase
379  {
380  protected:
381    typedef G1 Graph1;
382    typedef G2 Graph2;
383
384    const G1 &_g1;
385    const G2 &_g2;
386
387    MappingType _mapping_type;
388
389    typedef typename G1::template NodeMap<typename G2::Node> Mapping;
390    bool _local_mapping;
391    void *_mapping;
392    void createMapping()
393    {
394      _mapping = new Mapping(_g1);
395    }
396
397    typedef bits::vf2::AlwaysEq NodeEq;
398    NodeEq _node_eq;
399
400    Vf2WizardBase(const G1 &g1,const G2 &g2)
401      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) {}
402  };
403
404  template<class TR>
405  class Vf2Wizard : public TR
406  {
407    typedef TR Base;
408    typedef typename TR::Graph1 Graph1;
409    typedef typename TR::Graph2 Graph2;
410
411    typedef typename TR::Mapping Mapping;
412    typedef typename TR::NodeEq NodeEq;
413
414    using TR::_g1;
415    using TR::_g2;
416    using TR::_mapping_type;
417    using TR::_mapping;
418    using TR::_node_eq;
419
420  public:
421    ///Copy constructor
422    Vf2Wizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {}
423
424    ///Copy constructor
425    Vf2Wizard(const Base &b) : Base(b) {}
426
427
428    template<class T>
429    struct SetMappingBase : public Base {
430      typedef T Mapping;
431      SetMappingBase(const Base &b) : Base(b) {}
432    };
433
434    ///\brief \ref named-templ-param "Named parameter" for setting
435    ///the mapping.
436    ///
437    ///\ref named-templ-param "Named parameter" function for setting
438    ///the map that stores the found embedding.
439    template<class T>
440    Vf2Wizard< SetMappingBase<T> > mapping(const T &t)
441    {
442      Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
443      Base::_local_mapping = false;
444      return Vf2Wizard<SetMappingBase<T> >(*this);
445    }
446
447    template<class NE>
448    struct SetNodeEqBase : public Base {
449      typedef NE NodeEq;
450      NodeEq _node_eq;
451      SetNodeEqBase(const Base &b, const NE &node_eq)
452        : Base(b), _node_eq(node_eq) {}
453    };
454
455    ///\brief \ref named-templ-param "Named parameter" for setting
456    ///the node equivalence relation.
457    ///
458    ///\ref named-templ-param "Named parameter" function for setting
459    ///the equivalence relation between the nodes.
460    template<class T>
461    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq)
462    {
463      return Vf2Wizard<SetNodeEqBase<T> >(SetNodeEqBase<T>(*this,node_eq));
464    }
465
466    ///\brief \ref named-templ-param "Named parameter" for setting
467    ///the node labels.
468    ///
469    ///\ref named-templ-param "Named parameter" function for setting
470    ///the node labels defining equivalence relation between them.
471    template<class M1, class M2>
472    Vf2Wizard< SetNodeEqBase<bits::vf2::MapEq<M1,M2> > >
473    nodeLabels(const M1 &m1,const M2 &m2)
474    {
475      return nodeEq(bits::vf2::MapEq<M1,M2>(m1,m2));
476    }
477
478    Vf2Wizard<Base> &mappingType(MappingType m_type)
479    {
480      _mapping_type = m_type;
481      return *this;
482    }
483
484    Vf2Wizard<Base> &induced()
485    {
486      _mapping_type = INDUCED;
487      return *this;
488    }
489
490    Vf2Wizard<Base> &iso()
491    {
492      _mapping_type = ISOMORPH;
493      return *this;
494    }
495
496    bool run()
497    {
498      if(Base::_local_mapping)
499        Base::createMapping();
500
501      Vf2<Graph1, Graph2, Mapping, NodeEq >
502        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);
503
504      alg.mappingType(_mapping_type);
505
506      bool ret = alg.find();
507
508      if(Base::_local_mapping)
509        delete reinterpret_cast<Mapping*>(_mapping);
510
511      return ret;
512    }
513  };
514
515  template<class G1, class G2>
516  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2)
517  {
518    return Vf2Wizard<Vf2WizardBase<G1,G2> >(g1,g2);
519  }
520
521}
522
523#endif
Note: See TracBrowser for help on using the repository browser.