##Initial Population size (N), Time span (T), detection probabilities (probs), and buffer widths (buffer, in m)

N <- 1000
T <- 20
I <- 500

probs <- c("None"=0.04,"Lo"=0.25,"Med"=0.5,"Hi"=0.75)
buffer <- c("Five"=1.52,"Ten"=3.05,"Fifteen"=4.6,"Twenty"=6.1,"T.Five"=7.6)

###Initial matrices and placeholders used in all models

#list into which growth is written (written over with each new parameter set):one slot for each individual
veg.pop <- matrix(0,nrow=N,ncol=T)

#vector for the number of years of growth (written over with each new parameter set): fed from "detection" function
realgrow <- numeric(N) 

#overall results array: contains the year at which each individual reach the buffer edge (or not, if 0) - raw output
veg.results <- array(NA,dim=c(I,T,length(probs),length(buffer)),dimnames=list("IT"=1:I,"Year"=1:T,"Detect"=names(probs),"Buff"=names(buffer)))

#array containing the number of escapes
yr.pool <- array(NA,dim=c(I,T,length(probs),length(buffer)),dimnames=list("IT"=1:I,"Year"=1:T,"Detect"=names(probs),"Buff"=names(buffer)))

#for the detection function (below): vector into which lifespan draws is written (written over with each new parameter set)
y <- numeric(N)

###Functions used in the model or data summary

#makes a vector of final grow years - one limit for each individual in the population(written into realgrow)
detection <- function(popsize,pD){
				
				y <- rgeom(popsize,pD)	
				y[y>(T-1)] = T-1
				return(y+1) }

#applied to veg.pop, identifies the year at which the plant reaches the buffer edge (if at all)
selex <- function(x) ifelse(any(x==1),min(which(x==1)),0)	



####Scenario 1: Clonal expansion only

#calculates cumulative plant growth (based on random draws from a normal distribution) to the time contained in realgrow (drawn from a geometric distribution)
###growth per year (in m) taken from Matlaga, Schutte, and Davis 2012
###Gaussian distribution: mean = 0.158m/yr; sd = 0.113 m/yr 

growth <- function(time) round(cumsum(rnorm(time,mean=0.158,sd=0.113)),3)	

###***The model

for(t in 1:I){	#for each iteration

for(i in seq(along=probs)){	#cycle through all detection levels
	for(j in seq(along=buffer)){	#combined with all possible buffer widths

veg.pop <- veg.pop*0 #resets the matrix of individual growth periods

realgrow <- detection(popsize=N,pD=probs[i])	#makes a vector of years different individuals are allowed to grow

for(k in seq(along=realgrow))	veg.pop[k,1:realgrow[k]] <- growth(time=realgrow[k])
					#calculates cumulative growth for the available years (drawm from realgrow)

#makes veg.pop binary: 1=at buffer edge, 0 not at edge or dead

veg.pop[veg.pop<buffer[j]]=0
veg.pop[veg.pop>0]=1

yr.pool[t,,i,j] <-colSums(veg.pop) #calculates the number present at the buffer edge each year

veg.results[t,,i,j] <- tabulate(apply(veg.pop,1,selex),nbins=T) #writes the number reaching the buffer each year

}##end buffer parameter
	}##end detection parameter (and single iteration) 
		} ##end all iterations

#########Clonal Results
###Checking results for zeroes: exclude Buffer=20 and Detect=High
any(yr.pool[,,,"Twenty"]>0)
any(yr.pool[,,,"Fifteen"]>0)
any(yr.pool[,,,"Ten"]>0)
any(yr.pool[,,,"Five"]>0)
any(yr.pool[,,,"T.Five"]>0)

any(yr.pool[,,"Hi","Twenty"]>0)

