Lecture13 <- function()
{
 set.seed(1954)
# Base <- Part1()
# Part2(Base)
# Part3(Base)
# Part4()
# Part5(Base)
 Part6()


}

# ================================================================================

DoProject <- function(r,Start_val=500,Thresh=10,Nyears=100,Nsim=1000,SigP=0.2)
{
 WentExtinct <- 0
 for (Isim in 1:Nsim)
  {
   Nvalue <- Start_val
   Extinct <- F
   for (Iyr in 1:Nyears)
    {
     Nvalue <- r*Nvalue*exp(rnorm(1,0,SigP) - SigP^2/2)
     if (Nvalue < Thresh) Extinct <- T
    }
   if (Extinct==T) WentExtinct <- WentExtinct + 1 
  }
 return(WentExtinct/Nsim*100)

}

# ================================================================================
Part1 <- function()
{
 xx <- c(0.85, 0.9, 0.95, 1, 1.05,1.1, 1.15)
 yy <- rep(0,length=length(xx))
 for (II in 1:length(xx))
  yy[II] <- DoProject(xx[II])
 print(yy) 
 
 par(mfrow=c(2,2))
 plot(xx,yy,xlab="r",ylab="Probability (extinction)",type='b',lty=1,lwd=2,pch=16)   

 return(yy)

}

# ================================================================================

DoProject2 <- function(r,Start_val=500,Thresh=10,Nyears=100,Nsim=1000,SigP1=0.2,SigP2=0.4,ProbCat=0.1)
{
 WentExtinct <- 0
 for (Isim in 1:Nsim)
  {
   Nvalue <- Start_val
   Extinct <- F
   for (Iyr in 1:Nyears)
    {
     Error <-rnorm(1,0,1)
     if (runif(1,0,1) <= ProbCat)
      Error <- Error*SigP2 - SigP2^2/2
     else 
      Error <- Error*SigP1 - SigP1^2/2
     Nvalue <- r*Nvalue*exp(Error)
     if (Nvalue < Thresh) Extinct <- T
    }
   if (Extinct==T) WentExtinct <- WentExtinct + 1 
  }
 return(WentExtinct/Nsim*100)

}

# ================================================================================
Part2 <- function(Part1)
{
 xx <- c(0.85, 0.9, 0.95, 1, 1.05,1.1, 1.15)
 yy <- rep(0,length=length(xx))
 for (II in 1:length(xx))
  yy[II] <- DoProject2(xx[II])
 print(yy)
 
 par(mfrow=c(2,2))
 plot(xx,Part1,xlab="r",ylab="Probability (extinction)",type='b',lty=1,lwd=2,pch=16)   
 lines(xx,yy,lty=3,lwd=2)

 return(yy)

}

# ================================================================================

DoProject3 <- function(r,Start_val=500,Thresh=10,Nyears=100,Nsim=1000,SigP=0.2,Prow=0.717)
{
 WentExtinct <- 0
 for (Isim in 1:Nsim)
  {
   Nvalue <- Start_val
   Extinct <- F
   Error <- rnorm(1,0,SigP)
   for (Iyr in 1:Nyears)
    {
     Error <- Prow*Error + sqrt(1-Prow^2)*rnorm(1,0,SigP)     
     Nvalue <- r*Nvalue*exp(Error - SigP^2/2)
     if (Nvalue < Thresh) Extinct <- T
    }
   if (Extinct==T) WentExtinct <- WentExtinct + 1 
  }
 return(WentExtinct/Nsim*100)

}

# ================================================================================
Part3 <- function(Part1)
{
 xx <- c(0.85, 0.9, 0.95, 1, 1.05,1.1, 1.15)
 yy <- rep(0,length=length(xx))
 for (II in 1:length(xx))
  yy[II] <- DoProject3(xx[II])
 print(yy)
 
 par(mfrow=c(2,2))
 plot(xx,Part1,xlab="r",ylab="Probability (extinction)",type='b',lty=1,lwd=2,pch=16)   
 lines(xx,yy,lty=3,lwd=2)

 return(yy)

}

# ==============================================================================================================================

