// morpho/cdl/bgl_expansions/algorithms/ullmann.hpp header file//
// Copyright (c) 2003-2008 Vladimir J. Sykora 
// Copyright (c) 2007-2008 Vladimir J. Sykora and NCU Studies Ltd
// Modifications by Greg Landrum, January 2009
//
//*****************************************************************************
// Permission is hereby granted, free of charge, to any person or organization
// obtaining a copy of the software and accompanying documentation covered by
// this license (the "Software") to use, reproduce, display, distribute,
// execute, and transmit the Software, and to prepare derivative works of the
// Software, and to permit third-parties to whom the Software is furnished to
// do so, all subject to the following:
// 
// The copyright notices in the Software and this entire statement, including
// the above license grant, this restriction and the following disclaimer,
// must be included in all copies of the Software, in whole or in part, and
// all derivative works of the Software, unless such copies or derivative
// works are solely in the form of machine-executable object code generated by
// a source language processor.
// 
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT
// SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE
// FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE,
// ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.

//*****************************************************************************

//-----------------------------------------------------------------------------

// Algorithm due to JR Ullmann. "An Algorithm for Subgraph
// Isomorphism". Journal of the Association for
// Computing Machinery, Vol 23, No.1, January 1976, pp 31-42.

#ifndef MORPHO_CDL_BGL_EXP_ULLMANN_HPP
#define MORPHO_CDL_BGL_EXP_ULLMANN_HPP

// -------- boost
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/graph/graph_traits.hpp>
// -------- std
#include <utility>    // for std::pair<>
#include <algorithm>

namespace boost {
  namespace detail {

    template <class Graph, class UblasMatrix, class BackInsertionSequence, class EdgeLabeling>
    bool forward_checking(const Graph& g1, const Graph& g2,
                          UblasMatrix& M, size_t count, BackInsertionSequence& F, size_t num_vert_g1,
                          size_t num_vert_g2, EdgeLabeling& edge_labeling){
      typedef std::pair<typename graph_traits<Graph>::edge_descriptor,bool> edge_presence;
      typename BackInsertionSequence::iterator   fi, fi_end=F.end();
      for (size_t k=count+1; k<num_vert_g1; ++k) {
        for (size_t l=0; l<num_vert_g2; ++l) {
          if(!M(k,l)) continue;
          // check mapping:
          fi=F.begin();
          while(fi!=fi_end) {
            bool flag1(0), flag1_1(0), flag2(0), flag2_1(0);
            edge_presence ep1=edge(k,fi->first,g1);
            if(ep1.second) {
              flag1_1=true;
              edge_presence ep2=edge(l,fi->second,g2);
              if(ep2.second) flag1=edge_labeling(ep1.first,ep2.first);
            }
            if(flag1_1 && flag1) {
              M(k,l)=1;
              ++fi;
              continue;
            }
            edge_presence ep2=edge(l,fi->second,g2);
            if(ep2.second) {
              flag2_1=true;
              ep1=edge(k,fi->first,g1);
              if(ep1.second){
                flag2=edge_labeling(ep1.first,ep2.first);
              } else {
                // ring closed in main structure, not closed in query. This should pass
                flag2=true;
              }
            }
            if(flag2_1 && flag2) {  //if one edge exists, there must be a mapping
              M(k,l)=1;
            }
            else if ( !flag1_1 && !flag2_1 ) { // or both edges are not present
              M(k,l)=1;
            }
            else M(k,l)=0;         // if not, there's no mapping
            ++fi;
          }
        }
      }
      // TODO: change the data structure of the M matrix to sparse matrix. This wouldn't be neccessary
      size_t cero_row(0);
      for (size_t k=0; k<num_vert_g1; ++k) {
        for (size_t l=0; l<num_vert_g2; ++l) {
          if(M(k,l)) break;
          else ++cero_row;
        }
        if(cero_row==num_vert_g2) return false; // if there is a cero row
        cero_row=0;
      }
      return true;
    }

    template <class Graph, class EdgeLabeling, class UblasMatrix, class BackInsertionSequence>
    bool backtrack(const Graph& g1, const Graph& g2, size_t count,
                   const UblasMatrix& M, BackInsertionSequence& F, const size_t num_vert_g1,
                   const size_t num_vert_g2, EdgeLabeling& edge_labeling){
      if(count==num_vert_g1) return true;
      for (size_t i=0; i<num_vert_g2; ++i) {
        if(M(count,i)) {
          F.push_back(std::make_pair(count,i));
          UblasMatrix M_prime(M);
          for (size_t m=count+1; m<num_vert_g1; ++m) {
            M_prime(m,i)=0;
          }
          
          if (forward_checking(g1,g2,M_prime,count,F,num_vert_g1,num_vert_g2,edge_labeling)){
            if(backtrack(g1,g2,count+1,M_prime,F,num_vert_g1,num_vert_g2,edge_labeling)){
              return true;
            }
          }
          F.erase(std::remove(F.begin(),F.end(),std::make_pair(count,i)), F.end());
        }
      }
      return false;
    }

