1 /****************************************************************************/
2 /* */
3 /* This file is part of CONCORDE */
4 /* */
5 /* (c) Copyright 1995--1999 by David Applegate, Robert Bixby, */
6 /* Vasek Chvatal, and William Cook */
7 /* */
8 /* Permission is granted for academic research use. For other uses, */
9 /* contact the authors for licensing options. */
10 /* */
11 /* Use at your own risk. We make no guarantees about the */
12 /* correctness or usefulness of this code. */
13 /* */
14 /****************************************************************************/
15
16 /****************************************************************************/
17 /* */
18 /* NEAREST NEIGHBORS FOR X-NORMS AND JUNK-NORMS */
19 /* */
20 /* TSP CODE */
21 /* */
22 /* */
23 /* Written by: Applegate, Bixby, Chvatal, and Cook */
24 /* Date: March 2, 1995 */
25 /* */
26 /* */
27 /* EXPORTED FUNCTIONS: */
28 /* */
29 /* int CCedgegen_x_k_nearest (int ncount, int k, CCdatagroup *dat, */
30 /* double *wcoord, int wantlist, int *ecount, int **elist, */
31 /* int silent) */
32 /* RETURNS the k_nearest neighbor graph (for X-Norms) */
33 /* -ncount is the number of nodes */
34 /* -k is the number of nearest neighbors wanted */
35 /* -dat contains the info to generate edge lengths */
36 /* -wcoord are nodeweights for Held-Karp style edge lengths, using */
37 /* len[i,j] + wcoord[i] + wcoord[j] (wcoord can be NULL) */
38 /* -wantlist should be set to 0 if you don't want the edges */
39 /* -ecount returns the number of edges if wantlist is 1 */
40 /* -elist returns the edges in end1 end2 format if wantlist is 1 */
41 /* */
42 /* int CCedgegen_x_quadrant_k_nearest (int ncount, int k, */
43 /* CCdatagroup *dat, double *wcoord, int wantlist, int *ecount, */
44 /* int **elist, int silent) */
45 /* RETURNS the quadrant k_nearest_graph (for X-Norms) */
46 /* */
47 /* int CCedgegen_x_node_k_nearest (CCxnear *xn, int n, int k, */
48 /* int ncount, int *list) */
49 /* RETURNS the k nearest neighbors from node n (for X-Norms */
50 /* -xn is a structure built by a call to CCedgegen_xnear_build () */
51 /* -list returns the neighbors of n. The calling routine should */
52 /* be sure that list points to an array of length at least num. */
53 /* */
54 /* int CCedgegen_x_node_quadrant_k_nearest (CCxnear *xn, int n, int k, */
55 /* int ncount, int *list) */
56 /* RETURNS the quadrant k nearest to node n (for X-Norms) */
57 /* -xn is a structure built by a call to CCedgegen_xnear_build () */
58 /* -list returns the neighbors of n. The calling routine should */
59 /* be sure that list points to a sufficiently large array (4*num */
60 /* for D2_SIZE norms and 8*num for D3_SIZE norms) */
61 /* */
62 /* int CCedgegen_x_node_nearest (CCxnear *xn, int ncount, int ni, */
63 /* char *marks) */
64 /* RETURNS the nearest unmarked node to node n (as the return value) */
65 /* -marks is an array. The entries that are nonzero correspond to */
66 /* nodes that will not be looked at in the search. */
67 /* */
68 /* int CCedgegen_x_nearest_neighbor_tour (int ncount, int start, */
69 /* CCdatagroup *dat, int *outcycle, double *val) */
70 /* RETURNS a nearest neighbor tour, starting at node start. */
71 /* -outcycle will contain the tour if it is not NULL (the calling */
72 /* routine should be sure it points to an array of length at */
73 /* least ncount if it is not set to NULL) */
74 /* -val will return the length of the tour. */
75 /* */
76 /* int CCedgegen_junk_k_nearest (int ncount, int k, CCdatagroup *dat, */
77 /* double *wcoord, int wantlist, int *ecount, int **elist, */
78 /* int silent) */
79 /* RETURNS the k-nearest graph (for JUNK-Norms) */
80 /* -see CCedgegen_x_k_nearest (above) for the variables */
81 /* */
82 /* int CCedgegen_junk_node_k_nearest (CCdatagroup *dat, double *wcoord, */
83 /* int n, int k, int ncount, int *list) */
84 /* RETURNS the k nearest neighbors to node n (for JUNK-Norms) */
85 /* -list returns the neighbors of n. The calling routine should */
86 /* be sure that list points to an array of length at least num. */
87 /* */
88 /* int CCedgegen_junk_node_nearest (CCdatagroup *dat, double *wcoord, */
89 /* int ncount, int n, char *marks) */
90 /* RETURNS the nearest unmarked node to node n (as the return value) */
91 /* -marks is an array, the nodes with marks[i] nonzero are ignored. */
92 /* */
93 /* int CCedgegen_junk_nearest_neighbor_tour (int ncount, int start, */
94 /* CCdatagroup *dat, int *outcycle, double *val, int silent) */
95 /* RETURNS a nearest neighbor tour starting at node start. Note that */
96 /* this will be slow for large problems (it is a quadratic routine) */
97 /* -see the describtion of CCedgegen_x_nearest_neighbor_tour above */
98 /* */
99 /* int CCedgegen_junk_greedy_tour (int ncount, CCdatagroup *dat, */
100 /* int *outcycle, double *val, int ecount, int *elist, int silent) */
101 /* RETURNS a greedy tour using edges from elist. Disjoint segments */
102 /* from elist will be connected by nearest neighbors (which are slow */
103 /* for large problems). */
104 /* */
105 /* int CCedgegen_junk_qboruvka_tour (int ncount, CCdatagroup *dat, */
106 /* int *outcycle, double *val, int ecount, int *elist, int silent) */
107 /* RETURNS a quick-boruvka tour using edges from elist. Disjoint */
108 /* segments from elist will be connected by nearest neighbors (which */
109 /* are slow for large problems). See the description of */
110 /* CCkdtree_qboruvka_tour for a description of quick-boruvka. */
111 /* */
112 /* int CCedgegen_xnear_build (int ncount, CCdatagroup *dat, */
113 /* double *wcoord, CCxnear *xn) */
114 /* RETURNS the structure needed for calls to */
115 /* CCedgegen_x_node_k_nearest and */
116 /* CCedgegen_x_node_quadrant_k_nearest (the calling routine should */
117 /* be sure that xn points to such a structure). All this routine */
118 /* does is permute the data so that the x coordinates are in */
119 /* nonincreasing order. */
120 /* */
121 /* void CCedgegen_xnear_free (CCxnear *xn) */
122 /* FREES the CCxnear structure pointed to by xn */
123 /* */
124 /* NOTES: */
125 /* All routines other than CCedgegen_xnear_free and */
126 /* CCedgegen_x_node_nearest return 0 on success and 1 on failure */
127 /* (normally due to running out of memory). */
128 /* The X-Norm functions will also work for KD-Norms, but they are */
129 /* much slower than the KD-Norm functions. */
130 /* */
131 /****************************************************************************/
132
133 #include "machdefs.h"
134 #include "edgegen.h"
135 #include "util.h"
136 #include "macrorus.h"
137
138 #define BIGDOUBLE (1e30)
139 #define NEAR_HEAP_CUTOFF 100 /* When to switch from list to heap */
140
141 #define dtrunc(x) (((x)>0.0)?floor(x):ceil(x))
142
143 typedef struct shortedge {
144 double length;
145 int end;
146 } shortedge;
147
148 typedef struct intptr {
149 int this;
150 struct intptr *next;
151 } intptr;
152
153 typedef struct tabledat {
154 intptr **table;
155 CCptrworld intptr_world;
156 } tabledat;
157
158 static void
159 add_to_list_and_reset (int *list, int *lcount, shortedge *nearlist,
160 int nearnum, int *nodenames),
161 insert (int n, int m, shortedge *nearlist, CCdatagroup *dat,
162 double *wcoord),
163 x_quicksort (int *list, double *x, int l, int u);
164 static int
165 run_x_k_nearest (int ncount, int num, CCdatagroup *dat, double *wcoord,
166 int wantlist, int *ecount, int **elist, int doquad, int silent),
167 put_in_table (tabledat *td, int i, int j, int *added);
168
169
170
CC_PTRWORLD_LIST_ROUTINES(intptr,int,intptralloc,intptr_bulk_alloc,intptrfree,intptr_listadd,intptr_listfree)171 CC_PTRWORLD_LIST_ROUTINES (intptr, int, intptralloc, intptr_bulk_alloc,
172 intptrfree, intptr_listadd, intptr_listfree)
173 CC_PTRWORLD_LEAKS_ROUTINE (intptr, intptr_check_leaks, this, int)
174
175
176 int CCedgegen_x_k_nearest (int ncount, int num, CCdatagroup *dat,
177 double *wcoord, int wantlist, int *ecount, int **elist, int silent)
178 {
179 return run_x_k_nearest (ncount, num, dat, wcoord, wantlist, ecount,
180 elist, 0, silent);
181 }
182
CCedgegen_x_quadrant_k_nearest(int ncount,int num,CCdatagroup * dat,double * wcoord,int wantlist,int * ecount,int ** elist,int silent)183 int CCedgegen_x_quadrant_k_nearest (int ncount, int num, CCdatagroup *dat,
184 double *wcoord, int wantlist, int *ecount, int **elist, int silent)
185 {
186 return run_x_k_nearest (ncount, num, dat, wcoord, wantlist, ecount,
187 elist, 1, silent);
188 }
189
CCedgegen_junk_k_nearest(int ncount,int num,CCdatagroup * dat,double * wcoord,int wantlist,int * ecount,int ** elist,int silent)190 int CCedgegen_junk_k_nearest (int ncount, int num, CCdatagroup *dat,
191 double *wcoord, int wantlist, int *ecount, int **elist, int silent)
192 {
193 return run_x_k_nearest (ncount, num, dat, wcoord, wantlist, ecount,
194 elist, 0, silent);
195 }
196
197
run_x_k_nearest(int ncount,int num,CCdatagroup * dat,double * wcoord,int wantlist,int * ecount,int ** elist,int doquad,int silent)198 static int run_x_k_nearest (int ncount, int num, CCdatagroup *dat,
199 double *wcoord, int wantlist, int *ecount, int **elist, int doquad,
200 int silent)
201 {
202 int rval = 0;
203 int i, n;
204 intptr *ip, *ipnext;
205 int total, onlist;
206 int added, ntotal = 0;
207 int *list = (int *) NULL;
208 int goal, usex;
209 CCxnear xn;
210 int norm;
211 tabledat td;
212
213 td.table = (intptr **) NULL;
214 CCptrworld_init (&td.intptr_world);
215
216 xn.nodenames = (int *) NULL;
217 xn.invnames = (int *) NULL;
218 xn.w = (double *) NULL;
219 CCutil_init_datagroup (&xn.dat);
220
221 if (wantlist) {
222 *ecount = 0;
223 *elist = (int *) NULL;
224 }
225
226 CCutil_dat_getnorm (dat, &norm);
227 if ((norm & CC_NORM_BITS) == CC_X_NORM_TYPE ||
228 (norm & CC_NORM_BITS) == CC_KD_NORM_TYPE) {
229 usex = 1;
230 if (!silent) {
231 printf ("Using x-norm nearest code\n"); fflush (stdout);
232 }
233 } else {
234 usex = 0;
235 if (!silent) {
236 printf ("Using junk-norm nearest code\n"); fflush (stdout);
237 }
238 }
239
240 if (wcoord != (double *) NULL) {
241 for (i = 0; i < ncount; i++) {
242 if (wcoord[i] < -0.00000001) {
243 fprintf (stderr, "Cannot use CCxnear with negative weights\n");
244 return 1;
245 }
246 }
247 }
248
249 if (doquad && dat->z != (double *) NULL)
250 goal = 8 * num;
251 else if (doquad)
252 goal = 4 * num;
253 else
254 goal = num;
255
256 if (usex) {
257 if (CCedgegen_xnear_build (ncount, dat, wcoord, &xn)) {
258 fprintf (stderr, "build_nodes failed\n");
259 rval = 1;
260 goto CLEANUP;
261 }
262 }
263
264 td.table = CC_SAFE_MALLOC (ncount, intptr *);
265 if (!td.table) {
266 rval = 1;
267 goto CLEANUP;
268 }
269 for (i = 0; i < ncount; i++)
270 td.table[i] = (intptr *) NULL;
271 list = CC_SAFE_MALLOC (goal, int);
272 if (!list) {
273 rval = 1;
274 goto CLEANUP;
275 }
276
277 if (!usex && doquad) {
278 printf ("NOTE: Cannot run quadrant nearest with a JUNK norm.\n");
279 printf (" Running nearest instead.\n");
280 fflush (stdout);
281 }
282
283 for (n = 0; n < ncount; n++) {
284 if (usex) {
285 if (doquad) {
286 if (CCedgegen_x_node_quadrant_k_nearest (&xn, n, num, ncount,
287 list)) {
288 fprintf (stderr, "CCedgegen_x_node_quadrant_k_nearest failed\n");
289 rval = 1;
290 goto CLEANUP;
291 }
292 } else {
293 if (CCedgegen_x_node_k_nearest (&xn, n, num, ncount, list)) {
294 fprintf (stderr, "CCedgegen_x_node_k_nearest failed\n");
295 rval = 1;
296 goto CLEANUP;
297 }
298 }
299 } else {
300 if (CCedgegen_junk_node_k_nearest (dat, wcoord, n, num, ncount,
301 list)) {
302 fprintf (stderr, "junk_node_k_nearest_failed\n");
303 rval = 1;
304 goto CLEANUP;
305 }
306 }
307 for (i = 0; i < goal; i++) {
308 if (list[i] != -1) {
309 if (put_in_table (&td, n, list[i], &added)) {
310 fprintf (stderr, "put_in_table failed\n");
311 rval = 1;
312 goto CLEANUP;
313 } else {
314 ntotal += added;
315 }
316 }
317 }
318 if (!silent) {
319 if (n % 1000 == 999) {
320 printf (".");
321 fflush (stdout);
322 }
323 }
324 }
325
326 if (!silent) {
327 printf (" %d edges\n", ntotal); fflush (stdout);
328 }
329
330 if (wantlist) {
331 int j = 0;
332 *elist = CC_SAFE_MALLOC (2 * ntotal, int);
333 if (!(*elist)) {
334 rval = 1;
335 goto CLEANUP;
336 }
337 *ecount = ntotal;
338 for (i = 0; i < ncount; i++) {
339 for (ip = td.table[i]; ip; ip = ipnext) {
340 ipnext = ip->next;
341 (*elist)[j++] = i;
342 (*elist)[j++] = ip->this;
343 intptrfree (&td.intptr_world, ip);
344 }
345 td.table[i] = (intptr *) NULL;
346 }
347 } else {
348 for (i = 0; i < ncount; i++) {
349 intptr_listfree (&td.intptr_world, td.table[i]);
350 td.table[i] = (intptr *) NULL;
351 }
352 }
353
354 if (intptr_check_leaks (&td.intptr_world, &total, &onlist)) {
355 fprintf (stderr, "WARNING: %d outstanding intptrs in kdnear\n",
356 total - onlist);
357 }
358
359 CLEANUP:
360
361 CCptrworld_delete (&td.intptr_world);
362 CC_IFFREE (list, int);
363 CC_IFFREE (td.table, intptr *);
364 if (usex)
365 CCedgegen_xnear_free (&xn);
366
367 return rval;
368 }
369
CCedgegen_x_node_quadrant_k_nearest(CCxnear * xn,int ni,int nearnum,int ncount,int * list)370 int CCedgegen_x_node_quadrant_k_nearest (CCxnear *xn, int ni, int nearnum,
371 int ncount, int *list)
372 {
373 int i, j, ntotal = 0;
374 shortedge *nearlist = (shortedge *) NULL;
375 double scale;
376 int goal = (xn->dat.z == (double *) NULL ? 4 * nearnum : 8 * nearnum);
377 int n = xn->invnames[ni];
378 int norm;
379
380 nearlist = CC_SAFE_MALLOC (nearnum + 1, shortedge);
381 if (!nearlist)
382 return 1;
383 for (i = 0; i < nearnum; i++)
384 nearlist[i].length = BIGDOUBLE;
385 nearlist[nearnum].length = -BIGDOUBLE;
386
387 CCutil_dat_getnorm (&xn->dat, &norm);
388 if (norm == CC_GEOGRAPHIC) scale = CC_GEOGRAPHIC_SCALE;
389 else if (norm == CC_GEOM) scale = CC_GEOM_SCALE;
390 else if (norm == CC_ATT) scale = CC_ATT_SCALE;
391 else scale = 1.0;
392 if ((norm & CC_NORM_SIZE_BITS) == CC_D3_NORM_SIZE) {
393 double ny = xn->dat.y[n];
394 double nz = xn->dat.z[n];
395 for (j = n - 1; j >= 0 && dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
396 < nearlist[0].length; --j) {
397 if (xn->dat.y[j] <= ny && xn->dat.z[j] <= nz)
398 insert (n, j, nearlist, &(xn->dat), xn->w);
399 }
400 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
401 xn->nodenames);
402 for (j = n - 1; j >= 0 && dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
403 < nearlist[0].length; --j) {
404 if (xn->dat.y[j] <= ny && xn->dat.z[j] >= nz)
405 insert (n, j, nearlist, &(xn->dat), xn->w);
406 }
407 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
408 xn->nodenames);
409 for (j = n - 1; j >= 0 && dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
410 < nearlist[0].length; --j) {
411 if (xn->dat.y[j] >= ny && xn->dat.z[j] <= nz)
412 insert (n, j, nearlist, &(xn->dat), xn->w);
413 }
414 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
415 xn->nodenames);
416 for (j = n - 1; j >= 0 && dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
417 < nearlist[0].length; --j) {
418 if (xn->dat.y[j] >= ny && xn->dat.z[j] >= nz)
419 insert (n, j, nearlist, &(xn->dat), xn->w);
420 }
421 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
422 xn->nodenames);
423 for (j = n + 1; j < ncount &&
424 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
425 < nearlist[0].length; j++) {
426 if (xn->dat.y[j] <= ny && xn->dat.z[j] <= nz)
427 insert (n, j, nearlist, &(xn->dat), xn->w);
428 }
429 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
430 xn->nodenames);
431 for (j = n + 1; j < ncount &&
432 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
433 < nearlist[0].length; j++) {
434 if (xn->dat.y[j] <= ny && xn->dat.z[j] >= nz)
435 insert (n, j, nearlist, &(xn->dat), xn->w);
436 }
437 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
438 xn->nodenames);
439 for (j = n + 1; j < ncount &&
440 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
441 < nearlist[0].length; j++) {
442 if (xn->dat.y[j] >= ny && xn->dat.z[j] <= nz)
443 insert (n, j, nearlist, &(xn->dat), xn->w);
444 }
445 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
446 xn->nodenames);
447 for (j = n + 1; j < ncount &&
448 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
449 < nearlist[0].length; j++) {
450 if (xn->dat.y[j] >= ny && xn->dat.z[j] >= nz)
451 insert (n, j, nearlist, &(xn->dat), xn->w);
452 }
453 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
454 xn->nodenames);
455 } else {
456 double ny = xn->dat.y[n];
457 for (j = n - 1; j >= 0 &&
458 dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
459 < nearlist[0].length; --j) {
460 if (xn->dat.y[j] <= ny)
461 insert (n, j, nearlist, &(xn->dat), xn->w);
462 }
463 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
464 xn->nodenames);
465 for (j = n - 1; j >= 0 &&
466 dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
467 < nearlist[0].length; --j) {
468 if (xn->dat.y[j] >= ny)
469 insert (n, j, nearlist, &(xn->dat), xn->w);
470 }
471 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
472 xn->nodenames);
473 for (j = n + 1; j < ncount &&
474 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
475 < nearlist[0].length; j++) {
476 if (xn->dat.y[j] <= ny)
477 insert (n, j, nearlist, &(xn->dat), xn->w);
478 }
479 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
480 xn->nodenames);
481 for (j = n + 1; j < ncount &&
482 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
483 < nearlist[0].length; j++) {
484 if (xn->dat.y[j] >= ny)
485 insert (n, j, nearlist, &(xn->dat), xn->w);
486 }
487 add_to_list_and_reset (list, &ntotal, nearlist, nearnum,
488 xn->nodenames);
489 }
490
491 if (ntotal < goal) {
492 for (i = ntotal; i < goal; i++)
493 list[i] = -1;
494 }
495 CC_IFFREE (nearlist, shortedge);
496 return 0;
497 }
498
add_to_list_and_reset(int * list,int * lcount,shortedge * nearlist,int nearnum,int * nodenames)499 static void add_to_list_and_reset (int *list, int *lcount, shortedge *nearlist,
500 int nearnum, int *nodenames)
501 {
502 int i;
503
504 for (i = 0; i < nearnum; i++) {
505 if (nearlist[i].length < BIGDOUBLE) {
506 list[(*lcount)++] = nodenames[nearlist[i].end];
507 nearlist[i].length = BIGDOUBLE;
508 }
509 }
510 }
511
CCedgegen_x_node_k_nearest(CCxnear * xn,int ni,int nearnum,int ncount,int * list)512 int CCedgegen_x_node_k_nearest (CCxnear *xn, int ni, int nearnum, int ncount,
513 int *list)
514 {
515 int i, j, ntotal;
516 shortedge *nearlist = (shortedge *) NULL;
517 double scale;
518 int n = xn->invnames[ni];
519 int norm;
520
521 nearlist = CC_SAFE_MALLOC (nearnum + 1, shortedge);
522 if (!nearlist)
523 return 1;
524 for (i = 0; i < nearnum; i++)
525 nearlist[i].length = BIGDOUBLE;
526 nearlist[nearnum].length = -BIGDOUBLE;
527
528 CCutil_dat_getnorm (&xn->dat, &norm);
529 if (norm == CC_GEOGRAPHIC) scale = CC_GEOGRAPHIC_SCALE;
530 else if (norm == CC_GEOM) scale = CC_GEOM;
531 else if (norm == CC_ATT) scale = CC_ATT_SCALE;
532 else scale = 1.0;
533 for (j = n - 1; j >= 0 &&
534 dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
535 < nearlist[0].length; --j)
536 insert (n, j, nearlist, &(xn->dat), xn->w);
537 for (j = n + 1; j < ncount &&
538 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
539 < nearlist[0].length; j++)
540 insert (n, j, nearlist, &(xn->dat), xn->w);
541
542 ntotal = 0;
543 for (i = 0; i < nearnum; i++) {
544 if (nearlist[i].length < BIGDOUBLE)
545 list[ntotal++] = xn->nodenames[nearlist[i].end];
546 }
547 if (ntotal < nearnum) {
548 fprintf (stderr, "WARNING: There do not exist %d neighbors\n",
549 nearnum);
550 for (i = ntotal; i < nearnum; i++)
551 list[i] = -1;
552 return 1;
553 }
554
555 CC_IFFREE (nearlist, shortedge);
556 return 0;
557 }
558
CCedgegen_x_node_nearest(CCxnear * xn,int ncount,int ni,char * marks)559 int CCedgegen_x_node_nearest (CCxnear *xn, int ncount, int ni, char *marks)
560 {
561 int n = xn->invnames[ni];
562 int j, bestnode = 0;
563 double scale, thisdist, bestdist = BIGDOUBLE;
564 int norm;
565
566 CCutil_dat_getnorm (&xn->dat, &norm);
567 if (norm == CC_GEOGRAPHIC) scale = CC_GEOGRAPHIC_SCALE;
568 else if (norm == CC_GEOM) scale = CC_GEOM;
569 else if (norm == CC_ATT) scale = CC_ATT_SCALE;
570 else scale = 1.0;
571 for (j = n - 1; j >= 0 &&
572 dtrunc((xn->dat.x[n] - xn->dat.x[j]) * scale)
573 < bestdist; --j) {
574 if (!marks[xn->nodenames[j]]) {
575 thisdist = CCutil_dat_edgelen (n, j, &(xn->dat));
576 if (xn->w)
577 thisdist += (xn->w[n] + xn->w[j]);
578 if (thisdist < bestdist) {
579 bestdist = thisdist;
580 bestnode = j;
581 }
582 }
583 }
584 for (j = n + 1; j < ncount &&
585 dtrunc((xn->dat.x[j] - xn->dat.x[n]) * scale)
586 < bestdist; j++) {
587 if (!marks[xn->nodenames[j]]) {
588 thisdist = CCutil_dat_edgelen (n, j, &(xn->dat));
589 if (xn->w)
590 thisdist += (xn->w[n] + xn->w[j]);
591 if (thisdist < bestdist) {
592 bestdist = thisdist;
593 bestnode = j;
594 }
595 }
596 }
597 return xn->nodenames[bestnode];
598 }
599
insert(int n,int m,shortedge * nearlist,CCdatagroup * dat,double * wcoord)600 static void insert (int n, int m, shortedge *nearlist, CCdatagroup *dat,
601 double *wcoord)
602 {
603 int i;
604 int thisdist;
605
606 thisdist = CCutil_dat_edgelen (n, m, dat);
607
608 if (wcoord != (double *) NULL)
609 thisdist += (wcoord[n] + wcoord[m]);
610
611 if (thisdist < nearlist[0].length) {
612 for (i = 0; nearlist[i+1].length > thisdist; i++) {
613 nearlist[i].end = nearlist[i + 1].end;
614 nearlist[i].length = nearlist[i + 1].length;
615 }
616 nearlist[i].length = thisdist;
617 nearlist[i].end = m;
618 }
619 }
620
CCedgegen_junk_node_k_nearest(CCdatagroup * dat,double * wcoord,int n,int nearnum,int ncount,int * list)621 int CCedgegen_junk_node_k_nearest (CCdatagroup *dat, double *wcoord, int n,
622 int nearnum, int ncount, int *list)
623 {
624 int i, j, ntotal;
625 shortedge *nearlist = (shortedge *) NULL;
626
627 nearlist = CC_SAFE_MALLOC (nearnum + 1, shortedge);
628 if (!nearlist)
629 return 1;
630 for (i = 0; i < nearnum; i++)
631 nearlist[i].length = BIGDOUBLE;
632 nearlist[nearnum].length = -BIGDOUBLE;
633
634 for (j = n - 1; j >= 0; j--) {
635 insert (n, j, nearlist, dat, wcoord);
636 }
637 for (j = n + 1; j < ncount; j++) {
638 insert (n, j, nearlist, dat, wcoord);
639 }
640
641 ntotal = 0;
642 for (i = 0; i < nearnum; i++) {
643 if (nearlist[i].length < BIGDOUBLE)
644 list[ntotal++] = nearlist[i].end;
645 }
646 if (ntotal < nearnum) {
647 fprintf (stderr, "WARNING: There do not exist %d neighbors\n",
648 nearnum);
649 for (i = ntotal; i < nearnum; i++)
650 list[i] = -1;
651 return 1;
652 }
653
654 CC_IFFREE (nearlist, shortedge);
655 return 0;
656 }
657
CCedgegen_junk_node_nearest(CCdatagroup * dat,double * wcoord,int ncount,int n,char * marks)658 int CCedgegen_junk_node_nearest (CCdatagroup *dat, double *wcoord, int ncount,
659 int n, char *marks)
660 {
661 int j, bestnode = 0;
662 double thisdist, bestdist = BIGDOUBLE;
663
664
665 if (wcoord) {
666 for (j = n - 1; j >= 0; j--) {
667 if (!marks[j]) {
668 thisdist = CCutil_dat_edgelen (n, j, dat) +
669 (wcoord[n] + wcoord[j]);
670 if (thisdist < bestdist) {
671 bestdist = thisdist;
672 bestnode = j;
673 }
674 }
675 }
676 for (j = n + 1; j < ncount; j++) {
677 if (!marks[j]) {
678 thisdist = CCutil_dat_edgelen (n, j, dat) +
679 (wcoord[n] + wcoord[j]);
680 if (thisdist < bestdist) {
681 bestdist = thisdist;
682 bestnode = j;
683 }
684 }
685 }
686 } else {
687 for (j = n - 1; j >= 0; j--) {
688 if (!marks[j]) {
689 thisdist = CCutil_dat_edgelen (n, j, dat);
690 if (thisdist < bestdist) {
691 bestdist = thisdist;
692 bestnode = j;
693 }
694 }
695 }
696 for (j = n + 1; j < ncount; j++) {
697 if (!marks[j]) {
698 thisdist = CCutil_dat_edgelen (n, j, dat);
699 if (thisdist < bestdist) {
700 bestdist = thisdist;
701 bestnode = j;
702 }
703 }
704 }
705 }
706
707 return bestnode;
708 }
709
CCedgegen_x_nearest_neighbor_tour(int ncount,int start,CCdatagroup * dat,int * outcycle,double * val)710 int CCedgegen_x_nearest_neighbor_tour (int ncount, int start, CCdatagroup *dat,
711 int *outcycle, double *val)
712 {
713 double len;
714 int i, current, next;
715 CCxnear xn;
716 char *marks;
717 int norm;
718
719 /*
720 printf ("Grow nearest neighbor tour from node %d\n", start);
721 fflush (stdout);
722 */
723
724 if (ncount < 3) {
725 fprintf (stderr, "Cannot find tour in an %d node graph\n", ncount);
726 return 1;
727 }
728 CCutil_dat_getnorm (dat, &norm);
729 if ((norm & CC_NORM_BITS) != CC_X_NORM_TYPE &&
730 (norm & CC_NORM_BITS) != CC_KD_NORM_TYPE) {
731 fprintf (stderr, "Cannot run x_nearest with norm %d\n", norm);
732 return 1;
733 }
734
735 if (CCedgegen_xnear_build (ncount, dat, (double *) NULL, &xn)) {
736 fprintf (stderr, "Unable to build CCxnear\n");
737 return 1;
738 }
739
740 marks = CC_SAFE_MALLOC (ncount, char );
741 if (!marks) {
742 CCedgegen_xnear_free (&xn);
743 return 1;
744 }
745
746 for (i = 0; i < ncount; i++)
747 marks[i] = 0;
748
749
750 len = 0.0;
751 current = start;
752 if (outcycle != (int *) NULL)
753 outcycle[0] = start;
754
755 for (i = 1; i < ncount; i++) {
756 marks[current] = 1;
757 next = CCedgegen_x_node_nearest (&xn, ncount, current, marks);
758 if (outcycle != (int *) NULL)
759 outcycle [i] = next;
760 len += (double) CCutil_dat_edgelen (current, next, dat);
761 current = next;
762 }
763 len += (double) CCutil_dat_edgelen (current, start, dat);
764 *val = len;
765 CCedgegen_xnear_free (&xn);
766 CC_IFFREE (marks, char);
767 return 0;
768 }
769
CCedgegen_x_greedy_tour(int ncount,CCdatagroup * dat,int * outcycle,double * val,int ecount,int * elist,int silent)770 int CCedgegen_x_greedy_tour (int ncount, CCdatagroup *dat, int *outcycle,
771 double *val, int ecount, int *elist, int silent)
772 {
773 int rval;
774 int *perm = (int *) NULL;
775 int *elen = (int *) NULL;
776 int *tcyc = (int *) NULL;
777 int *degree = (int *) NULL;
778 int *tail = (int *) NULL;
779 char *marks = (char *) NULL;
780 int norm;
781 CCxnear xn;
782 int tcount = 0;
783 int x;
784 int y;
785 int z;
786 int i;
787 double len;
788
789 if (!silent) {
790 printf ("Grow a greedy tour \n");
791 fflush (stdout);
792 }
793
794 CCutil_dat_getnorm (dat, &norm);
795 if ((norm & CC_NORM_BITS) != CC_X_NORM_TYPE &&
796 (norm & CC_NORM_BITS) != CC_KD_NORM_TYPE) {
797 fprintf (stderr, "Cannot run x_nearest with norm %d\n", norm);
798 return 1;
799 }
800
801 if (CCedgegen_xnear_build (ncount, dat, (double *) NULL, &xn)) {
802 fprintf (stderr, "Unable to build CCxnear\n");
803 return 1;
804 }
805
806 perm = CC_SAFE_MALLOC (ecount, int);
807 elen = CC_SAFE_MALLOC (ecount, int);
808 if (perm == (int *) NULL ||
809 elen == (int *) NULL) {
810 fprintf (stderr, "Out of memory in CCedgegen_x_greedy_tour\n");
811 rval = 1; goto CLEANUP;
812 }
813 for (i=0; i<ecount; i++) {
814 perm[i] = i;
815 elen[i] = CCutil_dat_edgelen (elist[2*i], elist[2*i+1], dat);
816 }
817 CCutil_int_perm_quicksort (perm, elen, ecount);
818
819 if (outcycle) {
820 tcyc = CC_SAFE_MALLOC (2 * ncount, int);
821 if (!tcyc) {
822 rval = 1;
823 goto CLEANUP;
824 }
825 }
826
827 degree = CC_SAFE_MALLOC (ncount, int);
828 if (!degree) {
829 rval = 1;
830 goto CLEANUP;
831 }
832 tail = CC_SAFE_MALLOC (ncount, int);
833 if (!tail) {
834 rval = 1;
835 goto CLEANUP;
836 }
837
838 for (i = 0; i < ncount; i++) {
839 degree[i] = 0;
840 tail[i] = -1;
841 }
842
843 len = 0.0;
844 for (i = 0; i < ecount; i++) {
845 x = elist[2*perm[i]];
846 y = elist[2*perm[i]+1];
847 if (degree[x] < 2 && degree[y] < 2 && y != tail[x]) {
848 if (tcyc) {
849 tcyc[2*tcount] = x;
850 tcyc[2*tcount+1] = y;
851 }
852 tcount++;
853 len += elen[perm[i]];
854 (degree[x])++;
855 (degree[y])++;
856 if (tail[x] == -1) {
857 if (tail[y] == -1) {
858 tail[x] = y;
859 tail[y] = x;
860 } else {
861 tail[x] = tail[y];
862 tail[tail[y]] = x;
863 }
864 } else if (tail[y] == -1) {
865 tail[tail[x]] = y;
866 tail[y] = tail[x];
867 } else {
868 tail[tail[x]] = tail[y];
869 tail[tail[y]] = tail[x];
870 }
871 }
872 }
873 CC_FREE (perm, int);
874 CC_FREE (elen, int);
875 marks = CC_SAFE_MALLOC (ncount, char);
876 if (!marks) {
877 rval = 1; goto CLEANUP;
878 }
879 for (i=0; i<ncount; i++) {
880 if (degree[i] == 2) marks[i] = 1;
881 else marks[i] = 0;
882 }
883 for (x = 0; marks[x] == 1; x++);
884 if (tail[x] != -1) marks[tail[x]] = 1;
885 while (tcount < ncount-1) {
886 marks[x] = 1;
887 y = CCedgegen_x_node_nearest (&xn, ncount, x, marks);
888 if (tcyc) {
889 tcyc[2*tcount] = x;
890 tcyc[2*tcount+1] = y;
891 }
892 tcount++;
893 len += (double) CCutil_dat_edgelen (x, y, dat);
894 if (tail[y] == -1) z = y;
895 else z = tail[y];
896 if (degree[x]) marks[x] = 1;
897 (degree[x])++;
898 if (degree[y]) marks[y] = 1;
899 (degree[y])++;
900 if (tail[x] == -1) {
901 if (tail[y] == -1) {
902 tail[x] = y;
903 tail[y] = x;
904 } else {
905 tail[x] = tail[y];
906 tail[tail[y]] = x;
907 }
908 } else if (tail[y] == -1) {
909 tail[tail[x]] = y;
910 tail[y] = tail[x];
911 } else {
912 tail[tail[x]] = tail[y];
913 tail[tail[y]] = tail[x];
914 }
915 x = z;
916 }
917 for (x = 0; degree[x] != 1; x++);
918 for (y = x + 1; degree[y] != 1; y++);
919 if (tcyc) {
920 tcyc[2*tcount] = x;
921 tcyc[2*tcount+1] = y;
922 }
923 tcount++;
924 len += (double) CCutil_dat_edgelen (x, y, dat);
925 *val = len;
926 if (!silent) {
927 printf ("Length of Greedy Tour: %.2f\n", len);
928 fflush (stdout);
929 }
930
931 if (tcyc) {
932 int istour;
933 rval = CCutil_edge_to_cycle (ncount, tcyc, &istour, outcycle);
934 if (rval) {
935 fprintf (stderr, "CCutil_edge_to_cycle failed\n"); goto CLEANUP;
936 }
937 if (istour == 0) {
938 fprintf (stderr, "ERROR: greedy tour is not a tour\n");
939 rval = 1; goto CLEANUP;
940 }
941 }
942 rval = 0;
943
944 CLEANUP:
945 CCedgegen_xnear_free (&xn);
946 CC_IFFREE (tcyc, int);
947 CC_IFFREE (degree, int);
948 CC_IFFREE (tail, int);
949 CC_IFFREE (marks, char);
950 CC_IFFREE (perm, int);
951 CC_IFFREE (elen, int);
952 return rval;
953 }
954
CCedgegen_x_qboruvka_tour(int ncount,CCdatagroup * dat,int * outcycle,double * val,int ecount,int * elist,int silent)955 int CCedgegen_x_qboruvka_tour (int ncount, CCdatagroup *dat, int *outcycle,
956 double *val, int ecount, int *elist, int silent)
957 {
958 int rval;
959 int *perm = (int *) NULL;
960 int *elen = (int *) NULL;
961 int *deg = (int *) NULL;
962 int **adj = (int **) NULL;
963 int *adjspace = (int *) NULL;
964 int *tcyc = (int *) NULL;
965 int *degree = (int *) NULL;
966 int *tail = (int *) NULL;
967 char *marks = (char *) NULL;
968 int tcount = 0;
969 double len;
970 int norm;
971 CCxnear xn;
972 int x;
973 int y = 0;
974 int i;
975 int j;
976 int k;
977 int count;
978 int found;
979
980 if (!silent) {
981 printf ("Grow a Quick-Boruvka tour \n");
982 fflush (stdout);
983 }
984
985 CCutil_dat_getnorm (dat, &norm);
986 if ((norm & CC_NORM_BITS) != CC_X_NORM_TYPE &&
987 (norm & CC_NORM_BITS) != CC_KD_NORM_TYPE) {
988 fprintf (stderr, "Cannot run x_nearest with norm %d\n", norm);
989 return 1;
990 }
991
992 if (CCedgegen_xnear_build (ncount, dat, (double *) NULL, &xn)) {
993 fprintf (stderr, "Unable to build CCxnear\n");
994 return 1;
995 }
996
997 perm = CC_SAFE_MALLOC (ecount, int);
998 elen = CC_SAFE_MALLOC (ecount, int);
999 if (perm == (int *) NULL ||
1000 elen == (int *) NULL) {
1001 fprintf (stderr, "Out of memory in CCedgegen_x_greedy_tour\n");
1002 rval = 1; goto CLEANUP;
1003 }
1004 for (i=0; i<ecount; i++) {
1005 perm[i] = i;
1006 elen[i] = CCutil_dat_edgelen (elist[2*i], elist[2*i+1], dat);
1007 }
1008 CCutil_int_perm_quicksort (perm, elen, ecount);
1009
1010 CC_FREE (elen, int);
1011
1012 deg = CC_SAFE_MALLOC (ncount, int);
1013 adj = CC_SAFE_MALLOC (ncount, int *);
1014 adjspace = CC_SAFE_MALLOC (ecount*2, int);
1015 if (deg == (int *) NULL ||
1016 adj == (int **) NULL ||
1017 adjspace == (int *) NULL) {
1018 fprintf (stderr, "Out of memory in CCedgegen_x_greedy_tour\n");
1019 rval = 1; goto CLEANUP;
1020 }
1021 for (i=0; i<ncount; i++) {
1022 deg[i] = 0;
1023 }
1024
1025 for (i=0; i<ecount; i++) {
1026 deg[elist[2*i]]++;
1027 deg[elist[2*i+1]]++;
1028 }
1029 j = 0;
1030 for (i=0; i<ncount; i++) {
1031 adj[i] = &(adjspace[j]);
1032 j += deg[i];
1033 deg[i] = 0;
1034 }
1035
1036 for (i=0; i<ecount; i++) {
1037 j = elist[2*perm[i]];
1038 k = elist[2*perm[i]+1];
1039 adj[j][deg[j]++] = k;
1040 adj[k][deg[k]++] = j;
1041 }
1042 CC_FREE (perm, int);
1043
1044 if (outcycle) {
1045 tcyc = CC_SAFE_MALLOC (2 * ncount, int);
1046 if (!tcyc) {
1047 rval = 1;
1048 goto CLEANUP;
1049 }
1050 }
1051
1052 degree = CC_SAFE_MALLOC (ncount, int);
1053 if (!degree) {
1054 rval = 1;
1055 goto CLEANUP;
1056 }
1057 tail = CC_SAFE_MALLOC (ncount, int);
1058 if (!tail) {
1059 rval = 1;
1060 goto CLEANUP;
1061 }
1062 marks = CC_SAFE_MALLOC (ncount, char);
1063 if (!marks) {
1064 rval = 1;
1065 goto CLEANUP;
1066 }
1067
1068 for (i = 0; i < ncount; i++) {
1069 degree[i] = 0;
1070 tail[i] = -1;
1071 marks[i] = 0;
1072 }
1073
1074 len = 0.0;
1075 count = 1;
1076 while (count < ncount) {
1077 for (x = 0; x < ncount && count < ncount; x++) {
1078 if (degree[x] != 2) {
1079 found = 0;
1080 for (j=0; j<deg[x] && !found; j++) {
1081 y = adj[x][j];
1082 if (degree[y] != 2 && y != tail[x]) {
1083 found = 1;
1084 }
1085 }
1086 if (!found) {
1087 if (tail[x] != -1) {
1088 marks[tail[x]] = 1;
1089 y = CCedgegen_x_node_nearest (&xn, ncount, x, marks);
1090 marks[tail[x]] = 0;
1091 } else {
1092 y = CCedgegen_x_node_nearest (&xn, ncount, x, marks);
1093 }
1094 }
1095
1096 /* add (x, y) to the tour */
1097 if (degree[x] != 0)
1098 marks[x] = 1;
1099 if (degree[y] != 0)
1100 marks[y] = 1;
1101 len += (double) CCutil_dat_edgelen (x, y, dat);
1102 degree[x]++;
1103 degree[y]++;
1104 if (tcyc) {
1105 tcyc[tcount++] = x;
1106 tcyc[tcount++] = y;
1107 }
1108 if (tail[x] == -1) {
1109 if (tail[y] == -1) {
1110 tail[x] = y;
1111 tail[y] = x;
1112 } else {
1113 tail[x] = tail[y];
1114 tail[tail[y]] = x;
1115 }
1116 } else if (tail[y] == -1) {
1117 tail[tail[x]] = y;
1118 tail[y] = tail[x];
1119 } else {
1120 tail[tail[x]] = tail[y];
1121 tail[tail[y]] = tail[x];
1122 }
1123 if (count % 10000 == 9999) {
1124 printf (".");
1125 fflush (stdout);
1126 }
1127 count++;
1128 }
1129 }
1130 }
1131 for (x = 0; degree[x] != 1; x++);
1132 for (y = x + 1; degree[y] != 1; y++);
1133 if (tcyc) {
1134 tcyc[tcount++] = x;
1135 tcyc[tcount++] = y;
1136 }
1137 len += (double) CCutil_dat_edgelen (x, y, dat);
1138 *val = len;
1139 if (ncount >= 10000)
1140 printf ("\n");
1141 printf ("Length of Quick-Boruvka Tour: %.2f\n", len);
1142
1143 if (tcyc) {
1144 int istour;
1145 rval = CCutil_edge_to_cycle (ncount, tcyc, &istour, outcycle);
1146 if (rval) {
1147 fprintf (stderr, "CCutil_edge_to_cycle failed\n"); goto CLEANUP;
1148 }
1149 if (istour == 0) {
1150 fprintf (stderr, "ERROR: greedy tour is not a tour\n");
1151 rval = 1; goto CLEANUP;
1152 }
1153 }
1154 rval = 0;
1155
1156 CLEANUP:
1157 CCedgegen_xnear_free (&xn);
1158 CC_IFFREE (tcyc, int);
1159 CC_IFFREE (degree, int);
1160 CC_IFFREE (tail, int);
1161 CC_IFFREE (marks, char);
1162 CC_IFFREE (deg, int);
1163 CC_IFFREE (adj, int *);
1164 CC_IFFREE (adjspace, int);
1165 CC_IFFREE (perm, int);
1166 CC_IFFREE (elen, int);
1167 return rval;
1168 }
1169
CCedgegen_junk_nearest_neighbor_tour(int ncount,int start,CCdatagroup * dat,int * outcycle,double * val,int silent)1170 int CCedgegen_junk_nearest_neighbor_tour (int ncount, int start,
1171 CCdatagroup *dat, int *outcycle, double *val, int silent)
1172 {
1173 double len;
1174 int i, current, next;
1175 char *marks;
1176
1177 /*
1178 printf ("Grow nearest neighbor tour from node %d\n", start);
1179 */
1180 if (!silent) {
1181 printf ("This is a JUNK norm, so expect a quadratic running time\n");
1182 fflush (stdout);
1183 }
1184
1185 if (ncount < 3) {
1186 fprintf (stderr, "Cannot find tour in an %d node graph\n", ncount);
1187 return 1;
1188 }
1189
1190 marks = CC_SAFE_MALLOC (ncount, char );
1191 if (!marks) {
1192 return 1;
1193 }
1194
1195 for (i = 0; i < ncount; i++)
1196 marks[i] = 0;
1197
1198 len = 0.0;
1199 current = start;
1200 if (outcycle != (int *) NULL)
1201 outcycle[0] = start;
1202
1203 for (i = 1; i < ncount; i++) {
1204 marks[current] = 1;
1205 next = CCedgegen_junk_node_nearest (dat, (double *) NULL, ncount,
1206 current, marks);
1207 if (outcycle != (int *) NULL)
1208 outcycle [i] = next;
1209 len += (double) CCutil_dat_edgelen (current, next, dat);
1210 current = next;
1211 }
1212 len += (double) CCutil_dat_edgelen (current, start, dat);
1213 *val = len;
1214 CC_IFFREE (marks, char);
1215 return 0;
1216 }
1217
CCedgegen_junk_greedy_tour(int ncount,CCdatagroup * dat,int * outcycle,double * val,int ecount,int * elist,int silent)1218 int CCedgegen_junk_greedy_tour (int ncount, CCdatagroup *dat, int *outcycle,
1219 double *val, int ecount, int *elist, int silent)
1220 {
1221 int rval;
1222 int *perm = (int *) NULL;
1223 int *elen = (int *) NULL;
1224 int *tcyc = (int *) NULL;
1225 int *degree = (int *) NULL;
1226 int *tail = (int *) NULL;
1227 char *marks = (char *) NULL;
1228 int tcount = 0;
1229 int x;
1230 int y;
1231 int z;
1232 int i;
1233 double len;
1234
1235 if (!silent) {
1236 printf ("Grow a greedy tour \n");
1237 fflush (stdout);
1238 }
1239
1240 perm = CC_SAFE_MALLOC (ecount, int);
1241 elen = CC_SAFE_MALLOC (ecount, int);
1242 if (perm == (int *) NULL ||
1243 elen == (int *) NULL) {
1244 fprintf (stderr, "Out of memory in CCedgegen_junk_greedy_tour\n");
1245 rval = 1; goto CLEANUP;
1246 }
1247 for (i=0; i<ecount; i++) {
1248 perm[i] = i;
1249 elen[i] = CCutil_dat_edgelen (elist[2*i], elist[2*i+1], dat);
1250 }
1251 CCutil_int_perm_quicksort (perm, elen, ecount);
1252
1253 if (outcycle) {
1254 tcyc = CC_SAFE_MALLOC (2 * ncount, int);
1255 if (!tcyc) {
1256 rval = 1;
1257 goto CLEANUP;
1258 }
1259 }
1260
1261 degree = CC_SAFE_MALLOC (ncount, int);
1262 if (!degree) {
1263 rval = 1;
1264 goto CLEANUP;
1265 }
1266 tail = CC_SAFE_MALLOC (ncount, int);
1267 if (!tail) {
1268 rval = 1;
1269 goto CLEANUP;
1270 }
1271
1272 for (i = 0; i < ncount; i++) {
1273 degree[i] = 0;
1274 tail[i] = -1;
1275 }
1276
1277 len = 0.0;
1278 for (i = 0; i < ecount; i++) {
1279 x = elist[2*perm[i]];
1280 y = elist[2*perm[i]+1];
1281 if (degree[x] < 2 && degree[y] < 2 && y != tail[x]) {
1282 if (tcyc) {
1283 tcyc[2*tcount] = x;
1284 tcyc[2*tcount+1] = y;
1285 }
1286 tcount++;
1287 len += elen[perm[i]];
1288 (degree[x])++;
1289 (degree[y])++;
1290 if (tail[x] == -1) {
1291 if (tail[y] == -1) {
1292 tail[x] = y;
1293 tail[y] = x;
1294 } else {
1295 tail[x] = tail[y];
1296 tail[tail[y]] = x;
1297 }
1298 } else if (tail[y] == -1) {
1299 tail[tail[x]] = y;
1300 tail[y] = tail[x];
1301 } else {
1302 tail[tail[x]] = tail[y];
1303 tail[tail[y]] = tail[x];
1304 }
1305 }
1306 }
1307 CC_FREE (perm, int);
1308 CC_FREE (elen, int);
1309 marks = CC_SAFE_MALLOC (ncount, char);
1310 if (!marks) {
1311 rval = 1; goto CLEANUP;
1312 }
1313 for (i=0; i<ncount; i++) {
1314 if (degree[i] == 2) marks[i] = 1;
1315 else marks[i] = 0;
1316 }
1317 for (x = 0; marks[x] == 1; x++);
1318 if (tail[x] != -1) marks[tail[x]] = 1;
1319 while (tcount < ncount-1) {
1320 marks[x] = 1;
1321 y = CCedgegen_junk_node_nearest (dat, (double *) NULL, ncount, x,
1322 marks);
1323 if (tcyc) {
1324 tcyc[2*tcount] = x;
1325 tcyc[2*tcount+1] = y;
1326 }
1327 tcount++;
1328 len += (double) CCutil_dat_edgelen (x, y, dat);
1329 if (tail[y] == -1) z = y;
1330 else z = tail[y];
1331 if (degree[x]) marks[x] = 1;
1332 (degree[x])++;
1333 if (degree[y]) marks[y] = 1;
1334 (degree[y])++;
1335 if (tail[x] == -1) {
1336 if (tail[y] == -1) {
1337 tail[x] = y;
1338 tail[y] = x;
1339 } else {
1340 tail[x] = tail[y];
1341 tail[tail[y]] = x;
1342 }
1343 } else if (tail[y] == -1) {
1344 tail[tail[x]] = y;
1345 tail[y] = tail[x];
1346 } else {
1347 tail[tail[x]] = tail[y];
1348 tail[tail[y]] = tail[x];
1349 }
1350 x = z;
1351 }
1352 for (x = 0; degree[x] != 1; x++);
1353 for (y = x + 1; degree[y] != 1; y++);
1354 if (tcyc) {
1355 tcyc[2*tcount] = x;
1356 tcyc[2*tcount+1] = y;
1357 }
1358 tcount++;
1359 len += (double) CCutil_dat_edgelen (x, y, dat);
1360 *val = len;
1361 if (!silent) {
1362 printf ("Length of Greedy Tour: %.2f\n", len);
1363 fflush (stdout);
1364 }
1365
1366 if (tcyc) {
1367 int istour;
1368 rval = CCutil_edge_to_cycle (ncount, tcyc, &istour, outcycle);
1369 if (rval) {
1370 fprintf (stderr, "CCutil_edge_to_cycle failed\n"); goto CLEANUP;
1371 }
1372 if (istour == 0) {
1373 fprintf (stderr, "ERROR: greedy tour is not a tour\n");
1374 rval = 1; goto CLEANUP;
1375 }
1376 }
1377 rval = 0;
1378
1379 CLEANUP:
1380 CC_IFFREE (tcyc, int);
1381 CC_IFFREE (degree, int);
1382 CC_IFFREE (tail, int);
1383 CC_IFFREE (marks, char);
1384 CC_IFFREE (perm, int);
1385 CC_IFFREE (elen, int);
1386 return rval;
1387 }
1388
CCedgegen_junk_qboruvka_tour(int ncount,CCdatagroup * dat,int * outcycle,double * val,int ecount,int * elist,int silent)1389 int CCedgegen_junk_qboruvka_tour (int ncount, CCdatagroup *dat, int *outcycle,
1390 double *val, int ecount, int *elist, int silent)
1391 {
1392 int rval;
1393 int *perm = (int *) NULL;
1394 int *elen = (int *) NULL;
1395 int *deg = (int *) NULL;
1396 int **adj = (int **) NULL;
1397 int *adjspace = (int *) NULL;
1398 int *tcyc = (int *) NULL;
1399 int *degree = (int *) NULL;
1400 int *tail = (int *) NULL;
1401 char *marks = (char *) NULL;
1402 int tcount = 0;
1403 double len;
1404 int x;
1405 int y = 0;
1406 int i;
1407 int j;
1408 int k;
1409 int count;
1410 int found;
1411
1412 if (!silent) {
1413 printf ("Grow a Quick-Boruvka tour \n");
1414 fflush (stdout);
1415 }
1416
1417 perm = CC_SAFE_MALLOC (ecount, int);
1418 elen = CC_SAFE_MALLOC (ecount, int);
1419 if (perm == (int *) NULL ||
1420 elen == (int *) NULL) {
1421 fprintf (stderr, "Out of memory in CCedgegen_junk_greedy_tour\n");
1422 rval = 1; goto CLEANUP;
1423 }
1424 for (i=0; i<ecount; i++) {
1425 perm[i] = i;
1426 elen[i] = CCutil_dat_edgelen (elist[2*i], elist[2*i+1], dat);
1427 }
1428 CCutil_int_perm_quicksort (perm, elen, ecount);
1429
1430 CC_FREE (elen, int);
1431
1432 deg = CC_SAFE_MALLOC (ncount, int);
1433 adj = CC_SAFE_MALLOC (ncount, int *);
1434 adjspace = CC_SAFE_MALLOC (ecount*2, int);
1435 if (deg == (int *) NULL ||
1436 adj == (int **) NULL ||
1437 adjspace == (int *) NULL) {
1438 fprintf (stderr, "Out of memory in CCedgegen_junk_greedy_tour\n");
1439 rval = 1; goto CLEANUP;
1440 }
1441 for (i=0; i<ncount; i++) {
1442 deg[i] = 0;
1443 }
1444
1445 for (i=0; i<ecount; i++) {
1446 deg[elist[2*i]]++;
1447 deg[elist[2*i+1]]++;
1448 }
1449 j = 0;
1450 for (i=0; i<ncount; i++) {
1451 adj[i] = &(adjspace[j]);
1452 j += deg[i];
1453 deg[i] = 0;
1454 }
1455
1456 for (i=0; i<ecount; i++) {
1457 j = elist[2*perm[i]];
1458 k = elist[2*perm[i]+1];
1459 adj[j][deg[j]++] = k;
1460 adj[k][deg[k]++] = j;
1461 }
1462 CC_FREE (perm, int);
1463
1464 if (outcycle) {
1465 tcyc = CC_SAFE_MALLOC (2 * ncount, int);
1466 if (!tcyc) {
1467 rval = 1;
1468 goto CLEANUP;
1469 }
1470 }
1471
1472 degree = CC_SAFE_MALLOC (ncount, int);
1473 if (!degree) {
1474 rval = 1;
1475 goto CLEANUP;
1476 }
1477 tail = CC_SAFE_MALLOC (ncount, int);
1478 if (!tail) {
1479 rval = 1;
1480 goto CLEANUP;
1481 }
1482 marks = CC_SAFE_MALLOC (ncount, char);
1483 if (!marks) {
1484 rval = 1;
1485 goto CLEANUP;
1486 }
1487
1488 for (i = 0; i < ncount; i++) {
1489 degree[i] = 0;
1490 tail[i] = -1;
1491 marks[i] = 0;
1492 }
1493
1494 len = 0.0;
1495 count = 1;
1496 while (count < ncount) {
1497 for (x = 0; x < ncount && count < ncount; x++) {
1498 if (degree[x] != 2) {
1499 found = 0;
1500 for (j=0; j<deg[x] && !found; j++) {
1501 y = adj[x][j];
1502 if (degree[y] != 2 && y != tail[x]) {
1503 found = 1;
1504 }
1505 }
1506 if (!found) {
1507 if (tail[x] != -1) {
1508 marks[tail[x]] = 1;
1509 y = CCedgegen_junk_node_nearest (dat, (double *) NULL,
1510 ncount, x, marks);
1511 marks[tail[x]] = 0;
1512 } else {
1513 y = CCedgegen_junk_node_nearest (dat, (double *) NULL,
1514 ncount, x, marks);
1515 }
1516 }
1517
1518 /* add (x, y) to the tour */
1519 if (degree[x] != 0)
1520 marks[x] = 1;
1521 if (degree[y] != 0)
1522 marks[y] = 1;
1523 len += (double) CCutil_dat_edgelen (x, y, dat);
1524 degree[x]++;
1525 degree[y]++;
1526 if (tcyc) {
1527 tcyc[tcount++] = x;
1528 tcyc[tcount++] = y;
1529 }
1530 if (tail[x] == -1) {
1531 if (tail[y] == -1) {
1532 tail[x] = y;
1533 tail[y] = x;
1534 } else {
1535 tail[x] = tail[y];
1536 tail[tail[y]] = x;
1537 }
1538 } else if (tail[y] == -1) {
1539 tail[tail[x]] = y;
1540 tail[y] = tail[x];
1541 } else {
1542 tail[tail[x]] = tail[y];
1543 tail[tail[y]] = tail[x];
1544 }
1545 if (count % 10000 == 9999) {
1546 printf (".");
1547 fflush (stdout);
1548 }
1549 count++;
1550 }
1551 }
1552 }
1553 for (x = 0; degree[x] != 1; x++);
1554 for (y = x + 1; degree[y] != 1; y++);
1555 if (tcyc) {
1556 tcyc[tcount++] = x;
1557 tcyc[tcount++] = y;
1558 }
1559 len += (double) CCutil_dat_edgelen (x, y, dat);
1560 *val = len;
1561 if (ncount >= 10000)
1562 printf ("\n");
1563 printf ("Length of Quick-Boruvka Tour: %.2f\n", len);
1564
1565 if (tcyc) {
1566 int istour;
1567 rval = CCutil_edge_to_cycle (ncount, tcyc, &istour, outcycle);
1568 if (rval) {
1569 fprintf (stderr, "CCutil_edge_to_cycle failed\n"); goto CLEANUP;
1570 }
1571 if (istour == 0) {
1572 fprintf (stderr, "ERROR: greedy tour is not a tour\n");
1573 rval = 1; goto CLEANUP;
1574 }
1575 }
1576 rval = 0;
1577
1578 CLEANUP:
1579 CC_IFFREE (tcyc, int);
1580 CC_IFFREE (degree, int);
1581 CC_IFFREE (tail, int);
1582 CC_IFFREE (marks, char);
1583 CC_IFFREE (deg, int);
1584 CC_IFFREE (adj, int *);
1585 CC_IFFREE (adjspace, int);
1586 CC_IFFREE (perm, int);
1587 CC_IFFREE (elen, int);
1588 return rval;
1589 }
1590
put_in_table(tabledat * td,int i,int j,int * added)1591 static int put_in_table (tabledat *td, int i, int j, int *added)
1592 {
1593 intptr *ip;
1594
1595 if (j < i) {
1596 int temp;
1597 CC_SWAP(i, j, temp);
1598 }
1599
1600 for (ip = td->table[i]; ip; ip = ip->next) {
1601 if (ip->this == j) {
1602 *added = 0;
1603 return 0;
1604 }
1605 }
1606 if (intptr_listadd (&td->table[i], j, &td->intptr_world)) {
1607 *added = 0;
1608 return 1;
1609 }
1610 *added = 1;
1611 return 0;
1612 }
1613
CCedgegen_xnear_build(int ncount,CCdatagroup * dat,double * wcoord,CCxnear * xn)1614 int CCedgegen_xnear_build (int ncount, CCdatagroup *dat, double *wcoord,
1615 CCxnear *xn)
1616 {
1617 int i;
1618 int norm;
1619
1620 xn->nodenames = (int *) NULL;
1621 xn->invnames = (int *) NULL;
1622 xn->w = (double *) NULL;
1623 CCutil_init_datagroup (&xn->dat);
1624
1625 CCutil_dat_getnorm (dat, &norm);
1626 if (CCutil_dat_setnorm (&xn->dat, norm))
1627 return 1;
1628
1629 xn->nodenames = CC_SAFE_MALLOC (ncount, int);
1630 if (!xn->nodenames)
1631 return 1;
1632 for (i = 0; i < ncount; i++)
1633 xn->nodenames[i] = i;
1634 xn->dat.x = CC_SAFE_MALLOC (ncount, double);
1635 if (!xn->dat.x) {
1636 CC_FREE (xn->nodenames, int);
1637 return 1;
1638 }
1639 for (i = 0; i < ncount; i++)
1640 xn->dat.x[i] = dat->x[i];
1641
1642 for (i = 1; i < ncount && dat->x[i] >= dat->x[i - 1]; i++);
1643 if (i < ncount) {
1644 x_quicksort (xn->nodenames, xn->dat.x, 0, ncount - 1);
1645 }
1646
1647 xn->invnames = CC_SAFE_MALLOC (ncount, int);
1648 if (!xn->invnames) {
1649 CC_FREE (xn->nodenames, int);
1650 CCutil_freedatagroup (&(xn->dat));
1651 return 1;
1652 }
1653 for (i = 0; i < ncount; i++)
1654 xn->invnames[xn->nodenames[i]] = i;
1655
1656 xn->dat.y = CC_SAFE_MALLOC (ncount, double);
1657 if (!xn->dat.y) {
1658 CC_FREE (xn->nodenames, int);
1659 CC_FREE (xn->invnames, int);
1660 CCutil_freedatagroup (&(xn->dat));
1661 return 1;
1662 }
1663 for (i = 0; i < ncount; i++)
1664 xn->dat.y[i] = dat->y[xn->nodenames[i]];
1665 if (dat->z != (double *) NULL) {
1666 xn->dat.z = CC_SAFE_MALLOC (ncount, double);
1667 if (!xn->dat.z) {
1668 CC_FREE (xn->nodenames, int);
1669 CC_FREE (xn->invnames, int);
1670 CCutil_freedatagroup (&(xn->dat));
1671 return 1;
1672 }
1673 for (i = 0; i < ncount; i++)
1674 xn->dat.z[i] = dat->z[xn->nodenames[i]];
1675 }
1676 if (wcoord != (double *) NULL) {
1677 xn->w = CC_SAFE_MALLOC (ncount, double);
1678 if (!xn->w) {
1679 CC_FREE (xn->nodenames, int);
1680 CC_FREE (xn->invnames, int);
1681 CCutil_freedatagroup (&(xn->dat));
1682 return 1;
1683 }
1684 for (i = 0; i < ncount; i++)
1685 xn->w[i] = wcoord[xn->nodenames[i]];
1686 }
1687 return 0;
1688 }
1689
CCedgegen_xnear_free(CCxnear * xn)1690 void CCedgegen_xnear_free (CCxnear *xn)
1691 {
1692 CC_IFFREE (xn->nodenames, int);
1693 CC_IFFREE (xn->invnames, int);
1694 CC_IFFREE (xn->w, double);
1695 CCutil_freedatagroup (&(xn->dat));
1696 }
1697
x_quicksort(int * list,double * x,int l,int u)1698 static void x_quicksort (int *list, double *x, int l, int u)
1699 {
1700 int i, j, itemp;
1701 double t, dtemp;
1702
1703 if (l >= u)
1704 return;
1705
1706 CC_SWAP (x[l], x[(l+u)/2], dtemp);
1707 CC_SWAP (list[l], list[(l+u)/2], itemp);
1708
1709 i = l;
1710 j = u + 1;
1711 t = x[l];
1712
1713 while (1) {
1714 do i++; while (i <= u && x[i] < t);
1715 do j--; while (x[j] > t);
1716 if (j < i) break;
1717 CC_SWAP (x[i], x[j], dtemp);
1718 CC_SWAP (list[i], list[j], itemp);
1719 }
1720 CC_SWAP (x[l], x[j], dtemp);
1721 CC_SWAP (list[l], list[j], itemp);
1722 x_quicksort (list, x, l, j - 1);
1723 x_quicksort (list, x, i, u);
1724 }
1725