/**------------------------------------------------------------------------
* @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: Marcos Da Silva
*
*-----------------------------------------------------------------------*/

#include "SGDistMeasurer.h"
#include "ShockGraph.h"
#include "DAG.h"
#include "HelperFunctions.h"
#include "SGNode.h"
#include <float.h>
#include <stdio.h>
#include <ctype.h>
#include <assert.h>

#include "ApproxPoly.h"

#ifdef WIN32
#define M_PI 3.14159265359
#endif

#define MINDIST 0
#define MAXDIST 1

using namespace std;
using namespace dml;
using namespace leda;

//////////////////////////////////////////////////////////////////////////////
// SGDistMeasurer3

void SGDistMeasurer3::LogResults(const ModelFit& m1,
								 const ModelFit& m2, dml::POINTS g1Pts, 
								 dml::POINTS g2Pts, bool bReverseOrder) const
{
	static bool bFirstTime = true;
	static fstream log;
	
	//Log the results if debug mode. In any case, empty the log file
	if (bFirstTime)
	{
		bFirstTime = false;
		log.open("nodedist.m", ios::out | ios::trunc);
	}
	//else
	//	log.open("nodedist.m", ios::out | ios::app);

	int d0, dN;
	
	// If reverted, put them back in place before displaying
   	if (bReverseOrder)
   	{
   		g1Pts = pNode1->GetVelocityRadiusArray(d0, dN, false);
   		g2Pts = pNode2->GetVelocityRadiusArray(d0, dN, false);
   	}
		
   	log << "disp('" << pG1->GetDAGLbl() << "');\n";
   	log << "disp('" << pG2->GetDAGLbl() << "');\n";
   	
   	log << "disp('" << pNode1->GetNodeLbl()
   		<< " vs " << pNode2->GetNodeLbl()
   		<< " [DFS indices: " << pG1->GetNodeDFSIndex(v1) << " vs " << pG2->GetNodeDFSIndex(v2) << "]');\n";
   	log << "% SEARCH CODE: " << pG1->GetViewNumber() << pNode1->GetNodeLbl()
   	 	<< pG2->GetViewNumber() << pNode2->GetNodeLbl() << "\n";
   		
   	m1.Plot(log, g2Pts, pNode2->GetSegments());
   	
   	log << "disp('" << pG1->GetDAGLbl() << "');\n";
   	log << "disp('" << pG2->GetDAGLbl() << "');\n";
		
   	log << "disp('" << pNode2->GetNodeLbl()
   		<< " vs " << pNode1->GetNodeLbl()
   		<< " [DFS indices: " << pG2->GetNodeDFSIndex(v2) << " vs " << pG1->GetNodeDFSIndex(v1) << "]');\n";
   		
   	m2.Plot(log, g1Pts, pNode1->GetSegments());	
		
   	//log << "disp('Avg error: " << error << "');\n";
}

