R/pcfmulti.R

In spatstat/spatstat.core: Core Functionality of the 'spatstat' Family

#
#   pcfmulti.R
#
#   $Revision: 1.8 $   $Date: 2016/09/21 07:28:58 $
#
#   multitype pair correlation functions
#

pcfcross <- 
  function(X, i, j, ...,
         r=NULL, kernel="epanechnikov", bw=NULL, stoyan=0.15,
         correction = c("isotropic", "Ripley", "translate"),
         divisor=c("r","d"))
{
  verifyclass(X, "ppp")
  stopifnot(is.multitype(X))
  if(missing(correction))
    correction <- NULL
  divisor <- match.arg(divisor)
  ##
  marx <- marks(X)
  if(missing(i))
    i <- levels(marx)[1]
  if(missing(j))
    j <- levels(marx)[2]
  I <- (marx == i)
  J <- (marx == j)
  Iname <- paste("points with mark i =", i)
  Jname <- paste("points with mark j =", j)
  ##
  result <- pcfmulti(X, I, J, ...,
                     r=r, 
                     kernel=kernel, bw=bw, stoyan=stoyan,
                     correction=correction,
                     divisor=divisor,
                     Iname=Iname, Jname=Jname)
  ##
  iname <- make.parseable(paste(i))
  jname <- make.parseable(paste(j))
  result <-
    rebadge.fv(result,
               substitute(g[i,j](r),
                          list(i=iname,j=jname)),
               c("g", paste0("list", paren(paste(iname, jname, sep=",")))),
               new.yexp=substitute(g[list(i,j)](r),
                                   list(i=iname,j=jname)))
  return(result)
}

pcfdot <- 
function(X, i, ...,
         r=NULL, kernel="epanechnikov", bw=NULL, stoyan=0.15,
         correction = c("isotropic", "Ripley", "translate"),
         divisor=c("r", "d"))
{
  verifyclass(X, "ppp")
  stopifnot(is.multitype(X))
  if(missing(correction))
    correction <- NULL
  divisor <- match.arg(divisor)

  marx <- marks(X)
  if(missing(i))
    i <- levels(marx)[1]

  I <- (marx == i)
  J <- rep.int(TRUE, X$n)  # i.e. all points
  Iname <- paste("points with mark i =", i)
  Jname <- "points"
	
  result <- pcfmulti(X, I, J, ...,
                     r=r, kernel=kernel, bw=bw, stoyan=stoyan,
                     correction=correction,
                     divisor=divisor,
                     Iname=Iname, Jname=Jname)

  iname <- make.parseable(paste(i))
  result <-
    rebadge.fv(result,
               substitute(g[i ~ dot](r), list(i=iname)),
               c("g", paste0(iname, "~symbol(\"\\267\")")),
               new.yexp=substitute(g[i ~ symbol("\267")](r),
                 list(i=iname)))
  return(result)
}


