chiark / gitweb /
set threads_started
[moebius2.git] / 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 }