#
#	Kmulti.S		
#
#	Compute estimates of cross-type K functions
#	for multitype point patterns
#
#	$Revision: 5.22 $	$Date: 2008/06/09 17:49:58 $
#
#
# -------- functions ----------------------------------------
#	Kcross()	cross-type K function K_{ij}
#                       between types i and j
#
#	Kdot()          K_{i\bullet}
#                       between type i and all points regardless of type
#
#       Kmulti()        (generic)
#
#
# -------- standard arguments ------------------------------	
#	X		point pattern (of class 'ppp')
#				including 'marks' vector
#	r		distance values at which to compute K	
#
# -------- standard output ------------------------------
#      A data frame with columns named
#
#	r:		same as input
#
#	trans:		K function estimated by translation correction
#
#	iso:		K function estimated by Ripley isotropic correction
#
#	theo:		K function for Poisson ( = pi * r ^2 )
#
#	border:		K function estimated by border method
#			using standard formula (denominator = count of points)
#
#       bord.modif:	K function estimated by border method
#			using modified formula 
#			(denominator = area of eroded window
#
# ------------------------------------------------------------------------

"Lcross" <- function(X, i, j, ...) {
  if(missing(i)) i <- levels(marks(X))[1]
  if(missing(j)) j <- levels(marks(X))[2]
  K <- Kcross(X, i, j, ...)
  L <- eval.fv(sqrt(K/pi))
  # relabel the fv object
  attr(L, "ylab") <- substitute(Lcross[i,j](r), list(i=paste(i),j=paste(j)))
  nama <- names(L)
  desc <- attr(L, "desc")
  labl <- attr(L, "labl")
  TAGS <- c("theo", "un", "bord.modif", "border", "trans", "iso")
  DESC <- c("theoretical Poisson Lcross(r)",
             "uncorrected estimate of Lcross(r)",
             "modified border-corrected estimate of Lcross(r)",
             "border-corrected estimate of Lcross(r)",
             "translation-corrected estimate of Lcross(r)",
             "Ripley isotropic correction estimate of Lcross(r)")
  LABL <- c("Lpois(r)",
             "Lun(r)",
             "Lbord*(r)",
             "Lbord(r)",
             "Ltrans(r)",
             "Liso(r)")

  for(i in seq(length(TAGS)))
    if(!is.na(m <- match(TAGS[i], nama))) {
      desc[m] <- DESC[i]
      labl[m] <- LABL[i]
    }
  attr(L, "desc") <- desc
  attr(L, "labl") <- labl
  return(L)  
}

"Ldot" <- function(X, i, ...) {
  if(missing(i)) i <- levels(marks(X))[1]
  K <- Kdot(X, i, ...)
  L <- eval.fv(sqrt(K/pi))
  # relabel the fv object
  attr(L, "ylab") <- substitute(Ldot[i](r), list(i=paste(i)))
  nama <- names(L)
  desc <- attr(L, "desc")
  labl <- attr(L, "labl")
  TAGS <- c("theo", "un", "bord.modif", "border", "trans", "iso")
  DESC <- c("theoretical Poisson Ldot(r)",
             "uncorrected estimate of Ldot(r)",
             "modified border-corrected estimate of Ldot(r)",
             "border-corrected estimate of Ldot(r)",
             "translation-corrected estimate of Ldot(r)",
             "Ripley isotropic correction estimate of Ldot(r)")
  LABL <- c("Lpois(r)",
             "Lun(r)",
             "Lbord*(r)",
             "Lbord(r)",
             "Ltrans(r)",
             "Liso(r)")

  for(i in seq(length(TAGS)))
    if(!is.na(m <- match(TAGS[i], nama))) {
      desc[m] <- DESC[i]
      labl[m] <- LABL[i]
    }
  attr(L, "desc") <- desc
  attr(L, "labl") <- labl
  return(L)  
}


"Kcross" <- 
function(X, i, j, r=NULL, breaks=NULL,
         correction =c("border", "isotropic", "Ripley", "translate") , ...)
{
	verifyclass(X, "ppp")
	if(!is.marked(X))
		stop("point pattern has no marks (no component 'marks')")
        marx <- marks(X)
        if(missing(i))
          i <- levels(marx)[1]
        if(missing(j))
          j <- levels(marx)[2]
        I <- (marx == i)
        if(!any(I))
          stop(paste("No points have mark i =", i))

        if(i == j) {
          result <- Kest(X[I],
                         r=r, breaks=breaks, correction=correction, ...)
        } else {
          J <- (marx == j)
          if(!any(J))
            stop(paste("No points have mark j =", j))
          result <- Kmulti(X, I, J,
                           r=r, breaks=breaks, correction=correction, ...)
        }
        attr(result, "ylab") <- substitute(Kcross[i,j](r), list(i=paste(i),j=paste(j)))
        result
}

"Kdot" <- 
function(X, i, r=NULL, breaks=NULL,
         correction = c("border", "isotropic", "Ripley", "translate") , ...)
{
	verifyclass(X, "ppp")
	if(!is.marked(X))
		stop("point pattern has no marks (no component 'marks')")

        marx <- marks(X)
        if(missing(i))
          i <- levels(marx)[1]
        
	I <- (marx == i)
	J <- rep(TRUE, X$n)  # i.e. all points
	
	if(!any(I)) stop(paste("No points have mark i =", i))
	
	result <- Kmulti(X, I, J,
                         r=r, breaks=breaks, correction=correction, ...)
        attr(result, "ylab") <- substitute(Kdot[i](r), list(i=paste(i)))
        result
}


