# Functions.R ---> defines reproduction, dispersal, and climate velocity functions used in simulations
#Laplace dispersal kernel
k<-function(x,y,b) return(1/2*b*exp(-b*abs(x-y)))
# dispersal matrix
d<-array(,c(w,w))
#dispersal probabilities to point i from every point j
for(i in 1:w) d[i,]=k(world[i],world,b)
#Beverton-Holt recruitment
f<-function(n,R0,K) return(R0*n/(1+(R0-1)/K*n))
#Standard error removing NAs
stderr <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))
moveMPA <- function(MPA.current = MPA.current, displaced = displaced, mpa.yes=mpa.yes, mpa.no=mpa.no, world){
##### Move MPAS
# move MPA forward by amount
next_MPA = MPA.current[displaced:length(MPA.current)]
lost <- MPA.current[1:(displaced-1)]
# are there any MPAs in ?
#any(test[[i]]$next_MPA==1)
# IF FALSE, then need to figure out how many 0s were lost in
# move, and make sure to preserve interval of zeros == mpa.no,
# then fill in mpa.yes,mpa.no to length of world. But also this
# is only for when mpa.no exists on both sides of MPA_next. If
# exactly to the edge of one reserve is lost, should shift down
# to next if statement
if(any(next_MPA==1)==FALSE & any(lost==1)==FALSE){
# these are the intervals that are left behind as world moves
# forward
lost <- MPA.current[1:(displaced-1)]
# this is the last continuous chunk of numbers at the end of
#
length_last <- rep(tail(rle(lost)$value,1),
tail(rle(lost)$lengths,1))
# want to know how long is so can make sure to
# get correct interval
L_int <- length(length_last)
# how many more zeros do we need before we start with mpa.yes
# again?
zero_add <- sum(mpa.no==0) - L_int - sum(next_MPA==0)
# this is what needs to be appended to
new_MPA <- c(rep(0,zero_add),
rep(c(mpa.yes,mpa.no),
length.out=length(world)))
}else{
# IF FALSE (there are some 1s) then we only care about the
# last interval of the , is that protected or not?
last_step = tail(next_MPA,1)
# IF ==1
if(last_step ==1){
#then how many 1s are at the end of the section?
end_step = rep(tail(rle(next_MPA)$value,1),
tail(rle(next_MPA)$lengths,1))
# number of 1s at very end of lost interval
# length of
end_int = length(end_step)
# need to prepend sum(mpa.yes==1) - to beginning
prepend = rep(1, (sum(mpa.yes==1) - end_int))
# and then fill out with mpa.no,mpa.yes for length of world
fillOut <- rep(c(mpa.no, mpa.yes),
length = (length(world) - length(prepend)))
}else{
# IF ==0
end_step = rep(tail(rle(next_MPA)$value,1),
tail(rle(next_MPA)$lengths,1))
# number of 0s at very end of lost interval
# length of
end_int = length(end_step)
# need to prepend sum(mpa.no==0) - to beginning
prepend = rep(0, (sum(mpa.no==0) - end_int))
# and then fill out with mpa.no,mpa.yes for length of world
fillOut <- rep(c(mpa.yes, mpa.no),
length = (length(world) ))
}
new_MPA = c(prepend,fillOut)
}
MPA_finish = c(next_MPA,new_MPA)
MPA_finish = MPA_finish[1:length(world)]
# reduce to just the size of the world
return(MPA_finish)
}
m <- function(n, s, Fthresh = NA, Fharv = NA, mpa.yes = NA,
mpa.no = NA, MPA.current=NA,effort_re_allocate=NA){
# steps
# 1. Harvest (check for thresholds, harvesting, MPA coverage)
# 2. Patch moves (and MPAs are adjusted)
# 3. Individuals outside patch die
# 4. Individuals still alive (ie inside the patch) reproduce
# harvesting occurs first - check to see how should
# re-allocate effort
if(!is.na(effort_re_allocate) & !is.na(Fharv)){ # harvesting non-zero and effort reallocate
#total_catch = sum(n)*Fharv
total_catch = sum(n)*Fharv * 1.5 # increasing effort by 50% to account for effort-reallocation
available_total_pop = sum(n[which(MPA.current==0)])
# pop with no MPA coverage
available_fish = rep(0,length(n))
available_fish[MPA.current==0] <-
n[MPA.current==0]/available_total_pop
# available_total_pop
# proportion at each point
catch_in_space <- total_catch*available_fish
# allocate catch
next_gen = n-catch_in_space
next_gen[next_gen<0] = 0
}else{
if(!is.na(Fthresh)) { # if thresholds
next_gen = ifelse(n < Fthresh,
n, n - (n - Fthresh) * Fharv)
}
if(!is.na(Fharv) & is.na(Fthresh)) {
# if harvesting, no thresholds
next_gen = n*(1-Fharv)
}
if(is.na(Fharv) & is.na(Fthresh)) {next_gen = n}
# if no harvesting of any kind
# but put fish back if places that were harvested were in
# the MPA
next_gen[MPA.current == 1] <- n[MPA.current == 1]
}
# move the patch
# calculate how far the patch will move through the population
# (if speed !=0)
displaced = ifelse(s>0,s/step_size,1)
# assign population that will still be inside the patch to
# moved patch
next_n = next_gen[displaced:length(next_gen)]
# fill in newly existing patch with 0s
next_n = c(next_n,rep(0,length.out=(displaced-1)))
# move MPAs?
if(s > 0){MPA_finish = moveMPA(MPA.current, displaced,
mpa.yes,mpa.no,world)}else{MPA_finish= MPA.current}
# let patch reproduce
next_patch = vector(mode="numeric",length(world))
# keep individuals still in patch + those now in it due to
# move
next_patch[1:length(patch)] = next_n[1:length(patch)]
babies = next_patch*f_ind
n2 = babies %*% d *step_size
n2 = sapply(n2,f,R0,K)
MPA = MPA_finish
return(list(n2,MPA)) # removed a mysterious 'harv' from here
}
# wrapper function to run simulation for 6000 generations and save outcome from 2000 additional generations to take average
longRun <- function(s, mpa.yes, mpa.no, Fthresh, Fharv, init,
MPA.start, generations_total, generations_av,
effort_re_allocate=effort_allocate){
MPA.current <- MPA.start
burn_in <- generations_total - generations_av
for(t in 1:(burn_in)){
output = m(n=init, s = s, Fthresh=Fthresh,Fharv=Fharv,
mpa.yes = mpa.yes, mpa.no = mpa.no,
MPA.current = MPA.current, effort_re_allocate=effort_allocate)
init= output[[1]]
MPA.current = output[[2]]
}
# make dataframe for simulation average
pop <- rep(0,generations_av)
for(keep in 1:generations_av){
output = m(n=init, s = s, Fthresh=Fthresh,Fharv=Fharv,
mpa.yes = mpa.yes, mpa.no = mpa.no,
MPA.current = MPA.current, effort_re_allocate = effort_allocate)
init = output[[1]]
MPA.current = output[[2]]
pop[keep] = sum(output[[1]])
}
# take mean for equil_abundance
equil.pop = mean(pop)
equil.sd = sd(pop)
return(list(equil.pop,equil.sd))
}
# to introduce population to empty landscape, harvesting before
# adding speed treatment
# wrapper function to initialize the population, only returns results from final generation
startUp <- function(s, mpa.yes, mpa.no, Fthresh, Fharv, init,
MPA.start, burn_in,effort_re_allocate=effort_allocate){
MPA.current <- MPA.start
for(t in 1:(burn_in)){
output = m(n=init, s = s, Fthresh=Fthresh,Fharv=Fharv,
mpa.yes = mpa.yes, mpa.no = mpa.no,
MPA.current = MPA.current,
effort_re_allocate=effort_allocate)
init= output[[1]]
MPA.current = output[[2]]
}
return(list(init,MPA.current))
}
#------------------------------------------------------------------------#