'******************************************************************************
'Subject: Teuvo Kohonen's 1982 self organizing feature map
'         a toy engine for 2D color clustering
'Author : Sjoerd.J.Schaper
'Date   : 12-16-2005
'Code   : all QBasic's
'Keys   : [Esc] quit program
'         press [Space] to read the next vector set
'******************************************************************************
'This program is copyright (c) 2005 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

DECLARE FUNCTION invecs ()
'read vector file
DECLARE FUNCTION inlink ()
'initialize weights and normalize input
DECLARE SUB shuffle (nr)
'random update sequence
DECLARE FUNCTION bmu% (n)
'return best matching unit
DECLARE FUNCTION adjust% (m, n)
'adjust weights relative to bmu
DECLARE SUB paints (k)
'paint network cell

CONST MODE = 13, dW = 320, dH = 200
' set screenmode, width, height

CONST L0 = .9, L1 = .1             '  learning rate range
CONST R0 = 1, R1 = .1              '  neighborhood range
CONST MIN = .1                     '  small number
CONST MAX = 2147483647             '  big number

CONST IX = 16, IY = 16             '  rectangular Kohonen layer
CONST NO = IX * IY                 '  number of cells
CONST NI = 3                       '  input vector dimension
CONST NX = 14                      '  max. number of input vectors

DIM SHARED wei(NO - 1, NI - 1) AS SINGLE ' weight matrix
DIM SHARED vec(NX - 1, NI - 1) AS SINGLE ' observation vectors
DIM SHARED ic(NX - 1) AS INTEGER         ' vector indices
DIM SHARED M1, NV, ef
'iteration count, number of input vectors, scaling factor

OPEN ".\colors.lif" FOR INPUT AS 1
RANDOMIZE TIMER
paints 0

DO
   NV = invecs
   IF NV = 0 THEN EXIT DO
   ef = inlink

   FOR m = 1 TO M1
      LOCATE 2, 2: PRINT "m"; m;
      shuffle NV

      FOR n = 0 TO NV - 1
         LOCATE 4, 2: PRINT "n"; ic(n) + 1;
         t = adjust(m, n)

         IF t THEN
            IF t = 27 THEN EXIT DO
            GOTO skip
         END IF
      NEXT n
   NEXT m
   SLEEP
   WHILE INKEY$ <> "": WEND

skip:
LOOP

CLOSE
SYSTEM

FUNCTION adjust% (m, n)
STATIC ar, br, al, bl
IF n = 0 THEN
   IF m = 1 THEN
      x = IX - 1: y = IY - 1 '        linear decay parameters
      c0 = R0 * SQR(x * x + y * y) * .5
      c1 = R1 * 29.1
      ar = (c1 - c0) / (M1 - 1)
      br = c0 - ar
      al = (L1 - L0) / (M1 - 1)
      bl = L0 - al
   END IF
END IF
adjust = 0

   cr = ar * m + br '                 compute neighbourhood,
   cr = cr * cr
   bh = al * m + bl '                 learning rate factor
   ah = -bh / cr

   k = ic(n): km = bmu(n) '           best match
   y = km \ IX '                       coordinates
   x = km - y * IX

   FOR v = 0 TO IY - 1
      vi = v * IX
      dy = (v - y) * (v - y)

      FOR u = 0 TO IX - 1
         dt = (u - x) * (u - x) + dy

         IF dt < cr THEN
            r = vi + u
            ch = ah * dt + bh '       quadratic smoothing

            FOR s = 0 TO NI - 1
               d = ch * (vec(k, s) - wei(r, s))
               wei(r, s) = wei(r, s) + d ' update weight
            NEXT s
          
            paints r
         END IF
      NEXT u

      g$ = INKEY$
      IF g$ <> "" THEN
         adjust = ASC(g$): EXIT FOR
      END IF
   NEXT v
END FUNCTION

FUNCTION bmu% (n)
   k = ic(n): dmin = MAX
   FOR r = 0 TO NO - 1

      ds = 0
      FOR s = 0 TO NI - 1
         d = vec(k, s) - wei(r, s) '  error difference
         ds = ds + d * d
      NEXT s
      IF ds < dmin THEN
         dmin = ds: bmu = r '         winning weight
      END IF
   NEXT r
END FUNCTION

FUNCTION inlink
   d = 2 * MIN
   FOR r = 0 TO NO - 1
      FOR s = 0 TO NI - 1 '           small random weights
         wei(r, s) = (RND - .5) * d
      NEXT s
   NEXT r

   dmax = 0
   FOR r = 0 TO NV - 1 '              max. vector length
      d = 0
      FOR s = 0 TO NI - 1
         d = d + vec(r, s) * vec(r, s)
      NEXT s
      IF d > dmax THEN dmax = d
   NEXT r

   d = SQR(dmax)
   FOR r = 0 TO NV - 1 '              normalize
      FOR s = 0 TO NI - 1
         vec(r, s) = vec(r, s) / d
      NEXT s
   NEXT r

inlink = d
END FUNCTION

FUNCTION invecs
invecs = 0: r = 0
   DO
      IF EOF(1) THEN
         invecs = r: EXIT DO
      END IF
      LINE INPUT #1, g$: g$ = RTRIM$(g$)
      '
      SELECT CASE LEFT$(g$, 2)
      CASE "#C", "#D", "#L" '         comments
         IF r THEN
            invecs = r: EXIT DO '     go!
         END IF
      CASE "#N", "#P", "#R" '         iteration count
         M1 = VAL(MID$(g$, 3))
         IF M1 < 1 THEN M1 = 1
      CASE ELSE '                     read color vector
         cl& = VAL(g$)
         FOR s = 0 TO NI - 1
            vec(r, s) = cl& MOD 256
            cl& = cl& \ 256
         NEXT s: r = r + 1
      END SELECT
   LOOP
END FUNCTION

SUB paints (k)
STATIC fl
IF NOT fl THEN
   SCREEN MODE
   cl& = &H282828 '                   initialize gfx
   COLOR 255
   r = 1 + dW * .25 * .618
   s = r + dW * .75 - 1: fl = -1
   VIEW (r, 1)-(s, dH - 2)

   dx = IX - .5: dy = IY - .5
   WINDOW SCREEN (-.5, -.5)-(dx, dy)
   FOR r = 0 TO IY - 1 '              draw attribute grid
       FOR s = 0 TO IX - 1
          t = r * IX + s: PALETTE t, cl&
          LINE (s - .5, r - .5)-(s + .5, r + .5), t, BF
       NEXT s
   NEXT r
END IF

   pw& = 256: cl& = 0
   FOR t = NI - 1 TO 0 STEP -1
      n = CINT(ef * wei(k, t))
      IF n < 0 THEN n = 0
      IF n > 63 THEN n = 63
      cl& = cl& * pw& + n
   NEXT t

   PALETTE k, cl& '                   update cell
END SUB

SUB shuffle (nr)
   FOR r = 0 TO NV - 1
      ic(r) = r
   NEXT r
   FOR t = 1 TO nr '                  switch indices
      r = INT(RND * NV)
      s = INT(RND * NV)
      SWAP ic(r), ic(s)
   NEXT t
END SUB