#
# linearpcfmulti.R
#
# $Revision: 1.12 $ $Date: 2017/02/07 08:12:05 $
#
# pair correlation functions for multitype point pattern on linear network
#
#

linearpcfdot <- function(X, i, r=NULL, ..., correction="Ang") {
  if(!is.multitype(X, dfok=FALSE)) 
	stop("Point pattern must be multitype")
  marx <- marks(X)
  lev <- levels(marx)
  if(missing(i) || is.null(i)) i <- lev[1L] else
    if(!(i %in% lev)) stop(paste("i = ", i , "is not a valid mark"))  
  I <- (marx == i)
  J <- rep(TRUE, npoints(X))  # i.e. all points
  result <- linearpcfmulti(X, I, J,
                           r=r, correction=correction, ...)
  correction <- attr(result, "correction")
  type <- if(correction == "Ang") "L" else "net"
  result <- rebadge.as.dotfun(result, "g", type, i)
  return(result)
}

linearpcfcross <- function(X, i, j, r=NULL, ..., correction="Ang") {
  if(!is.multitype(X, dfok=FALSE)) 
	stop("Point pattern must be multitype")
  marx <- marks(X)
  lev <- levels(marx)
  if(missing(i) || is.null(i)) i <- lev[1L] else
    if(!(i %in% lev)) stop(paste("i = ", i , "is not a valid mark"))
  if(missing(j) || is.null(j)) j <- lev[2L] else
    if(!(j %in% lev)) stop(paste("j = ", j , "is not a valid mark"))
  #
  if(i == j) {
    result <- linearpcf(X[marx == i], r=r, correction=correction, ...)
  } else {
    I <- (marx == i)
    J <- (marx == j)
    result <- linearpcfmulti(X, I, J, r=r, correction=correction, ...)
  }
  # rebrand
  correction <- attr(result, "correction")
  type <- if(correction == "Ang") "L" else "net"
  result <- rebadge.as.crossfun(result, "g", type, i, j)
  return(result)
}

