ext/OGDF/src/energybased/StressMajorizationSimple.cpp

   1 /*
   2  * $Revision: 2565 $
   3  *
   4  * last checkin:
   5  *   $Author: gutwenger $
   6  *   $Date: 2012-07-07 17:14:54 +0200 (Sa, 07. Jul 2012) $
   7  ***************************************************************/
   8
   9 /** \file
  10  * \brief Implementation of stress-majorization layout algorithm.
  11  *
  12  * \author Karsten Klein
  13  *
  14  * \par License:
  15  * This file is part of the Open Graph Drawing Framework (OGDF).
  16  *
  17  * \par
  18  * Copyright (C)<br>
  19  * See README.txt in the root directory of the OGDF installation for details.
  20  *
  21  * \par
  22  * This program is free software; you can redistribute it and/or
  23  * modify it under the terms of the GNU General Public License
  24  * Version 2 or 3 as published by the Free Software Foundation;
  25  * see the file LICENSE.txt included in the packaging of this file
  26  * for details.
  27  *
  28  * \par
  29  * This program is distributed in the hope that it will be useful,
  30  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  31  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  32  * GNU General Public License for more details.
  33  *
  34  * \par
  35  * You should have received a copy of the GNU General Public
  36  * License along with this program; if not, write to the Free
  37  * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  38  * Boston, MA 02110-1301, USA.
  39  *
  40  * \see  http://www.gnu.org/copyleft/gpl.html
  41  ***************************************************************/
  42
  43 #include <ogdf/energybased/StressMajorizationSimple.h>
  44
  45 //For initial layouts
  46 #include <ogdf/energybased/FMMMLayout.h>
  47
  48 #include <ogdf/basic/SList.h>
  49
  50 //only debugging
  51 #include <ogdf/basic/simple_graph_alg.h>
  52
  53 namespace ogdf {
  54 const double  StressMajorization::startVal = DBL_MAX - 1.0;
  55 const double  StressMajorization::minVal = DBL_MIN;
  56 const double  StressMajorization::desMinLength = 0.0001;
  57
  58 void  StressMajorization::initialize(
  59         GraphAttributes& GA,
  60         const EdgeArray<double>& eLength,
  61         NodeArray< NodeArray<double> >& oLength,
  62         NodeArray< NodeArray<double> >& weights,
  63         double & maxDist,
  64         bool simpleBFS)
  65 {
  66         node v;
  67         const Graph &G = GA.constGraph();
  68
  69         m_prevEnergy =  startVal;
  70         m_prevLEnergy =  startVal;
  71
  72         // all edges straight-line
  73         GA.clearAllBends();
  74         if (!m_useLayout)
  75                 shufflePositions(GA);
  76
  77         //the shortest path lengths
  78         forall_nodes (v, G) oLength[v].init(G, DBL_MAX);
  79         forall_nodes (v, G) weights[v].init(G, 0.0);
  80
  81         //-------------------------------------
  82         //computes shortest path distances d_ij
  83         //-------------------------------------
  84         if (simpleBFS)
  85         {
  86                 //we use simply BFS n times
  87                 //TODO experimentally compare speed, also with bintree dijkstra
  88 #ifdef OGDF_DEBUG
  89                 double timeUsed;
  90                 usedTime(timeUsed);
  91 #endif
  92                 maxDist = allpairsspBFS(G, oLength, weights);
  93 #ifdef OGDF_DEBUG
  94                 timeUsed = usedTime(timeUsed);
  95                 cout << "\n******APSP BFS runtime: \n";
  96 #endif
  97         }
  98         else
  99         {
 100                 EdgeArray<double> adaptedLength(G);
 101                 adaptLengths(G, GA, eLength, adaptedLength);
 102                 //we use simply the BFM n times or Floyd instead, leading to cubic runtime
 103                 //TODO experimentally compare speed, also with bintree dijkstra
 104                 maxDist = allpairssp(G, adaptedLength, oLength, weights, DBL_MAX);
 105         }
 106
 107         if (m_radial)
 108         {
 109                 //TODO: Also allow different ways to compute centrality (closeness,...)
 110                 // and incorporate node sizes
 111                 //if centrality
 112                 computeRadii(G, oLength, maxDist);
 113         }
 114 //Todo: Here we could add distance values depending on given node sizes
 115 }//initialize
 116
 117
 118 //uses closeness to compute the radii
 119 void StressMajorization::computeRadii(const Graph& G, const NodeArray< NodeArray<double> >& distances,
 120         double diameter)
 121 {
 122         m_radii.init(G, 1.0);
 123         //compute a center (simple version)
 124         node v, w;
 125         double minMax = DBL_MAX;
 126         node center = 0;
 127         int numCentralNodes = 0;
 128         double maxCloseness = 0.0;
 129         double minCloseness = DBL_MAX;
 130
 131         // inverse of sum of shortest path distances
 132         NodeArray<double> closeness(G, 0.0);
 133
 134         //Compute closeness values and min/max
 135         forall_nodes(v, G)
 136         {
 137                 double maxDist = 0.0;
 138                 forall_nodes(w, G)
 139                 {
 140                         if (v != w) closeness[v] += distances[v][w];
 141                         if (distances[v][w] > maxDist) maxDist = distances[v][w];
 142                 }
 143                 if (maxDist < minMax)
 144                 {
 145                         minMax = maxDist;
 146                         center = v;
 147                 }
 148                 closeness[v] = (G.numberOfNodes()-1)/closeness[v]; //was 1/
 149
 150                 //Check for min/max closeness
 151                 if (DIsGreater(closeness[v], maxCloseness))
 152                 {
 153                         maxCloseness = closeness[v];
 154                         numCentralNodes = 1;
 155                 }
 156                 else if (DIsEqual(closeness[v], maxCloseness)) numCentralNodes++;
 157                 if (DIsGreater(minCloseness, closeness[v]))
 158                         minCloseness = closeness[v];
 159         }
 160
 161         //see paper by Brandes, Kenis, Wagner
 162         double coincideOfs = min(0.5, double(numCentralNodes)/double((G.numberOfNodes()-1)));
 163
 164         OGDF_ASSERT(center != 0);
 165         forall_nodes(v, G)
 166         {
 167                 //cout << "Diameter: "<<diameter<<"\n";
 168                 m_radii[v] = diameter/double(2.0) * (1 - (closeness[v]-minCloseness)/(maxCloseness-minCloseness+coincideOfs));
 169                 //cout << "Radius: "<<m_radii[v]<<"\n";
 170         }
 171
 172 }//computeRadii
 173
 174
 175 void StressMajorization::mainStep(GraphAttributes& GA,
 176         NodeArray< NodeArray<double> >& oLength,
 177         NodeArray< NodeArray<double> >& weights,
 178         const double maxDist)
 179 {
 180         const Graph &G = GA.constGraph();
 181         node v;
 182
 183         //preparation for upward constraints
 184         edge e;
 185         NodeArray< NodeArray<double> > upWeight(G);
 186         forall_nodes(v, G) upWeight[v].init(G, 0.0);
 187         //try to add upward constraints by looking at the adjacent edges
 188         forall_edges(e, G)
 189         {
 190                 upWeight[e->source()][e->target()] = 1.0;
 191                 upWeight[e->target()][e->source()] = -1.0;
 192         }
 193
 194 #ifdef OGDF_DEBUG
 195         NodeArray<int> nodeCounts(G, 0);
 196 #endif
 197         NodeArray<double> invPosNorm(G);
 198         NodeArray< NodeArray<double> > invDistb(G); //bij in paper Pich Brandes
 199
 200         //this value is needed quite often
 201         NodeArray<double> wijSum(G, 0.0); //sum of w_ij for each i over j!= i
 202
 203         forall_nodes(v, G)
 204         {
 205                 invDistb[v].init(G);
 206                 node w;
 207                 forall_nodes(w, G)
 208                 {
 209                         if (v != w) wijSum[v]+= weights[v][w];
 210                 }
 211         }
 212
 213         //--------------------------------------------------
 214         // Iterations can either be done until a predefined
 215         // number m_numSteps+m_itFac*|G| is reached or a break condition
 216         // is fullfilled
 217         //--------------------------------------------------
 218         //for iteration checking
 219         double weightsum = 1.0;
 220         if (m_upward) weightsum = 0.93;
 221
 222         OGDF_ASSERT(DIsGreater(m_numSteps+m_itFac*G.numberOfNodes(), 0.0));
 223         double kinv = weightsum/double(m_numSteps+m_itFac*G.numberOfNodes());
 224
 225         NodeArray<double> newX(G);
 226         NodeArray<double> newY(G);
 227         //weighting for different optimization criteria
 228         double t = 0.0;
 229
 230         for (double tt = 0; tt <= weightsum; tt+= kinv)
 231         {
 232                 double iterStressSum = 0.0;
 233                 if (m_radial || m_upward) t = tt;
 234                 else t = 0.0;
 235                 //compute inverse norm of position for radial distance constraints
 236                 forall_nodes(v, G)
 237                 {
 238                         double px = GA.x(v);
 239                         double py = GA.y(v);
 240                         //set values a_i as in paper pseudo code (simplified)
 241                         double tmp = px*px+py*py;
 242                         if (DIsGreater(tmp, 0.0))
 243                                 invPosNorm[v] = 1/double(sqrt(tmp));
 244                         else invPosNorm[v] = 0.0;
 245                         //cout << "invPosNorm: "<<invPosNorm[v]<<"\n";
 246
 247                         //set values b_ij as in paper pseudo code (simplified)
 248                         node w;
 249                         forall_nodes(w, G)
 250                         {
 251                                 double tmpij = (px-GA.x(w))*(px-GA.x(w))+(py-GA.y(w))*(py-GA.y(w));
 252                                 invDistb[v][w] = (DIsGreater(tmpij, 0.0) ? 1/double(sqrt(tmpij)) : 0.0);
 253                                 //cout << "invdistb: "<<invDistb[v][w]<<"\n";
 254                         }
 255
 256                 }
 257
 258                 forall_nodes(v, G)
 259                 {
 260                         //value corresponding to radial constraint
 261                         double radOfsX = 0.0;
 262                         double radOfsY = 0.0;
 263                         double tmpinvsq = 0.0;
 264
 265                         double upWeightSum = 0.0;
 266
 267                         if(m_upward)
 268                         {
 269                                 upWeightSum = 0.05+v->degree()/100.0;
 270                         }
 271
 272                         if(m_radial)
 273                         {
 274                                 double tmp = radius(v);
 275                                 OGDF_ASSERT(DIsGreater(tmp, 0.0))
 276                                         tmpinvsq = 1/(tmp*tmp);
 277                                 double tmppart = t*tmpinvsq*tmp*invPosNorm[v];
 278
 279                                 radOfsX = tmppart*GA.x(v);
 280                                 radOfsY = tmppart*GA.y(v);
 281                         }
 282
 283                         //upward constraints
 284                         double upOfs = 0.0; //only for y coordinate
 285
 286                         //sum over all other nodes
 287                         node w;
 288                         double stressSumX = 0.0;
 289                         double stressSumY = 0.0;
 290                         forall_nodes(w, G)
 291                         {
 292                                 if ( v != w )
 293                                 {
 294                                         stressSumX += weights[v][w]*(GA.x(w)+oLength[v][w]*(GA.x(v)-GA.x(w))*invDistb[v][w]);
 295                                         stressSumY += weights[v][w]*(GA.y(w)+oLength[v][w]*(GA.y(v)-GA.y(w))*invDistb[v][w]);
 296                                         //also add the influence of adjacent nodes that are not placed correctly
 297                                         if (m_upward)
 298                                         {
 299                                                 double val = upWeight[v][w];
 300                                                 //adjacent nodes influence each other
 301                                                 if (!DIsEqual(val, 0.0))
 302                                                 {
 303                                                         if (DIsGreater(val, 0.0)) //v is source
 304                                                         {
 305                                                                 if (GA.y(v) > GA.y(w)-1.0)
 306                                                                 {
 307                                                                         upOfs -= val*(GA.y(w)-1.0*invDistb[v][w]);
 308                                                                 }
 309                                                         }
 310 //                                                      else //v is target
 311 //                              {
 312 //                                      if (GA.y(w) > GA.y(v)-1.0)
 313 //                                      {
 314 //                                              upOfs -= val*(GA.y(w)+(fabs(GA.y(v)-GA.y(w)))*invDistb[v][w]);
 315 //                                      }
 316 //                              }
 317                                                 }
 318
 319                                         }
 320                                 }
 321                         }
 322                         upOfs *= upWeightSum;
 323                         if (m_radial || m_upward)
 324                         {
 325                                 stressSumX *= (1-t);
 326                                 stressSumY *= (1-t);
 327                         }
 328
 329                         //main fraction
 330                         newX[v] = (stressSumX+radOfsX)/double(((1-t)*wijSum[v]+t*tmpinvsq));
 331                         //experimental, only valid if disjoint constraints
 332                         newY[v] = (stressSumY+radOfsY+0.2*upOfs)/double(((1-t)*wijSum[v]+t*tmpinvsq+0.2*upWeightSum));
 333                         GA.x(v) = newX[v];
 334                         GA.y(v) = newY[v];
 335                         iterStressSum += stressSumX;
 336                         iterStressSum += stressSumY;
 337                 }
 338                 //Alternatively, we can do the update after the computation step
 339                 //            forall_nodes(v, G)
 340                 //            {
 341                 //GA.x(v) = newX[v];
 342                 //GA.y(v) = newY[v];
 343                 //            }
 344                 if (finished(iterStressSum))
 345                 {
 346                         cout << iterStressSum<<"\n";
 347                         break;
 348                 }
 349         }//outer loop: iterations or threshold
 350
 351 }//mainStep
 352
 353
 354 void StressMajorization::doCall(GraphAttributes& GA, const EdgeArray<double>& eLength, bool simpleBFS)
 355 {
 356         const Graph& G = GA.constGraph();
 357         double maxDist; //maximum distance between nodes
 358         NodeArray< NodeArray<double> > oLength(G);//first distance, then original length
 359         NodeArray< NodeArray<double> > weights(G);//standard weights as in MCGee,Kamada/Kawai
 360
 361         //only for debugging
 362         OGDF_ASSERT(isConnected(G));
 363
 364         //compute relevant values
 365         initialize(GA, eLength, oLength, weights, maxDist, simpleBFS);
 366
 367         //main loop with node movement
 368         mainStep(GA, oLength, weights, maxDist);
 369
 370         if (simpleBFS) scale(GA);
 371 }
 372
 373
 374 void StressMajorization::call(GraphAttributes& GA)
 375 {
 376         const Graph &G = GA.constGraph();
 377         if(G.numberOfEdges() < 1)
 378                 return;
 379
 380         EdgeArray<double> eLength(G);//, 1.0);is not used
 381         doCall(GA, eLength, true);
 382 }//call
 383
 384
 385 void  StressMajorization::call(GraphAttributes& GA,  const EdgeArray<double>& eLength)
 386 {
 387         const Graph &G = GA.constGraph();
 388         if(G.numberOfEdges() < 1)
 389                 return;
 390
 391         doCall(GA, eLength, false);
 392 }//call with edge lengths
 393
 394
 395 //changes given edge lengths (interpreted as weight factors)
 396 //according to additional parameters like node size etc.
 397 void  StressMajorization::adaptLengths(
 398         const Graph& G,
 399         const GraphAttributes& GA,
 400         const EdgeArray<double>& eLengths,
 401         EdgeArray<double>& adaptedLengths)
 402 {
 403         //we use the edge lengths as factor and try to respect
 404         //the node sizes such that each node has enough distance
 405         edge e;
 406         //adapt to node sizes
 407         forall_edges(e, G)
 408         {
 409                 double smax = max(GA.width(e->source()), GA.height(e->source()));
 410                 double tmax = max(GA.width(e->target()), GA.height(e->target()));
 411                 if (smax+tmax > 0.0)
 412                         adaptedLengths[e] = (1+eLengths[e])*((smax+tmax));///2.0);
 413                 else adaptedLengths[e] = 5.0*eLengths[e];
 414         }
 415 }//adaptLengths
 416
 417
 418 void  StressMajorization::shufflePositions(GraphAttributes& GA)
 419 {
 420         //random layout? FMMM? classical MDS? see Paper of Pich and Brandes
 421         //just hope that we have low sigma values (distance error)
 422         FMMMLayout fm;
 423         //fm.call(GA);
 424 }//shufflePositions
 425
 426
 427
 428 /**
 429  * Initialise the original estimates from nodes and edges.
 430  */
 431
 432 //we could speed this up by not using nested NodeArrays and
 433 //by not doing the fully symmetrical computation on undirected graphs
 434 //All Pairs Shortest Path Floyd, initializes the whole matrix
 435 //returns maximum distance. Does not detect negative cycles (lead to neg. values on diagonal)
 436 //threshold is the value for the distance of non-adjacent nodes, distance has to be
 437 //initialized with
 438 // The weight parameter here is just for the stress majorization
 439 // and directly set here for speedup
 440 double StressMajorization::allpairssp(
 441         const Graph& G,
 442         const EdgeArray<double>& eLengths,
 443         NodeArray< NodeArray<double> >& distance,
 444         NodeArray< NodeArray<double> >& weights,
 445         const double threshold)
 446 {
 447         node v;
 448         edge e;
 449         double maxDist = -threshold;
 450
 451         forall_nodes(v, G)
 452         {
 453                 distance[v][v] = 0.0f;
 454                 weights[v][v] = 0.0f;
 455         }
 456
 457         //TODO: Experimentally compare this against
 458         // all nodes and incident edges (memory access) on huge graphs
 459         forall_edges(e, G)
 460         {
 461                 distance[e->source()][e->target()] = eLengths[e];
 462                 distance[e->target()][e->source()] = eLengths[e];
 463         }
 464
 465 ///**
 466 // * And run the main loop of the algorithm.
 467 // */
 468         node u, w;
 469         forall_nodes(v, G)
 470         {
 471                 forall_nodes(u, G)
 472                 {
 473                         forall_nodes(w, G)
 474                         {
 475                                 if ((distance[u][v] < threshold) && (distance[v][w] < threshold))
 476                                 {
 477                                         distance[u][w] = min( distance[u][w], distance[u][v] + distance[v][w] );
 478                                         weights[u][w] = 1/(distance[u][w]*distance[u][w]);
 479                                         //distance[w][u] = distance[u][w]; //is done anyway afterwards
 480                                 }
 481                                 if (distance[u][w] < threshold)
 482                                         maxDist = max(maxDist,distance[u][w]);
 483                         }
 484                 }
 485         }
 486         //debug output
 487 #ifdef OGDF_DEBUG
 488         forall_nodes(v, G)
 489         {
 490                 if (distance[v][v] < 0.0) cerr << "\n###Error in shortest path computation###\n\n";
 491         }
 492         cout << "Maxdist: "<<maxDist<<"\n";
 493         forall_nodes(u, G)
 494         {
 495                 forall_nodes(w, G)
 496                 {
 497                         cout << "Distance " << u->index() << " -> "<<w->index()<<" "<<distance[u][w]<<"\n";
 498                 }
 499         }
 500 #endif
 501         return maxDist;
 502 }//allpairssp
 503
 504
 505 //the same without weights, i.e. all pairs shortest paths with BFS
 506 //Runs in time |V|²
 507 //for compatibility, distances are double
 508 // The weight parameter here is just for the stress majorization
 509 // and directly set here for speedup
 510 double  StressMajorization::allpairsspBFS(
 511         const Graph& G,
 512         NodeArray< NodeArray<double> >& distance,
 513         NodeArray< NodeArray<double> >& weights)
 514 {
 515         node v;
 516         double maxDist = 0;
 517
 518         forall_nodes(v, G)
 519         {
 520                 distance[v][v] = 0.0f;
 521         }
 522
 523         v = G.firstNode();
 524
 525         //start in each node once
 526         while (v != 0)
 527         {
 528                 //do a bfs
 529                 NodeArray<bool> mark(G, true);
 530                 SListPure<node> bfs;
 531                 bfs.pushBack(v);
 532                 mark[v] = false;
 533
 534                 while (!bfs.empty())
 535                 {
 536                         node w = bfs.popFrontRet();
 537                         edge e;
 538                         double d = distance[v][w]+1.0f;
 539                         forall_adj_edges(e,w)
 540                         {
 541                                 node u = e->opposite(w);
 542                                 if (mark[u])
 543                                 {
 544                                         mark[u] = false;
 545                                         bfs.pushBack(u);
 546                                         distance[v][u] = d;
 547                                         weights[v][u] = 1/(d*d);
 548                                         maxDist = max(maxDist,d);
 549                                 }
 550                         }
 551                 }//while
 552
 553                 v = v->succ();
 554         }//while
 555         //check for negative cycles
 556         forall_nodes(v, G)
 557         {
 558                 if (distance[v][v] < 0.0) cerr << "\n###Error in shortest path computation###\n\n";
 559         }
 560
 561 //debug output
 562 #ifdef OGDF_DEBUG
 563         node u, w;
 564         cout << "Maxdist: "<<maxDist<<"\n";
 565         forall_nodes(u, G)
 566         {
 567                 forall_nodes(w, G)
 568                 {
 569                         cout << "Distance " << u->index() << " -> "<<w->index()<<" "<<distance[u][w]<<"\n";
 570                 }
 571         }
 572 #endif
 573         return maxDist;
 574 }//allpairsspBFS
 575
 576
 577 void  StressMajorization::scale(GraphAttributes& GA)
 578 {
 579 //Simple version: Just scale to max needed
 580 //We run over all nodes, find the largest distance needed and scale
 581 //the edges uniformly
 582         node v;
 583         edge e;
 584         double maxFac = 0.0;
 585         forall_edges(e, GA.constGraph())
 586         {
 587                 double w1 = sqrt(GA.width(e->source())*GA.width(e->source())+
 588                                                  GA.height(e->source())*GA.height(e->source()));
 589                 double w2 = sqrt(GA.width(e->target())*GA.width(e->target())+
 590                                                  GA.height(e->target())*GA.height(e->target()));
 591                 w2 = (w1+w2)/2.0; //half length of both diagonals
 592                 double xdist = GA.x(e->source())-GA.x(e->target());
 593                 double ydist = GA.y(e->source())-GA.y(e->target());
 594                 double elength = sqrt(xdist*xdist+ydist*ydist);
 595                 w2 = m_distFactor * w2 / elength;//relative to edge length
 596                 if (w2 > maxFac)
 597                         maxFac = w2;
 598         }
 599
 600         if (maxFac > 0.0)
 601         {
 602                 forall_nodes(v, GA.constGraph())
 603                 {
 604                         GA.x(v) = GA.x(v)*maxFac;
 605                         GA.y(v) = GA.y(v)*maxFac;
 606                 }
 607 #ifdef OGDF_DEBUG
 608                 cout << "Scaled by factor "<<maxFac<<"\n";
 609 #endif
 610         }
 611 }//Scale
 612
 613 }//namespace