R/miscFunctions.R

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

#' Create Heterogeneity with Fractional Brownian Nnoise
#'
#' This function adds some heterogeneity to the temperature gradient with fractional Brownian noise.
#' See Keitt (2000), Spectral representation of neutral landscapes, Landscape Ecology 15.
#' @param L Length of the temperature gradient.
#' @param H Hurst exponent (should be between 0 and 1). It relates to the autocorrelation.
#' When H is near 1, this function has positive long-term positive autocorrelation and will look relatively smooth.
#' When H=0.5, this function is Brownian motion.
#' When H is near 0, then autocorrelation is negative and positive values will more often be followed by negative values (and vice versa).
#' @param cZero If TRUE, it will center the whole output so the mean is 0.
#'
#' @return A heterogeneous vector
#' @export
tempVarH <-  function(L,H,cZero=T){
  # random phases uniformly distributed on [0,2pi]
  phif <- runif(L)*2*pi

  # adjusted exponent for amplitudes
  betaH <- 1+2*H

  # uniformly distributed random numbers
  xf <- rnorm(L)
  # to form the amplitudes
  af <- 1/seq(1,L)^(betaH/2)*xf

  # complex coeffcients
  cf <- af*exp(1i*phif)

  #  real part of the inverse fourier transform
  tH <- Re(pracma::ifft(cf))

  # center it around zero?
  if(cZero){
    tH <- tH-mean(tH)
  }

  # multiply the output to increase the magnitude of the heterogeneity
  # add that to the the temperature gradient
  return(tH)
}


#' Probability Mass Function for Double Geometric Distribution
#'
#' Probability mass function for "double geometric" distribution.
#' @param x Distance from origin to landing spot. Can be integer or vector of integers.
#' @param q Probability of remaining in a given patch (and not continuing to move); see supplemental
#'
#' @return Probability of landing in location x.
#' @export
doubGeom<-function(x,q){
  return((q/(2-q)*(1-q)^abs(x)))
}

#' Adjust Reproduction for Tradeoff
#'
#' This function creates a constant to adjust the reproduction rate so that the area under the curve is roughly equal for all species.
#' @param sig Thermal tolerance width of a species
#' @param B Desired total integrated area of positive growth
#' @param lam Skewness in thermal tolerance
#' @param eps Precision of estimate
#' @param len Length of temperature vector. Higher values are more precise
#'
#' @return Scaling parameter for reproductive output.
#' @export
rAdjust<-function(sig,B,lam=-2.7,eps=1e-06,len=2^13){
  # This function creates a constant to adjust the reproduction rate so that the area under the curve is roughly equal for all species
  # sig: Thermal tolerance width of a species
  # B:   Desired total integrated area of positive growth
  # lam: Skewness in thermal tolerance
  # eps: Precision of estimate
  # len: Length of temperature vector. Higher values are more precise

  # Set up an extended version of a linear tempereature gradient
  temp <- seq(-100,100,length=len)

  # The actual optimal temperature is not important here, so we use the center of the temperature gradient
  z <- 20
  r <- exp(-(temp-z)^2/sig^2)*(1+pracma::erf(lam*(temp-z)/sig))-1

  bL <- -125; bH <- 125

  # Binary search for a value of ro such that exp(ro*r) integrates to B over all temperature values where exp(ro*r) is positive

  for(i in 1:500){
    bM <- (bL+bH)/2
    R <- exp(bM*r)
    G <- pracma::trapz(temp,(R-1)*(R>1))
    differ<- G-B
    if(abs(differ)<eps){
      break
    } else{
      if(differ>0){
        bH<-bM
      } else{
        bL<-bM
      }
    }
  }
  return(bM)
}

#' Adjust Thermal Tolerance Optimum
#'
#' The reproduction function in Urban et al. 2012 is useful for creating the shape of the reproduction rate over temperature.
#' However, the zi "optimal temperature" doesn't end up where we might expect it to be.
#' This function adjusts so that argmax of ri over temperature is  zi.
#' @param sig The thermal tolerance width of a species.
#' @param zo Optimal temperature of species.
#' @param lam Skewness in thermal tolerance.
#' @param L Length of temperature vector. Higher values are more precise.
#'
#' @return Adjusted thermal optimum
#' @export
zAdjust<-function(sig,zo,lam=-2.7,L=2^13){
  # Set up an extended version of a linear tempereature gradient
  temp <- seq(-100,100,length=L)

  # We need to calculate the difference between the expected optimal temperature and the actual optimal temperature
  # To do so, we begin with a baseline at zc=20
  zc <- 20

  # Calculate the baseline reproductive rate
  r<-exp(-(temp-zc)^2/sig^2)*(1+pracma::erf(lam*(temp-zc)/sig))-1

  # index for which temperature has the maximum reproductive output with the baseline
  iZ<-which.max(r)
  # index for baseline optimal temperature
  oZ<-which.min(abs(temp-zc))
  # index for desired optimal temperature
  tZ<-which.min(abs(temp-zo))

  # adjusted z to make optimal temperature in the right place
  z<-temp[tZ+oZ-iZ]

  return(z)
}

#' Convert Binned Populations into Individual Location Vector
#'
#' Convert vector of population sizes over \code{X} into a vector of the location for each individual.
#' @param v A vector of population sizes
#' @return The locations of each individual
#' @examples
#' unbin(c(1,3,10,5,2))
unbin <- function(v){
  L<-sum(v)
  ub<-matrix(0,L)
  j<-0
  for(i in which(v>0)){
    ub[(j+1):(j+v[i])]<-i
    j<-j+v[i]
  }
  return(c(ub))
}

#' Convert Individual Location Vector into Binned Populations
#'
#' Converts a vector of individual locations into a vector of population sizes over \code{X}.
#' @param ub The locations of each individual
#' @param L The length of the desired space vector
#' @return A vector of population sizes
#' @examples
#' unbin(c(1,3,10,5,2))
rebin <- function(ub,L){
  v<-1:L
  for(i in 1:L){
    v[i]<-length(which(ub==i))
  }
  return(v)
}


#' Caluclate Diversity Index
#'
#' Calculates the diversity of a community.
#' @param n \code{S}x\code{L}x\code{W} array of population sizes over species and space
#' @param type What scale of diversity ("alpha" or "gamma")
#' @param index Which diversity index to use ("invSimp" for inverse Simpson's; "giniSimp" for gini-Simpson's)
#'
#' @return The calculated diversity
#' @export
diversityIndex <- function(n,type="alpha",index="invSimp"){
  n<-flatComm(n)
  if(type=='alpha'){
    p <- t(n)/colSums(n)
    p2 <- rowSums(p^2)
    p2[is.nan(p2)] <- 1
    if(index=="invSimp"){
      D <- mean(1/p2)
    } else if(index=="giniSimp"){
      D <- 1- mean(p2)
    }
  } else if (type=='gamma'){
    p <- rowSums(n)/sum(n)
    D <- 1/sum(p^2)
    if(index=="invSimp"){
      D <- 1/sum(p^2)
    } else if(index=="giniSimp"){
      D <- 1- sum(p^2)
    }
  }
  return(D)
}

#' Combine Community into 1-D Space
#'
#' Take a community with subpopulations and flatten it into a single patch each.
#' @param n \code{S}x\code{L}x\code{W} array of population sizes over species, space, and subpatches
#'
#' @return \code{S}x\code{L} array of population sizes over species and space
#' @export
flatComm <- function(n){
  nf <- apply(n,c(1,2),sum)
  return(nf)
}