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 "kkutils.h"
16#include "closest.h"
17#include <stdlib.h>
18
19/*****************************************
20** This module contains functions that **
21** given a 1-D layout construct a graph **
22** where close nodes in this layout are **
23** adjacent **
24*****************************************/
25
26typedef struct {
27 /* this structure represents two nodes in the 1-D layout */
28 int left; /* the left node in the pair */
29 int right; /* the right node in the pair */
30 double dist; /* distance between the nodes in the layout */
31} Pair;
32
33#define LT(p,q) ((p).dist < (q).dist)
34#define EQ(p,q) ((p).dist == (q).dist)
35
36/*
37Pair(int v, int u) {left=v; right=u;}
38bool operator>(Pair other) {return dist>other.dist;}
39bool operator>=(Pair other) {return dist>=other.dist;}
40bool operator<(Pair other) {return dist<other.dist;}
41bool operator<=(Pair other) {return dist<=other.dist;}
42bool operator==(Pair other) {return dist==other.dist;}
43*/
44
45typedef struct {
46 Pair *data;
47 int max_size;
48 int top;
49} PairStack;
50
51static void initStack(PairStack * s, int n)
52{
53 s->data = N_GNEW(n, Pair);
54 s->max_size = n;
55 s->top = 0;
56}
57
58static void freeStack(PairStack * s)
59{
60 free(s->data);
61}
62
63#define push(s,x) { \
64 if (s->top>=s->max_size) { \
65 s->max_size *= 2; \
66 s->data = (Pair*) realloc(s->data, s->max_size*sizeof(Pair)); \
67 } \
68 s->data[s->top++] = x; \
69}
70
71#define pop(s,x) ((s->top==0) ? FALSE : (s->top--, x = s->data[s->top], TRUE))
72
73#define read_top(h,x) ((s->top==0) ? FALSE : (x = s->data[s->top-1], TRUE))
74
75#define sub(h,i) (h->data[i])
76
77/* An auxulliary data structure (a heap) for
78 * finding the closest pair in the layout
79 */
80typedef struct {
81 Pair *data;
82 int heapSize;
83 int maxSize;
84} PairHeap;
85
86#define left(i) (2*(i))
87#define right(i) (2*(i)+1)
88#define parent(i) ((i)/2)
89#define insideHeap(h,i) ((i)<h->heapSize)
90#define greaterPriority(h,i,j) \
91 (LT(h->data[i],h->data[j]) || ((EQ(h->data[i],h->data[j])) && (rand()%2)))
92
93#define exchange(h,i,j) {Pair temp; \
94 temp=h->data[i]; \
95 h->data[i]=h->data[j]; \
96 h->data[j]=temp; \
97}
98#define assign(h,i,j) {h->data[i]=h->data[j]}
99
100static void heapify(PairHeap * h, int i)
101{
102 int l, r, largest;
103 while (1) {
104 l = left(i);
105 r = right(i);
106 if (insideHeap(h, l) && greaterPriority(h, l, i))
107 largest = l;
108 else
109 largest = i;
110 if (insideHeap(h, r) && greaterPriority(h, r, largest))
111 largest = r;
112 if (largest == i)
113 break;
114
115 exchange(h, largest, i);
116 i = largest;
117 }
118}
119
120#ifdef UNUSED
121static void mkHeap(PairHeap * h, int size)
122{
123 h->data = N_GNEW(size, Pair);
124 h->maxSize = size;
125 h->heapSize = 0;
126}
127#endif
128
129static void freeHeap(PairHeap * h)
130{
131 free(h->data);
132}
133
134static void initHeap(PairHeap * h, double *place, int *ordering, int n)
135{
136 int i;
137 Pair edge;
138 int j;
139
140 h->heapSize = n - 1;
141#ifdef REDO
142 if (h->heapSize > h->maxSize) {
143 h->maxSize = h->heapSize;
144 h->data = (Pair *) realloc(h->data, h->maxSize * sizeof(Pair));
145 }
146#else
147 h->maxSize = h->heapSize;
148 h->data = N_GNEW(h->maxSize, Pair);
149#endif
150
151 for (i = 0; i < n - 1; i++) {
152 edge.left = ordering[i];
153 edge.right = ordering[i + 1];
154 edge.dist = place[ordering[i + 1]] - place[ordering[i]];
155 h->data[i] = edge;
156 }
157 for (j = (n - 1) / 2; j >= 0; j--) {
158 heapify(h, j);
159 }
160}
161
162static boolean extractMax(PairHeap * h, Pair * max)
163{
164 if (h->heapSize == 0)
165 return FALSE;
166
167 *max = h->data[0];
168 h->data[0] = h->data[h->heapSize - 1];
169 h->heapSize--;
170 heapify(h, 0);
171 return TRUE;
172}
173
174static void insert(PairHeap * h, Pair edge)
175{
176 int i = h->heapSize;
177 if (h->heapSize == h->maxSize) {
178 h->maxSize *= 2;
179 h->data = (Pair *) realloc(h->data, h->maxSize * sizeof(Pair));
180 }
181 h->heapSize++;
182 h->data[i] = edge;
183 while (i > 0 && greaterPriority(h, i, parent(i))) {
184 exchange(h, i, parent(i));
185 i = parent(i);
186 }
187}
188
189/*
190static bool
191isheap(PairHeap* h)
192{
193 int i,l,r;
194 for (i=0; i<h->heapSize; i++) {
195 l=left(i); r=right(i);
196 if (insideHeap(h,l) && greaterPriority(h,l,i))
197 return FALSE;
198 if (insideHeap(h,r) && greaterPriority(h,r,i))
199 return FALSE;
200 }
201 return TRUE;
202}
203*/
204
205static void
206find_closest_pairs(double *place, int n, int num_pairs,
207 PairStack * pairs_stack)
208{
209 /* Fill the stack 'pairs_stack' with 'num_pairs' closest pairs int the 1-D layout 'place' */
210 int i;
211 PairHeap heap;
212 int *left = N_GNEW(n, int);
213 int *right = N_GNEW(n, int);
214 Pair pair = { 0, 0 }, new_pair;
215
216 /* Order the nodes according to their place */
217 int *ordering = N_GNEW(n, int);
218 int *inv_ordering = N_GNEW(n, int);
219
220 for (i = 0; i < n; i++) {
221 ordering[i] = i;
222 }
223 quicksort_place(place, ordering, 0, n - 1);
224 for (i = 0; i < n; i++) {
225 inv_ordering[ordering[i]] = i;
226 }
227
228 /* Initialize heap with all consecutive pairs */
229 initHeap(&heap, place, ordering, n);
230
231 /* store the leftmost and rightmost neighbors of each node that were entered into heap */
232 for (i = 1; i < n; i++) {
233 left[ordering[i]] = ordering[i - 1];
234 }
235 for (i = 0; i < n - 1; i++) {
236 right[ordering[i]] = ordering[i + 1];
237 }
238
239 /* extract the 'num_pairs' closest pairs */
240 for (i = 0; i < num_pairs; i++) {
241 int left_index;
242 int right_index;
243 int neighbor;
244
245 if (!extractMax(&heap, &pair)) {
246 break; /* not enough pairs */
247 }
248 push(pairs_stack, pair);
249 /* insert to heap "descendant" pairs */
250 left_index = inv_ordering[pair.left];
251 right_index = inv_ordering[pair.right];
252 if (left_index > 0) {
253 neighbor = ordering[left_index - 1];
254 if (inv_ordering[right[neighbor]] < right_index) {
255 /* we have a new pair */
256 new_pair.left = neighbor;
257 new_pair.right = pair.right;
258 new_pair.dist = place[pair.right] - place[neighbor];
259 insert(&heap, new_pair);
260 right[neighbor] = pair.right;
261 left[pair.right] = neighbor;
262 }
263 }
264 if (right_index < n - 1) {
265 neighbor = ordering[right_index + 1];
266 if (inv_ordering[left[neighbor]] > left_index) {
267 /* we have a new pair */
268 new_pair.left = pair.left;
269 new_pair.right = neighbor;
270 new_pair.dist = place[neighbor] - place[pair.left];
271 insert(&heap, new_pair);
272 left[neighbor] = pair.left;
273 right[pair.left] = neighbor;
274 }
275 }
276 }
277 free(left);
278 free(right);
279 free(ordering);
280 free(inv_ordering);
281 freeHeap(&heap);
282}
283
284static void add_edge(vtx_data * graph, int u, int v)
285{
286 int i;
287 for (i = 0; i < graph[u].nedges; i++) {
288 if (graph[u].edges[i] == v) {
289 /* edge already exist */
290 return;
291 }
292 }
293 /* add the edge */
294 graph[u].edges[graph[u].nedges++] = v;
295 graph[v].edges[graph[v].nedges++] = u;
296 if (graph[0].ewgts != NULL) {
297 graph[u].ewgts[0]--;
298 graph[v].ewgts[0]--;
299 }
300}
301
302static void
303construct_graph(int n, PairStack * edges_stack, vtx_data ** New_graph)
304{
305 /* construct an unweighted graph using the edges 'edges_stack' */
306 int i;
307 vtx_data *new_graph;
308
309 /* first compute new degrees and nedges; */
310 int *degrees = N_GNEW(n, int);
311 int top = edges_stack->top;
312 int new_nedges = 2 * top + n;
313 Pair pair;
314 int *edges = N_GNEW(new_nedges, int);
315 float *weights = N_GNEW(new_nedges, float);
316
317 for (i = 0; i < n; i++) {
318 degrees[i] = 1; /* save place for the self loop */
319 }
320 for (i = 0; i < top; i++) {
321 pair = sub(edges_stack, i);
322 degrees[pair.left]++;
323 degrees[pair.right]++;
324 }
325
326 /* copy graph into new_graph: */
327 for (i = 0; i < new_nedges; i++) {
328 weights[i] = 1.0;
329 }
330
331 *New_graph = new_graph = N_GNEW(n, vtx_data);
332 for (i = 0; i < n; i++) {
333 new_graph[i].nedges = 1;
334 new_graph[i].ewgts = weights;
335#ifdef USE_STYLES
336 new_graph[i].styles = NULL;
337#endif
338 new_graph[i].edges = edges;
339 *edges = i; /* self loop for Lap */
340 *weights = 0; /* self loop weight for Lap */
341 weights += degrees[i];
342 edges += degrees[i]; /* reserve space for possible more edges */
343 }
344
345 free(degrees);
346
347 /* add all edges from stack */
348 while (pop(edges_stack, pair)) {
349 add_edge(new_graph, pair.left, pair.right);
350 }
351}
352
353void
354closest_pairs2graph(double *place, int n, int num_pairs, vtx_data ** graph)
355{
356 /* build a graph with with edges between the 'num_pairs' closest pairs in the 1-D space: 'place' */
357 PairStack pairs_stack;
358 initStack(&pairs_stack, num_pairs);
359 find_closest_pairs(place, n, num_pairs, &pairs_stack);
360 construct_graph(n, &pairs_stack, graph);
361 freeStack(&pairs_stack);
362}
363