/*!
	Use pG1, pG2, v1, v2, pNode1, pNode2
*/
double SGDistMeasurer3::CompDistance() const
{	
	/*if (pG1->GetViewNumber() == 79 && pG2->GetViewNumber() == 6 &&
		pG1->GetNodeDFSIndex(v1) == 1 && pG2->GetNodeDFSIndex(v2) == 1)
	{
		cout << "\n" << pNode1->GetNodeLbl() << "\n" << pNode2->GetNodeLbl() << "\n";
	}*/
		
	int n1Type = pG1->NodeType(v1);
	int n2Type = pG2->NodeType(v2);
		
	// If label is different, do not compare nodes	
	if (n1Type != n2Type)
		return MAXDIST;	
	else if (n1Type == ROOT || n2Type == ROOT)
		return n1Type == n2Type ? MINDIST:MAXDIST;
	else if (pNode1->m_shocks.GetSize() == 1 && pNode2->m_shocks.GetSize() == 1)
	{
		double r1 = pNode1->m_shocks.GetHead().radius;
		double r2 = pNode2->m_shocks.GetHead().radius;
		
		return 1 - MIN(r1, r2) / MAX(r1, r2);
	}
	else if (pNode1->m_shocks.GetSize() == 1 || pNode2->m_shocks.GetSize() == 1)
	{
		if (pNode1->m_shocks.GetSize() > 4 || pNode2->m_shocks.GetSize() > 4)
			return MAXDIST; // if it is too long, it's unlikely that they are the same
			
		double r1 = pNode1->m_shocks.AvgRadius();
		double r2 = pNode2->m_shocks.AvgRadius();
		double l1 = pNode1->m_shocks.Length();
		double l2 = pNode2->m_shocks.Length();
		
		double rdiff = 1 - MIN(r1, r2) / MAX(r1, r2);
		double lendiff = 1 - MIN(l1, l2) / MAX(l1, l2);
				
		return 0.85 * rdiff + 0.15 * lendiff;
	}
	
	ASSERT(pNode1->GetSegments().GetSize() > 0 && pNode2->GetSegments().GetSize() > 0);
	
	// We are dealing with same label, non-root nodes, with more that 1 shock point
	// so...
	leda_node par1 = pG1->GetFirstParent(v1);
	leda_node par2 = pG2->GetFirstParent(v2);
	int n1Dir, n2Dir;
	
	if (pG1->NodeType(par1) != ROOT && pG2->NodeType(par2) != ROOT)
	{
		n1Dir = pG1->GetBranchDir(v1, par1);
		n2Dir = pG2->GetBranchDir(v2, par2);
	
		// If the nodes' directions wrt their parents differ,
		// we can say they directions don't match
		if (n1Dir != n2Dir)
			return MAXDIST;
	}
	
	// Now, given that their directions match, we need to make
	// sure that we get them right. The may go from 0 to N or form N to 0.
	n1Dir = pNode1->GetBranchDir(ShockBranch::FORWARD);
	n2Dir = pNode2->GetBranchDir(ShockBranch::FORWARD);
	bool bReverseOrder = ( n1Dir != n2Dir );
	
	// what if dir == 0 but there is a slope? We should still reverse it.
	// THIS IS A TO DO
	
	int d0, dN;
	POINTS g1Pts = pNode1->GetVelocityRadiusArray(d0, dN, bReverseOrder);
	POINTS g2Pts = pNode2->GetVelocityRadiusArray(d0, dN, bReverseOrder);
	
	
// 	if(!( pNode2->GetSegments().GetTail().p1.y == (bReverseOrder ? g2Pts.GetHead().y:g2Pts.GetTail().y) ))
// 	{
// 		DBG_VAL(g1Pts.GetSize())
// 		DBG_VAL(g2Pts.GetSize())
// 		DBG_VAL(pNode1->GetShockCount())
// 		DBG_VAL(pNode2->GetShockCount())
// 		DBG_VAL(n1Dir)
// 		DBG_VAL(n2Dir)
// 		DBG_VAL(pNode2->GetSegments().GetTail().p0.y)
// 		DBG_VAL(pNode2->GetSegments().GetTail().p1.y)
// 		DBG_VAL(pNode2->GetSegments().GetHead().p0.y)
// 		DBG_VAL(pNode2->GetSegments().GetHead().p1.y)
// 		DBG_VAL(pNode1->GetSegments().GetTail().p0.y)
// 		DBG_VAL(pNode1->GetSegments().GetTail().p1.y)
// 		DBG_VAL(pNode1->GetSegments().GetHead().p0.y)
// 		DBG_VAL(pNode1->GetSegments().GetHead().p1.y)
// 		DBG_VAL(g2Pts.GetHead().y)
// 		DBG_VAL(g2Pts.GetTail().y)
// 		DBG_VAL(g1Pts.GetHead().y)
// 		DBG_VAL(g1Pts.GetTail().y)
// 		
// 		for (int k = 0; k < pNode2->GetShockCount(); k++)
// 			cerr << (*pNode2)[k].radius << ", ";
// 		cerr << endl;
// 	}
	
	ASSERT( pNode1->GetSegments().GetHead().p0.y == (bReverseOrder ? g1Pts.GetTail().y:g1Pts.GetHead().y) || g1Pts.GetSize() <= 3 );
	ASSERT( pNode1->GetSegments().GetTail().p1.y == (bReverseOrder ? g1Pts.GetHead().y:g1Pts.GetTail().y) || g1Pts.GetSize() <= 3 );
	ASSERT( pNode2->GetSegments().GetHead().p0.y == (bReverseOrder ? g2Pts.GetTail().y:g2Pts.GetHead().y) || g2Pts.GetSize() <= 3 );
	ASSERT( pNode2->GetSegments().GetTail().p1.y == (bReverseOrder ? g2Pts.GetHead().y:g2Pts.GetTail().y) || g2Pts.GetSize() <= 3 );
	
	// Let's fit the model of node one to the data of node 2 and viceversa
	ModelFit m1(2), m2(2); // 2 is max diff in lenght
	
	double e1 = m1.Fit(g1Pts, pNode2->GetSegments());
	
	if (e1 >= INFINITY)
		return MAXDIST;
	
	double e2 = m2.Fit(g2Pts, pNode1->GetSegments());
	
	if (e2 >= INFINITY)
		return MAXDIST;
		
	if (DAG::IsDbgMode())
		LogResults(m1, m2, g1Pts, g2Pts, bReverseOrder);
		
	/*double e, s;
	
	for (int i = 0; i < pNode2->GetSegments().GetSize(); i++)
	{
		m1.GetSegmentError(i, &e, &s);
		DBG_SHOW(e)
		DBG_SHOW(s)
	}*/
	
	// We have non-infinite errors, so we must normalize them
	double avge1 = (e1 / g1Pts.GetSize()) / pNode1->m_shocks.AvgRadius();
	double avge2 = (e2 / g2Pts.GetSize()) / pNode2->m_shocks.AvgRadius();
	double error = (avge1 + avge2) / 2.0;
	
	if (error >= MAXDIST)
		return MAXDIST;
	
	double minL = MIN(pNode1->m_shocks.Length(), pNode2->m_shocks.Length());
	double maxL = MAX(pNode1->m_shocks.Length(), pNode2->m_shocks.Length());
		
	return 0.85 * error + 0.15 * (1 - minL / maxL);
}

//////////////////////////////////////////////////////////////////////////////
// SGDistMeasurer1

