/**------------------------------------------------------------------------
 * @file
 *
 * Copyright 2001 by the University of Toronto (UT).
 *
 *                       All Rights Reserved
 * 
 * Permission to copy and modify this software and its documentation 
 * only for internal use in your organization is hereby granted, 
 * provided that this notice is retained thereon and on all copies.  
 * UT makes no representations as to the suitability and operability 
 * of this software for any purpose.  It is provided "as is" without 
 * express or implied warranty.
 *
 * UT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING
 * ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, IN NO EVENT
 * SHALL UT BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
 * DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
 * OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
 * TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
 * PERFORMANCE OF THIS SOFTWARE.
 *
 * No other rights, including for example, the right to redistribute 
 * this software and its documentation or the right to prepare 
 * derivative works, are granted unless specifically provided in a 
 * separate license agreement.
 *
 * Author: Diego Macrini based on Ali Shokoufandeh's original code
 * 
 * $Log: DAG.cpp,v $
 * Revision 1.6  2002/06/17 23:09:11  dmac
 * Matching algorithm updated.
 *
 * Revision 1.5  2002/06/14 23:13:03  dmac
 * Matching algorithm updated.
 *
 * Revision 1.4  2002/06/13 21:46:37  dmac
 * Matching algorithm updated.
 *
 * Revision 1.3  2002/06/10 18:37:33  dmac
 * New class that handles the node distance function added.
 *
 * Revision 1.2  2002/05/16 22:48:56  dmac
 * Several bugs fixed and a new DAG-derived class added
 *
 * Revision 1.1.1.1  2002/03/10 07:38:07  dmac
 * Imported sources
 *
 *
 *-----------------------------------------------------------------------*/

#include "DAG.h"

//#include <LEDA/graph/templates/mwb_matching.h>
#include <LEDA/graph/mwb_matching.h>

// Fix an error in the next LEDA header
//#ifdef WIN32
using namespace std;
//#endif

#include <LEDA/graph/scale_weights.h>
//#include <LEDA/graphwin.h>

#include <float.h>
#include "HelperFunctions.h"
#include "Exceptions.h"
#include "BasicUtils.h"

#ifdef WIN32
#define isnan  _isnan
#define finite _finite
#endif

#define ROOT_NODE_LEVEL			0
#define ROOT_NODE_DFS_INDEX		0

using namespace dml;
using namespace leda;

DAG::MatchParams DAG::s_matchParams;
bool DAG::s_bDbgMode = false;

LogFile logFile;

/*!
	@brief Visits all the node in DFS order and computes they derived values from the
	info in the graph.

	@return the last DFS index + 1. i.e. The corresponding next DFS index.
*/

int DAG::ComputeNodesInfo(leda::node v, int nLevel, int nDFSIndex)
{
	leda::node u;
	DAGNodePtr& ptrNode = GetNode(v);

	// Default level is -1. Loopy nodes should have their greatest height
	if (ptrNode->GetLevel() < nLevel)
		ptrNode->SetLevel(nLevel);		
	
	if (ptrNode->GetDFSIndex() == -1)
		ptrNode->SetDFSIndex(nDFSIndex++);

	nLevel++;
	
	forall_adj_nodes(u, v)
		nDFSIndex = ComputeNodesInfo(u, nLevel, nDFSIndex);

	return nDFSIndex;
}

void DAG::ResetNodesInfo()
{
	leda::node v;
	
	forall_nodes(v, *this)
	{
		GetNode(v)->SetLevel(-1);
		GetNode(v)->SetDFSIndex(-1);
	}
}

void DAG::ComputeTransiveClosure()
{
	leda::node u, v;
	int i, n = GetNodeCount();

	transClosMat.ReSize(n, n, true);

#if LEDA_VERSION == 43
	leda::graph g = TRANSITIVE_CLOSURE(*this);
	
	// Translate the transitive closure graph into a matrix form
	forall_nodes(v, g)
	{
		i = GetNodeDFSIndex(v);

		forall_adj_nodes(u, v)
			transClosMat[i][GetNodeDFSIndex(u)] = 1;
	}
#else
	LEDA_GRAPH<leda::node,leda::edge> g = TRANSITIVE_CLOSURE(*this);
	
	// Translate the transitive closure graph into a matrix form
	forall_nodes(v, g)
	{
		i = GetNodeDFSIndex(g[v]);

		forall_adj_nodes(u, v)
			transClosMat[i][GetNodeDFSIndex(g[u])] = 1;
	}
#endif	
}

void DAG::ComputeNodesMass()
{
	leda::node v;
	ASSERT(GetNodeCount() == transClosMat.GetSize());

	/*forall_nodes(v, *this)
		GetNode(v)->SetMass(transClosMat.RowSum(GetNodeDFSIndex(v)));*/
	
	SmartArray<int> ones(transClosMat.NRows());
	ones.Set(1);
	
	SmartArray<int> relMasses = transClosMat * (transClosMat * ones);
	
	forall_nodes(v, *this)
		GetNode(v)->SetMass(relMasses[GetNodeDFSIndex(v)]);
}

/*!
	\brief Compute all the derived values in the graph.

	The derived values are those which depend on other values and don't need
	to be saved. These values are in fact actually read and saved but only in 
	order to improve computation performance.
*/
void DAG::ComputeDerivedValues()
{
	SmartArray<leda::node> roots = GetAllRoots();
	int i, nDFSIndex = ROOT_NODE_DFS_INDEX;
	leda::node v;
	
	ASSERT(roots.GetSize() > 0); // There has to be at least one root
	
	ResetNodesInfo(); // Make sure the node info is clean
	
	for (i = 0; i < roots.GetSize(); i++)
		nDFSIndex = ComputeNodesInfo(roots[i], ROOT_NODE_LEVEL, nDFSIndex);
	
	ComputeTransiveClosure();
	ComputeNodesMass();
	
	nCumulativeMass = 0;
	nMaxBFactor = 0;
	dTotalTSVSum = 0;

	for (i = 0; i < roots.GetSize(); i++)
	{
		v = roots[i];
		ComputeTSVs(v);
		nCumulativeMass += GetNodeMass(v);
	}
	
	forall_nodes(v, *this)
	{
		GetNode(v)->ComputeDerivedValues();

		dTotalTSVSum += GetNodeTSVNorm(v);

		if (nMaxBFactor < outdeg(v))
			nMaxBFactor = outdeg(v);
	}
}

