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

/* -*- mode: C++; indent-tabs-mode: nil; -*-
 *
 * This file is a part of LEMON, a generic C++ optimization library.
 *
 * Copyright (C) 2015
 * EMAXA Kutato-fejleszto Kft. (EMAXA Research Ltd.)
 *
 * Permission to use, modify and distribute this software is granted
 * provided that this copyright notice appears in all copies. For
 * precise terms see the accompanying LICENSE file.
 *
 * This software is provided "AS IS" with no warranty of any kind,
 * express or implied, and with no claim as to its suitability for any
 * purpose.
 *
 */

#ifndef LEMON_VF2_H
#define LEMON_VF2_H

#include <lemon/core.h>
#include <lemon/concepts/graph.h>
#include <lemon/dfs.h>
#include <lemon/bfs.h>
#include <test/test_tools.h>

#include <vector>
#include <set>

namespace lemon
{
  namespace bits
  {
    namespace vf2
    {
      class AlwaysEq
      {
      public:
        template<class T1, class T2>
        bool operator()(T1, T2) const
        {
          return true;
        }
      };

      template<class M1, class M2>
      class MapEq
      {
        const M1 &_m1;
        const M2 &_m2;
      public:
        MapEq(const M1 &m1, const M2 &m2) : _m1(m1), _m2(m2) {}
        bool operator()(typename M1::Key k1, typename M2::Key k2) const
        {
          return _m1[k1] == _m2[k2];
        }
      };

      template <class G>
      class DfsLeaveOrder : public DfsVisitor<G>
      {
        const G &_g;
        std::vector<typename G::Node> &_order;
        int i;
      public:
        DfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
          : i(countNodes(g)), _g(g), _order(order)
        {}
        void leave(const typename G::Node &node)
        {
          _order[--i]=node;
        }
      };

      template <class G>
      class BfsLeaveOrder : public BfsVisitor<G>
      {
        int i;
        const G &_g;
        std::vector<typename G::Node> &_order;
      public:
        BfsLeaveOrder(const G &g, std::vector<typename G::Node> &order)
          : i(0), _g(g), _order(order)
        {}
        void process(const typename G::Node &node)
        {
          _order[i++]=node;
        }
      };
    }
  }

  enum MappingType {
    SUBGRAPH = 0,
    INDUCED = 1,
    ISOMORPH = 2
  };

  template<class G1, class G2, class I, class NEq = bits::vf2::AlwaysEq >
  class Vf2
  {
    //Current depth in the DFS tree.
    int _depth;
    //Functor with bool operator()(G1::Node,G2::Node), which returns 1
    //if and only if the 2 nodes are equivalent.
    NEq _nEq;

    typename G2::template NodeMap<int> _conn;
    //Current matching. We index it by the nodes of g1, and match[v] is
    //a node of g2.
    I &_match;
    //order[i] is the node of g1, for which we find a pair if depth=i
    std::vector<typename G1::Node> order;
    //currEdgeIts[i] is an edge iterator, witch is last used in the ith
    //depth to find a pair for order[i].
    std::vector<typename G2::IncEdgeIt> currEdgeIts;
    //The small graph.
    const G1 &_g1;
    //The big graph.
    const G2 &_g2;
    //lookup tables for cut the searchtree
    typename G1::template NodeMap<int> rNew1t,rInOut1t;

    MappingType _mapping_type;

    //cut the search tree
    template<MappingType MT>
    bool cut(const typename G1::Node n1,const typename G2::Node n2) const
    {
      int rNew2=0,rInOut2=0;
      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
        {
          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
          if(_conn[currNode]>0)
            ++rInOut2;
          else if(MT!=SUBGRAPH&&_conn[currNode]==0)
            ++rNew2;
        }
      switch(MT)
        {
        case INDUCED:
          return rInOut1t[n1]<=rInOut2&&rNew1t[n1]<=rNew2;
        case SUBGRAPH:
          return rInOut1t[n1]<=rInOut2;
        case ISOMORPH:
          return rInOut1t[n1]==rInOut2&&rNew1t[n1]==rNew2;
        default:
          return false;
        }
    }

