#
#    summary.ppm.R
#
#    summary() method for class "ppm"
#
#    $Revision: 1.37 $   $Date: 2010/08/08 02:05:05 $
#
#    summary.ppm()
#    print.summary.ppm()
#
summary.ppm <- function(object, ..., quick=FALSE) {
  verifyclass(object, "ppm")

  x <- object
  y <- list()
  
  #######  Extract main data components #########################

  QUAD <- object$Q
  DATA <- QUAD$data
  TREND <- x$trend

  INTERACT <- x$interaction
  if(is.null(INTERACT)) INTERACT <- Poisson()
  
  #######  Check version #########################

  mpl.ver <- versionstring.ppm(object)
  int.ver <- versionstring.interact(INTERACT)
  current <- versionstring.spatstat()

  virgin <- min(package_version(c(mpl.ver, int.ver)))

  y$antiquated <- antiquated <- (virgin <= package_version("1.5"))
  y$old        <- old        <- (virgin < majorminorversion(current))

  y$version    <- as.character(virgin)
                
  ####### Determine type of model ############################
  
  y$entries <- list()
  y$no.trend <- identical.formulae(TREND, NULL) || identical.formulae(TREND, ~1)
  y$trendvar <- trendvar <- variablesinformula(TREND)
  y$stationary <- y$no.trend || all(trendvar == "marks")

  y$poisson <- is.poisson.interact(INTERACT)

  y$marked <- is.marked.ppp(DATA)
  y$multitype <- is.multitype.ppp(DATA)
  y$marktype <- if(y$multitype) "multitype" else
                if(y$marked) "marked" else "unmarked"

  if(y$marked) y$entries$marks <- marks(DATA)

  y$name <- paste(
          if(y$stationary) "Stationary " else "Nonstationary ",
          if(y$poisson) {
            if(y$multitype) "multitype "
            else if(y$marked) "marked "
            else ""
          },
          INTERACT$name,
          sep="")

  y$method <- x$method

  y$problems <- x$problems

  ######  Fitting algorithm ########################################

  y$fitter <- x$fitter
  if(y$fitter %in% c("glm", "gam"))
    y$converged <- x$internal$glmfit$converged
        
  ######  Extract fitted model coefficients #########################

  y$entries$coef <- COEFS <- x$coef

  ###### Extract fitted interaction and summarise  #################
  
  FITIN <- fitin(x)
  y$interaction <- summary(FITIN)

  y$entries$Vnames <- Vnames <- x$internal$Vnames

  # Exit here if quick=TRUE
    
  if(is.logical(quick) && quick) {
    class(y) <- "summary.ppm"
    return(y)
  }
  

  ######  Does it have external covariates?  ####################

  if(!antiquated) {
    covars <- x$covariates
    y$has.covars <- !is.null(covars) && (length(covars) > 0)
    used <- (y$trendvar %in% names(covars))
    y$covars.used <- y$trendvar[used]
    y$uses.covars <- sum(used)
  } else {
    # Antiquated format
    # Interpret the function call instead
    callexpr <- parse(text=x$call)
    callargs <- names(as.list(callexpr[[1]]))
    # Data frame of covariates was called 'data' in versions up to 1.4-x
    y$has.covars <- !is.null(callargs) && !is.na(pmatch("data", callargs))
    # conservative guess
    y$uses.covars <- y$has.covars
  }
    
  ######  Arguments in call ####################################
  
  y$args <- x[c("call", "correction", "rbord")]
  
  #######  Main data components #########################

  y$entries <- append(list(quad=QUAD,
                           data=DATA,
                           interaction=INTERACT),
                      y$entries)

  ####### Summarise data ############################

  y$data <- summary.ppp(DATA, checkdup=FALSE)
  y$quad <- summary.quad(QUAD, checkdup=FALSE)

  if(is.character(quick) && (quick == "no prediction"))
    return(y)
  
  ######  Trend component #########################

  y$trend <- list()

  y$trend$name <- if(y$poisson) "Intensity" else "Trend"

  y$trend$formula <- if(y$no.trend) NULL else TREND

  if(y$poisson && y$no.trend) {
    # uniform Poisson process
    y$trend$value <- lambda <- exp(COEFS[[1]])
    y$trend$label <- switch(y$marktype,
                            unmarked="Uniform intensity",
                            multitype="Uniform intensity for each mark level",
                            marked="Uniform intensity in product space",
                            "")
  } else if(y$stationary) {
    # stationary
    switch(y$marktype,
           unmarked={
             # stationary non-poisson non-marked
             y$trend$label <- "First order term"
             y$trend$value <- c(beta=exp(COEFS[[1]]))
           },
           multitype={
             # stationary, multitype
             mrk <- marks(DATA)
             y$trend$label <-
               if(y$poisson) "Intensities" else "First order terms"
             # Use predict.ppm to evaluate the fitted intensities
             lev <- factor(levels(mrk), levels=levels(mrk))
             nlev <- length(lev)
             marx <- list(x=rep(0, nlev), y=rep(0, nlev), marks=lev)
             betas <- predict(x, locations=marx, type="trend")
             names(betas) <- paste("beta_", as.character(lev), sep="")
             y$trend$value <- betas
           },
           marked={
             # stationary, marked
             y$trend$label <- "Fitted intensity coefficients"
             y$trend$value <- blankcoefnames(COEFS)
           })
  } else {
    # not stationary 
    y$trend$label <- "Fitted coefficients for trend formula"
    # extract trend terms without trying to understand them much
    if(is.null(Vnames)) 
      trendbits <- COEFS
    else {
      agree <- outer(names(COEFS), Vnames, "==")
      whichbits <- apply(!agree, 1, all)
      trendbits <- COEFS[whichbits]
    }
    y$trend$value <- blankcoefnames(trendbits)
  }
  
  class(y) <- "summary.ppm"
  return(y)
}

