#' ---
#' title:  "Fishing herring"
#' author: "Rob J. de Boer"
#' output:
#'   pdf_document: default
#' ---

source("grind.R")

#' ## a.
#' This function is a parabola that is zero when *N=0* and when *N=K*.
#' Because we are only interested in the shape of *f(N)* we set *r=K=1*,
#' and we sketch the function *f(N)=rN(1-N/K)* as follows:
curve(N*(1-N),xname="N",lwd=2,col="red")
points(0,0,pch=1); points(1,0,pch=16)
#' which has a maximum of *rK/4=0.25* when *N=K/2=0.5*.
#' 

#' ## b.
#' The maximum harvest is *Q=rK/4*, 
#' which is obtained with a herring population of *H=K/2*.
#' 

#' ## c.
#' The function *f(N)=rN(1-N/K)-rK/4* was sketched for *r=K=1* as follows:
curve(N*(1-N)-0.25,xname="N",lwd=2,col="red")
lines(c(0,1),c(0,0))
points(0.5,0,pch=1)
#' which is negative except when *N=K/2=0.5*.
#' Catching less than that, e.g., *Q=0.2*, 
#' allows for two steady states where *f(N)=0*:
curve(N*(1-N)-0.2,xname="N",lwd=2,col="red")
lines(c(0,1),c(0,0))
points(0.275,0,pch=1); points(0.725,0,pch=16)
#' The upper steady state is stable because the slope of *f(N)* is negative,
#' and the lower steady state is unstable because the slope of *f(N)* is positive.
#' 

#' ## d.
#' Defining *Q* as the harvest, the model becomes *dH/dt=rH(1-H/K)-Q*,
#' which was defined in Grind as follows:
model <- function(t, state, parms) {
  state <- ifelse(state < 0, 0, state)
  with(as.list(c(state,parms)), {
    dH <- r*H*(1 - H/K) - Q  
    return(list(dH))  
  }) 
}  
#'The parameters 
p <- c(r=0.1,K=1,Q=0.025)
#' define a scaled carrying capacity of *H=K=1*,
#' and the growth rate *r=0.1* defines a growth rate of 10% per year.
#' The parameter *Q=rK/4=0.025* defines the optimal quota leading to *H=K/2=0.5*.
#' We study such a harvested population starting at *H=K/2=0.5* with and without noise:
s <- c(H=0.5)
z <- run(100,xlab="Time in years")
after="parms[\"K\"]<-abs(rnorm(1,mean=1,sd=0.25));state[1]<-ifelse(state[1]<0,0,state[1])"
z <- run(100,tstep=1/52,after=after,add=TRUE)
z <- run(100,tstep=1/52,after=after,add=TRUE)
z <- run(100,tstep=1/52,after=after,add=TRUE)
#' Thus, running the model for 100 years with a carrying capacity that changes weekly
#' due to changing weather conditions ultimately leads to extinction.
#'

#' ## e.
#' Plotting *f(N)=rH(1-H/K)-fH* for *r=K=1* and the optimal *f=r/2* gives
curve(N*(1-N)-N/2,xname="N",from=0,to=1,ylim=c(0,0.1),lwd=2,col="red")
points(0,0,pch=1); points(0.5,0,pch=16)
#' which reveals that *N=K/2* becomes a stable steady state because the slope
#' of *f(N)* in *N=K/2* is negative.
#' 

#' ## f.
#' The new model *dH/dt=rH(1-H/K)-fH* is defined in Grind as
model <- function(t, state, parms) {
  state <- ifelse(state < 0, 0, state)
  with(as.list(c(state,parms)), {
    dH <- r*H*(1 - H/K) - f*H  
    return(list(dH))  
  }) 
}
#' and for the optimal harvest we set *f=r/2=0.05*
p <- c(r=0.1,K=1,f=0.05)
#' Running the model under the same conditions as before,
#' i.e., with the same total harvest and the same weather fluctuations,
z <- run(100,xlab="Time in years")
z <- run(100,tstep=1/52,after=after,add=TRUE)
z <- run(100,tstep=1/52,after=after,add=TRUE)
z <- run(100,tstep=1/52,after=after,add=TRUE)
#' does not lead to extinction because the state *H=0.5* is stable
#' Do note that the herring population is always lower than *H=K/2=0.5*.
#' 
#' We conclude that it is much more sustainable to catch the optimal fraction
#' of the herring population, rather than the optimal number of herring, and 
#' that this does not have any economic cost, as on average the same amount
#' of herring can be harvested.