tspbgen.pl.in

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

#! @PERL@ -w
# @configure_input@
# Generate TSPLIB instances from any of the Bentley distributions.
# See "Fast Algorithms for Geometric Traveling Salesman Problems",
#	Jon Louis Bentley, ORSA Journal on Computing, Vol 4, No 4, Fall 1992.

# This program is Copyright 1997 David Neto.
# It is distributed with the package LK under the terms of the 
# GNU Public License, version 2 or later.
# See the file named COPYING for the details of the license.

# The input distributions are listed in Table 1 of Bentley's paper, 
# reproduced here:
#
# Num  Name     Description
#  0   uni      Uniform within the unit square (U[0,1]^2)
#  1   annulus  Uniform on a circle (a zero width annulus)
#  2   arith    x0= 0, 1, 4, 9, 16, ... (arithmetic differences), x1=0
#  3   ball     Uniform inside a circle
#  4   clusnorm Choose 10 points from U[0,1]^2, ethen put a Normal(0.05) at each
#  5   cubediam x0=x1=U[0,1]
#  6   cubeedge x0=U[0,1]; x1=0
#  7   corners  U[0,1] at (0,0), (2,0), (0,2) (2,2)
#  8   grid     Choose N points froma  square grid that contains about 
#                   1.3N points
#  9   normal   Each dimension independent from Normal(1)
# 10   spokes   N/2 points at (U[0,1],0.5); N/2 at (0.5,U[0,1])


# Include directory hack, to find GB_flip.
push(@INC,"../script");
# Get my Perl clone of the Stanford GraphBase random number generator.
require GB_flip; import GB_flip;


my $progname="tspbgen.pl";
my $pkgname = "@PACKAGE@";
my $pkgversion = "@VERSION@";
my $version_banner="$progname ($pkgname $pkgversion)";

# See Knuth's TAOCP, volume 1.
my $Pi = 3.1415926535897932384626433832795028841972 ;
my $scale = 1000000;	# Bounding box size.

my $d_default = 0;
my $seed_default = 42;
my $N_default = 1000;
my %d_map = ("uni",0,"annulus",1,"arith",2,"ball",3,"clusnorm",4,"cubediam",5,
"cubeedge",6,"corners",7,"grid",8,"normal",9,"spokes",10);
my @basename_map=("uni","annulus","arith","ball","clusnorm","cubediam",
"cubeedge","corners","grid","normal","spokes");

my $usage = <<EOT;
$version_banner
Generate TSPLIB instances from one of the 10 Bentley distributions.
$progname [options]
  where options can be:
    -d <d> --distnum <d>   : Distribution number <d>, 0 through 10, see -D
    -D <d> --distname <d>  : Name of the distribution, <d> is one of 
                             {uni,annulus,arith,ball,clusnorm,cubediam,
                             cubeedge,corners,grid,normal,spokes}
    -h     --help          : Print this message
    -n N   --number N      : Instance has N cities, $N_default
    -s n   --seed n        : Set the random seed to n, default is $seed_default
                             seed 0 forces the identity permutation
           --version       : Print a version info, then exit
EOT

# Set the defaults;
my $d = $d_default;
my $seed=$seed_default;
my $N = $N_default;

