tspgen.c

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

{int i;
for(i= 0;i<n;i++){
hull_vertex[i].x= 0;
hull_vertex[i].y= 0;
hull_vertex[i].next= hull_vertex[i].prev= hull_vertex+i;
}
}

/*:47*/
#line 967 "./tspgen.w"


/*:46*//*55:*/
#line 1087 "./tspgen.w"

component= new_arr_of(component_t,n+ne);
{int i;
for(i= 0;i<n;i++){
component[i].hull= hull_vertex+i;
component[i].hull_size= 1;
transform_init(&component[i].transform);
}
for(i= 0;i<n+ne;i++)
component[i].left= component[i].right= component[i].ancestor= NULL;
}

/*:55*//*62:*/
#line 1172 "./tspgen.w"

hull_sortbox= new_arr_of(hull_sortbox_t,n);

/*:62*//*87:*/
#line 1574 "./tspgen.w"

hull_buf= new_arr_of(point_t,n);

/*:87*//*96:*/
#line 1686 "./tspgen.w"

stack= new_arr_of(int,n+1);

/*:96*/
#line 262 "./tspgen.w"

/*31:*/
#line 748 "./tspgen.w"

for(nc= n;nc>1;nc--){
component_t*a,*b;
const edge_t this_edge= in_tree->edge[n-nc];
const double this_length= this_edge.cost;
if(try_forcing_isomorphic){
const int u= this_edge.city[0],v= this_edge.city[1];
/*40:*/
#line 859 "./tspgen.w"

{
component_t*d;
d= component+u;/*41:*/
#line 874 "./tspgen.w"

{
component_t*updating= d;
while(d->ancestor)d= d->ancestor;
for(;updating!=d;updating= updating->ancestor)
updating->ancestor= d;
}


/*:41*/
#line 862 "./tspgen.w"
a= d;
d= component+v;/*41:*/
#line 874 "./tspgen.w"

{
component_t*updating= d;
while(d->ancestor)d= d->ancestor;
for(;updating!=d;updating= updating->ancestor)
updating->ancestor= d;
}


/*:41*/
#line 863 "./tspgen.w"
b= d;
errorif(a==b,
"Input wasn't a tree.  It contained a cycle involving edge (%d,%d)",
u+1,v+1);
}

/*:40*/
#line 755 "./tspgen.w"

/*42:*/
#line 887 "./tspgen.w"

{component_t*c;int w;
c= a;w= u;/*43:*/
#line 903 "./tspgen.w"

{
int found_w= 0;
hull_vertex_t*first= c->hull,*here= first;
do{
if(hull_city(here)==w){
c->hull= here;
found_w= 1;
/*115:*/
#line 1966 "./tspgen.w"

#if TSPGEN_DEBUG >= 200
fprintf(stderr,"Found %d and set c->hull = %p (city %d)\n",
w,c->hull,hull_city(c->hull));
#endif

/*:115*/
#line 911 "./tspgen.w"

break;
}
here= here->next;
}while(here!=first);
if(!found_w){
const double old_join_length_bias= join_length_bias;
join_length_bias= 1.0;
/*61:*/
#line 1158 "./tspgen.w"

{hull_vertex_t*here= c->hull;
int i,index;
for(i= 0,here= c->hull;i<c->hull_size;i++,here= here->next){
hull_sortbox[i].vertex= here;
hull_sortbox[i].len= hull_edge_len(here);
}
dsort(hull_sortbox,(size_t)c->hull_size,sizeof(hull_sortbox_t),cmp_hull_sortbox);
index= floor(scale_01(join_length_bias,0,c->hull_size));
if(index==c->hull_size)index--;
c->hull= hull_sortbox[index].vertex;
}

/*:61*/
#line 919 "./tspgen.w"

join_length_bias= old_join_length_bias;
force_failures++;
/*116:*/
#line 1973 "./tspgen.w"

#if TSPGEN_DEBUG >= 200
fprintf(stderr,"Force failed to find city %d on hull\n",w);
#endif

/*:116*/
#line 922 "./tspgen.w"

}
}

/*:43*/
#line 889 "./tspgen.w"

c= b;w= v;/*43:*/
#line 903 "./tspgen.w"

{
int found_w= 0;
hull_vertex_t*first= c->hull,*here= first;
do{
if(hull_city(here)==w){
c->hull= here;
found_w= 1;
/*115:*/
#line 1966 "./tspgen.w"

#if TSPGEN_DEBUG >= 200
fprintf(stderr,"Found %d and set c->hull = %p (city %d)\n",
w,c->hull,hull_city(c->hull));
#endif

/*:115*/
#line 911 "./tspgen.w"

break;
}
here= here->next;
}while(here!=first);
if(!found_w){
const double old_join_length_bias= join_length_bias;
join_length_bias= 1.0;
/*61:*/
#line 1158 "./tspgen.w"

{hull_vertex_t*here= c->hull;
int i,index;
for(i= 0,here= c->hull;i<c->hull_size;i++,here= here->next){
hull_sortbox[i].vertex= here;
hull_sortbox[i].len= hull_edge_len(here);
}
dsort(hull_sortbox,(size_t)c->hull_size,sizeof(hull_sortbox_t),cmp_hull_sortbox);
index= floor(scale_01(join_length_bias,0,c->hull_size));
if(index==c->hull_size)index--;
c->hull= hull_sortbox[index].vertex;
}

/*:61*/
#line 919 "./tspgen.w"

join_length_bias= old_join_length_bias;
force_failures++;
/*116:*/
#line 1973 "./tspgen.w"

#if TSPGEN_DEBUG >= 200
fprintf(stderr,"Force failed to find city %d on hull\n",w);
#endif

/*:116*/
#line 922 "./tspgen.w"

}
}

/*:43*/
#line 890 "./tspgen.w"

}

/*:42*/
#line 756 "./tspgen.w"

}else{
/*35:*/
#line 819 "./tspgen.w"

{component_t temp;
int base= 2*(n-nc),r;
r= base+gb_prng_unif_rand(component_prng,nc);
temp= component[base];component[base]= component[r];component[r]= temp;
r= base+1+gb_prng_unif_rand(component_prng,nc-1);
temp= component[base+1];component[base+1]= component[r];component[r]= temp;
a= component+base;
b= component+base+1;
}

/*:35*/
#line 758 "./tspgen.w"

/*58:*/
#line 1115 "./tspgen.w"

if(join_length_bias<0){
/*59:*/
#line 1127 "./tspgen.w"

{long steps;
for(steps= gb_prng_unif_rand(hull_prng,a->hull_size);steps;steps--)
a->hull= a->hull->next;
for(steps= gb_prng_unif_rand(hull_prng,b->hull_size);steps;steps--)
b->hull= b->hull->next;
}

/*:59*/
#line 1117 "./tspgen.w"

}else{
/*60:*/
#line 1136 "./tspgen.w"

{component_t*c;
c= a;/*61:*/
#line 1158 "./tspgen.w"

{hull_vertex_t*here= c->hull;
int i,index;
for(i= 0,here= c->hull;i<c->hull_size;i++,here= here->next){
hull_sortbox[i].vertex= here;
hull_sortbox[i].len= hull_edge_len(here);
}
dsort(hull_sortbox,(size_t)c->hull_size,sizeof(hull_sortbox_t),cmp_hull_sortbox);
index= floor(scale_01(join_length_bias,0,c->hull_size));
if(index==c->hull_size)index--;
c->hull= hull_sortbox[index].vertex;
}

/*:61*/
#line 1138 "./tspgen.w"

c= b;/*61:*/
#line 1158 "./tspgen.w"

{hull_vertex_t*here= c->hull;
int i,index;
for(i= 0,here= c->hull;i<c->hull_size;i++,here= here->next){
hull_sortbox[i].vertex= here;
hull_sortbox[i].len= hull_edge_len(here);
}
dsort(hull_sortbox,(size_t)c->hull_size,sizeof(hull_sortbox_t),cmp_hull_sortbox);
index= floor(scale_01(join_length_bias,0,c->hull_size));
if(index==c->hull_size)index--;
c->hull= hull_sortbox[index].vertex;
}

/*:61*/
#line 1139 "./tspgen.w"

}

/*:60*/
#line 1119 "./tspgen.w"

}

/*:58*/
#line 759 "./tspgen.w"

}
/*75:*/
#line 1333 "./tspgen.w"

/*76:*/
#line 1364 "./tspgen.w"

rotate_hull_with_bias(a,LAY_FORWARD,packing_factor);
if(try_forcing_isomorphic){
const int to_prev_is_longer= 
hull_edge_len(b->hull->prev)>hull_edge_len(b->hull);
const int join_longer= join_length_bias>=0.5;
if(join_longer==to_prev_is_longer){
rotate_hull_with_bias(b,LAY_BACKWARD,packing_factor);
}else{
rotate_hull_with_bias(b,LAY_FORWARD,packing_factor);
}
}else{
b->hull= b->hull->next;
rotate_hull_with_bias(b,LAY_BACKWARD,packing_factor);
}
component_then_rotate(b,pi);
/*106:*/
#line 1849 "./tspgen.w"

#if TSPGEN_DEBUG >=120
{component_t*comp;
int is_forward= -9999;
fprintf(stderr," Component a\n");fflush(stderr);
comp= a;/*105:*/
#line 1835 "./tspgen.w"

#if TSPGEN_DEBUG >= 125
{int i;hull_vertex_t*here;
fprintf(stderr,"  Hull size %d (city %d), `bot-hull' point first.  is_forward=%d\n",
comp->hull_size,hull_city(comp->hull),is_forward);
for(i= 0,here= comp->hull;i<comp->hull_size;i++,here= here->next){
fprintf(stderr,"\t%10.5f  %10.5f (city %d)\n",here->x,here->y,
hull_city(here));
}
errorif(here!=comp->hull,"hull_size wrong");
}
#endif

/*:105*/
#line 1854 "./tspgen.w"

fprintf(stderr," Component b\n");fflush(stderr);
comp= b;/*105:*/
#line 1835 "./tspgen.w"

#if TSPGEN_DEBUG >= 125
{int i;hull_vertex_t*here;
fprintf(stderr,"  Hull size %d (city %d), `bot-hull' point first.  is_forward=%d\n",
comp->hull_size,hull_city(comp->hull),is_forward);
for(i= 0,here= comp->hull;i<comp->hull_size;i++,here= here->next){
fprintf(stderr,"\t%10.5f  %10.5f (city %d)\n",here->x,here->y,
hull_city(here));
}
errorif(here!=comp->hull,"hull_size wrong");
}
#endif

/*:105*/
#line 1856 "./tspgen.w"

}
#endif

/*:106*/
#line 1380 "./tspgen.w"


/*:76*/
#line 1334 "./tspgen.w"

{
component_t*c= next_component_slot();
hull_vertex_t*ah= a->hull,*bh= b->hull;
/*77:*/
#line 1383 "./tspgen.w"

component_then_translate(a,bh->x-ah->x,bh->y-ah->y+this_length);
/*108:*/
#line 1872 "./tspgen.w"

#if TSPGEN_DEBUG >= 100
{
hull_vertex_t*here,*ahere,*bhere;
const double d= euc2hull(ah,bh);
int fail= 0;
here= ah;


for(here= bh->next;here!=bh;here= here->next){
const double hd= euc2hull(ah,here);
if(d>hd)fail= 1,fprintf(stderr,"d = %f > hd = %f\n",d,hd);
fprintf(stderr,"a");
}

for(here= ah->next;here!=ah;here= here->next){
const double hd= euc2hull(bh,here);
if(d>hd)fail= 1,fprintf(stderr,"d = %f > hd = %f\n",d,hd);
fprintf(stderr,"b");
}


ahere= ah;
do{
bhere= bh;
do{
const double hd= euc2hull(ahere,bhere);
if((ahere!=ah||bhere!=bh)&&d>hd)
fail= 1,fprintf(stderr,"full d = %f > hd = %f\n",d,hd);
bhere= bhere->next;
fprintf(stderr,"c");
}while(bhere!=bh);
ahere= ahere->next;
}while(ahere!=ah);



if(d<this_length-1e-6)
fail= 1,fprintf(stderr,"d = %f < %f = this_length\n",d,this_length);
if(d>this_length+1e-6)
fail= 1,fprintf(stderr,"d = %f > %f = this_length\n",d,this_length);
fprintf(stderr,"\n");
errorif(fail,"distances mucked up!");
}
#endif

/*:108*/
#line 1385 "./tspgen.w"


/*:77*/
#line 1338 "./tspgen.w"

/*83:*/
#line 1508 "./tspgen.w"

c->left= a;
c->right= b;
a->ancestor= b->ancestor= c;
/*85:*/
#line 1537 "./tspgen.w"

/*86:*/
#line 1551 "./tspgen.w"

{hull_vertex_t*start,*here;
nh= 0;
here= start= a->hull;
do{
hull_buf[nh].x= here->x;
hull_buf[nh].y= here->y;
hull_buf[nh].home= here;
here= here->next;
nh++;
}while(here!=start);
here= start= b->hull;
do{
hull_buf[nh].x= here->x;
hull_buf[nh].y= here->y;
hull_buf[nh].home= here;
here= here->next;
nh++;
}while(here!=start);
}

/*:86*/
#line 1538 "./tspgen.w"

/*91:*/
#line 1593 "./tspgen.w"

{int i;
west= hull_buf[0];
for(i= 1;i<nh;i++)
if(hull_buf[i].x<west.x){
west= hull_buf[i];
hull_buf[i]= hull_buf[0];
hull_buf[0]= west;
}
}

/*:91*/
#line 1539 "./tspgen.w"

/*92:*/
#line 1610 "./tspgen.w"

{int i;
for(i= 1;i<nh;i++){
hull_buf[i].angle= atan2(hull_buf[i].y-west.y,hull_buf[i].x-west.x);
}
dsort(&hull_buf[1],(size_t)nh-1,sizeof(point_t),cmp_point);
}

/*:92*/
#line 1540 "./tspgen.w"

/*94:*/
#line 1660 "./tspgen.w"

{
int i;
int top= -1;
/*111:*/
#line 1941 "./tspgen.w"

#if TSPGEN_DEBUG >= 400
fprintf(stderr,"Graham scan: nh==%d\n",nh);fflush(stderr);
#endif

/*:111*/
#line 1664 "./tspgen.w"

for(i= 0;i<=nh;i++){
stack[++top]= (i<nh)?i:0;
/*112:*/
#line 1947 "./tspgen.w"

#if TSPGEN_DEBUG >= 500
fprintf(stderr," push   %d, top==%d\n",stack[top],top);fflush(stderr);
#endif


/*:112*/
#line 1667 "./tspgen.w"

while(top>=2&&stack[top-2]!=stack[top]
&&!left_turn(hull_buf[stack[top-2]],
hull_buf[stack[top-1]],
hull_buf[stack[top]])){
/*113:*/
#line 1954 "./tspgen.w"

#if TSPGEN_DEBUG>=500
fprintf(stderr," remove %d, top==%d\n",stack[top-1],top-1);fflush(stderr);
#endif

/*:113*/
#line 1672 "./tspgen.w"

stack[top-1]= stack[top];
top--;
}
}
/*114:*/
#line 1960 "./tspgen.w"

#if TSPGEN_DEBUG >=400
fprintf(stderr,"Done Graham scan: top==%d\n",top);fflush(stdout);
#endif

/*:114*/
#line 1677 "./tspgen.w"

/*98:*/
#line 1702 "./tspgen.w"

for(i= 0;i<top;i++){
hull_vertex_t*here= hull_buf[stack[i]].home;
here->x= hull_buf[stack[i]].x;
here->y= hull_buf[stack[i]].y;
here->prev= hull_buf[stack[(i-1+top)%top]].home;
here->next= hull_buf[stack[(i+1)%top]].home;
}
c->hull= hull_buf[stack[0]].home;
c->hull_size= top;

/*:98*/
#line 1678 "./tspgen.w"

}

/*:94*/
#line 1541 "./tspgen.w"


/*:85*/
#line 1512 "./tspgen.w"


/*:83*/
#line 1339 "./tspgen.w"

/*107:*/
#line 1861 "./tspgen.w"

#if TSPGEN_DEBUG >=100
fprintf(stderr,"New component\n");fflush(stderr);
{component_t*comp= c;
int is_forward= -9999;
/*105:*/
#line 1835 "./tspgen.w"

#if TSPGEN_DEBUG >= 125
{int i;hull_vertex_t*here;
fprintf(stderr,"  Hull size %d (city %d), `bot-hull' point first.  is_forward=%d\n",
comp->hull_size,hull_city(comp->hull),is_forward);
for(i= 0,here= comp->hull;i<comp->hull_size;i++,here= here->next){
fprintf(stderr,"\t%10.5f  %10.5f (city %d)\n",here->x,here->y,
hull_city(here));
}
errorif(here!=comp->hull,"hull_size wrong");
}
#endif

/*:105*/
#line 1866 "./tspgen.w"

}
fprintf(stderr,"\n\n\n");fflush(stderr);
#endif

/*:107*/
#line 1340 "./tspgen.w"

if(nc>2){/*78:*/
#line 1397 "./tspgen.w"

component_then_rotate(c,uniform_sample(angle_prng,-pi,pi));

/*:78*/
#line 1341 "./tspgen.w"
}
}

/*:75*/
#line 761 "./tspgen.w"

}

