R/timeSteps.R

In gabackus/rcomms: Simulate Randomized Metacommunities under Stochastic Climate Change

#' Simulate the Model
#'
#' This simulates a community over a discrete number of years.
#' @param n Initial population sizes. \code{S}x\code{L}x\code{W} array of population nonnegative integers.
#' @param P List of biotic variables. See \code{\link{commSetup}}.
#' @param X List of abiotic variables. See \code{\link{commSetup}}.
#' @param y Current year.
#' @param years How many time steps to run the model.
#' @param init Is the model being initialized? If TRUE, then Tau is set to 0 (no climate chage).
#' @param extInit If \code{TRUE}, the simulation stops running after extThresh time steps are completed without any extinctions. When attempting to initialize a stable community, consider setting extInit to T.
#' @param extThresh If \code{extInit==TRUE}, then the simulation stops once extInit time steps have passed without any extinctions.
#' @param manage A list with a bunch of management options. See \code{\link{manageSetup}} If left blank, the model generates a list where all management options are set to \code{FALSE}.
#'
#' @return A list with \code{n} (\code{S}x\code{L}x\code{W}array of final population sizes over space), \code{N} (\code{S}x\code{Y} array of population sizes over time), \code{temps} (\code{Y} vector of mean temperature over time)
#' @export
commSimulate <- function(n,P,X,y=1,years=100,init=F,extInit=F,extThresh=100,manage=NULL){
  # Make an all-FALSE management list if none is provided
  if(is.null(manage)){
    manage <- manageSetup(P,X)
  }
  if(init==T){
    Tau <- 0
  } else{
    Tau <- X$Tau
  }

  # First, we set up a matrix to save the total population size over time
  N <- matrix(0,P$S,years+1)
  # Record the total initial population size of each species across the whole ecosystem
  N[,1] <- apply(n,1,sum)

  # Temperature changes over time, so we need to adjust this over the course of the model
  temp2d0 <- X$temp2d

  # For output, we want to keep track of the average temperature over time, and we do that with temps
  temps <- seq(0,years)
  temps[1] <- mean(temp2d0+X$tempY[y])

  # Keep track of tempY and tau outside of X
  tempY <- X$tempY

  # Run the model for a number of time steps equal to 'years'
  for (i in 1:years){
    # Temperature changes before each time step
    X$temp2d=temp2d0+Tau*i+tempY[i+y]
    # save the mean temperature
    temps[i+1]=mean(X$temp2d)
    # Run the time step function for population adjustment after the change in temperature
    if(sum(N[,i])>0){
      n <- timeStep(n,P,X,N[,i],i,manage)
    }
    # Record the population size
    N[,i+1]<-apply(n,1,sum)
  }
  return(list(n=n,N=N,temps=temps))
}

#' Simulate One Time Step
#'
#' Cycle through each step of the model.
#' Each time step could be one "year" or one "generation", but ultimately it runs through each part of the life cycle.
#' The various management techniques can optionally be added to the model between any two of the required steps.
#' reproduction -> dispersal -> mortality -> competition
#' @param n Initial population sizes. \code{S}x\code{L}x\code{W} array of population nonnegative integers.
#' @param P List of biotic variables. See \code{\link{commSetup}}.
#' @param X List of abiotic variables. See \code{\link{commSetup}}.
#' @param N Population size over time
#' @param t Current time
#' @param manage A list with a bunch of management options. See \code{\link{manageSetup}}. If left blank, the model generates a list where all management options are set to \code{FALSE}.
#'
#' @return \code{S}x\code{L}x\code{W} array of population sizes at the end of the time step
#' @export
timeStep <- function(n,P,X,N,t,manage){
  n1 <- reproduce(n,X$L,X$W,P$S,P$z,P$sig,P$ro,P$lam,X$temp2d)

  # If assisted migration is occurring in this simulation
  if(manage$AM$AM==T){
    am <- assistMigrate(n1,N,X$L,X$W,X$temp1d,t,manage$AM)
    # These are the individuals (of all species) that will still be relocating normally
    n1 <- am$nD
    # These are the individuals that were relocated
    nR <- am$nR
    # Update the time since last relocation vector
    manage$AM$tLR <- am$tLR
    # Make note of the times when the species was relocated
    manage$AM$relTimes[am$SReloc,t] <- 1
  }

  # To save computation time, we don't need to use the dispersal function on species without any propagules
  # dS are the species with extant propagules
  dS<-which(rowSums(n1)>0)
  # Thus, as long as there is at least one species dispersing, go through with dispersal on dS species
  if(!pracma::isempty(dS)){
    # Slice the array so we are only using those that will be dispersing
    dispn <- n1[dS,,,drop=F]
    # Now run the disperse function on the slice
    dispn2 <- disperse(dispn,X$L,X$W,P$S,P$K[dS])
    # Preallocate the dispersed array
    n2 <- n1*0
    # Add the dispersed individuals to that array
    n2[dS,,] <- dispn2
  } else {
    # If there are no propagules at all, n2 is just n1
    n2<-n1
  }

  # The survive function simulates mortality of adults
  n3 <- n2+survive(n,X$L,X$W,P$S,P$M)

  # All inidividuals then compete
  # (At this point, adults and offspring have equal competitive ability. We could change the compete function if this should change.)
  n4 <- compete(n3,X$L,X$W,P$S,X$Q,P$A,P$compType,P$z,P$sig,P$ro,P$lam,X$temp2d)

  return(n4)
}