    template<MappingType MT>
    bool feas(const typename G1::Node n1,const typename G2::Node n2)
    {
      if(!_nEq(n1,n2))
        return 0;

      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
        {
          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
          if(_match[currNode]!=INVALID)
            --_conn[_match[currNode]];
        }
      bool isIso=1;
      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
        {
          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
          if(_conn[currNode]<-1)
            ++_conn[currNode];
          else if(MT!=SUBGRAPH&&_conn[currNode]==-1)
            {
              isIso=0;
              break;
            }
        }

      for(typename G1::IncEdgeIt e1(_g1,n1); e1!=INVALID; ++e1)
        {
          const typename G1::Node currNode=_g1.oppositeNode(n1,e1);
          if(_match[currNode]!=INVALID&&_conn[_match[currNode]]!=-1)
            {
              switch(MT)
                {
                case INDUCED:
                case ISOMORPH:
                  isIso=0;
                  break;
                case SUBGRAPH:
                  if(_conn[_match[currNode]]<-1)
                    isIso=0;
                  break;
                }
              _conn[_match[currNode]]=-1;
            }
        }
      return isIso&&cut<MT>(n1,n2);
    }

    void addPair(const typename G1::Node n1,const typename G2::Node n2)
    {
      _conn[n2]=-1;
      _match.set(n1,n2);
      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
        if(_conn[_g2.oppositeNode(n2,e2)]!=-1)
          ++_conn[_g2.oppositeNode(n2,e2)];
    }

    void subPair(const typename G1::Node n1,const typename G2::Node n2)
    {
      _conn[n2]=0;
      _match.set(n1,INVALID);
      for(typename G2::IncEdgeIt e2(_g2,n2); e2!=INVALID; ++e2)
        {
          const typename G2::Node currNode=_g2.oppositeNode(n2,e2);
          if(_conn[currNode]>0)
            --_conn[currNode];
          else if(_conn[currNode]==-1)
            ++_conn[n2];
        }
    }

    void setOrder()//we will find pairs for the nodes of g1 in this order
    {
      // bits::vf2::DfsLeaveOrder<G1> v(_g1,order);
      //   DfsVisit<G1,bits::vf2::DfsLeaveOrder<G1> >dfs(_g1, v);
      //   dfs.run();

      //it is more efficient in practice than DFS
      bits::vf2::BfsLeaveOrder<G1> v(_g1,order);
      BfsVisit<G1,bits::vf2::BfsLeaveOrder<G1> >bfs(_g1, v);
      bfs.run();
    }

  public:

    template<MappingType MT>
    bool extMatch()
    {
      while(_depth>=0)
        {
          //there are not nodes in g1, which has not pair in g2.
          if(_depth==static_cast<int>(order.size()))
            {
              --_depth;
              return true;
            }
          //the node of g2, which neighbours are the candidates for
          //the pair of order[_depth]
          typename G2::Node currPNode;
          if(currEdgeIts[_depth]==INVALID)
            {
              typename G1::IncEdgeIt fstMatchedE(_g1,order[_depth]);
              //if _match[order[_depth]]!=INVALID, we dont use
              //fstMatchedE
              if(_match[order[_depth]]==INVALID)
                for(; fstMatchedE!=INVALID &&
                      _match[_g1.oppositeNode(order[_depth],
                                              fstMatchedE)]==INVALID;
                    ++fstMatchedE) ; //find fstMatchedE
              if(fstMatchedE==INVALID||_match[order[_depth]]!=INVALID)
                {
                  //We did not find an covered neighbour, this means
                  //the graph is not connected(or _depth==0).  Every
                  //uncovered(and there are some other properties due
                  //to the spec. problem types) node of g2 is
                  //candidate.  We can read the iterator of the last
                  //tryed node from the match if it is not the first
                  //try(match[order[_depth]]!=INVALID)
                  typename G2::NodeIt n2(_g2);
                  //if its not the first try
                  if(_match[order[_depth]]!=INVALID)
                    {
                      n2=++typename G2::NodeIt(_g2,_match[order[_depth]]);
                      subPair(order[_depth],_match[order[_depth]]);
                    }
                  for(; n2!=INVALID; ++n2)
                    if(MT!=SUBGRAPH&&_conn[n2]==0)
                      {
                        if(feas<MT>(order[_depth],n2))
                          break;
                      }
                    else if(MT==SUBGRAPH&&_conn[n2]>=0)
                      if(feas<MT>(order[_depth],n2))
                        break;
                  // n2 is the next candidate
                  if(n2!=INVALID)
                    {
                      addPair(order[_depth],n2);
                      ++_depth;
                    }
                  else // there is no more candidate
                    --_depth;
                  continue;
                }
              else
                {
                  currPNode=_match[_g1.oppositeNode(order[_depth],fstMatchedE)];
                  currEdgeIts[_depth]=typename G2::IncEdgeIt(_g2,currPNode);
                }
            }
          else
            {
              currPNode=_g2.oppositeNode(_match[order[_depth]],
                                         currEdgeIts[_depth]);
              subPair(order[_depth],_match[order[_depth]]);
              ++currEdgeIts[_depth];
            }
          for(; currEdgeIts[_depth]!=INVALID; ++currEdgeIts[_depth])
            {
              const typename G2::Node currNode =
                _g2.oppositeNode(currPNode, currEdgeIts[_depth]);
              //if currNode is uncovered
              if(_conn[currNode]>0&&feas<MT>(order[_depth],currNode))
                {
                  addPair(order[_depth],currNode);
                  break;
                }
            }
          currEdgeIts[_depth]==INVALID?--_depth:++_depth;
        }
      return false;
    }

    //calc. the lookup table for cut the searchtree
    void setRNew1tRInOut1t()
    {
      typename G1::template NodeMap<int> tmp(_g1,0);
      for(unsigned int i=0; i<order.size(); ++i)
        {
          tmp[order[i]]=-1;
          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
            {
              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
              if(tmp[currNode]>0)
                ++rInOut1t[order[i]];
              else if(tmp[currNode]==0)
                ++rNew1t[order[i]];
            }
          for(typename G1::IncEdgeIt e1(_g1,order[i]); e1!=INVALID; ++e1)
            {
              const typename G1::Node currNode=_g1.oppositeNode(order[i],e1);
              if(tmp[currNode]!=-1)
                ++tmp[currNode];
            }
        }
    }
  public:
    Vf2(const G1 &g1, const G2 &g2,I &i, const NEq &nEq = NEq() ) :
      _nEq(nEq), _conn(g2,0), _match(i), order(countNodes(g1)),
      currEdgeIts(countNodes(g1),INVALID), _g1(g1), _g2(g2), rNew1t(g1,0),
      rInOut1t(g1,0), _mapping_type(SUBGRAPH)
    {
      _depth=0;
      setOrder();
      setRNew1tRInOut1t();
    }

    MappingType mappingType() const { return _mapping_type; }
    void mappingType(MappingType m_type) { _mapping_type = m_type; }

    bool find()
    {
      switch(_mapping_type)
        {
        case SUBGRAPH:
          return extMatch<SUBGRAPH>();
        case INDUCED:
          return extMatch<INDUCED>();
        case ISOMORPH:
          return extMatch<ISOMORPH>();
        default:
          return false;
        }
    }
  };

  template<class G1, class G2>
  class Vf2WizardBase
  {
  protected:
    typedef G1 Graph1;
    typedef G2 Graph2;

    const G1 &_g1;
    const G2 &_g2;

    MappingType _mapping_type;

    typedef typename G1::template NodeMap<typename G2::Node> Mapping;
    bool _local_mapping;
    void *_mapping;
    void createMapping()
    {
      _mapping = new Mapping(_g1);
    }

