# SCCS  %W% %G%
# format a set of numbers using C's "g" format
#  It is applied on an element by element basis, which is more
#  appropriate for rpart output than the standard Splus format()
#  command.
# For instance if x=(123, 1.23, .00123)
#	  format(x) = "123.00000", "1.23000", "0.00123"
#  but formatg does not add all of those zeros to the first two numbers
#
formatg <- function(x, digits = getOption("digits"),
                    format = paste("%.", digits, "g", sep=''))
{
    if (!is.numeric(x)) stop("x must be a numeric vector")

    n <- length(x)
    #
    # the resultant strings could be up to 8 characters longer,
    #   assume that digits =4,  -0.dddde+104 is a worst case, where
    #   dddd are the 4 significant digits.
    dummy  <- paste(rep(" ", digits+8), collapse='')
    .C("formatg", as.integer(n),
                  as.double(x),
                  rep(format,n),
                  out= rep(dummy, n))$out
}
# SCCS %W% %G%
# Make the nice labels used by print and summary
#   digits = obvious
#   minlength = 0 = don't abbrev factors
#               1 = use single letters
#               2+= the same arg as the "abbreviate" function
#   collapse = an oddly named argument
#              F = return a matrix with two columns, containing the labels of
#                    the left and right descendants of each node
#              T = return a vector of 1 column, with the label of the parent
#   pretty: for historical compatability
#               0   -> minlength =0
#              NULL -> minlength =1
#               T   -> minlength =4
#   ... = other args for abbreviate()
#
labels.rpart <- function(object, digits=4, minlength=1, pretty,
			      collapse=T, ...) {
    if (missing(minlength) && !missing(pretty)) {
	if (is.null(pretty)) minlength <-1
	else if (is.logical(pretty)) {
	    if (pretty) minlength <- 4
	    else        minlength <- 0
	    }
	else minlength <- 0
	}

    ff <- object$frame
    n  <- nrow(ff)
    is.leaf <- (ff$var == "<leaf>")
    whichrow <- !is.leaf
    vnames <- ff$var[whichrow]  #the variable names for the primary splits

    index <- cumsum(c(1, ff$ncompete + ff$nsurrogate + 1*(!is.leaf)))
    irow  <- index[c(whichrow,F)]     #we only care about the primary split
    ncat  <- object$splits[irow, 2]

    # Now to work: first create labels for the left and right splits,
    #  but not for leaves of course
    #
    lsplit <- rsplit <- vector(mode='character', length= length(irow))

    if (any(ncat <2)) {  # any continuous vars ?
	jrow <- irow[ncat <2]
	cutpoint <- formatg(object$splits[jrow,4], digits)
	temp1 <- (ifelse(ncat<0, "<", ">"))[ncat <2]
	temp2 <- (ifelse(ncat<0, ">", "<"))[ncat <2]
	lsplit[ncat<2] <- paste(temp1, cutpoint, sep='')
	rsplit[ncat<2] <- paste(temp2, cutpoint, sep='')
	}

    if (any(ncat >1)) { # any categorical variables ?
	xlevels <- attr(object, 'xlevels')
	#
	# jrow will be the row numbers of factors within lsplit and rsplit
	# crow the row number in "csplit"
	# and cindex the index on the "xlevels" list
	#
	jrow <- (seq(along=ncat))[ncat>1]
	crow <- object$splits[irow[ncat>1],4]    #row number in csplit
	cindex <- (match(vnames, names(xlevels)))[ncat >1]

	# Now, abbreviate the levels
	if (minlength ==1) {
	    if (any(ncat>52))
		warning(paste("More than 52 levels in a predicting factor,",
			      "truncated for printout"))
	    xlevels <- lapply(xlevels,
			       function(z) {
				   k <- length(z)
				   k <- pmin(1:k, 52)
				   c(letters, LETTERS)[k]
				   })
	    }
	else if (minlength >1)
	    xlevels <- lapply(xlevels, abbreviate, minlength, ...)

	# Now tuck in the labels
	# I'll let some other clever person vectorize this
	for (i in 1:(length(jrow))) {
	    j <- jrow[i]
	    splits <- object$csplit[crow[i],]
	    # splits will contain 1=left, 2=right, 3= neither
	    ltemp <- (1:length(splits))[splits== 1]
	    rtemp <- (1:length(splits))[splits== 3]
	    if (minlength==1) {
		lsplit[j] <- paste((xlevels[[cindex[i]]])[ltemp], collapse='')
		rsplit[j] <- paste((xlevels[[cindex[i]]])[rtemp], collapse='')
		}
	    else {
		lsplit[j] <-paste((xlevels[[cindex[i]]])[ltemp], collapse=',')
		rsplit[j] <-paste((xlevels[[cindex[i]]])[rtemp], collapse=',')
		}
	    }
	}

    if (!collapse) {  #called by no routines that I know of
	ltemp <- rtemp <- rep("<leaf>", n)
	ltemp[whichrow] <- lsplit
	rtemp[whichrow] <- rsplit
	return(cbind(ltemp, rtemp))
	}

    lsplit <- paste(ifelse(ncat<2, "", "="), lsplit, sep='')
    rsplit <- paste(ifelse(ncat<2, "", "="), rsplit, sep='')

    #
    # Now match them up to node numbers
    #   The output will have one label per row of object$frame, each
    #   corresponding the the line segement joining this node to its parent
    #
    varname <- (as.character(vnames))
    node <- as.numeric(row.names(ff))
    parent <- match(node %/% 2, node[whichrow])
    odd <- (as.logical(node %%2))

    labels <- vector('character', length=n)
    labels[odd] <- paste(varname[parent[odd]], rsplit[parent[odd]], sep="")
    labels[!odd]<- paste(varname[parent[!odd]],lsplit[parent[!odd]], sep="")
    labels[1] <- "root"
    labels
    }





# SCCS 02/18/97 @(#)meanvar.rpart.s	1.2

meanvar.rpart <- function(tree, xlab = "ave(y)", ylab = "ave(deviance)", ...)

{
	if(!inherits(tree, "rpart"))
		stop("Not legitimate rpart object")
	if(!tree$method=='anova')
		stop("Plot not useful for classification or poisson trees")
	frame <- tree$frame
	frame <- frame[frame$var == "<leaf>",  ]
	x <- frame$yval
	y <- frame$dev/frame$n
	label <- row.names(frame)
	plot(x, y, xlab = xlab, ylab = ylab, type = "n", ...)
	text(x, y, label)
	invisible(list(x = x, y = y, label = label))
}

meanvar <- function(x,...) UseMethod('meanvar')
# sccs @(#)model.frame.rpart.s	1.3 01/21/97
model.frame.rpart <- function(formula, ...)
{
	m <- formula$model
	if(!is.null(m))
		return(m)
	oc <- formula$call
	if(substring(deparse(oc[[1]]), 1, 7) == "predict") {
		m <- eval(oc$newdata)
		if(is.null(attr(m, "terms"))) {
			object <- eval(oc$object)
			m <- model.frame(object$terms, m, na.rpart)
		}
		return(m)
	}
	while(deparse(oc[[1]]) != "rpart") oc <- eval(oc[[2]])$call
	oc$subset <- names(formula$where)
	oc$method <- formula$method
	eval(oc)
}

#SCCS  @(#)na.rpart.s	1.5 12/13/99
na.rpart <- function(x){
    Terms <- attr(x, 'terms')
    if(!is.null(Terms)) yvar <- attr(Terms, "response") else yvar <- 0
    if (yvar==0) {
	xmiss <- is.na(x)
	keep <-  (xmiss %*% rep(1,ncol(xmiss))) < ncol(xmiss)
	}
    else {
	xmiss <- is.na(x[-yvar])
	ymiss <- is.na(x[[yvar]])
	if (is.matrix(ymiss))
	    keep <- ((xmiss %*% rep(1,ncol(xmiss))) < ncol(xmiss)) &
		    ((ymiss %*% rep(1,ncol(ymiss))) == 0 )
	else
	    keep <- ((xmiss %*% rep(1,ncol(xmiss))) < ncol(xmiss)) & !ymiss
	}
    if (all(keep)) x
    else {
	temp <- seq(keep)[!keep]
	names(temp) <- row.names(x)[!keep]
	#the methods for this group are all the same as for na.omit
	class(temp) <- c("na.rpart", "omit")
	structure(x[keep,], na.action=temp)
	}
    }
## submitted by Anantha Prasad 1/26/98

path.rpart <- function(tree, nodes, pretty = 0, print.it = TRUE)
{
        if(!inherits(tree, "rpart"))
                stop("Not legitimate tree")
        splits <- labels.rpart(tree, pretty = pretty)
        frame <- tree$frame
        n <- row.names(frame)
        node <- as.numeric(n)
        which <- descendants(node)      #ancestors are columns
        path <- list()
        if(missing(nodes)) {
                xy <- rpartco(tree)
                while(length(i <- identify(xy, n = 1, plot = FALSE)) > 0) {
                        path[[n[i]]] <- path.i <- splits[which[, i]]
                        if(print.it) {
                                cat("\n", "node number:", n[i], "\n")
                                cat(paste("  ", path.i), sep = "\n")
                        }
                }
        }
        else {

                if(length(nodes <- node.match(nodes, node)) == 0)
                        return(invisible())
                for(i in nodes)
                       { path[[n[i]]] <- path.i <- splits[which[, i]]
			if(print.it) {
                                cat("\n", "node number:", n[i], "\n")
                                cat(paste("  ", path.i), sep = "\n")
                                }
                       }
        }
        invisible(path)
}





#SCCS @(#)plot.rpart.s	1.6 02/24/00
plot.rpart <- function(tree, uniform=FALSE, branch=1, compress=FALSE,
			     nspace, margin=0, minbranch=.3, ...){
    if(!inherits(tree, "rpart"))
	    stop("Not an rpart object")

    if (compress & missing(nspace)) nspace <- branch
    if (!compress) nspace <- -1     #means no compression
    dev <- dev.cur()
    if (dev == 1) dev <- 2
    assign(paste(".rpart.parms", dev, sep = "."),
            list(uniform=uniform, branch=branch, nspace=nspace,
		 minbranch=minbranch), envir=.GlobalEnv)

    #define the plot region
    temp <- rpartco(tree)
    x <- temp$x
    y <- temp$y
    temp1 <- range(x) + diff(range(x))*c(-margin, margin)
    temp2 <- range(y) + diff(range(y))*c(-margin, margin)
    plot(temp1, temp2, type='n', axes=FALSE, xlab='', ylab='', ...)

    # Draw a series of horseshoes or V's, left son, up, down to right son
    #   NA's in the vector cause lines() to "lift the pen"
    node <- as.numeric(row.names(tree$frame))
    temp <- rpart.branch(x, y, node, branch)

    if (branch>0) text(x[1], y[1], '|')
    lines(c(temp$x), c(temp$y))
    invisible(list(x=x, y=y))
}