#' Reproductive Performance
#'
#' @export
bi <- function(z,sig,ro,lam,temp){
  op<-ro*(exp(-((temp-z)/sig)^2)*(1+pracma::erf(lam*(temp-z)/sig))-1)
  return(op)
}

#' All Individuals Reproduce
#'
#' The number of offspring born for each species in each location is a Poisson random variable with mean \code{r*n}.
#' The base reproductive rate is a skewed function, adjust such that \code{min(r)=0} and \code{max(r)=2}.
#' Each species will have a different reproductive rate depending on the temperature at that space.
#'
#' @export
reproduce <- function(n,L,W,S,z,sig,ro,lam,temp2d){
  # r is the "continuous" form of the birth rate
  r <- sapply(1:S, function(i) bi(z[i],sig[i],ro[i],lam,temp2d))
  # R turns it into a discrete form
  R <- exp(r)
  # This just turns it into the correct array format
  R <- aperm(array(R,c(L,W,S)),c(3,1,2))

  # Mean number of offspring
  rn <-c(R*n)

  # The number of offspring is a Poisson random variable with mean=R*n
  nr<-array(sapply(rn, function(x) rpois(1,x)),c(S,L,W))

  return(nr)
}

#' Individuals Disperse Outwards
#'
#' Each individual spreads throughout the spatial landscape with a random double geometric dispersal kernel determined by the species' mean dispersal distance, gami.
#' For each species in each location, disperse!
#'
#' @export
disperse <- function(n,L,W,S,K){
  # Si is the total number of species that are dispersing
  Si <- nrow(n)

  if(is.null(Si)){
    # When there is 1 species, this function gets confused, so we can fix it here
    n <- array(n,c(S,L,W))
    Si <- S
  }

  # Flatten the metapopulation so that all microhabitats in one patch are summed together
  # This makes n1 an SxL matrix
  n1 <- apply(n,c(1,2),sum)
  # This disperses the propagules with multinomial random vectors across the temperature gradient X
  n2 <- t(sapply(1:Si, function(j) disperseMulti(n1[j,],L,K[[j]])))
  # This distributes the propagules randomly into the microhabitats for each patch x
  n3 <- c(sapply(1:Si, function(i) t(sapply(1:L,function(j) rebin(sample(1:W,n2[i,j],replace=T),W)))))
  # This just reformats n3 so it is in the previous SxLxW form
  n4 <- aperm(array(n3,c(L,W,Si)),c(3,1,2))

  # And now it's ready to go
  return(n4)
}


#' Function Used in Dispersal Function
#'
#' Used in in the disperse function
#' To save computation time, we only use the mulinomial random number generator for patches where local \code{n} is positive.
#'
#' @export
disperseMulti <- function(n,L,K){
  # y is just there to mark which indices we are going to run the multinomial generator
  y <- which(n>0)
  # Run the multinomial random generator
  n1 <- sapply(y,function(x) rmultinom(1,n[x],K[x,]))

  # Now we add all of these vectors together
  if(length(y)>1){
    # (Assuming that propagules dispersed from more than one patch)
    n2 <- rowSums(n1)
  } else{
    # (Otherwise, no summation is necessary)
    n2 <- n1
  }
  # This just cuts off the edges that are removed from the model (since we have absorbing boundaries)
  n3 <- n2[2:(L+1)]
  return(n3)
}

