lk.w,v

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

Fixed the semantics of fwrite.

Revision 1.129  96/08/16  13:04:58  neto
Added fixincludes.

Revision 1.128  96/08/16  12:42:06  neto
Converted putchar to printf.  Otherwise, I'd never get a prototype
for SunOS' \_flusbuf.

Revision 1.127  96/08/15  14:45:51  neto
Enable  definition of length\_rcs\_id

Revision 1.126  96/08/15  13:29:00  neto
Make it pass -Wall

Revision 1.125  96/08/15  12:36:18  neto
No longer use the "upper" module.

Revision 1.124  96/08/14  15:16:42  neto
Fixed printing of command line.

Revision 1.123  96/08/14  14:45:14  neto
Really fixed it this time.

Revision 1.122  96/08/14  14:41:38  neto
Fixed cast of qsort.

Revision 1.121  96/08/14  14:37:27  neto
Fix the declaration of sort to match ANSI (use void instead of char).

Revision 1.120  96/08/14  13:48:56  neto
Add an option to specify sorting procedure.

Revision 1.119  96/08/02  14:38:42  neto
Update documentation and printing of command line to reflect new
construct implementation.
.\

Revision 1.118  96/08/02  14:10:52  neto
Update command line switches to match what construct offers.
Updated the construct() call.

Revision 1.117  96/07/30  16:29:32  neto
Flush stdout after each call to resource\_mark

Revision 1.116  96/07/29  17:08:19  neto
Fixed to compile

Revision 1.115  96/07/29  16:20:02  neto
Added *\_rcs\_id.
Made sure RCS log is activated within this file.

}

@@*The Lin-Kernighan program.
This program implements the Lin-Kernighan heuristic for the 
traveling salesman problem (TSP).  

See the famous 1973 paper {\sl 
An effective heuristic algorithm for the traveling salesman
        problem},
Operations Research,
    {\bf 21},
    pp.~498--516,
    1973,
  by Shen Lin and Brian Kernighan
for the first description and motivation of this algorithm.
A more modern source is the chapter by Johnson and McGeoch, 
    {``The traveling salesman problem:
            a case study''},
Chapter 8, 
    pp.~215--310, in
    {\sl Local Search in Combinatorial Optimization},
    {Emile Aarts and Jan Karel Lenstra} editors,
    {John Wiley \& Sons}, 1997,
    {Wiley Interscience series in discrete mathematics
        and optimziation}.
That chapter
describes the Lin-Kernighan  algorithm (LK from now on) in the
context of other local search algorithms for the TSP.  If you plan to implement
LK, you really ought to study that chapter closely.  When it comes becomes
available,
you should also read the implementation report by Johnson, Bentley, McGeoch
and Ostheimer.

(Then again, why not just work from this code base?)

@@ The outline of the main module is as follows:
@@c
#include <config.h>
#include "lkconfig.h"
@@<System headers@@>@@;
@@<Early module headers@@>@@;
@@<Module headers@@>@@;

@@<Module definitions@@>@@;
@@<Module variables@@>@@;
@@<Global variables@@>@@;
@@<Static prototypes@@>@@;
@@<Module subroutines@@>@@;
@@<Subroutines@@>@@;
const char *lk_rcs_id="$Id: lk.w,v 1.274 1998/10/16 20:42:04 neto Exp neto $";


@@ We will be using many routines from external libraries.  The interfaces
to those routines are described in the following headers.

@@<System headers@@>=
#include <stdio.h>
#ifndef __USE_MISC
#define __USE_MISC		/* Linux needs this to get the definition of |nrand48| */
#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

@@ The exported interface is contained in the \file{lk.h} header file,
which has the following form.

