att2upper.pl.in,v

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

head	1.1;
access;
symbols
	zero-five-zero:1.1
	zero-four-seventeen:1.1
	zero-four-ten:1.1
	zero-four-nine:1.1
	zero-four-eight:1.1
	zero-four-five:1.1;
locks
	neto:1.1; strict;
comment	@# @;


1.1
date	98.08.08.00.07.05;	author neto;	state Exp;
branches;
next	;


desc
@Convert from various formats to upper row format.
@


1.1
log
@Initial revision
@
text
@#! @@PERL@@ -w
# @@configure_input@@
# vi:ts=4 sw=4:

# att2upper.pl
# Convert TSPLIB ATT instances to upper row-wise distance matrix format.

# This file is in the public domain, and comes with no warranty.
# David Neto, November 22, 1997.

# $Id$
# $Log$

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

my($float_expr)="(-?\\d+\\.?\[0-9\]*(?:\[eE\]\[+-\]\\d+)?|\\.\[0-9\]+(?:\[eE\]\[+-\]\\d+)?)";
my($line);  #Input line
my(@@comments)=();
my(@@buf)=(); # lines buffered for output, pending opening of output stream.
my($n); # Number of vertices
my($type)=0;
my(%type_map)=("EUC_2D",1,"CEIL_2D",2,"ATT",3,"GEO",4,"EXPLICIT",5);
my($type_string);

# 1-based array of coordinates
my(@@coordx)=();
my(@@coordy)=();

my($usage) = <<EOT;
$version_banner
Convert TSPLIB instances into DIMACS matching format
$progname [options]
  -h --help       : Print this help and exit successfully
     --version    : Print version info and exit successfully
EOT

# 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/^--version$/) { print "$version_banner\n"; exit; }
    die "$progname: Unknown option $option\n$usage";
}




############################################
# Form the output.

HEADER: while($line=<>) {
	$_ = $line;
	if (m/^\s*NAME\s*:/) {
		push(@@buf,$line);
	} elsif (m/^\s*COMMENT\s*:/) {
		chop $line;
		push(@@buf,$line." | att2upper.pl\n");
	} elsif (m/^\s*TYPE\s*:\s*(.*)/) {
		$1 =~ m/^TSP/ || die "I know TSPLIB files of type TSP, not $1";
		push(@@buf,$line);
	} elsif (m/^(\s*EDGE_WEIGHT_TYPE\s*:\s*)([^\s]*)\s*/) {
		#print STDERR "edge weight type .$2.\n";
		$type_string = $2;
		$type = $type_map{$type_string};
		$type > 0
			|| die "TSP is of type $type_string, not ATT, CEIL_2D, EUC_2D, or GEO.  Not converted.\n";
	} elsif (m/^\s*DIMENSION\s*:\s*(.*)/) {
		push(@@buf,$line);
		#print STDERR "dimension .$1.\n";
		$n=0+$1;
	} elsif (m/^\s*NODE_COORD_SECTION\s*/) {
		$type==1 || $type==2 || $type==3 || $type==4
			|| die "Can't have NODE_COORD_SECTION in $type_string";
		last HEADER;
	} elsif (m/^\s*DISPLAY_DATA_TYPE\s*/) {
		;
	} else { die "Unrecognized line: $line"; }
}

push(@@buf, "EDGE_WEIGHT_TYPE: EXPLICIT\n");
push(@@buf, "EDGE_WEIGHT_FORMAT: UPPER_ROW\n");
push(@@buf, "EDGE_WEIGHT_SECTION\n");
print @@buf;

SWITCH: {
	&two_d, 	last SWITCH if ($type==1 || $type==2 || $type==3);
	&geo, 		last SWITCH if $type==4;
	&explicit, 	last SWITCH if $type==5;
	die "Unkown TSPLIB type";
}

print "EOF\n";
exit 0;

sub two_d {
	# Read the coordinates and store them into @@coordx and @@coordy.
	my($i)=0;
	my($j)=0;
	while (($i<$n) && ($line=<>) ) {
		if ( $line=~ m/^\s*(\d+)\s+$float_expr\s+$float_expr/o ) {
			my($k)=0+$1;
			$coordx[$k] = 0+$2;
			$coordy[$k] = 0+$3;
			$i++;
		}
	}
	$i==$n || die "Not enough vertices on input";

	# Now write out the matrix.
	for ( $i=1;$i<$n;$i++) {
		my($d);
		for ( $j=$i+1;$j<=$n;$j++) {
			my($dx)=$coordx[$i]-$coordx[$j];
			my($dy)=$coordy[$i]-$coordy[$j];
			my($sumsq)=$dx*$dx+$dy*$dy;
			ROUND: {
				$d=int(0.5+sqrt($sumsq)),		last ROUND if $type==1; #EUC_2D
				$d=ceil(sqrt $sumsq),			last ROUND if $type==2; #CEIL_2D
				$d=ceil(sqrt ($sumsq/10.0)),	last ROUND if $type==3; #ATT
			}
			printf "%.0f ", $d;
		}
		print "\n";
	}
}