double SGDistMeasurer1::CompDistance() const
{		
	double level1 = NodeDistanceLevel1();
	double level2 = 0;

	//	double level2 = NodeDistanceLevel2(v1,v2);
		
	/*if(pG1->GetNodeLbl(v1) == pG2->GetNodeLbl(v2)){
//	cout << pG1->GetNodeLbl(v1) << "<>" << pG2->GetNodeLbl(v2) << endl;
//	cout << level1 << ":" <<level2 << endl;
}*/
	
	//double value = w1 * (level1) + w2 * (level2);
	double value = level1;
	
	//cout << " s = " << value;
	
	ASSERT(value >= 0 && value <= 1);
	
	return value;
}

/* 
Takes two numbers which represent the relative weights 
of level1 vs level2 in the node distance calculation.
*/
void SGDistMeasurer1::SetWeights(double level1Weight, double level2Weight)
{
	ASSERT(level1Weight + level2Weight == 1);
	
	w1 = level1Weight;
	w2 = level2Weight;
}

double SGDistMeasurer1::AverageRadius(const SGNode *sb) const
{
	double min_rad = 999999;
	double max_rad = 0;
	double avg_rad = 0;
	int shockCount;
	
	shockCount = sb->GetShockCount();
	
	for (int j = 0; j < shockCount; j++)
    {
		double curr_rad = sb->m_shocks[j].radius;	    
		if (curr_rad > max_rad)
			max_rad = curr_rad;
		if (curr_rad < min_rad)
			min_rad = curr_rad;
		
		avg_rad += curr_rad;
    }
	
	if (shockCount != 0)
		avg_rad /= shockCount;
	return avg_rad;
}

double SGDistMeasurer1::Length(const SGNode *sb) const
{
	int shockCount = sb->GetShockCount();
	
	double distance = 0;
	double x1,x2,y1,y2;
	
	for(int i = 1; i < shockCount; i++)
    {
		x1 =  sb->m_shocks[i-1].xcoord;
		y1 =  sb->m_shocks[i-1].ycoord;
		
		x2 =  sb->m_shocks[i].xcoord;
		y2 =  sb->m_shocks[i].ycoord;
		
		distance += sqrt(pow(x2-x1,2.0) + pow(y2-y1,2.0));
    }
	
	return distance;
}

double SGDistMeasurer1::Elongation(const SGNode *sb) const
{
	double avg_radius_sb = AverageRadius(sb);
	double length = Length(sb);
	double elongation;
	
	
	if(length != 0)
	  {
	    // This just normalizes the value.
	    elongation = (avg_radius_sb/length)/(avg_radius_sb + length);
	  }
	else
	  elongation = avg_radius_sb;
	//	cout << "elongation = " << elongation << endl;
	return elongation;
		
}

double SGDistMeasurer1::Bend(const SGNode *u) const
{
	double x1,x2,y1,y2;
	int numSinModel = u->GetShockCount();
	
	if(numSinModel <= 3)
		return M_PI;
	
	x1 =  u->m_shocks[numSinModel/2].xcoord - u->m_shocks[0].xcoord;
	x2 =  u->m_shocks[numSinModel-1].xcoord - u->m_shocks[numSinModel/2].xcoord;
	
	y1 =  u->m_shocks[numSinModel/2].ycoord - u->m_shocks[0].ycoord;
	y2 =  u->m_shocks[numSinModel-1].ycoord - u->m_shocks[numSinModel/2].ycoord;
	
	
	double value = (x1*x2+y1*y2)/(sqrt(pow(x1,2.0)+pow(y1,2.0))*sqrt(pow(x2,2.0)+pow(y2,2.0)));
	double theta;
	
	if(fabs(value) <= 1)
		theta = acos(value);
	else
		theta = 0;
	
	return theta;
}


double SGDistMeasurer1::SweepContrast(const SGNode *node) const
{
	
	double min_rad = 999999;
	double max_rad = 0;
	double curr_rad = 0;
	double sweep_contrast = 0.0;
	
	int shockCount = node->GetShockCount();

	if (shockCount <= 0)
		return 0;
	
	double LengthNode = Length(node);
	
	for (int j = 0; j < shockCount; j++)
    {
		double curr_rad = node->m_shocks[j].radius;
	    
		if (curr_rad > max_rad)
			max_rad = curr_rad;

		if (curr_rad < min_rad)
			min_rad = curr_rad;
    }

	if (max_rad > 0.0)
		sweep_contrast = ((double) max_rad - (double) min_rad) / (double) max_rad;

	return(sweep_contrast);
}


double SGDistMeasurer1::Taper(const SGNode *u) const
{
	
	int numSinNode = u->GetShockCount();  
	double length = Length(u);
	double slope;
	double del_radius;

	if (numSinNode < 2)
		return 0.0;
	
	del_radius = fabs(u->m_shocks[0].radius - u->m_shocks[numSinNode - 1].radius);
	
	if (del_radius == 0.0)
	{
		slope = 0.0;
	}
	else
	{
		slope = del_radius / length;
	}
	//	cout << "del_radius = " << del_radius << " length = " << length << " slope = " << slope << endl;
	
	return(slope);
}

