graph.c

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