R/ca.R

#' Run cellular automata simulation.
#' 
#' @param x A landscape object.
#' @param model A valid object of class 'ca_model'. Defaults to musselbed. Valid
#'   values are: \code{musselbed}, \code{grazing}.
#' @param parms A list of parameters with one or several parameters containing a
#'   vector of parameter values. Those will be combined full-factorially using 
#'   \code{expand.grid()}. Will be checked against template parameters in 
#'   \code{model}. If not provided, template parameters will be used.
#' @param t_max Maximal number of timesteps. Model will be terminated even if 
#'   still in transient dynamics (i.e. model did not reach steady state)
#' @param t_eval Timespan of moving window that is evaluated for the end of 
#'   transient dynamics.
#' @param saveeach Timespan between timesteps at which a full snapshot of the 
#'   landscape is saved into the output of the simulation.
#' @param stopifsteady Binary parameter, defaults to FALSE. If TRUE, the
#'   function provided in parameter \code{steady} will be applied to test in
#'   each timestep if steady state is reached.
#' @param steady A function returning TRUE or FALSE, taking exactly the
#'   parameters \code{i}, \code{result}, \code{steadyparms}. By default the
#'   function returns TRUE if the difference in mean cover of the primary cell
#'   state (i.e. the first in the vector provided in \code{model$states}) over
#'   two subsequent timespans of length \code{steadyparms$t_eval} is smaller
#'   than \code{steadyparms$accept}.
#' @param steadyparms a list of parameters that are required by the function
#'   provided in \code{steady}.
#' @param seed An integer number serving as seed for random number generation.
#' @param plotting A binary variable. If TRUE, simulation is plotted into an animated gif.
#' @param filename A character string. Filename of animated gif (defaults to "modelrun.gif") which will be placed in current working directory.
#' 
#'   
#'   If not provided global seeds of R apply.
#' @param ... Parameters handed over to update function in \code{model$update}.
#'   
#' @return The output is returned as a list object of class \code{ca_result}, 
#'   containing a full timeseries of global and local cover as well as snapshots
#'   of the landscape.
#'   
#'   \describe{ 
#'   \item{\code{$model}}{The entire model object
#'   used to generate this simulation run, including the parameters at
#'   \code{$model$parms}}
#'   \item{\code{$time}}{Vector of timesteps.} 
#'   \item{\code{$issteady}}{Binary vector reporting for each timestep if the criterion for steady state was fulfilled. The criterion can be customized by adjusting parameter \code{steady}.} 
#'   \item{\code{$cover}}{A list of cover timeseries for each state of the model.} 
#'   \item{\code{$local}}{A list of local cover timeseries for each
#'   state of the model.} 
#'   \item{\code{$snaps}}{A vector of indices of saved snapshots.} 
#'   \item{\code{$landscapes}}{A list of landscape objects at each
#'   point in \code{$snaps}} \item{\code{$issteady}}{A binary vector of the
#'   returned values of function \code{steady}} 
#'   }
#'   
#' @details Runs iterations of the update function \code{model$update()} on the 
#'   initial landscape \code{x} until a \code{t_max} is reached. The function
#'   saves the full timeseries, i.e. a value for each timestep, of the global
#'   cover of each state as well as the average local cover of each state. The
#'   full landscape object of each timestep is stored in a list
#'   \code{result$landscapes} of the output object, but frequency of these
#'   snapshots can be altered by increasing the parameter \code{saveeach}.
#'   
#' @export
#' 
#' @examples 
#' 
#' # 1. run simulation and save a snapshot each 50 timesteps. plot timeseries and snapshots.
#' 
#  # create initial landscape
#' l <- init_landscape(c("+","0","-"), c(0.6,0.2,0.2), width = 100) 
#' p <- list(r = 0.4, d = 0.9, delta = 0.01)   # set parameters
#' r <- ca(l, model = musselbed, parms = p, t_max = 200)    # run simulation
#' plot(r)
#' 
#' par(mfrow= c(2,3))
#' sapply(c(0,25,50,100,150,200)+1, function(i) plot(r$landscapes[[i]]) )
#' 
#' # 2. run simulation and save full landsape at each timestep. create animated gif.
#' 
#' l <- init_landscape(c("1","0"), c(0.6,0.4), 100)
#' r <- ca(l, model = life, t_max = 400)
#' animate(r, "life01.gif")