double SGDistMeasurer1::Acceleration(const SGNode *u) const
{
	int numShocks = u->GetShockCount();
	double accelerationSum = 0;
	
	for(int i = 2; i < numShocks; i++)
    {
		// These are the X and Y coordinates 
		// in the shocks that will be used 
		// to approximate the derivative.
		double x1 =  u->m_shocks[i-2].xcoord;
		double x2 =  u->m_shocks[i-1].xcoord;
		double x3 =  u->m_shocks[i].xcoord;
		
		double y1 =  u->m_shocks[i-2].ycoord;
		double y2 =  u->m_shocks[i-1].ycoord;
		double y3 =  u->m_shocks[i].ycoord;
		
		double radius1 =  u->m_shocks[i-2].radius;
		double radius2 =  u->m_shocks[i-1].radius;
		double radius3 =  u->m_shocks[i].radius;
		
		double D1 = sqrt(pow(x2-x1,2.0) + pow(y2-y1,2.0));
		double D2 = sqrt(pow(x3-x2,2.0) + pow(y3-y2,2.0));
		
		double dT1 = radius2 - radius1;
		double dT2 = radius3 - radius2;
		
		double V1 = D1/dT1;
		double V2 = D2/dT2;
		
		double A = (V2 - V1)/dT2 - dT2;
		
		accelerationSum += A;
    }  
	return accelerationSum;
}

double SGDistMeasurer1::ElongationDiff() const
{
	double u_elongation = Elongation(pNode1);
	double v_elongation = Elongation(pNode2);
	
	double elongation_diff;
	
	if(fabs(v_elongation + u_elongation) != 0)
		elongation_diff = fabs(v_elongation - u_elongation)/fabs(v_elongation + u_elongation);
	else{
		elongation_diff =  fabs(v_elongation - u_elongation);
	}
	
	return elongation_diff;
}

double SGDistMeasurer1::TaperDiff() const
{
	double taper_u = Taper(pNode1);
	double taper_v = Taper(pNode2);
	
	double taper_max, taper_diff;
	
	if (taper_u >= taper_v) 
	{
		taper_max = taper_u;
	}
	else
	{
		taper_max = taper_v;
	}
	if (taper_max == 0)
	{
		taper_diff = 0;
	}
	else
	{
		taper_diff = fabs(taper_u - taper_v)/taper_max;
	}
	return(taper_diff);
}

double SGDistMeasurer1::SweepContrastDiff() const
{
	double u_sweep_contrast = SweepContrast(pNode1);
	double v_sweep_contrast = SweepContrast(pNode2);
	double max_sweep_contrast = 0.0;
	
	
	double sweep_contrast_diff = 0.0;
	
	if (u_sweep_contrast >= v_sweep_contrast)
	{
		max_sweep_contrast = u_sweep_contrast;
	}
	else
	{
		max_sweep_contrast = v_sweep_contrast;
	}
	
	if((u_sweep_contrast == 0.0) && (v_sweep_contrast == 0.0))
	{
		sweep_contrast_diff = 0;
	}
	else
	{
		sweep_contrast_diff = fabs(u_sweep_contrast - v_sweep_contrast)/max_sweep_contrast;
	}
	
	return (sweep_contrast_diff);
}

/* 
This takes two nodes and returns a number representing
the distance between these two nodes in terms of the
2nd derivative of the radius vs time "graph".
*/

double SGDistMeasurer1::AccelerationDiff() const
{
	double u_acc = Acceleration(pNode1);
	double v_acc = Acceleration(pNode2);
	
	double acc_diff;
	
	if(fabs(v_acc) + fabs(u_acc) != 0)
		acc_diff = fabs(v_acc -u_acc)/(fabs(v_acc)+fabs(u_acc));
	else
		acc_diff =  fabs(v_acc - u_acc);
	
	return acc_diff;
}

/* 
This takes two nodes and returns a number representing
the distance between these two nodes in terms of an
approximation of an angle that represents it's medial
axis.
*/

double SGDistMeasurer1::BendDiff() const
{
	double u_bend = Bend(pNode1);
	double v_bend = Bend(pNode2);
	
	double bend_diff = fabs(u_bend - v_bend)/M_PI;
	
	return bend_diff;
}

