R/makeTmatrix.R

In latticeDensity: Density Estimation and Nonparametric Regression on Irregular Regions

#' Create the transition matrix for the diffusion.
#' 
#' This function generates a transition matrix for 
#' the diffusion process on the lattice.
#' 
#' @param formLatticeOutput A formLatticeOutput object, returned 
#' by the functions formLattice or by the function editLattice.
#' @param M A smoothing parameter. It is the maximum probability 
#' that the random walk moves from the node in a single step. 
#' It is a maximum probability in the sense that this is the
#' movement probability for nodes not near a boundary.  Of course,
#' near a boundary movement will be constrained proportional to
#' how many neighbors the node has.  Thus if interior nodes have
#' eight neighbors, a node with only four neighbors will move half as 
#' often.
#' Since the number of steps k also determines smoothing,
#' M is usually left at 0.5. Note that values of M=1 or M=0 
#' can lead to pathological results. The paper of Barry and 
#' McIntyre (2011) shows the exact construction of the 
#' transition matrix.
#' @param sparse logical. If TRUE, then uses sparse matrix 
#' computations from packages spdep and spam. If FALSE, 
#' uses full matrix computations. The use of sparse matrices 
#' is almost always more efficient.
#' 
#' @return An NxN transition matrix, where N is the number of nodes.
#' @author Ronald P. Barry
#' @references 
#' Ronald P. Barry, Julie McIntyre. Estimating animal 
#' densities and home range in regions with irregular 
#' boundaries and holes: A lattice-based alternative 
#' to the kernel density estimator. Ecological Modelling 
#' 222 (2011) 1666-1672. <doi:10.1016/j.ecolmodel.2011.02.016>
#' @examples 
#' plot.new()
#' data(polygon1)
#' nodeFillingOutput <- nodeFilling(poly=polygon1, node_spacing=0.02)
#' formLatticeOutput <- formLattice(nodeFillingOutput)
#' Pointdata <- splancs::csr(polygon1,75)
#' Pointdata <- Pointdata[Pointdata[,1]<0.5,]
#' poly.area <- splancs::areapl(polygon1)
#' init_prob <- addObservations(formLatticeOutput, Pointdata)
#' T = makeTmatrix(formLatticeOutput, M = 0.5, sparse=TRUE)
#' p10 <- Tkp(T, 10, p=init_prob$init_prob)
#' head(cbind(init_prob$init_prob, p10))
#' 
#' @import utils
#' @import graphics
#' @import stats
#' @importFrom spdep nb2listw
#' @importFrom spdep nb2mat
#' @importFrom spatialreg as.spam.listw
#' @export
makeTmatrix <- 
function(formLatticeOutput,M = 0.5, sparse = TRUE){
  # 
  M <- as.numeric(M)
  if(class(formLatticeOutput)!="formLatticeOutput"){
       stop("Should be the output from the function formLattice")}
  latt <- formLatticeOutput$latt
  #
  if(sparse){
      #  sparse matrix computations
    neigh_matrix <- spatialreg::as.spam.listw(spdep::nb2listw(
                            latt,style="B",zero.policy=TRUE))
    rowTotals <- apply(neigh_matrix,MARGIN=1,sum)
    diags <- 1 - M*(rowTotals/max(rowTotals))
    T <- spam::diag(as.vector(diags))+ M*(neigh_matrix/max(rowTotals))
    } else {
    neigh_matrix <- spdep::nb2mat(latt,style="B",zero.policy=TRUE)
    rowTotals <- rowSums(neigh_matrix)
    diags <- 1 - M*(rowTotals/max(rowTotals))
    T <- diag(as.vector(diags))+ M*(neigh_matrix/max(rowTotals))
    }
  return(T)
}