'Program COLLECT.BAS by J. C. Sprott
'Collects Lyapunov exponents in filename$ for a variety of n, d, and s.
$ERROR ALL OFF
DEFEXT A-Z           'Use extended precision
CLS
filename$ = "LYAP.DAT"
nd& = 250           'Number of points to discard
ns& = 1000          'Number of points to save for calculating LE
nmax% = 128               'Maximum number of units
dmax% = 512              'Maximum number of inputs
smax = 128            'Maximum standard deviation of normal weights
gr% = -1            'Turn graphics on (-1) or off (0)
IF gr% then SCREEN 12 else SCREEN 2
DIM huge x(dmax%), xe(dmax%), w!(nmax%, dmax%), b!(nmax%), nl%(639)
cases& = 0: period& = 0: cycle& = 0: chaos& = 0: nlmax% = 0
FOR i% = 0 TO 639: nl%(i%) = 0: NEXT i%
WHILE q$ <> CHR$(27)    'Collect cases until <Esc> key is pressed
    n%=2^int(ran*log2(2*nmax%)): 'n% = 2
    d%=2^int(ran*log2(2*dmax%)): 'd% = 2
    s=2^int(ran*log2(2*smax)): 's = 8
    q$ = ""
    CLS
    IF gr% THEN
        LINE (320, 459)-(320, 0), , , &H4444
        LINE (0, 0)-(639, 459), , B
        LINE (329,9)-(631,221), , B, &H4444
    END IF
    LOCATE 2, 3: PRINT "N ="; n%;
    LOCATE 3, 3: PRINT "D ="; d%;
    LOCATE 4, 3: PRINT "S ="; s;
    IF gr% THEN LOCATE 30, 1: PRINT "-1";
    IF gr% THEN LOCATE 30, 29: PRINT "Largest Lyapunov Exponent";
    IF gr% THEN LOCATE 30, 80: PRINT "1";
    IF cases& > 0 THEN      'Bin data and print status
        SELECT CASE l
            CASE < -.01: period& = period& + 1
            CASE -.01 TO .01: cycle& = cycle& + 1
            CASE > .01: chaos& = chaos& + 1
        END SELECT
        IF l > -100 THEN i% = 320 * (l + 1) ELSE i% = 0
        IF i% < 0 THEN i% = 0
        IF i% > 639 THEN i% = 639
        nl%(i%) = nl%(i%) + 1
        IF nl%(i%) > nlmax% THEN nlmax% = nl%(i%)
        LOCATE 5, 3: PRINT "Cases:"; cases&;
        LOCATE 6, 3: PRINT "Chaotic:"; CINT(1000 * chaos& / cases&) / 10; "%";
        IF gr% THEN PSET (0, 459)
        IF gr% THEN FOR j% = 0 TO 639: LINE -(j%, 459 - 440 * nl%(j%) / nlmax%): NEXT j%
    END IF
    its& = 0
    btot = 0
    FOR i% = 1 TO n%    'Initialize weights
        b!(i%) = ran
        btot = btot + ABS(b!(i%))
    NEXT i%
    FOR i% = 1 to n%: b!(i%) = b!(i%) / btot: NEXT i%        'make b sum to 1
    FOR i% = 1 to n%
        FOR j% = 0 TO d%
            g = 0
            FOR k% = 1 TO 24: g = g + ran: NEXT k%    'g is normal distribution
            g = 0.707106781 * (g - 12)
            w!(i%, j%) = s * g
        NEXT j%
    NEXT i%
    FOR j% = 1 TO d%: x(j%) = 2 * ran - 1: NEXT j%      'Initialize x (-1,1)
    WHILE its& < nd& + ns&      'Iterate net eqns
        CALL iterated(x(), d%, n%, y)
        xp = 480 + 150 * x(d%): yp = 115 - 105 * y
        FOR i% = 1 TO d% - 1: x(i%) = x(i% + 1): NEXT i%
        x(d%) = y
        IF its& > nd& THEN IF gr% THEN PSET (xp, yp)
        CALL lyapunov(x(), d%, n%, its&, nd&, l)
        CALL sumle(x(), d%, n%, its&, nd&, sl)
        its& = its& + 1
        q$ = INKEY$
        IF q$ = CHR$(27) THEN END
        IF q$ = CHR$(32) THEN CLS     'Clear screen if <Spacebar> is pressed
    WEND
    cases& = cases& + 1
    OPEN filename$ FOR BINARY AS #1
    seek 1,lof(1)
    ni%=log2(n%): di%=log2(d%): si%=log2(s)
    PUT$ #1, chr$(ni%)+chr$(di%)+chr$(si%)+mks$(l)
    CLOSE #1
    OPEN "SUMLYAP.DAT" FOR BINARY AS #1
    seek 1,lof(1)
    PUT$ #1, chr$(ni%)+chr$(di%)+chr$(si%)+mks$(sl)
    CLOSE #1
WEND
END

SUB iterated (x(), d%, n%, y)   'Iterate equations
    'y = 1 - 1.4 * x(2) * x(2) + .3 * x(1): EXIT SUB    'Henon map test
    SHARED w!(), b!()
    y = 0
    FOR i% = 1 TO n%
        u = w!(i%, 0)
        FOR j% = 1 TO d%: u = u + w!(i%, j%) * x(j%): NEXT j%
        y = y + b!(i%) * phi(u)
    NEXT i%