# SCCS @(#)plotcp.s	1.1 02/08/98
# Contributed by B.D. Ripley 97/07/17
#
plotcp <- function(x, minline=TRUE, lty=3, col=1,
		   upper=c("size", "splits", "none"), ...)
{
  if(!inherits(x, "rpart")) stop("Not legitimate rpart object")
  upper <- match.arg(upper)
  p.rpart <- x$cptable
  if(ncol(p.rpart) < 5)
    stop("cptable does not contain cross-validation results")
  xstd <- p.rpart[, 5]
  xerror <- p.rpart[, 4]
  nsplit <- p.rpart[, 2]
  ns <- seq(along=nsplit)
  cp0 <- p.rpart[ ,1]
  cp <- sqrt(cp0 * c(Inf, cp0[-length(cp0)]))
  ylim <- c(min(xerror - xstd) - 0.1, max(xerror + xstd) + 0.1)
  plot(ns, xerror, axes = FALSE, xlab = "cp", ylab =
       "X-val Relative Error", ylim = ylim, type = "o", ...)
  box()
  axis(2, ...)
  segments(ns, xerror - xstd, ns, xerror + xstd)
  axis(1, at = ns, lab = as.character(signif(cp, 2)), ...)
  switch(upper,
	 size = {
           axis(3, at = ns, lab = as.character(nsplit+1), ...)
           mtext("size of tree", side=3, line=3)
	 },
	 splits = {
           axis(3, at = ns, lab = as.character(nsplit), ...)
           mtext("number of splits", side=3, line=3)
	 },)
  minpos <- min(seq(along=xerror)[xerror==min(xerror)])
  if(minline) abline(h=(xerror+xstd)[minpos], lty=lty, col=col)
  invisible()
}
# SCCS 02/07/00 @(#)post.rpart.s	1.12
#
post.rpart <- function(tree, title.,
		       filename=paste(deparse(substitute(tree)),".ps",sep=""),
		       digits=getOption("digits") - 3, pretty=TRUE,
		       use.n=TRUE,  horizontal=TRUE, ...)
{
    if(filename !=""){
	postscript(file = filename, horizontal=horizontal, ...)
	par(mar=c(2,2,4,2)+.1)
	on.exit(dev.off())
	}
    else {
	oldpar <- par(mar=c(2,2,4,2)+.1)
	on.exit(invisible(par(oldpar)))
	}

    plot(tree, uniform=TRUE, branch=.2, compress=TRUE, margin=.1)
    text(tree, all=TRUE, use.n=use.n, fancy=TRUE, digits=digits, pretty=pretty)
    method <- tree$method

    if(missing(title.)) {
        temp  <- attr(tree$terms,'variables')[2]
        title(paste("Endpoint =",temp),cex=.8)
    } else if (title. !="") title(title.,cex=.8)
}

## SCCS @(#)post.s	1.3 02/27/98
post <- function(tree, ...) UseMethod("post")

# SCCS @(#)pred.rpart.s	1.3 09/03/97
#
# Do Rpart predictions given a tree and a matrix of predictors
pred.rpart <- function(fit, x) {

    frame <- fit$frame
    nc <- frame[, c('ncompete', 'nsurrogate')]
    frame$index <- 1 + c(0, cumsum((frame$var != "<leaf>") +
                                       nc[[1]] + nc[[2]]))[-(nrow(frame)+1)]
    frame$index[frame$var == "<leaf>"] <- 0
    vnum <- match(dimnames(fit$split)[[1]], dimnames(x)[[2]])
    if (any(is.na(vnum))) stop("Tree has variables not found in new data")
    temp <- .C("pred_rpart",
		    as.integer(dim(x)),
		    as.integer(dim(frame)[1]),
		    as.integer(dim(fit$splits)),
		    as.integer(if(is.null(fit$csplit)) rep(0,2)
		               else dim(fit$csplit)),
		    as.integer(row.names(frame)),
		    as.integer(unlist(frame[,
			     c('n', 'ncompete', 'nsurrogate', 'index')])),
		    as.integer(vnum),
		    as.double(fit$splits),
		    as.integer(fit$csplit -2),
		    as.integer((fit$control)$usesurrogate),
		    as.double(x),
		    as.integer(is.na(x)),
		    where = integer(dim(x)[1]),
		    NAOK =TRUE)
    temp <- temp$where
    names(temp) <- dimnames(x)[[1]]
    temp
    }
## SCCS %W% %G%
predict.rpart <-
function(object, newdata = list(),
	 type = c("vector", "tree", "class", "matrix")) {
    if(!inherits(object, "rpart"))
	stop("Not legitimate tree")
    type <- match.arg(type)
    if(missing(newdata) & type == "tree")
	return(object)  #idiot proofing
    if(missing(newdata))
	where <- object$where
    else {
	if(is.null(attr(newdata, "terms"))) {
	    Terms <- delete.response(object$terms)
	    act <- (object$call)$na.action
	    if (is.null(act)) act<- na.rpart
	    newdata <- model.frame(Terms, newdata, na.action = act,
                                      xlev=attr(object, "xlevels"))
	    }
	where <- pred.rpart(object, rpart.matrix(newdata))
	}
    frame <- object$frame
    method <- object$method
    ylevels <- attr(object,'ylevels')
    if(type == "vector" || (type=='matrix' && is.null(frame$yval2))) {
	pred <- frame$yval[where]
	names(pred) <- names(where)
	}
    else if (type == 'matrix') {
	pred <- frame$yval2[where,]
	dimnames(pred) <- list(names(where), NULL)
	}
    else if(type == "class") {
	if(length(ylevels) == 0)
		stop("Type class is only appropriate for classification")
	pred <- factor(ylevels[frame$yval[where]], levels=ylevels)
	names(pred) <- names(where)
	}
    else stop("Invalid prediction for rpart objects")

    #Expand out the missing values in the result
    # But only if operating on the original dataset
    if (missing(newdata) && !is.null(object$na.action)) {
	    pred <- naresid(object$na.action, pred)
	    }
    pred
    }

#SCCS  %W% %G%
print.rpart <- function(x, minlength=0, spaces=2, cp,
               digits=getOption("digits"), ...) {
    if(!inherits(x, "rpart")) stop("Not legitimate rpart object")

    if (!missing(cp)) x <- prune.rpart(x, cp=cp)
    frame <- x$frame
    ylevel <- attr(x,'ylevels')
    node <- as.numeric(row.names(frame))
    depth <- tree.depth(node)
    indent <- paste(rep(" ", spaces * 32), collapse = "")
    #32 is the maximal depth
    if(length(node) > 1) {
        indent <- substring(indent, 1, spaces * seq(depth))
        indent <- paste(c("", indent[depth]), format(node), ")", sep = "")
        }
    else indent <- paste(format(node), ")", sep = "")

    tfun <- (x$functions)$print
    if (!is.null(tfun)) {
	if (is.null(frame$yval2))
		yval <- tfun(frame$yval,  ylevel, digits)
	else    yval <- tfun(frame$yval2,  ylevel, digits)
	}
    else yval <- format(signif(frame$yval, digits = digits))
    term <- rep(" ", length(depth))
    term[frame$var == "<leaf>"] <- "*"
    z <- labels(x, digits=digits, minlength=minlength, ...)
    n <- frame$n
    z <- paste(indent, z, n, format(signif(frame$dev, digits = digits)),
            yval, term)

    omit <- x$na.action
    if (length(omit))
    cat("n=", n[1], " (", naprint(omit), ")\n\n", sep="")
    else cat("n=", n[1], "\n\n")

    #This is stolen, unabashedly, from print.tree
    if (x$method=='class')
         cat("node), split, n, loss, yval, (yprob)\n")
    else cat("node), split, n, deviance, yval\n")
    cat("      * denotes terminal node\n\n")

    cat(z, sep = "\n")
    return(invisible(x))
    #end of the theft
    }
#SCCS  @(#)printcp.s	1.6 01/20/00
# print out the cptable, along with some summary of the tree
printcp <- function(x, digits=getOption("digits")-2) {
    if (!inherits(x, 'rpart')) stop ("Must be an rpart x")
    cat(switch(x$method,anova = "\nRegression tree:\n" ,
			class = "\nClassification tree:\n" ,
			poisson="\nRates regression tree:\n",
			exp = "\nSurvival regression tree:\n")
			)

    if(!is.null(cl <- x$call)) {
	dput(cl)
	cat("\n")
      }
    frame <- x$frame
    leaves <- frame$var == "<leaf>"
    used <- unique(frame$var[!leaves])

    if(!is.null(used)) {
		cat("Variables actually used in tree construction:\n")
		print(sort(as.character(used)), quote=FALSE)
		cat("\n")
	}


    cat("Root node error: ", format(frame$dev[1], digits=digits), '/',
			frame$n[1], ' = ',
		         format(frame$dev[1]/frame$n[1], digits=digits),
 			'\n\n', sep='')


    n <- x$frame$n
    omit <- x$na.action
    if (length(omit))
    cat("n=", n[1], " (", naprint(omit), ")\n\n", sep="")
    else cat("n=", n[1], "\n\n")

    print (x$cptable, digits=digits)
    invisible(x$cptable)
    }

#SCCS @(#)prune.rpart.s	1.8 02/27/98
prune.rpart <- function(tree, cp) {
     ff <- tree$frame
     id <- as.integer(row.names(ff))
     toss <- id[ff$complexity <= cp &  ff$var!='<leaf>'] #not a leaf
     if (length(toss)==0) return(tree)   #all the tree is retained

     newx <- snip.rpart(tree, toss)

     # Now cut down the CP table
     temp <- pmax(tree$cptable[,1], cp)
     keep <- match(unique(temp), temp)
     newx$cptable <- tree$cptable[keep,]
     newx$cptable[max(keep),1] <- cp

     newx
     }
