lk.c

Lin-Kernighan heuristic for the TSP and minimum weight perfect matching

case REP_ARRAY:printf("-r array ");break;
case REP_TWO_LEVEL:printf("-r two-level ");break;
case REP_TWO_LEVEL_DEBUG:printf("-r tld ");break;
case REP_SPLAY_0:printf("-r splay 0 ");break;
case REP_SPLAY_1:printf("-r splay 1 ");break;
case REP_SPLAY_2:printf("-r splay 2 ");break;
case REP_SPLAY_3:printf("-r splay 3 ");break;
default:errorif(1,"Bad representation == %d\n",representation);
}

if(should_sfc_reorder)printf("--sfc ");

if(max_generic_flips!=INT_MAX)printf("--maxdepth %d ",max_generic_flips);

if(lower_bound_name!=NULL)printf("-l %s %f ",lower_bound_name,lower_bound_value);
if(upper_bound_name!=NULL)printf("-u %s %f ",upper_bound_name,upper_bound_value);

if(filename){
if(filename[0]=='-')printf("-- %s",filename);
else printf("%s",filename);
}else printf("-");
printf("\n");
}

/*:58*/
#line 1617 "./lk.w"

printf("Start time: ");/*59:*/
#line 1731 "./lk.w"

{
#if HAVE_TIME_H && HAVE_TIME && HAVE_CTIME
time_t now= time(NULL);
printf("%s",ctime(&now));
fflush(stdout);
#endif
}

/*:59*/
#line 1618 "./lk.w"

}

/*:54*/
#line 693 "./lk.w"

/*61:*/
#line 1761 "./lk.w"

last_resource_mark= resource_mark("Reading the instance");
if(filename){
TSPLIB_in= fopen(filename,"r");
errorif(TSPLIB_in==NULL,"Couldn't open \"%s\" for reading",filename);
}else TSPLIB_in= stdin;
if(PostScript_filename){
FILE*prolog;
char buf[8192];
size_t countin,countout;
ps_out= fopen(PostScript_filename,"w");
errorif(ps_out==NULL,"Couldn't open \"%s\" for writing",filename);
prolog= fopen("prolog.ps","r");
errorif(prolog==NULL,"Couldn't open prolog.ps for reading");
while((countin= fread(buf,1,8192,prolog))>0){
char*p= buf;
while((countout= fwrite(p,1,countin,ps_out))<countin){
countin-= countout;p+= countout;
}
}
fclose(prolog);
}else ps_out= NULL;
tsp_instance= read_tsp_file(TSPLIB_in,ps_out,do_weighted_perfect_matching);
n= tsp_instance->n;

/*:61*//*117:*/
#line 2511 "./lk.w"

if(do_weighted_perfect_matching){
errorif(n%2,"Must have even number of vertices to have a perfect matching");
}

/*:117*/
#line 694 "./lk.w"

/*68:*/
#line 1855 "./lk.w"

if(should_sfc_reorder){
int i,n= tsp_instance->n;
coord_2d*new_coord;
errorif(tsp_instance->coord==NULL,
"Space filling curve reordering applies only to geometric instances.\n");
original_city_num= new_arr_of(int,n);
for(i= 0;i<n;i++)original_city_num[i]= i;
sort(original_city_num,(size_t)n,sizeof(int),cmp_sfc_Moore);
new_coord= new_arr_of(coord_2d,n);
for(i= 0;i<n;i++)new_coord[i]= tsp_instance->coord[original_city_num[i]];
#if 0
free_mem(tsp_instance->coord);mem_deduct(n*sizeof(coord_2d));
tsp_instance->coord= new_coord;
#else
for(i= 0;i<n;i++)tsp_instance->coord[i]= new_coord[i];
free_mem(new_coord);mem_deduct(n*sizeof(coord_2d));
#endif
#if defined(LK_SHOW_AFTER_SFC)
if(ps_out){
length_t len= cost(0,n-1);int i;for(i= 1;i<n;i++,len+= cost(i,i-1));
fprintf(ps_out,"(SFC(Moore) tour, len "length_t_spec") title\n",
length_t_pcast(len));
fprintf(ps_out,"(%s) comment\n",tsp_instance->comment);
for(i= 0;i<n;i++){
fprintf(ps_out,"%d %f %f sfcs\n",i+1,
tsp_instance->coord[i].x,tsp_instance->coord[i].y);
fprintf(ps_out,"(%d) %d label\n",i,i);
if(i>0)fprintf(ps_out,"%d %d edge\n",i-1,i);
}
fprintf(ps_out,"%d %d edge\n",n-1,0);
fprintf(ps_out,
"1 1 N { dup xs exch get exch ys exch get circle } for\n");
fprintf(ps_out,"showpage\n");
fflush(ps_out);
}
#endif
}