#' Adults Die or Survive
#'
#' Adults survival is stochastic.
#'
#' @export
survive <- function(n,L,W,S,M){
  # First, we flatten n so it is one long vector (SLWx1) (just like M)
  nv <- c(n)

  # We can save computation time if we skip over cases where all species have 0 or 1 mortality probabilities.
  if(all(M==1)){
    # When all M is 1, all adults die
    ns <- n*0
  } else if(all(M==0)){
    # When all M is 0, all adults live
    ns <- n
  } else{
    # To save computation time, we find out which patches have living adults with some probability of surviving
    # wMort are all of the patches that fit this
    wMort <- which(nv>0 & M<1)

    # nsv is the full vector form of surviving adults. Pre-allocated to 0 for all.
    nsv <- 0*nv

    # Assuming that at least one patch with adults that might survive, calculate stochastic survival
    if(!pracma::isempty(wMort)){
      nsi<-sapply(wMort, function(i) rbinom(1,nv[i],1-M[i]))
      nsv[wMort]<-nsi
    }
    # Convert the output into original SxLxW form
    ns<-array(nsv,c(S,L,W))
  }

  # ns is all of the surviving adults
  return(ns)
}


#' Individuals Compete over Space
#'
#' The density dependence in this model is roughly a Beverton-Holt model that includes both interspecific and intraspecific competition.
#' Each individual has a random chance of survival based on a variety of conditions.
#' Competition coefficients depend on interactions between each species and the temperature at the location at the time.
#' These can be thought of as temperature-varying Lotka-Volterra competition coefficients.
#' Probability of survival depends on competition coefficients, number of individuals of each different species at that location, and the quality of the habitat at that location.
#'
#' @export
compete <- function(n,L,W,S,Q,A,compType='lottery',z=NULL,sig=NULL,ro=NULL,lam=NULL,temp2d=NULL){
  # Competition works differently dependenting on whether it is temperature-dependent or pure lottery competition
  if(compType=="temp"){
    # Use the same reproduction temperature dependence to determine competitive pressure
    r <- sapply(1:S, function(i) bi(z[i],sig[i],ro[i],lam,temp2d))
    R <- aperm(array(exp(r),c(L,W,S)),c(3,1,2))
  } else if(compType=="lottery"){
    # All individuals are equal
    R <- 1
  }

  # Convert Q into an SxLxW array
  Qrep <- array(rep(Q,each=S),c(S,L,W))
  # QR determines habitat quality the species
  QR <- 1/(Qrep*R)
  # nR is used to determine species interactions
  nR <- R*n
  # This puts the species interactions together
  anR <- sapply(1:S, function(s) colSums(A[s,]*nR))
  anR <- aperm(array(anR,c(L,W,S)),c(3,1,2))

  # Convert this into survival probability
  p <- 1/(1+QR*anR)

  # Binomial random variables to see who survives
  nc <- (sapply(1:S,function(s) mapply(rbinom,1,c(n[s,,]),p[s,,])))
  # Converted into proper SxLxW form
  nc2 <- array(t(nc),c(S,L,W))

  return(nc2)
}

