lk.w,v

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

	extra_backtrack = 0;
	random_seed = RIGHT_NOW;
	
@@
@@<Global variables@@>=
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;

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

@@
@@<Exported definitions@@>=
#define REP_ARRAY		1	/* Tour representations */
#define REP_TWO_LEVEL	2
#define REP_SPLAY_0		3
#define REP_SPLAY_1		4
#define REP_SPLAY_2		5
#define REP_SPLAY_3		6
#define REP_TWO_LEVEL_DEBUG		7
@@#

#define CAND_NN			1	/* Candidate predicate types --- they may be ||||'d */
#define CAND_NQ			2
#define CAND_DEL		4
@@#


@@  The \type{--help} option is the easiest one.
@@<Parse this option@@>=
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);
}

@@ If the user tells us the options are finished, then we shouldn't process more
options...  (Sort of obvious, huh?).

@@<Parse this option@@>=
if ( strcmp(argv[r],"--")==0 ) {
	more_options = 0;
	continue;
}

@@ Here we find out whether we should turn off rounding of Euclidean distance
computations.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--no-round")==0 ) {
	noround = 1;
	continue;
}

@@ Here we find out whether we should perform extra backtracking.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--extra-backtrack")==0 ) {
	extra_backtrack = 1;
	continue;
}

@@ Here we find out whether we should not perform extra backtracking.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--no-extra-backtrack")==0 ) {
	extra_backtrack = 0;
	continue;
}

@@ Here we find out whether we must print the tour in addition to
printing its length.
@@<Parse this option@@>=
if ( strcmp(argv[r],"-p")==0 || strcmp(argv[r],"--print")==0 ) {
	should_show_tour = 1;
	continue;
}

@@ Here we find out whether we should only print the minimum spanning
tree (and not compute a tour).
@@<Parse this option@@>=
if ( strcmp(argv[r],"-M")==0 || strcmp(argv[r],"--mst-only")==0 ) {
	mst_only = 1;
	continue;
}

@@ Here we find out whether we should only compute a Held-Karp lower
bound for the TSP.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--held-karp")==0 ) {
	held_karp_only = 1;
	continue;
}

@@ Here we find out whether we should only compute a Held-Karp lower
bound for the TSP and print the Lagrange multipliers used.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--held-karp-lambda")==0 ) {
	held_karp_lambda_only = 1;
	continue;
}

@@ Here we find out whether we should try to find a cheap weighted perfect
matching rather than a short tour.  
If we are to find a matching, then the tour representation is irrelevant,
because there is no tour.

@@<Parse this option@@>=
if ( strcmp(argv[r],"-m")==0 || strcmp(argv[r],"--matching")==0 ) {
	do_weighted_perfect_matching = 1;
	continue;
}



@@ Here we find out whether we should print version information.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--version")==0 ) {
	should_show_version = 1;
	@@<Print a version banner@@>@@;
	exit(0);
	continue;
}


@@ Here we find out whether we should reorder the cities according to
my spacefilling curve.  I suspect this improves the performance in
a non-uniform memory architecture (NUMA) machine.  I have my reasons.
@@^non-uniform memory architecture@@>
@@^NUMA@@>

@@<Parse this option@@>=
if ( strcmp(argv[r],"--sfc")==0 ) {
	should_sfc_reorder = 1;
	continue;
}

@@
@@<Set the option defaults@@>=
should_sfc_reorder = 0;

@@
@@<Global variables@@>=
int should_sfc_reorder;


@@ The quiet option sets the verbose variable.
@@<Parse this option@@>=
if ( strcmp(argv[r],"-q")==0 || strcmp(argv[r],"--quiet")==0 ) {
	verbose = 0;
	continue;
}

@@ The verbose option has an optional numerical argument.  We must check for
its presence.
@@<Parse this option@@>=
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;
}


@@ But we need to create the predicate |is_number|.  We allow an optional
single
unary minus.  We are paranoid about |NULL| parameters.

@@<Module subroutines@@>=
static int
is_number(char *p) 
{
	if ( p == NULL ) return 0;	/* Be paranoid. */
	if ( *p == '-' ) p++;		/* Skip over a single unary minus. */
	if ( *p == 0 ) return 0;	/* Need at least one digit. */
	for (;*p;p++) {
		if (!isdigit(*p)) return 0;
	}
	return 1;
}

