#
#   envelope.R
#
#   computes simulation envelopes 
#
#   $Revision: 2.112 $  $Date: 2022/04/06 02:11:50 $
#

envelope <- function(Y, fun, ...) {
  UseMethod("envelope")
}

  # .................................................................
  #     A 'simulation recipe' contains the following variables
  #
  #  type = Type of simulation
  #           "csr": uniform Poisson process
  #           "rmh": simulated realisation of fitted Gibbs or Poisson model 
  #          "kppm": simulated realisation of fitted cluster model 
  #          "slrm": simulated realisation of spatial logistic regression
  #          "expr": result of evaluating a user-supplied expression
  #          "list": user-supplied list of point patterns
  #
  #  expr = expression that is repeatedly evaluated to generate simulations
  #
  #    envir = environment in which to evaluate the expression `expr'
  #
  #    'csr' = TRUE iff the model is (known to be) uniform Poisson
  #
  #    pois  = TRUE if model is known to be Poisson
  #  
  #  constraints = additional information about simulation (e.g. 'with fixed n')
  #
  # ...................................................................

simulrecipe <- function(type, expr, envir, csr, pois=csr, constraints="") {
  if(csr && !pois) warning("Internal error: csr=TRUE but pois=FALSE")
  out <- list(type=type,
              expr=expr,
              envir=envir,
              csr=csr,
              pois=pois,
              constraints=constraints)
  class(out) <- "simulrecipe"
  out
}


envelope.ppp <-
  function(Y, fun=Kest, nsim=99, nrank=1, ...,
           funargs=list(), funYargs=funargs,
           simulate=NULL, fix.n=FALSE, fix.marks=FALSE,
           verbose=TRUE, clipdata=TRUE, 
           transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
           alternative=c("two.sided", "less", "greater"),
           scale=NULL, clamp=FALSE, 
           savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
           VARIANCE=FALSE, nSD=2,
           Yname=NULL,
           maxnerr=nsim, rejectNA=FALSE, silent=FALSE,
           do.pwrong=FALSE,
           envir.simul=NULL) {
  cl <- short.deparse(sys.call())
  if(is.null(Yname)) Yname <- short.deparse(substitute(Y))
  if(is.null(fun)) fun <- Kest
  envir.user <- if(!is.null(envir.simul)) envir.simul else parent.frame()
  envir.here <- sys.frame(sys.nframe())

  ismarked <- is.marked(Y)
  ismulti  <- is.multitype(Y)
  fix.marks <- fix.marks && ismarked
  
  if(!is.null(simulate)) {
    # ...................................................
    # Simulations are determined by 'simulate' argument
    if(fix.n || fix.marks) 
      warning("fix.n and fix.marks were ignored, because 'simulate' was given")
    # Processing is deferred to envelopeEngine
    simrecipe <- simulate
    # Data pattern is argument Y
    X <- Y
  } else if(!fix.n && !fix.marks) {
    # ...................................................
    # Realisations of complete spatial randomness
    # will be generated by rpoispp 
    # Data pattern X is argument Y
    # Data pattern determines intensity of Poisson process
    X <- Y
    sY <- summary(Y, checkdup=FALSE)
    Yintens <- sY$intensity
    nY <- npoints(Y)
    Ywin <- Y$window
    Ymarx <- marks(Y)
    # expression that will be evaluated
    simexpr <- if(is.null(Ymarx)) {
        # unmarked point pattern
        expression(rpoispp(Yintens, win=Ywin))
      } else if(is.null(dim(Ymarx))) {
        # single column of marks
        expression({
          A <- rpoispp(Yintens, win=Ywin);
          j <- sample(nY, npoints(A), replace=TRUE);
          A %mark% Ymarx[j]
        })
      } else {
        # multiple columns of marks
        expression({
          A <- rpoispp(Yintens, win=Ywin);
          j <- sample(nY, npoints(A), replace=TRUE);
          A %mark% Ymarx[j, , drop=FALSE]
        })
      }
    dont.complain.about(Yintens, Ywin)
    # evaluate in THIS environment
    simrecipe <- simulrecipe(type = "csr",
                             expr = simexpr,
                             envir = envir.here,
                             csr   = TRUE,
                             pois  = TRUE)
  } else if(fix.marks) {
    # ...................................................
    # Data pattern is argument Y
    X <- Y
    # Realisations of binomial process
    # with fixed number of points and fixed marks
    # will be generated by runifpoint
    nY <- npoints(Y)
    Ywin <- as.owin(Y)
    Ymarx <- marks(Y)
    # expression that will be evaluated
    simexpr <- expression(runifpoint(nY, Ywin) %mark% Ymarx)
    # suppress warnings from code checkers
    dont.complain.about(nY, Ywin, Ymarx)
    # simulation constraints (explanatory string)
    constraints <- if(ismulti) "with fixed number of points of each type" else
                   "with fixed number of points and fixed marks"
    # evaluate in THIS environment
    simrecipe <- simulrecipe(type = "csr",
                             expr = simexpr,
                             envir = envir.here,
                             csr   = TRUE,
                             pois  = TRUE,
                             constraints = constraints)
  } else {
    # ...................................................
    # Data pattern is argument Y
    X <- Y
    # Realisations of binomial process
    # will be generated by runifpoint
    nY <- npoints(Y)
    Ywin <- as.owin(Y)
    Ymarx <- marks(Y)
    # expression that will be evaluated
    simexpr <- if(is.null(Ymarx)) {
      ## unmarked
      expression(runifpoint(nY, Ywin))
    } else if(is.null(dim(Ymarx))) {
      ## single column of marks
      expression({
        A <- runifpoint(nY, Ywin);
        j <- sample(nY, npoints(A), replace=TRUE);
        A %mark% Ymarx[j]
      })
    } else {
      ## multiple columns of marks
      expression({
        A <- runifpoint(nY, Ywin);
        j <- sample(nY, npoints(A), replace=TRUE);
        A %mark% Ymarx[j, ,drop=FALSE]
      })
    }
    dont.complain.about(nY, Ywin)
    # evaluate in THIS environment
    simrecipe <- simulrecipe(type = "csr",
                             expr = simexpr,
                             envir = envir.here,
                             csr   = TRUE,
                             pois  = TRUE,
                             constraints = "with fixed number of points")
  }
  
  envelopeEngine(X=X, fun=fun, simul=simrecipe,
                 nsim=nsim, nrank=nrank, ...,
                 funargs=funargs, funYargs=funYargs,
                 verbose=verbose, clipdata=clipdata,
                 transform=transform,
                 global=global, ginterval=ginterval, use.theory=use.theory,
                 alternative=alternative, scale=scale, clamp=clamp,
                 savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
                 VARIANCE=VARIANCE, nSD=nSD,
                 Yname=Yname,
                 maxnerr=maxnerr, rejectNA=rejectNA, silent=silent,
                 cl=cl, envir.user=envir.user, do.pwrong=do.pwrong)
}

envelope.ppm <- 
  function(Y, fun=Kest, nsim=99, nrank=1, ..., 
           funargs=list(), funYargs=funargs,
           simulate=NULL, fix.n=FALSE, fix.marks=FALSE,
           verbose=TRUE, clipdata=TRUE, 
           start=NULL,
           control=update(default.rmhcontrol(Y), nrep=nrep), nrep=1e5, 
           transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL, 
           alternative=c("two.sided", "less", "greater"),
           scale=NULL, clamp=FALSE, 
           savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
           VARIANCE=FALSE, nSD=2,
           Yname=NULL,
           maxnerr=nsim, rejectNA=FALSE, silent=FALSE,
           do.pwrong=FALSE,
           envir.simul=NULL) {
  cl <- short.deparse(sys.call())
  if(is.null(Yname)) Yname <- short.deparse(substitute(Y))
  if(is.null(fun)) fun <- Kest
  envir.user <- if(!is.null(envir.simul)) envir.simul else parent.frame()
  envir.here <- sys.frame(sys.nframe())

  # Extract data pattern X from fitted model Y
  X <- data.ppm(Y)
  
  if(is.null(simulate)) {
    # ...................................................
    # Simulated realisations of the fitted model Y
    # will be generated
    pois <- is.poisson(Y)
    csr <- is.stationary(Y) && pois
    type <- if(csr) "csr" else "rmh"
    # Set up parameters for rmh
    rmodel <- rmhmodel(Y, verbose=FALSE)
    if(is.null(start))
      start <- list(n.start=npoints(X))
    rstart <- rmhstart(start)
    rcontr <- rmhcontrol(control)
    if(fix.marks) {
      rcontr <- update(rcontr, fixall=TRUE, p=1, expand=1)
      nst <- if(is.multitype(X)) table(marks(X)) else npoints(X)
      rstart <- update(rstart, n.start=nst)
      constraints <- "with fixed number of points of each type"
    } else if(fix.n) {
      rcontr <- update(rcontr, p=1, expand=1)
      rstart <- update(rstart, n.start=X$n)
      constraints <- "with fixed number of points"
    } else constraints <- ""
    # pre-digest arguments
    rmhinfolist <- rmh(rmodel, rstart, rcontr, preponly=TRUE, verbose=FALSE)
    # expression that will be evaluated
    simexpr <- expression(rmhEngine(rmhinfolist, verbose=FALSE))
    dont.complain.about(rmhinfolist)
    # evaluate in THIS environment
    simrecipe <- simulrecipe(type  = type,
                             expr  = simexpr,
                             envir = envir.here,
                             csr   = csr,
                             pois  = pois,
                             constraints = constraints)
  } else {
    # ...................................................
    # Simulations are determined by 'simulate' argument
    # Processing is deferred to envelopeEngine
    simrecipe <- simulate
  }
  envelopeEngine(X=X, fun=fun, simul=simrecipe, 
                 nsim=nsim, nrank=nrank, ...,
                 funargs=funargs, funYargs=funYargs,
                 verbose=verbose, clipdata=clipdata,
                 transform=transform,
                 global=global, ginterval=ginterval, use.theory=use.theory,
                 alternative=alternative, scale=scale, clamp=clamp, 
                 savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
                 VARIANCE=VARIANCE, nSD=nSD,
                 Yname=Yname,
                 maxnerr=maxnerr, rejectNA=rejectNA, silent=silent,
                 cl=cl, envir.user=envir.user, do.pwrong=do.pwrong)
}