#' Assisted Migration Function
#'
#' Relocate populations based on the AM management strategy.
#' First check to see if any species need to be relocated (are they below the threshold and not within the cooldown period).
#' Then we pick which individuals and how many from the donor population depending on donor strategy.
#' Then some individuals do not survive assisted migration.
#' Then they are moved to a target location depending on how we define target strategy.
#'
#' @export
assistMigrate<-function(n,N,L,W,temp1d,t,AM){
  # Attach AM parameters
  tLR <- AM$tLR; targs <- AM$targs; eta <- AM$eta; tCD <- AM$tCD

  # Preallocate the output arrays
  # nR is the array of relocated individuals
  nR <- n*0
  # nD is the array of individuals that will disperse naturally instead of assisted migration
  nD <- n


  ##########################################################
  # DO WE DO ANY ASSISTED MIGRATION DURING THIS TIME STEP? #
  ##########################################################
  # Which species need to be relocated?
  # Must be a target species, not during a cooldown period, and population less than the threshold but greater than 0
  SReloc <- which(targs & t-tLR>tCD & N<eta & N>0)
  SRL <- length(SReloc)

  # We only need to go through this bit if there are going to be any relocations during this time step
  if(!pracma::isempty(SReloc)){

    # Attach AM parameters
    rho <- AM$rho; mu <- AM$mu; zEst <- AM$zEst; xLoc <- AM$xLoc; recRad <- AM$recRad; donor <- AM$donor; recipient <- AM$recipient; randPick <- AM$randPick
    # Preallocate the relocation array
    nXWReloc <- n[SReloc,,,drop=FALSE]*0

    # Because relocation is occuring, we can update the tLR (time of last relocation) vector
    tLR[SReloc] <- t

    # We don't need to worry about microhabitats here, so we can flatten out the population array a bit
    # This makes n1 an SxL matrix
    nf <- apply(n,c(1,2),sum)

    # To help with this, set up the populations for each SReloc species into a vector for all patches with microhabitats summed up (SReloc x L)
    nv <- sapply(SReloc, function(s) c(nf[s,]))

    NSProps <- colSums(nv)
    ##################################################
    # HOW MANY AND WHICH INDIVIDUALS DO WE RELOCATE? #
    ##################################################
    # Do we pick individuals randomly or just a take a particular amount?
    # NDon is the total number of donated propagules for each species in SReloc
    if(randPick){
      NDon <- sapply(1:SRL, function(s) rbinom(1,NSProps[s],rho[s]))
    } else{
      NDon <- ceiling(NSProps*rho[SReloc])
    }


    # Now we pick the actual individuals out of the metapopulations
    # We can pick them in different ways
    if(donor==1){
      # 1 is randomly picked throughout the entire range
      # We can do this with a multivariate hypergeometric random vector
      nDon <- t(sapply(1:SRL, function(s) rmvhyper(1,nv[,s],NDon[s])))
    }else if(donor==2){
      # 2 is picking individuals from the trailing edge
      # First, we preallocate the nDon
      nDon <- nv*0
      for(s in 1:SRL){
        # We convert the spatial distribution of local populations into a vector that just shows where each individual is
        nUnbin <- unbin(nv[,s])
        # Pick the first NDon[,s] on the trailing edge
        nUnbinDon <- nUnbin[1:NDon[s]]
        # and put it back into the regular format
        nDon[s,] <- rebin(nUnbinDon,length(nv[,s]))
      }
    }else if(donor==3){
      # 3 is picking individuals from the leading edge
      # First, we preallocate the nDon
      nDon <- nv*0
      for(s in 1:SRL){
        # We convert the spatial distribution of local populations into a vector that just shows where each individual is
        nUnbin <- unbin(nv[,s])
        # Pick the first NDon[,s] on the leading edge
        nUnbinDon <- nUnbin[(length(nUnbin)-NDon[s]+1):length(nUnbin)]
        # and put it back into the regular format
        nDon[s,] <- rebin(nUnbinDon,length(nv[,s]))
      }
    }

    # Remove the donated indivudals from the nv array
    nvDisp <- nv-t(nDon)
    # Preallocate an array for relocated individuals
    nvReloc <- nv*0

    ############################
    # WHO SURVIVES RELOCATION? #
    ############################
    # Total individuals surviving
    NDonS <- sapply(1:SRL, function(s) rbinom(1,NDon[s],mu[s]))

    #########################
    # WHERE DO WE RELOCATE? #
    #########################
    # Which patch is closest to the estimated thermal optimum + xLoc patches ahead
    locS <- sapply(SReloc, function(s) which.min(abs(zEst[s]-temp1d))+xLoc[s])
    # If locS is too big or too small, we need to fix that
    for(s in 1:SRL){
      if(locS[s]<(1+recRad[SReloc[s]])){
        loc[s] <- 1+recRad[SReloc[s]]
      }else if(locS[s]>(L-recRad[SReloc[s]]))
        locS[s] <- L-recRad[SReloc[s]]
    }

    # Identify the locations that receive inidividuals
    locSs <- sapply(1:SRL, function(s) (locS[s]-recRad[SReloc[s]]):(locS[s]+recRad[SReloc[s]]))

    # Relocate the individuals based on shape
    if(recipient==1){
      # 1 is a square shape
      for(s in 1:SRL){
        # Length of recipient location
        recSize<-(2*recRad[s]+1)
        # Evenly distribute individuals through the recipient location
        nvReloc[locSs[,s],s] <- nvReloc[locSs[,s],s]+floor(NDonS[s]/recSize)
        # There will probably be a few extras after all is even. Place these randomly.
        extras <- sample(locSs[,s],NDonS[s]%%recSize)
        nvReloc[extras,s] <- nvReloc[extras,s]+1

        # Now we redistribute these amongst the width of the ecosystem
        # Similarly to above
        for(x in locSs[,s]){
          nXWReloc[s,x,] <- floor(nvReloc[x,s]/W)
          extras <- sample(1:W,nvReloc[x,s]%%W)
          nXWReloc[s,x,extras] <- nXWReloc[s,x,extras]+1
        }
      }
    }
    # Now put the relocated and dispersing individuals into full arrays
    nR[SReloc,,] <- nXWReloc
    # The non-relocated individuals will disperse naturally
    # When running the disperse() function, all microhabitats are combined, so there's no reason to separate these into microhabitats
    nD[SReloc,,1] <- t(nvDisp)
  }
  # Ultimately, we want to output the relocated array, the dispersing array, the species that were relocated, and the updated time since last relocation
  output <- list(nR=nR,
                 nD=nD,
                 SReloc=SReloc,
                 tLR=tLR)
  return(output)
}