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

GraphLayout.C

Go to the documentation of this file.
00001 /***************************************************************************
00002  *cr
00003  *cr            (C) Copyright 1995-2019 The Board of Trustees of the
00004  *cr                        University of Illinois
00005  *cr                         All Rights Reserved
00006  *cr
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010 * RCS INFORMATION:
00011 *
00012 *      $RCSfile: GraphLayout.C,v $
00013 *      $Author: johns $      $Locker:  $               $State: Exp $
00014 *      $Revision: 1.8 $         $Date: 2020/10/15 16:07:31 $
00015 *
00016 ***************************************************************************/
00024 // Enable the fully optimized code path, although it and the less
00025 // optimized path have different strengths and weaknesses from a 
00026 // numerical perspective.
00027 #define VMDFROPTIMIZED 1
00028 
00029 // Fruchterman-Reingold spring force graph layout  
00030 #include <math.h>
00031 #include <stdio.h>
00032 #include <algorithm>
00033 #include "GraphLayout.h"
00034 #include "utilities.h"   // needed for sincosf() on MacOS X
00035 
00036 GraphLayout::GraphLayout(int numverts, int numedges) {
00037   pos_x.set_size(numverts);
00038   pos_y.set_size(numverts);
00039   dsp_x.set_size(numverts);
00040   dsp_y.set_size(numverts);
00041 
00042   if (numedges > 0) {
00043     edges_v0.set_size(numedges);
00044     edges_v1.set_size(numedges);
00045     weights.set_size(numedges);
00046   }
00047 
00048   weight_matrix = NULL;
00049 }
00050 
00051 
00052 GraphLayout::~GraphLayout() {
00053   pos_x.clear();
00054   pos_y.clear();
00055   dsp_x.clear();
00056   dsp_y.clear();
00057 
00058   edges_v0.clear();
00059   edges_v1.clear();
00060   weights.clear();
00061 }
00062 
00063 
00064 void GraphLayout::add_edge(int v0, int v1, float weight) {
00065   // If the two indices are the same, or if the weight is zero,
00066   // do nothing (no link)
00067   if (v0 == v1 || weight == 0.f) {
00068     return;
00069   }
00070 
00071   // auto-extend the vertex array if we get an edge with a larger
00072   // index than the current size.  This is easy-to-use but inefficient
00073   // compared to initializing the array sizes at construction time.
00074   long maxvert = std::max(v0, v1);
00075   if (maxvert > pos_x.num()) {
00076     long extra = maxvert - pos_x.num(); 
00077     pos_x.extend(extra);
00078     pos_y.extend(extra);
00079     dsp_x.extend(extra);
00080     dsp_y.extend(extra);
00081   }
00082 
00083   edges_v0.append(v0);
00084   edges_v1.append(v1);
00085   weights.append(weight);
00086 }
00087 
00088 
00089 void GraphLayout::add_weight_matrix(const float *weight_matrix) {
00090   this->weight_matrix = weight_matrix;
00091 }
00092 
00093 
00094 // place nodes randomly in a square centered at origin
00095 void GraphLayout::init_positions_box() {
00096   int numverts = pos_x.num();
00097   float rm_inv = 1.0f / float(RAND_MAX);
00098 
00099   srand(314159); // make reproducible
00100   int v;
00101   for (v=0; v<numverts; v++) {
00102     pos_x[v] = rm_inv * rand() - 0.5f;
00103     pos_y[v] = rm_inv * rand() - 0.5f;
00104   }
00105 
00106   // clear displacement vectors
00107   memset(&dsp_x[0], 0, numverts * sizeof(float));
00108   memset(&dsp_y[0], 0, numverts * sizeof(float));
00109 }
00110 
00111 
00112 // place nodes randomly in a circle centered at origin
00113 void GraphLayout::init_positions_circle() {
00114   int numverts = pos_x.num();
00115   float radscale = 0.5f;
00116   float angle = 0.0f;
00117   float delta_angle = 2.0f * VMD_PIF / numverts;
00118 
00119   int v;
00120   for (v=0; v<numverts; v++) {
00121     float sa, ca;
00122     sincosf(angle, &sa, &ca);
00123     pos_x[v] = radscale * ca;
00124     pos_y[v] = radscale * sa;
00125     angle += delta_angle;
00126   }
00127 }
00128 
00129 
00130 // A variant of the Fruchterman-Reingold layout algorithm, 
00131 // with modifications to add support for weights, implicit N^2 all-to-all
00132 // node/vertex connectivity, and to allow efficient CPU vectorization 
00133 // on modern processors.  The original Fruchterman-Reingold algorithm
00134 // is described in:
00135 //  Graph Drawing by Force-Directed Placement. 
00136 //  Fruchterman, T. M. J., and Reingold, E. M.,
00137 //  Software: Practice and Experience, 21(11), 1991. 
00138 //
00139 void GraphLayout::compute(int iters, float area, float kscale, 
00140                           float tempscale, float distance_epsilon) {
00141   int numverts = pos_x.num();
00142   int numedges = edges_v0.num();
00143 
00144 printf("GraphLayout::compute(iters: %d  verts: %d  edges: %d  area: %g  kscale: %g  tempscale: %g\n", 
00145        iters, numverts, numedges, area, kscale, tempscale);
00146 
00147   if (iters < 1) 
00148     return;
00149 
00150 //  float area = 1.0f;
00151 //  float kscale = 1.0e3f;
00152 //  float tempscale = 0.2f;
00153 //  float distance_epsilon = 1.0e-06f; // prevent division by zero
00154 
00155   float k2 = area / numverts;
00156   float k = sqrtf(k2) * kscale;
00157   float k_inv = 1.0f / k;
00158   float iters_inv = 1.0f / iters; 
00159 
00160   for (int i=0; i<iters; i++) {
00161 printf("iter: %d\r", i); fflush(stdout);
00162 
00163     // Temperature cool down; starts at 1, ends at 0
00164     // (other formulas can be used for the cooling)
00165     float temp = (1.f -  i * iters_inv);
00166 
00167     // Rather than a purely linear temperature cooldown, we use a
00168     // quadratic curve, but it might be useful to allow a user-defined
00169     // exponent here, since it's an insignificant cost.
00170     // Higher exponents yield a cooldown with more iterations having
00171     // a lower temperature (smaller max displacement per step).
00172     float tempsquared = (temp * temp) * tempscale;
00173 
00174     //
00175     // calc repulsive force: Fr(d) = -k^2/d
00176     //
00177     int v0, v1;
00178     for (v0=0; v0<numverts; v0++) {
00179       float posx = pos_x[v0];
00180       float posy = pos_y[v0];
00181       float dispx = 0.0f;
00182       float dispy = 0.0f;
00183 
00184       // avoid per-iteration if test for v0!=v1, or (dx!=0, dy!=0)
00185       // by splitting the loop into two runs...
00186       for (v1=0; v1<v0; v1++) {
00187         float dx = posx - pos_x[v1];
00188         float dy = posy - pos_y[v1];
00189         float dxy2 = dx*dx + dy*dy + distance_epsilon; 
00190 #if defined(VMDFROPTIMIZED)
00191         float repulsion = k2 / dxy2;
00192         dispx += dx * repulsion;
00193         dispy += dy * repulsion;
00194 #else
00195         float d_1 = 1.0f / sqrtf(dxy2);
00196         float repulsion = k2 * d_1;
00197         dispx += dx * d_1 * repulsion;
00198         dispy += dy * d_1 * repulsion;
00199 #endif
00200       }
00201 
00202       // avoid per-iteration if test for v0!=v1, or (dx!=0, dy!=0)
00203       // by splitting the loop into two runs...
00204       for (v1=v0+1; v1<numverts; v1++) {
00205         float dx = posx - pos_x[v1];
00206         float dy = posy - pos_y[v1];
00207         float dxy2 = dx*dx + dy*dy + distance_epsilon;
00208 #if defined(VMDFROPTIMIZED)
00209         float repulsion = k2 / dxy2;
00210         dispx += dx * repulsion;
00211         dispy += dy * repulsion;
00212 #else
00213         float d_1 = 1.0f / sqrtf(dxy2);
00214         float repulsion = k2 * d_1;
00215         dispx += dx * d_1 * repulsion;
00216         dispy += dy * d_1 * repulsion;
00217 #endif
00218       }
00219 
00220       // write new displacements back...
00221       dsp_x[v0] = dispx;
00222       dsp_y[v0] = dispy;
00223     }
00224 
00225     //
00226     // calc attractive forces, only for connected vertices
00227     //   Fa(d) = d^2/k * weight -- weight factor is optional
00228     //
00229     // if the edges list is empty, assume a fully-connected graph
00230     if (numedges == 0) {
00231       if (weight_matrix == NULL) {
00232         // N^2/2 vertex combinations in fully-connected graph, no weights
00233         for (v0=0; v0<numverts; v0++) {
00234           float posx = pos_x[v0];
00235           float posy = pos_y[v0];
00236           float dispx = dsp_x[v0];
00237           float dispy = dsp_y[v0];
00238 
00239           for (v1=0; v1<v0; v1++) {
00240             float dx = posx - pos_x[v1];
00241             float dy = posy - pos_y[v1];
00242             float dxy2 = dx*dx + dy*dy;
00243 #if defined(VMDFROPTIMIZED)
00244             float attraction = sqrtf(dxy2) * k_inv;
00245             float dxa = dx * attraction;
00246             float dya = dy * attraction;
00247 #else
00248             float attraction = dxy2 * k_inv;
00249             float d_1 = 1.0f / sqrtf(dxy2 + distance_epsilon);
00250             float dxa = dx * d_1 * attraction;
00251             float dya = dy * d_1 * attraction;
00252 #endif
00253 
00254             // memory scatter here is bad for performance...
00255             dispx     -= dxa;
00256             dispy     -= dya;
00257             dsp_x[v1] += dxa;
00258             dsp_y[v1] += dya;
00259           }
00260 
00261           // write new displacements back...
00262           dsp_x[v0] = dispx;
00263           dsp_y[v0] = dispy;
00264         }
00265       } else {
00266         // N^2/2 vertex combinations in fully-connected graph
00267         for (v0=0; v0<numverts; v0++) {
00268           float posx = pos_x[v0];
00269           float posy = pos_y[v0];
00270           float dispx = dsp_x[v0];
00271           float dispy = dsp_y[v0];
00272 
00273           for (v1=0; v1<v0; v1++) {
00274             // NxN weight matrix
00275             float weight = weight_matrix[v1*numverts + v0];
00276 
00277             float dx = posx - pos_x[v1];
00278             float dy = posy - pos_y[v1];
00279             float dxy2 = dx*dx + dy*dy;
00280 #if defined(VMDFROPTIMIZED)
00281             float attraction = sqrtf(dxy2) * k_inv * weight;
00282             float dxa = dx * attraction;
00283             float dya = dy * attraction;
00284 #else
00285             float attraction = dxy2 * k_inv * weight;
00286             float d_1 = 1.0f / sqrtf(dxy2 + distance_epsilon);
00287             float dxa = dx * d_1 * attraction;
00288             float dya = dy * d_1 * attraction;
00289 #endif
00290 
00291             // memory scatter here is bad for performance...
00292             dispx     -= dxa;
00293             dispy     -= dya;
00294             dsp_x[v1] += dxa;
00295             dsp_y[v1] += dya;
00296           }
00297 
00298           // write new displacements back...
00299           dsp_x[v0] = dispx;
00300           dsp_y[v0] = dispy;
00301         }
00302       }
00303     } else {
00304       int e;
00305       for (e=0; e<numedges; e++) {
00306         int v0 = edges_v0[e];
00307         int v1 = edges_v1[e];
00308         float weight = weights[e];
00309 
00310         float dx = pos_x[v0] - pos_x[v1];
00311         float dy = pos_y[v0] - pos_y[v1];
00312         float dxy2 = dx*dx + dy*dy;
00313 #if defined(VMDFROPTIMIZED)
00314         float attraction = sqrtf(dxy2) * k_inv * weight;
00315         float dxa = dx * attraction;
00316         float dya = dy * attraction;
00317 #else
00318         float attraction = dxy2 * k_inv * weight;
00319         float d_1 = 1.0f / sqrtf(dxy2 + distance_epsilon);
00320         float dxa = dx * d_1 * attraction;
00321         float dya = dy * d_1 * attraction;
00322 #endif
00323 
00324         // memory scatter here is bad for performance...
00325         dsp_x[v0] -= dxa;
00326         dsp_y[v0] -= dya;
00327         dsp_x[v1] += dxa;
00328         dsp_y[v1] += dya;
00329       }
00330     }
00331 
00332     //
00333     // Integration: update graph vertex positions with
00334     // the accumulated displacement vectors, with maximum
00335     // displacement limited by the current squared "temperature", 
00336     // per a very simple simulated annealing scheme.
00337     //
00338     for (v0=0; v0<numverts; v0++) {
00339       float dx = dsp_x[v0];
00340       float dy = dsp_y[v0];
00341 
00342       float dxy2 = dx*dx + dy*dy;
00343 
00344       // limit maximum displacement according to current 
00345       // temperature factor, by rescaling the displacement magnitude 
00346       // to match the squared temperature magnitude.
00347       if (dxy2 > tempsquared * tempsquared) {
00348         float d = sqrtf(dxy2);
00349         float rescale = tempsquared / d;
00350         dx *= rescale;
00351         dy *= rescale;
00352       }
00353 
00354       pos_x[v0] += dx;
00355       pos_y[v0] += dy;
00356     } 
00357   }
00358 printf("\n");
00359 
00360 }
00361 
00362 
00363 

Generated on Wed Dec 4 02:43:45 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002