# SCCS @(#)prune.s	1.2 02/12/98
# This should be part of Splus proper -- make prune a method
prune <- function(tree, ...)  UseMethod("prune")
#SCCS  @(#)residuals.rpart.s	1.7 02/11/00


residuals.rpart <- function(object, type)
    {

    if(!inherits(object, "rpart"))
	    stop("Not legitimate rpart object")

    if (object$method=='anova' || object$method=='class')
      { ## code taken directly from residuals.tree
        if(is.null(y <- object$y))
                y <- model.extract(model.frame(object), "response")
        frame <- object$frame

        if(is.null(ylevels <- attr(object, "ylevels"))) {
                tmpr <-y - frame$yval[object$where]    #       y <- unclass(y)
		#Expand out the missing values in the result
                if (!is.null(object$na.action)) 
               	   tmpr <- naresid(object$na.action, tmpr)
                return(tmpr)
	      }

        if(missing(type))
                type <- "usual"
        else if(is.na(match(type, c("usual", "pearson", "deviance"))))
                stop("Don't know about this type of residual")
        if(type == "usual")
                yhat <- frame$yval[object$where]
        else yhat <- frame$yprob[object$where,  ][cbind(seq(y), unclass(y))]
        r <- switch(type,
                usual = as.integer(y != yhat),
                # misclassification
                pearson = (1 - yhat)/yhat,
                # sum((obs-fitted)/fitted)
                deviance = -2 * log(yhat))
        names(r) <- names(y)
       }

    else {
	if(is.null(y <- object$y))
		y <- model.extract(model.frame(object), "response")
	lambdat  <- (object$frame$yval)[object$where] * y[,1]

	events <- y[,2]
	temp <- pmax(events, 1)
	r <- sign(events-lambdat) *
		  sqrt(-2*((events - lambdat) + events*log(lambdat/temp)))
	}

    #Expand out the missing values in the result
    if (!is.null(object$na.action)) 
	r <- naresid(object$na.action, r)

    r
    }
#SCCS %W% %G%
rpart.anova <- function(y, offset, parms, wt) {
    if (!is.null(offset)) y <- y-offset
    list(y=y, parms=0, numresp=1, numy=1,
	 summary= function(yval, dev, wt, ylevel, digits ) {
	     paste("  mean=", formatg(yval, digits),
		   ", MSE=" , formatg(dev/wt, digits),
		   sep='')
	     })
    }
#SCCS @(#)rpart.branch.s	1.2 01/25/97
#
# Compute the "branches" to be drawn for an rpart object
#
rpart.branch <- function(x, y, node, branch) {
    if (missing(branch)) {
	if (exists(parms <-paste(".rpart.parms", dev.cur(), sep="." ),
                   envir=.GlobalEnv)) {
#	    parms <- get(parms, frame=0)
            parms <- get(parms, envir=.GlobalEnv)
            branch <- parms$branch
	    }
	else branch <- 0
        }

    # Draw a series of horseshoes, left son, up, over, down to right son
    #   NA's in the vector cause lines() to "lift the pen"
    is.left <- (node%%2 ==0)        #left hand sons
    node.left <- node[is.left]
    parent <- match(node.left/2, node)
    sibling <- match(node.left+1, node)
    temp <- (x[sibling] - x[is.left])*(1-branch)/2
    xx <- rbind(x[is.left], x[is.left]+ temp,
                x[sibling]- temp, x[sibling], NA)
    yy <- rbind(y[is.left], y[parent], y[parent], y[sibling], NA)
    list(x=xx, y=yy)
    }
#SCCS %W% %G%
rpart.class <- function(y, offset, parms, wt) {
    if (!is.null(offset)) stop("No offset allowed in classification models")
    fy <- as.factor(y)
    y <- as.integer(fy)
    numclass <- max(y[!is.na(y)])
    counts <- tapply(wt, factor(y, levels=1:numclass), sum)
    counts <- ifelse(is.na(counts), 0, counts)   #in case of an empty class
    numresp <- 1+numclass
    if (missing(parms) || is.null(parms))
	parms <- c(counts/sum(counts), rep(1,numclass^2)-diag(numclass),1)
    else if (is.list(parms)) {
	if (is.null(parms$prior)) temp <- c(counts/sum(counts))
	else {
	    temp <- parms$prior
	    if (sum(temp) !=1) stop("Priors must sum to 1")
	    if (any(temp<0)) stop("Priors must be >= 0")
	    if (length(temp) != numclass) stop("Wrong length for priors")
	    }

	if (is.null(parms$loss)) temp2<- 1 - diag(numclass)
	else {
	    temp2 <- parms$loss
	    if (length(temp2) != numclass^2)
			    stop("Wrong length for loss matrix")
	    temp2 <- matrix(temp2, ncol=numclass)
	    if (any(diag(temp2) != 0))
			stop("Loss matrix must have zero on diagonals")
	    if (any(temp2<0))
			stop("Loss matrix cannot have negative elements")
	    if (any(apply(temp2,1,sum)==0))
			stop("Loss matrix has a row of zeros")
	    }

	if (is.null(parms$split)) temp3 <- 1
 	    else {
		temp3 <- pmatch(parms$split, c("gini", "information"))
		if (is.null(temp3)) stop("Invalid splitting rule")
		}
	parms <- c(temp, temp2, temp3)
	}
    else stop("Parameter argument must be a list")

    list(y=y, parms=parms, numresp=numclass+1, counts=counts,
	 ylevels= levels(fy), numy=1,
	 print = function(yval, ylevel, digits) {
	     if (is.null(ylevel))
		     temp <- as.character(yval[,1])
	     else    temp <- ylevel[yval[,1]]

	     nclass <- (ncol(yval) -1)/2
	     if (nclass <6) {
		 yprob <- format(yval[, 1+nclass + 1:nclass],
				 digits=digits, nsmall=digits)
		 temp <- paste(temp, ' (', yprob[,1], sep='')
		 for(i in 2:ncol(yprob))
		     temp  <- paste(temp, yprob[, i], sep=' ')
		 temp <- paste(temp, ")", sep="")
		 }
	     temp
	     },
	 summary= function(yval, dev, wt, ylevel, digits) {
	     nclass <- (ncol(yval)-1) /2
	     group <- yval[, 1]
	     counts <- yval[, 1+ (1:nclass)]
	     yprob  <- yval[, 1+nclass + 1:nclass]
	     if(!is.null(ylevel)) group <- ylevel[group]

	     temp1 <- formatg(counts, digits)
	     temp2 <- formatg(yprob, digits)
	     if (nclass >1) {
		 temp1 <- apply(matrix(temp1, ncol=nclass), 1,
				    paste, collapse=' ')
		 temp2 <- apply(matrix(temp2, ncol=nclass), 1,
				    paste, collapse=' ')
		 }
	     paste("  predicted class=", format(group, justify='left'),
		   "  expected loss=", formatg(dev/wt, digits),"\n",
		   "    class counts: ", temp1,"\n",
		   "   probabilities: ", temp2,
		   sep='')
	     })
    }

#SCCS %W% %G%
rpart.control <-
  function(minsplit=20, minbucket= round(minsplit/3), cp=.01,
	   maxcompete=4, maxsurrogate=5, usesurrogate=2, xval=10, 
	   surrogatestyle =0, maxdepth=30, ... ) {

	if (maxcompete<0) {
	    warning("The value of maxcompete supplied is <0; the value 0 was used instead")
	    maxcompete <-0
	    }
	if (any(xval<0)) {
	    warning("The value of xval supplied is <0; the value 0 was used instead")
	    xval <-0
	    }
	if (maxdepth > 30) stop("Maximum depth is 30")
	if (maxdepth < 1)  stop("Maximum depth must be at least 1")

	if (missing(minsplit) && !missing(minbucket)) minsplit <- minbucket*3

	# Because xval can be of length either 1 or n, and the C code
	#   refers to parameters by number, i.e., "opt[5]" in rpart.c, 
	#   the xval parameter should always be last on the list.
	list(minsplit=minsplit, minbucket=minbucket, cp=cp,
	     maxcompete=maxcompete, maxsurrogate=maxsurrogate,
	     usesurrogate=usesurrogate, 
	     surrogatestyle=surrogatestyle, maxdepth=maxdepth, xval=xval )
	}
#SCCS %W% %G%
# rescaled exponential splitting
rpart.exp <- function(y, offset, parms, wt) {

    if (!inherits(y, "Surv")) 
	   stop("Response must be a survival object - use the Surv() function")

    if (ncol(y) !=2) stop("Can't use (start, stop] data")
    n  <- nrow(y)
    if (any(y[,1]<=0)) stop("Observation time must be >0")
    if (sum(y[,2])==0) stop("No deaths in data set")
    #
    # Rescale time so that the event rate in every interval is = 1
    #
    n  <- nrow(y)
    ord <- order(y[,1])
    temp <- .C("rpartexp", as.integer(n),	
		              as.double(y[ord,]),
		              as.double(wt[ord]),
		              newy = double(n),
	                      double(n))
    newy <- double(n)
    newy[ord] <- temp$newy
    if (length(offset)==n)  newy <- newy * exp(offset)

    if (missing(parms)) parms <- c(shrink=1, method=1)
    else {
	parms <- as.list(parms)
	if (is.null(parms$method)) method <- 1
	else method <- pmatch(parms$method, c("deviance", "sqrt"))
	if (is.na(method)) stop("Invalid error method for Poisson")
	
	if (is.null(parms$shrink)) shrink <- 2-method
	else shrink <- parms$shrink
	if (!is.numeric(shrink) || shrink < 0) 
		stop("Invalid shrinkage value")
	parms <- c(shrink=shrink, method=method)
	}
    list(y=cbind(newy, y[,2]), parms=parms, numresp=2, numy=2,
	 summary= function(yval, dev, wt, ylevel, digits) {
	     paste("  events=", formatg(yval[,2]),
		",  estimated rate=" , formatg(yval[,1], digits),
		" , mean deviance=",formatg(dev/wt, digits),
		sep = "")
	     })
    }
