R/results.R

In mptools: RAMAS Metapop Tools

#' Extract simulation results from RAMAS Metapop .mp files
#' 
#' Extract population size simulation results (mean, sd, min and max), including 
#' expected minimum abundance (EMA) and its standard deviation, from a RAMAS 
#' Metapop .mp file.
#' 
#' @param mp A character string containing the path to a RAMAS Metapop .mp file 
#'   containing simulation results. Metapop .mp files are plain text files that
#'   store settings describing RAMAS metapopulation models, and the results of 
#'   simulating population dynamics according to those models. 
#' @return A \code{list} containing: \item{version}{The version of RAMAS Metapop 
#'   from which the file indicated by \code{mp} originated (if
#'   identifiable).}\item{results}{An array containing simulation results
#'   extracted from \code{file}. The number of rows is equal to the number of
#'   time steps in the simulation. The array has four columns, containing mean,
#'   sd, min and max of population size across iterations at each time step
#'   (i.e. each row), and the number of array slices is equal to the number of
#'   populations. The third dimension is named according to population names
#'   (numeric component only).} \item{iter_min}{A sorted vector of minimum
#'   abundance, across time steps, for each iteration.} \item{iter_max}{A sorted
#'   vector of maximum abundance, across time steps, for each iteration.}
#'   \item{iter_terminal}{A sorted vector of terminal abundance for each
#'   iteration.} \item{qe_thr}{The quasi-extinction threshold. When the total
#'   abundance is beneath \code{qe_thr}, the metapopulation is considered to be
#'   quasi-extinct.} \item{qe_prob}{The probability and cumulative probability
#'   of exceeding the quasi-extinction threshold (\code{qe_thr}) at each time
#'   step.} \item{EMA}{The mean minimum abundance (i.e. the mean, across
#'   iterations, of the minimum abundance for each simulation trajectory).}
#'   \item{SDMA}{The standard deviation of minimum abundance (i.e. the sd,
#'   across iterations, of the minimum abundance for each simulation
#'   trajectory).} \item{timestamp}{A POSIXlt object representing the date and
#'   time at which the simulation was completed.} \item{n_pops}{The number of
#'   populations in the simulation.} \item{duration}{The number of time steps in
#'   the simulation} \item{n_iters}{The number of iterations performed.}
#' @seealso \code{\link{meta}}
#' @note \code{mptools} has been tested with outputs generated by RAMAS Metapop
#' version 5, and may produce unexpected results for other versions. A warning
#' is issued if the user attempts to access files originating from other 
#' versions of RAMAS Metapop.
#' @references \itemize{
#'   \item{Akcakaya, H. R., Burgman, M. A., Kindvall, O., Wood, C. C., 
#'   Sjogren-Gulve, P., Hatfield, J. S., & McCarthy, M. A. (2004). \emph{Species
#'   Conservation and Management: Case Studies}. New York: Oxford University
#'   Press.}
#'   \item{\href{https://www.ramas.com/ramas.htm#metapop}{RAMAS Software}}
#'   \item{\href{https://www.ramas.com/gis-faq.htm}{RAMAS GIS and RAMAS Metapop
#'   Frequently Asked Questions}}
#' }
#' @importFrom utils read.csv read.table
#' @importFrom stats sd
#' @export
#' @examples 
#' mp <- system.file('example.mp', package='mptools')
#' res <- results(mp)
#' str(res)
#' 
#' # look at the simulation results for the first array slice (NB: this slice is
#' # all pops combined):
#' res$results[,, 1]
#' # equivalently, subset by name:
#' res$results[,, 'ALL']
#' res$results[,, 'Pop 190']
#' res$results[,, '240A24']
#' dimnames(res$results)[[3]] # population names
#' 
#' # return a matrix of mean population sizes, where columns represent
#' # populations and rows are time steps:
#' res$results[, 1, ] # or res$results[, 'mean', ]
#' 
#' # sd across iterations:
#' res$results[, 2, ] # or res$results[, 'sd', ]
#' 
#' # min pop sizes across iterations:
#' res$results[, 3, ] # or res$results[, 'min', ]
#' 
#' # max pop sizes across iterations:
#' res$results[, 4, ] # or res$results[, 'max', ] 
results <- function(mp) {
  message('Extracting simulation results from file:\n', mp)
  metapop <- check_mp(mp)
  ver <- gsub('.*\\(|\\).*', '', metapop[1])
  if(ver==metapop[1]) ver <- 'unknown'
  title <- metapop[2]
  comment <- paste(metapop[3:6], collapse='\n')
  metapop <- metapop[-(1:6)]
  if(!any(grepl('Simulation results', metapop))) {
    stop(sprintf('No simulation results found in %s.', mp), call.=FALSE)
  }
  pops <- metapop[39:(grep('^Migration$', metapop) - 1)]
  pop.names <- utils::read.csv(
    text=pops, stringsAsFactors=FALSE, header=FALSE)[, 1]
  sim.res <- metapop[grep(
    '^Simulation results', metapop):(grep('^Occupancy', metapop)-1)]
  n_iters <- as.numeric(sub('(\\d*).*', '\\1', sim.res[2]))
  res <- strsplit(sim.res[-(1:3)], ' ')
  pop.ind <- grep('Pop', res)  
  res <- do.call(rbind, res[-pop.ind])
  res <- apply(res, 2, as.numeric)
  res.t <- t(res)
  dim(res.t) <- c(4, diff(pop.ind)[1]-1, nrow(res)/(diff(pop.ind)[1]-1))
  res <- aperm(res.t, c(2, 1, 3))
  dimnames(res) <- list(NULL, 
                        c('mean', 'sd', 'min', 'max'),
                        c('ALL', pop.names))
  n_pops <- dim(res)[3] - 1
  n_steps <- nrow(res)
  minmaxterm <- metapop[grep('^Min.  Max.  Ter.$', 
                             metapop):(grep('Time to cross', metapop)-1)]
  minmaxterm <- strsplit(minmaxterm[-1], ' ')
  minmaxterm <- apply(do.call(rbind, minmaxterm), 2, as.numeric)
  if(length(minmaxterm)==3) dim(minmaxterm) <- c(1, 3)
  colnames(minmaxterm) <- c('min', 'max', 'terminal')
  EMA <- mean(minmaxterm[, 'min'])
  SDMA <- stats::sd(minmaxterm[, 'min'])
  n_manage <- as.numeric(gsub('\\D', '', metapop[grep('pop mgmnt', metapop)]))
  qext_thr <- as.numeric(metapop[grep('pop mgmnt', metapop) + n_manage + 1])
  expl_thr <- as.numeric(metapop[grep('pop mgmnt', metapop) + n_manage + 2])
  
  cross <- utils::read.table(text=metapop[
    grep('Time to cross', metapop) + 
      ceiling(seq_len(n_steps)/as.numeric(
        metapop[grep('pop mgmnt', metapop) + n_manage + 3]))],
    col.names=c('quasi_extinction', 'explosion'))
  qext_prob <- cross$quasi_extinction/n_iters
  qext_cumprob <- cumsum(qext_prob)
  out <- list(
    version=ver,
    title=title,
    comment=comment,
    results=res, 
    iter_min=minmaxterm[, 1],
    iter_max=minmaxterm[, 2],
    iter_terminal=minmaxterm[, 3],
    qe_thr=qext_thr,
    qe_prob=data.frame(
      time=seq_len(n_steps),
      prob=qext_prob, cumulative_prob=qext_cumprob),
    EMA=EMA, SDMA=SDMA,
    timestamp=sub('Simulation results ', '', 
                  grep('^Simulation results', metapop, value=TRUE)), 
    n_pops=n_pops, duration=n_steps, n_iters=n_iters)
  class(out) <- 'mpresults'
  out
}