@@	We've just used the |isdigit| macro which is defined in \file{ctype.h}.
@@<System headers@@>=
#include <ctype.h>

@@ While we're doing optional numerical arguments, 
let's do the \type{--iterate} switch as well.
@@<Parse this option@@>=
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;
}


@@ The maximum probe depth option has an optional numerical argument, which
goes into |max_generic_flips|.  
We must check this option's presence, as in the verbose option.

@@<Parse this option@@>=
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;
	if ( r+1<argc && is_number(argv[r+1]) ) max_generic_flips = atoi(argv[++r]);
	continue;
}

@@
@@<Global variables@@>=
int max_generic_flips;

@@
@@<Exported variables@@>=
extern int max_generic_flips;

@@ The default maximum probe depth should really be $\infty$ but |INT_MAX|,
the maximum representable signed integer should suffice.
@@<Set the option defaults@@>=
max_generic_flips=INT_MAX;

@@  We need the definition of |INT_MAX|
@@<System headers@@>=
#include <limits.h>

@@ The random seed option has an optional argument which might be
an integer and which might be the string \type{now}.
@@<Parse this option@@>=
if ( strcmp(argv[r],"--seed")==0 ) {
	random_seed = RIGHT_NOW;
	if ( r+1<argc ) {
		if ( is_number(argv[r+1]) ) random_seed = atol(argv[++r]);
		else if ( strcmp("now",argv[r+1])==0 ) r++; /* We already got |RIGHT_NOW|. */
	}
	continue;
}


@@ Lets tackle the starting tour option.  This has a mandatory string
argument chosen from a fixed list.  Argument string \type{random} has
an optional numeric parameter.
@@<Parse this option@@>=
if ( strcmp(argv[r],"-s")==0 || strcmp(argv[r],"--start")==0 ) {
	if ( r+1>=argc ) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr, "Need one of {canonical,greedy,greedydet,random [seed]}\n");
		exit(1);
	}
	r++;

	if ( strcmp(argv[r],"greedy")==0 ) construction_algorithm = CONSTRUCT_GREEDY_RANDOM;
	else if ( strcmp(argv[r],"greedydet")==0 ) construction_algorithm = CONSTRUCT_GREEDY;
	else if ( strcmp(argv[r],"canonical")==0 ) construction_algorithm = CONSTRUCT_CANONICAL;
	else if ( strcmp(argv[r],"random")==0 ) {
		construction_algorithm = CONSTRUCT_RANDOM;
		if ( r+1<argc && is_number(argv[r+1]) ) start_heuristic_param = atol(argv[++r]);
	} else {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr, "Need one of {canonical,greedy,greedydet,random [seed],best}\n");
		exit(1);
	}
	continue;
}

@@ Next comes the sort procedure option.  It also has a mandatory string
argument chosen from a fixed list.
@@<Parse this option@@>=
if ( strcmp(argv[r],"-S")==0 || strcmp(argv[r],"--sort")==0 ) {
	if ( r+1>=argc ) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr, "Need one of {dsort,qsort}\n");
		exit(1);
	}
	r++;

	if ( strcmp(argv[r],"qsort")==0 ) 
		sort = (void (*)(void *, size_t, size_t,int(*)(const void*,const void*)))qsort;
	else if ( strcmp(argv[r],"dsort")==0 ) sort = dsort;
	else {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr, "Need one of {dsort,qsort}\n");
		exit(1);
	}
	continue;
}

@@ We need to import the interface to |dsort|.
@@<Module headers@@>=
#include "dsort.h"

@@ The \type{--representation} switch also has a mandatory string
argument chosen from a fixed list.

