Revolve: Adaptive Dynamics toolkit for R

##' @title R* model
##'
##' Model implementing the "R*" model, most closely based on that
##' presented by "Huisman, J., and Weissing, F.J. (2001).  Biological
##' conditions for oscillations and chaos generated by multispecies
##' competition. Ecology 82, 2682-2695."
##'
##' @rdname rstar
##' @export rstar
##' @export
##' @importFrom methods setRefClass
rstar <- setRefClass("rstar",
                     fields=list(
                       r="numeric",
                       m="numeric",
                       D="numeric",
                       S="numeric",
                       mat="list",
                       k="integer"))

rstar$methods(initialize=function(matrices, S, r=1, m=0.25, D=0.25) {
  assert_scalar(r)
  assert_scalar(m)
  assert_scalar(D)
  if (!(length(S) %in% c(1, matrices$k))) {
    stop("S must be scalar or length k")
  }
  r <<- r
  m <<- m
  D <<- D
  S <<- rep(S, length.out=matrices$k)
  mat <<- matrices
  # Not sure about these
  k <<- mat$k
})

rstar$methods(K=function(...) mat$K(...))
rstar$methods(C=function(...) mat$C(...))
rstar$methods(p=function(x, R) {
  "Specific growth rate at a particular resource level R"
  r * R / (K(x) + R)
})
rstar$methods(pinv=function(x, g) {
  "Resource level that would be required for particular growth rate g"
  g * K(x) / (r - g)
})

rstar$methods(max_growth_rate=function(x, global=FALSE) {
  (if (global) rep(r, ncol(x)) else .self$min_p(x, S)) - m
})

rstar$methods(Rstar=function(x) {
  "Compute equilibrium level of each resource (i.e. level it will
reach when limiting.  This is the resource level where the growth rate
equals the mortality rate."
  pinv(x, m)
})

rstar$methods(carrying_capacity=function(x) {
  "Compute the carrying capacity -- the population size that a species
with trait x will reach when alone."
  colMins(D * (S - Rstar(x)) / (m * C(x)))
})

## Methods to run the model:

## First, define a bunch of derivatives:
rstar$methods(dRdt=function(x, N, R) {
  "Change in resources with respect to time"
  D * (S - R) - drop(C(x) %*% (N * .self$min_p(x, R)))
})

rstar$methods(dNdt=function(x, N, R) {
  "Change in density with respect to time (total, not per capita)"
  N * (.self$min_p(x, R) - m)
})
## I might need to change these so that they are self contained free
## functions stored as fields within the class to avoid the overhead
## of calling reference methods.
rstar$methods(derivs=function(t, ode_y, sys, ...) {
  R <- ode_y[mat$i.R]
  N <- ode_y[-mat$i.R]
  c(dRdt(sys$x, N, R),
    dNdt(sys$x, N, R))
})
rstar$methods(derivs_R=function(t, R, sys, ...) {
  dRdt(sys$x, sys$y, R)
})

## Dynamically run the model:
rstar$methods(run=function(sys, times) {
  ode_y <- .self$initial_conditions(sys$y)
  res <- lsoda_nolist(ode_y, times, .self$derivs, sys)[,-1,drop=FALSE]
  list(t=times, x=sys$x, y=res[,-mat$i.R], R=res[,mat$i.R])
})

rstar$methods(run_fixed_density=function(sys, times) {
  R0 <- .self$R0_from_sys(sys)
  res <- lsoda_nolist(R0, times, .self$derivs_R, sys)
  list(t=times, R=res[,-1])
})

## Equilibrium levels:
rstar$methods(equilibrium=function(sys, ...) {
  "Given traits and initial densities as a system (sys) compute equilibirum numerically."
  ode_y <- .self$initial_conditions(sys$y)
  res <- equilibrium_(.self$derivs, sys, ode_y, ...)
  modifyList(sys, list(y=res[-mat$i.R], R=res[mat$i.R]))
})

rstar$methods(equilibrium_R=function(sys, ...) {
  "Given traits 'x' and initial densities 'y', compute the
  equilibrium.  Return the equilibrium resource availability at the
  same time."
  R0 <- .self$R0_from_sys(sys)
  sys$R <- equilibrium_(.self$derivs_R, sys, R0, ...)
  sys
})

