tspgen.c

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

#define version(out) (fprintf(out,"%s (LK %s)\n",prog_name,VERSION_STRING) )  \

#define MAX_LINE_LENGTH (20000)  \

#define TSPLIB_NAME_PREFIX "c TSPLIB: NAME: "
#define TSPLIB_COMMENT_PREFIX "c TSPLIB: COMMENT: "
#define LK_LENGTH_T_PREFIX "c LK: length_t: " \

#define next_component_slot() (component+n+n-nc)  \

#define hull_city(hull_vertex_ptr) ((hull_vertex_ptr) -hull_vertex)  \

#define sqr(X) ((X) *(X) ) 
#define hull_edge_len(HP) (sqr((HP) ->x-(HP) ->next->x) +sqr((HP) ->y-(HP) ->next->y) )  \

#define LAY_FORWARD 1
#define LAY_BACKWARD 0 \

#define pi (3.1415926535897932384626433832795028841972) 
#define two_to_the_31 ((unsigned long) 0x80000000) 
#define scale_01(s,a,b) ((a) +((b) -(a) ) *(s) ) 
#define uniform_sample(g,a,b) scale_01((((double) gb_prng_next_rand(g) ) /two_to_the_31) ,(a) ,(b) )  \

#define uniform_sample_incl(g,a,b) scale_01(((double) gb_prng_next_rand(g) ) /(two_to_the_31-1) ,(a) ,(b) )  \
 \

/*5:*/
#line 234 "./tspgen.w"


const char*prog_name= "tspgen";
const char*tspgen_rcs_id= "$Id: tspgen.w,v 1.31 1998/12/05 22:36:13 neto Exp neto $";
#include <config.h>
#include <lkconfig.h>
/*7:*/
#line 363 "./tspgen.w"

#define _POSIX_C_SOURCE 2   
#if HAVE_UNISTD_H
#include <unistd.h>
#endif

/*:7*//*19:*/
#line 578 "./tspgen.w"

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <string.h>



/*:19*//*72:*/
#line 1283 "./tspgen.w"

#include <math.h>


/*:72*//*80:*/
#line 1448 "./tspgen.w"

#include <errno.h>

/*:80*/
#line 240 "./tspgen.w"

/*8:*/
#line 370 "./tspgen.w"

#define FIXINCLUDES_NEED_GETOPT
#include "fixincludes.h"
#undef FIXINCLUDES_NEED_GETOPT


/*:8*//*20:*/
#line 588 "./tspgen.w"

#include "error.h"
#include "memory.h"
#include "length.h"

/*:20*//*28:*/
#line 697 "./tspgen.w"

#include "dsort.h"

/*:28*//*39:*/
#line 851 "./tspgen.w"

#include "gb_flip.h"

/*:39*/
#line 241 "./tspgen.w"

/*67:*/
#line 1231 "./tspgen.w"

typedef struct{
double theta;
double v[2];
}rigid_motion_t;

/*:67*/
#line 242 "./tspgen.w"

/*45:*/
#line 951 "./tspgen.w"


typedef struct hull_vertex_s{
double x,y;
struct hull_vertex_s*next,*prev;
}hull_vertex_t;

/*:45*//*54:*/
#line 1075 "./tspgen.w"

typedef struct component_s{
struct component_s*left,*right,*ancestor;
hull_vertex_t*hull;
int hull_size;
rigid_motion_t transform;
}component_t;

/*:54*//*65:*/
#line 1184 "./tspgen.w"

typedef struct{
hull_vertex_t*vertex;
double len;
}hull_sortbox_t;

/*:65*//*90:*/
#line 1588 "./tspgen.w"

typedef struct{double x,y,angle;hull_vertex_t*home;}point_t;

/*:90*/
#line 243 "./tspgen.w"


#if !defined(TSPGEN_DEBUG)
#define TSPGEN_DEBUG 0
#endif

/*15:*/
#line 545 "./tspgen.w"

typedef struct{
int city[2];
length_t cost;
}edge_t;

typedef struct{
edge_t*edge;
int ne;
int num_filled;
}tree_t;

/*:15*/
#line 249 "./tspgen.w"

/*13:*/
#line 502 "./tspgen.w"

static long seed;
static double packing_factor,join_length_bias;
static int try_forcing_isomorphic,force_failures= 0;
static char*option_name= NULL;
static char*name= NULL;