/*:31*/
#line 263 "./tspgen.w"

/*100:*/
#line 1751 "./tspgen.w"

{char*name= NULL;
if(!name)name= option_name;
if(!name)name= in_name;
if(!name){
char constructed_name[1000];
sprintf(constructed_name,"%s%d",prog_name,n);
name= dup_string(constructed_name);
}
printf("NAME: %s\n",name);
}


/*:100*//*102:*/
#line 1776 "./tspgen.w"

if(in_comment){
char*newline;
for(newline= strchr(in_comment,'\n');newline;newline= strchr(in_comment,'\n'))*newline= ' ';
}
printf("COMMENT: %s | (LK %s) %s ",
in_comment?in_comment:"<unknown source>",
VERSION_STRING,prog_name);
if(try_forcing_isomorphic){
printf("-i (%d or %.1f%% fail) ",force_failures,(100.0*force_failures)/ne);
}
printf("-s %ld -j %g -p %g\n",seed,join_length_bias,packing_factor);

/*:102*//*103:*/
#line 1799 "./tspgen.w"

/*47:*/
#line 973 "./tspgen.w"

{int i;
for(i= 0;i<n;i++){
hull_vertex[i].x= 0;
hull_vertex[i].y= 0;
hull_vertex[i].next= hull_vertex[i].prev= hull_vertex+i;
}
}

