.packageName <- "randomForest"
MDSplot <- function(rf, fac, k=2, ...) {
  if(!inherits(rf, "randomForest")) {
    stop(deparse(substitute(rf)), " must be a randomForest object")
  }
  if(is.null(rf$proximity)) {
    stop(deparse(substitute(rf)), " does not contain a proximity matrix")
  }
  op <- par(pty="s")
  on.exit(par(op))
  Rver <- substring(paste(R.version[c("major", "minor")], collapse="."), 1, 3)
  cmdscale <- if (as.numeric(Rver) >= 1.9) stats:::cmdscale else mva:::cmdscale
  rf.mds <- cmdscale(1 - rf$proximity, eig=TRUE, k=k)
  colnames(rf.mds$points) <- paste("Dim", 1:k)
  nlevs <- length(levels(fac))
  if( require(RColorBrewer) && nlevs < 12)
    pal <- brewer.pal(nlevs,"Set1")
  else
    pal <- rainbow(nlevs)
  if(k <= 2) {
    plot(rf.mds$points, col=pal[as.numeric(fac)], pch=20, ...)
  } else {
    pairs(rf.mds$points, col=pal[as.numeric(fac)], pch=20, ...)
  }
  invisible(rf.mds)
}
combine <- function(...) {
   pad0 <- function(x, len) c(x, rep(0, len-length(x)))
   padm0 <- function(x, len) rbind(x, matrix(0, nrow=len-nrow(x),
                                             ncol=ncol(x)))
   rflist <- list(...)
   areForest <- sapply(rflist, function(x) inherits(x, "randomForest")) 
   if (any(!areForest)) stop("Argument must be a list of randomForest objects")
   ## Use the first component as a template
   rf <- rflist[[1]]
   classRF <- rf$type == "classification"
   trees <- sapply(rflist, function(x) x$ntree)
   ntree <- sum(trees)
   rf$ntree <- ntree
   nforest <- length(rflist)
   haveTest <- !any(sapply(rflist, function(x) is.null(x$test)))
   
   ## Combine the forest component, if any
   haveForest <- sapply(rflist, function(x) !is.null(x$forest))
   if (all(haveForest)) {
       nrnodes <- max(sapply(rflist, function(x) x$forest$nrnodes))
       rf$forest$nrnodes <- nrnodes
       rf$forest$ndbigtree <-
           unlist(sapply(rflist, function(x) x$forest$ndbigtree))
       rf$forest$nodestatus <-
           do.call("cbind", lapply(rflist, function(x)
                                   padm0(x$forest$nodestatus, nrnodes)))
       rf$forest $bestvar <-
           do.call("cbind",
                   lapply(rflist, function(x)
                          padm0(x$forest$bestvar, nrnodes)))
       rf$forest$xbestsplit <-
           do.call("cbind",
                   lapply(rflist, function(x)
                          padm0(x$forest$xbestsplit, nrnodes)))
       rf$forest$nodepred <-
           do.call("cbind", lapply(rflist, function(x)
                                   padm0(x$forest$nodepred, nrnodes)))
       tree.dim <- dim(rf$forest$treemap)
       rf$forest$treemap <-
           array(unlist(lapply(rflist, function(x) apply(x$forest$treemap, 2:3,
                                                         pad0, nrnodes))),
                 c(nrnodes, 2, ntree))
       rf$forest$ntree <- ntree
       if (classRF) rf$forest$cutoff <- rflist[[1]]$forest$cutoff
   } else {
       rf$forest <- NULL
   }
   
   if (classRF) {
       ## Combine the votes matrix: 
       rf$votes <- 0
       rf$oob.times <- 0
       areVotes <- all(sapply(rflist, function(x) any(x$votes > 1)))
       if (areVotes) {
           for(i in 1:nforest) {
               rf$oob.times <- rf$oob.times + rflist[[i]]$oob.times
               rf$votes <- rf$votes + rflist[[i]]$votes
           }
       } else {
           for(i in 1:nforest) {
               rf$oob.times <- rf$oob.times + rflist[[i]]$oob.times            
               rf$votes <- rf$votes + rflist[[i]]$votes * rflist[[i]]$oob.times
           }
           rf$votes <- rf$votes / rf$oob.times
       }
       rf$predicted <- factor(colnames(rf$votes)[max.col(rf$votes)],
                              levels=levels(rf$predicted))
       if(haveTest) {
           rf$test$votes <- 0
           if (any(rf$test$votes > 1)) {
               for(i in 1:nforest)
                   rf$test$votes <- rf$test$votes + rflist[[i]]$test$votes
           } else {
               for (i in 1:nforest)
                   rf$test$votes <- rf$test$votes +
                       rflist[[i]]$test$votes * rflist[[i]]$ntree
           }
           rf$test$predicted <-
               factor(colnames(rf$test$votes)[max.col(rf$test$votes)],
                      levels=levels(rf$test$predicted))
       }
   } else {
       rf$predicted <- 0
       for (i in 1:nforest) rf$predicted <- rf$predicted +
           rflist[[i]]$predicted * rflist[[i]]$ntree
       rf$predicted <- rf$predicted / ntree
       if (haveTest) {
           rf$test$predicted <- 0
           for (i in 1:nforest) rf$test$predicted <- rf$test$predicted +
               rflist[[i]]$test$predicted * rflist[[i]]$ntree
           rf$test$predicted <- rf$test$predicted / ntree
       }
   }
   
   ## If variable importance is in all of them, compute the average
   ## (weighted by the number of trees in each forest)
   have.imp <- !any(sapply(rflist, function(x) is.null(x$importance)))
   if(have.imp) {
       rf$importance <- rf$importanceSD <- 0
       for(i in 1:nforest) {
           rf$importance <- rf$importance +
               rflist[[i]]$importance * rflist[[i]]$ntree
           rf$importance <- rf$importance / ntree
           ## Do the same thing with SD of importance, though that's not
           ## exactly right...
           rf$importanceSD <- rf$importanceSD +
               rflist[[i]]$importanceSD^2 * rflist[[i]]$ntree
           rf$importance <- sqrt(rf$importance / ntree)
       }
       haveCaseImp <- !any(sapply(rflist, function(x)
                                  is.null(x$localImportance)))
       ## Average casewise importance
       if (haveCaseImp) {
           rf$localImportance <- 0
           for (i in 1:nforest) {
               rf$localImportance <- rf$localImportance +
                   rflist[[i]]$localImportance
           }
           rf$localImportance <- rf$localImportance / ntree
       }
   }
   
   ## If proximity is in all of them, compute the average
   ## (weighted by the number of trees in each forest)
   have.prox <- !any(sapply(rflist, function(x) is.null(x$proximity)))
   if(have.prox) {
       rf$proximity <- 0
       for(i in 1:nforest)
           rf$proximity <- rf$proximity + rflist[[i]]$proximity * rflist[[i]]$ntree
       rf$proximity <- rf$proximity / ntree
   }
   
   ## Set confusion matrix and error rates to NULL
   if(classRF) {
       rf$confusion <- NULL
       rf$err.rate <- NULL
       if(haveTest) {
           rf$test$confusion <- NULL
           rf$err.rate <- NULL
       }
   } else {
       rf$mse <- rf$rsq <- NULL
       if(haveTest) rf$test$mse <- rf$test$rsq <- NULL
   }   
   rf
}
getTree <- function(rfobj, k=1) {
  if (is.null(rfobj$forest)) {
    stop("No forest component in ", deparse(substitute(rfobj)))
  }
  if (k > rfobj$ntree) {
    stop("There are fewer than ", k, "trees in the forest")
  }
  tree <- cbind(rfobj$forest$treemap[,,k], rfobj$forest$bestvar[,k],
                rfobj$forest$xbestsplit[,k], rfobj$forest$nodestatus[,k],
                rfobj$forest$nodepred[,k])[1:rfobj$forest$ndbigtree[k],]
  colnames(tree) <- c("left daughter", "right daughter", "split var",
                      "split point", "status", "prediction")
  tree
}