/*:13*//*18:*/
#line 574 "./tspgen.w"

static tree_t*in_tree= NULL;

/*:18*//*22:*/
#line 625 "./tspgen.w"

static char*in_comment= NULL;
static char*in_name= NULL;

/*:22*//*25:*/
#line 656 "./tspgen.w"

static int n,ne;


/*:25*//*34:*/
#line 807 "./tspgen.w"

static int nc;

/*:34*//*38:*/
#line 846 "./tspgen.w"

static gb_prng_t*angle_prng;

/*:38*//*49:*/
#line 988 "./tspgen.w"

static hull_vertex_t*hull_vertex;

/*:49*//*57:*/
#line 1104 "./tspgen.w"

static component_t*component;


/*:57*//*64:*/
#line 1180 "./tspgen.w"

static hull_sortbox_t*hull_sortbox;

/*:64*//*89:*/
#line 1582 "./tspgen.w"

static point_t*hull_buf= NULL,west;
static int nh= 0;

/*:89*/
#line 250 "./tspgen.w"

/*11:*/
#line 452 "./tspgen.w"

static void usage(FILE*out,char**argv);
static void sample(FILE*out);


/*:11*/
#line 251 "./tspgen.w"

/*10:*/
#line 411 "./tspgen.w"

static void
usage(FILE*out,char**argv)
{
version(out);
fprintf(out,
"Generate TSPLIB instances from minimum spanning trees\n"
"\nCopyright (C) 1997, 1998 David M. Neto\n"
"LK comes with NO WARRANTY, to the extent permitted by law.\n"
"You may redistribute and/or modify copies of LK under the terms of the\n"
"GNU General Public License, version 2 or any later version.\n"
"For more information about these matters, see the file named COPYING.\n\n"
"Usage: %s [options] < <weighted-tree.edge> \n",
argv[0]);
fprintf(out,
" -h                    : Output this help and quit\n"
" --help                : Output this help and quit\n"
" -i                    : Try hard to make the output have a minimum spanning\n"
"                         tree that is isomorphic to the input weighted tree\n"
" -j <join-length-bias> : Bias choice of hull point based on associated\n"
"                         hull edge length.\n"
"                           negative means random vertex\n"
"                           0 to 1: use this as a fractional ordinal in\n"
"                             sorted hull edge lengths (read the source!)\n"
" -n <name>             : Use <name> as the instance name, overriding \n"
"                         the 'c TSPLIB: NAME: <name>', if any \n"
" -p <packing-factor>   : Non-negative <packing-factor> is a hint on how\n"
"                         tight to try to make the instance:\n"
"                           0 means try hard to stretch it loose,\n"
"                           large values means try to pack it tight\n"
"                           1 means join angles are uniform over legal range\n"
" -s <seed>             : Use integer <seed> as the random number seed\n"
" --version             : Print version number and exit successfully\n"
"\nThe weighted tree file is in DIMACS matching input format in 'edge' style.\n"
"SAMPLE BEGINS NEXT LINE\n");
sample(out);
fprintf(out,"SAMPLE ENDS PREVIOUS LINE\n");
}

/*:10*//*12:*/
#line 470 "./tspgen.w"

static void
sample(FILE*out)
{
fprintf(out,
"c Minimum spanning tree generated by LK 0.4.9\n"
"c TSPLIB: NAME: foobar5\n"
"c TSPLIB: COMMENT: Sample bogus instance for tspgen.\n"
"c LK: option: -M\n"
"c LK: length_t: double\n"
"c LENGTH: 21876.000000000000000\n"
"p edge 5 4\n"
"e 1 2 4328.000000000000000\n"
"e 3 1 10055.000000000000000\n"
"e 1 4 5760.000000000000000\n"
"e 5 2 1733.000000000000000\n");
}



/*:12*//*29:*/
#line 701 "./tspgen.w"

static int
cmp_edge(const void*ap,const void*bp)
{
const edge_t*a= (const edge_t*)ap,*b= (const edge_t*)bp;
if(a->cost<b->cost)return-1;
if(a->cost>b->cost)return 1;
return 0;
}


/*:29*//*50:*/
#line 995 "./tspgen.w"