#SCCS  %W% %G%
#
# This differs from tree.matrix in xlevels -- we don't keep NULLS in
#   the list for all of the non-categoricals
#
rpart.matrix <- function(frame)
    {
    if(!inherits(frame, "data.frame"))
	    return(as.matrix(frame))
    frame$"(weights)" <- NULL
    terms <- attr(frame, "terms")
    if(is.null(terms)) predictors <- names(frame)
    else {
	a <- attributes(terms)
	predictors <- as.character(a$variables)[-1] # R change
	removals <- NULL
	if((TT <- a$response) > 0) {
	    removals <- TT
	    frame[[predictors[TT]]] <- NULL
	    }
	if(!is.null(TT <- a$offset)) {
	    removals <- c(removals, TT)
	    frame[[predictors[TT]]] <- NULL
	    }
	if(!is.null(removals)) predictors <- predictors[ - removals]
        labels <- a$term.labels
	if(abs(length(labels)-length(predictors))>0)
	  predictors <- predictors[match(labels,predictors)]
	}

    factors <- sapply(frame, function(x) !is.null(levels(x)))
    characters <- sapply(frame, is.character)
    if(any(factors | characters)) {
	# change characters to factors
	for (preds in predictors[characters])
		frame[[preds]] <- as.factor(frame[[preds]])
        factors <- factors | characters
        column.levels <- lapply(frame[factors], levels)

	# Now make them numeric
	for (preds in predictors[factors])
	     frame[[preds]] <- as.numeric(frame[[preds]])
	x <- as.matrix(frame)
	attr(x, "column.levels") <- column.levels
	}
    else x <- as.matrix(frame[predictors])
    class(x) <- "rpart.matrix"
    x
    }


#SCCS %W% %G%
rpart.poisson <- function(y, offset, parms, wt) {
    if (is.matrix(y)) {
	if (ncol(y)!=2) stop("response must be a 2 column matrix or a vector")
	if (!is.null(offset)) y[,1] <- y[,1] + exp(offset)
	}
    else {
	if (is.null(offset)) y <- cbind(1,y)
	else  y <- cbind( exp(offset), y)
	}
    if (any(y[,1] <=0)) stop("Observation time must be >0")
    if (any(y[,2] <0))  stop("Number of events must be >=0")

    if (missing(parms)) parms <- c(shrink=1, method=1)
    else {
	parms <- as.list(parms)
	if (is.null(parms$method)) method <- 1
	else method <- pmatch(parms$method, c("deviance", "sqrt"))
	if (is.na(method)) stop("Invalid error method for Poisson")
	
	if (is.null(parms$shrink)) shrink <- 2- method
	else shrink <- parms$shrink

	if (!is.numeric(shrink) || shrink <0) 
		stop("Invalid shrinkage value")
	parms <- c(shrink=shrink, method=method)
	}

    list(y=y, parms=parms, numresp=2, numy=2,
	 summary= function(yval, dev, wt, ylevel, digits) {
	     paste("  events=", formatg(yval[,2]),
		",  estimated rate=" , formatg(yval[,1], digits),
		" , mean deviance=",formatg(dev/wt, digits),
		sep = "")
	     })
    }
# SCCS  @(#)rpart.s	1.31 04/23/01
#
#  The recursive partitioning function, for S
#
rpart <- function(formula, data=NULL, weights, subset,
		   na.action=na.rpart, method, model=FALSE, x=FALSE, y=TRUE,
		   parms, control, ...) {


    call <- match.call()
    if (is.data.frame(model)) {
	m <- model
	model <- FALSE
	}
    else {
	m <- match.call(expand=FALSE)
	m$model <- m$method <- m$control<- NULL
	m$x <- m$y <- m$parms <- m$... <- NULL
	m$na.action <- na.action
	m[[1]] <- as.name("model.frame.default")
	m <- eval(m, parent.frame())
	}
    Terms <- attr(m, "terms")
    if(any(attr(Terms, "order") > 1))
	stop("Trees cannot handle interaction terms")

    Y <- model.extract(m, "response")
    wt <- model.extract(m, "weights")
    if(length(wt)==0) wt <- rep(1.0, nrow(m))
    offset <- attr(Terms, "offset")
    X <- rpart.matrix(m)
    nobs <- nrow(X)

    if (missing(method)) {
	if (is.factor(Y) || is.character(Y))      method <- 'class'
        else if (is.Surv(Y))   method <- 'exp'
	else if (is.matrix(Y)) method<- 'poisson'
	else                   method<- 'anova'
	}

    if (is.list(method)) {
	# User written split methods
	mlist <- method
	method <- 'user'

	if (missing(parms)) init <- mlist$init(Y, offset,,wt)
	else                init <- mlist$init(Y, offset, parms, wt)

	method.int <- 4      #the fourth entry in func_table.h
if(F) {
	if (!is.list(mlist) || length(mlist) !=3)
		stop("User written methods must have 3 functions")
	if (is.null(mlist$init) || typeof(mlist$init) != 'closure')
		stop("User written method does not contain an init function")
	if (is.null(mlist$split) || typeof(mlist$split) != 'closure')
		stop("User written method does not contain a split function")
	if (is.null(mlist$eval) || typeof(mlist$eval) != 'closure')
		stop("User written method does not contain an eval function")

	user.eval <- mlist$eval
	user.split <- mlist$split

	numresp <- init$numresp
	numy <-  init$numy
	parms <- init$parms

	#
	# expr2 is an expression that will call the user "evaluation"
	#   function, and check that what comes back is valid
	# expr1 does the same for the user "split" function
	#
	# For speed in the C interface, yback, xback, and wback are
	#  fixed S vectors of a fixed size, and nback tells us how
	#  much of the vector is actually being used on this particular
	#  callback.
	#
	if (numy==1) {
	    expr2 <- quote({
		temp <- user.eval(yback[1:nback], wback[1:nback], parms)
		if (length(temp$label) != numresp)
			stop("User eval function returned invalid label")
		if (length(temp$deviance) !=1)
			stop("User eval function returned invalid deviance")
		as.numeric(as.vector(c(temp$deviance, temp$label)))
		})
	    expr1 <- quote({
		if (nback <0) { #categorical variable
		    n2 <- -1*nback
		    temp  <- user.split(yback[1:n2], wback[1:n2],
					xback[1:n2], parms, FALSE)
		    ncat <- length(unique(xback[1:n2]))
		    if (length(temp$goodness) != ncat-1 ||
			length(temp$direction) != ncat)
			    stop("Invalid return from categorical split fcn")
		    }

		else {
		    temp <- user.split(yback[1:nback], wback[1:nback],
				       xback[1:nback], parms, TRUE)
		    if (length(temp$goodness) != (nback-1))
			stop("User split function returned invalid goodness")
		    if (length(temp$direction) != (nback-1))
			stop("User split function returned invalid direction")
		    }
		as.numeric(as.vector(c(temp$goodness, temp$direction)))
		})
	    }
	else {
	    expr2 <- quote({
		tempy <- matrix(yback[1:(nback*numy)], ncol=numy)
		temp <- user.eval(tempy, wback[1:nback], parms)
		if (length(temp$label) != numresp)
			stop("User eval function returned invalid label")
		if (length(temp$deviance !=1))
			stop("User eval function returned invalid deviance")
		as.numeric(as.vector(c(temp$deviance, temp$label)))
		})
	    expr1 <- quote({
		if (nback <0) { #categorical variable
		    n2 <- -1*nback
		    tempy <- matrix(yback[1:(n2*numy)], ncol=numy)
		    temp  <- user.split(tempy, wback[1:n2], xback[1:n2],
					parms, FALSE)
		    ncat <- length(unique(xback[1:n2]))
		    if (length(temp$goodness) != ncat-1 ||
			length(temp$direction) != ncat)
			    stop("Invalid return from categorical split fcn")
		    }
		else {
		    tempy <- matrix(yback[1:(nback*numy)], ncol=numy)
		    temp <- user.split(tempy, wback[1:nback],xback[1:nback],
				       parms, TRUE)
		    if (length(temp$goodness) != (nback-1))
			stop("User split function returned invalid goodness")
		    if (length(temp$direction) != (nback-1))
			stop("User split function returned invalid direction")
		    }
		as.numeric(as.vector(c(temp$goodness, temp$direction)))
		})
	    }
	#
	# The vectors nback, wback, xback and yback will have their
	#  contents constantly re-inserted by C code.  It's one way to make
	#  things very fast.  It is dangerous to do this, so they
	#  are tossed into a separate frame to isolate them.  Evaluations of
	#  the above expressions occur in that frame.
	#
        rho <- new.env()
        assign("nback", integer(1), envir = rho)
        assign("wback", double(nobs), envir = rho)
        assign("xback", double(nobs), envir = rho)
        assign("yback", double(nobs), envir = rho)
        assign("user.eval", user.eval, envir = rho)
        assign("user.split", user.split, envir = rho)
        assign("numy", numy, envir = rho)
        assign("numresp", numresp, envir = rho)
        assign("parms", parms, envir = rho)
	.Call("init_rpcallback", rho, as.integer(numy),
	                         as.integer(numresp),
	                         expr1, expr2)
    }
        ## assign this to avoid garbage collection
        keep <- rpartcallback(mlist, nobs, init)
    }
    else {
	method.int <- pmatch(method, c("anova", "poisson", "class", "exp"))
	if (is.na(method.int)) stop("Invalid method")
	method <- c("anova", "poisson", "class", "exp")[method.int]
	if (method.int==4) method.int <- 2

	if (missing(parms))
	  init <- (get(paste("rpart", method, sep='.')))(Y,offset, ,wt)
	else
	  init <- (get(paste("rpart", method, sep='.')))(Y,offset, parms, wt)
	}

    Y <- init$y

    xlevels <- attr(X, "column.levels")
    cats <- rep(0,ncol(X))
    if(!is.null(xlevels)) {
	cats[match(names(xlevels), dimnames(X)[[2]])] <-
		  unlist(lapply(xlevels, length))
	}

    controls <- rpart.control(...)
    if (!missing(control)) controls[names(control)] <- control

    xval <- controls$xval
    if (is.null(xval) || (length(xval)==1 && xval==0) || method=='user') {
	xgroups <-0
	xval <- 0
	}
    else if (length(xval)==1) {
	# make random groups
        xgroups <- sample(rep(1:xval, length=nobs), nobs, replace=FALSE)
	}
    else if (length(xval) == nobs) {
	xgroups <- xval
	xval <- length(unique(xgroups))
	}
    else stop("Invalid value for xval")

    #
    # Have s_to_rp consider ordered categories as continuous
    #  A right-hand side variable that is a matrix forms a special case
    # for the code.
    #
    tfun <- function(x) {
	if (is.matrix(x)) rep(is.ordered(x), ncol(x))
	else is.ordered(x)
	}
    isord <- unlist(lapply(m[attr(Terms, 'term.labels')], tfun))
    rpfit <- .C("s_to_rp",
		    n = as.integer(nobs),
		    nvarx = as.integer(ncol(X)),
		    ncat = as.integer(cats* !isord),
		    method= as.integer(method.int),
		    as.double(unlist(controls)),
		    parms = as.double(init$parms),
		    as.integer(xval),
		    as.integer(xgroups),
		    as.double(t(init$y)),
		    as.double(X),
		    as.integer(!is.finite(X)),
		    error = character(1),
		    wt = as.double(wt),
		    as.integer(init$numy),
		    NAOK=TRUE )
    if (rpfit$n == -1)  stop(rpfit$error)

    # rpfit$newX[1:n] contains the final sorted order of the observations
    nodes <- rpfit$n          # total number of nodes
    nsplit<- rpfit$nvarx      # total number of splits, primary and surrogate
    numcp <- rpfit$method     # number of lines in cp table
    ncat  <- rpfit$ncat[1]    #total number of categorical splits
    numresp<- init$numresp    # length of the response vector

    if (nsplit==0) stop("No splits found")
    cpcol <- if (xval>0) 5 else 3
    if (ncat==0) catmat <- 0
    else         catmat <- matrix(integer(1), ncat, max(cats))

    rp    <- .C("s_to_rp2",
		       as.integer(nobs),
		       as.integer(nsplit),
		       as.integer(nodes),
		       as.integer(ncat),
		       as.integer(cats *!isord),
		       as.integer(max(cats)),
		       as.integer(xval),
		       which = integer(nobs),
		       cptable = matrix(double(numcp*cpcol), nrow=cpcol),
		       dsplit =  matrix(double(1),  nsplit,3),
		       isplit =  matrix(integer(1), nsplit,3),
		       csplit =  catmat,
		       dnode  =  matrix(double(1),  nodes, 3+numresp),
		       inode  =  matrix(integer(1), nodes, 6))
    tname <- c("<leaf>", dimnames(X)[[2]])

    if (cpcol==3) temp <- c("CP", "nsplit", "rel error")
    else          temp <- c("CP", "nsplit", "rel error", "xerror", "xstd")
    dimnames(rp$cptable) <- list(temp, 1:numcp)

    splits<- matrix(c(rp$isplit[,2:3], rp$dsplit), ncol=5,
		     dimnames=list(tname[rp$isplit[,1]+1],
			  c("count", "ncat", "improve", "index", "adj")))
    index <- rp$inode[,2]  #points to the first split for each node

    # Now, make ordered categories look like categories again (a printout
    #  choice)
    nadd <- sum(isord[rp$isplit[,1]])
    if (nadd >0) {
	newc <- matrix(integer(1), nadd, max(cats))
	cvar <- rp$isplit[,1]
	indx <- isord[cvar]		     # vector of T/F
	cdir <- splits[indx,2]               # which direction splits went
	ccut <- floor(splits[indx,4])        # cut point
	splits[indx,2] <- cats[cvar[indx]]   #Now, # of categories instead
	splits[indx,4] <- ncat + 1:nadd      # rows to contain the splits

	# Next 4 lines can be done without a loop, but become indecipherable
	for (i in 1:nadd) {
	    newc[i, 1:(cats[(cvar[indx])[i]])] <- -1*as.integer(cdir[i])
	    newc[i, 1:ccut[i]] <- as.integer(cdir[i])
	    }
	if (ncat==0) catmat <- newc
	else         catmat <- rbind(rp$csplit, newc)
	ncat <- ncat + nadd
	}
    else catmat <- rp$csplit

    temp <- ifelse(index==0, 1, index)
    svar <- ifelse(index==0, 0, rp$isplit[temp,1]) #var number for each node
    frame <- data.frame(row.names=rp$inode[,1],
			   var=  factor(svar, 0:ncol(X), tname),
			   n =   rp$inode[,5],
			   wt=   rp$dnode[,3],
			   dev=  rp$dnode[,1],
			   yval= rp$dnode[,4],
			   complexity=rp$dnode[,2],
			   ncompete  = pmax(0, rp$inode[,3]-1),
			   nsurrogate=rp$inode[,4])

    if (method.int ==3 ) {
        numclass <- init$numresp -1
        temp <- rp$dnode[,-(1:4)] %*% diag(init$parms[1:numclass]*
					   sum(init$counts)/init$counts)
        yprob <- matrix(temp /rp$dnode[,3] ,ncol=numclass)
        yval2 <- matrix(rp$dnode[, -(1:3)], ncol=numclass+1)
	frame$yval2 <- cbind(yval2, yprob)
	}
    else if (init$numy >1) frame$yval2 <- rp$dnode[,-(1:3)]

    if (is.null(init$summary))
	    stop("Initialization routine is missing the summary function")
    if (is.null(init$print))
	    functions <- list(summary=init$summary)
    else    functions <- list(summary=init$summary, print=init$print)
    if (!is.null(init$text)) functions <- c(functions, list(text=init$text))
    if (method=='user')	functions <- c(functions, mlist)

    ans <- list(frame = frame,
                where = structure(rp$which, names = row.names(m)),
                call=call, terms=Terms,
    		cptable =  t(rp$cptable),
		splits = splits,
		method = method,
		parms  = init$parms,
		control= controls,
		functions= functions)

    if (ncat>0) ans$csplit <- catmat +2
    if (model) {
	ans$model <- m
	if (missing(y)) y <- FALSE
	}
    if (y) ans$y <- Y
    if (x) {
	ans$x <- X
	ans$wt<- wt
	}
    ans$ordered <- isord
    if (!is.null(xlevels)) attr(ans, 'xlevels') <- xlevels
    if(method=='class') attr(ans, "ylevels") <- init$ylevels
    na.action <- attr(m, "na.action")
    if (length(na.action)) ans$na.action <- na.action
    class(ans) <- c("rpart")
    ans
    }
