dtrace <-
function(x,  bw, sd = "nrd0", kernel = "gaussian", ...)
UseMethod("dtrace")

dtrace.numeric <-
function(x, bw, sd = "nrd0", kernel = "gaussian", ...)
{
    if (!missing(bw)) sd <-  bw / sqrt(12)
    ans <- density(x, bw = sd, kernel = kernel[1], ...)
    structure(list(density=structure(list(ans[c("x", "y")]),
                     names = deparse(substitute(x))),
                   xlim = range(ans$x),
                   ylim = range(ans$y),
                   nobs = ans$n,
                   kernel = kernel[1],
                   sd = ans$bw,
                   bw = sqrt(12) * ans$bw),
              class = "dtrace")
}

dtrace.list <-
function(x, bw, sd = "nrd0", kernel = "gaussian", ...)
{
    nx = length(x)
    if (!missing(bw))
      sd =  bw / sqrt(12)
    sd <- rep(sd, length = nx)
    kernel = rep(kernel, length = nx)
    xnames <- names(x)
    if (is.null(xnames))
      xnames <- paste("Group", 1:nx)
    dens <- vector("list", nx)
    sds <- bws <- kns <- nobs <- NULL
    xlim <- NULL
    xlim <- NULL
    ylim <- 0
    for(i in 1:nx) {
        d <- dtrace(x[[i]], sd = sd[i], kernel = kernel[i], ...)
        dens[i] <- d$density[1]
        xlim <- range(xlim, d$xlim)
        ylim <- range(ylim, d$ylim)
        nobs <- c(nobs, d$nobs)
        kns <- c(kns, d$kernel)
        bws <- c(bws, d$bw)
        sds <- c(sds, d$sd)
    }
    names(dens) = xnames
    structure(list(density = dens, xlim = xlim, ylim = ylim,
                   bs = bws, sd = sds, kernel = kns),
              class="dtrace")
}

dtrace.formula <-
function(formula, bw, sd = "nrd0", kernel = "gaussian",
         data = NULL, subset, ...)
{
    if (missing(formula) || (length(formula) != 3)) 
      stop("formula missing or incorrect")
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval(m$data, parent.frame()))) 
      m$data <- as.data.frame(data)
    m$... <- NULL
    m[[1]] <- as.name("model.frame")
    mf <- eval(m, parent.frame())
    response <- attr(attr(mf, "terms"), "response")
    dtrace(split(mf[[response]], mf[-response]),
                 bw, sd, kernel, ...)
}

plot.dtrace <-
function(x, col = NULL, lty=1:length(x$density), lwd = NULL,
         main = "", xlab = "", ylab = "Density", log = "",
         add = FALSE, draw.legend = length(x$density) > 1 && !add,
         ...)
{
    if (!is.null(col)) col <- rep(col, length = length(x))
    if (!is.null(lty)) lty <- rep(lty, length = length(x))
    if (!is.null(lwd)) lwd <- rep(lwd, length = length(x))
    if (!add) {
        xlim = x$xlim
        ylim = x$ylim
        plot(NA, NA, xlim = xlim, ylim = ylim,
             xlab = "", ylab = "Density", log = "")
        abline(h=0, col="gray")
        title(main = main, xlab = xlab, ylab = ylab)
    }
    nx = length(x$density)
    for(i in 1:nx)
      lines(x$density[[i]], col = col[i], lty = lty[i], lwd = lwd[i])
    if (draw.legend) {
        legend.args <- list(...)
        args <- list(x = xlim[2],
                     y = ylim[2],
                     xjust = 1,
                     yjust = 1,
                     legend = names(x$density),
                     bty = "n")
        args$col <- col
        args$lty <- lty
        args$lwd <- lwd
        if (is.list(legend.args)) {
            fmls <- names(formals(legend))
            fmls <- fmls[match(names(legend.args), fmls)]
            if (length(fmls) > 0)
              args[fmls] <- legend.args[fmls]
        }
        do.call("legend", args)
    }
    invisible(x)
}

summary.dtrace <-
function(object, ...)
{
    structure(list(
                   xlim = object$xlim,
                   ylim = object$ylim,
                   names = names(object$density),
                   nobs = object$nobs,
                   kernel = object$kernel,
                   bw =object$bw,
                   sd = object$sd),
              class = "summary.dtrace")
}

print.summary.dtrace <-
function(x, digits = max(4, getOption("digits") - 4), ...)
{
    cat("\nDensity Trace:\n\n")
    ans <- cbind(format(x$nobs), x$kernel,
                 format(x$bw, digits = digits),
                 format(x$sd, digits = digits))
    dimnames(ans) <- list(x$names, c("nobs", "kernel", "bw", "sd"))
    print(ans, right = TRUE, print.gap = 2, quote = FALSE, ...)
    cat("\n")
    invisible(x)
}

print.dtrace <-
function(x, digits=max(4, getOption("digits") - 4), ...)
{
    print(summary(x), digits=digits, ...)
    invisible(x)
}

"[.dtrace" <-
function(x, i)
{
    xlim <- ylim <- NULL
    for(d in x$density) {
        xlim <- range(xlim, range(d$x))
        ylim <- range(ylim, range(d$y))
    }
    structure(list(density = x$density[i],
                   xlim = xlim,
                   ylim = ylim,
                   nobs = x$nobs[i],
                   kernel = x$kernel[i],
                   bw = x$bw[i],
                   sd = x$sd[i]),
              class="dtrace")
}

"[[.dtrace" <-
function(x, i)
{
    x$density[[i]]
}