/*!
	This function separates a subtree rooted at <v> and returs it
	in the <subg> parameter.

	The transitive closure matrix and the cumulative mass of the entire
	graph is retained in both subgraphs.

	\param v root of the subtree.
	\param bRecomputeTSVs controls whether the TSVs that remains in 
	the original graph should be recomputed.
*/
DAGPtr DAG::SplitSubGraph(leda::node v, bool bRecomputeTSVs)
{
	DAGPtr subg = CreateObject();

	// if v is not a root...
	if (indeg(v) > 0)
	{
		// Create a map for the parents of nodes with more 
		// than one parent
		NodeMap nodeParentMap;
		
		SplitSubGraph(v, subg, nodeParentMap, nil);
		
		// Now it is safe to eliminate those nodes in the 
		// subtre that have multiple parents
		forall_defined(v, nodeParentMap)
			del_node(v);

		subg->transClosMat = GetTransClosMat();
		subg->nCumulativeMass = nCumulativeMass;

		if (bRecomputeTSVs) 
		{
			SmartArray<leda::node> roots = GetAllRoots();
			for (int i = 0; i < roots.GetSize(); i++)
				ComputeTSVs(roots[i]);
		}
	}
	else
	{
		*subg = *this;
		del_all_nodes();
	}
		
	return subg;
}

/*!
	Recursive function called by SplitSubGraph(leda::node, bool).
	It recursively adds a subgraph to the graph pointed by <ptrDag>
	and then it deletes the subgraph from the original graph.
	
	Note: if the desired subtree contains nodes with multiple parents,
	we need special care. There are two cases, (a) more than one of the
	multiple parents of a node are also contained in the subtree; (b) only 
	one of the node's parents is in the subtree. Case (b) is not a problem, 
	but case (a) requires keeping track of the multiple parents and avoding 
	deleting these nodes before traversing all their insident edges.
*/
leda::node DAG::SplitSubGraph(leda::node v, DAGPtr& ptrDag, NodeMap& nodeParentMap, leda::edge inEdge)
{
	leda::node a, b;
	leda::edge e;
	double w;
	bool bMultiParents;
	
	ASSERT(inEdge == nil || target(inEdge) == v);
	
	if (nodeParentMap.defined(v))
	{
		bMultiParents = true; // it either has OR had multiple parents
		a = nodeParentMap[v];
	}
	else
	{
		bMultiParents = indeg(v) > 1; // it has multiple parents
		a = ptrDag->NewNode(inf(v));
		
		if (bMultiParents)
			nodeParentMap[v] = a;
	}

	forall_adj_edges(e, v)
	{
		w = inf(e);
		b = SplitSubGraph(target(e), ptrDag, nodeParentMap, e);
		ptrDag->NewEdge(a, b, w);
	}

	// NOTE: nodes with multiple parents will be eliminated at the end,
	// but the edge that brought us to the node must be eliminated now.
	if (!bMultiParents)
		del_node(v);
	else if (inEdge != nil)
		del_edge(inEdge);

	return a;
}

//! Create an edge between vertices <i> and <j> of weight <dist>.
void DAG::AddEdge(Matrix& adj, leda::edge e, int srcNodeIdx, int tgtNodeIdx) const
{
	double w = GetEdgeWeight(e);

	adj(srcNodeIdx, tgtNodeIdx) = w;
	adj(tgtNodeIdx, srcNodeIdx) = -w;
}

/*!
	\brief Computes eigen-sum over a given adjacency matrix.

	We are supposed to compute the $\delta(T_i) - 1$ largest absolute eigenvalues, 
	where $\delta(T_i)$ is the degree of the root node of the tree $T_i$.

	The justification for summing the $\delta(T_i) - 1$ largest absolute eigenvalues follows:

	a) the largest absolute eigenvalues are the most informative of the subgraph structure.

	b) by summing $\delta(T_i) - 1$ elements we normalize the sum according to the local
	   complexity of the subgraph root.

	According to the Overton and Womersley's Theorem, for the sum of the first 
	k eigenvalues of a symmetric matrix A, the following semidefinite programming
	characterization holds:

	\[\lambda_1(A) + \hdots + \lamda_k(A) = max A \dot U  s.t. trace(U) = k for 0 \le U \le I.\]

	However, for our particular case, we can't use this theorem because A is antisymmetric.
	Therefore, for this implementation, we compute the SVD to obtain the eigenvalues and then 
	all the values are summed up, instead of only the $\delta(T_i) - 1$ largest ones.

	\todo Find out the impact of summing all the eigenvalues instead of the $k$ largest ones.5s3l1

*/
double DAG::ComputeEigenSum(const Matrix& adj, int nVals)
{
	using namespace NEWMAT;

	int n = adj.Nrows();
	Matrix U(n, n);
	DiagonalMatrix D(n);
	double td = 0.0;
	//double eigenVals[n];

	//nVals *= 2; // Values will be duplicated, so step will be 2.
	//ASSERT(nVals <= n);

	SVD(adj, D, U);
	
	for (int i = 1; i <= nVals; i++)
		td += D(i, i);

	return td;
}

double DAG::ComputeLaplacian(const Matrix& adj)
{
	using namespace NEWMAT;

	int n = adj.Nrows(), i, j;
	SymmetricMatrix L(n);
	ColumnVector T(n);

	for (i = 1; i <= n; i++)
	{
		T(i) = 0;

		for (j = 1; j <= n; j++)
			T(i) += fabs(adj(i, j));
	}

	for (i = 1; i <= n; i++)
	{
		for (j = 1; j <= n; j++)
		{
			if (i == j)
				L(i, j) = 1;
			else if (adj(i, j) != 0)
				L(i, j) = -1 / (sqrt(T(i) * T(j)));
			else
				L(i, j) = 0;
		}
	}

	DiagonalMatrix D(n);

	Jacobi(L, D);

	SortAscending(D);

	return fabs(D(1,1));
}

/*!
	Computes all the nodes' TSV's.
*/
double DAG::ComputeTSVs(leda::node root)
{
	int j = 1, n = GetNodeMass(root);
	ASSERT(n > 0);
	
	NodeIndexMap loopyNodeMap;
	Matrix adj(n, n);
	adj = 0.0;

	return ComputeTSVs(root, j, adj, loopyNodeMap);
}

