R/localTrapCalculations.R

In nimbleSCR: Spatial Capture-Recapture (SCR) Methods Using 'nimble'

#' Local Trap Calculations
#'
#' These functions are deprecated, and they will be removed from a future release.  Utility functions to enable local trap calculations in SCR models.  See details section for more information.
#'
#' @details
#'
#' These functions are deprecated, and they will be removed from a future release.
#' 
#' The makeGrid function is used in advance of model building.  It creates and returns a list of two objects: a table (grid) corresponding to the discretized grid, where each row gives the x-coordinate, the y-coordinate, and the id number for a grid cell; and second, a function (makeID) to be used in the model code which operates on a discretized AC location, and returns the id number of the corresponding grid cell.
#' 
#' The findLocalTraps function operates on the grid object returned from makeGrid, and an array of the trap location coordinates, and the desired maximal exposure radius for caluclations (dmax).  It returns a array (localTraps) with number of rows equal to the number of grid cells.  The first element of each row gives the number of local traps within exposure radius to that grid cell.  The following elements of each row give the id numbers of those local traps.
#'
#' A visualization function (plotTraps) is also provided in the example code, which displaces the discretized grid (small black points), all trap locations (green circles), a specified grid cell location (specified by i) as a large X, and the local traps to that specified grid cell (red circles).
#'
#' The getNumLocalTraps function is used inside the model code.  It operates on an id for a grid cell, the localTraps array (generated by findLocalTraps), and the constant value LTD1.  This function returns the number of traps which are local to a specified grid cell.
#'
#' The getLocalTrapIndices function is used inside the model code.  It returns a vector containing the ids of the local traps to a particular grid cell.
#'
#' The calcLocalTrapDists function is used inside the model code.  It calculates the distances from an activity center, to the local traps relative to the grid cell nearest that activity center.
#'
#' The calcLocalTrapExposure function is specific to the detection probability calculations used in this example.  This function should be modified specifically to the detection function, exposure function, or otherwise calculations to be done only for the traps in the vicinity of individual activity center locations
#'
#' @name localTrapCalculations
#'
#' @param xmin Minimal value among all trap location x-coordinates.
#' @param xmax Maximal value among all trap location x-coordinates.
#' @param ymin Minimal value among all trap location y-coordinates.
#' @param ymax Maximal value among all trap location y-coordinates.
#' @param resolution Desired resolution (in both x and y directions) of discretized grid.
#' @param buffer Horizontal and vertical buffer for discretized grid, specifying how much it should extend (above, below, left, and right) of the maximal trap locations.
#' @param grid The grid object returned from the makeGrid function.
#' @param trapCoords An nTraps x 2 array giving giving the x- and y-coordinate locations of all traps.
#' @param dmax The maximal radius from an activity center for performing trap calculations (dmax).
#' @param idarg A grid id, returned from the makeID function inside model code.
#' @param nLocalTraps The number of local traps to all grid cells, which is given by the first column of the localTraps array.
#' @param LTD1arg The number of columns in the localTraps array.
#' @param MAXNUM The maximum number of local traps among all grid cells.  This is given by the (number of rows)-1 of the localTraps array.
#' @param localTraps The array returned from the findLocalTraps function.
#' @param n The number of local traps to a specified grid cell, as return.
#' @param localTrapInd The indices of the local traps to a grid cell, as returned by the getLocalTrapIndices function.
#' @param s A length-2 vector giving the activity center of an indiviual.
#' @param R The total number of traps.
#' @param d A vector of distances from an activity center to the local traps.
#' @param sigma Scale of decay for detection probability.
#' @param p0 Baseline detection probability.
#'
#' @author Daniel Turek
#'
#' @import nimble
#'
#' @examples
#'
#' \dontrun{
#' 
#' ## generate random trap locations
#' nTraps <- 200
#' traps_xmin <- 0
#' traps_ymin <- 0
#' traps_xmax <- 100
#' traps_ymax <- 200
#' set.seed(0)
#' traps_xCoords <- round(runif(nTraps, traps_xmin, traps_xmax))
#' traps_yCoords <- round(runif(nTraps, traps_ymin, traps_ymax))
#' trap_coords <- cbind(traps_xCoords, traps_yCoords)
#'  
#' ## buffer distance surrounding sides of rectangular discretization grid
#' ## which overlays trap locations
#' buffer <- 10
#'  
#' ## resolution of rectangular discretization grid
#' resolution <- 10
#'  
#' ## creates grid and makeID function,
#' ## for grid overlaying trap locations,
#' ## and to lookup nearest grid cell to any AC
#' makeGridReturn <- makeGrid(xmin = traps_xmin, xmax = traps_xmax,
#'                            ymin = traps_ymin, ymax = traps_ymax,
#'                            buffer = buffer,
#'                            resolution = resolution)
#'  
#' grid <- makeGridReturn$grid
#' makeID <- makeGridReturn$makeID
#'  
#' ## maximum radis within an individual AC to perform trap calculations,
#' dmax <- 30
#'  
#' ## n = localTraps[i,1] gives the number of local traps
#' ## localTraps[i, 2:(n+1)] gives the indices of the local traps
#' localTraps <- findLocalTraps(grid, trap_coords, dmax)
#'  
#' plotTraps <- function(i, grid, trap_coords, localTraps) {
#'     plot(grid[,1], grid[,2], pch = '.', cex=2)
#'     points(trap_coords[,1], trap_coords[,2], pch=20, col='forestgreen', cex=1)
#'     if(!missing(i)) {
#'         i <- max(i %% dim(grid)[1], 1)
#'         n <- localTraps[i,1]
#'         trapInd <- numeric(0)
#'         if(n > 0)  trapInd <- localTraps[i,2:(n+1)]
#'         theseTraps <- trap_coords[trapInd,, drop = FALSE]
#'         points(theseTraps[,1], theseTraps[,2], pch = 20, col = 'red', cex=1.5)
#'         points(grid[i,1], grid[i,2], pch = 'x', col = 'blue', cex=3)
#'     }
#' }
#'  
#' ## visualise some local traps
#' plotTraps(10,  grid, trap_coords, localTraps)
#' plotTraps(200, grid, trap_coords, localTraps)
#' plotTraps(380, grid, trap_coords, localTraps)
#'  
#' ## example model code
#' ## using local trap calculations
#' code <- nimbleCode({
#'     sigma ~ dunif(0, 100)
#'     p0 ~ dunif(0, 1)
#'     for(i in 1:N) {
#'         S[i,1] ~ dunif(0, xmax)
#'         S[i,2] ~ dunif(0, ymax)
#'         Sdiscrete[i,1] <- round(S[i,1]/res) * res
#'         Sdiscrete[i,2] <- round(S[i,2]/res) * res
#'         id[i] <- makeID( Sdiscrete[i,1:2] )
#'         nLocalTraps[i] <- getNumLocalTraps(id[i], localTraps[1:LTD1,1], LTD1)
#'         localTrapIndices[i,1:maxTraps] <-
#'             getLocalTrapIndices(maxTraps, localTraps[1:LTD1,1:LTD2], nLocalTraps[i], id[i])
#'         d[i, 1:maxTraps] <- calcLocalTrapDists(
#'             maxTraps, nLocalTraps[i], localTrapIndices[i,1:maxTraps],
#'             S[i,1:2], trap_coords[1:nTraps,1:2])
#'         g[i, 1:nTraps] <- calcLocalTrapExposure(
#'             nTraps, nLocalTraps[i], d[i,1:maxTraps], localTrapIndices[i,1:maxTraps], sigma, p0)
#'         y[i, 1:nTraps] ~ dbinom_vector(prob = g[i,1:nTraps], size = trials[1:nTraps])
#'     }
#' })
#'  
#' ## generate random detection data; completely random
#' N <- 100
#' set.seed(0)
#' y <- array(rbinom(N*nTraps, size=1, prob=0.8), c(N, nTraps))
#'  
#' ## generate AC location initial values
#' Sinit <- cbind(runif(N, traps_xmin, traps_xmax),
#'                runif(N, traps_ymin, traps_ymax))
#'  
#' constants <- list(N = N,
#'                   nTraps = nTraps,
#'                   trap_coords = trap_coords,
#'                   xmax = traps_xmax,
#'                   ymax = traps_ymax,
#'                   res = resolution,
#'                   localTraps = localTraps,
#'                   LTD1 = dim(localTraps)[1],
#'                   LTD2 = dim(localTraps)[2],
#'                   maxTraps = dim(localTraps)[2] - 1)
#'  
#' data <- list(y = y, trials = rep(1,nTraps))
#'  
#' inits <- list(sigma = 1,
#'               p0 = 0.5,
#'               S = Sinit)
#'  
#' ## create NIMBLE model object
#' Rmodel <- nimbleModel(code, constants, data, inits,
#'                       calculate = FALSE, check = FALSE)
#'  
#' ## use model object for MCMC, etc.
#'
#' }
#'
NULL