/* 
This function is the control function for the first level of
matching. 
*/
double SGDistMeasurer1::NodeDistanceLevel1() const
{ 
	int uType = pG1->NodeType(v1);
	int vType = pG2->NodeType(v2);
	
	if (uType == ROOT && vType == ROOT)
		return 0;
	else if(uType == ROOT || vType == ROOT)
		return 1;
		
	//cout << endl << pG1->GetNodeLbl(v1) << " <> " << pG2->GetNodeLbl(v2);
	
	if(uType == SHOCK_1 && vType == SHOCK_1)
    {
		if (pNode1->GetShockCount() <= 1 || pNode2->GetShockCount() <= 1)
		{
			DBG_MSG("Warning: Type 1 node with 1 or 0 shocks!")
			return 1;
		}

		double taper = TaperDiff();
		//		double acc = AccelerationDiff();
		double bend = BendDiff();
		double elongation = ElongationDiff();
		
		//		cout << "taper diff = " << taper << " bend diff = " << bend << " elongation diff = " << elongation << endl;
		
		
		double value = 0.5*taper + 0.2*bend + 0.3*elongation; // + 0.1*acc;
		
		if(value <= 1 && value >= 0)
			return value;
		else if(value < 0)
		{
			cerr << "Error in the distance calculation type 1 SHOCK!!" << value << endl;
			assert(false);
			return -1;
		}
		else
		{
			cerr << "Error in the normalization type 1 SHOCK!!" << value << endl;
			assert(false);
			return -1;
		}
    }
	else if(uType == SHOCK_3 && vType == SHOCK_3)
    {
		if (pNode1->GetShockCount() <= 1 || pNode2->GetShockCount() <= 1)
		{
			DBG_MSG("Warning: Type 3 node with 1 or 0 shocks!")
			return 1;
		}

		double bend = BendDiff();
		double elongation = ElongationDiff();
		
		if(bend <= 1 && bend >= 0)
		{
			double value = 0.8*bend + 0.2*elongation;
			
			if(value >= 0 && value <= 1)
				return value;
			else
			{
				cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
				assert(false);
				return -1;
			}
		}
		else if(bend < 0 || elongation < 0)
		{
			cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
			assert(false);
			return -1;
		}
		else
		{
			cerr << "Error in the normalization type 3 SHOCK!!" << bend << endl;
			assert(false);
			return -11;
		}
    }
	else if(uType == SHOCK_2 && vType == SHOCK_2)
    {
		if (pNode1->GetShockCount() <= 0 || pNode2->GetShockCount() <= 0)
		{
			DBG_MSG("Warning: Type 2 node with 0 shocks!")
			return 1;
		}

		double bend = BendDiff();
		double elongation = ElongationDiff();
		double contrast_sweep = SweepContrastDiff();
		
		//		cout << "cont sweep diff = " << contrast_sweep << " bend diff = " << bend << " elongation diff = " << elongation << endl;
		
		if(bend <= 1 && bend >= 0)
		{
			double value = 0.3 * bend + 0.2 * elongation + 0.5 * contrast_sweep;
			
			if(value >= 0 && value <= 1)
				return value;
			else
			{
				cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
				assert(false);
				return -1;
			}
		}
		else if(bend < 0 || elongation < 0)
		{
			cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
			assert(false);
			return -1;
		}
		else{
			cerr << "Error in the normalization type 2 SHOCK!!" << bend << endl;
			assert(false);
			return -1;
		}
    }
	else if(uType == SHOCK_4 && vType == SHOCK_4)  //updated to compare non-perfect blobs, using same criteria as type 2.
    {
		double bend = BendDiff();
		double elongation = ElongationDiff();
		
		if(bend <= 1 && bend >= 0)
		{
			double value = 0.8*bend + 0.2*elongation;
			
			if(value >= 0 && value <= 1)
				return value;
			else
			{
				cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
				assert(false);
				return -1;
			}
		}
		else if(bend < 0 || elongation < 0)
		{
			cerr << "Error in the distance calculation type 2 SHOCK!!" << bend << endl;
			assert(false);
			return -1;
		}
		else{
			cerr << "Error in the normalization type 2 SHOCK!!" << bend << endl;
			assert(false);
			return -1;
		}
    }
	else if((vType == SHOCK_1 && uType == SHOCK_3) || (vType == SHOCK_3 && uType == SHOCK_1))
    {    
		double taper = TaperDiff();
		double bend = BendDiff();
		double elongation = ElongationDiff();
		
		double value = 0.5*taper + 0.3*elongation +  0.2*bend;
		
		if(value <= 1)
			return value;
		else if(value < 0){
			cerr << "Error in comparison of TYPE 3 and 1!! " << value << endl;
			assert(false);
			return -1;
		}
		else
		{
			cerr << "Error in normalization of TYPE 3 and 1!! " << value << endl;
			assert(false);
			return -1;
		}
    }
	/*	else if((uType == SHOCK_3 && vType == SHOCK_2) || (uType == SHOCK_2 && vType == SHOCK_3))
    {
	double bend = BendDiff();
	double elongation = ElongationDiff();
	
	 if(bend <= 1 && bend >= 0)
	 {
	 double value = 0.8*bend + 0.2*elongation;
	 
	  if(value >= 0 && value <= 1)
	  return value;
	  else
	  {
	  cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
	  assert(false);
	  return -1;
	  }
	  }
	  else if(bend < 0 || elongation < 0)
	  {
	  cerr << "Error in the distance calculation type 3 SHOCK!!" << bend << endl;
	  assert(false);
	  return -1;
	  }
	  else
	  {
	  cerr << "Error in the normalization type 3 SHOCK!!" << bend << endl;
	  assert(false);
	  return 1;
	  }
	  }
	*/
	else{
		
		return 1;
	}
}

/* 
This function implements rule 1 of the grammar. It tries to quantify how
close the roots (a and b) are to each other in terms of the sub tree they
generated. 
*/