/*:47*/
#line 1800 "./tspgen.w"

printf("TYPE: TSP\n");
printf("DIMENSION: %d\n",n);
printf("EDGE_WEIGHT_TYPE: EUC_2D\n");
printf("NODE_COORD_SECTION\n");
errorif(nc!=1,"Must be only one component, but nc==%d.",nc);
{int city_num= 1;
write_component(stdout,next_component_slot()-1,&city_num);
errorif(city_num-1!=n,"Wrote %d cities instead of %d.",city_num,n);
}
printf("EOF\n");

/*:103*/
#line 264 "./tspgen.w"

/*17:*/
#line 565 "./tspgen.w"

if(in_tree){
free_mem(in_tree->edge);mem_deduct(in_tree->ne*sizeof(edge_t));
in_tree->ne= -1;
in_tree->num_filled= 0;
free_mem(in_tree);
}

/*:17*//*23:*/
#line 630 "./tspgen.w"

if(in_comment){mem_deduct(strlen(in_comment));free_mem(in_comment);}
if(in_name){mem_deduct(strlen(in_name));free_mem(in_name);}


/*:23*//*48:*/
#line 984 "./tspgen.w"

free_mem(hull_vertex);

/*:48*//*56:*/
#line 1100 "./tspgen.w"

free_mem(component);

/*:56*//*63:*/
#line 1176 "./tspgen.w"

free_mem(hull_sortbox);

/*:63*//*88:*/
#line 1578 "./tspgen.w"

free_mem(hull_buf);

/*:88*//*97:*/
#line 1690 "./tspgen.w"

free_mem(stack);


/*:97*//*101:*/
#line 1765 "./tspgen.w"

if(name){mem_deduct(strlen(name));free_mem(name);}


/*:101*/
#line 265 "./tspgen.w"

return 0;
}

/*:5*/