/*:68*/
#line 695 "./lk.w"

/*87:*/
#line 2152 "./lk.w"

if(ALWAYS_BUILD_DECLUSTER_STRUCTURES||mst_only){
mst= decluster_setup(n);
}

/*:87*//*100:*/
#line 2338 "./lk.w"

if(do_weighted_perfect_matching){
tour_setup= NULL;
}else{
switch(representation){
case REP_TWO_LEVEL:
{
const int n= tsp_instance->n;
twolevel_setup
(n,n<50?n:(n<1000?50:(n>100000?200:100)));
}
break;
case REP_TWO_LEVEL_DEBUG:
{
#if defined(TWOLEVEL_DEBUG)
const int n= tsp_instance->n;
twolevel_debug_setup
(n,n<50?n:(n<1000?50:(n>100000?200:100)));
#else
errorif(1,"Debugging of two-level isn't possible.  Recompile with -DTWOLEVEL_DEBUG");
#endif
}
break;
default:
tour_setup(tsp_instance->n);
}
}


/*:100*//*114:*/
#line 2487 "./lk.w"

if(!do_weighted_perfect_matching)jbmr_setup(n);

/*:114*//*118:*/
#line 2519 "./lk.w"

if(do_weighted_perfect_matching){
match_setup(n);
}

/*:118*/
#line 696 "./lk.w"

/*78:*/
#line 2028 "./lk.w"

if(E2_supports(tsp_instance)){
begin_data_structures_mark= 
last_resource_mark= resource_mark("Build the 2-d tree");
E2_create(tsp_instance);
if(ps_out){E2_postscript_show(ps_out);}
/*141:*/
#line 2815 "./lk.w"

#ifdef LK_CHECK_KDTREE
if(E2_supports(tsp_instance)){
int i,c,d,*done= new_arr_of(int,n),last;
length_t last_dist,next_dist;
for(c= 0;c<n;c++){
if(verbose)printf("%d ",c);fflush(stdout);
for(i= 0;i<n;i++){
done[i]= 0;
}
done[c]= 1;
last_dist= 0;
last= c;
for(i= 0;i<n-1;i++){
d= E2_nn(c);
if(0<=d&&d<n){
next_dist= cost(c,d);
errorif(last_dist>next_dist,
"NN for %d out of order: #%d %d cost="length_t_spec", "
"#%d %d cost="length_t_spec,
c,i-1,last,length_t_pcast(last_dist),
i,d,length_t_pcast(next_dist));

if(verbose>=500){
printf("nn(%d) #%d = %d dist "length_t_spec"\n",c,i,d,length_t_pcast(cost(c,d)));
}

last_dist= next_dist;
last= d;
}else{
printf("Invalid city %d returned on nn query %d from %d\n",
d,i,c);
}
errorif(done[d],"Shouldn't return %d on nn query at %d",d,c);
E2_hide(d);
done[d]= 1;
}
d= E2_nn(c);
errorif(d!=-1,"nn(%d) returned %d when all others hidden; should be -1.\n",
c,d);
E2_unhide_all();
}
free_mem(done);
if(verbose)printf("\nkd tree passed an integrity test on nearest-neighbour searching\n");
}
#endif



/*:141*/
#line 2034 "./lk.w"

}else{
begin_data_structures_mark
= resource_mark("Begin building data structures (but not k-d tree)");
}


/*:78*//*81:*/
#line 2061 "./lk.w"