ca <- function(x, model = grazing, parms = "default", 
               t_max = 200, saveeach = 1, 
               stopifsteady = FALSE, 
               steady = caspr::steady, 
               steadyparms = list(t_eval = 200, accept = 0.001),
               t_min = 0, # min number of iterations
               plotting = FALSE,
               filename = "modelrun",
               seed = NULL, ... )  {
  
  # checking for valid input
  ## parms
  
  if(class(parms) == "parms_timeseries") {
    parms_temp <- as.list(parms[1,]) 
  } else {
    parms_temp <- parms   # set temporary parameter object
  }
  
  if(parms_temp[1] != "default") { 
     if(!all(names(model$parms) %in% names(parms_temp))) {
        stop(paste("missing parameter(s): ", names(parms_temp)[all(names(model$parms) %in% names(parms_temp))] , "! please specify!"    ) )
    }
  } else {
    parms <- model$parms
    warning("you did not specify 'parms'! using default parameters of model!")
  }
  if(saveeach > t_max) { 
    warning(paste0("Parameter saveeach is larger than t_max: I will only save initial and final landscape! Increase t_max!") ) 
    saveeach <- t_max
  }
  #if(stopifsteady && identical(steadyparms$dim, x$dim)  ) {
  #  warning(paste0("steadyness criterion was defined for lattice of dimensions: ",steadyparms$dim[1], "x",steadyparms$dim[1], " and might be invalid! Please adjust parameter 'steadyparms'!" ))
  #}
  ## states
  if(!all(x$states %in% model$states)) stop("invalid cell states specified in landscape object!")
 #notworking!  if(!identical(model$states, x$states)) warning("one or several of the model's default cell states are not present in initial landscape object!") 
  ## is interaction matrix provided? if not defaults to four cell neighborhood. 
  if(is.null(model$interact)) {
    I <- matrix(c(0, 1, 0,  1, NA, 1, 0, 1, 0), ncol = 3, byrow = TRUE) 
  } else {
    I <- model$interact 
  }
  
    
  mapping(x$dim[1], x$dim[2], i_matrix = I )
  xstats <- summary(x)
  states <- levels(x$cells)
  
  # calculate timesteps for snapshots:
  snaps <- seq(0, as.integer(t_max), as.integer(saveeach))
  
  # initialise result object:
  
  result <- list()  # generate list object
  result$model <- model
  result$model$parms <- parms
  result$time <- seq(0, t_max) # add vector of realized timesteps
  #result$evaluate <- c(t_min, t_min+2*t_eval)+1
  result$issteady <- logical(length = t_max + 1 )
  #result$steadyval <- numeric( length = t_max + 1 )

  result$cover <- as.data.frame(t(xstats$cover))
  result$cover <- result$cover[rep(1, t_max+1),] # preallocating memory
  
  result$local <- as.data.frame(t(xstats$local))
  result$local <- result$local[rep(1, t_max+1),] # preallocating memory
  
  result$seed <- seed       # save seed 
  result$ini_landscape <- x # save initial landscape object at t=0
  
  result$snaps <- snaps
  result$landscapes <- list()
  result$landscapes <- lapply(1:length(snaps), function(i) x) # preallocating memory
  
  # --------------- start simulation -----------------
  

  # initialise simulation variables: 
  x_old <- x  # ghost matrix at t_i
  x_new <- x_old
  
  if(0 %in% snaps) {  
    result$landscapes[[match(0, snaps)]] <- x_new
  }
  
  i = 1  # iterator for simulation timesteps ([Alex] this should really start 
         #   at 1 but I'm afraid of breaking something. The way it is now 
         #   overwrites the initial state).
  if(!is.null(seed)) set.seed(seed)  # get seed from function call
 
  # starting iterations:
  while ( i <= t_min || 
          i <= t_max && 
          ( !stopifsteady || !steady(i, result, steadyparms) ) ) {
    
    # call update function:
    if(class(parms) == "parms_timeseries") parms_temp <- as.list(parms[i,]) 
    
    model$update(x_old, parms_temp, ...) -> x_new
    
    # save stats of new landscape
    
    xstats <- summary(x_new)
    result$cover[i+1,] <- xstats$cover
    result$local[i+1,] <- xstats$local

    # save landscape if snapshot
    
    if(i %in% snaps) {  
      result$landscapes[[match(i, snaps)]] <- x_new
    }

    result$issteady[i+1] <- steady(i, result, steadyparms)

    # replace ghost matrix for next iteration
    
    x_old <- x_new 
    i = i+1
    
  } 
  # ------------ end simulation -----------------
  
  #result$steady_state$issteady <- steadiness <= steady
  #result$steady_state$steadiness <- steadiness 
  class(result) <- "ca_result"
  
  if(plotting) ca_animate(result, filename = filename)
  return(result)
}



#' @export