static void
hull_then_rotate(hull_vertex_t*h,double theta)
{
const double ct= cos(theta),st= sin(theta);
const hull_vertex_t*start= h;
do{
const double x= h->x,y= h->y;
h->x= x*ct-y*st;
h->y= x*st+y*ct;
h= h->next;
}while(h!=start);
}

static void
hull_then_translate(hull_vertex_t*h,double x,double y)
{
const hull_vertex_t*start= h;
do{
h->x+= x;
h->y+= y;
h= h->next;
}while(h!=start);
}

/*:50*//*66:*/
#line 1191 "./tspgen.w"

static int
cmp_hull_sortbox(const void*ap,const void*bp)
{
const double
al= ((const hull_sortbox_t*)ap)->len,
bl= ((const hull_sortbox_t*)bp)->len;
if(al<bl)return-1;
if(al>bl)return 1;
return 0;
}

/*:66*//*68:*/
#line 1241 "./tspgen.w"

static void
transform_init(rigid_motion_t*t)
{
t->theta= t->v[0]= t->v[1]= 0.0;
}

/*:68*//*70:*/
#line 1260 "./tspgen.w"

static void
transform_then_translate(rigid_motion_t*t,double x,double y)
{
t->v[0]+= x;
t->v[1]+= y;
}

/*:70*//*71:*/
#line 1271 "./tspgen.w"

static void
transform_then_rotate(rigid_motion_t*t,double theta)
{
const double ct= cos(theta),st= sin(theta);
const double x= t->v[0],y= t->v[1];
t->v[0]= x*ct-y*st;
t->v[1]= x*st+y*ct;
t->theta+= theta;
}

/*:71*//*73:*/
#line 1290 "./tspgen.w"

static void
transform_then_transform(rigid_motion_t*t,rigid_motion_t*t2)

{
rigid_motion_t copy;
if(t==t2){copy= *t2;t2= &copy;}
transform_then_rotate(t,t2->theta);
transform_then_translate(t,t2->v[0],t2->v[1]);
}

/*:73*//*74:*/
#line 1304 "./tspgen.w"

static void
component_then_rotate(component_t*c,double theta)
{
transform_then_rotate(&(c->transform),theta);
hull_then_rotate(c->hull,theta);
}

static void
component_then_translate(component_t*c,double x,double y)
{
transform_then_translate(&(c->transform),x,y);
hull_then_translate(c->hull,x,y);
}


/*:74*//*79:*/
#line 1426 "./tspgen.w"

void rotate_hull_with_bias(component_t*comp,int is_foward,double packing_factor);
void
rotate_hull_with_bias(component_t*comp,int is_forward,double packing_factor)
{
hull_vertex_t*comph= comp->hull;
if(errno){perror("main: about to rotate comp");errorif(1,"Oops!");}
if(comp->hull_size>1){
double rot;
double alpha= -atan2(comph->y-comph->prev->y,comph->x-comph->prev->x);
double beta= -atan2(comph->next->y-comph->y,comph->next->x-comph->x);
errorif(errno==EDOM,"Sorry, my call to atan2 had a bad argument.");
errorif(errno==ERANGE,"Sorry, the answer to my atan2 call was out of range.");
while(alpha<beta)beta-= 2*pi;

/*82:*/
#line 1481 "./tspgen.w"

/*81:*/
#line 1452 "./tspgen.w"

errno= 0;


/*:81*/
#line 1482 "./tspgen.w"

{double y= 1-uniform_sample(angle_prng,0,1),t= exp(packing_factor*log(y));
switch(errno){
case 0:break;
case ERANGE:t= 1.0;/*81:*/
#line 1452 "./tspgen.w"

errno= 0;


/*:81*/
#line 1486 "./tspgen.w"
;break;
default:perror("main: about to sample alpha beta");errorif(1,"Oops!");
}
errorif(t<0.0||1.0<t,"alpha-beta mediation is %f, not in [0,1]",t);
if(!is_forward)t= 1-t;
rot= t*alpha+(1-t)*beta;
/*110:*/
#line 1931 "./tspgen.w"

#if TSPGEN_DEBUG >= 200
{double frac= (rot-beta)/(alpha-beta);
fprintf(stderr," angle %10.5f%% alpha and %10.5f%% beta %s\n",
100*frac,100*(1-frac),(is_forward?"forward":"reverse"));
errorif(frac<0.0||1.0<frac,"Frac %f out of [0,1]",frac);
}
#endif

/*:110*/
#line 1492 "./tspgen.w"

}




/*:82*/
#line 1441 "./tspgen.w"


component_then_rotate(comp,rot);
}
}

