graph.c

   1 /*
   2  * graph layout energy
   3  */
   4
   5 #include "mgraph.h"
   6 #include "minimise.h"
   7
   8 static int sqdistances[N][N];
   9
  10 static double alpha, beta, beta_prime;
  11
  12 static void breadth_first_search(int start, int sqdistances_r[N]) {
  13   int d[N], buffer[N], *buf_pop=buffer, *buf_push=buffer;
  14   int v,e, current, future, dfuture;
  15
  16   buf_push= buf_pop= buffer;
  17   FOR_VERTEX(v) d[v]= -1;
  18
  19   d[start]= 0;
  20   *buf_push++= start;
  21
  22   while (buf_pop < buf_push) {
  23     current= *buf_pop++;
  24     dfuture= d[current] + 1;
  25     FOR_VEDGE(current,e,future) {
  26       if (d[future] >= 0) continue; /* already found this one */
  27       d[future]= dfuture;
  28       *buf_push++= future;
  29     }
  30   }
  31   assert(buf_pop==buf_push);
  32   assert(buf_push <= buffer+sizeof(buffer)/sizeof(buffer[0]));
  33
  34   FOR_VERTEX(v) {
  35     assert(d[v] >= 0);
  36     sqdistances_r[v]= d[v] * d[v];
  37   }
  38 }
  39
  40 void graph_layout_prepare() {
  41   int v1;
  42
  43   FOR_VERTEX(v1)
  44     breadth_first_search(v1, sqdistances[v1]);
  45
  46   alpha= 2;
  47   beta= -log(10)/log(alpha);
  48   beta_prime= (1-beta)/2;
  49   printf("alpha=%g beta=%g beta'=%g\n", alpha,beta,beta_prime);
  50 }
  51
  52
  53 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
  54   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
  55    * along edges, and actual distance d_ij = |vi-vj|.
  56    *
  57    * We will also use the `vertex areas': for each vertex vi the
  58    * vertex area a_vi is the mean area of the incident triangles.
  59    * This is computed elsewhere.
  60    *
  61    * Energy contribution is proportional to
  62    *
  63    *               -4          2
  64    *    a  a   .  d   . [ (s/d)  - 1 ]
  65    *     vi vj
  66    *
  67    * (In practice we compute d^2+epsilon and use it for the
  68    *  divisions, to avoid division by zero.)
  69    */
  70   //static const double d2_epsilon= 1e-6;
  71
  72   //  double edge_weights[V6<<ESHIFT], vertex_distances[N],
  73   double total_cost=0;
  74   int v1,v2,e, nedges=0;
  75   double totaledgelength=0, meanedgelength, meanedgelength2;
  76
  77   FOR_EDGE(v1,e,v2) {
  78     totaledgelength += hypotD(v[v1], v[v2]);
  79     nedges++;
  80   }
  81
  82   meanedgelength= totaledgelength / nedges;
  83   meanedgelength2= meanedgelength * meanedgelength;
  84 //  printf("mean=%g mean^2=%g\n", meanedgelength, meanedgelength2);
  85
  86   FOR_VERTEX(v1) {
  87     FOR_VERTEX(v2) {
  88       if (v1 == v2) continue;
  89
  90       double d2= hypotD2(v[v1],v[v2]);
  91
  92       int dist2= sqdistances[v1][v2];
  93       assert(dist2>0);
  94
  95       double s2= dist2 * meanedgelength2;
  96
  97       /* energy = (d/s)^(1-beta)     where beta is -log\_{alpha}(10)
  98        * energy = ((d/s)^2) ^ (1-beta)/2
  99        *          let beta' = (1-beta)/2
 100        */
 101
 102       double cost= pow(d2/s2, beta_prime);
 103
 104       //double s2= s*s + d2_epsilon;
 105       //double sd2= s2 / d2;
 106       //double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
 107       //double cost_contrib= sd2;
 108
 109       //printf("layout %03x..%03x dist^2=%d s^2=%g d^2=%g "
 110                //" cost+=%g\n", v1,v2, dist2,
 111                //            s2,d2, cost);
 112       total_cost += cost;
 113     }
 114   }
 115   return total_cost;
 116 }