#' @rdname localTrapCalculations
#' @export
makeGrid <- function( xmin = 0,
                      ymin = 0,
                      xmax,
                      ymax,
                      resolution = 1,
                      buffer = 0) {
    warning('This function is deprecated, and will be removed from a future release.')
    makeVals <- function(min, max, buf, res) {
        unique(c(rev(seq(min, min-buf, by = -res)), seq(min, max+buf, by = res)))
    }
    xvals <- makeVals(xmin, xmax, buffer, resolution)
    yvals <- makeVals(ymin, ymax, buffer, resolution)
    grid <- expand.grid(xvals, yvals)
    colnames(grid) <- c('x', 'y')
    ## unique ids:
    mult <- diff(range(grid$y/resolution)) + 1
    ids <- grid$x/resolution * mult + grid$y/resolution
    offset <- 1 - min(ids)
    makeIDdef <- substitute(
        nimbleFunction(
            run = function(xy = double(1)) {
                id <- xy[1]/RES * MULT + xy[2]/RES + OFFSET
                returnType(double())
                return(id)
            }
        ),
        list(RES = resolution,
             MULT = mult,
             OFFSET = offset))
    makeID <- eval(makeIDdef)
    ids2 <- apply(grid, 1, function(xy) makeID(xy))
    sorted <- sort(ids2, index.return = TRUE)
    gridReordered <- grid[sorted$ix, ]
    gridReordered$id <- sorted$x
    return(list(grid = gridReordered, makeID = makeID))
}

