VRPH  1.0
src/VRP.cpp
Go to the documentation of this file.
00001 
00002 //                                                        //
00003 // This file is part of the VRPH software package for     //
00004 // generating solutions to vehicle routing problems.      //
00005 // VRPH was developed by Chris Groer (cgroer@gmail.com).  //
00006 //                                                        //
00007 // (c) Copyright 2010 Chris Groer.                        //
00008 // All Rights Reserved.  VRPH is licensed under the       //
00009 // Common Public License.  See LICENSE file for details.  //
00010 //                                                        //
00012 
00013 #include "VRPH.h"
00014 #define VRPH_MAX_CYCLES 500
00015 
00016 void VRPH_version()
00017 {
00018     printf("--------------------------------------------\n");
00019     printf("VRPH, version 1.0.1\nCopyright 2010 Chris Groer\nDistributed under Common Public License 1.0\n");
00020     printf("--------------------------------------------\n\n");
00021 }
00022 
00023 
00024 VRP::VRP(int n)
00025 {
00029 
00030     int i,j;
00031 
00032     num_nodes=n;
00033     num_original_nodes=n;
00034     total_demand=0;
00035     num_days=0;
00036     
00037     next_array = new int[n+2];
00038     pred_array = new int[n+2];
00039     route_num = new int[n+2];
00040     route = new VRPRoute[n+2];
00041     routed = new bool[n+2];
00042     best_sol_buff = new int[n+2];
00043     current_sol_buff = new int[n+2];
00044     search_space = new int[n+2];
00045     nodes = new VRPNode[n+2];    
00046     // Add an extra spot for the VRPH_DEPOT and for the dummy node
00047     
00048     symmetric=true;
00049     // Set to false only when we encounter FULL_MATRIX file
00050 
00051     forbid_tiny_moves=true;    
00052     // Default is to allow these moves
00053 
00054     d=NULL;
00055     // The distance matrix is allocated when the problem is loaded
00056     fixed=new bool*[n+2];
00057     fixed[0]=new bool[(n+2)*(n+2)];
00058     for(i=1;i<n+2;i++)
00059         fixed[i]=fixed[i-1]+(n+2);
00060     for(i=0;i<n+2;i++)
00061     {
00062         routed[i]=false;
00063         for(j=0;j<n+2;j++)
00064             fixed[i][j]=false;
00065     }
00066 
00067     // Set these to default values--they may change once
00068     // we read the file.
00069 
00070     min_vehicles=-1;
00071     has_service_times=false;
00072     max_route_length=VRP_INFINITY;
00073     orig_max_route_length=VRP_INFINITY;
00074     total_route_length=0.0;
00075     best_known=VRP_INFINITY;
00076     depot_normalized=false; 
00077     // Will be set to true if we shift nodes so VRPH_DEPOT is at origin
00078     
00079     // These are for record-to-record travel
00080     record = 0.0;
00081     deviation = VRPH_DEFAULT_DEVIATION;        
00082 
00083     // For keeping track of the statistics
00084 
00085     for(i=0;i<NUM_HEURISTICS;i++)
00086     {
00087         num_evaluations[i]=0;
00088         num_moves[i]=0;
00089 
00090     }
00091 
00092     total_service_time = 0.0;
00093 
00094     // Create the solution warehouse
00095     this->solution_wh=new VRPSolutionWarehouse(NUM_ELITE_SOLUTIONS,n);
00096 
00097     this->tabu_list=new VRPTabuList(MAX_VRPH_TABU_LIST_SIZE);
00098     this->route_wh=NULL;
00099 
00100     // Set this to true only if we have valid coordinates 
00101     // This is valid only for plotting the solution
00102     can_display=false;
00103     
00104     
00105 }
00106 
00107 VRP::VRP(int n, int ndays)
00108 {
00112 
00113     int i,j;
00114 
00115     num_nodes=n;
00116     num_original_nodes=n;
00117     total_demand=0;
00118     num_days=ndays;
00119     
00120     next_array = new int[n+2];
00121     pred_array = new int[n+2];
00122     route_num = new int[n+2];
00123     route = new VRPRoute[n+2];
00124     routed = new bool[n+2];
00125     best_sol_buff = new int[n+2];
00126     current_sol_buff = new int[n+2];
00127     search_space = new int[n+2];
00128     nodes = new VRPNode[n+2];    
00129     // Add an extra spot for the VRPH_DEPOT and for the dummy node
00130     
00131     forbid_tiny_moves=true;    
00132     // Default is to forbid these moves
00133 
00134     d=NULL;
00135     // The distance matrix is allocated when the problem is loaded
00136     fixed=new bool*[n+2];
00137     fixed[0]=new bool[(n+2)*(n+2)];
00138     for(i=1;i<n+2;i++)
00139         fixed[i]=fixed[i-1]+(n+2);
00140     for(i=0;i<n+2;i++)
00141     {
00142         routed[i]=false;
00143         for(j=0;j<n+2;j++)
00144             fixed[i][j]=false;
00145     }
00146 
00147     // Set these to default values--they may change once
00148     // we read the file.
00149 
00150     min_vehicles=-1;
00151     has_service_times=false;
00152     max_route_length=VRP_INFINITY;
00153     orig_max_route_length=VRP_INFINITY;
00154     total_route_length=0.0;
00155     best_known=VRP_INFINITY;
00156 
00157     depot_normalized=false; 
00158     // Will be set to true if we shift nodes so VRPH_DEPOT is at origin
00159     
00160     // These are for record-to-record travel
00161     record = 0.0;
00162     deviation = VRPH_DEFAULT_DEVIATION;        
00163 
00164     // For keeping track of the statistics
00165 
00166     for(i=0;i<NUM_HEURISTICS;i++)
00167     {
00168         num_evaluations[i]=0;
00169         num_moves[i]=0;
00170 
00171     }
00172 
00173     total_service_time = 0.0;
00174 
00175     // Create the solution warehouse
00176     this->solution_wh=new VRPSolutionWarehouse(NUM_ELITE_SOLUTIONS,n);
00177 
00178     this->tabu_list=new VRPTabuList(MAX_VRPH_TABU_LIST_SIZE);
00179     this->route_wh=NULL;
00180 
00181     // Now allocate d days worth of storage at each of the nodes
00182     for(i=0;i<=n+1;i++)
00183     {
00184         this->nodes[i].daily_demands=new int[ndays+1];
00185         this->nodes[i].daily_service_times=new double[ndays+1];
00186     }
00187 
00188     // Set this to true only if we have valid coordinates 
00189     // This is valid only for plotting the solution
00190     can_display=false;
00191     
00192 }
00193 
00194 VRP::~VRP()
00195 {
00199 
00200     delete [] this->best_sol_buff;
00201     delete [] this->current_sol_buff;
00202     delete [] this->d[0];
00203     delete [] this->d;
00204     delete [] this->fixed[0];
00205     delete [] this->fixed;
00206     delete [] this->next_array;
00207     delete [] this->search_space;
00208     delete [] this->nodes;
00209     delete [] this->pred_array;
00210     delete [] this->route;
00211     delete [] this->route_num;
00212     delete [] this->routed;
00213     delete this->solution_wh;
00214     delete this->tabu_list;
00215 
00216 }
00217 
00218 
00219 // Accessor functions for private data 
00220 int VRP::get_num_nodes()
00221 {
00225 
00226     return this->num_nodes;
00227 }
00228 
00229 double VRP::get_total_route_length()
00230 {
00234 
00235     return this->total_route_length;
00236 }
00237 
00238 double VRP::get_total_service_time()
00239 {
00243     
00244     return this->total_service_time;
00245 }
00246 
00247 double VRP::get_best_sol_buff(int *sol_buff)
00248 {
00254 
00255     memcpy(sol_buff,this->best_sol_buff,(this->num_nodes+1)*sizeof(int));
00256     return this->best_total_route_length;
00257 }
00258 
00259 double VRP::get_best_total_route_length()
00260 {
00264 
00265     return this->best_total_route_length;        
00266 }
00267 
00268 int VRP::get_total_number_of_routes()
00269 {
00273 
00274     return this->total_number_of_routes;
00275 }
00276 
00277 int VRP::get_num_original_nodes()
00278 {
00282 
00283     return this->num_original_nodes;
00284 }
00285 
00286 int VRP::get_num_days()
00287 {
00291 
00292     return this->num_days;
00293 }
00294 
00295 double VRP::get_best_known()
00296 {
00297     return this->best_known;
00298     
00299 }
00300 
00301 int VRP::get_max_veh_capacity()
00302 {   
00303     return this->max_veh_capacity;
00304 }
00305 
00306 double VRP::get_max_route_length()
00307 {
00308     return this->max_route_length;     
00309 }
00310 
00311 void VRP::set_best_total_route_length(double val)
00312 {
00313     this->best_total_route_length=val;
00314     return;
00315 
00316 }
00317 bool VRP::clone(VRP *W)
00318 { 
00322 
00323     this->balance_parameter=W->balance_parameter;
00324     this->best_known=W->best_known;
00325     this->best_total_route_length=W->best_total_route_length;
00326     // Copy the best solution buffer
00327     memcpy(this->best_sol_buff, W->best_sol_buff, (sizeof(int))*(W->num_nodes+2));
00328 
00329     // Copy the currnet solution buffer
00330     memcpy(this->current_sol_buff, W->current_sol_buff, (sizeof(int))*(W->num_nodes+2));
00331 
00332     this->coord_type=W->coord_type;
00333     this->d=W->d;                // OK to just copy the pointers here
00334     this->deviation=W->deviation;
00335     this->display_type=W->display_type;
00336     this->edge_weight_format=W->edge_weight_format;
00337     this->dummy_index=W->dummy_index;
00338     this->has_service_times=W->has_service_times;
00339     this->matrix_size=W->matrix_size;
00340     this->max_route_length=W->max_route_length;
00341     this->neighbor_list_size=W->neighbor_list_size;
00342     this->nodes = W->nodes; // OK to just copy the pointers here?
00343 
00344     this->num_nodes=W->num_nodes;
00345     this->total_route_length=W->total_route_length;
00346     this->orig_max_route_length=W->orig_max_route_length;
00347     this->orig_max_veh_capacity=W->orig_max_veh_capacity;
00348 
00349     this->problem_type=W->problem_type;
00350     this->record=W->record;
00351 
00352 
00353     this->search_size=W->search_size;
00354 
00355     this->temperature=W->temperature;
00356     this->total_number_of_routes=W->total_number_of_routes;
00357     this->total_service_time=W->total_service_time;
00358     this->max_veh_capacity=W->max_veh_capacity;
00359     this->violation=W->violation;
00360     this->edge_weight_type=W->edge_weight_type;
00361 
00362     // Assume that the current_solution has been sent to W->
00363     W->export_solution_buff(W->current_sol_buff);
00364     this->import_solution_buff(W->current_sol_buff);
00365     return true;    
00366 
00367 }
00368 
00369 
00370 void VRP::create_distance_matrix(int type)
00371 {    
00378     
00379 
00380     int i,j,k,n;
00381     n=this->num_nodes;
00382 
00383     if(type==VRPH_EXPLICIT)
00384     {
00385         // We have presumably already loaded in the distance matrix
00386         // But still add in the service times!!
00387         if(this->has_service_times)
00388         {
00389             for(i=0;i<=n+1;i++)
00390             {
00391                 for(j=0;j<=n+1;j++)
00392                 {
00393                     //printf("Adding service time of %f\n",this->nodes[j].service_time);
00394                     this->d[i][j]+= .5*this->nodes[j].service_time;
00395                     this->d[i][j]+= .5*this->nodes[i].service_time;
00396                 }
00397             }
00398         }
00399         // We have presumably already loaded in the distance matrix!
00400         return;
00401     }
00402 
00403     n= this->num_original_nodes;
00404     k=0;
00405 
00406     // Otherwise construct the matrix - we store the whole thing even though
00407     // it is symmetric as we found that this was quite a bit faster...
00408 
00409     for(i=0;i<=n+1;i++)
00410     {
00411         for(j=0;j<=n+1;j++)
00412             // CSG: Changed Nov. 20 - add in .5 service time from i and j
00413             this->d[i][j]=VRPDistance(type, this->nodes[i].x,this->nodes[i].y,this->nodes[j].x,
00414                 this->nodes[j].y) + 
00415                 0.5*this->nodes[i].service_time +
00416                 0.5*this->nodes[j].service_time; 
00417     }
00418 
00419 
00420     return;
00421 
00422 }
00423 
00424 void VRP::create_neighbor_lists(int nsize)
00425 {
00430 
00431     if(nsize>num_nodes )
00432     {
00433         fprintf(stderr,"Requested neighbor list size is greater than num_nodes!\n%d>%d\n",
00434             nsize,num_nodes);
00435         report_error("%s: Neighbor list error!!\n",__FUNCTION__);
00436     }
00437 
00438     if(nsize>MAX_NEIGHBORLIST_SIZE )
00439     {
00440         fprintf(stderr,"Requested neighbor list size is greater than MAX_NEIGHBORLIST_SIZE!\n%d>%d\n",
00441             nsize,MAX_NEIGHBORLIST_SIZE);
00442         report_error("%s: Neighbor list error!!\n",__FUNCTION__);
00443     }
00444 
00445 
00446     int i,j,n,b;
00447     int k;
00448     double dd;
00449     double max;
00450     int maxpos;
00451     VRPNeighborElement *NList;
00452 
00453     n= num_nodes;
00454 
00455     // Set the neighbor_list_size value
00456     neighbor_list_size=nsize;
00457 
00458     NList=new VRPNeighborElement[nsize];
00459 
00460     // First, do the neighbor_list for the VRPH_DEPOT
00461     max=0.0;
00462     // NEW
00463     maxpos=0;
00464 
00465     for(i=1;i<=nsize;i++)
00466     {
00467         NList[i-1].val=d[VRPH_DEPOT][i];
00468 
00469         NList[i-1].position = i;
00470 
00471         // keep track of the largest value
00472         if(NList[i-1].val > max)
00473         {
00474             max=NList[i-1].val;
00475             maxpos=i-1;
00476         }
00477 
00478     }    
00479 
00480     // Now go through the rest of the nodes.
00481     for(i=nsize+1;i<=n;i++)
00482     {
00483         if(d[VRPH_DEPOT][i]<max)
00484         {
00485             // Replace element in position maxpos
00486             NList[maxpos].val=d[VRPH_DEPOT][i];
00487             //nodes[i].depot_distance;
00488             NList[maxpos].position=i;
00489         }
00490         // Now find the new max
00491         max=0.0;
00492         for(b=0;b<nsize;b++)
00493         {
00494             if(NList[b].val>max)
00495             {
00496                 maxpos=b;
00497                 max=NList[b].val;
00498             }
00499         }
00500     }
00501 
00502     qsort(NList,nsize,sizeof(VRPNeighborElement), VRPNeighborCompare);
00503     i=0;
00504     for(b=0;b<nsize;b++)
00505     {
00506         nodes[VRPH_DEPOT].neighbor_list[b].val=NList[b].val;
00507         nodes[VRPH_DEPOT].neighbor_list[b].position=NList[b].position;
00508 
00509 #        if NEIGHBOR_DEBUG
00510         printf("(%d,%d,%f) \n",VRPH_DEPOT,nodes[VRPH_DEPOT].neighbor_list[b].position,
00511             nodes[VRPH_DEPOT].neighbor_list[b].val);
00512 #endif
00513 
00514 
00515     }
00516 
00517 #    if NEIGHBOR_DEBUG
00518     printf("\n\n");
00519 #endif
00520 
00521 
00522     k=0;
00523     // Done with NeighborList for VRPH_DEPOT.  Now do it for the rest of the nodes
00524     for(i=1;i<=n;i++)
00525     {
00526         // Loop through all the nodes and create their Neighbor Lists
00527         // First initialize the NList
00528         for(j=0;j<nsize;j++)
00529         {
00530             NList[j].position=VRP_INFINITY;
00531             NList[j].val=VRP_INFINITY;
00532         }
00533 
00534 
00535         // VRPH_DEPOT can be in NList, so set NList[0] to depot_distance
00536         NList[0].position=VRPH_DEPOT;
00537         NList[0].val=d[i][VRPH_DEPOT];
00538 
00539 
00540         for(j=1;j<i;j++)
00541         {
00542 #if NEIGHBOR_DEBUG
00543             printf("Loop 1; i=%d; in j loop(%d)\n",i,j);
00544 #endif
00545 
00546             dd=d[i][j];
00547 
00548             if(j<nsize)
00549             {
00550                 // Construct the original nsize long neighbor list
00551                 NList[j].val=dd;
00552                 NList[j].position=j;
00553                 // Update max size
00554                 if(NList[j].val>max)
00555                 {
00556                     max=NList[j].val;
00557                     maxpos=j;
00558 
00559                 }
00560             }
00561             else
00562             {
00563                 // j >= nsize
00564                 if(dd<max)
00565                 {
00566                     NList[maxpos].val=dd;
00567                     NList[maxpos].position=j;
00568                     // Now find the new max element;
00569                     max=0.0;
00570                     for(b=0;b<nsize;b++)
00571                     {
00572                         if(NList[b].val>max)
00573                         {
00574                             max=NList[b].val;
00575                             maxpos=b;
00576                         }
00577                     }
00578                 }
00579             }
00580 
00581         }
00582 
00583 
00584         for(j=i+1;j<=n;j++)
00585         {
00586 #if NEIGHBOR_DEBUG
00587             printf("Loop 2; i=%d; in j loop(%d)\n",i,j);
00588 #endif
00589 
00590             dd=d[i][j];
00591 
00592             if(j<=nsize)// was <
00593             {
00594                 // Construct the original nsize long neighbor list
00595                 NList[j-1].val=dd;
00596                 NList[j-1].position=j;
00597                 // Update max size
00598                 if(NList[j-1].val>max)
00599                 {
00600                     max=NList[j-1].val;
00601                     maxpos=j-1;
00602 
00603                 }
00604             }
00605             else
00606             {
00607                 if(dd<max)
00608                 {
00609                     // Replace maxpos with this one
00610                     NList[maxpos].val=dd;
00611                     NList[maxpos].position=j;
00612                     // Now find the new max element;
00613                     max=0.0;
00614                     for(b=0;b<nsize;b++)
00615                     {
00616                         if(NList[b].val>max)
00617                         {
00618                             max=NList[b].val;
00619                             maxpos=b;
00620                         }
00621                     }
00622                 }
00623             }
00624 
00625         }
00626 
00627 
00628 #if NEIGHBOR_DEBUG
00629         printf("Node %d neighbor list before sorting\n",i);
00630         for(b=0;b<nsize;b++)
00631         {
00632             printf("NList[%d]=%d,%f[%f]\n",b,NList[b].position,NList[b].val,d[i][NList[b].position]);
00633         }
00634 #endif
00635         // Now sort the NList and stick the resulting list into
00636         // the neighbor list for node i
00637         qsort (NList, nsize, sizeof(VRPNeighborElement), VRPNeighborCompare);
00638 
00639         for(b=0;b<nsize;b++)
00640         {
00641             nodes[i].neighbor_list[b].position=NList[b].position;
00642             nodes[i].neighbor_list[b].val=NList[b].val;
00643             if(i== nodes[i].neighbor_list[b].position)
00644             {
00645                 fprintf(stderr,"ERROR:: Node %d is in it's own neighbor list!!\n", i);
00646                 fprintf(stderr,"i=%d; %d[%d],%f[%f])\n",i,nodes[i].neighbor_list[b].position,
00647                     NList[b].position,nodes[i].neighbor_list[b].val,NList[b].val);
00648 
00649                 report_error("%s: Error creating neighbor lists\n",__FUNCTION__);
00650             }
00651 #if NEIGHBOR_DEBUG
00652 
00653             printf("NList[%d]=%d,%f[%f]",b,nodes[i].neighbor_list[b].position,nodes[i].neighbor_list[b].val,d[i][nodes[i].neighbor_list[b].position]);
00654 
00655 #endif
00656         }
00657 
00658         // Hackk to make the VRPH_DEPOT in everyone's list
00659         //nodes[i].neighbor_list[nsize-1].position=VRPH_DEPOT;
00660         //nodes[i].neighbor_list[nsize-1].val=d[i][VRPH_DEPOT];
00661 
00662     }
00663 
00664     delete [] NList;
00665 
00666     return;
00667 
00668 }
00669 
00670 
00671 bool VRP::check_feasibility(VRPViolation *VV)
00672 {
00678 
00679     int i;
00680     bool is_feasible=true; 
00681     // We will set this to false if we find a violation
00682 
00683     VV->capacity_violation=-VRP_INFINITY;
00684     VV->length_violation=-VRP_INFINITY;
00685 
00686     this->normalize_route_numbers();
00687 
00688     for(i=1;i<=this->total_number_of_routes;i++)
00689     {
00690         if(this->route[i].length>this->orig_max_route_length)
00691         {
00692             // Length violation
00693             if(this->route[i].length - this->orig_max_route_length > VV->length_violation)
00694                 VV->length_violation = this->route[i].length - this->orig_max_route_length ;
00695 
00696             is_feasible=false;
00697         }
00698 
00699         if(this->route[i].load>this->orig_max_veh_capacity)
00700         {
00701             // Load violation
00702             if(this->route[i].load - this->orig_max_veh_capacity > VV->capacity_violation)
00703                 VV->capacity_violation = this->route[i].load - this->orig_max_veh_capacity ;
00704 
00705             is_feasible=false;
00706 
00707         }
00708 
00709     }
00710 
00711     return is_feasible;
00712 
00713 }    
00714 
00715 
00716 
00717 
00718 void VRP::refresh_routes()
00719 {
00724 
00725 
00726     int i, n, cnt;
00727     double len=0;
00728     double rlen=0;
00729     double tot_len=0;
00730     int next_node;
00731     int current_node, current_route, route_start, current_start, current_end;
00732     int total_load = 0;
00733     int current_load = 0;
00734 
00735 
00736     n=num_nodes;
00737 
00738     i = 1;
00739     cnt = 0;
00740     route_start = -next_array[VRPH_DEPOT];
00741 
00742     current_node = route_start;
00743     current_route = route_num[current_node];
00744 
00745     current_start = route[current_route].start;
00746     current_end = route[current_route].end;
00747 
00748     total_load+=nodes[current_node].demand;
00749     current_load+=nodes[current_node].demand;
00750     len+=d[VRPH_DEPOT][current_node];
00751     rlen+=d[VRPH_DEPOT][current_node];
00752 
00753     while(route_start != 0 && i<num_nodes+1)
00754     {
00755         // When we finally get a route beginning at 0, this is the last route 
00756         // and there is no next route, so break out
00757         if(next_array[current_node]==current_node)
00758         {
00759             // We've entered a loop
00760             fprintf(stderr,"(2)Self loop found in next array(%d)\n",current_node);
00761             report_error("%s: Error in refresh_routes()\n",__FUNCTION__);
00762         }
00763 
00764         if(next_array[current_node]==0)
00765         {
00766             // We're at the end of the Solution!
00767 
00768             len+=d[current_node][VRPH_DEPOT];
00769             rlen+=d[current_node][VRPH_DEPOT];
00770             current_route=route_num[current_node];
00771             route[current_route].length = rlen;
00772             route[current_route].load = current_load;
00773             total_route_length=len;
00774 
00775             // We're done now
00776             return;
00777         }
00778 
00779         if(next_array[current_node]>0)
00780         {
00781             // Next node is somewhere in the middle of a route        
00782 
00783             next_node = next_array[current_node];
00784             len+=d[current_node][next_node];
00785             rlen+=d[current_node][next_node];
00786             current_node=next_node;
00787             total_load+=nodes[current_node].demand;
00788             current_load+=nodes[current_node].demand;
00789             cnt++;
00790         }
00791         else
00792         {
00793             // We must have a non-positive "next" node indicating the beginning of a new route
00794             //len+=nodes[current_node].depot_distance;
00795             len+=d[current_node][VRPH_DEPOT];
00796             //rlen+=nodes[current_node].depot_distance;
00797             rlen+=d[current_node][VRPH_DEPOT];
00798             tot_len+=rlen;
00799             current_route=route_num[current_node];
00800             current_end = route[current_route].end;
00801 
00802             route[current_route].length = rlen;
00803             route[current_route].load = current_load;
00804 
00805             i++;
00806             route_start = - (next_array[current_node]);    
00807             current_route = route_num[route_start];
00808             current_start = route[current_route].start;
00809             current_end = route[current_route].end;
00810             if(route_start != current_start)
00811             {
00812                 fprintf(stderr,"Route %d:  %d != %d\n",current_route, route_start, current_start);
00813                 report_error("%s: Error in refresh_routes()\n",__FUNCTION__);
00814             }
00815 
00816             current_node = route_start;
00817             total_load+=nodes[current_node].demand;
00818             // reset current_load to 0
00819             current_load=nodes[current_node].demand;
00820             len+=d[VRPH_DEPOT][current_node];
00821             rlen=d[VRPH_DEPOT][current_node];
00822             cnt++;
00823         }
00824     }
00825 
00826 
00827     return;
00828 
00829 
00830 }
00831 void VRP::create_pred_array()
00832 {
00836     int i,j;
00837 
00838     i=VRPH_DEPOT;
00839     j=next_array[i];
00840     while(j!=VRPH_DEPOT)
00841     {
00842         if(j>0)
00843             pred_array[j]=i;
00844         else
00845             pred_array[-j]=-i;
00846         
00847         i=VRPH_ABS(j);
00848         j=next_array[i];
00849 
00850     }
00851     // The VRPH_DEPOT
00852     pred_array[j]=-i;
00853     return;
00854 
00855     
00856 
00857 
00858 
00859     /*for(i=0;i<=num_nodes;i++)
00860     {
00861         j=next_array[i];
00862         if(j>0)
00863             pred_array[j] = i;
00864         else
00865             pred_array[-j] = -i;
00866     }
00867     return;*/
00868 }
00869 
00870 
00871 
00872 
00873 
00874 bool VRP::get_segment_info(int a, int b, struct VRPSegment *S)
00875 {
00885 
00886 
00887     if(a==b)
00888     {
00889         // Degenerate case...
00890         S->segment_start=b;
00891         S->segment_end=a;
00892         S->len=0;  //nodes[a].service_time??;
00893         S->load=nodes[a].demand;
00894         S->num_custs=1;
00895         return true;
00896 
00897     }
00898 
00899     S->len=0;    
00900     S->segment_start=a;
00901     S->segment_end=b;
00902     S->num_custs=0;        
00903 
00904 
00905     // Check special cases
00906     if(a==VRPH_DEPOT)
00907     {
00908         // Safe to assume we want the segment 0-i-j-...-k-b-...
00909 
00910         S->segment_start=route[route_num[b]].start;
00911         S->segment_end=b;
00912         S->len+=d[VRPH_DEPOT][S->segment_start];
00913 
00914     }
00915 
00916     if(b==VRPH_DEPOT)
00917     {
00918         // Safe to assume we want the segment ...a-i-j-...-k-0
00919 
00920         S->segment_start=a;
00921         S->segment_end=route[route_num[a]].end;
00922 
00923     }
00924 
00925     // Now calculate the length, load, and # customer on this segment
00926 
00927     int current_node = S->segment_start;
00928     int next_node;
00929 
00930     S->load=nodes[current_node].demand;
00931     if(current_node!=dummy_index)
00932         S->num_custs++;;
00933 
00934     while(current_node != S->segment_end)
00935     {
00936         next_node = VRPH_MAX(next_array[current_node],0);
00937         S->len+=d[current_node][next_node];
00938         current_node=next_node;        
00939         S->load+=nodes[next_node].demand;
00940         if(current_node!= dummy_index)
00941             S->num_custs += 1;
00942     }
00943 
00944     if(b==VRPH_DEPOT)
00945         S->len+=d[S->segment_end][VRPH_DEPOT];
00946 
00947 
00948 
00949     return true;
00950 
00951 }
00952 
00953 int VRP::get_string_end(int a, int len)
00954 {
00960 
00961     int ctr=1;
00962     int current_node;
00963 
00964     current_node=a;
00965     while(ctr<len)
00966     {
00967         current_node= next_array[current_node];
00968         if(current_node<0)
00969         {
00970             // The string of length len beginning at a is not contained in a single route
00971             return -1;
00972         }
00973         ctr++;
00974 
00975     }
00976 
00977     return current_node;
00978 }
00979 
00980 void VRP::reverse_route(int i)
00981 {
00987 
00988     
00989     int current_node, start_node, last_node, prev_route, next_route, temp;
00990     int orig_start, orig_end ;
00991 
00992     if(i<=0)
00993     {
00994         // Illegal!
00995         fprintf(stderr,"Reversing route  of negative index?? i=%d\n",i);
00996         report_error("%s: Can't reverse route!\n",__FUNCTION__);
00997     }
00998 
00999     orig_end = route[i].end;
01000     orig_start = route[i].start;
01001 
01002     start_node=orig_start;
01003     current_node=start_node;
01004 
01005     // We begin with the first node in the route - we will change its next_array[]
01006     // value at the very end
01007 
01008     temp = next_array[current_node];
01009     prev_route = -pred_array[current_node];
01010     // So next[prev_route] should equal start_node
01011     // We will change this so that next[prev_route]=-last_node
01012 
01013     pred_array[current_node] = temp;
01014     current_node = temp;
01015     while( (temp = next_array[current_node]) >0)
01016     {        
01017         next_array[current_node] = pred_array[current_node];
01018         pred_array[current_node] = temp;
01019         current_node = temp;
01020     }
01021 
01022     // current_node is now the end of the original route, so we can finish the reversal
01023 
01024     temp= pred_array[current_node];
01025     last_node = current_node;
01026     next_route = -next_array[last_node];
01027 
01028     route[i].end=orig_start;
01029     route[i].start=orig_end;
01030 
01031 #if REVERSE_DEBUG
01032     printf("start_node: %d; last_node: %d; prev_route: %d; next_route: %d\n",start_node,last_node, prev_route,next_route);
01033 #endif
01034 
01035     // Final set of updates
01036     next_array[prev_route]=-last_node;
01037     pred_array[next_route]=-start_node;
01038     next_array[start_node]=-next_route;
01039     pred_array[last_node]= -prev_route;
01040     next_array[last_node]=temp;
01041     next_array[prev_route]=-last_node;
01042 
01043     // Need to update length if asymmetric
01044     if(!this->symmetric)
01045     {
01046         // Update length as it possibly changed!
01047         this->refresh_routes();
01048     }
01049     
01050     return;
01051 }
01052 
01053 
01054 
01055 bool VRP::postsert_dummy(int i)
01056 {
01057     // This function inserts dummy node after i
01058 
01059     int pre_i, post_i, dummy;
01060     int start, end, start_i, end_i;
01061     int n, i_route;
01062     //double tu, uv, tv, iu, ju, ij, u_loss, i_gain, i_change, i_length, u_length;
01063 
01064     // Special cases first
01065 
01066     if(i<=VRPH_DEPOT || i>matrix_size)
01067         report_error("%s: input doesn't make sense\n",__FUNCTION__);
01068 
01069     n= num_nodes;
01070     dummy= dummy_index;
01071 
01072     i_route= route_num[i];
01073 
01074     start_i= route[i_route].start;
01075     end_i= route[i_route].end;
01076 
01077     start=start_i;
01078     if(end_i==i)
01079         // i is last in the route
01080         end=dummy;
01081     else 
01082         end=end_i;
01083 
01084     // pre_i is what used to be before i
01085     pre_i= pred_array[i];
01086     // post_i is what used to be after i
01087     post_i= next_array[i];
01088 
01089     next_array[i] = dummy;
01090     next_array[dummy] = post_i;
01091 
01092     // Now need to update pred_array as well!
01093 
01094     // u's predecessor is now i
01095     pred_array[dummy]=i;
01096 
01097     // The element who used to be after i is now preceded by u 
01098     if(post_i>=0)
01099         pred_array[post_i]=dummy;
01100     else
01101         // post_i 
01102         pred_array[-post_i]=-dummy;
01103 
01104     //start_array[dummy]=start;
01105     //end_array[dummy]=end;
01106 
01107     // Update i_route information
01108     route_num[dummy]=i_route;
01109     route[i_route].end=end;
01110     route[i_route].start=start;
01111 
01112     return true;
01113 }
01114 
01115 bool VRP::presert_dummy(int i)
01116 {
01117     // This function inserts a dummy node BEFORE node i in 
01118     // whatever route node i is currently in
01119 
01120     int pre_i, post_i;
01121     int start, end, start_i, end_i, dummy;
01122     int n, i_route;
01123 
01124 
01125     n= num_nodes ;
01126     dummy = dummy_index; // dummy node has index (n+1)
01127 
01128     if(i<=VRPH_DEPOT)
01129         report_error("%s: bad index\n",__FUNCTION__);
01130 
01131     i_route = route_num[i];
01132 
01133     start_i= route[i_route].start;
01134     end_i= route[i_route].end;
01135 
01136 
01137     // Get the start and end of i's route
01138     start=start_i;
01139     end=end_i;
01140 
01141     if(start==i)
01142         // i is first in the route, so dummy is now the first node
01143         start=dummy;
01144 
01145     // pre_i is what used to be before i
01146     pre_i= pred_array[i];
01147     // post_i is what used to be after i
01148     post_i= next_array[i];
01149 
01150 
01151     // dummy is now followed by i
01152     next_array[dummy]=i;
01153     pred_array[i]=dummy;
01154     pred_array[dummy]=pre_i;    
01155 
01156     // The element who used to be after i is now preceded by dummy
01157     if(pre_i>0)// was >=!!
01158         next_array[pre_i]=dummy;
01159     else
01160         // post_i 
01161         next_array[VRPH_ABS(pre_i)]=-dummy;
01162 
01163     // Update i_route information
01164     route_num[dummy]=i_route;
01165     route[i_route].end=end;
01166     route[i_route].start=start;
01167 
01168     // Add in the relevant Data fields for the dummy node!!
01169 
01170     return true;
01171 
01172 }
01173 bool VRP::remove_dummy()
01174 {
01175     int pre_d, post_d, d_route, dummy, d_start, d_end, n;
01176 
01177     n= num_nodes;
01178     dummy= dummy_index;
01179 
01180     post_d= next_array[dummy];
01181     pre_d= pred_array[dummy];
01182 
01183     if(post_d > dummy || post_d < -dummy || pre_d > dummy || pre_d < -dummy)
01184     {
01185         fprintf(stderr,"post_d= %d; pre_d=%d\n",post_d,pre_d);
01186         report_error("%s: invalid indices\n",__FUNCTION__);
01187     }
01188 
01189     d_route= route_num[dummy];
01190     d_start= route[d_route].start;
01191     d_end= route[d_route].end;
01192 
01193     if(d_start==dummy)
01194     {
01195         if(post_d<0)
01196         {
01197             report_error("%s: post_d error in removal\n",__FUNCTION__);
01198         }
01199 
01200         route[d_route].start=post_d;
01201     }
01202 
01203     if(d_end==dummy)
01204     {
01205         if(pre_d<0)
01206         {
01207             report_error("%s: pre_d error in removal\n",__FUNCTION__);
01208         }
01209         route[d_route].end=pre_d;
01210     }
01211 
01212     next_array[VRPH_ABS(pre_d)]=post_d;
01213 
01214     if(d_start==dummy)
01215         next_array[VRPH_ABS(pre_d)]=-post_d;
01216 
01217 
01218     pred_array[VRPH_ABS(post_d)]=pre_d;
01219 
01220     if(d_end==dummy)
01221         pred_array[VRPH_ABS(post_d)]=-pre_d;
01222 
01223 
01224     return true;
01225 }
01226 bool VRP::create_default_routes()
01227 {
01234 
01235 
01236     int i,n;
01237     bool is_feasible=true;
01238     // No violations yet...
01239     violation.capacity_violation = 0;
01240     violation.length_violation = 0;
01241 
01242 
01243     total_route_length=0;
01244 
01245     n = this->num_original_nodes;
01246     // First, create the initial routes:  VRPH_DEPOT - node - VRPH_DEPOT for all nodes
01247 
01248     routed[VRPH_DEPOT]=true;
01249 
01250     next_array[VRPH_DEPOT] = -1;
01251     for(i=1;i<=n;i++)
01252     {
01253         next_array[i] = -(i+1);
01254         total_route_length+= (d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]);
01255 
01256         route_num[i]=i;
01257         route[i].start=i;
01258         route[i].end=i;
01259         route[i].load= nodes[i].demand;
01260         route[i].length= d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT];
01261 
01262         
01263         // Check capacities
01264         if(route[i].load > max_veh_capacity)
01265             is_feasible=false;
01266         if(route[i].length > max_route_length)
01267             is_feasible=false;
01268 
01269         route[i].num_customers=1;
01270 
01271         routed[i]=true;    
01272 
01273     }
01274 
01275 
01276     next_array[n]=VRPH_DEPOT;
01277 
01278     route_num[VRPH_DEPOT]=0;
01279     
01280     // Now create the associated pred_array implied by the newly created next_array
01281     create_pred_array();
01282 
01283     total_number_of_routes=n;
01284 
01285 
01286 
01287 
01288 
01289     if(is_feasible == false)
01290     { 
01291         // We had infeasible routes - record the violation in the Solution and 
01292         // return false;
01293 
01294 
01295         for(i=1;i<=n;i++)
01296         {
01297             routed[i]=false;    
01298             // Set routed status to false since default routes are infeasible
01299 
01300             // Check capacities
01301             if(route[i].load > max_veh_capacity)
01302             {
01303                 // Update the worst violations
01304 
01305                 printf("Default routes load violation: %d > %d\n",route[i].load,
01306                     max_veh_capacity);
01307 
01308                 if( (route[i].load - max_veh_capacity) > 
01309                     violation.capacity_violation)
01310                 {
01311                     violation.capacity_violation = 
01312                         (route[i].load - max_veh_capacity);
01313                 }
01314 
01315             }
01316             if(route[i].length > max_route_length)
01317             {
01318                 // Update the worst violations
01319 
01320                 printf("Default routes length violation: %f > %f\n",route[i].length,
01321                     max_route_length);
01322 
01323                 if( (route[i].length -max_route_length ) >
01324                     violation.length_violation)
01325                 {
01326                     violation.length_violation = 
01327                         (route[i].length -max_route_length );
01328                 }
01329 
01330 
01331             }
01332 
01333 
01334         }
01335 
01336         return false;
01337     }
01338     else
01339     {
01340         // All routes were feasible
01341         for(i=1;i<=n;i++)
01342         {
01343             routed[i]=true;
01344         }
01345         return true;
01346     }
01347 
01348 
01349 }
01350 
01351 
01352 bool VRP::create_default_routes(int day)
01353 {
01360 
01361     int i,n;
01362     bool is_feasible=true;
01363     // No violations yet...
01364     violation.capacity_violation = 0;
01365     violation.length_violation = 0;
01366     total_route_length=0;
01367 
01368     n = this->num_original_nodes;
01369     this->num_nodes=n;
01370 
01371     // First, create the initial routes:  VRPH_DEPOT - node - VRPH_DEPOT for all nodes
01372 
01373     routed[VRPH_DEPOT]=true;
01374     next_array[VRPH_DEPOT] = -1;
01375     for(i=1;i<=n;i++)
01376     {
01377 
01378         next_array[i] = -(i+1);
01379         total_route_length+= (d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT]);
01380 
01381         route_num[i]=i;
01382         route[i].start=i;
01383         route[i].end=i;
01384         route[i].load= nodes[i].demand;
01385         route[i].length= d[VRPH_DEPOT][i] + d[i][VRPH_DEPOT];
01386         route[i].num_customers=1;
01387         routed[i]=true;    
01388 
01389     }
01390 
01391 
01392     next_array[n]=VRPH_DEPOT;
01393 
01394     route_num[VRPH_DEPOT]=0;
01395     
01396     // Now create the associated pred_array implied by the newly created next_array
01397     create_pred_array();
01398 
01399     total_number_of_routes=n;
01400 
01401     // Now eject the nodes that don't require service on this day
01402     for(i=1;i<=n;i++)
01403     {
01404         if(this->nodes[i].daily_demands[day]==-1)// indicates no service required
01405             this->eject_node(i);
01406     }
01407 
01408     // Now check for feasibility;
01409     this->normalize_route_numbers();
01410 
01411     // Check capacities
01412     for(i=1;i<=this->total_number_of_routes;i++)
01413     {
01414         if(route[i].load > max_veh_capacity)
01415             is_feasible=false;
01416         if(route[i].length > max_route_length)
01417             is_feasible=false;
01418     }
01419 
01420     if(is_feasible == false)
01421     { 
01422         // We had infeasible routes - record the violation in the Solution and 
01423         // return false;
01424 
01425 
01426         for(i=1;i<=n;i++)
01427         {
01428             routed[i]=false;    
01429             // Set routed status to false since default routes are infeasible
01430 
01431             // Check capacities
01432             if(route[i].load > max_veh_capacity)
01433             {
01434                 // Update the worst violations
01435 
01436                 printf("Default routes load violation: %d > %d\n",route[i].load,
01437                     max_veh_capacity);
01438 
01439                 if( (route[i].load - max_veh_capacity) > 
01440                     violation.capacity_violation)
01441                 {
01442                     violation.capacity_violation = 
01443                         (route[i].load - max_veh_capacity);
01444                 }
01445 
01446             }
01447             if(route[i].length > max_route_length)
01448             {
01449                 // Update the worst violations
01450 
01451                 printf("Default routes length violation: %f > %f\n",route[i].length,
01452                     max_route_length);
01453 
01454                 if( (route[i].length -max_route_length ) >
01455                     violation.length_violation)
01456                 {
01457                     violation.length_violation = 
01458                         (route[i].length -max_route_length );
01459                 }
01460             }
01461         }
01462         return false;
01463     }
01464     else
01465         return true;
01466     
01467 }
01468 
01469 
01470 int VRP::count_num_routes()
01471 {
01475 
01476     int current,next;
01477     int num=0;
01478 
01479     current=VRPH_DEPOT;
01480     next=-1;
01481 
01482     while(next!=VRPH_DEPOT)
01483     {
01484         next= next_array[current];
01485         if(next<0)
01486         {
01487             // We're in a new route
01488             num++;
01489             current=-next;
01490         }
01491         else
01492             current=next;
01493     }
01494 
01495     return num;
01496 
01497 
01498 }
01499 
01500 bool VRP::perturb()
01501 {
01505 
01506     int i, n, pre, post;
01507     int current,next;
01508     n= num_nodes;
01509     VRPNeighborElement *v;
01510     VRPMove M;
01511 
01512     Postsert Postsert;
01513     Presert  Presert;
01514 
01515     v=new VRPNeighborElement[n];
01516 
01517     current=VRPH_ABS(next_array[VRPH_DEPOT]);
01518     i=0;
01519     while(current!=VRPH_DEPOT)
01520     {
01521         pre= VRPH_MAX(VRPH_DEPOT,pred_array[current]);
01522         post= VRPH_MAX(VRPH_DEPOT, next_array[current]);
01523 
01524         // Define s[i]= d(pre(i),i) + d(i,next(i))-d(pre(i),next(i))
01525         v[i].val=((double) nodes[current].demand) / 
01526             (VRPH_EPSILON+d[pre][current]+d[current][post]-d[pre][post]);
01527         v[i].position=current;
01528         i++;
01529         next=VRPH_ABS(next_array[current]);
01530         current=next;
01531     }
01532 
01533 
01534     // Now sort the v[i]'s in ascending order of the val's
01535 
01536     qsort(v,n,sizeof(VRPNeighborElement), VRPNeighborCompare);
01537 
01538     int m=(int) (VRPH_MIN(30,n/5));
01539     
01540 
01541     // Now insert the first m nodes from the sorted list into new locations.
01542 
01543     int a,b,c,j,k,node1=0, node2=0,b_route,b_load,k_demand, jj, ll;
01544     double  best_savings,b_len,jk,kl,jl;
01545 
01546 
01547     for(j=0;j<m;j++)
01548     {
01549         best_savings=VRP_INFINITY;
01550 
01551         k=v[j].position;
01552         k_demand= nodes[k].demand;
01553 
01554         // Node k will be moved
01555 
01556         // Former neighboring nodes of k
01557         jj=VRPH_MAX(pred_array[k],0);
01558         ll=VRPH_MAX(next_array[k],0);
01559 
01560         jk= d[jj][k];
01561         kl= d[k][ll];
01562         jl= d[jj][ll];
01563 
01564 
01565         // Node k will be moved - look for a new location
01566         // OLD SITUATION:  jj-k-ll
01567         b=VRPH_ABS(next_array[VRPH_DEPOT]);
01568         while(b!=VRPH_DEPOT)
01569         {
01570             b_route= route_num[b];
01571             b_load= route[b].load;
01572             b_len= route[b].length;
01573 
01574             a=VRPH_MAX(pred_array[b],0);
01575             c=VRPH_MAX(next_array[b],0);
01576 
01577             if(a!=k && b!=k && c!=k)// This guarantees a new location
01578             {
01579 
01580                 if(a!=VRPH_DEPOT)
01581                 {
01582                     if(Postsert.evaluate(this,k,a,&M)==true && M.savings<best_savings)
01583                     {
01584                         best_savings=M.savings;
01585                         node1=a;node2=b;
01586                     }
01587                 }
01588                 else
01589                 {
01590 
01591                     if(Presert.evaluate(this,k,b,&M)==true && M.savings<best_savings)
01592                     {
01593                         best_savings=M.savings;
01594                         node1=a;node2=b;
01595                     }
01596 
01597 
01598 
01599                 }
01600 
01601                 // Now try the edge b-c
01602                 if(b!=VRPH_DEPOT)
01603                 {
01604                     if(Postsert.evaluate(this,k,b,&M)==true && M.savings<best_savings)
01605                     {
01606                         best_savings=M.savings;
01607                         node1=b;node2=c;
01608                     }
01609                 }
01610                 else
01611                 {
01612                     // b is the VRPH_DEPOT
01613                     if(Presert.evaluate(this,k,c,&M)==true && M.savings<best_savings)
01614                     {
01615                         best_savings=M.savings;
01616                         node1=b;node2=c;
01617                     }
01618 
01619 
01620                 }
01621 
01622             }
01623             // Advance b
01624             b=VRPH_ABS(next_array[b]);
01625 
01626         }
01627         if(best_savings!=VRP_INFINITY)
01628         {
01629             // Now insert k in between node1 and node2
01630             if(node1!=VRPH_DEPOT)
01631                 Postsert.move(this,k,node1);
01632             else
01633                 Presert.move(this,k,node2);
01634         }
01635 
01636     }
01637 
01638 
01639 
01640     delete [] v;
01641     return true;        
01642 }
01643 
01644 
01645 bool VRP::eject_node(int j)
01646 {
01651 
01652     if(j<=VRPH_DEPOT || routed[j]==false )
01653     {
01654         fprintf(stderr,"Tried to eject index %d\n",j);
01655         report_error("%s: VRPH_DEPOT or unrouted index used in eject node\n",__FUNCTION__);
01656     }
01657 
01658 
01659     int k=j;
01660 
01661     int c, e, k_route, c_route, e_route, k_start, k_end,flag;
01662     double change,ck,ke,ce;
01663 
01664     flag=0;
01665 
01666     // k is no longer routed
01667     routed[k]=false;
01668     
01669     // Get k's current location - between c and e
01670     
01671     c= pred_array[k];
01672     e= next_array[k];
01673 
01674     k_route= route_num[k];
01675     c_route= route_num[VRPH_ABS(c)];
01676     e_route= route_num[VRPH_ABS(e)];
01677 
01678     k_start= route[k_route].start;
01679     k_end= route[k_route].end;
01680 
01681     if(k_start==k && k_end!=k)
01682     {
01683         // c is in a different route
01684         route[k_route].start=VRPH_ABS(e);
01685         next_array[VRPH_ABS(c)]=-VRPH_ABS(e);
01686         pred_array[VRPH_ABS(e)]=-VRPH_ABS(c);
01687         
01688     }
01689 
01690     if(k_end==k && k_start!=k)
01691     {
01692         // e is in a different route
01693         route[k_route].end=c;
01694         next_array[VRPH_ABS(c)]=-VRPH_ABS(e);
01695         pred_array[VRPH_ABS(e)]=-VRPH_ABS(c);
01696     }
01697 
01698     if(k_end==k && k_start==k)
01699     {
01700         // k is its own route
01701         next_array[VRPH_ABS(c)]=-VRPH_ABS(e);
01702         pred_array[VRPH_ABS(e)]=-VRPH_ABS(c);
01703         flag=1;
01704     }
01705 
01706     if(k_start!=k && k_end!=k)
01707     {
01708         // k was interior--easy case
01709         next_array[c]=e;
01710         pred_array[e]=c;
01711 
01712     }
01713 
01714 
01715 
01716     // Now adjust the objective function value
01717     // and the route length and load, etc.
01718 
01719     // old route c-k-e
01720     // new route c-e
01721 
01722     if(e<0)
01723         e=VRPH_DEPOT;
01724     if(c<0)
01725         c=VRPH_DEPOT;
01726 
01727     ce= this->d[c][e];
01728     ck= this->d[c][k];
01729     ke= this->d[k][e];
01730 
01731     change=ce-(ck+ke);
01732     total_route_length+=change;
01733     route[k_route].length+=change;
01734     route[k_route].load-= nodes[k].demand;
01735 
01736     // Removing a node from route k_route
01737     route[k_route].num_customers--;
01738     num_nodes--;
01739 
01740     if(flag)
01741         // We removed a singleton route
01742         this->total_number_of_routes--;
01743     
01744     route_num[k]=-1;
01745     normalize_route_numbers();
01746     
01747     return true;
01748 }
01749 
01750 bool VRP::eject_route(int r, int *route_buff)
01751 {
01756 
01757     int current,cnt;
01758 
01759     // First copy the route into the buffer
01760     current=this->route[r].start;
01761     cnt=0;
01762     while(current>0)
01763     {
01764         route_buff[cnt]=current;
01765         cnt++;
01766         current=this->next_array[current];
01767     }
01768     route_buff[cnt]=-1;
01769     // Terminate with a -1
01770 
01771     
01772 
01773     // Now eject the nodes
01774     for(int i=0;i<cnt;i++)
01775         this->eject_node(route_buff[i]);
01776     
01777     
01778     return true;
01779 
01780 }
01781 bool VRP::check_move(VRPMove *M, int rules)
01782 {
01783 
01788 
01789 
01790     double savings;
01791 
01792     savings=M->savings;
01793 
01794 
01795     if(this->forbid_tiny_moves)
01796     {
01797         // See if it is a "meaningless" move
01798         if(savings>-VRPH_EPSILON && savings < VRPH_EPSILON)
01799             return false;
01800     }
01801 
01802 
01803 
01804     if( (rules & VRPH_FREE) == VRPH_FREE )
01805     {
01806         // use with care!
01807         return true;
01808     }
01809 
01810 
01811     
01812 
01813     if( (rules & VRPH_DOWNHILL) == VRPH_DOWNHILL)
01814     {
01815         if(savings<-VRPH_EPSILON )
01816             return true;
01817         else
01818             return false;
01819 
01820     }
01821 
01822     
01823 
01824 
01825     if( rules & VRPH_RECORD_TO_RECORD )
01826     {
01827         if(savings<=-VRPH_EPSILON)
01828             return true;
01829 
01830         // o/w savings is positive but must be less than deviation*record
01831 
01832 
01833 
01834         if(!has_service_times)
01835         {
01836             if(total_route_length+savings<= (1+deviation)*record)
01837                 return true;
01838             else
01839                 return false;
01840 
01841         }
01842         else
01843         {
01844             // We have service times so remove the service time from the
01845             // deviation calculation
01846 
01847             if( (total_route_length - total_service_time) + savings <=
01848                 ((1+deviation)*(record-total_service_time)))
01849                 return true;
01850             else
01851                 return false;
01852         }
01853     }
01854 
01855     if( rules & VRPH_SIMULATED_ANNEALING )
01856     {
01857         if(M->evaluated_savings==true)
01858             return true;    // We already checked the random acceptance
01859 
01860         // Otherwise, determine whether to accept or not.
01861         if( exp(- (M->savings / this->temperature)) > lcgrand(10) )
01862             return true;
01863         else
01864             return false;
01865 
01866     }
01867 
01868 
01869 
01870     report_error("%s: didn't return yet!\n",__FUNCTION__);
01871 
01872     return false;
01873     
01874 
01875 }
01876 
01877 
01878 bool VRP::is_feasible(VRPMove *M, int rules)
01879 {
01885 
01886     for(int i=0;i<M->num_affected_routes;i++)
01887     {
01888         if( (M->route_lens[i]>this->max_route_length) || (M->route_loads[i]>this->max_veh_capacity) )
01889             return false;
01890     }
01891     
01892     return true;
01893 
01894 }
01895 bool VRP::inject_node(int j)
01896 {
01901 
01902     if(j==VRPH_DEPOT)
01903         report_error("%s: Can't inject VRPH_DEPOT!!\n",__FUNCTION__);
01904 
01905     int edge[4];
01906     double costs[4];
01907     this->find_cheapest_insertion(j,edge,costs,VRPH_USE_NEIGHBOR_LIST); 
01908 
01909     this->insert_node(j,edge[0],edge[1]);
01910 
01911     return true;
01912 
01913 }
01914 
01915 bool VRP::insert_node(int j, int i, int k)
01916 {    
01921 
01922     if(routed[j] || !routed[i] || !routed[k])
01923     {
01924         fprintf(stderr,"insert_node(%d,%d,%d)\n",j,i,k);
01925         report_error("%s: Improper nodes used in insert_node\n",__FUNCTION__);
01926     }
01927 
01928     double increase;
01929     int r;
01930     routed[j]=true;
01931 
01932     if(i==k && k==VRPH_DEPOT)
01933     {
01934         
01935         int last_node=VRPH_ABS(pred_array[VRPH_DEPOT]);
01936         this->next_array[last_node]=-j;
01937         this->pred_array[j]=-last_node;
01938         this->next_array[j]=VRPH_DEPOT;
01939         this->pred_array[VRPH_DEPOT]=-j;
01940 
01941         increase=this->d[0][j]+this->d[j][0];
01942         this->total_number_of_routes++;
01943         this->route_num[j]=this->total_number_of_routes;
01944         this->route[total_number_of_routes].length=increase;
01945         this->route[total_number_of_routes].load=nodes[j].demand;
01946         this->route[total_number_of_routes].num_customers=1;
01947         this->route[total_number_of_routes].start=j;
01948         this->route[total_number_of_routes].end=j;
01949         this->num_nodes++;
01950         this->total_route_length+=increase;
01951 
01952         return true;
01953 
01954 
01955         
01956     }
01957 
01958     if(i!=VRPH_DEPOT && k!=VRPH_DEPOT)
01959     {
01960         
01961         // i-k is interior edge
01962         if(VRPH_MAX(VRPH_DEPOT,next_array[i]) != k || VRPH_MAX(VRPH_DEPOT,pred_array[k])!=i)
01963         {
01964             fprintf(stderr,"edge doesn't exist: next(i(%d)) != k(%d)\n",i,k);
01965             report_error("%s\n",__FUNCTION__);
01966         }
01967 
01968         increase=d[i][j]+d[j][k]-d[i][k];
01969         r=route_num[i];
01970 
01971 
01972         num_nodes++;
01973         next_array[i]=j;
01974         pred_array[j]=i;
01975         next_array[j]=k;
01976         pred_array[k]=j;
01977         route_num[j]=r;
01978         route[r].length+=increase;
01979         route[r].load+= nodes[j].demand;
01980         route[r].num_customers++;
01981         total_route_length+=increase;
01982 
01983         return true;
01984     }
01985 
01986     if(i==VRPH_DEPOT)
01987     {
01988 
01989         increase=d[i][j]+d[j][k]-d[i][k];
01990 
01991         r=route_num[k];
01992 
01993 
01994         int pre=VRPH_ABS(pred_array[k]);
01995 
01996         // pre-VRPH_DEPOT-j-k situation...
01997         num_nodes++;
01998 
01999         next_array[pre]=-j;
02000         pred_array[j]=-pre;
02001         next_array[j]=k;
02002         pred_array[k]=j;
02003         route_num[j]=r;
02004         route[r].start=j;
02005         route[r].length+=increase;
02006         route[r].load+= nodes[j].demand;
02007         route[r].num_customers++;
02008 
02009         total_route_length+=increase;
02010         return true;
02011 
02012     }
02013 
02014     if(k==VRPH_DEPOT)
02015     {
02016 
02017         increase=d[i][j]+d[j][k]-d[i][k];
02018         r=route_num[i];
02019 
02020         int post=VRPH_ABS(next_array[i]);
02021         // i-j-VRPH_DEPOT-post situation...
02022 
02023         num_nodes++;
02024 
02025         next_array[i]=j;
02026         pred_array[j]=i;
02027         next_array[j]=-post;
02028         pred_array[post]=-j;
02029         route_num[j]=r;
02030         route[r].length+=increase;
02031         route[r].load+= nodes[j].demand;
02032         route[r].end=j;
02033         route[r].num_customers++;
02034         total_route_length+=increase;
02035         return true;
02036 
02037     }
02038 
02039     return false;
02040 
02041 
02042 }
02043 
02044 void VRP::find_cheapest_insertion(int j, int *edge, double *costs, int rules)
02045 {
02059 
02060     int h,i,k,m, best_route, new_route, next_node;
02061     double ij,ik,jk,min_feasible_increase, increase, min_increase;
02062 
02063     this->normalize_route_numbers();
02064     best_route = -1;
02065     k=-1;
02066 
02067     // Set the increase to be a singleton route since this must be feasible
02068 
02069     min_feasible_increase = this->d[VRPH_DEPOT][j] + this->d[j][VRPH_DEPOT];
02070     min_increase = this->d[VRPH_DEPOT][j] + this->d[j][VRPH_DEPOT];
02071     edge[0]=edge[1]=edge[2]=edge[3]=VRPH_DEPOT;
02072 
02073 
02074     if(!(rules & VRPH_USE_NEIGHBOR_LIST))
02075     {
02076         i=VRPH_DEPOT;
02077         next_node=VRPH_ABS(next_array[i]);
02078         int cnt=0;
02079         while(next_node != VRPH_DEPOT)
02080         {
02081             cnt++;
02082             if(next_node > 0)
02083             {
02084                 k=next_node;
02085 
02086                 // Try to put j between i and k
02087 
02088                 ij= d[i][j];
02089                 ik= d[i][k];
02090                 jk= d[j][k];
02091                 increase=ij+jk-ik;
02092 
02093                 if(increase<min_increase)
02094                 {
02095                     min_increase=increase;
02096                     edge[2]=i;
02097                     edge[3]=k;
02098                 }
02099                 if(increase<min_feasible_increase)
02100                 {
02101                     // Check feasibility
02102                     new_route= route_num[k];
02103 
02104                     if( (route[new_route].length+increase <= max_route_length) &&
02105                         (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02106                     {
02107                         edge[0]=i;
02108                         edge[1]=k;
02109                         best_route=new_route;
02110                         min_feasible_increase=increase;
02111 
02112                     }
02113 
02114                 }
02115                 // Move on to the next node.
02116                 i=k;
02117 
02118             }
02119             else
02120             {
02121                 // next_node <= 0, implying the edges i-VRPH_DEPOT, VRPH_DEPOT-VRPH_ABS(next_node)
02122                 // First, consider i-VRPH_DEPOT
02123                 k=VRPH_DEPOT;
02124 
02125                 ij= d[i][j];
02126                 ik= d[i][k];
02127                 jk= d[j][k];
02128                 increase=ij+jk-ik;
02129 
02130                 if(increase<min_increase)
02131                 {
02132                     min_increase=increase;
02133                     edge[2]=i;
02134                     edge[3]=k;
02135                 }
02136 
02137                 if(increase<min_feasible_increase)
02138                 { 
02139                     // Check feasibility
02140                     new_route= route_num[i];
02141 
02142                     if( (route[new_route].length+increase <= max_route_length) &&
02143                         (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02144                     {
02145                         edge[0]=i;
02146                         edge[1]=k;
02147                         best_route=new_route;
02148                         min_feasible_increase=increase;
02149                     }
02150                 }
02151 
02152                 // Now consider the edge VRPH_DEPOT-VRPH_ABS(next_node)
02153                 i=VRPH_DEPOT;
02154                 k=VRPH_ABS(next_node);
02155 
02156                 ij= d[i][j];
02157                 ik= d[i][k];
02158                 jk= d[j][k];
02159                 increase=ij+jk-ik;
02160 
02161                 if(increase<min_increase)
02162                 {
02163                     min_increase=increase;
02164                     edge[2]=i;
02165                     edge[3]=k;
02166                 }
02167 
02168                 if(increase<min_feasible_increase)
02169                 {
02170                     new_route= route_num[k];
02171 
02172                     if( (route[new_route].length+increase <= max_route_length )&&
02173                         (route[new_route].load + nodes[j].demand <= max_veh_capacity ) )
02174                     {
02175                         edge[0]=i;
02176                         edge[1]=k;
02177                         best_route=new_route;
02178                         min_feasible_increase=increase;
02179                     }
02180 
02181                 }
02182                 i=k;
02183             }
02184 
02185 
02186             next_node= next_array[i];
02187 
02188         }
02189         // Consider the final edge
02190         k=VRPH_DEPOT;
02191         ij= d[i][j];
02192         ik= d[i][k];
02193         jk= d[j][k];
02194         increase=ij+jk-ik;
02195 
02196         if(increase<min_increase)
02197         {
02198             min_increase=increase;
02199             edge[2]=i;
02200             edge[3]=k;
02201         }
02202 
02203         if(increase<min_feasible_increase)
02204         { 
02205             // Check feasibility
02206             new_route= route_num[i];
02207 
02208             if( (route[new_route].length+increase <= max_route_length) &&
02209                 (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02210             {
02211                 edge[0]=i;
02212                 edge[1]=k;
02213                 best_route=new_route;
02214                 min_feasible_increase=increase;
02215             }
02216         }
02217 
02218         costs[0]=min_feasible_increase;
02219         costs[1]=min_increase;
02220         return;
02221 
02222     }
02223     else
02224     {
02225         // Only search the neighbor list for positions
02226         // We already computed the VRPH_DEPOT-j-VRPH_DEPOT position.
02227         for(m=0;m<neighbor_list_size;m++)
02228         {
02229             i=nodes[j].neighbor_list[m].position;
02230             if(routed[i])
02231             {
02232                 h=VRPH_MAX(VRPH_DEPOT,pred_array[i]);
02233                 increase=d[h][j]+d[j][i]-d[h][i];            
02234 
02235                 if(increase<min_increase)
02236                 {
02237                     min_increase=increase;
02238                     edge[2]=h;
02239                     edge[3]=i;
02240                 }
02241 
02242                 if(increase<min_feasible_increase)
02243                 { 
02244                     // Check feasibility
02245                     if(i!=VRPH_DEPOT)
02246                         new_route= route_num[i];
02247                     else
02248                         new_route= route_num[k];
02249 
02250 
02251                     if( (route[new_route].length+increase <= max_route_length) &&
02252                         (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02253                     {
02254                         edge[0]=h;
02255                         edge[1]=i;
02256                         best_route=new_route;
02257                         min_feasible_increase=increase;
02258                     }
02259                 }
02260 
02261                 k=VRPH_MAX(VRPH_DEPOT,next_array[i]);
02262                 increase=d[i][j]+d[j][k]-d[i][k];    
02263 
02264                 if(increase<min_increase)
02265                 {
02266                     min_increase=increase;
02267                     edge[2]=i;
02268                     edge[3]=k;
02269                 }
02270 
02271                 if(increase<min_feasible_increase)
02272                 { 
02273                     // Check feasibility
02274                     if(i!=VRPH_DEPOT)
02275                         new_route= route_num[i];
02276                     else
02277                         new_route= route_num[k];
02278 
02279 
02280                     if( (route[new_route].length+increase <= max_route_length) &&
02281                         (route[new_route].load + nodes[j].demand <= max_veh_capacity) )
02282                     {
02283                         edge[0]=i;
02284                         edge[1]=k;
02285                         best_route=new_route;
02286                         min_feasible_increase=increase;
02287                     }
02288                 }
02289             }
02290         }
02291 
02292         costs[0]=min_feasible_increase;
02293         costs[1]=min_increase;
02294         return;
02295     }
02296 }
02297 int VRP::inject_set(int num, int *nodelist, int rules, int attempts)
02298 {
02305 
02306     if(rules != VRPH_RANDOM_SEARCH && rules != VRPH_REGRET_SEARCH)
02307         report_error("%s: invalid rules\n",__FUNCTION__);
02308 
02309     int i,j,k;
02310     double best_obj=VRP_INFINITY;
02311     int *best_sol, *orderings,*best_ordering, *start_sol;
02312     int best_index=0;
02313 
02314     best_sol=orderings=best_ordering=start_sol=NULL;
02315 
02316     for(i=0;i<num;i++)
02317     {
02318         if(nodelist[i]==VRPH_DEPOT)
02319         {
02320             fprintf(stderr,"nodelist[%d] of %d=VRPH_DEPOT\n",i,num);
02321             report_error("%s: Cannot inject VRPH_DEPOT!\n",__FUNCTION__);
02322         }
02323     }
02324 
02325     best_sol=new int[3+(this->num_nodes)+num];
02326     start_sol=new int[3+(this->num_nodes)+num];
02327     this->export_solution_buff(start_sol);
02328     this->import_solution_buff(start_sol);
02329     
02330     if(rules==VRPH_RANDOM_SEARCH)
02331     {
02332         // Try "attempts" different random orderings to inject the nodes
02333         orderings=new int[num];
02334         best_ordering=new int[num];
02335         int edge[4];
02336         double costs[2];
02337         for(i=0;i<num;i++)
02338             orderings[i]=i;
02339         
02340         for(i=0;i<attempts;i++)
02341         {
02342             // Create a random permutation
02343             random_permutation(orderings,num);
02344             for(j=0;j<num;j++)
02345             {
02346                 if(nodelist[orderings[j]]==VRPH_DEPOT)
02347                 {
02348                     fprintf(stderr,"inject_set: Bizarre nodelist[%d]:\n",j);
02349                     for(k=0;k<num;k++)
02350                         fprintf(stderr,"(%d,%d) ",orderings[k],nodelist[orderings[k]]);
02351 
02352                     report_error("%s: Found VRPH_DEPOT in nodelist!\n",__FUNCTION__);
02353                 }
02354                 this->inject_node(nodelist[orderings[j]]);
02355             }
02356 
02357             // Look for a new best
02358             // Look for a new best
02359             if(this->total_route_length < best_obj && VRPH_ABS(this->total_route_length-best_obj)>VRPH_EPSILON)
02360             {    
02361                 best_index=i;
02362                 best_obj=this->total_route_length;
02363                 this->export_solution_buff(best_sol);
02364                 memcpy(best_ordering,orderings,num*sizeof(int));
02365             }
02366 
02367             // Now return to the original solution...
02368             for(j=0;j<num;j++)
02369                 this->eject_node(nodelist[orderings[j]]);
02370         }
02371         // Now revisit the best ordering found and try to improve it!
02372         for(j=0;j<num;j++)
02373         {
02374             this->find_cheapest_insertion(nodelist[best_ordering[j]],edge,costs,VRPH_USE_NEIGHBOR_LIST);
02375             this->insert_node(nodelist[best_ordering[j]],edge[0],edge[1]);
02376         }        
02377     }
02378 
02379     if(rules==VRPH_REGRET_SEARCH)
02380     {
02381         // Try "attempts" different random orderings to inject the nodes
02382         orderings=new int[num];
02383         best_ordering=new int[num];
02384 
02385         int *current_list;
02386         current_list=new int[num];
02387         int edge[4];
02388         double costs[2];
02389         for(i=0;i<num;i++)
02390             orderings[i]=i;
02391         
02392         for(i=0;i<attempts;i++)
02393         {
02394             // Create a random permutation
02395             random_permutation(orderings,num);
02396             for(j=0;j<num;j++)
02397                 current_list[j]=nodelist[orderings[j]];
02398 
02399             // It is unfortunately possible to enter some cycles here - this is a cheap
02400             // way of getting out
02401             int cycle_ctr=0;
02402             for(j=0;j<num;j++)
02403             {
02404                 cycle_ctr++;
02405                 if(cycle_ctr==VRPH_MAX_CYCLES)
02406                 {
02407                     fprintf(stderr,"Cycle encountered in REGRET SEARCH!\nReverting to original solution\n");
02408                     if(best_obj<VRP_INFINITY)
02409                     {
02410                         this->import_solution_buff(best_sol);
02411                         delete [] best_sol;
02412                         delete [] orderings;
02413                         delete [] best_ordering;
02414                         delete [] start_sol;
02415                         delete [] current_list;
02416                         return 0;
02417                     }
02418                     else
02419                         report_error("%s: Couldn't escape cycle in REGRET search!\n",__FUNCTION__);
02420                 }
02421                 
02422                 // Find the cheapest way to insert the present node
02423                 this->find_cheapest_insertion(current_list[j],edge,costs,VRPH_USE_NEIGHBOR_LIST);
02424                 // Now loop over the nodes we already injected and compute their ejection costs
02425                 double max_ejection_cost=-VRP_INFINITY;
02426                 double current_cost;
02427                 int node_to_eject=-1;
02428                 for(k=0;k<j;k++)
02429                 {
02430                     // Make sure the node we are testing isn't involved in the cheapest insertion
02431                     if(current_list[k]!=edge[0] && current_list[k]!=edge[1])
02432                     {
02433                         if(( current_cost = this->ejection_cost(current_list[k])) > max_ejection_cost)
02434                         {
02435                             max_ejection_cost=current_cost;
02436                             node_to_eject=current_list[k];
02437                             if(!routed[node_to_eject])
02438                             {
02439                                 fprintf(stderr,"%d NOT ROUTED!!\n",node_to_eject);
02440                                 report_error("%s: Error in inject_set\n",__FUNCTION__);
02441                             }
02442                         }
02443                     }
02444                 }
02445 
02446                 if(node_to_eject==-1 || costs[0]>= max_ejection_cost)
02447                 {
02448                     // Insert
02449                     this->insert_node(current_list[j],edge[0],edge[1]);
02450                 }
02451                 else
02452                 {
02453                     // Eject and Insert
02454                     this->eject_node(node_to_eject);
02455                     this->insert_node(current_list[j],edge[0],edge[1]);
02456                     current_list[j]=node_to_eject;
02457                     j--;
02458                 }        
02459             }
02460 
02461             // Look for a new best
02462             if(this->total_route_length < best_obj && VRPH_ABS(this->total_route_length-best_obj)>VRPH_EPSILON)
02463             {    
02464                 best_index=i;
02465                 best_obj=this->total_route_length;
02466                 this->export_solution_buff(best_sol);
02467                 memcpy(best_ordering,orderings,num*sizeof(int));
02468             }
02469 
02470             // Now return to the original solution...
02471             this->import_solution_buff(start_sol);
02472         }
02473         delete [] current_list;
02474     }
02475 
02476     this->import_solution_buff(best_sol);
02477 
02478     delete [] best_sol;
02479     delete [] orderings;
02480     delete [] best_ordering;
02481     delete [] start_sol;
02482 
02483 
02484     return best_index;
02485 }
02486 
02487 void VRP::eject_neighborhood(int j, int num, int *nodelist)
02488 {
02495 
02496     int i,k,cnt;
02497     int *ejected;
02498 
02499     ejected=new int[this->num_nodes+1];
02500     memset(ejected,0,(this->num_nodes+1)*sizeof(int));
02501 
02502     cnt=0;
02503     
02504     // Always eject j
02505     nodelist[cnt]=j;
02506     ejected[j]=1;
02507     cnt++;
02508 
02509     i=0;
02510     while(cnt<num)
02511     {
02512         // Find a node near j to eject
02513         i=(int)(lcgrand(17)*2*num);
02514         
02515         // Grab it from the neighborlist
02516         i=VRPH_MIN(MAX_NEIGHBORLIST_SIZE-1,i);
02517 
02518         // Get the index
02519         k=this->nodes[j].neighbor_list[i].position;
02520         if(ejected[k]==0)
02521         {
02522             if(k!=VRPH_DEPOT)
02523             {
02524                 nodelist[cnt]=k;
02525                 ejected[k]=1;
02526                 cnt++;
02527             }
02528         }
02529             
02530     }
02531 
02532     // We now have cnt nodes to eject in nodelist
02533     for(i=0;i<cnt;i++)
02534     {
02535         if(nodelist[i]==VRPH_DEPOT)
02536         {
02537             fprintf(stderr,"Trying to eject VRPH_DEPOT??? j=%d\n",j);
02538             report_error("%s: Error in eject_neighborhood\n",__FUNCTION__);
02539         }
02540 
02541         this->eject_node(nodelist[i]);
02542     }
02543 
02544 #if 0
02545     this->verify_routes("After ejecting neighborhood\n");
02546 #endif
02547     delete [] ejected;
02548 
02549     return;
02550 }
02551 
02552 
02553 
02554 
02555 void VRP::normalize_route_numbers()
02556 {
02562 
02563     int current, end, n, R, i;
02564     int current_route;
02565     int *indices;
02566     R= count_num_routes();
02567     n= this->num_original_nodes;
02568 
02569 
02570     // First check to see if we are already normalized
02571 
02572     // Make the arrays 1-based to avoid insanity...
02573     indices=new int[n+1];
02574 
02575     for(i=0;i<=n;i++)
02576         indices[i]=1;
02577 
02578     int ctr=0;
02579     i=VRPH_ABS(next_array[VRPH_DEPOT]);
02580     while(i!=VRPH_DEPOT)
02581     {
02582         if(route_num[i] <= R && routed[i]==true)
02583         {
02584             // route_num[i] is less than R, so we have to
02585             // avoid this index.
02586             indices[route_num[i]]=0;
02587         }
02588         else
02589             ctr++;
02590 
02591         i=VRPH_ABS(next_array[i]);
02592     }
02593     
02594     if(ctr==0)
02595     {
02596         delete [] indices;
02597         return;// route numbers are already normalized since no route_nums were >R.
02598     }
02599 
02600     // Now all of the available indices that are less than or equal to R
02601     // are marked with a 1. The unavailable indices are set to zero.
02602 
02603     int next_index=1;
02604     while(indices[next_index]==0)
02605         next_index++;
02606 
02607     // So next_index is the first index available for use.
02608 
02609     // Get the first route and change the information
02610     current=VRPH_ABS(next_array[VRPH_DEPOT]);
02611 
02612     while(current!=VRPH_DEPOT)
02613     {
02614         current_route= route_num[current];
02615         end= route[current_route].end;
02616 
02617         // Get the next available index
02618         while(indices[next_index]==0)
02619             next_index++;
02620 
02621         // Modify the current index if necessary
02622         if(current_route>R)
02623         {
02624             // Update the route_num
02625             i=VRPH_ABS(next_array[VRPH_DEPOT]);
02626             while(i!=VRPH_DEPOT)
02627             {
02628                 if(route_num[i]==current_route)
02629                     route_num[i]=next_index;
02630                 i=VRPH_ABS(next_array[i]);
02631 
02632             }
02633             // Update the route information
02634             route[next_index].start = route[current_route].start;
02635             route[next_index].end = route[current_route].end;
02636             route[next_index].length = route[current_route].length;
02637             route[next_index].load = route[current_route].load;
02638             route[next_index].num_customers = route[current_route].num_customers;
02639             
02640             next_index++;
02641 
02642 
02643         }
02644 
02645         // Set current to be the first node in the next route
02646         current=VRPH_ABS(next_array[end]);
02647 
02648 
02649     }
02650 
02651     delete [] indices;
02652 
02653     return;
02654 
02655 
02656 
02657 }
02658 
02659 
02660 bool VRP::create_search_neighborhood(int j, int rules)
02661 {
02666 
02667     int i,k, cnt;
02668     // Define the search space
02669     
02670 
02671     if( rules & VRPH_USE_NEIGHBOR_LIST )
02672     {
02673 
02674         // Search only those nodes that are in the neighbor list
02675         search_size=0;
02676         cnt=0;
02677         for(i=0;i<neighbor_list_size;i++)
02678         {
02679             // Consider node k
02680             k=nodes[j].neighbor_list[i].position;
02681             if( routed[k] == true)
02682             {
02683                 // The node is routed 
02684                 if(rules & VRPH_INTER_ROUTE_ONLY)
02685                 {
02686                     // Make sure k is in a different route
02687                     if(route_num[k]!=route_num[j])
02688                     {
02689                         search_space[cnt] = k;
02690                         cnt++;
02691                         
02692                     }
02693                 }
02694                 else
02695                 {
02696                     if(rules & VRPH_INTRA_ROUTE_ONLY)
02697                     {
02698                         // Make sure k is in the same route
02699 
02700                         if(route_num[k]==route_num[j])
02701                         {
02702                             // The node is in a different route - add to the search space
02703                             search_space[cnt] = k;
02704                             cnt++;
02705                             
02706                         }
02707                     }
02708                     else
02709                     {
02710                         // No intra/inter restrictions
02711                         search_space[cnt] = k;
02712                         cnt++;
02713                         
02714             
02715                     }
02716                 }
02717             }
02718         }
02719         search_size=cnt;
02720         goto randomize;
02721     }
02722 
02723 
02724     // Otherwise no neighborlist
02725     if( (rules & VRPH_INTRA_ROUTE_ONLY) )
02726     {        
02727         // Search_space is just the route itself
02728         search_space[0]=VRPH_DEPOT;
02729         search_space[1]=this->route[route_num[j]].start;
02730         for(i=2;i<=this->route[route_num[j]].num_customers;i++)
02731             search_space[i]=this->next_array[search_space[i-1]];
02732 
02733         search_size=this->route[route_num[j]].num_customers+1;//add 1 for depot
02734 
02735         goto randomize;
02736 
02737         /*
02738 
02739         search_size=0;
02740         i=0;
02741         search_space[i]=VRPH_DEPOT;
02742         i++;
02743         search_space[i]=route[route_num[j]].start;
02744         for(;;)
02745         {
02746             search_space[i]=next_array[search_space[i]];
02747             if(search_space[i]<=VRPH_DEPOT)
02748             {
02749                 search_size=i;
02750                 goto randomize;
02751             }
02752             i++;
02753         }*/
02754 
02755     }
02756 
02757     if( (rules & VRPH_INTER_ROUTE_ONLY) )
02758     {        
02759         // Only nodes in a different route        
02760         i=0;
02761         search_space[i]=VRPH_DEPOT;
02762         search_size=1;
02763         k=VRPH_ABS(next_array[search_space[i]]);    
02764         i++;
02765 
02766         for(;;)
02767         {
02768             if(k==VRPH_DEPOT)
02769                 // end of solution
02770                 goto randomize;
02771 
02772             if(route_num[k]!=route_num[j])
02773             {
02774                 search_space[i]=k;
02775                 i++;
02776                 search_size++;
02777             }
02778             k=VRPH_ABS(next_array[k]);
02779         }
02780 
02781     }
02782 
02783     
02784     // Search space will consist of all positions in routes moving in the
02785     // given direction
02786 
02787     if((rules & VRPH_FORWARD))
02788     {
02789         search_size=0;
02790 
02791         i=0;
02792         search_space[i]=VRPH_ABS(next_array[route[route_num[j]].end]);
02793         search_size++;
02794 
02795         // First node in next route in the solution
02796         for(;;)
02797         {
02798             search_space[i+1]=VRPH_ABS(next_array[search_space[i]]);    
02799             search_size++;
02800             if(search_space[i+1]==VRPH_DEPOT)
02801                 goto randomize;
02802             i++;
02803 
02804         }
02805 
02806     }
02807 
02808     if((rules & VRPH_BACKWARD))
02809     {
02810         search_size=0;
02811 
02812         i=0;
02813         search_space[i]=VRPH_ABS(pred_array[route[route_num[j]].start]);
02814         search_size++;
02815 
02816         // Last node in previous route in the solution
02817         for(;;)
02818         {
02819             search_space[i+1]=VRPH_ABS(pred_array[search_space[i]]);    
02820             search_size++;
02821             if(search_space[i+1]==VRPH_DEPOT)
02822                 goto randomize;
02823             i++;
02824 
02825         }
02826     
02827     }
02828 
02829     // No rules given - search space is set of all nodes
02830     search_size=0;
02831     i=0;
02832 
02833     search_space[i]=VRPH_ABS(next_array[VRPH_DEPOT]);
02834     search_size++;
02835 
02836     for(;;)
02837     {
02838         search_space[i+1]=VRPH_ABS(next_array[search_space[i]]);    
02839         search_size++;
02840         if(search_space[i+1]==VRPH_DEPOT)
02841             goto randomize;
02842         i++;
02843     }
02844 
02845 randomize:
02846 
02847     // Now permute it if the rules is VRPH_RANDOMIZED
02848     if(rules & VRPH_RANDOMIZED)
02849         random_permutation(search_space, search_size);
02850 
02851     return true;
02852     
02853 
02854 
02855 }
02856 
02857 
02858 double VRP::insertion_cost(int u, int a, int b)
02859 {
02864 
02865     // a-b is an existing edge
02866     // u is the node we are trying to insert between a and b
02867     // Returns the increased cost to the route containing a-b of
02868     // inserting u in between a and b
02869 
02870     if(a==b && b==VRPH_DEPOT)
02871         return d[VRPH_DEPOT][u]+d[u][VRPH_DEPOT];// assume feasible!!
02872     
02873     if(a==u || u==b || u==VRPH_DEPOT) 
02874         report_error("%s: overlap or VRPH_DEPOT found\n",__FUNCTION__);
02875 
02876     if(routed[a]==false || routed[b]==false)
02877     {
02878         fprintf(stderr,"%d,%d:  not routed!\n",a,b);
02879         report_error("%s: Unrouted nodes in insertion_cost\n",__FUNCTION__);
02880     }
02881 
02882     if( (a!=VRPH_DEPOT) && VRPH_MAX(next_array[a],VRPH_DEPOT)!=b)
02883     {
02884         fprintf(stderr,"%d,%d:  not an edge!\n",a,b);
02885         report_error("%s: Invalid nodes in insertion_cost\n",__FUNCTION__);
02886     }
02887 
02888     if( (a==VRPH_DEPOT) && VRPH_MAX(pred_array[b],VRPH_DEPOT)!=a)
02889     {
02890         fprintf(stderr,"%d,%d:  not an edge!\n",a,b);
02891         report_error("%s: Invalid nodes in insertion_cost\n",__FUNCTION__);
02892     }
02893 
02894 
02895     int new_route;
02896     if(a==VRPH_DEPOT)
02897         new_route=route_num[b];
02898     else
02899         new_route=route_num[a];
02900 
02901     if(nodes[u].demand+route[new_route].load>max_veh_capacity)
02902         return VRP_INFEASIBLE;
02903     
02904     double increase=d[a][u] + d[u][b] - d[a][b];
02905     if(route[new_route].length+increase>max_route_length)
02906         return VRP_INFEASIBLE;
02907 
02908     return increase;
02909 
02910 }
02911 
02912 
02913 
02914 double VRP::ejection_cost(int u)
02915 {
02920 
02921     if( u==VRPH_DEPOT)
02922         report_error("%s: Cannot eject the VRPH_DEPOT!\n",__FUNCTION__);
02923 
02924     if(!routed[u])
02925         return -VRP_INFINITY;
02926 
02927     // t-u-v --> t-v
02928 
02929     return this->d[VRPH_MAX(pred_array[u],VRPH_DEPOT)][u]+d[u][VRPH_MAX(next_array[u],VRPH_DEPOT)] -
02930         d[VRPH_MAX(pred_array[u],VRPH_DEPOT)][VRPH_MAX(next_array[u],VRPH_DEPOT)];
02931 
02932 
02933 }
02934 void VRP::clean_route(int r, int heuristics)
02935 {
02940 
02941     int i,j, r_start, r_end;
02942     double start_val, end_val, start_rlen, end_rlen;
02943 
02944     OnePointMove    OPM;
02945     TwoPointMove    TPM;
02946     TwoOpt            TO;
02947     ThreeOpt        ThreeO;
02948     OrOpt            OrO;
02949 
02950     int rules= VRPH_INTRA_ROUTE_ONLY+VRPH_DOWNHILL+VRPH_FIRST_ACCEPT+VRPH_SAVINGS_ONLY;
02951 
02952 
02953 start_improving:
02954 
02955     start_val = route[r].length;
02956 
02957 #if CLEAN_DEBUG
02958     printf("CLEAN::start_val=%f\n",start_val);
02959 #endif
02960 
02961     end_val = - VRP_INFINITY;
02962     r_start=route[r].start;
02963     r_end=route[r].end;
02964 
02965     if((heuristics & ONE_POINT_MOVE)==ONE_POINT_MOVE)
02966     {
02967         // Try the ONE_POINT_MOVE for all nodes in route r;
02968         r_start=route[r].start;
02969 
02970 opm_search:
02971         start_rlen=route[r].length;
02972         i=r_start;
02973         j=VRPH_MAX(next_array[i],0);
02974         while(i!=VRPH_DEPOT)
02975         {
02976             while(OPM.search(this,i,rules));
02977             i=j;
02978             j=VRPH_MAX(next_array[j],0);
02979         }
02980 
02981         end_rlen = route[r].length;
02982         if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
02983             goto opm_search;
02984 
02985 #if CLEAN_DEBUG
02986         printf("CLEAN::OPM end_rlen=%f\n",end_rlen);
02987 #endif
02988 
02989     }
02990 
02991     if((heuristics & TWO_POINT_MOVE)==TWO_POINT_MOVE)
02992     {
02993         // Try the TWO_POINT_MOVE for all nodes in route r;
02994 tpm_search:
02995         r_start=route[r].start;
02996         start_rlen=route[r].length;
02997         i=r_start;
02998         j=VRPH_MAX(next_array[i],0);
02999         while(i!=VRPH_DEPOT)
03000         {
03001             while(TPM.search(this,i,rules));
03002             i=VRPH_MAX(next_array[i],VRPH_DEPOT);
03003         }
03004 
03005         end_rlen = route[r].length;
03006         if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
03007             goto tpm_search;
03008 
03009 #if CLEAN_DEBUG
03010         printf("CLEAN::TPM end_rlen=%f\n",end_rlen);
03011 #endif
03012 
03013 
03014     }
03015 
03016     if((heuristics & TWO_OPT)==TWO_OPT)
03017     {
03018 
03019 to_search:
03020         r_start=route[r].start;
03021         start_rlen=route[r].length;
03022         i=r_start;
03023         j=VRPH_MAX(next_array[i],0);
03024         while(i!=VRPH_DEPOT)
03025         {
03026             while(TO.search(this,i,rules));
03027             i=VRPH_MAX(next_array[i],0);
03028         }
03029 
03030         end_rlen = route[r].length;
03031         if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
03032             goto to_search;
03033 
03034 #if CLEAN_DEBUG
03035         printf("CLEAN::TO end_rlen=%f\n",end_rlen);
03036 #endif
03037 
03038 
03039     }
03040 
03041     if((heuristics & OR_OPT)==OR_OPT)
03042     {
03043         // Use strings of length 2 and 3
03044 or_search:
03045         r_start=route[r].start;
03046         start_rlen=route[r].length;
03047         i=r_start;
03048         j=VRPH_MAX(next_array[i],0);
03049         while(i!=VRPH_DEPOT)
03050         {
03051             OrO.search(this,i,3,rules);
03052             i=VRPH_MAX(next_array[i],0);
03053         }
03054 
03055         i=r_start;
03056         j=VRPH_MAX(next_array[i],0);
03057         while(i!=VRPH_DEPOT)
03058         {
03059             OrO.search(this,i,2,rules);
03060             i=VRPH_MAX(next_array[i],0);
03061         }
03062         end_rlen = route[r].length;
03063         if( (end_rlen<start_rlen) && (VRPH_ABS(end_rlen-start_rlen)>VRPH_EPSILON))
03064             goto or_search;
03065 
03066 #if CLEAN_DEBUG
03067         printf("CLEAN::TO end_rlen=%f\n",end_rlen);
03068 #endif
03069 
03070     }
03071 
03072     if((heuristics & THREE_OPT)==THREE_OPT)
03073     {
03074         //show_route(r);
03075 
03076         while(ThreeO.route_search(this,r,rules))
03077 
03078             end_val = route[r].length;
03079 #if CLEAN_DEBUG
03080         printf("CLEAN::3O end_val=%f\n",end_val);
03081 #endif
03082 
03083 
03084     }
03085 
03086     end_val= route[r].length;
03087 
03088 #if CLEAN_DEBUG
03089     printf("CLEAN::final end_val=%f\n",end_val);
03090 #endif
03091 
03092 
03093     if(VRPH_ABS(start_val-end_val)>VRPH_EPSILON)
03094         goto start_improving;
03095     else
03096         return;
03097 
03098 }
03099 
03100 bool VRP::before(int a, int b)
03101 {
03107 
03108     int i;
03109 
03110     if(a==VRPH_DEPOT || b==VRPH_DEPOT)
03111         report_error("%s: before called with VRPH_DEPOT\n",__FUNCTION__);
03112     
03113     if(route_num[a]!=route_num[b])
03114     {
03115         fprintf(stderr,"Ordering error: before called with %d and %d not in the same route!\n",a,b);
03116         report_error("%s: differnet routes\n",__FUNCTION__);
03117     }
03118 
03119     if(next_array[a]==b)
03120         return true;
03121     if(next_array[b]==a)
03122         return false;
03123 
03124     i=a;
03125     while(i>0 && i!=b)
03126         i=next_array[i];
03127 
03128     // At the end of this loop, if i<=0, then we're at the end of a route
03129     // and haven't encountered b ==> must have a after b in the route
03130     // If i==b, then we know that b follows a
03131 
03132     if(i==b)
03133         return true;
03134     else
03135         return false;
03136 
03137 }
03138 
03139 
03140 void VRP::update(VRPMove *M)
03141 {
03145 
03146     if(M->num_affected_routes==0)    // Nothing to do!
03147         return;    
03148 
03149     int i;
03150 
03151     for(i=0; i < M->num_affected_routes; i++)
03152     {
03153         // Update length
03154         route[M->route_nums[i]].length = M->route_lens[i];
03155 
03156         // Update load
03157         route[M->route_nums[i]].load = M->route_loads[i];
03158 
03159         // Update # of customers
03160         route[M->route_nums[i]].num_customers = M->route_custs[i];
03161     }
03162 
03163     // Now update total_route_length
03164     total_route_length = M->new_total_route_length;
03165 
03166     // # of routes
03167     total_number_of_routes = M->total_number_of_routes; 
03168 
03169     return;
03170 
03171 }
03172 
03173 
03174 
03175 
03176 void VRP::compute_route_center(int r)
03177 {
03182 
03183     int current_node = route[r].start;
03184     double total_x=0;
03185     double total_y=0;
03186 
03187     while(current_node != VRPH_DEPOT)
03188     {
03189         total_x += nodes[current_node].x;
03190         total_y += nodes[current_node].y;
03191         current_node=VRPH_MAX(VRPH_DEPOT, next_array[current_node]);
03192     }
03193 
03194     route[r].x_center = total_x / ((double) route[r].num_customers);
03195     route[r].y_center = total_y / ((double) route[r].num_customers);    
03196 
03197 }
03198 void VRP::find_neighboring_routes()
03199 {
03205 
03206     int i,j;
03207     VRPNeighborElement **rd;
03208 
03209     // Compute the route centers
03210     normalize_route_numbers();
03211     for(i=1;i<=total_number_of_routes;i++)
03212     {
03213         compute_route_center(i);
03214     }
03215 
03216     // Create the matrix
03217     rd = new VRPNeighborElement *[total_number_of_routes + 1];
03218     rd[0]= new VRPNeighborElement[(total_number_of_routes +1) * (total_number_of_routes + 1)];
03219     for(i=1;i<total_number_of_routes+1;i++)
03220         rd[i]=rd[i-1] + (total_number_of_routes+1);
03221 
03222     // Fill the matrix with inter-route distances
03223     for(i=1;i<=this->total_number_of_routes;i++)
03224     {
03225         rd[i][0].position=VRP_INFINITY;
03226         rd[i][0].val=VRP_INFINITY;
03227         for(j=1;j<=this->total_number_of_routes;j++)
03228         {
03229             rd[i][j].position = j;
03230             rd[i][j].val = VRPDistance(VRPH_EUC_2D,route[i].x_center, route[i].y_center,
03231                 route[j].x_center,route[j].y_center);
03232         }
03233     }
03234 
03235     // Now sort each row in put the top MAX_NEIGHBORING_ROUTES in the route structs.
03236 
03237     
03238     for(i=1;i<=total_number_of_routes;i++)
03239         qsort(rd[i],total_number_of_routes+1,sizeof(VRPNeighborElement),VRPNeighborCompare);
03240     
03241     
03242 
03243     for(i=1;i<=total_number_of_routes;i++)
03244     {
03245         for(j=0;j<MAX_NEIGHBORING_ROUTES;j++)
03246         {
03247             route[i].neighboring_routes[j]=rd[i][j+1].position;
03248             // don't include the route itself!
03249         }
03250     }
03251 
03252     delete [] rd[0];
03253     delete [] rd;
03254 
03255     
03256 
03257     return;
03258 }
03259 
03260 void VRP::capture_best_solution()
03261 {
03265 
03266     if( (this->total_route_length < this->best_total_route_length) && 
03267         (VRPH_ABS(this->total_route_length - this->best_total_route_length) > VRPH_EPSILON) )
03268     {
03269         this->best_total_route_length=this->total_route_length;
03270         this->export_solution_buff(this->best_sol_buff);
03271         
03272     }
03273 
03274     
03275     if(this->total_route_length < this->solution_wh->worst_obj || 
03276         this->solution_wh->num_sols < this->solution_wh->max_size)
03277     {
03278         VRPSolution this_sol(this->num_nodes);
03279 
03280         this_sol.obj=this->total_route_length;
03281         this_sol.in_IP=false;
03282 
03283         // Export buffer
03284         this->export_canonical_solution_buff(this_sol.sol);
03285         
03286         this->solution_wh->add_sol(&this_sol, 0); 
03287         // We don't know any information to help us know where to insert
03288     }
03289 
03290     return;
03291 
03292     
03293 }
03294 
03295 void VRP::update_solution_wh()
03296 {
03300 
03301     VRPSolution this_sol(this->num_nodes);
03302 
03303     this_sol.obj=this->total_route_length;
03304     this_sol.in_IP=false;
03305 
03306     // Export buffer
03307     this->export_canonical_solution_buff(this_sol.sol);
03308 
03309     this->solution_wh->add_sol(&this_sol, 0); 
03310     // We don't know any information to help us know where to insert
03311     
03312 
03313     return ;
03314 
03315 
03316 }
03317 
03318 
03319 void VRP::update_route(int j, VRPRoute *R)
03320 {
03326 
03327     int i, current;
03328     double st=0;// for the service time
03329 
03330     R->x[0]=this->nodes[VRPH_DEPOT].x;
03331     R->y[0]=this->nodes[VRPH_DEPOT].y;
03332     R->start=route[j].start;
03333     R->end=route[j].end;
03334     R->num_customers=route[j].num_customers;
03335     R->load=route[j].load;
03336     R->length=route[j].length;
03337     R->obj_val=this->total_route_length-this->total_service_time;
03338 
03339 
03340     if(R->start < R->end)
03341     {
03342         current=R->start;
03343         // Just output the normal ordering
03344         R->ordering[0]=current;
03345         R->x[1]=this->nodes[current].x;
03346         R->y[1]=this->nodes[current].y;
03347         st+=this->nodes[current].service_time;
03348 
03349         for(i=1; i<R->num_customers; i++)
03350         {
03351             current=this->next_array[current];
03352             st+=this->nodes[current].service_time;
03353             R->ordering[i]=current;
03354             R->x[i+1]=this->nodes[current].x;
03355             R->y[i+1]=this->nodes[current].y;
03356         }
03357 
03358         R->total_service_time=st;
03359         
03360         return;
03361     }
03362 
03363     // Otherwise start > end - store the reverse
03364     current=R->end;
03365     R->ordering[0]=current;
03366     R->x[1]=this->nodes[current].x;
03367     R->y[1]=this->nodes[current].y;
03368     st+=this->nodes[current].service_time;
03369 
03370 
03371     for(i=1; i<R->num_customers; i++)
03372     {
03373         current=this->pred_array[current];
03374         st+=this->nodes[current].service_time;
03375         R->ordering[i]=current;
03376         R->x[i+1]=this->nodes[current].x;
03377         R->y[i+1]=this->nodes[current].y;
03378 
03379     }
03380 
03381     R->total_service_time=st;
03382 
03383     return;
03384 
03385 }
03386 
03387 
03388 
03389 double VRP::split(double p)
03390 {
03398 
03399     if(p>.5)
03400         report_error("%s: p must be less than .5\n");
03401 
03402     int i,j,k;
03403     double beta;
03404     struct double_int *thetas;
03405 
03406     thetas=new double_int[this->num_nodes];
03407 
03408     for(i=1;i<=this->num_nodes;i++)
03409     {
03410         thetas[i-1].k=i;
03411         thetas[i-1].d=this->nodes[i].theta;
03412     }
03413 
03414     // Sort the list of thetas
03415     qsort(thetas,this->num_nodes,sizeof(struct double_int),double_int_compare);
03416 
03417     // Select a random beta in (min_theta, max_theta ) and see how the split is
03418     // Stop when we get a satisfactory split
03419 
03420     
03421     for(;;)
03422     {
03423         beta = this->min_theta + (double)(lcgrand(10))*.5*(max_theta - min_theta);
03424 
03425         // Now find out how many nodes we get in each half
03426         k=0;
03427         for(j=0;j<this->num_nodes;j++)
03428         {
03429             // Count the nodes above the line
03430             if(this->nodes[j+1].y >= tan(beta)*(this->nodes[j+1].x))
03431                 k++;
03432         }
03433     
03434         if( (k>=(int)(p*(double)(this->num_nodes))) && (k<=(int)((1-p)*(double)(this->num_nodes))))
03435             break;    // The split is good
03436     }
03437 
03438     // Eject them - note that we assume to have begun with a full solution...
03439     for(i=1;i<=this->num_original_nodes;i++)
03440     {
03441         if(routed[i])
03442         {
03443             if(this->nodes[i].y >= tan(beta)*(this->nodes[i].x))
03444                 this->eject_node(i);
03445         }
03446     }
03447 
03448     
03449     delete [] thetas;
03450     return beta;
03451 }
03452 
03453 
03454 int VRP::split_routes(double p, int **ejected_routes, double *t)
03455 {
03466 
03467     if(p>.5)
03468         report_error("%s: p must be less than .5\n");
03469 
03470     int i,j,k;
03471     double beta;
03472     struct double_int *thetas;
03473 
03474     thetas=new double_int[this->num_nodes];
03475 
03476     for(i=1;i<=this->num_nodes;i++)
03477     {
03478         thetas[i-1].k=i;
03479         thetas[i-1].d=this->nodes[i].theta;
03480     }
03481 
03482     // Sort the list of thetas
03483     qsort(thetas,this->num_nodes,sizeof(struct double_int),double_int_compare);
03484 
03485     // Select a random beta in (min_theta, max_theta ) and see how the split is
03486     // Stop when we get a satisfactory split
03487 
03488     
03489     for(;;)
03490     {
03491         beta = this->min_theta + (double)(lcgrand(10))*.5*(max_theta - min_theta);
03492 
03493         // Now find out how many nodes we get in each half
03494         k=0;
03495         for(j=0;j<this->num_nodes;j++)
03496         {
03497             // Count the nodes above the line
03498             if(this->nodes[j+1].y >= tan(beta)*(this->nodes[j+1].x))
03499                 k++;
03500         }
03501     
03502         if( (k>=(int)(p*(double)(this->num_nodes))) && (k<=(int)((1-p)*(double)(this->num_nodes))))
03503         {
03504             break;    // The split is good
03505         }
03506     }
03507 
03508     // Eject the routes that do not have any nodes in the larger portion
03509 
03510     int num_ejected=0;
03511     for(i=1;i<=this->total_number_of_routes;i++)
03512     {
03513         int start=this->route[i].start;
03514         int current=start;
03515         bool will_eject=true;
03516 
03517         while(current!=VRPH_DEPOT)
03518         {
03519             
03520             if(this->nodes[current].y <= tan(beta)*(this->nodes[current].x))
03521             {
03522                 // We have at least one node inside the large part
03523                 will_eject=false;
03524                 break;
03525             }
03526 
03527             current=VRPH_MAX(VRPH_DEPOT,this->next_array[current]);
03528         }
03529         if(will_eject)
03530         {
03531             // Eject route i from the solution - place the ejected nodes in 
03532             // the ejected_routes[][] array
03533             this->eject_route(i, ejected_routes[num_ejected]);
03534             num_ejected++;
03535 
03536         }
03537         
03538     }
03539     
03540     delete [] thetas;
03541 
03542     *t=beta;
03543     return num_ejected;
03544 }
03545 
03546 
03547 void VRP::fix_edge(int start, int end)
03548 {
03556 
03557     this->fixed[start][end]=true;
03558     this->fixed[end][start]=true;
03559 
03560     // Handle the dummy_node
03561     if(start==VRPH_DEPOT)
03562     {
03563         this->fixed[dummy_index][start]=true;
03564         this->fixed[start][dummy_index]=true;
03565     }
03566 
03567     if(end==VRPH_DEPOT)
03568     {
03569         this->fixed[dummy_index][end]=true;
03570         this->fixed[end][dummy_index]=true;
03571     }
03572 
03573 
03574 }
03575 
03576 void VRP::unfix_edge(int start, int end)
03577 {
03581 
03582     if(this->fixed[start][end])
03583         report_error("%s: Edge %d-%d is not already fixed!\n",__FUNCTION__,start,end);
03584 
03585     this->fixed[start][end]=false;
03586     this->fixed[end][start]=false;
03587 
03588     // Handle the dummy_node
03589     if(start==VRPH_DEPOT)
03590     {
03591         this->fixed[dummy_index][start]=false;
03592         this->fixed[start][dummy_index]=false;
03593     }
03594 
03595     if(end==VRPH_DEPOT)
03596     {
03597         this->fixed[dummy_index][end]=false;
03598         this->fixed[end][dummy_index]=false;
03599     }
03600 
03601 
03602 }
03603 
03604 void VRP::unfix_all()
03605 {
03609 
03610     for(int i=0;i<=matrix_size;i++)
03611         for(int j=0;j<=matrix_size;j++)
03612             fixed[i][j]=false;
03613 }
03614 
03615 void VRP::fix_string(int *node_string, int k)
03616 {
03623 
03624 
03625     for(int i=0;i<k-1;i++)
03626         this->fix_edge(node_string[i], node_string[i+1]);
03627 
03628     
03629 }
03630 
03631 
03632 
03633 
03634 
03635 int VRP::read_fixed_edges(const char *filename)
03636 {
03649 
03650 
03651     int i, a, b, k;
03652     FILE *in;
03653 
03654     if( (in=fopen(filename,"r"))==NULL)
03655     {
03656         fprintf(stderr,"Error opening %s for reading\n",filename);
03657         report_error("%s: Bad file in read_fixed_edges\n",__FUNCTION__);
03658     }
03659 
03660     fscanf(in,"%d\n",&k);
03661 
03662     for(i=0;i<k;i++)
03663     {
03664         fscanf(in,"%d %d\n",&a, &b);
03665         if(a<0 || b<0 || a>this->num_original_nodes || b>this->num_original_nodes)
03666         {
03667             fprintf(stderr,"Tried to fix edge %d-%d\n",a,b);
03668             report_error("%s: Error in read_fixed_edges\n",__FUNCTION__);
03669         }
03670         this->fix_edge(a,b);
03671     }
03672 
03673     // Return the number of fixed edges read in
03674     return k;
03675 }
03676 void VRP::list_fixed_edges(int *fixed_list)
03677 {
03684 
03685     int pos=0;
03686     int current=VRPH_DEPOT;
03687     int next=VRPH_ABS(next_array[current]);
03688 
03689     while(next!=VRPH_DEPOT)
03690     {
03691         if(fixed[current][next])
03692         {
03693             fixed_list[pos]=current;
03694             fixed_list[pos+1]=next;
03695             pos+=2;
03696         }
03697 
03698         current=next;
03699         next=next_array[current];
03700 
03701         if(next<0)
03702         {
03703             if(fixed[current][VRPH_DEPOT])
03704             {
03705                 fixed_list[pos]=current;
03706                 fixed_list[pos+1]=VRPH_DEPOT;
03707                 pos+=2;
03708             }
03709 
03710             if(fixed[VRPH_DEPOT][-next])
03711             {
03712                 fixed_list[pos]=VRPH_DEPOT;
03713                 fixed_list[pos+1]=-next;
03714                 pos+=2;
03715             }
03716 
03717             next=-next;
03718         }
03719     }
03720 
03721 
03722 
03723 }
03724 void VRP::perturb_locations(double c)
03725 {
03732 
03733     int i, pre, post;
03734     double v,theta;
03735 
03736     // First export the solution buffer
03737     this->export_solution_buff(this->current_sol_buff);
03738 
03739     for(i=1;i<=this->num_nodes;i++)
03740     {
03741         pre= VRPH_MAX(VRPH_DEPOT,this->pred_array[i]);
03742         post=VRPH_MAX(VRPH_DEPOT,this->next_array[i]);
03743 
03744         v= (this->d[pre][i]+this->d[i][post]) - (this->nodes[i].service_time + this->nodes[post].service_time);
03745         v= v*c;
03746 
03747         // Now alter the location of node i by shifting v units in a random direction theta
03748         theta= VRPH_PI*lcgrand(8);
03749 
03750         this->nodes[i].x += v*cos(theta);
03751         this->nodes[i].y += v*sin(theta);
03752         
03753     }
03754 
03755     // Now re-create the distance matrix
03756 
03757     this->create_distance_matrix(this->edge_weight_type);
03758 
03759     // Now re-import the solution with the new distances
03760     this->import_solution_buff(this->current_sol_buff);
03761 
03762     return;
03763 
03764 }
03765 void VRP::add_route(int *route_buff)
03766 {
03772 
03773     // We will do this by exporting the current solution to a buffer, appending the 
03774     // new route to this buffer, and then importing the new resulting solution.
03775     // Probably not the fastest way, but we have many fields to update when adding a
03776     // new route!
03777 
03778     int *temp_buff;
03779 
03780     this->verify_routes("Before adding route\n");
03781 
03782     temp_buff=new int[this->num_original_nodes+2];
03783     this->export_solution_buff(temp_buff);
03784 
03785 
03786     int old_num=this->num_nodes;
03787     int i=0;
03788     while(route_buff[i]!=-1)
03789     {
03790         temp_buff[old_num+1+i]=route_buff[i];
03791         if(i==0)
03792             temp_buff[old_num+1+i]=-temp_buff[old_num+1+i];
03793 
03794         temp_buff[0]++;
03795         i++;
03796     }
03797     
03798     temp_buff[old_num+1+i]=VRPH_DEPOT;
03799     
03800     // Now import the new solution
03801     this->import_solution_buff(temp_buff);
03802     
03803     this->verify_routes("After adding route\n");
03804 
03805     delete [] temp_buff;
03806 
03807 }
03808 void VRP::append_route(int *sol_buff, int *route_buff)
03809 {
03816 
03817     int i,j,current_num;
03818 
03819     current_num=sol_buff[0]; // The # of nodes in the sol_buff solution
03820 
03821     // Now find out how many nodes are in the route_buff - ends in a -1;
03822 
03823     j=0;
03824     while(route_buff[j]!=-1)
03825         j++;
03826 
03827     // The route has j nodes - we assume they are all different from those in the solution!
03828     // Increment the # of nodes in sol_buff;
03829     sol_buff[0]+=j;
03830 
03831     sol_buff[current_num+1]=-route_buff[0];
03832     for(i=1;i<j;i++)
03833         sol_buff[current_num+1+i]=route_buff[i];
03834     
03835 
03836     // End at the VRPH_DEPOT;
03837     sol_buff[current_num+1+j]=VRPH_DEPOT;
03838         
03839     return;
03840 
03841 
03842 }
03843 
03844 int VRP::intersect_solutions(int *new_sol, int **route_list, int *sol1, int *sol2, int min_routes)
03845 {
03856 
03857     int i,j,k,m,num_routes;
03858     int *rnums;
03859 
03860     rnums=new int[this->num_original_nodes+1]; // Plenty big...
03861 
03862     // Find the routes in common
03863     j = this->find_common_routes(sol1, sol2, rnums);
03864     
03865     if(j==0)
03866     {
03867         // No routes in common - just put sol1 into new_sol
03868         memcpy(new_sol,sol1, (this->num_original_nodes+2)*sizeof(int));
03869         delete[] rnums;
03870         return 0;
03871     }
03872 
03873     // Otherwise we have some routes in common - copy them to the routes[] array
03874     // Import sol1 to do this
03875     this->import_solution_buff(sol1);
03876     num_routes=this->total_number_of_routes;
03877     k=0;
03878     for(i=0;i<j;i++)
03879     {
03880         // Copy route i from sol1 into the route_list[] array
03881         route_list[i][0]=this->route[rnums[i]].start;
03882 
03883         m=0;
03884         while(route_list[i][m]!=this->route[rnums[i]].end)
03885         {
03886             m++;
03887             route_list[i][m]=this->next_array[route_list[i][m-1]];
03888         }
03889         // Add a -1 to indicate the end of the route
03890         m++;
03891         route_list[i][m]=-1;
03892 
03893         k++;
03894         num_routes--;
03895         if(num_routes==min_routes)
03896             break;
03897     }
03898 
03899     // We actually copied k of the j routes in the intersection
03900     // k may be less than j
03901 
03902     // Now eject the routes from sol1 (already imported)
03903     int *junk;
03904     junk=new int[this->num_original_nodes];
03905     // won't use this
03906     for(i=0;i<k;i++)
03907         this->eject_route(rnums[i], junk);
03908 
03909     // Now export the solution after ejections to new_sol
03910     this->export_canonical_solution_buff(new_sol);
03911 
03912     // Now import the new (possibly partial) solution
03913     this->import_solution_buff(new_sol);
03914 
03915     delete [] rnums;
03916     delete [] junk;
03917 
03918     return k;
03919 
03920 }
03921 bool VRP::osman_insert(int k, double alpha)
03922 {
03930 
03931     int i,j,l,ll,m,mm, best_l, best_m;
03932     double savings, best_savings, ik,kj,ij;
03933     Postsert postsert;
03934     Presert presert;
03935     VRPMove M;    
03936 
03937     if(k==VRPH_DEPOT)
03938         return false;
03939 
03940     i=VRPH_MAX(this->pred_array[k],VRPH_DEPOT);
03941     j=VRPH_MAX(this->next_array[k],VRPH_DEPOT);
03942 
03943 
03944     ik=this->d[i][k];
03945     kj=this->d[k][j];
03946     ij=this->d[i][j];
03947 
03948     best_savings=VRP_INFINITY;
03949     best_l=-1;
03950     best_m=-1;
03951 
03952     l=VRPH_DEPOT;
03953     m=abs(this->next_array[l]);
03954     while(m!=VRPH_DEPOT)
03955     {
03956         if(m>0)
03957         {
03958             savings = (ik + kj - this->d[l][m]) - alpha*(this->d[l][k] + this->d[k][m] - ij);
03959             if(savings<best_savings && l!=i && m!=j)
03960             {
03961                 // Now check feasibility-put k before m
03962                 if(presert.evaluate(this,k,m,&M)==true)
03963                 {
03964                     best_savings=savings;
03965                     best_l=l;
03966                     best_m=m;
03967                 }
03968             }
03969         }
03970         else
03971         {
03972             // m<0 indicating the end of a route
03973             
03974             // consider the edge l-VRPH_DEPOT
03975             mm=VRPH_DEPOT;
03976             savings = (ik + kj - this->d[l][mm]) - alpha*(this->d[l][k] + this->d[k][mm] - ij);
03977             if(savings<best_savings && l!=i && mm!=j)
03978             {
03979                 // put k after l
03980                 if(postsert.evaluate(this,k,l,&M)==true)
03981                 {
03982                     best_savings=savings;
03983                     best_l=l;
03984                     best_m=mm;
03985                 }
03986             }
03987 
03988             // Now consider the edge VRPH_DEPOT-abs(m)
03989             ll=VRPH_DEPOT;
03990             m=abs(m);
03991             savings = (ik + kj - this->d[ll][m]) - alpha*(this->d[ll][k] + this->d[k][m] - ij);
03992             if(savings<best_savings && ll!=i && m!=j)
03993             {
03994                 // put k before m
03995                 if(presert.evaluate(this,k,m,&M)==true)
03996                 {
03997                     best_savings=savings;
03998                     best_l=ll;
03999                     best_m=m;
04000                 }
04001             }
04002         }
04003         // Now advance to the next edge
04004         l=m;
04005         m=this->next_array[l];
04006     }
04007 
04008     // We now have the best location to move node k by 
04009     // placing in between best_l and best_m;
04010 
04011     if(best_savings==VRP_INFINITY)
04012     {
04013         // No locations found
04014         return false;
04015     }
04016 
04017     if(best_l!=VRPH_DEPOT)
04018     {
04019         if(postsert.move(this,k,best_l)==false)
04020             report_error("%s: postsert.move is false\n",__FUNCTION__);
04021     }
04022     else
04023     {
04024         if(presert.move(this,k,best_m)==false)
04025             report_error("%s: presert.move is false\n",__FUNCTION__);
04026     }
04027 
04028     return true;
04029 
04030 }
04031 
04032 int VRP::osman_perturb(int num, double alpha)
04033 {
04039 
04040     int tot=0;
04041     int k;
04042     int num_attempts=0;
04043 
04044     while(tot<num)
04045     {
04046         // Select a random node to move
04047         k=VRPH_MAX((int)(this->num_original_nodes * lcgrand(10)),1);//to avoid inserting VRPH_DEPOT
04048 
04049         if(routed[k])
04050         {
04051 
04052             // Try to insert
04053             if(this->osman_insert(k,alpha)==true)
04054                 tot++;
04055 
04056             num_attempts++;
04057             if(num_attempts>2*this->num_original_nodes)
04058                 return tot;// Only moved tot nodes
04059         }
04060     }
04061 
04062     // We were able to move num nodes
04063 
04064     return num;
04065 
04066     
04067 }
04068 
04069 bool VRP::check_fixed_edges(const char *message)
04070 {
04076 
04077 
04078 
04079     int i,j;
04080 
04081     for(i=0;i<=this->num_original_nodes;i++)
04082     {
04083         for(j=0;j<=this->num_original_nodes;j++)
04084         {
04085             if(this->fixed[i][j])
04086             {
04087                 // Make sure i-j or j-i exists
04088                 if(i!=VRPH_DEPOT)
04089                 {
04090                     if(VRPH_MAX(this->next_array[i],VRPH_DEPOT)!=j && VRPH_MAX(this->pred_array[i],VRPH_DEPOT)!=j)
04091                     {
04092                         fprintf(stderr,"Fixed edge %d-%d not in solution!!",i,j);
04093                         fprintf(stderr,"%d-%d-%d\n",VRPH_MAX(this->pred_array[i],VRPH_DEPOT),i,
04094                             VRPH_MAX(this->next_array[i],VRPH_DEPOT));
04095                         fprintf(stderr,"%s", message);
04096 
04097                         if(this->fixed[j][i])
04098                             fprintf(stderr,"%d-%d also fixed\n",j,i);
04099                         else
04100                             fprintf(stderr,"%d-%d NOT fixed!\n",j,i);
04101 
04102                         return false;
04103                     }
04104                 }
04105 
04106                 // Make sure i-j or j-i exists
04107                 if(j!=VRPH_DEPOT)
04108                 {
04109                     if(VRPH_MAX(this->next_array[j],VRPH_DEPOT)!=i && VRPH_MAX(this->pred_array[j],VRPH_DEPOT)!=i)
04110                     {
04111                         fprintf(stderr,"Fixed edge %d-%d not in solution!!",i,j);
04112                         fprintf(stderr,"%d-%d-%d\n",VRPH_MAX(this->pred_array[j],VRPH_DEPOT),j,
04113                             VRPH_MAX(this->next_array[j],VRPH_DEPOT));
04114                         fprintf(stderr,"%s", message);
04115                         if(this->fixed[j][i])
04116                             fprintf(stderr,"%d-%d also fixed\n",j,i);
04117                         else
04118                             fprintf(stderr,"%d-%d NOT fixed!\n",j,i);
04119                         return false;
04120 
04121                     }
04122                 }
04123             }
04124     
04125         }
04126 
04127     }
04128 
04129     return true;
04130 }
04131 int VRP::find_common_routes(int *sol1, int *sol2, int *route_nums)
04132 {
04138 
04139     int *h11, *h12, *h21, *h22; // to store the hash values
04140     double *L1, *L2;
04141     int i, j, cnt, r1, r2;    // to store the # of routes.
04142     VRPRoute rte(num_nodes);
04143 
04144     // First, import sol1
04145     this->import_solution_buff(sol1);
04146     
04147     r1=this->total_number_of_routes;
04148     h11=new int[r1+1]; h12=new int[r1+1];// 1-based arrays
04149     L1=new double[r1+1];
04150     for(i=1;i<=r1;i++)
04151     {
04152         // Import route i into rte
04153         this->update_route(i,&rte);
04154         // Hash the route
04155         h11[i]=rte.hash(SALT_1);
04156         h12[i]=rte.hash(SALT_2);
04157         L1[i]=rte.length;
04158     }
04159 
04160     // Now import sol2
04161     this->import_solution_buff(sol2);
04162     
04163     r2=this->total_number_of_routes;
04164     h21=new int[r2+1]; h22=new int[r2+1];// 1-based arrays
04165     L2=new double[r2+1];
04166     for(i=1;i<=r2;i++)
04167     {
04168         this->update_route(i,&rte);
04169         h21[i]=rte.hash(SALT_1);
04170         h22[i]=rte.hash(SALT_2);
04171         L2[i]=rte.length;
04172     }
04173     
04174 
04175 
04176     // Now compare the hash values
04177     cnt=0;
04178     for(i=1;i<=r1;i++)
04179     {
04180         for(j=1;j<=r2;j++)
04181         {
04182             if(h11[i]==h21[j] && h12[i]==h22[j] && VRPH_ABS(L1[i]-L2[j])<VRPH_EPSILON)
04183             {
04184                 // route i in sol1 is the same as route j in sol2
04185                 route_nums[cnt]=i;
04186                 cnt++;
04187             }
04188         }
04189     }
04190 
04191 
04192     delete [] h11; delete [] h12; delete [] h21; delete [] h22;
04193     delete [] L1; delete [] L2;
04194     return cnt;
04195 
04196 }
04197 
04198 void VRP::set_daily_demands(int day)
04199 {
04205 
04206     int i;
04207 
04208     if(day>0)
04209     {
04210         for(i=0;i<=this->num_original_nodes;i++)
04211         {
04212             if(this->nodes[i].daily_demands[day]>=0)
04213                 this->nodes[i].demand=this->nodes[i].daily_demands[day];
04214             else
04215                 this->nodes[i].demand=-1;
04216 
04217         }
04218     }
04219     else
04220     {
04221         // Day is 0 - use the mean demand
04222         for(i=0;i<=this->num_original_nodes;i++)
04223         {
04224             int mean_demand=0;
04225             int k=0;
04226             for(int j=1;j<=this->num_days;j++)
04227             {
04228                 if(this->nodes[i].daily_demands[j]>=0)
04229                 {
04230                     mean_demand+=this->nodes[i].daily_demands[j];
04231                     k++;
04232                 }
04233             }
04234             if(k>0)
04235                 this->nodes[i].demand=(int)((double)mean_demand/(double)k);
04236             else
04237                 this->nodes[i].demand=-1;
04238 
04239         }    
04240 
04241     }
04242 }
04243 
04244 void VRP::set_daily_service_times(int day)
04245 {
04250 
04251     int i;
04252 
04253     for(i=0;i<=this->num_original_nodes;i++)
04254     {
04255         this->nodes[i].service_time=this->nodes[i].daily_service_times[day];
04256     }
04257 
04258     // We also have to recompute the distance matrix if the service times are not identical
04259     // across the days
04260     this->create_distance_matrix(this->edge_weight_type);
04261 }
04262 
04263 
04264 void VRP::update_arrival_times()
04265 {
04269 
04270     int routenum, next, current;
04271     double t;
04272 
04273     // Set the arrival times to -1 so that the only ones with positive 
04274     // arrival_time values are those that are actually visited
04275     for(int i=1;i<=this->num_original_nodes;i++)
04276         this->nodes[i].arrival_time=-1;
04277 
04278 
04279     for(int i=1;i<=this->num_original_nodes;i++)
04280     {
04281 
04282         if(routed[i])
04283         {
04284             // Get the route number
04285             routenum=this->route_num[i];
04286 
04287             // Start at the beginning of this route
04288             current=this->route[routenum].start;
04289             t=this->d[VRPH_DEPOT][current];
04290 
04291             while(current!=i)
04292             {
04293                 next=VRPH_MAX(VRPH_DEPOT, this->next_array[current]);
04294                 t+=this->d[current][next];
04295                 current=next;
04296 
04297             }
04298 
04299             // Now subtract the service time at t as this is included in the distance 
04300             // matrix
04301             t-=this->nodes[i].service_time;
04302             this->nodes[i].arrival_time=t;
04303         }
04304     }
04305 
04306 }
04307 
04308 bool VRP::check_tabu_status(VRPMove *M, int *old_sol)
04309 {
04317 
04318 
04319 #if VRPH_TABU_DEBUG
04320     printf("Checking tabu status\n");
04321 #endif
04322 
04323 
04324     int i,j;
04325 
04326     // We will always accept a move that reduces the # of routes
04327     if(M->num_affected_routes>1)
04328     {
04329         for(i=0;i<M->num_affected_routes;i++)
04330         {
04331             if(M->route_custs[i]==0)
04332                 return true;
04333         }
04334     }    
04335     
04336     VRPRoute r(this->num_nodes);
04337     int num_tabu_routes=0;
04338 
04339 #if VRPH_TABU_DEBUG
04340     printf("Checking tabu status of %d routes\n",M->num_affected_routes);
04341 #endif
04342 
04343     for(i=0;i<M->num_affected_routes;i++)
04344     {
04345         this->update_route(M->route_nums[i],&r);
04346         r.hash_val = r.hash(SALT_1);
04347         r.hash_val2 = r.hash(SALT_2);    
04348         
04349         for(j=0;j<this->tabu_list->num_entries;j++)
04350         {
04351             if(r.hash_val==this->tabu_list->hash_vals1[j] && r.hash_val2==this->tabu_list->hash_vals2[j])
04352                 // The move is tabu!
04353                 num_tabu_routes++;            
04354         }
04355         
04356     }
04357     
04358 
04359     if(num_tabu_routes>0)
04360     {
04361         // At least one tabu route
04362         this->import_solution_buff(old_sol);
04363 #if VRPH_TABU_DEBUG
04364                 printf("Move is Tabu! Reverting to old solution\n");
04365 #endif
04366         return false;
04367     }
04368 
04369 #if VRPH_TABU_DEBUG
04370     printf("Move is not Tabu.  Updating list of %d routes\n",M->num_affected_routes);
04371 #endif
04372 
04373     // The move is not tabu - update the tabu list 
04374     for(i=0;i<M->num_affected_routes;i++)
04375     {    
04376         this->update_route(M->route_nums[i],&r);
04377 
04378         this->tabu_list->update_list(&r);
04379     }
04380 
04381     return true;
04382 
04383 }
04384 
04385 
04386 void VRP::print_stats()
04387 {
04392 
04393     printf("                       (Moves,      Evaluations)\n"); 
04394     printf("     One Point Move:   (%010d, %010d)\n",
04395         this->num_moves[ONE_POINT_MOVE_INDEX], this->num_evaluations[ONE_POINT_MOVE_INDEX]);
04396     printf("     Two Point Move:   (%010d, %010d)\n",
04397         this->num_moves[TWO_POINT_MOVE_INDEX], this->num_evaluations[TWO_POINT_MOVE_INDEX]);
04398     printf("   Three Point Move:   (%010d, %010d)\n",
04399         this->num_moves[THREE_POINT_MOVE_INDEX], this->num_evaluations[THREE_POINT_MOVE_INDEX]);
04400     printf("       Two-opt Move:   (%010d, %010d)\n",
04401         this->num_moves[TWO_OPT_INDEX], this->num_evaluations[TWO_OPT_INDEX]);
04402     printf("     Three-opt Move:   (%010d, %010d)\n",
04403         this->num_moves[THREE_OPT_INDEX], this->num_evaluations[THREE_OPT_INDEX]);
04404     printf("        Or-opt Move:   (%010d, %010d)\n",
04405         this->num_moves[OR_OPT_INDEX], this->num_evaluations[OR_OPT_INDEX]);
04406     printf("Cross-Exchange Move:   (%010d, %010d)\n\n",
04407         this->num_moves[CROSS_EXCHANGE_INDEX], this->num_evaluations[CROSS_EXCHANGE_INDEX]);
04408 
04409     return;
04410 
04411 }
04412 
04413 
04414 void VRP::reset()
04415 {
04419 
04420     this->solution_wh->liquidate();
04421     for(int j=1;j<=this->num_original_nodes;j++)
04422         this->routed[j]=false;
04423 }
04424