    typedef bits::vf2::AlwaysEq NodeEq;
    NodeEq _node_eq;

    Vf2WizardBase(const G1 &g1,const G2 &g2)
      : _g1(g1), _g2(g2), _mapping_type(SUBGRAPH), _local_mapping(true) {}
  };

  template<class TR>
  class Vf2Wizard : public TR
  {
    typedef TR Base;
    typedef typename TR::Graph1 Graph1;
    typedef typename TR::Graph2 Graph2;

    typedef typename TR::Mapping Mapping;
    typedef typename TR::NodeEq NodeEq;

    using TR::_g1;
    using TR::_g2;
    using TR::_mapping_type;
    using TR::_mapping;
    using TR::_node_eq;

  public:
    ///Copy constructor
    Vf2Wizard(const Graph1 &g1,const Graph2 &g2) : Base(g1,g2) {}

    ///Copy constructor
    Vf2Wizard(const Base &b) : Base(b) {}


    template<class T>
    struct SetMappingBase : public Base {
      typedef T Mapping;
      SetMappingBase(const Base &b) : Base(b) {}
    };

    ///\brief \ref named-templ-param "Named parameter" for setting
    ///the mapping.
    ///
    ///\ref named-templ-param "Named parameter" function for setting
    ///the map that stores the found embedding.
    template<class T>
    Vf2Wizard< SetMappingBase<T> > mapping(const T &t)
    {
      Base::_mapping=reinterpret_cast<void*>(const_cast<T*>(&t));
      Base::_local_mapping = false;
      return Vf2Wizard<SetMappingBase<T> >(*this);
    }

    template<class NE>
    struct SetNodeEqBase : public Base {
      typedef NE NodeEq;
      NodeEq _node_eq;
      SetNodeEqBase(const Base &b, const NE &node_eq)
        : Base(b), _node_eq(node_eq) {}
    };

    ///\brief \ref named-templ-param "Named parameter" for setting
    ///the node equivalence relation.
    ///
    ///\ref named-templ-param "Named parameter" function for setting
    ///the equivalence relation between the nodes.
    template<class T>
    Vf2Wizard< SetNodeEqBase<T> > nodeEq(const T &node_eq)
    {
      return Vf2Wizard<SetNodeEqBase<T> >(SetNodeEqBase<T>(*this,node_eq));
    }

    ///\brief \ref named-templ-param "Named parameter" for setting
    ///the node labels.
    ///
    ///\ref named-templ-param "Named parameter" function for setting
    ///the node labels defining equivalence relation between them.
    template<class M1, class M2>
    Vf2Wizard< SetNodeEqBase<bits::vf2::MapEq<M1,M2> > >
    nodeLabels(const M1 &m1,const M2 &m2)
    {
      return nodeEq(bits::vf2::MapEq<M1,M2>(m1,m2));
    }

    Vf2Wizard<Base> &mappingType(MappingType m_type)
    {
      _mapping_type = m_type;
      return *this;
    }

    Vf2Wizard<Base> &induced()
    {
      _mapping_type = INDUCED;
      return *this;
    }

    Vf2Wizard<Base> &iso()
    {
      _mapping_type = ISOMORPH;
      return *this;
    }

    bool run()
    {
      if(Base::_local_mapping)
        Base::createMapping();

      Vf2<Graph1, Graph2, Mapping, NodeEq >
        alg(_g1, _g2, *reinterpret_cast<Mapping*>(_mapping), _node_eq);

      alg.mappingType(_mapping_type);

      bool ret = alg.find();

      if(Base::_local_mapping)
        delete reinterpret_cast<Mapping*>(_mapping);

      return ret;
    }
  };

  template<class G1, class G2>
  Vf2Wizard<Vf2WizardBase<G1,G2> > vf2(const G1 &g1, const G2 &g2)
  {
    return Vf2Wizard<Vf2WizardBase<G1,G2> >(g1,g2);
  }

}

#endif