'******************************************************************************
'Subject: Approximate solution to the Euclidean traveling salesperson problem
' with the Kohonen self-organizing feature map, improved by 2-opt.
'Author : Sjoerd.J.Schaper
'Date : 07-17-2007
'Code : all Q-Basic's, FreeBasic extendable
'******************************************************************************
'This program is copyright (c) 2007 by the author. It is made available as is,
'and no warranty - about the program, its performance, or its conformity to any
'specification - is given or implied. It may be used, modified, and distributed
'freely, as long as the original author is credited as part of the final work.
'******************************************************************************
REM $STATIC
DEFSNG A-E
DEFINT F-Z
TYPE vec
x AS SINGLE
y AS SINGLE
END TYPE
DECLARE FUNCTION infile% (dt)
'read tsp coordinates
'
Link to a sample input file in TSPLIB-format.
DECLARE FUNCTION normlz! (cr)
'normalize input vectors
DECLARE SUB netini (cr)
'initialize neural net
DECLARE FUNCTION match% (c() AS vec, k AS vec, n)
'return best matching unit
DECLARE SUB adapt (m)
'adjust weights relative to bmu
DECLARE FUNCTION oneto1% ()
'1-to-1 city-neuron mapping
DECLARE SUB enlist ()
'build 2-opt neighbour list
DECLARE SUB twoopt ()
'2-opt algorithm
DECLARE SUB revers (p(), r, s, t)
'reverse tour between r and s
DECLARE FUNCTION dlink! (r, s)
'link (r, s) length
DECLARE FUNCTION dtour! ()
'tour length
DECLARE SUB shuffle ()
'random update sequence
DECLARE SUB paints (c() AS vec, sw)
'draw vector net
CONST MODE = 12, dW = 640, dH = 480
' set screenmode, width, height
CONST MX = 50 ' learning loops/city
CONST L0 = .7, L1 = .05 ' learning rate range
CONST R0 = 6, R1 = 1 ' neighbourhood range
CONST MAX = 2147483647 ' large number
CONST STX = 20 ' 2-opt neighbourhood
CONST UB = 2400 ' problem bound (higher w/FB)
DIM SHARED city(UB) AS vec ' city coordinates
DIM SHARED cell(UB) AS vec ' synaptic weights
REDIM SHARED nb(1, 1) AS INTEGER ' 2-opt neighbours
REDIM SHARED ir(1) AS INTEGER ' random indices
REDIM SHARED im(1) AS INTEGER ' modulo indices
REDIM SHARED pt(1) AS INTEGER ' tour pointers
DIM SHARED NC, NN ' number of cities, neurons
RANDOMIZE TIMER
paints city(), 0
DO
PRINT : INPUT "name: ", Ni$
IF Ni$ = "" THEN EXIT DO
IF Ni$ <> "r" THEN
OPEN ".\" + Ni$ + ".tsp" FOR INPUT AS 1
m = infile(dt) ' TSPLIB-file
CLOSE : IF m = 0 THEN EXIT DO
di = normlz(1.05)
IF dt = 0 THEN dt = di * .749 * SQR(NC - 1)
CLS 2: enlist
END IF
'
netini .25: tim# = TIMER
LOCATE 2, 1: PRINT "N"; NC; " ";
FOR m = 1 TO MX
LOCATE 3, 1: PRINT "m"; m; " ";
paints cell(), 1
adapt m
IF INKEY$ = CHR$(27) THEN EXIT FOR
NEXT m
WHILE INKEY$ <> "": WEND
'
IF oneto1 THEN
paints city(), 0
twoopt
paints city(), 0
d = di * dtour
LOCATE 4, 1: PRINT "d"; d; " ";
d = (d - dt) / d ' excess
LOCATE 5, 1: PRINT "e"; d * 100; "% ";
ELSE
PRINT : PRINT "fail: 1-to-1";
END IF
LOCATE 6, 1: PRINT "t"; CSNG(TIMER - tim#); "s "
LOOP
SYSTEM
SUB adapt (m)
STATIC ar, br, al, bl
IF m = 1 THEN ' decay parameters
c0 = R0 + NN * .04
ar = MX * R1
br = ar / c0 - 1
al = (L1 - L0) / (MX - 1)
bl = L0 - al
'
i = CINT(c0): n = NN
IF i > n THEN n = i
REDIM im(-i TO n + i) AS INTEGER
FOR r = -i TO n + i ' indices mod NN
im(r) = (r + n) MOD n
NEXT r
END IF
DIM k AS vec
cr = ar / (m + br) ' compute neighbourhood,
bh = al * m + bl ' learning rate factor
ah = -bh / cr: cr = INT(cr)
'
shuffle
FOR n = 0 TO NC - 1
k = city(ir(n))
r = match(cell(), k, NN) ' best match
FOR i = -cr TO cr ' neighbouring neurons
s = im(r + i)
ch = ah * ABS(i) + bh
d = ch * (k.x - cell(s).x) ' update weights
cell(s).x = cell(s).x + d
d = ch * (k.y - cell(s).y)
cell(s).y = cell(s).y + d
NEXT i
NEXT n
END SUB
FUNCTION dlink! (r, s)
dx = city(pt(r)).x - city(pt(s)).x
dy = city(pt(r)).y - city(pt(s)).y
dlink = SQR(dx * dx + dy * dy)
END FUNCTION
FUNCTION dtour!
d = dlink(NC - 1, 0)
FOR r = 0 TO NC - 2
d = d + dlink(r, r + 1)
NEXT r
dtour = d
END FUNCTION
SUB enlist
REDIM nb(NC - 1, STX) AS INTEGER ' 2-opt neighbours
DIM ds(STX) AS SINGLE ' distances
LINE (0, 0)-(1, 1), 0, BF
FOR r = 0 TO NC - 1
LOCATE 2, 1: PRINT " "; r; " ";
a = city(r).x: b = city(r).y ' current city
IF r THEN
LINE -(a, b), 7
ELSE
PSET (a, b), 7
END IF
CIRCLE (a, b), .003, 12
FOR i = 0 TO STX
ds(i) = MAX
NEXT i
'
j = 0
FOR s = 0 TO NC - 1
dx = city(r).x - city(s).x
dy = city(r).y - city(s).y
d = dx * dx + dy * dy
FOR i = 0 TO j ' truncated list
IF d < ds(i) AND r <> s THEN
FOR k = j TO i + 1 STEP -1
ds(k) = ds(k - 1)
nb(r, k) = nb(r, k - 1)
NEXT k ' insert s
ds(i) = d: nb(r, i) = s
IF j < STX THEN j = j + 1
EXIT FOR
END IF
NEXT i
NEXT s
NEXT r
END SUB
FUNCTION infile% (dt)
NC = 0: NN = 6: dt = 0
DO
LINE INPUT #1, g$ ' TSPLIB-format
g$ = LTRIM$(RTRIM$(g$))
IF (g$ = "EOF") OR EOF(1) THEN EXIT DO
'
s = INSTR(g$, ":")
IF INSTR(g$, "DIMENSION") THEN
NC = VAL(MID$(g$, s + 1))
ELSEIF INSTR(g$, "SOMNETSIZE") THEN
NN = VAL(MID$(g$, s + 1)) ' w/two extra records
ELSEIF INSTR(g$, "TOURLENGTH") THEN
dt = VAL(MID$(g$, s + 1))
ELSEIF g$ = "NODE_COORD_SECTION" THEN
FOR r = 0 TO NC - 1
LINE INPUT #1, g$
g$ = LTRIM$(RTRIM$(g$))
IF (g$ = "EOF") OR EOF(1) THEN EXIT FOR
s = 1 + INSTR(g$, " ")
g$ = LTRIM$(MID$(g$, s))
s = INSTR(g$, " ")
city(r).x = VAL(MID$(g$, 1, s - 1))
city(r).y = VAL(MID$(g$, s + 1))
NEXT r
IF r < NC THEN NC = r
EXIT DO
END IF
'
IF NC > UB OR NN > UB THEN ' out of bounds
PRINT "too big tsp": NC = 0: EXIT DO
END IF
LOOP
infile = NC > 2
END FUNCTION
FUNCTION match% (c() AS vec, k AS vec, n)
dmin = MAX
FOR r = 0 TO n - 1
dx = c(r).x - k.x
dy = c(r).y - k.y
d = dx * dx + dy * dy ' error difference
IF d < dmin THEN
dmin = d: match = r ' winning unit
END IF
NEXT r
END FUNCTION
SUB netini (cr)
n = NC: IF NN > n THEN n = NN
REDIM ir(n - 1) AS INTEGER
REDIM pt(n - 1) AS INTEGER
FOR r = 0 TO n - 1 ' random, city pts
ir(r) = r: pt(r) = r
NEXT r
c = 8 * ATN(1) / NN
FOR r = 0 TO NN - 1 ' neural polygon
a = c * r
cell(r).x = .5 + cr * COS(a)
cell(r).y = .5 + cr * SIN(a)
NEXT r
END SUB
FUNCTION normlz! (cr)
DIM d, dl AS vec, dr AS vec
dl.x = MAX: dl.y = MAX
FOR r = 0 TO NC - 1 ' map size
IF city(r).x < dl.x THEN dl.x = city(r).x
IF city(r).x > dr.x THEN dr.x = city(r).x
IF city(r).y < dl.y THEN dl.y = city(r).y
IF city(r).y > dr.y THEN dr.y = city(r).y
NEXT r
'
c = dr.x - dl.x
dr.x = dl.x + c / 2 ' center
IF c > d THEN d = c
c = dr.y - dl.y
dr.y = dl.y + c / 2
IF c > d THEN d = c
'
d = d * cr ' scale
FOR r = 0 TO NC - 1 ' normalize
city(r).x = .5 + (city(r).x - dr.x) / d
city(r).y = .5 + (city(r).y - dr.y) / d
NEXT r
normlz = d
END FUNCTION
FUNCTION oneto1%
DIM fl(NC - 1) AS INTEGER
FOR s = 0 TO NC - 1
fl(s) = -1 ' free city-flags
NEXT s
'
N0 = NN: r = 0
WHILE r < NN
s = match(city(), cell(r), NC)
IF fl(s) THEN
fl(s) = 0
cell(r) = city(s) ' assign neuron r to nearest city
pt(r) = s: r = r + 1
ELSE
NN = NN - 1
FOR i = r TO NN - 1 ' delete free neuron
cell(i) = cell(i + 1)
NEXT i
END IF
WEND
'
shuffle
FOR n = 0 TO NC - 1
s = ir(n)
IF fl(s) THEN ' insert neuron before best match
r = match(cell(), city(s), NN)
FOR i = NN TO r + 1 STEP -1
cell(i) = cell(i - 1)
pt(i) = pt(i - 1)
NEXT i
cell(r) = city(s) ' assign
pt(r) = s: NN = NN + 1
'
i = (r - 1 + NN) MOD NN ' locally shortest path
t = r + 1: j = (r + 2) MOD NN
d0 = dlink(i, r) + dlink(t, j)
d1 = dlink(i, t) + dlink(r, j)
IF d1 < d0 THEN
SWAP cell(r), cell(t)
SWAP pt(r), pt(t)
END IF
END IF
NEXT n
'
oneto1 = NN = NC
NN = N0
END FUNCTION
SUB paints (c() AS vec, sw)
STATIC fl
IF NOT fl THEN
SCREEN MODE
r = 1 + dW * .25 * .618
s = r + dW * .75 - 1
VIEW (r, 1)-(s, dH - 2)
WINDOW (0, 0)-(1, 1)
fl = -1: EXIT SUB ' initialize only
END IF
n = NC: k = 12
IF sw THEN n = NN: k = 9
WAIT &H3DA, 8
LINE (0, 0)-(1, 1), 0, BF
a = c(pt(0)).x: b = c(pt(0)).y
PSET (a, b), 7
CIRCLE (a, b), .003, k
FOR r = 1 TO n - 1 ' draw connexions
a = c(pt(r)).x: b = c(pt(r)).y
LINE -(a, b), 7
CIRCLE (a, b), .003, k
NEXT r
'
IF sw THEN
FOR r = 0 TO NC - 1 ' draw cities
CIRCLE (city(r).x, city(r).y), .003, 12
NEXT r
END IF
END SUB
SUB revers (p(), r, s, t)
i = r: j = s
FOR k = 0 TO t \ 2
SWAP pt(i), pt(j)
p(pt(i)) = i: i = im(i + 1)
p(pt(j)) = j: j = im(j - 1)
NEXT k
END SUB
SUB shuffle
FOR r = 0 TO NC - 1 ' switch indices
s = INT(RND * NC)
SWAP ir(r), ir(s)
NEXT r
END SUB
SUB twoopt
n2 = NC * 2 - 1
REDIM im(-NC TO n2) AS INTEGER
DIM ip(NC - 1) AS INTEGER
FOR r = -NC TO n2 ' indices mod NC
im(r) = (r + NC) MOD NC
NEXT r
FOR r = 0 TO NC - 1
ip(pt(r)) = r ' inverse pointers
NEXT r
n2 = NC \ 2
DO
sw = 0: shuffle
FOR r = 0 TO NC - 1
t1 = ir(r): dx = 0
t2 = im(t1 + 1)
d0 = dlink(t1, t2) ' optimize link (t1, t2)
FOR s = 0 TO STX
t3 = ip(nb(pt(t2), s)) ' list entry
d1 = dlink(t2, t3)
IF d1 > d0 THEN EXIT FOR
t4 = im(t3 - 1)
'
FOR k = 0 TO 1
d = d0 - dlink(t1, t4)
d = d + dlink(t4, t3) - d1
IF d > dx THEN ' maximum gain
w = (t1 = t3) OR (t3 = t2)
w = w OR (t4 = t2) OR (t1 = t4)
IF NOT w THEN dx = d: u = t3: v = t4
END IF
IF k THEN EXIT FOR
t4 = t3: t3 = im(t4 + 1)
d1 = dlink(t2, t3)
NEXT k
NEXT s
'
IF dx > 1E-08 THEN
t3 = u: t4 = v
w = im(t4 - t2 - 1) ' segment length
IF w < n2 THEN
u = t2: v = t4
ELSE
w = im(t1 - t3 - 1)
u = t3: v = t1
END IF
revers ip(), u, v, w ' improving move
sw = -1
'
u = pt(t1): v = pt(t2)
LINE (city(u).x, city(u).y)-(city(v).x, city(v).y), 9
u = pt(t4): v = pt(t3)
LINE (city(u).x, city(u).y)-(city(v).x, city(v).y), 9
END IF
NEXT r
LOOP WHILE sw
END SUB