Part4 <- function()
{
 par(mfrow=c(2,2))
 xx <- seq(from=0,to=1000,length=101)
 Surp <- xx*(1-xx/1000)
 Surp2 <- xx*(1-xx/1000)*(1-0.5^(xx/100))
 Surp3 <- xx*(1-xx/1000)*(1-0.5^(xx/200))

 plot(xx,Surp,xlab="Population size",ylab="Surplus production",lty=1,type='l',lwd=2,yaxs="i",ylim=c(0,400),xaxs="i")
 lines(xx,Surp2,lty=2,lwd=2)
 lines(xx,Surp3,lty=3,lwd=2)
 legend("topright",legend=c("No depensation",expression(N[50]==100),expression(N[50]==200)),lty=c(1:3),lwd=2,bty="n")
    
}
# ================================================================================

DoProject5 <- function(r,Start_val=500,Thresh=10,Nyears=100,Nsim=1000,SigP=0.2,Carry=1000,N50=200)
{
 WentExtinct <- 0
 for (Isim in 1:Nsim)
  {
   Nvalue <- Start_val
   Extinct <- F
   for (Iyr in 1:Nyears)
    {
     Error <-exp(rnorm(1,0,SigP) - SigP^2/2)
     Nvalue <- (Nvalue+(r-1)*Nvalue*(1-Nvalue/Carry)*(1-0.5^(Nvalue/N50)))*Error
     if (Nvalue < 0.1) Nvalue <- 0.1
     if (Nvalue < Thresh) Extinct <- T
    }
   if (Extinct==T) WentExtinct <- WentExtinct + 1 
  }
 return(WentExtinct/Nsim*100)

}

# ================================================================================
Part5 <- function(Part1)
{
 xx <- c(0.85, 0.9, 0.95, 1, 1.05,1.1, 1.15)
 yy <- rep(0,length=length(xx))
 for (II in 1:length(xx))
  yy[II] <- DoProject5(xx[II])
 print(yy)
 
 par(mfrow=c(2,2))
 plot(xx,Part1,xlab="r",ylab="Probability (extinction)",type='b',lty=1,lwd=2,pch=16)   
 lines(xx,yy,lty=3,lwd=2)

 return(yy)

}

# ==============================================================================================================================

Lynx <- function(K,s24,s1,Ninit,ProbFemCub=0.87,ProbCat=0.1,Thresh=2,Nsim=100,Nyear=50)
{
 WentExtinct <- 0
 for (Isim in 1:Nsim)
  {
    # Initialize
    n <- Ninit
    Extinct <- F
    for (Iyr in 1:Nyear)
     {
      # Define year-1 survival
      if (runif(1,0,1) < ProbCat) SurvJ <- s1[2] else SurvJ <- s1[1]
      s <- c(SurvJ,s24)
      
      # Update the stage-structure
      Nnext <- rep(0,4)
      for (Istage in 1:4)
       Nnext[Istage] <- rbinom(1,n[Istage],s[Istage])
       
      # Generate the cub 
      Cubs <- if (n[4] > 0) rbinom(1,n[4],ProbFemCub) else 0
      
      # Probaby becoming a breeder
      if (Nnext[2]+Nnext[3] > 0) ProbBreed <- (K-Nnext[4])/(Nnext[2]+Nnext[3]) else ProbBreed <- 0.5
      if (ProbBreed < 0) ProbBreed <- 0
      if (ProbBreed > 1) ProbBreed <- 1
      Ntemp <- Nnext[2]+Nnext[3]
      NewBreed <- if (Ntemp > 0) rbinom(1,Ntemp,ProbBreed) else 0
      n[4] <- Nnext[4] + NewBreed
      n[3] <- Nnext[2] + Nnext[3] - NewBreed
      n[2] <- Nnext[1]
      n[1] <- Cubs

      # Check for extinction
      Ntot <- sum(n)
      if (Ntot <= Thresh) Extinct <- T
     }
    if (Extinct==T) WentExtinct <- WentExtinct + 1 
  }
 return(WentExtinct/Nsim*100)

}


# ==============================================================================================================================
Part6 <- function()
{
 xx <- c(1,2,3,4,5,6,7,8,9,10)
 yy <- rep(0,length=length(xx))
 for (II in 1:length(xx))
  yy[II] <- Lynx(xx[II],c(0.7,0.7,0.86),c(0.5,0.2),Ninit=c(4,2,0,5),Nsim=1000)
 print(yy)
# yy <- c(100,81.9,41.8,21.9,13.4,6.7,6.3,3.4,2.7,2.7)
 
 par(mfrow=c(2,2))
 plot(xx,yy,xlab="Number of territories",ylab="Probability (extinction)",type='b',lty=1,lwd=2,pch=16)   
    
}

# ==============================================================================================================================

Lecture13()