---

SGDistMeasurer.cpp

Go to the documentation of this file.

#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

// SGDistMeasurer3

void SGDistMeasurer3::LogResults(const ModelFit& m1,
        const ModelFit& m2, POINTS g1Pts, 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";
}

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);
                
        ASSERT( pNode1->GetSegments().GetHead().p0.y == (bReverseOrder ? g1Pts.GetTail().y:g1Pts.GetHead().y) );
        ASSERT( pNode1->GetSegments().GetTail().p1.y == (bReverseOrder ? g1Pts.GetHead().y:g1Pts.GetTail().y) );
        ASSERT( pNode2->GetSegments().GetHead().p0.y == (bReverseOrder ? g2Pts.GetTail().y:g2Pts.GetHead().y) );
        ASSERT( pNode2->GetSegments().GetTail().p1.y == (bReverseOrder ? g2Pts.GetHead().y:g2Pts.GetTail().y) );
        
        // 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;
        }
}

/* 
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;
}*/

---