envelope.kppm <-
  function(Y, fun=Kest, nsim=99, nrank=1, ..., 
           funargs=list(), funYargs=funargs,
           simulate=NULL, verbose=TRUE, clipdata=TRUE, 
           transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
           alternative=c("two.sided", "less", "greater"),
           scale=NULL, clamp=FALSE,
           savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
           VARIANCE=FALSE, nSD=2, Yname=NULL, 
           maxnerr=nsim, rejectNA=FALSE, silent=FALSE,
           do.pwrong=FALSE, envir.simul=NULL)
{
  cl <- short.deparse(sys.call())
  if(is.null(Yname)) Yname <- short.deparse(substitute(Y))
  if(is.null(fun)) fun <- Kest
  envir.user <- if(!is.null(envir.simul)) envir.simul else parent.frame()
  envir.here <- sys.frame(sys.nframe())
  
  # Extract data pattern X from fitted model Y
  X <- Y$X
  
  if(is.null(simulate)) {
    # Simulated realisations of the fitted model Y
    # will be generated using simulate.kppm
    kmodel <- Y
    # expression that will be evaluated
    simexpr <- expression(simulate(kmodel)[[1L]])
    dont.complain.about(kmodel)
    # evaluate in THIS environment
    simrecipe <- simulrecipe(type = "kppm",
                             expr = simexpr,
                             envir = envir.here,
                             csr   = FALSE,
                             pois  = FALSE)
  } else {
    # ...................................................
    # Simulations are determined by 'simulate' argument
    # Processing is deferred to envelopeEngine
    simrecipe <- simulate
  }
  envelopeEngine(X=X, fun=fun, simul=simrecipe, 
                 nsim=nsim, nrank=nrank, ...,
                 funargs=funargs, funYargs=funYargs,
                 verbose=verbose, clipdata=clipdata,
                 transform=transform,
                 global=global, ginterval=ginterval, use.theory=use.theory,
                 alternative=alternative, scale=scale, clamp=clamp,
                 savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
                 VARIANCE=VARIANCE, nSD=nSD,
                 Yname=Yname, 
                 maxnerr=maxnerr, rejectNA=rejectNA, silent=silent,
                 cl=cl, envir.user=envir.user, do.pwrong=do.pwrong)

}

envelope.slrm <-
  function(Y, fun=Kest, nsim=99, nrank=1, ..., 
           funargs=list(), funYargs=funargs,
           simulate=NULL, verbose=TRUE, clipdata=TRUE, 
           transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
           alternative=c("two.sided", "less", "greater"),
           scale=NULL, clamp=FALSE,
           savefuns=FALSE, savepatterns=FALSE, nsim2=nsim,
           VARIANCE=FALSE, nSD=2, Yname=NULL, 
           maxnerr=nsim, rejectNA=FALSE, silent=FALSE,
           do.pwrong=FALSE, envir.simul=NULL)
{
  cl <- short.deparse(sys.call())
  if(is.null(Yname)) Yname <- short.deparse(substitute(Y))
  if(is.null(fun)) fun <- Kest
  envir.user <- if(!is.null(envir.simul)) envir.simul else parent.frame()
  envir.here <- sys.frame(sys.nframe())
  
  # Extract data pattern X from fitted model Y
  X <- response(Y)
  
  if(is.null(simulate)) {
    # Simulated realisations of the fitted model Y
    # will be generated using simulate.slrm
    smodel <- Y
    # expression that will be evaluated
    simexpr <- expression(simulate(smodel)[[1L]])
    dont.complain.about(smodel)
    # evaluate in THIS environment
    simrecipe <- simulrecipe(type = "slrm",
                             expr = simexpr,
                             envir = envir.here,
                             csr   = FALSE,
                             pois  = FALSE)
  } else {
    # ...................................................
    # Simulations are determined by 'simulate' argument
    # Processing is deferred to envelopeEngine
    simrecipe <- simulate
  }
  envelopeEngine(X=X, fun=fun, simul=simrecipe, 
                 nsim=nsim, nrank=nrank, ...,
                 funargs=funargs, funYargs=funYargs,
                 verbose=verbose, clipdata=clipdata,
                 transform=transform,
                 global=global, ginterval=ginterval, use.theory=use.theory,
                 alternative=alternative, scale=scale, clamp=clamp,
                 savefuns=savefuns, savepatterns=savepatterns, nsim2=nsim2,
                 VARIANCE=VARIANCE, nSD=nSD,
                 Yname=Yname, 
                 maxnerr=maxnerr, rejectNA=rejectNA, silent=silent,
                 cl=cl, envir.user=envir.user, do.pwrong=do.pwrong)

}

## .................................................................
##   Engine for simulating and computing envelopes
## .................................................................
#
#  X is the data point pattern, which could be ppp, pp3, ppx etc
#  X determines the class of pattern expected from the simulations
#