sub geo {
	# Read the coordinates and store them into @@lat and @@long.
	# The coordinates will be latitude and longitude.
	my($i)=0;
	my($j)=0;
	my($q1,$q2,$q3);
	my($DEG2RAD)=3.1415926535897932384626433832795028841972/180.0;
#	my($pi) = 3.141592; # according to TSPLIB doc.
	my($RRR)=6378.388;  	# according to TSPLIB doc.
	my($deg,$mn);
	my(@@long)=();	# 1-based longitude array
	my(@@lat)=();	# 1-based latitude array
	my($deg_expr)="(\[-+\]?)(\\d+)(\\.\\d+)";
	while (($i<$n) && ($line=<>) ) {
		if ( $line=~ m/^\s*(\d+)\s+$deg_expr\s+$deg_expr/o ) {
			my($k)=$1;
			my($sgn1)=$2;
			my($deg1)=$3;
			my($min1)=$4;
			my($sgn2)=$5;
			my($deg2)=$6;
			my($min2)=$7;
			$sgn1=($sgn1 eq "-" ? -1 : 1);
			$sgn2=($sgn2 eq "-" ? -1 : 1);
			$min1 < 0.6 || die "$min1 >= 0.6";
			$min2 < 0.6 || die "$min2 >= 0.6";
			$lat[$k]=$DEG2RAD*$sgn1*(5.0/3.0*$min1+$deg1);
			$long[$k]=$DEG2RAD*$sgn2*(5.0/3.0*$min2+$deg2);
print "$2$3$4 to $lat[$k] $5$6$7 to $long[$k]\n";
			$i++;
		} else {
			print STDERR $line;
			die "Rejected: not in \d+ degrees degrees form.";
		}
	}
	$i==$n || die "Not enough vertices on input";

	# Now write out the matrix.
	for ( $i=1;$i<$n;$i++) {
		my($d);
		for ( $j=$i+1;$j<=$n;$j++) {
#			my($q1)=cos($long[$i]-$long[$j]);
#			my($q2)=cos($lat[$i]-$lat[$j]);
#			my($q3)=cos($lat[$i]+$lat[$j]);
#			printf "%.0f ", 
#				int(1.5+$RRR*arccos(0.5*((1.0+$q1)*$q2-(1.0-$q1)*$q3)));
			# See ftp://ftp.netcom.com/pub/hb/hbaker/FAQ-lat-long.txt
			# Thanks to Henry Baker.
			my($lat1)=$lat[$i];
			my($lat2)=$lat[$j];
			my($dlat)=$lat1-$lat2;
			my($dlong)=$long[$i]-$long[$j];
			my($sdlato2)=sin($dlat/2);
			my($sdlongo2)=sin($dlong/2);
			my($d)=2*arcsin(min(1,sqrt($sdlato2*$sdlato2+cos($lat1)*cos($lat2)*$sdlongo2*$sdlongo2)));
			$d >= 0 || die "$d < 0";
			my($greatdist)=int(0.5+$RRR*$d);
			printf " %.0f",$greatdist;
		}
		print "\n";
	}
}

sub explicit {
	die "Yo David, write this code, will ya!";
}


###############

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

# arccos(x) = arctan((sqrt(1-x^2))/x)
sub arccos {
	my($t)=shift;
	return atan2(sqrt(1-$t*$t),$t);
}

# arcsin(x) = arctan(x/(sqrt(1-x^2)))
sub arcsin {
	my($t)=shift;
	return atan2($t,sqrt(1-$t*$t));
}

# convert to nearest integer, rounding halfway cases to larger magnitude.
sub nint {
	my($x)=shift;
	if ( $x < 0 ) {
		return int($x-0.5);
	} else {
		return int($x+0.5);
	}
}
sub min {
	my($x)=shift;
	my($y)=shift;
	if ( $x < $y ) {return $x; } else {return $y;}
}
@