/*!
	step 1: for all adjacent nodes w of the given root node, if w has 
	in-degree grater that 1 (multiple parents) and has already been visited, 
	bHasLoopyChild is set to true so that a new adj matrix will be computed,
	in which v is the root and use it to compute v's TSV. Note that 
	if j is the root of the matrix (i = j = 1), we don't care whether v has
	loopy children or not. Otherwise we would enter in an endless loop.

	step 2: add the edges to the adj matrix. Edges immediately below 
	a loopy	node must be added to the adj matrix, but those below them
	do not because they have already been added.

	step 3: compute v's eigen sum if necessary.
*/
double DAG::ComputeTSVs(leda::node v, int& j, Matrix& adj, NodeIndexMap& loopyNodeMap)
{
	using namespace NEWMAT;

	leda::node w;
	leda::edge e;
	double s;
	int i = j;
	bool bDefined, bHasLoopyChild = false;
	
	TSV& tsv = GetNodeTSV(v);
	tsv.ReSize(outdeg(v), true);

	forall_adj_edges(e, v)
	{
		w = target(e);
		bDefined = indeg(w) > 1 ? loopyNodeMap.defined(w):false;

		if (bDefined)
		{
			bHasLoopyChild = true;
			AddEdge(adj, e, i, loopyNodeMap[w]);
			s = GetEigenLbl(w);
		}
		else
		{
			AddEdge(adj, e, i, ++j);

			if (indeg(w) > 1)
				loopyNodeMap[w] = j; // loopy node visited

			s = ComputeTSVs(w, j, adj, loopyNodeMap);
		}

		tsv.Add(s);
	}

	// If v has loopy children and it isn't the root of the adj matrix
	// we'll compute its TSV sum by creating a new adj mat rooted at v.
	if (bHasLoopyChild && i > 1) 
		s = ComputeTSVs(v);
	else
	{
		s = ComputeEigenSum(adj.SubMatrix(i, j, i, j), outdeg(v));
		//s = ComputeLaplacian(adj.SubMatrix(i, j, i, j));

		SetEigenLbl(v, s);
		tsv.Sort();
	}

	return s;
}

leda::node DAG::GetNode(int nIndex)					
{ 
	const leda::list<leda::node>& nodeList = all_nodes();
	leda::list_item it = nodeList.get_item(nIndex); 

	return nodeList[it];
}

//! Returns the first root in the DAG. i.e. the first node v that has indeg(v) = 0
leda::node DAG::GetFirstRootNode() const
{
	leda::node v;

	forall_nodes(v, *this)
		if (indeg(v) == 0)
			return v;
	
	return nil;
}

//! Returns the next root in the DAG. i.e. the first SUCCESSOR of node v that has indeg(v) = 0
leda::node DAG::GetNextRootNode(leda::node v) const
{
	do {
		v = succ_node(v);
	} while (v != nil && indeg(v) != 0);
	
	return v;
}

//! Returns a SmartArray with all the graph's root nodes.
SmartArray<leda::node> DAG::GetAllRoots() const
{
	SmartArray<leda::node> roots;
	leda::node v;

	forall_nodes(v, *this)
		if (indeg(v) == 0)
			roots.AddTail(v);
	
	return roots;
}

/*static*/
/*leda::node DAG::WeightNodeSimilarity(const double& sim, leda::node v) const
{
	// The cumulative mass goes from 2N - 1 for the star graph
	// to N(N+1)/2 for a path graph.
	// We can convert the relative mass of a node to a value between 0 and 1
	// by the cumulative distribution function for a continuous uniform distribution.
	// (x - a) / (b - a)
	
	double relMass = GetNodeMass(v) / double(nCumulativeMass); // relative mass for the node
	int n = GetNodeCount();
	
	double w = (relMass - 2 * n + 1) / ((n * n) / 2 - 1.5 * n + 1);
	
	return .9 * sim + .1 * w;
}*/

