00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef VRP_CONCORDE_INCLUDED
00014 #define VRP_CONCORDE_INCLUDED
00015
00016
00017
00018
00019
00020
00021 #include "VRP_ConcordeI.h"
00022
00023
00024 #define VRP_CONCORDE_DEBUG 0
00025 #define VRP_CONCORDE_SOLVER_TINY 1
00026
00027
00028
00029
00030 class ConcordeGraph {
00031 public:
00032 int m_nEdges;
00033 int m_nVerts;
00034 int * m_edgeList;
00035 int * m_edgeValue;
00036
00037
00038
00039
00040
00041
00042 public:
00043 ConcordeGraph() :
00044 m_nEdges (0),
00045 m_nVerts (0),
00046 m_edgeList (0),
00047 m_edgeValue (0)
00048
00049
00050 {};
00051
00052 void clear(){
00053 m_nEdges = 0;
00054 m_nVerts = 0;
00055 UTIL_DELARR(m_edgeList);
00056 UTIL_DELARR(m_edgeValue);
00057
00058
00059 }
00060
00061 void setEdgeValue(const double * edgeValue){
00062
00063 }
00064
00065 void init(const int nVerts,
00066 const int nEdges){
00067 m_nVerts = nVerts;
00068 m_nEdges = nEdges;
00069 if(nEdges > 0){
00070 m_edgeList = new int[2 * m_nEdges];
00071 m_edgeValue = new int[ m_nEdges];
00072
00073
00074
00075
00076 CoinAssertHint(m_edgeList && m_edgeValue,
00077
00078 "Error: Out of Memory");
00079 }
00080 }
00081
00082 ~ConcordeGraph() {
00083 clear();
00084 };
00085 };
00086
00087
00088
00089 class VRP_Concorde {
00090 public:
00091 const VRP_Instance * m_vrp;
00092 ConcordeGraph m_cgExpand;
00093 int * m_CCoptTour;
00094 double m_CCoptVal;
00095
00096 double * tmpDblNOrigE;
00097
00098 int * tmpIntNOrigE;
00099
00100 public:
00101 VRP_Concorde() :
00102 m_vrp (0),
00103 m_cgExpand ( ) ,
00104 m_CCoptTour (0),
00105 m_CCoptVal ( ) ,
00106 tmpDblNOrigE( ) ,
00107 tmpIntNOrigE( ) {}
00108
00109 ~VRP_Concorde() {
00110 UTIL_DELARR(m_CCoptTour);
00111 UTIL_DELARR(tmpDblNOrigE);
00112 UTIL_DELARR(tmpIntNOrigE);
00113 }
00114
00115 public:
00116
00118 void init(const VRP_Instance * vrp) {
00119 m_vrp = vrp;
00120 buildExpandedCompleteGraph();
00121 m_CCoptTour = new int[m_cgExpand.m_nVerts];
00122 CoinAssertHint(m_CCoptTour, "Error: Out of Memory");
00123
00124 const UtilGraphLib & graphLib = m_vrp->m_graphLib;
00125 tmpDblNOrigE = new double[graphLib.n_edges];
00126 tmpIntNOrigE = new int [graphLib.n_edges];
00127 CoinAssertHint(tmpDblNOrigE && tmpIntNOrigE,
00128 "Error: Out of Memory");
00129 }
00130
00131
00132 void buildExpandedCompleteGraph(){
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154 int i, j, edgeIndex;
00155 const UtilGraphLib & graphLib = m_vrp->m_graphLib;
00156 const int nRoutes = m_vrp->m_numRoutes;
00157 int nCustomers = graphLib.n_vertices - 1;
00158 int nVerts = nCustomers + nRoutes;
00159 int nEdges = UtilNumEdgesU(nVerts);
00160 int * edgeList = NULL;
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170 m_cgExpand.init(nVerts, nEdges);
00171 edgeIndex = 0;
00172
00173 edgeList = m_cgExpand.m_edgeList;
00174
00175
00176 for(i = 1; i < nVerts; i++){
00177 for(j = 0; j < i; j++){
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187 edgeList[2*edgeIndex ] = i;
00188 edgeList[2*edgeIndex+1] = j;
00189 CoinAssert(edgeIndex < nEdges);
00190 edgeIndex++;
00191
00192
00193 }
00194 }
00195
00196
00197
00198
00199 }
00200
00201 void setExpandedCost(const double * origGraphCost){
00202
00203
00204
00205
00206
00207
00208
00209 const UtilGraphLib & graphLib = m_vrp->m_graphLib;
00210 const int nOrigEdges = graphLib.n_edges;
00211 const int nOrigVerts = graphLib.n_vertices;
00212 const int nRoutes = m_vrp->m_numRoutes;
00213 const int nEdges = m_cgExpand.m_nEdges;
00214 const int nCustomers = nOrigVerts - 1;
00215 double * origGraphCostDbl = tmpDblNOrigE;
00216 int * origGraphCostInt = tmpIntNOrigE;
00217 int * edgeValue = m_cgExpand.m_edgeValue;
00218 int i;
00219
00220
00221
00222 #if VRP_CONCORDE_DEBUG > 0
00223 for(i = 0; i < nOrigEdges; i++){
00224 printf("origGraphCost[%d -> %s]: %10.5f\n",
00225 i, UtilEdgeToStr(i).c_str(), origGraphCost[i]);
00226 }
00227 #endif
00228 const double * minElement = min_element(origGraphCost,
00229 origGraphCost + nOrigEdges);
00230 const double offset = min(0.0, *minElement);
00231 #if VRP_CONCORDE_DEBUG > 0
00232 printf("minElement = %10.5f offset = %10.5f\n", *minElement, -offset);
00233 #endif
00234 memcpy(origGraphCostDbl, origGraphCost, nOrigEdges * sizeof(double));
00235 if(offset < 0.0){
00236 UtilAddOffsetArr(nOrigEdges, -offset, origGraphCostDbl);
00237 }
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 UtilScaleDblToIntArr(nOrigEdges,
00249 origGraphCostDbl,
00250 origGraphCostInt, 1.0e-3);
00251
00252
00253
00254 #if VRP_CONCORDE_DEBUG > 0
00255
00256 for(i = 0; i < nOrigEdges; i++){
00257 printf("origGraphCost[%d]: %10.5f -> %10d\n",
00258 i,
00259 origGraphCostDbl[i],
00260 origGraphCostInt[i]);
00261
00262 }
00263 #endif
00264
00265
00266
00267
00268
00269 int bigCost = 0;
00270 double bigCostDbl = 0.0;
00271 int maxEdgeLen = *max_element(origGraphCostInt,
00272 origGraphCostInt+nOrigEdges);
00273 bigCostDbl = ((double)maxEdgeLen+1.0) * (double)m_cgExpand.m_nVerts;
00274 if((256*bigCostDbl) > (double)CCutil_MAXINT) {
00275 printf("WARNING: Large edge lengths!!!!\n");
00276 bigCost = CCutil_MAXINT / 256;
00277 }
00278 else
00279 bigCost = (maxEdgeLen+1) * m_cgExpand.m_nVerts;
00280
00281 while(maxEdgeLen > bigCost){
00282 printf("maxEdgeLen=%d > bigCost=%d !!! SCALE BY 10\n",
00283 maxEdgeLen, bigCost);
00284
00285 for(i = 0; i < nOrigEdges; i++){
00286 origGraphCostInt[i] /= 10;
00287 }
00288 maxEdgeLen /= 10;
00289 }
00290 assert(maxEdgeLen < bigCost);
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302 int origDepot = 0;
00303 int origCost = 0;
00304 int u, v, d, d1, d2, origIndex, newIndex, custIndex, depot;
00305 for(u = 1; u < nOrigVerts; u++){
00306 for(v = 1; v < u; v++){
00307 origIndex = UtilIndexU(u,v);
00308 newIndex = UtilIndexU(u-1, v-1);
00309 #if VRP_CONCORDE_DEBUG > 0
00310 printf("origIndex: [%d -> %s] to newIndex: [%d -> %s]\n",
00311 origIndex, UtilEdgeToStr(origIndex).c_str(),
00312 newIndex, UtilEdgeToStr(newIndex).c_str());
00313 #endif
00314 CoinAssert(origGraphCostInt[origIndex] >= 0);
00315 CoinAssert(origGraphCostInt[origIndex] < bigCost);
00316 edgeValue[newIndex] = origGraphCostInt[origIndex];
00317 }
00318 }
00319 for(u = 1; u <= nCustomers; u++){
00320
00321 origIndex = UtilIndexU(u, origDepot);
00322 origCost = origGraphCostInt[origIndex];
00323
00324 custIndex = u - 1;
00325
00326 for(d = 0; d < nRoutes; d++){
00327 depot = nCustomers + d;
00328
00329
00330 newIndex = UtilIndexU(custIndex, depot);
00331 CoinAssert(newIndex < nEdges);
00332
00333 #if VRP_CONCORDE_DEBUG > 0
00334 printf("(u:%d,0) origIndex=%d origCost=%d d=%d cInd=%d nInd=%d\n",
00335 u, origIndex, origCost, d,
00336 UtilIndexU(custIndex, depot), newIndex);
00337 #endif
00338 edgeValue[newIndex] = origCost;
00339 }
00340 }
00341 for(d1 = 1; d1 < nRoutes; d1++){
00342 for(d2 = 0; d2 < d1; d2++){
00343 newIndex = UtilIndexU(nCustomers + d1,
00344 nCustomers + d2);
00345 edgeValue[newIndex] = bigCost;
00346 }
00347 }
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360 }
00361
00362 int solveTSP(vector<int> & vrpRouteInd,
00363 vector<double> & vrpRouteEls){
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 int returnCode = 0;
00410 int silent = 2;
00411 int success, foundtour, hit_timebound;
00412
00413 CCrandstate rstate;
00414 int seed = 1;
00415 CCutil_sprand(seed, &rstate);
00416
00417
00418
00419 int i, d1, d2;
00420 const UtilGraphLib & graphLib = m_vrp->m_graphLib;
00421 int nRoutes = m_vrp->m_numRoutes;
00422 int nOrigVerts = graphLib.n_vertices;
00423 int nCustomers = nOrigVerts - 1;
00424 double upbound = DecompInf;
00425 double optval = DecompInf;
00426
00427
00428
00429
00430
00431 int * lower = new int[m_cgExpand.m_nEdges];
00432 int * upper = new int[m_cgExpand.m_nEdges];
00433 int * xsol = new int[m_cgExpand.m_nEdges];
00434 for(i = 0; i < m_cgExpand.m_nEdges; i++){
00435 lower[i] = 0;
00436 upper[i] = 1;
00437 xsol[i] = 0;
00438 }
00439 for(d1 = 1; d1 < nRoutes; d1++){
00440 for(d2 = 0; d2 < d1; d2++){
00441 upper[UtilIndexU(nCustomers + d1, nCustomers + d2)] = 0;
00442 }
00443 }
00444 returnCode = CCtiny_bnc_msp (m_cgExpand.m_nVerts,
00445 m_cgExpand.m_nEdges,
00446 m_cgExpand.m_edgeList,
00447 m_cgExpand.m_edgeValue,
00448 -1,
00449 lower,
00450 upper,
00451 &upbound,
00452 CC_TINYTSP_MINIMIZE,
00453 &optval,
00454 xsol,
00455 1,
00456 1000);
00457
00458 assert(returnCode != CC_TINYTSP_ERROR);
00459 assert(returnCode != CC_TINYTSP_INFEASIBLE);
00460 bool useFullSolver = false;
00461 if(returnCode == CC_TINYTSP_SEARCHLIMITEXCEEDED){
00462 printf("Search limit exceeded - let's go to full solver\n");
00463 useFullSolver = true;
00464 }
00465
00466 delete [] lower;
00467 delete [] upper;
00468
00469 if(useFullSolver){
00470 CCdatagroup dataGroup;
00471 returnCode = CCutil_graph2dat_matrix(m_cgExpand.m_nVerts,
00472 m_cgExpand.m_nEdges,
00473 m_cgExpand.m_edgeList,
00474 m_cgExpand.m_edgeValue,
00475 -1, &dataGroup);
00476 char nameC[100];
00477 string name = UtilStringRandom(10);
00478 strcpy(nameC, name.c_str());
00479
00480 #if VRP_CONCORDE_DEBUG > 0
00481 printf("CONCORDE GRAPH:\n");
00482 printf(" nVerts=%d nEdges=%d\n",
00483 m_cgExpand.m_nVerts,
00484 m_cgExpand.m_nEdges);
00485 for(i = 0; i < m_cgExpand.m_nEdges; i++){
00486 printf("(%4d,%4d) val:%4d\n",
00487 m_cgExpand.m_edgeList[2*i],
00488 m_cgExpand.m_edgeList[2*i+1],
00489 m_cgExpand.m_edgeValue[i]);
00490 }
00491 #endif
00492
00493 returnCode = CCtsp_solve_dat(m_cgExpand.m_nVerts,
00494 &dataGroup,
00495 NULL,
00496 m_CCoptTour,
00497 NULL,
00498 &m_CCoptVal,
00499 &success,
00500 &foundtour,
00501 nameC,
00502 NULL,
00503 &hit_timebound,
00504 silent,
00505 &rstate);
00506 CoinAssert((success == 1) && (foundtour == 1));
00507 }
00508
00509 #if VRP_CONCORDE_DEBUG > 0
00510 printf("Optimal Route:\n");
00511 #endif
00512
00513
00514 if(!useFullSolver){
00515 vector<int> edgeListSol;
00516 for(i = 0; i < m_cgExpand.m_nEdges; i++){
00517 if(xsol[i]){
00518 #if VRP_CONCORDE_DEBUG > 0
00519 printf("(%d,%d) wt:%d\n",
00520 m_cgExpand.m_edgeList[2*i],
00521 m_cgExpand.m_edgeList[2*i+1],
00522 m_cgExpand.m_edgeValue[i]);
00523 #endif
00524 edgeListSol.push_back(m_cgExpand.m_edgeList[2*i]);
00525 edgeListSol.push_back(m_cgExpand.m_edgeList[2*i+1]);
00526 }
00527 }
00528 createVrpRouteFromTspEdgeList(edgeListSol,
00529 vrpRouteInd,
00530 vrpRouteEls);
00531 }
00532 else{
00533 createVrpRouteFromTspRoute(m_CCoptTour,
00534 m_cgExpand.m_nVerts,
00535 vrpRouteInd,
00536 vrpRouteEls);
00537 }
00538 return returnCode;
00539 }
00540
00541 void createVrpRouteFromTspEdgeList(vector<int> & tspEdgeList,
00542 vector<int> & vrpRouteInd,
00543 vector<double> & vrpRouteEls){
00544
00545 int i, u, v;
00546 int nVerts = m_vrp->m_graphLib.n_vertices;
00547 int nCustomers = m_vrp->m_graphLib.n_vertices - 1;
00548 int * depotEdges = tmpIntNOrigE;
00549
00550
00551
00552
00553 UtilFillN(depotEdges, nVerts, 0);
00554 #if VRP_CONCORDE_DEBUG > 0
00555 printf("Optimal Route:\n");
00556 #endif
00557 for(i = 0; i < m_cgExpand.m_nVerts; i++){
00558
00559 u = tspEdgeList[2*i];
00560 v = tspEdgeList[2*i+1];
00561 u = u >= nCustomers ? 0 : u+1;
00562 v = v >= nCustomers ? 0 : v+1;
00563 CoinAssert((u + v) > 0);
00564 if(u == 0){
00565 depotEdges[v]++;
00566 }
00567 else if(v == 0){
00568 depotEdges[u]++;
00569 }
00570 else{
00571 vrpRouteInd.push_back(UtilIndexU(u,v));
00572 #if VRP_CONCORDE_DEBUG > 0
00573 printf("(%d,%d) -> %d\n", u, v, UtilIndexU(u,v));
00574 #endif
00575 }
00576 }
00577 UtilFillN(vrpRouteEls, vrpRouteInd.size(), 1.0);
00578
00579 for(v = 1; v < nVerts; v++){
00580 if(depotEdges[v]){
00581 CoinAssert(depotEdges[v] >= 0);
00582 CoinAssert(depotEdges[v] <= 2);
00583 vrpRouteInd.push_back(UtilIndexU(v, 0));
00584 vrpRouteEls.push_back(static_cast<double>(depotEdges[v]));
00585 #if VRP_CONCORDE_DEBUG > 0
00586 printf("%d D(%d,%d) -> %d\n",
00587 depotEdges[v], v, 0, UtilIndexU(v,0));
00588 #endif
00589 }
00590 }
00591 }
00592
00593 void createVrpRouteFromTspRoute(const int * tspRoute,
00594 const int tspRouteLen,
00595 vector<int> & vrpRouteInd,
00596 vector<double> & vrpRouteEls){
00597
00598 int i, u, v;
00599 int nVerts = m_vrp->m_graphLib.n_vertices;
00600 int nCustomers = m_vrp->m_graphLib.n_vertices - 1;
00601 int * depotEdges = tmpIntNOrigE;
00602
00603
00604
00605
00606 UtilFillN(depotEdges, nVerts, 0);
00607 #if VRP_CONCORDE_DEBUG > 0
00608 printf("Optimal Route Nodes:");
00609 for(i = 0; i < tspRouteLen; i++){
00610 printf("%d ", tspRoute[i]);
00611 }
00612 printf("\nOptimal Route Edges:\n");
00613 #endif
00614 for(i = 0; i < tspRouteLen-1; i++){
00615 u = tspRoute[i] >= nCustomers ? 0 : tspRoute[i] +1;
00616 v = tspRoute[i+1] >= nCustomers ? 0 : tspRoute[i+1]+1;
00617 CoinAssert((u + v) > 0);
00618 if(u == 0){
00619 depotEdges[v]++;
00620 }
00621 else if(v == 0){
00622 depotEdges[u]++;
00623 }
00624 else{
00625 vrpRouteInd.push_back(UtilIndexU(u,v));
00626 }
00627 #if VRP_CONCORDE_DEBUG > 0
00628 printf("(%d,%d) -> %d\n", u, v, UtilIndexU(u,v));
00629 #endif
00630 }
00631 u = tspRoute[tspRouteLen-1] >= nCustomers ? 0
00632 : tspRoute[tspRouteLen-1]+1;
00633 v = tspRoute[0] >= nCustomers ? 0 : tspRoute[0] + 1;
00634 CoinAssert((u + v) > 0);
00635 if(u == 0){
00636 depotEdges[v]++;
00637 }
00638 else if(v == 0){
00639 depotEdges[u]++;
00640 }
00641 else{
00642 vrpRouteInd.push_back(UtilIndexU(u,v));
00643 }
00644 #if VRP_CONCORDE_DEBUG > 0
00645 printf("(%d,%d) -> %d\n", u, v, UtilIndexU(u,v));
00646 #endif
00647 UtilFillN(vrpRouteEls, vrpRouteInd.size(), 1.0);
00648
00649 for(v = 1; v < nVerts; v++){
00650 if(depotEdges[v]){
00651 CoinAssert(depotEdges[v] >= 0);
00652 CoinAssert(depotEdges[v] <= 2);
00653 vrpRouteInd.push_back(UtilIndexU(v, 0));
00654 vrpRouteEls.push_back(static_cast<double>(depotEdges[v]));
00655 #if VRP_CONCORDE_DEBUG > 0
00656 printf("%d D(%d,%d) -> %d\n",
00657 depotEdges[v], v, 0, UtilIndexU(v,0));
00658 #endif
00659 }
00660 }
00661 }
00662
00663 void createTSPLIBFile(const string fileName){
00664 ofstream os;
00665 int i, j, edgeIndex;
00666 const int nVerts = m_cgExpand.m_nVerts;
00667 const int * edgeValue = m_cgExpand.m_edgeValue;
00668
00669 printf("Open File %s\n", fileName.c_str());
00670
00671 UtilOpenFile(os, fileName);
00672 os << "NAME: " << fileName << "\n"
00673 << "TYPE: TSP\n"
00674 << "DIMENSION: " << nVerts << "\n"
00675 << "EDGE_WEIGHT_TYPE: EXPLICIT\n"
00676 << "EDGE_WEIGHT_FORMAT: LOWER_DIAG_ROW\n"
00677 << "EDGE_WEIGHT_SECTION\n";
00678 os << " 0\n";
00679 edgeIndex = 0;
00680 for(i = 1; i < nVerts; i++){
00681 for(j = 0; j < i; j++){
00682 os << edgeValue[edgeIndex] << " ";
00683 edgeIndex++;
00684 }
00685 os << " 0\n";
00686 }
00687 os << "EOF\n";
00688
00689 os.close();
00690 }
00691
00692
00693 };
00694
00695
00696 #endif