Main Page Namespace List Class Hierarchy Compound List File List Namespace Members Compound Members File Members Related Pages

# DAG.cpp

Go to the documentation of this file.

    00001
    00052 #include <iostream.h>
    00053 #include "DAG.h"
    00054 #include "HelperFunctions.h"
    00055 #include "Exceptions.h"
    00056 #include "Debug.h"
    00057 #include <LEDA/graphwin.h>
    00058 #include <float.h>
    00059
    00060 #ifdef WIN32
    00061 #define isnan  _isnan
    00062 #define finite _finite
    00063 #endif
    00064
    00065 #define ROOT_NODE_LEVEL                 0
    00066 #define ROOT_NODE_DFS_INDEX             0
    00067
    00068 DAG::MatchParams DAG::s_matchParams;
    00069 bool DAG::s_bDbgMode = false;
    00070
    00071 LogFile logFile;
    00072
    00080 int DAG::ComputeNodesInfo(leda_node v, int nLevel, int nDFSIndex)
    00081 {
    00082         leda_node u;
    00083         DAGNodePtr& ptrNode = GetNode(v);
    00084
    00085         // Default level is -1. Loopy nodes should have their greatest height
    00086         if (ptrNode->GetLevel() < nLevel)
    00087                 ptrNode->SetLevel(nLevel);
    00088
    00089         if (ptrNode->GetDFSIndex() == -1)
    00090                 ptrNode->SetDFSIndex(nDFSIndex++);
    00091
    00092         nLevel++;
    00093
    00094         forall_adj_nodes(u, v)
    00095                 nDFSIndex = ComputeNodesInfo(u, nLevel, nDFSIndex);
    00096
    00097         return nDFSIndex;
    00098 }
    00099
    00100 void DAG::ResetNodesInfo()
    00101 {
    00102         leda_node v;
    00103
    00104         forall_nodes(v, *this)
    00105         {
    00106                 GetNode(v)->SetLevel(-1);
    00107                 GetNode(v)->SetDFSIndex(-1);
    00108         }
    00109 }
    00110
    00111 void DAG::ComputeTransiveClosure()
    00112 {
    00113         leda_node u, v;
    00114         int i, n = GetNodeCount();
    00115
    00116         transClosMat.ReSize(n, n, true);
    00117
    00118         LEDA_GRAPH<leda_node,leda_edge> g = TRANSITIVE_CLOSURE(*this);
    00119
    00120         // Translate the transitive closure graph into a matrix form
    00121         forall_nodes(v, g)
    00122         {
    00123                 i = GetNodeDFSIndex(g[v]);
    00124
    00125                 forall_adj_nodes(u, v)
    00126                         transClosMat[i][GetNodeDFSIndex(g[u])] = 1;
    00127         }
    00128 }
    00129
    00130 void DAG::ComputeNodesMass()
    00131 {
    00132         leda_node v;
    00133         ASSERT(GetNodeCount() == transClosMat.GetSize());
    00134
    00135         /*forall_nodes(v, *this)
    00136                 GetNode(v)->SetMass(transClosMat.RowSum(GetNodeDFSIndex(v)));*/
    00137
    00138         SmartArray<int> ones(transClosMat.NRows());
    00139         ones.Set(1);
    00140
    00141         SmartArray<int> relMasses = transClosMat * (transClosMat * ones);
    00142
    00143         forall_nodes(v, *this)
    00144                 GetNode(v)->SetMass(relMasses[GetNodeDFSIndex(v)]);
    00145 }
    00146
    00154 void DAG::ComputeDerivedValues()
    00155 {
    00156         leda_node v = GetFirstRootNode();
    00157
    00158         ResetNodesInfo(); // Make sure the node info is clean
    00159         ComputeNodesInfo(v, ROOT_NODE_LEVEL, ROOT_NODE_DFS_INDEX);
    00160         ComputeTransiveClosure();
    00161         ComputeNodesMass();
    00162         ComputeTSVs(v);
    00163
    00164         nMaxBFactor = 0;
    00165         dTotalTSVSum = 0;
    00166         //nCumulativeMass = 0;
    00167         nCumulativeMass = GetNodeMass(v);
    00168
    00169         forall_nodes(v, *this)
    00170         {
    00171                 GetNode(v)->ComputeDerivedValues();
    00172
    00173                 dTotalTSVSum += GetNodeTSVNorm(v);
    00174                 //nCumulativeMass += GetNodeMass(v);
    00175
    00176                 if (nMaxBFactor < outdeg(v))
    00177                         nMaxBFactor = outdeg(v);
    00178         }
    00179 }
    00180
    00192 DAGPtr DAG::SplitSubGraph(leda_node v, bool bRecomputeTSVs)
    00193 {
    00194         DAGPtr subg = CreateObject();
    00195
    00196         SplitSubGraph(v, subg);
    00197
    00198         subg->transClosMat = GetTransClosMat();
    00199         subg->nCumulativeMass = nCumulativeMass;
    00200
    00201         if (bRecomputeTSVs)
    00202                 ComputeTSVs(GetFirstRootNode());
    00203
    00204         return subg;
    00205 }
    00206
    00218 leda_node DAG::SplitSubGraph(leda_node v, DAGPtr& ptrDag)
    00219 {
    00220         leda_node a, b;
    00221         leda_edge e;
    00222         double w;
    00223
    00224         a = ptrDag->NewNode(inf(v));
    00225
    00226         forall_adj_edges(e, v)
    00227         {
    00228                 w = inf(e);
    00229                 b = SplitSubGraph(target(e), ptrDag);
    00230                 ptrDag->NewEdge(a, b, w);
    00231         }
    00232
    00233         del_node(v);
    00234
    00235         return a;
    00236 }
    00237
    00239 void DAG::AddEdge(Matrix& adj, leda_edge e, int srcNodeIdx, int tgtNodeIdx) const
    00240 {
    00241         double w = GetEdgeWeight(e);
    00242
    00243         adj(srcNodeIdx, tgtNodeIdx) = w;
    00244         adj(tgtNodeIdx, srcNodeIdx) = -w;
    00245 }
    00246
    00273 double DAG::ComputeEigenSum(const Matrix& adj, int nVals)
    00274 {
    00275         using namespace NEWMAT;
    00276
    00277         int n = adj.Nrows();
    00278         Matrix U(n, n);
    00279         DiagonalMatrix D(n);
    00280         double td = 0.0;
    00281         double eigenVals[n];
    00282
    00283         //nVals *= 2; // Values will be duplicated, so step will be 2.
    00284         //ASSERT(nVals <= n);
    00285
    00286         SVD(adj, D, U);
    00287
    00288         for (int i = 1; i <= nVals; i++)
    00289                 td += D(i, i);
    00290
    00291         return td;
    00292 }
    00293
    00294 double DAG::ComputeLaplacian(const Matrix& adj)
    00295 {
    00296         using namespace NEWMAT;
    00297
    00298         int n = adj.Nrows(), i, j;
    00299         SymmetricMatrix L(n);
    00300         ColumnVector T(n);
    00301
    00302         for (i = 1; i <= n; i++)
    00303         {
    00304                 T(i) = 0;
    00305
    00306                 for (j = 1; j <= n; j++)
    00307                         T(i) += fabs(adj(i, j));
    00308         }
    00309
    00310         for (i = 1; i <= n; i++)
    00311         {
    00312                 for (j = 1; j <= n; j++)
    00313                 {
    00314                         if (i == j)
    00315                                 L(i, j) = 1;
    00316                         else if (adj(i, j) != 0)
    00317                                 L(i, j) = -1 / (sqrt(T(i) * T(j)));
    00318                         else
    00319                                 L(i, j) = 0;
    00320                 }
    00321         }
    00322
    00323         DiagonalMatrix D(n);
    00324
    00325         Jacobi(L, D);
    00326
    00327         SortAscending(D);
    00328
    00329         return fabs(D(1,1));
    00330 }
    00331
    00335 double DAG::ComputeTSVs(leda_node root)
    00336 {
    00337         int j = 1, n = GetNodeMass(root);
    00338         ASSERT(n > 0);
    00339
    00340         NodeIndexMap loopyNodeMap;
    00341         Matrix adj(n, n);
    00342         adj = 0.0;
    00343
    00344         return ComputeTSVs(root, j, adj, loopyNodeMap);
    00345 }
    00346
    00361 double DAG::ComputeTSVs(leda_node v, int& j, Matrix& adj, NodeIndexMap& loopyNodeMap)
    00362 {
    00363         using namespace NEWMAT;
    00364
    00365         leda_node w;
    00366         leda_edge e;
    00367         double s;
    00368         int i = j;
    00369         bool bDefined, bHasLoopyChild = false;
    00370
    00371         TSV& tsv = GetNodeTSV(v);
    00372         tsv.ReSize(outdeg(v), true);
    00373
    00374         forall_adj_edges(e, v)
    00375         {
    00376                 w = target(e);
    00377                 bDefined = indeg(w) > 1 ? loopyNodeMap.defined(w):false;
    00378
    00379                 if (bDefined)
    00380                 {
    00381                         bHasLoopyChild = true;
    00382                         AddEdge(adj, e, i, loopyNodeMap[w]);
    00383                         s = GetEigenLbl(w);
    00384                 }
    00385                 else
    00386                 {
    00387                         AddEdge(adj, e, i, ++j);
    00388
    00389                         if (indeg(w) > 1)
    00390                                 loopyNodeMap[w] = j; // loopy node visited
    00391
    00392                         s = ComputeTSVs(w, j, adj, loopyNodeMap);
    00393                 }
    00394
    00395                 tsv.Add(s);
    00396         }
    00397
    00398         // If v has loopy children and it isn't the root of the adj matrix
    00399         // we'll compute its TSV sum by creating a new adj mat rooted at v.
    00400         if (bHasLoopyChild && i > 1)
    00401                 s = ComputeTSVs(v);
    00402         else
    00403         {
    00404                 s = ComputeEigenSum(adj.SubMatrix(i, j, i, j), outdeg(v));
    00405                 //s = ComputeLaplacian(adj.SubMatrix(i, j, i, j));
    00406
    00407                 SetEigenLbl(v, s);
    00408                 tsv.Sort();
    00409         }
    00410
    00411         return s;
    00412 }
    00413
    00414 leda_node DAG::GetNode(int nIndex)
    00415 {
    00416         const leda_list<leda_node>& nodeList = all_nodes();
    00417         leda_list_item it = nodeList.get_item(nIndex);
    00418
    00419         return nodeList[it];
    00420 }
    00421
    00423 leda_node DAG::GetFirstRootNode() const
    00424 {
    00425         leda_node v;
    00426
    00427         forall_nodes(v, *this)
    00428                 if (indeg(v) == 0)
    00429                         return v;
    00430
    00431         return nil;
    00432 }
    00433
    00434 /*static*/
    00435 /*leda_node DAG::WeightNodeSimilarity(const double& sim, leda_node v) const
    00436 {
    00437         // The cumulative mass goes from 2N - 1 for the star graph
    00438         // to N(N+1)/2 for a path graph.
    00439         // We can convert the relative mass of a node to a value between 0 and 1
    00440         // by the cumulative distribution function for a continuous uniform distribution.
    00441         // (x - a) / (b - a)
    00442
    00443         double relMass = GetNodeMass(v) / double(nCumulativeMass); // relative mass for the node
    00444         int n = GetNodeCount();
    00445
    00446         double w = (relMass - 2 * n + 1) / ((n * n) / 2 - 1.5 * n + 1);
    00447
    00448         return .9 * sim + .1 * w;
    00449 }*/
    00450
    00463 /*static*/ leda_edge DAG::ComputeBipartiteGraph(LEDA_GRAPH<leda_node, double>& G, const DAG& g1, const DAG& g2,
    00464                                                                          SimMatrix& simMat, const DAGNodeMap& paired1,
    00465                                                                          const DAGNodeMap& paired2, MatchedNodePair* pMatchedPair /*= NULL*/)
    00466 {
    00467         //g1 is the query graph and g2 is the model graph
    00468         leda_node u, v, g1Node, g2Node;
    00469
    00470         leda_list<leda_node> A, B;
    00471         leda_list_item lia, lib;
    00472
    00473         // Create a graph with n1 + n2 nodes
    00474         forall_nodes(u, g1)
    00475                 if (pMatchedPair->GetQueryNode() != g1.GetNodeDFSIndex(u))
    00476                         A.push(G.new_node(u));
    00477
    00478         forall_nodes(u, g2)
    00479                 if (pMatchedPair->GetModelNode() != g2.GetNodeDFSIndex(u))
    00480                         B.push(G.new_node(u));
    00481
    00482         leda_node_array<leda_node>& nodes = G.node_data();
    00483
    00484         // Create a bipartite graph between A and B nodes.
    00485         // do not establish edges between nodes that have been matched.
    00486         for (lia = A.first(); lia != nil; lia = A.succ(lia))
    00487         {
    00488                 // pick node from A
    00489                 u = A[lia];
    00490                 g1Node = nodes[u];
    00491
    00492                 if (paired1.defined(g1.GetNode(g1Node)))
    00493                         continue;       // loopy node that has been matched
    00494
    00495
    00496                 for (lib = B.first(); lib != nil; lib = B.succ(lib))
    00497                 {
    00498                         // pick node from B
    00499                         v = B[lib];
    00500                         g2Node = nodes[v];
    00501
    00502                         if (paired2.defined(g2.GetNode(g2Node)))
    00503                                 continue;       // loopy node that has been matched
    00504
    00505                         // establish an edge in proper direction if needed
    00506                         //if (g1.AreNodesRelated(g1Node, g2, g2Node))
    00507                         if (simMat[g1.GetNodeDFSIndex(g1Node)][g2.GetNodeDFSIndex(g2Node)] > 0.0)
    00508                                 G.new_edge(u, v, 1.0);
    00509                 }
    00510         }
    00511
    00512         //DBG_MSG1("Number of edges in Bipartite graph: ", G.number_of_edges());
    00513
    00514         if (G.number_of_edges() == 0)
    00515                 return nil;
    00516
    00517         // Assign weights to edges
    00518         leda_edge_array<double> edgeCosts(G);
    00519         leda_edge e;
    00520         double w;
    00521         int i, j;
    00522
    00523         ASSERT(s_matchParams.dBreakSibRelPen >= 0 && s_matchParams.dBreakSibRelPen <= 1);
    00524
    00525         forall_edges(e, G)
    00526         {
    00527                 g1Node = nodes[G.source(e)];
    00528                 g2Node = nodes[G.target(e)];
    00529
    00530                 i = g1.GetNodeDFSIndex(g1Node);
    00531                 j = g2.GetNodeDFSIndex(g2Node);
    00532
    00533                 if (DAG::GetMatchParams().nPreserveAncestorRel)
    00534                 {
    00535                         if (!pMatchedPair->IsEmpty() && !pMatchedPair->AnsestorRelPreserved(i, j))
    00536                                 simMat[i][j] = 0;
    00537                 }
    00538
    00539                 edgeCosts[e] = simMat[i][j];
    00540         }
    00541
    00542         /*GraphWin gw(G);
    00543         gw.display();
    00544         gw.place_into_win();
    00545         gw.edit();*/
    00546
    00547         // Note from the LEDA manual about the 'arithmetic demand' of MWMCB_MATCHING:
    00548         // "This scaling function should be used for the algorithm that compute minimum of
    00549         // maximum weight assignments or maximum weighted matchings of maximum cardinality." .
    00550         // MWA_SCALE_WEIGHTS(G, edgeCosts); // if we use MWMCB_MATCHING
    00551         MWBM_SCALE_WEIGHTS(G, edgeCosts); // if we use MAX_WEIGHT_BIPARTITE_MATCHING
    00552
    00553         // Solve the max weight max card bipartite matching
    00554         //leda_list<leda_edge> tmch = MWMCB_MATCHING(G, A, B, edgeCosts);
    00555         leda_list<leda_edge> tmch = MAX_WEIGHT_BIPARTITE_MATCHING(G, edgeCosts);
    00556
    00557         // Find the largest weight edge from the matching
    00558         leda_list_item lie;
    00559         double    sim, curr_weight, max_weight = -9999999;
    00560         leda_edge curr_edge, max_edge = nil;
    00561
    00562         const double g1CMass = g1.nCumulativeMass;  //cumulative mass for query
    00563         const double g2CMass = g2.nCumulativeMass;  //cumulative mass for model
    00564
    00565         double g1RelMass, g2RelMass;
    00566
    00567         // Only for debug
    00568         leda_node a = nil, b = nil, c = nil, d = nil;
    00569         double x,y, maxsim = -1;
    00570         int m1, m2;
    00571         // end debug
    00572
    00573         const double smw = s_matchParams.dSimilMassWeight;
    00574         ASSERT(smw >=0 && smw <= 1);
    00575
    00576         /*
    00577                 Note that node masses are relative to the original graph the belong to
    00578                 so the relative mass should be according to the size of that graph or to the root
    00579                 of the graph the belong to.
    00580         */
    00581
    00582         forall_items(lie, tmch)
    00583         {
    00584                 curr_edge = tmch[lie];
    00585
    00586                 g1Node = nodes[G.source(curr_edge)];
    00587                 g2Node = nodes[G.target(curr_edge)];
    00588
    00589                 g1RelMass = g1.GetNodeMass(g1Node) / g1CMass;
    00590                 g2RelMass = g2.GetNodeMass(g2Node) / g2CMass;
    00591
    00592                 ASSERT(g1RelMass >= 0 && g1RelMass <= 1);
    00593                 ASSERT(g2RelMass >= 0 && g2RelMass <= 1);
    00594
    00595                 sim = edgeCosts[curr_edge];
    00596
    00597                 //curr_weight = s_smw * sim * g1RelMass + (1 - s_smw) * sim * g2RelMass;
    00598                 curr_weight = smw * sim + (1 - smw) * MAX(g1RelMass, g2RelMass);
    00599
    00600                 ASSERT(curr_weight >=0 && curr_weight <= 1)
    00601
    00602                 if (DAG::IsDbgMode() && sim > maxsim)
    00603                 {
    00604                         maxsim = sim;
    00605                         c = g1Node;
    00606                         d = g2Node;
    00607                 }
    00608
    00609                 if (curr_weight > max_weight)
    00610                 {
    00611                         max_weight = curr_weight;
    00612                         max_edge = curr_edge;
    00613
    00614                         // only for debug
    00615                         if (DAG::IsDbgMode())
    00616                         {
    00617                         a = g1Node;
    00618                         b = g2Node;
    00619                         x = edgeCosts[curr_edge];
    00620                         y = MAX(g1RelMass, g2RelMass);
    00621                         m1 = g1.GetNodeMass(g1Node);
    00622                         m2 = g2.GetNodeMass(g2Node);
    00623                 }
    00624                 }
    00625         }
    00626
    00627         // Debug stuff
    00628         if (DAG::IsDbgMode() && max_edge != nil)
    00629         {
    00630                 logFile.Print("\nQuery node: %s (%d)", (LPCTSTR) g1.GetNodeLbl(a), g1.GetNodeDFSIndex(a));
    00631                 logFile.Print("\nModel node: %s (%d)", (LPCTSTR) g2.GetNodeLbl(b), g2.GetNodeDFSIndex(b));
    00632
    00633         if (a != c || b != d)
    00634                 {
    00635                         logFile.Print("\nTrue max Q node: %s (%d)", (LPCTSTR) g1.GetNodeLbl(c), g1.GetNodeDFSIndex(c));
    00636                         logFile.Print("\nTrue max M node: %s (%d)", (LPCTSTR) g2.GetNodeLbl(d), g2.GetNodeDFSIndex(d));
    00637         }
    00638
    00639                 logFile << "\nM1: " << m1 << "\nM2: " << m2;
    00640                 logFile << "\nSim: " << x << "\nRel Mass: " << y << "\nWSim: " << max_weight << "\n\n";
    00641         }
    00642
    00643         if (max_edge != nil && max_weight >= s_matchParams.dMinNodeSim)
    00644                 G[max_edge] = max_weight;
    00645         else if (max_edge != nil)
    00646                 max_edge = nil;
    00647
    00648         return max_edge;
    00649 }
    00650
    00661 /*static*/
    00662 double DAG::Match(DAG& g1, DAG& g2, SimMatrix& simMat, DAGNodeMap& paired1,
    00663                                   DAGNodeMap& paired2, MatchedNodePair* pMatchedPair /*= NULL*/)
    00664 {
    00665         int n1 = g1.GetNodeCount();
    00666         int n2 = g2.GetNodeCount();
    00667
    00668         // Edges should be directed from a smaller to a larger set.
    00669         // but s_rwm and swm are set assuming g1 is the query!!!!! so...
    00670         /*if (n1 > n2)
    00671         {
    00672                 DBG_MSG("SWAPPING...");
    00673                 ASSERT(false); //transpose simMat
    00674                 if (pMatchedPair)
    00675                         pMatchedPair->SwapNodes();
    00676
    00677                 return Match(g2, g1, simMat, paired2, paired1, pMatchedPair);
    00678         }*/
    00679
    00680         // Check that both graphs are of non-zero size
    00681         if (n1 * n2 == 0)
    00682                 return 0;
    00683
    00684         // Check that both sizes are >1, since if one of them
    00685         // is of size one, it has been matched (unless very first call).
    00686         if (((n1 == 1) || (n2 == 1)) && paired1.size() > 0)
    00687                 return 0;
    00688
    00689         // Create a bipartite graph between A and B nodes and
    00690         // check that the number of edges is non-zero.
    00691         LEDA_GRAPH<leda_node, double> G; // graph of meta-nodes.
    00692
    00693         leda_edge max_edge = ComputeBipartiteGraph(G, g1, g2, simMat, paired1, paired2, pMatchedPair);
    00694
    00695         if (max_edge == nil)
    00696                 return 0;
    00697
    00698         // Take the endpoints of the max weight edge
    00699         leda_node maxG1Node, maxG2Node;
    00700         leda_node_array<leda_node>& nodes = G.node_data();
    00701
    00702         maxG1Node = nodes[G.source(max_edge)];
    00703         maxG2Node = nodes[G.target(max_edge)];
    00704
    00705         // Append the found nodes to the matching dict.
    00706         // mark one with another's global number.
    00707         DAGNodePtr maxG1NodePtr = g1.GetNode(maxG1Node);
    00708         DAGNodePtr maxG2NodePtr = g2.GetNode(maxG2Node);
    00709
    00710         paired1[maxG1NodePtr] = maxG2NodePtr;
    00711         paired2[maxG2NodePtr] = maxG1NodePtr;
    00712
    00713         //MatchedNodePair matNodes(maxG1Node, g1, maxG2Node, g2);
    00714
    00715         // Compute rooted subgraphs and remainders
    00716         // We want to leave the splitting node in both subgraphs
    00717   DAGPtr sg1 = g1.SplitSubGraph(maxG1Node);
    00718         DAGPtr sg2 = g2.SplitSubGraph(maxG2Node);
    00719
    00720         int idxG1Node = g1.GetNodeDFSIndex(maxG1Node);
    00721         int idxG2Node = g2.GetNodeDFSIndex(maxG2Node);
    00722
    00723         double max_weight = simMat[idxG1Node][idxG2Node];
    00724
    00726         //For display purposes, save the similarity assigned to the node
    00727         DAGNode* pNode = (DAGNode*)(const DAGNode*)maxG2NodePtr;
    00728         pNode->SetSimilarity(max_weight);
    00730
    00731         // Recurse on the root-matched subtree first and on the remainder graph
    00732         pMatchedPair->SetEmpty();
    00733         max_weight += Match(*sg1, *sg2, simMat, paired1, paired2, pMatchedPair);
    00734
    00735         pMatchedPair->SetNodes(idxG1Node, idxG2Node);
    00736         pMatchedPair->UpdateSimMat(simMat, 1 - s_matchParams.dBreakSibRelPen);
    00737         max_weight += Match(g1, g2, simMat, paired1, paired2, pMatchedPair);
    00738
    00739         return max_weight;
    00740 }
    00741
    00751 /*static*/
    00752 double DAG::Similarity(const DAG& query, const DAG& model,
    00753                 DAGNodeMap& nodeMap, bool bRecomputeTSVs)
    00754 {
    00755         DAGNodeMap nodeMap2;
    00756         leda_node q, m;
    00757         int i, j;
    00758
    00759         int n1 = query.GetNodeCount();
    00760         int n2 = model.GetNodeCount();
    00761         int min = MIN(n1, n2);
    00762
    00763         // Create a similarity matrinx and a correspondence matrix
    00764         SmartMatrix<double> simMat, dbgMat;
    00765
    00766         simMat.ReSize(n1, n2, true);
    00767
    00768         // Populate the similarity matrix
    00769         forall_nodes(q, query)
    00770         {
    00771                 forall_nodes(m, model)
    00772                 {
    00773                         i = query.GetNodeDFSIndex(q);
    00774                         j = model.GetNodeDFSIndex(m);
    00775
    00776                         simMat[i][j] = query.NodeSimilarity(q, model, m);
    00777                 }
    00778         }
    00779
    00780         if (DAG::IsDbgMode())
    00781         {
    00782                 static bool bLogFileCreated = false;
    00783
    00784                 if (!bLogFileCreated)
    00785                 {
    00786                         bLogFileCreated = true;
    00787                         logFile.open("matching.log", ios::out | ios::trunc);
    00788                 }
    00789
    00790         logFile.Print("\nQuery: %s\nModel: %s\n\n", (LPCTSTR) query.GetDAGLbl(), (LPCTSTR) model.GetDAGLbl());
    00791                 simMat.Print(logFile, 6, 3, true);
    00792                 logFile.Print("\n\n");
    00793                 s_matchParams.Print(logFile);
    00794                 dbgMat = simMat;
    00795         }
    00796
    00797         // Make sure nodeMap is empty
    00798         nodeMap.clear();
    00799
    00800         // DAG::Match() does not preserve the contents of its parameters so
    00801         // we need to make a copy of the graphs.
    00802         DAGPtr cpQuery = query.CreateObject();
    00803         DAGPtr cpModel = model.CreateObject();
    00804
    00805         MatchedNodePair matNodes(query, model); // graphs use to apply penalties
    00806
    00807         double dMatchedNodes = DAG::Match(*cpQuery = query, *cpModel = model,
    00808                 simMat, nodeMap, nodeMap2, &matNodes);
    00809
    00810         if (DAG::IsDbgMode())
    00811         {
    00812                 for (int i = 0; i < simMat.NRows(); i++)
    00813                         for (int j = 0; j < simMat.NCols(); j++)
    00814                                 if (simMat[i][j] - dbgMat[i][j] != 0.0)
    00815                                 {
    00816                                         if (simMat[i][j] != 0.0)
    00817                                                 dbgMat[i][j] = simMat[i][j] - dbgMat[i][j];
    00818                                 }
    00819                                 else
    00820                                         dbgMat[i][j] = 0.0;
    00821
    00822                 dbgMat.Print(logFile, 6, 3, true);
    00823         }
    00824
    00825         return (dMatchedNodes / n1 + dMatchedNodes / n2) / 2.0;
    00826         //return dMatchedNodes / n1;
    00827 }
    00828
    00836 /*static*/
    00837 double DAG::Distance(const DAG& image, const DAG& model,
    00838         DAGNodeMap& nodeMap, bool bRecomputeTSVs)
    00839 {
    00840         return 1.0 / DAG::Similarity(image, model, nodeMap, bRecomputeTSVs);
    00841 }
    00842
    00843
    00845 // Virtual functions
    00846
    00848 DAG& DAG::operator=(const DAG& rhs)
    00849 {
    00850         Clear();
    00851
    00852         DAG_BASE_CLASS::operator=(rhs);
    00853
    00854         nFileOffset = rhs.nFileOffset;
    00855         nDAGId      = rhs.nDAGId;
    00856         strGraphLbl = rhs.strGraphLbl;
    00857
    00858         // Derived values
    00859         nMaxBFactor   = rhs.nMaxBFactor;
    00860         dTotalTSVSum  = rhs.dTotalTSVSum;
    00861         nViewNumber   = rhs.nViewNumber;
    00862         strObjectName = rhs.strObjectName;
    00863         eigenVals         = rhs.eigenVals; // Experimental. It isn't necessary
    00864         transClosMat  = rhs.transClosMat;
    00865         nCumulativeMass = rhs.nCumulativeMass;
    00866
    00867         return *this;
    00868 }
    00869
    00871 void DAG::Clear()
    00872 {
    00873         DAG_BASE_CLASS::clear();
    00874
    00875         nFileOffset = 0;
    00876         nDAGId = 0;
    00877         strGraphLbl.Clear();
    00878
    00879         // Derived values
    00880         nMaxBFactor = 0;
    00881         dTotalTSVSum = 0;
    00882         nViewNumber = 0;
    00883         strObjectName.Clear();
    00884         transClosMat.Clear();
    00885         nCumulativeMass = 0;
    00886 }
    00887
    00889 double DAG::NodeSimilarity(leda_node u, const DAG& from, leda_node v) const
    00890 {
    00891         double dNodeSim  = 1 - NodeDistance(u, from, v);
    00892         double dEigenSim = NodeTSVSimilarity(u, from, v);
    00893         double dSimilarity;
    00894
    00895         if (dNodeSim < 0 || dNodeSim > 1)
    00896         {
    00897                 if (DAG::IsDbgMode())
    00898                 {
    00899                         cerr << "\nWrong similarity value (" << dNodeSim << ") when matching: \n"
    00900                                 << GetDAGLbl() << ":" << GetNodeLbl(u) << " against "
    00901                                 << from.GetDAGLbl() << ":" << from.GetNodeLbl(v) << endl;
    00902                 }
    00903
    00904                 dNodeSim = (dNodeSim < 0) ? 0:1;
    00905         }
    00906
    00907         ASSERT(!isnan(dNodeSim)); ASSERT(finite(dNodeSim));
    00908         ASSERT(!isnan(dEigenSim)); ASSERT(finite(dEigenSim));
    00909
    00910         ASSERT(dNodeSim >= 0 && dNodeSim <= 1);
    00911         ASSERT(dEigenSim >= 0 && dEigenSim <= 1);
    00912
    00913         if (dNodeSim > 0.0)
    00914                 dSimilarity = s_matchParams.dTSVSimWeight * dEigenSim + (1 - s_matchParams.dTSVSimWeight) * dNodeSim;
    00915         else
    00916                 dSimilarity = 0.0;
    00917
    00918         /*switch (s_matchParams.nSimFuncType)
    00919         {
    00920                 case EXP_NEG_SUM:
    00921                         dSimilarity = exp (-dNodeDistance - dEigenDistance);
    00922                         break;
    00923                 case INV_NODE_DIST:
    00924                         dSimilarity = 1.0 / (1 + dNodeDistance);
    00925                         break;
    00926                 case INV_NODE_DIST_AND_TSV:
    00927                         dSimilarity = 1.0 /
    00928                                 (1 + sqrt(dNodeDistance * dNodeDistance + dEigenDistance * dEigenDistance));
    00929                         break;
    00930                 case ID_NODE_DIST:
    00931                         dSimilarity = dNodeDistance;
    00932                         break;
    00933                 default:
    00934                         ASSERT(false);
    00935         }*/
    00936
    00937         ASSERT(!isnan(dSimilarity)); ASSERT(finite(dSimilarity));
    00938         ASSERT(dSimilarity >= 0); ASSERT(dSimilarity <= 1);
    00939
    00940         return dSimilarity;
    00941 }
    00942
    00943
    00945 double DAG::EigenDistance(leda_node u, const DAG& from, leda_node v) const
    00946 {
    00947         //return NodeTSVSimilarity(u, from, v);
    00948
    00949         return Norm2(GetNodeTSV(u) - from.GetNodeTSV(v));
    00950 }
    00951
    00952
    00957 void DAG::Print(ostream& os /*=cout*/) const
    00958 {
    00959         os << "\nDAG label: " << GetDAGLbl() << endl;
    00960         os << "Object: " << GetObjName() << ", view: " << GetViewNumber() << endl;
    00961
    00962         /*
    00963         os << "\n--------------- Vertex Set ---------------\n"; */
    00964
    00965         leda_node v;
    00966
    00967         forall_nodes(v, *this)
    00968         {
    00969                 os << "\nNode:" << index(v) << endl;
    00970                 GetNode(v)->Print(os);
    00971         }
    00972 }
    00973
    00975 void DAG::PrintTSVs(ostream& os /*=cout*/) const
    00976 {
    00977         os << "\nDAG's adj matrix eigen values:\n";
    00978
    00979         eigenVals.Print(os);
    00980
    00981         os << "\n\n--------------- TSV Set ---------------\n";
    00982
    00983         leda_node v;
    00984
    00985         forall_nodes(v, *this)
    00986         {
    00987                 os << "\nNode " << GetNodeLbl(v) << ": ";
    00988                 GetNodeTSV(v).Print(os);
    00989         }
    00990
    00991         os << "\n---------------------------------------\n";
    00992 }
    00993
    00995 istream& DAG::Read(istream& is, bool bOnlyDataForMatching /*= false*/)
    00996 {
    00997         int n, i;
    00998         leda_node v, s, t;
    00999         leda_map<long, leda_node> map;
    01000         double w;
    01001         String strClassName;
    01002
    01003         Clear();
    01004
    01005         // Save the offset where we are about to read from
    01006         nFileOffset = is.tellg();
    01007
    01008         strClassName.Read(is);
    01009
    01010         if (ClassName() != strClassName) {
    01011                 DBG_MSG1("Reading DAG from offset:", nFileOffset);
    01012                 THROW_EXCEPTION("Incorrect DAG class");
    01013         }
    01014
    01015         strGraphLbl.Read(is);
    01016
    01017         is.read((char*) &n, sizeof(n));
    01018
    01019         for(i = 0; i < n; i++)
    01020         {
    01021                 is.read((char*) &v, sizeof(v));
    01022                 map[(long)v] = NewNode(ReadNode(is));
    01023         }
    01024
    01025         is.read((char*) &n, sizeof(n));
    01026
    01027         for(i = 0; i < n; i++)
    01028         {
    01029                 is.read((char*) &s, sizeof(s));
    01030                 is.read((char*) &t, sizeof(t));
    01031                 is.read((char*) &w, sizeof(w));
    01032                 NewEdge(map[(long)s], map[(long)t], w);
    01033         }
    01034
    01035         // Read derived values
    01036         is.read((char*) &nMaxBFactor, sizeof(nMaxBFactor));
    01037         is.read((char*) &dTotalTSVSum, sizeof(dTotalTSVSum));
    01038         is.read((char*) &nViewNumber, sizeof(nViewNumber));
    01039         strObjectName.Read(is);
    01040         eigenVals.Read(is);  // Experimental. It isn't necessary
    01041
    01042         transClosMat.Read(is);
    01043         is.read((char*) &nCumulativeMass, sizeof(nCumulativeMass));
    01044         ASSERT(nCumulativeMass >= GetNodeCount());
    01045
    01046         return is;
    01047 }
    01048
    01050 ostream& DAG::Write(ostream& os) const
    01051 {
    01052         leda_node v, s, t;
    01053         leda_edge e;
    01054         double w;
    01055
    01056         // Save the offset where we are about to write to ???
    01057         //nFileOffset = os.tellp();
    01058
    01059         ClassName().Write(os);
    01060         strGraphLbl.Write(os);
    01061
    01062         int n = GetNodeCount();
    01063
    01064         os.write((char*) &n, sizeof(n));
    01065
    01066         forall_nodes(v, *this)
    01067         {
    01068                 os.write((char*) &v, sizeof(v));
    01069                 os << inf(v);
    01070         }
    01071
    01072         n = GetEdgeCount();
    01073
    01074         os.write((char*) &n, sizeof(n));
    01075
    01076         forall_edges(e, *this)
    01077         {
    01078                 s = source(e);
    01079                 t = target(e);
    01080                 w = inf(e);
    01081                 os.write((char*) &s, sizeof(s));
    01082                 os.write((char*) &t, sizeof(t));
    01083                 os.write((char*) &w, sizeof(w));
    01084         }
    01085
    01086         // Read derived values
    01087         os.write((char*) &nMaxBFactor, sizeof(nMaxBFactor));
    01088         os.write((char*) &dTotalTSVSum, sizeof(dTotalTSVSum));
    01089         os.write((char*) &nViewNumber, sizeof(nViewNumber));
    01090         strObjectName.Write(os);
    01091         eigenVals.Write(os); // Experimental. It isn't necessary
    01092
    01093         transClosMat.Write(os);
    01094         os.write((char*) &nCumulativeMass, sizeof(nCumulativeMass));
    01095
    01096         return os;
    01097 }
    01098
    01100 void DAG::DelEdge(Matrix& adj, int srcNodeIdx, int tgtNodeIdx) const
    01101 {
    01102         adj(srcNodeIdx, tgtNodeIdx) = 0.0;
    01103         adj(tgtNodeIdx, srcNodeIdx) = 0.0;
    01104 }
    01105
    01117 double DAG::NodeTSVSimilarity(leda_node u, const DAG& from, leda_node v) const
    01118 {
    01119         const TSV& tsv1 = GetNodeTSV(u);
    01120         const TSV& tsv2 = from.GetNodeTSV(v);
    01121         double n1 = GetNodeTSVNorm(u);
    01122         double n2 = from.GetNodeTSVNorm(v);
    01123
    01124         if (n1 == 0 && n2 == 0)
    01125                 return 1;
    01126
    01127         TSV diff = tsv1 - tsv2;
    01128         double diffNorm = diff.Norm2();
    01129
    01130         double max = MAX(n1, n2);
    01131
    01132         // max should be a normalization factor to ensure
    01133         // the the result would be between zero and 1. However,
    01134         // it is possible that diffNorm is grater than max, and so
    01135         // we need to handle this case somehow. Clearly, a better solution
    01136         // would be to use a better normalization factor.
    01137         if (diffNorm > max)
    01138                 return 0;
    01139
    01140         return 1 - diffNorm / max;
    01141 }
    01142
    01149 void DAG::DeleteSubDAG(leda_node v)
    01150 {
    01151         leda_node u;
    01152
    01153         forall_adj_nodes(u, v)
    01154                 DeleteSubDAG(u);
    01155
    01156         del_node(v);
    01157 }
    01158
    01162 void DAG::PrintAdjMatrix(ostream& os)
    01163 {
    01164         int j = 1, n = GetNodeCount();
    01165         NodeIndexMap loopyNodeMap;
    01166         Matrix adj(n, n);
    01167
    01168         adj = 0.0;
    01169         ComputeTSVs(GetFirstRootNode(), j, adj, loopyNodeMap);
    01170
    01171         os << "\nAdjacency matrix:\n";
    01172         os << setw(1) << setprecision(0) << adj << setprecision(6) << endl;
    01173 }
    01174
    01175 /*
    01176         Depth First Search
    01177 */
    01178 leda_node DAG::DFSGetNext(leda_node v, EdgeList& l, bool bVisitNode) const
    01179 {
    01180         leda_edge e;
    01181
    01182         ASSERT(v != nil);
    01183
    01184         if (bVisitNode)
    01185                 forall_adj_edges(e, v)
    01186                         l.append(e);
    01187
    01188         return l.empty() ? nil:target(l.pop_back());
    01189 }
    01190
    01191 /*
    01192         Breadth First Search
    01193 */
    01194 leda_node DAG::BFSGetNext(leda_node v, EdgeList& l, bool bVisitNode) const
    01195 {
    01196         leda_edge e;
    01197
    01198         ASSERT(v != nil);
    01199
    01200         if (bVisitNode)
    01201                 forall_adj_edges(e, v)
    01202                         l.append(e);
    01203
    01204         return l.empty() ? nil:target(l.pop_front());
    01205 }
    01206
    01207 void SimilarityGraph::Add(leda_node g1Node, leda_node g2Node, double dSimilarity)
    01208 {
    01209         new_edge(MapNode(g1Map, g1Node), MapNode(g2Map, g2Node), dSimilarity);
    01210 }
    01211
    01212 /*void SimilarityGraph::Add(leda_node g1Node, leda_node g2Node, double dSimilarity,
    01213                                                   leda_node g1MatchedNode, leda_node g2MatchedNode)
    01214 {
    01215         leda_node u = MapNode(g1Map, g1Node);
    01216         leda_node v = MapNode(g2Map, g2Node);
    01217
    01218         leda_edge e = new_edge(u, v, dSimilarity);
    01219
    01220         g1MatchedNodeMap[e] = g1MatchedNode;
    01221         g2MatchedNodeMap[e] = g2MatchedNode;
    01222 }*/
    01223
    01224 /*
    01225         Remaining node, deleting node
    01226 */
    01227 void SimilarityGraph::MergeG2Nodes(leda_node remNode, leda_node delNode)
    01228 {
    01229         leda_node r, d;
    01230         leda_edge e;
    01231
    01232         r = g2Map[remNode];
    01233         d = g2Map[delNode];
    01234
    01235         forall_in_edges(e, d)
    01236                 move_edge(e, source(e), r);
    01237
    01238         g2Map.undefine(delNode);
    01239         del_node(d);
    01240 }
    01241
    01242 double SimilarityGraph::ComputeBestAssignment(METHOD nMethod, int* pnMatched,
    01243                                                                                           NodeMap* pNodeMap /*= NULL*/)
    01244 {
    01245         leda_edge_array<double> edgeCosts(*this, 0.0);
    01246         leda_edge e;
    01247
    01248         forall_edges(e, *this)
    01249         {
    01250                 edgeCosts[e] = inf(e);
    01251                 operator[](e) = 1.0;
    01252         }
    01253
    01254         //MWA_SCALE_WEIGHTS(*this, edgeCosts);
    01255         MWBM_SCALE_WEIGHTS(*this, edgeCosts);
    01256
    01257         EdgeList l;
    01258
    01259         if (nMethod == MWMC)
    01260         {
    01261                 NodeList setA, setB;
    01262                 bool bBipartite = Is_Bipartite(*this, setA, setB);
    01263
    01264                 ASSERT(bBipartite);
    01265                 l = MWMCB_MATCHING(*this, setA, setB, edgeCosts);
    01266         }
    01267         else
    01268                 l = MAX_WEIGHT_BIPARTITE_MATCHING(*this, edgeCosts);
    01269
    01270         leda_list_item i;
    01271         leda_node u, v;
    01272         double s = 0;
    01273
    01274         forall_items(i, l)
    01275         {
    01276                 e = l[i];
    01277                 s += inf(e);
    01278
    01279                 if (pNodeMap)
    01280                 {
    01281                         if (g1MatchedNodeMap.defined(e))
    01282                         {
    01283                                 u = g1MatchedNodeMap[e];
    01284                                 v = g2MatchedNodeMap[e];
    01285                         }
    01286                         else
    01287                         {
    01288                                 u = inf(source(e));
    01289                                 v = inf(target(e));
    01290                         }
    01291
    01292                         (*pNodeMap)[u] = v;
    01293                 }
    01294         }
    01295
    01296         if (pnMatched)
    01297                 *pnMatched = l.size();
    01298
    01299         return s;
    01300 }
    01301
    01302 DAG::MatchParams::MatchParams()
    01303 {
    01304         memset(this, 0, sizeof(*this));
    01305 }
    01306
    01307 void DAG::MatchParams::Print(ostream& os) const
    01308 {
    01309         PRINT_OPEN(os, (int)nSimFuncType);
    01310         PRINT(os, nNodeDistFunc);
    01311         PRINT(os, dTSVSimWeight);
    01312         PRINT(os, dSimilMassWeight);
    01313         PRINT(os, dRelMassWeight);
    01314         PRINT(os, dMinNodeSim);
    01315         PRINT(os, dBreakSibRelPen);
    01316         PRINT(os, nPreserveAncestorRel);
    01317         PRINT(os, nUseNewVoteWeightFunc);
    01318         PRINT_CLOSE(os, nUseMOOVC);
    01319 }
    01320
    01321 /*static*/
    01322 NodeIdxArray MatchedNodePair::GetSiblingsVector(int node, const SmartMatrix<int>& tcm, const NodeIdxMatrix& parents)
    01323 {
    01324         const NodeIdxArray& vpar = parents[node];
    01325         NodeIdxArray sibs;
    01326         int parent;
    01327
    01328         // look for all the node tha can be reached from node's parents
    01329         // i.e., the node's siblings and their descendents
    01330
    01331         for (int i = 0; i < vpar.GetSize(); i++)
    01332         {
    01333                 parent = vpar[i];
    01334
    01335                 if (sibs.GetSize() == 0)
    01336                         sibs = tcm[parent];
    01337                 else
    01338                 {
    01339                         const NodeIdxArray& halfSibs = tcm[parent];
    01340                         ASSERT(sibs.GetSize() == halfSibs.GetSize());
    01341
    01342                         // Do sibs |= halfSibs;
    01343                         for (int j = 0; j < halfSibs.GetSize(); j++)
    01344                                 if (halfSibs[j] == 1 && sibs[j] != 1) // only modifies sibs if necessary
    01345                                         sibs[j] = 1;
    01346                 }
    01347
    01348                 // we consider that a node is a sibling of itself so that case is okay. But,
    01349                 // the node's parent isn't a sibling, so we must correct for this.
    01350                 sibs[parent] = 0;
    01351         }
    01352
    01353         if (sibs.Sum() <= 1)
    01354                 sibs.Clear();
    01355
    01356         return sibs;
    01357 }
    01358
    01377 void MatchedNodePair::UpdateSimMat(SimMatrix& simMat, const double& alpha) const
    01378 {
    01379         int rows = simMat.NRows(), cols = simMat.NCols();
    01380         ASSERT(alpha >= 0 && alpha <= 1);
    01381
    01382         if (alpha == 1) //nothing to do
    01383                 return;
    01384
    01385         SmartArray<int> sq, sm;
    01386         int i, j;
    01387
    01388         sq = GetSiblingsVector(q, *pQueryTCM, qparents);
    01389         sm = GetSiblingsVector(m, *pModelTCM, mparents);
    01390
    01391         rows = sq.GetSize();
    01392         cols = sm.GetSize();
    01393
    01394         // rows or cols will be zero when there are no siblings
    01395         // if either q or m has no siblings, there is no point in
    01396         // penalizing the siblings that do exist, because there is no
    01397         // sibling relations to preserve.
    01398         if (rows == 0 || cols == 0)
    01399                 return;
    01400
    01401         // Update the similarity matrix
    01402         for (i = 0; i < rows; i++)
    01403                 if (sq[i] == 1)                                         // if 'v' in Q is a sibling of 'q'
    01404                         for (j = 0; j < cols; j++)              // travers its matches from M
    01405                                 if (sm[j] == 0)                         // if 'w' in M isn't a sibling of 'm'
    01406                                         simMat[i][j] *= alpha;
    01407
    01408         for (j = 0; j < cols; j++)
    01409                 if (sm[j] == 1)                                         // if 'v' in Q is a sibling of 'q'
    01410                         for (i = 0; i < rows; i++)              // travers its matches from M
    01411                                 if (sq[i] == 0)                         // if 'w' in M isn't a sibling of 'm'
    01412                                         simMat[i][j] *= alpha;
    01413 }
    01414
    01415