program swarm c This is "Swarm explained" a fortran code to replicate 1995 results c showing the self organized formation of lines of ants trafic near the line c of a second order phase transition. This was the first demonstration of complex c trail structures emerging from the simplest (local, memoryless, homogeneus) model. c Figure 5 of Chialvo DR, Millonas MM. (1995) How swarms build cognitive maps. c In “The biology and technology of intelligent autonomous agents”. c Luc Steels (Ed.) NATO ASI Series, (144) 439-450. c Figure 2 in Rauch, Millonas and Chialvo, (1995) Pattern Formation and c Functionality in Swarm models, Physics Letters A 207, 185-193. c Comments to d-chialvo@northwestern.edu or dchialvo@gmail.com c-------------------------------------------------------------------------------- ! m: total number of posible positions ! mf: number of ants (In principle mf= 0.3 of m ) ! ig: array holding each ant position (between 1 and m) ! ip: array holding each ant orientacion, is an integer bewteen 1 y 8, ! for north ip=1; northeast ip=2; east, ip=3, southeast ip=4; etc until 8. ! g: Array holding the amount of pheromone at each position ! n: Accesory arrays for calculating first and second lattice neighbors ! u , w: Arrays ants position parameter (mm=32,l=mm,m=mm*mm) parameter (mf=341) dimension ig(1:mf) dimension ip(1:mf) dimension g(1:m) dimension n(1:8,1:m) dimension u(1:8) dimension w(1:8) c ---------------------------------Ants Physiology------------------------------ c ------- Antenna Sensory parameters (here is an array can be same for all) dimension beta(1:mf) dimension delta(1:mf) c---------Pheromone ! eta: Amount of pheromone each ant deposite in each step ! kappa: Pheromone evaporation rate eta=0.070d0 rkappa=0.98510d0 c--------Ants "Antenna" Directionality (Not crucial, only to make pretty pictures) c Opcion A): As in Chialvo and Millonas paper c factor0=1.0d0 front c factor1=1.0d0/2.0d0 right ad left front c factor2=1.0d0/4.0d0 right and left c factor3=1.0d0/12.0d0 right and left back c factor4=1.0d0/20.0d0 back c factor5=1.0d0/200.0d0 stay put c Opcion B): As in Rauch Chialvo and Millonas Phys Letter Paper -------- c factor0=1.0d0 c factor1=1.0d0 c factor2=1.0d0/8.0d0 c factor3=1.0d0/12.0d0 c factor4=1.0d0/50.0d0 c factor5=1.0d0/200.0d0 c------------------------------------------------------------------------------- llavor =12223 call dran_ini(llavor) c-------------------Save all ants final positions open(9,file='antspositions.txt',status='unknown') c Accesory arrays do 1 i=1,l do 1 j=1,l ij=(j-1)*l+i i1=i+1 if (i1.gt.l) i1=1 i2=i-1 if (i2.lt.1) i2=l j1=j+1 if (j1.gt.l) j1=1 j2=j-1 if (j2.lt.1) j2=l n(1,ij)=(j1-1)*l+i n(2,ij)=(j1-1)*l+i1 n(3,ij)=(j-1)*l+i1 n(4,ij)=(j2-1)*l+i1 n(5,ij)=(j2-1)*l+i n(6,ij)=(j2-1)*l+i2 n(7,ij)=(j-1)*l+i2 n(8,ij)=(j1-1)*l+i2 g(ij)=0.0 ! Set the pheromone field to zero to start 1 continue u(1)=0. u(2)=1./sqrt(2.) u(3)=1. u(4)=1./sqrt(2.) u(5)=0 u(6)=-1./sqrt(2.) u(7)=-1. u(8)=-1./sqrt(2.) w(1)=1. w(2)=1./sqrt(2.) w(3)=0. w(4)=-1./sqrt(2.) w(5)=-1. w(6)=-1./sqrt(2.) w(7)=0. w(8)=1./sqrt(2.) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c----------------------Main loops------------------------------ inte=1000 ! Total time to simulate for each beta,delta parameter ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc do 100 bbeta= 0.50d0,10.0d0,1.0d0 do 99 ddelta=0.0d0,0.50d0,0.10d0 write(*,*)bbeta,ddelta do 40 i=1,mf beta(i)=bbeta delta(i)=ddelta 40 continue c---Distribute all ants at random positionsand orientations do 30 i=1,mf ig(i)=m*dran_u()+1. ip(i)=8.*dran_u()+1. 30 continue do i=1,m g(i)=0.0 end do c Let the ants walk! do ite=1,inte call ritera(ite,inte,it,beta,delta,eta,rkappa, & ig,ip,g,n) enddo 99 continue 100 continue close(9) end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine ritera(ite,inte,it,beta,delta,eta,rkappa, & ig,ip,g,n) parameter (mm=32,l=mm,m=mm*mm) parameter (mf=341 ) dimension ig(1:mf) dimension ip(1:mf) dimension g(1:m) dimension n(1:8,1:m) dimension pro(1:8) dimension f(0:8) dimension kk(1:mf) dimension kp(1:mf) dimension beta(1:mf) dimension delta(1:mf) c Move the ants to the new position acording to the probability c dictated by the amount of pheromone in the neighbor sites, c plus the sensory parameters and c weighted by the antenna directionality c---------------Antenna directionality parameters f(0)=1.0 f(1)=1.0 !or 0.5 f(7)=f(1) f(2)=1./8. ! or 0.25 f(6)=f(2) f(3)=1./12. f(5)=f(3) f(4)=1./50. c-------Calculate probabilties based on beta and delta and pheromone field do 1000 k=1,mf prob=0.0d0 do 900 i=1,8 iss=n(i,ig(k)) pro(i)=1.+g(iss)/(1.0+delta(k)*g(iss)) pro(i)=pro(i)**beta(k) pro(i)=pro(i)*f(abs(ip(k)-i)) prob=prob+pro(i) 900 continue do 950 i=1,8 pro(i)=pro(i)/prob 950 continue c------------------Weight the probabilty with the antenna bias ur=dran_u() s2=pro(1)+pro(2) s3=s2+pro(3) s4=s3+pro(4) s5=s4+pro(5) s6=s5+pro(6) s7=s6+pro(7) s8=s7+pro(8) if(ur.le.pro(1)) then kk(k)=n(1,ig(k)) kp(k)=1 elseif(ur.le.s2) then kk(k)=n(2,ig(k)) kp(k)=2 elseif(ur.le.s3) then kk(k)=n(3,ig(k)) kp(k)=3 elseif(ur.le.s4) then kk(k)=n(4,ig(k)) kp(k)=4 elseif(ur.le.s5) then kk(k)=n(5,ig(k)) kp(k)=5 elseif(ur.le.s6) then kk(k)=n(6,ig(k)) kp(k)=6 elseif(ur.le.s7) then kk(k)=n(7,ig(k)) kp(k)=7 elseif(ur.le.s8) then kk(k)=n(8,ig(k)) kp(k)=8 endif 1000 continue do 1100 k=1,mf c Deposit the pheromone g(ig(k))=g(ig(k))+eta c Move the ants position and update directionality ig(k)=kk(k) ip(k)=kp(k) 1100 continue c Evaporate the pheromene do 1200 i=1,m if (g(i).gt.0.0) then g(i)=g(i)*rkappa else g(i)=0.0 endif 1200 continue c------------------------------------------Files----------------------------- if(ite.ge.inte-10)then ! in the last 10 iterations ! save (oversimposed)each ants positions do i=1,mf !1st translate m to square lattice jk=int(ig(i)/l)+1 ik=mod(ig(i),l) if (ik.eq.0)ik=l 13 format (f8.4,f8.4) fik=float(ik)/float(l)*0.95d0 fjk=float(jk)/float(l)*0.95d0 write(9,13)delta(i)+fik/10,beta(i)+fjk !Save in the format of Chialvo paper enddo end if c--------------------------------------------------------------------------- return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c///////////Random Numbers Subroutines from Imedea//////////////////////// subroutine dran_ini(iseed0) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) integer ix(ip) dimension g(0:m) data c0,c1,c2/2.515517,0.802853,0.010328/ data d1,d2,d3/1.432788,0.189269,0.001308/ c common /ixx/ ix common /icc/ ic common /gg/ g c dseed=iseed0 do 200 i=1,ip ix(i)=0 do 200 j=0,nbit-1 if(rand_XX(dseed).lt.0.5) ix(i)=ibset(ix(i),j) 200 continue ic=0 c pi=4.0d0*datan(1.0d0) do 1 i=m/2,m p=1.0-real(i+1)/(m+2) t=sqrt(-2.0*log(p)) x=t-(c0+t*(c1+c2*t))/(1.0+t*(d1+t*(d2+t*d3))) g(i)=x g(m-i)=-x 1 continue u2th=1.0-real(m+2)/m*sqrt(2.0/pi)*g(m)*exp(-g(m)*g(m)/2) u2th=nn1*sqrt(u2th) do 856 i=0,m 856 g(i)=g(i)/u2th return end subroutine dran_read(iunit) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(np=14) parameter(m=2**np) integer ix(ip) dimension g(0:m) common /ixx/ ix common /icc/ ic common /gg/ g read (iunit,*) ic read (iunit,*) (ix(i),i=1,ip) read (iunit,*) (g(i),i=0,m) return end subroutine dran_write(iunit) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(np=14) parameter(m=2**np) integer ix(ip) dimension g(0:m) common /ixx/ ix common /icc/ ic common /gg/ g write (iunit,*) ic write (iunit,*) (ix(i),i=1,ip) write (iunit,*) (g(i),i=0,m) return end function i_dran(n) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(iq=418) parameter(is=ip-iq) integer ix(ip) common /ixx/ ix common /icc/ ic ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif i_ran=ix(ic) if (n.gt.0) i_dran=mod(i_ran,n)+1 return end function dran_u() implicit double precision(a-h,o-z) parameter(ip=1279) parameter(iq=418) parameter(is=ip-iq) parameter (rmax=2147483647.0) integer ix(ip) common /ixx/ ix common /icc/ ic ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif dran_u=real(ix(ic))/rmax return end function dran_g() implicit double precision(a-h,o-z) parameter(ip=1279) parameter(iq=418) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) parameter(is=ip-iq) integer ix(ip) dimension g(0:m) common /ixx/ ix common /icc/ ic common /gg/ g ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif i=ishft(ix(ic),-np1) i2=iand(ix(ic),nn) dran_g=i2*g(i+1)+(nn1-i2)*g(i) return end subroutine dran_bm(x1,x2) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(iq=418) parameter(is=ip-iq) parameter (rmax=2147483647.0) integer ix(ip) common /ixx/ ix common /icc/ ic data pi2 /6.2831853072/ ic=ic+1 if (ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif u=pi2*real(ix(ic))/rmax ic=ic+1 if(ic.gt.ip) ic=1 if(ic.gt.iq) then ix(ic)=ieor(ix(ic),ix(ic-iq)) else ix(ic)=ieor(ix(ic),ix(ic+is)) endif v=(real(ix(ic))+0.5)/rmax v=sqrt(-2.0*log(v)) x1=cos(u)*v x2=sin(u)*v return end subroutine dran_gv(u,n) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(iq=418) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) parameter(is=ip-iq) dimension g(0:m) dimension u(n) dimension ix(ip) common /gg/ g common /ixx/ ix common /icc/ic n1=0 10 if (ic.lt.iq) then kmax=min(n-n1,iq-ic) do 99 k=1,kmax ic=ic+1 ix(ic)=ieor(ix(ic),ix(ic+is)) i=ishft(ix(ic),-np1) i2=iand(ix(ic),nn) u(n1+k)=i2*g(i+1)+(nn1-i2)*g(i) 99 continue else kmax=min(n-n1,ip-ic) do 98 k=1,kmax ic=ic+1 ix(ic)=ieor(ix(ic),ix(ic-iq)) i=ishft(ix(ic),-np1) i2=iand(ix(ic),nn) u(n1+k)=i2*g(i+1)+(nn1-i2)*g(i) 98 enddo endif if(ic.ge.ip) ic=0 n1=n1+kmax if (n1.lt.n) goto 10 return end subroutine dran_gvv(iv,u,n) implicit double precision(a-h,o-z) parameter(ip=1279) parameter(iq=418) parameter(np=14) parameter(nbit=31) parameter(m=2**np,np1=nbit-np,nn=2**np1-1,nn1=nn+1) parameter(is=ip-iq) dimension iv(n+ip) dimension g(0:m) dimension u(n) dimension ix(ip) common /gg/ g common /ixx/ ix common /icc/ic do 10 i=ic+1,ip 10 iv(i-ic)=ix(i) do 20 i=1,ic 20 iv(ip-ic+i)=ix(i) cCDEC$ INIT_DEP_FWD do 99 k=1,n ir=ieor(iv(k+is),iv(k)) i=ishft(ir,-np1) i2=iand(ir,nn) u(k)=i2*g(i+1)+(nn1-i2)*g(i) iv(k+ip)=ir 99 continue do 11 i=ic+1,ip 11 ix(i)=iv(n+i-ic) do 21 i=1,ic 21 ix(i)=iv(n+ip-ic+i) return end function rand_XX(dseed) double precision a,c,xm,rm,dseed,rand_XX parameter (xm=2.d0**32,rm=1.d0/xm,a=69069.d0,c=1.d0) dseed=mod(dseed*a+c,xm) rand_XX=dseed*rm return end