###Model Results
vegbuf<-ceiling(apply(veg.results,3:4,colMeans)) #####response variable: number of plants arriving at buffer edge each year
veg.ci.buf<-1.96*apply(veg.results,c(2:4),sd)
veg.pool <- ceiling(apply(yr.pool,3:4,colMeans))#####response variable: number of escapes over time
veg.ci <- 1.96*apply(yr.pool,c(2:4),sd) #####response variable: 95% confidence intervals around mean estimates

write.csv(vegbuf,"vegbuff.csv")
write.csv(veg.pool,"vegpop.csv")
write.csv(veg.ci,"vegci.csv")
write.csv(veg.ci.buf,"vegcibf.csv")
##############Scenario 2: Clonal Expansion plus Fragmentation

##Data for bankfull and scouring events are available in Appendix 2
##Import number of deposition events ("deposition") based on USGS data (frequency of bankfull events) and probability of deposition over a year
#(deposition events / 365 days in a year) 

s.events <-read.csv("scour.csv")
rownames(s.events) <- s.events[,1]
s.events <-s.events[,-1]
colnames(s.events) <- seq(1,20)

deposition <-read.csv("deps.csv")
rownames(deposition) <- deposition[,1]
deposition<-deposition[,-1]
colnames(deposition) <- seq(1,20)

#probability of floating, deposition, and establishment if you fragment 
establishment <- 0.02  #see article text (from Matlaga et al. 2012)
deposit <- 0.25	#see article text
float <- 0.7	#see Appendix 1

p.depfrag <- round(as.matrix(deposition/365)*float*deposit*establishment,4)

####Growth mirrors the clonal expansion model

growth <- function(time) round(cumsum(rnorm(time,mean=0.158,sd=0.113)),3)

####Accumulation of fragments from plants scoured at the edge: "vg" is "veg.pop" in each model run

fragmentation <- function(vg) ceiling(colMeans(t(t(scour)*colSums(vg))*prp.dep))

####Inital scramble for mixing up scouring events
iterat <- t(replicate(I,sample(1:T,T,replace=TRUE)))

####placeholder vector for number of fragments produced through time

frags <- numeric(T) 

#############The Model

for(t in 1:I){

scour <- as.matrix(s.events[,iterat[t,]])	#sets up the scouring events

prp.dep <- as.matrix(p.depfrag[,iterat[t,]])	#sets up the probability of a fragment successfully being deposited and established
		
####Veg only Model

for(i in seq(along=probs)){		#for all detection probabilities
	for(j in seq(along=buffer)){	#for all buffer widths

veg.pop <- veg.pop*0	#reset veg.pop to zero each time

realgrow <- detection(popsize=N,pD=probs[i])	#makes a vector of grow years lost to detection

for(k in seq(along=realgrow))	veg.pop[k,1:realgrow[k]] <- growth(time=realgrow[k])
					#calculates cumulative growth for the available years (minus years detected)

#makes veg.pop binary: 1=at buffer edge, 0 not at edge or dead

veg.pop[veg.pop<buffer[j]]=0
veg.pop[veg.pop>0]=1

frags <- fragmentation(veg.pop)	#calculates the number of fragments produced by plants at the edge (see function details above)

#determines what year, if any, plants reach the buffer edge or fragments establish, and writes it into the results matrix

veg.results[t,,i,j] <- tabulate(apply(veg.pop,1,selex),nbins=T) + frags  

#totals the number of plants each year plus the number of successful fragments

yr.pool[t,,i,j] <- colSums(veg.pop) + cumsum(frags) 

}##end buffer parameter
	}##end detection parameter (and single iteration) 
		} ##end all iterations


###Checking results for zeroes
any(yr.pool[,,,"T.Five"]>0)
any(yr.pool[,,,"Twenty"]>0)
any(yr.pool[,,,"Fifteen"]>0)
any(yr.pool[,,,"Ten"]>0)
any(yr.pool[,,,"Five"]>0)

any(yr.pool[,,"Hi","Twenty"]>0)


