# run.lake.sim.ADD creates a simulated dataset by sampling TP values of size n, and then simulating the SDE described in Vitense et al. (2017). The function outputs a data frame containing the sampled TP values, as well as starting and ending values for Chla, attracting state, and SAV.


Equilibria <- read.csv("TP_equilibria.csv",header=T) # steady state values for parameters in Table S1
# ********** WARNING: THESE STEADY STATE VALUES ARE ONLY RELEVANT FOR PARAMETERS IN TABLE S1!! IF USERS WISH TO USE DIFFERENT PARAMETERS, THEY NEED TO SUPPLY A NEW TABLE OF EQUILIBRIA AND CHANGE TP THRESHOLD VALUES IN LINES 96-97 IF STATE-DEPENDENT PERTURBATION PARAMETERS ARE USED!! ****************


TPval <- Equilibria$TP
Lower <- Equilibria$Lower   # Lower stable state at TPval (can be NA)
Middle <- Equilibria$Middle # Middle unstable state at TPval (can be NA)
Upper <- Equilibria$Upper   # Upper stable state at TPval (can be NA)

numSS <- apply(Equilibria, 1, function(z) 3-sum(is.na(z))) # vector of number of SS at each TPval
head(cbind(Equilibria,numSS))

### Left bifurcation point:     69.0
# 68.9 3.72212138       NA       NA
# 69.0 3.72642476 20.98533 23.46424

### Right bifurcation point:    303.3
# 303.3 10.506268 12.00103 66.65784
# 303.4        NA       NA 66.66819