if(ALWAYS_BUILD_DECLUSTER_STRUCTURES||mst_only||held_karp_only||held_karp_lambda_only){
last_resource_mark= resource_mark("Build a MST (decluster)");
mst_len= decluster_mst(tsp_instance,mst);
if(mst_only){
/*84:*/
#line 2112 "./lk.w"

if(ps_out&&tsp_instance->coord){
int i;
fprintf(ps_out,"(MST length "length_t_spec") title\n",
length_t_pcast(mst_len));
fprintf(ps_out,"(%s) comment\n",tsp_instance->comment);
for(i= 0;i<n-1;i++){
const int*city= mst->edge[i].city;
fprintf(ps_out,"%f %f %f %f ue\n",
tsp_instance->coord[city[0]].x[0],
tsp_instance->coord[city[0]].x[1],
tsp_instance->coord[city[1]].x[0],
tsp_instance->coord[city[1]].x[1]);
}
fprintf(ps_out,"showpage\ngrestore\n%%%%EOF");
fclose(ps_out);ps_out= NULL;
}


/*:84*/
#line 2066 "./lk.w"

/*82:*/
#line 2082 "./lk.w"

printf("c Minimum spanning tree generated by LK %s\n",VERSION_STRING);
printf("c TSPLIB: NAME: %s\n",tsp_instance->name);
printf("c TSPLIB: COMMENT: %s\n",tsp_instance->comment);
printf("c LK: version: %s\n",VERSION_STRING);
printf("c LK: length_t: %s\n",LENGTH_TYPE_STRING);
printf("c LK: option: -M\n");
if(noround)printf("c LK: option: --no-round\n");
printf("c LENGTH: "length_t_native_spec"\n",length_t_native_pcast(mst_len));
printf("p edge %d %d\n",mst->n+1,mst->n);
/*83:*/
#line 2099 "./lk.w"

{
int i;
const int m= mst->n;
for(i= 0;i<m;i++)
printf("e %d %d "length_t_native_spec"\n",
mst->edge[i].city[0]+1,
mst->edge[i].city[1]+1,
length_t_native_pcast(mst->edge[i].cost));
}

/*:83*/
#line 2092 "./lk.w"

exit(0);


/*:82*/
#line 2067 "./lk.w"

}
}


/*:81*//*85:*/
#line 2136 "./lk.w"

#if ALWAYS_BUILD_DECLUSTER_STRUCTURES
last_resource_mark= resource_mark("Preprocess the MST (decluster)");
decluster_preprocess(mst);
#endif

/*:85*//*90:*/
#line 2180 "./lk.w"

last_resource_mark= resource_mark("Build the adjacency structure");
errorif((candidate_expr&CAND_NN)&&(cand_nn_k<1),
"Neighbourhood bound must be positive, but is %d",cand_nn_k);
errorif((candidate_expr&CAND_NQ)&&(cand_nq_k<1),
"Neighbourhood quadrant bound must be positive, but is %d",cand_nq_k);
errorif(candidate_expr&CAND_DEL,"Candidate structure %d not supported",candidate_expr);
nn_build(
(candidate_expr&CAND_NN)?cand_nn_k:0,
(candidate_expr&CAND_NQ)?cand_nq_k:0,
(candidate_expr&CAND_DEL)?cand_del_d:0);

/*:90*/
#line 697 "./lk.w"