print.ca_result <- function(x) {
  
  obj <- deparse(substitute(x))
  
  cat(" \n ")
  cat("Model run of ", x$model$name, " over ", max(x$time),
      " timesteps. \n", sep = "")
  
  #if(x$steady_state$issteady == TRUE) {
  # cat("\n average cover reached steady state ('steadiness' <= ", 
  #    x$steady_state$steady ,"): \n", sep = "")
  #} else {
  #  cat("\n average cover did not reach steady state ('steadiness' = ", 
  #      round(x$steady_state$steadiness, 6) ,"): \n", sep = "")
  #}
  cat("\n average cover: \n")
  cat( "  ", formatC(names(summary(x)$mean_cover), width = 6, flag = " " ) , "\n")
  cat( "  ", formatC(round(summary(x)$mean_cover, digits = 3), width = 6, flag = " " ) )
  cat(" \n \n")
  cat(" for more details call 'summary(", obj,")' or 'plot(", obj,")'! \n", sep = "")
  
  cat( " access simulation results: \n")
  cat( "   '", obj,"$model$parms' : simulation parameters \n", sep = "")
  cat( "   '", obj,"$time' : a vector of timesteps \n", sep = "")
  cat( "   '", obj,"$cover' : a dataframe of the states' timeseries \n", sep = "")
  cat( "   '", obj,"$ini_landscape' : the initial landscape object \n", sep = "")
  cat( "   '", obj,"$snaps' : an index table of saved snapshots \n", sep = "")
  cat( "   '", obj,"$landscapes[[i]]' : extract snapshot from list of snapshots \n", sep = "")
  cat(" \n ")
  
}

#' plot method for 'ca_result' objects
#' 
#' @param x An output object obtained from running function \code{ca()}.
#' @param plotstates A logical vector of the same length as 
#'   \code{x$model$states}. TRUE activates plotting of other cell state populations. 
#' @param cols A vector of valid colors to replace the default color vector.
#' @param lwd Line width for timeseries plot.
#' @param ... Parameters handed to standard plot function, e.g. `bty`, `xlab`.
#'   See \link{par}.
#' @details This prompts a (series of) summary plot for the simulation run. By
#' default this is a timeseries of the primary cell state (i.e. the first state
#' given in \code{x$model$states}).
#' 
#' @export

plot.ca_result <- function(x, plotstates = c(TRUE, rep(FALSE, length(x$model$states)-1)), snapshots = FALSE, cols = x$model$cols , lwd = 1, ...) {
  
  if(snapshots) {
    
  }
  
  plot(NA,NA,  
       type = "l", 
       col = x$model$cols[1], 
       xlab = "time", ylab = paste("cover of", x$model$states[1]) , 
       xlim = c(1, max(x$time)), ylim = c(0,1), 
       ...)
  
  if(length(plotstates) == length(x$model$states)  ) {
    for(i in (1:length(plotstates))[plotstates]) {
      lines(x$time,x$cover[[i]], col = cols[i], lwd = lwd)
      
    }
    
  }
  
}



#' summary method for `ca_result`
#'
#' @param x 
#' 
#' @param length.eval Number of time steps at the end of the simulation to consider
#'   when computing mean and sd of covers
#' 
#' @return Returns a list \code{out} containing the model name, the final time,
#'   the mean cover  and standard deviation in the last 10 timesteps.
#'   
#' @note The mean and sd are supposed to depend on the model being in steady
#'   state dynamics, i.e. beyond transitory. Lacking a universal criterion for
#'   that, we report the final 10 timesteps.
#'   
#' @export

summary.ca_result <- function(x, length.eval = 10) {
  out <- list()
  class(out) <- "ca_summary"
  eval <- tail(seq_along(x$time), length.eval)
  out$name <- x$model$name
  out$time <- c(min(x$time), max(x$time))
  out$mean_cover <- colMeans(x$cover[eval,])
  out$sd_cover <- sapply(x$cover[eval,], sd)
  
  return(out)
}


#' Transfer output of ca() into a list of matrices. 
#'
#' @param x 
#'
#' @return a list of matrices.
#' @export
#'

as.list.ca_result <- function(x) {
  lapply(x$landscapes, as.matrix)
}


#' Transfer output of runca() into an array. 
#'
#' @param x 
#'
#' @return an array with dimensions 'width', 'height' and 'snaps'.
#' @export
#'

as.array.ca_result <- function(x) {
  width <- x$landscapes[[1]]$dim[1]
  height <- x$landscapes[[1]]$dim[2]
  snaps <- length(x$landscapes)
  array(unlist(lapply(x$landscapes, as.matrix)), dim = c(width, height, snaps), dimnames = c("width", "height", "snaps"))
}
fdschneider/caspr documentation built on May 16, 2019, 12:12 p.m.