envelopeEngine <-
  function(X, fun, simul,
           nsim=99, nrank=1, ..., funargs=list(), funYargs=funargs,
           verbose=TRUE, clipdata=TRUE, 
           transform=NULL, global=FALSE, ginterval=NULL, use.theory=NULL,
           alternative=c("two.sided", "less", "greater"),
           scale=NULL, clamp=FALSE,
           savefuns=FALSE, savepatterns=FALSE,
           saveresultof=NULL,
           weights=NULL,
           nsim2=nsim,
           VARIANCE=FALSE, nSD=2,
           Yname=NULL,
           maxnerr=nsim,
           rejectNA=FALSE,
           silent=FALSE,
           maxerr.action=c("fatal", "warn", "null"),
           internal=NULL, cl=NULL,
           envir.user=envir.user,
           expected.arg="r",
           do.pwrong=FALSE,
           foreignclass=NULL,
           collectrubbish=FALSE)
{
  #
  envir.here <- sys.frame(sys.nframe())

  alternative <- match.arg(alternative)
  maxerr.action <- match.arg(maxerr.action)  

  foreignclass <- as.character(foreignclass)
  if(length(foreignclass) != 0 && clipdata) {
    warning(paste("Ignoring clipdata=TRUE:",
                  "I don't know how to clip objects of class",
                  sQuote(paste(foreignclass, collapse=","))))
    clipdata <- FALSE
  }
  
  # ----------------------------------------------------------
  # Determine Simulation
  # ----------------------------------------------------------

  check.1.integer(nsim)
  stopifnot(nsim > 1)

  # Identify class of patterns to be simulated, from data pattern X
  Xclass <- if(is.ppp(X)) "ppp" else
            if(is.pp3(X)) "pp3" else
            if(is.ppx(X)) "ppx" else
            if(inherits(X, foreignclass)) foreignclass else
            stop("Unrecognised class of point pattern")
  Xobjectname <- paste("point pattern of class", sQuote(Xclass))

  # Option to use weighted average
  if(use.weights <- !is.null(weights)) {
    # weight can be either a numeric vector or a function
    if(is.numeric(weights)) {
      compute.weights <- FALSE
      weightfun <- NULL
    } else if(is.function(weights)) {
      compute.weights <- TRUE
      weightfun <- weights
      weights <- NULL  
    } else stop("weights should be either a function or a numeric vector")
  } else compute.weights <- FALSE
    
  # Undocumented option to generate patterns only.
  patterns.only <- identical(internal$eject, "patterns")

  # Undocumented option to evaluate 'something' for each pattern
  if(savevalues <- !is.null(saveresultof)) {
    stopifnot(is.function(saveresultof))
    SavedValues <- list()
  }

  # Identify type of simulation from argument 'simul'
  if(inherits(simul, "simulrecipe")) {
    # ..................................................
    # simulation recipe is given
    simtype <- simul$type
    simexpr <- simul$expr
    envir   <- simul$envir
    csr     <- simul$csr
    pois    <- simul$pois
    constraints <- simul$constraints
  } else {
    # ...................................................
    # simulation is specified by argument `simulate' to envelope()
    simulate <- simul
    # which should be an expression, or a list of point patterns,
    # or an envelope object, or a function to be applied to the data
    csr <- FALSE
    # override
    if(!is.null(icsr <- internal$csr)) csr <- icsr
    pois <- csr
    constraints <- ""
#    model <- NULL
    if(inherits(simulate, "envelope")) {
      # envelope object: see if it contains stored point patterns
      simpat <- attr(simulate, "simpatterns")
      if(!is.null(simpat))
        simulate <- simpat
      else
        stop(paste("The argument", sQuote("simulate"),
                   "is an envelope object but does not contain",
                   "any saved point patterns."))
    }
    if(is.expression(simulate)) {
      ## The user-supplied expression 'simulate' will be evaluated repeatedly
      simtype <- "expr"
      simexpr <- simulate
      envir <- envir.user
    } else if(is.function(simulate)) {
      ## User-supplied function 'simulate' will be repeatedly evaluated on X
      simtype <- "func"
      simexpr <- expression(simulate(X))
      envir <- envir.here
    } else if(is.list(simulate) &&
              all(sapply(simulate, inherits, what=Xclass))) {
      # The user-supplied list of point patterns will be used
      simtype <- "list"
      SimDataList <- simulate
      # expression that will be evaluated
      simexpr <- expression(SimDataList[[i+nerr]])
      dont.complain.about(SimDataList)
      envir <- envir.here
      # ensure that `i' is defined
      i <- 1L
      nerr <- 0L
      maxnerr <- min(length(SimDataList)-nsim, maxnerr)
      # any messages?
      if(!is.null(mess <- attr(simulate, "internal"))) {
        # determine whether these point patterns are realisations of CSR
        csr <- !is.null(mc <- mess$csr) && mc
      }
    } else stop(paste(sQuote("simulate"),
                      "should be an expression,",
                      "or a list of point patterns of the same kind as X"))
  }
  # -------------------------------------------------------------------
  # Determine clipping window
  # ------------------------------------------------------------------

  if(clipdata) {
    # Generate one realisation
    Xsim <- eval(simexpr, envir=envir)
    if(!inherits(Xsim, Xclass))
      switch(simtype,
             csr=stop(paste("Internal error:", Xobjectname, "not generated")),
             rmh=stop(paste("Internal error: rmh did not return an",
               Xobjectname)),
             kppm=stop(paste("Internal error: simulate.kppm did not return an",
               Xobjectname)),
             slrm=stop(paste("Internal error: simulate.slrm did not return an",
               Xobjectname)),
             expr=stop(paste("Evaluating the expression", sQuote("simulate"),
               "did not yield an", Xobjectname)),
             func=stop(paste("Evaluating the function", sQuote("simulate"),
               "did not yield an", Xobjectname)),
             list=stop(paste("Internal error: list entry was not an",
               Xobjectname)),
             stop(paste("Internal error:", Xobjectname, "not generated"))
             )
    # Extract window
    clipwin <- Xsim$window
    if(!is.subset.owin(clipwin, X$window))
      warning("Window containing simulated patterns is not a subset of data window")
  }
  
  # ------------------------------------------------------------------
  # Summary function to be applied 
  # ------------------------------------------------------------------

  if(is.null(fun))
    stop("Internal error: fun is NULL")

  # Name of function, for error messages
  fname <- if(is.name(substitute(fun))) short.deparse(substitute(fun)) else
  if(is.character(fun)) fun else "fun"
  fname <- sQuote(fname)

  # R function to apply
  if(is.character(fun)) {
    gotfun <- try(get(fun, mode="function"))
    if(inherits(gotfun, "try-error"))
      stop(paste("Could not find a function named", sQuote(fun)))
    fun <- gotfun
  } else if(!is.function(fun)) 
    stop(paste("unrecognised format for function", fname))
  fargs <- names(formals(fun))
  if(!any(c(expected.arg, "...") %in% fargs))
    stop(paste(fname, "should have",
               ngettext(length(expected.arg), "an argument", "arguments"),
               "named", commasep(sQuote(expected.arg)),
               "or a", sQuote("..."), "argument"))
  usecorrection <- any(c("correction", "...") %in% fargs)
  
  # ---------------------------------------------------------------------
  # validate other arguments
  if((nrank %% 1) != 0)
    stop("nrank must be an integer")
  if((nsim %% 1) != 0)
    stop("nsim must be an integer")
  stopifnot(nrank > 0 && nrank < nsim/2)

  rgiven <- any(expected.arg %in% names(list(...)))

  if(tran <- !is.null(transform)) {
    stopifnot(is.expression(transform))
    # prepare expressions to be evaluated each time 
    transform.funX    <- inject.expr("with(funX,.)",    transform)
    transform.funXsim <- inject.expr("with(funXsim,.)", transform)
    # .... old code using 'eval.fv' ......
    # transform.funX <- dotexpr.to.call(transform, "funX", "eval.fv")
    # transform.funXsim <- dotexpr.to.call(transform, "funXsim", "eval.fv")
    # 'transform.funX' and 'transform.funXsim' are unevaluated calls to eval.fv
  }
  if(!is.null(ginterval)) 
    stopifnot(is.numeric(ginterval) && length(ginterval) == 2)
    
  # ---------------------------------------------------------------------
  # Evaluate function for data pattern X
  # ------------------------------------------------------------------
  Xarg <- if(!clipdata) X else X[clipwin]
  corrx <- if(usecorrection) list(correction="best") else NULL
  dont.complain.about(Xarg)
  funX <- do.call(fun,
                  resolve.defaults(list(quote(Xarg)),
                                   list(...),
                                   funYargs,
                                   corrx))
                                     
  if(!inherits(funX, "fv"))
    stop(paste("The function", fname,
               "must return an object of class", sQuote("fv")))

  ## catch 'conservation' parameters
  conserveargs <- attr(funX, "conserve")
  if(!is.null(conserveargs) && !any(c("conserve", "...") %in% fargs))
    stop(paste("In this usage, the function", fname,
               "should have an argument named 'conserve' or '...'"))

  ## warn about 'dangerous' arguments
  if(!is.null(dang <- attr(funX, "dangerous")) &&
     any(uhoh <- dang %in% names(list(...)))) {
    nuh <- sum(uhoh)
    warning(paste("Envelope may be invalid;",
                  ngettext(nuh, "argument", "arguments"),
                  commasep(sQuote(dang[uhoh])),
                  ngettext(nuh, "appears", "appear"),
                  "to have been fixed."),
            call.=FALSE)
  }
  
  argname <- fvnames(funX, ".x")
  valname <- fvnames(funX, ".y")
  has.theo <- "theo" %in% fvnames(funX, "*")
  csr.theo <- csr && has.theo
  use.theory <- if(is.null(use.theory)) csr.theo else (use.theory && has.theo)
  
  if(tran) {
    # extract only the recommended value
    if(use.theory) 
      funX <- funX[, c(argname, valname, "theo")]
    else
      funX <- funX[, c(argname, valname)]
    # apply the transformation to it
    funX <- eval(transform.funX)
  }
    
  rvals <- funX[[argname]]
#  fX    <- funX[[valname]]

  # default domain over which to maximise
  alim <- attr(funX, "alim")
  if(global && is.null(ginterval))
    ginterval <- if(rgiven || is.null(alim)) range(rvals) else alim
  
  #--------------------------------------------------------------------
  # Determine number of simulations
  # ------------------------------------------------------------------
  #
  ## determine whether dual simulations are required
  ## (one set of simulations to calculate the theoretical mean,
  ##  another independent set of simulations to obtain the critical point.)
  dual <- (global && !use.theory && !VARIANCE)
  if(dual) {
    check.1.integer(nsim2)
    stopifnot(nsim2 >= 1)
  }
  Nsim <- if(!dual) nsim else (nsim + nsim2)

  # if taking data from a list of point patterns,
  # check there are enough of them
  if(simtype == "list" && Nsim > length(SimDataList))
    stop(paste("Number of simulations",
               paren(if(!dual)
                     paste(nsim) else
                     paste(nsim, "+", nsim2, "=", Nsim)
                     ),
               "exceeds number of point pattern datasets supplied"))

  # Undocumented secret exit
  # ------------------------------------------------------------------
  if(patterns.only) {
    # generate simulated realisations and return only these patterns
    if(verbose) {
      action <- if(simtype == "list") "Extracting" else "Generating"
      descrip <- switch(simtype,
                        csr = "simulations of CSR",
                        rmh = paste("simulated realisations of fitted",
                          if(pois) "Poisson" else "Gibbs",
                          "model"),
                        kppm = "simulated realisations of fitted cluster model",
                        slrm = "simulated realisations of spatial logistic regression model",
                        expr = "simulations by evaluating expression",
                        func = "simulations by evaluating function",
                        list = "point patterns from list",
                        "simulated realisations")
      if(!is.null(constraints) && nzchar(constraints))
        descrip <- paste(descrip, constraints)
      explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
                                     nsim, "to calculate envelopes")) else ""
      splat(action, Nsim, descrip, explan, "...")
    }
    XsimList <- list()
  # start simulation loop
    sstate <- list() 
    for(i in 1:Nsim) {
      if(verbose) sstate <- progressreport(i, Nsim, state=sstate)
      Xsim <- eval(simexpr, envir=envir)
      if(!inherits(Xsim, Xclass))
        switch(simtype,
               csr={
                 stop(paste("Internal error:", Xobjectname, "not generated"))
               },
               rmh={
                 stop(paste("Internal error: rmh did not return an",
                            Xobjectname))
               },
               kppm={
                 stop(paste("Internal error: simulate.kppm did not return an",
                            Xobjectname))
               },
               slrm={
                 stop(paste("Internal error: simulate.slrm did not return an",
                            Xobjectname))
               },
               expr={
                 stop(paste("Evaluating the expression", sQuote("simulate"),
                            "did not yield an", Xobjectname))
               },
               func={
                 stop(paste("Evaluating the function", sQuote("simulate"),
                            "did not yield an", Xobjectname))
               },
               list={
                 stop(paste("Internal error: list entry was not an",
                            Xobjectname))
               },
               stop(paste("Internal error:", Xobjectname, "not generated"))
               )
      XsimList[[i]] <- Xsim
    }
    if(verbose) {
      cat(paste("Done.\n"))
      flush.console()
    }
    attr(XsimList, "internal") <- list(csr=csr)
    return(XsimList)
  }
  
  # capture main decision parameters
  envelopeInfo <-  list(call=cl,
                        Yname=Yname,
                        valname=valname,
                        csr=csr,
                        csr.theo=csr.theo,
                        use.theory=use.theory,
                        pois=pois,
                        simtype=simtype,
                        constraints=constraints,
                        nrank=nrank,
                        nsim=nsim,
                        Nsim=Nsim,
                        global=global,
                        ginterval=ginterval,
                        dual=dual,
                        nsim2=nsim2,
                        VARIANCE=VARIANCE,
                        nSD=nSD,
                        alternative=alternative,
                        scale=scale,
                        clamp=clamp,
                        use.weights=use.weights,
                        do.pwrong=do.pwrong)

  # ----------------------------------------
  ######### SIMULATE #######################
  # ----------------------------------------

  if(verbose) {
    action <- if(simtype == "list") "Extracting" else "Generating"
    descrip <- switch(simtype,
                      csr = "simulations of CSR",
                      rmh = paste("simulated realisations of fitted",
                        if(pois) "Poisson" else "Gibbs",
                        "model"),
                      kppm = "simulated realisations of fitted cluster model",
                      slrm = "simulated realisations of fitted spatial logistic regression model",
                      expr = "simulations by evaluating expression",
                      func = "simulations by evaluating function",
                      list = "point patterns from list",
                      "simulated patterns")
    if(!is.null(constraints) && nzchar(constraints))
      descrip <- paste(descrip, constraints)
    explan <- if(dual) paren(paste(nsim2, "to estimate the mean and",
                                   nsim, "to calculate envelopes")) else ""
    splat(action, Nsim, descrip, explan, "...")
  }
  # determine whether simulated point patterns should be saved
  catchpatterns <- savepatterns && simtype != "list"
  Caughtpatterns <- list()
  # allocate space for computed function values
  nrvals <- length(rvals)
  simvals <- matrix(, nrow=nrvals, ncol=Nsim)
  # allocate space for weights to be computed
  if(compute.weights)
    weights <- numeric(Nsim)
  
  # inferred values of function argument 'r' or equivalent parameters
  if(identical(expected.arg, "r")) {
    # Kest, etc
    inferred.r.args <- list(r=rvals)
  } else if(identical(expected.arg, c("rmax", "nrval"))) {
    # K3est, etc
    inferred.r.args <- list(rmax=max(rvals), nrval=length(rvals))
  } else
  stop(paste("Don't know how to infer values of", commasep(expected.arg)))
    
  # arguments for function when applied to simulated patterns
  funargs <-
    resolve.defaults(funargs,
                     inferred.r.args,
                     list(...),
                     conserveargs,
                     if(usecorrection) list(correction="best") else NULL)

  # reject simulated pattern if function values are all NA (etc)
  rejectNA <- isTRUE(rejectNA)
  
  # start simulation loop
  nerr <- 0
  gaveup <- FALSE
  if(verbose) pstate <- list()
  for(i in 1:Nsim) {
    ## safely generate a random pattern and apply function
    success <- FALSE
    while(!success && !gaveup) {
      Xsim <- eval(simexpr, envir=envir)
      ## check valid point pattern
      if(!inherits(Xsim, Xclass))
        switch(simtype,
               csr=stop(paste("Internal error:", Xobjectname, "not generated")),
               rmh=stop(paste("Internal error: rmh did not return an",
                 Xobjectname)),
               kppm=stop(paste("Internal error:",
                 "simulate.kppm did not return an",
                 Xobjectname)),
               slrm=stop(paste("Internal error:",
                 "simulate.slrm did not return an",
                 Xobjectname)),
               expr=stop(paste("Evaluating the expression", sQuote("simulate"),
                 "did not yield an", Xobjectname)),
               func=stop(paste("Evaluating the function", sQuote("simulate"),
                 "did not yield an", Xobjectname)),
               list=stop(paste("Internal error: list entry was not an",
                 Xobjectname)),
               stop(paste("Internal error:", Xobjectname, "not generated"))
               )
      if(catchpatterns)
        Caughtpatterns[[i]] <- Xsim
      if(savevalues)
        SavedValues[[i]] <- saveresultof(Xsim)
      if(compute.weights) {
        wti <- weightfun(Xsim)
        if(!is.numeric(wti))
          stop("weightfun did not return a numeric value")
        if(length(wti) != 1L)
          stop("weightfun should return a single numeric value")
        weights[i] <- wti
      }
      ## apply function safely
      funXsim <- try(do.call(fun, c(list(Xsim), funargs)), silent=silent)

      success <-
        !inherits(funXsim, "try-error") &&
        inherits(funXsim, "fv") &&
        (!rejectNA || any(is.finite(funXsim[[valname]])))

      if(!success) {
        #' error in computing summary function
        nerr <- nerr + 1L 
        if(nerr > maxnerr) {
          gaveup <- TRUE
          errtype <- if(rejectNA) "fatal errors or NA function values"
          if(simtype == "list") {
            whinge <- paste("Exceeded maximum possible number of errors",
                          "when evaluating summary function:",
                          length(SimDataList), "patterns provided,",
                          nsim, "patterns required,",
                          nerr, ngettext(nerr, "pattern", "pattern"),
                          "rejected due to", errtype)
          } else {
            whinge <- paste("Exceeded maximum permissible number of",
                            errtype,
                            paren(paste("maxnerr =", maxnerr)),
                            "when evaluating summary function",
                            "for simulated point patterns")
          }
          switch(maxerr.action,
                 fatal = stop(whinge, call.=FALSE),
                 warn  = warning(whinge, call.=FALSE),
                 null  = {})
        } else if(!silent) cat("[retrying]\n")
      }

      #' ..... end of while(!success) ................
    }
    
    if(gaveup) break; # exit loop now
    
    ## sanity checks
    if(i == 1L) {
      if(!inherits(funXsim, "fv"))
        stop(paste("When applied to a simulated pattern, the function",
                   fname, "did not return an object of class",
                   sQuote("fv")))
      argname.sim <- fvnames(funXsim, ".x")
      valname.sim <- fvnames(funXsim, ".y")
      if(argname.sim != argname)
        stop(paste("The objects returned by", fname,
                   "when applied to a simulated pattern",
                   "and to the data pattern",
                   "are incompatible. They have different argument names",
                   sQuote(argname.sim), "and", sQuote(argname), 
                   "respectively"))
      if(valname.sim != valname)
        stop(paste("When", fname, "is applied to a simulated pattern",
                   "it provides an estimate named", sQuote(valname.sim), 
                   "whereas the estimate for the data pattern is named",
                   sQuote(valname),
                   ". Try using the argument", sQuote("correction"),
                   "to make them compatible"))
      rfunX    <- with(funX,    ".x")
      rfunXsim <- with(funXsim, ".x")
      if(!identical(rfunX, rfunXsim))
        stop(paste("When", fname, "is applied to a simulated pattern,",
                   "the values of the argument", sQuote(argname.sim),
                   "are different from those used for the data."))
    }

    if(tran) {
      # extract only the recommended value
      if(use.theory) 
        funXsim <- funXsim[, c(argname, valname, "theo")]
      else
        funXsim <- funXsim[, c(argname, valname)]
      # apply the transformation to it
      funXsim <- eval(transform.funXsim)
    }

    # extract the values for simulation i
    simvals.i <- funXsim[[valname]]
    if(length(simvals.i) != nrvals)
      stop("Vectors of function values have incompatible lengths")
      
    simvals[ , i] <- funXsim[[valname]]
    if(verbose)
      pstate <- progressreport(i, Nsim, state=pstate)
    
    if(collectrubbish) {
      rm(Xsim)
      rm(funXsim)
      gc()
    }
  }
  ##  end simulation loop
  
  if(verbose) {
    cat("\nDone.\n")
    flush.console()
  }

  # ...........................................................
  # save functions and/or patterns if so commanded

  if(!gaveup) {
    if(savefuns) {
      alldata <- cbind(rvals, simvals)
      simnames <- paste("sim", 1:Nsim, sep="")
      colnames(alldata) <- c("r", simnames)
      alldata <- as.data.frame(alldata)
      SimFuns <- fv(alldata,
                    argu="r",
                    ylab=attr(funX, "ylab"),
                    valu="sim1",
                    fmla= deparse(. ~ r),
                    alim=attr(funX, "alim"),
                    labl=names(alldata),
                    desc=c("distance argument r",
                           paste("Simulation ", 1:Nsim, sep="")),
                    fname=attr(funX, "fname"),
                    yexp=attr(funX, "yexp"),
                    unitname=unitname(funX))
      fvnames(SimFuns, ".") <- simnames
    } 
    if(savepatterns)
      SimPats <- if(simtype == "list") SimDataList else Caughtpatterns
  }
  
  ######### COMPUTE ENVELOPES #######################

  etype <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
  if(dual) {
    jsim <- 1:nsim
    jsim.mean <- nsim + 1:nsim2
  } else {
    jsim <- jsim.mean <- NULL
  }

  result <- envelope.matrix(simvals, funX=funX,
                            jsim=jsim, jsim.mean=jsim.mean,
                            type=etype, alternative=alternative,
                            scale=scale, clamp=clamp,
                            csr=csr, use.theory=use.theory,
                            nrank=nrank, ginterval=ginterval, nSD=nSD,
                            Yname=Yname, do.pwrong=do.pwrong,
                            weights=weights, gaveup=gaveup)

  ## tack on envelope information
  attr(result, "einfo") <- resolve.defaults(envelopeInfo,
                                            attr(result, "einfo"))

  if(!gaveup) {
    ## tack on functions and/or patterns if so commanded   
    if(savefuns)
      attr(result, "simfuns") <- SimFuns
    if(savepatterns) {
      attr(result, "simpatterns") <- SimPats
      attr(result, "datapattern") <- X
    }
    ## undocumented - tack on values of some other quantity
    if(savevalues) {
      attr(result, "simvalues") <- SavedValues
      attr(result, "datavalue") <- saveresultof(X)
    }
  }

  ## save function weights 
  if(use.weights)
    attr(result, "weights") <- weights
  return(result)
}


