%% 

%% Copyright 2007, 2008, 2009 Elsevier Ltd

%% 

%% This file is part of the 'Elsarticle Bundle'.

%% ---------------------------------------------

%% 

%% It may be distributed under the conditions of the LaTeX Project Public

%% License, either version 1.2 of this license or (at your option) any

%% later version.  The latest version of this license is in

%%    http://www.latex-project.org/lppl.txt

%% and version 1.2 or later is part of all distributions of LaTeX

%% version 1999/12/01 or later.

%% 

%% The list of all files belonging to the 'Elsarticle Bundle' is

%% given in the file `manifest.txt'.

%% 

%% Template article for Elsevier's document class `elsarticle'

%% with numbered style bibliographic references

%% SP 2008/03/01

\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing

%% \documentclass[authoryear,preprint,review,12pt]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn

%% for a journal layout:

%% \documentclass[final,1p,times]{elsarticle}

%% \documentclass[final,1p,times,twocolumn]{elsarticle}

%% \documentclass[final,3p,times]{elsarticle}

%% \documentclass[final,3p,times,twocolumn]{elsarticle}

%% \documentclass[final,5p,times]{elsarticle}

%% \documentclass[final,5p,times,twocolumn]{elsarticle}

%% For including figures, graphicx.sty has been loaded in

%% elsarticle.cls. If you prefer to use the old commands

%% please give \usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols

\usepackage{amssymb}

%% The amsthm package provides extended theorem environments

%% \usepackage{amsthm}

%% The lineno packages adds line numbers. Start line numbering with

%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on

%% for the whole article with \linenumbers.

%% \usepackage{lineno}

\usepackage{amsmath}

%% \usepackage[pdftex]{graphicx}

\usepackage{pgfplots}

\pgfplotsset{width=9cm}

\pgfplotsset{compat=1.8}

\usepackage{caption}

\usepackage{subcaption} 

\usepackage{algorithm}

\usepackage{algpseudocode}

\usepackage{tikz}

\usepackage{amsthm,amssymb}

\renewcommand{\qedsymbol}{\rule{0.7em}{0.7em}}

\newtheorem{theorem}{Theorem}[subsection]

\newtheorem{corollary}{Corollary}[theorem]

\newtheorem{claim}[theorem]{Claim}

\newtheorem{definition}{Definition}[subsection]

\newtheorem{notation}{Notation}[subsection]

\newtheorem{example}{Example}[subsection]

\usetikzlibrary{decorations.markings}

\let\oldproofname=\proofname

%% \renewcommand{\proofname}{\rm\bf{Proof:}}

\captionsetup{font=normalsize}

\journal{Discrete Applied Mathematics}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;

%% use the tnotetext command for theassociated footnote;

%% use the fnref command within \author or \address for footnotes;

%% use the fntext command for theassociated footnote;

%% use the corref command within \author for corresponding author footnotes;

%% use the cortext command for theassociated footnote;

%% use the ead command for the email address,

%% and the form \ead[url] for the home page:

%% \title{Title\tnoteref{label1}}

%% \tnotetext[label1]{}

%% \author{Name\corref{cor1}\fnref{label2}}

%% \ead{email address}

%% \ead[url]{home page}

%% \fntext[label2]{}

%% \cortext[cor1]{}

%% \address{Address\fnref{label3}}

%% \fntext[label3]{}

\title{VF2++ --- An Improved Subgraph Isomorphism Algorithm}

%% use optional labels to link authors explicitly to addresses:

%% \author[label1,label2]{}

%% \address[label1]{}

%% \address[label2]{}

