# Spatial Invasion Simulator v2
#      plant feedback lattice model with two perennial plants
#	
# W Stan Harpole, 2007
#
# parameters
	# Species
	sp1.lh  <- 1 		# life history: 0=annual, 1=perennial
	sp1.dis <- 1 		# dispersal: 0=local, 1= global
	sp1.fb  <- 9 		# feedback: 0-9= none-strong
	sp1.fb<- -sp1.fb 	# species 1 values less than zero
	sp2.lh  <- 1 		# 0=annual, 1=perennial
	sp2.dis <- 0 		# 0=local, 1= global
	sp2.fb  <- 0 		# 0-9= none-strong
		
	# Scenario: invasion: few invasives into mostly natives
	#			restoration: establish natives into mostly exotics
	scenario <- 0 	# 0=invasion, 1=restoration, 2=equal
	
	# Management
	mgt.herb <- 0   # 0= no, 1=  herbicide
	mgt.seed <- 0   # 0= no, 1= seed natives
	mgt.fb   <- 0   # 0= no, 1= reset feedback state to 0
	mgt.crit <- 0   # criteria: 0=patch, 1=relative landscape abundance
	mgt.lev  <- 0 	# level: 0-1= proportion of cells or abundance
	mgt.resp <- 0   # response: 0=none, 1=patch, small scale; 2=diffuse, broad scale 
	mgt.int  <- 0 	# intensity: 0-1 proportion of lattice that can be managed each time step
	
	# Model
	lat   <- 20 # nxn lattice dimension 10-100
	lat2  <- lat^2 # number of cells
	steps <- 100 # time steps to run 1-1000
	#speed <- # how fast to run the simulation
	edge <- 0 # boundary type: 0=absorbing, 1=wrapping

# initialize
	# species abundance
	x1<-c()
	x2<-c()
	x0<-c()

	# fill nxn grid with two species (1 or 2) randomly
	if(scenario==0){ # invasion, 5% species 1
		x<-matrix(rep(2,lat2),nrow=lat,ncol=lat)
		sp1.0<-sample(1:lat2,0.05*lat2)
		x[sp1.0]<-1
		state<-matrix(sp2.fb, nrow=lat, ncol=lat) # feedback state
	} else if(scenario==1){ #restoration, 5% species 2
		x<-matrix(sample(rep(1,lat2)),nrow=lat,ncol=lat)
		sp2.0<-sample(2:lat2,0.05*lat2)
		x[sp2.0]<-2		
		state<-matrix(sp1.fb, nrow=lat, ncol=lat)
	} else {	# equal initial abundance
		x<-matrix(sample(rep(1:2,lat2/2)),nrow=lat,ncol=lat)
		state<-matrix(0, nrow=lat, ncol=lat)
	}

	# create 5% empty sites
	empty<-sample(1:lat2,0.05*lat2)
	x[empty]<-0
	# maximum number of managment weed control applications per time step/"season" where each = 9 cells
	mgt.max<-round(mgt.int*lat2/9)

palette(c("white","orange","blue"))
# plot initial spatial config
par(mfrow=c(1,2), pty="s", mar=c(2,2,0,1) )
symbols( row(x), col(x), squares=rep(1,lat2), inches=F, fg="gray", bg=x+1 )