## Special equilibria:
rstar$methods(single_equilibrium=function(x) {
  "Given species traits \\code{x}, compute equilibrium system.  If multiple species are given, the output will be multiple separate equilibira, rather than an actual community."
  N <- carrying_capacity(x)
  R <- Rstar(x)

  extinct <- N < 0
  N[extinct] <- 0

  ## For 1 resource both R and N must be at equilibrium  
  if (nrow(R) > 1) {
    len.min <- colMins((S - R) / C(x))
    R <- S - C(x) * rep(len.min, each=nrow(R))
  }
  ## Fixup cases where species went extinct:
  R[,extinct] <- S

  sys(x, N, R=R)
})

rstar$methods(single_equilibrium_R=function(sys,
                force.numerical=FALSE) {
  "Solve for the equilibrium level of a resource, given a vector of
densities (and species traits).  Do do this, we solve dRdt == 0:
    $D (S - R) - c N p(R) == 0$
    $D * (S - R) - c N r R / (K + R) == 0$
 As a quadratic with respect to 'R':
    $-D * R^2 + (D (S - K) - c r N) * R + D K S == 0$"
  if (k == 1 && !force.numerical) {
    Kx <- K(sys$x)
    Cx <- C(sys$x)
    N <- sys$y
    ans <- quadratic_roots(-D, (D * (S - Kx) - Cx * r * N), D * Kx * S)
    ans <- ans[ans >= 0 & ans <= S]
    if (length(ans) != 1) {
      stop("Did not find a unique solution")
    }
    ans
  } else {
    equilibrium_R(sys, "nleqslv")$R
  }
})

## General fitness
##
## For fitness we want the *per-capita growth rate*; that is (min_p -
## m), rather than y * (min_p - m).  We can't just divide by y
## because we want it when y = 0 (rather than just being close).  For
## now, I think we can just fake it by repeating a bunch of the
## stuff.
##
## The issue that we have in terms of interface is that this doesn't
## really depend on {x,y} as such, but on 'R', which is itself
## determined by {x,y}.  So this function allows any R to be plugged
## in, and if none is given we compute the equilibrium resource
## density *holding density y constant*.
rstar$methods(fitness=function(x_new, x, y, R=NULL) {
  if (is.null(R)) {
    if (length(y) != 1) {
      stop("Only works on a single species at the moment")
    }
    R <- single_equilibrium_R(sys(x, y))
  }
  .self$min_p(x_new, drop(R)) - m
})

## Helper methods:
rstar$methods(initial_conditions=function(y0) {
  if (!is.na(mat$n) && length(y0) != mat$n) {
    stop("Invalid length initial conditions")
  }
  c(S, y0)
})
rstar$methods(R0_from_sys=function(sys) {
  if (is.null(sys$R)) S else sys$R
})
rstar$methods(min_p=function(x, R) {
  colMins(p(x, R))
})

## I'd like to be able to switch between different functions for K and
## C, partly so that we can test the behaviour of the model against
## results in the paper.  The simplest way of doing this would be to
## allow arbitrary functions to be passed in as K(x) and C(x).  Or if
## matrices are passed in then build simple functions around these.

## For a single resouce case, there is a single trivial non-constant
## function:

##' @export
##' @rdname rstar
##' @param k number of resources
##' @param x matrix of traits (for the \code{rstar_mat_} functions).
make_rstar_identity <- function(k) {
  ret <- function(x) {
    if (!(is.matrix(x) && nrow(x) == k)) {
      stop("x must be a matrix with ", k, " rows")
    }
    x
  }
  attr(ret, "n") <- k
  ret
}

##' @export
##' @rdname rstar
rstar_mat_1 <- make_rstar_identity(1L)
##' @export
##' @rdname rstar
rstar_mat_2 <- make_rstar_identity(2L)

##' @export
##' @rdname rstar
rstar_mat_2_tradeoff <- function(x) {
  rbind(x, 1-x, deparse.level=0)
}
attr(rstar_mat_2_tradeoff, "n") <- 1

##' @export
##' @rdname rstar
##' @param M coumn matrix with K or C matrix that will be shared by
##' all species.
make_rstar_mat_constant <- function(M) {
  if (!is.matrix(M) || ncol(M) != 1) {
    stop("M must be a single column matrix")
  }
  ret <- function(x) {
    matrix(M, nrow=nrow(M), ncol=ncol(x))
  }
  attr(ret, "n") <- 0
  ret
}