plot.envelope <- function(x, ..., main) {
  if(missing(main)) main <- short.deparse(substitute(x))
  shade.given <- ("shade" %in% names(list(...)))
  shade.implied <- !is.null(fvnames(x, ".s"))
  if(!(shade.given || shade.implied)) {
    # ensure x has default 'shade' attribute
    # (in case x was produced by an older version of spatstat)
    if(all(c("lo", "hi") %in% colnames(x)))
      fvnames(x, ".s") <- c("lo", "hi")
    else warning("Unable to determine shading for envelope")
  }
  NextMethod("plot", main=main)
}

print.envelope <- function(x, ...) {
  e <- attr(x, "einfo")
  g <- e$global
  csr <- e$csr
  pois <- e$pois
  if(is.null(pois)) pois <- csr
  simtype <- e$simtype
  constraints <- e$constraints
  nr <- e$nrank
  nsim <- e$nsim
  V <- e$VARIANCE
  fname <- flat.deparse(attr(x, "ylab"))
  type <- if(V) paste("Pointwise", e$nSD, "sigma") else
          if(g) "Simultaneous" else "Pointwise"
  splat(type, "critical envelopes for", fname,
        "\nand observed value for", sQuote(e$Yname))
  if(!is.null(valname <- e$valname) && waxlyrical('extras'))
    splat("Edge correction:", dQuote(valname))
  ## determine *actual* type of simulation
  descrip <-
    if(csr) "simulations of CSR"
    else if(!is.null(simtype)) {
      switch(simtype,
             csr="simulations of CSR",
             rmh=paste("simulations of fitted",
               if(pois) "Poisson" else "Gibbs",
               "model"),
             kppm="simulations of fitted cluster model",
             slrm="simulations of fitted spatial logistic regression model",
             expr="evaluations of user-supplied expression",
             func="evaluations of user-supplied function",
             list="point pattern datasets in user-supplied list",
             funs="columns of user-supplied data")
    } else "simulations of fitted model"
  if(!is.null(constraints) && nzchar(constraints))
    descrip <- paste(descrip, constraints)
  #  
  splat("Obtained from", nsim, descrip)
  #
  if(waxlyrical('extras')) {
    dual <- isTRUE(e$dual)
    usetheory <- isTRUE(e$use.theory)
    hownull <- if(usetheory) {
                 "(known exactly)"
               } else if(dual) {
                 paste("(estimated from a separate set of",
                       e$nsim2, "simulations)")
               } else NULL
    formodel <- if(csr) "for CSR" else NULL
    if(g) {
      splat("Envelope based on maximum deviation of", fname,
            "from null value", formodel, hownull)
    } else if(dual) {
      splat("Null value of", fname, formodel, hownull)
    }
    if(!is.null(attr(x, "simfuns"))) 
      splat("(All simulated function values are stored)")
    if(!is.null(attr(x, "simpatterns"))) 
      splat("(All simulated point patterns are stored)")
  }
  splat("Alternative:", e$alternative)
  if(!V && waxlyrical('extras')) {
    ## significance interpretation!
    alpha <- if(g) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
    splat("Significance level of",
          if(g) "simultaneous" else "pointwise",
          "Monte Carlo test:",
          paste0(if(g) nr else 2 * nr,
                 "/", nsim+1),
          "=", signif(alpha, 3))
  }
  if(waxlyrical('gory') && !is.null(pwrong <- attr(x, "pwrong"))) {
    splat("\t[Estimated significance level of pointwise excursions:",
          paste0("pwrong=", signif(pwrong, 3), "]"))
  }
  NextMethod("print")
}
                  