\author[egres,elte]{Alp{\'a}r J{\"u}ttner}

\ead{alpar@cs.elte.hu}

\author[elte]{P{\'e}ter Madarasi}

\ead{madarasip@caesar.elte.hu}

\address[egres]{MTA-ELTE Egerv{\'a}ry Research Group, Budapest, Hungary.}

\address[elte]{Department of Operations Research, E{\"o}tv{\"o}s Lor{\'a}nd University, Budapest, Hungary.}

\begin{abstract}

  This paper presents a largely improved version of the VF2 algorithm

  for the \emph{Subgraph Isomorphism Problem}. The improvements are

  twofold. Firstly, it is based on a new approach for determining the

  matching order of the nodes, and secondly, more efficient ---

  nevertheless easier to compute --- cutting rules are proposed. They

  together reduce the search space significantly.

  In addition to the usual \emph{Subgraph Isomorphism Problem}, the paper also

  presents specialized algorithms for the \emph{Induced Subgraph

    Isomorphism} and for the \emph{Graph Isomorphism Problems}.

  Finally, an extensive experimental evaluation is provided using a

  wide range of inputs, including both real-life biological and

  chemical datasets and standard randomly generated graph series. The

  results show major and consistent running time improvements over the

  other known methods.

  The C++ implementations of the algorithms are available open-source as 

  part of the LEMON graph and network optimization library.

\end{abstract}

\begin{keyword}

  Computational Biology, Subgraph Isomorphism Problem

%% keywords here, in the form: keyword \sep keyword

%% PACS codes here, in the form: \PACS code \sep code

%% MSC codes here, in the form: \MSC code \sep code

%% or \MSC[2008] code \sep code (2000 is the default)

\end{keyword}

\end{frontmatter}

%% \linenumbers

%% main text

\section{Introduction}

\label{sec:intro}

In the last decades, combinatorial structures, and especially graphs

have been considered with ever increasing interest, and applied to the

solution of several new and revised questions. The expressiveness,

the simplicity and the deep theoretical background

of graphs make it one of the most useful

modeling tool and appears constantly in several seemingly independent

fields, such as bioinformatics and chemistry.

Complex biological systems arise from the interaction and cooperation

of plenty of molecular components. Getting acquainted with the

structure of such systems at the molecular level is of primary

importance, since protein-protein interaction, DNA-protein

interaction, metabolic interaction, transcription factor binding,

neuronal networks, and hormone singling networks can be understood

this way.

Many chemical and biological structures can easily be modeled

as graphs, for instance, a molecular structure can be

considered as a graph, whose nodes correspond to atoms and whose

edges to chemical bonds. The similarity and dissimilarity of

objects corresponding to nodes are incorporated to the model

by \emph{node labels}. Understanding such networks basically

requires finding specific subgraphs, thus it calls for efficient

graph matching algorithms.

Other real-world fields related to some

variants of graph matching include pattern recognition

and machine vision~\cite{HorstBunkeApplications}, symbol recognition~\cite{CordellaVentoSymbolRecognition}, and face identification~\cite{JianzhuangYongFaceIdentification}.  \\

Subgraph and induced subgraph matching problems are known to be

NP-Complete~\cite{SubgraphNPC}, while the graph isomorphism problem is

one of the few problems in NP neither known to be in P nor

NP-Complete. Although polynomial-time isomorphism algorithms are known

for various graph classes, like trees and planar

graphs~\cite{PlanarGraphIso}, bounded valence

graphs~\cite{BondedDegGraphIso}, interval graphs~\cite{IntervalGraphIso}

or permutation graphs~\cite{PermGraphIso}. Furthermore, an FPT algorithm has also been presented for the colored hypergraph isomorphism problem in~\cite{ColoredHiperGraphIso}.

In the following, some algorithms which do not need any restrictions on the graphs are

summarized. Even though,

an overall polynomial behavior is not expectable from such an

alternative, they may often have good practical performance, in fact,

they might be the best choice in practice even on a graph class for which polynomial

algorithm is known.

The first practically usable approach was due to

\emph{Ullmann}~\cite{Ullmann}, which is a commonly used algorithm based on depth-first

search with a complex heuristic for reducing the

number of visited states. A major problem is its $\Theta(n^3)$ space

complexity, which makes it impractical in the case of big sparse

graphs.

In a recent paper, Ullmann~\cite{UllmannBit} presents an

improved version of this algorithm based on a bit-vector solution for

the binary Constraint Satisfaction Problem.

The \emph{Nauty} algorithm~\cite{Nauty} transforms the two graphs to

a canonical form before starting to look for an isomorphism. It has

been considered as one of the fastest graph isomorphism algorithms,

although graph categories were shown in which it takes exponentially

many steps. This algorithm handles only the graph isomorphism problem.

The \emph{LAD} algorithm~\cite{Lad} uses a depth-first search

strategy and formulates the matching as a Constraint Satisfaction

Problem to prune the search tree. The constraints are that the mapping

has to be invective and edge-preserving, hence it is possible to

handle new matching types as well.

The \emph{RI} algorithm~\cite{RI} and its variations are based on a

state space representation. After reordering the nodes of the graphs,

it uses some fast executable heuristic checks without using any

complex pruning rules. It seems to run really efficiently on graphs

coming from biology, and won the International Contest on Pattern

Search in Biological Databases~\cite{Content}.

Currently, the most commonly used algorithm is the

\emph{VF2}~\cite{VF2}, an improved version of \emph{VF}~\cite{VF}, which was

designed for solving pattern matching and computer vision problems,

and has been one of the best overall algorithms for more than a

decade. Although, it is not as fast as some of the new specialized algorithms, it is still widely used due to its simplicity and space efficiency. VF2 uses

a state space representation and checks specific conditions in each state

to prune the search tree.

Meanwhile, another variant called \emph{VF2~Plus}~\cite{VF2Plus} has

been published. It is considered to be as efficient as the RI

algorithm and has a strictly better behavior on large graphs.  The

main idea of VF2~Plus is to precompute a heuristic node order of the graph to be embedded, on which VF2 works more efficiently.

This paper introduces \emph{VF2++}, a new further improved algorithm

for the graph and (induced) subgraph isomorphism problems. It is based on

efficient cutting rules and determines a node order in which VF2 runs

significantly faster on practical inputs.

The rest of the paper is structured as

follows. Section~\ref{sec:ProbStat} defines the exact problems to be

solved, Section~\ref{sec:VF2Alg} provides a description of VF2. Based

on that, Section~\ref{sec:VF2ppAlg} introduces VF2++. Some technical

details necessary for an efficient implementation are discussed in

Section~\ref{sec:VF2ppImpl}. Finally, Section~\ref{sec:ExpRes}

provides a detailed experimental evaluation of VF2++ and its comparison

to the state-of-the-art algorithm.

It must also be mentioned that the C++ implementations of the

algorithms have been made available for evaluation and use under an

open-source license as a part of

LEMON~\cite{o11:_lemon_open_sourc_c_graph_templ_librar} graph library.

\section{Problem Statement}\label{sec:ProbStat}

This section provides a formal description of the problems to be

solved.

\subsection{Definitions}

Throughout the paper $G_{1}=(V_{1}, E_{1})$ and

$G_{2}=(V_{2}, E_{2})$ denote two undirected graphs.

\begin{definition}

$\mathcal{L}: (V_{1}\cup V_{2}) \longrightarrow K$ is a \textbf{node

    label function}, where K is an arbitrary set. The elements in K

  are the \textbf{node labels}. Two nodes, u and v are said to be

  \textbf{equivalent} if $\mathcal{L}(u)=\mathcal{L}(v)$.

\end{definition}

For the sake of simplicity, the graph, subgraph and induced subgraph

isomorphisms are defined in a more general way.

\begin{definition}\label{sec:ismorphic}

$G_{1}$ and $G_{2}$ are \textbf{isomorphic} (by the node label $\mathcal{L}$) if $\exists \mathfrak{m}:

  V_{1} \longrightarrow V_{2}$ bijection, for which the

  following is true:

\begin{center}

$\forall u\in{V_{1}} : \mathcal{L}(u)=\mathcal{L}(\mathfrak{m}(u))$ and\\

$\forall u,v\in{V_{1}} : (u,v)\in{E_{1}} \Leftrightarrow (\mathfrak{m}(u),\mathfrak{m}(v))\in{E_{2}}$

\end{center}

\end{definition}

\begin{definition}

$G_{1}$ is a \textbf{subgraph} of $G_{2}$ (by the node label $\mathcal{L}$) if $\exists \mathfrak{m}:

  V_{1}\longrightarrow V_{2}$ injection, for which the

  following is true:

\begin{center}

$\forall u\in{V_{1}} : \mathcal{L}(u)=\mathcal{L}(\mathfrak{m}(u))$ and\\

$\forall u,v \in{V_{1}} : (u,v)\in{E_{1}} \Rightarrow (\mathfrak{m}(u),\mathfrak{m}(v))\in E_{2}$

\end{center}

\end{definition}

\begin{definition} 

$G_{1}$ is an \textbf{induced subgraph} of $G_{2}$ (by the node label $\mathcal{L}$) if $\exists

  \mathfrak{m}: V_{1}\longrightarrow V_{2}$ injection, for which the

  following is true:

\begin{center}

$\forall u\in{V_{1}} : \mathcal{L}(u)=\mathcal{L}(\mathfrak{m}(u))$ and

$\forall u,v \in{V_{1}} : (u,v)\in{E_{1}} \Leftrightarrow

  (\mathfrak{m}(u),\mathfrak{m}(v))\in E_{2}$

\end{center}

\end{definition}

\subsection{Common problems}\label{sec:CommProb}

The focus of this paper is on the following problems appearing in many applications.

The \textbf{subgraph isomorphism problem} is the following: is

$G_{1}$ isomorphic to any subgraph of $G_{2}$ by a given node

label?

The \textbf{induced subgraph isomorphism problem} asks the same about the

existence of an induced subgraph.

The \textbf{graph isomorphism problem} can be defined as induced

subgraph matching problem where the sizes of the two graphs are equal.

In addition, one may either want to find a \textbf{single} embedding or \textbf{enumerate} all of them.

\section{The VF2 Algorithm}\label{sec:VF2Alg}

This algorithm is the basis of both the VF2++ and the VF2~Plus.  VF2

is able to handle all the variations mentioned in Section~\ref{sec:CommProb}.  Although it can also handle directed graphs,

for the sake of simplicity, only the undirected case is

discussed.

\subsection{Common notations}

\indent Assume $G_{1}$ is searched in $G_{2}$.  The following

definitions and notations is used throughout this paper.

\begin{definition}

An injection $\mathfrak{m} : D \longrightarrow V_2$ is called (partial) \textbf{mapping}, where $D\subseteq V_1$.

\end{definition}

\begin{notation}

$\mathfrak{D}(f)$ and $\mathfrak{R}(f)$ denote the domain and the range of a function $f$, respectively.

\end{notation}

\begin{definition}

Mapping $\mathfrak{m}$ \textbf{covers} a node $u\in V_1\cup V_2$ if $u\in \mathfrak{D}(\mathfrak{m})\cup \mathfrak{R}(\mathfrak{m})$.

\end{definition}

\begin{definition}

A mapping $\mathfrak{m}$ is $\mathbf{whole\ mapping}$ if $\mathfrak{m}$ covers all the

nodes of $V_{1}$, i.e. $\mathfrak{D}(\mathfrak{m})=V_1$.

\end{definition}

\begin{definition}

Let \textbf{extend}$(\mathfrak{m},(u,v))$ denote the function $f : \mathfrak{D}(\mathfrak{m})\cup\{u\}\longrightarrow\mathfrak{R}(\mathfrak{m})\cup\{v\}$, for which $\forall w\in \mathfrak{D}(\mathfrak{m}) : f(w)=\mathfrak{m}(w)$ and $f(u)=v$ holds, where $u\in V_1\setminus\mathfrak{D}(\mathfrak{m})$ and $v\in V_2\setminus\mathfrak{R}(\mathfrak{m})$; otherwise $extend(\mathfrak{m},(u,v))$ is undefined.

\end{definition}

\begin{notation}

Throughout the paper, $\mathbf{PT}$ denotes a generic problem type

which can be substituted by any of the $\mathbf{SUB}$, $\mathbf{IND}$

and $\mathbf{ISO}$ problems, which stand for the problems mentioned in Section~\ref{sec:CommProb} respectively.

\end{notation}

\begin{definition}

Let $\mathfrak{m}$ be a mapping. The \textbf{consistency function for } $\mathbf{PT}$ is a logical function, for which $\mathbf{Cons_{PT}}(\mathfrak{m})$ is true if and only if $\mathfrak{m}$ satisfies the requirements of $\mathbf{PT}$ considering the subgraphs of $G_{1}$ and $G_{2}$ induced by $\mathfrak{D}(\mathfrak{m})$ and $\mathfrak{R}(\mathfrak{m})$, respectively.

%$\mathbf{Cons_{PT}}(\mathfrak{m})$ is true if there exists a whole mapping $w$ satisfying the requirements of $PT$, for which $\mathfrak{m}$ is exactly $w$ restricted to $\mathfrak{D}(\mathfrak{m})$.

\end{definition}

\begin{definition} 

Let $\mathfrak{m}$ be a mapping. A logical function $\mathbf{Cut_{PT}}$ is a

\textbf{cutting function for } $\mathbf{PT}$ if the following

holds. $\mathbf{Cut_{PT}(\mathfrak{m})}$ is false if there exists a sequence of extend operations, which results in a whole mapping satisfying the requirements of $PT$.

\end{definition}

\begin{definition}

A mapping $\mathfrak{m}$ is said to be \textbf{consistent mapping by} $\mathbf{PT}$ if

  $Cons_{PT}(\mathfrak{m})$ is true.

\end{definition}

$Cons_{PT}$ and $Cut_{PT}$ will often be used in the following form.

\begin{notation}

Let $\mathbf{Cons_{PT}(p, \mathfrak{m})}:=Cons_{PT}(extend(\mathfrak{m},p))$, and

$\mathbf{Cut_{PT}(p, \mathfrak{m})}:=Cut_{PT}(extend(\mathfrak{m},p))$, where

$p\in{V_{1}\backslash\mathfrak{D}(\mathfrak{m}) \!\times\!V_{2}\backslash\mathfrak{R}(\mathfrak{m})}$.

\end{notation}

$Cons_{PT}$ will be used to check the consistency of the already

covered nodes, while $Cut_{PT}$ is for looking ahead to recognize if

no whole consistent mapping can contain the current mapping.

\subsection{Overview of the algorithm}

VF2 begins with an empty mapping and gradually extends it with respect to the consistency and cutting functions until a whole mapping is reached.

Algorithm~\ref{alg:VF2Pseu} is a high-level description of

the VF2 algorithm. Each state of the matching process can

be associated with a mapping $\mathfrak{m}$. The initial state

is associated with a mapping $\mathfrak{m}$, for which

$\mathfrak{D}(\mathfrak{m})=\emptyset$, i.e. it starts with an empty mapping.

\begin{algorithm}

\algtext*{EndIf}%ne nyomtasson end if-et

\algtext*{EndFor}%ne

\algtext*{EndProcedure}%ne nyomtasson ..

\caption{\hspace{0.5cm}$A\ high\ level\ description\ of\ VF2$}\label{alg:VF2Pseu}

\begin{algorithmic}[1]

\Procedure{VF2}{Mapping $\mathfrak{m}$, ProblemType $PT$}

  \If{$\mathfrak{m}$ covers

    $V_{1}$} \State Output($\mathfrak{m}$)

  \Else

  \State Compute the set $P_\mathfrak{m}$ of the candidate pairs for extending $\mathfrak{m}$ \ForAll{$p\in{P_\mathfrak{m}}$} \If{Cons$_{PT}$($p,\mathfrak{m}$) $\wedge$

    $\neg$Cut$_{PT}$($p,\mathfrak{m}$)}

    \State \textbf{call}

  VF2($extend(\mathfrak{m},p)$, $PT$) \EndIf \EndFor \EndIf \EndProcedure

\end{algorithmic}

\end{algorithm}

For the current mapping $\mathfrak{m}$, the algorithm computes $P_\mathfrak{m}$, the set of

candidate node pairs for extending the current mapping $\mathfrak{m}$.

For each pair $p$ in $P_\mathfrak{m}$, $Cons_{PT}(p,\mathfrak{m})$ and

$Cut_{PT}(p,\mathfrak{m})$ are evaluated. If the former is true and

the latter is false, the whole process is recursively applied to

$extend(\mathfrak{m},p)$. Otherwise, $extend(\mathfrak{m},p)$ is not consistent by $PT$, or it

can be proved that $\mathfrak{m}$ can not be extended to a whole mapping.

The correctness of the procedure follows from the claim below.

\begin{claim}\label{claim:consMapps}

Through consistent mappings, only consistent whole mappings can be

reached, and all the consistent whole mappings are reachable through

consistent mappings.

\end{claim}

Note that a mapping may be reached in exponentially many different ways, since the

order of extensions does not influence the nascent mapping.

However, one may make the following observations.

%\begin{claim}

%\label{claim:claimTotOrd}

%Let $\prec$ be an arbitrary total ordering relation on $V_{1}$.  If

%the algorithm ignores each $p=(u,v) \in P_\mathfrak{m}$, for which

%\begin{center}

%$\exists (\tilde{u},\tilde{v})\in P_\mathfrak{m}: \tilde{u} \prec u$,

%\end{center}

%then no mapping can be reached more than once, and each whole mapping %remains reachable.

%\end{claim}

\begin{definition}

A total order $(u_{\sigma(1)},u_{\sigma(2)},..,u_{\sigma(|V_{1}|)})$ of

$V_{1}$ is \textbf{matching order} if VF2 can cover $u_{\sigma(d)}$ on the $d$-th level for all $d\in\{1,..,|V_{1}|\}$.

\end{definition}

\begin{claim}

\label{claim:claimTotOrd}

If VF2 is prescribed to cover the nodes of $G_1$ according to a matching order, then no mapping can be reached more than once and each whole mapping remains reachable.

\end{claim}

Note that the cornerstone of the improvements to VF2 is to choose a proper

matching order.

\subsection{The candidate set}

\label{candidateComputingVF2}

Let $P_\mathfrak{m}$ be the set of the candidate pairs for inclusion in $\mathfrak{m}$.

\begin{notation}

Let $\mathbf{T_{1}(\mathfrak{m})}:=\{u \in V_{1}\backslash\mathfrak{D}(\mathfrak{m}) : \exists \tilde{u}\in{\mathfrak{D}(\mathfrak{m}): (u,\tilde{u})\in E_{1}}\}$, and

 $\mathbf{T_{2}(\mathfrak{m})} := \{v \in V_{2}\backslash\mathfrak{R}(\mathfrak{m}) : \exists\tilde{v}\in{\mathfrak{R}(\mathfrak{m}):(v,\tilde{v})\in E_{2}}\}$.

\end{notation}

The set $P_\mathfrak{m}$ contains the pairs of uncovered neighbours of covered

nodes, and if there is not such a node pair, all the pairs containing

two uncovered nodes are added. Formally, let

\[

 P_\mathfrak{m}\!=\!

  \begin{cases} 

   T_{1}(\mathfrak{m})\times T_{2}(\mathfrak{m})&\hspace{-0.15cm}\text{if }

   T_{1}(\mathfrak{m})\!\neq\!\emptyset\ \text{and }T_{2}(\mathfrak{m})\!\neq

   \emptyset,\\ (V_{1}\!\setminus\!\mathfrak{D}(\mathfrak{m}))\!\times\!(V_{2}\!\setminus\!\mathfrak{R}(\mathfrak{m}))

   &\hspace{-0.15cm}\text{otherwise}.

  \end{cases}

\]

\subsection{Consistency}

Let $p=(u,v)\in V_{1}\times V_{2}$, and suppose $\mathfrak{m}$ is a consistent mapping by

$PT$. $Cons_{PT}(p,\mathfrak{m})$ checks whether

adding pair $p$ into $\mathfrak{m}$ leads to a consistent mapping by $PT$.

For example, the consistency function of the induced subgraph isomorphism problem is the following.

\begin{notation}

Let $\mathbf{\Gamma_{1} (u)}:=\{\tilde{u}\in V_{1} :

(u,\tilde{u})\in E_{1}\}$, and $\mathbf{\Gamma_{2}

  (v)}:=\{\tilde{v}\in V_{2} : (v,\tilde{v})\in E_{2}\}$, where $u\in V_{1}$ and $v\in V_{2}$. That is, $\mathbf{\Gamma_{i} (w)}$ denotes the set of neighbours of node $w$ in $G_i$ $(i=1,2)$.

\end{notation}

\begin{claim}

$extend(\mathfrak{m},(u,v))$ is a consistent mapping by $IND$ if and only if $\mathfrak{m}$ is consistent and $(\forall \tilde{u}\in \mathfrak{D}(\mathfrak{m}): (u,\tilde{u})\in E_{1}

\Leftrightarrow (v,\mathfrak{m}(\tilde{u}))\in E_{2})$. 

\end{claim}

The following formulation gives an efficient way of calculating $Cons_{IND}$.

\begin{claim}

$Cons_{IND}((u,v),\mathfrak{m}):=Cons_{IND}(\mathfrak{m})\wedge\mathcal{L}(u)\!\!=\!\!\mathcal{L}(v)\wedge(\forall \tilde{v}\in \Gamma_{2}(v)\cap\mathfrak{R}(\mathfrak{m}):(u,\mathfrak{m}^{-1}(\tilde{v}))\in E_{1})\wedge

  (\forall \tilde{u}\in \Gamma_{1}(u)

  \cap \mathfrak{D}(\mathfrak{m}):(v,\mathfrak{m}(\tilde{u}))\in E_{2})$ is the consistency function for $IND$.

\end{claim}

\subsection{Cutting rules}

$Cut_{PT}(p,\mathfrak{m})$ is defined by a collection of efficiently

verifiable conditions. The requirement is that $Cut_{PT}(p,\mathfrak{m})$ can

be true only if it is impossible to extend $extend(\mathfrak{m},p)$ to a

whole mapping.

As an example, a cutting function of induced subgraph isomorphism problem is presented.

\begin{notation}

Let $\mathbf{\tilde{T}_{1}}(\mathfrak{m}):=(V_{1}\backslash

\mathfrak{D}(\mathfrak{m}))\backslash T_{1}(\mathfrak{m})$, and

\\ $\mathbf{\tilde{T}_{2}}(\mathfrak{m}):=(V_{2}\backslash

\mathfrak{R}(\mathfrak{m}))\backslash T_{2}(\mathfrak{m})$.

\end{notation}

\begin{claim}

$Cut_{IND}((u,v),\mathfrak{m}):= |\Gamma_{2} (v)\ \cap\ T_{2}(\mathfrak{m})| <

  |\Gamma_{1} (u)\ \cap\ T_{1}(\mathfrak{m})| \vee |\Gamma_{2}(v)\cap

  \tilde{T}_{2}(\mathfrak{m})| < |\Gamma_{1}(u)\cap

  \tilde{T}_{1}(\mathfrak{m})|$ is a cutting function for $IND$.

\end{claim}

\section{The VF2++ Algorithm}\label{sec:VF2ppAlg}

Although any matching order makes the search space of VF2 a

tree, its choice turns out to dramatically influence the number of

visited states. The goal is to determine an efficient one as quickly

as possible.

The main reason for the superiority of VF2++ over VF2 is twofold. Firstly,

taking into account the structure and the node labeling of the graph,

VF2++ determines a matching order in which most of the unfruitful

branches of the search space can be pruned immediately. Secondly,

introducing more efficient --- nevertheless still easier to compute

--- cutting rules reduces the chance of going astray even further.

In addition to the usual subgraph isomorphism problem, specialized versions

for induced subgraph and graph isomorphism problems have also been

designed.

Note that a weaker version of the cutting rules of VF2++ and an efficient

candidate set calculation method were described in~\cite{VF2Plus}.

It should be noted that all the methods described in this section are

extendable to handle directed graphs and edge labels as well.

The basic ideas and the detailed description of VF2++ are provided in

the following.\newline

The goal is to find a matching order in which the algorithm is able to

recognize inconsistency or prune the infeasible branches on the

highest levels and goes deep only if it is needed.

\begin{notation}

Let $\mathbf{Conn_{H}(u)}:=|\Gamma_{1}(u)\cap H|$, that is the

number of neighbours of u which are in H, where $u\in V_{1} $ and

$H\subseteq V_{1}$.

\end{notation}

The principal question is the following. Suppose a mapping $\mathfrak{m}$ is

given. For which node of $T_{1}(\mathfrak{m})$ is the hardest to find a

consistent pair in $G_{2}$? The more covered neighbors a node in

$T_{1}(\mathfrak{m})$ has --- i.e. the largest $Conn_{\mathfrak{D}(\mathfrak{m})}$ it has

---, the more rare-to-satisfy consistency constraints for its pair

are given.

Most of the graphs of biological and chemical structures are sparse, thus several nodes in

$T_{1}(\mathfrak{m})$ may have the same $Conn_{\mathfrak{D}(\mathfrak{m})}$, which makes

reasonable to define a secondary and a tertiary order between them.

The observation above proves itself to be as determining, that the

secondary ordering prefers nodes with the most uncovered neighbors

among which have the same $Conn_{\mathfrak{D}(\mathfrak{m})}$ to increase

$Conn_{\mathfrak{D}(\mathfrak{m})}$ of uncovered nodes as much, as possible.  The tertiary ordering prefers nodes having the rarest uncovered labels in $G_2$.

Note that the secondary ordering is the same as ordering by degrees,

which is a static data in front of the above used.

These rules can easily result in a matching order which contains the

nodes of a long path successively, whose nodes may have low $Conn$ and

is easily to match into $G_{2}$. To try to avoid that, a

Breadth-First-Search order is used, and on each of its levels, the

ordering procedure described above is applied.

\newline

In the following, examples are shown, demonstrating that VF2 may be

slow are, even though a matching can be found easily by using a proper

matching order.

\begin{example}

Suppose $G_{1}$ can be mapped into $G_{2}$ in many ways

without node labels. Let $u\in V_{1}$ and $v\in V_{2}$.

\newline

$\mathcal{L}(u):=black$

\newline

$\mathcal{L}(v):=black$

\newline

$\mathcal{L}(\tilde{u}):=red \ \forall \tilde{u}\in V_{1}\backslash

\{u\}$

\newline

$\mathcal{L}(\tilde{v}):=red \ \forall \tilde{v}\in V_{2}\backslash

\{v\}$

Now, any mapping by $\mathcal{L}$ must contain $(u,v)$, since

$u$ is black and no node in $V_{2}$ has a black label except

$v$. If unfortunately $u$ were the last node which will get covered,

VF2 would check only in the last steps, whether $u$ can be matched to

$v$.

However, had $u$ been the first matched node, $u$ would have been

matched immediately to $v$, so all the mappings would have been

precluded in which node labels can not correspond.

\end{example}

\begin{example}

Suppose there is no node label given, and $G_{1}$ is a small graph that can not be mapped into $G_{2}$ and $u\in V_{1}$.

\newline

Let $G'_{1}:=(V_{1}\cup

\{u'_{1},u'_{2},..,u'_{k}\},E_{1}\cup

\{(u,u'_{1}),(u'_{1},u'_{2}),..,(u'_{k-1},u'_{k})\})$, that is,

$G'_{1}$ is $G_{1}\cup \{ a\ k$ long path, which is disjoint

from $G_{1}$ and one of its starting points is connected to $u\in

V_{1}\}$.

If, unfortunately, the nodes of the path were the first $k$ nodes in the

matching order, the algorithm would iterate through all the possible $k$

long paths in $G_{2}$, and it would recognize that no path can be

extended to $G'_{1}$.

\newline

However, had it started by the matching of $G_{1}$, it would not

have matched any nodes of the path.

\end{example}

These examples may look artificial, but the same problems also appear

in real-world instances, even though in a less obvious way.

%\subsection{Preparations}

%\begin{claim}

%\label{claim:claimCoverFromLeft}

%The total ordering relation uniquely determines a node order, in which

%the nodes of $V_{1}$ will be covered by VF2. From the point of

%view of the matching procedure, this means, that always the same node

%of $G_{1}$ will be covered on the $d$-th level.

%\end{claim}

%\begin{definition}

%An order $(u_{\sigma(1)},u_{\sigma(2)},..,u_{\sigma(|V_{1}|)})$ of

%$V_{1}$ is \textbf{matching order} if there exists $\prec$ total

%ordering relation, s.t. the VF2 with $\prec$ on the d-th level finds

%pair for $u_{\sigma(d)}$ for all $d\in\{1,..,|V_{1}|\}$.

%\end{definition}

%\begin{claim}\label{claim:MOclaim}

%A total ordering is matching order iff the nodes of every component

%form an interval in the node sequence, and every node connects to a

%previous node in its component except the first node of each component.

%\end{claim}

%In summary, a total ordering always uniquely determines a matching

%order, and every matching order can be determined by a total ordering,

%however, more than one different total orderings may determine the

%same matching order.

\subsection{Matching order}

\begin{notation}

Let \textbf{F$_\mathcal{M}$(l)}$:=|\{v\in V_{2} :

l=\mathcal{L}(v)\}|-|\{u\in \mathcal{M} : l=\mathcal{L}(u)\}|$,

where $l$ is a label and $\mathcal{M}\subseteq V_{1}$.

\end{notation}

\begin{definition}Let $\mathbf{arg\ max}_{f}(S) :=\{u\in S : f(u)=max_{v\in S}\{f(v)\}\}$ and $\mathbf{arg\ min}_{f}(S) := arg\ max_{(-f)}(S)$, where $S$ is a finite set and $f:S\longrightarrow \mathbb{R}$.

\end{definition}

\begin{notation}

Let $deg(v)$ denote the degree of node $v$.

\end{notation}

\begin{algorithm}[H]

\algtext*{EndIf}

\algtext*{EndProcedure}

\algtext*{EndWhile}

\algtext*{EndFor}

\caption{\hspace{0.5cm}$The\ method\ of\ VF2++\ for\ determining\ the\ node\ order$}\label{alg:VF2PPPseu}

\begin{algorithmic}[1]

\Procedure{VF2++order}{} \State $\mathcal{M}$ := $\emptyset$

\Comment{matching order} \While{$V_{1}\backslash \mathcal{M}

  \neq\emptyset$} \State $r\in$ arg max$_{deg}$ (arg

min$_{F_\mathcal{M}\circ \mathcal{L}}(V_{1}\backslash

\mathcal{M})$)\label{alg:findMin} \State Compute $T$, a BFS tree with

root node $r$.  \For{$d=0,1,...,depth(T)$} \State $V_d$:=nodes of the

$d$-th level \State Process $V_d$ \Comment{See Algorithm~\ref{alg:VF2PPProcess1}} \EndFor

\EndWhile \EndProcedure

\end{algorithmic}

\end{algorithm}

\begin{algorithm}

\algtext*{EndIf}

\algtext*{EndProcedure}%ne nyomtasson ..

\algtext*{EndWhile}

\caption{\hspace{.5cm}$The\ method\ for\ processing\ a\ level\ of\ the\ BFS\ tree$}\label{alg:VF2PPProcess1}

\begin{algorithmic}[1]

\Procedure{VF2++ProcessLevel}{$V_{d}$} \While{$V_d\neq\emptyset$}

\State $m\in$ arg min$_{F_\mathcal{M}\circ\ \mathcal{L}}($ arg max$_{deg}($arg

max$_{Conn_{\mathcal{M}}}(V_{d})))$ \State $V_d:=V_d\backslash m$

\State Append node $m$ to the end of $\mathcal{M}$ \State Refresh

$F_\mathcal{M}$ \EndWhile \EndProcedure

\end{algorithmic}

\end{algorithm}

Algorithm~\ref{alg:VF2PPPseu} is a high-level description of the

matching order procedure of VF2++. It computes a BFS tree for each

component in ascending order of their rarest node labels and largest $deg$,

whose root vertex is the minimal node of its component. Algorithm~\ref{alg:VF2PPProcess1} is a method to process a level of the BFS tree, which appends the nodes of the current level in descending

lexicographic order by $(Conn_{\mathcal{M}},deg,-F_\mathcal{M})$ separately

to $\mathcal{M}$, and refreshes $F_\mathcal{M}$ immediately.

\begin{claim}

Algorithm~\ref{alg:VF2PPPseu} provides a matching order.

\end{claim}

\subsection{Cutting rules}

\label{VF2PPCuttingRules}

This section presents the cutting rules of VF2++, which are improved by using extra information coming from the node labels.

\begin{notation}

Let $\mathbf{\Gamma_{1}^{l}(u)}:=\{\tilde{u} : \mathcal{L}(\tilde{u})=l

\wedge \tilde{u}\in \Gamma_{1} (u)\}$ and

$\mathbf{\Gamma_{2}^{l}(v)}:=\{\tilde{v} : \mathcal{L}(\tilde{v})=l \wedge

\tilde{v}\in \Gamma_{2} (v)\}$, where $u\in V_{1}$, $v\in

V_{2}$ and $l$ is a label.

\end{notation}

\begin{claim}[Cutting function for ISO]

\[LabCut_{ISO}((u,v),\mathfrak{m}):=\bigvee_{l\ is\ label}|\Gamma_{2}^{l} (v) \cap T_{2}(\mathfrak{m})|\!\neq\!|\Gamma_{1}^{l}(u)\cap T_{1}(\mathfrak{m})|\  \vee\]\[\bigvee_{l\ is\ label} \newline |\Gamma_{2}^{l}(v)\cap \tilde{T}_{2}(\mathfrak{m})| \neq |\Gamma_{1}^{l}(u)\cap \tilde{T}_{1}(\mathfrak{m})|\] is a cutting function for ISO.

\end{claim}

\begin{claim}[Cutting function for IND]

\[LabCut_{IND}((u,v),\mathfrak{m}):=\bigvee_{l\ is\ label}|\Gamma_{2}^{l} (v) \cap T_{2}(\mathfrak{m})|\!<\!|\Gamma_{1}^{l}(u)\cap T_{1}(\mathfrak{m})|\ \vee\]\[\bigvee_{l\ is\ label} \newline |\Gamma_{2}^{l}(v)\cap \tilde{T}_{2}(\mathfrak{m})| < |\Gamma_{1}^{l}(u)\cap \tilde{T}_{1}(\mathfrak{m})|\] is a cutting function for IND.

\end{claim}

\begin{claim}[Cutting function for SUB]

\[LabCut_{SU\!B}((u,v),\mathfrak{m}):=\bigvee_{l\ is\ label}|\Gamma_{2}^{l} (v) \cap T_{2}(\mathfrak{m})|\!<\!|\Gamma_{1}^{l}(u)\cap T_{1}(\mathfrak{m})|\] is a cutting function for SUB.

\end{claim}

\section{Implementation details}\label{sec:VF2ppImpl}

This section provides a detailed summary of an efficient

implementation of VF2++.

\begin{notation}

Let $\Delta_1$ and $\Delta_2$ denote the largest degree in $G_1$ and $G_2$, respectively, and let $\Delta=\max\{\Delta_1,\Delta_2\}$.

\end{notation}

\subsection{Storing a mapping}

After fixing an arbitrary node order ($u_0, u_1, ..,

u_{|V_{1}|-1}$) of $G_{1}$, an array $M$ can be used to store

the current mapping in the following way.

\[

 M[i] =

  \begin{cases} 

   v & if\ (u_i,v)\ is\ in\ the\ mapping\\ INV\!ALI\!D &

   if\ no\ node\ has\ been\ mapped\ to\ u_i,

  \end{cases}

\]

where $i\in\{0,1, ..,|V_{1}|-1\}$, $v\in V_{2}$ and $INV\!ALI\!D$

means "no node".

\subsection{Avoiding the recurrence}

The recursion of Algorithm~\ref{alg:VF2Pseu} can be realized

as a \textit{while loop}, which has a loop counter $depth$ denoting the

current depth of the recursion. Fixing a matching order, let $M$

denote the array storing the current mapping. Observe that

$M$ is $INV\!ALI\!D$ from index $depth$+1 and not $INV\!ALI\!D$ before

$depth$. $M[depth]$ changes

while the state is being processed, but the property is held before

both stepping back to a predecessor state and exploring a successor

state.

The necessary part of the candidate set is easy to maintain or

compute by following

the steps described in Section~\ref{candidateComputingVF2}. A much faster method

has been designed for biological and sparse graphs, see the next

section for details.

\subsection{Calculating the candidates for a node}

The task is not to maintain the candidate set, but to generate the

candidate nodes in $G_{2}$ for a given node $u\in V_{1}$.  In

case of any of the three problem types and a mapping $\mathfrak{m}$, if a node $v\in

V_{2}$ is a potential pair of $u\in V_{1}$, then $\forall

u'\in \mathfrak{D}(\mathfrak{m}) : (u,u')\in

E_{1}\Rightarrow (v,\mathfrak{m}(u'))\in

E_{2}$. That is, each covered neighbor of $u$ has to be mapped to

a covered neighbor of $v$, i.e. selecting arbitrarily a covered neighbor $u'$ of $u$, all of the admissible candidates for $u$ are among the neighbors of $\mathfrak{m}(u')$.

Having said that, an algorithm running in $\Theta(\Delta_2)$ time is

describable if there exists a covered node in the component containing

$u$, and a linear one otherwise.

\subsection{Determining the node order}

This section describes how the node order preprocessing method of

VF2++ can efficiently be implemented.

For using lookup tables, the node labels are associated with the

numbers $\{0,1,..,|K|-1\}$, where $K$ is the set of the labels. It

enables $F_\mathcal{M}$ to be stored in an array. At first, the node order

$\mathcal{M}=\emptyset$, so $F_\mathcal{M}[i]$ is the number of nodes

in $V_{2}$ having label $i$, which is easy to compute in

$\Theta(|V_{2}|)$ steps.

Representing $\mathcal{M}\subseteq V_{1}$ as an array of

size $|V_{1}|$, both the computation of the BFS tree, and processing its levels by Algorithm~\ref{alg:VF2PPProcess1} can be done in-place by swapping nodes.

\subsection{Cutting rules}

Section~\ref{VF2PPCuttingRules} described the cutting rules

using the sets $T_{1}$, $T_{2}$, $\tilde T_{1}$

and $\tilde T_{2}$, which are dependent on the current mapping. The aim is to check the labeled cutting

rules of VF2++ in $\Theta(\Delta)$ time.

Firstly, suppose that these four sets are given a way, that

checking whether a node is in a certain set takes constant time,

e.g. they are given by their 0-1 characteristic vectors. Let $L$ be an

initially zero integer lookup table of size $|K|$. After incrementing

$L[\mathcal{L}(u')]$ for all $u'\in \Gamma_{1}(u) \cap T_{1}(\mathfrak{m})$ and

decrementing $L[\mathcal{L}(v')]$ for all $v'\in\Gamma_{2} (v) \cap

T_{2}(\mathfrak{m})$, the first part of the cutting rules can be checked in

$\Theta(\Delta)$ time by considering the proper signs of $L$. Setting $L$

to zero takes $\Theta(\Delta)$ time again, which makes it possible to use

the same table through the whole algorithm. The second part of the

cutting rules can be verified using the same method with $\tilde

T_{1}$ and $\tilde T_{2}$ instead of $T_{1}$ and

$T_{2}$. Thus, the overall time complexity is $\Theta(\Delta)$.

To maintain the sets $T_{1}$, $T_{2}$, $\tilde T_{1}$

and $\tilde T_{2}$, two other integer lookup tables storing the number of covered neighbours of the nodes of the two graphs can be used. This representation allows constant-time membership checking, furthermore it is maintainable in $\Theta(\Delta)$ time whenever a node pair is added or subtracted by incrementing

or decrementing the proper indices. A further improvement is that the

values of $L[\mathcal{L}(u')]$ in case of checking $u$ are dependent only on

$u$, i.e. on the current depth of the recursion, so for each $u\in V_{1}$, an array of pairs \textit{(label, number of such labels)} can store $L$. Note that these arrays are at most of size

$\Delta_1$ if pairs with non-appearing node labels are discarded.

Using similar techniques, the consistency function can be evaluated in

$\Theta(\Delta)$ steps, as well.

\section{Experimental results}\label{sec:ExpRes}

This section compares the performance of VF2++ and VF2~Plus. According to

our experience, both algorithms run faster than VF2 with orders of

magnitude, thus its inclusion was not reasonable.

The algorithms were implemented in C++ using the open-source LEMON

graph and network optimization

library~\cite{LEMON}\cite{o11:_lemon_open_sourc_c_graph_templ_librar}. The

tests were carried out on a Linux-based system with an Intel i7 X980

3.33 GHz CPU and 6 GB of RAM.

\subsection{Biological graphs}

The tests have been executed on a recent biological dataset created

for the International Contest on Pattern Search in Biological

Databases~\cite{Content}, which has been constructed of molecule,

protein and contact map graphs extracted from the Protein Data

Bank~\cite{ProteinDataBank}.

The molecule dataset contains small graphs with less than 100 nodes

and an average degree of less than 3. The protein dataset contains

graphs having 500-10 000 nodes and an average degree of 4, while the

contact map dataset contains graphs with 150-800 nodes and an average

degree of 20.  \\

In the following, both the induced subgraph and the graph

isomorphism problems will be examined.

This dataset provides graph pairs, between which all the induced subgraph isomorphisms have to be found. For the running times, please see Figure~\ref{fig:bioIND}.

In an other experiment, the nodes of each graph in the database had been

shuffled, and an isomorphism between the shuffled and the original

graph was searched. The running times are shown on Figure~\ref{fig:bioISO}.

\begin{figure}[H]

\vspace*{-2cm}

\hspace*{-1.5cm}

\begin{subfigure}[b]{0.55\textwidth}

\begin{figure}[H]

\begin{tikzpicture}[trim axis left, trim axis right]

\begin{axis}[title=Molecules IND,xlabel={$|V_2|$},ylabel={time (ms)},legend entries={VF2 Plus,VF2++},grid

=major,mark size=1.2pt, legend style={at={(0,1)},anchor=north

  west},scaled x ticks = false,x tick label style={/pgf/number

  format/1000 sep = \kern 0.08em},y tick label style={/pgf/number

  format/1000 sep = \kern 0.08em}]

%\addplot+[only marks] table {proteinsOrig.txt};

\addplot table {Orig/Molecules.32.txt}; \addplot[mark=triangle*,mark

  size=1.8pt,color=red] table {VF2PPLabel/Molecules.32.txt};

\end{axis}

\end{tikzpicture}

\caption{In the case of molecules, the algorithms have

  similar behaviour, but VF2++ is almost two times faster even on such

  small graphs.} \label{fig:INDMolecule}

\end{figure}

\end{subfigure}

\hspace*{1.5cm}

\begin{subfigure}[b]{0.55\textwidth}

\begin{figure}[H]

\begin{tikzpicture}[trim axis left, trim axis right]

\begin{axis}[title=Contact maps IND,xlabel={$|V_2|$},ylabel={time (ms)},legend entries={VF2 Plus,VF2++},grid

=major,mark size=1.2pt, legend style={at={(0,1)},anchor=north

  west},scaled x ticks = false,x tick label style={/pgf/number

  format/1000 sep = \kern 0.08em},y tick label style={/pgf/number

  format/1000 sep = \kern 0.08em}]

%\addplot+[only marks] table {proteinsOrig.txt};

\addplot table {Orig/ContactMaps.128.txt};

\addplot[mark=triangle*,mark size=1.8pt,color=red] table

        {VF2PPLabel/ContactMaps.128.txt};

\end{axis}

\end{tikzpicture}

\caption{On contact maps, VF2++ runs almost in constant time, while VF2

  Plus has a near-linear behavior.} \label{fig:INDContact}

\end{figure}

\end{subfigure}

\begin{center}

\vspace*{-0.5cm}

\begin{subfigure}[b]{0.55\textwidth}

\begin{figure}[H]

\begin{tikzpicture}[trim axis left, trim axis right]

  \begin{axis}[title=Proteins IND,xlabel={$|V_2|$},ylabel={time (ms)},legend entries={VF2 Plus,VF2++},grid

  =major,mark size=1.2pt, legend style={at={(0,1)},anchor=north

    west},scaled x ticks = false,x tick label style={/pgf/number

  format/1000 sep = \kern 0.08em},y tick label style={/pgf/number

  format/1000 sep = \kern 0.08em}] %\addplot+[only marks] table

    {proteinsOrig.txt}; \addplot[mark=*,mark size=1.2pt,color=blue]

    table {Orig/Proteins.256.txt}; \addplot[mark=triangle*,mark

      size=1.8pt,color=red] table {VF2PPLabel/Proteins.256.txt};

  \end{axis}

  \end{tikzpicture}

\caption{Both of the algorithms have linear behaviour on protein

  graphs. VF2++ is more than 10 times faster than VF2

  Plus.} \label{fig:INDProt}

\end{figure}

\end{subfigure}

\end{center}

\vspace*{-0.5cm}

\caption{\normalsize{Induced subgraph isomorphism problem on biological graphs}}\label{fig:bioIND}

\end{figure}

\begin{figure}[H]

\vspace*{-2cm}

\hspace*{-1.5cm}

\begin{subfigure}[b]{0.55\textwidth}

\begin{figure}[H]

\begin{tikzpicture}[trim axis left, trim axis right]

\begin{axis}[title=Molecules ISO,xlabel={$|V_2|$},ylabel={time (ms)},legend entries={VF2 Plus,VF2++},grid

=major,mark size=1.2pt, legend style={at={(0,1)},anchor=north

  west},scaled x ticks = false,x tick label style={/pgf/number

  format/1000 sep = \kern 0.08em},y tick label style={/pgf/number

  format/1000 sep = \kern 0.08em}]

%\addplot+[only marks] table {proteinsOrig.txt};

\addplot table {Orig/moleculesIso.txt}; \addplot[mark=triangle*,mark

  size=1.8pt,color=red] table {VF2PPLabel/moleculesIso.txt};

\end{axis}

\end{tikzpicture}

\caption{The results are close to each other on contact maps, but VF2++ seems to be slightly faster as the number of nodes increases.

  }\label{fig:ISOMolecule}

\end{figure}

\end{subfigure}

\hspace*{1.5cm}

\begin{subfigure}[b]{0.55\textwidth}

\begin{figure}[H]

\begin{tikzpicture}[trim axis left, trim axis right]

\begin{axis}[title=Contact maps ISO,xlabel={$|V_2|$},ylabel={time (ms)},legend entries={VF2 Plus,VF2++},grid

=major,mark size=1.2pt, legend style={at={(0,1)},anchor=north

  west},scaled x ticks = false,x tick label style={/pgf/number

  format/1000 sep = \kern 0.08em},y tick label style={/pgf/number

  format/1000 sep = \kern 0.08em}]

%\addplot+[only marks] table {proteinsOrig.txt};

\addplot table {Orig/contactMapsIso.txt}; \addplot[mark=triangle*,mark

  size=1.8pt,color=red] table {VF2PPLabel/contactMapsIso.txt};

\end{axis}

\end{tikzpicture}