#
#  effectfun.R
#
#   $Revision: 1.9 $ $Date: 2012/01/24 10:29:03 $
#

effectfun <- function(model, covname, ..., se.fit=FALSE) {
  stopifnot(is.ppm(model))
  dotargs <- list(...)
  # determine names of covariates involved
  intern.names <-
    if(is.marked.ppm(model)) c("x", "y", "marks") else c("x", "y")
  co <- model$covariates
  extern.names <- names(co)
  # find the relevant covariate 
  if(missing(covname))
    stop("covname must be provided")
  if(!(covname %in% c(intern.names, extern.names)))
    stop(paste("model does not have a covariate called", sQuote(covname)))
  # check that fixed values for all other covariates are provided 
  given.covs <- names(dotargs)
  if(any(uhoh <- !(extern.names %in% c(given.covs, covname)))) {
    nuh <- sum(uhoh)
    stop(paste(ngettext(nuh,
                        "A value for the covariate",
                        "Values for the covariates"),
               commasep(dQuote(extern.names[uhoh])),
               "must be provided (as",
               ngettext(nuh, "an argument", "arguments"),
               "to effectfun)"))
  }
  # establish type and range of covariate values
  N0 <- 256
  if(covname == "x") {
    covtype <- "real"
    W <- as.owin(data.ppm(model))
    Zr <- W$xrange
    Zvals <- seq(from=Zr[1], to=Zr[2], length.out=N0)
  } else if(covname == "y") {
    covtype <- "real"
    W <- as.owin(data.ppm(model))
    Zr <- W$yrange
    Zvals <- seq(from=Zr[1], to=Zr[2], length.out=N0)
  } else if(covname == "marks") {
    covtype <- "factor"
    Zvals <- levels(marks(data.ppm(model)))
  } else {
    # covariate is external
    if(is.data.frame(co)) {
      Z <- co$covname
      covtype <- typeof(Z)
      if(covtype == "double")
        covtype <- "real"
      switch(covtype,
             real={
               Zr <- range(Z)
               Zvals <- seq(from=Zr[1], to=Zr[2], length.out=N0)
             },
             integer={
               Zr <- range(Z)
               Zvals <- seq(from=Zr[1], to=Zr[2], by=ceiling((diff(Zr)+1)/N0))
             },
             factor={
               Zvals <- levels(Z)
             },
             logical={
               Zvals <- c(FALSE, TRUE)
             },
             stop(paste("Cannot handle covariate of type", dQuote(covtype)))
             )
    } else if(is.list(co)) {
      Z <- co[[covname]]
      # convert to image
      if(!is.im(Z))
        Z <- as.im(Z, W=as.owin(model))
      covtype <- Z$type
      switch(covtype,
             real={
               Zr <- summary(Z)$range
               Zvals <- seq(from=Zr[1], to=Zr[2], length.out=N0)
             },
             factor={
               Zvals <- levels(Z)
             },
             logical={
               Zvals <- c(FALSE, TRUE)
             },
             stop(paste("Cannot handle covariate of type", dQuote(covtype)))
             )
    } else stop("Unrecognised format for covariates in model")
  }
  # set up data frames of fake data for predict method
  # First set up default, constant value for each covariate
  N <- length(Zvals)
  fakeloc <- resolve.defaults(dotargs,
                              list(x=0, y=0))[c("x","y")]
  if(is.marked.ppm(model)) {
    lev <- levels(marks(data.ppm(model)))
    fakeloc$marks <- lev[1]
  }
  fakeloc <- lapply(fakeloc, function(x,N) { rep(x[1],N)}, N=N)
  fakecov <- lapply(dotargs, function(x,N) { rep(x[1],N)}, N=N)
  # Overwrite value for covariate of interest
  if(covname %in% intern.names)
    fakeloc[[covname]] <- Zvals
  else fakecov[[covname]] <- Zvals
  # convert to data frame
  fakeloc <- do.call("data.frame", fakeloc)
  fakecov <- if(length(fakecov) > 0) do.call("data.frame", fakecov) else NULL
  #
  # Now predict
  lambda <- predict(model, locations=fakeloc, covariates=fakecov)
  if(se.fit) {
    se <- predict(model, locations=fakeloc, covariates=fakecov, type="se")
    sedf <- data.frame(se =se,
                       hi = lambda + 2 * se,
                       lo = lambda - 2 * se)
  }
  #
  dfin <- if(!is.null(fakecov)) cbind(fakeloc, fakecov) else fakeloc 
  dfin <- dfin[covname]
  df <- cbind(dfin, data.frame(lambda=lambda))
  #
  if(covtype != "real") {
    result <- df
    if(se.fit) result <- cbind(result, sedf)
  } else {
    bc <- paren(covname)
    result <- fv(df, argu=covname, 
                 ylab=substitute(lambda(X), list(X=as.name(covname))),
                 labl=c(covname,
                   paste("hat(%s)", bc)),
                 valu="lambda", alim=Zr,
                 desc=c(paste("value of covariate", covname),
                   "fitted intensity"),
                 fname="lambda")
    if(se.fit) {
      result <- bind.fv(result, sedf,
                        labl=c(paste("se[%s]", bc),
                          paste("%s[hi]", bc),
                          paste("%s[lo]", bc)),
                        desc=c("standard error of fitted trend",
                          "upper limit of pointwise 95%% CI for trend",
                          "lower limit of pointwise 95%% CI for trend"))
      fvnames(result, ".") <- c("lambda", "hi", "lo")
      attr(result, "fmla") <- as.formula(paste(". ~ ", covname))
    }
  }
  return(result)
}