'Program BRAIN.BAS is an artificial neural network architecture like the brain
'(c) 1998 by J. C. Sprott <http://sprott.physics.wisc.edu/>
'Compile with PowerBASIC 3.5 <http://www.powerbasic.com/>
'March 19, 1998 version
defext a-z                              'Use extended (80-bit) precision
cls
d0=1e-8                                 'Perturbation for LE calculation
input"Number of neurons";n%: if n%=0 then n%=100
mmax%=(fre(-1)-17000-40*n%)/(6*n%)-1    'Check available memory
if n%*mmax%>16383! then mmax%=int(16383!/n%)    'Arrays must fit in 64K
print"Maximum number of inputs is";mmax%
input"Number of inputs to each neuron";m%: if m%=0 then m%=10
input"Scale factor on weights";s: if s=0 then s=2
n$="Random"
screen 12
option base 1                           'Start arrays from 1, not 0
dim x(n%),y(n%),xe(n%),ye(n%),inn%(n%,m%),w!(n%,m%)
randomize timer                         'Reseed the random number generator
open"brainle.dat" for output as #2
do
    cls
    locate 1,64: print n$;" Neighbors";
    t%=0
    for i%=1 to n%
        y(i%)=1-2*rnd                   'Initial conditions
        ye(i%)=y(i%)+d0/sqr(n%)         'Perturbed orbit
        for j%=1 to m%
            if n$="Random" then
                inn%(i%,j%)=rnd(1,n%)   'Input connections
            else
                inn%(i%,j%)=1+(2*n%-m%\2+i%+j%-1) mod n%
            end if
        next j%
        for j%=1 to m%
            w!(i%,j%)=0
            for k%=1 to 12              'Make Gaussian random deviates
                w!(i%,j%)=w!(i%,j%)+rnd
            next k%
            w!(i%,j%)=(w!(i%,j%)-6)/sqr(m%)    'Input weights
        next j%
    next i%
    lsum=0: iter&=0
    while q$<>chr$(27)                  'Iterate net until Esc key is pressed
    	incr iter%
        q$=inkey$
        if len(q$) then                 'A key was pressed
            select case ucase$(q$)
                case chr$(13): exit loop        'Display a new network
                case " ": d%=(d%+1) mod 5: cls  'Toggle display mode
                case chr$(0,72)
                    s=1.1##*s
                    locate 1,1: print using"s =##.###^^^^";s
                    lsum=0: iter&=0
                case chr$(0,80)
                    s=s/1.1##
                    locate 1,1: print using"s =##.###^^^^";s
                    lsum=0: iter&=0
                case "N"
                    if n$="Near" then n$="Random" else n$="Near"
                    exit loop
                case else               'Display menu of commands
                    cls
                    print "Commands:": print
                    print "<Enter>","Display a new network"
                    print "<Space Bar>","Toggle display mode"
                    print "<Up Arrow>","Increase s by 10%"
                    print "<Down Arrow>","Decrease s by 10%"
                    print "<Esc>","Exit program"
                    print "N","Toggle neighborhood (now ";n$;")"
                    do: q$=inkey$: loop until len(q$): cls: iterate
            end select
        end if
        for i%=1 to n%
            x(i%)=y(i%)                 'Update input vector
        next i%
        for i%=1 to n%                  'Iterate the network equations
            u=0
            for j%=1 to m%              'This is the innermost (slowest) loop
                jnn%=inn%(i%,j%)
                u=u+w!(i%,j%)*x(jnn%)
            next j%
            y(i%)=phi(s*u)              'Calculate new output vector
        next i%
        dsum=0
        for i%=1 to n%                  'Do everything again to get the LE
            xe(i%)=ye(i%)               'Update input vector
        next i%
        for i%=1 to n%                  'Iterate the network equations
            u=0
            for j%=1 to m%              'This is the innermost (slowest) loop
                jnn%=inn%(i%,j%)
                u=u+w!(i%,j%)*xe(jnn%)
            next j%
            ye(i%)=phi(s*u)             'Calculate new output vector
            dy=ye(i%)-y(i%)
            dsum=dsum+dy*dy
        next i%
        d1=sqr(dsum)
        for i%=1 to n%                  'Readjust ye for separation d0
            ye(i%)=y(i%)+(d0/d1)*(ye(i%)-y(i%))
        next i%
        lsum=lsum+log(d1/d0)
        incr iter&
        le=lsum/iter&                   'Here's the LE
        locate 30,1: print using"LE =##.####";le;
        if iter%=1250 then
        	print#2, le
            iter%=0
            exit loop
        end if
        t%=(t%+1) mod 640
        if d%=0 then                    'Plot EEG (6 random neurons)
            line(t%,0)-(t%,479),0
            for i%=1 to min(6,n%)
                line(t%-1,80*i%-40+30*x(i%))-(t%,80*i%-40+30*y(i%)),i%+9
            next i%
        elseif d%=1 then                'Plot all neurons with color code
            for i%=1 to 480
                pset(t%-1,i%-1),7.5*(y(int(i%*n%/480+1))+1)
            next i%
        elseif d%=2 then                'Plot Poincare section
            pset(320+320*x(1),240-240*x(2))
        elseif d%=3 then                'Plot average of all neurons
            line(t%,0)-(t%,479),0
            xav=0: yav=0
            for i%=1 to n%: xav=xav+x(i%): yav=yav+y(i%): next i%
            line(t%-1,240-240*xav/n%)-(t%,240-240*yav/n%)
        elseif d%=4 then                'Plot directed graph of network
            call dirgraph(inn%()): screen 12: d%=0
        end if
    wend
