1 /*===========================================================================*/
2 /* */
3 /* This file is part of a demonstration application for use with the */
4 /* SYMPHONY Branch, Cut, and Price Library. This application is a solver for */
5 /* Capacitated Network Routing Problems. */
6 /* */
7 /* (c) Copyright 2000-2007 Ted Ralphs. All Rights Reserved. */
8 /* */
9 /* This application was developed by Ted Ralphs (ted@lehigh.edu) */
10 /* */
11 /* This software is licensed under the Eclipse Public License. Please see */
12 /* accompanying file for terms. */
13 /* */
14 /*===========================================================================*/
15
16 /* system include files */
17 #include <math.h>
18 #include <memory.h>
19 #include <stddef.h>
20 #include <stdlib.h>
21
22 /* SYMPHONY include files */
23 #include "sym_constants.h"
24 #include "sym_macros.h"
25 #include "sym_qsort.h"
26
27 /* CNRP include files */
28 #include "network.h"
29
30 /*===========================================================================*\
31 * Create the solution graph using adjacency lists
32 \*===========================================================================*/
33
create_net(int * xind,double * xval,int edgenum,double etol,int * edges,double * demand,int vertnum)34 network *create_net(int *xind, double *xval, int edgenum, double etol,
35 int *edges, double *demand, int vertnum)
36 {
37 register edge *net_edges;
38 network *n;
39 vertex *verts;
40 int nv0, nv1;
41 elist *adjlist;
42 int i;
43 #ifdef DIRECTED_X_VARS
44 int total_edgenum = vertnum*(vertnum-1)/2;
45 elist *edge1;
46 char edge_exists;
47 #endif
48
49 /*------------------------------------------------------------------------*\
50 * Allocate the needed memory and set up the data structures
51 \*------------------------------------------------------------------------*/
52
53 n = (network *) calloc (1, sizeof(network));
54 n->vertnum = vertnum;
55 n->edgenum = edgenum;/*the number of edges is equal to the number
56 of nonzeros in the LP solution*/
57 n->verts = (vertex *) calloc(n->vertnum, sizeof(vertex));
58 n->adjlist = (elist *) calloc(2*n->edgenum, sizeof(elist));
59 n->edges = (edge *) calloc(n->edgenum, sizeof(edge));
60 net_edges = n->edges;
61 verts = n->verts;
62 adjlist = n->adjlist;
63 n->is_integral = TRUE;
64
65 /*------------------------------------------------------------------------*\
66 * set up the adjacency list
67 \*------------------------------------------------------------------------*/
68
69 for (i = 0; i < edgenum; i++, xval++, xind++){
70 if (*xval < etol) continue;
71 #ifdef DIRECTED_X_VARS
72 if (*xind < total_edgenum){
73 nv0 = edges[(*xind) << 1];
74 nv1 = edges[((*xind)<< 1) + 1];
75 }else{
76 nv0 = edges[(*xind-total_edgenum) << 1];
77 nv1 = edges[((*xind-total_edgenum)<< 1) + 1];
78 }
79 for (edge_exists = FALSE, edge1 = verts[nv0].first; edge1;
80 edge1 = edge1->next_edge){
81 if (edge1->other_end == nv1){
82 if (fabs(floor(*xval+.5) - *xval) > etol){
83 edge1->data->weight += *xval;
84 }else{
85 edge1->data->weight += floor(*xval+.5);
86 }
87 edge_exists = TRUE;
88 break;
89 }
90 }
91 if (edge_exists)
92 continue;
93 #else
94 nv0 = edges[(*xind) << 1];
95 nv1 = edges[((*xind)<< 1) + 1];
96 #endif
97
98 if (fabs(floor(*xval+.5) - *xval) > etol){
99 net_edges->weight += *xval;
100 }else{
101 net_edges->weight = floor(*xval+.5);
102 }
103 net_edges->v0 = nv0;
104 net_edges->v1 = nv1;
105 if (!verts[nv0].first){
106 verts[nv0].first = verts[nv0].last = adjlist;
107 verts[nv0].degree++;
108 }
109 else{
110 verts[nv0].last->next_edge = adjlist;
111 verts[nv0].last = adjlist;
112 verts[nv0].degree++;
113 }
114 adjlist->data = net_edges;
115 adjlist->other_end = nv1;
116 adjlist->other = verts + nv1;
117 adjlist++;
118 if (!verts[nv1].first){
119 verts[nv1].first = verts[nv1].last = adjlist;
120 verts[nv1].degree++;
121 }
122 else{
123 verts[nv1].last->next_edge = adjlist;
124 verts[nv1].last = adjlist;
125 verts[nv1].degree++;
126 }
127 adjlist->data = net_edges;
128 adjlist->other_end = nv0;
129 adjlist->other = verts + nv0;
130 adjlist++;
131
132 net_edges++;
133 }
134
135 for (i = 0; i < edgenum; i++){
136 if (fabs(floor(n->edges[i].weight + .5) - n->edges[i].weight) > etol){
137 n->is_integral = FALSE;
138 break;
139 }
140 }
141
142 /*set the demand for each node*/
143 for (i = 0; i < vertnum; i++){
144 verts[i].demand = demand[i];
145 verts[i].orignodenum = i;
146 }
147
148 return(n);
149 }
150
151 /*===========================================================================*\
152 * Create the solution graph using adjacency lists
153 \*===========================================================================*/
154
create_flow_net(int * xind,double * xval,int edgenum,double etol,int * edges,double * demand,int vertnum)155 network *create_flow_net(int *xind, double *xval, int edgenum, double etol,
156 int *edges, double *demand, int vertnum)
157 {
158 register edge *net_edges;
159 network *n;
160 vertex *verts;
161 int nv0, nv1;
162 elist *adjlist;
163 int i = 0;
164 #if defined(DIRECTED_X_VARS) || defined(ADD_FLOW_VARS)
165 int total_edgenum = vertnum*(vertnum-1)/2;
166 char edge_exists, flow_var;
167 elist *edge1;
168 #endif
169 #if defined(DIRECTED_X_VARS) && defined(ADD_FLOW_VARS)
170 int j;
171 char d_x_vars = TRUE;
172 #elif defined(ADD_FLOW_VARS)
173 int j;
174 char d_x_vars = FALSE;
175 #endif
176
177 /*------------------------------------------------------------------------*\
178 * Allocate the needed memory and set up the data structures
179 \*------------------------------------------------------------------------*/
180
181 n = (network *) calloc (1, sizeof(network));
182 n->vertnum = vertnum;
183 n->verts = (vertex *) calloc(vertnum, sizeof(vertex));
184 n->adjlist = (elist *) calloc(2*edgenum, sizeof(elist));
185 n->edges = (edge *) calloc(edgenum, sizeof(edge));
186 net_edges = n->edges;
187 verts = n->verts;
188 adjlist = n->adjlist;
189 n->is_integral = TRUE;
190 n->edgenum = 0;
191
192 /*------------------------------------------------------------------------*\
193 * set up the adjacency list
194 \*------------------------------------------------------------------------*/
195
196 for (i = 0; i < edgenum; i++){
197 #if !defined(DIRECTED_X_VARS) && defined(ADD_FLOW_VARS)
198 if (xind[i] < total_edgenum){
199 nv0 = edges[xind[i] << 1];
200 nv1 = edges[(xind[i] << 1) + 1];
201 flow_var = FALSE;
202 }else if (xind[i] < 2 * total_edgenum){
203 nv0 = edges[(xind[i] - total_edgenum) << 1];
204 nv1 = edges[((xind[i] - total_edgenum) << 1) + 1];
205 flow_var = TRUE;
206 }else if (xind[i] < 3 * total_edgenum){
207 nv0 = edges[((xind[i] - 2*total_edgenum) << 1) + 1];
208 nv1 = edges[(xind[i] - 2*total_edgenum) << 1];
209 flow_var = TRUE;
210 }else{
211 continue;
212 }
213
214 #elif defined(DIRECTED_X_VARS)
215 if (xind[i] < total_edgenum){
216 nv0 = edges[xind[i] << 1];
217 nv1 = edges[(xind[i] << 1) + 1];
218 flow_var = FALSE;
219 }else if (xind[i] < 2 * total_edgenum){
220 nv0 = edges[((xind[i] - total_edgenum) << 1) + 1];
221 nv1 = edges[(xind[i] - total_edgenum) << 1];
222 flow_var = FALSE;
223 }else if (xind[i] < 3 * total_edgenum){
224 nv0 = edges[(xind[i] - 2 * total_edgenum) << 1];
225 nv1 = edges[((xind[i] - 2 * total_edgenum) << 1) + 1];
226 flow_var = TRUE;
227 }else if (xind[i] < 4 * total_edgenum){
228 nv0 = edges[((xind[i] - 3 * total_edgenum) << 1) + 1];
229 nv1 = edges[(xind[i] - 3 * total_edgenum) << 1];
230 flow_var = TRUE;
231 }else{
232 continue;
233 }
234 #endif
235
236 #ifdef ADD_FLOW_VARS
237 if (flow_var){
238 for (edge_exists = FALSE, edge1 = verts[nv0].first; edge1;
239 edge1 = edge1->next_edge){
240 if (edge1->other_end == nv1){
241 if (nv0 < nv1){
242 edge1->data->flow1 = xval[i];
243 }else{
244 edge1->data->flow2 = xval[i];
245 }
246 edge_exists = TRUE;
247 break;
248 }
249 }
250 if (edge_exists){
251 continue;
252 }else{
253 if (nv0 < nv1){
254 net_edges->flow1 = xval[i];
255 }else{
256 net_edges->flow2 = xval[i];
257 }
258 }
259 }else{
260 #ifdef DIRECTED_X_VARS
261 for (edge_exists = FALSE, edge1 = verts[nv0].first; edge1;
262 edge1 = edge1->next_edge){
263 if (edge1->other_end == nv1){
264 if (nv0 < nv1){
265 edge1->data->weight += xval[i];
266 edge1->data->weight1 = xval[i];
267 }else{
268 edge1->data->weight += xval[i];
269 edge1->data->weight2 = xval[i];
270 }
271 edge_exists = TRUE;
272 break;
273 }
274 }
275 if (edge_exists){
276 continue;
277 }else{
278 if (nv0 < nv1){
279 net_edges->weight += xval[i];
280 net_edges->weight1 = xval[i];
281 }else{
282 net_edges->weight += xval[i];
283 net_edges->weight2 = xval[i];
284 }
285 }
286 #else
287 net_edges->weight = xval[i];
288 #endif
289 }
290 #else
291 net_edges->weight = xval[i];
292 #endif
293
294 net_edges->v0 = nv0 < nv1 ? nv0 : nv1;
295 net_edges->v1 = nv0 < nv1 ? nv1 : nv0;
296 if (!verts[nv0].first){
297 verts[nv0].first = verts[nv0].last = adjlist;
298 verts[nv0].degree++;
299 }
300 else{
301 verts[nv0].last->next_edge = adjlist;
302 verts[nv0].last = adjlist;
303 verts[nv0].degree++;
304 }
305 adjlist->data = net_edges;
306 adjlist->other_end = nv1;
307 adjlist->other = verts + nv1;
308 adjlist++;
309 if (!verts[nv1].first){
310 verts[nv1].first = verts[nv1].last = adjlist;
311 verts[nv1].degree++;
312 }
313 else{
314 verts[nv1].last->next_edge = adjlist;
315 verts[nv1].last = adjlist;
316 verts[nv1].degree++;
317 }
318 adjlist->data = net_edges;
319 adjlist->other_end = nv0;
320 adjlist->other = verts + nv0;
321 adjlist++;
322
323 net_edges++;
324 n->edgenum++;
325 }
326
327 for (i = 0; i < edgenum; i++){
328 if (fabs(floor(n->edges[i].weight + .5) - n->edges[i].weight) > etol){
329 n->is_integral = FALSE;
330 break;
331 }
332 }
333
334 /*set the demand for each node*/
335 if (demand){
336 for (i = 0; i < vertnum; i++){
337 verts[i].demand = demand[i];
338 verts[i].orignodenum = i;
339 }
340 }
341
342 return(n);
343 }
344
345 /*===========================================================================*/
346
347 /*===========================================================================*\
348 * Calculates the connected components of the solution graph after removing
349 * the depot. Each node is assigned the number of the component in which it
350 * resides. The number of nodes in each component is put in "compnodes", the
351 * total demand of all customers in the component is put in "compdemands", and
352 * the value of the cut induced by the component is put in "compcuts".
353 \*===========================================================================*/
354
connected(network * n,int * compnodes,double * compdemands,int * compmembers,double * compcuts,double * compdensity)355 int connected(network *n, int *compnodes, double *compdemands,
356 int *compmembers, double *compcuts, double *compdensity)
357 {
358 int cur_node = 0, cur_comp = 0;
359 int cur_member = 0;
360 vertex *verts = n->verts;
361 int vertnum = n->vertnum;
362 elist *cur_edge;
363 int *nodes_to_scan, num_nodes_to_scan = 0, i;
364 char *is_not_integral = NULL;
365
366 nodes_to_scan = (int *) calloc (vertnum, sizeof(int));
367 if (compdensity)
368 is_not_integral = (char *) calloc (vertnum, sizeof(char));
369
370 while (TRUE){
371 for (cur_node = 1; cur_node < vertnum; cur_node++)
372 if (!verts[cur_node].comp){/*look for a node that hasn't been assigned
373 to a component yet*/
374 break;
375 }
376
377 if (cur_node == n->vertnum) break;/*this indicates that all nodes have been
378 assigned to components*/
379
380 nodes_to_scan[num_nodes_to_scan++] = cur_node;/*add the first node to the
381 list of nodes to be scanned*/
382
383 if (compmembers) compmembers[++cur_member] = cur_node;
384
385 verts[cur_node].comp = ++cur_comp;/*add the first node into the new
386 component*/
387 compnodes[cur_comp] = 1;
388 verts[cur_node].comp = cur_comp;
389 compdemands[cur_comp] = verts[cur_node].demand;
390 if (compdensity){
391 compdensity[cur_comp] = verts[cur_node].degree;
392 if (verts[cur_node].degree != 2)
393 is_not_integral[cur_comp] = TRUE;
394 }
395 while(TRUE){/*continue to execute this loop until there are no more
396 nodes to scan if there is a node to scan, then add all of
397 its neighbors in to the current component and take it off
398 the list*/
399 for (cur_node = nodes_to_scan[--num_nodes_to_scan],
400 verts[cur_node].scanned = TRUE,
401 cur_edge = verts[cur_node].first, cur_comp = verts[cur_node].comp;
402 cur_edge; cur_edge = cur_edge->next_edge){
403 if (cur_edge->other_end){
404 if (!verts[cur_edge->other_end].comp){
405 verts[cur_edge->other_end].comp = cur_comp;
406 compnodes[cur_comp]++;
407 if (compmembers) compmembers[++cur_member] = cur_edge->other_end;
408 compdemands[cur_comp] += verts[cur_edge->other_end].demand;
409 if (compdensity){
410 compdensity[cur_comp] += verts[cur_edge->other_end].degree;
411 if (verts[cur_edge->other_end].degree != 2)
412 is_not_integral[cur_comp] = TRUE;
413 }
414 nodes_to_scan[num_nodes_to_scan++] = cur_edge->other_end;
415 }
416 }
417 else{/*if this node is connected to the depot, then
418 update the value of the cut*/
419 if (compcuts) compcuts[cur_comp] += cur_edge->data->weight;
420 if (compdensity) compdensity[cur_comp] += 1;
421 }
422 }
423 if (!num_nodes_to_scan) break;
424 }
425 }
426
427 if (compdensity){
428 for (i = 1; i <= cur_comp; i++){
429 if (is_not_integral[i])
430 compdensity[i] /= 2*(compnodes[i]+1);
431 else
432 compdensity[i] = MAXDOUBLE;
433 }
434 }
435
436 FREE(nodes_to_scan);
437 FREE(is_not_integral);
438
439 return(cur_comp);
440 }
441
442 /*===========================================================================*/
443
flow_connected(network * n,int * compnodes,double * compdemands,int * compmembers,double * compcuts,double * compdensity,double etol)444 int flow_connected(network *n, int *compnodes, double *compdemands,
445 int *compmembers, double *compcuts, double *compdensity,
446 double etol)
447 {
448 int cur_node = 0, cur_comp = 0;
449 int cur_member = 0;
450 vertex *verts = n->verts;
451 int vertnum = n->vertnum;
452 elist *cur_edge;
453 int *nodes_to_scan, num_nodes_to_scan = 0, i;
454 char *is_not_integral = NULL;
455
456 nodes_to_scan = (int *) calloc (vertnum, sizeof(int));
457 if (compdensity)
458 is_not_integral = (char *) calloc (vertnum, sizeof(char));
459
460 while (TRUE){
461 for (cur_node = 1; cur_node < vertnum; cur_node++)
462 if (!verts[cur_node].comp){/*look for a node that hasn't been assigned
463 to a component yet*/
464 break;
465 }
466
467 if (cur_node == n->vertnum) break;/*this indicates that all nodes have been
468 assigned to components*/
469
470 nodes_to_scan[num_nodes_to_scan++] = cur_node;/*add the first node to the
471 list of nodes to be scanned*/
472
473 if (compmembers) compmembers[++cur_member] = cur_node;
474
475 verts[cur_node].comp = ++cur_comp;/*add the first node into the new
476 component*/
477 compnodes[cur_comp] = 1;
478 verts[cur_node].comp = cur_comp;
479 compdemands[cur_comp] = verts[cur_node].demand;
480 if (compdensity){
481 compdensity[cur_comp] = verts[cur_node].degree;
482 if (verts[cur_node].degree != 2)
483 is_not_integral[cur_comp] = TRUE;
484 }
485 while(TRUE){/*continue to execute this loop until there are no more
486 nodes to scan if there is a node to scan, then add all of
487 its neighbors in to the current component and take it off
488 the list*/
489 for (cur_node = nodes_to_scan[--num_nodes_to_scan],
490 verts[cur_node].scanned = TRUE,
491 cur_edge = verts[cur_node].first, cur_comp = verts[cur_node].comp;
492 cur_edge; cur_edge = cur_edge->next_edge){
493 if (cur_edge->data->weight < etol) continue;
494 if (cur_edge->other_end){
495 if (!verts[cur_edge->other_end].comp){
496 verts[cur_edge->other_end].comp = cur_comp;
497 compnodes[cur_comp]++;
498 if (compmembers) compmembers[++cur_member]=cur_edge->other_end;
499 compdemands[cur_comp] += verts[cur_edge->other_end].demand;
500 if (compdensity){
501 compdensity[cur_comp] += verts[cur_edge->other_end].degree;
502 if (verts[cur_edge->other_end].degree != 2)
503 is_not_integral[cur_comp] = TRUE;
504 }
505 nodes_to_scan[num_nodes_to_scan++] = cur_edge->other_end;
506 }
507 }
508 else{/*if this node is connected to the depot, then
509 update the value of the cut*/
510 if (compcuts) compcuts[cur_comp] += cur_edge->data->weight;
511 if (compdensity) compdensity[cur_comp] += 1;
512 }
513 }
514 if (!num_nodes_to_scan) break;
515 }
516 }
517
518 if (compdensity){
519 for (i = 1; i <= cur_comp; i++){
520 if (is_not_integral[i])
521 compdensity[i] /= 2*(compnodes[i]+1);
522 else
523 compdensity[i] = MAXDOUBLE;
524 }
525 }
526
527 FREE(nodes_to_scan);
528 FREE(is_not_integral);
529
530 return(cur_comp);
531 }
532
533 /*===========================================================================*/
534
535 /*===========================================================================*\
536 * Free the memory associated with the solution graph data structures
537 \*===========================================================================*/
538
free_net(network * n)539 void free_net(network *n)
540 {
541 int i;
542 if (n){
543 if (n->adjlist) free ((char *) n->adjlist);
544 if (n->verts){
545 for (i = 0; i < n->vertnum; i++)
546 if (n->verts[i].orig_node_list)
547 free((char *)n->verts[i].orig_node_list);
548 free ((char *) n->verts);
549 }
550 if (n->edges) free((char *) n->edges);
551 free((char *) n);
552 }
553 }
554