grow <- function(x, ...) UseMethod("grow")

grow.default <- function(x, ...)
  stop("grow has not been implemented for this class of object")

grow.randomForest <- function(x, how.many, ...) {
  y <- update(x, ntree=how.many)
  combine(x, y)
}
importance <- function(x, ...)  UseMethod("importance")

importance.default <- function(x, ...)
    stop("No method implemented for this class of object")

importance.randomForest <- function(x, type=NULL, class=NULL, scale=TRUE,
                                    ...) {
    if (!inherits(x, "randomForest"))
        stop("x is not of class randomForest")
    classRF <- x$type != "regression"
    hasImp <- !is.null(dim(x$importance))
    hasType <- !is.null(type)
    allImp <- is.null(type) && hasImp
    if (hasType) {
        if (!(type %in% 1:2)) stop("Wrong type specified")
        if (type == 1 && !hasImp) stop("That measure was not computed")
        if (type == 2 && !is.null(class))
            stop("No class-specific measure for that type")
    }
    
    imp <- x$importance
    if (hasType && type == 2) {
        if (hasImp) imp <- imp[, ncol(imp)]
    } else {
        if (hasImp) {
            if (scale) {
                SD <- x$importanceSD
                imp[,-ncol(imp)] <-
                    imp[,-ncol(imp)] / ifelse(SD < .Machine$double.eps, 1, SD)
            }
            if (!allImp) {
                if (is.null(class)) {
                    ## The average decrease in accuracy measure:
                    imp <- imp[, ncol(imp)-1]
                } else {
                    if (classRF) {
                        whichCol <- match(class, colnames(x$importance))
                    } else {
                        whichCol <- 1
                    }
                    imp <- imp[, whichCol]
                }
                names(imp) <- rownames(x$importance)
            }
        }
    }
    imp
}
margin <- function(rf, observed) {
    if( !inherits(rf, "randomForest") ) {
        stop("margin defined for Random Forests")
    }
    if( is.null(rf$votes) ) {
        stop("margin is only defined if votes are present")
    }
    if( !is.factor(observed) ) {
        stop(deparse(substitute(observed)), " is not a factor")
    }
    augD <- rf$votes
    if( any(augD > 1) ) {
        augD <- sweep(augD, 1, rowSums(augD), "/")
    }
    augD <- data.frame(augD, observed)
    names(augD) <- c(dimnames(rf$votes)[[2]], "observed")
    nlev <- length(levels(observed))
    
    ans<- apply(augD, 1, function(x) { pos <- match(x[nlev+1], names(x));
                                       t1 <- as.numeric(x[pos]);
                                       t2 <- max(as.numeric(x[-c(pos, nlev+1)]));
                                       t1 - t2 }
                )
    names(ans) <- observed
    class(ans) <- "margin"
    ans
}

plot.margin <- function(x, sort=TRUE, ...) {
    if (sort) x <- sort(x)
    nF <- factor(names(x))
    nlevs <- length(levels(nF))
    if ( require(RColorBrewer) && nlevs < 12) {
        pal <- brewer.pal(nlevs,"Set1")
    } else {
        pal <- rainbow(nlevs)
    }
    plot.default(x, col=pal[as.numeric(nF)], pch=20, ... )
}
na.roughfix <- function(object, ...)
  UseMethod("na.roughfix")

na.roughfix.data.frame <- function(object, ...) {
##  n <- length(object)
##  vars <- seq(length = n)
  isfac <- sapply(object, is.factor)
  isnum <- sapply(object, is.numeric)
  hasNA <- which(sapply(object, function(x) any(is.na(x))))
  if (any(!(isfac | isnum)))
      stop("na.roughfix only works for numeric or factor")
  for (j in hasNA) {
    if (isfac[j]) {
      freq <- table(object[[j]])
      xmode <- names(freq)[max.col(rbind(freq))]
      object[[j]][is.na(object[[j]])] <- xmode
    } else {
      xmed <- median(object[[j]], na.rm=TRUE)
      object[[j]][is.na(object[[j]])] <- xmed
    }
  }
  object
}

na.roughfix.default <- function(object, ...) {
  if (!is.atomic(object)) 
    return(object)
  d <- dim(object)
  if (length(d) > 2) 
    stop("can't handle objects with more than two dimensions")
  if (all(!is.na(object)))
    return(object)
  if (!is.numeric(object))
    stop("roughfix can only deal with numeric data.")
  if (d == 2) {
      hasNA <- which(apply(object, 2, function(x) any(is.na(x))))
      for (j in hasNA) 
          object[is.na(object[, j]), j] <- median(object[, j], na.rm=TRUE)
  } else {
      object[is.na(object)] <- median(object, na.rm=TRUE)
  }
  object
}
partialPlot <- function(x, ...) UseMethod("partialPlot")

partialPlot.default <- function(x, ...)
    stop("partial dependence plot not implemented for this class of objects.\n")