    template <class Graph, class EdgeLabeling, class UblasMatrix, class DoubleBackInsertionSequence>
    void backtrack_all(const Graph& g1, const Graph& g2, size_t count,
                       const UblasMatrix& M, DoubleBackInsertionSequence& FF, const size_t num_vert_g1,
                       const size_t num_vert_g2, EdgeLabeling& edge_labeling) {
      if(count==num_vert_g1) return;
      DoubleBackInsertionSequence holdFF;
      holdFF.insert(holdFF.begin(),FF.begin(),FF.end());
      FF.clear();
      for (size_t i=0; i<num_vert_g2; ++i) {
        if(M(count,i)) {
          DoubleBackInsertionSequence tFF;

          UblasMatrix M_prime(M);
          for (size_t m=count+1; m<num_vert_g1; ++m){
            M_prime(m,i)=0;
          }

          if(holdFF.size()){
            for(typename DoubleBackInsertionSequence::const_iterator iter=holdFF.begin();
                iter!=holdFF.end();++iter){
              typename DoubleBackInsertionSequence::value_type F=*iter;
              F.push_back(std::make_pair(count,i));
              if (forward_checking(g1,g2,M_prime,count,F,num_vert_g1,num_vert_g2,edge_labeling)){
                tFF.push_back(F);
              }
            }
          } else {
            typename DoubleBackInsertionSequence::value_type F;
            F.push_back(std::make_pair(count,i));
            if (forward_checking(g1,g2,M_prime,count,F,num_vert_g1,num_vert_g2,edge_labeling)){
              tFF.push_back(F);
            }
          }
          backtrack_all(g1,g2,count+1,M_prime,tFF,num_vert_g1,num_vert_g2,edge_labeling);
          if(tFF.size()){
            for(typename DoubleBackInsertionSequence::const_iterator iter=tFF.begin();
                iter!=tFF.end();++iter){
              FF.push_back(*iter);
            }
          }
        }
      }
    }

    template <  class Graph
                , class VertexLabeling    // binary predicate
                ,  class UblasMatrix
                >
    void prepareM(const Graph& g1, const Graph& g2,
                  VertexLabeling& vertex_labeling,UblasMatrix &M){
      size_t rows(num_vertices(g1));
      size_t cols(num_vertices(g2));
      M.resize(rows,cols);
      // initialize the matrix:
      for (size_t i=0; i<rows; ++i) {
        for (size_t j=0; j<cols; ++j) {
          if(out_degree(j,g2)>=out_degree(i,g1) && vertex_labeling(i,j)) {
            M(i,j)=1;
          }
          else M(i,j)=0;
        }
      }
    }
  }  // namespace detail


  // test if g1 is a subgraph of g2. mapped vertices are returned in F
  // Mapping : first: g1 vertices, second : g2 vertices
  // O( num_vertices(g1)! num_vertices(g1) ^ 3 )
  // This function doesnt Work with filtered graphs!
  // The size of F doesn't indicate match!
  template <  class Graph
              , class VertexLabeling    // binary predicate
              , class EdgeLabeling      // binary predicate
              , class BackInsertionSequence   // contains std::pair<vertex_descriptor,vertex_descriptor>
              >
  bool ullmann(const Graph& g1, const Graph& g2,
               VertexLabeling& vertex_labeling, EdgeLabeling& edge_labeling, BackInsertionSequence& F) {
    typedef ::boost::numeric::ublas::matrix<int>   matrix_t;
    size_t rows(num_vertices(g1));
    size_t cols(num_vertices(g2));
    matrix_t M;
    detail::prepareM(g1,g2,vertex_labeling,M);
    size_t    count(0);
    return detail::backtrack(g1,g2,count,M,F,rows,cols,edge_labeling);
  }

  // test if g1 is a subgraph of g2.
  // F returns all mappings of g1 in g2. mapping in separate containers
  template <  class Graph
              , class VertexLabeling    // binary predicate
              , class EdgeLabeling      // binary predicate
              , class DoubleBackInsertionSequence   // contains a back insertion sequence
              >
  bool ullmann_all(const Graph& g1, const Graph& g2,
                   VertexLabeling& vertex_labeling, EdgeLabeling& edge_labeling, DoubleBackInsertionSequence& F) {
    typedef ::boost::numeric::ublas::matrix<int>   matrix_t;
    size_t rows(num_vertices(g1));
    size_t cols(num_vertices(g2));
    matrix_t M;
    detail::prepareM(g1,g2,vertex_labeling,M);
    size_t    count(0);
    detail::backtrack_all(g1,g2,count,M,F,rows,cols,edge_labeling);
    return !F.empty();
  }

}  // namespace boost


#endif // MORPHO_CDL_BGL_EXP_ULLMANN_HPP