#
#    QQ plot of smoothed residual field against model
#
#  qqplot.ppm()       QQ plot (including simulation)
#
#  $Revision: 1.5 $   $Date: 2005/05/24 21:11:26 $
#

qqplot.ppm <-
  function(fit, nsim=100, expr=NULL, ..., type="raw", style="mean",
           fast=TRUE, verbose=TRUE, plot.it=TRUE,
           dimyx=NULL, nrep=if(fast) 5e4 else 1e5,
           control=list(nrep=nrep, expand=default.expand(fit),
           periodic=TRUE), saveall=FALSE)
{
  verifyclass(fit, "ppm")

  if(fast) {
    oldnpixel <- spatstat.options("npixel")[[1]]
    if(is.null(dimyx)) 
      dimyx <- pmin(40, oldnpixel)
    spatstat.options(npixel=dimyx)
  } 
    
  ################   How to evaluate residuals ##########################
  
  # Quantiles of the residual field will be computed.
  residualfield <- function(fit, ...) {
    d <- diagnose.ppm(fit, which="smooth",
                      plot.it=FALSE, compute.cts=FALSE, compute.sd=FALSE, ...)
    return(d$smooth$Z$v)
  }

  # Data values
  dat <- residualfield(fit, type=type, ..., dimyx=dimyx)

  # How to refit the model properly!
  refit <- function(fit, pattern) {
    update.ppm(fit, Q=pattern)
  }
  
  ##################  How to perform simulations?  #######################

  simulate.from.fit <- is.null(expr)
  
  if(!simulate.from.fit) {
     # 'expr' will be evaluated 'nsim' times
    if(!is.expression(expr))
      stop("Argument \'expr\' should be an expression")
  } else{
    # We will simulate from the fitted model 'nsim' times
    # and refit the model to these simulations

    # prepare rmh arguments
    rcontrol <- rmhcontrol(control)
    rmodel   <- rmhmodel(fit, control=rcontrol, project=TRUE, verbose=verbose)
    rstart   <- rmhstart(n.start=data.ppm(fit)$n)

    # expression to be evaluated each time
    expr <- expression(
        refit(fit, 
              rmh.default(rmodel, rstart, rcontrol, verbose=FALSE)))
   }

  ######  Perform simulations
  if(verbose) cat(paste("Simulating", nsim, "realisations... "))
  for(i in 1:nsim) {
    ei <- eval(expr)
    fiti <-
      if(simulate.from.fit)
        ei
      else if(is.ppm(ei))
        ei
      else if(is.ppp(ei))
        refit(fit, ei)
      else
        stop("result of eval(expr) is not a ppm or ppp object")
    resi <- residualfield(fiti, type=type, ..., dimyx=dimyx)
    if(i == 1)
      sim <- array(, dim=c(dim(resi), nsim))
    sim[,,i] <- resi
    if(verbose) 
      cat(paste(i, ifelse(i %% 10 == 0, "\n", ", "), sep=""))
  }

  ############ Plot them
  switch(style,
         classical = {
           rr <- range(c(dat,sim))
           result <- qqplot(sim, dat, xlim=rr, ylim=rr, asp=1.0,
                            xlab="Quantiles of simulation",
                            ylab="Quantiles of data",plot.it=plot.it)
           title(sub=paste("Residuals:", type))
           abline(0,1, lty=2)
           result <- append(result,
                            list(data=dat,
                                 sim=sim,
                                 xlim=rr,
                                 ylim=rr,
                                 xlab="Quantiles of simulation",
                                 ylab="Quantiles of data",
                                 rtype=type,
                                 nsim=nsim,
                                 fit=fit,
                                 expr=
                                 if(simulate.from.fit) NULL else deparse(expr),
                                 simulate.from.fit=simulate.from.fit
                                 )
                            )
         },
         mean = {
           # compute quantiles corresponding to probabilities p[i]
           # separately in each realisation.
           if(verbose) cat("Calculating quantiles...")
           if(fast) {
             p <- ppoints(min(100,length(dat)), 3/8)
             qsim <- apply(sim, 3, quantile, probs=p, na.rm=TRUE)
           } else {
             qsim <- apply(sim, 3, sort, na.last=TRUE)
           }
           if(verbose) cat("averaging...")
           # sample mean of each quantile
           meanq <- apply(qsim, 1, mean, na.rm=TRUE)
           # et cetera
           varq <- apply(qsim, 1, var, na.rm=TRUE)
           sdq <- sqrt(varq)
           q.025 <- apply(qsim, 1, quantile, probs=0.025, na.rm=TRUE)
           q.975 <- apply(qsim, 1, quantile, probs=0.975, na.rm=TRUE)
  
           rr <- range(c(meanq,dat), na.rm=TRUE)

           dats <- if(fast) quantile(dat, probs=p, na.rm=TRUE) else sort(dat, na.last=TRUE)

           if(verbose) cat("..Done.\n")
           if(plot.it) {
           	plot(meanq, dats,
                     xlab="Mean quantile of simulations", ylab="data quantile",
                     xlim=rr, ylim=rr, asp=1.0)
                abline(0,1)
                lines(meanq, q.025, lty=2, col="red")
                lines(meanq, q.975, lty=2, col="red")
                title(sub=paste("Residuals:", type))
           }
           result <- list(x=meanq, y=dats, sdq=sdq,
                          q.025=q.025, q.975=q.975,
                          data=dat, sim=sim,
                          xlim=rr, ylim=rr,
                          xlab="Mean quantile of simulations",
                          ylab="data quantile",
                          rtype=type,
                          nsim=nsim,
                          fit=fit,
                          expr=if(simulate.from.fit) NULL else deparse(expr),
                          simulate.from.fit=simulate.from.fit)
         },
         stop("Unrecognised option for \"style\"")
       )

  # Throw out baggage if not wanted         
  if(!saveall) {
    result$fit <- summary(fit, quick=TRUE)
    result$sim <- NULL
  }
         
  # reset npixel
  if(fast)
    spatstat.options(npixel=oldnpixel)
  #
  class(result) <- c("qqppm", class(result))
  return(invisible(result))
}