double SGDistMeasurer1::rule1(leda_node a, leda_node b) const
{
	int aType = pG1->NodeType(a);
	int bType = pG2->NodeType(b);
	
	double num_3A = 0, num_3B = 0, num_4A = 0,num_4B = 0;
	
	leda_node w,s;
	
	if(aType == ROOT &&  bType == ROOT)
    {
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_3)
				num_3A++;
			else if(pG1->NodeType(w) == SHOCK_4)
				num_4A++;
		}
		
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_3)
				num_3B++;
			else if(pG2->NodeType(s) == SHOCK_4)
				num_4B++;
		}
		
		if(num_3A + num_4A == 0){
			cerr << "Tree has no children of the root" << endl;
			assert(false);
			return -1;
		}
		
		if(num_3B + num_4B == 0){
			cerr << "Tree has no children of the root" << endl;
			assert(false);
			return -1;
		}
		
		double score;
		
		if(num_3A + num_3B > 0)
			score = fabs(num_3A - num_3B)/(num_3A + num_3B);
		else if(num_3A + num_3B == 0)
		{
			score = 0;
		}
		else
		{
			cerr << "Tree has a negative number of children of the root node" << endl;
			assert(false);
			return -1;
		}
		
		if(num_4A + num_4B > 0)
			score += fabs(num_4A - num_4B)/(num_4A + num_4B);
		else if(num_4A + num_4B == 0)
		{
			score += 0;
		}
		else
		{
			cerr << "Tree has a negative number of children of the root node" << endl;
			assert(false);
			return -1;
		}
		
		return score/2;
    }
	else
    {
		return 1;
    }
}

/* 
This function tries to check how closely related both nodes
that were generated from rule 2 match. 
*/

double SGDistMeasurer1::rule2(leda_node a, leda_node b)  const
{
	leda_node w,s;
	int aType = pG1->NodeType(a);
	int bType = pG2->NodeType(b);
	
	int num_child_a = 0, num_child_b = 0;
	
	if(aType == SHOCK_4 && bType == SHOCK_4)
    {
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_1)
				num_child_a++;
		}
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_1)
				num_child_b++;
		}
		
		if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
			return 0;
		else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
			return 1;
    }
	else
    {
		return 1;
    }
	
	double sim = abs(num_child_a - num_child_b);
	
	
	// If this isn't true then it must go into death
	// and it's outdeg should be zero!! I check for 
	// this above.
	assert(num_child_a + num_child_b > 0);
	
	// normalizing
	sim = sim/(num_child_a +num_child_b);
	
	return sim;
}

/* 
This function tries to check how closely related both nodes
that were generated from rule 3 match. NOTE this function
assumes that the nodes being passed as parameters are of
type 3.
*/

double SGDistMeasurer1::rule3(leda_node a, leda_node b)  const
{
	
	leda_edge e = pG1->first_in_edge(a);
	leda_edge f = pG2->first_in_edge(b);
	
	int aParentType, bParentType;
	int num_child_a = 0, num_child_b = 0;
	
	int aType = pG1->NodeType(a);
	int bType = pG2->NodeType(b);  
	
	
	// If a node has been removed such that it
	//  no longer has any parent then we will
	// assume that it's parent is the ROOT.
	if(e == nil && f == nil)
    {
		aParentType = ROOT;
		bParentType = ROOT;   
    }
	else if (e == nil || f == nil)
    {
		return 1;
    }
	else
    {
		aParentType = pG1->NodeType(pG1->source(e));
		bParentType = pG2->NodeType(pG2->source(f));
    }
	
	leda_node w,s;
	
	// rule 3a
	if(aParentType == SHOCK_1 && bParentType == SHOCK_1 
		&& pG1->indeg(a) == 1 && pG2->indeg(b) == 1)
    {
		
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_1)
				num_child_a++;
		}
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_1)
				num_child_b++;
		}
		
		// If both go to death...
		if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
			return 0;
		else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
			return 1;
    }
	else if(aParentType == ROOT && bParentType == ROOT)
    {
		// Rule 3 b
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_1)
				num_child_a++;
		}
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_1)
				num_child_b++;
		}
		
		if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
			return 0;
		if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
			return 1;
    }
	else
    {
		return 1;
    }
	
	double sim = abs(num_child_a - num_child_b);
	
	// If this isn't true then it must go into death
	// and it's outdeg should be zero!! I check for 
	// this above.
	assert(num_child_a + num_child_b > 0);
	
	// normalizing
	sim = sim/(num_child_a + num_child_b);
	
	return sim;
}

/* 
This function tries to check how closely related both nodes
that were generated from rule 4 match. 
*/

double SGDistMeasurer1::rule4(leda_node a, leda_node b)  const
{
	int aType = pG1->NodeType(a);
	int bType = pG2->NodeType(b);
	int num_child_a = 0, num_child_b = 0;
	leda_node w,s;
	
	if(aType == SHOCK_1 && bType == SHOCK_1)
    {
		// If they go to death
		if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
			return 0;
		else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
			return 1;
		
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_1)
				num_child_a++;
		}
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_1)
				num_child_b++;
		}
    }
	else
    {
		return 1;
    }
	
	double sim = abs(num_child_a - num_child_b);
	
	// If this isn't true then it must go into death
	// and it's outdeg should be zero!! I check for 
	// this above.
	assert(num_child_a + num_child_b > 0);
	
	
	sim = sim/(num_child_a +num_child_b);
	
	return sim;
}


/* 
This function tries to check how closely related both nodes
that were generated from rule 5 match. 
*/

