caspr: Cellular Automata for Spatial Pressure in R

# core functions
# --------------

#' Create vector maps for cellular automata updating.
#' 
#' @param width An integer number giving the width of the landscape or a landscape
#'   object.
#' @param height An integer number. The height of the landscape. Defaults to
#'   value given in width.
#' @param boundary A character value representing the type of boundary
#'   conditions. Defaults to \code{"periodic"}. No other boundaries implemented,
#'   yet.
#' @param i_matrix A matrix defining the interaction matrix to be evaluated.
#'   Defaults to 4-cell neighborhood.
#'   
#' @return Creates a map of the landscape object to translate it into a vector
#'   with boundaries (\code{x_with_border}) and another one to back-translate it
#'   to a vector without boundaries  (\code{x_to_evaluate}) into the global
#'   environment. This will automatically be called within \code{ca()}.
#'   
#' @details This function and the objects generated by it are necessary for the
#'   function \code{\link{neighbors}} to vectorise the calculation of
#'   neighbors. The interaction matrix can be flexibly defined as a matrix that
#'   has precisely one cell with value \code{NA}, marking the focal cell. Any
#'   number larger than 0 marks a cell that is to be taken into account for
#'   counting. 
#'   
#' @export

mapping <- function(width, height = width, boundary = "periodic", i_matrix = matrix(c(0,1,0,1,NA,1,0,1,0), ncol = 3, byrow = TRUE)) {
  
  if("landscape" %in% class(width)) {
    mapping(width$dim[1], width$dim[2]) 
  } else {
    
    # derive helper vectors for counting: 
    # transformation vector for evaluation at the border of the grid
    # set evaluation matrix 
    X <- matrix(as.integer(1:(width*height)), ncol = width, byrow =TRUE)
    # setting the border of the evaluation matrix X
    X <- cbind(X[,width], X, X[,1] )  
    X <- rbind(X[height,], X, X[1,] ) 
    # transformation vector which adds the border to the lattice:
    assign("x_with_border", as.integer(t(X))  , envir = .GlobalEnv )
    
    # from the matrix X (with copied border cells), which cells are the actual cells (reverse transformation vector of the previous lines) 
    #x_to_evaluate <- sort(matrix(1:prod(dim(X)), ncol = dim(X)[2], byrow =TRUE)[-c(1, dim(X)[1]), -c(1,dim(X)[2])]  )    
    
    
    assign("x_to_evaluate", sort(matrix(1:prod(dim(X)), ncol = dim(X)[2], byrow =TRUE)[-c(1, dim(X)[1]), -c(1,dim(X)[2])]  )	, envir = .GlobalEnv )
    
    
    # defining the neighborhood which is to be evaluated	
    # set interaction matrix
    I <- i_matrix	
    # coordinates of neighbours in Interaction matrix I: 
    neighbours_in_I <- which(is.finite(abs(I)/abs(I)), arr.ind = TRUE)
    # coordinates relative to the evaluated cell (=  which(is.na(I) ) 
    relrow <- neighbours_in_I[,1]-which(is.na(I), arr.ind = TRUE)[1]
    relcol <- neighbours_in_I[,2]-which(is.na(I), arr.ind = TRUE)[2]
    
    # relative position of the four direct neighbours of a cell
    #interact <- (relrow * dim(X)[2] + relcol)
    
    assign("interact", relrow * dim(X)[2] + relcol, envir = .GlobalEnv )
  }
}


#' Count neighbors.
#' 
#' @param x A landscape object.
#' @param state A character value. The state to count. Defaults to the primary
#'   cell state, i.e. the first state given in \code{levels(x$cells)}.
#'   
#' @return Returns a vector of the counts in the neighborhood of each cell in
#'   the landscape. Neighborhood is specified by the global
#'   \code{\link{mapping}} objects, by default the 4-cell neighborhood.
#'   
#' @details This function is vectorised by using the mapping objects defined by
#'   a previous call of \code{\link{mapping}}.
#'   
#' @seealso For counting the neighborhood of \emph{one} target file, e.g. when
#'   doing sequential updating of randomly drawn cells, the function
#'   \code{\link{neighborsij}} is slightly faster. However, for such updating
#'   processes we recommend using C++ via \code{Rcpp} (see implementation of
#'   the predator prey model, \code{predprey.R}, for an example).
#'   
#' @export

neighbors <- function(x, state = levels(x$cells)[1]) {
  
  out <- integer(length = prod(x$dim))
  x_binary_with_border <- (x$cells %in% state)[x_with_border]
  for(k in interact) {
    out <- out + x_binary_with_border[x_to_evaluate+k]
  }
  return(out)  
}


#' Count neighbors of one cell.
#' 
#' @export
#' 
neighborsij <- function(x, state = levels(x$cells)[1], i = NULL, j = NULL, select = 1:prod(x$dim) ) {
  
  if(!is.null(i) & !is.null(j)) select <- as.integer(matrix(1:prod(x$dim), ncol = x$dim[1], byrow = TRUE)[j,i]) 
  
  x_binary_with_border <- (x$cells %in% state)[x_with_border]
  colSums(sapply(select, function(k) x_binary_with_border[x_to_evaluate[k]+interact] ))
  
}



# helper functions
# ----------------

#' @export
grayscale <- colorRampPalette(c("black", "white"), space = "rgb")