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 //  FOR_VERTEX(v) {
  40 //
  41 }
  42
  43 void graph_layout_prepare() {
  44   int v1;
  45
  46   FOR_VERTEX(v1)
  47     breadth_first_search(v1, sqdistances[v1]);
  48
  49   alpha= 2;
  50   beta= log(10)/log(alpha);
  51   beta_prime= (1-beta)/2;
  52   printf("alpha=%g beta=%g beta'=%g\n", alpha,beta,beta_prime);
  53 }
  54
  55
  56 double graph_layout_cost(const Vertices v) {
  57   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
  58    * along edges, and actual distance d_ij = |vi-vj|.
  59    *
  60    * We will also use the `vertex areas': for each vertex vi the
  61    * vertex area a_vi is the mean area of the incident triangles.
  62    * This is computed elsewhere.
  63    *
  64    * Energy contribution is proportional to
  65    *
  66    *               -4          2
  67    *    a  a   .  d   . [ (s/d)  - 1 ]
  68    *     vi vj
  69    *
  70    * (In practice we compute d^2+epsilon and use it for the
  71    *  divisions, to avoid division by zero.)
  72    */
  73   //static const double d2_epsilon= 1e-6;
  74
  75   double total_cost=0;
  76   int v1,v2,e, nedges=0;
  77   double totaledgelength=0, meanedgelength, meanedgelength2;
  78
  79   FOR_EDGE(v1,e,v2) {
  80     totaledgelength += hypotD(v[v1], v[v2]);
  81     nedges++;
  82   }
  83
  84   meanedgelength= totaledgelength / nedges;
  85   meanedgelength2= meanedgelength * meanedgelength;
  86 //  printf("mean=%g mean^2=%g\n", meanedgelength, meanedgelength2);
  87
  88   FOR_VERTEX(v1) {
  89     FOR_VERTEX(v2) {
  90       if (v1 == v2) continue;
  91
  92       double d2= hypotD2(v[v1],v[v2]);
  93
  94       int dist2= sqdistances[v1][v2];
  95       assert(dist2>0);
  96
  97       double s2= dist2 * meanedgelength2;
  98
  99       /* energy = (d/s)^(1-beta)     where beta is log\_{alpha}(10)
 100        * energy = ((d/s)^2) ^ (1-beta)/2
 101        *          let beta' = (1-beta)/2
 102        */
 103
 104       double cost= (vertex_mean_edge_lengths[v1]) * pow(d2/s2, beta_prime);
 105
 106       //printf("layout %03x..%03x dist^2=%d s^2=%g d^2=%g "
 107                //" cost+=%g\n", v1,v2, dist2,
 108                //            s2,d2, cost);
 109       total_cost += cost;
 110     }
 111   }
 112   return total_cost/meanedgelength;
 113 }