print.summary.ppm <- function(x, ...) {

  if(x$old)
    warning("Model was fitted by an older version of spatstat")
  
  if(is.null(x$args)) {
    # this is the quick version
    cat(paste(x$name, "\n"))
    return(invisible(NULL))
  }

  # otherwise - full details
  cat("Point process model\n")
  methodchosen <-
    if(is.null(x$method))
      "unspecified method"
    else 
      switch(x$method,
             mpl="maximum pseudolikelihood (Berman-Turner approximation)",
             ho="Huang-Ogata method (approximate maximum likelihood)",
             paste("unrecognised method", sQuote(x$method)))
  cat(paste("Fitting method:", methodchosen, "\n"))
  howfitted <- switch(x$fitter,
                      exact= "analytically",
                      gam  = "using gam()",
                      glm  = "using glm()",
                      ho   = NULL,
                      paste("using unrecognised fitter", sQuote(x$fitter)))
  if(!is.null(howfitted)) cat(paste("Model was fitted", howfitted, "\n"))
  if(x$fitter %in% c("glm", "gam")) {
    if(x$converged) cat("Algorithm converged\n")
    else cat("*** Algorithm did not converge ***\n")
  }

  cat("Call:\n")
  print(x$args$call)

  if(x$old) 
    cat(paste("** Executed by old spatstat version", x$version, " **\n"))
  
  cat(paste("Edge correction:", dQuote(x$args$correction)))
  if(x$args$rbord > 0)
    cat(paste("border correction distance r =", x$args$rbord,"\n"))

  cat("\n----------------------------------------------------\n")

  # print summary of quadrature scheme
  print(x$quad)
  
  cat("\n----------------------------------------------------\n")
  cat("FITTED MODEL:\n\n")

  # This bit is currently identical to print.ppm()
  # except for a bit more fanfare
  # and the inclusion of the 'gory details' bit
  
  notrend <-    x$no.trend
  stationary <- x$stationary
  poisson <-    x$poisson
  markeddata <- x$marked
  multitype  <- x$multitype
        
  markedpoisson <- poisson && markeddata

  # ----------- Print model type -------------------
        
  cat(x$name)
  cat("\n")

  if(markeddata) mrk <- x$entries$marks
  if(multitype) {
    cat("Possible marks: \n")
    cat(paste(levels(mrk)))
  }

  # ----- trend --------------------------

  cat(paste("\n\n ---- ", x$trend$name, ": ----\n\n", sep=""))

  if(!notrend) {
    cat("Trend formula: ")
    print(x$trend$formula)
    if(x$uses.covars) 
      cat(paste("Model depends on external",
                ngettext(length(x$covars.used), "covariate", "covariates"),
                commasep(sQuote(x$covars.used)), "\n"))
    else if(x$has.covars)
      cat("Model object contains external covariates\n")
  }
        
  cat(paste("\n", x$trend$label, ":\n", sep=""))
  
  tv <- x$trend$value
  if(!is.list(tv))
    print(tv)
  else 
    for(i in seq(tv))
      print(tv[[i]])
        
  # ---- Interaction ----------------------------

  if(!poisson) {
    cat("\n\n ---- Interaction: -----\n\n")
    print(x$interaction)
  }

  ####### Gory details ###################################
  cat("\n\n----------- gory details -----\n")
  COEFS <- x$entries$coef
  
  cat("\nFitted regular parameters (theta): \n")
  print(COEFS)

  cat("\nFitted exp(theta): \n")
  print(exp(unlist(COEFS)))

  ##### Warnings issued #######

  probs <- x$problems
  if(!is.null(probs) && is.list(probs) && (length(probs) > 0)) 
    lapply(probs,
           function(a) {
             if(is.list(a) && !is.null(p <- a$print))
               cat(paste("Problem:\n", p, "\n\n"))
           })
          
  return(invisible(NULL))
}

no.trend.ppm <- function(x) {
  summary.ppm(x, quick=TRUE)$no.trend
}

is.stationary.ppm <- function(x) {
  summary.ppm(x, quick=TRUE)$stationary
}

is.poisson.ppm <- function(x) {
  summary.ppm(x, quick=TRUE)$poisson
}

is.marked.ppm <- function(X, ...) {
  summary.ppm(X, quick=TRUE)$marked
}

is.multitype.ppm <- function(X, ...) {
  summary.ppm(X, quick=TRUE)$multitype
}

blankcoefnames <- function(x) {
  # remove name labels from ppm coefficients
  # First decide whether there are 'labels within labels'
  unlabelled <- unlist(lapply(x,
                              function(z) { is.null(names(z)) } ))
  if(all(unlabelled))
    value <- unlist(x)
  else {
    value <- list()
    for(i in seq(x))
      value[[i]] <- if(unlabelled[i]) unlist(x[i]) else x[[i]]
  }
  return(value)
}