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 /* the Vehicle Routing Problem and the Traveling Salesman Problem. */
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 /* VRP include files */
28 #include "network.h"
29
30 /*===========================================================================*\
31 * Create the solution graph using adjacency lists
32 \*===========================================================================*/
33
createnet(int * xind,double * xval,int edgenum,double etol,int * edges,int * demand,int vertnum)34 network *createnet(int *xind, double *xval, int edgenum, double etol,
35 int *edges, int *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
44 double *val_low, *val_high, val_aux;
45 int *ind_low, *ind_high, ind_aux;
46
47 /*------------------------------------------------------------------------*\
48 * Allocate the needed memory and set up the data structures
49 \*------------------------------------------------------------------------*/
50
51 n = (network *) calloc (1, sizeof(network));
52 n->vertnum = vertnum;
53 n->edgenum = edgenum;/*the number of edges is equal to the number
54 of nonzeros in the LP solution*/
55 n->verts = (vertex *) calloc(n->vertnum, sizeof(vertex));
56 n->adjlist = (elist *) calloc(2*n->edgenum, sizeof(elist));
57 n->edges = (edge *) calloc(n->edgenum, sizeof(edge));
58 net_edges = n->edges;
59 verts = n->verts;
60 adjlist = n->adjlist;
61 n->is_integral = TRUE;
62
63 qsort_di(xval, xind, edgenum);
64 /* qsort_di sorts the array in nondecreasing order;
65 now need to translate it into nonincreasing */
66
67 for (i = 0, val_low = xval, val_high = xval + edgenum - 1, ind_low = xind,
68 ind_high = xind + edgenum - 1; i < (edgenum/2); i++){
69 val_aux = *val_low;
70 *val_low++ = *val_high;
71 *val_high-- = val_aux;
72 ind_aux = *ind_low;
73 *ind_low++ = *ind_high;
74 *ind_high-- = ind_aux;
75 }
76
77
78 /*------------------------------------------------------------------------*\
79 * set up the adjacency list
80 \*------------------------------------------------------------------------*/
81
82 for (i = 0; i < edgenum; i++, xval++, xind++){
83 if (*xval < etol) continue;
84 if (fabs(floor(*xval+.5) - *xval) > etol){
85 n->is_integral = FALSE;
86 net_edges->weight = *xval;
87 }else{
88 net_edges->weight = floor(*xval+.5);
89 }
90 nv0 = net_edges->v0 = edges[(*xind) << 1];
91 nv1 = net_edges->v1 = edges[((*xind)<< 1) + 1];
92 if (!verts[nv0].first){
93 verts[nv0].first = verts[nv0].last = adjlist;
94 verts[nv0].degree++;
95 }
96 else{
97 verts[nv0].last->next_edge = adjlist;
98 verts[nv0].last = adjlist;
99 verts[nv0].degree++;
100 }
101 adjlist->data = net_edges;
102 adjlist->other_end = nv1;
103 adjlist->other = verts + nv1;
104 adjlist++;
105 if (!verts[nv1].first){
106 verts[nv1].first = verts[nv1].last = adjlist;
107 verts[nv1].degree++;
108 }
109 else{
110 verts[nv1].last->next_edge = adjlist;
111 verts[nv1].last = adjlist;
112 verts[nv1].degree++;
113 }
114 adjlist->data = net_edges;
115 adjlist->other_end = nv0;
116 adjlist->other = verts + nv0;
117 adjlist++;
118
119 net_edges++;
120 }
121
122 /*set the demand for each node*/
123 for (i = 0; i < vertnum; i++){
124 verts[i].demand = demand[i];
125 verts[i].orignodenum = i;
126 }
127
128 /*__BEGIN_EXPERIMENTAL_SECTION__*/
129 #if 0
130 /*allocate memory for existing nodes list and the binary tree used by
131 capforest*/
132 n->enodes = (vertex **) calloc (n->vertnum, sizeof(vertex *));
133 n->tnodes = (vertex **) calloc (n->vertnum, sizeof(vertex *));
134 #endif
135
136 /*___END_EXPERIMENTAL_SECTION___*/
137 return(n);
138 }
139
140 /*__BEGIN_EXPERIMENTAL_SECTION__*/
141 /*===========================================================================*/
142
createnet2(int * xind,double * xval,int edgenum,double etol,int * edges,int * demand,int vertnum,char * status)143 network *createnet2(int *xind, double *xval, int edgenum, double etol,
144 int *edges, int *demand, int vertnum, char *status)
145 {
146 register edge *net_edges;
147 network *n;
148 vertex *verts;
149 int nv0, nv1;
150 elist *adjlist;
151 int i;
152 char *stat = status;
153 int *sort_order;
154 double *tmp;
155
156 double *val_low, *val_high, val_aux;
157 int *ind_low, *ind_high, ind_aux;
158
159 /*------------------------------------------------------------------------*\
160 * Allocate the needed memory and set up the data structures
161 \*------------------------------------------------------------------------*/
162
163 n = (network *) calloc (1, sizeof(network));
164 n->vertnum = vertnum;
165 n->edgenum = edgenum;/*the number of edges is equal to the number
166 of nonzeros in the LP solution*/
167 n->verts = (vertex *) calloc(n->vertnum, sizeof(vertex));
168 n->adjlist = (elist *) calloc(2*n->edgenum, sizeof(elist));
169 n->edges = (edge *) calloc(n->edgenum, sizeof(edge));
170 net_edges = n->edges;
171 verts = n->verts;
172 adjlist = n->adjlist;
173 n->is_integral = TRUE;
174
175 tmp = (double *) malloc(edgenum*DSIZE);
176 sort_order = (int *) malloc(edgenum*ISIZE);
177 for (i = edgenum - 1; i >= 0; i--)
178 sort_order[edgenum - i -1] = i;
179 memcpy(tmp, xval, edgenum * DSIZE);
180
181 qsort_di(tmp, sort_order, edgenum);
182 qsort_ic(sort_order, status, edgenum);
183 FREE(tmp);
184 FREE(sort_order);
185 qsort_di(xval, xind, edgenum);
186 /* qsort_di sorts the array in nondecreasing order;
187 now need to translate it into nonincreasing */
188
189 for (i = 0, val_low = xval, val_high = xval + edgenum - 1, ind_low = xind,
190 ind_high = xind + edgenum - 1; i < (edgenum/2); i++){
191 val_aux = *val_low;
192 *val_low++ = *val_high;
193 *val_high-- = val_aux;
194 ind_aux = *ind_low;
195 *ind_low++ = *ind_high;
196 *ind_high-- = ind_aux;
197 }
198
199
200 /*------------------------------------------------------------------------*\
201 * set up the adjacency list
202 \*------------------------------------------------------------------------*/
203
204 for (i = 0; i < edgenum; i++, xval++, xind++, stat++){
205 if (*xval < etol) continue;
206 if (fabs(floor(*xval+.5) - *xval) > etol){
207 n->is_integral = FALSE;
208 net_edges->weight = *xval;
209 }else{
210 net_edges->weight = floor(*xval+.5);
211 }
212 #ifdef COMPILE_OUR_DECOMP
213 net_edges->status = *stat;
214 #endif
215 nv0 = net_edges->v0 = edges[(*xind) << 1];
216 nv1 = net_edges->v1 = edges[((*xind)<< 1) + 1];
217 if (!verts[nv0].first){
218 verts[nv0].first = verts[nv0].last = adjlist;
219 verts[nv0].degree++;
220 }
221 else{
222 verts[nv0].last->next_edge = adjlist;
223 verts[nv0].last = adjlist;
224 verts[nv0].degree++;
225 }
226 adjlist->data = net_edges;
227 adjlist->other_end = nv1;
228 adjlist->other = verts + nv1;
229 adjlist++;
230 if (!verts[nv1].first){
231 verts[nv1].first = verts[nv1].last = adjlist;
232 verts[nv1].degree++;
233 }
234 else{
235 verts[nv1].last->next_edge = adjlist;
236 verts[nv1].last = adjlist;
237 verts[nv1].degree++;
238 }
239 adjlist->data = net_edges;
240 adjlist->other_end = nv0;
241 adjlist->other = verts + nv0;
242 adjlist++;
243
244 net_edges++;
245 }
246
247 /*set the demand for each node*/
248 for (i = 0; i < vertnum; i++){
249 verts[i].demand = demand[i];
250 verts[i].orignodenum = i;
251 }
252
253 #if 0
254 /*allocate memory for existing nodes list and the binary tree used by
255 capforest*/
256 n->enodes = (vertex **) calloc (n->vertnum, sizeof(vertex *));
257 n->tnodes = (vertex **) calloc (n->vertnum, sizeof(vertex *));
258 #endif
259
260 return(n);
261 }
262
263 /*___END_EXPERIMENTAL_SECTION___*/
264 /*===========================================================================*/
265
266 /*===========================================================================*\
267 * Calculates the connected components of the solution graph after removing
268 * the depot. Each node is assigned the number of the component in which it
269 * resides. The number of nodes in each component is put in "compnodes", the
270 * total demand of all customers in the component is put in "compdemands", and
271 * the value of the cut induced by the component is put in "compcuts".
272 \*===========================================================================*/
273
connected(network * n,int * compnodes,int * compdemands,int * compmembers,double * compcuts,double * compdensity)274 int connected(network *n, int *compnodes, int *compdemands, int *compmembers,
275 double *compcuts, double *compdensity)
276 {
277 int cur_node = 0, cur_comp = 0;
278 int cur_member = 0;
279 vertex *verts = n->verts;
280 int vertnum = n->vertnum;
281 elist *cur_edge;
282 int *nodes_to_scan, num_nodes_to_scan = 0, i;
283 char *is_not_integral = NULL;
284
285 nodes_to_scan = (int *) calloc (vertnum, sizeof(int));
286 if (compdensity)
287 is_not_integral = (char *) calloc (vertnum, sizeof(char));
288
289 while (TRUE){
290 for (cur_node = 1; cur_node < vertnum; cur_node++)
291 if (!verts[cur_node].comp){/*look for a node that hasn't been assigned
292 to a component yet*/
293 break;
294 }
295
296 if (cur_node == n->vertnum) break;/*this indicates that all nodes have been
297 assigned to components*/
298
299 nodes_to_scan[num_nodes_to_scan++] = cur_node;/*add the first node to the
300 list of nodes to be scanned*/
301
302 if (compmembers) compmembers[++cur_member] = cur_node;
303
304 verts[cur_node].comp = ++cur_comp;/*add the first node into the new
305 component*/
306 compnodes[cur_comp] = 1;
307 verts[cur_node].comp = cur_comp;
308 compdemands[cur_comp] = verts[cur_node].demand;
309 if (compdensity){
310 compdensity[cur_comp] = verts[cur_node].degree;
311 if (verts[cur_node].degree != 2)
312 is_not_integral[cur_comp] = TRUE;
313 }
314 while(TRUE){/*continue to execute this loop until there are no more
315 nodes to scan if there is a node to scan, then add all of
316 its neighbors in to the current component and take it off
317 the list*/
318 for (cur_node = nodes_to_scan[--num_nodes_to_scan],
319 verts[cur_node].scanned = TRUE,
320 cur_edge = verts[cur_node].first, cur_comp = verts[cur_node].comp;
321 cur_edge; cur_edge = cur_edge->next_edge){
322 if (cur_edge->other_end){
323 if (!verts[cur_edge->other_end].comp){
324 verts[cur_edge->other_end].comp = cur_comp;
325 compnodes[cur_comp]++;
326 if (compmembers) compmembers[++cur_member] = cur_edge->other_end;
327 compdemands[cur_comp] += verts[cur_edge->other_end].demand;
328 if (compdensity){
329 compdensity[cur_comp] += verts[cur_edge->other_end].degree;
330 if (verts[cur_edge->other_end].degree != 2)
331 is_not_integral[cur_comp] = TRUE;
332 }
333 nodes_to_scan[num_nodes_to_scan++] = cur_edge->other_end;
334 }
335 }
336 else{/*if this node is connected to the depot, then
337 update the value of the cut*/
338 if (compcuts) compcuts[cur_comp] += cur_edge->data->weight;
339 if (compdensity) compdensity[cur_comp] += 1;
340 }
341 }
342 if (!num_nodes_to_scan) break;
343 }
344 }
345
346 if (compdensity){
347 for (i = 1; i <= cur_comp; i++){
348 if (is_not_integral[i])
349 compdensity[i] /= 2*(compnodes[i]+1);
350 else
351 compdensity[i] = MAXDOUBLE;
352 }
353 }
354
355 FREE(nodes_to_scan);
356 FREE(is_not_integral);
357
358 return(cur_comp);
359 }
360
361 /*===========================================================================*/
362
363 /*===========================================================================*\
364 * Free the memory associated with the solution graph data structures
365 \*===========================================================================*/
366
free_net(network * n)367 void free_net(network *n)
368 {
369 int i;
370 if (n){
371 if (n->adjlist) free ((char *) n->adjlist);
372 if (n->verts){
373 for (i = 0; i < n->vertnum; i++)
374 if (n->verts[i].orig_node_list)
375 free((char *)n->verts[i].orig_node_list);
376 free ((char *) n->verts);
377 }
378 if (n->edges) free((char *) n->edges);
379 /*__BEGIN_EXPERIMENTAL_SECTION__*/
380 #if 0
381 if (n->tnodes) free((char *)n->tnodes);
382 if (n->enodes) free((char *)n->enodes);
383 #endif
384 /*___END_EXPERIMENTAL_SECTION___*/
385 free((char *) n);
386 }
387 }
388