#
#  clftest.R
#
#  $Revision: 1.9 $  $Date: 2012/05/15 09:07:53 $
#
#  Monte Carlo tests for CSR (etc)
#

clf.test <- function(X, ..., rinterval=NULL, use.theo=FALSE) {
  Xname <- deparse(substitute(X))
  envelopeTest(X, ..., power=2,
               use.theo=use.theo, rinterval=rinterval, Xname=Xname)
}

mad.test <- function(X, ..., rinterval=NULL, use.theo=FALSE) {
  Xname <- deparse(substitute(X))
  envelopeTest(X, ..., power=Inf,
               use.theo=use.theo, rinterval=rinterval, Xname=Xname)
}

envelopeTest <- function(X, ...,
                         power=1, rinterval=NULL,
                         use.theo=FALSE, Xname=NULL, verbose=TRUE,
                         internal=NULL) {
  if(is.null(Xname)) Xname <- deparse(substitute(X))
  check.1.real(power)
  explain.ifnot(power >= 0)
  # compute or extract simulated functions
  if(use.theo) {
    # using theoretical function as reference.
    # ensure envelope object includes theoretical function.
    internal <- resolve.defaults(internal, list(csr=TRUE))
  }
  X <- envelope(X, ..., savefuns=TRUE, Yname=Xname,
                        internal=internal, verbose=verbose)
  Y <- attr(X, "simfuns")
  # extract values
  r   <- with(X, .x)
  obs <- with(X, .y)
  sim <- as.matrix(as.data.frame(Y))[, -1]
  nsim <- ncol(sim)
  # choose function as reference
  has.theo <- ("theo" %in% names(X))
  if(use.theo && !has.theo)
    warning("No theoretical function available; use.theo ignored")
  if(use.theo && has.theo) {
    reference <- with(X, theo)
    used.theo <- TRUE
  } else {
    # compute sample mean of simulations *and* observed 
    reference <- apply(cbind(sim, obs), 1, mean, na.rm=TRUE)
    used.theo <- FALSE
  }
  # determine interval of r values for computation
  if(!is.null(rinterval)) {
    stopifnot(is.numeric(rinterval))
    stopifnot(length(rinterval) == 2)
    stopifnot(rinterval[1] < rinterval[2])
    if(max(r) < rinterval[2]) {
      oldrinterval <- rinterval
      rinterval <- intersect.ranges(rinterval, range(r))
      if(verbose)
        warning(paste("The interval", prange(oldrinterval),
                      "is too large for the available data;",
                      "it has been trimmed to", prange(rinterval)))
    }
    ok <- (rinterval[1] <= r & r <= rinterval[2])
    obs <- obs[ok]
    sim <- sim[ok, ]
    reference <- reference[ok]
  } else rinterval <- range(r)

  # compute test statistic
  if(is.infinite(power)) {
    # MAD
    devdata <- max(abs(obs-reference))
    names(devdata) <- "mad"
    devsim <- apply(abs(sim-reference), 2, max)
    testname <- "Maximum absolute deviation test"
  } else {
    a <- diff(rinterval) * (if(used.theo) 1 else (nsim/(nsim - 1))^power)
    if(power == 2) {
      # Cramer-von Mises
      devdata <- a * mean((obs - reference)^2)
      names(devdata) <- "u"
      devsim <- a * colMeans((sim - reference)^2)
      testname <- "Cressie-Loosmore-Ford test"
    } else if(power == 1) {
      # integral absolute deviation
      devdata <- a * mean(abs(obs - reference))
      names(devdata) <- "L1"
      devsim <- a * colMeans(abs(sim - reference))
      testname <- "Integral absolute deviation test"
    } else {
      # general p
      devdata <- a * mean((abs(obs - reference)^power))
      names(devdata) <- "Lp"
      devsim <- a * colMeans((abs(sim - reference)^power))
      testname <- paste("Integrated", ordinal(power), "Power Deviation test")
    }
  }
  # compute rank
  datarank <- sum(devdata < devsim) + sum(devdata == devsim)/2 + 1
  pvalue <- datarank/(nsim+1)
  # bookkeeping
  statistic <- data.frame(devdata, rank=datarank)
  colnames(statistic)[1] <- names(devdata)
  e <- attr(X, "einfo")
  nullmodel <-
    if(identical(e$csr, TRUE)) "CSR" else 
    if(!is.null(e$simtype)) {
      switch(e$simtype,
             csr = "CSR",
             rmh = paste("fitted",
               if(identical(e$pois, TRUE)) "Poisson" else "Gibbs",
               "model"),
             kppm = "fitted cluster model",
             expr = "model simulated by evaluating expression",
             list = "model simulated by drawing patterns from a list",
             "unrecognised model")
    } else "unrecognised model"
  fname <- deparse(attr(X, "ylab"))
  testname <- c(paste(testname, "of", nullmodel),
                paste("Monte Carlo test based on", nsim, "simulations"),
                paste("Summary function:", fname),
                paste("Reference function:",
                      if(used.theo) "theoretical" else "sample mean"))
  result <- structure(list(statistic = statistic,
                           p.value = pvalue,
                           method = testname,
                           data.name = e$Yname),
                      class = "htest")
  return(result)
}