Contributor: PIERRE TOURIGNY
program citytour;
{modified anneal.pas in nrpas13.zip from the book Numerical Recipes}
{95-06-23 by Pierre Tourigny, pierre@panpan.synapse.net}
{solves the traveling salesman's problem by simulated annealing}
uses graph;
const
title = 'Tour of 48 US Capitals';
maxcity = 99; {maximum city index}
maxstep = 100; {maximum number of temperature steps}
tfactor = 0.9; {temperature step is 90% of previous step}
radius = 3956.67; {average radius of the earth in miles}
bgidir = 'c:\bp\bgi'; {directory of bgi drivers; change at need}
type
city = record name: string[80]; lat,long: real; gx,gy: integer end;
tdistances = array[0..maxcity,0..maxcity] of real;
var
cities: array[0..maxcity] of city;
tour: array[0..maxcity] of integer;
path,delta,temperature,lowerbound: real;
count,success,maxpath,maxsuccess,offcut,offlast,offsplice,
first,last,cutafter,cutbefore,spliceafter,splicebefore,
mx,my,grdriver,grmode: integer;
distance: ^tdistances;
procedure usage;
begin
writeln;
writeln('Usage: citytour ');
writeln;
writeln('Example: citytour us_cptl.dat us_cptl.sol');
writeln;
writeln('Data File Format:');
writeln(' at most 100 cities, one city per line;');
writeln(' line format: ,,'+
',');
writeln(' latitude and longitude: degrees, minutes and hemisphere');
writeln(' Example: Salt Lake City,40 45 n,111 58 w,Utah,US');
halt;
end;
function dmc2r(degree,minute: real; hemi: string): real;
{degrees, minutes, hemisphere to radians}
begin
if hemi[1] in ['s','S','w','W'] then
dmc2r := -(degree + minute/60.0)*pi/180.0
else dmc2r := (degree + minute/60.0)*pi/180.0;
end;
function gcd(i1,i2: integer): real;
{approximate great-circle distance in miles}
var
dlat,dlong: real;
begin
with cities[i1] do begin
dlat := abs(lat-cities[i2].lat);
if dlat> pi then dlat := 2*pi-dlat;
dlong := abs(long-cities[i2].long);
if dlong> pi then dlong := 2*pi-dlong;
gcd := radius*sqrt(sqr(dlat)+
abs(sqr(dlong)*cos(lat)*cos(cities[i2].lat)));
end;
end;
function split(var s : string; d : string) : string;
var
i,minpos,dpos : byte;
maxs : byte absolute s; maxd : byte absolute d;
begin
minpos := succ(maxs);
for i := 1 to maxd do begin
dpos := pos(d[i],s);
if (dpos> 0) and (dpos < minpos) then minpos := dpos; end; split := copy(s,1,pred(minpos)); s := copy(s,minpos+1,255); end; function fval(s: string): integer; var i,code: integer; begin val(s,i,code); fval := i end; procedure findlow; {find lower bound for optimal solution} {note: the lower bound is not the optimal length; it only means that the best tour cannot be shorter} var i,j: integer; mind1,mind2: real; begin lowerbound := 0.0; for i := 0 to count-1 do begin mind1 := 1.0e38; mind2 := mind1; for j := 0 to count-1 do begin if i = j then continue; if distance^[i,j] < mind1 then begin mind2 := mind1; mind1 := distance^[i,j]; end else if distance^[i,j] < mind2 then mind2 := distance^[i,j]; end; lowerbound := lowerbound+mind1+mind2; end; lowerbound := lowerbound/2; end; procedure drawtour; var i,j: integer; s,sn: string; begin cleardevice; outtext(title); setcolor(white); with cities[0] do moveto(gx,gy); j := 0; for i := 0 to count-1 do begin j := tour[j]; with cities[j] do lineto(gx,gy); end; setcolor(red); setfillstyle(solidfill,red); j := 0; for i := 0 to count-1 do begin with cities[j] do fillellipse(gx,gy,2,2); j := tour[j]; end; setcolor(white); str(path:8:0,sn); s := 'Path Length:'+sn; str(temperature:8:0,sn); s := s+' Temperature:'+sn; str(success:6,sn); s := s+' Successful Moves:'+sn; outtextxy(10,my-15,s); end; procedure init; var i,j: integer; data: text; s: string; gxscale,gxmax,gyscale,gymax: real; begin if paramcount 2 then usage;
randomize;
fillchar(cities,sizeof(cities),0);
fillchar(tour,sizeof(tour),0);
assign(data,paramstr(1));
{$I-} reset(data); {$I+}
if ioresult 0 then begin
writeln('Data file not found'); halt; end;
grdriver := detect;
initgraph(grdriver,grmode,bgidir);
mx := getmaxx; my := getmaxy;
gyscale := my*180.0/(25.0*pi); {map is 25 degrees tall}
gxscale := mx*180.0/(60.0*pi); {map is 60 degrees wide}
gymax := 50.0*pi/180.0; gxmax := 125.0*pi/180.0; {graph point 0,0}
count := 0;
while not eof(data) do begin
readln(data,s);
with cities[count] do begin
name := split(s,',');
lat := dmc2r(fval(split(s,' ')),fval(split(s,' ')),split(s,','));
long := dmc2r(fval(split(s,' ')),fval(split(s,' ')),split(s,','));
gx := trunc((gxmax+long)*gxscale);
gy := trunc((gymax-lat)*gyscale);
end;
inc(count);
end;
close(data);
if (count = 0) or (count> 100) then usage;
for i := 0 to count-1 do tour[i] := (i+1) mod count; {following city}
maxpath := 100*count; maxsuccess := 10*count; path := 0.0;
temperature := 1500.0 {miles};
new(distance); fillchar(distance^,sizeof(distance^),0);
for i := 0 to count-2 do for j := i+1 to count-1 do begin
distance^[i,j] := gcd(i,j); distance^[j,i] := distance^[i,j];
end;
findlow;
for i := 0 to (count-1) do path := path+distance^[i,(i+1) mod count];
drawtour;
end;
procedure findedges;
begin
offcut := random(count); {offset from 0}
offlast := random(count-2); {offset from first}
offsplice := random(count-offlast-2); {offset from cutbefore}
cutafter := 0;
while offcut> 0 do begin
cutafter := tour[cutafter];
dec(offcut);
end;
first := tour[cutafter]; {first city of segment}
last := first; {last city of segment}
while offlast> 0 do begin
last := tour[last];
dec(offlast);
end;
cutbefore := tour[last];
end;
procedure revcost; {difference in path length if segment reversed}
begin
delta := -distance^[cutafter,first]-distance^[last,cutbefore]
+distance^[cutafter,last]+distance^[first,cutbefore];
end;
procedure reverse;
var
i,next,after: integer;
begin
tour[cutafter] := last;
i := first; after := cutbefore;
repeat
next := tour[i];
tour[i] := after;
after := i;
i := next;
until i = cutbefore;
end;
procedure transcost; {difference due to transporting segment}
begin
delta :=
-distance^[cutafter,first]-distance^[last,cutbefore]
-distance^[spliceafter,splicebefore]+distance^[cutafter,cutbefore]
+distance^[spliceafter,first]+distance^[last,splicebefore];
end;
procedure transport;
begin
tour[cutafter] := cutbefore;
tour[spliceafter] := first;
tour[last] := splicebefore;
end;
function metropolis: boolean;
begin
if delta < 0 then metropolis := true else metropolis := random < exp(-delta/temperature); end; procedure anneal; var i,j,k: integer; begin for j := 1 to maxstep do begin success := 0; for k := 1 to maxpath do begin findedges; if random(2) = 1 then begin {reverse segment} revcost; if metropolis then begin inc(success); path := path+delta; reverse; end end else begin {transport segment somewhere else} spliceafter := cutbefore; while offsplice> 0 do begin
spliceafter := tour[spliceafter];
dec(offsplice);
end;
splicebefore := tour[spliceafter];
transcost;
if metropolis then begin
inc(success); path := path+delta; transport;
end;
end;
if success>= maxsuccess then break;
end;
drawtour;
temperature := temperature*tfactor;
if success = 0 then break;
end;
end;
procedure report;
var
i: integer; d,p: real; solution: text;
begin
assign(solution,paramstr(2));
rewrite(solution);
i := 0; p := 0.0;
writeln('City':20,'Distance (miles)':20,'Cumulative':20);
writeln(solution,'City':20,'Distance (miles)':20,'Cumulative':20);
writeln(solution);
repeat
d := distance^[i,tour[i]]; p := p + d;
writeln(cities[i].name:20,d:20:0,p:20:0);
writeln(solution,cities[i].name:20,d:20:0,p:20:0);
i := tour[i];
until i = 0;
writeln('Lower bound on optimal solution: ',lowerbound:8:0);
writeln(solution,'Lower bound on optimal solution: ',lowerbound:8:0);
close(solution);
end;
begin
init;
anneal;
closegraph;
report;
dispose(distance);
end.