#  SCCS %W% %G%
#This routine sets up the callback code for user-written split
#  routines in rpart
#
rpartcallback <- function(mlist, nobs, init)
{
	if (length(mlist) < 3)
		stop("User written methods must have 3 functions")
	if (is.null(mlist$init) || typeof(mlist$init) != 'closure')
		stop("User written method does not contain an init function")
	if (is.null(mlist$split) || typeof(mlist$split) != 'closure')
		stop("User written method does not contain a split function")
	if (is.null(mlist$eval) || typeof(mlist$eval) != 'closure')
		stop("User written method does not contain an eval function")

	user.eval <- mlist$eval
	user.split <- mlist$split

	numresp <- init$numresp
	numy <-  init$numy
	parms <- init$parms

	#
	# expr2 is an expression that will call the user "evaluation"
	#   function, and check that what comes back is valid
	# expr1 does the same for the user "split" function
	#
	# For speed in the C interface, yback, xback, and wback are
	#  fixed S vectors of a fixed size, and nback tells us how
	#  much of the vector is actually being used on this particular
	#  callback.
	#
	if (numy==1) {
	    expr2 <- quote({
		temp <- user.eval(yback[1:nback], wback[1:nback], parms)
		if (length(temp$label) != numresp)
			stop("User eval function returned invalid label")
		if (length(temp$deviance) !=1)
			stop("User eval function returned invalid deviance")
		as.numeric(as.vector(c(temp$deviance, temp$label)))
		})
	    expr1 <- quote({
		if (nback <0) { #categorical variable
		    n2 <- -1*nback
		    temp  <- user.split(yback[1:n2], wback[1:n2],
					xback[1:n2], parms, FALSE)
		    ncat <- length(unique(xback[1:n2]))
		    if (length(temp$goodness) != ncat-1 ||
			length(temp$direction) != ncat)
			    stop("Invalid return from categorical split fcn")
		    }

		else {
		    temp <- user.split(yback[1:nback], wback[1:nback],
				       xback[1:nback], parms, TRUE)
		    if (length(temp$goodness) != (nback-1))
			stop("User split function returned invalid goodness")
		    if (length(temp$direction) != (nback-1))
			stop("User split function returned invalid direction")
		    }
		as.numeric(as.vector(c(temp$goodness, temp$direction)))
		})
	    }
	else {
	    expr2 <- quote({
		tempy <- matrix(yback[1:(nback*numy)], ncol=numy)
		temp <- user.eval(tempy, wback[1:nback], parms)
		if (length(temp$label) != numresp)
			stop("User eval function returned invalid label")
		if (length(temp$deviance !=1))
			stop("User eval function returned invalid deviance")
		as.numeric(as.vector(c(temp$deviance, temp$label)))
		})
	    expr1 <- quote({
		if (nback <0) { #categorical variable
		    n2 <- -1*nback
		    tempy <- matrix(yback[1:(n2*numy)], ncol=numy)
		    temp  <- user.split(tempy, wback[1:n2], xback[1:n2],
					parms, FALSE)
		    ncat <- length(unique(xback[1:n2]))
		    if (length(temp$goodness) != ncat-1 ||
			length(temp$direction) != ncat)
			    stop("Invalid return from categorical split fcn")
		    }
		else {
		    tempy <- matrix(yback[1:(nback*numy)], ncol=numy)
		    temp <- user.split(tempy, wback[1:nback],xback[1:nback],
				       parms, TRUE)
		    if (length(temp$goodness) != (nback-1))
			stop("User split function returned invalid goodness")
		    if (length(temp$direction) != (nback-1))
			stop("User split function returned invalid direction")
		    }
		as.numeric(as.vector(c(temp$goodness, temp$direction)))
		})
	    }
	#
	# The vectors nback, wback, xback and yback will have their
	#  contents constantly re-inserted by C code.  It's one way to make
	#  things very fast.  It is dangerous to do this, so they
	#  are tossed into a separate frame to isolate them.  Evaluations of
	#  the above expressions occur in that frame.
	#
        rho <- new.env()
        assign("nback", integer(1), envir = rho)
        assign("wback", double(nobs), envir = rho)
        assign("xback", double(nobs), envir = rho)
        assign("yback", double(nobs), envir = rho)
        assign("user.eval", user.eval, envir = rho)
        assign("user.split", user.split, envir = rho)
        assign("numy", numy, envir = rho)
        assign("numresp", numresp, envir = rho)
        assign("parms", parms, envir = rho)
	.Call("init_rpcallback", rho, as.integer(numy),
	                         as.integer(numresp),
	                         expr1, expr2)
        list(expr1 = expr1, expr2 = expr2, rho = rho)
}
#SCCS @(#)rpartco.s	1.7 02/07/00
# Compute the x-y coordinates for a tree
rpartco <- function(tree, parms =  paste(".rpart.parms", dev.cur(), sep = "."))
    {

    frame <- tree$frame
    method <- tree$method
    node <- as.numeric(row.names(frame))
    depth <- tree.depth(node)
    is.leaf <- (frame$var == '<leaf>')
    if (exists(parms, envir=.GlobalEnv)) {
	parms <- get(parms, envir=.GlobalEnv)
	uniform <- parms$uniform
	nspace <-parms$nspace
	minbranch <- parms$minbranch
	}
    else {
	uniform <- FALSE
	nspace <- -1
	minbranch <- .3
        }

    if(uniform) y <- (1 + max(depth) -depth) / max(depth,4)
    else {                    #make y- (parent y) = change in deviance
	y <- dev <- frame$dev
        temp <- split(seq(node), depth)     #depth 0 nodes, then 1, then ...
        parent <- match(floor(node/2), node)
        sibling <- match(ifelse(node %% 2, node - 1, node + 1), node)

	# assign the depths
        for(i in temp[-1]) {
	    temp2 <- dev[parent[i]] - (dev[i] + dev[sibling[i]])
            y[i] <- y[parent[i]] - temp2
	    }
	#
	# For some problems, classification & loss matrices in particular
	#   the gain from a split may be 0.  This is ugly on the plot.
	# Hence the "fudge" factor of  .3* the average step
	#
	fudge <-  minbranch * diff(range(y)) / max(depth)
        for(i in temp[-1]) {
	    temp2 <- dev[parent[i]] - (dev[i] + dev[sibling[i]])
	    haskids <- !(is.leaf[i] & is.leaf[sibling[i]])
	    y[i] <- y[parent[i]] - ifelse(temp2<=fudge & haskids, fudge, temp2)
	    }
	y <- y / (max(y))
        }

    # Now compute the x coordinates, by spacing out the leaves and then
    #   filling in
    x   <-  double(length(node))         #allocate, then fill it in below
    x[is.leaf] <- seq(sum(is.leaf))      # leaves at 1, 2, 3, ....
    left.child <- match(node * 2, node)
    right.child <- match(node * 2 + 1, node)

    # temp is a list of non-is.leaf, by depth
    temp <- split(seq(node)[!is.leaf], depth[!is.leaf])
    for(i in rev(temp))
            x[i] <- 0.5 * (x[left.child[i]] + x[right.child[i]])

    if (nspace < 0) return(list(x=x, y=y))

    #
    # Now we get fancy, and try to do overlapping
    #
    #  The basic algorithm is, at each node:
    #      1: get the left & right edges, by depth, for the left and
    #           right sons, of the x-coordinate spacing.
    #      2: find the minimal free spacing.  If this is >0, slide the
    #           right hand son over to the left
    #      3: report the left & right extents of the new tree up to the
    #           parent
    #   A way to visualize steps 1 and 2 is to imagine, for a given node,
    #      that the left son, with all its descendants, is drawn on a
    #      slab of wood.  The left & right edges, per level, give the
    #      width of this board.  (The board is not a rectangle, it has
    #      'stair step' edges). Do the same for the right son.  Now
    #      insert some spacers, one per level, and slide right hand
    #      board over until they touch.  Glue the boards and spacer
    #      together at that point.
    #
    #  If a node has children, its 'space' is considered to extend left
    #    and right by the amount "nspace", which accounts for space
    #    used by the arcs from this node to its children.  For
    #    horseshoe connections nspace usually is 1.
    #
    #  To make it global for a recursive function, the x coordinate list
    #    is written into frame 0.
    #
    compress <- function(me, depth) {
        lson <- me +1
	x <- x
	if (is.leaf[lson]) left <- list(left=x[lson], right=x[lson],
						depth=depth+1, sons=lson)
        else               left <- compress(me+1, depth+1)

        rson <- me + 1 + length(left$sons)        #index of right son
	if (is.leaf[rson]) right<- list(left=x[rson], right=x[rson],
						depth=depth+1, sons=rson)
	else               right<- compress(rson, depth+1)

	maxd <- max(left$depth, right$depth) - depth
        mind <- min(left$depth, right$depth) - depth

	# Find the smallest distance between the two subtrees
	#   But only over depths that they have in common
	# 1 is a minimum distance allowed
	slide <- min(right$left[1:mind] - left$right[1:mind]) -1
	if (slide >0) { # slide the right hand node to the left
	    x[right$sons] <- x[right$sons] - slide;
	    x[me] <- (x[right$sons[1]] + x[left$sons[1]])/2
#	    assign("x", x)
            x <<- x
	    }
	else slide <- 0

	# report back
        if (left$depth > right$depth) {
	    templ <- left$left
            tempr <- left$right
            tempr[1:mind] <- pmax(tempr[1:mind], right$right -slide)
	    }
        else {
	    templ <- right$left  - slide
	    tempr <- right$right - slide
	    templ[1:mind] <- pmin(templ[1:mind], left$left)
	    }

	list(left = c(x[me]- nspace*(x[me] -x[lson]), templ),
	     right= c(x[me]- nspace*(x[me] -x[rson]), tempr),
	     depth= maxd+ depth, sons=c(me, left$sons, right$sons))
	}
#    assign('compress', compress)
#    assign('x', x)
#    assign('is.leaf', is.leaf)
#    assign('nspace', nspace)

#    temp <-
    compress(1, 1)
#    x <- get('x')
#    remove(c('compress', 'x', 'is.leaf', 'nspace'))
    list(x = x, y = y)
}