/*! 
	@brief Static helper function for DAG::Match(). 
	
	It creates a bipartite graph between g1 and g2's nodes,
	as long as they weren't previously matched, in which edges are directed 
	from a smaller to a larger set. 
	
	It the solution to the max weight max card bipartite matching and returns
	the largest weight edge from the matching.

	@return largest weight edge from the MWMCB matching.
*/
/*static*/ leda::edge DAG::ComputeBipartiteGraph(LEDA_GRAPH<leda::node, double>& G, const DAG& g1, const DAG& g2,
									 SimMatrix& simMat, const DAGNodeMap& paired1,
									 const DAGNodeMap& paired2, MatchedNodePair* pMatchedPair /*= NULL*/)
{
	//g1 is the query graph and g2 is the model graph
	leda::node u, v, g1Node, g2Node;

	leda::list<leda::node> A, B;
	leda::list_item lia, lib;

	// Create a graph with n1 + n2 nodes
	forall_nodes(u, g1)
		if (pMatchedPair->GetQueryNode() != g1.GetNodeDFSIndex(u))
			A.push(G.new_node(u));

	forall_nodes(u, g2)
		if (pMatchedPair->GetModelNode() != g2.GetNodeDFSIndex(u))
			B.push(G.new_node(u));

	leda::node_array<leda::node>& nodes = G.node_data();

	// Create a bipartite graph between A and B nodes.
	// do not establish edges between nodes that have been matched.
	for (lia = A.first(); lia != nil; lia = A.succ(lia))
	{ 
		// pick node from A
		u = A[lia];
		g1Node = nodes[u];

		if (paired1.defined(g1.GetNode(g1Node)))
			continue;	// loopy node that has been matched
		

		for (lib = B.first(); lib != nil; lib = B.succ(lib))
		{ 
			// pick node from B
			v = B[lib];
			g2Node = nodes[v];

			if (paired2.defined(g2.GetNode(g2Node)))
				continue;	// loopy node that has been matched

			// establish an edge in proper direction if needed
			//if (g1.AreNodesRelated(g1Node, g2, g2Node)) 
			if (simMat[g1.GetNodeDFSIndex(g1Node)][g2.GetNodeDFSIndex(g2Node)] > 0.0)
				G.new_edge(u, v, 1.0);
		}
	}
		
	//DBG_MSG1("Number of edges in Bipartite graph: ", G.number_of_edges());
	
	if (G.number_of_edges() == 0)
		return nil;

	// Assign weights to edges
	leda::edge_array<double> edgeCosts(G);
	leda::edge e;
	int i, j;

	ASSERT(s_matchParams.dBreakSibRelPen >= 0 && s_matchParams.dBreakSibRelPen <= 1);

	forall_edges(e, G)
	{
		g1Node = nodes[G.source(e)];
		g2Node = nodes[G.target(e)];

		i = g1.GetNodeDFSIndex(g1Node);
		j = g2.GetNodeDFSIndex(g2Node);
		
		if (DAG::GetMatchParams().nPreserveAncestorRel)
		{
			if (!pMatchedPair->IsEmpty() && !pMatchedPair->AnsestorRelPreserved(i, j))
				simMat[i][j] = 0;
		}
		
		edgeCosts[e] = simMat[i][j];
	}

	/*GraphWin gw(G);
	gw.display();
	gw.place_into_win();
	gw.edit();*/

	// Solve the max weight max card bipartite matching
	
	// Note from the LEDA manual about the 'arithmetic demand' of MWA_SCALE_WEIGHTS:
	// "This scaling function should be used for the algorithm that compute minimum of 
	// maximum weight assignments or maximum weighted matchings of maximum cardinality." 
	// while MWBM_SCALE_WEIGHTS "is appropriate for the maximum weight matching algorithm.".
	
	//MWA_SCALE_WEIGHTS(G, edgeCosts); // also see Notes.txt
	//leda::list<leda::edge> tmch = MWMCB_MATCHING(G, A, B, edgeCosts);
	
	std::cerr << "Please ignore this warnings [ ";
	//scale_weights(G, edgeCosts, 3.0);
	MWBM_SCALE_WEIGHTS(G, edgeCosts);
	
	leda::list<leda::edge> tmch = MAX_WEIGHT_BIPARTITE_MATCHING(G, edgeCosts);
	//leda::list<leda::edge> tmch = MAX_WEIGHT_BIPARTITE_MATCHING_T(G, edgeCosts);

	std::cerr << " ]" << std::endl;
	
	// Find the largest weight edge from the matching
	leda::list_item lie;
	double    sim, curr_weight, max_weight = -9999999;
	leda::edge curr_edge, max_edge = nil;

	const double g1CMass = g1.nCumulativeMass;  //cumulative mass for query
	const double g2CMass = g2.nCumulativeMass;  //cumulative mass for model

	double g1RelMass, g2RelMass;

	// Only for debug
	leda::node a = nil, b = nil, c = nil, d = nil;
	double x,y, maxsim = -1;
	int m1, m2;
	// end debug

	const double smw = s_matchParams.dSimilMassWeight;
	ASSERT(smw >=0 && smw <= 1);

	/*
		Note that node masses are relative to the original graph the belong to
		so the relative mass should be according to the size of that graph or to the root
		of the graph the belong to.
	*/

	forall_items(lie, tmch) 
	{
		curr_edge = tmch[lie];

		g1Node = nodes[G.source(curr_edge)];
		g2Node = nodes[G.target(curr_edge)];

		g1RelMass = g1.GetNodeMass(g1Node) / g1CMass;
		g2RelMass = g2.GetNodeMass(g2Node) / g2CMass;

		ASSERT(g1RelMass >= 0 && g1RelMass <= 1);
		ASSERT(g2RelMass >= 0 && g2RelMass <= 1);

		sim = edgeCosts[curr_edge];

		//curr_weight = s_smw * sim * g1RelMass + (1 - s_smw) * sim * g2RelMass;
		curr_weight = smw * sim + (1 - smw) * MAX(g1RelMass, g2RelMass);
		
		ASSERT(curr_weight >=0 && curr_weight <= 1)
		
		if (DAG::IsDbgMode() && sim > maxsim)
		{
			maxsim = sim;
			c = g1Node;
			d = g2Node;
		}

		if (curr_weight > max_weight) 
		{
			max_weight = curr_weight;
			max_edge = curr_edge;
			
			// only for debug
			if (DAG::IsDbgMode())
			{
    			a = g1Node;
    			b = g2Node;
    			x = edgeCosts[curr_edge];
    			y = MAX(g1RelMass, g2RelMass);
    			m1 = g1.GetNodeMass(g1Node);
    			m2 = g2.GetNodeMass(g2Node);
    		}
		}
	}

	// Debug stuff
	if (DAG::IsDbgMode() && max_edge != nil)
	{
		logFile.Print("\nQuery node: %s (%d)", (LPCTSTR) g1.GetNodeLbl(a), g1.GetNodeDFSIndex(a));
		logFile.Print("\nModel node: %s (%d)", (LPCTSTR) g2.GetNodeLbl(b), g2.GetNodeDFSIndex(b));
   		
    	if (a != c || b != d)
		{
			logFile.Print("\nTrue max Q node: %s (%d)", (LPCTSTR) g1.GetNodeLbl(c), g1.GetNodeDFSIndex(c));
			logFile.Print("\nTrue max M node: %s (%d)", (LPCTSTR) g2.GetNodeLbl(d), g2.GetNodeDFSIndex(d));
    	}
    	
   		logFile << "\nM1: " << m1 << "\nM2: " << m2;
   		logFile << "\nSim: " << x << "\nRel Mass: " << y << "\nWSim: " << max_weight << "\n\n";
	}

	/*if (max_edge != nil && max_weight >= s_matchParams.dMinNodeSim)
		G[max_edge] = max_weight;
	else if (max_edge != nil)
		max_edge = nil;*/

	return max_edge;
}