END SUB

SUB lyapunov (x(), d%, n%, its&, nd&, l)     'Get largest Lyapunov exponent
    SHARED xe()
    STATIC lsum, k%
    IF its& = 0 THEN
        xe(1) = x(1) + .00000001#
        FOR i% = 2 TO d%: xe(i%) = x(i%): NEXT i%
        lsum = 0: k% = 0
    END IF
    IF its& = nd& THEN lsum = 0     'Discard the first hundred iterates
    k% = (k% + 1) MOD 50
    CALL iterated(xe(), d%, n%, y)
    FOR i% = 1 TO d% - 1: xe(i%) = xe(i% + 1): NEXT i%
    xe(d%) = y
    dl2 = 0
    FOR i% = 1 TO d%: dlx = xe(i%) - x(i%): dl2 = dl2 + dlx * dlx: NEXT i%
    IF dl2 > 0 THEN
        df = 1D+16 * dl2
        rs = 1 / SQR(df)
        FOR i% = 1 TO d%: xe(i%) = x(i%) + rs * (xe(i%) - x(i%)): NEXT i%
        lsum = lsum + LOG(df)
    ELSE
        lsum = -2e30 * (its& - nd&)
    END IF
    IF k% = 1 AND its& > nd& THEN      'Don't print results too often
        l = .5 * lsum / (its& - nd&)
        LOCATE 2, 67: PRINT "I ="; its& - nd&;
        LOCATE 3, 67: PRINT USING "L =##.###"; l;
    END IF
END SUB

SUB sumle (x(), d%, n%, its&, nd&, sl)        'get sum of lypapunov exponents
        SHARED b!(), w!()
        STATIC lsum, k%
        IF its& = 0 THEN lsum = 0: k% = 0
        IF its& = nd& THEN lsum = 0            'discard iterates
        k% = (k% + 1) MOD 50
        df = 0
        FOR i% = 1 TO n%
            u = w!(i%, 0)
            FOR j% = 1 TO d%: u = u + w!(i%, j%) * x(j%): NEXT j%
            sechu = sech(u)        'hyperbolic secant of u
            df = df + b!(i%) * sechu * sechu * w!(i%, d%)
        NEXT i%
        df = abs(df)
        IF csng(df) > 0 THEN
            lsum = lsum + LOG(df)
        ELSE
            lsum = -1E+30 * (its& - nd&)
        END IF
        IF k% = 1 AND its& > nd& THEN    'don't print result too often
            sl = lsum / (its& - nd&)
            LOCATE 4, 67: PRINT USING "äL =###.###"; sl;
        END IF
END SUB

FUNCTION phi (u)        'Define squashing function (hyperbolic tangent)
    if abs(u) > 40 then
        phi = sgn(u)
    else
        phi = 2 / (1 + EXP(-2 * u)) - 1
    end if
END FUNCTION

FUNCTION sech (u)        'Hyperbolic secant of u
    if abs(u) > 40 then
        sech = 2 / EXP(abs(u))
    else
        sech = 2 / (EXP(u) + EXP(-u))
    end if
END FUNCTION

FUNCTION ran! STATIC    'Generate random number [0,1)
'$1000 if ran! isn't random over [0,1)
'Ref: W. A. Press and S. A. Teukolsky, Computers in Physics 6, 522 (1992)
'This routine is about 5X slower than the PowerBASIC RND function
'It reseeds itself with the hardware clock the first time it is called
IF idum& = 0 THEN      'Initialize on first call
    im1& = 2147483563
    im2& = 2147483399
    am! = 1 / im1&
    imm1& = im1& - 1
    ia1& = 40014
    ia2& = 40692
    iq1& = 53668
    iq2& = 52774
    ir1& = 12211
    ir2& = 3791
    ntab% = 32
    ndiv& = 1 + imm1& \ ntab%
    DIM iv&(ntab%)
    idum2& = 123456789
    FOR j% = 1 TO ntab%: iv&(j%) = 0: NEXT j%
    iy& = 0
    idum& = -100 * TIMER    'Reseed generator
    idum2& = idum&
    FOR j% = ntab% + 8 TO 1 STEP -1
        k& = idum& \ iq1&
        idum& = ia1& * (idum& - k& * iq1&) - k& * ir1&
        IF idum& < 0 THEN idum& = idum& + im1&
        IF j% <= ntab% THEN iv&(j%) = idum&
    NEXT j%
    iy& = iv&(1)
END IF
k& = idum& \ iq1&
idum& = ia1& * (idum& - k& * iq1&) - k& * ir1&
IF idum& < 0 THEN idum& = idum& + im1&
k& = idum2& \ iq2&
idum2& = ia2& * (idum2& - k& * iq2&) - k& * ir2&
IF idum2& < 0 THEN idum2& = idum2& + im2&
j% = 1 + iy& \ ndiv&
iy& = iv&(j%) - idum2&
iv&(j%) = idum&
IF iy& < 1 THEN iy& = iy& + imm1&
ran! = am! * iy&
END FUNCTION