## This function plots the approximate r-square for the different
## splits (assumes using anova method).

## SCCS @(#)rsq.rpart.s	1.6 08/28/97

rsq.rpart <- function(x) {

  if(!inherits(x,'rpart')) stop("Not legitimate rpart")

  p.rpart <- printcp(x)
  xstd <- p.rpart[,5]
  xerror <- p.rpart[,4]
  rel.error <- p.rpart[,3]
  nsplit <- p.rpart[,2]
  method <- x$method

  if(!method=='anova') cat("May not be applicable for this method\n")

  plot(nsplit, 1-rel.error, xlab='Number of Splits', ylab='R-square',
       ylim=c(0,1), type='o')
  par(new=TRUE)
  plot(nsplit, 1-xerror, type='o', ylim=c(0,1),lty=2, xlab=' ', ylab=' ')
  legend(0,1, c('Apparent','X Relative'), lty=1:2)


  ylim <- c(min(xerror-xstd) -.1, max(xerror + xstd) + .1)
  plot(nsplit, xerror, xlab='Number of Splits', ylab='X Relative Error',
       ylim=ylim, type='o')
  segments(nsplit, xerror - xstd, nsplit, xerror + xstd)
  invisible()

  }

# SCCS %W% %G%
#
#  Interactively snip off part of a tree
#

snip.rpart.mouse <- function(tree,
		      parms=paste(".rpart.parms", dev.cur(), sep = ".")) {
    xy <- rpartco(tree)
    toss <- NULL
    ff <- tree$frame
    if (exists(parms, envir=.GlobalEnv)) {
        parms <- get(parms, envir=.GlobalEnv)
	branch <- parms$branch
	}
    else branch <- 1

    node <- as.numeric(row.names(tree$frame))
    draw <- rpart.branch(xy$x,xy$y, node, branch)

    lastchoice <- 0
    while (length(choose <- identify(xy, n=1, plot=FALSE)) >0 ) {
	if (ff$var[choose] == '<leaf>') {
		cat("Terminal node -- try again\n")
		next
		}

	if (choose != lastchoice) {
	    # print out some info on the click
	    cat("node number:", node[choose], " n=", ff$n[choose], "\n")
	    cat("    response=", format(ff$yval[choose]))
	    if (is.null(ff$yval2)) cat ("\n")
	    else if (is.matrix(ff$yval2))
		  cat(" (", format(ff$yval2[choose,]), ")\n")
	    else  cat(" (", format(ff$yval2[choose]), ")\n")
	    cat("    Error (dev) = ", format(ff$dev[choose]), "\n")
	    lastchoice <- choose
	    }
	else {
	    # second click-- erase all of the descendants
	    #   (stolen from snip.tree)
	    id  <- node[choose]
	    id2 <- node
	    while (any(id2>1)) {
		id2 <- floor(id2/2)
		temp  <- (match(id2, id, nomatch=0) >0)
  	        id <- c(id, node[temp])
		id2[temp] <- 0
		}
	    temp <- match(id, node[ff$var != '<leaf>'], nomatch=0)
	    lines(c(draw$x[,temp]), c(draw$y[,temp]), col=0)
	    toss <- c(toss, node[choose])
	    }
	}
    toss
    }
