/**------------------------------------------------------------------------
 * @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
 *
 *-----------------------------------------------------------------------*/

typedef leda_d_array<leda_node, int> NodeIndexMap;


//! Create an edge between vertices <i> and <j> of weight specified by the edge label.
void DAG::AddEdge(Matrix& adj, leda_edge e, int srcNodeIdx, int tgtNodeIdx) const
{
	double w = GetEdgeWeight(e); // w == 1 in our case

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

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

	We are supposed to compute the $nVals = \delta(T_i)$ largest absolute eigenvalues, 
	where $\delta(T_i)$ is the outdegree of the root node.

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

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

	b) by summing $\delta(T_i)$ elements we normalize the sum according to the local
	   complexity of the subgraph root.
*/
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];

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

	return td;
}

/*!
	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));

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

	return s;
}