if(do_weighted_perfect_matching){
/*120:*/
#line 2531 "./lk.w"

last_resource_mark= resource_mark("Construct starting matching");
incumbent_len= match_construct(construction_algorithm,start_heuristic_param,random_seed);
if(verbose>=10)printf("Initial matching length: "length_t_spec"\n",
length_t_pcast(incumbent_len));
/*121:*/
#line 2540 "./lk.w"

if(ps_out&&tsp_instance->coord){
char heuristic_name[200];
switch(construction_algorithm){
case CONSTRUCT_CANONICAL:sprintf(heuristic_name,"Canonical");break;
case CONSTRUCT_GREEDY:sprintf(heuristic_name,"Deterministic Greedy");break;
case CONSTRUCT_GREEDY_RANDOM:sprintf(heuristic_name,"Randomized Greedy");break;
case CONSTRUCT_RANDOM:sprintf(heuristic_name,"Random %ld",start_heuristic_param);break;
default:sprintf(heuristic_name,"unknown--Bug!");
}
match_ps_out(ps_out,heuristic_name);
}



/*:121*/
#line 2536 "./lk.w"



/*:120*/
#line 699 "./lk.w"

}else{
/*102:*/
#line 2379 "./lk.w"

last_resource_mark= resource_mark("Construct starting tour");
tour= new_arr_of(int,n);
incumbent_len= 
construct(n,tour,construction_algorithm,start_heuristic_param,random_seed);
if(verbose>=10)printf("Initial tour length: "length_t_spec"\n",length_t_pcast(incumbent_len));
/*103:*/
#line 2389 "./lk.w"

if(ps_out&&tsp_instance->coord&&verbose>=50&&
!do_weighted_perfect_matching){
int i;
const char*tour_name;
switch(construction_algorithm){
case CONSTRUCT_CANONICAL:tour_name= "Canonical";break;
case CONSTRUCT_GREEDY_RANDOM:tour_name= "Randomized Greedy";break;
case CONSTRUCT_GREEDY:tour_name= "Deterministic Greedy";break;
case CONSTRUCT_RANDOM:tour_name= "Random";break;
default:errorif(1,"Unkown start tour parameter %d\n",construction_algorithm);
tour_name= "";
}
fprintf(ps_out,"(%s tour length "length_t_native_spec") title\n",
tour_name,
length_t_native_pcast(incumbent_len));
fprintf(ps_out,"(%s) comment\n",tsp_instance->comment);
for(i= 0;i<n;i++){
const int city= tour[i],next_city= tour[(i+1)%n];
fprintf(ps_out,"%f %f %f %f ue\n",
tsp_instance->coord[city].x[0],
tsp_instance->coord[city].x[1],
tsp_instance->coord[next_city].x[0],
tsp_instance->coord[next_city].x[1]);
}
fprintf(ps_out,"showpage\n");
}

/*:103*/
#line 2385 "./lk.w"


/*:102*//*110:*/
#line 2452 "./lk.w"

tour_set(tour);

/*:110*/
#line 701 "./lk.w"

/*66:*/
#line 1813 "./lk.w"

if(held_karp_only||held_karp_lambda_only){
length_t held_karp_bound;
ascend_setup(n);
held_karp_bound= ascend(n,incumbent_len);
printf("Held-Karp lower bound: "length_t_native_spec" ",
length_t_native_pcast(held_karp_bound));
if(upper_bound_value>0)
printf("(%.2f%% below upper bound %f) ",
100.0*((double)(upper_bound_value-held_karp_bound))/upper_bound_value,
upper_bound_value);
printf("(%.2f%% below incumbent "length_t_native_spec")\n",
100.0*((double)(incumbent_len-held_karp_bound))/incumbent_len,
length_t_native_pcast(incumbent_len));
if(held_karp_lambda_only){
int i;
double*const lambda= ascend_best_lambda();
errorif(lambda==NULL,"No Lagrange multilpier lambda vector!\n");
printf("Lagrange multipliers:\n");
for(i= 0;i<n;i++){
printf("%d %f\n",1+i,lambda[i]);
}
}
ascend_cleanup();
exit(0);
}

/*:66*/
#line 702 "./lk.w"

}
/*111:*/
#line 2462 "./lk.w"

{
prng_t*prng= prng_new(PRNG_DEFAULT,1998^random_seed);



last_resource_mark= resource_mark("Lin-Kernighan");
if(do_weighted_perfect_matching){
match_run(extra_backtrack?3:2,iterations,prng);
}else{
jbmr_run(iterations,prng);
}
prng_free(prng);
}

/*:111*/
#line 704 "./lk.w"

/*123:*/
#line 2564 "./lk.w"

last_resource_mark= resource_mark("The end");
if(verbose>=2){
int i;
if(verbose>=5)
for(i= 0;i<last_resource_mark;i++){
resource_report(stdout,i,i+1);
}
if(last_resource_mark>0){
resource_report(stdout,begin_data_structures_mark,last_resource_mark);
}
if(verbose>=10){
printf("End time: ");/*59:*/
#line 1731 "./lk.w"

{
#if HAVE_TIME_H && HAVE_TIME && HAVE_CTIME
time_t now= time(NULL);
printf("%s",ctime(&now));
fflush(stdout);
#endif
}

/*:59*/
#line 2576 "./lk.w"

mem_report();
}
/*125:*/
#line 2596 "./lk.w"

#if HAVE_GETHOSTNAME
{
#if !defined(MAXHOSTNAMELEN)
#define MAXHOSTNAMELEN 256
#endif 
char buf[MAXHOSTNAMELEN+1];
gethostname(buf,MAXHOSTNAMELEN);
buf[MAXHOSTNAMELEN]= 0;
printf("Machine: %s\n",buf);
}
#endif 


/*:125*/
#line 2579 "./lk.w"

}


/*:123*/
#line 705 "./lk.w"

/*135:*/
#line 2679 "./lk.w"