summary.envelope <- function(object, ...) {
  e <- attr(object, "einfo")
  g <- e$global
  V <- e$VARIANCE
  nr <- e$nrank
  nsim <- e$nsim
  csr <- e$csr
  pois <- e$pois
  if(is.null(pois)) pois <- csr
  simtype <- e$simtype
  constraints <- e$constraints
  alternative <- e$alternative
  use.theory <- e$use.theory
  has.theo <- "theo" %in% fvnames(object, "*")
  csr.theo <- csr && has.theo
  use.theory <- if(is.null(use.theory)) csr.theo else (use.theory && has.theo)
  fname <- deparse(attr(object, "ylab"))
  type <- if(V) paste("Pointwise", e$nSD, "sigma") else
          if(g) "Simultaneous" else "Pointwise"
  splat(type, "critical envelopes for", fname, 
      "\nand observed value for", sQuote(e$Yname))
  # determine *actual* type of simulation
  descrip <-
    if(csr) "simulations of CSR"
    else if(!is.null(simtype)) {
      switch(simtype,
             csr="simulations of CSR",
             rmh=paste("simulations of fitted",
               if(pois) "Poisson" else "Gibbs",
               "model"),
             kppm="simulations of fitted cluster model",
             slrm="simulations of fitted spatial logistic regression model",
             expr="evaluations of user-supplied expression",
             func="evaluations of user-supplied function",
             list="point pattern datasets in user-supplied list",
             funs="columns of user-supplied data",
             "simulated point patterns")
    } else "simulations of fitted model"
  if(!is.null(constraints) && nzchar(constraints))
    descrip <- paste(descrip, constraints)
  #  
  splat("Obtained from", nsim, descrip)
  #
  if(waxlyrical('extras')) {
    if(!is.null(e$dual) && e$dual) 
      splat("Theoretical (i.e. null) mean value of", fname,
            "estimated from a separate set of",
            e$nsim2, "simulations")
    if(!is.null(attr(object, "simfuns")))
      splat("(All", nsim, "simulated function values",
            "are stored in attr(,", dQuote("simfuns"), ") )")
    if(!is.null(attr(object, "simpatterns")))
      splat("(All", nsim, "simulated point patterns",
            "are stored in attr(,", dQuote("simpatterns"), ") )")
  }
  #
  splat("Alternative:", alternative)
  if(V) {
    # nSD envelopes
    splat(switch(alternative,
                 two.sided = "Envelopes",
                 "Critical boundary"),
          "computed as sample mean",
          switch(alternative,
                 two.sided="plus/minus",
                 less="minus",
                 greater="plus"),
          e$nSD, "sample standard deviations")
  } else {
    # critical envelopes
    lo.ord <- if(nr == 1L) "minimum" else paste(ordinal(nr), "smallest")
    hi.ord <- if(nr == 1L) "maximum" else paste(ordinal(nr), "largest")
    if(g) 
      splat(switch(alternative,
                   two.sided = "Envelopes",
                   "Critical boundary"),
            "computed as",
            if(use.theory) "theoretical curve" else "mean of simulations",
            switch(alternative,
                   two.sided="plus/minus",
                   less="minus",
                   greater="plus"),
            hi.ord,
            "simulated value of maximum", 
            switch(alternative,
                   two.sided="absolute",
                   less="negative",
                   greater="positive"),
            "deviation")
    else {
      if(alternative != "less")
        splat("Upper envelope: pointwise", hi.ord, "of simulated curves")
      if(alternative != "greater")
        splat("Lower envelope: pointwise", lo.ord, "of simulated curves")
    }
    symmetric <- (alternative == "two.sided") && !g
    alpha <- if(!symmetric) { nr/(nsim+1) } else { 2 * nr/(nsim+1) }
    splat("Significance level of Monte Carlo test:",
          paste0(if(!symmetric) nr else 2 * nr,
                 "/", nsim+1),
          "=", alpha)
  } 
  splat("Data:", e$Yname)
  return(invisible(NULL))
}
  

# envelope.matrix

# core functionality to compute envelope values

# theory = funX[["theo"]]
# observed = fX

