#Example 3 - Spatial host parasitoid model

Z<-30
H<- matrix (NA, nrow=Z, ncol=Z) 
P<- matrix (NA, nrow=Z, ncol=Z) 

G<- 100 

Beddington<-function(fNt,fPt,fr,falpha,fK,fc){
NtPlus1<-fNt*exp((fr*(1-(fNt/fK)))-(falpha*fPt))
PtPlus1<-fc*fNt*(1-(exp(-falpha*fPt)))
answer<-c(NtPlus1,PtPlus1)
answer
}

random_asign_to_neighbours<-function(to_Add_f, Number_moving, neighbours){
cell<-rep(NA,8)
breaks<- sort(sample(1: Number_moving, 7, replace=TRUE))
for(i in 2:7) {cell[i] <-breaks[i]- breaks[i-1]}
cell[1]<-breaks[1]
cell[8]<- Number_moving - breaks[7]
for(x in 1:8){
to_Add_f[neighbours[x,1],neighbours [x,2]]<- to_Add_f[neighbours[x,1],neighbours [x,2]] + cell[x]
} 
to_Add_f
} # random_asign_to_neighbours function




H[1:Z,1:Z]<-100
P[1:Z,1:Z]<-10
alpha<-0.0025
r<-2.1
c<-1
K<-1000 
proportion_to_migrate<-0.1
neighbours<-matrix(NA, nrow=8, ncol=2)
par(mfrow=c(3,2))


save_runs<-c(10,30,60,100) 


hold_H<-matrix(NA, nrow=Z*Z, ncol=length(save_runs))# 4 runs
hold_P<- matrix(NA, nrow=Z*Z, ncol=length(save_runs))
count<-0

for(g in 1:G) { 
	for(ba in 1:Z){ 
		for(bb in 1:Z){
			answer<-Beddington(H[ba,bb],P[ba,bb],r,alpha,K,c)
			H[ba,bb]<-round(answer[1])
			P[ba,bb]<-round(answer[2])
		} # for bb loop
	} # for ba loop


tomove_Host<-round(H* proportion_to_migrate)
tomove_Parasitoid<-round(P* proportion_to_migrate)
H<- round(H-tomove_Host) 
P<- round(P-tomove_Parasitoid) 
to_add_Host<-matrix(0, nrow=Z, ncol=Z)
to_add_Parasitoid<-matrix(0, nrow=Z, ncol=Z)

for(ba in 1:Z){
	for(bb in 1:Z){  
		neighbours[1,1:2]<- c(ba-1, bb-1) 
		neighbours[2,1:2]<- c(ba-1, bb)
		neighbours[3,1:2]<- c(ba-1, bb+1)
		neighbours[4,1:2]<- c(ba, bb-1)
		neighbours[5,1:2]<- c(ba, bb+1)
		neighbours[6,1:2]<- c(ba+1,bb-1)
		neighbours[7,1:2]<- c(ba+1,bb)
		neighbours[8,1:2]<- c(ba+1,bb+1)
		neighbours[which(neighbours==0)]<- Z 
		neighbours[which(neighbours==Z+1)]<-1 
to_add_Host <- random_asign_to_neighbours (to_add_Host, H[ba,bb], neighbours)
to_add_Parasitoid <- random_asign_to_neighbours (to_add_Parasitoid, P[ba,bb], neighbours) 
	} #end bb loop
} # end ba loop


H<-H+ to_add_Host
P<- P+ to_add_Parasitoid

if(g %in% save_runs){
count<-count+1
hold_H[,count]<- unlist(expand.grid(H))
hold_P[,count]<- unlist(expand.grid(P))
} #end if g in save list
}   # for g loop




plot_expandgrid<-function(X,Z){ 
mcol<-max(X)
start<-0
for(j in 1:Z){
	for(k in 1:Z){
	     start<-start+1
             points(j,k, col=round(10*X[start]/mcol),pch=15)
             } # for k loop
        } # for j loop
} # end plot_generation function


getcoloursrainbow<-function(x){
purplef<-220
f<-unlist(as.list(rainbow(256))) 
mn<-0
mx<-max(x,na.rm=TRUE)
cor<-(x-mn)/(mx-mn)
ans<-f[(1+ceiling(purplef*(cor)))]
ans
}


par(mfrow=c(2,2))

point_colours<- getcoloursrainbow(hold_H[,1])
point_colours2<-matrix(point_colours,nrow=Z)


plot(NULL, xlim=c(1,Z),ylim=c(1,Z), axes=F, xlab="",ylab="", main="Host")
for(j in 1:Z){
	for(k in 1:Z){
points(j,k, col= point_colours2[j,k],pch=15)
}}



point_colours3<- getcoloursrainbow(hold_P[,1])
point_colours4<-matrix(point_colours3,nrow=Z)
plot(NULL, xlim=c(1,Z), ylim=c(1,Z), axes=F, xlab="" ,ylab="", main="Parasitoid")
for(j in 1:Z){
	for(k in 1:Z){
points(j,k, col= point_colours4[j,k],pch=15)
}}


point_colours<- getcoloursrainbow(hold_H[,4])
point_colours2<-matrix(point_colours,nrow=Z)
plot(NULL, xlim=c(1,Z),ylim=c(1,Z), axes=F, xlab="",ylab="", main="Host")
for(j in 1:Z){
	for(k in 1:Z){
points(j,k, col= point_colours2[j,k],pch=15)
}}


point_colours3<- getcoloursrainbow(hold_P[,4])
point_colours4<-matrix(point_colours3,nrow=Z)
plot(NULL, xlim=c(1,Z), ylim=c(1,Z), axes=F, xlab="" ,ylab="", main="Parasitoid")
for(j in 1:Z){
	for(k in 1:Z){
points(j,k, col= point_colours4[j,k],pch=15)
}}



centreX<- -7
centreY<- 45
LX<-8
lowY<- centreY-2
highY<- centreY+2
increment<-0.08
line_colour<- getcoloursrainbow(1:220)
for(i in 1:220){
X<- centreX-LX+(i*increment)
lines (c(X,X),c(lowY,highY), col=line_colour[i], xpd=NA)
}
text(xpd=NA, -24, centreY, "Population = 0")
text(xpd=NA, 15, centreY, "Maximum population")

text(centreX, centreY+6,paste("After",save_runs[1],"generations", sep=" "),xpd=NA)
text(centreX, -5,paste("After",save_runs[4],"generations", sep=" "),xpd=NA)