plot.qqppm <- function(x, ..., limits=TRUE) {
  stopifnot(inherits(x, "qqppm"))
  default.type <- if(length(x$x) > 150) "l" else "p"
  myplot <- function(object,
                     xlab = object$xlab, ylab = object$ylab,
                     xlim = object$xlim, ylim = object$ylim,
                     asp = 1,
                     type = default.type,
                     ..., limits=TRUE) {
    plot(object$x, object$y, xlab = xlab, ylab = ylab,
         xlim = xlim, ylim = ylim, asp = asp, type = type, ...)
    abline(0, 1)
    if(limits) {
      if(!is.null(object$q.025))
        lines(object$x, object$q.025, lty = 2, col="red")
      if(!is.null(object$q.975))
        lines(object$x, object$q.975, lty = 2, col="red")
    }
    title(sub=paste("Residuals:", object$rtype))
  }
  myplot(x, ..., limits=limits)
  return(invisible(x))
}

  
print.qqppm <- function(x, ...) {
  stopifnot(inherits(x, "qqppm"))
  cat(paste("Q-Q plot of point process residuals ",
            "of type \`", x$rtype, "\'\n",
            "based on ", x$nsim, " simulations\n",
            sep=""))
  if(x$simulate.from.fit) {
    fit  <- x$fit
    sumfit <- if(is.ppm(fit)) summary(fit, quick=TRUE)
              else if(inherits(fit, "summary.ppm")) fit
              else list(name="(unrecognised format)")
    cat(paste("\nSimulations from fitted model:",
              sumfit$name, "\n"))
  } else {
    cat("Simulations obtained by evaluating the following expression:\n")
    print(x$expr)
  } 
  invisible(NULL)
}