#
#   pcfinhom.R
#
#   $Revision: 1.6 $   $Date: 2011/07/23 02:37:00 $
#
#   inhomogeneous pair correlation function of point pattern 
#
#

pcfinhom <- function(X, lambda=NULL, ..., r=NULL,
                     kernel="epanechnikov", bw=NULL, stoyan=0.15,
                     correction=c("translate", "Ripley"),
                     renormalise=TRUE,
                     normpower=1,
                     reciplambda=NULL, 
                     sigma=NULL, varcov=NULL)
{
  verifyclass(X, "ppp")
  r.override <- !is.null(r)

  win <- X$window
  area <- area.owin(win)
  npts <- npoints(X)
  
  correction.given <- !missing(correction)
  correction <- pickoption("correction", correction,
                           c(isotropic="isotropic",
                             Ripley="isotropic",
                             translate="translate",
                             best="best"),
                           multi=TRUE)

  correction <- implemented.for.K(correction, win$type, correction.given)
  
  if(is.null(bw) && kernel=="epanechnikov") {
    # Stoyan & Stoyan 1995, eq (15.16), page 285
    h <- stoyan /sqrt(npts/area)
    # conversion to standard deviation
    bw <- h/sqrt(5)
  }

  ########## intensity values #########################

  if(missing(lambda) && is.null(reciplambda)) {
    # No intensity data provided
    # Estimate density by leave-one-out kernel smoothing
    lambda <- density(X, ..., sigma=sigma, varcov=varcov,
                      at="points", leaveoneout=TRUE)
    lambda <- as.numeric(lambda)
    reciplambda <- 1/lambda
  } else if(!is.null(reciplambda)) {
    # 1/lambda values provided
    if(is.im(reciplambda)) 
      reciplambda <- safelookup(reciplambda, X)
    else if(is.function(reciplambda))
      reciplambda <- reciplambda(X$x, X$y)
    else if(is.numeric(reciplambda) && is.vector(as.numeric(reciplambda)))
      check.nvector(reciplambda, npts)
    else stop(paste(sQuote("reciplambda"),
                    "should be a vector, a pixel image, or a function"))
  } else {
    # lambda values provided
    if(is.im(lambda)) 
      lambda <- safelookup(lambda, X)
    else if(is.ppm(lambda))
      lambda <- predict(lambda, locations=X, type="trend")
    else if(is.function(lambda)) 
      lambda <- lambda(X$x, X$y)
    else if(is.numeric(lambda) && is.vector(as.numeric(lambda)))
      check.nvector(lambda, npts)
    else stop(paste(sQuote("lambda"),
                    "should be a vector, a pixel image, or a function"))
    # evaluate reciprocal
    reciplambda <- 1/lambda
  }
  
  # renormalise
  if(renormalise) {
    check.1.real(normpower)
    stopifnot(normpower %in% 1:2)
    renorm.factor <- (area/sum(reciplambda))^normpower
  } 
  
  ########## r values ############################
  # handle arguments r and breaks 

  rmaxdefault <- rmax.rule("K", win, lambda)        
  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))

  ########## smoothing parameters for pcf ############################  
  # arguments for 'density'

  from <- 0
  to <- max(r)
  nr <- length(r)
  
  denargs <- resolve.defaults(list(kernel=kernel, bw=bw),
                              list(...),
                              list(n=nr, from=from, to=to))
  
  #################################################
  
  # compute pairwise distances
  
  close <- closepairs(X, max(r))
  dIJ <- close$d
  I <- close$i
  J <- close$j
  XI <- ppp(close$xi, close$yi, window=win, check=FALSE)
  wIJ <- reciplambda[I] * reciplambda[J]

  # initialise fv object
  
  df <- data.frame(r=r, theo=rep(1,length(r)))
  out <- fv(df, "r",
            substitute(g[inhom](r), NULL), "theo", ,
            alim,
            c("r","{%s^{Pois}}(r)"),
            c("distance argument r", "theoretical Poisson %s"),
            fname="g[inhom]")

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

  if(any(correction=="translate")) {
    # translation correction
    XJ <- ppp(close$xj, close$yj, window=win, check=FALSE)
    edgewt <- edge.Trans(XI, XJ, paired=TRUE)
    gT <- sewpcf(dIJ, edgewt * wIJ, denargs, area)$g
    if(renormalise) gT <- gT * renorm.factor
    out <- bind.fv(out,
                   data.frame(trans=gT),
                   "hat(%s^{Trans})(r)",
                   "translation-corrected estimate of %s",
                   "trans")
  }
  if(any(correction=="isotropic")) {
    # Ripley isotropic correction
    edgewt <- edge.Ripley(XI, matrix(dIJ, ncol=1))
    gR <- sewpcf(dIJ, edgewt * wIJ, denargs, area)$g
    if(renormalise) gR <- gR * renorm.factor
    out <- bind.fv(out,
                   data.frame(iso=gR),
                   "hat(%s^{Ripley})(r)",
                   "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?
  nama2 <- names(out)
  corrxns <- rev(nama2[nama2 != "r"])

  # default is to display them all
  attr(out, "fmla") <- deparse(as.formula(paste(
                       "cbind(",
                        paste(corrxns, collapse=","),
                        ") ~ r")))
  unitname(out) <- unitname(X)
  return(out)
}