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