pcfmulti <- function(X, I, J, ...,
                     r=NULL, 
                     kernel="epanechnikov", bw=NULL, stoyan=0.15,
                     correction=c("translate", "Ripley"),
                     divisor=c("r","d"),
                     Iname="points satisfying condition I",
                     Jname="points satisfying condition J")
{
  verifyclass(X, "ppp")
#  r.override <- !is.null(r)
  divisor <- match.arg(divisor)

  win <- X$window
  areaW <- area(win)
  npts <- npoints(X)
  
  correction.given <- !missing(correction) && !is.null(correction)
  if(is.null(correction))
    correction <- c("translate", "Ripley")
  correction <- pickoption("correction", correction,
                           c(isotropic="isotropic",
                             Ripley="isotropic",
                             trans="translate",
                             translate="translate",
                             translation="translate",
                             best="best"),
                           multi=TRUE)

  correction <- implemented.for.K(correction, win$type, correction.given)
  
  ## .......... indices I and J .............................
  
  I <- ppsubset(X, I)
  J <- ppsubset(X, J)
  if(is.null(I) || is.null(J))
    stop("I and J must be valid subset indices")

  nI <- sum(I)
  nJ <- sum(J)
  if(nI == 0) stop(paste("There are no", Iname))
  if(nJ == 0) stop(paste("There are no", Jname))

  XI <- X[I]
  XJ <- X[J]

#  lambdaI <- nI/areaW
  lambdaJ <- nJ/areaW
  nIJ <- sum(I & J)
  lambdaIJarea <- (nI * nJ - nIJ)/areaW
  
  ## ...........  kernel bandwidth and support .........................
  
  if(is.null(bw) && kernel=="epanechnikov") {
    # Stoyan & Stoyan 1995, eq (15.16), page 285
    h <- stoyan /sqrt(lambdaJ)
    hmax <- h
    # conversion to standard deviation
    bw <- h/sqrt(5)
  } else if(is.numeric(bw)) {
    # standard deviation of kernel specified
    # upper bound on half-width
    hmax <- 3 * bw
  } else {
    # data-dependent bandwidth selection: guess upper bound on half-width
    hmax <- 2 * stoyan /sqrt(lambdaJ)
  }


########## r values ############################
  # handle argument r 

  rmaxdefault <- rmax.rule("K", win, lambdaJ)
  breaks <- handle.r.b.args(r, NULL, win, rmaxdefault=rmaxdefault)
  if(!(breaks$even))
    stop("r values must be evenly spaced")
  # extract r values
  r <- breaks$r
  rmax <- breaks$max
  # recommended range of r values for plotting
  alim <- c(0, min(rmax, rmaxdefault))

  # initialise fv object
  
  df <- data.frame(r=r, theo=rep.int(1,length(r)))
  fname <- c("g", "list(I,J)")
  yexp <- quote(g[list(I,J)](r))
  out <- fv(df, "r",
            quote(g[I,J](r)), "theo", ,
            alim,
            c("r", makefvlabel(NULL, NULL, fname, "Pois")),
            c("distance argument r", "theoretical Poisson %s"),
            fname=fname,
            yexp=yexp)
  
  ########## smoothing parameters for pcf ############################  
  # arguments for 'density'

  denargs <- resolve.defaults(list(kernel=kernel, bw=bw),
                              list(...),
                              list(n=length(r), from=0, to=rmax))
  
  #################################################
  
  ## compute pairwise distances
  
  ## identify close pairs of points
  what <- if(any(correction == "translate")) "all" else "ijd"
  close <- crosspairs(XI, XJ, rmax+hmax, what=what)
  ## map (i,j) to original serial numbers in X
  orig <- seq_len(npts)
  imap <- orig[I]
  jmap <- orig[J]
  iX <- imap[close$i]
  jX <- jmap[close$j]
  ## eliminate any identical pairs
  if(nIJ > 0) {
    ok <- (iX != jX)
    if(!all(ok))
      close <- as.list(as.data.frame(close)[ok, , drop=FALSE])
  }
  ## extract information for these pairs (relative to orderings of XI, XJ)
  dclose <- close$d
  icloseI  <- close$i
#  jcloseJ  <- close$j

  ###### compute #######

  if(any(correction=="translate")) {
    # translation correction
    edgewt <- edge.Trans(dx=close$dx, dy=close$dy, W=win, paired=TRUE)
    gT <- sewpcf(dclose, edgewt, denargs, lambdaIJarea, divisor)$g
    out <- bind.fv(out,
                   data.frame(trans=gT),
                   makefvlabel(NULL, "hat", fname, "Trans"),
                   "translation-corrected estimate of %s",
                   "trans")
  }
  if(any(correction=="isotropic")) {
    # Ripley isotropic correction
    edgewt <- edge.Ripley(XI[icloseI], matrix(dclose, ncol=1))
    gR <- sewpcf(dclose, edgewt, denargs, lambdaIJarea, divisor)$g
    out <- bind.fv(out,
                   data.frame(iso=gR),
                   makefvlabel(NULL, "hat", fname, "Ripley"),
                   "isotropic-corrected estimate of %s",
                   "iso")
  }
  
  ## sanity check
  if(is.null(out)) {
    warning("Nothing computed - no edge corrections chosen")
    return(NULL)
  }
  
  # which corrections have been computed?
  corrxns <- rev(setdiff(names(out), "r"))

  # default is to display them all
  formula(out) <- . ~ r
  fvnames(out, ".") <- corrxns

  # 
  unitname(out) <- unitname(X)
  return(out)
}