lk.c

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

#define DEFAULT_VERBOSE 5
#define XNORM(V) (tsp_instance->xmax==tsp_instance->xmin?0.0: \
((V) -tsp_instance->xmin) /(tsp_instance->xmax-tsp_instance->xmin) /(1+10*DBL_EPSILON) ) ;
#define YNORM(V) (tsp_instance->ymax==tsp_instance->ymin?0.0: \
((V) -tsp_instance->ymin) /(tsp_instance->ymax-tsp_instance->ymin) /(1+10*DBL_EPSILON) ) ; \

#define sgn(value) ((value) <0?-1:((value) >0?1:0) )  \
 \

#define ALWAYS_BUILD_DECLUSTER_STRUCTURES  \
(LK_BUILD_DECLUSTER_STRUCTURES \
||JBMR_DECLUSTER_IN_ELIGIBILITY_TEST \
||JBMR_DECLUSTER_IN_GREEDY)  \

#define my_abs(A) ((A) <0?-(A) :(A) ) 
#define within_epsilon(A,B) (my_abs(((A) -(B) ) ) <0.5)  \

/*2:*/
#line 623 "./lk.w"

#include <config.h>
#include "lkconfig.h"
/*3:*/
#line 642 "./lk.w"

#include <stdio.h>
#ifndef __USE_MISC
#define __USE_MISC  
#endif
#include <stdlib.h>
#include <stddef.h>
#include <string.h>
#include <math.h>
#if defined(OS_HAS_BROKEN_HEADERS)
#define FIXINCLUDES_USE_RCS_ID
#define FIXINCLUDES_NEED_GETHOSTNAME
#define FIXINCLUDES_NEED_TIME_STUFF
#include "fixincludes.h"
#undef FIXINCLUDES_USE_RCS_ID
#undef FIXINCLUDES_NEED_GETHOSTNAME
#undef FIXINCLUDES_NEED_TIME_STUFF
#endif

/*:3*//*35:*/
#line 1278 "./lk.w"

#include <ctype.h>

/*:35*//*41:*/
#line 1320 "./lk.w"

#include <limits.h>

/*:41*//*60:*/
#line 1741 "./lk.w"

#if HAVE_TIME_H
#include <time.h>
#endif

/*:60*//*73:*/
#line 1982 "./lk.w"

#include <float.h>


/*:73*//*124:*/
#line 2584 "./lk.w"

#if HAVE_SYS_PARAM_H
#include <sys/param.h>
#endif
#if HAVE_UNISTD_H
#include <unistd.h> 
#endif

/*:124*/
#line 626 "./lk.w"

/*63:*/
#line 1793 "./lk.w"

#define LENGTH_USE_RCS_ID
#include "length.h"
#undef LENGTH_USE_RCS_ID
#include "read.h"

/*:63*//*104:*/
#line 2421 "./lk.w"

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

/*:104*//*113:*/
#line 2483 "./lk.w"

#include "prng.h"

/*:113*/
#line 627 "./lk.w"

/*5:*/
#line 674 "./lk.w"

#include "lk.h"


/*:5*//*45:*/
#line 1385 "./lk.w"

#include "dsort.h"

/*:45*//*57:*/
#line 1655 "./lk.w"

#include "pool.h"
#include "dict.h"
#include "kdtree.h"


/*:57*//*67:*/
#line 1841 "./lk.w"

#include "ascend.h"

/*:67*//*75:*/
#line 2007 "./lk.w"

#include "resource.h"

/*:75*//*89:*/
#line 2171 "./lk.w"

#include "declevel.h"
#include "decluster.h"

/*:89*//*91:*/
#line 2193 "./lk.w"

#include "nn.h"

/*:91*//*99:*/
#line 2321 "./lk.w"

#include "array.h"
#include "twolevel.h"

/*:99*//*105:*/
#line 2426 "./lk.w"

#include "construct.h"

/*:105*//*112:*/
#line 2479 "./lk.w"

#include "jbmr.h"

/*:112*//*116:*/
#line 2505 "./lk.w"

#include "match.h"