double SGDistMeasurer1::rule5_6(leda_node a, leda_node b)  const
{
	
	int aType = pG1->NodeType(a);
	int bType = pG2->NodeType(b);
	int num_2child_a = 0; 
	int num_2child_b = 0;
	int num_3child_a = 0; 
	int num_3child_b = 0;
	leda_node w,s;
	
	if(aType == SHOCK_1 && bType == SHOCK_1)
    {
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_2)
				num_2child_a++;
			else if (pG1->NodeType(w) == SHOCK_3)
				num_3child_a++;
		}
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_2)
				num_2child_b++;
			else if (pG2->NodeType(s) == SHOCK_3)
				num_3child_b++;
		}    
		
		/* This checks that a type 1 doesn't have more than
		1 type 2 child. */
		assert(num_2child_a <= 1 && num_2child_b <= 1);
		
		/* This checks that a type 1 doesn't have more than
		1 type 3 child. */
		assert(num_3child_a <= 1 && num_3child_b <= 1);
		
		if (pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
			return 0;
		else if (pG1->outdeg(a) == 0 || pG2->outdeg(b) == 0)
			return 1;
		
		// rule 5
		if(num_2child_a == 1 && num_2child_b == 1){
			return 0;
		}
		else if(num_3child_a == 1 && num_3child_b == 1)
		{
			// rule 6a
			leda_node t,r;
			
			forall_adj_nodes(t, a)
			{
				forall_adj_nodes(r, b)
				{
					if(pG1->indeg(t) == 2 && pG2->indeg(r) == 2)
						return 0;
					else if(pG1->indeg(t) == 1 && pG2->indeg(r) == 1)
						return 0;
					else{
						return 1;
					}
				}
			}
		}
		else{
			// Nodes have different types of children
			return 1;
		}
    }
	else
	{
		return 1;
	}

	return 0;
}

/* 
This function is the control function for the second level of
matching. 
*/
double SGDistMeasurer1::NodeDistanceLevel2(leda_node a, leda_node b) const
{
	
	int aType = pG1->NodeType(a);
	int bType = pG2->NodeType(b);
	
	if (aType == ROOT && bType == ROOT)
		return rule1(a,b);
	else if(aType == ROOT || bType == ROOT)
		return 1;
	
	if(pG1->NodeRole(a) == DUMMY && pG2->NodeRole(b) == DUMMY)
    {
		return 0;
    }
	else if(pG1->NodeRole(a) == DUMMY || pG2->NodeRole(b) == DUMMY)
    {
		return 1;
    }
	
	int aParentType = pG1->NodeType(pG1->source(pG1->first_in_edge(a)));
	int bParentType = pG2->NodeType(pG2->source(pG2->first_in_edge(b)));
	
	leda_node parent_a = pG1->source(pG1->first_in_edge(a));
	leda_node parent_b = pG2->source(pG2->first_in_edge(b));
	
	if(aType == SHOCK_1 && bType == SHOCK_1)
    {
		double parent_score;
		
		// IF BOTH EMERGED FROM RULE 2
		if(aParentType == SHOCK_4 && bParentType == SHOCK_4)
		{
			// This will return a normalized number which
			// is the difference in siblings for each of
			// a and b
			
			parent_score = rule2(parent_a, parent_b);
		}
		// IF BOTH EMERGED FROM RULE 3
		else if(aParentType == SHOCK_3 && bParentType == SHOCK_3)
		{
			parent_score = rule3(parent_a, parent_b);
		}
		// IF BOTH EMERGED FROM RULE 4
		else if(aParentType == SHOCK_1 && bParentType == SHOCK_1)
		{
			parent_score = rule4(parent_a, parent_b);
		}
		else
		{
			// Nodes didn't emerge from the same production rule.
			parent_score = 1;
		}
		
		
		
		leda_node w,s;
		bool rule4a = false, rule4b = false;
		
		forall_adj_nodes(w, a)
		{
			if (pG1->NodeType(w) == SHOCK_1){
				rule4a = true;
				break;
			}
		}
		
		forall_adj_nodes(s, b)
		{
			if (pG2->NodeType(s) == SHOCK_1){
				rule4b = true;
				break;
			} 
		}    
		
		// DESCENDANTS SIMILARITY
		
		double child_score;  
		
		if(rule4a && rule4b)
		{
			child_score = rule4(a, b);
		}
		else if(!rule4a && !rule4b)
		{
			child_score = rule5_6(a, b);
		}
		else if(pG1->indeg(a) == 0 && pG2->indeg(b) == 0){
			// If both go to death
			child_score = 0;
		}
		else
		{
			child_score = 1;
		}
		return (parent_score + child_score)/2;
    }
	else if(aType == SHOCK_2 && bType == SHOCK_2)
    {
		// ANCESTOR SIMILARITY
		if(aParentType == SHOCK_1 && bParentType == SHOCK_1 && pG1->outdeg(a) == 0 && pG2->outdeg(b) == 0)
			return 0;
		else
		{
			cerr << "A type 2 shock doesn't go to death or have type 1 parents!!" << endl;
			assert(false);
			return -1;
		}
    }
	else if(aType == SHOCK_3 && bType == SHOCK_3)
    {
		double child_score = rule3(a, b);    
		double ancestor_score;
		
		if(aParentType == ROOT && bParentType == ROOT)
			ancestor_score = 0;
		else if(pG1->indeg(a) == 2 && pG2->indeg(b) == 2)
			ancestor_score = 0;
		else if(pG1->indeg(a) == 1 && pG2->indeg(b) == 1)
		{
			ancestor_score = 0;
		}
		else
			ancestor_score = 1;
		
		return (child_score + ancestor_score)/2;
    }
	else if(aType == SHOCK_4 && bType == SHOCK_4)
    {
		// Descendants
		double child_score = rule2(a, b);
		// parents
		double parent_score = rule1(parent_a,parent_b);
		
		return (child_score + parent_score)/2;
    }
	else{
		return 1;
	}
}