if(do_weighted_perfect_matching){
match_validate(&validate_len,&double_validate_len,&ordered_double_len,&raw_len);
}else{
int i,c,cn;
double*lengths= new_arr_of(double,n);
double*raw_lengths= new_arr_of(double,n);
length_t*length_t_lengths= new_arr_of(length_t,n);
validate_len= 0;
double_validate_len= ordered_double_len= raw_len= 0.0;
for(i= 0,c= 0;i<n;i++){
errorif(c==0&&i>0,"Not a tour");
cn= tour_next(c);
length_t_lengths[i]= cost(c,cn);
double_validate_len+= (double)length_t_lengths[i];
lengths[i]= (double)length_t_lengths[i];
if(tsp_instance->edge_weight_type==EUC_2D||
tsp_instance->edge_weight_type==CEIL_2D){
raw_lengths[i]= cost_from_euc2d_raw(c,cn);
}else{
raw_lengths[i]= lengths[i];
}
c= cn;
}
sort(lengths,(unsigned)n,sizeof(double),lk_double_cmp);
sort(raw_lengths,(unsigned)n,sizeof(double),lk_double_cmp);
sort(length_t_lengths,(unsigned)n,sizeof(length_t),lk_length_t_cmp);
for(i= 0;i<n;i++){
ordered_double_len+= lengths[i];
raw_len+= raw_lengths[i];
validate_len+= length_t_lengths[i];
}
free_mem(lengths);
free_mem(raw_lengths);
free_mem(length_t_lengths);
mem_deduct(n*(sizeof(double)+sizeof(double)+sizeof(length_t)));
errorif(c!=0,"Not a tour");
}

/*:135*//*136:*/
#line 2719 "./lk.w"

if(verbose>=2){
printf("Instance name: %s\n",tsp_instance->name);
printf("Instance comment: %s\n",tsp_instance->comment);
}
if(should_show_tour){
if(do_weighted_perfect_matching){
match_show(stdout);
}else{
int i,c;
printf("Tour:\n");
for(i= 0,c= 0;i<n;i++,c= tour_next(c)){
printf("%d ",(original_city_num?original_city_num[c]:c)+1);

if((i%10)==9||i==n-1)printf("\n");
}
}
}
if(verbose>=2)printf("Length: ");
if(verbose)printf(noround?"%f\n":"%.0f\n",(double)validate_len);
if(verbose>=10){
printf("\tincumbent_len       == "length_t_spec"\n"
"\tvalidate_len        == "length_t_spec"\n"
"\tdouble_validate_len == %f\n"
"\tordered_double_len  == %f\n"
"\traw_len             == %f\n"
"\tdiscrepancy == (incumbent_len-ordered_double_len) == %-10g\n",
length_t_pcast(incumbent_len),
length_t_pcast(validate_len),
double_validate_len,
ordered_double_len,
raw_len,
((double)incumbent_len)-ordered_double_len);
}
errorif(!within_epsilon(((double)incumbent_len),ordered_double_len),
"%s mistaken about improvement",
do_weighted_perfect_matching?"match_run()":"jbmr_run()");

/*:136*//*138:*/
#line 2764 "./lk.w"

if(ps_out){
if(do_weighted_perfect_matching){
match_ps_out(ps_out,(const char*)"LK-opt");
}else{
int i,c,cn;
fprintf(ps_out,"%%Here's the final tour\n");
fprintf(ps_out,"(LK opt, tour len "length_t_native_spec") title\n",
length_t_native_pcast(incumbent_len));
fprintf(ps_out,"(%s) comment\n",tsp_instance->comment);
for(i= 0,c= 0;i<n;i++,c= cn){
cn= tour_next(c);
errorif(c==0&&i>0,"Not a tour");
fprintf(ps_out,"%f x %f y %f x %f y rawedge\n",
tsp_instance->coord[cn].x[0],
tsp_instance->coord[cn].x[1],
tsp_instance->coord[c].x[0],
tsp_instance->coord[c].x[1]);
}
fprintf(ps_out,"showpage\n");fflush(ps_out);
}
fprintf(ps_out,"end\ngrestore\n%%EOF\n");fclose(ps_out);
}


/*:138*/
#line 706 "./lk.w"

/*71:*/
#line 1907 "./lk.w"

free_mem(original_city_num);

/*:71*//*126:*/
#line 2613 "./lk.w"

lk_cleanup();

/*:126*/
#line 707 "./lk.w"

return 0;
}

/*:6*//*139:*/
#line 2790 "./lk.w"

int
lk_double_cmp(const void*a,const void*b)
{
const double da= *((const double*)a),db= *((const double*)b);
if(da<db)return-1;
if(da>db)return+1;
return 0;
}

int
lk_length_t_cmp(const void*a,const void*b)
{
const length_t da= *((const length_t*)a),db= *((const length_t*)b);
if(da<db)return-1;
if(da>db)return+1;
return 0;
}

/*:139*/
#line 635 "./lk.w"

const char*lk_rcs_id= "$Id: lk.w,v 1.275 1999/02/19 17:57:26 neto Exp neto $";


/*:2*/