/*!
	@brief The main matching algorithm for various DAG's.

	It is a static function that receives two DAG's as input 
	and computes a value representing their similarity, as well as 
	a list of corresponding nodes in the two DAG's.  Uses two
	node functions, <node_similarity> and <nodes_are_related>.

	@see NodeSimilarity AreNodesRelated ComputeBipartiteGraph
*/
/*static*/
double DAG::Match(DAG& g1, DAG& g2, SimMatrix& simMat, DAGNodeMap& paired1,
				  DAGNodeMap& paired2, MatchedNodePair* pMatchedPair /*= NULL*/)
{
	int n1 = g1.GetNodeCount();
	int n2 = g2.GetNodeCount();

	// Edges should be directed from a smaller to a larger set.
	// but s_rwm and swm are set assuming g1 is the query!!!!! so...
	/*if (n1 > n2)
	{
		DBG_MSG("SWAPPING...");
		ASSERT(false); //transpose simMat
		if (pMatchedPair)
			pMatchedPair->SwapNodes();

		return Match(g2, g1, simMat, paired2, paired1, pMatchedPair);
	}*/

	// Check that both graphs are of non-zero size
	if (n1 * n2 == 0)
		return 0;

	// Check that both sizes are >1, since if one of them
	// is of size one, it has been matched (unless very first call).
	if (((n1 == 1) || (n2 == 1)) && paired1.size() > 0)
		return 0;

	// Create a bipartite graph between A and B nodes and 
	// check that the number of edges is non-zero.
	LEDA_GRAPH<leda::node, double> G; // graph of meta-nodes.

	leda::edge max_edge = ComputeBipartiteGraph(G, g1, g2, simMat, paired1, paired2, pMatchedPair);

	if (max_edge == nil)
		return 0;

	// Take the endpoints of the max weight edge
	leda::node maxG1Node, maxG2Node;
	leda::node_array<leda::node>& nodes = G.node_data();

	maxG1Node = nodes[G.source(max_edge)];
	maxG2Node = nodes[G.target(max_edge)];
	
	// Append the found nodes to the matching dict.
	// mark one with another's global number.
	DAGNodePtr maxG1NodePtr = g1.GetNode(maxG1Node);
	DAGNodePtr maxG2NodePtr = g2.GetNode(maxG2Node);

	paired1[maxG1NodePtr] = maxG2NodePtr;
	paired2[maxG2NodePtr] = maxG1NodePtr;

	//MatchedNodePair matNodes(maxG1Node, g1, maxG2Node, g2);

	// Get the DFS index BEFORE spliting the graphs
	int idxG1Node = g1.GetNodeDFSIndex(maxG1Node);
	int idxG2Node = g2.GetNodeDFSIndex(maxG2Node);
	
	// Compute rooted subgraphs and remainders
	// We want to leave the splitting node in both subgraphs
  DAGPtr sg1 = g1.SplitSubGraph(maxG1Node);
	DAGPtr sg2 = g2.SplitSubGraph(maxG2Node);
	
	double max_weight = simMat[idxG1Node][idxG2Node];
	
	////////////////////////////////////////////////////////////////
	//For display purposes, save the similarity assigned to the node
	DAGNode* pNode = (DAGNode*)(const DAGNode*)maxG2NodePtr;
	pNode->SetSimilarity(max_weight);
	////////////////////////////////////////////////////////////////

	// Recurse on the root-matched subtree first...
	pMatchedPair->SetEmpty();
	max_weight += Match(*sg1, *sg2, simMat, paired1, paired2, pMatchedPair);

	// ...set the info needed to update the similarity matrix ...
	pMatchedPair->SetNodes(idxG1Node, idxG2Node);
	
	if (g1.GetNodeCount() > 0 && g2.GetNodeCount() > 0)
	{
		// ...update sim matrix to avoid breaking hierarchical constraints...
		pMatchedPair->UpdateSimMat(simMat, 1 - s_matchParams.dBreakSibRelPen);

		// ...and then recurse on the complement graph
		max_weight += Match(g1, g2, simMat, paired1, paired2, pMatchedPair);
	}

	return max_weight;
}

/*!
	\brief Returns the similarity between two dags.

	The similarity is computed by matching image's and model's nodes and 
	then normalizing the number of matched nodes according to the order of
	the model graph.

	/see NodeSimilarity AreNodesRelated ComputeBipartiteGraph Match
*/ 
/*static*/
double DAG::Similarity(const DAG& query, const DAG& model, 
		DAGNodeMap& nodeMap, bool bRecomputeTSVs)
{

// 	See who's being matched:
// 	DBG_VAL(model.GetObjName())
// 	DBG_VAL(model.GetViewNumber())
	

	DAGNodeMap nodeMap2;
	leda::node q, m;
	int i, j;

	int n1 = query.GetNodeCount();
	int n2 = model.GetNodeCount();
	int min = MIN(n1, n2);

	// Create a similarity matrinx and a correspondence matrix
	SmartMatrix<double> simMat, dbgMat;

	simMat.ReSize(n1, n2, true);

	// Populate the similarity matrix
	forall_nodes(q, query)
	{
		forall_nodes(m, model)
		{
			i = query.GetNodeDFSIndex(q);
			j = model.GetNodeDFSIndex(m);
			
			simMat[i][j] = query.NodeSimilarity(q, model, m);
		}
	}
	
	if (DAG::IsDbgMode())
	{
		static bool bLogFileCreated = false;
		
		if (!bLogFileCreated)
		{
			bLogFileCreated = true;
			logFile.open("matching.log", std::ios::out | std::ios::trunc);
		}
		
    	logFile.Print("\nQuery: %s\nModel: %s\n\n", (LPCTSTR) query.GetDAGLbl(), (LPCTSTR) model.GetDAGLbl());
		simMat.Print(logFile, 6, 3, true);
		logFile.Print("\n\n");
		s_matchParams.Print(logFile);
		dbgMat = simMat;
	}

	// Make sure nodeMap is empty
	nodeMap.clear(); 

	// DAG::Match() does not preserve the contents of its parameters so
	// we need to make a copy of the graphs.
	DAGPtr cpQuery = query.CreateObject();
	DAGPtr cpModel = model.CreateObject();
	
	MatchedNodePair matNodes(query, model); // graphs use to apply penalties
	
	double dMatchedNodes = DAG::Match(*cpQuery = query, *cpModel = model, 
		simMat, nodeMap, nodeMap2, &matNodes);
		
	if (DAG::IsDbgMode())
	{		
		for (int i = 0; i < simMat.NRows(); i++)
			for (int j = 0; j < simMat.NCols(); j++)
				if (simMat[i][j] - dbgMat[i][j] != 0.0)
				{
					if (simMat[i][j] != 0.0)
						dbgMat[i][j] = simMat[i][j] - dbgMat[i][j];
				}
				else
					dbgMat[i][j] = 0.0;

		dbgMat.Print(logFile, 6, 3, true);
	}

	return (dMatchedNodes / n1 + dMatchedNodes / n2) / 2.0;
	//return dMatchedNodes / n1;
}