#' @rdname localTrapCalculations
#' @export
findLocalTraps <- function( grid,
                            trapCoords,
                            dmax) {
    warning('This function is deprecated, and will be removed from a future release.')
    trtrapsBool <- apply(grid, 1, function(row) {
        apply(trapCoords, 1, function(tp) {
            sqrt(sum((row[1:2]-tp)^2)) <= dmax
        })
    })
    trapsBool <- t(trtrapsBool)
    trapsInd <- apply(trapsBool, 1, which)
    numsTraps <- sapply(trapsInd, length)
    localTraps <- array(as.numeric(NA), c(dim(grid)[1], max(numsTraps)+1))
    for(i in seq_along(trapsInd)) {
        n <- numsTraps[i]
        localTraps[i,1] <- n
        if(n > 0)    localTraps[i, 2:(n+1)] <- trapsInd[[i]]
    }
    localTraps
}

#' @rdname localTrapCalculations
#' @export
getNumLocalTraps <- nimbleFunction(
    run = function( idarg = double(),
                    nLocalTraps = double(1),
                    LTD1arg = double()) {
        stop('This function is deprecated, and will be removed from a future release.')
        if(idarg < 1)       {   return(0)   }
        if(idarg > LTD1arg) {   return(0)   }
        n <- nLocalTraps[idarg]
        returnType(double())
        return(n)
    })

#' @rdname localTrapCalculations
#' @export
getLocalTrapIndices <- nimbleFunction(
    run = function( MAXNUM = double(),
                    localTraps = double(2),
                    n = double(),
                    idarg = double()) {
        stop('This function is deprecated, and will be removed from a future release.')
        indices <- numeric(MAXNUM, 0)
        if(n > 0) {
            indices[1:n] <- localTraps[idarg, 2:(n+1)]
        }
        returnType(double(1))
        return(indices)
    })

#' @rdname localTrapCalculations
#' @export
calcLocalTrapDists <- nimbleFunction(
    run = function( MAXNUM = double(),
                    n = double(),
                    localTrapInd = double(1),
                    s = double(1),
                    trapCoords = double(2)) {
        stop('This function is deprecated, and will be removed from a future release.')
        Ds <- numeric(MAXNUM, 0)
        if(n > 0) {
            Ds[1:n] <- sqrt((s[1] - trapCoords[localTrapInd[1:n],1])^2 + 
                                (s[2] - trapCoords[localTrapInd[1:n],2])^2)
        }
        returnType(double(1))
        return(Ds)
    })

#' @rdname localTrapCalculations
#' @export
calcLocalTrapExposure <- nimbleFunction(
    run = function( R = double(),
                    n = double(),
                    d = double(1),
                    localTrapInd = double(1),
                    sigma = double(),
                    p0 = double()) {
        stop('This function is deprecated, and will be removed from a future release.')
        g <- numeric(R, 0.00000000000001)      ## small value
        if(n > 0) {
            g[localTrapInd[1:n]] <- p0 * exp(-d[1:n]/sigma)
        }
        returnType(double(1))
        return(g)
    })