# Parse the command line options
while ($#ARGV >= 0 && $ARGV[0] =~ m/^-/) {
    my $option = $_ = shift(@ARGV);
    if (m/^-h$/ || m/^--help$/) { print $usage; exit; }
    if (m/^-s$/ || m/^--seed$/) { # mandatory seed parameter
		if ( $#ARGV >= 0 ) {$seed=0xffffffff & int(shift(@ARGV)); next;}
		else { die "$progname: option $option requires a seed argument\n";}
	}
    if (m/^-d$/ || m/^--distnum$/) { # mandatory number parameter
		if ( $#ARGV >= 0 ) {
			my($raw_d)=shift(@ARGV);
			$d=int($raw_d);
			(0 <= $d && $d <= 10) 
				|| die "$progname: option $option needs integer value 0..10, given $raw_d";
			next;
		} else { die "$progname: option $option requires a seed argument\n";}
	}
    if (m/^-D$/ || m/^--distname$/) { # mandatory name keyword parameter
		if ( $#ARGV >= 0 ) {
			my($name)=shift(@ARGV);
			if ($name eq "uni") { $d=0; }
			else { 
				$d=$d_map{$name} 
				|| die "$progname: option $option needs distribution name argument, given $name";
			}
			next;
		} else { die "$progname: option $option requires a distribution name argument\n";}
	}
    if (m/^-n$/ || m/^--number$/) { # mandatory number parameter
		if ( $#ARGV >= 0 ) {
			my($raw_N)=shift(@ARGV);
			$N=int($raw_N);
			(3 <= $N)
				|| die "$progname: option $option needs integer value at least 3, given $raw_N";
			next;
		} else { die "$progname: option $option requires an integer argument\n";}
	}
    if (m/^--version$/) { print "$version_banner\n"; exit; }
    die "$progname: Unknown option $option\n$usage";
}


############################################
# Form of the output.
gb_init_rand($seed);

print "NAME: $basename_map[$d].$seed.$N\n";
print "TYPE: TSP\n";
print "COMMENT: Bentley $basename_map[$d], seed $seed, n=$N, by $version_banner\n";
print "DIMENSION: $N\n";
print "EDGE_WEIGHT_TYPE: EUC_2D\n";
print "NODE_COORD_SECTION\n";

# Plotting the output to the following provides an eyeball check on unif01.
# Use gnuplot to view the output.
#open(FOO,">foo.gpl");
#for $i (0..20000) {print FOO "$i ",&unif01,"\n";}
#close(FOO);

#print "option d is $d\n";
#print "option N is $N\n";
SWITCH: {
	&uni,      last SWITCH if $d==0;
	&annulus,  last SWITCH if $d==1;
	&arith,    last SWITCH if $d==2;
	&ball,     last SWITCH if $d==3;
	&clusnorm, last SWITCH if $d==4;
	&cubediam, last SWITCH if $d==5;
	&cubeedge, last SWITCH if $d==6;
	&corners,  last SWITCH if $d==7;
	&grid,     last SWITCH if $d==8;
	&normal,   last SWITCH if $d==9;
	&spokes,   last SWITCH if $d==10;
}

print "EOF\n";

exit 0;


############################################
# Random sampling.

sub unif01 { # Uniform sample over [0,1], to at least 53 bits precision.
	my ($a,$b,$quo);
	# Assert an ordering.
	$a = &gb_unif_rand(0x40000000);  # 30 bit random number
	$b = &gb_unif_rand(0x40000000);  # 30 bit random number
	$quo = (1<<30);
	$quo *= (1<<30);
	$quo -= 1;  # This likely has no effect, but theoretically forces 
		# the sampling interval to be closed at 1.
	return ($a*(1<<30) + $b)/$quo;
}


my $norm_saved=0;
my $norm_have_saved=0;
sub norm { # One argument, the standard deviation
	# See Knuth's TAOCP, vol2, sec 3.4.1, algorithm P.
	my($stddev)=shift;
	my($v1,$v2,$s);
	if ( $norm_have_saved ) { $norm_have_saved = 0; return $norm_saved*$stddev; }
	else {
		do { ($v1,$v2) = (2*&unif01-1,2*&unif01-1); 
			$s = $v1*$v1 + $v2*$v2;
		} while ( $s >= 1 );
		$norm_saved = $v2*sqrt(-2*log($s)/$s); $norm_have_saved=1;
		return $v1*sqrt(-2*log($s)/$s)*$stddev;
	}
}


############################################
# Generating individual distributions

sub putline { # three arguments: i, x, y
	# I put this is one spot so that I can easily and globally control output precision.
	my($i,$x,$y)=@_;
	printf "%d %d %d\n",$i,(int($x+0.5)),(int($y+0.5));
}

sub ceil {
	my($x,$ix)=shift;
	$ix = int($x);
	if ( $ix == $x ) {return $ix};
	return $ix+1;
}

sub uni {
	my $i;
	for $i (1..$N) {
		putline($i,gb_unif_rand($scale+1),gb_unif_rand($scale+1));
	}
}

sub annulus {
	my($angle);
	for my $i (1..$N) {
		$angle =  &unif01 * 2 * $Pi;
		putline($i,$scale*cos($angle),$scale*sin($angle));
	}
}

sub arith {
	for my $i (1..$N) {
		putline($i,($i-1)*($i-1),0);
	}
}

sub ball {
	for my $i (1..$N) {
		my($v1,$v2,$s);
		do { ($v1,$v2) = (2*&unif01-1,2*&unif01-1); 
			$s = $v1*$v1 + $v2*$v2;
		} while ( $s >= 1 );
		putline($i,$scale*$v1,$scale*$v2);
	}
}

sub clusnorm {
	my(@x,@y,$xn,$yn);
	for ( $i=0; $i < 10; $i++ ) {
		$x[$i]=&unif01; 
		$y[$i]=&unif01;
	}
	for $i (1..$N) {  # Lousy myity...
		$xn = &norm(0.05);  # Force ordering.
		$yn = &norm(0.05);
		putline($i, $scale*($x[$i%10] + $xn), $scale*($y[$i%10] + $yn));
	}
}

sub cubediam {
	my($v);
	for $i (1..$N) {
		$v=&unif01;  # Perl is not referentially transparent....
		putline($i,$scale*$v,$scale*$v);
	}
}

sub cubeedge {
	my($v);
	for $i (1..$N) {
		putline($i,$scale*&unif01,0);
	}
}

sub corners {
	my(@x,@y);
	@x=(0,2,0,2);
	@y=(0,0,2,2);
	for $i (1..$N) {  # Lousy myity...
		putline($i,$scale*($x[$i%4]+&unif01),$scale*($y[$i%4]+&unif01));
	}
}

sub grid {
	my($total,$side,$omit_num,%omit,$x,$y);
	$side= ceil(sqrt(1.3*$N));
	my $my_scale = int($scale/$side);
	$total = $side*$side;
	$omit_num=$total-$N;
	#print "side=$side total=$total omit_num=$omit_num\n";
	%omit=();
	for ( $i=0; $i<$omit_num; $i++) {
		do {
			$x=int(&unif01*$side);
			$y=int(&unif01*$side);
			#printf "%5d ($x,$y)\n",$i;
		} while( $x >= $side || $y >= $side || $omit{"$x#$y"} );
		$omit{"$x#$y"}=1;
		#print "   omit($x,$y)\n";
	}
	$i=0;
	for ( $x=0; $x < $side; $x++  ) {
		for ( $y=0; $y < $side; $y++  ) {
			do { $i++;	putline($i,$my_scale*$x,$my_scale*$y); } unless $omit{"$x#$y"};
		}
	}
	$i==$N || die "$progname: grid generated $i cities but should have generated $i\n";
}

sub normal {
	for $i (1..$N) { 
		putline($i,$scale*&norm(1),$scale*&norm(1));
	}
}

sub spokes {
	for $i (1..$N) { # Lousy myity...
		if ( $i % 2 ) {
			putline($i,$scale*&unif01,$scale*0.5);
		} else {
			putline($i,$scale*0.5,$scale*&unif01);
		}
	}
}