envelope.matrix <- function(Y, ...,
                            rvals=NULL, observed=NULL, theory=NULL, 
                            funX=NULL,
                            nsim=NULL, nsim2=NULL,
                            jsim=NULL, jsim.mean=NULL,
                            type=c("pointwise", "global", "variance"),
                            alternative=c("two.sided", "less", "greater"),
                            scale = NULL, clamp=FALSE,
                            csr=FALSE, use.theory = csr, 
                            nrank=1, ginterval=NULL, nSD=2,
                            savefuns=FALSE,
                            check=TRUE,
                            Yname=NULL,
                            do.pwrong=FALSE,
                            weights=NULL,
                            precomputed=NULL,
                            gaveup=FALSE) {
  if(is.null(Yname))
    Yname <- short.deparse(substitute(Y))

  type <- match.arg(type)
  alternative <- match.arg(alternative)

  if(!is.null(funX))
    stopifnot(is.fv(funX))

  pwrong <- NULL
  use.weights <- !is.null(weights)
  cheat <- !is.null(precomputed)

  if(is.null(rvals) && is.null(observed) && !is.null(funX)) {
    ## assume funX is summary function for observed data
    rvals <- with(funX, .x)
    observed <- with(funX, .y)
    theory <- if(use.theory) (theory %orifnull% funX[["theo"]]) else NULL
    if(check) stopifnot(nrow(funX) == nrow(Y)) 
  } else if(check) {
    ## validate vectors of data
    if(is.null(rvals)) stop("rvals must be supplied")
    if(is.null(observed)) stop("observed must be supplied")
    stopifnot(length(rvals) == nrow(Y))
    stopifnot(length(observed) == length(rvals))
  }

  use.theory <- use.theory && !is.null(theory)
  if(use.theory && check) stopifnot(length(theory) == length(rvals))

  simvals <- Y
  fX <- observed

  atr <- if(!is.null(funX)) attributes(funX) else
         list(alim=range(rvals),
              ylab=quote(f(r)),
              yexp=quote(f(r)),
              fname="f")

  fname <- atr$fname

  NAvector <- rep(NA_real_, length(rvals))
  
  if(!cheat) {
    ## ................   standard calculation .....................
    ## validate weights
    if(use.weights && !gaveup) 
      check.nvector(weights, ncol(simvals), 
                    things="simulated functions", naok=TRUE, vname="weights")

    ## determine numbers of columns used
    Ncol <- if(!gaveup) ncol(simvals) else Inf
    if(Ncol < 2)
      stop("Need at least 2 columns of function values")

    ## all columns are used unless 'nsim' or 'jsim' given.
    if(!(is.null(nsim) && is.null(jsim))) {
      if(is.null(jsim)) {
        jsim <- 1:nsim
      } else if(is.null(nsim)) {
        nsim <- length(jsim)
      } else stopifnot(length(jsim) == nsim)
      if(nsim > Ncol)
        stop(paste(nsim, "simulations are not available; only",
                   Ncol, "columns provided"))
    }
    
    ## nsim2 or jsim.mean may be given, and imply dual calculation
    if(!(is.null(nsim2) && is.null(jsim.mean))) {
      if(is.null(jsim.mean)) {
        jsim.mean <- setdiff(seq_len(Ncol), jsim)[1:nsim2]
      } else if(is.null(nsim2)) {
        nsim2 <- length(jsim.mean)
      } else stopifnot(length(jsim.mean) == nsim2)
      if(nsim + nsim2 > Ncol)
        stop(paste(nsim, "+", nsim2, "=", nsim+nsim2, 
                   "simulations are not available; only",
                   Ncol, "columns provided"))
      if(length(intersect(jsim, jsim.mean)))
        warning("Internal warning: Indices in jsim and jsim.mean overlap")
    }
      
    restrict.columns <- !is.null(jsim)
    dual <- !is.null(jsim.mean)

  } else {
    ## ................ precomputed values ..................
    ## validate weights
    if(use.weights) 
      check.nvector(weights, nsim,
                    things="simulations", naok=TRUE, vname="weights")
    restrict.columns <- FALSE
    dual <- FALSE
  }

  shadenames <- NULL
  nsim.mean <- NULL
  
  switch(type,
         pointwise = {
           ## ....... POINTWISE ENVELOPES ...............................
           if(gaveup) {
             lo <- hi <- NAvector
           } else if(cheat) {
             stopifnot(checkfields(precomputed, c("lo", "hi")))
             lo <- precomputed$lo
             hi <- precomputed$hi
           } else {
             simvals[is.infinite(simvals)] <- NA
             if(restrict.columns) {
               simvals <- simvals[,jsim]
               if(use.weights) weights <- weights[jsim]
             }
             nsim <- ncol(simvals)
             if(nrank == 1L) {
               lohi <- apply(simvals, 1L, range)
             } else {
               lohi <- apply(simvals, 1L,
#                             function(x, n) { sort(x)[n] },
                             orderstats,
                             k=c(nrank, nsim-nrank+1L))
             }
             lo <- lohi[1L,]
             hi <- lohi[2L,]
           }
           lo.name <- "lower pointwise envelope of %s from simulations"
           hi.name <- "upper pointwise envelope of %s from simulations"
           ##
           if(!gaveup)
             switch(alternative,
                    two.sided = { },
                    less = {
                      hi <- rep.int(Inf, length(hi))
                      hi.name <- "infinite upper limit"
                    },
                    greater = {
                      lo <- rep.int(-Inf, length(lo))
                      lo.name <- "infinite lower limit"
                    })
           if(use.theory) {
             results <- data.frame(r=rvals,
                                   obs=fX,
                                   theo=theory,
                                   lo=lo,
                                   hi=hi)
           } else {
             m <- if(gaveup) NAvector else
                  if(cheat) precomputed$mmean else 
                  if(!use.weights) apply(simvals, 1L, mean, na.rm=TRUE) else
                  apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
             results <- data.frame(r=rvals,
                                   obs=fX,
                                   mmean=m,
                                   lo=lo,
                                   hi=hi)
           }
           shadenames <- c("lo", "hi")
           if(do.pwrong) {
             ## estimate the p-value for the 'wrong test'
             if(gaveup) {
               pwrong <- NA_real_
             } else if(cheat) {
               pwrong <- precomputed$pwrong
               do.pwrong <- !is.null(pwrong) && !badprobability(pwrong, FALSE)
             } else {
               dataranks <- t(apply(simvals, 1, rank, ties.method="random"))
               upper.signif <- (dataranks <= nrank)
               lower.signif <- (dataranks >= nsim-nrank+1L)
               is.signif <- switch(alternative,
                                   less = lower.signif,
                                   greater = upper.signif,
                                   two.sided = lower.signif | upper.signif)
               is.signif.somewhere <- matcolany(is.signif)
               pwrong <- sum(is.signif.somewhere)/nsim
             }
           }
         },
         global = {
           ## ..... SIMULTANEOUS ENVELOPES ..........................
           if(gaveup) {
             lo <- hi <- reference <- NAvector
           } else if(cheat) {
             ## ... use precomputed values ..
             stopifnot(checkfields(precomputed, c("lo", "hi")))
             lo <- precomputed$lo
             hi <- precomputed$hi
             if(use.theory) {
               reference <- theory
             } else {
               stopifnot(checkfields(precomputed, "mmean"))
               reference <- precomputed$mmean
             }
             domain <- rep.int(TRUE, length(rvals))
           } else {
             ## ... normal case: compute envelopes from simulations
             if(!is.null(ginterval)) {
               domain <- (rvals >= ginterval[1L]) & (rvals <= ginterval[2L])
               funX <- funX[domain, ]
               simvals <- simvals[domain, ]
             } else domain <- rep.int(TRUE, length(rvals))
             simvals[is.infinite(simvals)] <- NA
             if(use.theory) {
               reference <- theory[domain]
               if(restrict.columns) {
                 simvals <- simvals[, jsim]
                 if(use.weights) weights <- weights[jsim]
               }
             } else if(dual) {
               # Estimate the mean from one set of columns
               # Form envelopes from another set of columns
               simvals.mean <- simvals[, jsim.mean]
               # mmean <-
               reference <- 
                 if(!use.weights) apply(simvals.mean, 1L, mean, na.rm=TRUE) else
                 apply(simvals.mean, 1L, weighted.mean, w=weights[jsim.mean],
                       na.rm=TRUE)
               nsim.mean <- ncol(simvals.mean)
               # retain only columns used for envelope
               simvals <- simvals[, jsim]
             } else {
               # Compute the mean and envelopes using the same data
               if(restrict.columns) {
                 simvals <- simvals[, jsim]
                 if(use.weights) weights <- weights[jsim]
               }
               # mmean <-
               reference <- 
                 if(!use.weights) apply(simvals, 1L, mean, na.rm=TRUE) else
                 apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
             }
             nsim <- ncol(simvals)
             # compute deviations
             deviations <- sweep(simvals, 1L, reference)
             deviations <-
               switch(alternative,
                      two.sided = abs(deviations),
                      greater = if(clamp) pmax(0, deviations) else deviations,
                      less = if(clamp) pmax(0, -deviations) else (-deviations))
             deviations <- matrix(deviations,
                                  nrow=nrow(simvals), ncol=ncol(simvals))
             ## rescale ?
             sc <- 1
             if(!is.null(scale)) {
               stopifnot(is.function(scale))
               sc <- scale(rvals)
               sname <- "scale(r)"
               ans <- check.nvector(sc, length(rvals), things="values of r",
                                    fatal=FALSE, vname=sname)
               if(!ans)
                 stop(attr(ans, "whinge"), call.=FALSE)
               if(any(bad <- (sc <= 0))) {
                 ## issue a warning unless this only happens at r=0
                 if(any(bad[rvals > 0]))
                   warning(paste("Some values of", sname,
                                 "were negative or zero:",
                                 "scale was reset to 1 for these values"),
                           call.=FALSE)
                 sc[bad] <- 1
               }
               deviations <- sweep(deviations, 1L, sc, "/")
             }
             ## compute max (scaled) deviations
             suprema <- apply(deviations, 2L, max, na.rm=TRUE)
             # ranked deviations
             dmax <- sort(suprema)[nsim-nrank+1L]
             # simultaneous bands
             lo <- reference - sc * dmax
             hi <- reference + sc * dmax
           }

           lo.name <- "lower critical boundary for %s"
           hi.name <- "upper critical boundary for %s"

           if(!gaveup)
             switch(alternative,
                    two.sided = { },
                    less = {
                      hi <- rep.int(Inf, length(hi))
                      hi.name <- "infinite upper boundary"
                    },
                    greater = {
                      lo <- rep.int(-Inf, length(lo))
                      lo.name <- "infinite lower boundary"
                    })

           if(use.theory)
             results <- data.frame(r=rvals[domain],
                                   obs=fX[domain],
                                   theo=reference,
                                   lo=lo,
                                   hi=hi)
           else
             results <- data.frame(r=rvals[domain],
                                   obs=fX[domain],
                                   mmean=reference,
                                   lo=lo,
                                   hi=hi)
           shadenames <- c("lo", "hi")
           if(do.pwrong)
             warning(paste("Argument", sQuote("do.pwrong=TRUE"), "ignored;",
                           "it is not relevant to global envelopes"))
         },
         variance={
           ## ....... POINTWISE MEAN, VARIANCE etc ......................
           if(gaveup) {
             Ef <- varf <- NAvector
           } else if(cheat) {
             # .... use precomputed values ....
             stopifnot(checkfields(precomputed, c("Ef", "varf")))
             Ef   <- precomputed$Ef
             varf <- precomputed$varf
           } else {
             ## .... normal case: compute from simulations
             simvals[is.infinite(simvals)] <- NA
             if(restrict.columns) {
               simvals <- simvals[, jsim]
               if(use.weights) weights <- weights[jsim]
             }
             nsim <- ncol(simvals)
             if(!use.weights) {
               Ef   <- apply(simvals, 1L, mean, na.rm=TRUE)
               varf <- apply(simvals, 1L, var,  na.rm=TRUE)
             } else {
               Ef   <- apply(simvals, 1L, weighted.mean, w=weights, na.rm=TRUE)
               varf <- apply(simvals, 1L, weighted.var,  w=weights, na.rm=TRUE)
             }
           }
           if(gaveup) {
             sd <- stdres <- lo <- hi <- loCI <- hiCI <- NAvector
           } else {
             ## derived quantities
             sd <- sqrt(varf)
             stdres <- (fX-Ef)/sd
             stdres[!is.finite(stdres)] <- NA
             ## critical limits
             lo <- Ef - nSD * sd
             hi <- Ef + nSD * sd
             ## confidence interval 
             loCI <- Ef - nSD * sd/sqrt(nsim)
             hiCI <- Ef + nSD * sd/sqrt(nsim)
           }
           lo.name <- paste("lower", nSD, "sigma critical limit for %s")
           hi.name <- paste("upper", nSD, "sigma critical limit for %s")
           loCI.name <- paste("lower", nSD, "sigma confidence bound",
                              "for mean of simulated %s")
           hiCI.name <- paste("upper", nSD, "sigma confidence bound",
                                "for mean of simulated %s")
           ##
           if(!gaveup)
             switch(alternative,
                    two.sided = { },
                    less = {
                      hi <- hiCI <- rep.int(Inf, length(hi))
                      hi.name <- "infinite upper boundary"
                      hiCI.name <- "infinite upper confidence limit"
                    },
                    greater = {
                      lo <- loCI <- rep.int(-Inf, length(lo))
                      lo.name <- "infinite lower boundary"
                      loCI.name <- "infinite lower confidence limit"
                    })
           ## put together
           if(use.theory) {
             results <- data.frame(r=rvals,
                                   obs=fX,
                                   theo=theory,
                                   lo=lo,
                                   hi=hi)
             shadenames <- c("lo", "hi")
             morestuff <- data.frame(mmean=Ef,
                                     var=varf,
                                     res=fX-Ef,
                                     stdres=stdres,
                                     loCI=loCI,
                                     hiCI=hiCI)
             loCIlabel <-
               if(alternative == "greater" && !gaveup) "-infinity" else
               makefvlabel(NULL, NULL, fname, "loCI")
             hiCIlabel <-
               if(alternative == "less" && !gaveup) "infinity" else 
               makefvlabel(NULL, NULL, fname, "hiCI")
             mslabl <- c(makefvlabel(NULL, "bar", fname),
                         makefvlabel("var", "hat", fname),
                         makefvlabel("res", "hat", fname),
                         makefvlabel("stdres", "hat", fname),
                         loCIlabel,
                         hiCIlabel)
             wted <- if(use.weights) "weighted " else NULL
             msdesc <- c(paste0(wted, "sample mean of %s from simulations"),
                         paste0(wted, "sample variance of %s from simulations"),
                         "raw residual",
                         "standardised residual",
                         loCI.name, hiCI.name)
           } else {
             results <- data.frame(r=rvals,
                                   obs=fX,
                                   mmean=Ef,
                                   lo=lo,
                                   hi=hi)
             shadenames <- c("lo", "hi")
             morestuff <- data.frame(var=varf,
                                     res=fX-Ef,
                                     stdres=stdres,
                                     loCI=loCI,
                                     hiCI=hiCI)
             loCIlabel <-
               if(alternative == "greater" && !gaveup) "-infinity" else
               makefvlabel(NULL, NULL, fname, "loCI")
             hiCIlabel <-
               if(alternative == "less" && !gaveup) "infinity" else 
               makefvlabel(NULL, NULL, fname, "hiCI")
             mslabl <- c(makefvlabel("var", "hat", fname),
                         makefvlabel("res", "hat", fname),
                         makefvlabel("stdres", "hat", fname),
                         loCIlabel,
                         hiCIlabel)
             msdesc <- c(paste0(if(use.weights) "weighted " else NULL,
                                "sample variance of %s from simulations"),
                         "raw residual",
                         "standardised residual",
                         loCI.name, hiCI.name)
           }
           if(do.pwrong) {
             ## estimate the p-value for the 'wrong test'
             if(gaveup) {
               pwrong <- NA_real_
             } else if(cheat) {
               pwrong <- precomputed$pwrong
               do.pwrong <- !is.null(pwrong) && !badprobability(pwrong, FALSE)
             } else {
               upper.signif <- (simvals > hi)
               lower.signif <- (simvals < lo)
               is.signif <- switch(alternative,
                                   less = lower.signif,
                                   greater = upper.signif,
                                   two.sided = lower.signif | upper.signif)
#               is.signif.somewhere <- apply(is.signif, 2, any)
               is.signif.somewhere <- matcolany(is.signif)
               pwrong <- sum(is.signif.somewhere)/nsim
             }
           }
         }
         )

  ############  WRAP UP #########################

  if(use.theory) {
    # reference is computed curve `theo'
    reflabl <- makefvlabel(NULL, NULL, fname, "theo")
    refdesc <- paste0("theoretical value of %s", if(csr) " for CSR" else NULL)
  } else {
    # reference is sample mean of simulations
    reflabl <- makefvlabel(NULL, "bar", fname)
    refdesc <- paste0(if(use.weights) "weighted " else NULL,
                      "sample mean of %s from simulations")
  }

  lolabl <- if(alternative == "greater" && !gaveup) "-infinity" else
             makefvlabel(NULL, "hat", fname, "lo")
  hilabl <- if(alternative == "less"&& !gaveup) "infinity" else
             makefvlabel(NULL, "hat", fname, "hi")

  result <- fv(results,
               argu="r",
               ylab=atr$ylab,
               valu="obs",
               fmla= deparse(. ~ r),
               alim=intersect.ranges(atr$alim, range(results$r)),
               labl=c("r",
                 makefvlabel(NULL, "hat", fname, "obs"),
                 reflabl,
                 lolabl,
                 hilabl),
               desc=c("distance argument r",
                 "observed value of %s for data pattern",
                 refdesc, lo.name, hi.name),
               fname=atr$fname,
               yexp =atr$yexp)

  # columns to be plotted by default
  dotty <- c("obs", if(use.theory) "theo" else "mmean", "hi", "lo")

  if(type == "variance") {
    # add more stuff
    result <- bind.fv(result, morestuff, mslabl, msdesc)
    if(use.theory) dotty <- c(dotty, "mmean")
  }

  fvnames(result, ".") <- dotty
  fvnames(result, ".s") <- shadenames

  unitname(result) <- unitname(funX)
  class(result) <- c("envelope", class(result))

  # tack on envelope information
  attr(result, "einfo") <- list(global = (type =="global"),
                                ginterval = ginterval,
                                alternative=alternative,
                                scale = scale,
                                clamp = clamp,
                                csr = csr,
                                use.theory = use.theory,
                                csr.theo = csr && use.theory,
                                simtype = "funs",
                                constraints = "",
                                nrank = nrank,
                                nsim = nsim,
                                VARIANCE = (type == "variance"),
                                nSD = nSD,
                                valname = NULL,
                                dual = dual,
                                nsim = nsim,
                                nsim2 = nsim.mean,
                                Yname = Yname,
                                do.pwrong=do.pwrong,
                                use.weights=use.weights,
                                gaveup = gaveup)

  # tack on saved functions
  if(savefuns && !gaveup) {
    nSim <- ncol(Y)
    alldata <- cbind(rvals, Y)
    simnames <- paste("sim", 1:nSim, sep="")
    colnames(alldata) <- c("r", simnames)
    alldata <- as.data.frame(alldata)
    SimFuns <- fv(alldata,
                  argu="r",
                  ylab=atr$ylab,
                  valu="sim1",
                  fmla= deparse(. ~ r),
                  alim=atr$alim,
                  labl=names(alldata),
                  desc=c("distance argument r",
                          paste("Simulation ", 1:nSim, sep="")),
                  unitname=unitname(funX))
    fvnames(SimFuns, ".") <- simnames
    attr(result, "simfuns") <- SimFuns
  }
  if(do.pwrong)
    attr(result, "pwrong") <- pwrong
  if(use.weights)
    attr(result, "weights") <- weights
  return(result)
}