/*!
	\brief Returns the match distance between the two dags. 
	
	This function is simply the inverse of the Similarity function

	\see Similarity
*/
/*static*/
double DAG::Distance(const DAG& image, const DAG& model, 
	DAGNodeMap& nodeMap, bool bRecomputeTSVs)
{
	return 1.0 / DAG::Similarity(image, model, nodeMap, bRecomputeTSVs);
}


/////////////////////////////////////////////////////////////////////////////////////
// Virtual functions

//! Assignment operator for the DAG class. Must be overwritten in the derived classes.
DAG& DAG::operator=(const DAG& rhs)
{
	Clear();

	DAG_BASE_CLASS::operator=(rhs);

	nFileOffset = rhs.nFileOffset;
	nDAGId      = rhs.nDAGId;
	strGraphLbl = rhs.strGraphLbl;

	// Derived values
	nMaxBFactor   = rhs.nMaxBFactor;
	dTotalTSVSum  = rhs.dTotalTSVSum;
	nViewNumber   = rhs.nViewNumber;
	strObjectName = rhs.strObjectName;
	eigenVals	  = rhs.eigenVals; // Experimental. It isn't necessary
	transClosMat  = rhs.transClosMat;
	nCumulativeMass = rhs.nCumulativeMass;

	return *this;
}

//! It sets to zero all the DAG's member variables
void DAG::Clear()
{ 
	DAG_BASE_CLASS::clear(); 

	nFileOffset = 0;
	nDAGId = 0;
	strGraphLbl.Clear();

	// Derived values
	nMaxBFactor = 0;
	dTotalTSVSum = 0;
	nViewNumber = 0;
	strObjectName.Clear();
	transClosMat.Clear();
	nCumulativeMass = 0;
}

/// Similarity function for two nodes.
double DAG::NodeSimilarity(leda::node u, const DAG& from, leda::node v) const
{
	double dNodeSim  = 1 - NodeDistance(u, from, v); 
	double dEigenSim = NodeTSVSimilarity(u, from, v);
	double dSimilarity;
	
	if (dNodeSim < 0 || dNodeSim > 1)
	{
		if (DAG::IsDbgMode())
		{
			std::cerr << "\nWrong similarity value (" << dNodeSim << ") when matching: \n"
				<< GetDAGLbl() << ":" << GetNodeLbl(u) << " against "
				<< from.GetDAGLbl() << ":" << from.GetNodeLbl(v) << std::endl;
		}
			
		dNodeSim = (dNodeSim < 0) ? 0:1;
	}

	ASSERT(!isnan(dNodeSim)); ASSERT(finite(dNodeSim));
	ASSERT(!isnan(dEigenSim)); ASSERT(finite(dEigenSim));

	ASSERT(dNodeSim >= 0 && dNodeSim <= 1);
	ASSERT(dEigenSim >= 0 && dEigenSim <= 1);

	if (dNodeSim > 0.0)
		dSimilarity = s_matchParams.dTSVSimWeight * dEigenSim + (1 - s_matchParams.dTSVSimWeight) * dNodeSim;
	else
		dSimilarity = 0.0;
	
	/*switch (s_matchParams.nSimFuncType)
	{
		case EXP_NEG_SUM: 
			dSimilarity = exp (-dNodeDistance - dEigenDistance);
			break;
		case INV_NODE_DIST: 
			dSimilarity = 1.0 / (1 + dNodeDistance);
			break;
		case INV_NODE_DIST_AND_TSV:
			dSimilarity = 1.0 / 
				(1 + sqrt(dNodeDistance * dNodeDistance + dEigenDistance * dEigenDistance));
			break;
		case ID_NODE_DIST: 
			dSimilarity = dNodeDistance;
			break;
		default:
			ASSERT(false);
	}*/

	ASSERT(!isnan(dSimilarity)); ASSERT(finite(dSimilarity));
	ASSERT(dSimilarity >= 0); ASSERT(dSimilarity <= 1);

	return dSimilarity;
}


/// Eigen-value based distance function for two nodes.
double DAG::EigenDistance(leda::node u, const DAG& from, leda::node v) const
{
	//return NodeTSVSimilarity(u, from, v);

	return Norm2(GetNodeTSV(u) - from.GetNodeTSV(v));
}


/*! 
	Print out vertices, their labels and eigenlabels, and the
	adjacency structure of the graph.
*/
void DAG::Print(std::ostream& os /*=cout*/) const
{
	os << "\nDAG label: " << GetDAGLbl() << std::endl;
	os << "Object: " << GetObjName() << ", view: " << GetViewNumber() << std::endl;

	/*
	os << "\n--------------- Vertex Set ---------------\n"; */ 

	leda::node v;

	forall_nodes(v, *this)
	{
		os << "\nNode:" << index(v) << std::endl; 
		GetNode(v)->Print(os);
	}
}

//! Prints the set of all TSVs in the graph.
void DAG::PrintTSVs(std::ostream& os /*=cout*/) const
{
	os << "\nDAG's adj matrix eigen values:\n";
	
	eigenVals.Print(os);

	os << "\n\n--------------- TSV Set ---------------\n";
	
	leda::node v;

	forall_nodes(v, *this)
	{
		os << "\nNode " << GetNodeLbl(v) << ": ";
		GetNodeTSV(v).Print(os);
	}

	os << "\n---------------------------------------\n";
}

//! Reads a DAG from the specified stream.
std::istream& DAG::Read(std::istream& is, bool bOnlyDataForMatching /*= false*/)
{
	int n, i;
	leda::node v, s, t;
	leda::map<long, leda::node> map;
	double w;
	String strClassName;

	Clear();

	// Save the offset where we are about to read from
	nFileOffset = is.tellg();

	strClassName.Read(is);

	if (ClassName() != strClassName) {
		DBG_MSG1("Reading DAG from offset:", nFileOffset);
		THROW_EXCEPTION("Incorrect DAG class");
	}

	strGraphLbl.Read(is);

	is.read((char*) &n, sizeof(n));

	for(i = 0; i < n; i++)
	{
		is.read((char*) &v, sizeof(v));
		map[(long)v] = NewNode(ReadNode(is));
	}

	is.read((char*) &n, sizeof(n));

	for(i = 0; i < n; i++)
	{
		is.read((char*) &s, sizeof(s));
		is.read((char*) &t, sizeof(t));
		is.read((char*) &w, sizeof(w));
		NewEdge(map[(long)s], map[(long)t], w);
	}
	
	// Read derived values
	is.read((char*) &nMaxBFactor, sizeof(nMaxBFactor));
	is.read((char*) &dTotalTSVSum, sizeof(dTotalTSVSum));
	is.read((char*) &nViewNumber, sizeof(nViewNumber));
	strObjectName.Read(is);
	eigenVals.Read(is);  // Experimental. It isn't necessary

	transClosMat.Read(is);
	is.read((char*) &nCumulativeMass, sizeof(nCumulativeMass));
	ASSERT(nCumulativeMass >= GetNodeCount());

	return is;
}