@@<Parse this option@@>=
if ( strcmp(argv[r],"-r")==0 || strcmp(argv[r],"--representation")==0 ) {
	if ( r+1>=argc ) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr, "Need one of {array,splay [level],two-level,tld}\n");
		exit(1);
	}
	r++;
	if ( strcmp(argv[r],"array")==0 ) {
		representation = REP_ARRAY;
	} else if ( strcmp(argv[r],"two-level")==0 ) {
		representation = REP_TWO_LEVEL;
	} else if ( strcmp(argv[r],"tld")==0 ) {
		representation = REP_TWO_LEVEL_DEBUG;
	} else if ( strcmp(argv[r],"splay")==0 ) {
		int level=0;
		if ( r+1<argc && is_number(argv[r+1]) ) level = atoi(argv[++r]);
		if ( level<0 || level>3 ) {
			@@<Print the command line; split at |r|@@>@@;
			fprintf(stderr,"Splay level must be 0, 1, 2, or 3\n");
			exit(1);
		}
		representation = REP_SPLAY_0 + level;
	} else {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr, "Need one of {array,splay [level],two-level,tld}\n");
		exit(1);
	}
	continue;
}

@@ The PostScript option  needs a filename.
@@<Parse this option@@>=
if ( strcmp(argv[r],"-P")==0 || strcmp(argv[r],"--postscript")==0 ) {
	if (r+1>=argc) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr,"Need a file name\n");
		exit(1);
	}
	r++;
	if ( postscript_filecount ) {
		fprintf(stderr,
			"Warning: %s already specified as PostScript output file; "
			"%s overrides\n", 
			PostScript_filename, argv[r]);
	}
	postscript_filecount++;
	PostScript_filename = argv[r];
	continue;
}

@@ The lower bound option needs a name and a floating-point number.
In most UNIX shells, names with spaces can be specified if quoted.

@@<Parse this option@@>=
if ( strcmp(argv[r],"-l")==0 || strcmp(argv[r],"--lower-bound")==0 ) {
	if (r+1>=argc) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr,"Need a lower bound name\n");
		exit(1);
	}
	if (r+2>=argc) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr,"Need a lower bound on length\n");
		exit(1);
	}
	r++;
	lower_bound_name=dup_string(argv[r]);
	r++;
	lower_bound_value=atof(argv[r]);
	continue;
}

@@ The upper bound option is very much like the lower bound option.

@@<Parse this option@@>=
if ( strcmp(argv[r],"-u")==0 || strcmp(argv[r],"--upper-bound")==0 ) {
	if (r+1>=argc) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr,"Need an upper bound name\n");
		exit(1);
	}
	if (r+2>=argc) {
		@@<Print the command line; split at |r|@@>@@;
		fprintf(stderr,"Need an upper bound on length\n");
		exit(1);
	}
	r++;
	upper_bound_name=dup_string(argv[r]);
	r++;
	upper_bound_value=atof(argv[r]);
	continue;
}

@@ The Held-Karp upper bounding option is 
valid only for the TSP.

@@<Check option prerequisites@@>=
errorif( (held_karp_only || held_karp_lambda_only) && do_weighted_perfect_matching,
	"Held-Karp lower bound valid only for the TSP, not perfect matching\n");


@@ Now we get to the most involved option: the candidate list predicate 
expression.  Fortunately, it's at worst just the disjunction of a number of
basic predicates.

@@<Parse this option@@>=
if ( strcmp(argv[r],"-c")==0 || strcmp(argv[r],"--candidate")==0 ) {
	int numeric_param;
	candidate_expr = cand_nn_k = cand_nq_k = cand_del_d = 0;
	r++;
	do {
		if ( r>=argc ||
				(strcmp(argv[r],"nn")!=0 && strcmp(argv[r],"nq")!=0 &&
				strcmp(argv[r],"del")!=0) 
		   ) {
			@@<Print the command line; split at |r|@@>@@;
			fprintf(stderr,"Need one of {nn,nq,del}\n");
			exit(1);
		}
		if ( r+1>=argc || !is_number(argv[r+1]) ) {
			r++;
			@@<Print the command line; split at |r|@@>@@;
			fprintf(stderr,"Need a numeric parameter\n");
			exit(1);
		} else numeric_param = atoi(argv[r+1]);
		if ( strcmp(argv[r],"nn")==0 ) {
			candidate_expr |= CAND_NN;
			if ( cand_nn_k < numeric_param ) cand_nn_k = numeric_param;
		} else if ( strcmp(argv[r],"nq")==0 ) {
			candidate_expr |= CAND_NQ;
			if ( cand_nq_k <