# Epidemics
#
# Probability of catching the disease
#   exp(alfa+beta*neighbours)/1+exp(alfa+beta*neighbours)
#
# Probability of being cured
#   exp(delta+gama*neighbours)/1+exp(delta+gama*neighbours)
map_matrix(0,ncol=30,nrow=50)
ney_map
# Parameters
alfa_10.
beta_-5.5
delta_-10.
gama_5.9
# Indices, odd & even
odx_seq(1,50,2)
ody_seq(1,30,2)
evex_seq(2,50,2)
evey_seq(2,30,2)
# Initialisation
map[1,20:21]_1
map[2,19:22]_1
map[3,20:21]_1
map[4,21]_1
image(map,col=heat.colors(2))
delay(3)
# Propagation by Swendson-Wang
for (t in 1:10000){
# Neighbours, inside then borders, then corners
  ney[2:49,2:29]_map[1:48,2:29]+map[3:50,2:29]+map[2:49,1:28]+map[2:49,3:30]
  ney[1,2:29]_map[2,2:29]+map[1,1:28]+map[1,3:30]
  ney[50,2:29]_map[49,2:29]+map[50,1:28]+map[50,3:30]
  ney[2:49,1]_map[2:49,2]+map[1:48,1]+map[3:50,1]
  ney[2:49,30]_map[2:49,29]+map[1:48,30]+map[3:50,30]
  ney[1,1]_map[1,2]+map[2,1]
  ney[1,30]_map[1,29]+map[2,30]
  ney[50,1]_map[50,2]+map[49,1]
  ney[50,30]_map[50,29]+map[49,30]
# Moves
  alea_matrix(runif(50*30),ncol=30)
  map[odx,ody]_map[odx,ody]*(alea[odx,ody]<(1/(1+exp(delta+gama*ney[odx,ody]))))+
     (1-map[odx,ody])*(alea[odx,ody]<(1/(1+exp(alfa+beta*ney[odx,ody]))))*(ney[odx,ody]>0)
  map[evex,evey]_map[evex,evey]*(alea[evex,evey]<(1/(1+exp(delta+gama*ney[evex,evey]))))+
     (1-map[evex,evey])*(alea[evex,evey]<(1/(1+exp(alfa+beta*ney[evex,evey]))))*(ney[evex,evey]>0)
  map[evex,ody]_map[evex,ody]*(alea[evex,ody]<(1/(1+exp(delta+gama*ney[evex,ody]))))+
     (1-map[evex,ody])*(alea[evex,ody]<(1/(1+exp(alfa+beta*ney[evex,ody]))))*(ney[evex,ody]>0)
  map[odx,evey]_map[odx,evey]*(alea[odx,evey]<(1/(1+exp(delta+gama*ney[odx,evey]))))+
     (1-map[odx,evey])*(alea[odx,evey]<(1/(1+exp(alfa+beta*ney[odx,evey]))))*(ney[odx,evey]>0)
# Plot
  image(map,col=heat.colors(2))
  delay(3)
}