1/* $Id$ $Revision$ */
2/* vim:set shiftwidth=4 ts=8: */
3
4/*************************************************************************
5 * Copyright (c) 2011 AT&T Intellectual Property
6 * All rights reserved. This program and the accompanying materials
7 * are made available under the terms of the Eclipse Public License v1.0
8 * which accompanies this distribution, and is available at
9 * http://www.eclipse.org/legal/epl-v10.html
10 *
11 * Contributors: See CVS logs. Details at http://www.graphviz.org/
12 *************************************************************************/
13
14
15#include "bfs.h"
16#include "dijkstra.h"
17#include "kkutils.h"
18#include <stdlib.h>
19#include <math.h>
20
21int common_neighbors(vtx_data * graph, int v, int u, int *v_vector)
22{
23 /* count number of common neighbors of 'v' and 'u' */
24 int neighbor;
25 int num_shared_neighbors = 0;
26 int j;
27 for (j = 1; j < graph[u].nedges; j++) {
28 neighbor = graph[u].edges[j];
29 if (v_vector[neighbor] > 0) {
30 /* a shared neighobr */
31 num_shared_neighbors++;
32 }
33 }
34 return num_shared_neighbors;
35}
36void fill_neighbors_vec_unweighted(vtx_data * graph, int vtx, int *vtx_vec)
37{
38 /* a node is NOT a neighbor of itself! */
39 /* unlike the other version of this function */
40 int j;
41 for (j = 1; j < graph[vtx].nedges; j++) {
42 vtx_vec[graph[vtx].edges[j]] = 1;
43 }
44}
45
46void empty_neighbors_vec(vtx_data * graph, int vtx, int *vtx_vec)
47{
48 int j;
49 /* a node is NOT a neighbor of itself! */
50 /* unlike the other version ofthis function */
51 for (j = 1; j < graph[vtx].nedges; j++) {
52 vtx_vec[graph[vtx].edges[j]] = 0;
53 }
54}
55
56/* compute_apsp_dijkstra:
57 * Assumes the graph has weights
58 */
59static DistType **compute_apsp_dijkstra(vtx_data * graph, int n)
60{
61 int i;
62 DistType *storage;
63 DistType **dij;
64
65 storage = N_GNEW(n * n, DistType);
66 dij = N_GNEW(n, DistType *);
67 for (i = 0; i < n; i++)
68 dij[i] = storage + i * n;
69
70 for (i = 0; i < n; i++) {
71 dijkstra(i, graph, n, dij[i]);
72 }
73 return dij;
74}
75
76static DistType **compute_apsp_simple(vtx_data * graph, int n)
77{
78 /* compute all pairs shortest path */
79 /* for unweighted graph */
80 int i;
81 DistType *storage = N_GNEW(n * n, int);
82 DistType **dij;
83 Queue Q;
84
85 dij = N_GNEW(n, DistType *);
86 for (i = 0; i < n; i++) {
87 dij[i] = storage + i * n;
88 }
89 mkQueue(&Q, n);
90 for (i = 0; i < n; i++) {
91 bfs(i, graph, n, dij[i], &Q);
92 }
93 freeQueue(&Q);
94 return dij;
95}
96
97DistType **compute_apsp(vtx_data * graph, int n)
98{
99 if (graph->ewgts)
100 return compute_apsp_dijkstra(graph, n);
101 else
102 return compute_apsp_simple(graph, n);
103}
104
105DistType **compute_apsp_artifical_weights(vtx_data * graph, int n)
106{
107 DistType **Dij;
108 /* compute all-pairs-shortest-path-length while weighting the graph */
109 /* so high-degree nodes are distantly located */
110
111 float *old_weights = graph[0].ewgts;
112
113 compute_new_weights(graph, n);
114 Dij = compute_apsp_dijkstra(graph, n);
115 restore_old_weights(graph, n, old_weights);
116 return Dij;
117}
118
119
120/**********************/
121/* */
122/* Quick Sort */
123/* */
124/**********************/
125
126static void
127split_by_place(double *place, int *nodes, int first, int last, int *middle)
128{
129 unsigned int splitter=((unsigned int)rand()|((unsigned int)rand())<<16)%(unsigned int)(last-first+1)+(unsigned int)first;
130 int val;
131 double place_val;
132 int left = first + 1;
133 int right = last;
134 int temp;
135
136 val = nodes[splitter];
137 nodes[splitter] = nodes[first];
138 nodes[first] = val;
139 place_val = place[val];
140
141 while (left < right) {
142 while (left < right && place[nodes[left]] <= place_val)
143 left++;
144 /* use here ">" and not ">=" to enable robustness
145 * by ensuring that ALL equal values move to the same side
146 */
147 while (left < right && place[nodes[right]] > place_val)
148 right--;
149 if (left < right) {
150 temp = nodes[left];
151 nodes[left] = nodes[right];
152 nodes[right] = temp;
153 left++;
154 right--; /* (1) */
155
156 }
157 }
158 /* at this point either, left==right (meeting), or
159 * left=right+1 (because of (1))
160 * we have to decide to which part the meeting point (or left) belongs.
161 *
162 * notice that always left>first, because of its initialization
163 */
164 if (place[nodes[left]] > place_val)
165 left = left - 1;
166 *middle = left;
167 nodes[first] = nodes[left];
168 nodes[left] = val;
169}
170
171double distance_kD(double **coords, int dim, int i, int j)
172{
173 /* compute a k-D Euclidean distance between 'coords[*][i]' and 'coords[*][j]' */
174 double sum = 0;
175 int k;
176 for (k = 0; k < dim; k++) {
177 sum +=
178 (coords[k][i] - coords[k][j]) * (coords[k][i] - coords[k][j]);
179 }
180 return sqrt(sum);
181}
182
183static float* fvals;
184static int
185fcmpf (int* ip1, int* ip2)
186{
187 float d1 = fvals[*ip1];
188 float d2 = fvals[*ip2];
189 if (d1 < d2) return -1;
190 else if (d1 > d2) return 1;
191 else return 0;
192}
193
194void quicksort_placef(float *place, int *ordering, int first, int last)
195{
196 if (first < last) {
197 fvals = place;
198 qsort(ordering+first, last-first+1, sizeof(ordering[0]), (qsort_cmpf)fcmpf);
199 }
200}
201
202static int
203sorted_place(double * place, int * ordering, int first,int last)
204{
205 int i, isSorted = 1;
206 for (i=first+1; i<=last && isSorted; i++) {
207 if (place[ordering[i-1]]>place[ordering[i]]) {
208 isSorted = 0;
209 }
210 }
211 return isSorted;
212}
213
214/* quicksort_place:
215 * For now, we keep the current implementation for stability, but
216 * we should consider replacing this with an implementation similar to
217 * quicksort_placef above.
218 */
219void quicksort_place(double *place, int *ordering, int first, int last)
220{
221 if (first < last) {
222 int middle;
223#ifdef __cplusplus
224 split_by_place(place, ordering, first, last, middle);
225#else
226 split_by_place(place, ordering, first, last, &middle);
227#endif
228 quicksort_place(place, ordering, first, middle - 1);
229 quicksort_place(place, ordering, middle + 1, last);
230 /* Checking for "already sorted" dramatically improves running time
231 * and robustness (against uneven recursion) when not all values are
232 * distinct (thefore we expect emerging chunks of equal values)
233 * it never increased running time even when values were distinct
234 */
235 if (!sorted_place(place,ordering,first,middle-1))
236 quicksort_place(place,ordering,first,middle-1);
237 if (!sorted_place(place,ordering,middle+1,last))
238 quicksort_place(place,ordering,middle+1,last);
239 }
240}
241
242void compute_new_weights(vtx_data * graph, int n)
243{
244 /* Reweight graph so that high degree nodes will be separated */
245
246 int i, j;
247 int nedges = 0;
248 float *weights;
249 int *vtx_vec = N_GNEW(n, int);
250 int deg_i, deg_j, neighbor;
251
252 for (i = 0; i < n; i++) {
253 nedges += graph[i].nedges;
254 }
255 weights = N_GNEW(nedges, float);
256
257 for (i = 0; i < n; i++) {
258 vtx_vec[i] = 0;
259 }
260
261 for (i = 0; i < n; i++) {
262 graph[i].ewgts = weights;
263 fill_neighbors_vec_unweighted(graph, i, vtx_vec);
264 deg_i = graph[i].nedges - 1;
265 for (j = 1; j <= deg_i; j++) {
266 neighbor = graph[i].edges[j];
267 deg_j = graph[neighbor].nedges - 1;
268 weights[j] =
269 (float) (deg_i + deg_j -
270 2 * common_neighbors(graph, i, neighbor,
271 vtx_vec));
272 }
273 empty_neighbors_vec(graph, i, vtx_vec);
274 weights += graph[i].nedges;
275 }
276 free(vtx_vec);
277}
278
279void restore_old_weights(vtx_data * graph, int n, float *old_weights)
280{
281 int i;
282 free(graph[0].ewgts);
283 graph[0].ewgts = NULL;
284 if (old_weights != NULL) {
285 for (i = 0; i < n; i++) {
286 graph[i].ewgts = old_weights;
287 old_weights += graph[i].nedges;
288 }
289 }
290}
291