#
#
#     markcorr.R
#
#     $Revision: 1.22 $ $Date: 2007/10/24 09:41:15 $
#
#    Estimate the mark correlation function
#
#
# ------------------------------------------------------------------------

"markcorr"<-
function(X, f = function(m1, m2) { m1 * m2 }, r=NULL, 
         correction=c("isotropic", "Ripley", "translate"),
         method="density", ...)
{
	verifyclass(X, "ppp")
        stopifnot(is.marked(X))

        fmul <- function(m1, m2) { m1 * m2 }
        fequ <- function(m1, m2) { m1 == m2 }
        if(missing(f)) {
          if(is.multitype(X)) {
            f <- fequ
            ftype <- "equality"
          } else {
            f <- fmul
            ftype <- "multiplication"
          }
        } else {
          if(!is.function(f))
            stop("Argument f must be a function")
          same <- function(f, g) {
            environment(g) <- environment(f)
            identical(f,g)
          }
          ftype <-
            if(same(f, fmul)) "multiplication" else
            if(same(f, fequ)) "equality" else "unknown"
          if(ftype == "multiplication" && is.multitype(X))
            stop(paste("Inappropriate choice of function f;",
                       "point pattern is multitype;",
                       "types cannot be multiplied."))
        }

        ##
        
        npoints <- X$n
        W <- X$window

        # r values 
        rmaxdefault <- rmax.rule("K", W, npoints/area.owin(W))
        breaks <- handle.r.b.args(r, NULL, W, rmaxdefault=rmaxdefault)
        r <- breaks$r
        rmax <- breaks$max
        
        if(length(method) > 1)
          stop("Select only one method, please")
        if(method=="density" && !breaks$even)
          stop(paste("Evenly spaced r values are required if method=",
                     sQuote("density"), sep=""))
        
        # available selection of edge corrections depends on window
        if(W$type == "mask") {
          iso <- (correction == "isotropic") | (correction == "Ripley")
          if(any(iso)) {
            if(!missing(correction))
              warning("Isotropic correction not implemented for binary masks")
            correction <- correction[!iso]
          }
        }
         
        # this will be the output data frame
        result <- data.frame(r=r, theo= rep(1,length(r)))
        desc <- c("distance argument r", "theoretical Poisson m(r)=1")
        alim <- c(0, min(rmax, rmaxdefault))
        result <- fv(result,
                     "r", substitute(m(r), NULL), "theo", , alim, c("r","mpois(r)"), desc)


        # find close pairs of points
        close <- closepairs(X, rmax)
        dIJ <- close$d
        I   <- close$i
        J   <- close$j
        XI <- ppp(close$xi, close$yi, window=W, check=FALSE)

        # apply f to marks of close pairs of points
        #
        marx <- marks(X, dfok=FALSE)
        ff <- f(marx[I], marx[J])

        if(any(is.na(ff)))
          stop("function f returned some NA values")
        if(is.logical(ff))
          ff <- as.numeric(ff)
        else if(!is.numeric(ff))
          stop("function f did not return numeric values")

        # apply f to every possible pair of marks, and average
        Ef <- switch(ftype,
                     multiplication = {
                       mean(marx)^2
                     },
                     equality = {
                       mtable <- table(marx)
                       sum(mtable^2)/sum(mtable)^2
                     },
                     unknown = {
                       mean(outer(marx, marx, f))
                     },
                     stop("Internal error: invalid ftype"))

        #### Compute estimates ##############
        
        if(any(correction == "translate")) {
          # translation correction
          XJ <- ppp(close$xj, close$yj, window=W, check=FALSE)
          edgewt <- edge.Trans(XI, XJ, paired=TRUE)
          # get smoothed estimate of mark covariance
          Mtrans <- mkcor(dIJ, ff, edgewt, Ef, r, method, ...)
          result <- bind.fv(result,
                            data.frame(trans=Mtrans), "mtrans(r)",
                            "translation-corrected estimate of m(r)",
                            "trans")
        }
        if(any(correction == "isotropic" | correction == "Ripley")) {
          # Ripley isotropic correction
          edgewt <- edge.Ripley(XI, matrix(dIJ, ncol=1))
          # get smoothed estimate of mark covariance
          Miso <- mkcor(dIJ, ff, edgewt, Ef, r, method, ...)
          result <- bind.fv(result,
                            data.frame(iso=Miso), "miso(r)",
                            "Ripley isotropic correction estimate of m(r)",
                            "iso")
        }
        # which corrections have been computed?
        nama2 <- names(result)
        corrxns <- rev(nama2[nama2 != "r"])

        # default is to display them all
        attr(result, "fmla") <- deparse(as.formula(paste(
                       "cbind(",
                        paste(corrxns, collapse=","),
                        ") ~ r")))
        unitname(result) <- unitname(X)
        return(result)
}
	