for(k in 1:steps){		 # k time steps 
	
	# random order of updates each time step
	y<-sample(1:lat2)
	
	# management effort (number of weed control applications)
	mgt.ct<-0 # initialize each time step/"season"

	
		for(i in 1:lat2){	   # single turnover of population, asynchronous updating
		    yi<-y[i]     # update random cell 
		    xi<-x[yi]    # occupant of random cell
		   	
		    
		    # Boundaries
			if(edge==0){			# absorbing edges
		    	temp<-cbind(rep(0,lat), x[,1:lat], rep(0,lat) )
		    	wrap<-rbind(rep(0,lat+2), temp[1:lat,], rep(0,lat+2))
		    } else {				 # wrapping edges
		    	temp<-cbind(x[,lat], x[,1:lat], x[,1])
		    	wrap<-rbind(temp[lat,], temp[1:lat,], temp[1,])
			}

 		 	# nearest neighbors
		    m<- row(x)[yi]
		    n<- col(x)[yi]
		    nearest<- c( wrap[m,n], wrap[m,n+1], wrap[m,n+2], wrap[m+1,n], wrap[m+1,n+1],
		           		wrap[m+1,n+2], wrap[m+2,n], wrap[m+2,n+1], wrap[m+2,n+2] ) # include current cell
		    if( length(nearest[nearest>0])==0 ){ next } # skip if not at least one neighboring plant		  
		  	# management/ weed control application
		  	if(mgt.resp>0 & mgt.ct<mgt.max){
		  			# herbicide patch if all spp1 but allow previous residents to still contribute seed for colonization
		  			if( length(nearest[nearest==1])/9 >= mgt.lev ){ 
		  					if(mgt.herb==1){ wrap[ c(m,m+1,m+2), c(n,n+1,n+2) ] <-0 }  # herbicide the patch
		   					if(mgt.seed==1){ wrap[ c(m,m+1,m+2), c(n,n+1,n+2) ] <-2 }  # herbicide and reseed patch w/ spp2
		   						# reseeding alone not supported yet
		   					x<-wrap[2:(lat+1),2:(lat+1)] # update x
		   					mgt.ct=mgt.ct+1 # increment managment count
		   			}
		   	}
		   	
		   	# feedback effect on survival and establishment
		   	s1<- s2<-1
		    if(state[yi]>0){ s1<- (1-abs(state[yi])/10)   
		    } else { 
		    	s2<- (1-abs(state[yi])/10)  
		    }
			
			 # empty cell, colonization, proportion of local or global neighborhood
		    if(xi==0 ){  
		    		if(sp1.dis==0) { #local
		    			c1<- ( length(nearest[nearest==1])/9 )
		    		} else { #global
		    			c1<- ( length(x[x==1])/lat2 )
		    		}
		    		if(sp2.dis==0) { 
		           		c2<- ( length(nearest[nearest==2])/length(nearest[nearest>0]) ) 
		    		} else { 
		    			c2<- ( length(x[x==2])/length(x[x>0]) )
		    		}
		    		# standardize probabilities of colonization weighted by feedback effect
		    		p1<- c1*s1/(c1*s1+c2*s2)
		    		p2<- 1-p1
		    		if(runif(1) < p1){ x[yi]<-1 } else { x[yi]<-2 }
 
		  	}
		  	
		  	# species 1 occupies cell, 5% mortality + additional feedback effect
		  	if(xi==1){     
		           if(runif(1) < (0.05 + (1-s1)) ){ x[yi]<- 0 }
		           if(state[yi] > sp1.fb){ state[yi]<- state[yi] - 1 } # decrement cell state in favor of species 1
		  	} 
		  	
		  	#species 2
		  	if(xi==2){  
		  		   if(runif(1) < (0.05 + (1-s2)) ){ x[yi]<- 0 }  # survive?
		           if(state[yi] < sp2.fb){state[yi]<- state[yi] + 1} # increment cell state in favor of species 2
		  	}
		}

x1[k]<-length(x[x==1])
x2[k]<-length(x[x==2])
x0[k]<-length(x[x==0]) #empty cells

# spatial plot
symbols( row(x), col(x), squares=rep(1,lat2), inches=F, fg="gray", bg=x+1, add=T)
#hist(state,breaks=c(-5.5,-4.5,-3.5,-2.5,-1.5,-.5,0.5,1.5,2.5,3.5,4.5,5.5), ylim=c(0,lat2) )
}



# abundance vs time
plot(1:steps,1:steps*(lat2/steps),type="n")
lines(1:steps,x1,col=2)
lines(1:steps,x2,col=3)

palette("default")