"Kmulti"<-
function(X, I, J, r=NULL, breaks=NULL,
         correction = c("border", "isotropic", "Ripley", "translate") , ...)
{
	verifyclass(X, "ppp")

	npoints <- X$n
	x <- X$x
	y <- X$y
        W <- X$window
	area <- area.owin(W)

        # available selection of edge corrections depends on window
        if(W$type == "mask") {
           iso <- (correction == "isotropic") | (correction == "Ripley")
           if(all(iso))
             stop("Isotropic correction not implemented for binary masks")
           if(any(iso)) {
             if(!missing(correction))
               warning("Isotropic correction not implemented for binary masks")
             correction <- correction[!iso]
           }
        }
         
	if(!is.logical(I) || !is.logical(J))
		stop("I and J must be logical vectors")
	if(length(I) != npoints || length(J) != npoints)
	     stop("The length of I and J must equal \
 the number of points in the pattern")
	
	if(!any(I)) stop("no points satisfy I")
	if(!any(J)) stop("no points satisfy J")
		
	nI <- sum(I)
	nJ <- sum(J)
	lambdaI <- nI/area
	lambdaJ <- nJ/area

        # r values 
        rmaxdefault <- rmax.rule("K", W, lambdaJ)
        breaks <- handle.r.b.args(r, breaks, W, rmaxdefault=rmaxdefault)
        r <- breaks$r
        rmax <- breaks$max
        
        # recommended range of r values
        alim <- c(0, min(rmax, rmaxdefault))
        
        # this will be the output data frame
        # It will be given more columns later
        K <- data.frame(r=r, theo= pi * r^2)
        desc <- c("distance argument r", "theoretical Poisson K(r)")
        K <- fv(K, "r", substitute(Kmulti(r), NULL),
                "theo", , alim, c("r","Kpois(r)"), desc)

# find close pairs of points
        XI <- X[I]
        XJ <- X[J]
        close <- crosspairs(XI, XJ, max(r))
# close$i and close$j are serial numbers in XI and XJ respectively;        
# map them to original serial numbers in X
        orig <- seq(npoints)
        imap <- orig[I]
        jmap <- orig[J]
        iX <- imap[close$i]
        jX <- jmap[close$j]
# eliminate any identical pairs
        if(any(I & J)) {
          ok <- (iX != jX)
          if(!all(ok)) {
            close$i  <- close$i[ok]
            close$j  <- close$j[ok]
            close$xi <- close$xi[ok]
            close$yi <- close$yi[ok]
            close$xj <- close$xj[ok]
            close$yj <- close$yj[ok]
            close$dx <- close$dx[ok]
            close$dy <- close$dy[ok]
            close$d  <- close$d[ok]
          }
        }
# extract information for these pairs (relative to orderings of XI, XJ)
        dcloseIJ <- close$d
        icloseI  <- close$i
        jcloseJ  <- close$j
        
# Compute estimates by each of the selected edge corrections.
        
        if(any(correction == "border" | correction == "bord.modif")) {
          # border method
          # distance to boundary from each point of type I
          bI <- bdist.points(XI)
          # distance to boundary from first element of each (i, j) pair
          bcloseI <- bI[icloseI]
          # apply reduced sample algorithm
          RS <- Kount(dcloseIJ, bcloseI, bI, breaks)
          if(any(correction == "bord.modif")) {
            denom.area <- eroded.areas(W, r)
            Kbm <- RS$numerator/(denom.area * nI * nJ)
            K <- bind.fv(K, data.frame(bord.modif=Kbm), "Kbord*(r)",
                         "modified border-corrected estimate of K(r)",
                         "bord.modif")
          }
          if(any(correction == "border")) {
            Kb <- RS$numerator/(lambdaJ * RS$denom.count)
            K <- bind.fv(K, data.frame(border=Kb), "Kbord(r)",
                         "border-corrected estimate of K(r)",
                         "border")
          }
        }
        if(any(correction == "translate")) {
          # translation correction
            edgewt <- edge.Trans(XI[icloseI], XJ[jcloseJ], paired=TRUE)
            wh <- whist(dcloseIJ, breaks$val, edgewt)
            Ktrans <- cumsum(wh)/(lambdaI * lambdaJ * area)
            rmax <- diameter(W)/2
            Ktrans[r >= rmax] <- NA
            K <- bind.fv(K, data.frame(trans=Ktrans), "Ktrans(r)",
                         "translation-corrected estimate of K(r)",
                         "trans")
        }
        if(any(correction == "isotropic" | correction == "Ripley")) {
          # Ripley isotropic correction
            edgewt <- edge.Ripley(XI[icloseI], matrix(dcloseIJ, ncol=1))
            wh <- whist(dcloseIJ, breaks$val, edgewt)
            Kiso <- cumsum(wh)/(lambdaI * lambdaJ * area)
            rmax <- diameter(W)/2
            Kiso[r >= rmax] <- NA
            K <- bind.fv(K, data.frame(iso=Kiso), "Kiso(r)",
                         "Ripley isotropic correction estimate of K(r)",
                         "iso")
        }
        # default is to display them all
        attr(K, "fmla") <- . ~ r
        unitname(K) <- unitname(X)
        return(K)
          
}