https://github.com/cran/spatstat
Raw File
Tip revision: 3aca716ce2576a0dab83f08052acd47afed8ee6a authored by Adrian Baddeley on 29 February 2012, 00:00:00 UTC
version 1.25-4
Tip revision: 3aca716
dgs.R
#
#
#    dgs.R
#
#    $Revision: 1.5 $	$Date: 2012/01/18 10:06:43 $
#
#    Diggle-Gates-Stibbard process
#
#
# -------------------------------------------------------------------
#	

DiggleGatesStibbard <- local({

  # .......... auxiliary functions ................
  dgsTerms <- function(X, Y, idX, idY, rho) {
    stopifnot(is.numeric(rho))
    # sort in increasing order of x coordinate
    oX <- order(X$x)
    oY <- order(Y$x)
    Xsort <- X[oX]
    Ysort <- Y[oY]
    idXsort <- idX[oX]
    idYsort <- idY[oY]
    nX <- npoints(X)
    nY <- npoints(Y)
    # call C routine
    out <- .C("Ediggatsti",
            nnsource = as.integer(nX),
            xsource  = as.double(Xsort$x),
            ysource  = as.double(Xsort$y),
            idsource = as.integer(idXsort),
            nntarget = as.integer(nY),
            xtarget  = as.double(Ysort$x),
            ytarget  = as.double(Ysort$y),
            idtarget = as.integer(idYsort),
            rrho     = as.double(rho),
            values   = as.double(double(nX)),
            PACKAGE  = "spatstat")
    answer <- integer(nX)
    answer[oX] <- out$values
    return(answer)
  }

  # ...... template object ......................
  BlankDGS <- 
    list(
         name   = "Diggle-Gates-Stibbard process",
         creator = "DiggleGatesStibbard",
         family  = "pairwise.family",  # evaluated later
         pot    = function(d, par) {
           rho <- par$rho
           v <- log(sin((pi/2) * d/rho)^2)
           v[ d > par$rho ] <- 0
           attr(v, "IsOffset") <- TRUE
           v
         },
         par    = list(rho = NULL),  # to be filled in later
         parnames = "interaction range", 
         init   = function(self) {
           rho <- self$par$rho
           if(!is.numeric(rho) || length(rho) != 1 || rho <= 0)
             stop("interaction range rho must be a positive number")
         },
         update = NULL,       # default OK
         print = NULL,        # default OK
         interpret =  function(coeffs, self) {
           return(NULL)
         },
         valid = function(coeffs, self) {
           return(TRUE)
         },
         project = function(coeffs, self) {
           return(NULL)
         },
         irange = function(self, coeffs=NA, epsilon=0, ...) {
           rho <- self$par$rho
           return(rho)
         },
         version=NULL, # evaluated later
         # fast evaluation is available for the border correction only
         can.do.fast=function(X,correction,par) {
           return(all(correction %in% c("border", "none")))
         },
         fasteval=function(X,U,EqualPairs,pairpot,potpars,correction, ...) {
           # fast evaluator for DiggleGatesStibbard interaction
           if(!all(correction %in% c("border", "none")))
             return(NULL)
           if(spatstat.options("fasteval") == "test")
             message("Using fast eval for DiggleGatesStibbard")
           rho <- potpars$rho
           idX <- seq_len(npoints(X))
           idU <- rep(-1, npoints(U))
           idU[EqualPairs[,2]] <- EqualPairs[,1]
           v <- dgsTerms(U, X, idU, idX, rho)
           v <- matrix(v, ncol=1)
           attr(v, "IsOffset") <- TRUE
           return(v)
         }
         )
  class(BlankDGS) <- "interact"

  DiggleGatesStibbard <- function(rho) {
    instantiate.interact(BlankDGS, list(rho = rho))
  }

  DiggleGatesStibbard
})
back to top