PROGRAM genetic
DIM s$(1000),T(400,2),Eselect(1000)
CALL initial(s$(),L,L2,T(,),npop,nrecombine,nmutation,ngeneration)
FOR igeneration = 1 to ngeneration
    LET ntot = npop
    FOR iswap = 1 to nrecombine
        CALL recombine(s$(),L2,npop,ntot)
    NEXT iswap
    FOR i = 1 to nmutation
        CALL mutate(s$(),L2,npop,ntot)
    NEXT i
    CALL selection(s$(),L,L2,T(,),npop,ntot,Eselect())
NEXT igeneration
CALL showoutput(s$(),npop,Eselect())
END

SUB initial(s$(),L,L2,T(,),npop,nrecombine,nmutation,ngeneration)
    RANDOMIZE
    INPUT prompt "string (lattice) size = ": L
    LET L2 = L*L
    INPUT prompt "population number = ": npop
    ! L is linear dimension of phenotype
    INPUT prompt "number of recombinations per generation = ": nrecombine
    INPUT prompt "number of mutations per generation = ": nmutation
    INPUT prompt "number of generations = ": ngeneration
    ! create random population of genotypes
    FOR ipop = 1 to npop
        LET s$(ipop) = ""
        FOR i = 1 to L2
            IF rnd > 0.5 then
               LET s$(ipop) = s$(ipop) & "1"
            ELSE
               LET s$(ipop) = s$(ipop) & "0"
            END IF
        NEXT i
    NEXT ipop
    ! create random bond network of Tij's
    FOR i = 1 to L2
        FOR j = 1 to 2
            IF rnd > 0.5 then
               LET T(i,j) = 1
            ELSE
               LET T(i,j) = -1
            END IF
        NEXT j
    NEXT i
END SUB

SUB recombine(s$(),L2,npop,ntot)
    ! choose two strings (genotypes) to recombine
    LET i = int(npop*rnd) + 1
    DO
       LET j = int(npop*rnd) + 1
    LOOP until i <> j
    LET size = int(rnd*(L2/2)) + 1
    LET pos = int(rnd*L2) + 1
    LET s1$ = s$(i)
    LET s2$ = s$(j)
    IF pos + size <= L2 then
       LET s$(ntot+1) = s1$[1:pos-1]&s2$[pos:pos+size]&s1$[pos+size+1:L2]
       LET s$(ntot+2) = s2$[1:pos-1]&s1$[pos:pos+size]&s2$[pos+size+1:L2]
    ELSE                          ! apply periodic eboundarey conditions
       LET pbc = pos + size - L2
       LET s$(ntot+1) = s2$[1:pbc] & s1$[pbc+1:pos-1] & s2$[pos:L2]
       LET s$(ntot+2) = s1$[1:pbc] & s2$[pbc+1:pos-1] & s1$[pos:L2]
    END IF
    LET ntot = ntot + 2
END SUB

SUB mutate(s$(),L2,npop,ntot)
    LET i = int(rnd*npop) + 1
    LET pos = int(rnd*L2) + 1
    LET c$ = s$(i)[pos:pos]
    IF c$ = "1" then
       LET c$ = "0"
    ELSE
       LET c$ = "1"
    END IF
    LET s$(ntot + 1) = s$(i)[1:pos-1] & c$ & s$(i)[pos+1:L2]
    LET ntot = ntot + 1
END SUB

SUB convert(a$,L,s(,))
    ! coverts strings of 0's and 1's to 2D array of spins
    ! that is, genotype to phenotype
    FOR i = 1 to L
        FOR j = 1 to L
            LET n = (j-1)*L + i
            IF a$[n:n] = "1" then
               LET s(i,j) = 1
            ELSE
               LET s(i,j) = -1
            END IF
        NEXT j
    NEXT i
END SUB

SUB energy(L,s(,),T(,),E)
    LET E = 0
    FOR i = 1 to L
        LET ip = i + 1
        IF ip > L then LET ip = 1
        FOR j = 1 to L
            LET jp = j + 1
            IF jp > L then LET jp = 1
            LET n = (j-1)*L + i
            LET E = E - T(n,1)*s(i,j)*s(ip,j) - T(n,2)*s(i,j)*s(i,jp)
        NEXT j
    NEXT i
END SUB

SUB selection(s$(),L,L2,T(,),npop,ntot,Eselect())
    DIM s(30,30),save$(1000),Elist(0:1000)
    LET Esum = 0
    LET Elist(0) = 0
    FOR i = 1 to ntot
        CALL convert(s$(i),L,s(,))
        CALL energy(L,s(,),T(,),E)
        LET Esum = Esum -  E + 2*L2    ! contribution to Esum > 0
        LET Elist(i) = Esum
    NEXT i
    MAT save$ = s$
    ! select new population
    FOR ipop = 1 to npop
        LET E = Esum*rnd
        LET i = 0
        DO
           LET i = i + 1
        LOOP until E < Elist(i)   ! choose according to energy
        LET s$(ipop) = save$(i)
        LET Eselect(ipop) = Elist(i-1) - Elist(i) + 2*L2
    NEXT ipop
END SUB

SUB showoutput(s$(),npop,Eselect())
    FOR ipop = 1 to npop
        PRINT ipop,s$(ipop),Eselect(ipop)
    NEXT ipop
    PRINT
END SUB