partialPlot.randomForest <-
    function (x, pred.data, x.var, which.class, add = FALSE,
              n.pt = min(length(unique(pred.data[, xname])), 51), rug = TRUE,
              xlab=deparse(substitute(x.var)), ylab="",
              main=paste("Partial Dependence on", deparse(substitute(x.var))),
              ...) 
{
    classRF <- x$type != "regression"
    if (is.null(x$forest)) 
        stop("The randomForest object must contain the forest.\n")
    xname <- if (is.name(substitute(x.var))) deparse(substitute(x.var)) else eval(substitute(x.var))
    xv <- pred.data[, xname]
    n <- nrow(pred.data)
    if (classRF) {
        if (missing(which.class)) {
            focus <- 1
        }
        else {
            focus <- charmatch(which.class, colnames(x$votes))
            if (is.na(focus)) 
                stop(which.class, "is not one of the class labels.")
        }
    }
    if (is.factor(xv) && !is.ordered(xv)) {
        x.pt <- levels(xv)
        y.pt <- numeric(length(x.pt))
        for (i in seq(along = x.pt)) {
            x.data <- pred.data
            x.data[, xname] <- factor(rep(x.pt[i], n), levels = x.pt)
            if (classRF) {
                pr <- predict(x, x.data, type = "prob")
                y.pt[i] <- mean(log(ifelse(pr[, focus] > 0,
                                           pr[, focus], 1)) -
                                rowMeans(log(ifelse(pr > 0, pr, 1))))
            } else {
                y.pt[i] <- mean(predict(x, x.data))
            }
        }
        if (add) {
            points(1:length(x.pt), y.pt, type = "h", lwd = 2, 
                   ...)
        } else {
            barplot(1:length(x.pt), y.pt, col="blue", xlab = xlab, 
                    ylab = ylab, main=main, ...)
        }
    } else {
        if (is.ordered(xv)) 
            xv <- as.numeric(xv)
        x.pt <- seq(min(xv), max(xv), length = n.pt)
        y.pt <- numeric(length(x.pt))
        for (i in seq(along = x.pt)) {
            x.data <- pred.data
            x.data[, xname] <- rep(x.pt[i], n)
            if (classRF) {
                pr <- predict(x, x.data, type = "prob")
                y.pt[i] <- mean(log(ifelse(pr[, focus] == 0, 1, pr[, focus]))
                                - rowMeans(log(ifelse(pr == 0, 1, pr))))
            } else {
                y.pt[i] <- mean(predict(x, x.data))
            }
        }
        if (add) {
            lines(x.pt, y.pt, ...)
        } else {
            plot(x.pt, y.pt, type = "l", 
                 xlab=xlab, ylab=ylab, main = main, ...)
        }
        if (rug) {
            if (n.pt > 10) {
                rug(quantile(xv, seq(0.1, 0.9, by = 0.1)), side = 1)
            } else {
                rug(unique(xv, side = 1))
            }
        }
    }
    invisible(list(x = x.pt, y = y.pt))
}
plot.randomForest <- function(x, type="l", main=deparse(substitute(x)), ...) {
  if(x$type == "unsupervised")
    stop("No plot for unsupervised randomForest.")
  test <- !(is.null(x$test$mse) || is.null(x$test$err.rate))
  if(x$type == "regression") {
    err <- x$mse
    if(test) err <- cbind(err, x$test$mse)
  } else {
    err <- x$err.rate
    if(test) err <- cbind(err, x$test$err.rate)
  }
  if(test) {
    colnames(err) <- c("OOB", "Test")
    matplot(1:x$ntree, err, type = type, xlab="trees", ylab="Error",
            main=main, ...)
  } else {
    matplot(1:x$ntree, err, type = type, xlab="trees", ylab="Error",
            main=main, ...)
  }
  invisible(err)
}

  
"predict.randomForest" <-
  function (object, newdata, type = "response", norm.votes = TRUE,
            predict.all=FALSE, proximity = FALSE, nodes=FALSE, ...) 
{
  if (!inherits(object, "randomForest")) 
    stop("object not of class randomForest")
  if (is.null(object$forest)) stop("No forest component in the object")
  out.type <- charmatch(tolower(type),
                        c("response", "prob", "vote", "class"))
  if (is.na(out.type)) 
    stop("type must be one of 'response', 'prob', 'vote'")
  if (out.type == 4) out.type <- 1
  if (out.type != 1 && object$type == "regression")
    error("'prob' or 'vote' not meaningful for regression")
  if (out.type == 2) 
    norm.votes <- TRUE
  if (missing(newdata)) {
    if (object$type == "regression") return(object$predicted)
    if (proximity & is.null(object$proximity))
      warning("cannot return proximity without new data if random forest object does not already have proximity")
    if (out.type == 1) {
      if (proximity) {
        return(list(pred = object$predicted,
                    proximity = object$proximity))
      } else return(object$predicted)
    }
    if (norm.votes) { 
      t1 <- t(apply(object$votes, 1, function(x) { x/sum(x) }))
      if(proximity) return(list(pred = t1, proximity = object$proximity))
      else return(t1)
    } else {
      if(proximity) return(list(pred = object$votes, proximity = object$proximity))
      else return(object$votes)
    }
  }

  if (object$type == "unsupervised") 
    stop("Can't predict unsupervised forest.")

  if (inherits(object, "randomForest.formula")) {
    newdata <- as.data.frame(newdata)
    rn <- row.names(newdata)
    Terms <- delete.response(object$terms)
    x <- model.frame(Terms, newdata, na.action = na.omit)
    keep <- match(row.names(x), rn)
  } else {
    if (is.null(dim(newdata))) 
      dim(newdata) <- c(1, length(newdata))
    x <- newdata
    if (nrow(x) == 0) 
      stop("newdata has 0 rows")
    if (any(is.na(x))) 
      stop("missing values in newdata")
    keep <- 1:nrow(x)
    rn <- rownames(x)
  }
  if (is.data.frame(x)) {
    for(i in seq(along=ncol(x))) {
      if(is.ordered(x[[i]])) x[[i]] <- as.numeric(x[[i]])
    }
    cat.new <- sapply(x, function(x) if (is.factor(x) && !is.ordered(x)) 
                      length(levels(x))
    else 1)
    if (length(cat.new) != length(object$forest$ncat))
        stop("Number of variables in newdata does not match the model.")
    if (!all(object$forest$ncat == cat.new)) 
      stop("Type of predictors in new data do not match that of the training data.")
  }
  vname <- if (is.null(dim(object$importance))) {
      names(object$importance)
  } else {
      rownames(object$importance)
  }
  if (any(colnames(x) != vname)) stop("names of predictor variables do not match")
  mdim <- ncol(x)
  ntest <- nrow(x)
  ntree <- object$forest$ntree
  maxcat <- object$forest$maxcat
  nclass <- object$forest$nclass
  nrnodes <- object$forest$nrnodes
  show.error <- 0
  x <- t(data.matrix(x))

  if (predict.all) {
    if (object$type == "regression") {
      treepred <- double(ntest * ntree)
    } else {
      treepred <- integer(ntest * ntree)
    }
  } else {
    treepred <- numeric(ntest)
  }
  
  if (proximity) {
    proxmatrix <- matrix(0, ntest, ntest)
  } else {
    proxmatrix <- numeric(1)
  }

  nodexts <- if (nodes) integer(ntest*ntree) else integer(ntest)
  
  if(object$type == "regression") {
    keepIndex <- "ypred"
    if (predict.all) keepIndex <- c(keepIndex, "treepred")
    if (proximity) keepIndex <- c(keepIndex, "proximity")
    ans <- .C("runrforest",
              as.double(x),
              ypred = double(ntest),
              as.integer(mdim),
              as.integer(ntest),
              as.integer(ntree),
              as.integer(object$forest$ndbigtree),
              as.integer(aperm(object$forest$treemap, c(2, 1, 3))),
              as.integer(object$forest$nodestatus),
              as.integer(object$forest$nrnodes),
              as.double(object$forest$xbestsplit),
              as.double(object$forest$nodepred),
              as.integer(object$forest$bestvar),
              as.integer(object$forest$ncat),
              as.integer(predict.all),
              treepred = as.double(treepred),
              as.integer(proximity),
              proximity = as.double(proxmatrix),
              DUP=FALSE,
              PACKAGE = "randomForest")[keepIndex]
    ## Apply bias correction if needed.
    if (!is.null(object$coefs)) {
      yhat <- object$coefs[1] + object$coefs[2] * ans$ypred
    } else {
      yhat <- ans$ypred
    }
    if (predict.all) {
      treepred <- matrix(ans$treepred, length(keep),
                         dimnames=list(rn[keep], NULL))
    }
    if (!proximity) {
      res <- if (predict.all)
        list(aggregate=yhat, individual=treepred) else yhat
    } else {
      res <- list(predicted = yhat, proximity = structure(ans$proximity,
                                     dim=c(ntest, ntest),
                                     dimnames=list(rn, rn)))
    }
  } else {
    countts <- matrix(0, ntest, nclass)
    t1 <- .C("runforest",
             mdim = as.integer(mdim),
             ntest = as.integer(ntest), 
             nclass = as.integer(object$forest$nclass),
             maxcat = as.integer(maxcat), 
             nrnodes = as.integer(nrnodes),
             jbt = as.integer(ntree),
             xts = as.double(x),
             xbestsplit = as.double(object$forest$xbestsplit), 
             pid = as.double(object$forest$pid),
             cutoff = as.double(object$forest$cutoff),
             countts = as.double(countts),
             treemap = as.integer(aperm(object$forest$treemap, 
               c(2, 1, 3))),
             nodestatus = as.integer(object$forest$nodestatus), 
             cat = as.integer(object$forest$ncat),
             cbestsplit = as.integer(numeric(maxcat * nrnodes)),
             nodepred = as.integer(object$forest$nodepred), 
             treepred = as.integer(treepred),
             jet = as.integer(numeric(ntest)), 
             bestvar = as.integer(object$forest$bestvar),
             nodexts = nodexts,
             ndbigtree = as.integer(object$forest$ndbigtree), 
             predict.all = as.integer(predict.all),
             prox = as.integer(proximity),
             proxmatrix = as.double(proxmatrix),
             nodes = as.integer(nodes),
             DUP=TRUE,
             PACKAGE = "randomForest")
    if (out.type > 1) {
      out.class.votes <- t(matrix(t1$countts, nr = nclass, nc = ntest))
      if (norm.votes) 
        out.class.votes <-
          sweep(out.class.votes, 1, rowSums(out.class.votes), "/")
      z <- matrix(NA, ntest, nclass, dimnames = list(rn, levels(object$predicted)))
      z[keep, ] <- out.class.votes
      res <- z
    } else {
      out.class <- factor(rep(NA, length(rn)),
                          levels=1:length(object$classes),
                          labels=object$classes)
      out.class[keep] <- object$classes[t1$jet]
      names(out.class[keep]) <- rn[keep]
      res <- out.class
    }
    if (predict.all) {
      treepred <- matrix(object$classes[t1$treepred],
                         nrow=length(keep), dimnames=list(rn[keep], NULL))
      res <- list(aggregate=res, individual=treepred)
    }
    if(proximity)
      res <- list(predicted = res, proximity = structure(t1$proxmatrix,
                                dim = c(ntest, ntest),
                                dimnames = list(rn[keep], rn[keep])))
    if (nodes) attr(res, "nodes") <- matrix(t1$nodexts, ntest, ntree,
                                            dimnames=list(rn[keep], 1:ntree))
  }
  res
}
"print.randomForest" <-
function(x, ...) {
  cat("\nCall:\n", deparse(x$call), "\n")
  cat("               Type of random forest: ", x$type, "\n", sep="")
  cat("                     Number of trees: ", x$ntree, "\n",sep="")
  cat("No. of variables tried at each split: ", x$mtry, "\n\n", sep="")
  if(x$type == "classification") {
    if(!is.null(x$confusion)) {
      cat("        OOB estimate of  error rate: ",
          round(x$err.rate[x$ntree, "OOB"]*100, dig=2), "%\n", sep="")
      cat("Confusion matrix:\n")
      print(x$confusion)
      if(!is.null(x$test$err.rate)) {
        cat("                Test set error rate: ",
            round(x$test$err.rate[x$ntree, "Test"]*100, dig=2), "%\n",
            sep="")
        cat("Confusion matrix:\n")
        print(x$test$confusion)
      }
    }
  }
  if(x$type == "regression") {
    if(!is.null(x$mse)) {
      cat("          Mean of squared residuals: ", x$mse[length(x$mse)],
          "\n", sep="")
      cat("                    % Var explained: ",
          round(100*x$rsq[length(x$rsq)], dig=2), "\n", sep="")
      if(!is.null(x$test$mse)) {
        cat("                       Test set MSE: ",
            round(x$test$mse[length(x$test$mse)], dig=2), "\n", sep="")
        cat("                    % Var explained: ",
            round(100*x$test$rsq[length(x$test$rsq)], dig=2), "\n", sep="")
      }      
    }
    if (!is.null(x$coefs)) {
      cat("  Bias correction applied:\n")
      cat("  Intercept: ", x$coefs[1], "\n")
      cat("      Slope: ", x$coefs[2], "\n")
    }
  }
}
"randomForest" <-
function(x, ...)
  UseMethod("randomForest")
