'Program CHAOSFIT.BAS fits a chaotic model to experimental data
'(c)1997 by J. C. Sprott (http://sprott.physics.wisc.edu/)
SCREEN 12
jmax% = 50  'Maximum time to correlate
ip% = 50    'Number of time points to plot
kmax% = 6   'Number of coefficients
dmax% = 2   'Number of variables (embedding dimension)
DIM cd(jmax%), cx(jmax%), y(dmax%), a(kmax%), am(kmax%)
OPEN "CHAOSFIT.DAT" FOR INPUT AS #1     'Get length of input data record
imax% = 0
WHILE NOT EOF(1)
    imax% = imax% + 1
    INPUT #1, d
WEND
CLOSE #1
imax% = imax% - 1
DIM d(imax%), x(imax%)
OPEN "CHAOSFIT.DAT" FOR INPUT AS #1     'Read input data
dav = 0
FOR i% = 0 TO imax%
    INPUT #1, d(i%)
    dav = dav + d(i%)
NEXT i%
dav = dav / (imax% + 1)     'Average value of input data
CLOSE #1
FOR j% = 0 TO jmax%         'Calculate autocorrelation of input data
    cd(j%) = 0
    FOR i% = 0 TO imax% - j%
        cd(j%) = cd(j%) + (d(i%) - dav) * (d(i% + j%) - dav)
    NEXT i%
    IF j% = 0 THEN cd0 = cd(0)
    cd(j%) = cd(j%) / cd0
NEXT j%
cd0 = SQR(cd0 / (imax% + 1))    'RMS value of input data
FOR i% = 0 TO imax%: d(i%) = (d(i%) - dav) / cd0: NEXT i%'Normalize to cd0
RANDOMIZE TIMER
rmin = 1E+38: eps = 4: rbest = rmin
WHILE INKEY$ <> CHR$(27)        'Loop until the <Esc> key is pressed
    FOR k% = 1 TO kmax%: a(k%) = am(k%) + eps * (RND - .5): NEXT k%
    FOR j% = 1 TO dmax%: y(j%) = 0: NEXT j% 'Zero initial conditions
    xav = 0
    FOR i% = 1 TO 1000 + imax%  'Generate model equation data
        IF ABS(y(1)) > 100 THEN EXIT FOR    'unbounded
        y(0) = a(1) + a(2) * y(1) + a(3) * y(1) * y(1) + a(4) * y(1) * y(2) + a(5) * y(2) + a(6) * y(2) * y(2)
        IF y(0) = y(1) THEN EXIT FOR        'fixed point
        FOR j% = dmax% TO 1 STEP -1: y(j%) = y(j% - 1): NEXT j%
            IF i% > 999 THEN x(i% - 1000) = y(1): xav = xav + y(1)
    NEXT i%
    IF i% = imax% + 1001 THEN   'Calculate autocorrelation of model data
        xav = xav / (imax% + 1)
        FOR i% = 0 TO imax%: x(i%) = x(i%) - xav: NEXT i%   'Subtract average
        FOR j% = 0 TO jmax%
            cx(j%) = 0
            FOR i% = 0 TO imax% - j%
                cx(j%) = cx(j%) + x(i%) * x(i% + j%)
            NEXT i%
            IF j% = 0 THEN cx0 = cx(0)
            cx(j%) = cx(j%) / cx0
        NEXT j%
        cx0 = SQR(cx0 / (imax% + 1))    'RMS value of model data
        FOR i% = 0 TO imax%: x(i%) = x(i%) / cx0: NEXT i%   'Normalize to cx0
        r = 0
        FOR j% = 1 TO jmax%     'Calculate sum of squared errors (r)
            rsqrt = cd(j%) - cx(j%)
            r = r + rsqrt * rsqrt
        NEXT j%
        IF r <= rmin AND r > 0 THEN     'Test quality of solution
            rmin = r
            IF rmin < rbest THEN        'A better solution was found
                BEEP
                rbest = rmin
                CLS                     'Graph it
                LINE (0, 100)-(639, 100), 2
                PSET (0, 100 - 40 * d(0))
                FOR i% = 1 TO ip%
                    LINE -(640 * i% / ip%, 100 - 40 * d(i%)), 15
                NEXT i%
                LINE (0, 240)-(639, 240), 2
                PSET (0, 240 - 40 * x(0))
                FOR i% = 1 TO ip%
                    LINE -(640 * i% / ip%, 240 - 40 * x(i%)), 12
                NEXT i%
                LINE (0, 400)-(639, 400), 2
                PSET (0, 300)
                FOR j% = 1 TO jmax%
                    LINE -(640 * j% / jmax%, 400 - 100 * cd(j%)), 15
                NEXT j%
                PSET (0, 300)
                FOR j% = 1 TO jmax%
                    LINE -(640 * j% / jmax%, 400 - 100 * cx(j%)), 12
                NEXT j%
                LOCATE 1, 1
                FOR k% = 1 TO kmax%: PRINT a(k%), ; : NEXT k%
                LINE (0, 0)-(639, 479), , B
            END IF
            FOR k% = 1 TO kmax%: am(k%) = a(k%): NEXT k%    'Move coefficients
            IF eps < 1 THEN eps = 4 * eps
            OPEN "CHAOSFIT.FIT" FOR OUTPUT AS #1            'Save model data
            FOR i% = 0 TO imax%: PRINT #1, x(i%): NEXT i%
            CLOSE #1
        ELSE
            eps = eps / 1.001   'Shrink the search neighborhood (gradually)
            IF eps < .000001 THEN   'Neighborhood too small, start over
                rmin = 1E+38
                eps = 4
                FOR k% = 1 TO kmax%: am(k%) = 0: NEXT k%
            END IF
        END IF
    END IF
WEND