---

DAG.cpp

Go to the documentation of this file.

#include <iostream.h>
#include "DAG.h"
#include "HelperFunctions.h"
#include "Exceptions.h"
#include "Debug.h"
#include <LEDA/graphwin.h>
#include <float.h>

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

#define ROOT_NODE_LEVEL                 0
#define ROOT_NODE_DFS_INDEX             0

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

LogFile logFile;

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);

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

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)]);
}

void DAG::ComputeDerivedValues()
{
        leda_node v = GetFirstRootNode();

        ResetNodesInfo(); // Make sure the node info is clean
        ComputeNodesInfo(v, ROOT_NODE_LEVEL, ROOT_NODE_DFS_INDEX);
        ComputeTransiveClosure();
        ComputeNodesMass();
        ComputeTSVs(v);

        nMaxBFactor = 0;
        dTotalTSVSum = 0;
        //nCumulativeMass = 0;
        nCumulativeMass = GetNodeMass(v);

        forall_nodes(v, *this)
        {
                GetNode(v)->ComputeDerivedValues();

                dTotalTSVSum += GetNodeTSVNorm(v);
                //nCumulativeMass += GetNodeMass(v);

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

DAGPtr DAG::SplitSubGraph(leda_node v, bool bRecomputeTSVs)
{
        DAGPtr subg = CreateObject();

        SplitSubGraph(v, subg);

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

        if (bRecomputeTSVs)
                ComputeTSVs(GetFirstRootNode());

        return subg;
}

leda_node DAG::SplitSubGraph(leda_node v, DAGPtr& ptrDag)
{
        leda_node a, b;
        leda_edge e;
        double w;

        a = ptrDag->NewNode(inf(v));

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

        del_node(v);

        return a;
}

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

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));
}

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);
}

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];
}

leda_node DAG::GetFirstRootNode() const
{
        leda_node v;

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

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

/*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;
        double w;
        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();*/

        // Note from the LEDA manual about the 'arithmetic demand' of MWMCB_MATCHING:
        // "This scaling function should be used for the algorithm that compute minimum of 
        // maximum weight assignments or maximum weighted matchings of maximum cardinality." .
        // MWA_SCALE_WEIGHTS(G, edgeCosts); // if we use MWMCB_MATCHING
        MWBM_SCALE_WEIGHTS(G, edgeCosts); // if we use MAX_WEIGHT_BIPARTITE_MATCHING

        // Solve the max weight max card bipartite matching
        //leda_list<leda_edge> tmch = MWMCB_MATCHING(G, A, B, edgeCosts);
        leda_list<leda_edge> tmch = MAX_WEIGHT_BIPARTITE_MATCHING(G, edgeCosts);

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

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

        // 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);
  
        int idxG1Node = g1.GetNodeDFSIndex(maxG1Node);
        int idxG2Node = g2.GetNodeDFSIndex(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 and on the remainder graph
        pMatchedPair->SetEmpty();
        max_weight += Match(*sg1, *sg2, simMat, paired1, paired2, pMatchedPair);
        
        pMatchedPair->SetNodes(idxG1Node, idxG2Node);
        pMatchedPair->UpdateSimMat(simMat, 1 - s_matchParams.dBreakSibRelPen);
        max_weight += Match(g1, g2, simMat, paired1, paired2, pMatchedPair);

        return max_weight;
}

/*static*/
double DAG::Similarity(const DAG& query, const DAG& model, 
                DAGNodeMap& nodeMap, bool bRecomputeTSVs)
{
        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", ios::out | 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;
}

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

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

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

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())
                {
                        cerr << "\nWrong similarity value (" << dNodeSim << ") when matching: \n"
                                << GetDAGLbl() << ":" << GetNodeLbl(u) << " against "
                                << from.GetDAGLbl() << ":" << from.GetNodeLbl(v) << 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;
}


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));
}


void DAG::Print(ostream& os /*=cout*/) const
{
        os << "\nDAG label: " << GetDAGLbl() << endl;
        os << "Object: " << GetObjName() << ", view: " << GetViewNumber() << endl;

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

        leda_node v;

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

void DAG::PrintTSVs(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";
}

istream& DAG::Read(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;
}

ostream& DAG::Write(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;
}

void DAG::DelEdge(Matrix& adj, int srcNodeIdx, int tgtNodeIdx) const
{
        adj(srcNodeIdx, tgtNodeIdx) = 0.0;
        adj(tgtNodeIdx, srcNodeIdx) = 0.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;
}

void DAG::DeleteSubDAG(leda_node v)
{
        leda_node u;

        forall_adj_nodes(u, v)
                DeleteSubDAG(u);

        del_node(v);
}

void DAG::PrintAdjMatrix(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 << setw(1) << setprecision(0) << adj << setprecision(6) << 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;
        }

        //MWA_SCALE_WEIGHTS(*this, edgeCosts);
        MWBM_SCALE_WEIGHTS(*this, edgeCosts);

        EdgeList l;

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

                ASSERT(bBipartite);
                l = MWMCB_MATCHING(*this, setA, setB, edgeCosts);
        }
        else
                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(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;
}

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

---