envelope.envelope <- function(Y, fun=NULL, ...,
                              transform=NULL, global=FALSE, VARIANCE=FALSE) {

  Yname <- short.deparse(substitute(Y))

  stopifnot(inherits(Y, "envelope"))
  Yorig <- Y

  aargh <- list(...)

  X  <- attr(Y, "datapattern")
  sf <- attr(Y, "simfuns")
  sp <- attr(Y, "simpatterns")
  wt <- attr(Y, "weights")
  einfo <- attr(Y, "einfo")

  csr <- aargh$internal$csr %orifnull% einfo$csr

  if(is.null(fun) && is.null(sf)) {
    # No simulated functions - must compute them from simulated patterns
    if(is.null(sp))
      stop(paste("Cannot compute envelope:",
                 "Y does not contain simulated functions",
                 "(was not generated with savefuns=TRUE)",
                 "and does not contain simulated patterns",
                 "(was not generated with savepatterns=TRUE)"))
    # set default fun=Kest
    fun <- Kest
  }
  
  if(!is.null(fun)) {
    # apply new function 
    # point patterns are required
    if(is.null(sp))
      stop(paste("Object Y does not contain simulated point patterns",
                 "(attribute", dQuote("simpatterns"), ");",
                 "cannot apply a new", sQuote("fun")))
    if(is.null(X))
      stop(paste("Cannot apply a new", sQuote("fun"),
                 "; object Y generated by an older version of spatstat"))
    ## send signal if simulations were CSR
    internal <- aargh$internal
    if(csr) {
        if(is.null(internal)) internal <- list()
        internal$csr <- TRUE
    }
    ## compute new envelope
    result <- do.call(envelope,
                      resolve.defaults(list(Y=quote(X), fun=fun, simulate=sp),
                                       aargh,
                                       list(transform=transform,
                                            global=global,
                                            VARIANCE=VARIANCE,
                                            internal=internal,
                                            Yname=Yname,
                                            nsim=einfo$nsim,
                                            nsim2=einfo$nsim2,
                                            weights=wt),
                                       .StripNull=TRUE))
  } else {
    #' compute new envelope with existing simulated functions
    if(is.null(sf)) 
      stop(paste("Y does not contain a", dQuote("simfuns"), "attribute",
                 "(it was not generated with savefuns=TRUE)"))
    
    if(!is.null(transform)) {
      # Apply transformation to Y and sf
      stopifnot(is.expression(transform))
      ##      cc <- dotexpr.to.call(transform, "Y", "eval.fv")
      cc <- inject.expr("with(Y, .)", transform)
      Y <- eval(cc)
      ##      cc <- dotexpr.to.call(transform, "sf", "eval.fv")
      cc <- inject.expr("with(sf, .)", transform)
      sf <- eval(cc)
    }

    #' catch discrepancy between domains of observed and simulated functions
    if(nrow(sf) != nrow(Y)) {
      rrsim <- sf[[fvnames(sf, ".x")]]
      rrobs <- Y[[fvnames(Y, ".x")]]
      ra <- intersect.ranges(range(rrsim), range(rrobs))
      delta <- min(mean(diff(rrsim)), mean(diff(rrobs)))/2
      oksim <- (rrsim >= ra[1] - delta) & (rrsim <= ra[2] + delta)
      okobs <- (rrobs >= ra[1] - delta) & (rrobs <= ra[2] + delta)
      if(sum(oksim) != sum(okobs))
        stop("Internal error: Unable to reconcile the domains",
             "of the observed and simulated functions", call.=FALSE)
      if(mean(abs(rrsim[oksim] - rrobs[okobs])) > delta)
        stop("Internal error: Unable to reconcile the r values",
             "of the observed and simulated functions", call.=FALSE)
      sf <- sf[oksim, ,drop=FALSE]
      Y  <- Y[okobs,  ,drop=FALSE]
    }
    
    # extract simulated function values
    df <- as.data.frame(sf)
    rname <- fvnames(sf, ".x")
    df <- df[, (names(df) != rname)]

    # interface with 'envelope.matrix'
    etype <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
    dfm <- as.matrix(df)
    dont.complain.about(dfm)
    result <- do.call(envelope.matrix,
                      resolve.defaults(list(Y=quote(dfm)),
                                       aargh,
                                       list(type=etype,
                                            csr=csr,
                                            funX=Y, 
                                            Yname=Yname,
                                            weights=wt),
                                       .StripNull=TRUE))
  }

  if(!is.null(transform)) {
    # post-process labels
    labl <- attr(result, "labl")
    dnames <- colnames(result)
    dnames <- dnames[dnames %in% fvnames(result, ".")]
    # expand "."
    ud <- as.call(lapply(c("cbind", dnames), as.name))
    dont.complain.about(ud)
    expandtransform <- eval(substitute(substitute(tr, list(.=ud)),
                                       list(tr=transform[[1L]])))
    # compute new labels 
    attr(result, "fname") <- attr(Yorig, "fname")
    mathlabl <- as.character(fvlegend(result, expandtransform))
    # match labels to columns
    evars <- all.vars(expandtransform)
    used.dotnames <- evars[evars %in% dnames]
    mathmap <- match(colnames(result), used.dotnames)
    okmath <- !is.na(mathmap)
    # update appropriate labels
    labl[okmath] <- mathlabl[mathmap[okmath]]
    attr(result, "labl") <- labl
  }
  
  # Tack on envelope info
  copyacross <- c("Yname", "csr.theo", "use.theory", "simtype", "constraints")
  attr(result, "einfo")[copyacross] <- attr(Yorig, "einfo")[copyacross]
  attr(result, "einfo")$csr <- csr
  # Save data
  
  return(result)
}

