/* Richard A. DeVenezia
* www.devenezia.com
*
* Original posted to SAS-L Jan. 9, 2004
*/
/*
* swanson_d wrote:
* Subject: Traveling Salesman Problem
* Newsgroups: comp.soft-sys.sas
* Date: 2004-01-08 04:02:52 PST
*
* I’m looking for a macro similar to the TSP.
*
* I have a list of n cities along with their geographic locations, i.e.,
* their (x,y) coordinates. I want to find the route I should take in which I
* always travel to the nearest un-visited city. Does anybody have a macro
* that does it?
*
* I don’t know if it makes any difference, but I want to measure distances
* with the L1 distance function instead of the usual L2 Euclidean distance
* function. That is, I want the distance from (x0,y0) to (x1,y1) to be ABS
* (x0-x1) + ABS(y0-y1) instead of SQRT((x0-x1)**2 + (y0-y1)**2). The reason
* is that roads in the midwest form a grid.
*/
/* This sample program generates M cities (at random grid points)
* in a NxN grid. A graph of the salesman path is plotted
*
* An initial city is selected and the path (along the grid)
* to the first nearest (L1 distance) city is selected and a line drawn
* to show the path. Furthermore, the lines between cities are
* gradient colored to allow an easier visualation. The number of colors
* are restricted to 254 due to limitations of SAS/Graph.
*
* Note: if the spiral parameter is set to non-zerp, then the cities are
* laid out in the fashion of a spiral, with the spiral value effecting the
* direction and tightness of winding.
*
* This program should not be construed as trying to find a suitable global
* path for minimizing total distance.
*
* A challenge for the inclined: Adjust the program to select the next N cities
* having a minimal path from the current city. This would involve an exhaustive
* search amongst the C(N_remaining,N) possible paths. Laying out the paths
* amongst the N chosen is yet another challenge.
*/
%macro L1_Salesman (
nCities=1520
, nGrid=400
, seed=0
, r0=0, g0=0, b0=255
, r1=0, g1=255,b1=0
, r2=255,g2=0 ,b2=0
, spiral=0
);
data travels ;
* generate city locations;
array x[&nCities] _temporary_;
array y[&nCities] _temporary_;
do citynum = 1 to &nCities;
%if %sysevalf(&spiral = 0) %then %do;
x[citynum] = int(&nGrid * ranuni(&seed)); * 0..&nGrid-1;
y[citynum] = int(&nGrid * ranuni(0)) ; * 0..&nGrid-1;
%end;
%else %do;
theta = citynum / &nCities * Constant('PI') * &spiral ;
r = &nGrid/2 * citynum / &nCities;
x[citynum] = &nGrid/2 + r * cos (theta);
y[citynum] = &nGrid/2 + r * sin (theta);
%end;
end;
* select an initial city;
city0 = 1 + int(&nCities * ranuni(0));
city0 = 1;
x0 = x[city0];
y0 = y[city0];
nLeft = &nCities-1;
* compute distances to untraveled cities
* until all cities traveled to;
do until (nLeft = 0);
mindistance = 1e9;
* determine nearest city;
do i = 1 to &nCities;
if x[i] = . then continue;
if i = city0 then continue;
distance = abs(x[i]-x0) + abs(y[i]-y0);
if distance < minDistance then do;
* record information about the current nearest city;
minDistance = distance;
minI = i;
end;
end;
city1 = minI;
x1 = x[city1];
y1 = y[city1];
distance = minDistance;
output;
* mark city as visited;
x[city0] = .;
y[city0] = .;
* shift values for next distance determination;
city0 = city1;
x0 = x1;
y0 = y1;
nLeft + (-1);
end;
stop;
keep city0 city1 x0 y0 x1 y1 distance;
run;
* create an annotation dataset to draw the path
* traveled by the salesman;
data anno;
set travels nobs=ncities end=end;
length function $8 color $8;
retain xsys ysys '2';
* color gradient path lines for easier eyeball tracking;
* only 255 colors allowed in a sas gplot (news to me!);
if _n_ < ncities/2 then do;
factor = _n_ / ncities * 2; * map [0,0.5] to [0,1];
factor = int(factor*127) / 127; * restrict to 127 levels of alpha;
r = &r0 + factor * (&r1-&r0);
g = &g0 + factor * (&g1-&g0);
b = &b0 + factor * (&b1-&b0);
end;
else do;
factor = _n_ / ncities * 2 - 1; * map [.5,1] to [0,1];
factor = int(factor*127) / 127; * restrict to 127 levels of alpha;
r = &r1 + factor * (&r2-&r1);
g = &g1 + factor * (&g2-&g1);
b = &b1 + factor * (&b2-&b1);
end;
if r<0 then r=0; else if r>255 then r=255;
if g<0 then g=0; else if g>255 then g=255;
if b<0 then b=0; else if b>255 then b=255;
color = "CX" || put (r,hex2.) || put (g,hex2.) || put (b,hex2.);
if priorcolor ne color then ncolors+1;
if ncolors > 254 then color = priorcolor;
priorcolor = color;
retain priorcolor;
if _n_ = 1 then do;
function = 'MOVE'; x = x0; y = y0; output;
function = 'LABEL'; text='S'; size=3; output;
end;
if (0 <= x0 < &nGrid)
and (0 <= x1 < &nGrid)
and (0 <= y0 < &nGrid)
and (0 <= y1 < &nGrid)
then do;
function = 'MOVE'; x = x0; y = y0; output;
function = 'DRAW'; x = x1; y = y1; output;
end;
if end
and (0 <= x1 < &nGrid)
and (0 <= y1 < &nGrid)
then do;
function = 'MOVE'; x = x1; y = y1; output;
function = 'LABEL'; text='E'; size=3; output;
end;
run;
title1 h=2.5pct f="Arial/bold" "Travels of an L1 salesman";
%if %sysevalf(&spiral=0)
%then footnote1 j=r h=2pct "Seed=&SEED, nCities=&nCities, nGrid=&nGrid";
%else footnote1 j=r h=2pct "Seed=&SEED, nCities=&nCities, nGrid=&nGrid, spiral=&spiral";
;
symbol1 v=circle ;
axis1 order=0 to %eval(&nGrid-1)
label=none value=none minor=none major=none ;
proc gplot data=travels anno=anno;
plot y0*x0 / haxis=axis1 vaxis=axis1;
run;
quit;
%mend;
options mprint;
ods listing;
goptions targetdevice=pdf ftext='Arial' htext=12pt hsize=8in vsize=8in;
goptions targetdevice=png ftext='Arial' htext=4pct hsize=1in vsize=1in device=png gsfname=tsp;
filename tsp "\\Extreme\samples\tsp.png";
%L1_Salesman (seed=1, nCities=200, nGrid=100, spiral=0);
filename tsp "\\Extreme\samples\tsp-6.2.png";
%L1_Salesman (seed=1, nCities=200, nGrid=100, spiral=-6.2);
goptions device=win gsfname=;