@@(lk.h@@>=
extern const char *lk_rcs_id;
@@<Exported type definitions@@>@@;
@@<Exported definitions@@>@@;
@@<Exported variables@@>@@;
@@<Exported subroutines@@>@@;

@@ To ensure consistency between the interface and the implementation,
we import our own interface.

@@<Module headers@@>=
#include "lk.h"


@@*The main module.
The program parses the command line, reads an instance,
builds support data structures, constructs a starting tour,
optimized that tour with Lin-Kernighan, then prints the results.

If the system supports resource measurement, then resource usage
is reported over the phases of execution.

@@<Subroutines@@>=
int
lk_main(int argc, char **argv) 
{
	@@<|main| variables@@>@@#
	@@<Do basic initialization@@>@@;
	@@<Parse the command line@@>@@;
	@@<Print a banner@@>@@;
	@@<Read the TSP instance@@>@@;
	@@<Possibly reorder the cities@@>@@;
	@@<Allocate the space for this instance@@>@@;
	@@<Build the data structures@@>@@;
	if ( do_weighted_perfect_matching ) {
		@@<Construct a starting matching@@>@@;
	} else {
		@@<Construct a starting tour@@>@@;
		@@<Possibly produce only a Held-Karp lower bound@@>@@;
	}
	@@<Run the Lin-Kernighan algorithm@@>@@;
	@@<Stop the timers and print interval times@@>@@;
	@@<Validate and print the result@@>@@;
	@@<Free the allocated space@@>@@;
	return 0;
}

@@ Just to make it clear that this is a library and not a program, 
|lk_main| is exported as a library function.
@@<Exported subroutines@@>=
int lk_main(int argc, char **argv);


@@*Command-line options.
We'll be experimenting with this program, so we need to put in lots of
control knobs, better known as command-line options.  

This program follows the
GNU conventions.  Each option has both a short name consisting of a 
minus sign followed by a single letter, and a long name, 
consisting of two minus
signs  followed by a string of alphanumeric characters.  Parameters to
an option follow whitespace after the option name.  Options are terminated
by `{\tt --}', thus allowing filename parameters to the program to
begin with a leading dash while not being in conflict with command-line
option names or their parameters.

In this program, parameters are processed from left to right, so later
options may override earlier options.

Option \type{-h} or \type{--help} prints out help for command-line 
options, and then quits the program.

Option \type{--version} makes LK print out version information
and then exit successfully.


Option \type{-v} or \type{--verbose} turns on the ``verbose'' mode of 
output.  It has an optional numeric argument, where 0 means no verbose
reporting, and higher values lead to more verbiage.  If no number
is provided, then a default value of 100 is used.  If this option is not
used at all, then a default value of |DEFAULT_VERBOSE| is used.

Option \type{-P} or \type{--postscript} has a mandatory filename argument.
Debugging PostScript output is placed in that file.

Option \type{-q} or \type{--quiet} is synonymous with \type{--verbose 0}.

Option \type{-i} or \type{--iterate} tells us that Lin-Kernighan should
be iterated in the sense of Johnson'90.  An optional numeric
parameter specifies the number of iterations. 
If no numeric parameter is given, a value of 20 is used.  

Option \type{-r} or \type{--representation} has a mandatory string argument
specifying which oriented tour representation we should use.  
Allowable arguments are \type{array} (the default), 
\type{splay}, 
\type{two-level}, 
\type{tld}, 
and \type{segment}.   
Only \type{array} and \type{two-level} are supported for now.
Unimplemented suboption \type{splay} has an optional numeric argument that must be either
0, 1, 2, or 3, corresponding to each of the four levels of faithfullness
to the splaying discipline (lower numbers mean less splaying).  See
Fredman \etal.~for the details on this.
Option \type{tld} is ``two-level-debug''; it checks the two-level tree
implementation against the array implementation; the program must have
been compiled
with symbol |TWOLEVEL_DEBUG| defined.


Option \type{-c} or \type{--candidate} has a mandatory expression 
as a parameter.  For each city $i$ in the tour, the expresion specifies
a predicate that defines which other cities may appear on $i$'s candidate list.
The basic predicates are as follows: \type{nn $k$} is satisfied by the $k$
nearest neighbours of city $i$; \type{nq $k$} is satisfied by the $k$
nearest neighbours of $i$ in each quadrant surrounding $i$;
\type{del $d$} is satisfied by any city that is at most $d$ steps away 
from city $i$ in the Delauney triangulation of all the cities.
(The last two options are only applicable to two-dimensional Euclidean 
instances.)  Basic predicates may be combined by the keyword \type{or}.
Delauney triangulation is not yet supported.

Option \type{-s} or \type{--start} has a mandatory string argument specifying
what the starting tour should be.  Allowable arguments are 
\type{random} (random tour),
\type{canonical} (tour $1,2,3,\ldots,n$),
\type{greedy} (ranomized greedy, or multi-fragment heuristic---the default)
\type{greedydet} (deterministic greedy)
%\type{ai} (arbitrary insertion), 
%and \type{best} (best of the above---the default).
Suboption \type{random} takes an optional long integer seed.
If the points are to be sorted by a spacefilling curve, 
by specifying option \type{--sfc}, then the canonical tour is the cities
in space-filling curve order.

Option \type{-S} or \type{--sort} has a mandatory string argument specifying
what the sorting procedure should be.  Allowable arguments are
\type{qsort} (the system C library's implementation of Quicksort),
or
\type{dsort} (a determinstic sorting function---Bentley and McIlroy's 
implementation of Quicksort, see INSERT REFERENCE).
It so happens that some system library's implementations of |qsort|
are not deterministic; running times of this program have varied
by up to a factor of 2 because of it, and they sometimes return different
length tours.  Using \type{-S dsort} ensures consistent results.

Option \type{--seed} takes an optional argument, an integer, or 
the string \type{now}.   If no argument is given, then \type{now} is
assumed.  Argument \type{now} is interpreted as the number of seconds 
since midnight January 1, 1970.  That is, |now==time(NULL)| under any
POSIX-compliant system.  
(However, if |time()| is not supported on the system, then 
a constant of 902598118 is used.  That's the current time as I write.)
The seed is used to seed generators
for various randomized subalgorithms.
@@^POSIX@@>

@@<Module definitions@@>=
#if HAVE_TIME
#define RIGHT_NOW ((long)time(NULL))
#else
#define RIGHT_NOW (902598118L)
#endif

@@ Option \type{-p} or \type{--print} tells us to print the LK-optimal
tour (or matching)
when
we're done.  The default is to not print the tour (matching), because I'm
usually
interested in the far less bulky output of only run times and tour lengths.

Option \type{-M} or \type{--mst-only} tells us to 
compute the minimum spanning tree, output it, 
and then exit.  In particular, a tour is not computed
when this option is used.  This option is useful for cloning an input.
See the program \file{tspgen}.  For post-processing, this option is
probably most useful with the \type{--quiet} so that human-oriented
noise is removed.

Option \type{--no-round} tells us to avoid rounding computed distances.  
The TSPLIB file format specifies that instances of type |EUC_2D| use
rounded distances, \ie, 
$$\hbox{\it cost}(u,v) = \left\lfloor 0.5 + \sqrt{(u_x-v_x)^2+(u_y-v_y)^2} 
\right\rfloor.$$
However, David Johnson's runs on Euclidean inputs do not do any rounding
(I had a pleasant conversation with him at FOCS'96 in Burlington, Vermont).
Using this option allows us to override the TSPLIB definition and just use
the ordinary Euclidean metric:
$$\hbox{\it cost}(u,v) = \sqrt{(u_x-v_x)^2+(u_y-v_y)^2}.$$
This is useful only when the |length_t| type is either an exact type or 
a floating point type.

Option \type{--sfc} asks us to reorder the cities according to my
space-filling curve.  This only makes sense for geometric instances,
and is currently only implemented for 2-dimensional Euclidean instances.
This uses extra memory to remember the mapping.
However, I expect this to use the cache hierarchy more effectively.
This is a research hypothesis.

Option \type{--maxdepth} takes an optional integer argument, $d$.
It asks us to limit the probe depth to at most $d$ generic flips, \ie\
$d$ flips deep in the greedy portion of the search for an improving
sequential $\lambda$-change.
Some researchers speed up their implementations of Lin-Kernighan
by specifiying that probes should stop at 50 or 100 flips after the
backtracking depth.
I believe
this is unnatural.    (Or I used to, rather$\ldots$)
Instead I think limiting the depth due to the
inherent topology of the instance is preferable.
This is a research hypothesis.
If this option is not specified, then no limit is put on the probing depth.
If \type{--maxdepth} is specified without giving a value for $d$, then
a default of $d=50$ is used.

Option \type{-l} or \type{--lower-bound}
takes two mandatory arguments, a string and a number.
The number is taken as a lower bound on the tour length, useful for
keeping track of milestones.  The string is a name for the bound
given by the user.  In verbose output the milestones are labelled with
the name of the bound. 
For example, the user might indicate a known optimal tour length by
specifying \type{-l optimal 14379}. 
(Optimal tours for \type{lin105} have
length 14379.) 
Then the output might say something like: 
\type{2.5\% above optimal after 1.5 user seconds}.
Simliarly, \type{-l Held-Karp 702.3} might be used to get messages like
\type{2.5\% above Held-Karp after 5.5 user seconds}.
If a bound is specified, milestones are printed when the verbose level 
is 25 or more.

Option \type{-u} or \type{--upper-bound}
is similar to \type{-l}, 
taking two mandatory arguments, a string and a number.
The number is taken as an upper bound on the tour length,
and the name is taken as its description, \eg, \type{optimal}, or
\type{best-known}.  

Option \type{--held-karp} specifies that we should compute the
Held-Karp lower bound on the tour length, and then exit.
See module \module{ASCEND} for a description of the Held-Karp lower bound.

Option \type{--held-karp-lambda} specifies that we should compute
the Held-Karp lower bound on tour length and print out a representation
of the Lagrange multipliers with which we achieve that bound.

Option \type{-m} or \type{--matching} tells us to tackle the Minimum
Weighted Perfect Matching problem rather than the Traveling Salesman
Problem.  It affects several of the other options: the tour representation
is irrelevant; the start kind (\type{canonical}, \type{greedy},
\type{random}) specifies algorithms constructing weighted
perfect matchings that correspond to similar algorithms for
constructing tours.

Option \type{-} specifies that the TSPLIB input will be found on the standard
input stream.  

@@ 
Options \type{--extra-backtrack} and \type{--no-extra-backtrack}
require a bit more discussion.

Johnson and McGeoch say the following about backtracking (page
257):\emphpar{Note that for some choices of $t_1$ through $t_5$ there
will be two possible ways of generating a corresponding 3-Opt move 
(two potential candidates for $t_6$ that yield legal 3-Opt moves), and
both are considered.  For other choices, there may be no way to generate
a legal 3-Opt move, and in these cases we attempt to
find cities $t_6$, $t_7$, and $t_8$ that together with $t_1$ through
$t_5$ yield a legal 4-Opt move that meets the one-tree restriction
(with $t_7$  restricted to those cities on the
neighbor list for $t_6$).  The first such move found is used to produce the
starting tour.}

Until 1998/8/6 I had missed the sense of the last sentence and backtracked
over all normal possibilites for $t_6$, $t_7$ and $t_8$. 
This potentially does a bunch more work than the implementation described
by Johnson and McGeoch.  
Since I want my code to be comparable,
I've now made their behaviour the default.  Their behaviour can
be specified explicitly
by option \type{--no-extra-backtrack}.  The old standard behaviour
of my code, more backtracking, is specified by option
\type{--extra-backtrack}.

A fine point needs some attention.
By my case analysis, it appears that some relative orderings of
the first five $t$ cities can lead to a legal 3-change, and no legal
3-change, depending on the choice of $t_6$.
For example,  with relative tour ordering 12435 we can place 6 immediately
before or immediately after 5.  The first case, 124365,  is a legal
3-change for the TSP.
The second case, 124356, is not, but can be extended to legal 4-change
TSP moves: 12874356, 12784356, 12437856, and 12438756.

So, in my humble opinion, predicating the split behaviour on only
the relative positions of $t_1$ through $t_5$ doesn't quite make sense:
the legality of the available 3-Opt moves isn't completely determined
by then.  So I'll take ``their behaviour'' to mean taking only the
first choice of $t_7$ and $t_8$, given $t_1$ through $t_5$ {\it and} $t_6$.


@@ Anything that remains on the command line after options are processed will
be treated as a filename from which our TSPLIB input is taken.  It is an
error to specify more than one file (standard input counts as a file).
If no files are specified, then the standard input is used.

@@ We scan the command line arguments from left to right.

@@<Parse the command line@@>=
@@<Set the option defaults@@>@@;
{ 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 ) {
			@@<Parse this option@@>@@;
			@@<Flag illegal option@@>@@;
		} else {
			if ( filecount ) {
				@@<Print the command line; split at |r|@@>@@;
				fprintf(stderr,"Only one input file allowed\n");
				exit(1);
			}
			if ( !(more_options && strcmp(argv[r],"-")==0) ) 
				filename = argv[r];
				/* |else| case means use |stdin| */
			filecount++;
		}
	}
}
@@<Check option prerequisites@@>@@;


@@
@@<Module variables@@>=
static char *filename;
static int mst_only;
static int do_weighted_perfect_matching;
static int held_karp_only;
static int held_karp_lambda_only;

@@
@@d DEFAULT_VERBOSE 5
@@<Set the option defaults@@>=
	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; /* Dummy. */
	upper_bound_name = NULL;
	upper_bound_value=0.0; /* Dummy. */
	sort = (void (*)(void *, size_t, size_t, int(*)(const void*,const void*)))qsort;
	noround = 0;