//! Writes the DAG to the specified stream.
std::ostream& DAG::Write(std::ostream& os) const
{
	leda::node v, s, t;
	leda::edge e;
	double w;

	// Save the offset where we are about to write to ???
	//nFileOffset = os.tellp();

	ClassName().Write(os);
	strGraphLbl.Write(os);

	int n = GetNodeCount();
	
	os.write((char*) &n, sizeof(n));

	forall_nodes(v, *this)
	{
		os.write((char*) &v, sizeof(v));
		os << inf(v);
	}

	n = GetEdgeCount();

	os.write((char*) &n, sizeof(n));

	forall_edges(e, *this)
	{
		s = source(e);
		t = target(e);
		w = inf(e);
		os.write((char*) &s, sizeof(s));
		os.write((char*) &t, sizeof(t));
		os.write((char*) &w, sizeof(w));
	}

	// Read derived values
	os.write((char*) &nMaxBFactor, sizeof(nMaxBFactor));
	os.write((char*) &dTotalTSVSum, sizeof(dTotalTSVSum));
	os.write((char*) &nViewNumber, sizeof(nViewNumber));
	strObjectName.Write(os);
	eigenVals.Write(os); // Experimental. It isn't necessary

	transClosMat.Write(os);
	os.write((char*) &nCumulativeMass, sizeof(nCumulativeMass));

	return os;
}

//! Delete and edge between vertices <i> and <j>.
void DAG::DelEdge(Matrix& adj, int srcNodeIdx, int tgtNodeIdx) const
{
	adj(srcNodeIdx, tgtNodeIdx) = 0.0;
	adj(tgtNodeIdx, srcNodeIdx) = 0.0;
}

/*! 
	@brief Computes the TSV similarity between two nodes

	The formula used is: 
		
	\[ \fraq{|tsv_1 - tsv_2|}{\mbox{max} (|tsv_1|, (|tsv_2|)} \]
	
	If both terms are zero, the similarity returned is 1.
	If the numerator is grater than or equal to the numerator,
	the TSV similarity is 0.
*/
double DAG::NodeTSVSimilarity(leda::node u, const DAG& from, leda::node v) const
{
	const TSV& tsv1 = GetNodeTSV(u);
	const TSV& tsv2 = from.GetNodeTSV(v);	
	double n1 = GetNodeTSVNorm(u);
	double n2 = from.GetNodeTSVNorm(v);

	if (n1 == 0 && n2 == 0)
		return 1;

	TSV diff = tsv1 - tsv2;
	double diffNorm = diff.Norm2();

	double max = MAX(n1, n2);

	// max should be a normalization factor to ensure
	// the the result would be between zero and 1. However, 
	// it is possible that diffNorm is grater than max, and so
	// we need to handle this case somehow. Clearly, a better solution
	// would be to use a better normalization factor.
	if (diffNorm > max)
		return 0;

	return 1 - diffNorm / max;
}

/*! 
	@brief Deletes the subtree rooted at <v>

	\todo This function does not support some kinds of DAGs. More specifically,
	DAGs with loops are not allowed. This must be fixed in the future.
*/
void DAG::DeleteSubDAG(leda::node v)
{
	leda::node u;

	forall_adj_nodes(u, v)
		DeleteSubDAG(u);

	del_node(v);
}

/*! 
	@brief Prints the DAG's adjacency matrix
	
	FIXME It only works where there is one single root node
*/
void DAG::PrintAdjMatrix(std::ostream& os)
{
	int j = 1, n = GetNodeCount();
	NodeIndexMap loopyNodeMap;
	Matrix adj(n, n);

	adj = 0.0;
	ComputeTSVs(GetFirstRootNode(), j, adj, loopyNodeMap);

	os << "\nAdjacency matrix:\n";
	os << std::setw(1) << std::setprecision(0) << adj << std::setprecision(6) << std::endl;
}

/*
	Depth First Search
*/
leda::node DAG::DFSGetNext(leda::node v, EdgeList& l, bool bVisitNode) const
{
	leda::edge e;

	ASSERT(v != nil);

	if (bVisitNode)
		forall_adj_edges(e, v)
			l.append(e);

	return l.empty() ? nil:target(l.pop_back());
}

/*
	Breadth First Search
*/
leda::node DAG::BFSGetNext(leda::node v, EdgeList& l, bool bVisitNode) const
{
	leda::edge e;

	ASSERT(v != nil);

	if (bVisitNode)
		forall_adj_edges(e, v)
			l.append(e);

	return l.empty() ? nil:target(l.pop_front());
}

void SimilarityGraph::Add(leda::node g1Node, leda::node g2Node, double dSimilarity)
{
	new_edge(MapNode(g1Map, g1Node), MapNode(g2Map, g2Node), dSimilarity);
}

/*void SimilarityGraph::Add(leda::node g1Node, leda::node g2Node, double dSimilarity,
						  leda::node g1MatchedNode, leda::node g2MatchedNode)
{
	leda::node u = MapNode(g1Map, g1Node);
	leda::node v = MapNode(g2Map, g2Node);

	leda::edge e = new_edge(u, v, dSimilarity);

	g1MatchedNodeMap[e] = g1MatchedNode;
	g2MatchedNodeMap[e] = g2MatchedNode;
}*/

/*
	Remaining node, deleting node
*/
void SimilarityGraph::MergeG2Nodes(leda::node remNode, leda::node delNode)
{
	leda::node r, d;
	leda::edge e;

	r = g2Map[remNode];
	d = g2Map[delNode];

	forall_in_edges(e, d)
		move_edge(e, source(e), r);

	g2Map.undefine(delNode);
	del_node(d);
}