/*:116*/
#line 628 "./lk.w"


/*8:*/
#line 822 "./lk.w"

#if HAVE_TIME
#define RIGHT_NOW ((long)time(NULL))
#else
#define RIGHT_NOW (902598118L)
#endif

/*:8*/
#line 630 "./lk.w"

/*13:*/
#line 999 "./lk.w"

static char*filename;
static int mst_only;
static int do_weighted_perfect_matching;
static int held_karp_only;
static int held_karp_lambda_only;

/*:13*//*62:*/
#line 1787 "./lk.w"

static int n;
static FILE*TSPLIB_in,*ps_out;

/*:62*//*76:*/
#line 2014 "./lk.w"

static int last_resource_mark;

/*:76*//*106:*/
#line 2430 "./lk.w"

static int*tour;

/*:106*/
#line 631 "./lk.w"

/*15:*/
#line 1035 "./lk.w"

int verbose,iterations,should_show_tour,should_show_version;
int representation,construction_algorithm;
long start_heuristic_param;
int candidate_expr,cand_nn_k,cand_nq_k,cand_del_d;
char*PostScript_filename,*lower_bound_name,*upper_bound_name;
void(*sort)(void*a,size_t n,size_t es,int(*cmp)(const void*,const void*));
int noround;
double upper_bound_value,lower_bound_value;
int extra_backtrack;
long random_seed;

/*:15*//*31:*/
#line 1239 "./lk.w"

int should_sfc_reorder;


/*:31*//*38:*/
#line 1307 "./lk.w"

int max_generic_flips;

/*:38*//*64:*/
#line 1801 "./lk.w"

tsp_instance_t*tsp_instance;

/*:64*//*69:*/
#line 1898 "./lk.w"

int*original_city_num= NULL;

/*:69*//*79:*/
#line 2042 "./lk.w"

int begin_data_structures_mark;

/*:79*//*88:*/
#line 2161 "./lk.w"

decluster_tree_t*mst;
length_t mst_len;

/*:88*//*96:*/
#line 2260 "./lk.w"

int(*tour_next)(int)= NULL;
int(*tour_prev)(int)= NULL;
int(*tour_between)(int,int,int)= NULL;
void(*tour_flip)(int,int,int,int)= NULL;
void(*tour_set)(int const*)= NULL;
void(*tour_setup)(int n)= NULL;
void(*tour_cleanup)(void)= NULL;

/*:96*//*108:*/
#line 2442 "./lk.w"

length_t incumbent_len;

/*:108*/
#line 632 "./lk.w"

/*128:*/
#line 2626 "./lk.w"

static void lk_cleanup(void);

/*:128*/
#line 633 "./lk.w"

/*34:*/
#line 1264 "./lk.w"

static int
is_number(char*p)
{
if(p==NULL)return 0;
if(*p=='-')p++;
if(*p==0)return 0;
for(;*p;p++){
if(!isdigit(*p))return 0;
}
return 1;
}

/*:34*//*72:*/
#line 1941 "./lk.w"