linearpcfmulti <- function(X, I, J, r=NULL, ..., correction="Ang") {
  stopifnot(inherits(X, "lpp"))
  correction <- pickoption("correction", correction,
                           c(none="none",
                             Ang="Ang",
                             best="Ang"),
                           multi=FALSE)
  
  # extract info about pattern
  np <- npoints(X)
  lengthL <- volume(domain(X))
  # validate I, J
  if(!is.logical(I) || !is.logical(J))
    stop("I and J must be logical vectors")
  if(length(I) != np || length(J) != np)
    stop(paste("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)
  nIandJ <- sum(I & J)
#  lambdaI <- nI/lengthL
#  lambdaJ <- nJ/lengthL
  # compute pcf
  denom <- (nI * nJ - nIandJ)/lengthL
  g <- linearPCFmultiEngine(X, I, J, r=r, denom=denom, correction=correction, ...)
  # set appropriate y axis label
  correction <- attr(g, "correction")
  type <- if(correction == "Ang") "L" else "net"
  g <- rebadge.as.crossfun(g, "g", type, "I", "J")
  attr(g, "correction") <- correction
  return(g)
}

# ................ inhomogeneous ............................

linearpcfdot.inhom <- function(X, i, lambdaI, lambdadot,
                             r=NULL, ..., correction="Ang", normalise=TRUE) {
  if(!is.multitype(X, dfok=FALSE)) 
	stop("Point pattern must be multitype")
  marx <- marks(X)
  lev <- levels(marx)
  if(missing(i)) i <- lev[1L] else
    if(!(i %in% lev)) stop(paste("i = ", i , "is not a valid mark"))  
  I <- (marx == i)
  J <- rep(TRUE, npoints(X))  # i.e. all points
  # compute
  result <- linearpcfmulti.inhom(X, I, J, lambdaI, lambdadot, 
                               r=r, correction=correction, normalise=normalise,
                               ...)
  correction <- attr(result, "correction")
  type <- if(correction == "Ang") "L, inhom" else "net, inhom"
  result <- rebadge.as.dotfun(result, "g", type, i)
  return(result)
}

linearpcfcross.inhom <- function(X, i, j, lambdaI, lambdaJ,
                               r=NULL, ...,
                               correction="Ang", normalise=TRUE) {
  if(!is.multitype(X, dfok=FALSE)) 
	stop("Point pattern must be multitype")
  marx <- marks(X)
  lev <- levels(marx)
  if(missing(i)) i <- lev[1L] else
    if(!(i %in% lev)) stop(paste("i = ", i , "is not a valid mark"))
  if(missing(j)) j <- lev[2L] else
    if(!(j %in% lev)) stop(paste("j = ", j , "is not a valid mark"))
  #
  if(i == j) {
    I <- (marx == i)
    result <- linearpcfinhom(X[I], lambda=lambdaI, r=r,
                           correction=correction, normalise=normalise, ...)
  } else {
    I <- (marx == i)
    J <- (marx == j)
    result <- linearpcfmulti.inhom(X, I, J, lambdaI, lambdaJ,
                                 r=r, correction=correction,
                                 normalise=normalise, ...)
  }
  # rebrand
  correction <- attr(result, "correction")
  type <- if(correction == "Ang") "L, inhom" else "net, inhom"
  result <- rebadge.as.crossfun(result, "g", type, i, j)
  return(result)
}

linearpcfmulti.inhom <- function(X, I, J, lambdaI, lambdaJ,
                               r=NULL, ...,
                               correction="Ang",
                               normalise=TRUE) {
  stopifnot(inherits(X, "lpp"))
  correction <- pickoption("correction", correction,
                           c(none="none",
                             Ang="Ang",
                             best="Ang"),
                           multi=FALSE)
  
  # extract info about pattern
  np <- npoints(X)
  lengthL <- volume(domain(X))
  # validate I, J
  if(!is.logical(I) || !is.logical(J))
    stop("I and J must be logical vectors")
  if(length(I) != np || length(J) != np)
    stop(paste("The length of I and J must equal",
               "the number of points in the pattern"))
	
  if(!any(I)) stop("no points satisfy I")

  # validate lambda vectors
  lambdaI <- getlambda.lpp(lambdaI, X, subset=I, ...)
  lambdaJ <- getlambda.lpp(lambdaJ, X, subset=J, ...)

  # compute pcf
  weightsIJ <- outer(1/lambdaI, 1/lambdaJ, "*")
  denom <- if(!normalise) lengthL else sum(1/lambdaI) 
  g <- linearPCFmultiEngine(X, I, J, r=r,
                            reweight=weightsIJ, denom=denom,
                            correction=correction, ...)
  # set appropriate y axis label
  correction <- attr(g, "correction")
  type <- if(correction == "Ang") "L, inhom" else "net, inhom"
  g <- rebadge.as.crossfun(g, "g", type, "I", "J")
  attr(g, "correction") <- correction
  attr(g, "dangerous") <- union(attr(lambdaI, "dangerous"),
                                attr(lambdaJ, "dangerous"))
  return(g)
}

# .............. internal ...............................

linearPCFmultiEngine <- function(X, I, J, ..., r=NULL, reweight=NULL, denom=1,
                          correction="Ang", showworking=FALSE) {
  # ensure distance information is present
  X <- as.lpp(X, sparse=FALSE)
  # extract info about pattern
  np <- npoints(X)
  # extract linear network
  L <- domain(X)
  W <- Window(L)
  # determine r values
  rmaxdefault <- 0.98 * boundingradius(L)
  breaks <- handle.r.b.args(r, NULL, W, rmaxdefault=rmaxdefault)
  r <- breaks$r
  rmax <- breaks$max
  #
  if(correction == "Ang") {
    fname <- c("g", "list(L, I, J)")
    ylab <- quote(g[L,I,J](r))
  } else {
    fname <- c("g", "list(net, I, J)")
    ylab <- quote(g[net,I,J](r))
  }
  #
   if(np < 2) {
    # no pairs to count: return zero function
    zeroes <- rep(0, length(r))
    df <- data.frame(r = r, est = zeroes)
    g <- fv(df, "r", ylab,
            "est", . ~ r, c(0, rmax),
            c("r", makefvlabel(NULL, "hat", fname)), 
            c("distance argument r", "estimated %s"),
            fname = fname)
    unitname(g) <- unitname(X)
    attr(g, "correction") <- correction
    return(g)
  }
  #
  nI <- sum(I)
  nJ <- sum(J)
  whichI <- which(I)
  whichJ <- which(J)
  clash <- I & J
  has.clash <- any(clash)
  # compute pairwise distances
  if(exists("crossdist.lpp")) {
    DIJ <- crossdist(X[I], X[J], check=FALSE)
    if(has.clash) {
      # exclude pairs of identical points from consideration
      Iclash <- which(clash[I])
      Jclash <- which(clash[J])
      DIJ[cbind(Iclash,Jclash)] <- Inf
    }
  } else {
    D <- pairdist(X)
    diag(D) <- Inf
    DIJ <- D[I, J]
  }
  #---  compile into pair correlation function ---
  if(correction == "none" && is.null(reweight)) {
    # no weights (Okabe-Yamada)
    g <- compilepcf(DIJ, r, denom=denom, check=FALSE, fname=fname)
    g <- rebadge.as.crossfun(g, "g", "net", "I", "J")    
    unitname(g) <- unitname(X)
    attr(g, "correction") <- correction
    return(g)
  }
  if(correction == "none")
     edgewt <- 1
  else {
     # inverse m weights (Ang's correction)
     # determine tolerance
     toler <- default.linnet.tolerance(L)
     # compute m[i,j]
     m <- matrix(1, nI, nJ)
     XI <- X[I]
     if(!has.clash) {
       for(k in seq_len(nJ)) {
         j <- whichJ[k]
         m[,k] <- countends(L, XI, DIJ[, k], toler=toler)
       }
     } else {
       # don't count identical pairs
       for(k in seq_len(nJ)) {
         j <- whichJ[k]
         inotj <- (whichI != j)
         m[inotj, k] <- countends(L, XI[inotj], DIJ[inotj, k], toler=toler)
       }
     }
     edgewt <- 1/m
  }
  # compute pcf
  wt <- if(!is.null(reweight)) edgewt * reweight else edgewt
  g <- compilepcf(DIJ, r, weights=wt, denom=denom, check=FALSE, ...,
                  fname=fname)
  ## rebadge and tweak
  g <- rebadge.as.crossfun(g, "g", "L", "I", "J")
  fname <- attr(g, "fname")
  # tack on theoretical value
  g <- bind.fv(g, data.frame(theo=rep(1,length(r))),
               makefvlabel(NULL, NULL, fname, "pois"),
               "theoretical Poisson %s")
  unitname(g) <- unitname(X)
  fvnames(g, ".") <- rev(fvnames(g, "."))
  # show working
  if(showworking)
    attr(g, "working") <- list(DIJ=DIJ, wt=wt)
  attr(g, "correction") <- correction
  return(g)
}