#' Print a metapop results object
#'
#'Print summary information about a RAMAS Metapop simulation.
#'
#' @param x An mpresults object
#' @param ... ignored
#' @details The print method for \code{mpresults} objects shows only summary
#'   information. To see the full structure of \code{mpresults} object \code{x},
#'   see \code{str(x)}. Individual elements can be accessed as for a normal 
#'   \code{list}.
#' @seealso \code{\link{results}}
#' @method print mpresults
#' @export
print.mpresults <- function(x, ...) {
  stamp <- sprintf('Simulation completed at %s', x$timestamp)
  n <- sprintf('%s populations; %s time steps; %s iterations',
               x$n_pops, x$duration, x$n_iters) 
  x$title <- gsub('^\\s+|\\s+$|(?<=\\s)\\s+', '', x$title, perl=TRUE)
  if(x$title=='') x$title <- 'Untitled model'
  x$comment <- gsub('^\\s+|\\s+$|(?<=\\s)\\s+', '', x$comment, perl=TRUE)
  cat('results:\n---\n')
  pre <- paste(strwrap(c(x$title, strsplit(x$comment, '\n')[[1]]), width=60), 
               collapse='\n')
  cat(pre, '\n', '---', '\n', sep='')
  cat(stamp, '\n', n, '\n\n', sep='')
  print.default(
    c(EMA=x$EMA, `Extinction risk`=mean(x$iter_min==0)),
    print.gap=2L)
  invisible(x)
}