#SCCS  @(#)snip.rpart.s	1.8 01/20/00
#
#  This routine "throws away" branches
#
snip.rpart <- function(x, toss) {
    if (!inherits(x, 'rpart')) stop("Not an rpart object")

    if (missing(toss) || length(toss)==0) {
        toss <- snip.rpart.mouse(x)
	if (length(toss)==0) return(x)
	}

    where <- x$where
    ff   <- x$frame
    id    <- as.numeric(row.names(ff))
    index <- ff$index
    ff.n  <- length(id)

    toss <- unique(toss)
    toss.idx <- match(toss, id, nomatch=0) #the rows of the named nodes
    if (any(toss.idx ==0)) {
	warning(paste("Nodes", toss[toss.idx==0], "are not in this tree"))
	toss <- toss[toss.indx>0]
        toss.idx <- toss.idx[toss.idx>0]
        }

    if (any(toss==1))  {
	# a special case that causes grief later
	warning("Can't prune away the root node and still have a tree!")
        return(NULL)
	}

    # Now add all of the descendants of the selected nodes
    #   We do this be finding all node's parents.
    #        (Division by 2 gives the parent of any node.)
    #   At each step we make id2 <- parent(id2), and augment 'toss' with
    #     found children.  The loop should take <  log_2(maxdepth)/2 steps
    id2 <- id
    while (any(id2>1)) {
	id2 <- floor(id2/2)
	xx <- (match(id2, toss, nomatch=0) >0)
	toss <- c(toss, id[xx])
        id2[xx] <- 0
	}

    # Now "toss" contains all of the nodes that should not be splits
    temp <- match(floor(toss/2) , toss, nomatch=0)  #which are leaves?
    newleaf <- match(toss[temp==0], id)             # row numbers, leaves
    keepit <- (1:ff.n)[is.na(match(id,toss))]  # row numbers to be let be

    # Compute the parent row for each row in the splits structure
    #  Then "thin out" the splits and csplit components
    n.split <- rep((1:ff.n), ff$ncompete + ff$nsurrogate+ 1*(ff$var!='<leaf>'))
    split <- x$splits[match(n.split, keepit, nomatch=0) >0, ,drop=FALSE]
    temp <- split[,2] >1      #which rows point to categoricals?
    if (any(temp)) {
        x$csplit <- x$csplit[split[temp,4], , drop=FALSE]
	split[temp,4] <- 1
        if(is.matrix(x$csplit)) split[temp,4] <- 1:nrow(x$csplit)
	}
    else x$csplit <- NULL
    x$splits <- split

    # Thin out unneeded rows in the frame component
    ff$ncompete[newleaf] <- ff$nsurrogate[newleaf] <- 0
    ff$var[newleaf]     <- "<leaf>"
##    ff$splits[newleaf,] <- ""
    x$frame <- ff[sort(c(keepit, newleaf)),]

    # Now do the 'parents' loop one more time, to fix up the "where"
    #   vector
    # This pass requires log_2(depth) iterations
    #
    id2 <- id[x$where]         #the list of old leaf nodes
    id3 <- id[sort(c(keepit, newleaf))]
    temp <- match(id2, id3, nomatch=0)
    while (any(temp==0)) {
	id2[temp==0] <- floor(id2[temp==0]/2)
	temp <- match(id2, id3, nomatch=0)
	}
    x$where <- match(id2, id3)

    x
    }
#SCCS  %W% %G%
summary.rpart <- function(object, cp=0, digits=getOption("digits"), file,  ...)
{
    if(!inherits(object, "rpart")) stop("Not legitimate rpart object")

     if (!missing(file)) {
	  sink(file)
	  on.exit(sink())
	  }

    if(!is.null(object$call)) {
        cat("Call:\n")
        dput(object$call)
        }

     omit <- object$na.action
     n <- object$frame$n
     if (length(omit))
      cat("  n=", n[1], " (", naprint(omit), ")\n\n", sep="")
     else cat("  n=", n[1], "\n\n")

     print(object$cptable, digits=digits)
     ff <- object$frame
     ylevel <- attr(object,'ylevels')
     id <- as.integer(row.names(ff))
     parent.id <- ifelse(id==1,1, floor(id/2))
     parent.cp <- ff$complexity[match(parent.id, id)]
     rows <- (1:length(id))[parent.cp > cp]
     if (length(rows)>0) rows <- rows[order(id[rows])]
     else rows <- 1
     is.leaf <- (ff$var=='<leaf>')
     index <- cumsum(c(1, ff$ncompete + ff$nsurrogate + 1*(!is.leaf)))

     sname <- dimnames(object$splits)[[1]]
     cuts <- vector(mode='character', length=nrow(object$splits))
     temp <- object$splits[ ,2]
     for (i in 1:length(cuts)) {
	 if (temp[i] == -1)
	     cuts[i] <-paste("<", format(signif(object$splits[i,4], digits=digits)))
	 else if (temp[i] ==1)
	     cuts[i] <-paste("<", format(signif(object$splits[i,4], digits=digits)))
	 else cuts[i]<- paste("splits as ",
	     paste(c("L", "-", "R")[object$csplit[object$splits[i,4], 1:temp[i]]],
                   collapse='', sep=''), collapse='')
	 }
    # S-PLUS 4.0 can't handle null vectors here
    if(any(temp<2)) cuts[temp<2 ] <- format(cuts[temp<2],justify="left")
    cuts <- paste(cuts, ifelse(temp >=2, ",",
			 ifelse(temp==1, " to the right,", " to the left, ")),
			 sep = '')

    if (is.null(ff$yval2))
	    tprint <- object$functions$summary(ff$yval[rows], ff$dev[rows],
					  ff$wt[rows], ylevel, digits)
    else
	    tprint <- object$functions$summary(ff$yval2[rows,], ff$dev[rows],
					  ff$wt[rows], ylevel, digits)

    for (ii in 1:length(rows)) {
	i <- rows[ii]
	nn <- ff$n[i]
	twt <- ff$wt[i]
	cat("\nNode number ", id[i], ": ", nn, " observations", sep='')
	if (ff$complexity[i] < cp || is.leaf[i]) cat("\n")
	else cat(",    complexity param=",
		       format(signif(ff$complexity[i], digits)), "\n", sep="")

	cat(tprint[ii], "\n")
	if (ff$complexity[i] > cp && !is.leaf[i] ){
	    sons <- 2*id[i] + c(0,1)
	    sons.n <- ff$n[match(sons, id)]
	    cat("  left son=", sons[1], " (", sons.n[1], " obs)",
		" right son=", sons[2], " (", sons.n[2], " obs)", sep='')
	    j <- nn - (sons.n[1] + sons.n[2])
	    if (j>1) cat(", ", j, " observations remain\n", sep='')
	    else if (j==1) cat(", 1 observation remains\n")
	    else     cat("\n")
	    cat("  Primary splits:\n")
	    j <- seq(index[i], length=1+ff$ncompete[i])
	    if (all(nchar(cuts[j]) < 25))
		  temp <- format(cuts[j], justify="left")
	    else  temp <- cuts[j]
	    cat(paste("      ", format(sname[j], justify="left"), " ", temp,
		      " improve=", format(signif(object$splits[j,3], digits)),
		      ", (", nn - object$splits[j,1], " missing)", sep=''),
		      sep="\n")
	    if (ff$nsurrogate[i] >0) {
		cat("  Surrogate splits:\n")
		j <- seq(1 +index[i] + ff$ncompete[i], length=ff$nsurrogate[i])
		agree <- object$splits[j,3]
		adj   <- object$splits[j,5]
		if (all(nchar(cuts[j]) < 25))
		      temp <- format(cuts[j], justify="left")
		else  temp <- cuts[j]
		cat(paste("      ", format(sname[j], justify="left"), " ",
		      temp,
		      " agree=", format(round(agree, 3)),
                      ", adj=" , format(round(adj, 3)),
		      ", (", object$splits[j,1], " split)", sep=''),
		      sep="\n")
		}
	    }
	}
    cat("\n")
    invisible(object)
    }
# SCCS %W% %G%
# This is a modification of text.tree.
# Fancy option has been added in (to mimic post.tree)
#