loop until q$=chr$(27)
end

FUNCTION phi(u)
    'Returns hyperbolic tangent of u (squashing function)
    phi=1-2/(exp(2*u)+1)
END FUNCTION

SUB dirgraph (inn%())
    'Makes a directed graph of n% nodes each with m% inputs in the array inn%
    n%=ubound(inn%,1): m%=ubound(inn%,2)
    if n%>2020 then exit sub            'Maximum number of nodes is 2020
    dim x%(2020), y%(2020), xm%(2020), ym%(2020)
    indx%=0
    for j%=38 to 0 step -1: for i%=51 to 0 step -1
        incr indx%
        if indx%>2020 then exit for
        xm%(indx%)=14+12*i%: ym%(indx%)=12+12*j%
    next i%,j%
    for i%=1 to 2020                    'Shuffle nodes
        io%=1+int(2020*rnd)
        swap xm%(i%),xm%(io%)
        swap ym%(i%),ym%(io%)
    next i%
    screen 11
    q$=""
    maxpixels&=0
    temp!=1
    mintemp!=1/n%                       'Minimum temperature allowed
    xtmax!=-638/log(mintemp!)
    cf!=exp(-1/xtmax!)                  'Cooling factor
    maxpix&=2232+29*n%                  'Maximum number of pixels
    maxpixels&=-4000000
    minpix&=0
    do
        suml&=0
        cls
        q$=ucase$(inkey$)
        if q$=chr$(27) or q$=" " then exit loop
        if q$="H" then
            temp!=2*temp!
            if temp!>1 then temp!=1 else decr maxpixels&
        end if
        line(0,0)-(639,479),,b
        xt%=1+xtmax!*log(temp!/mintemp!)        'Temperature marker
        preset(xt%,479): preset(xt%-1,479)
        yt%=1+478*(maxpixels&-minpix&)/(maxpix&-minpix&)
        preset(639,yt%): preset(639,yt%-1)
        call arraycopy(xm%(),x%()): call arraycopy(ym%(),y%())
        if q$<>"B" then
            for i%=1 to n%
                if rnd<temp! then
                    io%=1+int(2020*rnd)
                    swap x%(i%),x%(io%)
                    swap y%(i%),y%(io%)
                end if
            next i%
        end if
        temp!=cf!*temp!                 'Simulated annealing schedule
        if temp!<mintemp! then temp!=mintemp!
        for i%=1 to n%: pset(x%(i%),y%(i%)): circle(x%(i%),y%(i%)),5: next i%
        for i%=1 to n%: for j%=1 to m%
            x1%=x%(inn%(i%,j%)): y1%=y%(inn%(i%,j%))
            x2%=x%(i%): y2%=y%(i%)
            dx%=x2%-x1%: dy%=y2%-y1%
            dmax%=max%(abs(dx%),abs(dy%))
            dmin%=min%(abs(dx%),abs(dy%))
            suml&=suml&+dmax%-2
            if 2*dmin%>dmax% and 8*dmin%<5*dmax% then suml&=suml&+2    'Correct
            if 1.143*dmin%>dmax% then suml&=suml&+2                    'diagonals
            'suml&=suml&+dmax%\128      'Penalize long connections
            'if dmax% then suml&=suml&+128\dmax%    'and overly short ones too
            line(x1%,y1%)-(x2%,y2%)
            if q$="B" then
                pset(xt%,479): pset(xt%-1,479)
                pset(639,yt%): pset(639,yt%-1)
                locate 1,1: print file2$;
                l%=sqr(dx%*dx%+dy%*dy%)
                if l% then              'Put arrows on connections
                    xo%=(x1%+x2%)/2: yo%=(y1%+y2%)/2
                    dx%=4*dx%/l%: dy%=4*dy%/l%
                    line(xo%,yo%)-(xo%-2*dx%+dy%,yo%-2*dy%-dx%)
                    line(xo%,yo%)-(xo%-2*dx%-dy%,yo%-2*dy%+dx%)
                else                    'Put loop back to node
                    circle(x1%+9,y1%),9
                end if
            end if
        next j%,i%
        if q$="B" then delay 3: iterate
        pix&=pixels&-suml&
        if minpix&=0 then minpix&=pix&
        if pix&>=maxpixels& then
            call arraycopy(x%(),xm%()): call arraycopy(y%(),ym%())
            temp!=2*temp!               'Add heat if still settling
            if temp!>1 then temp!=1
            if pix&>maxpixels& then     'A better solution was found
                maxpixels&=pix&
                delay 1                 'Admire the solution for a second
            end if
        end if
    loop