/*:79*//*93:*/
#line 1623 "./tspgen.w"

#define abs(A) ((A)<0 ? -(A) : (A))
static int
cmp_point(const void*ap,const void*bp)
{
const point_t a= *(const point_t*)ap,b= *(const point_t*)bp;
if(a.angle<b.angle)return-1;
if(a.angle>b.angle)return 1;
if(a.x<b.x)return-1;
if(a.x>b.x)return 1;
if(abs(a.y-west.y)<abs(b.y-west.y))return-1;
if(abs(a.y-west.y)>abs(b.y-west.y))return 1;
return 0;
}

/*:93*//*99:*/
#line 1731 "./tspgen.w"

static int
left_turn(point_t u,point_t v,point_t w)
{
const double
ad= (u.x-w.x)*(v.y-w.y),
bc= (u.y-w.y)*(v.x-w.x);
return ad>bc;
}


/*:99*//*104:*/
#line 1813 "./tspgen.w"

static void
write_component(FILE*file,component_t*c,int*city_num)
{
if(c->left){
transform_then_transform(&(c->left->transform),&(c->transform));
write_component(file,c->left,city_num);
}
if(c->left==NULL&&c->right==NULL){
printf("%10d %.16g %.16g\n",
*city_num,
c->transform.v[0],c->transform.v[1]);
*city_num+= 1;
}
if(c->right){
transform_then_transform(&(c->right->transform),&(c->transform));
write_component(file,c->right,city_num);
}
}

/*:104*//*109:*/
#line 1919 "./tspgen.w"

#if TSPGEN_DEBUG >= 100
static double
euc2hull(hull_vertex_t*ah,hull_vertex_t*bh)
{
const double dx= ah->x-bh->x,dy= ah->y-bh->y;
return sqrt(dx*dx+dy*dy);
}
#endif


/*:109*/
#line 252 "./tspgen.w"