text.rpart <-  function(x, splits = TRUE, label = "yval", FUN = text, all=FALSE,
		        pretty = NULL, digits = getOption("digits") - 3,
                        use.n=FALSE, fancy=FALSE, fwidth=.8, fheight =.8, ...)
{
	if(!inherits(x, "rpart")) stop("Not legitimate rpart")
	frame <- x$frame
	col <- names(frame)
	method <- x$method
	ylevels <- attr(x,'ylevels')
	if(!is.null(ylevels <- attr(x, "ylevels")))
		col <- c(col, ylevels)
	if(is.na(match(label, col)))
		stop("Label must be a column label of the frame component of the tree"
			)
	cxy <- par("cxy") #character width and height
	if(!is.null(srt <- list(...)$srt) && srt == 90)
		cxy <- rev(cxy)
	xy <- rpartco(x)



        node <- as.numeric(row.names(x$frame))
        is.left <- (node%%2 ==0)        #left hand sons
        node.left <- node[is.left]
	parent <- match(node.left/2, node)

#Put left splits at the parent node

	if(splits) {
		left.child <- match(2 * node, node)
		right.child <- match(node * 2 + 1, node)
		rows <- labels(x, pretty = pretty)

		if(fancy) {
		  ## put split labels on branches instead of nodes

		  xytmp <- rpart.branch(x=xy$x,y=xy$y,node=node)
		  leftptx <- (xytmp$x[2,]+xytmp$x[1,])/2
		  leftpty <- (xytmp$y[2,]+xytmp$y[1,])/2
		  rightptx <- (xytmp$x[3,]+xytmp$x[4,])/2
		  rightpty <- (xytmp$y[3,]+xytmp$y[4,])/2

		  FUN(leftptx,leftpty+.52*cxy[2],
		      rows[left.child[!is.na(left.child)]],...)
		  FUN(rightptx,rightpty-.52*cxy[2],
		      rows[right.child[!is.na(right.child)]],...)
		}

		else FUN(xy$x, xy$y + 0.5 * cxy[2], rows[left.child], ...)
	}
	leaves <- if(all) rep(T, nrow(frame)) else frame$var == "<leaf>"
        if(method=='class') {
	    nclass <- (ncol(frame$yval2)-1) /2
	    ycount <- frame$yval2[, 1 + 1:nclass]
	    yprob  <- frame$yval2[, 1+nclass + 1:nclass]

            if (label=='yval') stat <- ylevels[frame$yval[leaves]]
	    else  if(!is.na(lev <- match(label, ylevels)))
		stat <- format(signif(yprob[leaves, lev],
				      digits = digits))
            else if(label=='yprob'){
                sub <- matrix(c(1:length(frame$yval),frame$yval),
		              nrow=length(frame$yval))
	        stat <- format(signif(yprob[sub][leaves],
				      digits = digits))
				      }
            else stat <- format(signif(frame[leaves,label],digits=digits))
            if(use.n)
		  stat <- paste(stat,'\n','(',
			   apply(ycount[leaves,], 1, paste, collapse='/'),
      			     ')', sep='')
	      }
	else if(method=='anova') {
	    stat <- format(signif(frame[leaves, label], digits =digits))
            if(use.n) stat <-
	       paste(stat,'\n',' n=',frame$n[leaves], sep='')
	   }
	else if(method=='poisson'|method=='exp')
	  {
	   stat <- format(signif(frame[leaves, label], digits = digits))
	   if(use.n)
             { stat <-
	     paste(stat,'\n', frame$yval2[leaves,2],'/',frame$n[leaves], sep="")
	       }

	 }

        oval <- function(middlex,middley,a,b) {

             theta <- seq(0,2*pi,pi/30)
	     newx <- middlex + a*cos(theta)
	     newy <- middley + b*sin(theta)

	     polygon(newx,newy,border=TRUE,col=0)
#	     polygon(newx,newy,border=T)
	   }

        rectangle <- function(middlex, middley,a,b) {

	  newx <- middlex + c(a,a,-a,-a)
	  newy <- middley + c(b,-b,-b,b)

	  polygon(newx,newy,border=TRUE,col=0)
#	  polygon(newx,newy,border=T)
          }

        if(fancy) {

	        ## find maximum length of stat
	        maxlen <- max(string.bounding.box(stat)$columns) + 1
	        maxht <- max(string.bounding.box(stat)$rows) +1

		if(fwidth<1)  a.length <- fwidth*cxy[1]*maxlen
		else a.length <- fwidth*cxy[1]

		if(fheight<1) b.length <- fheight*cxy[2]*maxht
		else b.length <- fheight*cxy[2]

	        ### create ovals and rectangles here
		## sqrt(2) creates the smallest oval that fits around the
		## best fitting rectangle
		for(i in parent) oval(xy$x[i],xy$y[i],
  			          a=sqrt(2)*a.length/2, b=sqrt(2)*b.length/2)
		child <- match(node[frame$var=="<leaf>"],node)
		for(i in child) rectangle(xy$x[i],xy$y[i],
				  a=a.length/2,b=b.length/2)
	      }

#if FUN=text then adj=1 puts the split label to the left of the
#    split rather than centered
#Allow labels at all or just leaf nodes

	## stick values on nodes
	if(fancy) FUN(xy$x[leaves], xy$y[leaves] + .5 * cxy[2], stat, ...)
	else FUN(xy$x[leaves], xy$y[leaves] - 0.5 * cxy[2], stat, adj=.5, ...)

	invisible()
}
# SCCS
#
#  Get a set of cross-validated predictions
xpred.rpart <- function(fit, xval=10, cp) {
    if (!inherits(fit, 'rpart')) stop("Invalid fit object")

    method <- fit$method
    method.int <- pmatch(method, c("anova", "poisson", "class", "user", "exp"))
    if (method.int==5) method.int <- 2
    Terms <- fit$terms

    Y <- fit$y
    X <- fit$x
    wt<- fit$wt
    if (is.null(Y) || is.null(X)) {
	m <- fit$model
	if (is.null(m)) {
	    m <-fit$call[match(c("", 'formula', 'data', 'weights', 'subset',
					 'na.action'),
				names(fit$call), nomatch=0)]
	    if (is.null(m$na.action)) m$na.action<- na.rpart
	    m[[1]] <- as.name("model.frame.default")
	    m <- eval(m, parent.frame())
	    }
	if (is.null(X)) X <- rpart.matrix(m)
	if (is.null(wt)) wt <- model.extract(m, "weights")
	if (is.null(Y)) {
	    yflag <- TRUE
	    Y <- model.extract(m, "response")
            offset <- attr(Terms, "offset")
	    if (method != user) {
		init <- (get(paste("rpart", method, sep='.')))(Y,offset, NULL)
		Y <- as.matrix(init$y)
		numy <- ncol(Y)
		}
	    }
	else {
	    yflag <- FALSE
	    Y <- as.matrix(Y)
	    numy <- ncol(Y)
	    offset <- 0
	    }
	}

    nobs <- nrow(X)
    if (length(wt)==0) wt <- rep(1.0, nobs)

    cats <- rep(0, ncol(X))
    xlevels <- attr(fit, "xlevels")
    if (!is.null(xlevels)){
        cats[match(names(xlevels), dimnames(X)[[2]])] <-
		 unlist(lapply(xlevels, length))
        }

    controls <- fit$control
    if (missing(cp)) {
	cp<- fit$cptable[,1]
	cp <- sqrt(cp * c(10, cp[-length(cp)]))
	cp[1] <- (1+fit$cptable[1,1])/2
	}
    ncp <- length(cp)

    if (length(xval)==1) {
	# make random groups
	xgroups <- sample(rep(1:xval, length=nobs), nobs, replace=FALSE)
	}
    else if (length(xval) == nrow(Y)) {
	xgroups <- xval
	xval <- length(unique(xgroups))
	}
    else stop("Invalid value for xval")

    parms <- as.double(fit$parms)
    if (method=='user') {
	mlist <- fit$functions

        if (!is.list(mlist) || length(mlist) <3)
		stop("User written methods must have 3 functions")
	if (is.null(mlist$init) || typeof(mlist$init) != 'closure')
		stop("User written method does not contain an init function")

	# I need to call init to find out numresp, even though I
	#   only intend to use the first element of the response vector
	#  (For the "built-in" routines numresp is hardcoded into C).
	# If yflag==F, the "transformed" Y is already in hand, so we don't
	#   need to replace it with init$y
	if (length(parms)==0) init <- mlist$init(Y, offset,,wt)
	else                  init <- mlist$init(Y, offset, parms, wt)

        keep <- rpartcallback(mlist, nobs, init)

        if(F) {
	numresp <- init$numresp
	numy <-  init$numy
	if (yflag)  Y <- init$y

	if (is.null(mlist$split) || class(mlist$split) != 'function')
		stop("User written method does not contain a split function")
	if (is.null(mlist$eval) || class(mlist$eval) != 'function')
		stop("User written method does not contain an eval function")

	user.eval <- mlist$eval
	user.split <- mlist$split

	numresp <- init$numresp
	numy <-  init$numy
	parms <- init$parms

	#
	# expr2 is an expression that will call the user "evaluation"
	#   function, and check that what comes back is valid
	# expr1 does the same for the user "split" function
	#
	# For speed in the C interface, yback, xback, and wback are
	#  fixed S vectors of a fixed size, and nback tells us how
	#  much of the vector is actually being used on this particular
	#  callback.
	#
	if (numy==1) {
	    expr2 <- Quote({
		temp <- user.eval(yback[1:nback], wback[1:nback], parms)
		if (length(temp$label) != numresp)
			stop("User eval function returned invalid label")
		if (length(temp$deviance) !=1)
			stop("User eval function returned invalid deviance")
		as.numeric(as.vector(c(temp$deviance, temp$label)))
		})
	    expr1 <- Quote({
		if (nback <0) { #categorical variable
		    n2 <- -1*nback
		    temp  <- user.split(yback[1:n2], wback[1:n2],
					xback[1:n2], parms, F)
		    ncat <- length(unique(xback[1:n2]))
		    if (length(temp$goodness) != ncat-1 ||
			length(temp$direction) != ncat)
			    stop("Invalid return from categorical split fcn")
		    }

		else {
		    temp <- user.split(yback[1:nback], wback[1:nback],
				       xback[1:nback], parms, T)
		    if (length(temp$goodness) != (nback-1))
			stop("User split function returned invalid goodness")
		    if (length(temp$direction) != (nback-1))
			stop("User split function returned invalid direction")
		    }
		as.numeric(as.vector(c(temp$goodness, temp$direction)))
		})
	    }
	else {
	    expr2 <- Quote({
		tempy <- matrix(yback[1:(nback*numy)], ncol=numy)
		temp <- user.eval(tempy, wback[1:nback], parms)
		if (length(temp$label) != numresp)
			stop("User eval function returned invalid label")
		if (length(temp$deviance !=1))
			stop("User eval function returned invalid deviance")
		as.numeric(as.vector(c(temp$deviance, temp$label)))
		})
	    expr1 <- Quote({
		if (nback <0) { #categorical variable
		    n2 <- -1*nback
		    tempy <- matrix(yback[1:(n2*numy)], ncol=numy)
		    temp  <- user.split(tempy, wback[1:n2], xback[1:n2],
					parms, F)
		    ncat <- length(unique(xback[1:n2]))
		    if (length(temp$goodness) != ncat-1 ||
			length(temp$direction) != ncat)
			    stop("Invalid return from categorical split fcn")
		    }
		else {
		    tempy <- matrix(yback[1:(nback*numy)], ncol=numy)
		    temp <- user.split(tempy, wback[1:nback],xback[1:nback],
				       parms, T)
		    if (length(temp$goodness) != (nback-1))
			stop("User split function returned invalid goodness")
		    if (length(temp$direction) != (nback-1))
			stop("User split function returned invalid direction")
		    }
		as.numeric(as.vector(c(temp$goodness, temp$direction)))
		})
	    }
	#
	# The vectors nback, wback, xback and yback will have their
	#  contents constantly re-inserted by C code.  It's one way to make
	#  things very fast.  It is dangerous to do this, so they
	#  are tossed into a separate frame to isolate them.  Evaluations of
	#  the above expressions occur in that frame.
	#
	eframe <- new.frame(list(nback = integer(1),
				 wback = double(nobs),
				 xback = double(nobs),
				 yback = double(nobs),
				 user.eval =  user.eval,
				 user.split = user.split,
				 numy = numy,
				 numresp = numresp,
				 parms = parms), protect=T)
	.Call("init_rpcallback", eframe, as.integer(numy),
	                         as.integer(numresp),
	                         expr1, expr2)
	}
    }

    rpfit <- .C("s_xpred",
		    n = as.integer(nobs),
		    nvarx = as.integer(ncol(X)),
		    ncat = as.integer(cats * !fit$ordered),
		    method= as.integer(method.int),
		    as.double(unlist(controls)),
		    parms = as.double(parms),
		    as.integer(xval),
		    as.integer(xgroups),
		    as.double(t(Y)),
		    as.double(X),
		    as.integer(is.na(X)),
		    pred = double(ncp* nobs),
		    as.integer(ncp),
		    as.double(cp * fit$frame[1,"dev"]),
		    error = character(1),
		    wt = as.double(wt),
		    as.integer(numy),
		    NAOK=TRUE )
    if (rpfit$n == -1)  stop(rpfit$error)

    matrix(rpfit$pred, ncol=ncp, byrow=TRUE,
		dimnames=list(dimnames(X)[[1]], format(cp)) )
    }
.First.lib <- function(lib, pkg) library.dynam("rpart", pkg, lib)

is.Surv <- function(x) inherits(x, "Surv")

tree.depth <- function (nodes)
{
    depth <- floor(log(nodes, base = 2) + 1e-7)
    as.vector(depth - min(depth))
}

string.bounding.box <- function(s)
{
    s2 <- strsplit(s, "\n")
    rows <- sapply(s2, length)
    columns <- sapply(s2, function(x) max(nchar(x)))
    list(columns=columns, rows=rows)
}