mkcor <- function(d, ff, wt, Ef, rvals, method="smrep", ..., nwtsteps=500) {
  d <- as.vector(d)
  ff <- as.vector(ff)
  wt <- as.vector(wt)
  switch(method,
         density={
           fw <- ff * wt
           sum.fw <- sum(fw)
           sum.wt <- sum(wt)
           # smooth estimate of kappa_f
           est <- density(d, weights=fw/sum.fw,
                          from=min(rvals), to=max(rvals), n=length(rvals),
                          ...)$y
           numerator <- est * sum.fw
           # smooth estimate of kappa_1
           est0 <- density(d, weights=wt/sum.wt, 
                          from=min(rvals), to=max(rvals), n=length(rvals),
                          ...)$y
           denominator <- est0 * Ef * sum.wt
           result <- numerator/denominator
         },
         sm={
           # This is slow!
           require(sm)
           # smooth estimate of kappa_f
           fw <- ff * wt
           est <- sm.density(d, weights=fw,
                             eval.points=rvals,
                             display="none", nbins=0, ...)$estimate
           numerator <- est * sum(fw)/sum(est)
           # smooth estimate of kappa_1
           est0 <- sm.density(d, weights=wt,
                              eval.points=rvals,
                              display="none", nbins=0, ...)$estimate
           denominator <- est0 * (sum(wt)/ sum(est0)) * Ef
           result <- numerator/denominator
         },
         smrep={
           require(sm)

           hstuff <- resolve.defaults(list(...), list(hmult=1, h.weights=NA))
           if(hstuff$hmult == 1 && all(is.na(hstuff$h.weights)))
             warning("default smoothing parameter may be inappropriate")
           
           # use replication to effect the weights (it's faster)
           nw <- round(nwtsteps * wt/max(wt))
           drep.w <- rep(d, nw)
           fw <- ff * wt
           nfw <- round(nwtsteps * fw/max(fw))
           drep.fw <- rep(d, nfw)

           # smooth estimate of kappa_f
           est <- sm.density(drep.fw,
                             eval.points=rvals,
                             display="none", ...)$estimate
           numerator <- est * sum(fw)/sum(est)
           # smooth estimate of kappa_1
           est0 <- sm.density(drep.w,
                              eval.points=rvals,
                              display="none", ...)$estimate
           denominator <- est0 * (sum(wt)/ sum(est0)) * Ef
           result <- numerator/denominator
         },
         loess = {
           # set up data frame
           df <- data.frame(d=d, ff=ff, wt=wt)
           # fit curve to numerator using loess
           fitobj <- loess(ff ~ d, data=df, weights=wt, ...)
           # evaluate fitted curve at desired r values
           Eff <- predict(fitobj, newdata=data.frame(d=rvals))
           # normalise:
           # denominator is the sample mean of all ff[i,j],
           # an estimate of E(ff(M1,M2)) for M1,M2 independent marks
           result <- Eff/Ef
         },
         )
  return(result)
}