1 /*===========================================================================*
2 * This file is part of a solver for the Vehicle Routing Problem *
3 * developed using the BiCePS Linear Integer Solver (BLIS). *
4 * *
5 * This solver is distributed under the Eclipse Public License as part of *
6 * the COIN-OR repository (http://www.coin-or.org). *
7 * *
8 * Authors: Yan Xu, Lehigh University *
9 * Ted Ralphs, Lehigh University *
10 * *
11 * Copyright (C) 2007 Yan Xu and Ted Ralphs. *
12 * All Rights Reserved. *
13 *===========================================================================*/
14
15 #include "BcpsObjectPool.h"
16 #include "BlisConstraint.h"
17 #include "BlisHelp.h"
18
19 #include "VrpModel.h"
20 #include "VrpCutGenerator.h"
21 #include "VrpMacros.h"
22 #include "VrpParams.h"
23
24 #define DELETE_POWER 3
25 #define DELETE_AND 0x07
26
27 /*===========================================================================*/
28
VrpCutGenerator(VrpModel * vrp,int vertnum)29 VrpCutGenerator::VrpCutGenerator(VrpModel *vrp, int vertnum)
30 {
31 model_ = vrp;
32 if (vertnum){
33 ref_ = new int[model_->vertnum_];
34 cutVal_ = new double[model_->vertnum_];
35 cutList_ = new char[((model_->vertnum_ >> DELETE_POWER) + 1)*
36 (model_->VrpPar_->entry(VrpParams::maxNumCutsInShrink) + 1)];
37 inSet_ = new char[model_->vertnum_];
38 }else{
39 ref_ = 0;
40 cutVal_ = 0;
41 cutList_ = 0;
42 inSet_ = 0;
43 }
44 SRANDOM(1);
45 setName("VRP");
46 }
47
48 /*===========================================================================*/
49
50 // Return if need resolve LP immediately.
51 // New constraints are stored in BcpsConstraintPool
52 bool
generateConstraints(BcpsConstraintPool & conPool)53 VrpCutGenerator::generateConstraints(BcpsConstraintPool &conPool)
54 {
55 int vertnum = model_->vertnum_;
56 int rcnt, cur_bins = 0, i, k, max_node;
57 double cur_slack = 0.0, node_cut, max_node_cut;
58 int cut_size = (vertnum >> DELETE_POWER) + 1, num_cuts = 0;
59 elist *cur_edge = NULL;
60 VrpParams *par = model_->VrpPar_;
61 int which_connected_routine = par->entry(VrpParams::whichConnectedRoutine);
62 bool do_greedy = par->entry(VrpParams::doGreedy);
63 double etol = model_->etol_;
64
65 int *demand = model_->demand_;
66 elist *cur_edge1 = NULL, *cur_edge2 = NULL;
67 int node1 = 0, node2 = 0;
68 int total_demand = demand[0], num_trials = 0;
69 int type, rhs, capacity = model_->capacity_;
70 VrpNetwork *n = model_->n_;
71
72 // Get dense solution
73 const double *denseSol = model_->getLpSolution();
74 // Transform it to a sparse vector.
75 CoinPackedVector *sol = model_->getSolution(denseSol);
76 model_->createNet(sol);
77
78 if (n->isIntegral_){
79 /* if the network is integral, check for connectivity */
80 n->connected();
81 delete sol;
82 return connectivityCuts(conPool) ? true: false;
83 }
84
85 #ifdef DO_TSP_CUTS
86 if (par->tspProb){
87 delete sol;
88 return tspCuts(model_, conPool) ? true : false;
89 }
90 #endif
91
92 vertex *verts = n->verts_;
93
94 if (which_connected_routine == BOTH) which_connected_routine = CONNECTED;
95
96 int *compnodes_copy = new int[vertnum + 1];
97 int *compnodes = n->compNodes_;
98 double *compcuts = n->compCuts_;
99 int *compdemands = n->compDemands_;
100
101 do{
102 /*------------------------------------------------------------------*\
103 * Get the connected components of the solution graph without the
104 * depot and see if the number of components is more than one
105 \*------------------------------------------------------------------*/
106 rcnt = which_connected_routine == BICONNECTED ?
107 n->biconnected() : n->connected();
108
109 /* copy the arrays as they will be needed later */
110 if (!which_connected_routine && do_greedy){
111 compnodes_copy = (int *) memcpy((char *)compnodes_copy,
112 (char*)compnodes,
113 (vertnum + 1)*sizeof(int));
114 compnodes = compnodes_copy;
115 }
116 if (rcnt > 1){
117 /*---------------------------------------------------------------*\
118 * If the number of components is more then one, then check each
119 * component to see if it violates a capacity constraint
120 \*---------------------------------------------------------------*/
121
122 char **coef_list = new char *[rcnt];
123 memset(coef_list, 0, rcnt*sizeof(char *));
124 coef_list[0] = new char[rcnt * cut_size];
125 memset(coef_list[0], 0, rcnt*cut_size*sizeof(char));
126 for(i = 1; i < rcnt; i++)
127 coef_list[i] = coef_list[0]+i*cut_size;
128
129 for(i = 1; i < vertnum; i++)
130 (coef_list[(verts[i].comp)-1][i >> DELETE_POWER]) |=
131 (1 << (i & DELETE_AND));
132
133 for (i = 0; i < rcnt; i++){
134 if (compnodes[i+1] < 2) continue;
135 /*check ith component to see if it violates a constraint*/
136 if (which_connected_routine == BOTH &&
137 which_connected_routine == BICONNECTED && compcuts[i+1]==0)
138 continue;
139 if (compcuts[i+1] < 2*BINS(compdemands[i+1], capacity)-etol){
140 /*the constraint is violated so impose it*/
141 type = (compnodes[i+1] < vertnum/2 ?
142 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
143 rhs = (type == SUBTOUR_ELIM_SIDE ?
144 RHS(compnodes[i+1],compdemands[i+1],
145 capacity): 2*BINS(compdemands[i+1],capacity));
146 num_cuts += addVrpCut(conPool, coef_list[i], rhs, type);
147 }
148 else{/*if the constraint is not violated, then try generating a
149 violated constraint by deleting customers that don't
150 change the number of trucks required by the customers in
151 the component but decrease the value of the cut*/
152 cur_bins = BINS(compdemands[i+1], capacity);/*the current
153 number of trucks
154 required*/
155 /*current slack in the constraint*/
156 cur_slack = (compcuts[i+1] - 2*cur_bins);
157 while (compnodes[i+1]){/*while there are still nodes in the
158 component*/
159 for (max_node = 0, max_node_cut = 0, k = 1;
160 k < vertnum; k++){
161 if (verts[k].comp == i+1){
162 if (BINS(compdemands[i+1]-verts[k].demand, capacity)
163 == cur_bins){
164 /*if the number of trucks doesn't decrease upon
165 deleting this customer*/
166 for (node_cut = 0, cur_edge = verts[k].first;
167 cur_edge; cur_edge = cur_edge->next_edge){
168 node_cut += (cur_edge->other_end ?
169 -cur_edge->data->weight :
170 cur_edge->data->weight);
171 }
172 if (node_cut > max_node_cut){/*check whether the
173 value of the cut
174 decrease is the best
175 seen so far*/
176 max_node = k;
177 max_node_cut = node_cut;
178 }
179 }
180 }
181 }
182 if (!max_node){
183 break;
184 }
185 /*delete the customer that exhibited the greatest
186 decrease in cut value*/
187 compnodes[i+1]--;
188 compdemands[i+1] -= verts[max_node].demand;
189 compcuts[i+1] -= max_node_cut;
190 cur_slack -= max_node_cut;
191 verts[max_node].comp = 0;
192 coef_list[i][max_node >> DELETE_POWER] ^=
193 (1 << (max_node & DELETE_AND));
194 if (cur_slack < -etol){/*if the cut is now violated, impose
195 it*/
196 type = (compnodes[i+1] < vertnum/2 ?
197 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
198 rhs = (type == SUBTOUR_ELIM_SIDE ?
199 RHS(compnodes[i+1], compdemands[i+1],
200 capacity): 2*cur_bins);
201 num_cuts += addVrpCut(conPool, coef_list[i], rhs, type);
202 break;
203 }
204 }
205 }
206 }
207 delete [] coef_list[0];
208 delete [] coef_list;
209 }
210 which_connected_routine++;
211 }while(!num_cuts && which_connected_routine == BOTH &&
212 which_connected_routine < 2);
213
214 compnodes = n->compNodes_;
215
216 if (!do_greedy){
217 delete sol;
218 return num_cuts ? true : false;
219 }
220
221 if (num_cuts < 10 && do_greedy){
222 int numroutes = model_->numroutes_;
223 char *coef = new char[cut_size];
224 for (cur_edge=verts[0].first; cur_edge; cur_edge=cur_edge->next_edge){
225 for (cur_edge1 = cur_edge->other->first; cur_edge1;
226 cur_edge1 = cur_edge1->next_edge){
227 if (cur_edge1->data->weight + cur_edge->data->weight < 1 - etol)
228 continue;
229 node1 = cur_edge->other_end;
230 node2 = cur_edge1->other_end;
231 for (cur_edge2 = verts[node2].first; cur_edge2;
232 cur_edge2 = cur_edge2->next_edge){
233 if (!(cur_edge2->other_end) && node2){
234 if ((BINS(total_demand - demand[node1] - demand[node2],
235 capacity) > numroutes -1) &&
236 (cur_edge1->data->weight + cur_edge->data->weight +
237 cur_edge2->data->weight>2+etol)){
238 type = SUBTOUR_ELIM_ACROSS;
239 rhs =2*BINS(total_demand - demand[node1] -
240 demand[node2],capacity);
241 memset(coef, 0, cut_size);
242 for (i = 1; i <vertnum ; i++)
243 if ((i != node1) && (i != node2))
244 (coef[i >> DELETE_POWER]) |= (1 << (i&DELETE_AND));
245 num_cuts += addVrpCut(conPool, coef, rhs, type);
246 }
247 break;
248 }
249 }
250 }
251 }
252 delete [] coef;
253 }
254
255 //n->compNodes_ = compnodes_copy;
256
257 if (num_cuts < 10 && do_greedy){
258
259 memcpy((char *)n->newDemand_, (char *)demand, vertnum*sizeof(int));
260
261 n->reduce_graph(model_->etol_);
262 if (n->numComps_ > 1){
263 num_cuts += greedyShrinking1(
264 model_, par->entry(VrpParams::maxNumCutsInShrink), conPool);
265 }else{
266 num_cuts += greedyShrinking1One(
267 model_, par->entry(VrpParams::maxNumCutsInShrink), conPool);
268 }
269 }
270
271 if (num_cuts < 10 && do_greedy){
272 if (par->entry(VrpParams::doExtraInRoot)){
273 num_trials =
274 //level ? par->entry(VrpParams::greedyNumTrials):
275 2*par->entry(VrpParams::greedyNumTrials);
276 }else{
277 num_trials = par->entry(VrpParams::greedyNumTrials);
278 }
279 if (n->numComps_){
280 num_cuts += greedyShrinking6(
281 model_,par->entry(VrpParams::maxNumCutsInShrink),
282 num_cuts ? num_trials : 2 * num_trials, 10.5, conPool);
283 }else{
284 num_cuts += greedyShrinking6One(
285 model_, par->entry(VrpParams::maxNumCutsInShrink),
286 num_cuts ? num_trials : 2 * num_trials, 10.5, conPool);
287 }
288 }
289
290 n->compNodes_ = compnodes;
291
292 delete[] compnodes_copy;
293 delete sol;
294
295 return num_cuts ? true : false;
296 }
297
298 /*===========================================================================*/
299
300 int
connectivityCuts(BcpsConstraintPool & conPool)301 VrpCutGenerator::connectivityCuts(BcpsConstraintPool &conPool)
302 {
303 int vertnum = model_->vertnum_;
304 elist *cur_route_start;
305 edge *edge_data;
306 int weight = 0, reduced_weight, *route;
307 int cur_vert = 0, prev_vert, cust_num = 0, cur_route, rcnt;
308 int i, reduced_cust_num, vert1, vert2, type, rhs;
309 int cut_size = (vertnum >> DELETE_POWER) +1;
310 int capacity = model_->capacity_;
311 VrpNetwork *n = model_->n_;
312 vertex *verts = n->verts_;
313 double *compcuts = n->compCuts_;
314 int *compnodes = n->compNodes_;
315 int *compdemands = n->compDemands_;
316 double etol = model_->etol_;
317 int num_cuts = 0;
318
319 if (!n->isIntegral_) return 0;
320
321 /* This is a flag to tell the cut generator that the network has not been
322 constructed until the next call to userFeasibleSolution();*/
323 n->isIntegral_ = false;
324
325 /*get the components of the solution graph without the depot to check if the
326 graph is connected or not*/
327 /* rcnt = n->connected(); */ /* This was previously executed */
328 rcnt = n->numComps_;
329 char **coef_list = new char *[rcnt];
330 memset(coef_list, 0, rcnt*sizeof(char *));
331 coef_list[0] = new char[rcnt * cut_size];
332 memset(coef_list[0], 0, rcnt*cut_size*sizeof(char));
333 for(i = 1; i < rcnt; i++)
334 coef_list[i] = coef_list[0]+i*cut_size;
335 for(i = 1; i < rcnt; i++)
336 coef_list[i] = coef_list[0]+i*cut_size;
337
338 for(i = 1; i < vertnum; i++)
339 (coef_list[(verts[i].comp)-1][i >> DELETE_POWER]) |=
340 (1 << (i & DELETE_AND));
341
342 /*------------------------------------------------------------------------*\
343 | For each component check to see if the cut it induces is nonzero -- each |
344 | component's cut value must be either 0 or 2 since we have integrality |
345 \*------------------------------------------------------------------------*/
346
347 for (i = 0; i < rcnt; i++){
348 if (compcuts[i+1] < etol){/*if the cut value is zero, the graph is
349 disconnected and we have a violated cut*/
350 type = (compnodes[i+1] < vertnum/2 ?
351 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
352 rhs = (type == SUBTOUR_ELIM_SIDE ?
353 RHS(compnodes[i+1], compdemands[i+1], capacity) :
354 2*BINS(compdemands[i+1], capacity));
355 num_cuts += addVrpCut(conPool, coef_list[i], rhs, type);
356 }
357 }
358
359 delete [] coef_list[0];
360 delete [] coef_list;
361
362 /*-------------------------------------------------------------------------*\
363 | if the graph is connected, check each route to see if it obeys the |
364 | capacity constraints |
365 \*-------------------------------------------------------------------------*/
366
367 int numroutes = model_->numroutes_;
368 route = new int[vertnum];
369 char *coef = new char[cut_size];
370 for (cur_route_start = verts[0].first, cur_route = 0,
371 edge_data = cur_route_start->data; cur_route < numroutes;
372 cur_route++){
373 edge_data = cur_route_start->data;
374 edge_data->scanned = true;
375 cur_vert = edge_data->v1;
376 prev_vert = weight = cust_num = 0;
377
378 memset(coef, 0, cut_size*sizeof(char));
379
380 route[0] = cur_vert;
381 while (cur_vert){
382 /*keep tracing around the route and whenever the addition
383 of the next customer causes a violation, impose the
384 constraint induced
385 by the set of customers seen so far on the route*/
386 coef[cur_vert >> DELETE_POWER]|=(1 << (cur_vert & DELETE_AND));
387 cust_num++;
388 if ((weight += verts[cur_vert].demand) > capacity){
389 type = (cust_num < vertnum/2 ?
390 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
391 rhs = (type ==SUBTOUR_ELIM_SIDE ? RHS(cust_num, weight, capacity):
392 2*BINS(weight, capacity));
393 num_cuts += addVrpCut(conPool, coef, rhs, type);
394 vert1 = route[0];
395 reduced_weight = weight;
396 reduced_cust_num = cust_num;
397 while (true){
398 if ((reduced_weight -= verts[vert1].demand) > capacity){
399 reduced_cust_num--;
400 coef[vert1 >> DELETE_POWER] &=
401 ~(1 << (vert1 & DELETE_AND));
402 type = (reduced_cust_num < vertnum/2 ?
403 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
404 rhs = (type ==SUBTOUR_ELIM_SIDE ?
405 RHS(reduced_cust_num, reduced_weight, capacity):
406 2*BINS(reduced_weight, capacity));
407 num_cuts += addVrpCut(conPool, coef, rhs, type);
408 vert1 = route[vert1];
409 }else{
410 break;
411 }
412 }
413 vert2 = route[0];
414 while (vert2 != vert1){
415 coef[vert2 >> DELETE_POWER] |= (1 << (vert2 & DELETE_AND));
416 vert2 = route[vert2];
417 }
418 }
419 if (verts[cur_vert].first->other_end != prev_vert){
420 prev_vert = cur_vert;
421 edge_data = verts[cur_vert].first->data;
422 cur_vert = verts[cur_vert].first->other_end;
423 }
424 else{
425 prev_vert = cur_vert;
426 edge_data = verts[cur_vert].last->data; /*This statement could
427 possibly be taken out to
428 speed things up a bit*/
429 cur_vert = verts[cur_vert].last->other_end;
430 }
431 route[prev_vert] = cur_vert;
432 }
433 edge_data->scanned = true;
434
435 while (cur_route_start->data->scanned){/*find the next edge leading out of
436 the depot which has not yet been
437 traversed to start the next
438 route*/
439 if (!(cur_route_start = cur_route_start->next_edge)) break;
440 }
441 }
442
443 delete [] route;
444 delete [] coef;
445
446 for (cur_route_start = verts[0].first; cur_route_start;
447 cur_route_start = cur_route_start->next_edge)
448 cur_route_start->data->scanned = false;
449
450 // return if need resolve LP immediately.
451 return num_cuts;
452 }
453
454 /*===========================================================================*/
455
456 int
greedyShrinking1(VrpModel * m,int max_shrink_cuts,BcpsConstraintPool & conPool)457 VrpCutGenerator::greedyShrinking1(VrpModel *m,
458 int max_shrink_cuts,
459 BcpsConstraintPool &conPool)
460 {
461 VrpNetwork *n = m->n_;
462 double set_cut_val, set_demand;
463 vertex *verts = n->verts_;
464 elist *e;
465 int shrink_cuts = 0, i, j, k;
466 char *pt, *cutpt;
467 int *ipt;
468 double *dpt;
469 int vertnum = n->vertnum_;
470 int truck_cap = m->capacity_, type;
471
472 int max_vert = 0, set_size, begin = 1, cur_comp, end = 1, other_end;
473 double maxval, weight;
474 vertex *cur_nodept;
475 int *compmembers = n->compMembers_;
476 int *compnodes = n->compNodes_;
477 int *demand = n->newDemand_;
478 double etol = m->etol_;
479
480 int rhs;
481 int size = (vertnum >> DELETE_POWER) + 1;
482 char *coef = new char[size];
483 memset(coef, 0, size);
484 memset(cutList_, 0, size * (max_shrink_cuts + 1));
485
486 *inSet_ = 0;
487
488 for (i = 1; i < vertnum; i++){
489 if (verts[compmembers[i]].deleted) compmembers[i] = 0;
490 ref_[compmembers[i]] = i;
491 }
492 *ref_ = 0;
493 /* ref_ is a reference array for compmembers: gives a place
494 in which a vertex is listed in compmembers */
495
496 for (cur_comp = 1; cur_comp <= n->numComps_;
497 begin += compnodes[cur_comp], cur_comp++){ /* for every component */
498 for (i = begin, end = begin + compnodes[cur_comp]; i < end; i++){
499 if (compmembers[i] == 0) continue;
500 /* for every node as a starting one */
501 /*initialize the data structures */
502 memset(inSet_ + begin, 0, compnodes[cur_comp] * sizeof(char));
503 memset(cutVal_ + begin, 0, compnodes[cur_comp] * sizeof(double));
504 inSet_[i] = 1;
505 set_size = 1 + verts[compmembers[i]].orig_node_list_size;
506 set_cut_val = 0;
507 for (e = verts[compmembers[i]].first; e; e = e->next_edge){
508 if (e->other_end)
509 cutVal_[ref_[e->other_end]] = e->data->weight;
510 set_cut_val += e->data->weight;
511 }
512 set_demand = demand[compmembers[i]];
513
514 while(true){
515 if (set_cut_val < 2*(ceil(set_demand/truck_cap)) - etol &&
516 set_size > 2){
517 memset(coef, 0, size*sizeof(char));
518 for (j = begin, ipt = compmembers + begin; j < end; j++, ipt++){
519 if (inSet_[j]){
520 cur_nodept = verts + (*ipt);
521 if (cur_nodept->orig_node_list_size)
522 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
523 (coef[(cur_nodept->orig_node_list)[k] >>
524 DELETE_POWER]) |=
525 (1 << ((cur_nodept->orig_node_list)[k] &
526 DELETE_AND));
527 (coef[(*ipt) >> DELETE_POWER]) |=
528 (1 << ((*ipt) & DELETE_AND));
529 }
530 }
531 type = (set_size < vertnum/2 ?
532 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
533 rhs = (type == SUBTOUR_ELIM_SIDE ?
534 RHS((int)set_size,(int)set_demand,
535 (int)truck_cap):
536 2*BINS((int)set_demand, (int)truck_cap));
537 for (k = 0, cutpt = cutList_; k < shrink_cuts; k++,
538 cutpt += size)
539 if (!memcmp(coef, cutpt, size*sizeof(char)))
540 break;/* same cuts */
541 if (k >= shrink_cuts){
542 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
543 memcpy(cutpt, coef, size);
544 }
545 if (shrink_cuts > max_shrink_cuts){
546 delete [] coef;
547 return(shrink_cuts);
548 }
549 }
550 for (maxval = -1, pt = inSet_ + begin, dpt = cutVal_ + begin,
551 j = begin; j < end; pt++, dpt++, j++){
552 if (!(*pt) && *dpt > maxval){
553 maxval = cutVal_[j];
554 max_vert = j;
555 }
556 }
557 if (maxval > 0){ /* add the vertex to the set */
558 inSet_[max_vert]=1;
559 set_size += 1+ verts[compmembers[max_vert]].orig_node_list_size;
560 set_demand += demand[compmembers[max_vert]];
561 cutVal_[max_vert] = 0;
562 for (e=verts[compmembers[max_vert]].first; e; e = e->next_edge){
563 other_end = ref_[e->other_end];
564 weight = e->data->weight;
565 set_cut_val += (inSet_[other_end]) ? (-weight):weight;
566 cutVal_[other_end] += (inSet_[other_end]) ? 0 : weight;
567
568 }
569 }
570 else{ /* can't add anything to the set */
571 break;
572 }
573 }
574 }
575 }
576
577 delete [] coef;
578
579 return(shrink_cuts);
580 }
581
582 /*===========================================================================*/
583
584 int
greedyShrinking6(VrpModel * m,int max_shrink_cuts,int trial_num,double prob,BcpsConstraintPool & conPool)585 VrpCutGenerator::greedyShrinking6(VrpModel *m,
586 int max_shrink_cuts,
587 int trial_num,
588 double prob,
589 BcpsConstraintPool &conPool)
590 {
591 VrpNetwork *n = m->n_;
592 double set_cut_val, set_demand;
593 vertex *verts = n->verts_;
594 elist *e;
595 int i, j, k, shrink_cuts = 0;
596 char *pt, *cutpt;
597 double *dpt;
598 int vertnum = n->vertnum_, type;
599 int truck_cap = m->capacity_;
600
601 int max_vert = 0, set_size, begin = 1, cur_comp, end = 1, num_trials;
602 double maxval;
603 double denominator=pow(2.0,31.0)-1.0;
604 double r, q;
605
606 int other_end;
607 double weight;
608 int *ipt;
609 vertex *cur_nodept;
610 int *compmembers = n->compMembers_;
611 int *compnodes = n->compNodes_;
612 int *demand = n->newDemand_;
613 double etol = m->etol_;
614
615 int rhs;
616 int size = (vertnum >> DELETE_POWER) + 1;
617 char *coef = new char[size];
618 memset(coef, 0, size);
619 memset(cutList_, 0, size * (max_shrink_cuts +1));
620
621
622 *inSet_=0;
623
624 for(i = 1; i < vertnum; i++){
625 if (verts[compmembers[i]].deleted) compmembers[i] = 0;
626 ref_[compmembers[i]] = i;
627 }
628 *ref_ = 0;
629
630 /* ref_ is a reference array for compmembers: gives a place
631 in which a vertex is listed in compmembers */
632
633 for (cur_comp = 1; cur_comp <= n->numComps_; begin += compnodes[cur_comp],
634 cur_comp++){
635 /* for every component */
636 if (compnodes[cur_comp] <= 7) continue;
637
638 for (num_trials = 0; num_trials < trial_num * compnodes[cur_comp];
639 num_trials++){
640 end = begin + compnodes[cur_comp];
641 /*initialize the data structures */
642 memset(inSet_ + begin, 0, compnodes[cur_comp] * sizeof(char));
643 memset(cutVal_+ begin, 0, compnodes[cur_comp] * sizeof(double));
644
645 set_cut_val = 0;
646 set_size = 0;
647 set_demand = 0;
648 for (i = begin; i < end; i++ ){
649 if (compmembers[i] == 0) continue;
650 r = (RANDOM()/denominator);
651 q = (prob/compnodes[cur_comp]);
652 if (r < q){
653 inSet_[i] = 1;
654 set_size += 1 + verts[compmembers[i]].orig_node_list_size;
655 set_demand += demand[compmembers[i]];
656 for (e = verts[compmembers[i]].first; e; e = e-> next_edge){
657 other_end = ref_[e->other_end];
658 weight = e->data->weight;
659 set_cut_val += (inSet_[other_end]) ? (-weight) : weight;
660 cutVal_[other_end] += (inSet_[other_end]) ? 0 : weight;
661 }
662 }
663 }
664 while(set_size){
665 if (set_cut_val < 2*(ceil(set_demand/truck_cap)) - etol &&
666 set_size > 2){
667 memset(coef, 0, size*sizeof(char));
668 for (j = begin, ipt = compmembers + begin; j < end; j++, ipt++)
669 if (inSet_[j]){
670 cur_nodept = verts + (*ipt);
671 if (cur_nodept->orig_node_list_size)
672 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
673 (coef[(cur_nodept->orig_node_list)[k] >>
674 DELETE_POWER]) |=
675 (1 << ((cur_nodept->orig_node_list)[k] &
676 DELETE_AND));
677 (coef[(*ipt) >> DELETE_POWER]) |= (1 << ((*ipt) &
678 DELETE_AND));
679 }
680 type = (set_size < vertnum/2 ?
681 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
682 rhs = (type == SUBTOUR_ELIM_SIDE ?
683 RHS((int)set_size,(int)set_demand, (int)truck_cap):
684 2*BINS((int)set_demand, (int)truck_cap));
685 for (k = 0, cutpt = cutList_; k < shrink_cuts; k++,
686 cutpt += size)
687 if (!memcmp(coef, cutpt, size*sizeof(char))) break;
688 if ( k >= shrink_cuts){
689 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
690 memcpy(cutpt, coef, size);
691 }
692
693 if ( shrink_cuts > max_shrink_cuts){
694 delete [] coef;
695 return(shrink_cuts);
696 }
697 }
698 for (maxval = -1, pt = inSet_+begin, dpt = cutVal_+begin,
699 j = begin; j < end; pt++, dpt++, j++){
700 if (!(*pt) && *dpt > maxval){
701 maxval = cutVal_[j];
702 max_vert = j;
703 }
704 }
705 if (maxval > 0){ /* add the vertex to the set */
706 inSet_[max_vert]=1;
707 set_size+=1+ verts[compmembers[max_vert]].orig_node_list_size;
708 set_demand+=demand[compmembers[max_vert]];
709 cutVal_[max_vert]=0;
710 for (e = verts[compmembers[max_vert]].first; e;
711 e = e->next_edge){
712 other_end = ref_[e->other_end];
713 weight = e->data->weight;
714 set_cut_val += (inSet_[other_end]) ? (-weight) : weight;
715 cutVal_[other_end]+=(inSet_[other_end]) ? 0 : weight;
716 }
717 }
718 else{ /* can't add anything to the set */
719 break;
720 }
721 }
722 }
723 }
724
725 delete [] coef;
726 return shrink_cuts ? true : false;
727 }
728
729 /*===========================================================================*/
730
731 int
greedyShrinking1One(VrpModel * m,int max_shrink_cuts,BcpsConstraintPool & conPool)732 VrpCutGenerator::greedyShrinking1One(VrpModel *m,
733 int max_shrink_cuts,
734 BcpsConstraintPool &conPool)
735 {
736 VrpNetwork *n = m->n_;
737 double set_cut_val, set_demand;
738 vertex *verts = n->verts_;
739 elist *e;
740 int i, j, k, shrink_cuts = 0;
741 char *pt, *cutpt;
742 double *dpt;
743 int vertnum = n->vertnum_, type;
744 int truck_cap = m->capacity_;
745 int max_vert = 0;
746 int set_size;
747 /* int flag=0; */
748
749 double complement_demand, total_demand = verts[0].demand;
750 double complement_cut_val;
751 int complement_size;
752 double maxval;
753 int other_end;
754 double weight;
755 vertex *cur_nodept;
756 int *demand = n->newDemand_;
757 double etol = m->etol_;
758
759 int rhs;
760 int size = (vertnum >> DELETE_POWER) + 1;
761 char *coef = new char[size];
762 memset(coef, 0, size);
763 memset(cutList_, 0, size * (max_shrink_cuts + 1));
764
765 for (i = 1; i < vertnum; i++ ){
766 if (verts[i].deleted) continue;/* for every node as a starting one */
767 /*initialize the data structures */
768 memset(inSet_, 0, vertnum*sizeof(char));
769 memset(cutVal_, 0,vertnum* sizeof(double));
770 inSet_[i] = 1;
771 set_size = 1 + verts[i].orig_node_list_size;
772 set_cut_val = 0;
773 for (e= verts[i].first; e; e = e-> next_edge){
774 weight = e->data->weight;
775 cutVal_[e->other_end] = weight;
776 set_cut_val += weight;
777 }
778 set_demand = demand[i];
779
780 while(true){
781 if (set_cut_val < 2*(ceil(set_demand/truck_cap)) - etol &&
782 set_size > 2){
783 memset(coef, 0, size*sizeof(char));
784 /* printf("%d :", i); */
785 /* printf("%d ", j); */
786 for (j = 1; j < vertnum; j++)
787 if (inSet_[j]){
788 cur_nodept = verts + j;
789 if (cur_nodept->orig_node_list_size)
790 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
791 (coef[(cur_nodept->orig_node_list)[k] >>
792 DELETE_POWER]) |=
793 (1 << ((cur_nodept->orig_node_list)[k] &
794 DELETE_AND));
795 (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
796 }
797 /* printf("%f ", set_demand);
798 printf("%f \n",set_cut_val);*/
799 type = (set_size < vertnum/2 ?
800 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
801 rhs = (type == SUBTOUR_ELIM_SIDE ?
802 RHS((int)set_size,(int)set_demand,(int)truck_cap):
803 2*BINS((int)set_demand,(int)truck_cap));
804 for (k = 0, cutpt = cutList_; k < shrink_cuts; k++,
805 cutpt += size)
806 if (!memcmp(coef, cutpt, size*sizeof(char)))
807 break; /* same cuts */
808 if ( k >= shrink_cuts){
809 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
810 memcpy(cutpt, coef, size);
811 }
812
813 if ( shrink_cuts > max_shrink_cuts){
814 delete [] coef;
815 return(shrink_cuts);
816 }
817 }
818 /* check the complement */
819
820 complement_demand = total_demand - set_demand;
821 complement_cut_val =set_cut_val- 2*(*cutVal_) + 2*m->numroutes_;
822 complement_size = vertnum - 1 - set_size;
823 if (complement_cut_val< 2*(ceil(complement_demand/truck_cap))-etol &&
824 complement_size > 2){
825 memset(coef, 0, size*sizeof(char));
826 for (j = 1; j < vertnum; j++)
827 if (!(inSet_[j]) && !(verts[j].deleted)){
828 cur_nodept = verts + j;
829 if (cur_nodept->orig_node_list_size)
830 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
831 (coef[(cur_nodept->orig_node_list)[k] >>
832 DELETE_POWER]) |=
833 (1 << ((cur_nodept->orig_node_list)[k] &
834 DELETE_AND));
835 (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
836 }
837 type = (complement_size < vertnum/2 ?
838 SUBTOUR_ELIM_SIDE : SUBTOUR_ELIM_ACROSS);
839 rhs = (type == SUBTOUR_ELIM_SIDE ?
840 RHS((int)complement_size,(int)complement_demand,
841 (int)truck_cap):
842 2*BINS((int)complement_demand,(int)truck_cap));
843 for (k=0, cutpt = cutList_; k < shrink_cuts; k++,
844 cutpt += size)
845 if (!memcmp(coef, cutpt, size*sizeof(char))) break;
846 if ( k >= shrink_cuts){
847 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
848 memcpy(cutpt, coef, size);
849 }
850
851 if (shrink_cuts > max_shrink_cuts){
852 delete [] coef;
853 return(shrink_cuts);
854 }
855 }
856
857 for (maxval = -1, pt = inSet_, dpt = cutVal_,pt++, dpt++,
858 j = 1; j < vertnum; pt++, dpt++, j++){
859 if (!(*pt) && *dpt > maxval){
860 maxval = cutVal_[j];
861 max_vert = j;
862 }
863 }
864 if (maxval > 0){ /* add the vertex to the set */
865 inSet_[max_vert] = 1;
866 set_size += 1 + verts[max_vert].orig_node_list_size ;
867 set_demand += demand[max_vert];
868 cutVal_[max_vert] = 0;
869 for (e = verts[max_vert].first; e; e = e-> next_edge){
870 other_end = e->other_end;
871 weight = e->data->weight;
872 set_cut_val += (inSet_[other_end]) ? (-weight): weight;
873 cutVal_[other_end] += weight;
874
875 }
876 }
877 else{ /* can't add anything to the set */
878 break;
879 }
880 }
881 }
882
883 delete [] coef;
884 return(shrink_cuts);
885 }
886
887 /*===========================================================================*/
888
889 int
greedyShrinking6One(VrpModel * m,int max_shrink_cuts,int trial_num,double prob,BcpsConstraintPool & conPool)890 VrpCutGenerator::greedyShrinking6One(VrpModel *m,
891 int max_shrink_cuts,
892 int trial_num,
893 double prob,
894 BcpsConstraintPool &conPool)
895 {
896 VrpNetwork *n = m->n_;
897 double set_cut_val, set_demand;
898 vertex *verts=n->verts_;
899 elist *e;
900 int i, j, k, shrink_cuts = 0;
901 char *pt, *cutpt;
902 double *dpt;
903 int vertnum = n->vertnum_, type;
904 int truck_cap = m->capacity_;
905
906 int max_vert = 0, set_size, begin = 1, end = 1, num_trials;
907 double maxval, r, q;
908 double denominator = pow(2.0, 31.0) - 1.0;
909
910 int other_end;
911 double weight;
912
913 double complement_demand, total_demand = verts[0].demand;
914 double complement_cut_val;
915 int complement_size;
916 vertex *cur_nodept;
917 /* int flag=0;*/
918 int *demand = n->newDemand_;
919 double etol = m->etol_;
920
921 int rhs;
922 int size = (vertnum >> DELETE_POWER) + 1;
923 char *coef = new char[size];
924 memset(coef, 0, size);
925 memset(cutList_, 0, size * (max_shrink_cuts +1));
926
927 *inSet_ = 0;
928
929 for (num_trials = 0; num_trials < trial_num*vertnum ; num_trials++){
930
931 /*initialize the data structures */
932 memset(inSet_, 0, vertnum*sizeof(char));
933 memset(cutVal_, 0,vertnum* sizeof(double));
934
935 set_cut_val = 0;
936 set_size = 0;
937 set_demand = 0;
938 for (i = 1 ; i < vertnum; i++ ){
939 if (verts[i].deleted) continue;
940 r = (RANDOM()/denominator);
941 q = (prob/vertnum);
942 if (r < q){
943 inSet_[i] = 1;
944 set_size += 1 + verts[i].orig_node_list_size;
945 set_demand += demand[i];
946 for (e = verts[i].first; e; e = e-> next_edge){
947 other_end = e->other_end;
948 weight = e->data->weight;
949 set_cut_val += (inSet_[other_end]) ? (-weight) : weight;
950 cutVal_[other_end] += (inSet_[other_end]) ? 0 : weight;
951 }
952 }
953 }
954 while(set_size){
955 if (set_cut_val < 2*(ceil(set_demand/truck_cap)) - etol &&
956 set_size > 2){
957 memset(coef, 0, size*sizeof(char));
958 for (j = 1; j < vertnum; j++ )
959 if (inSet_[j]){
960 cur_nodept = verts + j;
961 if (cur_nodept->orig_node_list_size)
962 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
963 (coef[(cur_nodept->orig_node_list)[k] >>
964 DELETE_POWER]) |=
965 (1 << ((cur_nodept->orig_node_list)[k] &
966 DELETE_AND));
967 (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
968 }
969 type = (set_size < vertnum/2 ?
970 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
971 rhs = (type == SUBTOUR_ELIM_SIDE ?
972 RHS((int)set_size, (int)set_demand, (int)truck_cap):
973 2*BINS((int)set_demand, (int)truck_cap));
974 for (k = 0, cutpt = cutList_; k < shrink_cuts; k++,
975 cutpt += size)
976 if (!memcmp(coef, cutpt, size*sizeof(char))) break;
977 if ( k >= shrink_cuts){
978 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
979 memcpy(cutpt, coef, size);
980 }
981 if (shrink_cuts > max_shrink_cuts){
982 delete [] coef;
983 return(shrink_cuts);
984 }
985 }
986
987 /* check the complement */
988
989 complement_demand = total_demand - set_demand;
990 complement_cut_val = set_cut_val - 2*(*cutVal_) + 2*m->numroutes_;
991 complement_size = vertnum - 1 - set_size;
992 if (complement_cut_val< 2*(ceil(complement_demand/truck_cap))-etol &&
993 complement_size > 2){
994 memset(coef, 0, size*sizeof(char));
995 for (j = 1; j < vertnum; j++)
996 if (!(inSet_[j])&& !(verts[j].deleted)){
997 cur_nodept = verts + j;
998 if (cur_nodept->orig_node_list_size)
999 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1000 (coef[(cur_nodept->orig_node_list)[k] >>
1001 DELETE_POWER]) |=
1002 (1 << ((cur_nodept->orig_node_list)[k] &
1003 DELETE_AND));
1004 (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
1005 }
1006 type = (complement_size < vertnum/2 ?
1007 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
1008 rhs = (type == SUBTOUR_ELIM_SIDE ?
1009 RHS((int)complement_size,(int)complement_demand,
1010 (int)truck_cap):
1011 2*BINS((int)complement_demand,(int)truck_cap));
1012 for (k = 0, cutpt = cutList_; k < shrink_cuts; k++,
1013 cutpt += size)
1014 if (!memcmp(coef, cutpt, size*sizeof(char))) break;
1015 if ( k >= shrink_cuts){
1016 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
1017 memcpy(cutpt, coef, size);
1018 }
1019
1020 if (shrink_cuts > max_shrink_cuts){
1021 delete [] coef;
1022 return(shrink_cuts);
1023 }
1024 }
1025
1026 for (maxval = -1, pt = inSet_ + begin, dpt = cutVal_ + begin,
1027 j = begin; j < end; pt++, dpt++, j++){
1028 if (!(*pt) && *dpt > maxval){
1029 maxval = cutVal_[j];
1030 max_vert = j;
1031 }
1032 }
1033 if (maxval > 0){ /* add the vertex to the set */
1034 inSet_[max_vert] = 1;
1035 set_size += 1 + verts[max_vert].orig_node_list_size ;
1036 set_demand += demand[max_vert];
1037 cutVal_[max_vert] = 0;
1038 for (e = verts[max_vert].first; e; e = e-> next_edge){
1039 other_end = e->other_end;
1040 weight = e->data->weight;
1041 set_cut_val += (inSet_[other_end]) ? (-weight) : weight;
1042 cutVal_[other_end] += weight;
1043 }
1044 }
1045 else{ /* can't add anything to the set */
1046 break;
1047 }
1048 }
1049 }
1050
1051 delete [] coef;
1052
1053 return(shrink_cuts);
1054 }
1055
1056 /*===========================================================================*/
1057
1058 int
greedyShrinking2One(VrpModel * m,int max_shrink_cuts,BcpsConstraintPool & conPool)1059 VrpCutGenerator::greedyShrinking2One(VrpModel *m,
1060 int max_shrink_cuts,
1061 BcpsConstraintPool &conPool)
1062 {
1063 VrpNetwork *n = m->n_;
1064 double set_cut_val, set_demand;
1065 vertex *verts = n->verts_;
1066 elist *e, *cur_edge1, *cur_edge2;
1067 int j, k, shrink_cuts = 0;
1068 char *pt;
1069 double *dpt;
1070 int vertnum = n->vertnum_, type;
1071 int truck_cap = m->capacity_;
1072
1073 int max_vert = 0, set_size, begin = 1, end = 1;
1074 double maxval;
1075
1076 int other_end;
1077 double weight;
1078
1079 double complement_demand, total_demand = verts[0].demand;
1080 double complement_cut_val;
1081 int complement_size;
1082 vertex *cur_nodept;
1083 int *demand = n->newDemand_;
1084 double etol = m->etol_;
1085
1086 int rhs;
1087 int size = (vertnum >> DELETE_POWER) + 1;
1088 char *coef = new char[size];
1089 memset(coef, 0, size);
1090
1091 *inSet_=0;
1092
1093 for (cur_edge1 = verts[0].first; cur_edge1;
1094 cur_edge1 = cur_edge1->next_edge){
1095 for (cur_edge2 = cur_edge1->next_edge; cur_edge2;
1096 cur_edge2 = cur_edge2->next_edge){
1097
1098 /*initialize the data structures */
1099 memset(inSet_, 0, vertnum*sizeof(char));
1100 memset(cutVal_, 0,vertnum* sizeof(double));
1101
1102 set_cut_val = 2;
1103 set_size = 2 + cur_edge1->other->orig_node_list_size +
1104 cur_edge2->other->orig_node_list_size;
1105 set_demand = demand[cur_edge1->other_end] +
1106 demand[cur_edge2->other_end];
1107 inSet_[cur_edge1->other_end] = 1;
1108
1109 for (e = verts[cur_edge1->other_end].first; e; e = e-> next_edge){
1110 cutVal_[e->other_end] += e->data->weight;
1111 }
1112
1113 inSet_[cur_edge2->other_end] = 1;
1114 for (e = verts[cur_edge2->other_end].first; e; e = e-> next_edge){
1115 other_end = e->other_end;
1116 weight = e->data->weight;
1117 set_cut_val += (inSet_[other_end]) ? (-weight) : weight;
1118 cutVal_[other_end] += (inSet_[other_end]) ? 0 : weight;
1119 }
1120 while(set_size){
1121 if (set_cut_val < 2*(ceil(set_demand/truck_cap)) - etol &&
1122 set_size > 2){
1123 memset(coef, 0, size*sizeof(char));
1124 for (j = 1; j < vertnum; j++ )
1125 if (inSet_[j]){
1126 cur_nodept = verts + j;
1127 if (cur_nodept->orig_node_list_size)
1128 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1129 (coef[(cur_nodept->orig_node_list)[k] >>
1130 DELETE_POWER]) |=
1131 (1 << ((cur_nodept->orig_node_list)[k] &
1132 DELETE_AND));
1133 (coef[j >> DELETE_POWER]) |= (1 << (j & DELETE_AND));
1134 }
1135 type = (set_size < vertnum/2 ?
1136 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
1137 rhs = (type == SUBTOUR_ELIM_SIDE ?
1138 RHS((int)set_size, (int)set_demand, (int)truck_cap):
1139 2*BINS((int)set_demand, (int)truck_cap));
1140 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
1141 }
1142
1143 /* check the complement */
1144
1145 complement_demand = total_demand - set_demand;
1146 complement_cut_val = set_cut_val - 2*(*cutVal_) + 2*m->numroutes_;
1147 complement_size = vertnum - 1 - set_size;
1148 if (complement_cut_val<2*(ceil(complement_demand/truck_cap))-etol&&
1149 complement_size > 2){
1150 memset(coef, 0, size*sizeof(char));
1151 for (j = 1; j < vertnum; j++)
1152 if (!inSet_[j]){
1153 cur_nodept=verts + j;
1154 if (cur_nodept->orig_node_list_size)
1155 for (k = 0; k < cur_nodept->orig_node_list_size; k++)
1156 (coef[(cur_nodept->orig_node_list)[k] >>
1157 DELETE_POWER]) |=
1158 (1 << ((cur_nodept->orig_node_list)[k] &
1159 DELETE_AND));
1160 (coef[j>> DELETE_POWER]) |= (1 << ( j & DELETE_AND));
1161 }
1162 type = (complement_size < vertnum/2 ?
1163 SUBTOUR_ELIM_SIDE:SUBTOUR_ELIM_ACROSS);
1164 rhs = (type == SUBTOUR_ELIM_SIDE ?
1165 RHS((int)complement_size,(int)complement_demand,
1166 (int)truck_cap):
1167 2*BINS((int)complement_demand,(int)truck_cap));
1168 shrink_cuts += addVrpCut(conPool, coef, rhs, type);
1169 }
1170
1171 for (maxval = -1, pt = inSet_+begin, dpt = cutVal_+begin,
1172 j = begin; j < end; pt++, dpt++, j++){
1173 if (!(*pt) && *dpt > maxval){
1174 maxval = cutVal_[j];
1175 max_vert = j;
1176 }
1177 }
1178 if (maxval > 0){ /* add the vertex to the set */
1179 inSet_[max_vert] = 1;
1180 set_size += 1 + verts[max_vert].orig_node_list_size ;
1181 set_demand += demand[max_vert];
1182 cutVal_[max_vert] = 0;
1183 for (e = verts[max_vert].first; e; e = e-> next_edge){
1184 other_end = e->other_end;
1185 weight = e->data->weight;
1186 set_cut_val += (inSet_[other_end]) ? (-weight) : weight;
1187 cutVal_[other_end] += weight;
1188 }
1189 }
1190 else{ /* can't add anything to the set */
1191 break;
1192 }
1193 }
1194 }
1195 }
1196
1197 delete [] coef;
1198
1199 return(shrink_cuts);
1200 }
1201
1202 /*===========================================================================*/
1203
1204 int
addVrpCut(BcpsConstraintPool & conPool,char * coef,int rhs,int type)1205 VrpCutGenerator::addVrpCut(BcpsConstraintPool &conPool,
1206 char *coef,
1207 int rhs,
1208 int type)
1209 {
1210 int i, nzcnt = 0, nzcnt_side = 0, nzcnt_across = 0;
1211 int v0, v1;
1212 std::vector<VrpVariable *>edges = model_->getEdgeList();
1213 int *matind = NULL, *matind_across, *matind_side;
1214 double *matval = NULL;
1215 int edgenum = model_->edgenum_;
1216 double infinity = model_->solver()->getInfinity();
1217 char sense = 'A';
1218
1219 switch (type){
1220 /*-------------------------------------------------------------------*\
1221 * The subtour elimination constraints are stored as a vector of
1222 * bits indicating which side of the cut each customer is on
1223 \*-------------------------------------------------------------------*/
1224
1225 case SUBTOUR_ELIM:
1226 matind_side = new int[edgenum];
1227 matind_across = new int[edgenum];
1228 for (i = 0, nzcnt = 0; i < edgenum; i++){
1229 v0 = edges[i]->getv0();
1230 v1 = edges[i]->getv1();
1231 if (coef[v0 >> DELETE_POWER] &
1232 (1 << (v0 & DELETE_AND)) &&
1233 (coef[v1 >> DELETE_POWER]) &
1234 (1 << (v1 & DELETE_AND))){
1235 matind_side[nzcnt_side++] = i;
1236 }else if ((coef[v0 >> DELETE_POWER] >>
1237 (v0 & DELETE_AND) & 1) ^
1238 (coef[v1 >> DELETE_POWER] >>
1239 (v1 & DELETE_AND) & 1)){
1240 matind_across[nzcnt_across++] = i;
1241 }
1242 }
1243 type = nzcnt_side < nzcnt_across ? SUBTOUR_ELIM_SIDE :
1244 SUBTOUR_ELIM_ACROSS;
1245 switch (type){
1246 case SUBTOUR_ELIM_SIDE:
1247 nzcnt = nzcnt_side;
1248 matind = matind_side;
1249 rhs = 0; /*RHS(compnodes[i+1],compdemands[i+1], capacity)*/
1250 sense = 'L';
1251 delete [] matind_across;
1252 break;
1253
1254 case SUBTOUR_ELIM_ACROSS:
1255 nzcnt = nzcnt_across;
1256 matind = matind_across;
1257 rhs = 0; /*2*BINS(compdemands[i+1], capacity)*/
1258 sense = 'G';
1259 delete [] matind_side;
1260 break;
1261 }
1262
1263 break;
1264
1265
1266 case SUBTOUR_ELIM_SIDE:
1267 matind = new int[edgenum];
1268 for (i = 0, nzcnt = 0; i < edgenum; i++){
1269 v0 = edges[i]->getv0();
1270 v1 = edges[i]->getv1();
1271 if (coef[v0 >> DELETE_POWER] &
1272 (1 << (v0 & DELETE_AND)) &&
1273 (coef[v1 >> DELETE_POWER]) &
1274 (1 << (v1 & DELETE_AND))){
1275 matind[nzcnt++] = i;
1276 }
1277 }
1278 sense = 'L';
1279 break;
1280
1281 case SUBTOUR_ELIM_ACROSS:
1282 matind = new int[edgenum];
1283 for (i = 0, nzcnt = 0; i < edgenum; i++){
1284 v0 = edges[i]->getv0();
1285 v1 = edges[i]->getv1();
1286 if ((coef[v0 >> DELETE_POWER] >>
1287 (v0 & DELETE_AND) & 1) ^
1288 (coef[v1 >> DELETE_POWER] >>
1289 (v1 & DELETE_AND) & 1)){
1290 matind[nzcnt++] = i;
1291 }
1292 }
1293 sense = 'G';
1294 break;
1295
1296 default:
1297 break;
1298
1299 }
1300
1301 matval = new double[nzcnt];
1302 for (i = nzcnt-1; i >= 0; i--)
1303 matval[i] = 1;
1304
1305 BlisConstraint *blisCon = NULL;
1306
1307 if (sense == 'L'){
1308 blisCon = new BlisConstraint(-infinity, rhs,
1309 -infinity, rhs,
1310 nzcnt, matind, matval);
1311 blisCon->setValidRegion(BcpsValidGlobal);
1312 }
1313 else if (sense == 'G'){
1314 blisCon = new BlisConstraint(rhs, infinity,
1315 rhs, infinity,
1316 nzcnt, matind, matval);
1317 blisCon->setValidRegion(BcpsValidGlobal);
1318 }
1319 else{
1320 return 0;
1321 }
1322
1323 conPool.addConstraint(blisCon);
1324
1325 delete [] matind;
1326 delete [] matval;
1327
1328 return 1;
1329 }
1330
1331 /*===========================================================================*/
1332
1333 #ifdef DO_TSP_CUTS
1334 int
tspCuts(VrpModel * m,BcpsConstraintPool & conPool)1335 VrpCutGenerator::tspCuts(VrpModel *m, BcpsConstraintPool &conPool)
1336 {
1337 VrpNetwork *n = m->n_;
1338 VrpParams *par = model_->VrpPar_;
1339 int edgenum = n->edgenum_;
1340 int vertnum = m->vertnum_;
1341 int verbosity = par->verbosity;
1342
1343 edge *edges = n->edges_;
1344 CCtsp_lpcut_in *tsp_cuts = NULL;
1345 int *tsp_edgelist = new int[2*edgenum];
1346 double *tsp_x = new double[edgenum];
1347 int i, cutnum = 0, cuts_added = 0, rval, seed;
1348 CCrandstate rstate;
1349 CCtsp_cutselect *sel = new CCtsp_cutselect;
1350 CCtsp_tighten_info *stats = new CCtsp_tighten_info;
1351 memset (sel, 0, sizeof(CCtsp_cutselect));
1352 memset (stats, 0, sizeof(CCtsp_tighten_info));
1353 CCtsp_lpgraph g;
1354
1355 sel->connect = 1;
1356 if (par->whichTspCuts & SUBTOUR){
1357 sel->segments = 1;
1358 sel->exactsubtour = 1;
1359 }
1360 if (par->whichTspCuts & BLOSSOM){
1361 sel->fastblossom = 1;
1362 sel->ghfastblossom = 1;
1363 sel->exactblossom = 0;
1364 }
1365 if (par->whichTspCuts & COMB){
1366 sel->blockcombs = 1;
1367 sel->growcombs = 0;
1368 sel->prclique = 0;
1369 }
1370
1371 for (i = 0; i < edgenum; i++, edges++){
1372 tsp_edgelist[i << 1] = edges->v0;
1373 tsp_edgelist[(i << 1) + 1] = edges->v1;
1374 tsp_x[i] = edges->weight;
1375 }
1376
1377 CCtsp_init_lpgraph_struct (&g);
1378 CCtsp_build_lpgraph (&g, vertnum, edgenum, tsp_edgelist, (int *) NULL);
1379 CCtsp_build_lpadj (&g, 0, edgenum);
1380
1381 if (sel->connect){
1382 rval = CCtsp_connect_cuts(&tsp_cuts, &cutnum, vertnum, edgenum,
1383 tsp_edgelist, tsp_x);
1384 if (rval) {
1385 fprintf(stderr, "CCtsp_connect_cuts failed\n");
1386 printf("CCtsp_connect_cuts failed\n");
1387 rval = 1;
1388 }
1389 if (verbosity > 3)
1390 printf("Found %2d connect cuts\n", cutnum);
1391 if (!rval && cutnum > 0){
1392 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1393 if (cuts_added){
1394 if (verbosity > 3)
1395 printf("%i connect cuts added\n", cuts_added);
1396 goto CLEANUP;
1397 }
1398 }
1399 }
1400
1401 if (sel->segments){
1402 rval = CCtsp_segment_cuts(&tsp_cuts, &cutnum, vertnum, edgenum,
1403 tsp_edgelist, tsp_x);
1404 if (rval) {
1405 fprintf(stderr, "CCtsp_segment_cuts failed\n");
1406 printf("CCtsp_segment_cuts failed\n");
1407 rval = 1;
1408 }
1409 if (verbosity > 3)
1410 printf("Found %2d segment cuts\n", cutnum);
1411 if (!rval && cutnum > 0){
1412 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1413 if (cuts_added){
1414 if (verbosity > 3)
1415 printf("%i segment cuts added\n", cuts_added);
1416 goto CLEANUP;
1417 }
1418 }
1419 }
1420
1421 if (sel->fastblossom){
1422 rval = CCtsp_fastblossom(&tsp_cuts, &cutnum, vertnum, edgenum,
1423 tsp_edgelist, tsp_x);
1424 if (rval) {
1425 fprintf(stderr, "CCtsp_fastblossom failed\n");
1426 printf("CCtsp_fastblossom failed\n");
1427 rval = 1;
1428 }
1429 if (verbosity > 3)
1430 printf("Found %2d fastblossom cuts\n", cutnum);
1431 if (!rval && cutnum > 0){
1432 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1433 if (cuts_added){
1434 if (verbosity > 3)
1435 printf("%i fastblossom cuts added\n", cuts_added);
1436 goto CLEANUP;
1437 }
1438 }
1439 }
1440
1441 if (sel->ghfastblossom){
1442 rval = CCtsp_ghfastblossom(&tsp_cuts, &cutnum, vertnum, edgenum,
1443 tsp_edgelist, tsp_x);
1444 if (rval) {
1445 fprintf(stderr, "CCtsp_ghfastblossom failed\n");
1446 printf("CCtsp_ghfastblossom failed\n");
1447 rval = 1;
1448 }
1449 if (verbosity > 3)
1450 printf("Found %2d ghfastblossom cuts\n", cutnum);
1451 if (!rval && cutnum > 0){
1452 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1453 if (cuts_added){
1454 if (verbosity > 3)
1455 printf("%i ghfastblossom cuts added\n", cuts_added);
1456 goto CLEANUP;
1457 }
1458 }
1459 }
1460
1461 if (sel->blockcombs){
1462 rval = CCtsp_block_combs(&tsp_cuts, &cutnum, vertnum, edgenum,
1463 tsp_edgelist, tsp_x, true);
1464 if (rval) {
1465 fprintf(stderr, "CCtsp_block_combs failed\n");
1466 printf("CCtsp_block_combs failed\n");
1467 rval = 1;
1468 }
1469 if (verbosity > 3)
1470 printf("Found %2d block combs\n", cutnum);
1471 if (!rval && cutnum > 0){
1472 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1473 if (cuts_added){
1474 if (verbosity > 3)
1475 printf("%i block combs added\n", cuts_added);
1476 goto CLEANUP;
1477 }
1478 }
1479 }
1480
1481 if (sel->growcombs){
1482 rval = CCtsp_edge_comb_grower(&tsp_cuts, &cutnum, vertnum,
1483 edgenum, tsp_edgelist, tsp_x, stats);
1484 if (rval) {
1485 fprintf(stderr, "CCtsp_edge_comb_grower failed\n");
1486 printf("CCtsp_edge_comb_grower failed\n");
1487 rval = 1;
1488 }
1489 if (verbosity > 3)
1490 printf("Found %2d grown combs\n", cutnum);
1491 if (!rval && cutnum > 0){
1492 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1493 if (cuts_added){
1494 if (verbosity > 3)
1495 printf("%i grown combs added\n", cuts_added);
1496 goto CLEANUP;
1497 }
1498 }
1499 }
1500
1501 if (sel->prclique){
1502 rval = CCtsp_pr_cliquetree(&tsp_cuts, &cutnum, vertnum,
1503 edgenum, tsp_edgelist, tsp_x, stats);
1504 if (rval) {
1505 fprintf(stderr, "CCtsp_pr_cliquetree failed\n");
1506 printf("CCtsp_pr_cliquetree failed\n");
1507 rval = 1;
1508 }
1509 if (verbosity > 3)
1510 printf("Found %2d PR cliquetrees\n", cutnum);
1511 if (!rval && cutnum > 0){
1512 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1513 if (cuts_added){
1514 if (verbosity > 3)
1515 printf("%i PR cliquetrees added\n", cuts_added);
1516 goto CLEANUP;
1517 }
1518 }
1519 }
1520
1521 if (sel->exactsubtour){
1522 rval = CCtsp_exact_subtours(&tsp_cuts, &cutnum, vertnum,
1523 edgenum, tsp_edgelist, tsp_x);
1524 if (rval) {
1525 fprintf(stderr, "CCtsp_exact_subtours failed\n");
1526 printf("CCtsp_exact_subtours failed\n");
1527 rval = 1;
1528 }
1529 if (verbosity > 3)
1530 printf("Found %2d exact subtours\n", cutnum);
1531 if (!rval && cutnum > 0){
1532 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1533 if (cuts_added){
1534 if (verbosity > 3)
1535 printf("%i exactsubtours added\n", cuts_added);
1536 goto CLEANUP;
1537 }
1538 }
1539 }
1540
1541 if (sel->exactblossom){
1542 seed = (int) CCutil_real_zeit ();
1543 CCutil_sprand(seed, &rstate);
1544 rval = CCtsp_exactblossom(&tsp_cuts, &cutnum, vertnum, edgenum,
1545 tsp_edgelist, tsp_x, &rstate);
1546 if (rval) {
1547 fprintf(stderr, "CCtsp_exactblossom failed\n");
1548 printf("CCtsp_exactblossom failed\n");
1549 rval = 1;
1550 }
1551 if (verbosity > 3)
1552 printf("Found %2d exactblossoms\n", cutnum);
1553 if (!rval && cutnum > 0){
1554 cuts_added += addTspCuts(m, conPool, &tsp_cuts, &g);
1555 if (cuts_added){
1556 if (verbosity > 3)
1557 printf("%i exact blossoms added\n", cuts_added);
1558 goto CLEANUP;
1559 }
1560 }
1561 }
1562
1563 CLEANUP:
1564
1565 delete [] stats;
1566 delete [] tsp_edgelist;
1567 delete [] tsp_x;
1568 delete [] sel;
1569
1570 return(cuts_added);
1571
1572 }
1573
1574 /*===========================================================================*/
1575
1576 int
addTspCuts(VrpModel * m,BcpsConstraintPool & conPool,CCtsp_lpcut_in ** tsp_cuts,CCtsp_lpgraph * g)1577 VrpCutGenerator::addTspCuts(VrpModel *m, BcpsConstraintPool &conPool,
1578 CCtsp_lpcut_in **tsp_cuts, CCtsp_lpgraph *g)
1579 {
1580 #if 1
1581 int clique_size = (m->vertnum_ >> DELETE_POWER) + 1;
1582 char *clique_array, *clique_set;
1583 int i, j, k, size, cliquecount, val, *matind, nzcnt;
1584 double *matval, rhs;
1585 BlisConstraint *blisCon = NULL;
1586 int num_cuts = 0;
1587 CCtsp_lpcut_in *tsp_cut, *tsp_cut_next;
1588 int v0, v1, jj, edgenum = m->edgenum_;;
1589 std::vector<VrpVariable *>edges = model_->getEdgeList();
1590 double infinity = model_->solver()->getInfinity();
1591
1592 for (tsp_cut = *tsp_cuts; tsp_cut; tsp_cut = tsp_cut->next){
1593 cliquecount = tsp_cut->cliquecount;
1594 size = cliquecount * clique_size;
1595 rhs = (cliquecount == 1 ? 0.0 : -((double)cliquecount)/2.0 + 1.0);
1596 clique_set = clique_array = new char[size];
1597 memset(clique_array, 0, size);
1598
1599 for (i = 0; i < cliquecount; i++, clique_set += clique_size){
1600 for(j = 0; j < tsp_cut->cliques[i].segcount; j++){
1601 for(k = tsp_cut->cliques[i].nodes[j].lo;
1602 k <= tsp_cut->cliques[i].nodes[j].hi; k++){
1603 rhs++;
1604 clique_set[k >> DELETE_POWER] |= (1 << (k & DELETE_AND));
1605 }
1606 }
1607 /*For each tooth, we want to add |T|-1 to the rhs so we have to
1608 subtract off the one here. It subtracts one for the handle too
1609 but that is compensated for above*/
1610 rhs--;
1611 }
1612 matind = new int[cliquecount*edgenum];
1613 matval = new double[cliquecount*edgenum];
1614 for (nzcnt = 0, i = 0; i < edgenum; i++){
1615 v0 = edges[i]->getv0();
1616 v1 = edges[i]->getv1();
1617 val = 0;
1618 for (jj = 0; jj < cliquecount; jj++){
1619 clique_set = clique_array + clique_size * jj;
1620 if (clique_set[v0 >> DELETE_POWER] &
1621 (1 << (v0 & DELETE_AND)) &&
1622 (clique_set[v1 >> DELETE_POWER]) &
1623 (1 << (v1 & DELETE_AND))){
1624 val += 1;
1625 }
1626 }
1627 if (val){
1628 matind[nzcnt] = i;
1629 matval[nzcnt++] = val;
1630 }
1631 }
1632
1633 if (tsp_cut->sense == 'L'){
1634 blisCon = new BlisConstraint(-infinity, (double) rhs,
1635 -infinity, (double) rhs, nzcnt,
1636 matind, matval);
1637 }else{
1638 blisCon = new BlisConstraint((double) rhs, infinity,
1639 (double) rhs, infinity, nzcnt,
1640 matind, matval);
1641 }
1642
1643 conPool.addConstraint(blisCon);
1644 num_cuts++;
1645
1646 delete [] matind;
1647 delete [] matval;
1648 delete [] clique_array;
1649 }
1650
1651 #else
1652
1653 int nzlist, nzcnt;
1654 CCtsp_lpcut new_cut;
1655 int rval = 0, rhs, i, *matind;
1656 double *matval;
1657 BlisConstraint *blisCon = NULL;
1658 int num_cuts = 0;
1659 CCtsp_lpcut_in *tsp_cut;
1660 double infinity = model_->solver()->getInfinity();
1661
1662 for (tsp_cut = *tsp_cuts; tsp_cut; tsp_cut = tsp_cut->next){
1663
1664 CCtsp_init_lpcut (&new_cut);
1665
1666 new_cut.rhs = tsp_cut->rhs;
1667 new_cut.sense = tsp_cut->sense;
1668 new_cut.branch = tsp_cut->branch;
1669
1670 nzlist = CCtsp_lpcut_in_nzlist (g, tsp_cut);
1671
1672 rval = CCtsp_copy_skeleton (&tsp_cut->skel, &new_cut.skel);
1673
1674 rhs = new_cut.rhs;
1675 for (i=0; i<new_cut.modcount; i++) {
1676 rhs += 2*(((int) new_cut.mods[i].mult) - 128);
1677 }
1678
1679 nzcnt = 0;
1680 for (i=nzlist; i != -1; i = g->edges[i].coefnext) {
1681 if (g->edges[i].coef) nzcnt++;
1682 }
1683
1684 if (nzcnt != 0) {
1685 matind = new int[nzcnt];
1686 matval = new double[nzcnt];
1687 for (nzcnt = 0; nzlist != -1; nzlist = i) {
1688 i = g->edges[nzlist].coefnext;
1689 g->edges[nzlist].coefnext = -2;
1690 if (g->edges[nzlist].coef) {
1691 matind[nzcnt] = nzlist;
1692 matval[nzcnt] = g->edges[nzlist].coef;
1693 g->edges[nzlist].coef = 0;
1694 nzcnt++;
1695 }
1696 }
1697 if (new_cut.sense == 'L'){
1698 blisCon = new BlisConstraint(-infinity, (double) rhs,
1699 -infinity, (double) rhs, nzcnt,
1700 matind, matval);
1701 }else{
1702 blisCon = new BlisConstraint((double) rhs, infinity,
1703 (double) rhs, infinity, nzcnt,
1704 matind, matval);
1705 }
1706
1707 num_cuts++;
1708
1709 blisCon->setValidRegion(BcpsValidGlobal);
1710
1711 delete [] matind;
1712 delete [] matval;
1713 } else {
1714 printf ("WARNING: Adding an empty cut\n");
1715 }
1716 }
1717
1718 #endif
1719
1720 return(num_cuts);
1721
1722 }
1723 #endif
1724