!
! Neural network for pattern recognition using Ising spins on square lattice.
! Version with some bells. No cheating.
!   H. Nakanishi
!
program net
   option nolet
   randomize
   dim spinp(1,1,1),spin(1,1),jij(1,1)
   call initialize(nmax,imax,spinp(,,),spin(,),jij(,),nattempts,t,dt,pdamage,step)
   mat redim spinp(0 to imax,nmax,nmax)
   mat redim spin(nmax,nmax),jij(nmax^2,nmax^2)
   width=min(1/(imax+1),0.5)
   clear
   for i=1 to imax
     open #i: screen i*width,(i+1)*width,0,width
   next i
   open #(imax+1): screen 0,width,0,width
   open #(imax+2): screen 0.41,0.9,0.51,1
   open #(imax+3): screen 0.1,0.45,0.7,0.98
   for i=1 to imax
     window #i
     call setpattern(nmax,imax,i,spinp(,,),width)
   next i
   window #(imax+1)
   call setpattern(nmax,imax,0,spinp(,,),width)
   call getmin(kmin,nmax,imax,spinp(,,),sum0)
   call init_energy(nmax,imax,spinp(,,),spin(,),jij(,),energy,pdamage,damage)
   energy0=energy
   window #(imax+2)
   call calculate(nmax,imax,nattempts,nmoves,ntotal,spin(,),jij(,),energy,t,dt,step)
   window #(imax+3)
   set background color "white"
   clear
   set color "black"
   print "energy=";energy
   print "energy0=";energy0
   print "pdamage=";pdamage
   print "damaged bonds=";damage
   print "init hamd=";sum0;"(";kmin;")"
   print "temp=";t
   print "dt=";dt
   print "tot sweeps=";ntotal
   print "tot accepted=";nmoves
   get key z
   close #(imax+2)
   close #(imax+3)
   for i=1 to imax+1
     close #i
   next i
end

! initialize variables
! nmax = side of lattice
sub initialize(nmax,imax,spinp(,,),spin(,),jij(,),nattempts,t,dt,pdamage,step)
   input prompt "lattice size nmax => ": nmax ! lattice size is nmax by nmax
   input prompt "number of patterns in memory => ": imax
   input prompt "zero temp [1] or finite temp [2]? => ": anneal
   if anneal=1 then
      t=0
   else
      input prompt "temperature (in units of J/k) => ": t
      input prompt "temperature increment unit => ": dt
   end if
   input prompt "damage probability => ": pdamage
   input prompt "plot after this many lattice sweeps => ": nattempts
   input prompt "step through manually? Yes [1] or No [2] => ": step
end sub

sub setpattern(nmax,imax,i,spinp(,,),width)
   set window 0,nmax+2,0,nmax+2
   set background color "white"
   clear
   set color "black"
   plot 1,1;1,nmax+1;nmax+1,nmax+1;nmax+1,1;1,1
   set color mix(45) 2/3, 2/3, 2/3
   for j=1 to nmax
     for k=1 to nmax
       spinp(i,j,k)=-1
     next k
   next j
   do
     get point x,y
     if x<1 or x>nmax+1 or y<1 or y>nmax+1 then exit do
     x1=int(x)
     y1=int(y)
     spinp(i,x1,y1)=-spinp(i,x1,y1)
     if spinp(i,x1,y1)=1 then
       if i=0 then
         set color 45
       else
         set color "black"
       end if
     else
       set color "white"
     end if
     box area x1,x1+1,y1,y1+1
   loop
end sub

sub calculate(nmax,imax,nattempts,nmoves,ntotal,spin(,),jij(,),energy,t,dt,step)
   dim spin0(1,1)
   mat redim spin0(nmax,nmax)
   set window 0,nmax+2,0,nmax+2
   set background color "white"
   clear
   set color "black"
   box lines 1,nmax+1,1,nmax+1
   ntotal=0
   nmoves=0
   do 
     for i=1 to nmax
       for j=1 to nmax
         spin0(i,j)=spin(i,j)
       next j
     next i
     na = 0
     for i=1 to nattempts
     for nx=1 to nmax
     for ny=1 to nmax
       oldspin=spin(nx,ny)
       spin(nx,ny)=-spin(nx,ny)
       oldenergy=energy
       call update_energy(nmax,spin(,),jij(,),energy,nx,ny)
       if energy > oldenergy then
         if t=0 or exp((oldenergy-energy)/t)<rnd then
            spin(nx,ny)=oldspin
            energy=oldenergy
         else
           na=na+1
         end if
       else
         na=na+1
       end if
     next ny
     next nx
     next i
     ntotal = ntotal + nattempts
     nmoves = nmoves + na
     set color "white"
     call display(spin0(,),nmax)
     set color "red"
     call display(spin(,),nmax)
     if key input or step=1 then
       get key z
       c$=chr$(z)
       if c$="q" then
         exit do
       else if c$="p" then
         get key z
       else if c$="d" then
         dt=dt*2
       else if c$="h" then
         dt=dt*0.5
       else if c$="+" then
         t=t+dt
       else if c$="-" then
         t=t-dt
       end if
     end if
     if t=0 and na=0 then exit do
   loop
end sub

sub display(spin(,),nmax)
   for i=1 to nmax
     for j=1 to nmax
       if spin(i,j)=1 then
         box area i,i+1,j,j+1
       end if
     next j
   next i
end sub

sub init_energy(nmax,imax,spinp(,,),spin(,),jij(,),energy,pdamage,damage)
   for i=1 to nmax^2
     yi=int((i-1)/nmax)+1
     xi=i-(yi-1)*nmax
     for j=i to nmax^2
      jij(i,j)=0
      if i<>j then
       yj=int((j-1)/nmax)+1
       xj=j-(yj-1)*nmax
       for k=1 to imax
         jij(i,j)=jij(i,j)+spinp(k,xi,yi)*spinp(k,xj,yj)
       next k
      end if
     next j
   next i
   for j=1 to nmax^2
     for i=j+1 to nmax^2
       jij(i,j)=jij(j,i)
     next i
   next j
   damage=0
   if pdamage<>0 then
     for i=1 to nmax^2
       for j=1 to nmax^2
         if rnd<pdamage then
           jij(i,j)=0
           damage=damage+1
         end if
       next j
     next i
   end if
   energy=0
   for i=1 to nmax
     for j=1 to nmax
       spin(i,j)=spinp(0,i,j)
     next j
   next i
   for i=1 to nmax^2
     for j=i+1 to nmax^2
       yi=int((i-1)/nmax)+1
       xi=i-(yi-1)*nmax 
       yj=int((j-1)/nmax)+1
       xj=j-(yj-1)*nmax
       energy=energy-jij(i,j)*spin(xi,yi)*spin(xj,yj)
     next j
   next i
end sub

sub update_energy(nmax,spin(,),jij(,),energy,nx,ny)
   i=(ny-1)*nmax+nx
   for j=1 to nmax^2
     yj=int((j-1)/nmax)+1
     xj=j-(yj-1)*nmax
     energy=energy-2*jij(i,j)*spin(nx,ny)*spin(xj,yj)
   next j
end sub

sub getmin(k,nmax,imax,spinp(,,),sum0)
   k=0
   sum0=1e9
   for i=1 to imax
     sum=0
     for l=1 to nmax
       for m=1 to nmax
         sum=sum+(spinp(0,l,m)-spinp(i,l,m))^2
       next m
     next l
     if sum<sum0 then
       k=i
       sum0=sum
     end if
   next i
end sub