"randomForest.default" <-
    function(x, y=NULL,  xtest=NULL, ytest=NULL, addclass=0, ntree=500,
             mtry=ifelse(!is.null(y) && !is.factor(y),
             max(floor(ncol(x)/3), 1), 
             floor(sqrt(ncol(x)))), replace=TRUE, classwt=NULL, cutoff,
             sampsize = if (replace) nrow(x) else ceiling(.632*nrow(x)),
             nodesize = if (!is.null(y) && !is.factor(y)) 5 else 1, 
             importance=FALSE, localImp=FALSE,
             proximity=FALSE, oob.prox=proximity,
             outscale=FALSE, norm.votes=TRUE, do.trace=FALSE,
             keep.forest=is.null(xtest), corr.bias=FALSE, ...)
{
    classRF <- is.null(y) || is.factor(y)
    if (!classRF && length(unique(y)) <= 5) {
        warning("The response has five or fewer unique values.  Are you sure you want to do regression?")
    }
    if (classRF && !is.null(y) && addclass==0 && length(unique(y)) < 2) {
        stop("Need at least two classes to do classification.")
    }
    n <- nrow(x)
    p <- ncol(x)
    if (n == 0) stop("data (x) has 0 rows")
    x.row.names <- rownames(x)
    x.col.names <- if (is.null(colnames(x))) 1:ncol(x) else colnames(x)
    
    ## overcome R's lazy evaluation:
    keep.forest <- keep.forest
    
    testdat <- !is.null(xtest)
    if (testdat) {
        if (ncol(x) != ncol(xtest))
            stop("x and xtest must have same number of columns") 
        ntest <- nrow(xtest)
        xts.row.names <- rownames(xtest)
    }
    
    if(mtry > p) {
        mtry <- p
        warning("mtry can not be larger than number of predictors.  Reset to equal to number of predictors")
    }
    if (!is.null(y)) {
        if (length(y) != n) stop("length of response must be the same as predictors")
        addclass <- 0
    } else {
        if (addclass == 0) addclass <- 1
        y <- factor(c(rep(1, n), rep(2, n)))
        x <- rbind(x, x)
        keep.forest <- FALSE
    }
  if (!is.element(addclass,0:2))
      stop("addclass can only take on values 0, 1, or 2")
  
  if (any(is.na(x))) stop("NA not permitted in predictors")
  if (testdat && any(is.na(xtest))) stop("NA not permitted in xtest")
  if (any(is.na(y))) stop("NA not permitted in response")
  if (!is.null(ytest) && any(is.na(ytest))) stop("NA not permitted in ytest")

    if (is.data.frame(x)) {
        ncat <- sapply(x, function(x) if(is.factor(x) && !is.ordered(x))
                       length(levels(x)) else 1)
        x <- data.matrix(x)
        if(testdat) {
            if(!is.data.frame(xtest))
                stop("xtest must be data frame if x is")
            ncatts <- sapply(xtest, function(x) if(is.factor(x) &&
                                                   !is.ordered(x))
                             length(levels(x)) else 1)
            if(!all(ncat == ncatts))
                stop("columns of x and xtest must be the same type")
            xtest <- data.matrix(xtest)
        }
    } else {
        ncat <- rep(1, p)
    }
    maxcat <- max(ncat)
    if (maxcat > 32)
        stop("Can not handle categorical predictors with more than 32 categories.")
    
    if (classRF) {
        nclass <- length(levels(y))
        if (!is.null(ytest)) {
            if (!is.factor(ytest)) stop("ytest must be a factor")
            if (!all(levels(y) == levels(ytest)))
                stop("y and ytest must have the same levels")
        }
        if (missing(cutoff)) {
            cutoff <- rep(1 / nclass, nclass)
        } else {
            if (sum(cutoff) > 1 || sum(cutoff) < 0 || !all(cutoff > 0) ||
                length(cutoff) != nclass) {
                stop("Incorrect cutoff specified.")
            }
            if (!is.null(names(cutoff))) {
                if (!all(names(cutoff) %in% levels(y))) {
                    stop("Wrong name(s) for cutoff")
                }
                cutoff <- cutoff[levels(y)]
            }
        }
        if (!is.null(classwt)) {
            if (length(classwt) != nclass)
                stop("length of classwt not equal to number of classes")
            ## If classwt has names, match to class labels.
            if (!is.null(names(classwt))) {
                if (!all(names(cutoff) %in% levels(y))) {
                    stop("Wrong name(s) for cutoff")
                }
                classwt <- classwt[levels(y)]
            }
            if (any(classwt <= 0)) stop("classwt must be positive")
            ipi <- 1
        } else {
            classwt <- rep(1, nclass)
            ipi <- 0
        }
    } else addclass <- 0  
    
    if(outscale) {
        outscale <- 1
        proximity <- TRUE
        outlier <- rep(0, n)
    } else {
        outlier <- 0
        outscale <- 0
    }
    
    if(proximity) {
        prox <- matrix(0.0, n, n)
        proxts <- if (testdat) matrix(ntest, ntest + n) else double(1)
    } else {
        prox <- proxts <- double(1)
    }

    if (localImp) {
        importance <- TRUE
        impmat <- matrix(0, p, n)
    } else impmat <- double(1)
    
    if (importance) {
        if (classRF) {
            impout <- matrix(0.0, p, nclass + 2)
            impSD <- matrix(0.0, p, nclass + 1)
        } else {
            impout <- matrix(0.0, p, 2)
            impSD <- double(p)
            names(impSD) <- x.col.names
        }
    } else {
        impout <- double(p)
        impSD <- double(1)
    }
    
    nsample <- if (addclass == 0) n else 2*n
    Stratify <- length(sampsize) > 1
    if ((!Stratify) && sampsize > nrow(x)) stop("sampsize too large")
    if (Stratify && (!classRF)) stop("sampsize should be of length one")
    if (classRF) {
        if (Stratify) {
            nsum <- sum(sampsize)
            if (length(sampsize) > nlevels(y))
                stop("sampsize has too many elements.")
            if (any(sampsize <= 0) || nsum == 0)
                stop("Bad sampsize specification")
            ## If sampsize has names, match to class labels.
            if (!is.null(names(sampsize))) {
                sampsize <- sampsize[levels(y)]
            }
            if (any(sampsize > table(y)))
              stop("sampsize can not be larger than class frequency")
        } else {
            nsum <- sampsize
        }
        nrnodes <- 2 * trunc(nsum / nodesize) + 1
    } else {
        ## For regression trees, need to do this to get maximal trees.
        nrnodes <- 2 * trunc(sampsize/max(1, nodesize - 3)) + 1
    }
    

    x <- t(x)
    storage.mode(x) <- "double"
    if (testdat) {
        xtest <- t(xtest)
        storage.mode(xtest) <- "double"
        if (is.null(ytest)) {
            ytest <- labelts <- 0
        } else {
            labelts <- TRUE
        }
    } else {
        xtest <- double(1)
        ytest <- double(1)
        ntest <- 1
        labelts <- FALSE
    }
    nt <- if (keep.forest) ntree else 1

    if (classRF) {
        error.test <- if (labelts) double((nclass+1) * ntree) else double(1)
        rfout <- .C("classRF",
                    x = x,
                    xdim = as.integer(c(p, n)),
                    y = as.integer(y),
                    nclass = as.integer(nclass),
                    ncat = as.integer(ncat), 
                    maxcat = as.integer(maxcat),
                    sampsize = as.integer(sampsize),
                    Options = as.integer(c(addclass,
                    importance,
                    localImp,
                    proximity,
                    oob.prox,
                    outscale,
                    do.trace,
                    keep.forest,
                    replace,
                    Stratify)),
                    ntree = as.integer(ntree),
                    mtry = as.integer(mtry),
                    ipi = as.integer(ipi),
                    classwt = as.double(classwt),
                    cutoff = as.double(cutoff),
                    nodesize = as.integer(nodesize),
                    outlier = as.double(outlier),
                    outcl = integer(nsample),
                    counttr = integer(nclass * nsample),
                    prox = prox,
                    impout = impout,
                    impSD = impSD,
                    impmat = impmat,
                    nrnodes = as.integer(nrnodes),
                    ndbigtree = integer(ntree),
                    nodestatus = integer(nt * nrnodes),
                    bestvar = integer(nt * nrnodes),
                    treemap = integer(nt * 2 * nrnodes),
                    nodepred = integer(nt * nrnodes),
                    xbestsplit = double(nt * nrnodes),
                    pid = double(max(2, nclass)),
                    errtr = double((nclass+1) * ntree),
                    testdat = as.integer(testdat),
                    xts = as.double(xtest),
                    clts = as.integer(ytest),
                    nts = as.integer(ntest),
                    countts = double(nclass * ntest),
                    outclts = as.integer(numeric(ntest)),
                    labelts = as.integer(labelts),
                    proxts = proxts,
                    errts = error.test,
                    DUP=FALSE,
                    PACKAGE="randomForest")[-1]
        if (addclass == 0) {
            if (keep.forest) {
        ## deal with the random forest outputs
                max.nodes <- max(rfout$ndbigtree)
                treemap <- array(rfout$treemap, dim = c(2, nrnodes, ntree))
                treemap <- aperm(treemap, c(2,1,3))[1:max.nodes, , ,drop=FALSE]
            }
            ## Turn the predicted class into a factor like y.
            out.class <- factor(rfout$outcl, levels=1:nclass,
                                label=levels(y))
            names(out.class) <- x.row.names
            con <- table(observed = y,
                         predicted = out.class)[levels(y), levels(y)]
            con <- cbind(con, class.error = 1 - diag(con)/rowSums(con))
        }
        out.votes <- t(matrix(rfout$counttr, nclass, nsample))[1:n, ]
        oob.times <- rowSums(out.votes)
        if(norm.votes) 
            out.votes <- t(apply(out.votes, 1, function(x) x/sum(x)))
        dimnames(out.votes) <- list(x.row.names, levels(y))
        if(testdat) {
            out.class.ts <- factor(rfout$outclts, levels=1:nclass,
                                   label=levels(y))
            names(out.class.ts) <- xts.row.names
            out.votes.ts <- t(matrix(rfout$countts, nclass, ntest))
            dimnames(out.votes.ts) <- list(xts.row.names, levels(y))
            if (norm.votes)
                out.votes.ts <- t(apply(out.votes.ts, 1,
                                        function(x) x/sum(x)))
            if (labelts) {
                testcon <- table(observed = ytest,
                                 predicted = out.class.ts)[levels(y), levels(y)]
                testcon <- cbind(testcon,
                                 class.error = 1 - diag(testcon)/rowSums(testcon))
            }
        }
        out <- list(call = match.call(),
                    type = ifelse(addclass == 0, "classification",
                    "unsupervised"),
                    predicted = if (addclass == 0) out.class else NULL,
                    err.rate = if (addclass == 0) t(matrix(rfout$errtr,
                                   nclass+1,
                                   ntree, dimnames=list(c("OOB", levels(y)),
                                          NULL))) else NULL, 
                    confusion = if(addclass == 0) con else NULL,
                    votes = out.votes,
                    oob.times = oob.times,
                    classes = levels(y),
                    importance = if (importance) 
                    matrix(rfout$impout, p, nclass+2,
                           dimnames = list(x.col.names,
                           c(levels(y), "MeanDecreaseAccuracy",
                             "MeanDecreaseGini")))
                    else structure(rfout$impout, names=x.col.names),
                    importanceSD = if (importance)
                    matrix(rfout$impSD, p, nclass + 1,
                           dimnames = list(x.col.names,
                           c(levels(y), "MeanDecreaseAccuracy")))
                    else NULL,
                    localImportance = if (localImp)
                    matrix(rfout$impmat, p, n,
                           dimnames = list(x.col.names,x.row.names)) else NULL,
                    proximity = if (proximity) matrix(rfout$prox, n, n,
                    dimnames = list(x.row.names, x.row.names)) else NULL,
                    outlier = if (outscale) rfout$outlier else NULL,
                    ntree = ntree,
                    mtry = mtry,
                    forest = if (addclass > 0 || !keep.forest) NULL else {
                        list(ndbigtree = rfout$ndbigtree, 
                             nodestatus = matrix(rfout$nodestatus,
                             nc = ntree)[1:max.nodes,],
                             bestvar = matrix(rfout$bestvar, nc = ntree)[1:max.nodes,],
                             treemap = treemap,
                             nodepred = matrix(rfout$nodepred,
                             nc = ntree)[1:max.nodes,],
                             xbestsplit = matrix(rfout$xbestsplit,
                             nc = ntree)[1:max.nodes,],
                             pid = rfout$pid, cutoff = cutoff, ncat = ncat, maxcat = maxcat, 
                             nrnodes = max.nodes, ntree = ntree, nclass = nclass)
                    },
                    test = if(!testdat) NULL else list(
                    predicted = out.class.ts,
                    err.rate = if (labelts) t(matrix(rfout$errts, nclass+1,
                    ntree,
                dimnames=list(c("Test", levels(y)), NULL))) else NULL,
                    confusion = if (labelts) testcon else NULL,
                    votes = out.votes.ts,
                    proximity = if(proximity) matrix(rfout$proxts, nrow=ntest,
                    dimnames = list(xts.row.names, c(xts.row.names,
                    x.row.names))) else NULL))
    } else {
        rfout <- .C("regRF",
                    x,
                    as.double(y),
                    as.integer(n),
                    as.integer(p),
                    as.integer(sampsize),
                    as.integer(nodesize),
                    as.integer(nrnodes),
                    as.integer(ntree),
                    as.integer(mtry),
                    as.integer(c(importance, localImp)),
                    as.integer(ncat),
                    as.integer(do.trace),
                    as.integer(proximity),
                    as.integer(oob.prox),
                    as.integer(corr.bias),
                    ypred = double(n),
                    impout = impout,
                    impmat = impmat,
                    impSD = impSD,
                    prox = prox,
                    ndbigtree = integer(ntree),
                    nodestatus = integer(nt * nrnodes),
                    treemap = integer(nt * 2 * nrnodes),
                    nodepred = double(nt * nrnodes),
                    bestvar = integer(nt * nrnodes),
                    xbestsplit = double(nt * nrnodes),
                    mse = double(ntree),
                    keepf = as.integer(keep.forest),
                    replace = as.integer(replace),
                    testdat = as.integer(testdat),
                    xts = xtest,
                    ntest = as.integer(ntest),
                    yts = as.double(ytest),
                    labelts = as.integer(labelts),
                    ytestpred = double(ntest),
                    proxts = proxts,
                    msets = double(if (labelts) ntree else 1),
                    coef = double(2),
                    oob.times = integer(n),
                    DUP=FALSE,
                    PACKAGE="randomForest")[c(16:27, 35:39)]
        ## Format the forest component, if present.
        if (keep.forest) {
            max.nodes <- max(rfout$ndbigtree)
            
            rfout$nodestatus <-
              matrix(rfout$nodestatus, ncol = ntree)[1:max.nodes,,drop=FALSE]
            rfout$bestvar <-
              matrix(rfout$bestvar, ncol = ntree)[1:max.nodes,,drop=FALSE]
            rfout$nodepred <-
              matrix(rfout$nodepred, ncol = ntree)[1:max.nodes,,drop=FALSE]
            rfout$xbestsplit <-
              matrix(rfout$xbestsplit, ncol = ntree)[1:max.nodes,,drop=FALSE]
            rfout$treemap <- aperm(array(rfout$treemap,
                                         dim = c(2, nrnodes, ntree)),
                                   c(2, 1, 3))[1:max.nodes, , ,drop=FALSE]
        }
        
        out <- list(call = match.call(),
                    type = "regression",
                    predicted = structure(rfout$ypred, names=x.row.names),
                    mse = rfout$mse,
                    rsq = 1 - rfout$mse / (var(y) * (n-1) / n),
                    oob.times = rfout$oob.times,
                    importance = if (importance) matrix(rfout$impout, p, 2,
                    dimnames=list(x.col.names,
                                  c("%IncMSE","IncNodePurity"))) else
                        structure(rfout$impout, names=x.col.names),
                    importanceSD=if (importance) rfout$impSD else NULL,
                    localImportance = if (importance)
                    matrix(rfout$impmat, p, n, dimnames=list(x.col.names,
                                               x.row.names)) else NULL,
                    proximity = if (proximity) matrix(rfout$prox, n, n,
                    dimnames = list(x.row.names, x.row.names)) else NULL,
                    ntree = ntree,
                    mtry = mtry,
                    forest = if (keep.forest)
                    c(rfout[c("ndbigtree", "nodestatus", "treemap",
                              "nodepred", "bestvar", "xbestsplit")],
                      list(ncat = ncat), list(nrnodes=max.nodes),
                      list(ntree=ntree)) else NULL,
                    coefs = if (corr.bias) rfout$coef else NULL,
                    test = if(testdat) {
                        list(predicted = structure(rfout$ytestpred,
                             names=xts.row.names),
                             mse = if(labelts) rfout$msets else NULL,
                             rsq = if(labelts) 1 - rfout$msets /
                                        (var(ytest)*(n-1)/n) else NULL,
                             proximity = if (proximity) 
                             matrix(rfout$proxts / ntree, nrow = ntest,
                                    dimnames = list(xts.row.names,
                                    c(xts.row.names,
                                    x.row.names))) else NULL)
                    } else NULL)
    }
    class(out) <- "randomForest"
    return(out)
}
"randomForest.formula" <-
    function(x, data = NULL, ..., subset, na.action = na.fail) {	
### formula interface for randomForest.
### code gratefully stolen from svm.formula (package e1071).
###
    if (!inherits(x, "formula"))
        stop("method is only for formula objects")
    call <- match.call()
    m <- match.call(expand = FALSE)
    names(m)[2] <- "formula"
    if (is.matrix(eval(m$data, parent.frame())))
        m$data <- as.data.frame(data)
    m$... <- NULL
    m$na.action <- na.action
    m[[1]] <- as.name("model.frame")
    m <- eval(m, parent.frame())
    Terms <- attr(m, "terms")
    attr(Terms, "intercept") <- 0
    y <- model.response(m)
    if(!is.null(y)) m <- m[,-1]
    for (i in seq(along=ncol(m))) {
        if(is.ordered(m[[i]])) m[[i]] <- as.numeric(m[[i]])
    }
    ret <- randomForest.default(m, y, ...)
    ret$terms <- Terms
    cl <- match.call()
    cl[[1]] <- as.name("randomForest")
    ret$call <- cl

    if (!is.null(attr(m, "na.action"))) 
        ret$na.action <- attr(m, "na.action")
    class(ret) <- c("randomForest.formula", "randomForest")
    return(ret)
}
rfImpute <- function(x, ...)
    UseMethod("rfImpute")