int main(int argc,char**argv)
{
/*37:*/
#line 842 "./tspgen.w"

gb_prng_t*component_prng,*hull_prng;

/*:37*//*95:*/
#line 1682 "./tspgen.w"

int*stack;

/*:95*/
#line 257 "./tspgen.w"

/*6:*/
#line 324 "./tspgen.w"

/*9:*/
#line 394 "./tspgen.w"

{int i;
for(i= 1;i<argc;i++){
if(0==strcmp(argv[i],"--help")){
usage(stdout,argv);
exit(0);
}else if(0==strcmp(argv[i],"--version")){
version(stdout);
exit(0);
}
}
}


/*:9*/
#line 325 "./tspgen.w"

seed= -42L;
packing_factor= 100;
join_length_bias= 1;
try_forcing_isomorphic= 0;
force_failures= 0;
option_name= NULL;
while(1){
extern char*optarg;
const int opt= getopt(argc,argv,"s:p:j:hin:");
if(opt==EOF)break;
switch(opt){
case's':seed= atol(optarg);break;
case'n':option_name= optarg;break;
case'p':
packing_factor= atof(optarg);
errorif(packing_factor<0,
"Packing factor should be non-negative, but given %f",packing_factor);
break;
case'j':
join_length_bias= atof(optarg);
errorif(join_length_bias>1,
"Join length bias should be no more than 1, but given %f",
join_length_bias);
break;
case'i':try_forcing_isomorphic= 1;break;
case'h':usage(stdout,argv);exit(0);break;
case':':errorif(1,"some option is missing an argument");break;
case'?':usage(stderr,argv);errorif(1,"Unrecognized option");break;
default:errorif(1,"getopt returned character 0%o",opt);
}
}

/*:6*/
#line 258 "./tspgen.w"

/*36:*/
#line 831 "./tspgen.w"

component_prng= gb_prng_new(seed);
hull_prng= gb_prng_new(seed^42);
angle_prng= gb_prng_new(seed^(42*42));
errorif(
component_prng==NULL
||hull_prng==NULL
||angle_prng==NULL
,"Couldn't allocate three random number generators");

/*:36*/
#line 259 "./tspgen.w"

/*14:*/
#line 520 "./tspgen.w"

/*81:*/
#line 1452 "./tspgen.w"

errno= 0;


/*:81*/
#line 521 "./tspgen.w"

{int problem_line_read= 0;
char line[MAX_LINE_LENGTH];
while(fgets(line,MAX_LINE_LENGTH,stdin)){
switch(line[0]){
case'c':/*21:*/
#line 606 "./tspgen.w"

{char*p,parsed_suffix[MAX_LINE_LENGTH];
if((p= strstr(line,TSPLIB_COMMENT_PREFIX))!=NULL){
if(in_comment){mem_deduct(strlen(in_comment));free_mem(in_comment);}
in_comment= dup_string(p+strlen(TSPLIB_COMMENT_PREFIX));
}else if(1==sscanf(line,TSPLIB_NAME_PREFIX"%s",parsed_suffix)){

if(in_name){mem_deduct(strlen(in_name));free_mem(in_name);}
in_name= dup_string(parsed_suffix);
}else if(1==sscanf(line,LK_LENGTH_T_PREFIX"%s",parsed_suffix)){
if(strcmp(LENGTH_TYPE_STRING,parsed_suffix)){
fprintf(stderr,
"tspgen:Warning: tspgen uses %s for lengths, input used %s\n",
LENGTH_TYPE_STRING,p+strlen(LK_LENGTH_T_PREFIX));
}
}
}

/*:21*/
#line 526 "./tspgen.w"
break;
case'p':/*24:*/
#line 642 "./tspgen.w"

errorif(problem_line_read,
"A problem line has already been read, but I now see a second one: %s",
line);
{int num_read;
num_read= sscanf(line,"p edge %d %d ",&n,&ne);
errorif(num_read!=2,"Couldn't understand the problem line."
"  It should be of form: p edge <n> <n-1>");
errorif(n<2,"n, the number vertices, is %d but must be at least 2.",n);
errorif(ne!=n-1,"The number of edges is %d but should be %d.",ne,n-1);
/*16:*/
#line 558 "./tspgen.w"

in_tree= new_of(tree_t);
in_tree->ne= ne;
in_tree->num_filled= 0;
in_tree->edge= new_arr_of(edge_t,ne);

/*:16*/
#line 652 "./tspgen.w"

}

/*:24*/
#line 527 "./tspgen.w"
break;
case'e':/*26:*/
#line 661 "./tspgen.w"

errorif(in_tree==NULL,"We got an edge (e) line before a problem (p) line");
errorif(in_tree->num_filled>=in_tree->ne,
"Too many edges specified on input.  Got this one as first extra line: %s",
line);
{
int num_read;
int u= 0,v= 0;
length_t len;
num_read= sscanf(line,"e %d %d "length_t_native_inspec,
&u,&v,length_t_native_incast(&len));
errorif(num_read!=(2+length_t_native_incount),
"Couldn't understand an edge line. Parsed %d items; line is: %s",
num_read,line);
errorif(u<1||u>n,"First vertex %d is out of bounds in: %s",u,line);
errorif(v<1||v>n,"Second vertex %d is out of bounds in: %s",v,line);
u--;
v--;
in_tree->edge[in_tree->num_filled].city[0]= u;
in_tree->edge[in_tree->num_filled].city[1]= v;
in_tree->edge[in_tree->num_filled].cost= len;
in_tree->num_filled++;
}

/*:26*/
#line 528 "./tspgen.w"
break;
default:errorif(1,"Don't know how to parse input line:\n%s",line);
}
}
errorif(in_tree==NULL,"No tree found on the input");
errorif(in_tree->num_filled!=in_tree->ne,
"Read only %d edges when %d were expected",
in_tree->num_filled,in_tree->ne);
}
/*81:*/
#line 1452 "./tspgen.w"

errno= 0;


/*:81*/
#line 537 "./tspgen.w"


/*:14*/
#line 260 "./tspgen.w"

/*27:*/
#line 693 "./tspgen.w"

dsort(in_tree->edge,(size_t)ne,sizeof(edge_t),cmp_edge);

/*:27*/
#line 261 "./tspgen.w"

/*46:*/
#line 965 "./tspgen.w"

hull_vertex= new_arr_of(hull_vertex_t,n);
/*47:*/
#line 973 "./tspgen.w"