int cmp_sfc_Moore(const void*voida,const void*voidb);
int
cmp_sfc_Moore(const void*voida,const void*voidb)
{
int ai= *(const int*)voida,bi= *(const int*)voidb;
double a[2],b[2];
int ax,ay,bx,by,aq,bq,lastq= 12;
int order[13][4]= 
{{0,3,1,2},{0,1,3,2},{0,3,1,2},{2,3,1,0},
{0,3,1,2},{2,3,1,0},{0,1,3,2},{2,1,3,0},
{0,1,3,2},{2,3,1,0},{2,1,3,0},{2,1,3,0},{3,2,0,1}};
int rewrite[13][4]= 
{{1,3,2,0},{4,1,7,6},{8,9,2,0},{9,0,5,10},
{1,9,2,0},{9,0,5,11},{2,1,11,6},{7,10,8,5},
{2,1,7,6},{9,4,5,10},{7,10,6,3},{7,10,6,5},{7,10,2,0}};

a[0]= XNORM(tsp_instance->coord[ai].x[0]);
b[0]= XNORM(tsp_instance->coord[bi].x[0]);
a[1]= YNORM(tsp_instance->coord[ai].x[1]);
b[1]= YNORM(tsp_instance->coord[bi].x[1]);



errorif(a[0]<0||a[0]>=1,"cmp_sfc_Moore: a[0] out of range [0,1): %f",a[0]);
errorif(a[1]<0||a[1]>=1,"cmp_sfc_Moore: a[1] out of range [0,1): %f",a[1]);
errorif(b[0]<0||b[0]>=1,"cmp_sfc_Moore: b[0] out of range [0,1): %f",b[0]);
errorif(b[1]<0||b[1]>=1,"cmp_sfc_Moore: b[1] out of range [0,1): %f",b[1]);
while(1){
if(a[0]==b[0]&&a[1]==b[1])return 0;
ax= a[0]>=0.5;ay= a[1]>=0.5;aq= (ax<<1)+ay;
bx= b[0]>=0.5;by= b[1]>=0.5;bq= (bx<<1)+by;
if(aq!=bq)return order[lastq][aq]-order[lastq][bq];
lastq= rewrite[lastq][aq];
a[0]= (2*a[0])-ax;a[1]= (2*a[1])-ay;
b[0]= (2*b[0])-bx;b[1]= (2*b[1])-by;
}
}

/*:72*//*127:*/
#line 2619 "./lk.w"

static void lk_cleanup(void)
{
/*86:*/
#line 2143 "./lk.w"

#if ALWAYS_BUILD_DECLUSTER_STRUCTURES
if(verbose>=50)printf("Cleaning up MST\n");
if(mst!=NULL)decluster_cleanup_tree(mst);
if(verbose>=50)printf("Cleaning up decluster data structures\n");
decluster_cleanup();
#endif

/*:86*//*92:*/
#line 2200 "./lk.w"

if(verbose>=50)printf("Cleaning up nn structure\n");
nn_cleanup();


/*:92*//*101:*/
#line 2369 "./lk.w"

if(!do_weighted_perfect_matching){
if(verbose>=50)printf("Cleaning up tour data structure\n");
if(tour_cleanup!=NULL)(*tour_cleanup)();
}

/*:101*//*107:*/
#line 2434 "./lk.w"

if(!do_weighted_perfect_matching){
if(verbose>=50)printf("Cleaning initial tour\n");
free_mem(tour);
}

/*:107*//*115:*/
#line 2491 "./lk.w"

if(!do_weighted_perfect_matching){
if(verbose>=50)printf("Cleaning up jbmr structures\n");
jbmr_cleanup();
}

/*:115*//*119:*/
#line 2525 "./lk.w"

if(do_weighted_perfect_matching){
match_cleanup();
}

/*:119*//*131:*/
#line 2652 "./lk.w"

resource_cleanup();


/*:131*/
#line 2622 "./lk.w"

}

/*:127*/
#line 634 "./lk.w"

/*6:*/
#line 686 "./lk.w"