//////////////////////////////////////////////////////////////////////////////
// SGDistMeasurer2

#define MAX_DIST       100
#define DIFF_ROLE_DIST 100
#define DIFF_TYPE_DIST 1

/*void DbgNodeSimilarity(const ShockGraph& sg1, leda_node u, const ShockGraph& sg2, leda_node v,
double scale, double dr, double dv, 
double branchSim, double tsvSim, double nodeSim)
{
static bool bFirstTime = true;
static fstream dbg;

 if (bFirstTime)
 {
 bFirstTime = false;
 dbg.open("nodesim.log", ios::trunc | ios::out);
 }
 
  dbg << sg1.GetDAGLbl() << " vs " << sg2.GetDAGLbl() << endl
  << sg1.GetNodeLbl(u) << " vs " << sg2.GetNodeLbl(v) << endl
  << "role " << sg1.NodeRole(u) << " vs role " << sg2.NodeRole(v) << endl
  << "scale = " << scale << ", dr = " << dr << ", dv = " << dv << endl
  << "Branch similarity = " << branchSim << endl
  << "TSV similarity = " << tsvSim << endl
  << "Node similarity = " << nodeSim << endl
  << endl << flush;
}*/

/*bool SGDistMeasurer1::AreNodesRelated() const
{
const ShockGraph& sgFrom = dynamic_cast<const ShockGraph&>(from);

 NODE_ROLE ru = NodeRole(u), rv = sgFrom.NodeRole(v);
 
  if (ru == UNK_ROLE || rv == UNK_ROLE)
  return NodeType(u) == sgFrom.NodeType(v);		
  else if (ru == DUMMY || rv == DUMMY)
  return false;
  else
  return ru == rv;
}*/

double SGDistMeasurer2::CompDistance() const
{
	ASSERT(false);
	return 0;
}

/*double SGDistMeasurer2::CompDistance()
{
double da = 0, dv = 0, scale = 0, branchSim;
const ShockBranch& b1 = pG1->GetSGNode(v1)->m_shocks;
const ShockBranch& b2 = pG2->GetSGNode(v2)->m_shocks;

 int uType = pG1->NodeType(v1);
 
  if (pG1->NodeRole(v1) != BIRTH)
  {
  if (uType == SHOCK_4)
  {
  branchSim = 0;
  }
  else
  {
  //scale = ComputeScale(u, sgFrom, v);
  
	  da = b1.CompareAreas(b2, 1); //use scale!!!
	  dv = b1.CompareVelocities(b2, 1); //use scale!!!
	  
	   const double w = .7;
	   branchSim = w * da + (1 - w) * dv;
	   }
	   }
	   else
	   {
	   double n1, n2;
	   
		n1 = pG1->NeighborsCount(v1, SHOCK_3);
		n2 = pG2->NeighborsCount(v2, SHOCK_3);
		
		 branchSim = SIMILARITY(n1, n2, 0.5);
		 
		  n1 = pG1->NeighborsCount(v1, SHOCK_4);
		  n2 = pG2->NeighborsCount(v2, SHOCK_4);
		  
		   branchSim *= SIMILARITY(n1, n2, 0.5);
		   }
		   
			double tsvSim = pG1->NodeTSVSimilarity(v1, *pG2, v2);
			double nodeSim = (1 - s_dTSVSimWeight) * branchSim + s_dTSVSimWeight * tsvSim;
			
			 return nodeSim;
			 }
			 
			  int SGDistMeasurer2::NeighborsCount(leda_node v, int nNodeType) const
			  {
			  int nCount = 0;
			  leda_edge e;
			  leda_node w;
			  
			   forall_in_edges(e, v)
			   if (pG1->NodeType(pG1->source(e)) == nNodeType)
			   nCount++;
			   
				forall_adj_nodes(w, v)
				if (NodeType(w) == nNodeType)
				nCount++;
				
				 return nCount;
}*/

// Computes the difference in scale between two branches by using their parents info.
/*double SGDistMeasurer2::ComputeScale(leda_node u, const ShockGraph& from, leda_node v) const
{
double r1, r2;
leda_node uparent = source(first_in_edge(u));
leda_node vparent = from.source(from.first_in_edge(v));

 // Root node shouldn't be analysed!
 ASSERT(uparent != nil && vparent != nil);
 
  if (NodeType(uparent) == ROOT)
  r1 = GetSGNode(u)->m_shocks.AvgRadius();
  else
  r1 = GetSGNode(uparent)->m_shocks.AvgRadius();
  
	  if (from.NodeType(vparent) == ROOT)
	  r2 = from.GetSGNode(v)->m_shocks.AvgRadius();
	  else
	  r2 = from.GetSGNode(vparent)->m_shocks.AvgRadius();
	  
	   return  (r1 == 0 || r2 == 0) ? 1:r1 / r2;
}*/