##distribution of years reached the buffer edge
fragbuf<-ceiling(apply(veg.results,c(2:4),mean))###Results: number of plants reaching the buffer edge each year
frag.ci.buf<-1.96*apply(veg.results,c(2:4),sd)
##number of escape over time
frag.pool <- ceiling(apply(yr.pool,c(2:4),mean))###Results:Accumulation of escapes over time
frag.ci <- 1.96*apply(yr.pool,c(2:4),sd)###95% confidence intervals around means of the different iterations

write.csv(fragbuf,"fragbuf.csv")
write.csv(frag.pool,"fragpop.csv")
write.csv(frag.ci,"fragci.csv")
write.csv(frag.ci.buf,"fragcibf.csv")

##############Scenario 3: Rhizome Movement plus Clonal Expansion

#####Prep code (run before model)

###Number of potential rhizome fragments left in the field after harvest

RF <- 500

###Year of rhizome harvest
H.yr <- 5

#rainfall events from 1992-2012 from 8 central IL weather stations (for frequency of year 1 movement events - see "rain" below)
rainfall <- read.csv("storm events 2.csv") 

#Long distance dispersal distance (in m) taken from the flume experiment (to populate "long", see below)
distances <- read.csv("distance.csv")$Dist.m

#probabilities and distance maximums (in m) for short distance rhizome movement (from the rainfall simulator experiment - see Appendix 2)
parms <- list(small=c("prob"=0.16,"hi"=0.1),large=c("prob"=0.08,"hi"=0.16))

#calculates plant growth through time, through the year provided from realgrow (fed from the "detection" function, see below)

growth.frag <- function(time,v,run){
rhiz.move[v,run] + round(cumsum(rnorm(time,mean=0.158,sd=0.113)),3)
#adding rhizome movement to the cumulative sum of growth over the remaining available years (from realgrow)
}

###Quantifying movement in Harvest year

##determines the number of rainfall events for each iteration (I) from "
rain <- sample(rainfall[,"Storms"],I,replace=TRUE) + 1

S.movt <- function(x,r) sum(runif(x,min=0.01,max=r)) #Short distance movement
L.movt<- function(r) rnbinom(RF,size=rain[r],prob=round(runif(1,min=0.995,max=0.999),3))#Long distance movement

#placeholder array - filled with movement parameters
blerg <-array(NA,dim=c(RF,I,2),dimnames=list("RF"=1:RF,"I"=1:I,"par"=c("max.move","move.prob")))

size <- 1*(replicate(I,rlnorm(RF,meanlog=0.21,sdlog=1.08)) < 1.1) +1 #generates rhizome sizes based on field collected samples (see Appendix 2)

blerg[,,"max.move"][which(size==1)] = parms[[1]]["hi"]
blerg[,,"max.move"][which(size==2)] = parms[[2]]["hi"]

blerg[,,"move.prob"][which(size==1)] = parms[[1]]["prob"]
blerg[,,"move.prob"][which(size==2)] = parms[[2]]["prob"]

#calculating short distance movement
num.move <- ceiling(t(t(blerg[,,"move.prob"])*rain))

move <-apply(num.move,c(1:2),S.movt,r=blerg[,,"max.move"])

#calculating long distance movement
long <- sapply(1:length(rain),L.movt)

long[which(long>0)] <- as.numeric(lapply(sapply(long[which(long>0)],sample,replace=TRUE,x=distances),sum))

#calculating the number of long distance movements that carried rhizomes over the buffer edge.  These individuals
#were considered outside the monitoring zone, and their survival mirrored the "No detection" plants

esc.p <- array(NA,c(RF,I,length(buffer)),dimnames=list("RF"=1:RF,"I"=1:I,"Buff"=names(buffer)))
esc.p[1:RF,1:I,]=long

for(r in seq(along=buffer)) esc.p[,,r] <- 1*(esc.p[,,r]>buffer[r])
esc <- apply(esc.p,3,colSums)

#adding short and long distance movement
rhiz.move <- move + long

#proportion of year one escapes surviving to T

prp.surv<-function(p)(1-(p*0.2))^(T-H.yr)