rfImpute.formula <- function(x, data, ..., subset) {
    if (!inherits(x, "formula"))
        stop("method is only for formula objects")
    call <- match.call()
    m <- match.call(expand = FALSE)
    names(m)[2] <- "formula"
    if (is.matrix(eval(m$data, parent.frame())))
        m$data <- as.data.frame(data)
    m$... <- NULL
    m$na.action <- as.name("na.pass")
    m[[1]] <- as.name("model.frame")
    m <- eval(m, parent.frame())
    Terms <- attr(m, "terms")
    attr(Terms, "intercept") <- 0
    y <- model.response(m)
    if (!is.null(y)) m <- m[,-1]
    for (i in seq(along=ncol(m))) {
        if(is.ordered(m[[i]])) m[[i]] <- as.numeric(m[[i]])
    }
    ret <- rfImpute.default(m, y, ...)
    names(ret)[1] <- deparse(as.list(x)[[2]])
    ret
}

rfImpute.default <- function(x, y, iter=5, ntree=300, ...) {
    if (any(is.na(y))) stop("Can't have NAs in", deparse(substitute(y)))
    if (!any(is.na(x))) stop("No NAs found in ", deparse(substitute(x)))
    xf <- na.roughfix(x)
    hasNA <- which(apply(x, 2, function(x) any(is.na(x))))
    if (is.data.frame(x)) {
        isfac <- sapply(x, is.factor)
    } else {
        isfac <- rep(FALSE, ncol(x))
    }
    
    for (i in 1:iter) {
        prox <- randomForest(xf, y, ntree=ntree, ..., do.trace=ntree,
                             proximity=TRUE)$proximity
        for (j in hasNA) {
            miss <- which(is.na(x[, j]))
            if (isfac[j]) {
                lvl <- levels(x[[j]])
                catprox <- apply(prox[-miss, miss, drop=FALSE], 2,
                                 function(v) lvl[which.max(tapply(v, x[[j]][-miss], mean))])
                xf[miss, j] <- catprox
            } else {
                sumprox <- colSums(prox[-miss, miss, drop=FALSE])
                xf[miss, j] <- (prox[miss, -miss, drop=FALSE] %*% xf[,j][-miss]) / sumprox
            }
            NULL
        }
    }
    xf <- cbind(y, xf)
    names(xf)[1] <- deparse(substitute(y))
    xf
}
rfNews <- function() {
    newsfile <- file.path(system.file(package="randomForest"), "NEWS")
    file.show(newsfile)
}
treesize <- function(x, terminal=TRUE) {
  if(!inherits(x, "randomForest"))
    stop("This function only works for objects of class `randomForest'")
  if(is.null(x$forest)) stop("The object must contain the forest component")
  if(terminal) return((x$forest$ndbigtree+1)/2) else return(x$forest$ndbigtree)
}
tuneRF <- function(x, y, mtryStart=if(is.factor(y)) floor(sqrt(ncol(x))) else
                   floor(ncol(x)/3), ntreeTry=50, stepFactor=2,
                   improve=0.05, trace=TRUE, plot=TRUE, doBest=FALSE, ...) {
  if (improve < 0) stop ("improve must be non-negative.")
  classRF <- is.factor(y)
  errorOld <- if (classRF) {
    randomForest(x, y, mtry=mtryStart, ntree=ntreeTry,
                 keep.forest=FALSE, ...)$err.rate[ntreeTry,1]
  } else {
    randomForest(x, y, mtry=mtryStart, ntree=ntreeTry,
                 keep.forest=FALSE, ...)$mse[ntreeTry]
  }
  if (trace) {
    cat("mtry =", mtryStart, " OOB error =",
        if (classRF) paste(100*round(errorOld, 4), "%", sep="") else
        errorOld, "\n")
  }

  oobError <- list()
  oobError[[as.character(mtryStart)]] <- errorOld  
  
  for (direction in c("left", "right")) {
    if (trace) cat("Searching", direction, "...\n")
    Improve <- 1.1*improve
    mtryBest <- mtryStart
    mtryCur <- mtryStart
    while (Improve >= improve) {
      mtryOld <- mtryCur
      mtryCur <- if (direction == "left") {
        max(1, ceiling(mtryCur / stepFactor))
      } else {
        min(ncol(x), floor(mtryCur * stepFactor))
      }
      if (mtryCur == mtryOld) break
      errorCur <- if (classRF) {
        randomForest(x, y, mtry=mtryCur, ntree=ntreeTry,
                     keep.forest=FALSE, ...)$err.rate[ntreeTry,"OOB"]
      } else {
        randomForest(x, y, mtry=mtryCur, ntree=ntreeTry,
                     keep.forest=FALSE, ...)$mse[ntreeTry]
      }
      if (trace) {
        cat("mtry =",mtryCur, "\tOOB error =",
            if (classRF) paste(100*round(errorCur, 4), "%", sep="") else
            errorCur, "\n")
      }
      oobError[[as.character(mtryCur)]] <- errorCur
      Improve <- 1 - errorCur/errorOld
      cat(Improve, improve, "\n")
      if (Improve > improve) {
        errorOld <- errorCur
        mtryBest <- mtryCur
      }
    }
  }
  res <- unlist(oobError)[order(as.numeric(names(oobError)))]
  res <- cbind(mtry=as.numeric(names(res)), OOBError=res)

  if (plot) {
    plot(res, xlab=expression(m[try]), ylab="OOB Error", type="o", log="x",
         xaxt="n")
    axis(1, at=res[,"mtry"])
  }

  if (doBest) 
    res <- randomForest(x, y, mtry=res[which.min(res[,2]), 1], ...)
  
  res
}
varImpPlot <- function(x, sort=TRUE,
                         n.var=min(30, if(is.null(dim(x$importance)))
                           length(x$importance) else nrow(x$importance)),
                         class = NULL, scale=TRUE, xlab="Importance", ylab="",
                         main=deparse(substitute(x)), ...) {
    if (!inherits(x, "randomForest"))
        stop("This function only works for objects of class `randomForest'")
    ## If only impurity-based measures exists, or only class-specific
    ## measure is requested, then only one plot to draw.
    if (is.null(dim(x$importance)) || !is.null(class)) {
        if (is.null(class)) {
            imp <- x$importance
        } else {
            imp <- x$importance[,class]
            if (scale) {
                SD <- x$importanceSD[,class]
                imp <- ifelse(SD < .Machine$double.eps, 0, imp / SD)
            }
        }
        if (sort) {
            ord <- order(imp, decreasing=TRUE)[1:n.var]
            imp <- imp[ord]
            dotchart(rev(imp), xlab=xlab, ylab=ylab, main=main, ...)
        } else {
            dotchart(imp, xlab="Importance", ylab="", main=main, ...)
        }
    } else {
        ## Extract the last two columns.
        impmat <- x$importance[,(ncol(x$importance)-1):ncol(x$importance)]
        if (scale) {
            SD <- if (is.null(dim(x$importanceSD))) x$importanceSD else
                   x$importanceSD[,ncol(x$importanceSD)]
            impmat[,1] <- ifelse(SD < .Machine$double.eps, 0, impmat[,1] / SD)
        }
        imp <- vector(2, mode="list")
        names(imp) <- colnames(impmat)
        op <- par(mfrow=c(1, 2), mar=c(4, 5, 4, 1), mgp=c(2, .8, 0),
                  oma=c(0,0,2,0))
        on.exit(par(op))
        
        for (i in 1:2) {
            if(sort) {
                ord <- order(impmat[,i], decreasing=TRUE)[1:n.var]
                imp[[i]] <- impmat[ord, i]
                maximp <- max(imp[[i]])
                dotchart(rev(imp[[i]]), xlab=xlab,
                         ylab=ylab, main=names(imp)[i], xlim=c(0, maximp),
                         ...)
            } else {
                imp[[i]] <- impmat[, i]
                dotchart(imp[[i]], xlab=xlab, ylab=ylab,
                         main=colnames(imp)[i], ...)
            }
        }
        mtext(outer=TRUE, side=3, text=main, cex=1.1)
    }
    invisible(imp)
}
varUsed <- function(x, by.tree=FALSE, count=TRUE) {
    if (!inherits(x, "randomForest"))
        stop(deparse(substitute(x)), "is not a randomForest object")
    if (is.null(x$forest))
        stop(deparse(substitute(x)), "does not contain forest")
    
    p <- length(x$forest$ncat)  # Total number of variables.
    if (count) {
        if (by.tree) {
            v <- apply(x$forest$bestvar, 2, function(x) {
                xx <- numeric(p)
                y <- table(x[x>0])
                xx[as.numeric(names(y))] <- y
                xx
            })
        } else {
            v <- numeric(p)
            vv <- table(x$forest$bestvar[x$forest$bestvar > 0])
            v[as.numeric(names(vv))] <- vv
        }
    } else {
        v <- apply(x$forest$bestvar, 2, function(x) sort(unique(x[x>0])))
        if(!by.tree) v <- sort(unique(unlist(v)))
    }
    v
}