run.lake.sim.ADD  <- function(alpha=14, p=8, r=.63, s=.17, pi=270, c=1/200,mulogTP=4.5,sigmalogTP=1,n=130, sigma=c(.2,.2),dt=.001,numdays=30,plotpoints=TRUE, TPdist="lnorm"){
  
  # mulogTP is mean log(TP) and sigmalogTP is SD of log(TP) for sampling from the log-normal distribution to simulate a population of lakes
  # n is number of lakes
  # sigma is standard deviation associated with daily perturbations to Chla, which can be state-dependent (first value is for clear state, second value for turbid)
  
  if (TPdist=="lnorm"){
    # Sample n lakes TP_i values from log normal distribution based on empirical data
    TP0 <- rlnorm(n,meanlog=mulogTP,sdlog=sigmalogTP)
  }
  
  if (TPdist=="unif"){
    # sample from uniform distribution on log scale, then exponentiate. So here the log is UNIFORM
    TP.tmp <- runif(n,0,6.4)
    TP0 <- exp(TP.tmp)
  }
  
  ntemp <- length(TP0[TP0>600])
  TP0[TP0<1] <- 1; TP0[TP0>600] <- runif(ntemp,200,600) # restrict TP to be in [1,600]
  TP0 <- round(TP0,digits=1) # round to one decimal place (accuracy of numerically computed steady states)
  
  # Get current basin of attraction and Chla steady state for t=0 state variables
  SS.bin0 <- numeric(n) ; SS.Chla0 <- numeric(n); Chla0 <- numeric(n)
  MiddleVals <- numeric(n)
  
  for (i in 1:n){
    TP.ind <- which(TPval==TP0[i]) # find TP in the steady state matrix
    Equilibria[TP.ind,] # print the row in the SS matrix
    MiddleVals[i] <- Equilibria[TP.ind,]$Middle
    
    # Clear basin
    # If only one SS at TPval=TP0 and Lower SS exists, it's in clear basin
    if (numSS[TP.ind]==1 & is.na(Lower[TP.ind])==FALSE){
      SS.bin0[i] <- 0  # classify as clear at t=0
      SS.Chla0[i] <- Lower[TP.ind]  # Chla steady state at t=0
      Chla0[i] <- SS.Chla0[i]
    }
    
    # Turbid basin
    # If only one SS at TPval=TP0 and Lower SS exists, it's in turbid basin
    if (numSS[TP.ind]==1 & is.na(Upper[TP.ind])==FALSE){
      SS.bin0[i] <- 1  # classify as turbid at t=0
      SS.Chla0[i] <- Upper[TP.ind]  # Chla steady state at t=0
      Chla0[i] <- SS.Chla0[i] 
    }
    
    # Bistable region 
    # If three SS at TPval=TP0 can be in either basin
    if (numSS[TP.ind]==3){
      SS.bin0[i] <- sample(c(0,1),size=1, prob=c(.6, .4))  # sample from two deterministic SSs (higher probability of being in clear state to reflect empirical distribution)
      if (SS.bin0[i] == 1){ SS.Chla0[i] <- Upper[TP.ind] ; SS.bin0[i] <- 1}
      if (SS.bin0[i] == 0){ SS.Chla0[i] <- Lower[TP.ind] ; SS.bin0[i] <- 0}
      Chla0[i] <- SS.Chla0[i] 
    }
    
  }
  
  
  # Simulate Chla SDE for numdays
  timesteps <- seq(0,numdays,dt) # discretize time over 'numdays' days
  num.inter <- length(timesteps)  # number of time points (including t=0)  
  sqrtdt <- sqrt(dt)       # SD of the stochastic piece
  
  sig <- NULL
  Chla.old <- Chla0
  for (i in 2:num.inter){ 
    
    # This is code to allow perturbation strength to differ in each basin of attraction. Simulations for the paper were run with equal perturbation strength for each basin.
    sig <- as.numeric(Chla.old > MiddleVals)    
    sig[sig==1] <- sigma[2]
    sig[sig==0] <- sigma[1]
    sig[TP0>303.3] <- sigma[2] # above upper TP threshold (this threshold changes with different model parameters)
    sig[TP0<69.0] <- sigma[1]  # below lower TP threshold (this threshold changes with different model parameters)
    
    
    Chla.new <- Chla.old + dt*Chla.old*((r*TP0/(TP0+pi))*(s*(Chla.old^p+alpha^p)/( alpha^p +s*(Chla.old^p+alpha^p)))-c*Chla.old) + (Chla.old^(1))*sig*sqrtdt*rnorm(n)

    Chla.new[Chla.new<.2] <- .2 # don't let Chla go below .2 (.66 is min in DNR data)
    Chla.old <- Chla.new
  }
  
  
  # Get starting and ending SAV (negative sigmoidal relationship to Chla)
  V0 <- (alpha^p)/(Chla0^p+alpha^p)
  V1 <- (alpha^p)/(Chla.new^p+alpha^p)
  
  # Get ending basin of attraction and Chla steady state
  SS.bin1 <- numeric(n) ; SS.Chla1 <- numeric(n)
  for (i in 1:n){
    TP.ind <- which(TPval==TP0[i]) # find TP value in the steady state matrix
    Equilibria[TP.ind,] # print the row in the SS matrix
    
    # Clear basin
    # If only one SS at TPval=TP0 and Lower SS exists, it's in clear basin
    if (numSS[TP.ind]==1 & is.na(Lower[TP.ind])==FALSE){
      SS.bin1[i] <- 0  # classify as clear
      SS.Chla1[i] <- Lower[TP.ind]  # Chla steady state
    }
    # Turbid basin
    # If only one SS at TPval=TP0 and Lower SS exists, in turbid basin
    if (numSS[TP.ind]==1 & is.na(Upper[TP.ind])==FALSE){
      SS.bin1[i] <- 1  # classify as turbid 
      SS.Chla1[i] <- Upper[TP.ind]  # Chla steady state
    }
    
    # Bistable region 
    # If three SS at TPval=TP0 can be in turbid or clear basin
    if (numSS[TP.ind]==3){  
      # Classify basin after perturbation
      if (Chla.new[i] > Middle[TP.ind]){ SS.Chla1[i] <- Upper[TP.ind] ; SS.bin1[i] <- 1}
      if (Chla.new[i] < Middle[TP.ind]){ SS.Chla1[i] <- Lower[TP.ind] ; SS.bin1[i] <- 0}
    }
    
  }
  
  # If desired, plot starting and ending Chla values for each lake on the bifurcation diagram
  if(plotpoints==TRUE){
    # Make bifurcation diagram
    
    Ps <- TPval
    
    plot(Ps,Lower,xlim = range(Ps),ylim=c(0,90),type="l",xlab="TP",ylab="Chla steady state",las=1,main=paste("Start-End Chla levels after ",numdays," days: sigma=",sigma),lwd=2, col="royalblue4")
    lines(Ps,Middle,type="l",lty=2, lwd=2, col="gray60")
    lines(Ps,Upper,type="l", lwd=2, col="seagreen")
    legend("bottomright",pch=c(1,16,NA,NA),lty=c(NA,NA,1,2),c("Start","End","Stable","Unstable"),bty="n",cex=.9,lwd=2)
    
    # Plot starting and ending points
    SS0.col <- rep("royalblue4",n); SS0.col[SS.bin0==1] <- "seagreen"
    SS1.col <- rep("royalblue4",n); SS1.col[SS.bin1==1] <- "seagreen"
    
    points(TP0,Chla0, col=SS0.col)  # starting points (open circles)
    points(TP0,Chla.new,pch=16, col=SS1.col) # ending points (filled circles)
    for (i in 1:length(TP0)){
      lines(c(TP0[i],TP0[i]),c(Chla0[i], Chla.new[i])) # add lines to connect them
    }
    
    plot(log(Ps),log(Lower),xlim = c(0,7),ylim=c(-4,5),type="l",xlab="log(TP)",ylab="log(Chla) steady state",las=1,main="log(Chla) Steady State vs. log(TP)",lwd=2, col="royalblue4")
    lines(log(Ps),log(Middle),type="l",lty=2, lwd=2, col="gray60")
    lines(log(Ps),log(Upper),type="l", lwd=2, col="seagreen")
    legend("bottomright",pch=c(1,16,NA,NA),lty=c(NA,NA,1,2),c("Start","End","Stable","Unstable"),bty="n",cex=.9, lwd=2)
    points(log(TP0),log(Chla0), col=SS0.col)  # starting points (open circles), t=0
    points(log(TP0),log(Chla.new),pch=16, col=SS1.col) # ending points (filled circles), t=1
    for (i in 1:length(TP0)){
      lines(c(log(TP0[i]),log(TP0[i])),c(log(Chla0[i]), log(Chla.new[i]))) # add lines to connect them
    }
  }
  
  SAV1 <- round(V1,digits=10)
  SAV1[SAV1==1] <- .9999999999
  SAV1[SAV1==0] <- .0000000001
  
  SummaryStates <- data.frame(TP=TP0,SAV0=round(V0,digits=2),TrueBasin0=SS.bin0,TrueSS_Chla0=round(SS.Chla0,digits=2),Chla1=round(Chla.new,digits=2),SAV1=SAV1,TrueBasin1=SS.bin1,TrueSS_Chla1=round(SS.Chla1,digits=2))
  # TP = lake i's TP
  # SAV0 = starting SAV
  # TrueBasin0 = classified starting steady state (0=clear, 1=turbid)
  # TrueSS_Chla0 = value of Chla at the starting steady state
  # Chla1 = ending Chla concentration
  # SAV1 = ending SAV
  # TrueBasin1 = classified ending steady state (0=clear, 1=turbid)
  # TrueSS_Chla1 = ending Chla steady state
  
  return(SummaryStates=SummaryStates)
  
}