#list into which growth is written (written over with each new parameter set):one slot for each individual
frag.pop <- matrix(0,nrow=N,ncol=T)

#vector for the number of years of growth (written over with each new parameter set): fed from "detection" function
frag.grow <- numeric(RF) 


#######End required pre-run

#####*****The Model

for(t in 1:I){

for(i in seq(along=probs)){
	for(j in seq(along=buffer)){
		
veg.pop <- veg.pop*0

realgrow <- detection(popsize=N,pD=probs[i])	#makes a vector of grow years lost to detection
###**those rhizomes carried over the buffer edge in year 1 retain the same survivorship as the "No detection" plants, regardless of detection level
###so not included

for(k in seq(along=realgrow))	veg.pop[k,1:realgrow[k]] <- growth(time=realgrow[k])

#calculates cumulative growth for the available years (lifespan) plus movement in year 1

#determines what year, if any, the plant reaches the buffer edge, and write it into the results matrix

#makes veg.pop binary: 1=at buffer edge, 0 not at edge or dead

veg.pop[veg.pop<buffer[j]]=0
veg.pop[veg.pop>0]=1


####Fragments added and moved in Harvest year (H.yr above)

frag.grow <- detection(popsize=RF-esc[t,j],pD=(probs[i]*0.9))	#makes a vector of grow years lost to detection - excluding those that move past the buffer
frag.grow[frag.grow > (T-H.yr)] <- T-H.yr #years of growth available for fragments
frag.end <-frag.grow+H.yr	#inserting growth over the correct time span (see below growth equation)
###**those rhizomes carried over the buffer edge in year 1 retain the same survivorship as the "No detection" plants, regardless of detection level
###so not included

for(k in seq(along=frag.grow))	frag.pop[k,H.yr:frag.end[k]] <- c(rhiz.move[k,t],growth.frag(time=frag.grow[k],v=k,run=t))

#calculates cumulative growth for the available years (lifespan) plus movement in year 1

#determines what year, if any, the plant reaches the buffer edge, and write it into the results matrix

#makes veg.pop binary: 1=at buffer edge, 0 not at edge or dead

frag.pop[frag.pop<buffer[j]]=0
frag.pop[frag.pop>0]=1

##combining original and harvest populations

yr.pool[t,,i,j] <-colSums(veg.pop) + colSums(frag.pop) 
yr.pool[t,H.yr:T,i,j]<-yr.pool[t,H.yr:T,i,j] + ceiling(prp.surv(probs[i])*esc[t,j])

veg.results[t,,i,j] <- tabulate(apply(veg.pop,1,selex),nbins=T) + tabulate(apply(frag.pop,1,selex),nbins=T)

veg.results[t,H.yr:T,i,j] <- veg.results[t,H.yr:T,i,j] + esc[t,j]

}##end buffer parameter
	}##end detection parameter (and single iteration) 
		} ##end all iterations

###Results

###Checking results for zeroes
any(yr.pool[,,,"T.Five"]>0)
any(yr.pool[,,,"Twenty"]>0)
any(yr.pool[,,,"Fifteen"]>0)
any(yr.pool[,,,"Ten"]>0)
any(yr.pool[,,,"Five"]>0)

any(yr.pool[,,"Hi","Twenty"]>0)

mvtbuf<-ceiling(apply(veg.results,3:4,colMeans))##distribution of plants reaching the edge each year
mvt.ci.buf<-1.96*apply(veg.results,c(2:4),sd)
mvt.pool <- ceiling(apply(yr.pool,c(2:4),mean))##accumulation of plant escapes over time
mvt.ci <- 1.96*apply(yr.pool,c(2:4),sd)## 95% confidence intervals around mean estimates of escapes

write.csv(mvtbuf,"mvtbuf.csv")
write.csv(mvt.pool,"mvtpop.csv")
write.csv(mvt.ci,"mvtci.csv")
write.csv(mvt.ci.buf,"mvtcibf.csv")