source:Blis/trunk/Blis/examples/VRP/VrpModel.cpp@1001

Last change on this file since 1001 was 1001, checked in by yanxu, 6 years ago

File size: 29.6 KB
Line
1/*===========================================================================*
2 * This file is part of a solver for the Vehicle Routing Problem             *
3 * developed using the BiCePS Linear Integer Solver (BLIS).                  *
4 *                                                                           *
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.                                 *
13 *===========================================================================*/
14
15#include <vector>
16
17#include "BlisConstraint.h"
18#include "BlisTreeNode.h"
19
20#include "VrpConstants.h"
21#include "VrpModel.h"
22#include "VrpSolution.h"
23#include "VrpVariable.h"
24
25//#############################################################################
26
27static int VrpRoundToNearest(double x)
28{
29    return (int)(x + 0.5);
30}
31
32//#############################################################################
33
34/** For parallel code, only the master calls this function.
35 *  1) Read in the instance data
36 *  2) Set colMatrix_, varLB_, varUB_, conLB_, conUB
37 *     numCols_, numRows_
38 *  3) Set objCoef_ and objSense_
39 *  4) Set colType_ ('C', 'I', or 'B')
40 *  5) Create variables and constraints
41 *  6) Set numCoreVariables_ and numCoreConstraints_
42 */
43void
45{
46
47    int msgLevel = AlpsPar_->entry(AlpsParams::msgLevel);
48
49   static char keywords[KEY_NUM][22] = {
50      "NAME",
51      "NAME:",                 /* This section lists the names of the */
52      "TYPE",                  /* possible fields in the data file    */
53      "TYPE:",
54      "COMMENT",
55      "COMMENT:",
56      "DIMENSION",
57      "DIMENSION:",
58      "CAPACITY",
59      "CAPACITY:",
60      "EDGE_WEIGHT_TYPE",
61      "EDGE_WEIGHT_TYPE:",
62      "EDGE_WEIGHT_FORMAT",
63      "EDGE_WEIGHT_FORMAT:",
64      "DISPLAY_DATA_TYPE",
65      "DISPLAY_DATA_TYPE:",
66      "EDGE_WEIGHT_SECTION",
67      "EDGE_WEIGHT_SECTION:",
68      "DISPLAY_DATA_SECTION",
69      "DISPLAY_DATA_SECTION:",
70      "NODE_COORD_SECTION",
71      "NODE_COORD_SECTION:",
72      "NODE_COORD_TYPE",
73      "NODE_COORD_TYPE:",
74      "DEPOT_SECTION",
75      "DEPOT_SECTION:",
76      "CAPACITY_VOL",
77      "CAPACITY_VOL:",
78      "DEMAND_SECTION",
79      "DEMAND_SECTION:",
80      "TIME_WINDOW_SECTION",
81      "TIME_WINDOW_SECTION:",
82      "STANDTIME_SECTION",
83      "STANDTIME_SECTION:",
84      "PICKUP_SECTION",
85      "PICKUP_SECTION:",
86      "EOF",
87      "EOF.",
88      "NUMBER_OF_TRUCKS",
89      "NUMBER_OF_TRUCKS:",
90      "",
91      "",
92      "NO_MORE_TYPE"
93   };
94
95#define NCTYPE_NUM 3
96
97   static char nctypes[NCTYPE_NUM][14] = {
98      "TWOD_COORDS",
99      "THREED_COORDS",     /*This section lists the possible node*/
100      "NO_COORDS"          /*coordinate data types               */
101   };
102
103#define WTYPE_NUM 10
104
105   static char wtypes[WTYPE_NUM][9] = {
106      "EXPLICIT",
107      "EUC_2D",            /*This is a list of the possible data types for */
108      "EUC_3D",            /*edge weights                                  */
109      "MAX_2D",
110      "MAX_3D",
111      "MAN_2D",
112      "MAN_3D",
113      "CEIL_2D",
114      "GEO",
115      "ATT"
116   };
117
118#define WFORMAT_NUM 9
119
120   static char wformats[WFORMAT_NUM][20] = {
121      "UPPER_ROW",
122      "LOWER_ROW",          /*This is a list of the possible formats that*/
123      "UPPER_DIAG_ROW",     /*the edge weight matrix could be given in   */
124      "LOWER_DIAG_ROW",     /*if it is given explicitly                  */
125      "UPPER_COL",
126      "LOWER_COL",
127      "UPPER_DIAG_COL",
128      "LOWER_DIAG_COL",
129      "FULL_MATRIX"
130   };
131
132#define DTYPE_NUM 3
133
134   static char dtypes[DTYPE_NUM][14] = {
135      "COORD_DISPLAY",
136      "TWOD_DISPLAY",     /*This is a list of the various display data*/
137      "NO_DISPLAY"        /*types                                     */
138   };
139
140   char line[LENGTH], line1[LENGTH], key[30], tmp[80];
141   int wformat=-1, dtype=-1, nctype=-1;
142   double fdummy;
143   int i, j = 0;
144   int l, m;
145   FILE *f;
146   int node;
147   double deg, min, coord_x, coord_y, coord_z;
148   double x, y;
149   int capacity_vol = false;
150   int k;
151   numroutes_ = VrpPar_->entry(VrpParams::numRoutes);
152
153   if (!strcmp(dataFile, "")){
154      printf("\nVrp I/O: No problem data file specified\n\n");
155      exit(1);
156   }
157
158   if ((f = fopen(dataFile, "r")) == NULL){
159      fprintf(stderr, "Vrp I/O: file '%s' can't be opened\n", dataFile);
160      exit(1);
161   }
162
163   /*This loop reads in the next line of the data file and compares it
164     to the list of possible keywords to determine what data will follow.
165     It then reads the data into the appropriate field and iterates */
166
167   while(NULL != fgets( line1, LENGTH, f)){
168      strcpy(key,"");
169      sscanf(line1,"%s",key); /*read in next keyword*/
170
171      for (k = 0; k < KEY_NUM; k++) /*which data field comes next?*/
172         if (strcmp(keywords[k], key) == 0) break;
173
174      if (k == KEY_NUM){
175         continue;
176         fprintf(stderr, "Unknown keyword! bye.\n");
177         exit(1); /*error check for acceptable data field*/
178      }
179
180      k >>= 1; /* This is a bit shift operation that divides k by 2    */
181      /* since in the list of keywords, there are two possible*/
182      /* formats for the keyword                              */
183
184      if (strchr(line1,':')){
185         strcpy(line, strchr(line1, ':')+1);
186      }
187
188      switch (k){
189
190       case 0: /* NAME */
191         if (!sscanf(line, "%s", name_))
192            fprintf(stderr, "\nVrp I/O: error reading NAME\n\n");
193
194         if (msgLevel > 0) {
195             printf("PROBLEM NAME: \t\t%s\n", name_);
196         }
197         break;
198       case 1 : /*TYPE*/
199         sscanf(line, "%s", tmp);
200         if (strcmp("CVRP", tmp) != 0){
201           if (strcmp("TSP", tmp) == 0){
202              VrpPar_->setEntry(VrpParams::tspProb, true);
203           }else{
204              fprintf(stderr, "This is not a recognized file type!\n");
205              exit(1);
206           }
207         }
208         if (msgLevel > 0) {
209             printf("TYPE: \t\t\t%s\n", tmp);
210         }
211       case 2 : /*COMMENT*/
212#if 0
213         if (!strncpy(tmp, line, 80))
214            fprintf(stderr, "\nVrp I/O: error reading COMMENT\n\n");
215         printf("DESCRIPTION: \t\t%s\n", tmp);
216#endif
217         break;
218       case 3 : /* DIMENSION */
219         if (!sscanf(line, "%i", &k)){
220            fprintf(stderr, "Vrp I/O: error reading DIMENSION\n\n");
221            exit(1);
222         }
223         vertnum_ = (int) k;
224         edgenum_ = (int) vertnum_ * (vertnum_ - 1)/2;
225         if (msgLevel > 0) {
226             printf("DIMENSION: \t\t%i\n", k);
227         }
228         break;
229       case 4 : /*CAPACITY*/
230         if (!sscanf(line, "%i", &k)){
231            fprintf(stderr, "Vrp I/O: error reading CAPACITY\n\n");
232            exit(1);
233         }
234         capacity_ = (int) k;
235         break;
236       case 5 : /* EDGE_WEIGHT_TYPE */
237         sscanf(line, "%s", tmp);
238         for (wtype_ = 0; wtype_ < WTYPE_NUM; (wtype_)++)
239            if (strcmp(wtypes[wtype_], tmp) == 0) break;
240         if (wtype_ == WTYPE_NUM) {
241            fprintf(stderr, "Unknown weight type : %s !!!\n", tmp);
242            exit(1);
243         }
244         break;
245       case 6 : /* EDGE_WEIGHT_FORMAT */
246         sscanf(line, "%s", tmp);
247         for (wformat = 0; wformat < WFORMAT_NUM; wformat++)
248            if (strcmp(wformats[wformat], tmp) == 0) break;
249         if (wformat == WFORMAT_NUM) {
250            fprintf(stderr, "Unknown weight type : %s !!!\n", tmp);
251            exit(1);
252         }
253         break;
254       case 7 : /* DISPLAY_DATA_TYPE */
255         sscanf(line, "%s", tmp);
256         for (dtype = 0; dtype < DTYPE_NUM; dtype++)
257            if (strcmp(dtypes[dtype], tmp) == 0) break;
258         if (dtype == DTYPE_NUM) {
259            fprintf(stderr, "Unknown display type : %s !!!\n", tmp);
260            exit(1);
261         }
262         break;
263       case 8: /* EDGE_WEIGHT_SECTION */
264         /*------------------------break if not EXPLICIT -*/
265         if (wtype_ != _EXPLICIT) break;
266         switch (wformat){
267          case 1 : /* LOWER_ROW */
268          case 4 : /* UPPER_COL */
269          case 3 : /* LOWER_DIAG_ROW */
270          case 6 : /* UPPER_DIAG_COL */
271            for (i = 0; i < vertnum_; i++){
272               for (j = 0; j < i; j++){
273                  if (!fscanf(f,"%lf", &fdummy)){
274                     fprintf(stderr, "Not enough data -- DIMENSION or "
275                             "EDGE_WEIGHT_TYPE declared wrong\n");
276                     exit(1);
277                  } else {
278                     edges_.push_back(new VrpVariable(i, j, (int) fdummy));
279                  }
280               }
281               if ((wformat==3 || wformat==6) &&
282                   !fscanf(f,"%lf", &fdummy)){
283                  fprintf(stderr, "Not enough data -- DIMENSION or "
284                          "EDGE_WEIGHT_TYPE declared wrong\n");
285                  exit(1);
286               }
287            }
288            if (fscanf(f,"%lf", &fdummy)){
289               fprintf(stderr, "Too much data -- DIMENSION or "
290                       "EDGE_WEIGHT_TYPE declared wrong\n");
291               exit(1);
292            }
293            break;
294         case 0 : /* UPPER_ROW */
295         case 2 : /* UPPER_DIAG_ROW */
296         {
297
298             // Deem not to be able solve large problem
299             int ** lowDiag = new int* [vertnum_];
300             for (i = 0; i < vertnum_; ++i) {
301                 lowDiag[i] = new int [vertnum_];
302             }
303             for (i = 0; i < vertnum_; i++){
304                 if (wformat == 2) {
305                     if (!fscanf(f,"%lf", &fdummy)){
306                         fprintf(stderr, "Not enough data -- DIMENSION or "
307                                 "EDGE_WEIGHT_TYPE declared wrong");
308                         exit(1);
309                     }
310                 }
311                 for (j= i + 1; j < vertnum_; j++){
312                     if (!fscanf(f,"%lf", &fdummy)){
313                         fprintf(stderr, "Not enough data -- DIMENSION or "
314                                 "EDGE_WEIGHT_TYPE declared wrong");
315                         exit(1);
316                     } else {
317                         // Create lower diag
318                         lowDiag[j][i] = (int) fdummy;
319                     }
320                 }
321             }
322             if (fscanf(f,"%lf", &fdummy)){
323                 fprintf(stderr, "Too much data -- DIMENSION or "
324                         "EDGE_WEIGHT_TYPE declared wrong\n");
325                 exit(1);
326             }
327             // Create edges
328             for (i = 1, k = 0; i < vertnum_; i++){
329                 for (j = 0; j < i; j++){
330                     VrpVariable *aVar = new VrpVariable(i, j, lowDiag[i][j]);
331                     edges_.push_back(aVar);
332                 }
333             }
334             // Free lowDiag
335             for (i = 0; i < vertnum_; ++i) {
336                 delete [] lowDiag[i];
337             }
338             delete [] lowDiag;
339             break;
340         }
341         case 5 : /* LOWER_COL */
342         case 7 : /* LOWER_DIAG_COL */
343            for (i = 0; i < vertnum_; i++){
344               if (wformat==7)
345                  if (!fscanf(f,"%lf", &fdummy)){
346                     fprintf(stderr, "Not enough data -- DIMENSION or "
347                             "EDGE_WEIGHT_TYPE declared wrong");
348                     exit(1);
349                  }
350               for (j= i + 1; j < vertnum_; j++){
351                  if (!fscanf(f,"%lf", &fdummy)){
352                     fprintf(stderr, "Not enough data -- DIMENSION or "
353                             "EDGE_WEIGHT_TYPE declared wrong");
354                     exit(1);
355                  } else {
356                          edges_.push_back(new VrpVariable(i, j, (int) fdummy));
357                  }
358               }
359            }
360            if (fscanf(f,"%lf", &fdummy)){
361               fprintf(stderr, "Too much data -- DIMENSION or "
362                       "EDGE_WEIGHT_TYPE declared wrong\n");
363               exit(1);
364            }
365            break;
366          case 8 : /* FULL_MATRIX */
367            for (i = 0; i < vertnum_; i++){
368                for (j = 0; j < i; j++){
369                    if(!fscanf(f,"%lf", &fdummy)){
370                        fprintf(stderr, "Not enough data -- DIMENSION or "
371                                "EDGE_WEIGHT_TYPE declared wrong");
372                        exit(1);
373                    }
374                    edges_.push_back(new VrpVariable(i, j, (int) fdummy));
375                    //edges_[index(i, j)] = new VrpVariable(i, j, (int) fdummy);
376                }
377               for (j = i; j < vertnum_; j++){
378                   if(!fscanf(f,"%lf", &fdummy)){
379                       fprintf(stderr, "Not enough data -- DIMENSION or "
380                               "EDGE_WEIGHT_TYPE declared wrong");
381                       exit(1);
382                   }
383               }
384#if 0 // BUG: can read upper diag
385               for (j = 0; j <= i; j++)
386                  if(!fscanf(f,"%lf", &fdummy)){
387                     fprintf(stderr, "Not enough data -- DIMENSION or "
388                             "EDGE_WEIGHT_TYPE declared wrong");
389                     exit(1);
390                  }
391               for (j = i + 1; j < vertnum_; j++){
392                  if(!fscanf(f,"%lf", &fdummy)){
393                     fprintf(stderr, "Not enough data -- DIMENSION or "
394                             "EDGE_WEIGHT_TYPE declared wrong");
395                     exit(1);
396                  }
397                  edges_.push_back(new VrpVariable(i, j, (int) fdummy));
398                  //edges_[index(i, j)] = new VrpVariable(i, j, (int) fdummy);
399               }
400#endif
401
402            }
403            if (fscanf(f,"%lf", &fdummy)){
404               fprintf(stderr, "Too much data -- DIMENSION or "
405                       "EDGE_WEIGHT_TYPE declared wrong\n");
406               exit(1);
407            }
408            break;
409         }
410         break;
411       case 9 : /* DISPLAY_DATA_SECTION */
412         /*--------------------- break if NO_DISPLAY -*/
413         if (dtype != 1){
414            fprintf(stderr, "DISPLAY_DATA_SECTION exists"
415                    "but not TWOD_DISPLAY!\n");
416            exit(1);
417         }
418         /* posx, posy -*/
419         posx_ = new int[vertnum_];
420         posy_ = new int[vertnum_];
421         for (i = 0; i < vertnum_; i++){
422            if ((k = fscanf(f,"%i%lf%lf", &node, &x, &y)) != 3){
423               fprintf(stderr, "\nVrp I/O: error reading DISPLAY_DATA\n");
424               break;
425            }
426            posx_[node-1] = (int)(x + 0.5);
427            posy_[node-1] = (int)(y + 0.5);
428         }
429         if (fscanf(f,"%lf", &fdummy)){
430            fprintf(stderr, "\nVrp I/O: too much display data\n");
431            break;
432         }
433         break;
434       case 10 : /* NODE_COORD_SECTION */
435         if (nctype == -1) nctype = 0;  /*if not given: TWOD_COORDS*/
436         if (dtype == -1 && ((wtype_ == _EUC_2D) || /*display type*/
437                             (wtype_ == _MAX_2D) ||      /*not defd yet*/
438                             (wtype_ == _MAN_2D)   ))/*&& can disp.*/
439            dtype = 0;                              /* COORD_DISPLAY */
440         if (dtype == 0){
441            posx_ = new int[vertnum_];
442            posy_ = new int[vertnum_];
443         }
444         coordx_ = new double[vertnum_];
445         coordy_ = new double[vertnum_];
446         if (nctype == 1)
447            coordz_ = new double[vertnum_];
448         for (i=0; i<vertnum_; i++){
449            if (nctype == 0)         /* TWOD_COORDS */
450               if (fscanf(f,"%i%lf%lf", &node, &coord_x, &coord_y) != 3){
451                  fprintf(stderr, "\nVrp I/O: error reading NODE_COORD\n\n");
452                  exit(1);
453               }
454            if (nctype == 1)         /* THREED_COORDS */
455               if (fscanf(f,"%i%lf%lf%lf", &node, &coord_x, &coord_y,
456                          &coord_z) != 4){
457                  fprintf(stderr, "\nVrp I/O: error reading NODE_COORD\n\n");
458                  exit(1);
459               }
460            coordx_[node-1] = coord_x;
461            coordy_[node-1] = coord_y;
462            /*since position is an integer and coord is a double, I must
463              round off here if dtype is EXPLICIT*/
464            if (dtype == 0){
465               posx_[node-1] = (int)coord_x;
466               posy_[node-1] = (int)coord_y;
467            }
468            if (nctype == 1) coordz_[node-1] = coord_z;
469            if (wtype_ == _GEO){ /* GEO */
470               /*--- latitude & longitude for node ------------*/
471                deg = (int)(coordx_[node-1]);
472                //deg = VrpRoundToNearest(coordx_[node-1]);
473               min = coordx_[node-1] - deg;
474               coordx_[node-1] = MY_PI * (deg + 5.0*min/3.0 ) / 180.0;
475               deg = (int)(coordy_[node-1]);
476               //deg = floor(coordy_[node-1]);
477               //deg = VrpRoundToNearest(coordy_[node-1]);
478               min = coordy_[node-1] - deg;
479               coordy_[node-1] = MY_PI * (deg + 5.0*min/3.0 ) / 180.0;
480            }
481         }
482         if (fscanf(f,"%i%lf%lf%lf", &node, &coord_x, &coord_y, &coord_z)){
483            fprintf(stderr, "\nVrp I/O: too much data in NODE_COORD\n\n");
484            exit(1);
485         }
486         break;
487       case 11: /* NODE_COORD_TYPE */
488         sscanf(line, "%s", tmp);
489         for (nctype = 0; nctype < NCTYPE_NUM; nctype++)
490            if (strcmp(nctypes[nctype], tmp) == 0) break;
491         if (nctype == NCTYPE_NUM) {
492            fprintf(stderr, "Unknown node_coord_type : %s !!!\n", tmp);
493            exit(1);
494         }
495         break;
496       case 12: /*DEPOT_SECTION*/
497         fscanf(f, "%i", &k);
498         if (k != 1){
499            fprintf(stderr, "Error in data: depot must be node 1");
500            exit(1);
501         }
502         depot_ = k - 1;
503         while (-1 != k) fscanf(f, "%i", &k);
504         break;
505       case 13: /*CAPACITY_VOL*/
506         sscanf(line, "%i", &k);
507         capacity_vol = true;
508         break;
509       case 14: /*DEMAND_SECTION*/
510         demand_ = new int[vertnum_];
511         for (i = 0; i < vertnum_; i++){
512            if (capacity_vol){
513               if (fscanf(f, "%i%i%i", &k, &l, &m) != 3){
514                  fprintf(stderr,"\nVrp I/O: error reading DEMAND_SECTION\n\n");
515                  exit(1);
516               }
517            }
518            else if (fscanf(f, "%i%i", &k, &l) != 2){
519               fprintf(stderr, "\nVrp I/O: error reading DEMAND_SECTION\n\n");
520               exit(1);
521            }
522            demand_[k-1] = l;
523            demand_[0] += l;
524         }
525         if (fscanf(f, "%i%i", &k, &l)){
526            fprintf(stderr, "\nVrp I/O: too much data in DEMAND_SECTION\n\n");
527            exit(1);
528         }
529         break;
530       case 15: /*TIME_WINDOW_SECTION*/ /*These sections are not used*/
531         while (fscanf(f, "%d %*d:%*d %*d:%*d", &k));
532         break;
533       case 16: /*STANDTIME_SECTION*/
534         while (fscanf(f, "%d%*d", &k));
535         break;
536       case 17: /*PICKUP_SECTION*/
537         while (fscanf(f, "%d%*d%*d", &k));
538         break;
539       case 18: /*  EOF */
540         break;
541       case 19: /*  NUMBER_OF_TRUCKS  */
542         if (!sscanf(line, "%i", &k)){
543            fprintf(stderr, "Vrp I/O: error reading NO_OF_TRUCKS\n\n");
544            exit(1);
545         }
546         numroutes_ = (int) k;
547       default:
548         break;
549      }
550   }
551
552   if (f != stdin)
553      fclose(f);
554
555   /*calculate all distances explicitly and then use distance type EXPLICIT*/
556   if (wtype_ != _EXPLICIT){
557      for (i = 1, k = 0; i < vertnum_; i++){
558         for (j = 0; j < i; j++){
559             VrpVariable *aVar = new VrpVariable(i, j, computeCost(i, j));
560             edges_.push_back(aVar);
561         }
562      }
563      wtype_ = _EXPLICIT;
564   }
565
566   if (VrpPar_->entry(VrpParams::tspProb)){
567      capacity_ = vertnum_;
568      numroutes_ = 1;
569      demand_ = new int[vertnum_];
570      demand_[0] = vertnum_;
571      for (i = vertnum_ - 1; i > 0; i--)
572         demand_[i] = 1;
573   }
574
575   //-------------------------------------------------------
576   // 2) Set colMatrix_, varLB_, varUB_, conLB_, conUB
577   //    numCols_, numRows_
578   // 3) Set objCoef_ and objSense_
579   // 4) Set colType_ ('C', 'I', or 'B')
580   //-------------------------------------------------------
581
582   setModelData();
583
584   //-------------------------------------------------------
585   // 5) Create variables and constraints
586   // 6) Set numCoreVariables_ and numCoreConstraints_
587   //-------------------------------------------------------
588
589   for (k = 0; k < numCols_; ++k) {
590       variables_.push_back(edges_[k]);
591   }
592
593#if 0
594   std::cout << "numCols = " << numCols_
595             << "; variables_.size() = " << variables_.size() << std::endl;
596#endif
597
598   for (k = 0; k < numRows_; ++k) {
599       BlisConstraint *con = new BlisConstraint(conLB_[k],
600                                                conUB_[k],
601                                                conLB_[k],
602                                                conUB_[k]);
603       con->setObjectIndex(k);
604       con->setRepType(BCPS_CORE);
605       con->setStatus(BCPS_NONREMOVALBE);
606       constraints_.push_back(con);
607       con = NULL;
608   }
609
610   // Set all objects as core
611   numCoreVariables_ = numCols_;
612   numCoreConstraints_ = numRows_;
613
614   // Allocate space for network for later use
615   n_ = new VrpNetwork(edgenum_, vertnum_);
616
617   VrpCutGenerator *cg = new VrpCutGenerator(this, vertnum_);
618
619   cg->setStrategy(BlisCutStrategyPeriodic);
620   cg->setCutGenerationFreq(1);
621
623}
624
625//#############################################################################
626
627CoinPackedVector *
628VrpModel::getSolution()
629{
630   // Can get LP solution information from solver;
631   int varnum = solver()->getNumCols();
632   const double *sol = solver()->getColSolution();
633   std::vector<VrpVariable *>vars = getEdgeList();
634   double etol = etol_;
635   int *indices = new int[varnum];
636   double *values = new double[varnum]; /* n */
637   int i, cnt = 0;
638
639   assert(varnum == edgenum_);
640
641   for (i = 0; i < varnum; i++){
642      if (sol[i] > etol || sol[i] < -etol){
643         indices[cnt] = vars[i]->getIndex();
644         values[cnt++] = sol[i];
645      }
646   }
647
648   return(new CoinPackedVector(varnum, cnt, indices, values, false));
649}
650
651//#############################################################################
652
653void VrpModel::createNet(CoinPackedVector *vec)
654{
655   n_->createNet(vec, demand_, getEdgeList(), etol_, vertnum_);
656}
657
658//#############################################################################
659
660BlisSolution *
661VrpModel::userFeasibleSolution(bool &userFeasible)
662{
663    CoinPackedVector *sol = getSolution();
664    VrpSolution *vrpSol = NULL;
665
666    userFeasible = true;
667
668    createNet(sol);
669
670    if (!n_->isIntegral_){
671        userFeasible = false;
672    }
673    else {
674        int rcnt = n_->connected();
675        int i;
676        for (i = 0; i < rcnt; i++){
677            if (n_->compCuts_[i+1] < 2 - etol_){
678                userFeasible = false;
679                break;
680            }
681            else if (n_->compDemands_[i+1] > capacity_){
682                userFeasible = false;
683                break;
684            }
685        }
686    }
687
688    if (userFeasible) {
689        // Create a VRP solution
690        vrpSol = new VrpSolution(getNumCols(),
691                                 getLpSolution(),
692                                 getLpObjValue() * objSense_,
693                                 this);
694
696    }
697
698    delete sol;
699
700    return vrpSol;
701}
702
703//#############################################################################
704
705int
706VrpModel::computeCost(int v0, int v1){
707
708   if (wtype_ == _EXPLICIT){
709      return((int)edges_[index(v0, v1)]->getObjCoef());
710   }
711
712   double q1, q2, q3, dx, dy, dz;
713   int cost = 0;
714
715   if (wtype_ == _GEO){
716      q1 = cos( coordy_[v0] - coordy_[v1] );
717      q2 = cos( coordx_[v0] - coordx_[v1] );
718      q3 = cos( coordx_[v0] + coordx_[v1] );
719      //cost = VrpRoundToNearest(RRR*acos(0.5*((1.0+q1)*q2-(1.0-q1)*q3))+1.0);
720      cost = (int) (RRR*acos(0.5*((1.0+q1)*q2-(1.0-q1)*q3))+1.0);
721   }else{
722      dx = coordx_[v0] - coordx_[v1];
723      dy = coordy_[v0] - coordy_[v1];
724      switch (wtype_){
725       case _EUC_2D : cost = (int) floor( sqrt( dx*dx + dy*dy ) + 0.5);
726         break;
727       case _EUC_3D : dz = coordz_[v0] - coordz_[v1];
728         cost = (int) floor( sqrt( dx*dx + dy*dy + dz*dz) + 0.5);
729         break;
730       case _MAX_2D : cost = (int) fabs(dx);
731         if (cost < fabs(dy)) cost = (int) fabs(dy);
732         break;
733       case _MAX_3D : dz = coordz_[v0] - coordz_[v1];
734         cost = (int) fabs(dx);
735         if (cost < fabs(dy)) cost = (int) fabs(dy);
736         if (cost < fabs(dz)) cost = (int) fabs(dz);
737         break;
738       case _MAN_2D : cost = (int) floor( dx+dy+0.5 );
739         break;
740       case _MAN_3D : dz = coordz_[v0] - coordz_[v1];
741         cost = (int) floor( dx+dy+dz+0.5 );
742         break;
743       case _CEIL_2D : cost = (int)ceil( sqrt( dx*dx + dy*dy ) + 0.5);
744         break;
745       case _ATT     : cost = (int)( sqrt( (dx*dx + dy*dy ) / 10 ) + 1);
746         break;
747      }
748   }
749   return(cost);
750}
751
752//#############################################################################
753
754void
755VrpModel::setModelData()
756{
757   CoinBigIndex i, j, numNonzeros = 0, numInt = 0;
758   int size;
759
760   int* varIndices = 0;
761   double* varValues = 0;
762
763   double* collb = new double [edgenum_];
764   double* colub = new double [edgenum_];
765   double* conUpper = new double[vertnum_];
766   double* conLower = new double[vertnum_];
767   double* objCoef = new double[edgenum_];
768
769   for (i = 0; i < edgenum_; i++){
770      objCoef[i] = edges_[i]->getObjCoef();
771   }
772
773   conLower[0] = conUpper[0] = 2*numroutes_;
774   for (i = 1; i < vertnum_; i++){
775      conUpper[i] = conLower[i] = 2.0;
776   }
777
778   char* colType = new char[edgenum_];
779
780   // Get number of integer and number of nonzero for memory allocation.
781   for (i = 0; i < edgenum_; ++i) {
782       numNonzeros += edges_[i]->getSize();
783       colType[i] = edges_[i]->getIntType();
784   }
785
786   int* indices = new int[numNonzeros];
787   double* values = new double[numNonzeros];
788   CoinBigIndex * start = new CoinBigIndex [edgenum_+1];
789
790   // Get collb, colub, obj, and matrix from variables
791   for (numInt = 0, numNonzeros = 0, i = 0; i < edgenum_; ++i) {
792      collb[i] = edges_[i]->getLbHard();
793      colub[i] = edges_[i]->getUbHard();
794
795      start[i] = numNonzeros;
796      varValues = edges_[i]->getValues();
797      varIndices = edges_[i]->getIndices();
798
799      size = edges_[i]->getSize();
800      for (j = 0; j < size; ++j, ++numNonzeros){
801         indices[numNonzeros] = varIndices[j];
802         values[numNonzeros] = varValues[j];
803      }
804   }
805   start[i] = numNonzeros;
806
807   // Set input data
808   setColMatrix(new CoinPackedMatrix(true, vertnum_, edgenum_, numNonzeros,
809                                     values, indices, start, 0));
810
811   setNumCons(vertnum_);
812   setNumVars(edgenum_);
813
814   setConLb(conLower);
815   setConUb(conUpper);
816
817   setVarLb(collb);
818   setVarUb(colub);
819
820   setObjCoef(objCoef);
821
822   setColType(colType);
823
824   delete [] start;
825   delete [] indices;
826   delete [] values;
827}
828
829//#############################################################################
830
832void
833VrpModel::readParameters(const int argnum, const char * const * arglist)
834{
836    int msgLevel = AlpsPar_->entry(AlpsParams::msgLevel);
837    if (msgLevel > 0) {
838        std::cout << "Reading in ALPS parameters ..." << std::endl;
839        std::cout << "Reading in BLIS parameters ..." << std::endl;
840        std::cout << "Reading in VRP parameters ..." << std::endl;
841    }
844}
845
846//#############################################################################
847
848/** Register knowledge. */
849void
850VrpModel::registerKnowledge() {
851    // Register model, solution, and tree node
852    assert(broker_);
853    broker_->registerClass(AlpsKnowledgeTypeModel, new VrpModel);
854    if (broker_->getMsgLevel() > 5) {
855        std::cout << "BLIS: Register Alps model." << std::endl;
856    }
857
858    broker_->registerClass(AlpsKnowledgeTypeNode, new BlisTreeNode(this));
859    if (broker_->getMsgLevel() > 5) {
860        std::cout << "BLIS: Register Alps node." << std::endl;
861    }
862
863    broker_->registerClass(AlpsKnowledgeTypeSolution, new VrpSolution);
864    if (broker_->getMsgLevel() > 5) {
865        std::cout << "BLIS: Register Alps solution." << std::endl;
866    }
867
868    broker_->registerClass(BcpsKnowledgeTypeConstraint, new BlisConstraint);
869    if (broker_->getMsgLevel() > 5) {
870        std::cout << "BLIS: Register Bcps constraint." << std::endl;
871    }
872
873    broker_->registerClass(BcpsKnowledgeTypeVariable, new VrpVariable);
874    if (broker_->getMsgLevel() > 5) {
875        std::cout << "BLIS: Register Bcps variable." << std::endl;
876    }
877}
878
879//#############################################################################
880
881AlpsReturnStatus
882VrpModel::encodeVrp(AlpsEncoded *encoded) const
883{
884    AlpsReturnStatus status = AlpsReturnStatusOk;
885
886    //encoded->writeRep(name_, 100); // No funtion
887    encoded->writeRep(vertnum_);
888    encoded->writeRep(edgenum_);
889    encoded->writeRep(numroutes_);
890    encoded->writeRep(depot_);
891    encoded->writeRep(capacity_);
892    encoded->writeRep(wtype_);
893
894    encoded->writeRep(demand_, vertnum_);
895
896    if (posx_) {
897        encoded->writeRep(posx_, vertnum_);
898    }
899    else {
900        encoded->writeRep(posx_, 0);
901    }
902    if (posy_) {
903        encoded->writeRep(posy_, vertnum_);
904    }
905    else {
906        encoded->writeRep(posy_, 0);
907    }
908
909    if (coordx_) {
910        encoded->writeRep(coordx_, vertnum_);
911    }
912    else {
913        encoded->writeRep(coordx_, 0);
914    }
915    if (coordy_) {
916        encoded->writeRep(coordy_, vertnum_);
917    }
918    else {
919        encoded->writeRep(coordy_, 0);
920    }
921    if (coordz_) {
922        encoded->writeRep(coordz_, vertnum_);
923    }
924    else {
925        encoded->writeRep(coordz_, 0);
926    }
927
928    encoded->writeRep(etol_);
929
930    VrpPar_->pack(*encoded);
931
932    return status;
933}
934
935//#############################################################################
936
937AlpsReturnStatus
938VrpModel::decodeVrp(AlpsEncoded &encoded)
939{
940    AlpsReturnStatus status = AlpsReturnStatusOk;
941    int k, tempInt = 0;
942
950
952    assert(tempInt == vertnum_);
953
955    assert(tempInt == vertnum_ || tempInt == 0);
957    assert(tempInt == vertnum_ || tempInt == 0);
958
960    assert(tempInt == vertnum_ || tempInt == 0);
962    assert(tempInt == vertnum_ || tempInt == 0);
964    assert(tempInt == vertnum_ || tempInt == 0);
965
967
968    VrpPar_->unpack(encoded);
969
970    // Allocate space for network for later use
971    n_ = new VrpNetwork(edgenum_, vertnum_);
972
973    // Insert variables to edges_.
974    for (k = 0; k < edgenum_; ++k) {
975        edges_.push_back(dynamic_cast<VrpVariable *>(variables_[k]));
976    }
977
978    // Add a cut generator
979    VrpCutGenerator *cg = new VrpCutGenerator(this, vertnum_);
980    cg->setStrategy(BlisCutStrategyPeriodic);
981    cg->setCutGenerationFreq(1); // Generate cuts at every node
983
984    return status;
985}
986
987//#############################################################################
988
989AlpsEncoded*
990VrpModel::encode() const
991{
992    AlpsReturnStatus status = AlpsReturnStatusOk;
993
994    // NOTE: "AlpsKnowledgeTypeModel" is the type name.
995    AlpsEncoded* encoded = new AlpsEncoded(AlpsKnowledgeTypeModel);
996
997    status = encodeAlps(encoded);
998    status = encodeBcps(encoded);
999    status = encodeBlis(encoded);
1000    status = encodeVrp(encoded);
1001
1002    return encoded;
1003}
1004
1005//#############################################################################
1006
1007void
1008VrpModel::decodeToSelf(AlpsEncoded& encoded)
1009{
1010    AlpsReturnStatus status = AlpsReturnStatusOk;
1011
1012    status = decodeAlps(encoded);
1013    status = decodeBcps(encoded);
1014    status = decodeBlis(encoded);
1015    status = decodeVrp(encoded);
1016}
1017
1018//#############################################################################
1019
Note: See TracBrowser for help on using the repository browser.