##' @export
##' @rdname rstar
##' @param K Matrix of Monod constants
##' @param C Matrix of consumption constants
rstar_matrices_fixed <- function(K, C) {
  if (!is.matrix(K) || !is.matrix(C)) {
    stop("K and C must both be matrices")
  }
  if (!identical(dim(K), dim(C))) {
    stop("K and C must have the same dimensions")
  }
  # Note that this does not use make_rstar_mat_constant because
  # these matrices are constant with respect to both the number of
  # resources *and* the number of species.  This exists only for
  # checking against the paper.
  k <- nrow(K)
  matrices <- list(K=function(x) K, C=function(x) C,
                   n=ncol(K), k=k,
                   i.K=integer(0), i.C=integer(0),
                   i.R=seq_len(k))
}


##' @export
##' @rdname rstar
rstar_matrices <- function(K, C) {
  if (is.matrix(K)) {
    rstar_matrices(make_rstar_mat_constant(K), C)
  } else if (is.matrix(C)) {
    rstar_matrices(K, make_rstar_mat_constant(C))
  } else {
    n.K <- get_attribute(K, "n")
    n.C <- get_attribute(C, "n")
    i.K <- seq_len(n.K)
    i.C <- seq_len(n.C) + n.K
    # Number of resouces, computed with dummy vector
    k <- nrow(K(matrix(rep(0.5, n.K), n.K, 1)))
    if (nrow(C(matrix(rep(0.5, n.C), n.C, 1))) != k) {
      stop("K and C disagree on k")
    }
    list(K=function(x) K(x[i.K,,drop=FALSE]),
         C=function(x) C(x[i.C,,drop=FALSE]),
         n=NA, k=k, i.K=i.K, i.C=i.C, i.R=seq_len(k))
  }
}
##' Plot Rstar ZNGIs
##'
##' Will change
##' @title Plot Rstar ZNGIs
##' @param m Rstar model, made by \code{make_rstar}
##' @param sys Model state (x is used)
##' @param xlim X axis limits (resource 1)
##' @param ylim Y axis limits (resource 2)
##' @param col Vector of colours along species
##' @author Rich FitzJohn
##' @export
rstar_plot <- function(m, sys, xlim=c(0, 1), ylim=c(0, 1),
                       col=seq_len(ncol(x))) {
  if (m$k != 2) {
    stop("Can only plot 2 resource case at the moment")
  }
  x <- sys$x
  n <- ncol(x)
  col <- rep(col, length.out=n)

  rs <- m$Rstar(x)
  plot(NA, xlim=xlim, ylim=ylim, xlab="R1", ylab="R2")
  abline(v=rs[1,], lty=3, col=col)
  abline(h=rs[2,], lty=3, col=col)
  segments(rs[1,], rs[2,], rs[1,], par("usr")[4], col=col)
  segments(rs[1,], rs[2,], par("usr")[2], rs[2,], col=col)

  R.eq <- apply(rs, 1, max)
  # This only makes sense for some cases?
  points(R.eq[1], R.eq[2], pch=19)

  C <- m$C(x)
  for (i in seq_len(n)) {
    abcline(R.eq[1], R.eq[2], C[2,i] / C[1,i], col=col[i], lty=2)
  }
}

##' @export
##' @rdname rstar_plot
##' @param eps Density at which species are considered extinct
##' @param t_max Maximum time to run simulation until
##' @param t_len Number of time intervals to run simulation over
##' @param S Resource supply rates (and initial resource levels)
##' @param col_died Colour of communities that go extinct.
rstar_trajectory <- function(m, sys,
                             col=seq_len(ncol(sys$x)), eps=1e-6,
                             t_max=100, t_len=101, S=NULL,
                             col_died="grey") {
  check_sys(sys)
  if (ncol(sys$x) > 2) {
    stop("Only working for up to two species so far")
  }
  
  t <- seq(0, t_max, length.out=t_len)

  if (!is.null(S)) {
    if (length(S) != m$k)
      stop("Invalid length S")
    oS <- m$S
    m$S <- S
    on.exit(m$S <- oS)
  } else {
    S <- m$S
  }

  eq <- m$equilibrium(sys)
  tr <- m$run(sys, t)
  survived <- eq$y > eps

  if (sum(survived) == 1) {
    col <- col[which(survived)]
  } else if (sum(survived) == 2) {
    col <- mix(col[1], col[2])
  } else {
    col <- col_died
  }

  lines(tr$R[,1], tr$R[,2], lty=2, col=col)
  points(c(S[1], eq$R[1]),
         c(S[2], eq$R[2]), pch=19, col=col, cex=.5)
}