pool.envelope <- local({

  pool.envelope <- function(..., savefuns=FALSE, savepatterns=FALSE) {
    Yname <- short.deparse(sys.call())
    if(nchar(Yname) > 60) Yname <- paste(substr(Yname, 1L, 40L), "[..]")
    Elist <- unname(list(...))
    nE <-  length(Elist)
    if(nE == 0) return(NULL)
    #' ........ validate envelopes .....................
    #' All arguments must be envelopes
    notenv <- !unlist(lapply(Elist, inherits, what="envelope"))
    if(any(notenv)) {
      n <- sum(notenv)
      why <- paste(ngettext(n, "Argument", "Arguments"),
                   commasep(which(notenv)),
                   ngettext(n, "does not", "do not"),
                   "belong to the class",
                   dQuote("envelope"))
      stop(why)
    }
    ## Only one envelope?
    if(nE == 1)
      return(Elist[[1L]])
    ## envelopes must be compatible
    ok <- do.call(compatible, Elist)
    if(!ok)
      stop("Envelopes are not compatible")
    ## ... reconcile parameters in different envelopes .......
    eilist <- lapply(Elist, attr, which="einfo")
    nrank  <- resolveEinfo(eilist, "nrank", 1)
    global    <- resolveEinfo(eilist, "global",   FALSE)
    ginterval    <- resolveEinfo(eilist, "ginterval", NULL, atomic=FALSE)
    VARIANCE  <- resolveEinfo(eilist, "VARIANCE", FALSE)
    alternative      <- resolveEinfo(eilist, "alternative", FALSE)
    scale <- resolveEinfo(eilist, "scale", NULL, atomic=FALSE)
    clamp <- resolveEinfo(eilist, "clamp", FALSE)
    resolveEinfo(eilist, "simtype",  "funs",
                 "Envelopes were generated using different types of simulation")
    resolveEinfo(eilist, "constraints",  "",
            "Envelopes were generated using different types of conditioning")
    resolveEinfo(eilist, "csr.theo", FALSE, NULL)
    csr         <- resolveEinfo(eilist, "csr", FALSE, NULL)
    use.weights <- resolveEinfo(eilist, "use.weights" , FALSE,
     "Weights were used in some, but not all, envelopes: they will be ignored")
    use.theory <- resolveEinfo(eilist, "use.theory", csr, NULL)
    ##
    weights <-
      if(use.weights) unlist(lapply(Elist, attr, which="weights")) else NULL
    type <- if(global) "global" else if(VARIANCE) "variance" else "pointwise"
    
    ## ........ validate saved functions .....................
    if(savefuns || !VARIANCE) {
      ## Individual simulated functions are required
      SFlist <- lapply(Elist, attr, which="simfuns")
      isnul <- unlist(lapply(SFlist, is.null))
      if(any(isnul)) {
        n <- sum(isnul)
        comply <- if(!VARIANCE) "compute the envelope:" else
                  "save the simulated functions:"
        why <- paste("Cannot", comply,
                     ngettext(n, "argument", "arguments"),
                     commasep(which(isnul)),
                     ngettext(n, "does not", "do not"),
                     "contain a", dQuote("simfuns"), "attribute",
                     "(not generated with savefuns=TRUE)")
        stop(why)
      }
      ## Simulated functions must be the same function
      fnames <- unique(lapply(SFlist, attr, which="fname"))
      if(length(fnames) > 1L) {
        fnames <- unlist(lapply(fnames, flatfname))
        stop(paste("Envelope objects contain values",
                   "of different functions:",
                   commasep(sQuote(fnames))))
      }
      ## vectors of r values must be identical
      rlist <- lapply(SFlist, getrvals)
      rvals <- rlist[[1L]]
      samer <- unlist(lapply(rlist, identical, y=rvals))
      if(!all(samer))
        stop(paste("Simulated function values are not compatible",
                   "(different values of function argument)"))
      ## Extract function values and assemble into one matrix
      matlist <- lapply(SFlist, getdotvals)
      SFmatrix <- do.call(cbind, matlist)
    }
    ## compute pooled envelope
    switch(type,
           pointwise = {
             result <- envelope(SFmatrix, funX=Elist[[1L]],
                                type=type, alternative=alternative,
                                clamp=clamp,
                                nrank=nrank,
                                csr=csr, use.theory=use.theory,
                                Yname=Yname, weights=weights,
                                savefuns=savefuns)
           },
           global = {
             simfunmatrix <- if(is.null(ginterval)) SFmatrix else {
               ## savefuns have not yet been clipped to ginterval
               ## while envelope data have been clipped.
               domain <- (rvals >= ginterval[1L]) & (rvals <= ginterval[2L])
               SFmatrix[domain, , drop=FALSE]
             }
             result <- envelope(simfunmatrix, funX=Elist[[1L]],
                                type=type, alternative=alternative,
                                scale=scale, clamp=clamp,
                                csr=csr, use.theory=use.theory,
                                nrank=nrank,
                                ginterval=ginterval,
                                Yname=Yname, weights=weights,
                                savefuns=savefuns)
           },
           variance = {
             ## Pool sample means and variances
             nsims <- unlist(lapply(eilist, getElement, name="nsim"))
             mmeans <- lapply(Elist, getElement, name="mmean")
             vars   <- lapply(Elist, getElement, name="var")
             mmeans <- matrix(unlist(mmeans), ncol=nE)
             vars   <- matrix(unlist(vars),   ncol=nE)
             if(!use.weights) {
               w.mean <- nsims
               d.mean <- sum(nsims)
               w.var  <- nsims - 1
               d.var  <- sum(nsims) - 1
             } else {
               weightlist <- lapply(Elist, attr, which="weights")
               w.mean <- unlist(lapply(weightlist, sum))
               d.mean <- sum(w.mean)
               ssw <- unlist(lapply(weightlist, meansqfrac))
               ##  meansqfrac :  function(x) {sum((x/sum(x))^2)}))
               w.var  <- w.mean * (1 - ssw)
               d.var <-  d.mean * (1 - sum(ssw))
             }
             poolmmean <- as.numeric(mmeans %*% matrix(w.mean/d.mean, ncol=1L))
             within <- vars %*% matrix(w.var, ncol=1L)
             between <- ((mmeans - poolmmean[])^2) %*% matrix(w.mean, ncol=1L)
             poolvar <- as.numeric((within + between)/d.var)
             ## feed precomputed data to envelope.matrix
             pc <- list(Ef=poolmmean[],
                        varf=poolvar[])
             nsim <- sum(nsims)
             result <- envelope.matrix(NULL, funX=Elist[[1L]],
                                       type=type, alternative=alternative,
                                       csr=csr, Yname=Yname,
                                       weights=weights,
                                       savefuns=savefuns,
                                       nsim=nsim,
                                       precomputed=pc)
           })
  
    ## Copy envelope info that is not handled by envelope.matrix
    copyacross <- c("Yname", "csr.theo", "use.theory", "simtype", "constraints")
    attr(result, "einfo")[copyacross] <- attr(Elist[[1L]], "einfo")[copyacross]
  
    ## ..............saved patterns .....................
    if(savepatterns) {
      SPlist <- lapply(Elist, attr, which="simpatterns")
      isnul <- unlist(lapply(SPlist, is.null))
      if(any(isnul)) {
        n <- sum(isnul)
        why <- paste("Cannot save the simulated patterns:",
                     ngettext(n, "argument", "arguments"),
                     commasep(which(isnul)),
                     ngettext(n, "does not", "do not"),
                     "contain a", dQuote("simpatterns"), "attribute",
                     "(not generated with savepatterns=TRUE)")
        warning(why)
      } else {
        attr(result, "simpatterns") <- Reduce(append, SPlist)
      }
    }
    ## ..............saved summary functions ................
    if(savefuns) {
      alldata <- cbind(rvals, SFmatrix)
      Nsim <- ncol(SFmatrix)
      simnames <- paste0("sim", 1:Nsim)
      colnames(alldata) <- c("r", simnames)
      alldata <- as.data.frame(alldata)
      SFtemplate <- SFlist[[1L]]
      SimFuns <- fv(alldata,
                    argu="r",
                    ylab=attr(SFtemplate, "ylab"),
                    valu="sim1",
                    fmla= deparse(. ~ r),
                    alim=attr(SFtemplate, "alim"),
                    labl=names(alldata),
                    desc=c("distance argument r",
                      paste("Simulation ", 1:Nsim, sep="")),
                    fname=attr(SFtemplate, "fname"),
                    yexp=attr(SFtemplate, "yexp"),
                    unitname=unitname(SFtemplate))
      fvnames(SimFuns, ".") <- simnames
      attr(result, "simfuns") <- SimFuns
    } 
  
    dotnames   <- lapply(Elist, fvnames, a=".")
    dn <- dotnames[[1L]]
    if(all(unlist(lapply(dotnames, identical, y=dn))))
      fvnames(result, ".") <- dn
    
    shadenames <- lapply(Elist, fvnames, a=".s")
    sh <- shadenames[[1L]]
    if(all(unlist(lapply(shadenames, identical, y=sh))))
      fvnames(result, ".s") <- sh
  
    return(result)
  }

  getrvals <- function(z) { as.matrix(z)[, fvnames(z, ".x")] }
  
  getdotvals <- function(z) { as.matrix(z)[, fvnames(z, "."), drop=FALSE] }

  meansqfrac <- function(x) {sum((x/sum(x))^2)} 

  pool.envelope
})

# resolve matching entries in different envelope objects
#   x is a list of envelope info objects

resolveEinfo <- function(x, what, fallback, warn, atomic=TRUE) {
  if(atomic) {
    y <- unique(unlist(lapply(x, getElement, name=what)))
    if(length(y) == 1L)
      return(y)
  } else {
    y <- unique(lapply(x, getElement, name=what))
    if(length(y) == 1L)
      return(y[[1L]])
  }
  if(missing(warn))
    warn <- paste("Envelopes were generated using different values",
                  "of argument", paste(sQuote(what), ";", sep=""),
                  "reverting to default value")
  if(!is.null(warn))
    warning(warn, call.=FALSE)
  return(fallback)
}

as.data.frame.envelope <- function(x, ..., simfuns=FALSE) {
  if(simfuns && !is.null(sf <- attr(x, "simfuns"))) {
    # tack on the simulated functions as well
    y <- as.data.frame(bind.fv(x, sf, clip=TRUE))
    return(y)
  } 
  NextMethod("as.data.frame")
}