double SimilarityGraph::ComputeBestAssignment(METHOD nMethod, int* pnMatched, 
											  NodeMap* pNodeMap /*= NULL*/)
{
	leda::edge_array<double> edgeCosts(*this, 0.0);
	leda::edge e;

	forall_edges(e, *this)
	{
		edgeCosts[e] = inf(e);
		operator[](e) = 1.0;
	}

	EdgeList l;

	if (nMethod == MWMC)
	{
		MWA_SCALE_WEIGHTS(*this, edgeCosts);
		
		NodeList setA, setB;
		bool bBipartite = Is_Bipartite(*this, setA, setB);

		ASSERT(bBipartite);
		l = MWMCB_MATCHING(*this, setA, setB, edgeCosts);
	}
	else
	{
		MWBM_SCALE_WEIGHTS(*this, edgeCosts);
		l = MAX_WEIGHT_BIPARTITE_MATCHING(*this, edgeCosts);
	}

	leda::list_item i;
	leda::node u, v;
	double s = 0;

	forall_items(i, l) 
	{
		e = l[i];
		s += inf(e);

		if (pNodeMap)
		{
			if (g1MatchedNodeMap.defined(e))
			{
				u = g1MatchedNodeMap[e];
				v = g2MatchedNodeMap[e];
			}
			else
			{
				u = inf(source(e));
				v = inf(target(e));
			}

			(*pNodeMap)[u] = v;
		}
	}

	if (pnMatched)
		*pnMatched = l.size();

	return s;
}

DAG::MatchParams::MatchParams()
{
	memset(this, 0, sizeof(*this));
}

void DAG::MatchParams::Print(std::ostream& os) const
{
	PRINT_OPEN(os, (int)nSimFuncType);
	PRINT(os, nNodeDistFunc);
	PRINT(os, dTSVSimWeight);
	PRINT(os, dSimilMassWeight);
	PRINT(os, dRelMassWeight);
	PRINT(os, dMinNodeSim);
	PRINT(os, dBreakSibRelPen);
	PRINT(os, nPreserveAncestorRel);
	PRINT(os, nUseNewVoteWeightFunc);
	PRINT_CLOSE(os, nUseMOOVC);
}

/*static*/
NodeIdxArray MatchedNodePair::GetSiblingsVector(int node, const SmartMatrix<int>& tcm, const NodeIdxMatrix& parents)
{
	const NodeIdxArray& vpar = parents[node];
	NodeIdxArray sibs;
	int parent;
		
	// look for all the node tha can be reached from node's parents
	// i.e., the node's siblings and their descendents
	
	for (int i = 0; i < vpar.GetSize(); i++)
	{	
		parent = vpar[i];
		
		if (sibs.GetSize() == 0)
			sibs = tcm[parent];
		else
		{
			const NodeIdxArray& halfSibs = tcm[parent];
			ASSERT(sibs.GetSize() == halfSibs.GetSize());
			
			// Do sibs |= halfSibs;
			for (int j = 0; j < halfSibs.GetSize(); j++)
				if (halfSibs[j] == 1 && sibs[j] != 1) // only modifies sibs if necessary
					sibs[j] = 1;
		}
	
		// we consider that a node is a sibling of itself so that case is okay. But,
		// the node's parent isn't a sibling, so we must correct for this.	
		sibs[parent] = 0;
	}
	
	if (sibs.Sum() <= 1)
		sibs.Clear();
		
	return sibs;
}

/*!
	@brief Updates the similarity matrix with the penalty terms for breaking
	sibling relations and ancestor relations.
	
	Let S_{Qv} the vector of sibling of node v in graph Q, which are computed from
	adjacency matrix and the trans. clos. mat. of the graph.
	
	S_{Qv} = OR(T_Q(A_Q(:,v), :))
	S_{Mw} = OR(T_M(A_M(:,w), :))
	
	The vector of penalties for nodes in M mathing non-siblings of v in Q can be computed
	by: S'_{Qv} = penalty * (1 - S_{Qv}). Similarly for S'_{Mw} = penalty * (1 - S_{Mw}).
	
	Then, we proceed to update the sim matrix by
	
	simMat(S_{Qv}, :) = simMat(S_{Qv}, :) .* S'_{Mw};
	simMat(S_{Mw}, :) = simMat(S_{Mw}, :) .* S'_{Qv};
*/
void MatchedNodePair::UpdateSimMat(SimMatrix& simMat, const double& alpha) const
{
	int rows = simMat.NRows(), cols = simMat.NCols();
	ASSERT(alpha >= 0 && alpha <= 1);
	
	if (alpha == 1) //nothing to do
		return;
	
	SmartArray<int> sq, sm;	
	int i, j;
	
	sq = GetSiblingsVector(q, *pQueryTCM, qparents);
	sm = GetSiblingsVector(m, *pModelTCM, mparents);
	
	rows = sq.GetSize();
	cols = sm.GetSize();
	
	// rows or cols will be zero when there are no siblings
	// if either q or m has no siblings, there is no point in
	// penalizing the siblings that do exist, because there is no
	// sibling relations to preserve.
	if (rows == 0 || cols == 0)
		return;
					
	// Update the similarity matrix
	for (i = 0; i < rows; i++)
		if (sq[i] == 1) 					// if 'v' in Q is a sibling of 'q'
			for (j = 0; j < cols; j++)		// travers its matches from M
				if (sm[j] == 0) 			// if 'w' in M isn't a sibling of 'm'
					simMat[i][j] *= alpha;
	
	for (j = 0; j < cols; j++)
		if (sm[j] == 1) 					// if 'v' in Q is a sibling of 'q'
			for (i = 0; i < rows; i++)		// travers its matches from M
				if (sq[i] == 0) 			// if 'w' in M isn't a sibling of 'm'
					simMat[i][j] *= alpha;
}

/*!
 	@brief Builds a sort of reduced adjacency matrix for the graph.

 	Needs to save the parents before the DAGs are splitted and
 	links on loopy nodes are lost.
 */
void MatchedNodePair::SetParents(const DAG& g, NodeIdxMatrix& parents)
{
 	leda::node v;
 	leda::edge e;
 	int i, j;

 	parents.ReSize(g.GetNodeCount());

 	forall_nodes(v, g)
 	{
 		i = g.GetNodeDFSIndex(v);	
 		ASSERT(parents[i].GetSize() == 0);
 		parents[i].ReSize(g.indeg(v));

 		j = 0;
 		forall_in_edges(e, v)
 			parents[i][j++] = g.GetNodeDFSIndex(g.source(e));
 	}
}