int
lk_main(int argc,char**argv)
{
/*137:*/
#line 2758 "./lk.w"

length_t validate_len;
double double_validate_len,ordered_double_len,raw_len;


/*:137*/
#line 690 "./lk.w"

/*74:*/
#line 2000 "./lk.w"

mem_usage_reset();
resource_setup(50);

error_postcleanup_stats= resource_abnormal_exit_output;

/*:74*//*129:*/
#line 2634 "./lk.w"

error_cleanup= lk_cleanup;

/*:129*//*130:*/
#line 2642 "./lk.w"

error_postcleanup_stats= mem_report;

/*:130*/
#line 691 "./lk.w"

/*12:*/
#line 974 "./lk.w"

/*14:*/
#line 1008 "./lk.w"

verbose= DEFAULT_VERBOSE;
iterations= 1;
should_show_tour= 0;
should_show_version= 0;
mst_only= 0;
held_karp_only= 0;
held_karp_lambda_only= 0;
do_weighted_perfect_matching= 0;
representation= REP_ARRAY;
candidate_expr= CAND_NN;
cand_nn_k= 20;
cand_nq_k= 5;
cand_del_d= 3;
construction_algorithm= CONSTRUCT_GREEDY_RANDOM;
start_heuristic_param= 42L;
PostScript_filename= NULL;
lower_bound_name= NULL;
lower_bound_value= 1.0;
upper_bound_name= NULL;
upper_bound_value= 0.0;
sort= (void(*)(void*,size_t,size_t,int(*)(const void*,const void*)))qsort;
noround= 0;
extra_backtrack= 0;
random_seed= RIGHT_NOW;

/*:14*//*30:*/
#line 1235 "./lk.w"

should_sfc_reorder= 0;

/*:30*//*40:*/
#line 1316 "./lk.w"

max_generic_flips= INT_MAX;

/*:40*/
#line 975 "./lk.w"

{int r,filecount= 0,postscript_filecount= 0,more_options= 1;
filename= NULL;
for(r= 1;r<argc;r++){
if(more_options&&argv[r][0]=='-'&&argv[r][1]!=0){
/*18:*/
#line 1078 "./lk.w"

if(strcmp(argv[r],"-h")==0||strcmp(argv[r],"--help")==0){
printf("Usage: %s [options] [filename]\n",argv[0]);
printf(" -                          TSPLIB input on stdin\n");
printf(" --                         End options\n");
printf(" -l --lower-bound <name> <length> \n"
"                            Give lower bound, enable milestones\n");
printf(" -c --candidate <p>         Specify candidate cities\n");
printf("      <p> ::= <bp> | <bp> or <p>\n");
printf("      <bp> ::= nn <k> | nq <k>\n");
printf("        nn <k> is k nearest neighbours\n");
printf("        nq <k> is k nearest neighbours in each quadrant\n");
printf("        (nq 10 gives a typical city 40 neighbours)\n");
printf(" --held-karp                Compute approx Held-Karp TSP lower bound,\n");
printf("                            then exit successfully.\n");
printf(" --held-karp-lambda         Compute approx Held-Karp TSP lower bound,\n");
printf("                            print the best Lagrange multipliers found,\n");
printf("                            then exit successfully.\n");
printf(" --extra-backtrack          TSP: Try more possibilities for t7 and t8\n");
printf("                            when t1..t6 is not a legal 3-change\n");
printf("                            Matching: Backtrack to t6, not t8\n");
printf(" --no-extra-backtrack       Turn off --extra-backtrack\n");
printf("                            This is default; compatible with JBMR\n");
printf(" -h --help                  Give this help, then exit\n");
printf(" -i --iterate [n]           Iterated LK n times (default is 1, n default is 20)\n");
printf(" -m --matching              Find cheap weighted perfect matchings\n"
"                            rather than short tours\n");
printf(" -M --mst-only              Print the a minimum spanning tree\n");
printf("                            and then exit successfully\n");
printf(" --maxdepth [d]             Limit probe depth to d generic flips\n");
printf(" --no-round                 Don't round distance computations\n");
printf(" -p --print                 Print the LK-optimal tour or matching\n");
printf(" -P --postscript <file>     Generate PostScript output and write to <file>\n");
printf(" -q --quiet                 Same as --verbose 0\n");
printf(" -r --representation <rep>  Specify tour representation (default is array)\n");
printf("      <rep> ::= array | two-level | tld\n");
printf("                            (tld is two-level debugging mode.)\n");
printf("                            (Option -r has no effect if finding a matching)\n");
printf(" -s --start <kind>          Specify staring tour (matching) algorithm\n");
printf("      <kind> ::= canonical | greedy | greedydet | random [seed]\n"
"                            canonical is 1,2,3,...,n;\n"
"                            random is random order;\n"
"                            greedy is 2/3-1/3 randomized Greedy tour;\n"
"                            greedydet is pure (deterministic) Greedy tour;\n"
"                            Default is `greedy'.\n");
printf(" --seed [<s>]                 Random number seed\n");
printf("      <s> ::= <integer> | now \n");
printf("      Default is now, which means seconds since midnight Jan 1, 1970\n");
printf(" --sfc                      Reorder cities by Moore's space filling curve\n");
printf(" -S --sort <kind>           Specify sorting procedure (default is qsort)\n");
printf("      <kind> ::= qsort (from system library) | dsort (Bentley&McIlroy qsort)\n");
printf(" -u --upper-bound <name> <length> \n"
"                            Give upper bound, required for Held-Karp\n");
printf(" -v --verbose [n]           Set verbose level (default is 5, n default is 100)\n");
printf(" --version                  Print LK version number, then exit successfully\n");
exit(0);
}

/*:18*//*19:*/
#line 1139 "./lk.w"

if(strcmp(argv[r],"--")==0){
more_options= 0;
continue;
}

/*:19*//*20:*/
#line 1147 "./lk.w"

if(strcmp(argv[r],"--no-round")==0){
noround= 1;
continue;
}

/*:20*//*21:*/
#line 1154 "./lk.w"

if(strcmp(argv[r],"--extra-backtrack")==0){
extra_backtrack= 1;
continue;
}

/*:21*//*22:*/
#line 1161 "./lk.w"

if(strcmp(argv[r],"--no-extra-backtrack")==0){
extra_backtrack= 0;
continue;
}

/*:22*//*23:*/
#line 1169 "./lk.w"

if(strcmp(argv[r],"-p")==0||strcmp(argv[r],"--print")==0){
should_show_tour= 1;
continue;
}

/*:23*//*24:*/
#line 1177 "./lk.w"

if(strcmp(argv[r],"-M")==0||strcmp(argv[r],"--mst-only")==0){
mst_only= 1;
continue;
}

/*:24*//*25:*/
#line 1185 "./lk.w"

if(strcmp(argv[r],"--held-karp")==0){
held_karp_only= 1;
continue;
}

/*:25*//*26:*/
#line 1193 "./lk.w"

if(strcmp(argv[r],"--held-karp-lambda")==0){
held_karp_lambda_only= 1;
continue;
}

/*:26*//*27:*/
#line 1204 "./lk.w"

if(strcmp(argv[r],"-m")==0||strcmp(argv[r],"--matching")==0){
do_weighted_perfect_matching= 1;
continue;
}



/*:27*//*28:*/
#line 1213 "./lk.w"

if(strcmp(argv[r],"--version")==0){
should_show_version= 1;
/*55:*/
#line 1631 "./lk.w"

printf("LK %s",VERSION_STRING);
#if JBMR_DECLUSTER_IN_ELIGIBILITY_TEST
printf("de");
# if JBMR_DECLUSTER_IN_GREEDY
printf("g");
# endif
#endif
printf("\n");

/*:55*/
#line 1216 "./lk.w"

exit(0);
continue;
}


/*:28*//*29:*/
#line 1228 "./lk.w"

if(strcmp(argv[r],"--sfc")==0){
should_sfc_reorder= 1;
continue;
}

/*:29*//*32:*/
#line 1244 "./lk.w"

if(strcmp(argv[r],"-q")==0||strcmp(argv[r],"--quiet")==0){
verbose= 0;
continue;
}

/*:32*//*33:*/
#line 1252 "./lk.w"

if(strcmp(argv[r],"-v")==0||strcmp(argv[r],"--verbose")==0){
verbose= 100;
if(r+1<argc&&is_number(argv[r+1]))verbose= atoi(argv[++r]);
continue;
}


/*:33*//*36:*/
#line 1283 "./lk.w"

if(strcmp(argv[r],"-i")==0||strcmp(argv[r],"--iterate")==0){
iterations= 20;
if(r+1<argc&&is_number(argv[r+1]))iterations= atoi(argv[++r]);
continue;
}


/*:36*//*37:*/
#line 1295 "./lk.w"

if(strcmp(argv[r],"--maxdepth")==0){
#if !defined(JBMR_LIMIT_PROBE_DEPTH)
errorif(1,
"Option --maxdepth requires JBMR_LIMIT_PROBE_DEPTH compilation flag.\n");
#endif
max_generic_flips= 50;