END SUB

SUB ArrayCopy (x%(), y%())
    'Copies the integer array x%() to the array y%(i)
    'Only works for arrays with < 16376 elements
    '(ArrayLen% must fit in one string segment of 32750 bytes)
    i1% = LBOUND(x%)
    Elements% = UBOUND(x%) - i1% + 1
    SourceSeg& = VARSEG(x%(i1%))
    SourceOfs& = VARPTR(x%(i1%))
    DestSeg& = VARSEG(y%(i1%))
    DestOfs& = VARPTR(y%(i1%))
    SourceAbs& = SourceSeg& * 16 + SourceOfs&
    DestAbs& = DestSeg& * 16 + DestOfs&
    ArrayLen% = Elements% * 2           'Integers have 2 bytes per element
    DEF SEG = 0
    POKE$ DestAbs&, PEEK$(SourceAbs&, ArrayLen%)
    DEF SEG
END SUB

FUNCTION pixels&
    'Counts the pixels illuminated on the screen (screen mode 11 only)
    LOCAL n1??, n2??
    if PbvScrnMode <> 11 then exit function    'Not a supported screen mode
    asm        push ds
    asm        push si
    asm        mov ax, &HA000
    asm        mov ds, ax
    asm        xor ax, ax
    asm        mov bx, ax
    asm        mov si, ax
    asm        mov cx, &H4B00
    screen_count:
    asm        push cx
    asm        mov cx, &H10
    asm        mov dx, word ptr[si]
    asm        inc si
    asm        inc si
    byte_count:
    asm        shl dx, 1
    asm        jnc around_inc
    asm        add ax, 1
    asm        jnc around_inc
    asm        inc bx
    around_inc:
    asm        loop byte_count
    asm        pop cx
    asm        loop screen_count
    asm        mov n1??, bx
    asm        mov n2??, ax
    asm        pop si
    asm        pop ds
    pixels& = 65536 * n1?? + n2??
END FUNCTION